LBNETrackingAction.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
// LBNETrackingAction.cc
// $Id: LBNETrackingAction.cc,v 1.1.1.1.2.2 2013/10/28 12:06:50 lebrun Exp $
//----------------------------------------------------------------------


#include "G4TrackingManager.hh"
#include "G4Track.hh"
#include "G4Trajectory.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4EventManager.hh"

#include "LBNETrackInformation.hh"
#include "LBNETrackingAction.hh"
#include "LBNERunManager.hh"
#include "LBNEAnalysis.hh"
#include "LBNETrajectory.hh"
#include "LBNEPrimaryGeneratorAction.hh"
#include "LBNEEventAction.hh"
#include "LBNEQuickPiToNu.hh"
#include "LBNESteppingAction.hh"


//--------------------------------------------------------------------------------------
LBNETrackingAction::LBNETrackingAction()
{
  pRunManager=(LBNERunManager*)LBNERunManager::GetRunManager();
  NPGA=(LBNEPrimaryGeneratorAction*)pRunManager->GetUserPrimaryGeneratorAction();
 // No message, fpTrackingManager pointer not set yet ..
  doImportanceWeightSelection = true; // Per G4Numi tradition.. 
}
//--------------------------------------------------------------------------------------
LBNETrackingAction::~LBNETrackingAction()
{
   if(fpTrackingManager->GetVerboseLevel() > 0)
   {
      std::cout << "LBNETrackingAction Destructor Called." << std::endl;
   }
}    
void LBNETrackingAction::OpenHadronAtVertex() const { // Pseudo const.. 
   if (fOutHadronAtVectex.is_open()) {
     std::cerr << " From LBNETrackingAction::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("./hadronFluxPreTrackingASCII.txt");
   if ((clusterNum != 0) || (procNum != 0)) {
//     const char *userEnv = getenv("USER"); // assume it always works
//     std::string aUserStr(userEnv); 
     std::ostringstream fNStrStr; 
     fNStrStr << "./hadronFluxAtPreTrackingASCII_" 
              << clusterNum << "_" << procNum << ".txt";
     fName = fNStrStr.str();
   }
   fOutHadronAtVectex.open(fName.c_str());
   fOutHadronAtVectex << " Evt X Y Z Px Py Pz Id " << std::endl;
   std::cout << " LBNETrackingAction::OpenHadronAtVertex, Opening file " << 
        fName << std::endl;
   
} 

//--------------------------------------------------------------------------------------
void LBNETrackingAction::PreUserTrackingAction(const G4Track* aTrack)
{

   const LBNESteppingAction* LBNEUStep=
    reinterpret_cast<const LBNESteppingAction*>(pRunManager->GetUserSteppingAction());
   LBNEUStep->ResetFlagParticleGotToHAFront();
   if (fOutHadronAtVectex.is_open()) {
      int myPid = aTrack->GetParticleDefinition()->GetPDGEncoding();
      if ((std::abs(myPid) == 211) || (myPid == 2212) || 
          (std::abs(myPid) == 321) || (std::abs(myPid) == 311)) {     
        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 > 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 << " LBNETrackingAction::PreUserTrackingAction, track ID " << aTrack->GetTrackID() 
//             << " Track Ptr " << (void *) aTrack   << " User info " << oldTrInfo << std::endl;
//   if ((oldTrInfo != 0) && (oldTrInfo->GetPFFlag() == 1)) {
//     std::cerr << " ............. skip, job done in LBNEPerfectFocusProcess " << std::endl;
//     return;
//   }
   // 
   //
   // This was a track whose Pt has been set to zero, and redefine. No further information is needed. 
   //

//   G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
   if(fpTrackingManager->GetVerboseLevel() > 1)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno << ": LBNETrackingAction::PreUserTrackingAction() Called." << std::endl;
      if(fpTrackingManager->GetVerboseLevel() > 2)
      {
         LBNETrackInformation* trackinfo=new LBNETrackInformation(); 
         trackinfo -> Print(aTrack);
         delete trackinfo;
      }
   }
   // Debugging the infinite memory blowup.. 
//   if (evtno == 32135) {
//     std::cerr << " At track  " << aTrack->GetTrackID() << " At " << aTrack->GetPosition() 
//               << " P " << aTrack->GetMomentumDirection() << " E " << aTrack->GetTotalEnergy() 
//             << " name " << aTrack->GetParticleDefinition()->GetParticleName() <<  std::endl;
//   }
   
    const LBNESteppingAction *pStep = dynamic_cast<const LBNESteppingAction*>(pRunManager->GetUserSteppingAction());
    pStep->ResetNumSteps();
   
      LBNETrackingAction::SetPreLBNETrackInformation(aTrack);
      
    if(fpTrackingManager->GetVerboseLevel() > 2)  std::cerr << " .... ..... about to store trajectory " << std::endl;
     
      //Store tracks in trajectory container
      fpTrackingManager->SetStoreTrajectory(true);
      fpTrackingManager->SetTrajectory(new LBNETrajectory(aTrack));
      
      if(fpTrackingManager->GetVerboseLevel() > 2)  std::cerr << " .... ..... did store trajectory " << std::endl;
      LBNETrackingAction::AnalyzeIfNeutrino(aTrack);
      if(fpTrackingManager->GetVerboseLevel() > 2)  std::cerr << " .... ..... did Analyze if neutrino " << std::endl;
   
   
   
  if(fpTrackingManager->GetVerboseLevel() > 1)
  {
     std::cout << "         LBNETrackingAction::PreUserTrackingAction() - done." << std::endl;
  }

}
//--------------------------------------------------------------------------------------
void LBNETrackingAction::PostUserTrackingAction(const G4Track* aTrack)
{
   if(fpTrackingManager->GetVerboseLevel() > 1)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno << ": LBNETrackingAction::PostUserTrackingAction() Called." << std::endl;
      if(fpTrackingManager->GetVerboseLevel() > 2)
      {
         LBNETrackInformation* trackinfo=new LBNETrackInformation(); 
         trackinfo -> Print(aTrack);
         delete trackinfo;
      }
   }
// 
// specific checks for impotance weight. 
//
//   if (std::abs(aTrack->GetParticleDefinition()->GetPDGEncoding()) == 211) { 
//    G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
//    LBNETrackInformation *trInfo = reinterpret_cast<LBNETrackInformation*> (aTrack->GetUserInformation());  
//    if (trInfo != 0) std::cerr << "LBNETrackingAction::PostUserTrackingAction, Event " << evtno 
//               << " Track ID " << aTrack->GetTrackID() << " Imp Weight " << trInfo->GetNImpWt() << std::endl;
//    else  std::cerr << "LBNETrackingAction::PostUserTrackingAction, Event " << evtno 
//               << " Track ID " << aTrack->GetTrackID() << " has no weight " << std::endl;       
//    
//   }
   LBNETrackingAction::SetPostLBNETrackInformation(aTrack);


    if(fpTrackingManager->GetVerboseLevel() > 1)
    {
       std::cout << "         LBNETrackingAction::PostUserTrackingAction() - done." << std::endl;
    }


}





//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
void LBNETrackingAction::SetPreLBNETrackInformation(const G4Track* aTrack)
{

   if(fpTrackingManager->GetVerboseLevel() > 3)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno 
                << ": NumiTrackingAction::SetPreNumiTrackInformation() Called." 
                << std::endl;
   }

   //set tgen (and weight for fluka nad mars input) 
   if (aTrack->GetTrackID()==1) 
   {   
      LBNETrackInformation* info = new LBNETrackInformation();
      LBNERunManager* theRunManager=dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
      if (theRunManager->GetUseFlukaInput() || theRunManager->GetUseMarsInput())
      {
         info->SetTgen(NPGA->GetTgen());
         info->SetNImpWt(NPGA->GetWeight());
         if (doImportanceWeightSelection) info->SetNImpWt(NPGA->GetWeight()); // set weight of the new track equal to parent weight
         else info->SetNImpWt(1.0);
//         std::cerr << " LBNETrackingAction::SetPreLBNETrackInformation, setting imp weight " << NPGA->GetWeight() << std::endl;
      }
      else
      {
         info->SetTgen(1);
         info->SetNImpWt(1.);
//         std::cerr << " LBNETrackingAction::SetPreLBNETrackInformation, setting imp weight to 1 " << std::endl;
      }
      G4Track* theTrack = (G4Track*)aTrack;
      theTrack->SetUserInformation(info);
   }
   

}

