LBNESteppingAction.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
// LBNESteppingAction.cc
// $Id: LBNESteppingAction.cc,v 1.1.1.1.2.16 2013/12/02 20:24:52 lebrun Exp $
//----------------------------------------------------------------------

//C++
#include <string>

#include "LBNESteppingAction.hh"

#include "G4ProcessManager.hh"
#include "G4SteppingManager.hh"
#include "G4Track.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4TrajectoryContainer.hh"
#include "G4RunManager.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Event.hh"
#include "G4VTouchable.hh"
#include "G4TouchableHistory.hh"
#include "G4EventManager.hh"
#include "G4FieldManager.hh"
#include "G4TrackingManager.hh"

#include "LBNETrackInformation.hh"
#include "LBNETrajectory.hh"
#include "LBNETrackInformation.hh"
#include "LBNEAnalysis.hh"
#include "LBNEEventAction.hh"
#include "LBNERunManager.hh"
#include "LBNEMagneticField.hh"
#include "LBNFMagneticFieldPolyconeHorn.hh"
#include "LBNEVolumePlacements.hh"
#include "LBNEQuickPiToNu.hh"
#include "LBNEDetectorConstruction.hh"
#include "LBNFDeDxMap.hh"
#include "TFile.h"
#include "TTree.h"

//----------------------------------------------------------------------
LBNESteppingAction::LBNESteppingAction()
{
  pRunManager=(LBNERunManager*)LBNERunManager::GetRunManager();
 fKillTrackingThreshold = 0.050*CLHEP::GeV;
 if (pRunManager->GetVerboseLevel() > 0) {
   std::cerr << " LBNESteppingAction::LBNESteppingAction called ... " << std::endl;
 }
 pMessenger = new LBNESteppingActionMessenger(this); 
 fEvtIdPrevious = -1;
 doGenerateHadronFluxFromTargetASCII = false;
 fStudyGeantinoMode=G4String("none");
 fKeyVolumeForOutput=G4String("blank");
 fKeyVolumeForOutputTo=G4String("blank");
 fNConsecutivSmallSteps=0;
 fNumTracksKilledAsStuck=0;
 fNumStepsCurrentTrack=0;
 doStudyParticleThroughHorns=false;
 GenerateMuonSculptedAbsorberFlux = false;
 fEnergyDepInGraphite = 0.;
 fEnergyDepInArgonGasHorn1 = 0;
 fEnergyDepInArgonGasHorn2 = 0;
 fGenerateVoxelsZScale = 25.0;
 fGenerateVoxelsZOffset = 5000.0;
 fGenerateVoxelsZMax = 25000.0;
 fGenerateVoxelsRScale = 5.0;
 fGenerateVoxelsRMax = 2000.0;
 fGenerateVoxelsIRKMax = fGenerateVoxelsRScale*std::log(fGenerateVoxelsRMax);
 fOutMuonSculptedAbsorberFilename="MuonFluxAtSculptedAbsorber";
 // 
 // April 2017 study for MARS, cmpare meson rate.. 
 //
 bool doSergeiStriganovApr2017 = false;
 if (doSergeiStriganovApr2017) { 
   std::string dirNameTmp("/dune/data/users/lebrun/LBNFMARSApril2017/");
   std::string fNameTgt(dirNameTmp); fNameTgt += std::string("DumpOptimizedEngineered_Tgt.txt");
   std::ofstream *fOutTargetSergei = new std::ofstream(fNameTgt.c_str());
   fOutPtrsForMarsCmpApr2017.push_back(fOutTargetSergei);
   //
   std::string fNameHornA(dirNameTmp); fNameHornA += std::string("DumpOptimizedEngineered_HornA.txt");
   std::ofstream *fOutHornASergei = new std::ofstream(fNameHornA.c_str());
   fOutPtrsForMarsCmpApr2017.push_back(fOutHornASergei);
   //
   std::string fNameHornB(dirNameTmp); fNameHornB += std::string("DumpOptimizedEngineered_HornB.txt");
   std::ofstream *fOutHornBSergei = new std::ofstream(fNameHornB.c_str());
   fOutPtrsForMarsCmpApr2017.push_back(fOutHornBSergei);
   //
   std::string fNameHornC(dirNameTmp); fNameHornC += std::string("DumpOptimizedEngineered_HornC.txt");
   std::ofstream *fOutHornCSergei = new std::ofstream(fNameHornC.c_str());
   fOutPtrsForMarsCmpApr2017.push_back(fOutHornCSergei);
   for (size_t k=0; k!= fOutPtrsForMarsCmpApr2017.size(); k++) {
     (*fOutPtrsForMarsCmpApr2017[k]) << " evt trId pdg x y z px py pz ek " << std::endl;
   }
 }
}
//----------------------------------------------------------------------
LBNESteppingAction::~LBNESteppingAction()
{
 if (fOutStudy.is_open()) {
    if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) {
     std::cerr << " Dumping in for RZVoxels, num Voxels = " << fRZVoxelsData.size() << std::endl;
     for(std::map<int, int>::const_iterator itMRZ=fRZVoxelsData.begin(); 
              itMRZ != fRZVoxelsData.end(); itMRZ++) {
       int ik = itMRZ->first;
       int iz = ik/fGenerateVoxelsIRKMax;
       int ir = ik - iz*fGenerateVoxelsIRKMax;
       const double z = fGenerateVoxelsZScale*iz - fGenerateVoxelsZOffset; //see in GenerateRZVoxel..
       const double r = std::exp(static_cast<double>(ir)/fGenerateVoxelsRScale);
       fOutStudy << " " << ik << " " << r << " " << z << " " << itMRZ->second << std::endl;
     }
   }
   fOutStudy.close();
 }
 if (fOutStepStudy.is_open()) fOutStepStudy.close();
 if(makeSteppingTuple)  {
   steppingTupleFile->cd();
   steppingTuple->Write();
   steppingTupleFile->Close();
 } 
 
 delete pMessenger;
 // delete 
 //
   for (size_t k=0; k!= fOutPtrsForMarsCmpApr2017.size(); k++) {
     fOutPtrsForMarsCmpApr2017[k]->close();
     delete fOutPtrsForMarsCmpApr2017[k];
   }
   
}

