Distinguish detectors for readout

BenjaminW · December 7, 2022, 9:24am

Hi,

I want to make use out of multiple detectors (in total something between 5 and 15) and at the moment I’m building up the geometry such that a detector is always placed as a daughter object within another object which I’m actually placing multiple times.

Now I wonder, as I basically have only one(?) detector instance, how can I distinguish these detectors in terms of readout resp. energy deposition within them?

My code looks like this atm:

// Outer transportation box
G4Box box_dects_out = new G4Box(“box_dects_out”, 353/2mm, 25/2mm, 250/2mm);
G4LogicalVolume *box_detectors_out = new G4LogicalVolume(box_dects_out, Pur, “box_detectors_out”);
G4VPhysicalVolume box_detectors_out_phys = new G4PVPlacement(0, G4ThreeVector(0.m, -304mm, 0.m), box_detectors_out, “box_detectors_out”, logicWorld, false, 1, true);
// Inner box
G4Box box_dects_in = new G4Box(“box_dects_in”, 350/2mm, 22/2mm, 247/2mm);
G4LogicalVolume *box_detectors_in = new G4LogicalVolume(box_dects_in, worldMat, “box_detectors_in”);
G4VPhysicalVolume *box_detectors_in_phys = new G4PVPlacement(0, G4ThreeVector(0.*m, 0.*mm, 0.*m), box_detectors_in, “box_detectors_in”, box_detectors_out, false, 1, true);

// Set up detector case
G4Box detector_case = new G4Box(“detector_case”, 52/2mm, 9/2mm, 22/2mm);
G4LogicalVolume *logDetector_case = new G4LogicalVolume(detector_case, worldMat, “logDetector_case”);
G4VPhysicalVolume *physDetector_case1 = new G4PVPlacement(0, G4ThreeVector(0.mm, -0mm, 0.*mm), logDetector_case, “logDetector_case1”, box_detectors_in, false, 1, true);
G4VPhysicalVolume *physDetector_case2 = new G4PVPlacement(0, G4ThreeVector(100.mm, -0mm, 0.*mm), logDetector_case, “logDetector_case2”, box_detectors_in, false, 1, true);
G4VPhysicalVolume *physDetector_case3 = new G4PVPlacement(0, G4ThreeVector(-100.mm, -0mm, 0.*mm), logDetector_case, “logDetector_case3”, box_detectors_in, false, 1, true);

// Set up outer 0.5 mm Al filter
G4Box *al_filter_out = new G4Box(“al_filter_out”, (4.7/2+0.5)*mm, (0.5/2+0.5)*mm, (4.7/2+0.5)*mm);
G4LogicalVolume *logAl_filter_out = new G4LogicalVolume(al_filter_out, Al, “logAl_filter_out”);
G4VPhysicalVolume *physAl_filter_out1 = new G4PVPlacement(0, G4ThreeVector(-15.mm, -0mm, 0.*mm), logAl_filter_out, “logAl_filter_out1”, logDetector_case, false, 1, true);
G4VPhysicalVolume *physAl_filter_out2 = new G4PVPlacement(0, G4ThreeVector(15.mm, -0mm, 0.*mm), logAl_filter_out, “logAl_filter_out2”, logDetector_case, false, 1, true);
// Set up inner 0.5 mm Al filter
G4Box *al_filter_in = new G4Box(“al_filter_in”, (4.7/2+0.5/2)*mm, (0.5/2+0.5/2)*mm, (4.7/2+0.5/2)*mm);
G4LogicalVolume *logAl_filter_in = new G4LogicalVolume(al_filter_in, worldMat, “logAl_filter_in”);
G4VPhysicalVolume *physAl_filter_in = new G4PVPlacement(0, G4ThreeVector(0.mm, -0mm, 0.*mm), logAl_filter_in, “logAl_filter_in”, logAl_filter_out, false, 1, true);

// Set up detectors
G4Box solidDetector = new G4Box(“solidDetector”, 4.7/2mm, 0.5/2mm, 4.7/2mm);
logicDetector = new G4LogicalVolume(solidDetector, BeO, “logicDetector”);
G4VPhysicalVolume *physDetector1 = new G4PVPlacement(0, G4ThreeVector(0.mm, 0mm, 0.*mm), logicDetector, “physDetector”, logAl_filter_in, false, 1, true);

Is this a meaningful way to go or should I do it differently?

I believe this is discussed here: Identification of unique physical volumes with IDs ,yet, I’m not quite sure that I understand it correctly. It seems I need to get the data out of G4StepPoint with a pointer to the physical volumes. Correct? But how exactly?

mkelsey · December 7, 2022, 7:10pm

