Hands-on 9

Code for download: session9_start.tar.gz

Start code: session8_start with updated geometry, more collections for the calorimeter, and added ScreenSD. Bounds for calorimeter histogram have been changed also. Prints for hit have been removed, in order to run statistics without being annoyed by long outputs.

Exercise 9a (physics)

1) Energy deposit depending on particle species**

Two calorimeter hits have been added in the hits collection to

  • collect deposit from charged particles in all layers
  • collect deposit from neutral particles

Add histograms for these new quantities, and fill these historams. In the run action, you will have to add the creation of these histograms as, for example:

// histogram 10 
analysisManager->CreateH1("AllCharged", "Charged Edep in all layers", 150, 0., 1500.);
// histogram 11 
analysisManager->CreateH1("AllNeutral", "Neutral Edep in all layers" , 100, 0., 100.);

Run 1000 protons in batch mode (much recommended !). You will see that energy deposit from neutral particle is quite marginal. Explain why neutral particles deposit so few compared to charged ones.

2) Energy deposit depending on cut:

Two other calorimeter hits have been added to

  • collect deposit from the primary track
  • collect deposit from other particles

Add histograms for these new quantities, and fill these historams. In the run action, you will have to add the creation of these histograms as, for example:

// histogram 12 
analysisManager->CreateH1("EdepPrimary", "Edep \[MeV\] by primary in calorimeter", 150, 0., 1500);
// histogram 13 
analysisManager->CreateH1("EdepOthers", "Edep \[MeV\] by non-primary in calorimeter", 150, 0., 1500);

Make two productions, running in batch mode 1000 protons with cuts of

  • 1 mm (everywhere)
  • 1 km (everywhere).

Compare the histogram for the primary and for the other particles. Explain what happens.

  • Hint: the calorimeter volume has its own cuts defined via the associated region. Find the appropriate command in the /run command directory.

The effect can be made more spectacular shooting EM particles, as most of the secondary production by EM particles is sensitive to the production cuts.

Prepare a runPositron.mac file in order to run in batch mode 1000 positrons (/gun/particle e+) of 1.2 GeV (so that they hit the calorimeter - verify interactively if needed). Run several jobs with cuts:

  • 1 mm
  • 1 cm
  • 10 cm
  • 1 m
  • 1 km

Compare again the histograms of energy deposit.

Compare also the histogram per layer. Observe that these histograms are not much different as long as the cut value is not large compared to the layer dimensions.

Exercise 9b

In the start code, a new volume has been added, a thin screen, far after the calorimeter. For this reason, the world volume has been made larger than in the previous exercices, as can be seen in the detector construction:

// Make world larger for this example:
hx += 5*m;
hy += 5*m;
hz += 5*m;

This screen will be used for counting particles that exit from the calorimeter. It will be made sensitive, but we will not create hits and hits collections, but just store in the ntuple the quantities we are interested in.

To this screen logical volume, attach a sensitive detector that you will have to complete. In the ProcessHits() of this sensitive detector, you will fill a Root ntuple:

// Store hit in the ntuple
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
analysisManager->FillNtupleIColumn(fNtupleId, 0, ID);
analysisManager->FillNtupleIColumn(fNtupleId, 1, pdgCode);
analysisManager->FillNtupleDColumn(fNtupleId, 2, Ekin/MeV);
analysisManager->FillNtupleDColumn(fNtupleId, 3, localPosition.x()/cm);
analysisManager->FillNtupleDColumn(fNtupleId, 4, localPosition.y()/cm);
analysisManager->FillNtupleDColumn(fNtupleId, 5, time/ns);

where ID is the track ID, pdgCode its PDG code, Ekin the kinetic energy, localPosition is the G4ThreeVector of local coordinate in the frame of the screen (see ChamberSD to see how to get these coordinates, for example), time is the hit time, and creatorModelID telling which model has created this track. Note that this model is internal to a physics process, and this ID is defined only in a few cases (too sad…). The creator model can be obtained as:

G4int creatorModelID = track->GetCreatorModelID(); 

The full list of model ID is printed at the beginning of the job (and indices change depending on the physics list). We will run quite statistics (10000 protons) several times and here the MT is quite useful. A useful Root class is the TChain:

root> TChain myChain("Screen");   
root> myChain.Add("ED.root");

then you can use this chain myChain as if it was a ntuple like Screen, eg:

root> myChain.Draw("Ekin");

to see Ekin for a given particle, eg. proton:

`root> myChain.Draw("Ekin","PdgCode==2212");

Change of neutron behavior with FTFP_BERT versus FTFP_BERT_HP:

When the job starts, there is a long print out of the physics list configuration. For the hadronics, this starts with the lines:

                  HADRONIC PROCESSES SUMMARY (verbose level 1)

Study the differences between these prints for neutrons when using FTFP_BERT and FTFP_BERT but with the high precision neutron physics. In particular, pay attention the differences for energies below 20 MeV.

Run 10k protons with FTFP_BERT and 10k protons FTFP_BERT_HP. Compare for neutrons (pdgCode == 2112) the difference in of the Ekin and time spectra of ntuple Screen.

NB : the long prints at the beginning of the job contain lots of information. As for HP versus non HP, you may look at the differences when using -say- FTFP_BERT versus FTFP_BERT_PEN, the print starting with:

phot:   for  gamma    SubType= 12  BuildTable= 0
      LambdaPrime table from 200 keV to 10 TeV in 154 bins

Note in particular the differences about Rayl, Deexcitation.

Exercise 9c (geometry persistency):

Implement the newly added command-line option [-g gdmlFile] in the main program: when the GDML file is provided export the geometry in GDML.

  • Hint: The world volume can be accessed in main in this way:
    G4LogicalVolume* world = G4TransportationManager::GetTransportationManager() ->  GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume());

Solution: session9_solution.tar.gz