test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
StackingAction.cc
Go to the documentation of this file.
6 
8 
9 #include "G4VProcess.hh"
10 #include "G4EmProcessSubType.hh"
11 #include "G4LogicalVolumeStore.hh"
12 #include "G4RegionStore.hh"
13 #include "Randomize.hh"
14 #include "G4SystemOfUnits.hh"
15 #include "G4VSolid.hh"
16 #include "G4TransportationManager.hh"
17 
19  const CMSSteppingVerbose* sv)
20  : trackAction(trka),steppingVerbose(sv)
21 {
22  trackNeutrino = p.getParameter<bool>("TrackNeutrino");
23  killHeavy = p.getParameter<bool>("KillHeavy");
24  killGamma = p.getParameter<bool>("KillGamma");
25  kmaxGamma = p.getParameter<double>("GammaThreshold")*MeV;
26  kmaxIon = p.getParameter<double>("IonThreshold")*MeV;
27  kmaxProton = p.getParameter<double>("ProtonThreshold")*MeV;
28  kmaxNeutron = p.getParameter<double>("NeutronThreshold")*MeV;
29  killDeltaRay = p.getParameter<bool>("KillDeltaRay");
30  limitEnergyForVacuum = p.getParameter<double>("CriticalEnergyForVacuum")*MeV;
31  maxTrackTime = p.getParameter<double>("MaxTrackTime")*ns;
32  maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
33  maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
34  deadRegionNames= p.getParameter<std::vector<std::string> >("DeadRegions");
36  p.getUntrackedParameter<bool>("SaveAllPrimaryDecayProductsAndConversions",
37  true);
39  p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInTracker",
40  false);
42  p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo",
43  false);
45  p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon",
46  false);
48  p.getUntrackedParameter<bool>("SaveFirstLevelSecondary",false);
49  killInCalo = false;
50  killInCaloEfH = false;
51 
52  // Russian Roulette
53  regionEcal = 0;
54  regionHcal = 0;
55  regionMuonIron = 0;
56  regionPreShower= 0;
57  regionCastor = 0;
58  regionWorld = 0;
59 
60  gRusRoEnerLim = p.getParameter<double>("RusRoGammaEnergyLimit")*MeV;
61  nRusRoEnerLim = p.getParameter<double>("RusRoNeutronEnergyLimit")*MeV;
62 
63  gRusRoEcal = p.getParameter<double>("RusRoEcalGamma");
64  gRusRoHcal = p.getParameter<double>("RusRoHcalGamma");
65  gRusRoMuonIron = p.getParameter<double>("RusRoMuonIronGamma");
66  gRusRoPreShower = p.getParameter<double>("RusRoPreShowerGamma");
67  gRusRoCastor = p.getParameter<double>("RusRoCastorGamma");
68  gRusRoWorld = p.getParameter<double>("RusRoWorldGamma");
69 
70  nRusRoEcal = p.getParameter<double>("RusRoEcalNeutron");
71  nRusRoHcal = p.getParameter<double>("RusRoHcalNeutron");
72  nRusRoMuonIron = p.getParameter<double>("RusRoMuonIronNeutron");
73  nRusRoPreShower = p.getParameter<double>("RusRoPreShowerNeutron");
74  nRusRoCastor = p.getParameter<double>("RusRoCastorNeutron");
75  nRusRoWorld = p.getParameter<double>("RusRoWorldNeutron");
76 
77  gRRactive = false;
78  nRRactive = false;
79  if(gRusRoEnerLim > 0.0 &&
80  (gRusRoEcal < 1.0 || gRusRoHcal < 1.0 ||
81  gRusRoMuonIron < 1.0 || gRusRoPreShower < 1.0 || gRusRoCastor < 1.0 ||
82  gRusRoWorld < 1.0)) { gRRactive = true; }
83  if(nRusRoEnerLim > 0.0 &&
84  (nRusRoEcal < 1.0 || nRusRoHcal < 1.0 ||
85  nRusRoMuonIron < 1.0 || nRusRoPreShower < 1.0 || nRusRoCastor < 1.0 ||
86  nRusRoWorld < 1.0)) { nRRactive = true; }
87 
88  if ( p.exists("TestKillingOptions") ) {
89  killInCalo = (p.getParameter<edm::ParameterSet>("TestKillingOptions"))
90  .getParameter<bool>("KillInCalo");
91  killInCaloEfH = (p.getParameter<edm::ParameterSet>("TestKillingOptions"))
92  .getParameter<bool>("KillInCaloEfH");
93  edm::LogWarning("SimG4CoreApplication")
94  << " *** Activating special test killing options in StackingAction \n"
95  << " *** Kill secondaries in Calorimetetrs volume = " << killInCalo << "\n"
96  << " *** Kill electromagnetic secondaries from hadrons in Calorimeters volume= "
97  << killInCaloEfH;
98 
99  }
100 
101  initPointer();
102  newTA = new NewTrackAction();
103 
104  edm::LogInfo("SimG4CoreApplication") << "StackingAction initiated with"
105  << " flag for saving decay products in "
106  << " Tracker: " << savePDandCinTracker
107  << " in Calo: " << savePDandCinCalo
108  << " in Muon: " << savePDandCinMuon
109  << " everywhere: " << savePDandCinAll
110  << "\n saveFirstSecondary"
111  << ": " << saveFirstSecondary
112  << " Tracking neutrino flag: "
113  << trackNeutrino
114  << " Kill Delta Ray flag: "
115  << killDeltaRay
116  << " Kill hadrons/ions flag: "
117  << killHeavy;
118 
119 
120  if(killHeavy) {
121  edm::LogInfo("SimG4CoreApplication") << "StackingAction kill protons below "
122  << kmaxProton/MeV <<" MeV, neutrons below "
123  << kmaxNeutron/MeV << " MeV and ions"
124  << " below " << kmaxIon/MeV << " MeV";
125  }
126  killExtra = killDeltaRay || killHeavy || killInCalo || killInCaloEfH;
127 
128  edm::LogInfo("SimG4CoreApplication") << "StackingAction kill tracks with "
129  << "time larger than " << maxTrackTime/ns
130  << " ns ";
131  numberTimes = maxTimeNames.size();
132  if(0 < numberTimes) {
133  for (unsigned int i=0; i<numberTimes; ++i) {
134  edm::LogInfo("SimG4CoreApplication") << "StackingAction MaxTrackTime for "
135  << maxTimeNames[i] << " is "
136  << maxTrackTimes[i] << " ns ";
137  maxTrackTimes[i] *= ns;
138  }
139  }
140  if(limitEnergyForVacuum > 0.0) {
141  edm::LogInfo("SimG4CoreApplication")
142  << "StackingAction LowDensity regions - kill if E < "
143  << limitEnergyForVacuum/MeV << " MeV";
144  printRegions(lowdensRegions,"LowDensity");
145  }
146  if(deadRegions.size() > 0.0) {
147  edm::LogInfo("SimG4CoreApplication")
148  << "StackingAction Dead regions - kill all secondaries ";
149  printRegions(deadRegions, "Dead");
150  }
151  if(gRRactive) {
152  edm::LogInfo("SimG4CoreApplication")
153  << "StackingAction: "
154  << "Russian Roulette for gamma Elimit(MeV)= "
155  << gRusRoEnerLim/MeV << "\n"
156  << " ECAL Prob= " << gRusRoEcal << "\n"
157  << " HCAL Prob= " << gRusRoHcal << "\n"
158  << " MuonIron Prob= " << gRusRoMuonIron << "\n"
159  << " PreShower Prob= " << gRusRoPreShower << "\n"
160  << " CASTOR Prob= " << gRusRoCastor << "\n"
161  << " World Prob= " << gRusRoWorld;
162  }
163  if(nRRactive) {
164  edm::LogInfo("SimG4CoreApplication")
165  << "StackingAction: "
166  << "Russian Roulette for neutron Elimit(MeV)= "
167  << nRusRoEnerLim/MeV << "\n"
168  << " ECAL Prob= " << nRusRoEcal << "\n"
169  << " HCAL Prob= " << nRusRoHcal << "\n"
170  << " MuonIron Prob= " << nRusRoMuonIron << "\n"
171  << " PreShower Prob= " << nRusRoPreShower << "\n"
172  << " CASTOR Prob= " << nRusRoCastor << "\n"
173  << " World Prob= " << nRusRoWorld;
174  }
175 
176  if(savePDandCinTracker) {
177  edm::LogInfo("SimG4CoreApplication") << "StackingAction Tracker regions: ";
178  printRegions(trackerRegions,"Tracker");
179  }
180  if(savePDandCinCalo) {
181  edm::LogInfo("SimG4CoreApplication") << "StackingAction Calo regions: ";
182  printRegions(caloRegions, "Calo");
183  }
184  if(savePDandCinMuon) {
185  edm::LogInfo("SimG4CoreApplication") << "StackingAction Muon regions: ";
186  printRegions(muonRegions,"Muon");
187  }
188  worldSolid = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume()->GetSolid();
189 }
190 
192 {
193  delete newTA;
194 }
195 
196 G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track * aTrack)
197 {
198  // G4 interface part
199  G4ClassificationOfNewTrack classification = fUrgent;
200  int pdg = aTrack->GetDefinition()->GetPDGEncoding();
201  int abspdg = std::abs(pdg);
202 
203  // primary
204  if (aTrack->GetCreatorProcess()==nullptr || aTrack->GetParentID()==0) {
205  if (!trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18)) {
206  classification = fKill;
207  } else if (worldSolid->Inside(aTrack->GetPosition()) == kOutside) {
208  classification = fKill;
209  } else {
210  newTA->primary(aTrack);
211  }
212  } else {
213 
214  // secondary
215  const G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
216  // definetly killed tracks
217  if (aTrack->GetTrackStatus() == fStopAndKill) { classification = fKill; }
218  else if (!trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18))
219  { classification = fKill; }
220  else if (isItOutOfTimeWindow(reg, aTrack)) { classification = fKill; }
221 
222  // potentially good for tracking
223  else {
224  double ke = aTrack->GetKineticEnergy();
225 
226  // kill tracks in specific regions
227  if (classification != fKill && isThisRegion(reg, deadRegions)) {
228  classification = fKill;
229  }
230  if (ke <= limitEnergyForVacuum && isThisRegion(reg, lowdensRegions) ) {
231  classification = fKill;
232 
233  } else {
234  // very low-energy gamma
235  if (pdg == 22 && killGamma && ke < kmaxGamma) { classification = fKill; }
236 
237  // specific track killing - not for production
238  if(killExtra && classification != fKill) {
239  if (killHeavy && classification != fKill) {
240  if (((pdg/1000000000 == 1) && (((pdg/10000)%100) > 0) &&
241  (((pdg/10)%100) > 0) && (ke<kmaxIon)) ||
242  ((pdg == 2212) && (ke < kmaxProton)) ||
243  ((pdg == 2112) && (ke < kmaxNeutron))) { classification = fKill; }
244  }
245 
246  if (killDeltaRay && classification != fKill
247  && aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation) {
248  classification = fKill;
249  }
250  if (killInCalo && classification != fKill
251  && isThisRegion(reg, caloRegions)) {
252  classification = fKill;
253  }
254  if (killInCaloEfH && classification != fKill) {
255  int pdgMother = std::abs(trackAction->geant4Track()->GetDefinition()->GetPDGEncoding());
256  if ((pdg == 22 || abspdg == 11) && pdgMother != 11 && pdgMother != 22 &&
257  isThisRegion(reg,caloRegions)) {
258  classification = fKill;
259  }
260  }
261  }
262 
263  // Russian roulette && MC truth
264  if(classification != fKill) {
265  const G4Track* mother = trackAction->geant4Track();
266  int flag = 0;
267  if(savePDandCinAll) {
268  flag = isItPrimaryDecayProductOrConversion(aTrack, *mother);
269  } else {
273  flag = isItPrimaryDecayProductOrConversion(aTrack, *mother);
274  }
275  }
276  if (saveFirstSecondary && 0 == flag) {
277  flag = isItFromPrimary(*mother, flag);
278  }
279 
280  // Russian roulette
281  if(2112 == pdg || 22 == pdg) {
282  double currentWeight = aTrack->GetWeight();
283 
284  if(1.0 >= currentWeight) {
285  double prob = 1.0;
286  double elim = 0.0;
287 
288  // neutron
289  if(nRRactive && pdg == 2112) {
290  elim = nRusRoEnerLim;
291  if(reg == regionEcal) { prob = nRusRoEcal; }
292  else if(reg == regionHcal) { prob = nRusRoHcal; }
293  else if(reg == regionMuonIron) { prob = nRusRoMuonIron; }
294  else if(reg == regionPreShower) { prob = nRusRoPreShower; }
295  else if(reg == regionCastor) { prob = nRusRoCastor; }
296  else if(reg == regionWorld) { prob = nRusRoWorld; }
297 
298  // gamma
299  } else if(gRRactive && pdg == 22) {
300  elim = gRusRoEnerLim;
301  if(reg == regionEcal || reg == regionPreShower) {
302  if(rrApplicable(aTrack, *mother)) {
303  if(reg == regionEcal) { prob = gRusRoEcal; }
304  else { prob = gRusRoPreShower; }
305  }
306  } else {
307  if(reg == regionHcal) { prob = gRusRoHcal; }
308  else if(reg == regionMuonIron) { prob = gRusRoMuonIron; }
309  else if(reg == regionCastor) { prob = gRusRoCastor; }
310  else if(reg == regionWorld) { prob = gRusRoWorld; }
311  }
312  }
313  if(prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
314  if(G4UniformRand() < prob) {
315  const_cast<G4Track*>(aTrack)->SetWeight(currentWeight/prob);
316  } else {
317  classification = fKill;
318  }
319  }
320  }
321  }
322  if (classification != fKill) {
323  newTA->secondary(aTrack, *mother, flag);
324  }
325  LogDebug("SimG4CoreApplication") << "StackingAction:Classify Track "
326  << aTrack->GetTrackID() << " Parent "
327  << aTrack->GetParentID() << " Type "
328  << aTrack->GetDefinition()->GetParticleName()
329  << " K.E. " << aTrack->GetKineticEnergy()/MeV
330  << " MeV from process/subprocess "
331  << aTrack->GetCreatorProcess()->GetProcessType()
332  << "|"
333  <<aTrack->GetCreatorProcess()->GetProcessSubType()
334  << " as " << classification << " Flag " << flag;
335  }
336  }
337  }
338  }
339  if(nullptr != steppingVerbose) {
340  steppingVerbose->StackFilled(aTrack, (classification == fKill));
341  }
342  return classification;
343 }
344 
346 
348 
350 {
351  // Russian roulette
352  const G4RegionStore * rs = G4RegionStore::GetInstance();
353  std::vector<G4Region*>::const_iterator rcite;
354  std::vector<std::string>::const_iterator rd;
355 
356  for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
357  const G4Region* reg = (*rcite);
358  G4String rname = reg->GetName();
359  if ((gRusRoEcal < 1.0 || nRusRoEcal < 1.0) && rname == "EcalRegion") {
360  regionEcal = reg;
361  }
362  if ((gRusRoHcal < 1.0 || nRusRoHcal < 1.0) && rname == "HcalRegion") {
363  regionHcal = reg;
364  }
365  if ((gRusRoMuonIron < 1.0 || nRusRoMuonIron < 1.0) && rname == "MuonIron") {
366  regionMuonIron = reg;
367  }
368  if ((gRusRoPreShower < 1.0 || nRusRoPreShower < 1.0) && rname == "PreshowerRegion") {
369  regionPreShower = reg;
370  }
371  if ((gRusRoCastor < 1.0 || nRusRoCastor < 1.0) && rname == "CastorRegion") {
372  regionCastor = reg;
373  }
374  if ((nRusRoWorld < 1.0 || nRusRoWorld < 1.0) && rname == "DefaultRegionForTheWorld") {
375  regionWorld = reg;
376  }
377 
378  // time limits
379  unsigned int num = maxTimeNames.size();
380  maxTimeRegions.resize(num, 0);
381  if (num > 0) {
382  for (unsigned int i=0; i<num; i++) {
383  if (rname == (G4String)(maxTimeNames[i])) {
384  maxTimeRegions[i] = reg;
385  break;
386  }
387  }
388  }
389  //
390  if(savePDandCinTracker &&
391  (rname == "BeamPipe" || rname == "BeamPipeVacuum"
392  || rname == "TrackerPixelSensRegion"
393  || rname == "TrackerPixelDeadRegion"
394  || rname == "TrackerDeadRegion" || rname == "TrackerSensRegion")) {
395  trackerRegions.push_back(reg);
396  }
397  if(savePDandCinCalo &&
398  (rname == "HcalRegion" || rname == "EcalRegion"
399  || rname == "PreshowerSensRegion" || rname == "PreshowerRegion")) {
400  caloRegions.push_back(reg);
401  }
402  if(savePDandCinMuon &&
403  (rname == "MuonChamber" || rname == "MuonSensitive_RPC"
404  || rname == "MuonIron" || rname == "Muon"
405  || rname == "MuonSensitive_DT-CSC") ) {
406  muonRegions.push_back(reg);
407  }
408  if(rname == "BeamPipeOutside" || rname == "BeamPipeVacuum") {
409  lowdensRegions.push_back(reg);
410  }
411  for (rd = deadRegionNames.begin(); rd != deadRegionNames.end(); ++rd) {
412  if(rname == (G4String)(*rd)) {
413  deadRegions.push_back(reg);
414  }
415  }
416  }
417 }
418 
419 bool StackingAction::isThisRegion(const G4Region* reg,
420  std::vector<const G4Region*>& regions) const
421 {
422  bool flag = false;
423  unsigned int nRegions = regions.size();
424  for(unsigned int i=0; i<nRegions; ++i) {
425  if(reg == regions[i]) {
426  flag = true;
427  break;
428  }
429  }
430  return flag;
431 }
432 
434  const G4Track & mother) const
435 {
436  int flag = 0;
437  const TrackInformation & motherInfo(extractor(mother));
438  // Check whether mother is a primary
439  if (motherInfo.isPrimary()) {
440  if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) {
441  flag = 1;
442  } else if (aTrack->GetCreatorProcess()->GetProcessSubType()==fGammaConversion) {
443  flag = 2;
444  }
445  }
446  return flag;
447 }
448 
449 bool StackingAction::rrApplicable(const G4Track * aTrack,
450  const G4Track & mother) const
451 {
452  const TrackInformation & motherInfo(extractor(mother));
453 
454  // Check whether mother is gamma, e+, e-
455  int genID = motherInfo.genParticlePID();
456  return (22 == genID || 11 == genID || -11 == genID) ? false : true;
457 }
458 
459 int StackingAction::isItFromPrimary(const G4Track & mother, int flagIn) const
460 {
461  int flag = flagIn;
462  if (flag != 1) {
463  const TrackInformation & motherInfo(extractor(mother));
464  if (motherInfo.isPrimary()) { flag = 3; }
465  }
466  return flag;
467 }
468 
469  bool StackingAction::isItOutOfTimeWindow(const G4Region* reg, const G4Track* aTrack) const
470 {
471  double tofM = maxTrackTime;
472  for (unsigned int i=0; i<numberTimes; ++i) {
473  if (reg == maxTimeRegions[i]) {
474  tofM = maxTrackTimes[i];
475  break;
476  }
477  }
478  return (aTrack->GetGlobalTime() > tofM) ? true : false;
479 }
480 
481 void StackingAction::printRegions(const std::vector<const G4Region*>& reg,
482  const std::string& word) const
483 {
484  for (unsigned int i=0; i<reg.size(); ++i) {
485  edm::LogInfo("SimG4CoreApplication") << " StackingAction: " << word
486  << "Region " << i
487  << ". " << reg[i]->GetName();
488  }
489 }
#define LogDebug(id)
void primary(const G4Track *aSecondary) const
double nRusRoPreShower
T getParameter(std::string const &) const
NewTrackAction * newTA
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track *aTrack) final
StackingAction(const TrackingAction *, const edm::ParameterSet &ps, const CMSSteppingVerbose *)
std::vector< const G4Region * > lowdensRegions
void StackFilled(const G4Track *, bool isKilled) const
std::vector< double > maxTrackTimes
TrackInformationExtractor extractor
bool exists(std::string const &parameterName) const
checks if a parameter exists
double gRusRoPreShower
double nRusRoEnerLim
const G4Region * regionCastor
const G4Region * regionMuonIron
bool rrApplicable(const G4Track *, const G4Track &) const
int isItPrimaryDecayProductOrConversion(const G4Track *, const G4Track &) const
virtual ~StackingAction()
std::vector< const G4Region * > deadRegions
const double MeV
const CMSSteppingVerbose * steppingVerbose
std::vector< const G4Region * > muonRegions
std::vector< std::string > maxTimeNames
const G4Region * regionEcal
int genParticlePID() const
unsigned int numberTimes
void printRegions(const std::vector< const G4Region * > &reg, const std::string &word) const
void secondary(const G4Track *aSecondary, const G4Track &mother, int) const
double nRusRoMuonIron
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
G4VSolid * worldSolid
bool isThisRegion(const G4Region *, std::vector< const G4Region * > &) const
int ke
bool isPrimary() const
bool savePDandCinTracker
std::vector< std::string > deadRegionNames
double gRusRoEnerLim
const G4Region * regionHcal
bool isItOutOfTimeWindow(const G4Region *, const G4Track *) const
std::vector< const G4Region * > maxTimeRegions
const G4Region * regionPreShower
const TrackingAction * trackAction
int isItFromPrimary(const G4Track &, int) const
double gRusRoMuonIron
std::vector< const G4Region * > trackerRegions
std::vector< const G4Region * > caloRegions
double limitEnergyForVacuum
const G4Region * regionWorld
const G4Track * geant4Track() const