When you do G4PVPlacement(), assign a different copy number to each instance (look at $G4INCLUDE/G4PVPlacement.hh to see which constructor argument is the copy number). In your SteppingAction, you’ll have access to the G4Touchable for the actual detector that was hit, and you can query that to get the different copy numbers corresponding to the different mother volume placements.

BenjaminW · December 8, 2022, 11:54am

Thank you, but as far as I understand it correctly, I’m using G4PVPlacement() only once. But as being a daughter object it will be placed multiple times nevertheless.
How can I give all these detectors a different copy number when I use the command only once?

Sorry if this sounds silly or if I don’t understand it, I’m trying to grasp the (detector/readout) concept in general and also the way to go.

mkelsey · December 8, 2022, 3:25pm

You must be using G4PVPlacement multiple times, as you described yourself,

When you put the detector in its mother, that’s only one placement. You probably give that a copy number of 0 or 1 (that’s what we do in our experiment).
But then you have a whole separate placement of the mother. Each of those placements should be given its own unique copy number (1, 2, 3, 4, etc.).

In your SteppingAction, you can use the G4Touchable for the post step point to get both the placement of the detector, but also the placement of the mother volume. That’s where you can find the unique identifier for which of your real detectors was hit.

BenjaminW · December 8, 2022, 3:27pm

Omg, yeah, sure, now I (finally) got it! Thanks a lot, that makes totally sense
I’ll try to get this working, thank you!

Ok, I guess it is somehow working yet I have to figure out how to use the information resp. which I have to use exactly.I noticed that only

const G4VTouchable *post_touchable = aStep->GetPostStepPoint()->GetTouchable();
detector_num = post_touchable->GetCopyNumber(2);
G4cout << "Detector number (post)(2): " << detector_num << G4endl;

shows varying numbers where other commands like for GetPreStepPoint() or even GetCopyNumber(1), GetCopyNumber(0) for GetPostStepPoint() do not show any different numbers.
Also

G4int motherCopyNo = touchable->GetReplicaNumber(1);
G4cout << "Mother volume (pre): " << motherCopyNo << G4endl;

prints out 0 permanently.

G4int motherCopyNo = touchable->GetCopyNumber(1);
G4cout << "Mother volume (pre): " << motherCopyNo << G4endl;		

G4int copyNo = touchable->GetCopyNumber();
G4cout << "copyNo: " << copyNo << G4endl;

is also always 0.
It could also be everything valid (I guess it is), just a bit difficult to figure out whether I’m using the right commands when it’s not quite clear what they should or could print
I guess I understand the concept and so on but it is not so easy with only getting numbers of {0,…,4} to determine what’s going on. Especially when it’s mostly or only 0.

Which information do I need resp. how do I exactly access the information I need?
Is there also a way to get a char/label or the like (like physDetector1, physDetector2,…) instead of only an integer?

loydms · December 8, 2022, 8:01pm

When you give a parameter to GetReplicaNumber(), you are specifying which level of parent to reach (I believe, mkselsey knows more and has answered this in this link).

If you don’t give it any arguments, it will give you the copy number.

So if you have a single mother object with say 10 detectors inside of it, then calling GetReplicaNumber() will give you the detector number (0-9). But calling GetReplicaNumber(1) will give you the number of the mother object, which will always be the same.

I have done parameterization, and GetRplicaNumber() without arguments returns a unique parameter index

mkelsey · December 8, 2022, 9:52pm

Right, @loydms ! What you describe is correct, when you have made multiple placements of the “lowest level” objects. What @BenjaminW has is slightly more complex.

Create one LV “detector”. Place that one LV inside a single mother volume “holder”. There’s only one “detector” PV, so it will always have the same copy number, let’s say 0.
Now take the mother “holder” volume, and do 10 placements of it somewhere in your geometry. You give each of those “holder” PV’s a unique copy number, let’s say 1, 2, 3, … 9, 10.

Run an event, you get it hit in one of those detectors. When you call touchable->GetReplicaNumber() you will always get back 0, because that’s the lowest level “detector” object. But if you call touchable->GetReplicaNumber(1), in this case (which is different from your case), you will get one of the unique “holder” PVs, and get back 1 or 2 or 7 or something.

BenjaminW · December 9, 2022, 8:21am

Thanks a lot for your input/help! I guess I got it, probably only have to play a bit around to determine the correct “levels/instances” of these multiple mother-relationships I have within my code.

Can I suggest improvements or the like for G4 somewhere? I believe changing the data type from like integer to a string shouldn’t be difficult. It would probably help a lot to identify and ensure all the volumes
My geometry is already rather simple I would say, won’t imagine the real large geometries of the large detector facilities how they do that?

edit: I can assign arbitrary numbers for CopyNo, right? So I could use e.g. 100, 101, 102,… for one volume level, 200, 201, 202 for the daughter volume or the like?
Or is there maybe a better approach?

Btw, how do I select the hits which have hit only the detector? At the moment I track all particles and they pass through very different volumes and thus also the CopyNo is not constant.

maire · December 9, 2022, 9:47am

Probably, it will not help you; but, in case :
Geant4 → User support → FAQ -->Tracks and steps

BenjaminW · December 9, 2022, 10:11am

thank you, I’m reading here and there actually quite a lot but my progress is rather slow as I have to put all the pieces together. When one understands the situation it is quite easy to use a manual but when learning reading a manual is much harder though all information is there

However, regarding my question whether I can select only the particles which hit a detector: As far as I understand the logical detector is declared as a (sensitive) detector:

MySensitiveDetector *sensDet = new MySensitiveDetector(SensitiveDetector);
logicDetector->SetSensitiveDetector(sensDet);

so I’d assume I’m tracking already only particles which have hit the detector but I’m definitely tracking a lot of particles which don’t hit the detector. The class MySensitiveDetector starts as follows:

G4Track *track = aStep->GetTrack();
track->SetTrackStatus(fStopAndKill);			

G4StepPoint *preStepPoint = aStep->GetPreStepPoint();		
G4ThreeVector posPhoton = preStepPoint->GetPosition();	
const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();

does this help someone? Do I have to change something or so? I’m also fine with an example

weller · December 9, 2022, 11:41am

you should only track particles that are inside the logical volume which you defined as the sensitive detector. the correct location for that would be in

G4bool mySensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory* touchables) {
...

This is more or less the UserSteppingAction that is only valid in the sensitive detector’s volume.

edit: example

G4bool fluxSD::ProcessHits(G4Step* aStep, G4TouchableHistory* /*touchables*/)
{
    G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
    if (preStepPoint->GetStepStatus() == fGeomBoundary ) {     // inbound particles only
        if(aStep->GetTrack()->GetParticleDefinition() == G4Gamma::Gamma()) {
            auto hit = new fluxHit();
            hit->setEnergy(preStepPoint->GetKineticEnergy());
            hit->setWeight(preStepPoint->GetWeight());
            hit->setCopyNumber(preStepPoint->GetTouchableHandle()->GetCopyNumber(fDepth));
            fHitsCollection->insert(hit);
//            hit->Print();
        }
    }

  return true;
}

BenjaminW · December 9, 2022, 1:39pm

Thank you, I’m just not sure how to use hit and fHitsCollection and how to proceed from there on (e.g. storing deposited energy in a specific detector in a histogram?).

What I have is this:

#include “detector.hh”
MySensitiveDetector::MySensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{}

MySensitiveDetector::~MySensitiveDetector()
{}

G4bool MySensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist)
{
G4Track *track = aStep->GetTrack();
track->SetTrackStatus(fStopAndKill);

G4StepPoint *preStepPoint = aStep->GetPreStepPoint(); // particle enters detector
G4StepPoint *postStepPoint = aStep->GetPostStepPoint(); // leaves detector
G4ThreeVector posPhoton = preStepPoint->GetPosition();

const G4VTouchable *post_touchable = aStep->GetPostStepPoint()->GetTouchable();
const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();

G4VPhysicalVolume *physVol = touchable->GetVolume();
G4cout << "Volume: " << physVol << G4endl;
G4ThreeVector posDetector = physVol->GetTranslation();

G4int evt = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
G4double partEdep = aStep->GetTotalEnergyDeposit();
G4double energy = preStepPoint->GetTotalEnergy();

//Convert from MeV to keV
partEdep *= 1000;
energy *= 1000;

// Create histograms
G4AnalysisManager *man = G4AnalysisManager::Instance();
man->FillNtupleIColumn(0, evt);
man->FillNtupleDColumn(1, partEdep);
man->FillNtupleDColumn(2, energy);
man->FillNtupleDColumn(3, posDetector[0]);
man->FillNtupleDColumn(4, posDetector[1]);
man->FillNtupleDColumn(5, posDetector[2]);
man->AddNtupleRow(0);

return true;
}

I’m working on posDetector[]… from my basic knowledge I don’t see why I am also receiving information of particles which didn’t hit the detector.
Any idea how I can get the information/hit of a specific detector from here on?
I guess what I could do is something like this:

G4Box *detector_case = new G4Box("detector_case", 52/2*mm, 9/2*mm, 22/2*mm);
G4LogicalVolume *logDetector_case = new G4LogicalVolume(detector_case, worldMat, "logDetector_case");
G4VPhysicalVolume *physDetector_case1 = new G4PVPlacement(0, G4ThreeVector(0.*mm, 0*mm, 0.*mm), logDetector_case, "logDetector_case1", box_detectors_in, false, 200, true);	
G4VPhysicalVolume *physDetector_case2 = new G4PVPlacement(0, G4ThreeVector(100.*mm, 0*mm, 0.*mm), logDetector_case, "logDetector_case2", box_detectors_in, false, 210, true);	
G4VPhysicalVolume *physDetector_case3 = new G4PVPlacement(0, G4ThreeVector(-100.*mm, 0*mm, 0.*mm), logDetector_case, "logDetector_case3", box_detectors_in, false, 220, true);

