LBNEAnalysis Class Reference

#include <LBNEAnalysis.hh>

Public Member Functions
	LBNEAnalysis ()
	~LBNEAnalysis ()
G4bool	CreateOutput ()
G4bool	CreateDk2NuOutput ()
void	CloseOutput ()
void	CloseDk2NuOutput ()
void	FillNeutrinoNtuple (const G4Track &track, const std::vector< G4VTrajectory * > &nuHistory)
void	FillTrackingNtuple (const G4Track &track, LBNETrajectory *currTrajectory)
void	TrackThroughGeometry (const LBNETrajectory *TrackTrajectory)
void	fillDkMeta ()
void	FillAlcoveTrackingPlaneData (const G4Step &aStep)
void	FillTrackingPlaneData (const G4Step &aStep)
void	FillTrackingPlaneH1Data (const G4Step &aStep)
void	FillTrackingPlaneH2Data (const G4Step &aStep)
void	FillTrackingPlaneDPData (const G4Step &aStep)
void	FillTargetOutputData (const G4Step &aStep)
LBNETrajectory *	GetParentTrajectory (G4int parentID)
LBNETrajectory *	GetTrajectory (G4int trackID)
void	CalcLocationWeights (const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu, bool useRealisticNearDetectorVolume)
void	ResetEvent ()
void	FillEvent ()
void	SetCount (G4int count)
G4int	GetCount ()
void	SetEntry (G4int entry)
G4int	GetEntry ()
G4double	GetDistanceInVolume (LBNETrajectory *wanted_traj, G4String wanted_vol)

Static Public Member Functions
static LBNEAnalysis *	getInstance ()

Private Member Functions
void	setDetectorPositions ()

Private Attributes
bool	fDk2NuDetectorFileRead
G4double	x
G4double	y
G4double	z
int	fTrackID [kMaxP]
int	fParticlePDG [kMaxP]
double	fNImpWt [kMaxP]
double	fParticleX [kMaxP]
double	fParticleY [kMaxP]
double	fParticleZ [kMaxP]
double	fParticleEnergy [kMaxP]
double	fParticleMass [kMaxP]
double	fParticleDX [kMaxP]
double	fParticleDY [kMaxP]
double	fParticleDZ [kMaxP]
int	fNParticles
int	fH1TrackID
int	fH1ParticlePDG
double	fH1NImpWt
int	fH1DParticlePDG
int	fH1DTrackID
double	fH1ParticleX
double	fH1ParticleY
double	fH1ParticleZ
double	fH1ParticleEnergy
double	fH1ParticleMass
double	fH1ParticleDX
double	fH1ParticleDY
double	fH1ParticleDZ
double	fH1ParticlePX
double	fH1ParticlePY
double	fH1ParticlePZ
double	fH1ParticlePXPZ
double	fH1ParticlePYPZ
double	fH1PProductionX
double	fH1PProductionY
double	fH1PProductionZ
double	fH1PProductionDX
double	fH1PProductionDY
double	fH1PProductionDZ
std::string	fH1CProcess
G4ThreeVector	NuMomentum
int	fH1NParticles
int	fH2TrackID
int	fH2ParticlePDG
double	fH2NImpWt
double	fH2ParticleX
double	fH2ParticleY
double	fH2ParticleZ
double	fH2ParticleEnergy
double	fH2ParticleMass
double	fH2ParticleDX
double	fH2ParticleDY
double	fH2ParticleDZ
double	fH2ParticlePX
double	fH2ParticlePY
double	fH2ParticlePZ
double	fH2ParticlePXPZ
double	fH2ParticlePYPZ
double	fH2PProductionX
double	fH2PProductionY
double	fH2PProductionZ
double	fH2PProductionDX
double	fH2PProductionDY
double	fH2PProductionDZ
G4String	fH2CProcess
int	fH2CNProcess
int	fH2NParticles
int	fDPTrackID
int	fDPParentID
int	fDPParticlePDG
double	fDPNImpWt
double	fDPParticleX
double	fDPParticleY
double	fDPParticleZ
double	fDPParticleEnergy
double	fDPParticleMass
double	fDPParticleDX
double	fDPParticleDY
double	fDPParticleDZ
double	fDPParticlePX
double	fDPParticlePY
double	fDPParticlePZ
double	fDPParticlePXPZ
double	fDPParticlePYPZ
double	fDPPProductionX
double	fDPPProductionY
double	fDPPProductionZ
double	fDPPProductionDX
double	fDPPProductionDY
double	fDPPProductionDZ
int	fDPNParticles
int	fTTrackID
int	fTParentID
double	fTParticleX
double	fTParticleY
double	fTParticleZ
int	fTParticlePDG
double	fTParticleEnergy
double	fTParticlePX
double	fTParticlePY
double	fTParticlePZ
int	fTNParticles
G4double	noProtons
int	fMuRunNo
int	fMuEvtNo
int	fMuNParticles
int	fMuTrackID [kMaxP]
int	fMuParentID [kMaxP]
int	fMuPDG [kMaxP]
double	fMuNimpWt [kMaxP]
double	fMuX [kMaxP]
double	fMuY [kMaxP]
double	fMuZ [kMaxP]
double	fMuT [kMaxP]
double	fMuStartX [kMaxP]
double	fMuStartY [kMaxP]
double	fMuStartZ [kMaxP]
double	fMuStartE [kMaxP]
double	fMuMass [kMaxP]
double	fMuEnergy [kMaxP]
double	fMuPX [kMaxP]
double	fMuPY [kMaxP]
double	fMuPZ [kMaxP]
double	fMuTheta [kMaxP]
double	fMudEdx [kMaxP]
double	fMudEdx_ion [kMaxP]
double	fMuDE [kMaxP]
double	fMuDEion [kMaxP]
int	fMuNSteps [kMaxP]
int	fMuParPDG [kMaxP]
double	fMuParX [kMaxP]
double	fMuParY [kMaxP]
double	fMuParZ [kMaxP]
double	fMuParE [kMaxP]
double	fMuParPX [kMaxP]
double	fMuParPY [kMaxP]
double	fMuParPZ [kMaxP]
TTree *	fHorn1TrackingTree
TTree *	fHorn2TrackingTree
TTree *	fDPTrackingTree
TTree *	fTargetOutputTree
TTree *	fAlcoveTrackingTree
char	nuNtupleFileName [1024]
char	nuNtupleFileNameDK2Nu [1024]
std::map< int, int >	code
TFile *	fOutFile
TTree *	fOutTree
TFile *	nuNtuple
TFile *	fOutFileDk2Nu
TTree *	fOutTreeDk2Nu
TTree *	fOutTreeDk2NuMeta
TTree *	fTrackingTree
bsim::Dk2Nu *	fDk2Nu
bsim::DkMeta *	fDkMeta
LBNEDataNtp_t *	fLBNEOutNtpData
LBNEDataNtp_t *	fTrackingPlaneData
G4int	fcount
G4int	fentry
G4int	nGenAbs = 3
G4int	nVolAbs = 4
std::vector< G4double >	fXdet_near
std::vector< G4double >	fYdet_near
std::vector< G4double >	fZdet_near
std::vector< G4double >	fXdet_far
std::vector< G4double >	fYdet_far
std::vector< G4double >	fZdet_far
std::vector< G4String >	fDetNameNear
std::vector< G4String >	fDetNameFar
std::vector< bsim::Traj >	vec_traj
int	tar_trackID =-1
int	dk_trackID
int	nVintTot
int	nVdblTot
G4double	dist_IC1 [3]
G4double	dist_IC2 [3]
G4double	dist_DPIP [3]
G4double	dist_DVOL [3]
vstring_t	VolVdblName
vstring_t	VolAbsName
vstring_t	VolVintName
vstring_t	GenAbsName
bool	doPPFx
std::ofstream	fOutDBGDk2nu_
int	numEntFillNTu_
int	numEntFillNTuPzR_
int	numEntFillNTuDk2nuPzR_

Static Private Attributes
static LBNEAnalysis *	instance = 0

Detailed Description

Definition at line 35 of file LBNEAnalysis.hh.

Constructor & Destructor Documentation

LBNEAnalysis::LBNEAnalysis

(

)

Definition at line 57 of file LBNEAnalysis.cc.

                           :
    fDk2NuDetectorFileRead(false),
    fHorn1TrackingTree(0),
    fHorn2TrackingTree(0),
    fDPTrackingTree(0),
    fTargetOutputTree(0),
    fOutFile(0),
    fOutTree(0),
    nuNtuple(0),
    fOutFileDk2Nu(0),
    fOutTreeDk2Nu(0),
    fOutTreeDk2NuMeta(0),
    fTrackingTree(0),
    fDk2Nu(0),
    fDkMeta(0),
    fLBNEOutNtpData(0),
    fTrackingPlaneData(0)
{
  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetVerboseLevel() > 0) 
  {
     std::cout << "LBNEAnalysis Constructor Called." << std::endl;
  }

#ifdef G4ANALYSIS_USE
#endif

  
  //g4data = new data_t();
  fLBNEOutNtpData = new LBNEDataNtp_t();
  fTrackingPlaneData = new LBNEDataNtp_t();  // Not sure what this does?


  fcount = 0;
  fentry = 0;
  
  //
  //need code map defined but initialized
  //or get a seg fault don't know why
  //
  /*code[-13]   = 10;
  code[13]    = 11;
  code[111]   = 23;
  code[211]   = 13;
  code[-211]  = 14;
  code[130]   = 12;
  code[321]   = 15;
  code[-321]  = 16;
  code[2112]  = 8;
  code[2212]  = 1;
  code[-2212] = 2;
  code[310]   = 19;
  code[3122]  = 17;
  code[3222]  = 21;
  code[3212]  = 22;
  code[3112]  = 20;*/
  
  setDetectorPositions();
  
}

LBNEAnalysis::~LBNEAnalysis

(

)

Definition at line 120 of file LBNEAnalysis.cc.

{ 
  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetVerboseLevel() > 0) {
      std::cout << "LBNEAnalysis Destructor Called." << std::endl;
   }

#ifdef G4ANALYSIS_USE
  // delete things
#endif
}

Member Function Documentation

void LBNEAnalysis::CalcLocationWeights	(	const bsim::DkMeta *	dkmeta,
		bsim::Dk2Nu *	dk2nu,
		bool	useRealisticNearDetectorVolume
	)

Definition at line 2325 of file LBNEAnalysis.cc.

{       

  size_t nloc = dkmeta->location.size();        
  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
    // skip calculation for random location ... should already be filled        
    const std::string rkey = "random decay";    
    if ( dkmeta->location[iloc].name == rkey ) {        
      if ( iloc != 0 ) {        
        std::cerr << "calcLocationWeights \"" << rkey << "\""   
                  << " isn't the 0-th entry" << std::endl;      
        assert(0);      
      } 
      if ( dk2nu->nuray.size() != 1 ) { 
        std::cerr << "calcLocationWeights \"" << rkey << "\""   
                  << " nuenergy[" << iloc << "] not filled" << std::endl;       
        assert(0);      
      } 
      continue; 
    }   
    double xdet = dkmeta->location[iloc].x;
    double ydet = dkmeta->location[iloc].y;
    double zdet = dkmeta->location[iloc].z;

    if(useRealisticNearDetectorVolume) {
      // randomize location within a volume
      // for now, assume a 2 m x 2m x 6 m near detector
      xdet += 200*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
      ydet += 200*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
      zdet += 600*(CLHEP::HepRandom::getTheEngine()->flat()-0.5);
    }

    TVector3 xyzDet(xdet,       
                    ydet,       
                    zdet);  // position to evaluate     
    double enu_xy = 0;  // give a default value 
    double wgt_xy = 0;  // give a default value 
    int status = bsim::calcEnuWgt(dk2nu,xyzDet,enu_xy,wgt_xy);  
    if ( status != 0 ) {        
      std::cerr << "bsim::calcEnuWgt returned " << status << " for " 
                << dkmeta->location[iloc].name << std::endl;
    }   
    // with the recalculated energy compute the momentum components     
    TVector3 xyzDk(dk2nu->decay.vx,dk2nu->decay.vy,dk2nu->decay.vz);  // origin of decay
    TVector3 p3 = enu_xy * (xyzDet - xyzDk).Unit();
    bsim::NuRay anuray(p3.x(), p3.y(), p3.z(), enu_xy, wgt_xy);
    dk2nu->nuray.push_back(anuray);
  }
}       

void LBNEAnalysis::CloseDk2NuOutput

(

)

Definition at line 591 of file LBNEAnalysis.cc.

{  
  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
  if (pRunManager->GetCreateDk2NuOutput()) {
        fOutFileDk2Nu->cd();
        std::cerr << " Writing the Dk2Nu Ntuple ... " << std::endl; 
        fOutTreeDk2Nu->Write();
        fOutTreeDk2NuMeta->Write();
        std::cerr << " Written the Dk2Nu Ntuple and meta data.. ... " << std::endl; 
        
        fOutFileDk2Nu->Close();
        delete fOutFileDk2Nu;
        //}
  }

}

void LBNEAnalysis::CloseOutput

(

)

Definition at line 522 of file LBNEAnalysis.cc.

{  
  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
  if (pRunManager->GetCreateOutput()) {
     //if(fOutFile && fOutFile->IsOpen())
     //{
        fOutFile->cd();
        fOutTree->Write();

        if (pRunManager->GetCreateTrkPlaneOutput()){
        //---tracking plane data 
           std::cout << "   Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
           fTrackingTree->Write();
        }


          
        if(pRunManager->GetCreateTrkPlaneH1Output()){
        //---tracking plane data 
           std::cout << "   Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
           if (fHorn1TrackingTree != 0) fHorn1TrackingTree->Write();
        }

        if (pRunManager->GetCreateTrkPlaneH2Output()){
        //---tracking plane data 
           std::cout << "   Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
          if (fHorn2TrackingTree != 0) fHorn2TrackingTree->Write();
        }
        if
          (pRunManager->GetCreateTrkPlaneDPOutput()){
        //---tracking plane data 
           std::cout << "   Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
           if (fDPTrackingTree != 0) fDPTrackingTree->Write();
        }

        if (pRunManager->GetCreateTargetOutput()){
        //---tracking plane data 
           std::cout << "   Writing out Tracking Plane Ntuple to " << fOutFile -> GetName() << std::endl;
           if (fTargetOutputTree != 0) fTargetOutputTree->Write();
        }

        if (pRunManager->GetCreateAlcoveTrackingOutput())
        {
           std::cout << "  Writeing out sculpted absorber tracking plane to "<<fOutFile->GetName()<<std::endl;
           if (fAlcoveTrackingTree != 0 ) fAlcoveTrackingTree->Write();

        }


        fOutFile->Close();


/*
        if(fOutFile->IsOpen())
        {
           std::cout << "   PROBLEM: Failed to close Output Ntuple " << fOutFile -> GetName() << std::endl;
        }
        else 
        {*/
           std::cout << "   Sucessfully closed Output Ntuple " << fOutFile -> GetName() << std::endl;
/*      }
*/
        delete fOutFile;
        //}
  }

}

G4bool LBNEAnalysis::CreateDk2NuOutput

(

)

Definition at line 334 of file LBNEAnalysis.cc.

{

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

   if (pRunManager->GetCreateDk2NuOutput()) 
   {
      G4String spaces = "      ";
      std::cout << "    LBNEAnalysis::CreateOutput() - Creating output ntuple..." << std::endl; 

      
      {
         //std::cout << spaces << "Creating G4LBNEData Ntuple..." << std::endl;
         //return LBNEAnalysis::CreateG4LBNEDataNtp();

         pRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
         fDk2Nu  = new bsim::Dk2Nu;
         fDkMeta = new bsim::DkMeta;
         snprintf(nuNtupleFileName, 1023, "%s_%05d.dk2nu.root",
           (pRunManager->GetOutputDk2NuFileName()).c_str(),pRunManager->GetCurrentRun()->GetRunID());
         fOutFileDk2Nu = new TFile(nuNtupleFileName,"RECREATE","root ntuple");
         std::cerr << "Creating neutrino ntuple, Dk2Nu : " << nuNtupleFileName << std::endl;
         
         fOutTreeDk2Nu = new TTree("dk2nuTree", "g4lbne neutrino ntuple, dk2Nu format");
         fOutTreeDk2Nu->Branch("dk2nu", "bsim::Dk2Nu", &fDk2Nu, 32000, 99);
         fOutTreeDk2NuMeta  = new TTree("dkmetaTree", "g4lbne neutrino ntuple, dk2Nu format, metadata");
         fOutTreeDk2NuMeta->Branch("dkmeta", "bsim::DkMeta", &fDkMeta, 32000,99);
         std::cerr << "Done Creating neutrino ntuple, Dk2Nu : " << nuNtupleFileName << std::endl;

         return true;
         
      }
      
      return false;
      
   }
   
   return false;
   
   
}

G4bool LBNEAnalysis::CreateOutput

(

)

Definition at line 140 of file LBNEAnalysis.cc.

