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 
15 #include<algorithm>
16 
17 //#define DebugLog
18 
20 
21  trackNeutrino = p.getParameter<bool>("TrackNeutrino");
22  killHeavy = p.getParameter<bool>("KillHeavy");
23  kmaxIon = p.getParameter<double>("IonThreshold")*MeV;
24  kmaxProton = p.getParameter<double>("ProtonThreshold")*MeV;
25  kmaxNeutron = p.getParameter<double>("NeutronThreshold")*MeV;
26  killDeltaRay = p.getParameter<bool>("KillDeltaRay");
27  maxTrackTime = p.getParameter<double>("MaxTrackTime")*ns;
28  maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
29  for(unsigned int i=0; i<maxTrackTimes.size(); ++i) { maxTrackTimes[i] *= ns; }
30  maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
31  savePDandCinTracker = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInTracker",false);
32  savePDandCinCalo = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo",false);
33  savePDandCinMuon = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon",false);
34  saveFirstSecondary = p.getUntrackedParameter<bool>("SaveFirstLevelSecondary",false);
35  killInCalo = false;
36  killInCaloEfH = false;
37 
38  // Russian Roulette
39  regionEcal = 0;
40  regionHcal = 0;
41  regionQuad = 0;
42  regionMuonIron = 0;
43  regionPreShower= 0;
44  regionCastor = 0;
46 
47  nRusRoEnerLim = p.getParameter<double>("RusRoNeutronEnergyLimit")*MeV;
48  pRusRoEnerLim = p.getParameter<double>("RusRoProtonEnergyLimit")*MeV;
49 
50  nRusRoEcal = p.getParameter<double>("RusRoEcalNeutron");
51  nRusRoHcal = p.getParameter<double>("RusRoHcalNeutron");
52  nRusRoQuad = p.getParameter<double>("RusRoQuadNeutron");
53  nRusRoMuonIron = p.getParameter<double>("RusRoMuonIronNeutron");
54  nRusRoPreShower = p.getParameter<double>("RusRoPreShowerNeutron");
55  nRusRoCastor = p.getParameter<double>("RusRoCastorNeutron");
56  nRusRoBeam = p.getParameter<double>("RusRoBeamPipeOutNeutron");
57  nRusRoWorld = p.getParameter<double>("RusRoWorldNeutron");
58 
59  pRusRoEcal = p.getParameter<double>("RusRoEcalProton");
60  pRusRoHcal = p.getParameter<double>("RusRoHcalProton");
61  pRusRoQuad = p.getParameter<double>("RusRoQuadProton");
62  pRusRoMuonIron = p.getParameter<double>("RusRoMuonIronProton");
63  pRusRoPreShower = p.getParameter<double>("RusRoPreShowerProton");
64  pRusRoCastor = p.getParameter<double>("RusRoCastorProton");
65  pRusRoBeam = p.getParameter<double>("RusRoBeamPipeOutProton");
66  pRusRoWorld = p.getParameter<double>("RusRoWorldProton");
67 
68  nRRactive = false;
69  pRRactive = false;
70  if(nRusRoEcal < 1.0 || nRusRoHcal < 1.0 || nRusRoQuad < 1.0 ||
71  nRusRoMuonIron < 1.0 || nRusRoPreShower < 1.0 || nRusRoCastor < 1.0 ||
72  nRusRoBeam < 1.0 || nRusRoWorld < 1.0) { nRRactive = true; }
73  if(pRusRoEcal < 1.0 || pRusRoHcal < 1.0 || pRusRoQuad < 1.0 ||
74  pRusRoMuonIron < 1.0 || pRusRoPreShower < 1.0 || pRusRoCastor < 1.0 ||
75  pRusRoBeam < 1.0 || pRusRoWorld < 1.0) { pRRactive = true; }
76 
77  if ( p.exists("TestKillingOptions") ) {
78 
79  killInCalo = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCalo");
80  killInCaloEfH = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCaloEfH");
81  edm::LogWarning("SimG4CoreApplication") << " *** Activating special test killing options in StackingAction \n"
82  << " *** Kill secondaries in Calorimetetrs volume = " << killInCalo << "\n"
83  << " *** Kill electromagnetic secondaries from hadrons in Calorimeters volume = " << killInCaloEfH;
84 
85  }
86 
87  edm::LogInfo("SimG4CoreApplication") << "StackingAction initiated with"
88  << " flag for saving decay products in "
89  << " Tracker: " << savePDandCinTracker
90  << " in Calo: " << savePDandCinCalo
91  << " in Muon: " << savePDandCinMuon
92  << "\n saveFirstSecondary"
93  << ": " << saveFirstSecondary
94  << " Flag for tracking neutrino: "
95  << trackNeutrino << " Killing Flag "
96  << killHeavy << " protons below "
97  << kmaxProton <<" MeV, neutrons below "
98  << kmaxNeutron << " MeV and ions"
99  << " below " << kmaxIon << " MeV\n"
100  << " kill tracks with "
101  << "time larger than " << maxTrackTime
102  << " ns and kill Delta Ray flag set to "
103  << killDeltaRay;
104  for (unsigned int i=0; i<maxTrackTimes.size(); i++) {
105  maxTrackTimes[i] *= ns;
106  edm::LogInfo("SimG4CoreApplication") << "StackingAction::MaxTrackTime for "
107  << maxTimeNames[i] << " is "
108  << maxTrackTimes[i];
109  }
110  if(nRRactive) {
111  edm::LogInfo("SimG4CoreApplication")
112  << "StackingAction: "
113  << "Russian Roulette for neutron Elimit(MeV)= "
114  << nRusRoEnerLim/MeV << "\n"
115  << " ECAL Prob= " << nRusRoEcal << "\n"
116  << " HCAL Prob= " << nRusRoHcal << "\n"
117  << " QUAD Prob= " << nRusRoQuad << "\n"
118  << " MuonIron Prob= " << nRusRoMuonIron << "\n"
119  << " PreShower Prob= " << nRusRoPreShower << "\n"
120  << " CASTOR Prob= " << nRusRoCastor << "\n"
121  << " BeamPipeOut Prob= " << nRusRoBeam << "\n"
122  << " World Prob= " << nRusRoWorld;
123  }
124  if(pRRactive) {
125  edm::LogInfo("SimG4CoreApplication")
126  << "StackingAction: "
127  << "Russian Roulette for proton Elimit(MeV)= "
128  << pRusRoEnerLim/MeV << "\n"
129  << " ECAL Prob= " << pRusRoEcal << "\n"
130  << " HCAL Prob= " << pRusRoHcal << "\n"
131  << " QUAD Prob= " << pRusRoQuad << "\n"
132  << " MuonIron Prob= " << pRusRoMuonIron << "\n"
133  << " PreShower Prob= " << pRusRoPreShower << "\n"
134  << " CASTOR Prob= " << pRusRoCastor << "\n"
135  << " BeamPipeOut Prob= " << pRusRoBeam << "\n"
136  << " World Prob= " << pRusRoWorld;
137  }
138  initPointer();
139 }
140 
142 
143 G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track * aTrack) {
144 
145  // G4 interface part
146  G4ClassificationOfNewTrack classification = fUrgent;
147  int flag = 0;
148 
149  NewTrackAction newTA;
150  if (aTrack->GetCreatorProcess()==0 || aTrack->GetParentID()==0) {
151  /*
152  std::cout << "StackingAction: primary weight= "
153  << aTrack->GetWeight() << " "
154  << aTrack->GetDefinition()->GetParticleName()
155  << " " << aTrack->GetKineticEnergy()
156  << " Id=" << aTrack->GetTrackID()
157  << " trackInfo " << aTrack->GetUserInformation()
158  << std::endl;
159  */
160  newTA.primary(aTrack);
161  /*
162  if (!trackNeutrino) {
163  int pdg = std::abs(aTrack->GetDefinition()->GetPDGEncoding());
164  if (pdg == 12 || pdg == 14 || pdg == 16 || pdg == 18) {
165  classification = fKill;
166  }
167  }
168  */
169  } else if (aTrack->GetTouchable() == 0) {
170  edm::LogError("SimG4CoreApplication")
171  << "StackingAction: no touchable for track " << aTrack->GetTrackID()
172  << " from " << aTrack->GetParentID()
173  << " with PDG code " << aTrack->GetDefinition()->GetParticleName();
174  classification = fKill;
175  } else {
176  const G4Track * mother = CurrentG4Track::track();
177  int pdg = aTrack->GetDefinition()->GetPDGEncoding();
178  if ((savePDandCinTracker && isThisVolume(aTrack->GetTouchable(),tracker))||
179  (savePDandCinCalo && isThisVolume(aTrack->GetTouchable(),calo)) ||
180  (savePDandCinMuon && isThisVolume(aTrack->GetTouchable(),muon)))
181  flag = isItPrimaryDecayProductOrConversion(aTrack, *mother);
182  if (saveFirstSecondary) flag = isItFromPrimary(*mother, flag);
183  newTA.secondary(aTrack, *mother, flag);
184 
185  if (aTrack->GetTrackStatus() == fStopAndKill) classification = fKill;
186  if (killHeavy && classification != fKill) {
187  double ke = aTrack->GetKineticEnergy()/MeV;
188  if (((pdg/1000000000 == 1) && (((pdg/10000)%100) > 0) &&
189  (((pdg/10)%100) > 0) && (ke<kmaxIon)) ||
190  ((pdg == 2212) && (ke < kmaxProton)) ||
191  ((pdg == 2112) && (ke < kmaxNeutron))) classification = fKill;
192  }
193  if (!trackNeutrino && classification != fKill) {
194  if (pdg == 12 || pdg == 14 || pdg == 16 || pdg == 18)
195  classification = fKill;
196  }
197  if (isItLongLived(aTrack)) classification = fKill;
198  if (killDeltaRay) {
199  if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
200  aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation)
201  classification = fKill;
202  }
203  if (killInCalo && classification != fKill) {
204  if ( isThisVolume(aTrack->GetTouchable(),calo)) {
205  classification = fKill;
206  }
207  }
208  if (killInCaloEfH && classification != fKill) {
209  int pdgMother = mother->GetDefinition()->GetPDGEncoding();
210  if ( (pdg == 22 || std::abs(pdg) == 11) &&
211  (std::abs(pdgMother) < 11 || std::abs(pdgMother) > 17) &&
212  pdgMother != 22 ) {
213  if ( isThisVolume(aTrack->GetTouchable(),calo)) {
214  classification = fKill;
215  }
216  }
217  }
218  // Russian roulette
219  if(classification != fKill && (2112 == pdg || 2212 == pdg)) {
220  double currentWeight = aTrack->GetWeight();
221  if(1.0 >= currentWeight) {
222  double prob = 1.0;
223  double elim = 0.0;
224  if(nRRactive && pdg == 2112) {
225  elim = nRusRoEnerLim;
226  G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
227  if(reg == regionEcal) { prob = nRusRoEcal; }
228  else if(reg == regionHcal) { prob = nRusRoHcal; }
229  else if(reg == regionQuad) { prob = nRusRoQuad; }
230  else if(reg == regionMuonIron) { prob = nRusRoMuonIron; }
231  else if(reg == regionPreShower) { prob = nRusRoPreShower; }
232  else if(reg == regionCastor) { prob = nRusRoCastor; }
233  else if(reg == regionBeamPipeOut) { prob = nRusRoBeam; }
234  else if(reg == regionWorld) { prob = nRusRoWorld; }
235  } else if(pRRactive && pdg == 2212) {
236  elim = pRusRoEnerLim;
237  G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
238  if(reg == regionEcal) { prob = pRusRoEcal; }
239  else if(reg == regionHcal) { prob = pRusRoHcal; }
240  else if(reg == regionQuad) { prob = pRusRoQuad; }
241  else if(reg == regionMuonIron) { prob = pRusRoMuonIron; }
242  else if(reg == regionPreShower) { prob = pRusRoPreShower; }
243  else if(reg == regionCastor) { prob = pRusRoCastor; }
244  else if(reg == regionBeamPipeOut) { prob = pRusRoBeam; }
245  else if(reg == regionWorld) { prob = pRusRoWorld; }
246  }
247  if(prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
248  if(G4UniformRand() < prob) {
249  const_cast<G4Track*>(aTrack)->SetWeight(currentWeight/prob);
250  } else {
251  classification = fKill;
252  }
253  }
254  }
255  }
256  /*
257  double wt2 = aTrack->GetWeight();
258  if(wt2 != 1.0) {
259  G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
260  std::cout << "StackingAction: weight= " << wt2 << " "
261  << aTrack->GetDefinition()->GetParticleName()
262  << " " << aTrack->GetKineticEnergy()
263  << " Id=" << aTrack->GetTrackID()
264  << " IdP=" << aTrack->GetParentID();
265  const G4VProcess* pr = aTrack->GetCreatorProcess();
266  if(pr) std::cout << " from " << pr->GetProcessName();
267  if(reg) std::cout << " in " << reg->GetName()
268  << " trackInfo " << aTrack->GetUserInformation();
269  std::cout << std::endl;
270  }
271  */
272 #ifdef DebugLog
273  LogDebug("SimG4CoreApplication") << "StackingAction:Classify Track "
274  << aTrack->GetTrackID() << " Parent "
275  << aTrack->GetParentID() << " Type "
276  << aTrack->GetDefinition()->GetParticleName()
277  << " K.E. " << aTrack->GetKineticEnergy()/MeV
278  << " MeV from process/subprocess "
279  << aTrack->GetCreatorProcess()->GetProcessType() << "|"
280  << aTrack->GetCreatorProcess()->GetProcessSubType()
281  << " as " << classification << " Flag " << flag;
282 #endif
283  }
284  return classification;
285 }
286 
288 
290 
292 
293  const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
294  if (lvs) {
295  std::vector<G4LogicalVolume*>::const_iterator lvcite;
296  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
297  if (savePDandCinTracker) {
298  if (strcmp("Tracker",(*lvcite)->GetName().c_str()) == 0) tracker.push_back(*lvcite);
299  if (strcmp("BEAM",(*lvcite)->GetName().substr(0,4).c_str()) == 0) tracker.push_back(*lvcite);
300  }
302  if (strcmp("CALO",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
303  if (strcmp("VCAL",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
304  }
305  if (savePDandCinMuon) {
306  if (strcmp("MUON",(*lvcite)->GetName().c_str()) == 0) muon.push_back(*lvcite);
307  }
308  }
309  edm::LogInfo("SimG4CoreApplication") << "# of LV for Tracker "
310  << tracker.size() << " for Calo "
311  << calo.size() << " for Muon "
312  << muon.size();
313  for (unsigned int i=0; i<tracker.size(); ++i)
314  edm::LogInfo("SimG4CoreApplication") << "Tracker vol " << i << " name "
315  << tracker[i]->GetName();
316  for (unsigned int i=0; i<calo.size(); ++i)
317  edm::LogInfo("SimG4CoreApplication") << "Calorimeter vol " <<i <<" name "
318  << calo[i]->GetName();
319  for (unsigned int i=0; i<muon.size(); ++i)
320  edm::LogInfo("SimG4CoreApplication") << "Muon vol " << i << " name "
321  << muon[i]->GetName();
322  }
323 
324  const G4RegionStore * rs = G4RegionStore::GetInstance();
325  unsigned int num = maxTimeNames.size();
326  if (num > 0) {
327  std::vector<double> tofs;
328  if (rs) {
329  std::vector<G4Region*>::const_iterator rcite;
330  for (rcite = rs->begin(); rcite != rs->end(); rcite++) {
331  if ((nRusRoEcal < 1.0 || pRusRoEcal < 1.0) &&
332  (*rcite)->GetName() == "EcalRegion") { regionEcal = (*rcite); }
333  if ((nRusRoHcal < 1.0 || pRusRoHcal < 1.0) &&
334  (*rcite)->GetName() == "HcalRegion") { regionHcal = (*rcite); }
335  if ((nRusRoQuad < 1.0 || pRusRoQuad < 1.0) &&
336  (*rcite)->GetName() == "QuadRegion") { regionQuad = (*rcite); }
337  if ((nRusRoMuonIron < 1.0 || pRusRoMuonIron < 1.0) &&
338  (*rcite)->GetName() == "MuonIron") { regionMuonIron = (*rcite); }
339  if ((nRusRoPreShower < 1.0 || pRusRoPreShower < 1.0) &&
340  (*rcite)->GetName() == "PreshowerRegion") { regionPreShower = (*rcite); }
341  if ((nRusRoCastor < 1.0 || pRusRoCastor < 1.0) &&
342  (*rcite)->GetName() == "CastorRegion") { regionCastor = (*rcite); }
343  if ((nRusRoBeam < 1.0 || pRusRoBeam < 1.0) &&
344  (*rcite)->GetName() == "BeamPipeOutsideRegion") { regionBeamPipeOut = (*rcite); }
345  if ((nRusRoWorld < 1.0 || pRusRoWorld < 1.0) &&
346  (*rcite)->GetName() == "DefaultRegionForTheWorld") { regionWorld = (*rcite); }
347  for (unsigned int i=0; i<num; i++) {
348  if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
349  maxTimeRegions.push_back(*rcite);
350  tofs.push_back(maxTrackTimes[i]);
351  break;
352  }
353  }
354  if (tofs.size() == num) break;
355  }
356  }
357  for (unsigned int i=0; i<tofs.size(); i++) {
358  maxTrackTimes[i] = tofs[i];
359  G4String name = "Unknown";
360  if (maxTimeRegions[i]) name = maxTimeRegions[i]->GetName();
361  edm::LogInfo("SimG4CoreApplication") << name << " with pointer "
362  << maxTimeRegions[i]<<" KE cut off "
363  << maxTrackTimes[i];
364  }
365  }
366 }
367 
368 bool StackingAction::isThisVolume(const G4VTouchable* touch,
369  std::vector<G4LogicalVolume*> & lvs) const {
370 
371  bool flag = false;
372  if (lvs.size() > 0 && touch !=0) {
373  int level = ((touch->GetHistoryDepth())+1);
374  if (level >= 3) {
375  unsigned int ii = (unsigned int)(level - 3);
376  flag = (std::count(lvs.begin(),lvs.end(),(touch->GetVolume(ii)->GetLogicalVolume())) != 0);
377  }
378  }
379  return flag;
380 }
381 
383  const G4Track & mother) const {
384 
385  int flag = 0;
386  TrackInformationExtractor extractor;
387  const TrackInformation & motherInfo(extractor(mother));
388  // Check whether mother is a primary
389  if (motherInfo.isPrimary()) {
390  if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) flag = 1;
391  else if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
392  aTrack->GetCreatorProcess()->GetProcessSubType() == fGammaConversion) flag = 2;
393  }
394  return flag;
395 }
396 
397 int StackingAction::isItFromPrimary(const G4Track & mother, int flagIn) const {
398 
399  int flag = flagIn;
400  if (flag != 1) {
401  TrackInformationExtractor extractor;
402  const TrackInformation & motherInfo(extractor(mother));
403  if (motherInfo.isPrimary()) flag = 3;
404  }
405  return flag;
406 }
407 
408 bool StackingAction::isItLongLived(const G4Track * aTrack) const {
409 
410  bool flag = false;
411  double time = (aTrack->GetGlobalTime())/nanosecond;
412  double tofM = maxTrackTime;
413  if (maxTimeRegions.size() > 0) {
414  G4Region* reg = aTrack->GetTouchable()->GetVolume(0)->GetLogicalVolume()->GetRegion();
415  for (unsigned int i=0; i<maxTimeRegions.size(); i++) {
416  if (reg == maxTimeRegions[i]) {
417  tofM = maxTrackTimes[i];
418  break;
419  }
420  }
421  }
422  if (time > tofM) flag = true;
423  return flag;
424 }
425 
426 
#define LogDebug(id)
void primary(const G4Track *aSecondary) const
double nRusRoPreShower
std::vector< G4LogicalVolume * > tracker
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
std::vector< G4Region * > maxTimeRegions
long int flag
Definition: mlp_lapack.h:47
std::vector< double > maxTrackTimes
G4Region * regionPreShower
bool exists(std::string const &parameterName) const
checks if a parameter exists
#define abs(x)
Definition: mlp_lapack.h:159
G4Region * regionHcal
double nRusRoEnerLim
double pRusRoEnerLim
G4Region * regionMuonIron
int ii
Definition: cuy.py:588
G4Region * regionEcal
int isItPrimaryDecayProductOrConversion(const G4Track *, const G4Track &) const
virtual ~StackingAction()
virtual G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track *aTrack)
StackingAction(const edm::ParameterSet &ps)
std::vector< std::string > maxTimeNames
virtual void PrepareNewEvent()
double pRusRoPreShower
void secondary(const G4Track *aSecondary, const G4Track &mother, int) const
double nRusRoMuonIron
G4Region * regionCastor
G4Region * regionWorld
std::vector< G4LogicalVolume * > calo
long long int num
Definition: procUtils.cc:71
int ke
bool isPrimary() const
bool savePDandCinTracker
G4Region * regionBeamPipeOut
std::vector< G4LogicalVolume * > muon
bool isThisVolume(const G4VTouchable *, std::vector< G4LogicalVolume * > &) const
G4Region * regionQuad
virtual void NewStage()
tuple level
Definition: testEve_cfg.py:34
int isItFromPrimary(const G4Track &, int) const
double pRusRoMuonIron
bool isItLongLived(const G4Track *) const
static const G4Track * track()