protoana::truepionXsection Class Reference

Inheritance diagram for protoana::truepionXsection:

Public Member Functions
	truepionXsection (fhicl::ParameterSet const &pset)
virtual	~truepionXsection ()
void	beginJob ()
void	endJob ()
void	beginRun (const art::Run &run)
void	analyze (const art::Event &evt)
void	reset ()
Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDAnalyzer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept
	Analyzer (fhicl::ParameterSet const &pset)
template<typename Config >
	Analyzer (Table< Config > const &config)
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Observer
	~Observer () noexcept
	Observer (Observer const &)=delete
	Observer (Observer &&)=delete
Observer &	operator= (Observer const &)=delete
Observer &	operator= (Observer &&)=delete
void	registerProducts (ProductDescriptions &, ModuleDescription const &)
void	fillDescriptions (ModuleDescription const &)
fhicl::ParameterSetID	selectorConfig () const
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
ProtoDUNETruthUtils	truthUtil
	below are some addition to be able to use protodune utitilites More...
const art::InputTag	fBeamModuleLabel
std::string	fCalorimetryTag
std::string	fTrackerTag
std::string	fShowerTag
std::string	fPFParticleTag
std::string	fGeneratorTag
TTree *	fEventTree
Int_t	run
Int_t	subrun
Int_t	event
TH1D *	h_DEUniform
TH1D *	h_DXUniform
TH1D *	h_DEDXUniform
TH1D *	h_DeltaE
TH1D *	h_UniformDistances
TH1D *	hInteractingKE
TH1D *	hInteractingKEEl
TH1D *	hInteractingKEElDep
TH1D *	hInteractingKEInel
TH1D *	hIncidentKE
TH1D *	hCrossSection
TH1D *	hCrossSectionEl
TH1D *	hCrossSectionInel
TH1D *	hKEAtTPCFF
TH1D *	hInitialKE
TH1D *	hInitialPz
TH2D *	hXZ
TH2D *	hYZ
TH2D *	hXZPre
TH2D *	hYZPre
TH2D *	hdEVsdX
TH2D *	hdEVsKE
TH1D *	processtot
TH1D *	processEl
TH1D *	processInel
double	trueVtxX
double	trueVtxY
double	trueVtxZ
double	trueEndX
double	trueEndY
double	trueEndZ
double	finalKE
std::vector< std::string >	G4Process
double	minX = -360.0
double	maxX = 360.0
double	minY =0.0
double	maxY = 600.0
double	minZ = 0.0
double	maxZ = 695.0
int	evtsPreTPC = 0
int	evtsInTheMiddle = 0
int	evtsPostTPC = 0
int	throughgoing = 0
int	interactingInTPC = 0
std::string	fHitsModuleLabel
std::string	fTrackModuleLabel
std::string	fCalorimetryModuleLabel
bool	fSaveCaloInfo
bool	fSaveTrackInfo

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >
using	ModuleType = EDAnalyzer
Protected Member Functions inherited from art::Observer
std::string const &	processName () const
bool	wantAllEvents () const noexcept
bool	wantEvent (ScheduleID id, Event const &e) const
Handle< TriggerResults >	getTriggerResults (Event const &e) const
	Observer (fhicl::ParameterSet const &config)
	Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 96 of file truepionXsection_module.cc.

Constructor & Destructor Documentation

protoana::truepionXsection::truepionXsection ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 190 of file truepionXsection_module.cc.

                                                                  :
    EDAnalyzer(pset),
    ///new lines 
    //dataUtil(pset.get<fhicl::ParameterSet>("DataUtils")),
    fBeamModuleLabel(pset.get< art::InputTag >("BeamModuleLabel")),
    //fTrackerTag(pset.get<std::string>("TrackerTag")),
    //fShowerTag(pset.get<std::string>("ShowerTag")),
    //fPFParticleTag(pset.get<std::string>("PFParticleTag")),
    fGeneratorTag(pset.get<std::string>("GeneratorTag")),
    ///new lines end here

    fHitsModuleLabel          (pset.get< std::string >("HitsModuleLabel","")         ),
    fTrackModuleLabel         (pset.get< std::string >("TrackModuleLabel","")        ),
    fCalorimetryModuleLabel   (pset.get< std::string >("CalorimetryModuleLabel","")  ),
    fSaveCaloInfo             (pset.get< bool>("SaveCaloInfo",false)),
    fSaveTrackInfo            (pset.get< bool>("SaveTrackInfo",false))
  {
    if (fSaveTrackInfo == false) fSaveCaloInfo = false;
  }

