Reference Manual

Go to the documentation of this file.
00001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00002 
00003 #include "SimG4Core/GFlash/interface/GflashHadronShowerModel.h"
00004 #include "SimG4Core/Application/interface/SteppingAction.h"
00005 
00006 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
00007 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
00008 #include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h"
00009 #include "SimGeneral/GFlash/interface/GflashKaonMinusShowerProfile.h"
00010 #include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h"
00011 #include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h"
00012 #include "SimGeneral/GFlash/interface/GflashNameSpace.h"
00013 #include "SimGeneral/GFlash/interface/GflashHistogram.h"
00014 #include "SimGeneral/GFlash/interface/GflashHit.h"
00015 
00016 #include "G4FastSimulationManager.hh"
00017 #include "G4TransportationManager.hh"
00018 #include "G4TouchableHandle.hh"
00019 #include "G4VSensitiveDetector.hh"
00020 #include "G4VPhysicalVolume.hh"
00021 #include "G4EventManager.hh"
00022 
00023 #include "G4PionMinus.hh"
00024 #include "G4PionPlus.hh"
00025 #include "G4KaonMinus.hh"
00026 #include "G4KaonPlus.hh"
00027 #include "G4Proton.hh"
00028 #include "G4AntiProton.hh"
00029 #include "G4VProcess.hh"
00030 
00031 #include <vector>
00032 
00033 GflashHadronShowerModel::GflashHadronShowerModel(G4String modelName, G4Region* envelope, edm::ParameterSet parSet)
00034   : G4VFastSimulationModel(modelName, envelope), theParSet(parSet)
00035 {
00036   theWatcherOn = parSet.getParameter<bool>("watcherOn");
00037 
00038   theProfile = new GflashHadronShowerProfile(parSet);
00039   thePiKProfile = new GflashPiKShowerProfile(parSet);
00040   theKaonPlusProfile = new GflashKaonPlusShowerProfile(parSet);
00041   theKaonMinusProfile = new GflashKaonMinusShowerProfile(parSet);
00042   theProtonProfile = new GflashProtonShowerProfile(parSet);
00043   theAntiProtonProfile = new GflashAntiProtonShowerProfile(parSet);
00044   theHisto = GflashHistogram::instance();
00045 
00046   theGflashStep = new G4Step();
00047   theGflashTouchableHandle = new G4TouchableHistory();
00048   theGflashNavigator = new G4Navigator();
00049 
00050 }
00051 
00052 GflashHadronShowerModel::~GflashHadronShowerModel()
00053 {
00054   if(theProfile) delete theProfile;
00055   if(theGflashStep) delete theGflashStep;
00056 }
00057 
00058 G4bool GflashHadronShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
00059 {
00060   return 
00061     &particleType == G4PionMinus::PionMinusDefinition() ||
00062     &particleType == G4PionPlus::PionPlusDefinition() ||
00063     &particleType == G4KaonMinus::KaonMinusDefinition() ||
00064     &particleType == G4KaonPlus::KaonPlusDefinition() ||
00065     &particleType == G4AntiProton::AntiProtonDefinition() ||
00066     &particleType == G4Proton::ProtonDefinition() ;
00067 }
00068 
00069 G4bool GflashHadronShowerModel::ModelTrigger(const G4FastTrack& fastTrack)
00070 {
00071   // ModelTrigger returns false for Gflash Hadronic Shower Model if it is not
00072   // tested from the corresponding wrapper process, GflashHadronWrapperProcess. 
00073   // Temporarily the track status is set to fPostponeToNextEvent at the wrapper
00074   // process before ModelTrigger is really tested for the second time through 
00075   // PostStepGPIL of enviaged parameterization processes.  The better 
00076   // implmentation may be using via G4VUserTrackInformation of each track, which
00077   // requires to modify a geant source code of stepping (G4SteppingManager2)
00078 
00079   G4bool trigger = false;
00080 
00081   // mininum energy cutoff to parameterize
00082   if (fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV < Gflash::energyCutOff) return trigger;
00083 
00084   // check whether this is called from the normal GPIL or the wrapper process GPIL
00085   if(fastTrack.GetPrimaryTrack()->GetTrackStatus() == fPostponeToNextEvent ) {
00086 
00087     // Shower pameterization start at the first inelastic interaction point
00088     G4bool isInelastic  = isFirstInelasticInteraction(fastTrack);
00089     
00090     // Other conditions
00091     if(isInelastic) {
00092       trigger = (!excludeDetectorRegion(fastTrack));
00093     }
00094   }
00095 
00096   return trigger;
00097 
00098 }
00099 
00100 void GflashHadronShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
00101 {
00102   // kill the particle
00103   fastStep.KillPrimaryTrack();
00104   fastStep.ProposePrimaryTrackPathLength(0.0);
00105 
00106   // parameterize energy depostion by the particle type
00107   G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
00108   
00109   theProfile = thePiKProfile;
00110   if(particleType == G4KaonMinus::KaonMinusDefinition()) theProfile = theKaonMinusProfile;
00111   else if(particleType == G4KaonPlus::KaonPlusDefinition()) theProfile = theKaonPlusProfile;
00112   else if(particleType == G4AntiProton::AntiProtonDefinition()) theProfile = theAntiProtonProfile;
00113   else if(particleType == G4Proton::ProtonDefinition()) theProfile = theProtonProfile;
00114 
00115   //input variables for GflashHadronShowerProfile
00116   G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV;
00117   G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
00118   G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
00119   G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition()/cm;
00120   G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum()/GeV;
00121   G4int showerType = Gflash::findShowerType(position);
00122 
00123   theProfile->initialize(showerType,energy,globalTime,charge,position,momentum);
00124   theProfile->loadParameters();
00125   theProfile->hadronicParameterization();
00126 
00127   // make hits
00128   makeHits(fastTrack);
00129 
00130 }
00131 
00132 void GflashHadronShowerModel::makeHits(const G4FastTrack& fastTrack) {
00133 
00134   std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
00135 
00136   theGflashStep->SetTrack(const_cast<G4Track*>(fastTrack.GetPrimaryTrack()));
00137   theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(const_cast<G4VProcess*>
00138     (fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
00139   theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
00140 
00141   std::vector<GflashHit>::const_iterator spotIter    = gflashHitList.begin();
00142   std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
00143 
00144   for( ; spotIter != spotIterEnd; spotIter++){
00145 
00146     theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(spotIter->getPosition(),G4ThreeVector(0,0,0),
00147                                                                   theGflashTouchableHandle, false);
00148     updateGflashStep(spotIter->getPosition(),spotIter->getTime());
00149 
00150     // if there is a watcher defined in a job and the flag is turned on
00151     if(theWatcherOn) {
00152       theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy());
00153       SteppingAction* userSteppingAction = (SteppingAction*) G4EventManager::GetEventManager()->GetUserSteppingAction();
00154       userSteppingAction->m_g4StepSignal(theGflashStep);
00155     }
00156 
00157     G4VPhysicalVolume* aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
00158     if( aCurrentVolume == 0 ) continue;
00159 
00160     G4LogicalVolume* lv = aCurrentVolume->GetLogicalVolume();
00161     if(lv->GetRegion()->GetName() != "CaloRegion") continue;
00162           
00163     theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
00164     G4VSensitiveDetector* aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
00165 
00166     if( aSensitive == 0 ) continue;
00167     
00168     G4String nameCalor = aCurrentVolume->GetName();
00169     nameCalor.assign(nameCalor,0,2);
00170     G4double samplingWeight = 1.0; 
00171     if(nameCalor == "HB" ) {
00172       samplingWeight = Gflash::scaleSensitiveHB;
00173     }
00174     else if(nameCalor=="HE" || nameCalor == "HT") {
00175       samplingWeight = Gflash::scaleSensitiveHE;
00176     }
00177     theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy()*samplingWeight);
00178 
00179     aSensitive->Hit(theGflashStep);
00180 
00181   }
00182 }
00183 
00184 void GflashHadronShowerModel::updateGflashStep(G4ThreeVector spotPosition, G4double timeGlobal)
00185 {
00186   theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
00187   theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
00188   theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
00189   theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
00190 }
00191 
00192 G4bool GflashHadronShowerModel::isFirstInelasticInteraction(const G4FastTrack& fastTrack)
00193 {
00194   G4bool isFirst=false;
00195 
00196   G4StepPoint* preStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint();
00197   G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
00198 
00199   G4String procName = postStep->GetProcessDefinedStep()->GetProcessName();  
00200   G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
00201 
00202   //@@@ this part is still temporary and the cut for the variable ratio should be optimized later
00203 
00204   if((particleType == G4PionPlus::PionPlusDefinition() && procName == "WrappedPionPlusInelastic") || 
00205      (particleType == G4PionMinus::PionMinusDefinition() && procName == "WrappedPionMinusInelastic") ||
00206      (particleType == G4KaonPlus::KaonPlusDefinition() && procName == "WrappedKaonPlusInelastic") ||
00207      (particleType == G4KaonMinus::KaonMinusDefinition() && procName == "WrappedKaonMinusInelastic") ||
00208      (particleType == G4AntiProton::AntiProtonDefinition() && procName == "WrappedAntiProtonInelastic") ||
00209      (particleType == G4Proton::ProtonDefinition() && procName == "WrappedProtonInelastic")) {
00210 
00211     //skip to the second interaction if the first inelastic is a quasi-elastic like interaction
00212     //@@@ the cut may be optimized later
00213 
00214     const G4TrackVector* fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary();
00215     G4double leadingEnergy = 0.0;
00216 
00217     //loop over 'all' secondaries including those produced by continuous processes.
00218     //@@@may require an additional condition only for hadron interaction with the process name,
00219     //but it will not change the result anyway
00220 
00221     for (unsigned int isec = 0 ; isec < fSecondaryVector->size() ; isec++) {
00222       G4Track* fSecondaryTrack = (*fSecondaryVector)[isec];
00223       G4double secondaryEnergy = fSecondaryTrack->GetKineticEnergy();
00224 
00225       if(secondaryEnergy > leadingEnergy ) {
00226         leadingEnergy = secondaryEnergy;
00227       }
00228     }
00229 
00230     if((preStep->GetTotalEnergy()!=0) && 
00231        (leadingEnergy/preStep->GetTotalEnergy() < Gflash::QuasiElasticLike)) isFirst = true;
00232 
00233     //Fill debugging histograms and check information on secondaries -
00234     //remove after final implimentation
00235 
00236     if(theHisto->getStoreFlag()) {
00237       theHisto->preStepPosition->Fill(preStep->GetPosition().getRho()/cm);
00238       theHisto->postStepPosition->Fill(postStep->GetPosition().getRho()/cm);
00239       theHisto->deltaStep->Fill((postStep->GetPosition() - preStep->GetPosition()).getRho()/cm);
00240       theHisto->kineticEnergy->Fill(fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV);
00241       theHisto->energyLoss->Fill(fabs(fastTrack.GetPrimaryTrack()->GetStep()->GetDeltaEnergy()/GeV));
00242       theHisto->energyRatio->Fill(leadingEnergy/preStep->GetTotalEnergy());
00243     }
00244 
00245  }
00246   return isFirst;
00247 } 
00248 
00249 G4bool GflashHadronShowerModel::excludeDetectorRegion(const G4FastTrack& fastTrack)
00250 {
00251   G4bool isExcluded=false;
00252   int verbosity = theParSet.getUntrackedParameter<int>("Verbosity");
00253   
00254   //exclude regions where geometry are complicated
00255   //+- one supermodule around the EB/EE boundary: 1.479 +- 0.0174*5
00256   G4double eta =   fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity() ;
00257   if(std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566) {
00258     if(verbosity>0) {
00259        edm::LogInfo("SimGeneralGFlash") << "GflashHadronShowerModel: excluding region of eta = " << eta;
00260     }
00261     return true;  
00262   }
00263   else {
00264     G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
00265 
00266     Gflash::CalorimeterNumber kCalor = Gflash::getCalorimeterNumber(postStep->GetPosition()/cm);
00267     G4double distOut = 9999.0;
00268 
00269     //exclude the region where the shower starting point is inside the preshower
00270     if( std::fabs(eta) > Gflash::EtaMin[Gflash::kENCA] &&
00271         std::fabs((postStep->GetPosition()).getZ()/cm) < Gflash::Zmin[Gflash::kENCA]) {
00272       return true;
00273     }
00274 
00275     //<---the shower starting point is always inside envelopes
00276     //@@@exclude the region where the shower starting point is too close to the end of
00277     //the hadronic envelopes (may need to be optimized further!)
00278     //@@@if we extend parameterization including Magnet/HO, we need to change this strategy
00279     if(kCalor == Gflash::kHB) {
00280       distOut =  Gflash::Rmax[Gflash::kHB] - postStep->GetPosition().getRho()/cm;
00281       if (distOut < Gflash::MinDistanceToOut ) isExcluded = true;
00282     }
00283     else if(kCalor == Gflash::kHE) {
00284       distOut =  Gflash::Zmax[Gflash::kHE] - std::fabs(postStep->GetPosition().getZ()/cm);
00285       if (distOut < Gflash::MinDistanceToOut ) isExcluded = true;
00286     }
00287 
00288     //@@@remove this print statement later
00289     if(isExcluded && verbosity > 0) {
00290       std::cout << "GflashHadronShowerModel: skipping kCalor = " << kCalor << 
00291         " DistanceToOut " << distOut << " from (" <<  (postStep->GetPosition()).getRho()/cm << 
00292         ":" << (postStep->GetPosition()).getZ()/cm << ") of KE = " << fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV << std::endl;
00293     }
00294   }
00295 
00296   return isExcluded;
00297 }
/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4Core/GFlash/src/GflashHadronShowerModel.cc
