How to find out the breakdown of hit energies in Monte Carlo into particle types.
ECAL example
These notes describe what I did for the ECAL (as presented at the Lyon meeting). The situation will be more complicated for the HCAL because of the need to gang cells together. I have not yet summoned up the strength to tackle this problem.
- First, when you run Mokka you need to make sure include the following in your steering file:
/Mokka/init/lcioDetailedShowerMode true
This will increase the size of the output files significantly (a factor of 2-3 maybe).
- If you use digitisation code, you may need to make some changes. I was using my own SimpleEcalDigi processor for the ECAL. This was converting SimCalorimeterHits into CalorimeterHits. But in order to pass the true hit composition information on to subsequent analysis processors, I had to change it to write out SimCalorimeterHits. The modified versions are here:
SimpleEcalDigi.cc ;
SimpleEcalDigi.h.
When this processor adds noise in to the hits, it ascribes it a PDG code 0.
- To access the hit composition information in your analysis job, use the getEnergyCont method of SimCalorimeterHit. Something like this:
.....
Ehit=hit->getEnergy();
int NC=hit->getNMCContributions();
// Ee should give the energy associated with e+ or e-
// Emuon the energy associated with muons, pions or K+/K-
// Ep the energy associated with (anti)protons
// Eother the rest (noise, hyperons, deuterons, alpha-particles etc.)
Ee=0; Emeson=0; Ep=0; Eother=0;
for( int ii=0; ii<NC ; ii++) {
if(abs(hit->getPDGCont(ii))==11) Ee+=hit->getEnergyCont(ii);
else if(abs(hit->getPDGCont(ii))==13 || abs(hit->getPDGCont(ii))==211
|| abs(hit->getPDGCont(ii))==321) Emeson+=hit->getEnergyCont(ii);
else if(abs(hit->getPDGCont(ii))==2212) Ep+=hit->getEnergyCont(ii);
else {Eother+=hit->getEnergyCont(ii);}
}
.....
You can obviously customise this as you please to look at separate groups of PDG codes.
- I hope other users will add to this page if they try to use this information.
--
DavidWard - 22-Oct-2009
AHCAL example
Here you can find an example processor for the AHCAL case. The processor creates ROOT file with histograms of the energy sums for the different components, as well as of the transverse profile. Note that for this you need simulated hits from the drift chambers, in order to fit the track (i.e. shower axis).
When testing it on a Mokka file, don't forget to create the Mokka file using
/Mokka/init/lcioDetailedShowerMode true
as David said, and don't generate more than 10000 events, since the file size increases fast (and 10000 events is more than enough for most purposes).
AnaGen.hh ,
AnaGen.cc
In case you discover problems/bugs, or you have any suggestions of improvements, please let me know, and feel free to make the corresponding changes.
--
AngelaLucaciTimoce - 11-Jan-2010
- Transverse profile for 15 GeV pion, with MC components: