LBNEStackingAction Class Reference

#include <LBNEStackingAction.hh>

Inheritance diagram for LBNEStackingAction:

Public Member Functions
	LBNEStackingAction ()
virtual	~LBNEStackingAction ()
virtual G4ClassificationOfNewTrack	ClassifyNewTrack (const G4Track *aTrack)
virtual void	NewStage ()
virtual void	PrepareNewEvent ()
void	KillEMParticles (G4ClassificationOfNewTrack &classification, const G4Track *aTrack)
void	KillltZeroPzParticles (G4ClassificationOfNewTrack &classification, const G4Track *aTrack)
void	KillThresholdParticles (G4ClassificationOfNewTrack &classification, const G4Track *aTrack)
void	KillUnimportantParticles (G4ClassificationOfNewTrack &classification, const G4Track *aTrack)
void	UnKillMuonAlcoveParticles (G4ClassificationOfNewTrack &classification, const G4Track *aTrack)
void	OpenHadronAtVertex () const
void	SetDoImportanceWeightSelection (bool t) const
bool	GetDoImportanceWeightSelection () const
void	SetStackingKillingThreshold (double v) const
double	GetStackingKillingThreshold () const

Private Attributes
LBNERunManager *	pRunManager
double	fStackingKillingThreshold
std::ofstream	fOutHadronAtVectex
std::ofstream	fOutHadronAtVectexCl
bool	doImportanceWeightSelection

Detailed Description

Definition at line 15 of file LBNEStackingAction.hh.

Constructor & Destructor Documentation

LBNEStackingAction::LBNEStackingAction

(

)

Definition at line 24 of file LBNEStackingAction.cc.

{ 
  pRunManager=(LBNERunManager*)LBNERunManager::GetRunManager();
  std::cout << "LBNEStackingAction Constructor Called." << std::endl;
  doImportanceWeightSelection = true; // Per G4Numi tradition.. 
  fStackingKillingThreshold = 0.050*CLHEP::GeV;
}

LBNEStackingAction::~LBNEStackingAction ( )

virtual

Definition at line 33 of file LBNEStackingAction.cc.

{
   int verboseLevel = G4EventManager::GetEventManager()->GetTrackingManager()->GetVerboseLevel();
   if(verboseLevel > 0)
   {
      std::cout << "LBNEStackingAction Destructor Called." << std::endl;
   }
}

Member Function Documentation

G4ClassificationOfNewTrack LBNEStackingAction::ClassifyNewTrack ( const G4Track *  aTrack )

virtual

Definition at line 85 of file LBNEStackingAction.cc.