{

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

   if (pRunManager->GetCreateOutput()) 
   {
      G4String spaces = "      ";
      std::cout << "    LBNEAnalysis::CreateOutput() - Creating output ntuple..." << std::endl; 

      
      {
         //std::cout << spaces << "Creating G4LBNEData Ntuple..." << std::endl;
         //return LBNEAnalysis::CreateG4LBNEDataNtp();

         snprintf(nuNtupleFileName, 1023, "%s_%03d.root",
           (pRunManager->GetOutputNtpFileName()).c_str(),pRunManager->GetCurrentRun()->GetRunID());   //April 29, 2015->Fixed the dk2nu extension (line 176)
         fOutFile = new TFile(nuNtupleFileName,"RECREATE","root ntuple");                             //which was giving the error due to conflict in name.  
         std::cerr << "Creating neutrino ntuple: " << nuNtupleFileName << std::endl;
         //fOutTree = new TTree("nudata","g4numi Neutrino ntuple");
         //fOutTree->Branch("data","data_t",&g4data,32000,1);
         fOutTree = new TTree("nudata","g4lbne Neutrino ntuple");
         fOutTree->Branch("data","LBNEDataNtp_t",&fLBNEOutNtpData,32000,1);

         if (pRunManager->GetCreateTrkPlaneOutput()) {
           // create tree for tracking planes
           std::cout << spaces << "Creating Tracking Plane Ntuple..." << std::endl;
           fTrackingTree = new TTree("trackingdata", "Particles in Tracking Plane");
           fTrackingTree->Branch("num", &fNParticles, "fNParticles/I");
           fTrackingTree->Branch("trackID", &fTrackID, "fTrackID[fNParticles]/I]");
           fTrackingTree->Branch("pdgCode", &fParticlePDG, "fParticlePDG[fNParticles]/I");
           fTrackingTree->Branch("ImpWeight", &fNImpWt, "fNImpWt[fNParticles]/D");
           fTrackingTree->Branch("positionX", &fParticleX, "fParticleX[fNParticles]/D");
           fTrackingTree->Branch("positionY", &fParticleY, "fParticleY[fNParticles]/D");
           fTrackingTree->Branch("positionZ", &fParticleZ, "fParticleZ[fNParticles]/D");
           fTrackingTree->Branch("momentumX", &fParticleDX, "fParticleDX[fNParticles]/D");
           fTrackingTree->Branch("momentumY", &fParticleDY, "fParticleDY[fNParticles]/D");
           fTrackingTree->Branch("momentumZ", &fParticleDZ, "fParticleDZ[fNParticles]/D");
           fTrackingTree->Branch("mass", &fParticleMass, "fParticleMass[fNParticles]/D");
           fTrackingTree->Branch("energy", &fParticleEnergy, "fParticleEnergy[fNParticles]/D");
         }
 if (pRunManager->GetCreateTrkPlaneH1Output()) {
           // create tree for Horn 1 tracking planes Amit Bashyal   
           std::cout << spaces << "Creating Horn1 Tracking Plane Ntuple..." << std::endl;
           fHorn1TrackingTree = new TTree("H1trackingdata", "Particles in Horn1 Tracking Plane");
           fHorn1TrackingTree->Branch("nH1P", &fH1NParticles, "fH1NParticles/I");
           fHorn1TrackingTree->Branch("trackID", &fH1TrackID, "fH1TrackID/I]");
           fHorn1TrackingTree->Branch("pdgCode", &fH1ParticlePDG, "fH1ParticlePDG/I");
           fHorn1TrackingTree->Branch("H1ImpWeight", &fH1NImpWt, "fH1NImpWt/D");
           fHorn1TrackingTree->Branch("H1positionX", &fH1ParticleX, "fH1ParticleX/D");
           fHorn1TrackingTree->Branch("H1positionY", &fH1ParticleY, "fH1ParticleY/D");
           fHorn1TrackingTree->Branch("H1positionZ", &fH1ParticleZ, "fH1ParticleZ/D");
           fHorn1TrackingTree->Branch("H1momentumdirX", &fH1ParticleDX, "fH1ParticleDX/D");
           fHorn1TrackingTree->Branch("H1momentumdirY", &fH1ParticleDY, "fH1ParticleDY/D");
           fHorn1TrackingTree->Branch("H1momentumdirZ", &fH1ParticleDZ, "fH1ParticleDZ/D");
           fHorn1TrackingTree->Branch("mass", &fH1ParticleMass, "fH1ParticleMass/D");
           fHorn1TrackingTree->Branch("energy", &fH1ParticleEnergy, "fH1ParticleEnergy/D");
           fHorn1TrackingTree->Branch("H1momentumX",&fH1ParticlePX,"fH1ParticlePX/D");
           fHorn1TrackingTree->Branch("H1momentumY",&fH1ParticlePY,"fH1ParticlePY/D");
           fHorn1TrackingTree->Branch("H1momentumZ",&fH1ParticlePZ,"fH1ParticlePZ/D");
           fHorn1TrackingTree->Branch("H1PXPZ",&fH1ParticlePXPZ,"fH1ParticlePXPZ/D");
           fHorn1TrackingTree->Branch("H1PYPZ",&fH1ParticlePYPZ,"fH1ParticlePYPZ/D");
           fHorn1TrackingTree->Branch("H1PProductionX",&fH1PProductionX,"fH1PProductionX/D");
           fHorn1TrackingTree->Branch("H1PProductionY",&fH1PProductionY,"fH1PProductionY/D");
           fHorn1TrackingTree->Branch("H1PProductionZ",&fH1PProductionZ,"fH1PProductionZ/D");   
           fHorn1TrackingTree->Branch("H1PmomentumdirX",&fH1PProductionDX,"fH1PProductionDX/D");
           fHorn1TrackingTree->Branch("H1PmomentmdirY",&fH1PProductionDY,"fH1PProductionDY/D");
           fHorn1TrackingTree->Branch("H1PmomentumdirZ",&fH1PProductionDZ,"fH1PProductionDZ/D");
           }
                
                
          if (pRunManager->GetCreateTrkPlaneDPOutput()) {
          // Create Planes for Decay Pipe End Tracking Plane Amit Bashyal
           fDPTrackingTree = new TTree("DecayPipetrackingdata", "Particles Making to the Decay Pipe");
            std::cerr << spaces << "Creating DecaypipeOutput Ntuples"<<std::endl;
           fDPTrackingTree->Branch("DPImpWeight", &fDPNImpWt, "fDPNImpWt/D");
           fDPTrackingTree->Branch("DPPProductionX",&fDPPProductionX,"fDPPProductionX/D");
           fDPTrackingTree->Branch("DPPProductionY",&fDPPProductionY,"fDPPProductionY/D");
           fDPTrackingTree->Branch("DPPProductionZ",&fDPPProductionZ,"fDPPProductionZ/D");
           fDPTrackingTree->Branch("DPPtrackID", &fDPTrackID, "fDPTrackID/I]");
           fDPTrackingTree->Branch("DPParentID", &fDPParentID, "fDPParentID/I]");
           fDPTrackingTree->Branch("DPenergy", &fDPParticleEnergy, "fDPParticleEnergy/D");
           fDPTrackingTree->Branch("DPpositionX", &fDPParticleX, "fDPParticleX/D");
           fDPTrackingTree->Branch("DPpositionY", &fDPParticleY, "fDPParticleY/D");
           fDPTrackingTree->Branch("DPpositionZ", &fDPParticleZ, "fDPParticleZ/D");
           fDPTrackingTree->Branch("DPmomentumX",&fDPParticlePX,"fDPParticlePX/D");
           fDPTrackingTree->Branch("DPmomentumY",&fDPParticlePY,"fDPParticlePY/D");
           fDPTrackingTree->Branch("DPmomentumZ",&fDPParticlePZ,"fDPParticlePZ/D");
           fDPTrackingTree->Branch("DPpdgCode", &fDPParticlePDG, "fDPParticlePDG/I");
           }
            if (pRunManager->GetCreateTargetOutput()) {
            std::cout<<spaces<<"Creating TargetOutPut Ntuples"<<std::endl;
           fTargetOutputTree = new TTree("TargetOutputData", "Particles Produced in the Target");
           fTargetOutputTree->Branch("TtrackID", &fTTrackID, "fTTrackID/I]");
           fTargetOutputTree->Branch("TpdgCode", &fTParticlePDG, "fTParticlePDG/I");
           fTargetOutputTree->Branch("TpositionX", &fTParticleX, "fTParticleX/D");
           fTargetOutputTree->Branch("TpositionY", &fTParticleY, "fTParticleY/D");
           fTargetOutputTree->Branch("TpositionZ", &fTParticleZ, "fTParticleZ/D");
           fTargetOutputTree->Branch("Tenergy", &fTParticleEnergy, "fTParticleEnergy/D");
           fTargetOutputTree->Branch("TmomentumX",&fTParticlePX,"fTParticlePX/D");
           fTargetOutputTree->Branch("TmomentumY",&fTParticlePY,"fTParticlePY/D");
           fTargetOutputTree->Branch("TmomentumZ",&fTParticlePZ,"fTParticlePZ/D");   

         }
         
     //    fOutTree->Branch("horn2data","LBNEDataNtp_t",&fLBNEOutNtpData,32000,1);
         if  (pRunManager->GetCreateTrkPlaneH2Output()) {
           // create tree for Horn 2 tracking planes
           std::cout << spaces << "Creating Horn2 Tracking Plane Ntuple..." << std::endl;
           fHorn2TrackingTree = new TTree("H2trackingdata", "Particles in Horn2 Tracking Plane");
           fHorn2TrackingTree->Branch("nH2P", &fH2NParticles, "fH2NParticles/I");
           fHorn2TrackingTree->Branch("trackID", &fH2TrackID, "fH2TrackID/I]");
           fHorn2TrackingTree->Branch("pdgCode", &fH2ParticlePDG, "fH2ParticlePDG/I");
           fHorn2TrackingTree->Branch("H2ImpWeight", &fH2NImpWt, "fH2NImpWt/D");
           fHorn2TrackingTree->Branch("H2positionX", &fH2ParticleX, "fH2ParticleX/D");
           fHorn2TrackingTree->Branch("H2positionY", &fH2ParticleY, "fH2ParticleY/D");
           fHorn2TrackingTree->Branch("H2positionZ", &fH2ParticleZ, "fH2ParticleZ/D");
           fHorn2TrackingTree->Branch("H2momentumdirX", &fH2ParticleDX, "fH2ParticleDX/D");
           fHorn2TrackingTree->Branch("H2momentumdirY", &fH2ParticleDY, "fH2ParticleDY/D");
           fHorn2TrackingTree->Branch("H2momentumdirZ", &fH2ParticleDZ, "fH2ParticleDZ/D");
           fHorn2TrackingTree->Branch("mass", &fH2ParticleMass, "fH2ParticleMass/D");
           fHorn2TrackingTree->Branch("energy", &fH2ParticleEnergy, "fH2ParticleEnergy/D");
           fHorn2TrackingTree->Branch("H2momentumX",&fH2ParticlePX,"fH2ParticlePX/D");
           fHorn2TrackingTree->Branch("H2momentumY",&fH2ParticlePY,"fH2ParticlePY/D");
           fHorn2TrackingTree->Branch("H2momentumZ",&fH2ParticlePZ,"fH2ParticlePZ/D");
           fHorn2TrackingTree->Branch("H2PXPZ",&fH2ParticlePXPZ,"fH2ParticlePXPZ/D");
           fHorn2TrackingTree->Branch("H2PYPZ",&fH2ParticlePYPZ,"fH2ParticlePYPZ/D");
           fHorn2TrackingTree->Branch("H2PProductionX",&fH2PProductionX,"fH2PProductionX/D");
           fHorn2TrackingTree->Branch("H2PProductionY",&fH2PProductionY,"fH2PProductionY/D");
           fHorn2TrackingTree->Branch("H2PProductionZ",&fH2PProductionZ,"fH2PProductionZ/D");
           fHorn2TrackingTree->Branch("H2PmomentumdirX",&fH2PProductionDX,"fH2PProductionDX/D");
           fHorn2TrackingTree->Branch("H2PmomentumdirY",&fH2PProductionDY,"fH2PProductionDY/D");
           fHorn2TrackingTree->Branch("H2PmomentumdirZ",&fH2PProductionDZ,"fH2PProductionDZ/D");
           fHorn2TrackingTree->Branch("H2CNProcess",&fH2CNProcess,"fH2CNProcess/I");

         }

         if (pRunManager->GetCreateAlcoveTrackingOutput()){
           std::cout << spaces << "Creating Sculpted Absorber Alcove Tracking Plane Ntuple..." << std::endl;
           fAlcoveTrackingTree = new TTree("AlcoveTracks","Alcove Tracking");
            
           fAlcoveTrackingTree->Branch("run",&fMuRunNo,"run/I");
           fAlcoveTrackingTree->Branch("event",&fMuEvtNo,"event/I");
           fAlcoveTrackingTree->Branch("np",&fMuNParticles,"np/I");
           fAlcoveTrackingTree->Branch("ID",&fMuTrackID,"ID[np]/I");
           fAlcoveTrackingTree->Branch("ParID",&fMuParentID,"ParID[np]/I");
           fAlcoveTrackingTree->Branch("PDG",&fMuPDG,"PDG[np]/I");
           fAlcoveTrackingTree->Branch("impwt",&fMuNimpWt,"impwt[np]/D");
           fAlcoveTrackingTree->Branch("x",&fMuX,"x[np]/D");
           fAlcoveTrackingTree->Branch("y",&fMuY,"y[np]/D");
           fAlcoveTrackingTree->Branch("z",&fMuZ,"z[np]/D");
           fAlcoveTrackingTree->Branch("startx",&fMuStartX,"startx[np]/D");
           fAlcoveTrackingTree->Branch("starty",&fMuStartY,"starty[np]/D");
           fAlcoveTrackingTree->Branch("startz",&fMuStartZ,"startz[np]/D");
           fAlcoveTrackingTree->Branch("startE",&fMuStartE,"startE[np]/D");
           fAlcoveTrackingTree->Branch("px",&fMuPX,"px[np]/D");
           fAlcoveTrackingTree->Branch("py",&fMuPY,"py[np]/D");
           fAlcoveTrackingTree->Branch("pz",&fMuPZ,"pz[np]/D");
           fAlcoveTrackingTree->Branch("m",&fMuMass,"m[np]/D");
           fAlcoveTrackingTree->Branch("E",&fMuEnergy,"E[np]/D");
           fAlcoveTrackingTree->Branch("dEdx",&fMudEdx,"dEdx[np]/D");
           fAlcoveTrackingTree->Branch("dEdx_ion",&fMudEdx_ion,"dEdx_ion[np]/D");
           fAlcoveTrackingTree->Branch("cosTheta",&fMuTheta,"cosTheta[np]/D");
           fAlcoveTrackingTree->Branch("edep",&fMuDE,"edep[np]/D");
           fAlcoveTrackingTree->Branch("edep_ion",&fMuDEion,"edep_ion[np]/D");
           fAlcoveTrackingTree->Branch("nsteps",&fMuNSteps,"nsteps[np]/I");
           fAlcoveTrackingTree->Branch("parPDG",&fMuParPDG,"parPDG[np]/I");
           fAlcoveTrackingTree->Branch("parE",&fMuParE,"parE[np]/D");
           fAlcoveTrackingTree->Branch("parX",&fMuParX,"parX[np]/D");
           fAlcoveTrackingTree->Branch("parY",&fMuParY,"parY[np]/D");
           fAlcoveTrackingTree->Branch("parZ",&fMuParZ,"parZ[np]/D");
           fAlcoveTrackingTree->Branch("parPX",&fMuParPX,"parPX[np]/D");
           fAlcoveTrackingTree->Branch("parPY",&fMuParPY,"parPY[np]/D");
           fAlcoveTrackingTree->Branch("parPZ",&fMuParPZ,"parPZ[np]/D");


           fMuNParticles = 0;//ensure that we start with no particles
         }//Sculpted absorber tracking plane

         return true;
         
      }
      
      return false;
      
   }
   
   return false;
   
   
}

void LBNEAnalysis::FillAlcoveTrackingPlaneData

(

const G4Step &

aStep

)

Definition at line 2130 of file LBNEAnalysis.cc.