protoana::truepionXsection::~truepionXsection ( )

virtual

Definition at line 211 of file truepionXsection_module.cc.

                                     {
  }

Member Function Documentation

void protoana::truepionXsection::analyze

(

const art::Event &

evt

)

Definition at line 289 of file truepionXsection_module.cc.

                                                    {
    reset();  
    bool verbose=false;
    run = evt.run();
    subrun = evt.subRun();
    event = evt.id().event();
    art::ServiceHandle<cheat::BackTrackerService> bt_serv;
    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
    const sim::ParticleList& plist=pi_serv->ParticleList();
    // sim::ParticleList::const_iterator itPartp=plist.begin();
    bool keepInteraction=false;

    //for(size_t iPart=0;(iPart<plist.size()) && (plist.begin()!=plist.end());++iPart){  //particle loop begins here
    auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
    const simb::MCParticle* particle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);

    //const simb::MCParticle* particle=(itPartp++)->second;
    if(particle!=0x0){
      std::cout<<"particle pdg code is "<<particle->PdgCode()<<std::endl;
      if(particle->PdgCode()==211 && particle->Process()=="primary"){
        std::cout<<"found a pi+"<<std::endl;
        art::ServiceHandle<geo::Geometry> geom;
        simb::MCTrajectory truetraj=particle->Trajectory();
        geo::View_t view = geom->View(2);
        auto simIDE_prim=bt_serv->TrackIdToSimIDEs_Ps(particle->TrackId(),view);
        std::map<double, sim::IDE> orderedSimIDE;
        for (auto& ide : simIDE_prim) orderedSimIDE[ide->z]= *ide;
        std::string interactionLabel="";
        double mass=particle->Mass();
        if(verbose) std::cout<<"mass "<<mass<<"\n";
        //Store the kinetic energy and momentum on z at WC4. Just for cross check 
        auto inTPCPoint  = truetraj.begin(); 
        auto Momentum0   = inTPCPoint->second;
        double KE0 = 1000*(TMath::Sqrt(Momentum0.X()*Momentum0.X() + Momentum0.Y()*Momentum0.Y() + Momentum0.Z()*Momentum0.Z() + mass*mass ) - mass); //I want this in MeV
        hInitialKE->Fill(KE0);
        hInitialPz->Fill(1000*Momentum0.Z());
        if(verbose) std::cout<<KE0;

        //--------------------------------------------------------
        // Identify the first trajectory point inside the TPC
        // Loop From First TrajPoint --> First Point in TPC 
        // Stop when you get into the TPC
        for ( auto t = truetraj.begin(); t!= std::prev(truetraj.end()); t++)
          {
            auto pos = t->first;
            if (pos.Z() < minZ) continue;
            else if (pos.X() < minX || pos.X() > maxX ) continue;
            else if (pos.Y() < minY || pos.Y() > maxY ) continue;
            else {
              inTPCPoint = t;
              break;
            }
          }// End search for first point in TPC
        evtsPreTPC++;
        hXZPre->Fill((inTPCPoint->first).Z(), (inTPCPoint->first).X());
        hYZPre->Fill((inTPCPoint->first).Z(), (inTPCPoint->first).Y());
        if (inTPCPoint !=truetraj.begin()){
        evtsPostTPC++;     
        hXZ   ->Fill((inTPCPoint->first).Z(), (inTPCPoint->first).X());
        hYZ   ->Fill((inTPCPoint->first).Z(), (inTPCPoint->first).Y());

        // Identify the last interesting trajectory point in TPC
        auto finTPCPoint = std::prev(truetraj.end()); 
        // The last point is a bit more complicated:
        // if there's no interaction, then it is simply the last point in the TPC
        // if there's one or more interaction points, it's the first interaction point deemed interesting (no coulomb)
        // Take the interaction Map... check if there's something there
        auto thisTracjectoryProcessMap =  truetraj.TrajectoryProcesses();
    
        int el=0;
        int inel=0;
        int totint=0;
        if (thisTracjectoryProcessMap.size())
          {
            for(auto const& couple: thisTracjectoryProcessMap) 
              { 
                totint++;
                // I'm not interested in the CoulombScattering, discard this case
                if(!verbose) std::cout<<(truetraj.KeyToProcess(couple.second))<<" Position "<<((truetraj.at(couple.first)).first).Z()<<"  Momentum "<<((truetraj.at(couple.first)).second).Z()<<std::endl;
                if ((truetraj.KeyToProcess(couple.second)).find("CoulombScat")!= std::string::npos) continue;
              
                // Let's check if the interaction is in the the TPC
                auto     interactionPos4D =  (truetraj.at(couple.first)).first ;        
                if      (interactionPos4D.Z() <  minZ || interactionPos4D.Z() > maxZ ) continue;
                else if (interactionPos4D.X() <  minX || interactionPos4D.X() > maxX ) continue;
                else if (interactionPos4D.Y() <  minY || interactionPos4D.Y() > maxY ) continue;
                if ((truetraj.KeyToProcess(couple.second)).find("hadElastic")!= std::string::npos) el++;
                if ((truetraj.KeyToProcess(couple.second)).find("pi+InElastic")!= std::string::npos) inel++;

                // If we made it here, then this is the first interesting interaction in the TPC
                // Our job is done!!! Great! Store the interaction label and the iterator for the final point
                interactionLabel = truetraj.KeyToProcess(couple.second);
                std::cout<<"interaction Label "<<interactionLabel<<" EndProcess "<<particle->EndProcess()<<std::endl;
                finTPCPoint = truetraj.begin() + couple.first; 
                keepInteraction=true;
                interactingInTPC++;
                break; //commented for now
              }// Loop on interaction points
          } // If there are G4 interactions 

        std::cout<<"this trjectory size "<<thisTracjectoryProcessMap.size()<<std::endl;
        processtot->Fill(totint);
        processEl->Fill(el);
        processInel->Fill(inel);
  

        sim::ParticleList::const_iterator itPart=plist.begin();
  
        // If I didn't find anything interesting in the intereaction map, let's loop back!
        if (!keepInteraction)
          {
            //Loop on the daughters 
            for(size_t iPart1=0;(iPart1<plist.size()) && (plist.begin()!=plist.end());++iPart1)
              {
                const simb::MCParticle* pPart=(itPart++)->second;
                //   art::Ptr<simb::MCTruth> const& mcDaught=pi_serv->ParticleToMCTruth_P(pPart)
                //We keep only the dauthers of the primary not coming from elastic or inelastic scattering
                if (pPart->Mother()  != 1 ) continue;
                //  std::cout<<"pPart mother "<<pPart->Mother()<<std::endl;
                if ((pPart->Process()).find("astic")!= std::string::npos) continue;
                if ((pPart->Process()).find("CoulombScat")!= std::string::npos) continue;
                //Is the daughter born inside the TPC? If yes, store the process which created it 
                simb::MCTrajectory trueDaugthTraj = pPart->Trajectory();              
                if (trueDaugthTraj.begin()->first.Z() < minZ || trueDaugthTraj.begin()->first.Z() > maxZ) continue;
                else if (trueDaugthTraj.begin()->first.X() <   minX || trueDaugthTraj.begin()->first.X() > maxX ) continue;
                else if (trueDaugthTraj.begin()->first.Y() <   minY || trueDaugthTraj.begin()->first.Y() > maxY ) continue;
                else {
                  interactionLabel = pPart->Process();
            
                  break;
                }
              }        
            for ( auto t = std::prev(truetraj.end()); t!= truetraj.begin(); t--)
              {
                auto pos = t->first;
            
                if (pos.Z() > maxZ) continue;
                else if (pos.X() <   minX || pos.X() > maxX ) continue;
                else if (pos.Y() <   minY || pos.Y() > maxY ) continue;
                else {
                  finTPCPoint = t;
                  break;
                }
              }
          }
        //}//new bracket added here

        if (finTPCPoint != inTPCPoint){
          auto posFin = finTPCPoint->first;
          auto posIni = inTPCPoint->first;
          //Let's record what the initial and final points are.
          trueVtxX = posIni.X();
          trueVtxY = posIni.Y();
          trueVtxZ = posIni.Z();
          trueEndX = posFin.X();
          trueEndY = posFin.Y();
          trueEndZ = posFin.Z();
          // std::cout<<"initial x, y and z && final x, y and z "<<trueVtxX<<" "<<trueVtxY<<" "<<trueVtxZ<<" "<<trueEndX<<" "<<trueEndY<<" "<<trueEndZ<<std::endl;
          auto totLength = distance(posFin.X(), posFin.Y(), posFin.Z(),posIni.X(), posIni.Y(), posIni.Z() );
          // std::cout<<"totalLength "<<totLength<<std::endl;
      
          // We want to chop up the points between the first and last uniformly
          // and ordered by Z
          // Order them in order of increasing Z
          std::map<double, TVector3> orderedUniformTrjPts;
          // We want the first and uniform point to coincide with the 
          // the first and last points we just found 
          auto positionVector0 = (inTPCPoint ->first).Vect(); 
          auto positionVector1 = (finTPCPoint->first).Vect(); 
          orderedUniformTrjPts[positionVector0.Z()] = positionVector0;
          orderedUniformTrjPts[positionVector1.Z()] = positionVector1;
          // const double trackPitch = 0.47;
          const double trackPitch = 0.4792;
          // I do have space for at least one extra point, so let's put it there!
          // Calculate how many extra points I need to put between the new first point and the second TrajPoint
          int    nPts            = (int) (totLength/trackPitch);
          for (int iPt = 1; iPt <= nPts; iPt++ )
            {
              auto newPoint = positionVector0 + iPt*(trackPitch/totLength) * (positionVector1 - positionVector0);
              orderedUniformTrjPts[newPoint.Z()] = newPoint;
            }
          // If the distance between the last point and the second to last is less then 0.235
          // eliminate the second to last point
          auto lastPt         = (orderedUniformTrjPts.rbegin())->second;
          auto secondtoLastPt = (std::next(orderedUniformTrjPts.rbegin()))->second;
          double lastDist = distance(lastPt.X(),lastPt.Y(),lastPt.Z(),secondtoLastPt.X(),secondtoLastPt.Y(),secondtoLastPt.Z());
          if (lastDist < 0.240)
            {
              orderedUniformTrjPts.erase((std::next(orderedUniformTrjPts.rbegin()))->first );
            } 
    
          // Calculate the initial kinetic energy
          auto initialMom =     inTPCPoint->second;
          double initialKE = 1000*(TMath::Sqrt(initialMom.X()*initialMom.X() + initialMom.Y()*initialMom.Y() + initialMom.Z()*initialMom.Z() + mass*mass ) - mass); 
          hKEAtTPCFF->Fill(initialKE);
          double kineticEnergy = initialKE;
          auto old_it = orderedUniformTrjPts.begin();
          for (auto it = std::next(orderedUniformTrjPts.begin()); it != orderedUniformTrjPts.end(); it++, old_it++ )
            {
          
              if (verbose)  std::cout << it->first<<" : " << (it->second).Z() << std::endl ;
              auto oldPos        = old_it->second;
              auto currentPos    =     it->second;
        
              double uniformDist =  (currentPos - oldPos).Mag();
              h_UniformDistances->Fill(uniformDist);
          
              //Calculate the energy deposited in this slice          
              auto old_iter = orderedSimIDE.begin();
              double currentDepEnergy = 0.;
              for ( auto iter= orderedSimIDE.begin(); iter!= orderedSimIDE.end(); iter++,old_iter++)
                {
                  auto currentIde = iter->second;
                  // std::cout<<"Z position of the trajectory hits "<<currentIde.z<<std::endl;
                  if ( currentIde.z < oldPos.Z()) continue;
                  if ( currentIde.z > currentPos.Z()) continue;
                  currentDepEnergy += currentIde.energy;
                }// Determing which simIDE is within the current slice
              // avoid overfilling super tiny energy depositions
              if (currentDepEnergy/uniformDist < 0.1 ) continue;
              //Calculate the current kinetic energy
              kineticEnergy -= currentDepEnergy;
              hdEVsdX->Fill(currentDepEnergy,(currentPos.Z()-oldPos.Z()) );
              hdEVsKE->Fill(currentDepEnergy,kineticEnergy);
              hIncidentKE->Fill(kineticEnergy);
              h_DEUniform->Fill(currentDepEnergy);
              h_DXUniform->Fill(uniformDist);
              h_DEDXUniform->Fill(currentDepEnergy/uniformDist);
          
            }// Loop on OrderedPoints

          if(interactionLabel.find("Inelastic")!= std::string::npos )//added pi+Inelastic instead of Inelastic
            //if(interactionLabel=="pi+Inelastic")
            // if(particle->EndProcess()=="pi+Inelastic")
            {
              // std::cout<<"Interaction Label: "<<interactionLabel<<"  EndProcess "<<particle->EndProcess(); 
              hInteractingKE->Fill(kineticEnergy);
              hInteractingKEInel->Fill(kineticEnergy);
              //std::cout<<"inside pi+Inelastic loop and ke "<<kineticEnergy<<std::endl;
            }
      
          //Fill the Elastic and Total Interacting with the last point
          if ( interactionLabel.find("Elastic")!= std::string::npos ) //added hadElastic instead of Elastic
            //   if(interactionLabel=="hadElastic")
            // if(particle->EndProcess()=="hadElastic")
            {
              //  std::cout<<"Interaction Label: "<<interactionLabel<<"  EndProcess "<<particle->EndProcess(); 
              h_DeltaE ->Fill(kineticEnergy -  1000*((finTPCPoint->second).E() - mass) );
              hInteractingKEElDep->Fill(kineticEnergy);
              hInteractingKE->Fill(kineticEnergy);
              //        auto MomentumF = finTPCPoint->second;
              //        double KEF = 1000*(TMath::Sqrt(MomentumF.X()*MomentumF.X() + MomentumF.Y()*MomentumF.Y() + MomentumF.Z()*MomentumF.Z() + mass*mass ) - mass); //I want this in MeV
              //        hInteractingKEEl->Fill(KEF);
              hInteractingKEEl->Fill(kineticEnergy);
              //        std::cout<<"inside hadElastic Loop  and KEF "<<KEF<<std::endl;
            }
          finalKE = kineticEnergy;
          keepInteraction = false;
          if (!interactionLabel.size())
            {
              throughgoing++;
              G4Process.push_back("throughgoing");
            }else
            {
              G4Process.push_back(interactionLabel);
            }



          ////***********************new segment of the code ends here******************************************************************************//
        }
        }//if finTPC!=inTPC
      }//if primary
    } //if particle 
    fEventTree->Fill();
  } // end of analyze function