//----------------------------------------------------------------------
void LBNESteppingAction::UserSteppingAction(const G4Step * theStep)
{
 // 
 // Temporary timing study: kill the track outside the target region,
 //
// { 
//   G4StepPoint *prePtr = theStep->GetPostStepPoint();
//   G4ThreeVector posT = prePtr->GetPosition();
//   const double rrSq = posT[0]*posT[0] + posT[1]*posT[1];
//   G4Track *theTrack = theStep->GetTrack();
//   if (rrSq > 310.) {
//      if (posT[2] < 2500.) theTrack->SetTrackStatus(fStopAndKill);
//   } else {
//     if (posT[2] > 2500.) theTrack->SetTrackStatus(fStopAndKill);
//   }
// }
 // Killing tracks that are stuck on a boundary crossing 
 // This should not happen if the geometry has no overlap and is strictly 
 // correct.  For the nominal case, except for the hadron absorver, 
 // no overlaps have been detected at construction stage. 
 // However, the tests are "weak", i.e., do not explore all corners. 
 // Let us protect our-self here again infinite loop 
 //
 if (theStep->GetStepLength() < 1.0e-15) {   
   fNConsecutivSmallSteps++;
   if (fNConsecutivSmallSteps > 5000) {
       G4StepPoint *prePtr = theStep->GetPreStepPoint();
       G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();       
        std::cerr << " LBNESteppingAction, too man consecutive small steps  " << std::endl
          << ".. from " << volPre->GetName() 
                     << " detected at " <<  prePtr->GetPosition() <<
                      " last step length " << theStep->GetStepLength() << std::endl;
        G4Track* aTrack= theStep->GetTrack();                 
        std::cerr << " ...  for track  " << aTrack->GetTrackID() << " At " << aTrack->GetPosition() 
               << " P " << aTrack->GetMomentumDirection() << " E " << aTrack->GetTotalEnergy() 
               << " name " << aTrack->GetParticleDefinition()->GetParticleName() <<  std::endl;
        aTrack->SetTrackStatus(fStopAndKill);      
        std::cerr << " Track killed .... " << std::endl;
        fNumTracksKilledAsStuck++;
   }
   
 } else fNConsecutivSmallSteps=0;
 // Last check on memory explosion: No more than 500 mgB usage. Check than our step counter (can't trust G4 anylonger) 
 // gets too large.. 
 // G4 10 seems be a bit more demanding.  Relax these limits. 
  fNumStepsCurrentTrack++;
  if (fNumStepsCurrentTrack > 200000) {
      std::cerr << " LBNESteppingAction... Anomalously large number of steps.. kill this event.. " << std::endl;
      G4Exception("LBNESteppingAction", " ", EventMustBeAborted, " Too many steps.. ");
      return; 
//    const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
//    runAct->CheckMemoryUsage(1000000);
  }
  //
  // June 2017: we were asked to provide dedx in the horns. 
  //
  const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
  if (runAct->GetDoComputeEDepInHorns()) { 
     LBNFDeDxMap *aMapDedX = runAct->GetG4MARSDeDxMap();
     if (aMapDedX != 0) aMapDedX->fill(theStep);
     G4StepPoint *postPtr = theStep->GetPreStepPoint();
     // We don't track past the horns, not interested in neutrino spectrum in this slow tedious mode.. 
     // Assume extended chase. In fact, I am asked to do it only for horn A 
     if (postPtr->GetPosition()[2] > 3000.)  theStep->GetTrack()->SetTrackStatus(fStopAndKill);
  }   
  //
  // Code bloat, April 2017.. 
  //
    if (fOutPtrsForMarsCmpApr2017.size() > 0) this->CheckHadronsMarsCmpApr2017(theStep);

   // Debugging the bad neutrons in 4.9.6.p02 
//   G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();

//   if (evtno == 32135) {
//     const G4Track* aTrack= theStep->GetTrack(); 
//     if (aTrack->GetTrackID() == 79) {
//        G4StepPoint *prePtr = theStep->GetPreStepPoint();
//        G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();    
//        std::cout << " LBNESteppingAction::StudyPropagation, bad track 79 " << std::endl
//        << ".. from " << volPre->GetName() 
//                   << " detected at " <<  prePtr->GetPosition() <<
//                    " length " << theStep->GetStepLength() << std::endl;
// Found :                    
//      LBNESteppingAction::StudyPropagation, bad track 79 
//.. from Horn1IOTransCont detected at (31.0918,-3.71985,109.441) length 1.77636e-15
// LBNESteppingAction::StudyPropagation, bad track 79 
//.. from Horn1IOTransInnerPart3 detected at (31.0918,-3.71985,109.441) length 0
// LBNESteppingAction::StudyPropagation, bad track 79 
//.. from Horn1IOTransCont detected at (31.0918,-3.71985,109.441) length 1.77636e-15
// For ever... 
//
//     }
//   }

   int verboseLevel = 
   G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager()->GetverboseLevel();
   if(verboseLevel > 3)
   {
      G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
      std::cout << "Event " << evtno << ": LBNESteppingAction::UserSteppingAction() Called." << std::endl;
   }
   if((fStudyGeantinoMode.find("Mag") != std::string::npos) 
       && (pRunManager->GetCurrentEvent()->GetEventID() < 5)) 
        this->dumpStepCheckVolumeAndFields(theStep); 

   LBNESteppingAction::KillNonNuThresholdParticles(theStep);
   
   if (doGenerateHadronFluxFromTargetASCII) this->FillHadronFluxFromTargetASCII(theStep);
   
   LBNESteppingAction::CheckDecay(theStep);
   LBNESteppingAction::SetMomentumInfoForParticle(theStep);
   if (LBNEQuickPiToNuVect::Instance()->doIt()) StudyPionsThroughHorn2(theStep);
   if (fOutStudy.is_open()) {
    if(doStudyParticleThroughHorns) {
//       std::cerr << " Stepping .... doStudyParticleThroughHorns.. And quit  " << std::endl; exit(2);
      StudyParticleThroughHorns(theStep);
    }
    else {

      if (fStudyGeantinoMode.find("Absorb") != std::string::npos) StudyAbsorption(theStep);
      if (fStudyGeantinoMode.find("Propa") != std::string::npos) StudyPropagation(theStep);
      if (fStudyGeantinoMode.find("PropCO") != std::string::npos) StudyCheckOverlap(theStep);
      if (fStudyGeantinoMode.find("MagTilt") != std::string::npos) StudyCheckMagneticTilts(theStep);
      if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) GenerateVolumeCrossingsRZMap(theStep);
    }
   }
   if (makeSteppingTuple) 
     StudyAbsorption(theStep);
   if (GenerateMuonSculptedAbsorberFlux) this->StudyMuonSculptedAbsorberFlux(theStep);

   // Add energy loss if requested.. 
   const LBNERunAction *runUserAction =  reinterpret_cast<const LBNERunAction *>(G4RunManager::GetRunManager()->GetUserRunAction());
   if (runUserAction->GetDoComputeEDepInGraphite()) {
     const std::string tgt("Target");
     const std::string tg = theStep->GetTrack()->GetMaterial()->GetName();
     if (tg == tgt) fEnergyDepInGraphite += theStep->GetTotalEnergyDeposit(); // Not true for all targets!!!! Must have the right material name. 
   }
   if (runUserAction->GetDoComputeEDepInArgonGas()) {
     const G4StepPoint *prePt = theStep->GetPreStepPoint();
     if ((prePt != 0) && (prePt->GetPhysicalVolume() != 0)) {
         const G4LogicalVolume *aVol = prePt->GetPhysicalVolume()->GetLogicalVolume();
         std::string theMatName(aVol->GetMaterial()->GetName());
         if (theMatName == std::string("Argon")) {
           std::string theVolName(aVol->GetName());
           if (theVolName.find("orn1") != std::string::npos) {
              fEnergyDepInArgonGasHorn1 += theStep->GetTotalEnergyDeposit();
//            std::cerr << " fEnergyDepInArgonGasHorn1 now " << fEnergyDepInArgonGasHorn1 << std::endl;
           } else if (theVolName.find("orn2") != std::string::npos) {
              fEnergyDepInArgonGasHorn2 += theStep->GetTotalEnergyDeposit();
//            std::cerr << " fEnergyDepInArgonGasHorn2 now " << fEnergyDepInArgonGasHorn2 << std::endl;
           } else {
             std::string msg(" Unknow volume that contains Argon, volName "); msg+=theVolName;
             G4Exception("LBNESteppingaction", "", FatalException, msg.c_str());
           }
         }           
      }
   }

  // 
  // Adding checks for Perfect Focsuing..
  //
//  const G4Track *theTrack = theStep->GetTrack();
//  const G4DynamicParticle *pDef = theTrack->GetDynamicParticle();
//  const int pId = pDef->GetDefinition()->GetPDGEncoding();
//  if (((std::abs(pId) == 211) || (std::abs(pId) == 321)) && (theStep->GetPreStepPoint()->GetPosition()[2] < 2000.*CLHEP::mm)) {
//      const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(theTrack->GetUserInformation());
//      std::cerr << " Pion or Kaon, Id " << theTrack->GetTrackID() << "  at " << theStep->GetPreStepPoint()->GetPosition() 
//                << " Momentum " << theTrack->GetMomentum() << " Flag " << oldTrInfo->GetPFFlag() << std::endl;
//  }

  if (pRunManager->GetCreateTrkPlaneOutput()){
    //---tracking planes
    TrkPlnLogical =
         ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
          ->GetUserDetectorConstruction())->GetMuonDetectorLogical();

   G4StepPoint *postStep = theStep->GetPostStepPoint();
   if(postStep->GetPhysicalVolume()){
        G4LogicalVolume *postStepVolume =
          postStep->GetPhysicalVolume()->GetLogicalVolume();
        if(postStepVolume){
          //  std::cout << postStepVolume->GetName() << std::endl;        
          if(postStepVolume == TrkPlnLogical){
            CheckInTrackingDetectorPlane(theStep);
            }
        }
   }
  }

if(pRunManager->GetCreateTargetOutput()){
          
            CheckInTargetOutput(theStep);
            }
 /*       
if(pRunManager->GetCreateTrkPlaneDPOutput()){

           CheckInTrackingDetectorDPPlane(theStep);
}   
*/

// To stepping method to save the particle info that makes through the tracking plane at the end of Horn1 and horn 2. Amit Bashyal
 if (pRunManager->GetCreateTrkPlaneH1Output()){
    //---tracking planes
    TrkPlnH1Logical =
         ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
          ->GetUserDetectorConstruction())->GetHorn1PlaneLogical();

   G4StepPoint *postStep = theStep->GetPostStepPoint();
   if(postStep->GetPhysicalVolume()){
        G4LogicalVolume *postStepVolume =
          postStep->GetPhysicalVolume()->GetLogicalVolume();
        if(postStepVolume){
         // std::cout << postStepVolume->GetName() << std::endl;          
          if(postStepVolume == TrkPlnH1Logical){
          
            CheckInTrackingDetectorH1Plane(theStep);
            }
        }
   }
}

if (pRunManager->GetCreateTrkPlaneH2Output()){
    //---tracking planes
    TrkPlnH2Logical =
         ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
          ->GetUserDetectorConstruction())->GetHorn2PlaneLogical();

   G4StepPoint *postStep = theStep->GetPostStepPoint();
   if(postStep->GetPhysicalVolume()){
        G4LogicalVolume *postStepVolume =
          postStep->GetPhysicalVolume()->GetLogicalVolume();
        if(postStepVolume){
          //  std::cout << postStepVolume->GetName() << std::endl;        
          if(postStepVolume == TrkPlnH2Logical){
            CheckInTrackingDetectorH2Plane(theStep);
            }
        }
   }
}
//Although the tracking plane at the end (or start) of Decay Pipe has been created and the method to save the particle info is set, the Hadron monitor is likely to
// be in this position.
if (pRunManager->GetCreateTrkPlaneDPOutput()){
    //---tracking planes
    TrkPlnDPLogical =
         ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
          ->GetUserDetectorConstruction())->GetDecayPipePlaneLogical();

   G4StepPoint *postStep = theStep->GetPostStepPoint();
   if((postStep != 0) && (postStep->GetPhysicalVolume() != 0)){
        G4LogicalVolume *postStepVolume =
          postStep->GetPhysicalVolume()->GetLogicalVolume();
        if(postStepVolume){
          //  std::cout << postStepVolume->GetName() << std::endl;        
          if(postStepVolume == TrkPlnDPLogical){
            CheckInTrackingDetectorDPPlane(theStep);
            }
        }
   }
}

