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.
Two calorimeter hits have been added in the hits collection to
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.
Two other calorimeter hits have been added to
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
Compare the histogram for the primary and for the other particles. Explain what happens.
/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:
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.
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);
analysisManager->AddNtupleRow(fNtupleId);
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");
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.