void protoana::truepionXsection::beginJob ( )

virtual

Reimplemented from art::EDAnalyzer.

Definition at line 226 of file truepionXsection_module.cc.

                                 {
    std::cout<<"job begin..."<<std::endl;
    art::ServiceHandle<art::TFileService> tfs;
    fEventTree = tfs->make<TTree>("Event", "Event Tree from Reco");
    fEventTree->Branch("event", &event,"event/I");
    fEventTree->Branch("run", &run,"run/I");
    fEventTree->Branch("subrun", &subrun,"surbrun/I");

    //  h_DE   = tfs->make<TH1D>("h_DE","h_DE; Energy Deposited [MeV]",200, 0,100);   
    // h_DX   = tfs->make<TH1D>("h_DX","h_DX; Distance between points  [cm]",400, 0,20);   
    // h_DEDX = tfs->make<TH1D>("h_DEDEX","h_DEDX; dE/dX [MeV/cm]",500, 0,50);   
    h_DEUniform   = tfs->make<TH1D>("h_DEUniform","h_DE; Energy Deposited [MeV]",200, 0,100);   
    h_DXUniform   = tfs->make<TH1D>("h_DXUniform","h_DX; Distance between points  [cm]",400, 0,20);   
    h_DEDXUniform = tfs->make<TH1D>("h_DEDEXUniform","h_DEDX; dE/dX [MeV/cm]",500, 0,50);   
    h_DeltaE = tfs->make<TH1D>("h_DeltaE","h_DeltaE; dEDep - TrjDE [MeV/cm]",500, -1000,1000);   
    // h_SimIDEDist= tfs->make<TH1D>("h_SimIDEDist","h_SimIDEDist; h_SimIDEDist [cm]",1000, 0,10);   
  
    h_UniformDistances = tfs->make<TH1D>("h_UniformDistances","h_UniformDistances; Distance between uniform points  [cm]",500, 0,5);   
    hInitialPz     = tfs->make<TH1D>("hInitialPz"    , "Initial Pz [MeV/c]"    , 42, -100, 2000);
    hInitialKE     = tfs->make<TH1D>("hInitialKE"    , "Initial Kinetic Energy [MeV]"    , 42, -100, 2000);
    hKEAtTPCFF     = tfs->make<TH1D>("hKEAtTPCFF"    , "Kinetic Energy @ TPC FF [MeV]"   , 42, -100, 2000);
    hIncidentKE        = tfs->make<TH1D>("hIncidentKE"   , "Incident Kinetic Energy [MeV]"   , 42, -100, 2000); 

    processtot=tfs->make<TH1D>("processtot","Total number of interaction for each particle;number of interaction;number of entries",50,0,50);
    processEl=tfs->make<TH1D>("processEl","Total number of Elastic interaction for each particle;number of interaction;number of entries",50,0,50);
    processInel=tfs->make<TH1D>("processInel","Total number of InElastic interaction for each particle;number of interaction;number of entries",50,0,50);
    hInteractingKE     = tfs->make<TH1D>("hInteractingKE", "Interacting Kinetic Energy [MeV]", 42, -100, 2000);
    hInteractingKEEl   = tfs->make<TH1D>("hInteractingKEEl", "Elastic Interacting Kinetic Energy [MeV]", 42, -100, 2000); 
    hInteractingKEElDep   = tfs->make<TH1D>("hInteractingKEElDep", "Dep Elastic Interacting Kinetic Energy [MeV]", 42, -100, 2000); 
    hInteractingKEInel = tfs->make<TH1D>("hInteractingKEInel", "Inelastic Interacting Kinetic Energy [MeV]", 42, -100, 2000);
    hCrossSection     = tfs->make<TH1D>("hCrossSection"     , "Cross-Section [barn]"             , 42, -100, 2000);
    hCrossSectionEl   = tfs->make<TH1D>("hCrossSectionEl"   , "Elastic Cross-Section [barn]"     , 42, -100, 2000);
    hCrossSectionInel = tfs->make<TH1D>("hCrossSectionInel" , "Inelastic Cross-Section [barn]"   , 42, -100, 2000);
    hXZ     = tfs->make<TH2D>("hXZ"     , "hXZ"  , 895, -200, 695  , 720, -360, 360);  
    hYZ     = tfs->make<TH2D>("hYZ"     , "hYZ"  , 895, -200, 695  , 600, 0, 600); 
    hXZPre  = tfs->make<TH2D>("hXZPre"  , "hXZPre", 200, -100, 100 , 200, -100, 100); 
    hYZPre  = tfs->make<TH2D>("hYZPre"  , "hYZPre", 200, -100, 100 , 200, 300, 500); 
    hdEVsdX  = tfs->make<TH2D>("hdEVsdX"  , "hdEVsdX" , 504, -1, 50, 1100, -10, 100); 
    hdEVsKE  = tfs->make<TH2D>("hdEVsKE"  , "hdEVsKE" , 504, -1, 50, 220,  -10, 1000); 
  
    fEventTree->Branch("trueVtxX" ,&trueVtxX ,"trueVtxX/D");
    fEventTree->Branch("trueVtxY" ,&trueVtxY ,"trueVtxY/D");
    fEventTree->Branch("trueVtxZ" ,&trueVtxZ ,"trueVtxZ/D");
    fEventTree->Branch("trueEndX" ,&trueEndX ,"trueEndX/D");
    fEventTree->Branch("trueEndY" ,&trueEndY ,"trueEndY/D");
    fEventTree->Branch("trueEndZ" ,&trueEndZ ,"trueEndZ/D");
    fEventTree->Branch("finalKE"  ,&finalKE  ,"finalKE/D" );
    fEventTree->Branch("G4Process",&G4Process);


  }

