LBNEStackingAction.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
// LBNEStackingAction.cc
// $Id: LBNEStackingAction.cc,v 1.2.2.2 2013/09/05 12:32:50 lebrun Exp $
//----------------------------------------------------------------------

#include "LBNEStackingAction.hh"
#include "G4SDManager.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4HCofThisEvent.hh"
#include "G4Track.hh"
#include "G4TrackStatus.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4ios.hh"
#include "LBNEImpWeight.hh"
#include "LBNETrackInformation.hh"
#include "LBNETrajectory.hh"
#include "LBNEAnalysis.hh"
#include "LBNERunManager.hh"
#include "LBNEQuickPiToNu.hh"

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

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

void LBNEStackingAction::OpenHadronAtVertex() const { // 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;
   
} 

//----------------------------------------------------------------------------------
G4ClassificationOfNewTrack LBNEStackingAction::ClassifyNewTrack(const G4Track * aTrack)
{
   //
   // 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;
}
//----------------------------------------------------------------------------------
void LBNEStackingAction::NewStage() 
{
   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::PrepareNewEvent()
{
   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::UnKillMuonAlcoveParticles(G4ClassificationOfNewTrack &classification, const G4Track* aTrack)
{
  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;
  }

}



//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------
void LBNEStackingAction::KillEMParticles(G4ClassificationOfNewTrack &classification, 
                                         const G4Track * aTrack)
{
    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)
{
   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)

{

    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::KillOutOfWorldParticles(G4ClassificationOfNewTrack& classification, 
                                                 const G4Track * aTrack)
{

   //check if track is inside world (some neutral particles make huge jumps from horns (field part) and
   // end up decaying in ~infinity which occasionaly causes g4numi to crash
   
   int verboseLevel = G4TrackingManager::GetVerboseLevel()();
   
  G4ThreeVector position=aTrack->GetPosition();
   if ((classification != fKill) &&
       ((sqrt(position[0]*position[0]+position[1]*position[1])>LBNEData->GetRockRadius())||
        position[2]>LBNEData->GetRockHalfLen()))
   {
      if (verboseLevel > 2)
      {
         G4cout <<"LBNEStackingAction: Killing Out of World Particle" <<G4endl;
         G4cout << "   Particle type: "<<aTrack->GetDefinition()->GetParticleName()
                << " ; Parent ID: "<<aTrack->GetParentID()
                << " ;  Kinetic Energy = "<<aTrack->GetKineticEnergy()/GeV<<" GeV"<<G4endl;
         G4cout << "   Position (m) = "<<position/m<<G4endl;
      }
      classification =fKill;
   }
   
}
*/
//----------------------------------------------------------------------------------
void LBNEStackingAction::KillUnimportantParticles(G4ClassificationOfNewTrack& classification, 
                                                  const G4Track * aTrack)
{
      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);
    }
  } 
}