{

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
  if (pRunManager->GetCreateAlcoveTrackingOutput())
  {
    G4Track *track = aStep.GetTrack();
    LBNETrackInformation *info =
            (LBNETrackInformation*)(track->GetUserInformation());
    G4StepPoint *point = aStep.GetPostStepPoint();
    G4ThreeVector mom = point->GetMomentum();
    G4ThreeVector pos = point->GetPosition();

    int theTrackID = aStep.GetTrack()->GetTrackID();
    int index = -1;


  //First look for this track
    for (int i = 0 ; i < fMuNParticles; i++)
    {
      if (fMuTrackID[i] == theTrackID){
        index = i;
      }
    }

    if (index > -1){//Found track

      fMuNSteps[index]++;
      fMuDE[index] += aStep.GetTotalEnergyDeposit()/CLHEP::MeV;
      fMuDEion[index] += (aStep.GetTotalEnergyDeposit()- aStep.GetNonIonizingEnergyDeposit())/CLHEP::MeV;
      fMudEdx[index] += aStep.GetStepLength()/CLHEP::cm;
      fMudEdx_ion[index] += aStep.GetStepLength()/CLHEP::cm;
    }else{//Track not already saved

      index = fMuNParticles;
      if (index >= kMaxP)
      {
        G4cout << "Alcove Tracking: Max Number of Particles Reached. Not saving ... " <<G4endl;
        return;
      }
      fMuNParticles++;
      fMuEvtNo = pRunManager->GetCurrentEvent()->GetEventID();
      fMuRunNo = pRunManager->GetCurrentRun()->GetRunID();

      G4ThreeVector vpos = track->GetVertexPosition();
      if(info != 0){
        fMuNimpWt[index] = info->GetNImpWt();
      } 
      
      G4ParticleDefinition* theParticle = track->GetDefinition();

      fMuTrackID[index] = theTrackID;
      fMuParentID[index] = track->GetParentID();
      
      fMuPDG[index] = theParticle->GetPDGEncoding();
      fMuMass[index] = theParticle->GetPDGMass()/CLHEP::GeV;

      fMuT[index] = point->GetGlobalTime()/CLHEP::ns;
      fMuPX[index] = mom.x()/CLHEP::GeV;
      fMuPY[index] = mom.y()/CLHEP::GeV;
      fMuPZ[index] = mom.z()/CLHEP::GeV;
      fMuX[index] = pos.x()/CLHEP::cm;
      fMuY[index] = pos.y()/CLHEP::cm;
      fMuZ[index] = pos.z()/CLHEP::cm;
      fMuTheta[index] = mom.z()/sqrt(mom.z()*mom.z()+mom.y()*mom.y()+mom.x()*mom.x());
      fMuEnergy[index] = point->GetKineticEnergy()/CLHEP::GeV;
      fMuStartX[index] = vpos.x()/CLHEP::cm;
      fMuStartY[index] = vpos.y()/CLHEP::cm;
      fMuStartZ[index] = vpos.z()/CLHEP::cm;
      fMuStartE[index] = track->GetVertexKineticEnergy()/CLHEP::GeV;

      fMuDE[index] = aStep.GetTotalEnergyDeposit()/CLHEP::MeV;
      fMuDEion[index] = fMuDE[index] - aStep.GetNonIonizingEnergyDeposit()/CLHEP::MeV;
      fMudEdx[index] = aStep.GetStepLength()/CLHEP::cm;
      fMudEdx_ion[index] = aStep.GetStepLength()/CLHEP::cm;
      fMuNSteps[index] = 1;     
      //Parent info 
      LBNETrajectory* parTraj = GetParentTrajectory(fMuParentID[index]);

      if (parTraj){
        int npt = parTraj->GetPointEntries();
        fMuParPDG[index] = parTraj->GetPDGEncoding();
        fMuParX[index] = parTraj->GetPoint(npt-1)->GetPosition().x()/CLHEP::cm;
        fMuParY[index] = parTraj->GetPoint(npt-1)->GetPosition().y()/CLHEP::cm;
        fMuParZ[index] = parTraj->GetPoint(npt-1)->GetPosition().z()/CLHEP::cm;
        fMuParPX[index] = parTraj->GetMomentum(npt-1).x()/CLHEP::GeV;
        fMuParPY[index] = parTraj->GetMomentum(npt-1).y()/CLHEP::GeV;
        fMuParPZ[index] = parTraj->GetMomentum(npt-1).z()/CLHEP::GeV;
        fMuParE[index] = sqrt(fMuParPX[index]*fMuParPX[index]+fMuParPY[index]*fMuParPY[index]+fMuParZ[index]*fMuParZ[index]  + parTraj->GetMass()*parTraj->GetMass() / (CLHEP::GeV*CLHEP::GeV)  )-parTraj->GetMass()/ CLHEP::GeV; 
      }else{
        fMuParPDG[index] = 0;
        fMuParE[index] = -1;
      }
    }
  }//end if sculpted muon alcove plane is on

}

void LBNEAnalysis::fillDkMeta

(

)

Definition at line 377 of file LBNEAnalysis.cc.

                              {

   VersionAndContext* versionContext = VersionAndContext::Instance();

   LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
   if (!theRunManager->GetCreateDk2NuOutput()) return;
   
   fDkMeta->job = theRunManager->GetCurrentRun()->GetRunID();
   fDkMeta->pots = theRunManager->GetNumberOfEventsLBNE();
   fDkMeta->beamsim = std::string("g4lbnf-");
   fDkMeta->beamsim += versionContext->GetGitDescription();
   fDkMeta->beamsim += std::string (", builddir=") + 
                        versionContext->GetG4lbnfDir();
   fDkMeta->beamsim += std::string (", on") + 
                        versionContext->GetBuildTimeStamp();   
   fDkMeta->physics = G4Version;
   fDkMeta->physcuts = theRunManager->GetPhysicsListName();
   fDkMeta->tgtcfg = "le00zdu";
   
   {
     std::ostringstream oStrStr; 
     oStrStr << " --Kine Cut " << theRunManager->GetLBNESteppingManager()->GetKillTrackingThreshold()/CLHEP::GeV <<"CLHEP::GeV";
     fDkMeta->physcuts += oStrStr.str();
   }
   
   LBNEVolumePlacements *plManager = LBNEVolumePlacements::Instance();
   {
     //if (plManager->GetUse1p2MW()) fDkMeta->tgtcfg += std::string("1p2MW");
     //else fDkMeta->tgtcfg += std::string("700kW");
     std::ostringstream oStrStr;
     oStrStr.precision(4);
     oStrStr << "-length" << (plManager->GetTargetSLengthGraphite())/CLHEP::cm;
     oStrStr << "-Pos" << (plManager->GetTargetSLengthGraphite() - 
                           plManager->GetTargetLengthIntoHorn() + 25.3*CLHEP::mm)/CLHEP::cm;
     oStrStr << "cm";
     fDkMeta->tgtcfg = "le00zdu";//versionContext->GetMacroFileName();
   }
   {
     std::ostringstream oStrStr;
     oStrStr.precision(4);
     const LBNEDetectorConstruction *pDet = 
            static_cast<const LBNEDetectorConstruction*> (theRunManager->GetUserDetectorConstruction());
     oStrStr << "Current-" << pDet->GetHornCurrent()/(CLHEP::ampere*1000.); //in kAmp
    // fDkMeta->horncfg = oStrStr.str();
     oStrStr << "Current-" << pDet->GetHornCurrent()/(CLHEP::ampere*1000.); //in kAmp    
   // fDkMeta->horncfg = std::string("230i");
   fDkMeta->horncfg = pDet->GetHornCurrent()/(CLHEP::ampere*1000.);
   }
  // fDkMeta->dkvolcfg = plManager->GetDecayPipeGas();
  std::cout<<"GetDecayPipeGas() "<<plManager->GetDecayPipeGas()<<std::endl;
  fDkMeta->dkvolcfg = std::string("helium");
   {
     std::ostringstream oStrStr;
     oStrStr.precision(4);
     oStrStr << "length-" << plManager->GetDecayPipeLength()/CLHEP::cm << "cm";
    // fDkMeta->dkvolcfg += oStrStr.str();
   }
   //
   // Beam parameters 
   //
   const LBNEPrimaryGeneratorAction* NPGA=
      static_cast<const LBNEPrimaryGeneratorAction*> (theRunManager->GetUserPrimaryGeneratorAction());
   fDkMeta->beam0x = NPGA->GetBeamOffsetX();
   fDkMeta->beam0y = NPGA->GetBeamOffsetY();
   fDkMeta->beam0z = NPGA->GetBeamOffsetZ();
   fDkMeta->beamhwidth = NPGA->GetBeamSigmaX();
   fDkMeta->beamvwidth = NPGA->GetBeamSigmaY();
   fDkMeta->beamdxdz = 0.;
   fDkMeta->beamdydz = NPGA->GetBeamAngleTheta();
   
      //A. Bashyal stuff for the ppfx:
  G4String name_gen[] = {"parent","granparent","greatgranparent"};
  G4String name_vol[] = {"PHorn1IC","PHorn2IC","DPIP","DVOL"}; //First we try to verify if this will work on ppfx.
  for(G4int ii=0;ii<nGenAbs;ii++){
    GenAbsName.push_back(name_gen[ii]);
  }
  for(G4int ii=0;ii<nVolAbs;ii++){
    VolAbsName.push_back(name_vol[ii]);
  }
//  char namevar[50];
  for(int ii=0; ii<nVolAbs; ii++){
    for(int jj=0; jj<nGenAbs; jj++){
//      sprintf(namevar,"Material_%s_%s",VolAbsName[ii].c_str(),GenAbsName[jj].c_str());
//      VolVdblName.push_back(G4String(namevar));
//       recoding this in C++, potential overwrite with fix C-style arrays. 
         std::ostringstream nameVarStrStr; nameVarStrStr << "Material_" << VolAbsName[ii] << "_" << GenAbsName[jj];
         const std::string nameVarStr(nameVarStrStr.str());
         VolVdblName.push_back(nameVarStr); 
    }
  }
  nVdblTot = int(VolVdblName.size());
  
  //For vint branch
  nVintTot = 2;
  VolVintName.clear(); 
  G4String name_vint[] = {"Index_Tar_In_Ancestry","playlistID"};
  for(G4int ii=0;ii<nVintTot;ii++){
    VolVintName.push_back(name_vint[ii]);
  }
  fDkMeta->vintnames.clear();
  for(int ii=0;ii<nVintTot;ii++){
    (fDkMeta->vintnames).push_back(VolVintName[ii]);
   std::cout<< " LBNEAnalysis:: nVintnames are : "<<VolVintName[ii]<<std::endl;
  }
  fDkMeta->vdblnames.clear();
  for(int ii=0;ii<nVdblTot;ii++){
    (fDkMeta->vdblnames).push_back(VolVdblName[ii]);
    std::cout<<" LBNEAnalysis: nVdblTot are : "<<VolVdblName[ii]<<std::endl;
  }
   
   // Detector positions  
   if (!fDk2NuDetectorFileRead) {
     fDkMeta->location.clear();
     fDk2NuDetectorFileRead=true;
     std::string aFileName(theRunManager->GetDetectorLocationFileName());
     if (aFileName == G4String("?")) {
       char *G4WORKDIRC = getenv("G4LBNE_DIR");
       if (G4WORKDIRC == 0) {
         G4WORKDIRC = getenv("G4LBNE_DIR");
         if (G4WORKDIRC == 0) {
          G4Exception("LBNEAnalysis::LBNEAnalysis::fillDkMeta", "InvalidSetup", 
                    FatalErrorInArgument, 
                    "Environmental variable G4LBNE_DIR or G4LBNE not set, can't find default Det. loc. file ");
         }
       }
       aFileName = std::string(G4WORKDIRC) + std::string("/locations/etc/LBNElocations.txt");
       std::cerr << "  ... And this file name is  ... " << aFileName << std::endl;
     }
     // Is this open? 
     std::ifstream fStr;
     fStr.open(aFileName.c_str());
     if (!fStr.is_open()) {
       std::cerr << " LBNEAnalysis::LBNEAnalysis::fillDkMeta, with DK2NU set to " << std::string(getenv("G4LBNE_DIR")) << std::endl;
       std::cerr << "  can't open detector location file ... " << aFileName << std::endl;
       G4Exception("LBNEAnalysis::LBNEAnalysis::fillDkMeta", "InvalidSetup", 
                    FatalErrorInArgument, "Invalid Detector FileName (not found or corrupted)");
     } else {
       fStr.close();
     }
     bsim::readWeightLocations(aFileName, fDkMeta);
   }
   fOutTreeDk2NuMeta->Fill();
}

void LBNEAnalysis::FillEvent

(

)

Definition at line 2304 of file LBNEAnalysis.cc.

{
  if(fNParticles>0){
    fTrackingTree->Fill();
}

 if(fTNParticles>0){
    fTargetOutputTree->Fill();
}


  for (int i = 0 ; i < fMuNParticles; i++)
  {
    fMudEdx[i] = fMuDE[i]/fMudEdx[i]; 
    fMudEdx_ion[i] = fMuDEion[i]/fMudEdx_ion[i]; 
  }
  fAlcoveTrackingTree->Fill();


}

void LBNEAnalysis::FillNeutrinoNtuple	(	const G4Track &	track,
		const std::vector< G4VTrajectory * > &	nuHistory
	)

Definition at line 634 of file LBNEAnalysis.cc.