{
   //
   // v410.2:  Rare bug somewhere, we have track with bad parameters, we do not want them.. 
   //
   G4ThreeVector posTmp = aTrack->GetPosition();
   G4ThreeVector momTmp = aTrack->GetMomentum();
   for (size_t k=0; k != 3; k++) {
     if (std::isnan(posTmp[k])) {
       std::cerr << " LBNEStackingAction::ClassifyNewTrack, position has NaN, rejected " << std::endl;
       return fKill; 
     }
     if (std::isnan(momTmp[k])) {
       std::cerr << " LBNEStackingAction::ClassifyNewTrack, momentum has NaN, rejected " << std::endl;
       return fKill; 
     }
   }
   if (std::isnan(aTrack->GetTotalEnergy()) || std::isnan(aTrack->GetKineticEnergy())) {
       std::cerr << " LBNEStackingAction::ClassifyNewTrack, Energy is NaN, rejected " << std::endl;
       return fKill; 
   }
   
   if (fOutHadronAtVectex.is_open()) {
      int myPid = aTrack->GetParticleDefinition()->GetPDGEncoding();
      if ((std::abs(myPid) == 211) || (myPid == 2212) || 
          (std::abs(myPid) == 321) || (std::abs(myPid) == 311) ||
          (std::abs(myPid) == 130) || (std::abs(myPid) == 310)) {     
        double pTotSq = 0; for (int k=0; k != 3; k++) pTotSq += aTrack->GetMomentum()[k]*aTrack->GetMomentum()[k];
        double pTot = std::sqrt(pTotSq)/CLHEP::GeV;
//        if ((myPid == 2212) || ((pTot > fStackingKillingThreshold) && (pTot < 15.))) { // take all the proton, to check we are on target, and the relevant pions 
        if ((myPid == 2212) || ((pTot > 0.001) && (pTot < 50.))) { // take all the proton, to check we are on target, and the relevant pions 
          int evtId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
          fOutHadronAtVectex << " " << evtId;
          for (int k=0; k != 3; k++) fOutHadronAtVectex << " " << aTrack->GetPosition()[k];
          for (int k=0; k != 3; k++) fOutHadronAtVectex << " " << (aTrack->GetMomentum()[k])/CLHEP::GeV;
          fOutHadronAtVectex  << " " << myPid << std::endl; 
        } 
      
     }
   }
//   const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(aTrack->GetUserInformation());
//   std::cerr << " G4ClassificationOfNewTrack LBNEStackingAction::ClassifyNewTrack, ptr track " 
//             << (void*) aTrack << " info " <<  (void*) oldTrInfo << std::endl;
//   if (oldTrInfo != 0) std::cerr << " ... PFFlag " << oldTrInfo->GetPFFlag() << std::endl;

  G4ClassificationOfNewTrack classification = fUrgent;
  G4ParticleDefinition * particleType = aTrack->GetDefinition();

   int verboseLevel = G4EventManager::GetEventManager()->GetTrackingManager()->GetVerboseLevel();
  if(verboseLevel > 0)
  {
     G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
     std::cout << "Event " << evtno << ": LBNEStackingAction::ClassifyNewTrack for Particle " 
               << "\"" << particleType->GetParticleName() << "\"" << " Called." << std::endl;
  }
   if (LBNEQuickPiToNuVect::Instance()->doIt() && (std::abs(particleType->GetPDGEncoding()) == 211)) {
      int sign = (particleType->GetPDGEncoding() > 0) ? 1 : -1;
      LBNEQuickPiToNuVect::Instance()->addPion(aTrack->GetTrackID(), sign, aTrack->GetTotalEnergy(), aTrack->GetMomentum());
   }
   const LBNERunAction *runUserAction =  reinterpret_cast<const LBNERunAction *>(G4RunManager::GetRunManager()->GetUserRunAction());
   // 
   // If we are in the target, and tally the energy deposition, skip the re-classification done below..  
   // 
   if (runUserAction->GetDoComputeEDepInGraphite() && (aTrack->GetVolume() != 0)) {
     const std::string vName = aTrack->GetVolume()->GetName();
     if(vName.find("Target") != std::string::npos) return classification;
   }
   if (runUserAction->GetDoComputeEDepInArgonGas()) {
     if (aTrack->GetPosition()[2] < 9900.)  return classification; // rough  Z- cut
   }
   // 
   // Now do the reclassification, for CPU speed optimization.. 
   //
   LBNEStackingAction::KillEMParticles(classification, aTrack);
     //LBNEStackingAction::KillltZeroPzParticles(classification, aTrack);
   LBNEStackingAction::KillThresholdParticles(classification, aTrack);
//   LBNEStackingAction::KillOutOfWorldParticles(classification, aTrack); Obsolete in G4 
   LBNEStackingAction::KillUnimportantParticles(classification, aTrack);
//    
   if (fOutHadronAtVectexCl.is_open() && (classification == fUrgent)) {
      int myPid = aTrack->GetParticleDefinition()->GetPDGEncoding();
      if (std::abs(myPid) == 211) {     
        double pTotSq = 0; for (int k=0; k != 3; k++) pTotSq += aTrack->GetMomentum()[k]*aTrack->GetMomentum()[k];
        double pTot = std::sqrt(pTotSq)/CLHEP::GeV;
        if ((myPid == 2212) || ((pTot > fStackingKillingThreshold) && (pTot < 15.))) { // take all the proton, to check we are on target, and the relevant pions 
          int evtId = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
          fOutHadronAtVectexCl << " " << evtId;
          for (int k=0; k != 3; k++) fOutHadronAtVectexCl << " " << aTrack->GetPosition()[k];
          for (int k=0; k != 3; k++) fOutHadronAtVectexCl << " " << (aTrack->GetMomentum()[k])/CLHEP::GeV;
          fOutHadronAtVectexCl  << " " << myPid << std::endl; 
        } 
      
     }
   }

   LBNEStackingAction::UnKillMuonAlcoveParticles(classification,aTrack);
// 
// specific checks for impotance weight. 
//
//   if (std::abs(aTrack->GetParticleDefinition()->GetPDGEncoding()) == 211) { 
//      LBNETrackInformation *trInfo = reinterpret_cast<LBNETrackInformation*> (aTrack->GetUserInformation());  
//      if (trInfo != 0) std::cerr << "LBNEStackingAction::classifyNewTrack, Track ID " << aTrack->GetTrackID() 
//              << " Imp Weight " << trInfo->GetNImpWt() << std::endl;
//      else std::cerr << "LBNEStackingAction::classifyNewTrack, Track ID " << aTrack->GetTrackID() 
//              << " has no weight ... " << std::endl;
//   }
      
  return classification;
}

bool LBNEStackingAction::GetDoImportanceWeightSelection ( ) const

inline

Definition at line 52 of file LBNEStackingAction.hh.

{return doImportanceWeightSelection; }

double LBNEStackingAction::GetStackingKillingThreshold ( ) const

inline

