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