{

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
   const LBNERunAction *runUserAction =  reinterpret_cast<const LBNERunAction *>(pRunManager->GetUserRunAction());
   
   if (runUserAction->GetDoComputeEDepInHorns()) return;
 
 if (pRunManager->GetVerboseLevel() > 3) 
   { G4cout << "LBNEAnalysis::FillNeutrinoNtuple() called." << G4endl;}
   
  if ((!pRunManager->GetCreateOutput()) && (!pRunManager->GetCreateDk2NuOutput()))  return;

   if (pRunManager->GetCreateOutput()) fLBNEOutNtpData -> Clear();
   

   LBNEPrimaryGeneratorAction *NPGA = (LBNEPrimaryGeneratorAction*)(pRunManager)->GetUserPrimaryGeneratorAction();
   G4int evtno                      = pRunManager->GetCurrentEvent()->GetEventID();
    
  
   //Neutrino vertex position and momentum
   G4ThreeVector pos = track.GetPosition()/CLHEP::mm; 
   bool validTrack = true;
   if (std::isnan(pos.x()) || std::isnan(pos.y()) || std::isnan(pos.z())) {
     std::cerr << " LBNEAnalysis::FillNeutrinoNtuple  Bogus track position for track ID " 
               << track.GetTrackID() << " parentID " <<  track.GetParentID() << std::endl;
        validTrack = false;
        return;       
   }
   x = pos.x();
   y = pos.y();
   z = pos.z();
   NuMomentum = track.GetMomentum();
   G4int parentID           = track.GetParentID();

   if (pRunManager->GetCreateOutput()) fLBNEOutNtpData->ptrkid = parentID;
   
   if (parentID == 0)
   { std::cout << "Event no: " << evtno << " LBNEAnalysis::FillNeutrinoNtuple - "
               << "PROBLEM: Direct Nu Parent has track ID = 0." << std::endl;
      return; //I have to make some changes so that neutrinos in fluka/mars ntuples don't crash
   }
   
   // Debugging anomalous abosrption here... 
   if(LBNEQuickPiToNuVect::Instance()->doIt()) {
      LBNEQuickPiToNuVect::Instance()->blessPion(parentID, track.GetTotalEnergy(), track.GetMomentum());
   }

   LBNETrajectory* NuParentTrack     = GetParentTrajectory(parentID);
   G4int point_no                    = NuParentTrack->GetPointEntries();
   G4ThreeVector ParentMomentumFinal = NuParentTrack->GetMomentum(point_no-1);
   G4ThreeVector vertex_r            = NuParentTrack->GetPoint(point_no-1)->GetPosition(); //Should be the same as Neutrino vertex
   G4String parent_name              = NuParentTrack->GetParticleName();
   G4double Parent_mass              = NuParentTrack->GetMass();
   G4double gamma                    = sqrt(ParentMomentumFinal*ParentMomentumFinal+Parent_mass*Parent_mass)/Parent_mass; 
   G4double Parent_energy            = gamma*Parent_mass;
   G4ThreeVector beta_vec            = ParentMomentumFinal/Parent_energy;
   G4double partial                  = gamma*(beta_vec*NuMomentum);
   G4double enuzr                    = gamma*(track.GetTotalEnergy())-partial; //neutrino energy in parent rest frame
   G4double enuzrInGeV               = enuzr/CLHEP::GeV; //neutrino energy in parent rest frame, in CLHEP::GeV 
   //fill histograms, ntuples,...
   if (!validTrack) {
   
     std::cerr << " More info about bogus track... Parent name " << parent_name << "  Mom " 
               << ParentMomentumFinal.x() << " / " << ParentMomentumFinal.y() << " / " << ParentMomentumFinal.z() 
               <<  std::endl;
     
   }
   if (pRunManager->GetCreateOutput()) {
     fLBNEOutNtpData->protonN    = NPGA->GetProtonNumber();
   
     fLBNEOutNtpData->run        = pRunManager->GetCurrentRun()->GetRunID();
     fLBNEOutNtpData->evtno      = pRunManager->GetCurrentEvent()->GetEventID();
     fLBNEOutNtpData->beamHWidth = NPGA->GetBeamSigmaX()/CLHEP::cm;
     fLBNEOutNtpData->beamVWidth = NPGA->GetBeamSigmaY()/CLHEP::cm;
     fLBNEOutNtpData->beamX      = NPGA->GetBeamOffsetX()/CLHEP::cm;
     fLBNEOutNtpData->beamY      = NPGA->GetBeamOffsetY()/CLHEP::cm;
   }
   //G4int particleID = track.GetParentID();
   G4ThreeVector protonOrigin = NPGA->GetProtonOrigin();
   if (pRunManager->GetCreateOutput()) {
     fLBNEOutNtpData->protonX            = protonOrigin[0];
     fLBNEOutNtpData->protonY            = protonOrigin[1];
     fLBNEOutNtpData->protonZ            = protonOrigin[2];
   }
   
   G4ThreeVector protonMomentum = NPGA->GetProtonMomentum();
   if (pRunManager->GetCreateOutput()) {
     fLBNEOutNtpData->protonPx             = protonMomentum[0];
     fLBNEOutNtpData->protonPy             = protonMomentum[1];
     fLBNEOutNtpData->protonPz             = protonMomentum[2];
   }
//   std::cerr << " protonOrigin " << protonOrigin << " momentum " << protonMomentum << std::endl;
   
   LBNEVolumePlacements *volDB = LBNEVolumePlacements::Instance();
   if (pRunManager->GetCreateOutput()) {
     fLBNEOutNtpData->nuTarZ      = volDB->GetTargetLengthIntoHorn(); // A better info that the somewhat meaningless Z0 
   }
   const LBNEDetectorConstruction *detDB = 
      dynamic_cast<const LBNEDetectorConstruction*>(pRunManager->GetUserDetectorConstruction());
   //other info
   // Neutrino origin:
   // 3 From muon decay
   // 1 From particle from target
   // 2 From scraping
   //check if nu is from muon decay or from a particle from target, otherwise Norig = 2
   G4int Norig = 2;
   if ((parent_name=="mu+") || (parent_name=="mu-")) Norig = 3;
   G4String firstvolname = NuParentTrack->GetPreStepVolumeName(0);
   if (firstvolname.contains("Baffle") || firstvolname.contains("TGT")) Norig = 1;
   
   if (pRunManager->GetCreateOutput()) {
      fLBNEOutNtpData->hornCurrent = detDB->GetHornCurrent()/CLHEP::ampere/1000.;
   
     // Random decay - these neutrinos rarely hit any of the detectors
     fLBNEOutNtpData->Ndxdz   = NuMomentum[0]/NuMomentum[2];
     fLBNEOutNtpData->Ndydz   = NuMomentum[1]/NuMomentum[2];
     fLBNEOutNtpData->Npz     = NuMomentum[2]/CLHEP::GeV;
     fLBNEOutNtpData->Nenergy = track.GetTotalEnergy()/CLHEP::GeV;
   
   
     fLBNEOutNtpData->Norig  = Norig;
     fLBNEOutNtpData->Ndecay = NuParentTrack->GetDecayCode();
   
     G4ParticleDefinition * particleType = track.GetDefinition();
     G4int ntype   = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(particleType->GetParticleName()));
     fLBNEOutNtpData->Ntype = ntype;
     fLBNEOutNtpData->Vx    = x/CLHEP::cm;
     fLBNEOutNtpData->Vy    = y/CLHEP::cm;
     fLBNEOutNtpData->Vz    = z/CLHEP::cm;
     fLBNEOutNtpData->pdPx  = ParentMomentumFinal[0]/CLHEP::GeV;
     fLBNEOutNtpData->pdPy  = ParentMomentumFinal[1]/CLHEP::GeV;
     fLBNEOutNtpData->pdPz  = ParentMomentumFinal[2]/CLHEP::GeV;
   
     G4ThreeVector ParentMomentumProduction = NuParentTrack->GetMomentum(0);
     fLBNEOutNtpData->ppdxdz = ParentMomentumProduction[0]/ParentMomentumProduction[2];
     fLBNEOutNtpData->ppdydz = ParentMomentumProduction[1]/ParentMomentumProduction[2];
     fLBNEOutNtpData->pppz   = ParentMomentumProduction[2]/CLHEP::GeV; 
   
     G4double parentp = sqrt(ParentMomentumProduction*ParentMomentumProduction);
   
     fLBNEOutNtpData->ppenergy = sqrt((parentp*parentp+Parent_mass*Parent_mass))/CLHEP::GeV;
   
     fLBNEOutNtpData->ppmedium = 0.; //this is still empty
   
     fLBNEOutNtpData->ptype = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(parent_name));
   
     G4ThreeVector production_vertex = NuParentTrack->GetPoint(0)->GetPosition(); 
     fLBNEOutNtpData->ppvx = production_vertex[0]/CLHEP::cm;
     fLBNEOutNtpData->ppvy = production_vertex[1]/CLHEP::cm;
     fLBNEOutNtpData->ppvz = production_vertex[2]/CLHEP::cm;
   
   //if nu parent is a muon then find muon parent info
     if ((parent_name=="mu+" || parent_name=="mu-") && NuParentTrack->GetParentID()!=0)
       {
         G4int mupar                   = NuParentTrack->GetParentID();
         LBNETrajectory* MuParentTrack = GetParentTrajectory(mupar);
         G4int nopoint_mupar           = MuParentTrack->GetPointEntries();
         G4ThreeVector muparp          = MuParentTrack->GetMomentum(nopoint_mupar-1);
         G4double muparm               = MuParentTrack->GetMass();
         fLBNEOutNtpData->muparpx               = muparp[0]/CLHEP::GeV; // vector of hadron parent of muon
         fLBNEOutNtpData->muparpy               = muparp[1]/CLHEP::GeV; // 
         fLBNEOutNtpData->muparpz               = muparp[2]/CLHEP::GeV;
         fLBNEOutNtpData->mupare                = (sqrt(muparp*muparp+muparm*muparm))/CLHEP::GeV;
      }
      else
      {
         fLBNEOutNtpData->muparpx = -999999.;  
         fLBNEOutNtpData->muparpy = -999999.;
         fLBNEOutNtpData->muparpz = -999999.;
         fLBNEOutNtpData->mupare = -999999.;
      }
      
     fLBNEOutNtpData->Necm               = enuzr/CLHEP::GeV; // Neutrino energy in parent rest frame
     LBNETrackInformation* info = (LBNETrackInformation*)(track.GetUserInformation());
     fLBNEOutNtpData->Nimpwt             = info->GetNImpWt();  // Importance weight
     fLBNEOutNtpData->tgen               = info->GetTgen()-1;
   
     fLBNEOutNtpData->xpoint = 0.;  // x, y, z of parent at user selected vol
     fLBNEOutNtpData->xpoint = 0.;
     fLBNEOutNtpData->xpoint = 0.;
   
     {
       G4ThreeVector ParticlePosition = NPGA->GetParticlePosition();
       fLBNEOutNtpData->tvx = ParticlePosition[0]/CLHEP::cm;
       fLBNEOutNtpData->tvy = ParticlePosition[1]/CLHEP::cm;
       fLBNEOutNtpData->tvz = ParticlePosition[2]/CLHEP::cm;
       G4ThreeVector ParticleMomentum = NPGA->GetParticleMomentum();
       fLBNEOutNtpData->tpx = ParticleMomentum[0]/CLHEP::GeV;
       fLBNEOutNtpData->tpy = ParticleMomentum[1]/CLHEP::GeV;
       fLBNEOutNtpData->tpz = ParticleMomentum[2]/CLHEP::GeV;
       fLBNEOutNtpData->tptype = NPGA->GetParticleType();
     }   
   //end find particle exiting target

  } // Partial sequence of info fill for fLBNEOutNtpData  
  
  //
  // Dk2Nu filling
  // 
  // 
  if (pRunManager->GetCreateDk2NuOutput()) {
     fDk2Nu->nuray.clear();
     fDk2Nu->ancestor.clear();
     fDk2Nu->vint.clear();
     LBNERunManager* theRunManager = dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
     fDk2Nu->job = theRunManager->GetCurrentRun()->GetRunID();
     fDk2Nu->potnum = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
     LBNETrackInformation* info = (LBNETrackInformation*)(track.GetUserInformation());
     
//     bsim::NuRay myNu(NuMomentum[0]/CLHEP::GeV, NuMomentum[1]/CLHEP::GeV, NuMomentum[2]/CLHEP::GeV, 
//                      track.GetTotalEnergy()/CLHEP::GeV, info->GetNImpWt());
     bsim::NuRay myNu(NuMomentum[0]/CLHEP::GeV, NuMomentum[1]/CLHEP::GeV, NuMomentum[2]/CLHEP::GeV, 
                      track.GetTotalEnergy()/CLHEP::GeV, 1); // Sept 24, we do not place the importance weight here.. 
//     std::cerr << " NuMomentum, Pz  in CLHEP::GeV " << NuMomentum[2]/CLHEP::GeV << " in nuray " <<    myNu.pz << std::endl;      
     fDk2Nu->nuray.push_back(myNu);
     fDk2Nu->decay.norig=Norig; // Same convention as in G4LBNE old Ntuple, i..e, as above
     fDk2Nu->decay.ndecay=NuParentTrack->GetDecayCode(); // See details in LBNESteppingAction::CheckDecay
     fDk2Nu->decay.ntype=track.GetDefinition()->GetPDGEncoding();
  
     fDk2Nu->decay.vx = x/CLHEP::cm;
     fDk2Nu->decay.vy = y/CLHEP::cm;
     fDk2Nu->decay.vz = z/CLHEP::cm;
     fDk2Nu->decay.pdpx = ParentMomentumFinal[0]/CLHEP::GeV;
     fDk2Nu->decay.pdpy = ParentMomentumFinal[1]/CLHEP::GeV;
     fDk2Nu->decay.pdpz = ParentMomentumFinal[2]/CLHEP::GeV;
     G4ThreeVector ParentMomentumProduction = NuParentTrack->GetMomentum(0);
     fDk2Nu->decay.ppdxdz = ParentMomentumProduction[0]/ParentMomentumProduction[2];
     fDk2Nu->decay.ppdydz = ParentMomentumProduction[1]/ParentMomentumProduction[2];
     fDk2Nu->decay.pppz   = ParentMomentumProduction[2]/CLHEP::GeV; 
     if (std::isnan(fDk2Nu->decay.pdpy) || (std::abs(fDk2Nu->decay.pdpy) < 1.0e-20) ||  std::isinf(fDk2Nu->decay.pdpy) || 
        std::isnan(fDk2Nu->decay.pdpx) || (std::abs(fDk2Nu->decay.pdpx) < 1.0e-20)  ||  std::isinf(fDk2Nu->decay.pdpx)){
          std::cerr << " Anomalous value for the transverse momentum of parent of neutrino " << 
           " Px " <<  fDk2Nu->decay.pdpy << "Py " << fDk2Nu->decay.pdpy << " Pz " << fDk2Nu->decay.pdpz << std::endl;
     }
     G4double parentp = sqrt(ParentMomentumProduction*ParentMomentumProduction);   
     fDk2Nu->decay.ppenergy = sqrt((parentp*parentp+Parent_mass*Parent_mass))/CLHEP::GeV;
     fDk2Nu->decay.ppmedium = NuParentTrack->GetMaterialNumber1rst();
     fDk2Nu->decay.ptype = NuParentTrack->GetPDGEncoding();
     if ((parent_name=="mu+" || parent_name=="mu-") && NuParentTrack->GetParentID()!=0)
       {
         G4int mupar                   = NuParentTrack->GetParentID();
         LBNETrajectory* MuParentTrack = GetParentTrajectory(mupar);
         G4int nopoint_mupar           = MuParentTrack->GetPointEntries();
         G4ThreeVector muparp          = MuParentTrack->GetMomentum(nopoint_mupar-1);
         G4double muparm               = MuParentTrack->GetMass();
         fDk2Nu->decay.muparpx               = muparp[0]/CLHEP::GeV; // vector of hadron parent of muon
         fDk2Nu->decay.muparpy               = muparp[1]/CLHEP::GeV; // 
         fDk2Nu->decay.muparpz               = muparp[2]/CLHEP::GeV;
         fDk2Nu->decay.mupare                = (sqrt(muparp*muparp+muparm*muparm))/CLHEP::GeV;
      }
      else
      {
       fDk2Nu->decay.muparpx = -9999999.; 
       fDk2Nu->decay.muparpy = -9999999.; 
       fDk2Nu->decay.muparpz = -9999999.;
       fDk2Nu->decay.mupare =  -9999999.; 
     }
     fDk2Nu->decay.necm = enuzrInGeV; // Now in CLHEP::GeV... 
     fDk2Nu->decay.nimpwt = info->GetNImpWt(); // duplicate info, to ease analysis. Per Robert Hatcher directive, April 24 2014.
     
     //For PPfx stuff..... A. Bashyal

    //So here  we write down the tgtexit ntuples for dk2nu. A. Bashyal 12/22/2015
    bsim::TgtExit tgt;
    G4int tar_pdg = 0;
      int Tptype = 99;
      //LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(PParentTrack->GetParticleName()));
      
      //int ptrkid = PParentTrack->GetTrackID();
      G4ThreeVector TParticleMomentum = G4ThreeVector(-999999,-999999,-999999);
      G4ThreeVector TParticlePosition = G4ThreeVector(-999999,-999999,-999999);
    
    //G4int tar_trackID = -1;
    //G4bool findtarget = false;
    //Trying to use an alternative method to find the target exit location for neutrino parents and see if it works:
     size_t ntraj = history.size();
      std::vector<G4VTrajectory*>TempHistory(history);
      std::vector<int>catchercontainer;
      std::vector<int>trackcontainer;
    for (size_t iA=0; iA != ntraj; iA++) {
    //Trying to see if I can inlcude the neutrino info in the trajectory
       LBNETrajectory *tgtraj = dynamic_cast<LBNETrajectory*>(TempHistory.at(iA)); 
       //As of june 20, 2016, it will record any particle exiting the target...this will be handled by ppfx later
       int Points = tgtraj->GetPointEntries();
       if((tgtraj->GetMaterialName(0)=="Target")&&(tgtraj->GetMaterialName(Points-1) != "Target")){
       int transit = Points-1;
       while (transit != 0){
       G4String volmat = tgtraj->GetMaterialName(transit);
       if(volmat.contains("Target")){
       catchercontainer.push_back(transit);
       trackcontainer.push_back(iA);
      // std::cout<<" looking for the loop..."<<std::endl;
       break;
       
       }
 
       else --transit;
      // std::cout<< "transit is "<<transit<<std::endl;
       }
       //Now we know the last point at which the parent of the neutrino saw the target.
       //if(transit<Points-1)std::cout<<" position and momentum "<<tgtraj->GetPoint(transit)->GetPosition()/CLHEP::cm<<std::endl;
 
       }
       }
      if(trackcontainer.size() != 0){
      int trackno = trackcontainer.at(0);
      int pointno = catchercontainer.at(0);
      LBNETrajectory *tgtraj = dynamic_cast<LBNETrajectory*>(TempHistory.at(trackno));
      TParticleMomentum = tgtraj->GetMomentum(pointno)/CLHEP::GeV;
      TParticlePosition = tgtraj->GetPoint(pointno)->GetPosition()/CLHEP::cm;
      Tptype = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(tgtraj->GetParticleName()));
      }
      tgt.tvx = TParticlePosition[0];
      tgt.tvy = TParticlePosition[1];
      tgt.tvz = TParticlePosition[2];
      tgt.tpx = TParticleMomentum[0];
      tgt.tpy = TParticleMomentum[1];
      tgt.tpz = TParticleMomentum[2];
      tgt.tptype = Tptype;
     // std::cout<<"The tar_pdg is "<<Tptype<<" and momentum is "<<TParticleMomentum[2]/CLHEP::GeV<<" "<<tgt.tvz<<std::endl;
     //std::cout<<"tgtexit : tvz "<<tgt.tvz<<std::endl;
      //tgt.tgen = info->GetTgen()-1;
      fDk2Nu->tgtexit = tgt;
        
    // ancestry info. 
     // start by counting the depth of the tree.
     size_t numAncestors = 0; 
     G4int trackIDTmp = track.GetParentID();
     LBNETrajectory *tmpTraj = GetTrajectory(trackIDTmp);
     const LBNEVolumePlacements *aPlacementHandler = LBNEVolumePlacements::Instance();
     std::vector<G4String> Horn2ICList = aPlacementHandler->GetHorn2ICList();
     std::vector<G4String>Horn1ICList = aPlacementHandler->GetHorn1ICList();     
     while (trackIDTmp > 0) {
       numAncestors++;
       trackIDTmp = tmpTraj->GetParentID();
       if(trackIDTmp > 0) tmpTraj = GetTrajectory(trackIDTmp);
     }
     std::vector<LBNETrajectory *> trajs(numAncestors, 0);
     // Revert order, starting with the proton...
     size_t nAnces = numAncestors;
     trackIDTmp = track.GetParentID();
     tmpTraj = GetTrajectory(trackIDTmp);
     while (trackIDTmp > 0) {
       nAnces--;
       trajs[nAnces] = tmpTraj;
       trackIDTmp = tmpTraj->GetParentID();
       if(trackIDTmp > 0) tmpTraj = GetTrajectory(trackIDTmp);
       }
     // Now fill.. 
     fDk2Nu->ancestor.clear();
     fDk2Nu->vint.clear();
     //A. Bashyal...few more additions for ppfx.
     size_t nmax = 0;
     int idx_tar_in_chain = -1;
     size_t ntrajectory = history.size();
      std::vector<G4VTrajectory*>TmpHistory(history);
      std::reverse(TmpHistory.begin(),TmpHistory.end());
    // for (size_t iA=0; iA != numAncestors; iA++) {
    //Trying to see if I can inlcude the neutrino info in the trajectory
       for (size_t iA=0;iA<ntrajectory;++iA){
       bsim::Ancestor a;
       //LBNETrajectory *t = trajs[iA];
       LBNETrajectory *t = dynamic_cast<LBNETrajectory*>(TmpHistory.at(iA));
       a.pdg = t->GetPDGEncoding();
       //if(t->GetPDGEncoding()==14)std::cout<<"Muon neutrino found "<<std::endl;
       G4ThreeVector x1rst = t->GetPoint(0)->GetPosition();
       a.startx = x1rst[0]/CLHEP::cm; 
       a.starty = x1rst[1]/CLHEP::cm; 
       a.startz = x1rst[2]/CLHEP::cm;
       a.startt = t->GetTimeStart(); // ? Units? 
       G4ThreeVector p1rst = t->GetMomentum(0);
       a.startpx = p1rst[0]/CLHEP::GeV; 
       a.startpy = p1rst[1]/CLHEP::GeV; 
       a.startpz = p1rst[2]/CLHEP::GeV;
       G4int np = t->GetPointEntries();
       G4ThreeVector pLast = t->GetMomentum(np-1);
       a.stoppx = pLast[0]/CLHEP::GeV;
       a.stoppy = pLast[1]/CLHEP::GeV;
       a.stoppz = pLast[2]/CLHEP::GeV;
       G4ThreeVector pol = t->GetPolarization();
       a.polx = pol[0];
       a.poly = pol[1];
       a.polz = pol[2];
       // We do not fill pprodpx, y,z, duplicate info 
       //a.pprodpx = 0.; a.pprodpy = 0.; a.pprodpz = 0.;
       // prodpxyz ntuples are filled as per the requirement in ppfx.....A. Bashyal
        //As per the NumiAnalysis, the prodpx saves the momentum of the 
       //particle before it produces the secondaries. I am assuming this is 
       //equilavent to the particle momentum before the particle decays into 
       //whatever it has to.
       G4ThreeVector momatProd = t->GetParentMomentumAtThisProduction();
       a.pprodpx = momatProd[0]/CLHEP::GeV;
       a.pprodpy = momatProd[1]/CLHEP::GeV;
       a.pprodpz = momatProd[2]/CLHEP::GeV;
       //std::cout<<"pprodpz "<<momatProd/CLHEP::GeV<<std::endl;
       //if((t->GetPDGEncoding() != 211) && (momatProd[2]/CLHEP::GeV == 0)std::cout<<" 0 momentum pion "<<t->GetProcessName()<<std::endl;
       a.nucleus = t->GetPDGNucleus();
       a.proc = t->GetProcessName();
       //a.ivol = t->GetVolName1rst();
       G4String volname = t->GetVolName1rst();
       for(size_t j = 0;j != Horn1ICList.size();j++)
       {
        if(volname==Horn1ICList.at(j))volname="Horn1IC";
       }
       for(size_t j = 0;j!=Horn2ICList.size();j++){
       if(volname==Horn2ICList.at(j))volname="Horn2IC";
       }
       //if(volname=="TargetNoSplitSegmentLeft"||volname=="TargetNoSplitSegmentRight")volname="TargetNoSplitSegment"; 
       //below changes are made because of the naming convention in the ppfx-ral
        if(volname=="TargetUpstrDownstrSegmentRight"||volname=="TargetUpstrDownstrSegmentLeft"||volname=="RAL-Target"){
       volname="TargetNoSplitSegment";    //A. bashyal at leaset upto v3r5c, reference design      
       }
       a.ivol = volname;
      // std::cout<<"volname "<<volname<<std::endl;
       a.imat = t->GetMaterialName1rst();
       //More ppfx stuff.....A. Bashyal
       //for(int ap = 0;ap<np-1;ap++){if(t->GetMaterialName(ap).contains("Target"))idx_tar_in_chain = int(iA);}
       //Okay I think I did this part wrong. Should only look at the end point of the trajectory to determine this.
       if(t->GetMaterialName(np-1).contains("Target"))idx_tar_in_chain = int (iA);
       fDk2Nu->ancestor.push_back(a);
     }
     // 
       fDk2Nu->vint.push_back(idx_tar_in_chain); //ppfx stuff A. Bashyal...gotta see this one more time..
       //Just for the sake of consistency with g4numi I am going to feed in a dummy variable. A. Bashyal
       fDk2Nu->vint.push_back(-1); //-1 since we are not shifting the target (see TargetAttenuationReweighter.cpp
         //Few stuff on ppfx distance travelled through out volume of interest A. Bashyal
        fDk2Nu->vdbl.clear();


        const int Nanc = 3; //parent, granparent and greatgranparent for ancestry
        const double fact_Al = (26.98/2.7);
        const double fact_steel316 = (52.73/8.0);
        const double fact_concrete = (33.63/2.03);
        const double fact_water = (18.01/1.0);
        const double fact_iron = (207.19/11.35);
        const double fact_helium = (4.003/0.000145);//This from G4numi
   
        for(int ii = 0;ii<Nanc;ii++){
          dist_IC1[ii] = -1*fact_Al;
          dist_IC2[ii] = -1*fact_Al;
          dist_DPIP[ii] = -1*fact_steel316;
          dist_DVOL[ii] = -1*fact_helium;
         }
         
        if(history.size()!=0){
        
          for(int ii = 0;ii<Nanc;ii++){
          dist_IC1[ii] = 0;
          dist_IC2[ii] = 0;
          dist_DPIP[ii] =0;
          dist_DVOL[ii] = 0;
         }      
        
        }

       LBNETrajectory* temp_traj;
       //we need two different kinds of initializations of these things
        //for (size_t ii = 0;ii<std::min(size_t(3),trajs.size());++ii){ //This basically gives parent at 0, greatgranparent at 1 and so on.
        for(int ii = 1; ii<4;ii++){
        if(history.size()-ii<1)continue;
        temp_traj = dynamic_cast<LBNETrajectory*>(TmpHistory.at(ii)); //parsing in the distance travelled info

        double horn1IC = 0;
        double horn2IC = 0;
        //implementation...yeah!!
        for(size_t i = 0;i != Horn1ICList.size();i++)horn1IC += GetDistanceInVolume(temp_traj,Horn1ICList.at(i));
        for(size_t i = 0;i != Horn2ICList.size();i++)horn2IC += GetDistanceInVolume(temp_traj,Horn2ICList.at(i));
        
        //std::cout<<"Total IC1 and IC2 distance "<<horn1IC/fact_Al<<"   "<<horn2IC/fact_Al<<std::endl;
        
        if(horn1IC != 0.)dist_IC1[ii]=horn1IC/fact_Al;
        
        if(horn2IC != 0.)dist_IC2[ii]=horn2IC/fact_Al;

         dist_DPIP[ii] = GetDistanceInVolume(temp_traj,"DecayPipeWall")/fact_steel316;
         dist_DVOL[ii] = GetDistanceInVolume(temp_traj,"DecayPipeVolume")/fact_helium;
         
    }
        
        for(int ii=0;ii<3;ii++){
        fDk2Nu->vdbl.push_back(dist_IC1[ii]);
                                        
        fDk2Nu->vdbl.push_back(dist_IC2[ii]);

        fDk2Nu->vdbl.push_back(dist_DPIP[ii]);

        fDk2Nu->vdbl.push_back(dist_DVOL[ii]);
                                }
        
     
         
         
        
        //if(dist_DVOL[ii] != 0)std::cout<<dist_DVOL[ii]<<" and level is "<<ii<<std::endl;
        
        
        //Now we write the traj branches....A. Bashyal December 17 2015
                bsim::Traj tmp_traj; //Declaration outside the loop??
        double trkx[10];
        double trky[10];
        double trkz[10];
        double trkpx[10];
        double trkpy[10];
        double trkpz[10];

      for (int ii = 0;ii<10;++ii)
      {
      //bsim::Traj tmp_traj; //Declaration outside the loop??
       trkx[ii]=-99999;
       trky[ii]=-99999;
       trkz[ii]=-99999;
       trkpx[ii]=-99999;
       trkpy[ii]=-99999;
       trkpz[ii]=-99999;
       }
       G4int npoint          = NuParentTrack->GetPointEntries();
       //std::cout<<"Number of iterations are :"<<npoint<<std::endl;
       G4ThreeVector ParentMomentum;
       G4ThreeVector ParentPosition;
      // if(h2entryindex>2){std::cout<<" For Particles in horn2 "<<h2entryindex<<" "<<NuParentTrack->GetPreStepVolumeName(h2entryindex)<<
        // h2entryindex+1<<" "<<NuParentTrack->GetPreStepVolumeName(h2entryindex+1)<<
        // " "<<h2entryindex-1<<NuParentTrack->GetPreStepVolumeName(h2entryindex-1)<<std::endl;}
         std::vector<int> h2exitindex;
         size_t points = 0;
         points = size_t(npoint);
         h2exitindex.clear();
         for(size_t ii = 0;ii != points-1;ii++){
         int iii = int(ii);
           G4String h2volname = NuParentTrack->GetPreStepVolumeName(iii);
         if(h2volname.contains("Horn2")&&(NuParentTrack->GetMaterialName(iii).contains("Aluminum"))){
            h2exitindex.push_back(iii);
         
         
         }
         }
         int h2exit = -1;
         int h2entry = -1;
        if(!h2exitindex.empty())h2exit = h2exitindex.back();
        if(!h2exitindex.empty())h2entry = h2exitindex.front();
        //if(h2exit>2&&h2exit<npoint){std::cout<<" For Particles in horn2 "<<h2exit<<" "<<NuParentTrack->GetPreStepVolumeName(h2exit)<<
        // h2exit+1<<" "<<NuParentTrack->GetPreStepVolumeName(h2exit+1)<<
        // " "<<h2exit-1<<NuParentTrack->GetPreStepVolumeName(h2exit-1)<<std::endl;}
      for(G4int ii = 0; ii<npoint-1;ii++){
         ParentMomentum = NuParentTrack->GetMomentum(ii);
         ParentPosition = NuParentTrack->GetPoint(ii)->GetPosition();
         G4String postvolname = " ";
         G4String prevolname = NuParentTrack->GetPreStepVolumeName(ii);
         if(ii<npoint-2) postvolname = NuParentTrack->GetPreStepVolumeName(ii+1);

        if((prevolname.contains("TargetNoSplitSegment")||prevolname.contains("TargetFinHorizontal"))&&ii==0){
        trkx[0] = ParentPosition[0]/CLHEP::cm;
        trky[0] = ParentPosition[1]/CLHEP::cm;
        trkz[0] = ParentPosition[2]/CLHEP::cm;
        trkpx[0] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[0] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[0] = ParentMomentum[2]/CLHEP::GeV;
          
    }
    if(prevolname.contains("TargetNoSplitM1")&&postvolname.contains("TargetHallAndHorn1")){
    //  if(ii==Points){ //I meant to keep it exiting the target but requires more detailed arguments
                        //std::cout<<"trkx1 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;  
                        trkx[1] = ParentPosition[0]/CLHEP::cm;
                        trky[1] = ParentPosition[1]/CLHEP::cm;
                        trkz[1] = ParentPosition[2]/CLHEP::cm;
                        trkpx[1] = ParentMomentum[0]/CLHEP::GeV;
                        trkpy[1] = ParentMomentum[1]/CLHEP::GeV;
                        trkpz[1] = ParentMomentum[2]/CLHEP::GeV;
                
                }
                //if(prevolname.contains("Target"))std::cout<<" Exiting target Prevolname is "<<prevolname<<" and postvol name is "<<
                //postvolname<<" at xyz = "<<ParentPosition[0]/CLHEP::cm<<" "<<ParentPosition[1]/CLHEP::cm<<" "<<
                //ParentPosition[2]/CLHEP::cm<<std::endl;
        if(prevolname.contains("TargetHallAndHorn1") && postvolname.contains("Horn1PolyM1"))
          //  && !(postvolname.contains("TargetHall")))
        {
      //  std::cout<<"trkx2 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;
        trkx[2] = ParentPosition[0]/CLHEP::cm;
        trky[2] = ParentPosition[1]/CLHEP::cm;
        trkz[2] = ParentPosition[2]/CLHEP::cm;
        trkpx[2] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[2] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[2] = ParentMomentum[2]/CLHEP::GeV;
        }
        //if(postvolname.contains("Horn2"))std::cout<<" Postvolname contains Horn2 "<< postvolname<< " and prevol name is "<<
        //      prevolname<<" at z = "<<ParentPosition[2]/CLHEP::cm<<std::endl;
        if(prevolname.contains("Horn1PolyM1")&& postvolname.contains("TargetHallAndHorn1"))
        {     
        //std::cout<<"trkx3 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<" and "<<ParentPosition[2]/CLHEP::cm<<std::endl;  
        trkx[3] = ParentPosition[0]/CLHEP::cm;
        trky[3] = ParentPosition[1]/CLHEP::cm;
        trkz[3] = ParentPosition[2]/CLHEP::cm;
        trkpx[3] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[3] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[3] = ParentMomentum[2]/CLHEP::GeV;
        }
        if(ii==h2entry)
        { //need to make sure we dont include the particle exiting horn2 condition.
        //std::cout<<"trkx4 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;  
        trkx[4] = ParentPosition[0]/CLHEP::cm;
        trky[4] = ParentPosition[1]/CLHEP::cm;
        trkz[4] = ParentPosition[2]/CLHEP::cm;
        trkpx[4] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[4] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[4] = ParentMomentum[2]/CLHEP::GeV;
        //std::cout<<" Entering from "<<prevolname<<" Going to "<<postvolname<<" at "<<
        //ParentPosition/CLHEP::cm<<std::endl;
        }
        if(ii==h2exit)
        {
        //std::cout<<"trkx5 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;  
        trkx[5] = ParentPosition[0]/CLHEP::cm;
        trky[5] = ParentPosition[1]/CLHEP::cm;
        trkz[5] = ParentPosition[2]/CLHEP::cm;
        trkpx[5] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[5] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[5] = ParentMomentum[2]/CLHEP::GeV;
        //std::cout<<"Entering from"<<prevolname<<" Going to "<<postvolname<<" at "<<ParentPosition/CLHEP::cm<<std::endl;
        }
        if(prevolname.contains("DecayPipeHall")&&postvolname.contains("DecayPipeVolume"))
        {
        //std::cout<<"trkx6 "<<ParentPosition[0]/CLHEP::cm<<" from "<<prevolname<<" to "<<postvolname<<std::endl;  
        trkx[6] = ParentPosition[0]/CLHEP::cm;
        trky[6] = ParentPosition[1]/CLHEP::cm;
        trkz[6] = ParentPosition[2]/CLHEP::cm;
        trkpx[6] = ParentMomentum[0]/CLHEP::GeV;
        trkpy[6] = ParentMomentum[1]/CLHEP::GeV;
        trkpz[6] = ParentMomentum[2]/CLHEP::GeV;
        }

        const double NomH1InnerRad = aPlacementHandler->GetHorn1NeckInnerRadius()/CLHEP::cm;
        const double NomH1NeckZloc = aPlacementHandler->GetHorn1NeckZPosition()/CLHEP::cm;
        double OptH1InnerRad = 0.0;
        double OptH1NeckZLoc = 0.0;
        if(aPlacementHandler->GetUseHorn1Polycone()||aPlacementHandler->GetNumberOfHornsPolycone() >0){
        OptH1InnerRad = 
        aPlacementHandler->GetHorn1PolyInnerConductorRadius(1)/CLHEP::cm;
        OptH1NeckZLoc = aPlacementHandler->GetHorn1PolyInnerConductorZCoord(3)/CLHEP::cm; //This way the target interactions are not ignored
        } //1 for the inner neck radius...see DUNECDR_Optimized.mac in macros.
        double H1NeckInnerRad = 0.0;
        double H1NeckZLoc = 0.0;
        if(NomH1InnerRad != 0){
        H1NeckInnerRad = NomH1InnerRad;
        H1NeckZLoc = NomH1NeckZloc;
        }
        else {
        H1NeckInnerRad = OptH1InnerRad;
        H1NeckZLoc = OptH1NeckZLoc;
        
        
        }
        //std::cout<<H1NeckZLoc<< " is neck NLocation"<<std::endl;
        //const double H1NeckInnerRad = (aPlacementHandler->GetHorn1NeckInnerRadius())/CLHEP::cm;
        const double H2NeckInnerRad = (aPlacementHandler->GetHorn2NeckInnerRadius())/CLHEP::cm;
        //std::cout<<"Horn1neck and horn2 neck inner rad are" << H1NeckInnerRad <<" "<<H2NeckInnerRad<<std::endl;
//Exclude the particles through the neck of the horns.
        if((sqrt(ParentPosition[0]*ParentPosition[0]+ParentPosition[1]*ParentPosition[1])<H1NeckInnerRad)&&ParentPosition[2]>H1NeckZLoc)
        {
        trkx[2]=99999;
        trky[2]=99999;
        trkz[2]=99999;
        }
        if(sqrt(ParentPosition[0]*ParentPosition[0]+ParentPosition[1]*ParentPosition[1])<H2NeckInnerRad)
        {
        trkx[4]=99999;
        trky[4]=99999;
        trkz[4]=99999;
        }
        }
            
  ParentMomentum = NuParentTrack->GetMomentum(npoint-1);
  ParentPosition = (NuParentTrack->GetPoint(npoint-1)->GetPosition()/CLHEP::m)*CLHEP::m;
  trkx[7] = ParentPosition[0]/CLHEP::cm;
  trky[7] = ParentPosition[1]/CLHEP::cm;
  trkz[7] = ParentPosition[2]/CLHEP::cm;
  trkpx[7] = ParentMomentum[0]/CLHEP::GeV;
  trkpy[7] = ParentMomentum[1]/CLHEP::GeV;
  trkpz[7] = ParentMomentum[2]/CLHEP::GeV;


     fDk2Nu->traj.clear();
     
    for(G4int ii = 0;ii<10;++ii)
        {
        //bsim::Traj tmp_traj;
        tmp_traj.trkx = trkx[ii];
        tmp_traj.trky = trky[ii];
        tmp_traj.trkz = trkz[ii];
        tmp_traj.trkpx = trkpx[ii];
        tmp_traj.trkpy = trkpy[ii];
        tmp_traj.trkpz = trkpz[ii];
        fDk2Nu->traj.push_back(tmp_traj);
        }
      
     // Last step : compute the so-called "location weights" for 3 detectors. 
     //
     
     
     CalcLocationWeights(fDkMeta, fDk2Nu,pRunManager->GetUseRealisticNearDetectorVolume());  
     // 
     // Check the weights..
     //
     for (size_t kLoc = 0; kLoc != fDkMeta->location.size(); kLoc++) {
       if ((!std::isnan(fDk2Nu->nuray[kLoc].E)) && (!std::isnan(fDk2Nu->nuray[kLoc].wgt))) continue;
       std::cerr << " LBNEAnalysis::FillNeutrinoNtuple Problem with detector weight calculation for location "
                << kLoc << std::endl;
                
       std::cerr << " Decay Vertex : X = " << fDk2Nu->decay.vx << " Y = " 
                 << fDk2Nu->decay.vy << " Z = " << fDk2Nu->decay.vz << std::endl;
       std::cerr << " Original track position " << pos.x() << " / " << pos.y() << " / " << pos.z() << std::endl;         
       std::cerr << " Decay Momentum : " << fDk2Nu->decay.pdpx 
                 << " / " <<  fDk2Nu->decay.pdpy << " / " << fDk2Nu->decay.pdpy << std::endl;
       std::cerr << " Neutrino from G4 decay momentum " << NuMomentum[0]/CLHEP::GeV
                << "/ " <<  NuMomentum[1]/CLHEP::GeV << " / " << NuMomentum[2]/CLHEP::GeV << std::endl;
       std::cerr << " decay type " << fDk2Nu->decay.ptype << " parent name " << parent_name <<std::endl;
                       
     }