Definition at line 54 of file LBNEStackingAction.hh.

{ return fStackingKillingThreshold;}

void LBNEStackingAction::KillEMParticles	(	G4ClassificationOfNewTrack &	classification,
		const G4Track *	aTrack
	)

Definition at line 244 of file LBNEStackingAction.cc.

{
    const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
    if (runAct->GetDoComputeEDepInHorns() || runAct->GetDoComputeEDepInGraphite ()
                                          || runAct->GetDoComputeEDepInArgonGas()) return;
   G4ParticleDefinition * particleType = aTrack->GetDefinition();
   // Discard Gammas, Electrons, ...
   if ( (particleType==G4Gamma::GammaDefinition()       ||
        particleType==G4Electron::ElectronDefinition()  ||
        particleType==G4Positron::PositronDefinition() ) &&
        (classification != fKill) )
   {classification = fKill;}
   
}

void LBNEStackingAction::KillltZeroPzParticles	(	G4ClassificationOfNewTrack &	classification,
		const G4Track *	aTrack
	)

Definition at line 261 of file LBNEStackingAction.cc.

{
   const LBNERunAction *runUserAction =  reinterpret_cast<const LBNERunAction *>(G4RunManager::GetRunManager()->GetUserRunAction());
   if (runUserAction->GetDoComputeEDepInHorns()) return;
   if (runUserAction->GetDoComputeEDepInGraphite() && (aTrack->GetMaterial() != 0)) {
     const std::string tgt("Target");
     const std::string tg = aTrack->GetMaterial()->GetName();
     if (tg == tgt) return; // Don't loose any track, computing E-Loss...  
   }
  //Discard particles with pz<0
   G4ThreeVector momentum=aTrack->GetMomentumDirection();
   if (momentum[2]<0 && (classification != fKill)) {classification = fKill;}
   
}

void LBNEStackingAction::KillThresholdParticles	(	G4ClassificationOfNewTrack &	classification,
		const G4Track *	aTrack
	)

Definition at line 278 of file LBNEStackingAction.cc.

{

    const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
    if (runAct->GetDoComputeEDepInHorns() || runAct->GetDoComputeEDepInGraphite ()
                                          || runAct->GetDoComputeEDepInArgonGas()) return;
    G4ParticleDefinition * particleType = aTrack->GetDefinition();
   //Discard particles with kinetic energy < KillTracking Threshold.GeV (that are not neutrinos)
   if ((particleType!=G4NeutrinoE::NeutrinoEDefinition())&&
       (particleType!=G4NeutrinoMu::NeutrinoMuDefinition())&&
       (particleType!=G4NeutrinoTau::NeutrinoTauDefinition())&&
       (particleType!=G4AntiNeutrinoE::AntiNeutrinoEDefinition())&&
       (particleType!=G4AntiNeutrinoMu::AntiNeutrinoMuDefinition())&&
       (particleType!=G4AntiNeutrinoTau::AntiNeutrinoTauDefinition()))
   {
      G4double energy = aTrack->GetKineticEnergy();
//
// We hard code this cut, as the LBNEDataInput obsolete.. Could become a data card 
//      
      if ((energy < fStackingKillingThreshold ) && (classification != fKill))
          {classification = fKill;} 
   }
   
   
}

void LBNEStackingAction::KillUnimportantParticles	(	G4ClassificationOfNewTrack &	classification,
		const G4Track *	aTrack
	)

Definition at line 337 of file LBNEStackingAction.cc.

{
      const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
    if (runAct->GetDoComputeEDepInHorns() || runAct->GetDoComputeEDepInGraphite ()
                                          || runAct->GetDoComputeEDepInArgonGas()) return;
 //If importance weighting might be on:
 if ((classification != fKill))
  {
    LBNETrackInformation* trackInfo=(LBNETrackInformation*)(aTrack->GetUserInformation());  
    if (doImportanceWeightSelection) { 
      if (trackInfo!=0) 
        {
           G4double Nimpweight=LBNEImpWeight::CalculateImpWeight(aTrack);
//         std::cerr << " LBNEStackingAction::KillUnimportantParticles, E " 
//                   << aTrack->GetTotalEnergy() << " Imp weight " <<  Nimpweight << std::endl;
          if(Nimpweight==0.) { classification = fKill; } 
           else               { trackInfo->SetNImpWt(Nimpweight); }
        }
    } else {
    // We do not change the classification.. 
      if (trackInfo!=0) trackInfo->SetNImpWt(1.0);
    }
  } 
}

void LBNEStackingAction::NewStage ( )

virtual

Definition at line 195 of file LBNEStackingAction.cc.

{
   int verboseLevel = G4EventManager::GetEventManager()->GetTrackingManager()->GetVerboseLevel();
   if(verboseLevel > 1)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno << ": LBNEStackingAction::NewStage Called." << std::endl;
   }


  // stackManager->ReClassify();
  //  return;
}