void protoana::truepionXsection::beginRun

(

const art::Run &

run

)

Definition at line 280 of file truepionXsection_module.cc.

                                              {
    mf::LogInfo("truepionXsection")<<"begin run..."<<std::endl;
  }

void protoana::truepionXsection::endJob ( )

virtual

Reimplemented from art::EDAnalyzer.

Definition at line 567 of file truepionXsection_module.cc.

                               {     

    std::cout<<"-------------------------------------------"<<std::endl;
    std::cout<<"True Events pre-TPC .............. "<<evtsPreTPC<<std::endl;
    std::cout<<"True Events pre-TPC .............. "<<evtsInTheMiddle<<std::endl;
    std::cout<<"True Events post-TPC ............. "<<evtsPostTPC<<std::endl;
    std::cout<<"True Throughgoing    ............. "<<throughgoing<<std::endl;
    std::cout<<"True interactingInTPC ............ "<<interactingInTPC<<std::endl;
    std::cout<<"-------------------------------------------"<<std::endl;
    float rho = 1396; //kg/m^3
    float molar_mass = 39.95; //g/mol
    float g_per_kg = 1000; 
    float avogadro = 6.022e+23; //number/mol
    float number_density = rho*g_per_kg/molar_mass*avogadro;
    //float slab_width = 0.0047;//in m
    float slab_width = 0.004792;//in m
    // Calculate the Cross Section
    // ###################################################################
    // #### Looping over the exiting bins to extract the cross-section ###
    // ###################################################################
    for( int iBin = 1; iBin <= hInteractingKE->GetNbinsX(); ++iBin )
      {
        // ### If an incident bin is equal to zero then skip that bin ###
        if( hIncidentKE->GetBinContent(iBin) == 0 )continue; //Temporary fix to ensure that no Infinities are propagated to pad
   
        // ### Cross-section = (Exit Bins / Incident Bins) * (1/Density) * (1/slab width) ###
        float TempCrossSection = (hInteractingKE->GetBinContent(iBin)/hIncidentKE->GetBinContent(iBin)) * (1/number_density) * (1/slab_width);
   
        float elCrossSection   = ((hInteractingKEEl->GetBinContent(iBin)/hIncidentKE->GetBinContent(iBin)) * (1/number_density) * (1/slab_width) ) * (1/1e-28);
        float inelCrossSection = ((hInteractingKEInel->GetBinContent(iBin)/hIncidentKE->GetBinContent(iBin)) * (1/number_density) * (1/slab_width) ) * (1/1e-28);
        // ### Convert this into Barns ###
        float crossSection = TempCrossSection * (1/1e-28); 
        
        // ### Putting the value on the plot
        hCrossSection    ->SetBinContent(iBin,crossSection);
        hCrossSectionEl  ->SetBinContent(iBin,elCrossSection);
        hCrossSectionInel->SetBinContent(iBin,inelCrossSection);
        // ###########################################################
        // ### Calculating the error on the numerator of the ratio ###
        // ###########################################################
        float denomError = pow(hIncidentKE->GetBinContent(iBin),0.5);
        float denom = hIncidentKE->GetBinContent(iBin);
        if(denom == 0) continue; 
        float term2 = denomError/denom;
      
        float numError = pow(hInteractingKE->GetBinContent(iBin),0.5);
        float num = hInteractingKE->GetBinContent(iBin);
        float numErrorEl = pow(hInteractingKEEl->GetBinContent(iBin),0.5);
        float numEl = hInteractingKEEl->GetBinContent(iBin);
        float numErrorInel = pow(hInteractingKEInel->GetBinContent(iBin),0.5);
        float numInel = hInteractingKEInel->GetBinContent(iBin);
        // ### Putting in a protection against dividing by zero ###   
        if(num != 0){
          float term1 = numError/num;
          float totalError = (TempCrossSection) * (pow( ( (term1*term1) + (term2*term2) ),0.5)) * (1/number_density) * (1/slab_width) * (1/1e-28) *(1e26);
          hCrossSection->SetBinError(iBin,totalError);
        }
      
        if(numEl != 0){
          float term1 = numErrorEl/numEl;
          float totalError = (elCrossSection) * (pow( ( (term1*term1) + (term2*term2) ),0.5)) * (1/number_density) * (1/slab_width)  *(1e26);
          hCrossSectionEl->SetBinError(iBin,totalError);
        }
      
        if(numInel != 0){
          float term1 = numErrorInel/numInel;
          float totalError = (inelCrossSection) * (pow( ( (term1*term1) + (term2*term2) ),0.5)) * (1/number_density) * (1/slab_width)  *(1e26);
          hCrossSectionInel->SetBinError(iBin,totalError);
        }
      }//<---End iBin Loop


  }