//Sculpted absorber alcove tracking plane
  if (pRunManager->GetCreateAlcoveTrackingOutput()){
    CheckInAlcoveTrackingPlane(theStep);
  }

}


//More PPFX stuff...yea!!! A. Bashyal
//----------------------------------------------------------------------------------------------------------------
void LBNESteppingAction::SetMomentumInfoForParticle(const G4Step * theStep)
{
G4Track * aTrack = theStep->GetTrack();
if(theStep->GetPreStepPoint()->GetProcessDefinedStep()==NULL)return;
    G4ThreeVector mom = G4ThreeVector(0.,0.,0.);
 //Basically imported ditto from NumiSteppingAction for the sake of dk2nu ntuples..A Bashyal
        G4int nSecAtRest = fpSteppingManager->GetfN2ndariesAtRestDoIt();
        G4int nSecAlong  = fpSteppingManager->GetfN2ndariesAlongStepDoIt();
        G4int nSecPost   = fpSteppingManager->GetfN2ndariesPostStepDoIt();
        G4int nSecTotal  = nSecAtRest+nSecAlong+nSecPost;
        G4TrackVector* secVec = fpSteppingManager->GetfSecondary();
        //std::cout<<(*secVec).size()<<" "<<nSecTotal<<std::endl;
        //G4bool isShower = false;
         LBNETrackInformation* trackInformation = (LBNETrackInformation*)(aTrack->GetUserInformation());
         if(nSecTotal==0)return;
  //if(nSecTotal>0) {
     
      for(size_t lp1=(*secVec).size()-nSecTotal; lp1<(*secVec).size(); ++lp1) {
          G4Track* secTrack = (*secVec)[lp1];
        //  std::cout<<"trackid for secondary tracks "<<(*secVec)[lp1]->GetTrackID()<<std::endl;
          G4ParticleDefinition* def = secTrack->GetDefinition();
          //std::cout<<"LBNESteppingAction... "<<def<<std::endl;
              // Showering particles (e+/-,gamma) are not interesting for neutrino production,
              // skip them
          if (def == G4Electron::Electron() || def == G4Positron::Positron() ||
              def == G4Gamma::Gamma())continue;
             mom = theStep->GetPreStepPoint()->GetMomentum();
              
          }    
         // if(!isShower) mom = theStep->GetPostStepPoint()->GetMomentum();
          
         // std::cout<<"LBNESteppingAction:: "<<mom<<std::endl;
          
        //  LBNETrackInformation* trackInformation = (LBNETrackInformation*)(aTrack->GetUserInformation());
         // std::cout<<"I am here SteppingAction"<<std::endl;
          if (trackInformation != 0) {
         // std::cout<<" Now inside the if statement LBNESteppingaction"<<std::endl;
         //std::cout<<"LBNESteppingAction:: "<<mom<<std::endl;
             // LBNETrackInformation* trackInformation = new LBNETrackInformation;
           //  assert(trackInformation != NULL);
              trackInformation->SetParentMomentumAtThisProduction(mom);
            //if(mom != G4ThreeVector(0.,0.,0.))  std::cout<<"checking the function "<<mom<<std::endl;
              aTrack->SetUserInformation(trackInformation);
         } else {
        // std::cout<<"And the else statement.."<<std::endl;
              LBNETrackInformation* newtrackInformation = new LBNETrackInformation();
                  //= dynamic_cast<LBNETrackInformation*>(trackInformation->GetUserInformation());
              newtrackInformation->SetParentMomentumAtThisProduction(mom);
              aTrack->SetUserInformation(newtrackInformation);
              
          }
        
        //}
    // }
  
}

//----------------------------------------------------------------------
//----------------------------------------------------------------------
void LBNESteppingAction::KillNonNuThresholdParticles(const G4Step * theStep)
{

   const LBNERunAction *runUserAction =  reinterpret_cast<const LBNERunAction*>(G4RunManager::GetRunManager()->GetUserRunAction());
   if (runUserAction->GetDoComputeEDepInHorns()) return; 
   if (runUserAction->GetDoComputeEDepInGraphite()) { 
     const std::string tgt("Target");
     const std::string tg = theStep->GetTrack()->GetMaterial()->GetName();
     if (tgt == tg) return; // we let them go until they escape the target.. 
   }
   if (runUserAction->GetDoComputeEDepInArgonGas()) {
      const G4StepPoint *prePt = theStep->GetPreStepPoint();
      if (prePt->GetPhysicalVolume() !=0 && (prePt->GetPosition()[2] > 9900.)) {
         std::string  vName = theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
         G4Track * theTrack = theStep->GetTrack();
         if (vName.find("Tunnel") != std::string::npos)  theTrack->SetTrackStatus(fStopAndKill); // We don't produce neutrino, 
                     // just compute energy loss in Argon, to compute ionization rate...
         return;              
      }
   }

   if (pRunManager->GetCreateAlcoveTrackingOutput()){
     if(theStep->GetPreStepPoint()->GetPhysicalVolume() != NULL &&
        theStep->GetPostStepPoint()->GetPhysicalVolume()!=NULL) {

       std::string preStepPointName  =
         theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
       std::string postStepPointName  =
         theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
       if (preStepPointName=="HadronAbsorberSculpBackWall" || postStepPointName=="HadronAbsorberSculpBackWall" ||
           preStepPointName=="HadronAbsorberSculpAirBuffer" || postStepPointName=="HadronAbsorberSculpAirBuffer" ||
           preStepPointName=="MuonTrackingVol" || postStepPointName=="MuonTrackingVol" )
         return;//Don't kill tracks in back wall, muon alcove or tracking plane

     }//if step volumes exist
   }//if saving alcove output

   G4Track * theTrack = theStep->GetTrack();
   G4ParticleDefinition * particleDefinition = theTrack->GetDefinition();
   //
   //kill the track if it is below the kill tracking threshold and it is NOT a neutrino
   //
   if ( (theTrack->GetKineticEnergy() < fKillTrackingThreshold) &&
       (particleDefinition != G4NeutrinoE::NeutrinoEDefinition())&&
       (particleDefinition != G4NeutrinoMu::NeutrinoMuDefinition()) &&
       (particleDefinition != G4NeutrinoTau::NeutrinoTauDefinition()) &&
       (particleDefinition != G4AntiNeutrinoE::AntiNeutrinoEDefinition()) &&
       (particleDefinition != G4AntiNeutrinoMu::AntiNeutrinoMuDefinition()) &&
       (particleDefinition != G4AntiNeutrinoTau::AntiNeutrinoTauDefinition()))
   {
      theTrack->SetTrackStatus(fStopAndKill);
   }
   
}


//----------------------------------------------------------------------
void LBNESteppingAction::CheckDecay(const G4Step * theStep)
{
   G4Track * theTrack = theStep->GetTrack();
   G4ParticleDefinition * particleDefinition = theTrack->GetDefinition();
  

  // Check if the Pi+, Pi-, K+, K-, K0L, mu+ or mu- decayed and set Ndecay code:
  // 1  K0L -> nu_e pi- e+
  // 2  K0L -> anti_nu_e pi+ e-
  // 3  K0L -> nu_mu pi- mu+
  // 4  K0L -> anti_nu_mu pi+ mu-
  // 5  K+  -> nu_mu mu+
  // 6  K+  -> nu_e pi0 e+
  // 7  K+  -> nu_mu pi0 mu+
  // 8  K-  -> anti_nu_mu mu-
  // 9  K-  -> anti_nu_e pi0 e-
  // 10 K-  -> anti_nu_mu pi0 mu-
  // 11 mu+ -> anti_nu_mu nu_e e+
  // 12 mu- -> nu_mu anti_nu_e e-
  // 13 pi+ -> nu_mu mu+
  // 14 pi- -> anti_nu_mu mu-

  if (theStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL)
  {
     G4int decay_code=0;
     if (theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()=="Decay")
     {
        G4int nSecAtRest = fpSteppingManager->GetfN2ndariesAtRestDoIt();
        G4int nSecAlong  = fpSteppingManager->GetfN2ndariesAlongStepDoIt();
        G4int nSecPost   = fpSteppingManager->GetfN2ndariesPostStepDoIt();
        G4int nSecTotal  = nSecAtRest+nSecAlong+nSecPost;
        G4TrackVector* secVec = fpSteppingManager->GetfSecondary();
        
        if (particleDefinition==G4PionPlus::PionPlusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
                 decay_code=13;
           }
        }
        if (particleDefinition==G4PionMinus::PionMinusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
                 decay_code=14;
           }
        }
        if (particleDefinition==G4KaonPlus::KaonPlusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
              {if (nSecTotal==2) decay_code=5;
                if (nSecTotal==3) decay_code=7;}
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_e")
                 decay_code=6;
           }
        }
        if (particleDefinition==G4KaonMinus::KaonMinusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
              {if (nSecTotal==2) decay_code=8;
                 if (nSecTotal==3) decay_code=10;}
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_e")
                 decay_code=9;
           }
        }
        if (particleDefinition==G4KaonZeroLong::KaonZeroLongDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
            {
               if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_e")
                  decay_code=1;
               if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_e")
                  decay_code=2;
               if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
                  decay_code=3;
               if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
                  decay_code=4;     
            }
        }
        if (particleDefinition==G4MuonPlus::MuonPlusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
                 decay_code=11; 
           }
        }
        if (particleDefinition==G4MuonMinus::MuonMinusDefinition())
        {
           for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
           {
              if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
                 decay_code=12; 
           }
        }
        
        LBNETrackInformation* oldinfo=(LBNETrackInformation*)(theTrack->GetUserInformation()); 
        if (oldinfo!=0)
        {
           oldinfo->SetDecayCode(decay_code);                                                      
          theTrack->SetUserInformation(oldinfo); 
        }
        else
        {
           LBNETrackInformation* newinfo=new LBNETrackInformation(); 
          newinfo->SetDecayCode(decay_code);                                                       
          theTrack->SetUserInformation(newinfo); 
        }
     }
  }
}