//--------------------------------------------------------------------------------------
void LBNETrackingAction::SetPostLBNETrackInformation(const G4Track* aTrack)
{
   if(fpTrackingManager->GetVerboseLevel() > 3)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno 
                << ": NumiTrackingAction::SetPostNumiTrackInformation() Called." 
                << std::endl;
   }
   
   // set tgen(secondary) = tgen(parent)+1
   LBNETrackInformation* info = (LBNETrackInformation*)(aTrack->GetUserInformation());
   if (info!=0) 
   {
      G4int tgen               = info->GetTgen();
      G4double nimpwt          = info->GetNImpWt();
      G4TrackVector* SecVector = fpTrackingManager->GimmeSecondaries();
      for (size_t ii=0;ii<(*SecVector).size();ii++)
      {
         LBNETrackInformation* newinfo = new LBNETrackInformation();
         newinfo->SetTgen(tgen+1); // set generation of daughter particles
         if (doImportanceWeightSelection) newinfo->SetNImpWt(nimpwt); // set weight of the new track equal to parent weight
         else newinfo->SetNImpWt(1.0);
         newinfo->SetPFFlag(info->GetPFFlag()); // set weight of the new track equal to parent weight
        (*SecVector)[ii]->SetUserInformation(newinfo);
        // The above algorithm makes sense, but, in most cases, is a waste of time, as it assign a weight 
        // to electrons, ions, etc,, that were collected while the track moved though the beam line
        // Performance could be improved.. 
        // 
//         if (doImportanceWeightSelection &&  
//           ((*SecVector)[ii]->GetDynamicParticle()->GetParticleDefinition()->GetPDGEncoding() == 211)) {
//         std::cerr << " LBNETrackingAction::SetPostLBNETrackInformation, setting imp weight " 
//                   << nimpwt << " E " << (*SecVector)[ii]->GetTotalEnergy() 
//                   << " name " << (*SecVector)[ii]->GetDynamicParticle()->GetParticleDefinition()->GetParticleName() <<  std::endl;
//         }
      }
   }
}