void protoana::truepionXsection::reset

(

)

Definition at line 645 of file truepionXsection_module.cc.

                              {
    run = -99999;
    subrun = -99999;
    event = -99999;
 
    trueVtxX = -999.;
    trueVtxY = -999.;
    trueVtxZ = -999.;
    trueEndX = -999.;
    trueEndY = -999.;
    trueEndZ = -999.;
    trueEndZ = -999.;
    finalKE  = -999.;
    G4Process.clear(); 


  }

Member Data Documentation

Int_t protoana::truepionXsection::event

private

Definition at line 131 of file truepionXsection_module.cc.

int protoana::truepionXsection::evtsInTheMiddle = 0

private

Definition at line 178 of file truepionXsection_module.cc.

int protoana::truepionXsection::evtsPostTPC = 0

private

Definition at line 179 of file truepionXsection_module.cc.

int protoana::truepionXsection::evtsPreTPC = 0

private

Definition at line 177 of file truepionXsection_module.cc.

const art::InputTag protoana::truepionXsection::fBeamModuleLabel

private

Definition at line 119 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fCalorimetryModuleLabel

private

Definition at line 184 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fCalorimetryTag

private

Definition at line 120 of file truepionXsection_module.cc.

TTree* protoana::truepionXsection::fEventTree