//----------------------------------------------------------------------
void LBNESteppingAction::CheckInHornEndPlane(const G4Step * theStep)
{

   G4Track * theTrack = theStep->GetTrack();
   
   if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

   std::string preStepPointName  = 
     theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
   std::string postStepPointName  = 
     theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();

   if (((preStepPointName == G4String("TargetHallAndHorn1")) && (postStepPointName == G4String("Tunnel"))) ||
       ((preStepPointName == G4String("Horn2Hall")) && (postStepPointName == G4String("Tunnel")))) 
   {
      theTrack->SetTrackStatus(fStopAndKill);
      G4VTrajectory* currTrajectory = 0;
      currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());

      if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;

      LBNETrajectory* currLBNETrajectory = 0;
      currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);

      if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;

      LBNETrajectory* newtrajectory = 0;
      newtrajectory = new LBNETrajectory(*currLBNETrajectory);
      if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
            
      newtrajectory ->AppendStep(theStep);

      LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      analysis->FillTrackingNtuple(*theTrack, newtrajectory);
      
      delete newtrajectory;
   }


}

//----------------------------------------------------------------------
void LBNESteppingAction::CheckInAlcoveTrackingPlane(const G4Step * theStep)
{

   //G4Track * theTrack = theStep->GetTrack();
   
   if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL || theStep->GetPostStepPoint()->GetPhysicalVolume()==NULL) return;

   G4int pdg = abs(theStep->GetTrack()->GetDefinition()->GetPDGEncoding());
   if (pdg == 12 || pdg == 14 || pdg == 16) return;//Don't save neutrino info

   std::string preStepPointName  = 
     theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
   std::string postStepPointName  = 
     theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();

   if ((preStepPointName == G4String("MuonTrackingVol")) || (postStepPointName == G4String("MuonTrackingVol")))      
   {

      G4VTrajectory* currTrajectory = 0;
      currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());

      if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;

      LBNETrajectory* currLBNETrajectory = 0;
      currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);

      if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;

      LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      analysis->FillAlcoveTrackingPlaneData(*theStep);
      
   }


}

//----------------------------------------------------------------------
void LBNESteppingAction::CheckInTgtEndPlane(const G4Step * theStep)
{

   G4Track * theTrack = theStep->GetTrack();
   
   if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

   std::string preStepPointName  = theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName();

   std::stringstream endName;
   endName << "Horn1TopLevelDown";

   if (preStepPointName.find(endName.str()) != std::string::npos)
   {
      theTrack->SetTrackStatus(fStopAndKill);
      G4VTrajectory* currTrajectory = 0;
      currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());

      if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;

      LBNETrajectory* currLBNETrajectory = 0;
      currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);

      if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;

      LBNETrajectory* newtrajectory = 0;
      newtrajectory = new LBNETrajectory(*currLBNETrajectory);
      if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
            
      newtrajectory ->AppendStep(theStep);

      LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      analysis->FillTrackingNtuple(*theTrack, newtrajectory);
      
      delete newtrajectory;
   }


}
void LBNESteppingAction::OpenAscii(const char *fname) {

   std::cerr << " LBNESteppingAction::OpenAscii Start with filename " << std::string(fname) << std::endl;
   fOutStudy.open(fname);
   if (!fOutStudy.is_open()) {
     std::string descr("Can not open output file "); 
     descr += std::string(fname);
     G4Exception("LBNESteppingAction::OpenAscii", "I/O Problem ", FatalException, descr.c_str());
   }
   if (fStudyGeantinoMode.find("Absorb") != std::string::npos) { 
     fOutStudy << 
    " id x y z ILDecayChan ILHorn1Neck ILHorn2Entr ILNCDecayChan ILNCHorn1Neck ILNCHorn2Entr ILAlHorn2Entr" << std::endl;
   } else if (fStudyGeantinoMode.find("Propa") != std::string::npos) {
     fOutStudy << " id x y z xo yo zo zPost step matPre matPost " << std::endl;
   } else if (fStudyGeantinoMode.find("PropCO") != std::string::npos) {
     fOutStudy << " id x y z xo yo zo step matPre matPost " << std::endl;
   } else if (fStudyGeantinoMode.find("MagTilt") != std::string::npos) {
     fOutStudy << " id x y z xp yp zp  " << std::endl;
   } else if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) {
     fOutStudy << " ik r z nCr  " << std::endl;
   } else if (doStudyParticleThroughHorns) {
     fOutStudy << " evt tr id sNum x y z px py pz w " << std::endl;
   }
   fOutStudy.flush();
   std::cerr << " LBNESteppingAction::OpenAscii " << std::string(fname) << std::endl;
}

void LBNESteppingAction::OpenNtuple(const char *fname) {

   std::cout << " LBNESteppingAction::OpenNtuple Start with filename " << std::string(fname) << std::endl;
   steppingTupleFile = new TFile(fname,"RECREATE");

   if (steppingTupleFile->IsZombie()) {
     std::string descr("Can not open output file "); 
     descr += std::string(fname);
     G4Exception("LBNESteppingAction::OpenNtuple", "I/O Problem ", FatalException, descr.c_str());
   }
   makeSteppingTuple = true;
   steppingTuple = new TTree("steppingTuple","steppingTuple");
   steppingTuple->SetDirectory(steppingTupleFile);
   steppingTuple->Branch("x",&steppingTuple_x,"x/D");
   steppingTuple->Branch("y",&steppingTuple_y,"y/D");
   steppingTuple->Branch("z",&steppingTuple_z,"z/D");
   steppingTuple->Branch("density",&steppingTuple_density,"density/D");

   std::cout << " LBNESteppingAction::OpenNtuple " << std::string(fname) << std::endl;
}