and then demanding somewhere in G4bool MySensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist):

const G4VTouchable *post_touchable = aStep->GetPostStepPoint()->GetTouchable();
G4int detector_num0 = post_touchable->GetReplicaNumber(0);
G4int detector_num1 = post_touchable->GetReplicaNumber(1);
G4int detector_num2 = post_touchable->GetReplicaNumber(2);
…

if((detector_num1 == 200) || (detector_num1 == 200) || …)
{
…store in physDetector1_histo
}

if(detector_num1 == 210) || (detector_num1 == 210) || …)
{
…store in physDetector2_histo
}
…

but this doesn’t look like good Geant4 code… if there were 10,000 detectors I doubt anyone would write 10,000 if clauses ?

loydms · December 9, 2022, 1:47pm

Ah I figured it was something more complicated but I wasn’t quite understanding how he had done his parameterization.

You learn something new everyday

weller · December 9, 2022, 3:47pm

not quite sure what you want to achieve, but you could define the copy or replica number in such a way that the copy/replica number is the histogram ID that you want to fill.
also look into 2d/3d histograms, if you want to have multiple detectors (1 dimension of the histogram) in multiple mother volumes (another dimension of the histogram), and want to store e.g. the energy (third dimension)

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/managers.html#creating-histograms

mkelsey · December 9, 2022, 7:25pm

@BenjaminW , what @weller wrote is the right approach if you’re filling histograms directly. You can make your storage really compact by using a TH2D (or whatever the G4Analysis equivalent is), putting the detector ID (copy number) on one axis and energy on the other. If you do that, you probably want the IDs to be sequential integers from 1 to N, and document how they relate to your geometry.

You could also use arbitrary detector IDs (we use a digit mapping so the top detector in the fourth tower of the second column is ID “241”), and store separate 1D histograms in a std::map<G4int, TH1D*>.

Either way, you don’t need a massive chain of if statements, just use the detector ID as an index/argument.

BenjaminW · December 12, 2022, 10:33am

Thank you! Would this mean, when I e.g. assign CopyNo = 210 to physDetector1 that for G4int CreateH1() I would use 210 as one of the parameters (which would this be?)?

I’m not sure about the 2D/3D histograms, in my case I’m placing the detectors rather arbitrary and they don’t represent any meaningful/larger geometry (like a big circular detector or the like), so I guess this doesn’t add information?

BenjaminW · December 12, 2022, 10:36am

Thank you as well! Guess I’ll go with the second advice.
I think I got the idea/concept, just have to figure out what this means as code

weller · December 12, 2022, 10:59am

i think the histograms are incrementally numbered starting from either 0 or a value that you define:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Analysis/managers.html?highlight=setfirsthistoid

The start histogram identifier value can be changed either with the SetFirstHistoId(G4int) method, which applies the new value to all histogram types, or with the SetFirstHNId(G4int) , where N = 1, 2, 3 methods, which apply the new value only to the relevant histogram type.

where the id is not a parameter that can be set. @mkelsey circumvents this with the digit mapping.

you would then follow this numbering with your copynumber, and not the other way around.

alternative: write a mapping from copynumber (could be "1,2,3,… " or “(0,1), (0,2)…(N,M)” or more dimensions as required to histoID = 0,1,2,…, N*M

BenjaminW · December 13, 2022, 7:40am

Thanks a lot, again! Though I didn’t totally grasp this, yet, I’ll try to realize it.
Will probably be back in a few days

BenjaminW · January 12, 2023, 2:19pm

Finally an update from my side. For a longer break and reading everything up again, I’d say I understood all your input (finally).
I just realize that I’m facing another issue which troubled me last year but couldn’t realize in which way.

So, the situation is the following:

I have detectors with a different amount of filters around them. Actually, one detector has e.g. 2 filters around and the other only one filter around. This is given for each LV.
Hence, when I try to identify the mother volumes of the detector via GetCopyNumber() like GetCopyNumber(5), in one case I’ll get the mother instance of the detector indeed. This is the case when there is only one filter around the detector.
But when the detector has got two filters around, with GetCopyNumber(5) I’ll get the second filter and with GetCopyNumber(7) I’ll get the level I need to identify the detector.

Any idea how I could handle this in a smart way?

Is there a command to receive the amount of available levels resp. to which number GetCopyNumber() can go?