//--------------------------------------------------------------------------------------
void LBNETrackingAction::AnalyzeIfNeutrino(const G4Track* aTrack)
{
   if(fpTrackingManager->GetVerboseLevel() > 3)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno 
                << ": NumiTrackingAction::AnalyzeIfNeutrino() Called." 
                << std::endl;
   }

  //if a particle is a neutrino then analyse and store in ntuple
  G4ParticleDefinition * particleDefinition = aTrack->GetDefinition();
  if ((particleDefinition == G4NeutrinoE::NeutrinoEDefinition())||
      (particleDefinition == G4NeutrinoMu::NeutrinoMuDefinition()) ||
      (particleDefinition == G4NeutrinoTau::NeutrinoTauDefinition()) ||
      (particleDefinition == G4AntiNeutrinoE::AntiNeutrinoEDefinition()) ||
      (particleDefinition == G4AntiNeutrinoMu::AntiNeutrinoMuDefinition()) ||
      (particleDefinition == G4AntiNeutrinoTau::AntiNeutrinoTauDefinition()))
    {
      
      //Following part might be temporary....implemented to better understand..
      //..the particles through volume. A. Bashyal
                const G4Event* event = pRunManager->GetCurrentEvent();
        G4TrajectoryContainer* trajectories = event->GetTrajectoryContainer();
            // Make a temporary map to make it easier to go up the trajectory stack
        std::map<int,G4VTrajectory*> trajectoryMap;
        for (std::size_t i = 0; i < trajectories->size(); ++i) {
            G4VTrajectory* traj = (*trajectories)[i];
            trajectoryMap[traj->GetTrackID()] = traj;
        }

        std::vector<G4VTrajectory*> nuHistory;
        G4VTrajectory* neutrino = fpTrackingManager->GimmeTrajectory();
        nuHistory.push_back(neutrino);
        int parentId = aTrack->GetParentID();
        while (parentId  != 0) {
            G4VTrajectory* parentTraj = trajectoryMap[parentId];
            if (!parentTraj) {
                G4cerr << "Invalid trajectory object" << G4endl;
                continue;
            }
            nuHistory.push_back(parentTraj);
            LBNETrajectory* lbneTrajectory = dynamic_cast<LBNETrajectory*>(parentTraj);
            if (!lbneTrajectory) continue;
            G4String startVol = lbneTrajectory->GetPreStepVolumeName(0);
            G4String stopVol = lbneTrajectory->GetPreStepVolumeName(lbneTrajectory->GetPointEntries()-1);
            parentId = parentTraj->GetParentID();
        }
      LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      analysis->FillNeutrinoNtuple(*aTrack,nuHistory);
     } 
}