void LBNESteppingAction::StudyAbsorption(const G4Step * theStep) {
//
//make sure we are dealing with a geantino, or a mu, to include lengthening of step due to curling in 
// B Field
//
   G4Track * theTrack = theStep->GetTrack();
   if (((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) && (
        ((theTrack->GetParticleDefinition()->GetParticleName()).find("mu+") == std::string::npos ))) return;

   

        
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
/*
   if ( theTrack->GetNextVolume() == 0 ) {
       fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
        for (size_t k=0; k!=3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
        fOutStudy << " " << totalAbsDecayChan << " " <<  totalAbsHorn1Neck 
                  << " " << totalAbsHorn2Entr;
        fOutStudy << " " << waterAbsDecayChan << " " <<  waterAbsHorn1Neck 
                  << " " << waterAbsHorn2Entr << " " <<  alumAbsHorn2Entr << std::endl;
                  return;
   }              
*/
   //
   // I set the position of the geantino production vertex at Z=0.;
   //
//   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   if (fEvtIdPrevious  != pRunManager->GetCurrentEvent()->GetEventID() ) { 
//     std::cerr << " Evt id " << 
//           pRunManager->GetCurrentEvent()->GetEventID() <<
//            " Starting point, z = " << prePtr->GetPosition()[2] << std::endl;
     totalAbsDecayChan= 0.;
     totalAbsHorn1Neck=0.;
     totalAbsHorn2Entr=0.;
     waterAbsDecayChan= 0.;
     waterAbsHorn1Neck=0.;
     waterAbsHorn2Entr=0.;
     alumAbsHorn2Entr=0.;
     goneThroughHorn1Neck = false;
     goneThroughHorn2Entr = false;
     fEvtIdPrevious = pRunManager->GetCurrentEvent()->GetEventID();
     return;
   } 

   if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
   const double ll = theStep->GetStepLength();
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
/*
   if (postPtr == NULL) {
       fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
        for (size_t k=0; k!=3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
        fOutStudy << " " << totalAbsDecayChan << " " <<  totalAbsHorn1Neck 
                  << " " << totalAbsHorn2Entr;
        fOutStudy << " " << waterAbsDecayChan << " " <<  waterAbsHorn1Neck 
                  << " " << waterAbsHorn2Entr << " " <<  alumAbsHorn2Entr << std::endl;
                  return;
   }              
*/
  if (postPtr == NULL) return;
   G4VPhysicalVolume *physVol = postPtr->GetPhysicalVolume();
   std::string vName(physVol->GetName());
   G4Material *material = theTrack->GetMaterial();
    
      if(makeSteppingTuple) {
      steppingTuple_x = postPtr->GetPosition()[0];
     steppingTuple_y = postPtr->GetPosition()[1];
     steppingTuple_z = postPtr->GetPosition()[2];
     steppingTuple_density = material->GetDensity();
     steppingTuple->Fill();
      }

   if (pRunManager->GetCurrentEvent()->GetEventID() < -3) {
      const double r = std::sqrt(postPtr->GetPosition()[0]*postPtr->GetPosition()[0] + 
                                 postPtr->GetPosition()[1]*postPtr->GetPosition()[1]); 
      const double t = r/std::abs(postPtr->GetPosition()[2] + 515.25); // ZOrigin = -515.25
      // debugging only valid for Zorigin of -515.
      std::cerr << " r = " << r << " theta " << t <<  " z = " << postPtr->GetPosition()[2] << 
              " In " << vName << " material " << material->GetName()
              << " InterLength " << material->GetNuclearInterLength() << std::endl;  
   } 
   if (postPtr->GetPosition()[2] > 500.) goneThroughHorn1Neck=true; // approximate... 
   if (postPtr->GetPosition()[2] > 6360.) goneThroughHorn2Entr=true; //truly approximate. 
   if (ll < 1.0e-10) return; 

   totalAbsDecayChan += ll/material->GetNuclearInterLength(); 
   if (vName.find("WaterLayer") != std::string::npos) waterAbsDecayChan += ll/material->GetNuclearInterLength(); 
   if (!goneThroughHorn1Neck) {
     totalAbsHorn1Neck += ll/material->GetNuclearInterLength(); 
     if (vName.find("WaterLayer") != std::string::npos) waterAbsHorn1Neck += ll/material->GetNuclearInterLength(); 
   }
   if (!goneThroughHorn2Entr) {
     totalAbsHorn2Entr += ll/material->GetNuclearInterLength();
//     if (theTrack->GetTrackLength() < (6000.0*CLHEP::mm)) 
//       std::cerr << " trackLength = " << theTrack->GetTrackLength() << " Z = " << postPtr->GetPosition()[2] << 
//         " Abs L " << totalAbsHorn2Entr << std::endl;
     
     if (vName.find("WaterLayer") != std::string::npos) {
        waterAbsHorn2Entr += ll/material->GetNuclearInterLength(); 
     } else {
       if ((vName.find("Horn1") != std::string::npos) && (material->GetName().find("Alumin") != std::string::npos))
         alumAbsHorn2Entr += ll/material->GetNuclearInterLength(); 
      }
   }
   if (vName.find("DecayPipe") !=  std::string::npos) {
        fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
        for (size_t k=0; k!=3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
        fOutStudy << " " << totalAbsDecayChan << " " <<  totalAbsHorn1Neck 
                  << " " << totalAbsHorn2Entr;
        fOutStudy << " " << waterAbsDecayChan << " " <<  waterAbsHorn1Neck 
                  << " " << waterAbsHorn2Entr << " " <<  alumAbsHorn2Entr << std::endl;
        theTrack->SetTrackStatus(fStopAndKill);
        return;
   } 
    
}
void LBNESteppingAction::StudyPropagation(const G4Step * theStep) {
//
//make sure we are dealing with a geantino... 
//
   G4Track * theTrack = theStep->GetTrack();
   if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (postPtr == 0) return;
   if (theStep->GetStepLength() < 0.1*CLHEP::mm) return;
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());

   if ((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) || 
      (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) {
       std::cout << " LBNESteppingAction::StudyPropagation, critical volume " << std::endl 
          << ".. from " << volNamePre << " to " << volNamePost 
                     << " detected at " <<  prePtr->GetPosition() << " going to " << postPtr->GetPosition() << std::endl;
       const G4VTouchable *preHist = prePtr->GetTouchable();
       const G4VTouchable *postHist = postPtr->GetTouchable();
       const G4int nDepthPre = preHist->GetHistoryDepth();
       const G4int nDepthPost = postHist->GetHistoryDepth();
       std::cerr << " ....  History depth for pre point " <<  nDepthPre << " Post " << nDepthPost << std::endl;
       for (int k=0; k!=  std::max(nDepthPre, nDepthPost); k++) {
         if ((k <nDepthPre) && (k <nDepthPost))
             std::cerr << " ............. Translation at depth, pre ... " << preHist->GetTranslation(k) 
                       << " Post " <<  postHist->GetTranslation(k) <<   std::endl;
         else if (k <nDepthPre)                
             std::cerr << " ..............Translation at depth, pre ... " << preHist->GetTranslation(k) << std::endl;
         else  if (k <nDepthPost)     
             std::cerr << " ..............Translation at depth, post ... " << postHist->GetTranslation(k) << std::endl;
       }
       std::cerr << " ------------------------------------------- " << std::endl;            
   }   
   if (volPre->GetMaterial()->GetName() == volPost->GetMaterial()->GetName()) return;
   fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
   for (int k=0; k != 3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
   for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
   fOutStudy << " " << postPtr->GetPosition()[2];
   fOutStudy << " " << theStep->GetStepLength();
   fOutStudy << " " << volPre->GetMaterial()->GetName();
   fOutStudy << " " << volPost->GetMaterial()->GetName();
   fOutStudy << std::endl;
   G4String vName(volPre->GetName());
//   if (vName.find("DecayPipe") !=  std::string::npos) {
//        theTrack->SetTrackStatus(fStopAndKill);
//    }
}
void LBNESteppingAction::StudyCheckOverlap(const G4Step * theStep) {
//
//make sure we are dealing with a geantino... 
//
   G4Track * theTrack = theStep->GetTrack();
   // Skip this cut for X-checking the muon-genatinos 
//   if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
   if( theTrack->GetNextVolume() == 0 ) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " Out of world with track length  = " << theTrack->GetTrackLength()  << std::endl; 
       return;
   }
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (postPtr == 0) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " No Post Point, Last step at Z = " << prePtr->GetPosition()[2] << " r " <<            
              std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl; 
 
       return;
   }
   if (postPtr->GetPhysicalVolume() == 0) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " No Post Point Volume Last step at Z = " << prePtr->GetPosition()[2] << " r " <<            
              std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl; 
   
    return;
   }
   if (prePtr->GetPhysicalVolume() == 0) return;
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   if (volPre == 0) return;
   if (volPost == 0) return;
  std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
   std::string matNamePre = volPre->GetMaterial()->GetName();
   std::string matNamePost = volPost->GetMaterial()->GetName();
   //
   if (pRunManager->GetCurrentEvent()->GetEventID() < 3) {
     std::cout << " at Z = " << prePtr->GetPosition()[2] << " r " << 
           std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) 
           << ", " << volNamePre  << " to " << postPtr->GetPosition() << ", " 
           << volNamePost << " mat from " << matNamePre << " to " << matNamePost;
           if (std::abs(theTrack->GetParticleDefinition()->GetPDGEncoding()) == 13) {
             std::cout << " Ek " << theTrack->GetKineticEnergy()/CLHEP::GeV << std::endl;
           } else { std::cout << std::endl; }
   /* These debugging statement may cause a fatal exception, because A and/or Z not defined for mixtures.. 
           if (matNamePost.find("Targ") != std::string::npos) {
             std::cout << " Target material found, Z " << volPost->GetMaterial()->GetZ() << " A " 
                      << volPost->GetMaterial()->GetA() << std::endl;
           }  
           if (matNamePost.find("DecayPipe") != std::string::npos) {
             std::cout << " Decay Pipe material found, Z " << volPost->GetMaterial()->GetZ() << " A " 
                      << volPost->GetMaterial()->GetA()/mole << " density " << volPost->GetMaterial()->GetDensity()/(g/cm3) << std::endl;
           }
    */  
    }
//   if (((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) || 
//      (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) &&
//      ( (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos) || 
//      (volNamePre.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos))) {
   if ( (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos) &&
        (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos)) {
     fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
     for (int k=0; k != 3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
     for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
     fOutStudy << " " << theStep->GetStepLength();
     fOutStudy << " " << volPre->GetMaterial()->GetName();
     fOutStudy << " " << volPost->GetMaterial()->GetName();
     fOutStudy << std::endl;
  }
}

