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"
20 : trackAction(trka),steppingVerbose(sv)
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 ||
88 if ( p.
exists(
"TestKillingOptions") ) {
90 .getParameter<
bool>(
"KillInCalo");
92 .getParameter<
bool>(
"KillInCaloEfH");
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= "
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: "
114 <<
" Kill Delta Ray flag: "
116 <<
" Kill hadrons/ions flag: "
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";
128 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction kill tracks with "
129 <<
"time larger than " << maxTrackTime/ns
134 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction MaxTrackTime for "
135 << maxTimeNames[
i] <<
" is "
136 << maxTrackTimes[
i] <<
" ns ";
137 maxTrackTimes[
i] *= ns;
140 if(limitEnergyForVacuum > 0.0) {
142 <<
"StackingAction LowDensity regions - kill if E < "
143 << limitEnergyForVacuum/
MeV <<
" MeV";
148 <<
"StackingAction Dead regions - kill all secondaries ";
153 <<
"StackingAction: "
154 <<
"Russian Roulette for gamma Elimit(MeV)= "
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"
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"
176 if(savePDandCinTracker) {
177 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction Tracker regions: ";
180 if(savePDandCinCalo) {
181 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction Calo regions: ";
184 if(savePDandCinMuon) {
185 edm::LogInfo(
"SimG4CoreApplication") <<
"StackingAction Muon regions: ";
188 worldSolid = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()->GetLogicalVolume()->GetSolid();
199 G4ClassificationOfNewTrack classification = fUrgent;
200 int pdg = aTrack->GetDefinition()->GetPDGEncoding();
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;
215 const G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
217 if (aTrack->GetTrackStatus() == fStopAndKill) { classification = fKill; }
218 else if (!
trackNeutrino && (abspdg == 12 || abspdg == 14 || abspdg == 16 || abspdg == 18))
219 { classification = fKill; }
224 double ke = aTrack->GetKineticEnergy();
228 classification = fKill;
231 classification = fKill;
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)) ||
243 ((pdg == 2112) && (ke <
kmaxNeutron))) { classification = fKill; }
247 && aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation) {
248 classification = fKill;
252 classification = fKill;
256 if ((pdg == 22 || abspdg == 11) && pdgMother != 11 && pdgMother != 22 &&
258 classification = fKill;
264 if(classification != fKill) {
281 if(2112 == pdg || 22 == pdg) {
282 double currentWeight = aTrack->GetWeight();
284 if(1.0 >= currentWeight) {
313 if(prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
314 if(G4UniformRand() < prob) {
315 const_cast<G4Track*
>(aTrack)->SetWeight(currentWeight/prob);
317 classification = fKill;
322 if (classification != fKill) {
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()
333 <<aTrack->GetCreatorProcess()->GetProcessSubType()
334 <<
" as " << classification <<
" Flag " << flag;
342 return classification;
352 const G4RegionStore * rs = G4RegionStore::GetInstance();
353 std::vector<G4Region*>::const_iterator rcite;
354 std::vector<std::string>::const_iterator rd;
356 for (rcite = rs->begin(); rcite != rs->end(); ++rcite) {
357 const G4Region* reg = (*rcite);
358 G4String rname = reg->GetName();
382 for (
unsigned int i=0;
i<
num;
i++) {
391 (rname ==
"BeamPipe" || rname ==
"BeamPipeVacuum"
392 || rname ==
"TrackerPixelSensRegion"
393 || rname ==
"TrackerPixelDeadRegion"
394 || rname ==
"TrackerDeadRegion" || rname ==
"TrackerSensRegion")) {
398 (rname ==
"HcalRegion" || rname ==
"EcalRegion"
399 || rname ==
"PreshowerSensRegion" || rname ==
"PreshowerRegion")) {
403 (rname ==
"MuonChamber" || rname ==
"MuonSensitive_RPC"
404 || rname ==
"MuonIron" || rname ==
"Muon"
405 || rname ==
"MuonSensitive_DT-CSC") ) {
408 if(rname ==
"BeamPipeOutside" || rname ==
"BeamPipeVacuum") {
412 if(rname == (G4String)(*rd)) {
420 std::vector<const G4Region*>& regions)
const
423 unsigned int nRegions = regions.size();
424 for(
unsigned int i=0;
i<nRegions; ++
i) {
425 if(reg == regions[
i]) {
434 const G4Track & mother)
const
440 if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) {
442 }
else if (aTrack->GetCreatorProcess()->GetProcessSubType()==fGammaConversion) {
450 const G4Track & mother)
const
456 return (22 == genID || 11 == genID || -11 == genID) ?
false :
true;
464 if (motherInfo.
isPrimary()) { flag = 3; }
478 return (aTrack->GetGlobalTime() > tofM) ?
true :
false;
484 for (
unsigned int i=0;
i<reg.size(); ++
i) {
485 edm::LogInfo(
"SimG4CoreApplication") <<
" StackingAction: " << word
487 <<
". " << reg[
i]->GetName();
void primary(const G4Track *aSecondary) const
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
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 ¶meterName) const
checks if a parameter exists
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 CMSSteppingVerbose * steppingVerbose
std::vector< const G4Region * > muonRegions
std::vector< std::string > maxTimeNames
const G4Region * regionEcal
void printRegions(const std::vector< const G4Region * > ®, const std::string &word) const
void secondary(const G4Track *aSecondary, const G4Track &mother, int) const
Abs< T >::type abs(const T &t)
bool isThisRegion(const G4Region *, std::vector< const G4Region * > &) const
std::vector< std::string > deadRegionNames
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
std::vector< const G4Region * > trackerRegions
std::vector< const G4Region * > caloRegions
double limitEnergyForVacuum
const G4Region * regionWorld
const G4Track * geant4Track() const