//     std::cerr << " Weight and energy for E [G4]  = " << track.GetTotalEnergy()/CLHEP::GeV << " CLHEP::GeV " << std::endl;
//     for (size_t kDet = 0; kDet != 3; kDet++ ) { 
//       std::cerr << " E [CLHEP::GeV] " << fDk2Nu->nuray[kDet].E << " Weight " << fDk2Nu->nuray[kDet].wgt << std::endl; 
//     } 
  } // on filling the Dk2nu Ntuple. 
  
   unsigned int ndets = fXdet_near.size();
   for(unsigned int idet = 0; idet < ndets; ++idet)
   if (pRunManager->GetCreateOutput()) {
      fLBNEOutNtpData->NdxdzNear[idet] = (x-fXdet_near[idet])/(z-fZdet_near[idet]);
      fLBNEOutNtpData->NdydzNear[idet] = (y-fYdet_near[idet])/(z-fZdet_near[idet]);
      
      LBNENuWeight nuwgh;
      G4double nu_wght;
      G4double nu_energy;
      std::vector<double> r_det;
      r_det.push_back(fXdet_near[idet]/CLHEP::cm);
      r_det.push_back(fYdet_near[idet]/CLHEP::cm);
      r_det.push_back(fZdet_near[idet]/CLHEP::cm);
      nuwgh.GetWeight(fLBNEOutNtpData, r_det,nu_wght,nu_energy);
      fLBNEOutNtpData->NenergyN[idet] = nu_energy; //in CLHEP::GeV
      fLBNEOutNtpData->NWtNear[idet]  = nu_wght;
   }
   //end near det

   //
   //Far Detector
   
   if(fXdet_far.size() != fYdet_far.size() || 
      fXdet_far.size() != fYdet_far.size() ||
      fYdet_far.size() != fZdet_far.size())
   {
      G4cout << "LBNEAnalysis::FillNeutrinoNtuple - "
             << "Far Detector vectors are not the same size." << G4endl; 
   }