void LBNESteppingAction::StudyCheckMagneticTilts(const G4Step * theStep) {
//
//make sure we are dealing with a geantino... 
//
   G4Track * theTrack = theStep->GetTrack();
   if ((theTrack->GetParticleDefinition()->GetParticleName()).find("mu") == std::string::npos) return;
   if( theTrack->GetNextVolume() == 0 ) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " Out of world with track length  = " << theTrack->GetTrackLength()  << std::endl; 
       return;
   }
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
   if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
     std::cout << " at Z = " << prePtr->GetPosition()[2] << " r " << 
           std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (postPtr == 0) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " No Post Point, Last step at Z = " << prePtr->GetPosition()[2] << " r " <<            
              std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl; 
 
       return;
   }
   if (postPtr->GetPhysicalVolume() == 0) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " No Post Point Volume Last step at Z = " << prePtr->GetPosition()[2] << " r " <<            
              std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
                 prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl; 
   
    return;
   }
   if (prePtr->GetPhysicalVolume() == 0) return;
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   if (volPre == 0) return;
   if (volPost == 0) return;
  std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
//   if (((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) || 
//      (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) &&
//      ( (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos) || 
//      (volNamePre.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos))) {
   if ( (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos) &&
        (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos)) {
     fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
     for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
     for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetMomentumDirection()[k];
     fOutStudy << std::endl;
  }
}

void LBNESteppingAction::GenerateVolumeCrossingsRZMap(const G4Step * theStep) {
//
//make sure we are dealing with a geantino... 
//
   G4Track * theTrack = theStep->GetTrack();
   if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
   if( theTrack->GetNextVolume() == 0 ) {
        if (pRunManager->GetCurrentEvent()->GetEventID() < 3) 
           std::cout << " Out of world with track length  = " << theTrack->GetTrackLength()  << std::endl; 
       return;
   }
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (postPtr == 0) return;
   if (postPtr->GetPhysicalVolume() == 0) return;
   if (prePtr->GetPhysicalVolume() == 0) return;
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   if (volPre == 0) return;
   if (volPost == 0) return;
   double zPost = postPtr->GetPosition()[2];
   double rPost = std::sqrt(postPtr->GetPosition()[0]*postPtr->GetPosition()[0] + 
                                  postPtr->GetPosition()[1]*postPtr->GetPosition()[1]);
   zPost = std::min(zPost, fGenerateVoxelsZMax);
   zPost = std::max(zPost, fGenerateVoxelsZScale);
   rPost = std::min(rPost, fGenerateVoxelsRMax );
   rPost = std::max(rPost, 1.);
   double lnR = std::log(rPost);
   int iz = static_cast<int>((zPost + fGenerateVoxelsZOffset)/fGenerateVoxelsZScale);
   int ir = static_cast<int>(lnR*fGenerateVoxelsRScale);
   int ik = iz*fGenerateVoxelsIRKMax + ir;
   std::map<int, int>::iterator itMRZ =  fRZVoxelsData.find(ik);
   if(itMRZ ==  fRZVoxelsData.end()) 
     fRZVoxelsData.insert(std::pair<int, int>(ik, 1));
   else 
     itMRZ->second += 1;
}

void LBNESteppingAction::dumpStepCheckVolumeAndFields(const G4Step * theStep) {

  int idEvt = pRunManager->GetCurrentEvent()->GetEventID();
  if (fEvtIdPrevious != idEvt) {
      if (fOutStepStudy.is_open()) fOutStepStudy.close();
      std::ostringstream fNameStrStr; fNameStrStr << "./StepSudiesMagn_Evt" << idEvt << ".txt";
      std::string fNameStr(fNameStrStr.str());
      fOutStepStudy.open(fNameStr.c_str());
      fOutStepStudy << " x y z xp yp dpt bphi bz radIC radOC factR vName " << std::endl;
      fEvtIdPrevious = idEvt;
  }
  G4StepPoint* prePtr = theStep->GetPreStepPoint();
  if (theStep->GetPreStepPoint()->GetPhysicalVolume() == 0) return;
  if (prePtr == 0) return;
  G4Track * theTrack = theStep->GetTrack();
  if ( theTrack->GetNextVolume() == 0 ) {
    if (fOutStepStudy.is_open()) fOutStepStudy.close();
    return;
  }
  G4StepPoint* postPtr = theStep->GetPostStepPoint();
  if (postPtr == 0) {
    if (fOutStepStudy.is_open()) fOutStepStudy.close();
    return;
  } 
  for (size_t k=0; k !=3; k++) fOutStepStudy << " " << postPtr->GetPosition()[k];
  G4ThreeVector postVec = postPtr->GetMomentum();
  G4ThreeVector preVec = prePtr->GetMomentum();
  fOutStepStudy << " " << postVec[0]/postVec[2];
  fOutStepStudy << " " << postVec[1]/postVec[2];
  const double ptPre = std::sqrt(preVec[0]*preVec[0] + preVec[1]*preVec[1])/CLHEP::GeV;
  const double ptPost = std::sqrt(postVec[0]*postVec[0] + postVec[1]*postVec[1])/CLHEP::GeV;
  fOutStepStudy << " " << ptPost - ptPre;
  // Use our our utility to navigate to find the right field
  const G4Field *aField = 0;
  const LBNEMagneticFieldHorn *myField=0;
  const LBNFMagneticFieldPolyconeHorn *myFieldPH=0;
  double effInnerRad = 0; double effOuterRad = 0; double effSkinD = 0.;
  LBNEVolumePlacements *aPlacementHandler= LBNEVolumePlacements::Instance();
  if ((aPlacementHandler->GetNumberOfHornsPolycone() == 0) &&
      (postPtr->GetPosition()[2] > -500.) && (postPtr->GetPosition()[2] < 5000.)) {
    const LBNEVolumePlacementData *pl = 
        aPlacementHandler->Find("FieldHorn1", "Horn1PolyM1", "DetectorConstruction");
    aField = pl->fCurrent->GetFieldManager()->GetDetectorField();
    if (aField != 0) {
      myField = dynamic_cast<const LBNEMagneticFieldHorn *> (aField); 
      if (myField != 0) {
         effInnerRad = myField->GetEffectiveInnerRad(); 
         effOuterRad = myField->GetEffectiveOuterRad();
         effSkinD = myField->GetEffectiveSkinDepthFact();
      }
    } 
  } else if ((aPlacementHandler->GetNumberOfHornsPolycone() == 0) && (postPtr->GetPosition()[2] < 13000.)) {
    const LBNEVolumePlacementData *pl = 
        aPlacementHandler->Find("FieldHorn2", "Horn2Hall", "DetectorConstruction");
    aField = pl->fCurrent->GetFieldManager()->GetDetectorField();
    myField = dynamic_cast<const LBNEMagneticFieldHorn *> (aField); 
    if (myField != 0) { 
      effInnerRad =myField->GetEffectiveInnerRad(); 
      effOuterRad = myField->GetEffectiveOuterRad();
       effSkinD = myField->GetEffectiveSkinDepthFact();
    }
  } else if ((aPlacementHandler->GetNumberOfHornsPolycone() != 0) && (postPtr->GetPosition()[2] < 22000.)) {
    G4FieldManager *aFieldMgr = postPtr->GetPhysicalVolume()->GetLogicalVolume()->GetFieldManager();
    if (aFieldMgr != 0) {
      aField = aFieldMgr->GetDetectorField();
       myFieldPH = dynamic_cast<const LBNFMagneticFieldPolyconeHorn *> (aField); 
       if (myFieldPH != 0) {
         effInnerRad =myFieldPH->GetEffectiveInnerRad(); 
         effOuterRad = myFieldPH->GetEffectiveOuterRad();
         effSkinD = myFieldPH->GetEffectiveSkinDepthFact();
       }
    }
  }
   
  if (aField == 0) {
   fOutStepStudy << " 0. 0. 0. 0. 0. "; 
  } else {
   double pos[3]; for (size_t k=0; k!= 3; ++k) pos[k] = postPtr->GetPosition()[k];
   double bf[3];
   if (myField != 0) myField->GetFieldValue(pos, &bf[0]);
   if (myFieldPH != 0) {
         myFieldPH->GetFieldValue(pos, &bf[0]); // We repeat the checks on effective radii.. 
         effInnerRad =myFieldPH->GetEffectiveInnerRad(); 
         effOuterRad = myFieldPH->GetEffectiveOuterRad();
         effSkinD = myFieldPH->GetEffectiveSkinDepthFact();
  } 
  if((myField == 0) && (myFieldPH == 0)) aField->GetFieldValue(pos, &bf[0]);
   fOutStepStudy << " " << std::sqrt(bf[0]*bf[0] + bf[1]*bf[1])/CLHEP::tesla;
   fOutStepStudy << " " << bf[2]/CLHEP::tesla;
   fOutStepStudy << " " << effInnerRad << " " << effOuterRad << " " << effSkinD; 
  } 
  fOutStepStudy << " " << postPtr->GetPhysicalVolume()->GetLogicalVolume()->GetName() << std::endl;
}