private

Definition at line 128 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fGeneratorTag

private

Definition at line 124 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fHitsModuleLabel

private

Definition at line 182 of file truepionXsection_module.cc.

double protoana::truepionXsection::finalKE

private

Definition at line 168 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fPFParticleTag

private

Definition at line 123 of file truepionXsection_module.cc.

bool protoana::truepionXsection::fSaveCaloInfo

private

Definition at line 185 of file truepionXsection_module.cc.

bool protoana::truepionXsection::fSaveTrackInfo

private

Definition at line 186 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fShowerTag

private

Definition at line 122 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fTrackerTag

private

Definition at line 121 of file truepionXsection_module.cc.

std::string protoana::truepionXsection::fTrackModuleLabel

private

Definition at line 183 of file truepionXsection_module.cc.

std::vector<std::string> protoana::truepionXsection::G4Process

private

Definition at line 169 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::h_DEDXUniform

private

Definition at line 137 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::h_DeltaE

private

Definition at line 138 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::h_DEUniform

private

Definition at line 135 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::h_DXUniform

private

Definition at line 136 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::h_UniformDistances

private

Definition at line 140 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hCrossSection

private

Definition at line 146 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hCrossSectionEl

private

Definition at line 147 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hCrossSectionInel

private

Definition at line 148 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hdEVsdX

