9 #include "G4VProcess.hh"
10 #include "G4EmProcessSubType.hh"
11 #include "G4LogicalVolumeStore.hh"
12 #include "G4RegionStore.hh"
13 #include "Randomize.hh"
29 for(
unsigned int i=0;
i<maxTrackTimes.size(); ++
i) { maxTrackTimes[
i] *= ns; }
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; }
77 if ( p.
exists(
"TestKillingOptions") ) {
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;
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: "
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 "
104 for (
unsigned int i=0;
i<maxTrackTimes.size();
i++) {
105 maxTrackTimes[
i] *= ns;
106 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction::MaxTrackTime for "
112 <<
"StackingAction: "
113 <<
"Russian Roulette for neutron Elimit(MeV)= "
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"
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"
146 G4ClassificationOfNewTrack classification = fUrgent;
150 if (aTrack->GetCreatorProcess()==0 || aTrack->GetParentID()==0) {
169 }
else if (aTrack->GetTouchable() == 0) {
171 <<
"StackingAction: no touchable for track " << aTrack->GetTrackID()
172 <<
" from " << aTrack->GetParentID()
173 <<
" with PDG code " << aTrack->GetDefinition()->GetParticleName();
174 classification = fKill;
177 int pdg = aTrack->GetDefinition()->GetPDGEncoding();
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)) ||
191 ((pdg == 2112) && (ke <
kmaxNeutron))) classification = fKill;
194 if (pdg == 12 || pdg == 14 || pdg == 16 || pdg == 18)
195 classification = fKill;
199 if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
200 aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation)
201 classification = fKill;
205 classification = fKill;
209 int pdgMother = mother->GetDefinition()->GetPDGEncoding();
210 if ( (pdg == 22 ||
std::abs(pdg) == 11) &&
214 classification = fKill;
219 if(classification != fKill && (2112 == pdg || 2212 == pdg)) {
220 double currentWeight = aTrack->GetWeight();
221 if(1.0 >= currentWeight) {
226 G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
237 G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
247 if(prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
248 if(G4UniformRand() < prob) {
249 const_cast<G4Track*
>(aTrack)->SetWeight(currentWeight/prob);
251 classification = fKill;
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;
284 return classification;
293 const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
295 std::vector<G4LogicalVolume*>::const_iterator lvcite;
296 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
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);
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);
306 if (strcmp(
"MUON",(*lvcite)->GetName().c_str()) == 0)
muon.push_back(*lvcite);
309 edm::LogInfo(
"SimG4CoreApplication") <<
"# of LV for Tracker "
310 <<
tracker.size() <<
" for Calo "
311 <<
calo.size() <<
" for Muon "
314 edm::LogInfo(
"SimG4CoreApplication") <<
"Tracker vol " <<
i <<
" name "
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();
324 const G4RegionStore * rs = G4RegionStore::GetInstance();
327 std::vector<double> tofs;
329 std::vector<G4Region*>::const_iterator rcite;
330 for (rcite = rs->begin(); rcite != rs->end(); rcite++) {
332 (*rcite)->GetName() ==
"EcalRegion") {
regionEcal = (*rcite); }
334 (*rcite)->GetName() ==
"HcalRegion") {
regionHcal = (*rcite); }
336 (*rcite)->GetName() ==
"QuadRegion") {
regionQuad = (*rcite); }
340 (*rcite)->GetName() ==
"PreshowerRegion") {
regionPreShower = (*rcite); }
342 (*rcite)->GetName() ==
"CastorRegion") {
regionCastor = (*rcite); }
344 (*rcite)->GetName() ==
"BeamPipeOutsideRegion") {
regionBeamPipeOut = (*rcite); }
346 (*rcite)->GetName() ==
"DefaultRegionForTheWorld") {
regionWorld = (*rcite); }
347 for (
unsigned int i=0;
i<
num;
i++) {
354 if (tofs.size() ==
num)
break;
357 for (
unsigned int i=0;
i<tofs.size();
i++) {
359 G4String
name =
"Unknown";
361 edm::LogInfo(
"SimG4CoreApplication") << name <<
" with pointer "
369 std::vector<G4LogicalVolume*> & lvs)
const {
372 if (lvs.size() > 0 && touch !=0) {
373 int level = ((touch->GetHistoryDepth())+1);
375 unsigned int ii = (
unsigned int)(level - 3);
376 flag = (
std::count(lvs.begin(),lvs.end(),(touch->GetVolume(ii)->GetLogicalVolume())) != 0);
383 const G4Track & mother)
const {
390 if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) flag = 1;
391 else if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
392 aTrack->GetCreatorProcess()->GetProcessSubType() == fGammaConversion) flag = 2;
411 double time = (aTrack->GetGlobalTime())/nanosecond;
414 G4Region* reg = aTrack->GetTouchable()->GetVolume(0)->GetLogicalVolume()->GetRegion();
422 if (time > tofM) flag =
true;
void primary(const G4Track *aSecondary) const
std::vector< G4LogicalVolume * > tracker
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
std::vector< G4Region * > maxTimeRegions
std::vector< double > maxTrackTimes
G4Region * regionPreShower
bool exists(std::string const ¶meterName) const
checks if a parameter exists
G4Region * regionMuonIron
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()
void secondary(const G4Track *aSecondary, const G4Track &mother, int) const
std::vector< G4LogicalVolume * > calo
G4Region * regionBeamPipeOut
std::vector< G4LogicalVolume * > muon
bool isThisVolume(const G4VTouchable *, std::vector< G4LogicalVolume * > &) const
int isItFromPrimary(const G4Track &, int) const
bool isItLongLived(const G4Track *) const
static const G4Track * track()