void LBNESteppingAction::StudyParticleThroughHorns(const G4Step* theStep) {

   G4Track * theTrack = theStep->GetTrack();
   if (theTrack == 0) return;
   //
   // Simplified version for animation with muons...
   //
   /*
   if (theTrack->GetPosition()[2] > 15000.) return;
   G4ThreeVector direction = theTrack->GetMomentumDirection();
   const double eTot = theTrack->GetTotalEnergy();
   const double mass= theTrack->GetParticleDefinition()->GetPDGMass();
   const double pTot = std::sqrt(eTot*eTot - mass*mass);
   fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
   fOutStudy << " " << theTrack->GetParticleDefinition()->GetPDGEncoding();
   const int stepNum = theTrack->GetCurrentStepNumber();
   fOutStudy << " " << stepNum;
   if (stepNum < 2) {
     for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetVertexPosition()[k];
     for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*theTrack->GetVertexMomentumDirection()[k];
     fOutStudy <<  std::endl;   
     fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID(); 
     fOutStudy << " " << theTrack->GetParticleDefinition()->GetPDGEncoding();
     fOutStudy << " " << stepNum;
   }
   for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetPosition()[k];
   for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*direction[k];
   fOutStudy <<  std::endl;   
   return;
   */
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (prePtr->GetPhysicalVolume() == 0) return;
   if (postPtr == 0) return;
   if (postPtr->GetPhysicalVolume() == 0) return;
   //
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
   bool ofInterest = false;
   // Entering the corridor 
   if ((volNamePost == std::string("Tunnel")) && 
       (volNamePre.find("Horn1") != std::string::npos )) ofInterest=true;
   if ((volNamePre.find("Horn2") != std::string::npos) && 
       (volNamePost == std::string("Tunnel"))) ofInterest=true;
   if ((volNamePre == std::string("Tunnel")) && 
       (volNamePost == std::string("DecayPipeHall"))) ofInterest=true;
   if (!ofInterest) return;
   const double eTot = theTrack->GetTotalEnergy();
   const double mass= theTrack->GetParticleDefinition()->GetPDGMass();
   const double pTot = std::sqrt(eTot*eTot - mass*mass);
   if (pTot < 0.1*CLHEP::GeV) return;
   if (pTot > 50.*CLHEP::GeV) return;
   int pdgInfo = theTrack->GetParticleDefinition()->GetPDGEncoding();
   // check only the muons for now. 
   if (std::abs(pdgInfo) != 211) return;
   fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
   fOutStudy << " " <<  theTrack->GetTrackID(); 
   const int stepNum = theTrack->GetCurrentStepNumber();
   fOutStudy  << " " << pdgInfo << " " << stepNum;
   for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
   G4ThreeVector direction = theTrack->GetMomentumDirection();
   for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*direction[k];
   LBNETrackInformation* info = reinterpret_cast<LBNETrackInformation*>(theTrack->GetUserInformation());
   fOutStudy << " "  <<  info->GetNImpWt();  // Importance weight
   fOutStudy << std::endl;

}
void LBNESteppingAction::StudyPionsThroughHorn2 (const G4Step* theStep) {
  const G4Track * theTrack = theStep->GetTrack();
  if (theTrack == 0) return;
  int pID = theTrack->GetParticleDefinition()->GetPDGEncoding();
  if (std::abs(pID) != 211) return;  
  G4StepPoint* prePtr = theStep->GetPreStepPoint();
  if (prePtr == 0) return;
  if (prePtr->GetPhysicalVolume() == 0) return;
  G4StepPoint* postPtr = theStep->GetPostStepPoint();
  if (postPtr == 0) return;
  if (postPtr->GetPhysicalVolume() == 0) return;
  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
  std::string volNamePost(volPost->GetName());
  std::string volNamePre(volPre->GetName());
  if ((volNamePre == std::string("Horn2Part3Ring")) && 
       (volNamePost.find("Horn2") != std::string::npos)) {
      LBNEQuickPiToNuVect* ptr = LBNEQuickPiToNuVect::Instance();
      int sign = (pID > 0) ? 1 : -1;
      ptr->RingHangerPion(theTrack->GetTrackID(), sign, theTrack->GetMomentum(), postPtr->GetPosition()); 
   }
}

void LBNESteppingAction::CheckInTrackingDetectorPlane(const G4Step *theStep)
{
  //  std::cout << "I am here" << std::endl;

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
    theStep->GetPostStepPoint()->GetPhysicalVolume()){
    G4LogicalVolume *preStepVolume =
    theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    G4LogicalVolume *postStepVolume =
    theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    if(postStepVolume == TrkPlnLogical){
      if(preStepVolume != TrkPlnLogical){
        G4VTrajectory* currTrajectory = 0;
        currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());

        if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;

        LBNETrajectory* currLBNETrajectory = 0;
        currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;

        LBNETrajectory* newtrajectory = 0;
        newtrajectory = new LBNETrajectory(*currLBNETrajectory);
        if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;

        newtrajectory ->AppendStep(theStep);

        LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
        //analysis->FillTrackingNtuple(*theTrack, newtrajectory);
        analysis->FillTrackingPlaneData(*theStep);

        delete newtrajectory;
      }
    }
  }
}


void LBNESteppingAction::CheckInTargetOutput(const G4Step *theStep)
{

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
std::string preStepPointName =
theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();
//std::string postStepPointName = 
//theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
std::stringstream endName, posvolname;
endName<<"Target";
posvolname<<"TargetHallandHorn1";
if(preStepPointName.find(endName.str())!=std::string::npos)
//if(preStepPointName==postStepPointName) return;
//if(preStepPointName!=postStepPointName)
{
//theTrack->SetTrackStatus(fStopAndKill);
G4VTrajectory* currTrajectory = 0;
      currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());

      if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;

      LBNETrajectory* currLBNETrajectory = 0;
      currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);

      if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;

      LBNETrajectory* newtrajectory = 0;
      newtrajectory = new LBNETrajectory(*currLBNETrajectory);
      if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
            
      newtrajectory ->AppendStep(theStep);

      LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      analysis->FillTargetOutputData(*theStep);
      
      delete newtrajectory;
}
}       


// Amit Bashyal
// Does this still work ? Paul Lebrun, May 26 2016. 
void LBNESteppingAction::CheckInTrackingDetectorH1Plane(const G4Step *theStep)
{
  //  std::cout << "I am here" << std::endl;
//  G4Track * theTrack = theStep->GetTrack();

  if ((theStep->GetPreStepPoint() == 0) || (theStep->GetPostStepPoint() == 0)) return;
  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
    theStep->GetPostStepPoint()->GetPhysicalVolume()){
    G4LogicalVolume *preStepVolume =
    theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    G4LogicalVolume *postStepVolume =
   theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
   if(postStepVolume == TrkPlnH1Logical){
      if(preStepVolume != TrkPlnH1Logical){
        LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
      //  analysis->FillTrackingPlaneH1Data(*theTrack, newtrajectory);
        // std::cout<<"Filling H1Trackingplane Ntuples"<<std::endl;
        analysis->FillTrackingPlaneH1Data(*theStep);
     }
    }
  }

}
//Amit Bashyal
void LBNESteppingAction::CheckInTrackingDetectorH2Plane(const G4Step *theStep)
{
  //  std::cout << "I am here" << std::endl;
  if ((theStep->GetPreStepPoint() == 0) || (theStep->GetPostStepPoint() == 0)) return;
  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
    theStep->GetPostStepPoint()->GetPhysicalVolume()){
    G4LogicalVolume *preStepVolume =
    theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    G4LogicalVolume *postStepVolume =
    theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    if(postStepVolume == TrkPlnH2Logical){
      if(preStepVolume != TrkPlnH2Logical){
        LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
        analysis->FillTrackingPlaneH2Data(*theStep);
      }
    }
  }
}
//Amit Bashyal
void LBNESteppingAction::CheckInTrackingDetectorDPPlane(const G4Step *theStep)
{
  //  std::cout << "I am here" << std::endl;
  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;

  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
    theStep->GetPostStepPoint()->GetPhysicalVolume()){
    G4LogicalVolume *preStepVolume =
    theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    G4LogicalVolume *postStepVolume =
    theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
    if(postStepVolume == TrkPlnDPLogical){
      if(preStepVolume != TrkPlnDPLogical){
        LBNEAnalysis* analysis = LBNEAnalysis::getInstance();
        analysis->FillTrackingPlaneDPData(*theStep);
      }
    }
  }
}
void LBNESteppingAction::InitiateHadronFluxFromTargetASCII() const {
   doGenerateHadronFluxFromTargetASCII=true;
   // 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("./hadronFluxFromTargetASCII.txt");
   if ((clusterNum != 0) || (procNum != 0)) {
//     const char *userEnv = getenv("USER"); // assume it always works
//     std::string aUserStr(userEnv); 
     std::ostringstream fNStrStr; 
     fNStrStr << "./hadronFluxFromTargetASCII_" 
              << clusterNum << "_" << procNum << ".txt";
     fName = fNStrStr.str();
   }
   if (fStrHadronFluxFromTargetASCII.is_open()) fStrHadronFluxFromTargetASCII.close();
   fStrHadronFluxFromTargetASCII.open(fName.c_str());
   fStrHadronFluxFromTargetASCII << " evt tr id x y z px py pz w " << std::endl;
   std::cerr << " LBNESteppingAction::InitiateHadronFluxFromTargetASCII, Opening file " << 
        fName << std::endl;
   if (doGenerateHadronFluxFromTargetASCII) 
     std::cerr << " O.K., ready for HadronFluxFromTargetASCII.... in... " << this << std::endl;
}
void LBNESteppingAction::FillHadronFluxFromTargetASCII(const G4Step *theStep) const {

   G4Track * theTrack = theStep->GetTrack();
   if (theTrack->GetPosition()[2] > 3500.) return; // Hopefully, the target is not longer.. 
//   std::cerr << " LBNESteppingAction::FillHadronFluxFromTargetASCII At z  " << theTrack->GetPosition()[2] << std::endl;
   G4StepPoint* prePtr = theStep->GetPreStepPoint();
   if (prePtr == 0) return;
   G4StepPoint* postPtr = theStep->GetPostStepPoint();
   if (postPtr == 0) return;
   if (postPtr->GetPhysicalVolume() == 0) return;
   if (prePtr->GetPhysicalVolume() == 0) return;
   G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
   std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
//   std::cerr << " LBNESteppingAction::FillHadronFluxFromTargetASCII volNamePre " 
//      << volNamePre << " Post " << volNamePost << std::endl;
//   if (((volNamePost.find("TargetNoSplitHelium") != std::string::npos) || 
//        (volNamePost.find("TargetNoSplitHeCont") != std::string::npos)) && 
//       (volNamePre.find("TargetUsptreamSimpleCyl") != std::string::npos)) {
     if ((volNamePre.find("TargetNoSplitHelium") != std::string::npos) && 
         (volNamePost.find("TargetNoSplitHeContainer") != std::string::npos)) {
     double pTotSq = 0; for (int k=0; k != 3; k++) pTotSq += theTrack->GetMomentum()[k]*theTrack->GetMomentum()[k];
     if (pTotSq < 2500) return;
     G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
     fStrHadronFluxFromTargetASCII << " " << evtno << " " << theTrack->GetTrackID() << " " 
                                   <<  theTrack->GetParticleDefinition()->GetPDGEncoding();
    
     for (int k=0; k != 3; k++) fStrHadronFluxFromTargetASCII << " " << postPtr->GetPosition()[k];
     for (int k=0; k != 3; k++) fStrHadronFluxFromTargetASCII << " " << (theTrack->GetMomentum()[k])/CLHEP::GeV;
     LBNETrackInformation* info = (LBNETrackInformation*)(theTrack->GetUserInformation());
     fStrHadronFluxFromTargetASCII << " " << info->GetNImpWt(); 
     fStrHadronFluxFromTargetASCII << std::endl;
     //
     // To speed up the detection of missing pion (with respect to FLUKA)
     //
//     theTrack->SetTrackStatus(fStopAndKill);
   }
}