void LBNEStackingAction::OpenHadronAtVertex

(

)

const

Definition at line 42 of file LBNEStackingAction.cc.

                                                  { // Pseudo const.. 
   if (fOutHadronAtVectex.is_open()) {
     std::cerr << " From LBNEStackingAction::OpenHadronAtVertex, file already opened, fatal... " << std::endl;
     exit(2);
   }
   // Set a file name that can be unique.. Particularly, tailored for the use on FermiGrid
   // But not mandated. 
   int clusterNum = 0;
   int procNum =0;
   const char *clusterEnv = getenv("CLUSTER");
   if (clusterEnv != 0) clusterNum = atoi(clusterEnv);
   const char *procEnv = getenv("PROCESS");
   if (procEnv != 0) procNum = atoi(procEnv);
   std::string fName("./hadronFluxStackingASCII.txt");
   if ((clusterNum != 0) || (procNum != 0)) {
//     const char *userEnv = getenv("USER"); // assume it always works
//     std::string aUserStr(userEnv); 
     std::ostringstream fNStrStr; 
     fNStrStr << "./hadronFluxStackingASCII_" 
              << clusterNum << "_" << procNum << ".txt";
     fName = fNStrStr.str();
   }
   fOutHadronAtVectex.open(fName.c_str());
   fOutHadronAtVectex << " Evt X Y Z Px Py Pz Id " << std::endl;
   std::cout << " LBNEStackingAction::OpenHadronAtVertex, Opening file " << 
        fName << std::endl;
   fName = std::string("./hadronFluxStackingASCIICl.txt");
   if ((clusterNum != 0) || (procNum != 0)) {
//     const char *userEnv = getenv("USER"); // assume it always works
//     std::string aUserStr(userEnv); 
     std::ostringstream fNStrStr; 
     fNStrStr << "./hadronFluxStackingASCIICl_" 
              << clusterNum << "_" << procNum << ".txt";
     fName = fNStrStr.str();
   }
   fOutHadronAtVectexCl.open(fName.c_str());
   fOutHadronAtVectexCl << " Evt X Y Z Px Py Pz Id " << std::endl;
   std::cout << " LBNEStackingAction::OpenHadronAtVertex, Opening file " << 
        fName << std::endl;
   
} 

void LBNEStackingAction::PrepareNewEvent ( )

virtual

Definition at line 209 of file LBNEStackingAction.cc.

{
   int verboseLevel = G4EventManager::GetEventManager()->GetTrackingManager()->GetVerboseLevel();
   if(verboseLevel > 1)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno << ": LBNEStackingAction::PrepareNewEvent Called." << std::endl;
   } 

}

void LBNEStackingAction::SetDoImportanceWeightSelection ( bool  t ) const

inline

Definition at line 51 of file LBNEStackingAction.hh.

{ doImportanceWeightSelection = t; } // Pseudo Const.. 

void LBNEStackingAction::SetStackingKillingThreshold ( double  v ) const

inline

Definition at line 53 of file LBNEStackingAction.hh.

{ fStackingKillingThreshold = v;} //Pseudo const 

void LBNEStackingAction::UnKillMuonAlcoveParticles	(	G4ClassificationOfNewTrack &	classification,
		const G4Track *	aTrack
	)

Definition at line 220 of file LBNEStackingAction.cc.

{
  if (classification != fKill) return;

  //
  LBNERunManager* runMan = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
  if (!runMan->GetCreateAlcoveTrackingOutput()) return;
  G4String str = "";
  if (aTrack->GetVolume()->GetLogicalVolume())
    str = aTrack->GetVolume()->GetLogicalVolume()->GetName();
  if (str == "HadronAbsorberSculpAirBuffer" || str =="MuonTrackingBox" || str == "HadronAbsorberSculpBackWall")
  {
    classification = fUrgent;
  }

}

Member Data Documentation

bool LBNEStackingAction::doImportanceWeightSelection

mutableprivate

Definition at line 48 of file LBNEStackingAction.hh.

std::ofstream LBNEStackingAction::fOutHadronAtVectex

mutableprivate

Definition at line 46 of file LBNEStackingAction.hh.

std::ofstream LBNEStackingAction::fOutHadronAtVectexCl

mutableprivate

Definition at line 47 of file LBNEStackingAction.hh.

double LBNEStackingAction::fStackingKillingThreshold

mutableprivate

Definition at line 42 of file LBNEStackingAction.hh.

LBNERunManager* LBNEStackingAction::pRunManager

private

Definition at line 41 of file LBNEStackingAction.hh.

---

The documentation for this class was generated from the following files:

• g4lbne/include/LBNEStackingAction.hh
• g4lbne/src/LBNEStackingAction.cc