Go to the documentation of this file.00001 #include "SimG4Core/Application/interface/StackingAction.h"
00002 #include "SimG4Core/Notification/interface/CurrentG4Track.h"
00003 #include "SimG4Core/Notification/interface/NewTrackAction.h"
00004 #include "SimG4Core/Notification/interface/TrackInformation.h"
00005 #include "SimG4Core/Notification/interface/TrackInformationExtractor.h"
00006
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008
00009 #include "G4VProcess.hh"
00010 #include "G4EmProcessSubType.hh"
00011 #include "G4LogicalVolumeStore.hh"
00012 #include "G4RegionStore.hh"
00013 #include "Randomize.hh"
00014
00015 #include<algorithm>
00016
00017
00018
00019 StackingAction::StackingAction(const edm::ParameterSet & p) {
00020
00021 trackNeutrino = p.getParameter<bool>("TrackNeutrino");
00022 killHeavy = p.getParameter<bool>("KillHeavy");
00023 kmaxIon = p.getParameter<double>("IonThreshold")*MeV;
00024 kmaxProton = p.getParameter<double>("ProtonThreshold")*MeV;
00025 kmaxNeutron = p.getParameter<double>("NeutronThreshold")*MeV;
00026 killDeltaRay = p.getParameter<bool>("KillDeltaRay");
00027 maxTrackTime = p.getParameter<double>("MaxTrackTime")*ns;
00028 maxTrackTimes = p.getParameter<std::vector<double> >("MaxTrackTimes");
00029 for(unsigned int i=0; i<maxTrackTimes.size(); ++i) { maxTrackTimes[i] *= ns; }
00030 maxTimeNames = p.getParameter<std::vector<std::string> >("MaxTimeNames");
00031 savePDandCinTracker = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInTracker",false);
00032 savePDandCinCalo = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInCalo",false);
00033 savePDandCinMuon = p.getUntrackedParameter<bool>("SavePrimaryDecayProductsAndConversionsInMuon",false);
00034 saveFirstSecondary = p.getUntrackedParameter<bool>("SaveFirstLevelSecondary",false);
00035 killInCalo = false;
00036 killInCaloEfH = false;
00037
00038
00039 regionEcal = 0;
00040 regionHcal = 0;
00041 regionQuad = 0;
00042 regionMuonIron = 0;
00043 regionPreShower= 0;
00044 regionCastor = 0;
00045 regionBeamPipeOut = 0;
00046
00047 nRusRoEnerLim = p.getParameter<double>("RusRoNeutronEnergyLimit")*MeV;
00048 pRusRoEnerLim = p.getParameter<double>("RusRoProtonEnergyLimit")*MeV;
00049
00050 nRusRoEcal = p.getParameter<double>("RusRoEcalNeutron");
00051 nRusRoHcal = p.getParameter<double>("RusRoHcalNeutron");
00052 nRusRoQuad = p.getParameter<double>("RusRoQuadNeutron");
00053 nRusRoMuonIron = p.getParameter<double>("RusRoMuonIronNeutron");
00054 nRusRoPreShower = p.getParameter<double>("RusRoPreShowerNeutron");
00055 nRusRoCastor = p.getParameter<double>("RusRoCastorNeutron");
00056 nRusRoBeam = p.getParameter<double>("RusRoBeamPipeOutNeutron");
00057 nRusRoWorld = p.getParameter<double>("RusRoWorldNeutron");
00058
00059 pRusRoEcal = p.getParameter<double>("RusRoEcalProton");
00060 pRusRoHcal = p.getParameter<double>("RusRoHcalProton");
00061 pRusRoQuad = p.getParameter<double>("RusRoQuadProton");
00062 pRusRoMuonIron = p.getParameter<double>("RusRoMuonIronProton");
00063 pRusRoPreShower = p.getParameter<double>("RusRoPreShowerProton");
00064 pRusRoCastor = p.getParameter<double>("RusRoCastorProton");
00065 pRusRoBeam = p.getParameter<double>("RusRoBeamPipeOutProton");
00066 pRusRoWorld = p.getParameter<double>("RusRoWorldProton");
00067
00068 nRRactive = false;
00069 pRRactive = false;
00070 if(nRusRoEcal < 1.0 || nRusRoHcal < 1.0 || nRusRoQuad < 1.0 ||
00071 nRusRoMuonIron < 1.0 || nRusRoPreShower < 1.0 || nRusRoCastor < 1.0 ||
00072 nRusRoBeam < 1.0 || nRusRoWorld < 1.0) { nRRactive = true; }
00073 if(pRusRoEcal < 1.0 || pRusRoHcal < 1.0 || pRusRoQuad < 1.0 ||
00074 pRusRoMuonIron < 1.0 || pRusRoPreShower < 1.0 || pRusRoCastor < 1.0 ||
00075 pRusRoBeam < 1.0 || pRusRoWorld < 1.0) { pRRactive = true; }
00076
00077 if ( p.exists("TestKillingOptions") ) {
00078
00079 killInCalo = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCalo");
00080 killInCaloEfH = (p.getParameter<edm::ParameterSet>("TestKillingOptions")).getParameter<bool>("KillInCaloEfH");
00081 edm::LogWarning("SimG4CoreApplication") << " *** Activating special test killing options in StackingAction \n"
00082 << " *** Kill secondaries in Calorimetetrs volume = " << killInCalo << "\n"
00083 << " *** Kill electromagnetic secondaries from hadrons in Calorimeters volume = " << killInCaloEfH;
00084
00085 }
00086
00087 edm::LogInfo("SimG4CoreApplication") << "StackingAction initiated with"
00088 << " flag for saving decay products in "
00089 << " Tracker: " << savePDandCinTracker
00090 << " in Calo: " << savePDandCinCalo
00091 << " in Muon: " << savePDandCinMuon
00092 << "\n saveFirstSecondary"
00093 << ": " << saveFirstSecondary
00094 << " Flag for tracking neutrino: "
00095 << trackNeutrino << " Killing Flag "
00096 << killHeavy << " protons below "
00097 << kmaxProton <<" MeV, neutrons below "
00098 << kmaxNeutron << " MeV and ions"
00099 << " below " << kmaxIon << " MeV\n"
00100 << " kill tracks with "
00101 << "time larger than " << maxTrackTime
00102 << " ns and kill Delta Ray flag set to "
00103 << killDeltaRay;
00104 for (unsigned int i=0; i<maxTrackTimes.size(); i++) {
00105 maxTrackTimes[i] *= ns;
00106 edm::LogInfo("SimG4CoreApplication") << "StackingAction::MaxTrackTime for "
00107 << maxTimeNames[i] << " is "
00108 << maxTrackTimes[i];
00109 }
00110 if(nRRactive) {
00111 edm::LogInfo("SimG4CoreApplication")
00112 << "StackingAction: "
00113 << "Russian Roulette for neutron Elimit(MeV)= "
00114 << nRusRoEnerLim/MeV << "\n"
00115 << " ECAL Prob= " << nRusRoEcal << "\n"
00116 << " HCAL Prob= " << nRusRoHcal << "\n"
00117 << " QUAD Prob= " << nRusRoQuad << "\n"
00118 << " MuonIron Prob= " << nRusRoMuonIron << "\n"
00119 << " PreShower Prob= " << nRusRoPreShower << "\n"
00120 << " CASTOR Prob= " << nRusRoCastor << "\n"
00121 << " BeamPipeOut Prob= " << nRusRoBeam << "\n"
00122 << " World Prob= " << nRusRoWorld;
00123 }
00124 if(pRRactive) {
00125 edm::LogInfo("SimG4CoreApplication")
00126 << "StackingAction: "
00127 << "Russian Roulette for proton Elimit(MeV)= "
00128 << pRusRoEnerLim/MeV << "\n"
00129 << " ECAL Prob= " << pRusRoEcal << "\n"
00130 << " HCAL Prob= " << pRusRoHcal << "\n"
00131 << " QUAD Prob= " << pRusRoQuad << "\n"
00132 << " MuonIron Prob= " << pRusRoMuonIron << "\n"
00133 << " PreShower Prob= " << pRusRoPreShower << "\n"
00134 << " CASTOR Prob= " << pRusRoCastor << "\n"
00135 << " BeamPipeOut Prob= " << pRusRoBeam << "\n"
00136 << " World Prob= " << pRusRoWorld;
00137 }
00138 initPointer();
00139 }
00140
00141 StackingAction::~StackingAction() {}
00142
00143 G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack(const G4Track * aTrack) {
00144
00145
00146 G4ClassificationOfNewTrack classification = fUrgent;
00147 int flag = 0;
00148
00149 NewTrackAction newTA;
00150 if (aTrack->GetCreatorProcess()==0 || aTrack->GetParentID()==0) {
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 newTA.primary(aTrack);
00161
00162
00163
00164
00165
00166
00167
00168
00169 } else if (aTrack->GetTouchable() == 0) {
00170 edm::LogError("SimG4CoreApplication")
00171 << "StackingAction: no touchable for track " << aTrack->GetTrackID()
00172 << " from " << aTrack->GetParentID()
00173 << " with PDG code " << aTrack->GetDefinition()->GetParticleName();
00174 classification = fKill;
00175 } else {
00176 const G4Track * mother = CurrentG4Track::track();
00177 int pdg = aTrack->GetDefinition()->GetPDGEncoding();
00178 if ((savePDandCinTracker && isThisVolume(aTrack->GetTouchable(),tracker))||
00179 (savePDandCinCalo && isThisVolume(aTrack->GetTouchable(),calo)) ||
00180 (savePDandCinMuon && isThisVolume(aTrack->GetTouchable(),muon)))
00181 flag = isItPrimaryDecayProductOrConversion(aTrack, *mother);
00182 if (saveFirstSecondary) flag = isItFromPrimary(*mother, flag);
00183 newTA.secondary(aTrack, *mother, flag);
00184
00185 if (aTrack->GetTrackStatus() == fStopAndKill) classification = fKill;
00186 if (killHeavy && classification != fKill) {
00187 double ke = aTrack->GetKineticEnergy()/MeV;
00188 if (((pdg/1000000000 == 1) && (((pdg/10000)%100) > 0) &&
00189 (((pdg/10)%100) > 0) && (ke<kmaxIon)) ||
00190 ((pdg == 2212) && (ke < kmaxProton)) ||
00191 ((pdg == 2112) && (ke < kmaxNeutron))) classification = fKill;
00192 }
00193 if (!trackNeutrino && classification != fKill) {
00194 if (pdg == 12 || pdg == 14 || pdg == 16 || pdg == 18)
00195 classification = fKill;
00196 }
00197 if (isItLongLived(aTrack)) classification = fKill;
00198 if (killDeltaRay) {
00199 if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
00200 aTrack->GetCreatorProcess()->GetProcessSubType() == fIonisation)
00201 classification = fKill;
00202 }
00203 if (killInCalo && classification != fKill) {
00204 if ( isThisVolume(aTrack->GetTouchable(),calo)) {
00205 classification = fKill;
00206 }
00207 }
00208 if (killInCaloEfH && classification != fKill) {
00209 int pdgMother = mother->GetDefinition()->GetPDGEncoding();
00210 if ( (pdg == 22 || std::abs(pdg) == 11) &&
00211 (std::abs(pdgMother) < 11 || std::abs(pdgMother) > 17) &&
00212 pdgMother != 22 ) {
00213 if ( isThisVolume(aTrack->GetTouchable(),calo)) {
00214 classification = fKill;
00215 }
00216 }
00217 }
00218
00219 if(classification != fKill && (2112 == pdg || 2212 == pdg)) {
00220 double currentWeight = aTrack->GetWeight();
00221 if(1.0 >= currentWeight) {
00222 double prob = 1.0;
00223 double elim = 0.0;
00224 if(nRRactive && pdg == 2112) {
00225 elim = nRusRoEnerLim;
00226 G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
00227 if(reg == regionEcal) { prob = nRusRoEcal; }
00228 else if(reg == regionHcal) { prob = nRusRoHcal; }
00229 else if(reg == regionQuad) { prob = nRusRoQuad; }
00230 else if(reg == regionMuonIron) { prob = nRusRoMuonIron; }
00231 else if(reg == regionPreShower) { prob = nRusRoPreShower; }
00232 else if(reg == regionCastor) { prob = nRusRoCastor; }
00233 else if(reg == regionBeamPipeOut) { prob = nRusRoBeam; }
00234 else if(reg == regionWorld) { prob = nRusRoWorld; }
00235 } else if(pRRactive && pdg == 2212) {
00236 elim = pRusRoEnerLim;
00237 G4Region* reg = aTrack->GetVolume()->GetLogicalVolume()->GetRegion();
00238 if(reg == regionEcal) { prob = pRusRoEcal; }
00239 else if(reg == regionHcal) { prob = pRusRoHcal; }
00240 else if(reg == regionQuad) { prob = pRusRoQuad; }
00241 else if(reg == regionMuonIron) { prob = pRusRoMuonIron; }
00242 else if(reg == regionPreShower) { prob = pRusRoPreShower; }
00243 else if(reg == regionCastor) { prob = pRusRoCastor; }
00244 else if(reg == regionBeamPipeOut) { prob = pRusRoBeam; }
00245 else if(reg == regionWorld) { prob = pRusRoWorld; }
00246 }
00247 if(prob < 1.0 && aTrack->GetKineticEnergy() < elim) {
00248 if(G4UniformRand() < prob) {
00249 const_cast<G4Track*>(aTrack)->SetWeight(currentWeight/prob);
00250 } else {
00251 classification = fKill;
00252 }
00253 }
00254 }
00255 }
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 #ifdef DebugLog
00273 LogDebug("SimG4CoreApplication") << "StackingAction:Classify Track "
00274 << aTrack->GetTrackID() << " Parent "
00275 << aTrack->GetParentID() << " Type "
00276 << aTrack->GetDefinition()->GetParticleName()
00277 << " K.E. " << aTrack->GetKineticEnergy()/MeV
00278 << " MeV from process/subprocess "
00279 << aTrack->GetCreatorProcess()->GetProcessType() << "|"
00280 << aTrack->GetCreatorProcess()->GetProcessSubType()
00281 << " as " << classification << " Flag " << flag;
00282 #endif
00283 }
00284 return classification;
00285 }
00286
00287 void StackingAction::NewStage() {}
00288
00289 void StackingAction::PrepareNewEvent() {}
00290
00291 void StackingAction::initPointer() {
00292
00293 const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00294 if (lvs) {
00295 std::vector<G4LogicalVolume*>::const_iterator lvcite;
00296 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
00297 if (savePDandCinTracker) {
00298 if (strcmp("Tracker",(*lvcite)->GetName().c_str()) == 0) tracker.push_back(*lvcite);
00299 if (strcmp("BEAM",(*lvcite)->GetName().substr(0,4).c_str()) == 0) tracker.push_back(*lvcite);
00300 }
00301 if (savePDandCinCalo || killInCalo || killInCaloEfH ) {
00302 if (strcmp("CALO",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
00303 if (strcmp("VCAL",(*lvcite)->GetName().c_str()) == 0) calo.push_back(*lvcite);
00304 }
00305 if (savePDandCinMuon) {
00306 if (strcmp("MUON",(*lvcite)->GetName().c_str()) == 0) muon.push_back(*lvcite);
00307 }
00308 }
00309 edm::LogInfo("SimG4CoreApplication") << "# of LV for Tracker "
00310 << tracker.size() << " for Calo "
00311 << calo.size() << " for Muon "
00312 << muon.size();
00313 for (unsigned int i=0; i<tracker.size(); ++i)
00314 edm::LogInfo("SimG4CoreApplication") << "Tracker vol " << i << " name "
00315 << tracker[i]->GetName();
00316 for (unsigned int i=0; i<calo.size(); ++i)
00317 edm::LogInfo("SimG4CoreApplication") << "Calorimeter vol " <<i <<" name "
00318 << calo[i]->GetName();
00319 for (unsigned int i=0; i<muon.size(); ++i)
00320 edm::LogInfo("SimG4CoreApplication") << "Muon vol " << i << " name "
00321 << muon[i]->GetName();
00322 }
00323
00324 const G4RegionStore * rs = G4RegionStore::GetInstance();
00325 unsigned int num = maxTimeNames.size();
00326 if (num > 0) {
00327 std::vector<double> tofs;
00328 if (rs) {
00329 std::vector<G4Region*>::const_iterator rcite;
00330 for (rcite = rs->begin(); rcite != rs->end(); rcite++) {
00331 if ((nRusRoEcal < 1.0 || pRusRoEcal < 1.0) &&
00332 (*rcite)->GetName() == "EcalRegion") { regionEcal = (*rcite); }
00333 if ((nRusRoHcal < 1.0 || pRusRoHcal < 1.0) &&
00334 (*rcite)->GetName() == "HcalRegion") { regionHcal = (*rcite); }
00335 if ((nRusRoQuad < 1.0 || pRusRoQuad < 1.0) &&
00336 (*rcite)->GetName() == "QuadRegion") { regionQuad = (*rcite); }
00337 if ((nRusRoMuonIron < 1.0 || pRusRoMuonIron < 1.0) &&
00338 (*rcite)->GetName() == "MuonIron") { regionMuonIron = (*rcite); }
00339 if ((nRusRoPreShower < 1.0 || pRusRoPreShower < 1.0) &&
00340 (*rcite)->GetName() == "PreshowerRegion") { regionPreShower = (*rcite); }
00341 if ((nRusRoCastor < 1.0 || pRusRoCastor < 1.0) &&
00342 (*rcite)->GetName() == "CastorRegion") { regionCastor = (*rcite); }
00343 if ((nRusRoBeam < 1.0 || pRusRoBeam < 1.0) &&
00344 (*rcite)->GetName() == "BeamPipeOutsideRegion") { regionBeamPipeOut = (*rcite); }
00345 if ((nRusRoWorld < 1.0 || pRusRoWorld < 1.0) &&
00346 (*rcite)->GetName() == "DefaultRegionForTheWorld") { regionWorld = (*rcite); }
00347 for (unsigned int i=0; i<num; i++) {
00348 if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
00349 maxTimeRegions.push_back(*rcite);
00350 tofs.push_back(maxTrackTimes[i]);
00351 break;
00352 }
00353 }
00354 if (tofs.size() == num) break;
00355 }
00356 }
00357 for (unsigned int i=0; i<tofs.size(); i++) {
00358 maxTrackTimes[i] = tofs[i];
00359 G4String name = "Unknown";
00360 if (maxTimeRegions[i]) name = maxTimeRegions[i]->GetName();
00361 edm::LogInfo("SimG4CoreApplication") << name << " with pointer "
00362 << maxTimeRegions[i]<<" KE cut off "
00363 << maxTrackTimes[i];
00364 }
00365 }
00366 }
00367
00368 bool StackingAction::isThisVolume(const G4VTouchable* touch,
00369 std::vector<G4LogicalVolume*> & lvs) const {
00370
00371 bool flag = false;
00372 if (lvs.size() > 0 && touch !=0) {
00373 int level = ((touch->GetHistoryDepth())+1);
00374 if (level >= 3) {
00375 unsigned int ii = (unsigned int)(level - 3);
00376 flag = (std::count(lvs.begin(),lvs.end(),(touch->GetVolume(ii)->GetLogicalVolume())) != 0);
00377 }
00378 }
00379 return flag;
00380 }
00381
00382 int StackingAction::isItPrimaryDecayProductOrConversion(const G4Track * aTrack,
00383 const G4Track & mother) const {
00384
00385 int flag = 0;
00386 TrackInformationExtractor extractor;
00387 const TrackInformation & motherInfo(extractor(mother));
00388
00389 if (motherInfo.isPrimary()) {
00390 if (aTrack->GetCreatorProcess()->GetProcessType() == fDecay) flag = 1;
00391 else if (aTrack->GetCreatorProcess()->GetProcessType() == fElectromagnetic &&
00392 aTrack->GetCreatorProcess()->GetProcessSubType() == fGammaConversion) flag = 2;
00393 }
00394 return flag;
00395 }
00396
00397 int StackingAction::isItFromPrimary(const G4Track & mother, int flagIn) const {
00398
00399 int flag = flagIn;
00400 if (flag != 1) {
00401 TrackInformationExtractor extractor;
00402 const TrackInformation & motherInfo(extractor(mother));
00403 if (motherInfo.isPrimary()) flag = 3;
00404 }
00405 return flag;
00406 }
00407
00408 bool StackingAction::isItLongLived(const G4Track * aTrack) const {
00409
00410 bool flag = false;
00411 double time = (aTrack->GetGlobalTime())/nanosecond;
00412 double tofM = maxTrackTime;
00413 if (maxTimeRegions.size() > 0) {
00414 G4Region* reg = aTrack->GetTouchable()->GetVolume(0)->GetLogicalVolume()->GetRegion();
00415 for (unsigned int i=0; i<maxTimeRegions.size(); i++) {
00416 if (reg == maxTimeRegions[i]) {
00417 tofM = maxTrackTimes[i];
00418 break;
00419 }
00420 }
00421 }
00422 if (time > tofM) flag = true;
00423 return flag;
00424 }
00425
00426
1.7.3