//   const int ntypeBefWeight = fLBNEOutNtpData->Ntype;  Old statement, lost history... Paul Lebrun, Feb. 2015.     
   
   ndets = fXdet_far.size();
   if (pRunManager->GetCreateOutput()) {
     for(unsigned int idet = 0; idet < ndets; ++idet)
     {
        fLBNEOutNtpData->NdxdzFar[idet] = (x-fXdet_far[idet])/(z-fZdet_far[idet]);
        fLBNEOutNtpData->NdydzFar[idet] = (y-fYdet_far[idet])/(z-fZdet_far[idet]);
      
        LBNENuWeight nuwgh;
        G4double nu_wght;
        G4double nu_energy;
        std::vector<double> r_det;
        r_det.push_back(fXdet_far[idet]/CLHEP::cm);
        r_det.push_back(fYdet_far[idet]/CLHEP::cm);
        r_det.push_back(fZdet_far[idet]/CLHEP::cm);
        nuwgh.GetWeight(fLBNEOutNtpData, r_det,nu_wght,nu_energy);
        fLBNEOutNtpData->NenergyF[idet] = nu_energy; //in CLHEP::GeV
        fLBNEOutNtpData->NWtFar[idet]   = nu_wght;
     }
    }
    
    
    //Horn info starts here Amit Bashyal
if(pRunManager->GetCreateOutput()){

   G4int numberOfPoints            = NuParentTrack->GetPointEntries();
  // std::cout<<NuParentTrack->GetPDGEncoding()<<std::endl;
   G4ThreeVector ParticlePosition1 = NuParentTrack->GetPoint(0)->GetPosition();
     // std::cout<<ParticlePosition1/CLHEP::cm<<std::endl;
   for (G4int ii=numberOfPoints-1; ii > -1; --ii)
   {
      
      G4String prevolname = NuParentTrack->GetPreStepVolumeName(ii);
      
      G4String postvolname = "";
      if(ii < numberOfPoints-1) postvolname = NuParentTrack->GetPreStepVolumeName(ii+1);
      //std::cout<<postvolname<<std::endl;
      G4ThreeVector ParticleMomentum = NuParentTrack->GetMomentum(ii);  
      //std::cout<<ParticleMomentum<<"   "<<"eventid"<<fLBNEOutNtpData->evtno<<"   "<<"trackid "<<NuParentTrack->GetTrackID()<<std::endl;            
      G4ThreeVector ParticlePosition = NuParentTrack->GetPoint(ii)->GetPosition();
      if (((numberOfPoints - ii) < 3) || (ii < 2))
      {
      //++count;
      //if(count > 1)
//      std::cout << "eventid: " << fLBNEOutNtpData->evtno << " trackid: " << trackID << " ii = " << ii 
//      << " particle = " << TrackTrajectory->GetParticleName()
//      << " prestepvolname = " << prevolname << " postvolname = " << postvolname 
//      << " steplength = " << steplength/mm << " mm" << std::endl;
      }
    
           //
      //particle created in the target
      // For the neutrino ancestors passing through Horn 1 end plane Amit Bashyal
      if (prevolname.contains("Horn1TrackingPlane"))
      {
         fLBNEOutNtpData -> h1posx = ParticlePosition[0]/CLHEP::cm;
         fLBNEOutNtpData -> h1posy = ParticlePosition[1]/CLHEP::cm;
         fLBNEOutNtpData -> h1posz = ParticlePosition[2]/CLHEP::cm;
         fLBNEOutNtpData -> h1momx = ParticleMomentum[0]/CLHEP::GeV;
         fLBNEOutNtpData -> h1momy = ParticleMomentum[1]/CLHEP::GeV;
         fLBNEOutNtpData -> h1momz = ParticleMomentum[2]/CLHEP::GeV;
         fLBNEOutNtpData ->h1trackid = NuParentTrack->GetTrackID();
      }
      
      if(prevolname.contains("Horn2TrackingPlane"))
      { //For the Neutrino Ancestors passing through Horn 2 end plane Amit Bashyal
         fLBNEOutNtpData -> h2posx = ParticlePosition[0]/CLHEP::cm;
         fLBNEOutNtpData -> h2posy = ParticlePosition[1]/CLHEP::cm;
         fLBNEOutNtpData -> h2posz = ParticlePosition[2]/CLHEP::cm;
         fLBNEOutNtpData -> h2momx = ParticleMomentum[0]/CLHEP::GeV;
         fLBNEOutNtpData -> h2momy = ParticleMomentum[1]/CLHEP::GeV;
         fLBNEOutNtpData -> h2momz = ParticleMomentum[2]/CLHEP::GeV;
         fLBNEOutNtpData ->h2trackid = NuParentTrack->GetTrackID();
      }
      
      
      }
      }
      
   G4int trackID                   = track.GetParentID();
   LBNETrajectory* TrackTrajectory = GetTrajectory(trackID);         
   

   while (trackID > 0)
   {
      LBNEAnalysis::TrackThroughGeometry(TrackTrajectory);
      
      trackID = TrackTrajectory->GetParentID();
      if(trackID > 0) TrackTrajectory = GetTrajectory(trackID);      

//       std::cerr << " Fill Ntuple, Trajectory access for track ancestry " << trackID 
//            << " From trajectory " << TrackTrajectory->GetParticleDefinition()->GetParticleName() 
//            << " Parent Track ID " <<   TrackTrajectory->GetParentID()
//            << " num Points " << TrackTrajectory->GetPointEntries() << std::endl;
                

   }//end while trackID > 0
   //end tracking through geometry

  //
  //Done. Fill Tree.
  //
  if (pRunManager->GetCreateOutput()) fOutTree->Fill();  
  if (pRunManager->GetCreateDk2NuOutput()) {
      fOutTreeDk2Nu->Fill(); 
  } 
  //
  // Write to ascii file
  //
  if (pRunManager->GetCreateASCIIOutput())
  {
     std::string asciiFileName = pRunManager->GetOutputASCIIFileName();
     std::ofstream asciiFile(asciiFileName.c_str(), std::ios::app);
     if(asciiFile.is_open())
     {
        asciiFile << fLBNEOutNtpData->Ntype<< " " << fLBNEOutNtpData->Nenergy << " " 
                  << fLBNEOutNtpData->NenergyN[0] << " " << fLBNEOutNtpData->NWtNear[0];
        asciiFile << " " << fLBNEOutNtpData->NenergyF[0] << " " << fLBNEOutNtpData->NWtFar[0] 
                  <<" "<<fLBNEOutNtpData->Nimpwt<< G4endl; 
        asciiFile.close();
     }
  }

}

void LBNEAnalysis::FillTargetOutputData

(

const G4Step &

aStep

)

Definition at line 2091 of file LBNEAnalysis.cc.

                                                         {
  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
 if (pRunManager->GetCreateTargetOutput()){

  if(fTNParticles < HkMaxP){
    G4Track *track = step.GetTrack();
    // We are ignoring the Importance weight here... Why.. ?  Paul L. G. Lebrun. 
   // LBNETrajectory *ptrack = (LBNETrajectory*)(track->GetTrajectory());
//    LBNETrackInformation *info =
//    (LBNETrackInformation*)(track->GetUserInformation());
//    if(info != 0){
  //    fTNImpWt = info->GetNImpWt();
   // }
    G4String prematerial = step.GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();
    G4String posmaterial = step.GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();  
    G4StepPoint *point = step.GetPostStepPoint();
    G4ThreeVector pos = point->GetPosition();
    G4ThreeVector mom = point->GetMomentum();
    fTParticlePX = mom.x()/CLHEP::GeV;
    fTParticlePY = mom.y()/CLHEP::GeV;
    fTParticlePZ = mom.z()/CLHEP::GeV;
    fTParticleX = pos.x()/CLHEP::cm;
    fTParticleY = pos.y()/CLHEP::cm;
    fTParticleZ = pos.z()/CLHEP::cm;
    G4ParticleDefinition *def = track->GetDefinition();
    fTParticlePDG = def->GetPDGEncoding();
    fTParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
    fTTrackID = track->GetTrackID();
    fTParentID = track->GetParentID();
    fTNParticles++;
  }
 // fTargetOutputTree->Fill();


 }// end if track planes is on

}

void LBNEAnalysis::FillTrackingNtuple	(	const G4Track &	track,
		LBNETrajectory *	currTrajectory
	)

Definition at line 1517 of file LBNEAnalysis.cc.

{
   
   LBNERunManager* pRunManager = 
       dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());
       
   if (!pRunManager->GetCreateOutput()) return;

   fLBNEOutNtpData -> Clear();

   LBNEPrimaryGeneratorAction *NPGA = 
      (LBNEPrimaryGeneratorAction*)(pRunManager)->GetUserPrimaryGeneratorAction();




   fLBNEOutNtpData->protonN    = NPGA->GetProtonNumber();
   
   fLBNEOutNtpData->run        = pRunManager->GetCurrentRun()->GetRunID();
   fLBNEOutNtpData->evtno      = pRunManager->GetCurrentEvent()->GetEventID();
   fLBNEOutNtpData->beamHWidth = NPGA->GetBeamSigmaX()/CLHEP::cm;
   fLBNEOutNtpData->beamVWidth = NPGA->GetBeamSigmaY()/CLHEP::cm;
   fLBNEOutNtpData->beamX      = NPGA->GetBeamOffsetX()/CLHEP::cm;
   fLBNEOutNtpData->beamY      = NPGA->GetBeamOffsetY()/CLHEP::cm;
   
   //G4int particleID = track.GetParentID();
   G4ThreeVector protonOrigin = NPGA->GetProtonOrigin();
   fLBNEOutNtpData->protonX            = protonOrigin[0];
   fLBNEOutNtpData->protonY            = protonOrigin[1];
   fLBNEOutNtpData->protonZ            = protonOrigin[2];
   
   G4ThreeVector protonMomentum = NPGA->GetProtonMomentum();
   fLBNEOutNtpData->protonPx             = protonMomentum[0];
   fLBNEOutNtpData->protonPy             = protonMomentum[1];
   fLBNEOutNtpData->protonPz             = protonMomentum[2];
   
   LBNEVolumePlacements *volDB = LBNEVolumePlacements::Instance();
   
   fLBNEOutNtpData->nuTarZ      = volDB->GetTargetLengthIntoHorn(); // A better info that the somewhat meaningless Z0 
   const LBNEDetectorConstruction *detDB = 
      dynamic_cast<const LBNEDetectorConstruction*>(pRunManager->GetUserDetectorConstruction());
   fLBNEOutNtpData->hornCurrent = detDB->GetHornCurrent()/CLHEP::ampere/1000.;

   G4int trackID                   = track.GetTrackID();
   G4String currentTrackVolName    = track.GetVolume() -> GetName();
   LBNETrajectory* TrackTrajectory = currTrajectory;

   //G4int eventID                   =  G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
   //std::cout << "Event: " << eventID << " TrackID: " << trackID << " curr vol = " << currentTrackVolName << std::endl;


   while (trackID > 0)
   {
      
      LBNEAnalysis::TrackThroughGeometry(TrackTrajectory);

      trackID = TrackTrajectory->GetParentID();
      if(trackID > 0) TrackTrajectory = GetTrajectory(trackID);

   }
  
   fOutTree->Fill();  
   
}

void LBNEAnalysis::FillTrackingPlaneData

(

const G4Step &

aStep

)

Definition at line 1853 of file LBNEAnalysis.cc.

                                                          {

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetCreateTrkPlaneOutput()){

  if(fNParticles < kMaxP){
    G4Track *track = step.GetTrack();
    LBNETrackInformation *info =
    (LBNETrackInformation*)(track->GetUserInformation());
    if(info != 0){
      fNImpWt[fNParticles] = info->GetNImpWt();
    }
    G4StepPoint *point = step.GetPostStepPoint();
    G4ThreeVector pos = point->GetPosition();
    G4ThreeVector dir = point->GetMomentumDirection();
    G4ParticleDefinition *def = track->GetDefinition();
    fParticleX[fNParticles] = pos.x()/CLHEP::m;
    fParticleY[fNParticles] = pos.y()/CLHEP::m;
    fParticleZ[fNParticles] = pos.z()/CLHEP::m;
    fParticleDX[fNParticles] = dir.x();
    fParticleDY[fNParticles] = dir.y();
    fParticleDZ[fNParticles] = dir.z();
    fParticlePDG[fNParticles] = def->GetPDGEncoding();
    fParticleMass[fNParticles] = def->GetPDGMass();
    fParticleEnergy[fNParticles] = point->GetKineticEnergy()/CLHEP::GeV;
    fTrackID[fNParticles] = track->GetTrackID();
    fNParticles++;
  }
  // std::cout << "I stored muon data" << std::endl;

 }// end if track planes is on

}

void LBNEAnalysis::FillTrackingPlaneDPData

(

const G4Step &

aStep

)

Definition at line 2039 of file LBNEAnalysis.cc.

                                                            {

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetCreateTrkPlaneDPOutput()){

    G4Track *track = step.GetTrack();
   // LBNETrajectory *ptrack = (LBNETrajectory*)(track->GetTrajectory());
    LBNETrackInformation *info =
    (LBNETrackInformation*)(track->GetUserInformation());
    if(info != 0){
      fDPNImpWt = info->GetNImpWt();
    }
    G4StepPoint *point = step.GetPostStepPoint();
    G4ThreeVector pos = point->GetPosition();
    G4ThreeVector mom = point->GetMomentum();
    const G4ThreeVector ppoint = track->GetVertexPosition();
    const G4ThreeVector pmom = track->GetVertexMomentumDirection();
    fDPPProductionDX = pmom.x();
    fDPPProductionDY = pmom.y();
    fDPPProductionDZ = pmom.z();
    fDPPProductionX = ppoint.x()/CLHEP::cm;
    fDPPProductionY = ppoint.y()/CLHEP::cm;
    fDPPProductionZ = ppoint.z()/CLHEP::cm;
    fDPParticlePX = mom.x()/CLHEP::GeV;
    fDPParticlePY = mom.y()/CLHEP::GeV;
    fDPParticlePZ = mom.z()/CLHEP::GeV;
    fDPParticlePXPZ = mom[0]/mom[2];
    fDPParticlePYPZ = mom[1]/mom[2]; 
    G4ThreeVector dir = point->GetMomentumDirection();
    G4ParticleDefinition *def = track->GetDefinition();
    fDPParticleX = pos.x()/CLHEP::cm;
    fDPParticleY = pos.y()/CLHEP::cm;
    fDPParticleZ = pos.z()/CLHEP::cm;
    fDPParticleDX = dir.x();
    fDPParticleDY = dir.y();
    fDPParticleDZ = dir.z();
    fDPParticlePDG = def->GetPDGEncoding();
    fDPParticleMass = def->GetPDGMass();
    fDPParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
    fDPTrackID = track->GetTrackID();
    fDPParentID = track->GetParentID();
    fDPNParticles++;
    if (fDPTrackingTree!= 0) fDPTrackingTree->Fill();

 }// end if track planes is on

}

void LBNEAnalysis::FillTrackingPlaneH1Data

(

const G4Step &

aStep

)

