Go to the documentation of this file.00001 #include "SimG4Core/Application/interface/SteppingAction.h"
00002 #include "SimG4Core/Application/interface/EventAction.h"
00003
00004 #include "G4LogicalVolumeStore.hh"
00005 #include "G4ParticleTable.hh"
00006 #include "G4PhysicalVolumeStore.hh"
00007 #include "G4RegionStore.hh"
00008 #include "G4Track.hh"
00009 #include "G4UnitsTable.hh"
00010
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012
00013 SteppingAction::SteppingAction(EventAction* e,const edm::ParameterSet & p)
00014 : eventAction_(e), initialized(false), tracker(0), calo(0) {
00015
00016 killBeamPipe = (p.getParameter<bool>("KillBeamPipe"));
00017 theCriticalEnergyForVacuum = (p.getParameter<double>("CriticalEnergyForVacuum")*MeV);
00018 theCriticalDensity = (p.getParameter<double>("CriticalDensity")*g/cm3);
00019 maxTrackTime = (p.getParameter<double>("MaxTrackTime")*ns);
00020 maxTrackTimes = (p.getParameter<std::vector<double> >("MaxTrackTimes"));
00021 maxTimeNames = (p.getParameter<std::vector<std::string> >("MaxTimeNames"));
00022 ekinMins = (p.getParameter<std::vector<double> >("EkinThresholds"));
00023 ekinNames = (p.getParameter<std::vector<std::string> >("EkinNames"));
00024 ekinParticles = (p.getParameter<std::vector<std::string> >("EkinParticles"));
00025 verbose = (p.getUntrackedParameter<int>("Verbosity",0));
00026
00027 edm::LogInfo("SimG4CoreApplication") << "SteppingAction:: KillBeamPipe = "
00028 << killBeamPipe << " CriticalDensity = "
00029 << theCriticalDensity << " g/cc;"
00030 << " CriticalEnergyForVacuum = "
00031 << theCriticalEnergyForVacuum << " Mev;"
00032 << " MaxTrackTime = " << maxTrackTime
00033 << " ns";
00034 for (unsigned int i=0; i<maxTrackTimes.size(); i++) {
00035 maxTrackTimes[i] *= ns;
00036 edm::LogInfo("SimG4CoreApplication") << "SteppingAction::MaxTrackTime for "
00037 << maxTimeNames[i] << " is "
00038 << maxTrackTimes[i];
00039 }
00040 edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Kill following "
00041 << ekinParticles.size()
00042 << " particles in " << ekinNames.size()
00043 << " volumes";
00044 for (unsigned int i=0; i<ekinParticles.size(); i++) {
00045 ekinMins[i] *= GeV;
00046 edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Particle " << i
00047 << " " << ekinParticles[i]
00048 << " (Threshold = " << ekinMins[i]
00049 << " MeV)";
00050 }
00051 for (unsigned int i=0; i<ekinNames.size(); i++)
00052 edm::LogInfo("SimG4CoreApplication") << "SteppingAction::Volume[" << i
00053 << "] = " << ekinNames[i];
00054 }
00055
00056 SteppingAction::~SteppingAction() {}
00057
00058 void SteppingAction::UserSteppingAction(const G4Step * aStep) {
00059 if (!initialized) initialized = initPointer();
00060 m_g4StepSignal(aStep);
00061
00062 if (killBeamPipe) {
00063 catchLowEnergyInVacuumHere(aStep);
00064 catchLowEnergyInVacuumNext(aStep);
00065 }
00066
00067 if (aStep->GetPostStepPoint()->GetPhysicalVolume() != 0) {
00068 bool ok = catchLongLived(aStep);
00069 if (ok) ok = catchLongLived(aStep);
00070
00071 ok = (isThisVolume(aStep->GetPreStepPoint()->GetTouchable(),tracker)&&
00072 isThisVolume(aStep->GetPostStepPoint()->GetTouchable(),calo));
00073 if (ok) {
00074
00075 math::XYZVectorD pos((aStep->GetPreStepPoint()->GetPosition()).x(),
00076 (aStep->GetPreStepPoint()->GetPosition()).y(),
00077 (aStep->GetPreStepPoint()->GetPosition()).z());
00078
00079 math::XYZTLorentzVectorD mom((aStep->GetPreStepPoint()->GetMomentum()).x(),
00080 (aStep->GetPreStepPoint()->GetMomentum()).y(),
00081 (aStep->GetPreStepPoint()->GetMomentum()).z(),
00082 aStep->GetPreStepPoint()->GetTotalEnergy());
00083
00084 uint32_t id = aStep->GetTrack()->GetTrackID();
00085
00086 std::pair<math::XYZVectorD,math::XYZTLorentzVectorD> p(pos,mom);
00087 eventAction_->addTkCaloStateInfo(id,p);
00088 }
00089 }
00090 }
00091
00092 void SteppingAction::catchLowEnergyInVacuumHere(const G4Step * aStep) {
00093 G4Track * theTrack = aStep->GetTrack();
00094 double theKenergy = theTrack->GetKineticEnergy();
00095 if (theTrack->GetVolume()!=0) {
00096 double density = theTrack->GetVolume()->GetLogicalVolume()->GetMaterial()->GetDensity();
00097 if (theKenergy <= theCriticalEnergyForVacuum && theKenergy > 0.0 &&
00098 density <= theCriticalDensity && theTrack->GetDefinition()->GetPDGCharge() != 0 &&
00099 theTrack->GetTrackStatus() != fStopAndKill) {
00100 if (verbose>1)
00101 edm::LogInfo("SimG4CoreApplication")
00102 << " SteppingAction: LoopCatchSteppingAction:catchLowEnergyInVacuumHere: "
00103 << " Track from " << theTrack->GetDefinition()->GetParticleName()
00104 << " of kinetic energy " << theKenergy/MeV << " MeV "
00105 << " killed in " << theTrack->GetVolume()->GetLogicalVolume()->GetName()
00106 << " of density " << density/(g/cm3) << " g/cm3" ;
00107 theTrack->SetTrackStatus(fStopAndKill);
00108 }
00109 }
00110 }
00111
00112 void SteppingAction::catchLowEnergyInVacuumNext(const G4Step * aStep) {
00113 G4Track * theTrack = aStep->GetTrack();
00114 double theKenergy = theTrack->GetKineticEnergy();
00115 if (theTrack->GetNextVolume()) {
00116 double density = theTrack->GetNextVolume()->GetLogicalVolume()->GetMaterial()->GetDensity();
00117 if (theKenergy <= theCriticalEnergyForVacuum && theKenergy > 0.0 &&
00118 density <= theCriticalDensity && theTrack->GetDefinition()->GetPDGCharge() != 0 &&
00119 theTrack->GetTrackStatus() != fStopAndKill) {
00120 if (verbose>1)
00121 edm::LogInfo("SimG4CoreApplication")
00122 << " SteppingAction: LoopCatchSteppingAction::catchLowEnergyInVacuumNext: "
00123 << " Track from " << theTrack->GetDefinition()->GetParticleName()
00124 << " of kinetic energy " << theKenergy/MeV << " MeV "
00125 << " stopped in " << theTrack->GetVolume()->GetLogicalVolume()->GetName()
00126 << " before going into "<< theTrack->GetNextVolume()->GetLogicalVolume()->GetName()
00127 << " of density " << density/(g/cm3) << " g/cm3" ;
00128 theTrack->SetTrackStatus(fStopButAlive);
00129 }
00130 }
00131 }
00132
00133 bool SteppingAction::catchLongLived(const G4Step * aStep) {
00134
00135 bool flag = true;
00136 double time = (aStep->GetPostStepPoint()->GetGlobalTime())/nanosecond;
00137 double tofM = maxTrackTime;
00138 if (maxTimeRegions.size() > 0) {
00139 G4Region* reg = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion();
00140 for (unsigned int i=0; i<maxTimeRegions.size(); i++) {
00141 if (reg == maxTimeRegions[i]) {
00142 tofM = maxTrackTimes[i];
00143 break;
00144 }
00145 }
00146 }
00147 if (time > tofM) {
00148 killTrack(aStep);
00149 flag = false;
00150 }
00151 return flag;
00152 }
00153
00154 bool SteppingAction::killLowEnergy(const G4Step * aStep) {
00155
00156 bool ok = true;
00157 if (ekinVolumes.size() > 0) {
00158 bool flag = false;
00159 G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
00160 for (unsigned int i=0; i<ekinVolumes.size(); i++) {
00161 if (lv == ekinVolumes[i]) {
00162 flag = true;
00163 break;
00164 }
00165 }
00166 if (flag) {
00167 G4Track * track = aStep->GetTrack();
00168 double ekin = track->GetKineticEnergy();
00169 double ekinM = 0;
00170 int pCode = track->GetDefinition()->GetPDGEncoding();
00171 for (unsigned int i=0; i<ekinPDG.size(); i++) {
00172 if (pCode == ekinPDG[i]) {
00173 ekinM = ekinMins[i];
00174 break;
00175 }
00176 }
00177 if (ekin > ekinM) {
00178 killTrack(aStep);
00179 ok = false;
00180 }
00181 }
00182 }
00183 return ok;
00184 }
00185
00186 bool SteppingAction::initPointer() {
00187
00188 bool flag = true;
00189 const G4PhysicalVolumeStore * pvs = G4PhysicalVolumeStore::GetInstance();
00190 if (pvs) {
00191 std::vector<G4VPhysicalVolume*>::const_iterator pvcite;
00192 for (pvcite = pvs->begin(); pvcite != pvs->end(); pvcite++) {
00193 if ((*pvcite)->GetName() == "Tracker") tracker = (*pvcite);
00194 if ((*pvcite)->GetName() == "CALO") calo = (*pvcite);
00195 if (tracker && calo) break;
00196 }
00197 edm::LogInfo("SimG4CoreApplication") << "Pointer for Tracker " << tracker
00198 << " and for Calo " << calo;
00199 if (tracker) LogDebug("SimG4CoreApplication") << "Tracker vol name "
00200 << tracker->GetName();
00201 if (calo) LogDebug("SimG4CoreApplication") << "Calorimeter vol name "
00202 << calo->GetName();
00203 } else {
00204 flag = false;
00205 }
00206
00207 const G4LogicalVolumeStore * lvs = G4LogicalVolumeStore::GetInstance();
00208 unsigned int num = ekinNames.size();
00209 if (num > 0) {
00210 if (lvs) {
00211 std::vector<G4LogicalVolume*>::const_iterator lvcite;
00212 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
00213 for (unsigned int i=0; i<num; i++) {
00214 if ((*lvcite)->GetName() == (G4String)(ekinNames[i])) {
00215 ekinVolumes.push_back(*lvcite);
00216 break;
00217 }
00218 }
00219 if (ekinVolumes.size() == num) break;
00220 }
00221 }
00222 if (ekinVolumes.size() != num) flag = false;
00223 for (unsigned int i=0; i<ekinVolumes.size(); i++) {
00224 edm::LogInfo("SimG4CoreApplication") << ekinVolumes[i]->GetName()
00225 <<" with pointer " <<ekinVolumes[i];
00226 }
00227 }
00228
00229 G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
00230 G4String particleName;
00231 for (unsigned int i=0; i<ekinParticles.size(); i++) {
00232 int pdg = theParticleTable->FindParticle(particleName=ekinParticles[i])->GetPDGEncoding();
00233 ekinPDG.push_back(pdg);
00234 edm::LogInfo("SimG4CoreApplication") << "Particle " << ekinParticles[i]
00235 << " with code " << ekinPDG[i]
00236 << " and KE cut off " << ekinMins[i];
00237 }
00238 if (!flag) edm::LogInfo("SimG4CoreApplication") << "SteppingAction fails to"
00239 << " initialize some the "
00240 << "LV pointers correctly";
00241
00242 const G4RegionStore * rs = G4RegionStore::GetInstance();
00243 num = maxTimeNames.size();
00244 if (num > 0) {
00245 std::vector<double> tofs;
00246 if (rs) {
00247 std::vector<G4Region*>::const_iterator rcite;
00248 for (rcite = rs->begin(); rcite != rs->end(); rcite++) {
00249 for (unsigned int i=0; i<num; i++) {
00250 if ((*rcite)->GetName() == (G4String)(maxTimeNames[i])) {
00251 maxTimeRegions.push_back(*rcite);
00252 tofs.push_back(maxTrackTimes[i]);
00253 break;
00254 }
00255 }
00256 if (tofs.size() == num) break;
00257 }
00258 }
00259 for (unsigned int i=0; i<tofs.size(); i++) {
00260 maxTrackTimes[i] = tofs[i];
00261 G4String name = "Unknown";
00262 if (maxTimeRegions[i]) name = maxTimeRegions[i]->GetName();
00263 edm::LogInfo("SimG4CoreApplication") << name << " with pointer "
00264 << maxTimeRegions[i]<<" KE cut off "
00265 << maxTrackTimes[i];
00266 }
00267 if (tofs.size() != num)
00268 edm::LogInfo("SimG4CoreApplication") << "SteppingAction fails to "
00269 << "initialize some the region "
00270 << "pointers correctly";
00271 }
00272 return true;
00273 }
00274
00275 bool SteppingAction::isThisVolume(const G4VTouchable* touch,
00276 G4VPhysicalVolume* pv) {
00277
00278 int level = ((touch->GetHistoryDepth())+1);
00279 if (level > 0 && level >= 3) {
00280 int ii = level - 3;
00281 return (touch->GetVolume(ii) == pv);
00282 }
00283 return false;
00284 }
00285
00286 void SteppingAction::killTrack(const G4Step * aStep) {
00287
00288 aStep->GetTrack()->SetTrackStatus(fStopAndKill);
00289 G4TrackVector tv = *(aStep->GetSecondary());
00290 for (unsigned int kk=0; kk<tv.size(); kk++) {
00291 if (tv[kk]->GetVolume() == aStep->GetPreStepPoint()->GetPhysicalVolume())
00292 tv[kk]->SetTrackStatus(fStopAndKill);
00293 }
00294 LogDebug("SimG4CoreApplication") << "SteppingAction: Kills track "
00295 << aStep->GetTrack()->GetTrackID() << " ("
00296 << aStep->GetTrack()->GetDefinition()->GetParticleName()
00297 << ") at " << aStep->GetPostStepPoint()->GetGlobalTime()/nanosecond
00298 << " ns in " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName()
00299 << " from region " << aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetRegion()->GetName();
00300 }
1.7.3