private

Definition at line 156 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hdEVsKE

private

Definition at line 157 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hIncidentKE

private

Definition at line 145 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInitialKE

private

Definition at line 150 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInitialPz

private

Definition at line 151 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInteractingKE

private

Definition at line 141 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInteractingKEEl

private

Definition at line 142 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInteractingKEElDep

private

Definition at line 143 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hInteractingKEInel

private

Definition at line 144 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::hKEAtTPCFF

private

Definition at line 149 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hXZ

private

Definition at line 152 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hXZPre

private

Definition at line 154 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hYZ

private

Definition at line 153 of file truepionXsection_module.cc.

TH2D* protoana::truepionXsection::hYZPre

private

Definition at line 155 of file truepionXsection_module.cc.

int protoana::truepionXsection::interactingInTPC = 0

private

Definition at line 181 of file truepionXsection_module.cc.

double protoana::truepionXsection::maxX = 360.0

private

Definition at line 172 of file truepionXsection_module.cc.

double protoana::truepionXsection::maxY = 600.0

private

Definition at line 174 of file truepionXsection_module.cc.

double protoana::truepionXsection::maxZ = 695.0

private

Definition at line 176 of file truepionXsection_module.cc.

double protoana::truepionXsection::minX = -360.0

private

Definition at line 171 of file truepionXsection_module.cc.