Definition at line 1890 of file LBNEAnalysis.cc.

                                                            {

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetCreateTrkPlaneH1Output()){

    G4Track *track = step.GetTrack();

    
    LBNETrackInformation *info =
    (LBNETrackInformation*)(track->GetUserInformation());
    if(info != 0){
      fH1NImpWt = info->GetNImpWt();
    }
    G4ParticleDefinition *def = track->GetDefinition();
    fH1ParticlePDG = def->GetPDGEncoding();

    G4StepPoint *point = step.GetPostStepPoint();
    G4ThreeVector pos = point->GetPosition();
    G4ThreeVector mom = point->GetMomentum();

    const G4ThreeVector ppoint = track->GetVertexPosition();
    const G4ThreeVector pmom = track->GetVertexMomentumDirection();
    fH1PProductionDX = pmom.x();
    fH1PProductionDY = pmom.y();
    fH1PProductionDZ = pmom.z();
    fH1PProductionX = ppoint.x()/CLHEP::cm;
    fH1PProductionY = ppoint.y()/CLHEP::cm;
    fH1PProductionZ = ppoint.z()/CLHEP::cm;
    fH1ParticlePX = mom.x()/CLHEP::GeV;
    fH1ParticlePY = mom.y()/CLHEP::GeV;
    fH1ParticlePZ = mom.z()/CLHEP::GeV;
    fH1ParticlePXPZ = mom[0]/mom[2];
    fH1ParticlePYPZ = mom[1]/mom[2]; 
    G4ThreeVector dir = point->GetMomentumDirection();
    fH1ParticleX = pos.x()/CLHEP::cm;
    fH1ParticleY = pos.y()/CLHEP::cm;
    fH1ParticleZ = pos.z()/CLHEP::cm;
    fH1ParticleDX = dir.x();
    fH1ParticleDY = dir.y();
    fH1ParticleDZ = dir.z();
    fH1ParticleMass = def->GetPDGMass();
    fH1ParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
    fH1TrackID = track->GetTrackID();
    
       
    
    fH1NParticles++;
    if (fHorn1TrackingTree != 0) fHorn1TrackingTree->Fill();
 // }
  // std::cout << "I stored muon data" << std::endl;

 }// end if track planes is on

}

void LBNEAnalysis::FillTrackingPlaneH2Data

(

const G4Step &

aStep

)

Definition at line 1947 of file LBNEAnalysis.cc.

                                                            {

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

 if (pRunManager->GetCreateTrkPlaneH2Output()){

//  if(fH2NParticles < HkMaxP){
    G4Track *track = step.GetTrack();
    LBNETrackInformation *info =
    (LBNETrackInformation*)(track->GetUserInformation());
    if(info != 0){
      fH2NImpWt = info->GetNImpWt();
    }
    
        if(track->GetCreatorProcess()!=0){
    
    fH2CProcess = track->GetCreatorProcess()->GetProcessName();
    if( fH2CProcess.contains("alphaInelastic"))fH2CNProcess = 1;
    else if(fH2CProcess.contains("AntiLambdaInelastic"))fH2CNProcess = 2;
    else if(fH2CProcess.contains("AntiNeutronInelastic"))fH2CNProcess = 3;
    else if(fH2CProcess.contains("AntiProtonInelastic"))fH2CNProcess = 4;
    else if(fH2CProcess.contains("AntiSigmaMinusInelastic"))fH2CNProcess = 5;
    else if(fH2CProcess.contains("AntiSigmaPlusInelastic"))fH2CNProcess = 6;
    else if(fH2CProcess.contains("AntiXiMinusInelastic"))fH2CNProcess = 7;
    else if(fH2CProcess.contains("Decay"))fH2CNProcess = 8;
    else if(fH2CProcess.contains("hadElastic"))fH2CNProcess = 9;
    else if(fH2CProcess.contains("KaonMinusInelastic"))fH2CNProcess = 10;
    else if(fH2CProcess.contains("KaonZeroLInelastic"))fH2CNProcess = 11;
    else if(fH2CProcess.contains("KaonZeroSInelastic"))fH2CNProcess = 12;
    else if(fH2CProcess.contains("LambdaInelastic"))fH2CNProcess = 13;
    else if(fH2CProcess.contains("NeutronInelastic"))fH2CNProcess = 14;
    else if(fH2CProcess.contains("PionMinusInelastic"))fH2CNProcess = 15;
    else if(fH2CProcess.contains("PionPlusInelastic"))fH2CNProcess = 16;
    else if(fH2CProcess.contains("ProtonInelastic"))fH2CNProcess = 17;
    else if(fH2CProcess.contains("SigmaMinusInelastic"))fH2CNProcess = 18;
    else if(fH2CProcess.contains("SigmaPlusInelastic"))fH2CNProcess = 19;
    else if(fH2CProcess.contains("XiMinusInelastic"))fH2CNProcess = 20;
    else if(fH2CProcess.contains("XiZeroInelastic"))fH2CNProcess = 21;
    else if(fH2CProcess.contains("tInelastic"))fH2CNProcess = 22;

    else fH2CNProcess = 25;
    
    }
    else{
    fH2CProcess = "NA";
    fH2CNProcess = 0;
    }
    
    
    const G4ThreeVector ppoint = track->GetVertexPosition();
     const G4ThreeVector pmom = track->GetVertexMomentumDirection();
    fH2PProductionDX = pmom.x();
    fH2PProductionDY = pmom.y();
    fH2PProductionDZ = pmom.z();
    G4StepPoint *point = step.GetPostStepPoint();
    G4ThreeVector pos = point->GetPosition();
    G4ThreeVector dir = point->GetMomentumDirection();
    G4ParticleDefinition *def = track->GetDefinition();
    G4ThreeVector mom = point->GetMomentum();
    fH2PProductionX = ppoint.x()/CLHEP::cm;
    fH2PProductionY = ppoint.y()/CLHEP::cm;
    fH2PProductionZ = ppoint.z()/CLHEP::cm;
    fH2ParticlePX = mom.x()/CLHEP::GeV;
    fH2ParticlePY = mom.y()/CLHEP::GeV;
    fH2ParticlePZ = mom.z()/CLHEP::GeV;
    fH2ParticlePXPZ = mom[0]/mom[2];
    fH2ParticlePYPZ = mom[1]/mom[2]; 
    fH2ParticleX = pos.x()/CLHEP::cm;
    fH2ParticleY = pos.y()/CLHEP::cm;
    fH2ParticleZ = pos.z()/CLHEP::cm;
    fH2ParticleDX = dir.x();
    fH2ParticleDY = dir.y();
    fH2ParticleDZ = dir.z();
    fH2ParticlePDG = def->GetPDGEncoding();
    fH2ParticleMass = def->GetPDGMass();
    fH2ParticleEnergy = point->GetTotalEnergy()/CLHEP::GeV;
    fH2TrackID = track->GetTrackID();
    
  //      std::cout<<fH2CProcess<<"   "<<fH2ParticlePDG<<"  "<<ppoint<<"   "<<mom<<std::endl;

    
    fH2NParticles++;
    if (fHorn2TrackingTree != 0) fHorn2TrackingTree->Fill();

  //}

 }// end if track planes is on

}

G4int LBNEAnalysis::GetCount

(

)

Definition at line 617 of file LBNEAnalysis.cc.

{
  return fcount;
}

G4double LBNEAnalysis::GetDistanceInVolume	(	LBNETrajectory *	wanted_traj,
		G4String	wanted_vol
	)

Definition at line 2230 of file LBNEAnalysis.cc.

                                                                                          {
  double dist_vol = 0;
  if(wanted_traj==0)return -1.;
  
   //   G4double fact_Al = (26.98/2.7);
//      G4double fact_steel316 = (52.73/8.0);
//      G4double fact_concrete = (33.63/2.03);
//      G4double fact_water = (18.01/1.0);
//      G4double fact_iron = (207.19/11.35);
//      G4double fact_helium = (4.003/0.000145);//This from G4numi
  
  G4ThreeVector ParticlePos;
  G4int npoints = wanted_traj->GetPointEntries();
  
  G4ThreeVector tmp_ipos,tmp_fpos;
  G4ThreeVector tmp_disp;
  
  G4double tmp_dist = 0.0;
  G4bool enter_vol = false;
  G4bool exit_vol  = false;
  for(G4int ii=0; ii<npoints; ++ii){ 
    ParticlePos = (wanted_traj->GetPoint(ii)->GetPosition()/CLHEP::m)*CLHEP::m;
    G4String postvol = "";
    G4String prevol  = wanted_traj->GetPreStepVolumeName(ii);
    if(ii<npoints-1) postvol = wanted_traj->GetPreStepVolumeName(ii+1);

    G4bool vol_in  = (!(prevol.contains(wanted_vol)) && (postvol.contains(wanted_vol)) ) || ( ii==0 && prevol.contains(wanted_vol));
    G4bool vol_out = (prevol.contains(wanted_vol)) && (!postvol.contains(wanted_vol));
    G4bool vol_through = (prevol.contains(wanted_vol))&&(postvol.contains(wanted_vol));
    //Get rid of the vol through and instead create an array of volume of interests and create new ntuples.
   // if(prevol.contains("Water"))vol_through = false; //We dont want water to be part of Inner Conductor.
   // if(vol_in) std::cout<<"Particle entered at "<<prevol<<" and exited at "<<postvol<<std::endl;
   //   if((vol_through)&&(prevol != "DecayPipeVolume_P")) std::cout<<"Particle through the volume "<<prevol<<" and "<<postvol<<std::endl;
    if(vol_in){
      enter_vol = true;
      exit_vol  = false;
      tmp_ipos = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
      tmp_dist = 0.0;
    }
    if(enter_vol &&vol_through&& !exit_vol){
      tmp_fpos = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
      tmp_disp = tmp_fpos - tmp_ipos;
      tmp_dist += tmp_disp.mag();
      tmp_ipos = tmp_fpos;
    }
    if(enter_vol && vol_out){
      tmp_fpos  = G4ThreeVector(ParticlePos[0]/CLHEP::cm,ParticlePos[1]/CLHEP::cm,ParticlePos[2]/CLHEP::cm);
      tmp_disp  = tmp_fpos - tmp_ipos;
      tmp_dist += tmp_disp.mag();
      tmp_ipos  = tmp_fpos;
      exit_vol  = true;
      enter_vol = false;
      dist_vol += tmp_dist;
    }

  }
        

  return dist_vol;
  
}

G4int LBNEAnalysis::GetEntry

(

)

Definition at line 627 of file LBNEAnalysis.cc.

{
  return fentry;
}

LBNEAnalysis * LBNEAnalysis::getInstance ( )

static

Definition at line 134 of file LBNEAnalysis.cc.

{
  if (instance == 0) instance = new LBNEAnalysis;
  return instance;
}

LBNETrajectory * LBNEAnalysis::GetParentTrajectory

(

G4int

parentID

)

Definition at line 1736 of file LBNEAnalysis.cc.

{

  LBNERunManager *pRunManager=
    dynamic_cast<LBNERunManager*>(G4RunManager::GetRunManager());

  if (pRunManager->GetVerboseLevel() == 10) {
   G4cout << "LBNEAnalysis::GetParentTrajectory() called." << G4endl;
  }
      
  G4TrajectoryContainer* container = 
    G4RunManager::GetRunManager()->GetCurrentEvent()->GetTrajectoryContainer();
  if(container==0) return 0;

  TrajectoryVector* vect = container->GetVector();
  G4VTrajectory* tr;
  G4int ii = 0; 
  while (ii<G4int(vect->size())){  
    tr = (*vect)[ii]; 
    LBNETrajectory* tr1 = (LBNETrajectory*)(tr);  
    if(tr1->GetTrackID() == parentID) return tr1; 
    ++ii; 
  }

  return 0;
}

LBNETrajectory * LBNEAnalysis::GetTrajectory

(

G4int

trackID

)

Definition at line 1764 of file LBNEAnalysis.cc.

{
   G4TrajectoryContainer* container = G4RunManager::GetRunManager()->GetCurrentEvent()->GetTrajectoryContainer();
   
   if(container==0) 
   {
      G4cout << "LBNEAnalysis::GetTrajectory - PROBLEM: No Trajectory Container for track ID = " << trackID << endl;
      return 0;
   }   

   TrajectoryVector* vect = container->GetVector();
   G4VTrajectory* tr;
   G4int ii = 0; 
   while (ii < G4int(vect->size()))
   {  
      tr = (*vect)[ii]; 
      LBNETrajectory* tr1 = (LBNETrajectory*)(tr);  
      if(tr1->GetTrackID() == trackID) return tr1; 
      ++ii; 
   }

   G4cout << "LBNEAnalysis::GetTrajectory - PROBLEM: Failed to find track with ID = " << trackID << endl;
   
   return 0;
}

void LBNEAnalysis::ResetEvent

(

)

Definition at line 2292 of file LBNEAnalysis.cc.

{
  fH1NParticles = 0;
  fH2NParticles = 0;
  fNParticles = 0;
  fDPNParticles = 0;
  fTNParticles = 0;
  fMuNParticles = 0;
  //  std::cout << "Resetting event" << std::endl;
}

void LBNEAnalysis::SetCount

(

G4int

count

)

Definition at line 612 of file LBNEAnalysis.cc.

{
  fcount = count;
}

void LBNEAnalysis::setDetectorPositions ( )

private

Definition at line 1789 of file LBNEAnalysis.cc.

                                        {  

  // Info passed along in the Ntuple files, old or new (Dk2Nu), so we place these data statements in one place. 
  
  const int nNear=5;
  const int nFar=3;
  G4double xdetNear[]    = { 0.       ,  0.      ,    11.50,     0.         ,     25.84   };
  G4double ydetNear[]    = { 0.       ,  0.      ,    -3.08,     0.         ,     78.42   };
  G4double zdetNear[]    = {574.     ,1036.49   ,   1000.97,   1030.99     ,    745.25   };
  fDetNameNear.clear();
  fDetNameNear.push_back(std::string("LBNE"));
//  fDetNameNear.push_back(std::string("Minos"));
//  fDetNameNear.push_back(std::string("Nova"));
//  fDetNameNear.push_back(std::string("Minerva"));
//  fDetNameNear.push_back(std::string("MiniBooNE"));
  G4double xdetFar[]     = { 0.       , 0.        , 11040.   };
  G4double ydetFar[]     = { 0.       , 0.        , -4200.   };
  G4double zdetFar[]     = { 1297000. , 735340.  , 810450.   };
  fDetNameFar.clear();
  fDetNameFar.push_back(std::string("LBNE"));
//  fDetNameFar.push_back(std::string("Minos"));
//  fDetNameFar.push_back(std::string("Nova"));
  
//  get(xdetNear[0], "FluxAreaX0near");
//  get(ydetNear[0], "FluxAreaY0near");
//  get(zdetNear[0], "FluxAreaZ0near");
  
//  get(xdetFar[0], "FluxAreaX0far");
//  get(ydetFar[0], "FluxAreaY0far");
//  get(zdetFar[0], "FluxAreaZ0far");
   fXdet_near.clear();
   fYdet_near.clear();
   fZdet_near.clear();
//   
   fXdet_near.resize(nNear);
   fYdet_near.resize(nNear);
   fZdet_near.resize(nNear);

   for(G4int ii=0;ii<nNear;ii++){
     fXdet_near[ii] = xdetNear[ii]*CLHEP::m;
     fYdet_near[ii] = ydetNear[ii]*CLHEP::m;
     fZdet_near[ii] = zdetNear[ii]*CLHEP::m;
   }
   //Far Detector
   // Neutrino data at different points
   // need neutrino parent info to be filled in fLBNEOutNtpData by this point
   // 
   fXdet_far.clear();
   fYdet_far.clear();
   fZdet_far.clear();
   fXdet_far.resize(nFar);
   fYdet_far.resize(nFar);
   fZdet_far.resize(nFar);
   
  for(G4int ii=0;ii<nFar;ii++){
    fXdet_far[ii] = xdetFar[ii]*CLHEP::m;
    fYdet_far[ii] = ydetFar[ii]*CLHEP::m;
    fZdet_far[ii] = zdetFar[ii]*CLHEP::m;
  }
   

}

void LBNEAnalysis::SetEntry

(

G4int

entry

)

Definition at line 622 of file LBNEAnalysis.cc.

{
  fentry = entry;
}

void LBNEAnalysis::TrackThroughGeometry

(

const LBNETrajectory *

TrackTrajectory

)

Definition at line 1583 of file LBNEAnalysis.cc.