void LBNESteppingAction::StudyMuonSculptedAbsorberFlux(const G4Step* aStep) {
// June 16 2015:  add the capability to generate 
  if (!fOutMuonSculptedAbsorber.is_open()){
    //fOutMuonSculptedAbsorber.open("./MuonFluxAtSculptedAbsorber.txt");
    G4cout <<"Opening muon file!!!!! "<< fOutMuonSculptedAbsorberFilename<<G4endl;
    fOutMuonSculptedAbsorber.open(fOutMuonSculptedAbsorberFilename);
    fOutMuonSculptedAbsorber << " evt tr id slice x y z px py pz ek w g2 " << std::endl;     
  }
  const G4Track* aTrack = aStep->GetTrack();
  int pdg = aTrack->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
  if (std::abs(pdg) != 13) return;
  const G4StepPoint *prePt = aStep->GetPreStepPoint(); 
  if (prePt == 0) return;
  if (prePt->GetPhysicalVolume() == 0) return;
  if (prePt->GetPosition()[2] < 200000.) return; // Rough cut!. 
  const G4StepPoint *postPt = aStep->GetPostStepPoint(); 
  if (postPt->GetPhysicalVolume() == 0) return;  
   G4LogicalVolume *volPost = postPt->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePt->GetPhysicalVolume()->GetLogicalVolume();
   std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
   //
   // debugging ... 
   // 
//   if ((volNamePre.find("HadronAbsorberSculpBackWall") == 0) || (volNamePre.find("HadronAbsorberSculpEndIronAC") == 0)) {
//     std::cerr << " Got the end of HA, vol name Pre " <<  volNamePre << " .. And quit " << std::endl; exit(2);
//   }
   int aSlice=-100;
   bool ofInterestLBNE = ((volNamePre == std::string("Tunnel")) && (volNamePost == std::string("HadronAbsorberTop"))) ||
                         ((volNamePre == std::string("AH_Muon_alkair")) && (volNamePost == std::string("detLogical")));
                         
   bool ofInterestLBNF = ((volNamePre == std::string("Tunnel")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) ||
                     ((volNamePost == std::string("Tunnel")) && (volNamePre == std::string("HadronAbsorberSculpTop"))) ||
                     (volNamePre.find("HadronAbsorberSculp") != std::string::npos);
   if ((!ofInterestLBNF) && (!ofInterestLBNE)) return; 
   if (volNamePre == std::string("Tunnel")) { 
     aSlice = -2;
     fThisParticleGotToHAFront = true;
   }     
   if (ofInterestLBNE) {
     if (volNamePre == std::string("AH_Muon_alkair")) aSlice = 30;  
   } else { 
     if ((volNamePre == std::string("HadronAbsorberSculpAibufferDiag")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) {
       aSlice = -1;
     } else if ((volNamePre == std::string("HadronAbsorberSculpSpoilerAir")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) {
       aSlice = 0;
     } else if (volNamePost == std::string("HadronAbsorberSculpTop")) {
       size_t iPosU = volNamePre.find("-");
       if (iPosU != std::string::npos) {
         iPosU++;
         std::string tSliceStr = volNamePre.substr(iPosU, std::string::npos);
         aSlice = atoi(tSliceStr.c_str());
       }
       if  (volNamePre == std::string("HadronAbsorberSculpBackWall"))   aSlice = 30; 
       if  (volNamePre == std::string("HadronAbsorberSculpEndIronAC"))  aSlice = 29;
     }
   }
   if (aSlice < -5) return;
   G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
   fOutMuonSculptedAbsorber << " " << evtno << " " << aTrack->GetTrackID() << " " << pdg << " " << aSlice;
   for (size_t k=0; k!=3; k++) fOutMuonSculptedAbsorber << " " << prePt->GetPosition()[k];
   for (size_t k=0; k!=3; k++) fOutMuonSculptedAbsorber << " " << aTrack->GetMomentum()[k]/CLHEP::GeV;
   LBNETrackInformation* info = (LBNETrackInformation*)(aTrack->GetUserInformation());
   fOutMuonSculptedAbsorber << " " << aTrack->GetDynamicParticle()->GetKineticEnergy()/CLHEP::GeV << " " 
                            << info->GetNImpWt();
   if (fThisParticleGotToHAFront) fOutMuonSculptedAbsorber << " 1";
   else fOutMuonSculptedAbsorber << " 0";
   fOutMuonSculptedAbsorber << std::endl;
}
void LBNESteppingAction::CheckHadronsMarsCmpApr2017(const G4Step* aStep) {

   size_t kSelOut = 99;
   const G4StepPoint *prePt = aStep->GetPreStepPoint(); 
   if (prePt == 0) return;
   if (prePt->GetPhysicalVolume() == 0) return;
   const G4StepPoint *postPt = aStep->GetPostStepPoint(); 
   if (postPt == 0) return;
   if (postPt->GetPhysicalVolume() == 0) return;
   G4LogicalVolume *volPost = postPt->GetPhysicalVolume()->GetLogicalVolume();
   G4LogicalVolume *volPre = prePt->GetPhysicalVolume()->GetLogicalVolume();
   std::string volNamePost(volPost->GetName());
   std::string volNamePre(volPre->GetName());
   if ((volNamePre == std::string("TargetNoSplitHeContainer")) && 
       (volNamePost == std::string("TargetNoSplitM1"))) kSelOut = 0;
   else if ((volNamePre == std::string("Horn1PolyM1")) && 
       (volNamePost == std::string("TargetHallAndHorn1"))) kSelOut = 1;
   else if ((volNamePre == std::string("LBNFSimpleHorn2Container")) && 
       (volNamePost == std::string("Tunnel"))) kSelOut = 2;
   else if ((volNamePre == std::string("LBNFSimpleHorn3Container")) && 
       (volNamePost == std::string("Tunnel"))) kSelOut = 3;
    
   if (kSelOut == 99) return;  
   G4Track* aTrack = aStep->GetTrack();
  int pdg = aTrack->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
   bool isMeson = (std::abs(pdg) == 211) || (std::abs(pdg) == 321);
   if (!isMeson) return;
   G4ThreeVector xyz = postPt->GetPosition();
   G4ThreeVector pMom = postPt->GetMomentum();
   G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
   (*fOutPtrsForMarsCmpApr2017[kSelOut]) << " " << evtno << " " << aTrack->GetTrackID() << " "  << pdg << " " 
                                         << xyz[0]/10. << " " << xyz[1]/10. << " " << xyz[2]/10. << " " 
                                         << pMom[0]/1000. << " " << pMom[1]/1000. << " " << pMom[2]/1000. << " "
                                         << postPt->GetKineticEnergy()/1000. << std::endl;
   if (kSelOut == 3) aTrack->SetTrackStatus(fStopAndKill);
}