double protoana::truepionXsection::minY =0.0

private

Definition at line 173 of file truepionXsection_module.cc.

double protoana::truepionXsection::minZ = 0.0

private

Definition at line 175 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::processEl

private

Definition at line 159 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::processInel

private

Definition at line 160 of file truepionXsection_module.cc.

TH1D* protoana::truepionXsection::processtot

private

Definition at line 158 of file truepionXsection_module.cc.

Int_t protoana::truepionXsection::run

private

Definition at line 129 of file truepionXsection_module.cc.

Int_t protoana::truepionXsection::subrun

private

Definition at line 130 of file truepionXsection_module.cc.

int protoana::truepionXsection::throughgoing = 0

private

Definition at line 180 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueEndX

private

Definition at line 165 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueEndY

private

Definition at line 166 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueEndZ

private

Definition at line 167 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueVtxX

private

Definition at line 162 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueVtxY

private

Definition at line 163 of file truepionXsection_module.cc.

double protoana::truepionXsection::trueVtxZ

private

Definition at line 164 of file truepionXsection_module.cc.

ProtoDUNETruthUtils protoana::truepionXsection::truthUtil

private

below are some addition to be able to use protodune utitilites

Definition at line 116 of file truepionXsection_module.cc.

---

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

• protoduneana/protoduneana/singlephase/truepionXsection/truepionXsection_module.cc