{
   G4int trackID = TrackTrajectory-> GetTrackID();

   G4int numberOfPoints            = TrackTrajectory->GetPointEntries();
   
   for (G4int ii=numberOfPoints-1; ii > -1; --ii)
   {
      
      G4String prevolname = TrackTrajectory->GetPreStepVolumeName(ii);
      G4String postvolname = "";
      if(ii < numberOfPoints-1) postvolname = TrackTrajectory->GetPreStepVolumeName(ii+1);
      
      G4ThreeVector ParticleMomentum = TrackTrajectory->GetMomentum(ii);              
      G4ThreeVector ParticlePosition = TrackTrajectory->GetPoint(ii)->GetPosition();
      
      if (((numberOfPoints - ii) < 3) || (ii < 2))
      {
      //++count;
//      G4double steplength = TrackTrajectory->GetStepLength(ii);              
      //if(count > 1)
//      std::cout << "eventid: " << fLBNEOutNtpData->evtno << " trackid: " << trackID << " ii = " << ii 
//      << " particle = " << TrackTrajectory->GetParticleName()
//      << " prestepvolname = " << prevolname << " postvolname = " << postvolname 
//      << " steplength = " << steplength/mm << " mm" << std::endl;
      }
    
      TrackPoint_t TrkPt;
      TrkPt.type = LBNEParticleCode::AsInt(LBNEParticleCode::StringToEnum(TrackTrajectory->GetParticleName()));
      TrkPt.x     = ParticlePosition[0]/CLHEP::cm;
      TrkPt.y     = ParticlePosition[1]/CLHEP::cm;
      TrkPt.z     = ParticlePosition[2]/CLHEP::cm;
      TrkPt.px    = ParticleMomentum[0]/CLHEP::GeV;
      TrkPt.py    = ParticleMomentum[1]/CLHEP::GeV;
      TrkPt.pz    = ParticleMomentum[2]/CLHEP::GeV;
      TrkPt.trkid = trackID;
      TrkPt.impwt = TrackTrajectory->GetNImpWt();
      TrkPt.gen   = (TrackTrajectory->GetTgen())-1;
      
      //
      //particle created in the target
      //
      if (prevolname.contains("TargetFin") && ii==0)
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Target"), TrkPt);
      }
      //
      //particle as exits target
      //
      if ((prevolname.contains("Target")) && 
           postvolname.contains("Horn1") && (!postvolname.contains("Target")))
      {
         /*if(fLBNEOutNtpData->evtno == 1450 || fLBNEOutNtpData->evtno == 1679)
           {
           std::cout << "EVENT ID: " << fLBNEOutNtpData->evtno << std::endl;
           std::cout << "   LBNEAnalysis particle exiting target..." << std::endl;
           std::cout << "      track id = " <<  TrkPt.trkid << std::endl;
           std::cout << "      type     = " <<  TrkPt.type << std::endl;
           std::cout << "      x pos    = " <<  TrkPt.x << std::endl;
           std::cout << "      y pos    = " <<  TrkPt.y << std::endl;
           std::cout << "      z pos    = " <<  TrkPt.z << std::endl;
           std::cout << "      px mom   = " <<  TrkPt.px << std::endl;
           std::cout << "      py mom   = " <<  TrkPt.py << std::endl;
           std::cout << "      pz mom   = " <<  TrkPt.pz << std::endl;
           std::cout << "   LBNEAnalysis Nu Parent found exiting target..." << std::endl;
           std::cout << "      track id = " <<  tptrkid << std::endl;
           std::cout << "      type     = " <<  fLBNEOutNtpData->tptype << std::endl;
           std::cout << "      x pos    = " <<  fLBNEOutNtpData->tvx << std::endl;
           std::cout << "      y pos    = " <<  fLBNEOutNtpData->tvy << std::endl;
           std::cout << "      z pos    = " <<  fLBNEOutNtpData->tvz << std::endl;
           std::cout << "      px mom   = " <<  fLBNEOutNtpData->tpx << std::endl;
           std::cout << "      py mom   = " <<  fLBNEOutNtpData->tpy << std::endl;
           std::cout << "      pz mom   = " <<  fLBNEOutNtpData->tpz << std::endl;
           }*/
         
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("TargetExit"), TrkPt);
      }
      //
      //particle at plane of end of tgt. This only makes sense 
      //if only running with the target and target hall constructed.
      //otherwise this will
      //
      if (prevolname.contains("TgtEndPlane"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("TargetEndPlane"), TrkPt);
      }


//      if(LBNEData->GetSimulation() == "Target Tracking") continue;
      //
      // The following volume names need updating.....
     // 

      //particle enters horn 1
      if ((prevolname.contains("Horn1TopLevelDownstr")) && postvolname.contains("Horn1TopLevelDownstr"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1Enter"), TrkPt);
      }
      //particle at neck of horn 1 the last entry will be just before it leaves the neck
      //if (prevolname.contains("PH01-02"))
      if (prevolname.contains("Horn1DownstrPart1"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1NeckPlane"), TrkPt);
      }
      //particle exists horn 1
      if ((prevolname.contains("Horn1")) && postvolname.contains("Tunnel"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1Exit"), TrkPt);
      }
      //particle at end plane of horn 1
      if (prevolname.contains("PH01EndPlane"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn1EndPlane"), TrkPt);
      }
      //particle enters horn 2
      if ((prevolname.contains("Tunnel")) && postvolname.contains("Horn2"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2Enter"), TrkPt);
      }
      //particle at neck of horn 2 the last entry will be just before it leaves the neck
      if (prevolname.contains("Horn2Part2"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2NeckPlane"), TrkPt);
      }
      //particle exits horn 2
      if ((prevolname.contains("Horn2")) && postvolname.contains("Tunnel"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2Exit"), TrkPt);
      }
      //particle at end plane of horn 2
      if (prevolname.contains("PH02EndPlane"))
      {
         fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("Horn2EndPlane"), TrkPt);
      }
      //Particle enters DP
      if ( (prevolname.contains("DecayPipeUsptrWindow")) &&  
           (postvolname.contains("DecayPipeVolume"))) {
           
           fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("DPipeEnter"), TrkPt);
      }
      //Particle exits DP
      if ( (prevolname.contains("DecayPipe"))  && (postvolname.contains("Tunnel")))
         
      {
       fLBNEOutNtpData -> AddTrkPoint(TrkPoint::StringToEnum("DPipeExit"), TrkPt);
      }
      
   }//end loop over npts in Track
   
}

Member Data Documentation

std::map<int, int> LBNEAnalysis::code

private

Definition at line 255 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::dist_DPIP[3]

private

Definition at line 299 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::dist_DVOL[3]

private

Definition at line 300 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::dist_IC1[3]

private

Definition at line 297 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::dist_IC2[3]

private

Definition at line 298 of file LBNEAnalysis.hh.

int LBNEAnalysis::dk_trackID

private

Definition at line 294 of file LBNEAnalysis.hh.

bool LBNEAnalysis::doPPFx

private

Definition at line 305 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fAlcoveTrackingTree

private

Definition at line 246 of file LBNEAnalysis.hh.

G4int LBNEAnalysis::fcount

private

Definition at line 278 of file LBNEAnalysis.hh.

std::vector<G4String> LBNEAnalysis::fDetNameFar

private

Definition at line 290 of file LBNEAnalysis.hh.

std::vector<G4String> LBNEAnalysis::fDetNameNear

private

Definition at line 289 of file LBNEAnalysis.hh.

bsim::Dk2Nu* LBNEAnalysis::fDk2Nu

private

Definition at line 270 of file LBNEAnalysis.hh.

bool LBNEAnalysis::fDk2NuDetectorFileRead

private

Definition at line 82 of file LBNEAnalysis.hh.

bsim::DkMeta* LBNEAnalysis::fDkMeta

private

Definition at line 271 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPNImpWt

private

Definition at line 165 of file LBNEAnalysis.hh.

int LBNEAnalysis::fDPNParticles

private

Definition at line 185 of file LBNEAnalysis.hh.

int LBNEAnalysis::fDPParentID

private

Definition at line 163 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleDX

private

Definition at line 171 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleDY

private

Definition at line 172 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleDZ

private

Definition at line 173 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleEnergy

private

Definition at line 169 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleMass

private

Definition at line 170 of file LBNEAnalysis.hh.

int LBNEAnalysis::fDPParticlePDG

private

Definition at line 164 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticlePX

private

Definition at line 174 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticlePXPZ

private

Definition at line 177 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticlePY

private

Definition at line 175 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticlePYPZ

private

Definition at line 178 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticlePZ

private

Definition at line 176 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleX

private

Definition at line 166 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleY

private

Definition at line 167 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPParticleZ

private

Definition at line 168 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionDX

private

Definition at line 182 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionDY

private

Definition at line 183 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionDZ

private

Definition at line 184 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionX

private

Definition at line 179 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionY

private

Definition at line 180 of file LBNEAnalysis.hh.

double LBNEAnalysis::fDPPProductionZ

private

Definition at line 181 of file LBNEAnalysis.hh.

int LBNEAnalysis::fDPTrackID

private

Definition at line 162 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fDPTrackingTree

private

Definition at line 244 of file LBNEAnalysis.hh.

G4int LBNEAnalysis::fentry

private

Definition at line 279 of file LBNEAnalysis.hh.

std::string LBNEAnalysis::fH1CProcess

private

Definition at line 127 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH1DParticlePDG

private

Definition at line 106 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH1DTrackID

private

Definition at line 107 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1NImpWt

private

Definition at line 105 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH1NParticles

private

Definition at line 129 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleDX

private

Definition at line 113 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleDY

private

Definition at line 114 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleDZ

private

Definition at line 115 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleEnergy

private

Definition at line 111 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleMass

private

Definition at line 112 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH1ParticlePDG

private

Definition at line 104 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticlePX

private

Definition at line 116 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticlePXPZ

private

Definition at line 119 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticlePY

private

Definition at line 117 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticlePYPZ

private

Definition at line 120 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticlePZ

private

Definition at line 118 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleX

private

Definition at line 108 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleY

private

Definition at line 109 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1ParticleZ

private

Definition at line 110 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionDX

private

Definition at line 124 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionDY

private

Definition at line 125 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionDZ

private

Definition at line 126 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionX

private

Definition at line 121 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionY

private

Definition at line 122 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH1PProductionZ

private

Definition at line 123 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH1TrackID

private

Definition at line 103 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH2CNProcess

private

Definition at line 157 of file LBNEAnalysis.hh.

G4String LBNEAnalysis::fH2CProcess

private

Definition at line 156 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2NImpWt

private

Definition at line 136 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH2NParticles

private

Definition at line 158 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleDX

private

Definition at line 142 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleDY

private

Definition at line 143 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleDZ

private

Definition at line 144 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleEnergy

private

Definition at line 140 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleMass

private

Definition at line 141 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH2ParticlePDG

private

Definition at line 135 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticlePX

private

Definition at line 145 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticlePXPZ

private

Definition at line 148 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticlePY

private

Definition at line 146 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticlePYPZ

private

Definition at line 149 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticlePZ

private

Definition at line 147 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleX

private

Definition at line 137 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleY

private

Definition at line 138 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2ParticleZ

private

Definition at line 139 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionDX

private

Definition at line 153 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionDY

private

Definition at line 154 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionDZ

private

Definition at line 155 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionX

private

Definition at line 150 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionY

private

Definition at line 151 of file LBNEAnalysis.hh.

double LBNEAnalysis::fH2PProductionZ

private

Definition at line 152 of file LBNEAnalysis.hh.

int LBNEAnalysis::fH2TrackID

private

Definition at line 134 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fHorn1TrackingTree

private

Definition at line 242 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fHorn2TrackingTree

private

Definition at line 243 of file LBNEAnalysis.hh.

LBNEDataNtp_t* LBNEAnalysis::fLBNEOutNtpData

private

Definition at line 273 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuDE[kMaxP]

private

Definition at line 227 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMudEdx[kMaxP]

private

Definition at line 225 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMudEdx_ion[kMaxP]

private

Definition at line 226 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuDEion[kMaxP]

private

Definition at line 228 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuEnergy[kMaxP]

private

Definition at line 220 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuEvtNo

private

Definition at line 205 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuMass[kMaxP]

private

Definition at line 219 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuNimpWt[kMaxP]

private

Definition at line 210 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuNParticles

private

Definition at line 206 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuNSteps[kMaxP]

private

Definition at line 230 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParE[kMaxP]

private

Definition at line 235 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuParentID[kMaxP]

private

Definition at line 208 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuParPDG[kMaxP]

private

Definition at line 231 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParPX[kMaxP]

private

Definition at line 236 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParPY[kMaxP]

private

Definition at line 237 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParPZ[kMaxP]

private

Definition at line 238 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParX[kMaxP]

private

Definition at line 232 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParY[kMaxP]

private

Definition at line 233 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuParZ[kMaxP]

private

Definition at line 234 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuPDG[kMaxP]

private

Definition at line 209 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuPX[kMaxP]

private

Definition at line 221 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuPY[kMaxP]

private

Definition at line 222 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuPZ[kMaxP]

private

Definition at line 223 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuRunNo

private

Definition at line 204 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuStartE[kMaxP]

private

Definition at line 218 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuStartX[kMaxP]

private

Definition at line 215 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuStartY[kMaxP]

private

Definition at line 216 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuStartZ[kMaxP]

private

Definition at line 217 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuT[kMaxP]

private

Definition at line 214 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuTheta[kMaxP]

private

Definition at line 224 of file LBNEAnalysis.hh.

int LBNEAnalysis::fMuTrackID[kMaxP]

private

Definition at line 207 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuX[kMaxP]

private

Definition at line 211 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuY[kMaxP]

private

Definition at line 212 of file LBNEAnalysis.hh.

double LBNEAnalysis::fMuZ[kMaxP]

private

Definition at line 213 of file LBNEAnalysis.hh.

double LBNEAnalysis::fNImpWt[kMaxP]

private

Definition at line 92 of file LBNEAnalysis.hh.

int LBNEAnalysis::fNParticles

private

Definition at line 101 of file LBNEAnalysis.hh.

std::ofstream LBNEAnalysis::fOutDBGDk2nu_

private

Definition at line 306 of file LBNEAnalysis.hh.

TFile* LBNEAnalysis::fOutFile

private

Definition at line 258 of file LBNEAnalysis.hh.

TFile* LBNEAnalysis::fOutFileDk2Nu

private

Definition at line 263 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fOutTree

private

Definition at line 259 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fOutTreeDk2Nu

private

Definition at line 264 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fOutTreeDk2NuMeta

private

Definition at line 265 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleDX[kMaxP]

private

Definition at line 98 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleDY[kMaxP]

private

Definition at line 99 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleDZ[kMaxP]

private

Definition at line 100 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleEnergy[kMaxP]

private

Definition at line 96 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleMass[kMaxP]

private

Definition at line 97 of file LBNEAnalysis.hh.

int LBNEAnalysis::fParticlePDG[kMaxP]

private

Definition at line 91 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleX[kMaxP]

private

Definition at line 93 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleY[kMaxP]

private

Definition at line 94 of file LBNEAnalysis.hh.

double LBNEAnalysis::fParticleZ[kMaxP]

private

Definition at line 95 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fTargetOutputTree

private

Definition at line 245 of file LBNEAnalysis.hh.

int LBNEAnalysis::fTNParticles

private

Definition at line 198 of file LBNEAnalysis.hh.

int LBNEAnalysis::fTParentID

private

Definition at line 189 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticleEnergy

private

Definition at line 194 of file LBNEAnalysis.hh.

int LBNEAnalysis::fTParticlePDG

private

Definition at line 193 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticlePX

private

Definition at line 195 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticlePY

private

Definition at line 196 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticlePZ

private

Definition at line 197 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticleX

private

Definition at line 190 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticleY

private

Definition at line 191 of file LBNEAnalysis.hh.

double LBNEAnalysis::fTParticleZ

private

Definition at line 192 of file LBNEAnalysis.hh.

int LBNEAnalysis::fTrackID[kMaxP]

private

Definition at line 90 of file LBNEAnalysis.hh.

LBNEDataNtp_t* LBNEAnalysis::fTrackingPlaneData

private

Definition at line 274 of file LBNEAnalysis.hh.

TTree* LBNEAnalysis::fTrackingTree

private

Definition at line 267 of file LBNEAnalysis.hh.

int LBNEAnalysis::fTTrackID

private

Definition at line 188 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fXdet_far

private

Definition at line 286 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fXdet_near

private

Definition at line 283 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fYdet_far

private

Definition at line 287 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fYdet_near

private

Definition at line 284 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fZdet_far

private

Definition at line 288 of file LBNEAnalysis.hh.

std::vector<G4double> LBNEAnalysis::fZdet_near

private

Definition at line 285 of file LBNEAnalysis.hh.

vstring_t LBNEAnalysis::GenAbsName

private

Definition at line 304 of file LBNEAnalysis.hh.

LBNEAnalysis * LBNEAnalysis::instance = 0

staticprivate

Definition at line 80 of file LBNEAnalysis.hh.

G4int LBNEAnalysis::nGenAbs = 3

private

Definition at line 280 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::noProtons

private

Definition at line 200 of file LBNEAnalysis.hh.

int LBNEAnalysis::numEntFillNTu_

private

Definition at line 308 of file LBNEAnalysis.hh.

int LBNEAnalysis::numEntFillNTuDk2nuPzR_

private

Definition at line 310 of file LBNEAnalysis.hh.

int LBNEAnalysis::numEntFillNTuPzR_

private

Definition at line 309 of file LBNEAnalysis.hh.

G4ThreeVector LBNEAnalysis::NuMomentum

private

Definition at line 128 of file LBNEAnalysis.hh.

TFile* LBNEAnalysis::nuNtuple

private

Definition at line 261 of file LBNEAnalysis.hh.

char LBNEAnalysis::nuNtupleFileName[1024]

private

Definition at line 250 of file LBNEAnalysis.hh.

char LBNEAnalysis::nuNtupleFileNameDK2Nu[1024]

private

Definition at line 251 of file LBNEAnalysis.hh.

int LBNEAnalysis::nVdblTot

private

Definition at line 296 of file LBNEAnalysis.hh.

int LBNEAnalysis::nVintTot

private

Definition at line 295 of file LBNEAnalysis.hh.

G4int LBNEAnalysis::nVolAbs = 4

private

Definition at line 281 of file LBNEAnalysis.hh.

int LBNEAnalysis::tar_trackID =-1

private

Definition at line 293 of file LBNEAnalysis.hh.

std::vector<bsim::Traj> LBNEAnalysis::vec_traj

private

Definition at line 292 of file LBNEAnalysis.hh.

vstring_t LBNEAnalysis::VolAbsName

private

Definition at line 302 of file LBNEAnalysis.hh.

vstring_t LBNEAnalysis::VolVdblName

private

Definition at line 301 of file LBNEAnalysis.hh.

vstring_t LBNEAnalysis::VolVintName

private

Definition at line 303 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::x

private

Definition at line 84 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::y

private

Definition at line 85 of file LBNEAnalysis.hh.

G4double LBNEAnalysis::z

private

Definition at line 86 of file LBNEAnalysis.hh.

---

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

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