PionSIManager Class Reference

#include <PionInteractionSystematic.hxx>

Public Member Functions
	PionSIManager ()
	ctor
virtual	~PionSIManager ()
	dtor
virtual Bool_t	InVOI (SubDetId_h det, Float_t *pos) const
virtual Bool_t	InVOI1 (Float_t *pos) const
virtual Bool_t	InVOI2 (Float_t *pos) const
virtual Bool_t	InVOIext (Float_t *pos) const
SteppingResult	StepBetweenPoints (SubDetId_h det, Int_t charge, TLorentzVector initPos, TVector3 initMom, TLorentzVector finalPos, Double_t stepLength, TGeoManager *geom) const
PionInteractionSystematic *	ComputePionWeightInfo (const AnaEventB &event, SubDetId_h det) const
PionInteractionSystematic *	ComputePionWeightInfo (const AnaEventB &event, const SelectionBase &sel, Int_t branch) const

Detailed Description

Combine everything that may be analysis/selection dependent int a simple class, so that can extended

Definition at line 132 of file PionInteractionSystematic.hxx.

Member Function Documentation

ComputePionWeightInfo()

PionInteractionSystematic * PionSIManager::ComputePionWeightInfo	(	const AnaEventB &	event,
		SubDetId_h	det
	)		const

Calculates the information needed to compute an event weight, as well as a weight correcting Geant4 to Data. Returns a PionInteractionSystematic object containing all the needed information The detector can be provided externally (not only through the selection since for the moment the info is filled prior to any cuts sequence applied) to give more flexibility

Definition at line 767 of file PionInteractionSystematic.cxx.

                                                                                                      {

  // Get the geometry from the GeometryManager
  TGeoManager* pionSIGeom = ND::hgman().GeoManager();

  //Only load in the geometry if it has not already been loaded by the InputConverter. Use in this case the default file
  if (!pionSIGeom){
    ND::hgman().LoadGeometry();  // Use default file
    pionSIGeom = ND::hgman().GeoManager();
  }

  AnaTrueParticleB* allTrajInTPCFGD[NMAXTRUEPARTICLES];
  int nTraj = anaUtils::GetAllTrajInBunch(event, allTrajInTPCFGD);
  if((UInt_t) nTraj>NMAXTRUEPARTICLES) nTraj = NMAXTRUEPARTICLES;

  //Declare the object the result will be written into.
  //If the geometry has not already been loaded, this will load the geometry.
  PionInteractionSystematic* result = new PionInteractionSystematic();
  //Now, go through and select only those trajectories that start in the volume we're concerned with.
  //FGD1FV plus dead material between that and TPC gas.  
  //It is a trajectory of interest if it starts in the volume of interest, as defined by InVOI, or its parent ends
  //in the volume of interest.
  std::vector<AnaTrueParticleB*> trajOfInterest;
  for (Int_t it = 0; it < nTraj; it++){
    
    AnaTrueParticleB* track = allTrajInTPCFGD[it];

    //If the trajectory starts in the volume of interest, save it.
    //One contiguous VOI, just use the one function.
    if(InVOI(det, track->Position) ){
      trajOfInterest.push_back(track);    
    }
    else{      
      //Otherwise, find its parent, and if the parent ends in the volume of interest,
      //save this trajectory.
      for (Int_t jt = 0; jt < nTraj; jt++){
        AnaTrueParticleB* track2 = allTrajInTPCFGD[jt];
        if(track->ParentID == track2->ID){
          if(InVOI(det, track2->PositionEnd) )
            trajOfInterest.push_back(track);
          break;
        }
      }
    }
  }

  //So at this point we have a vector of all true trajectories that start in the volume of interest and those
  //whose parents end in the volume of interest.  (This is to guard against unsaved trajectories resulting
  //in an apparent child being created outside the volume of interest, but corresponding to the interaction
  //that ended the trajectory in the volume of interest.)

  //We want to go through this list of trajectories, picking out pion trajectories and taking all trajectories with
  //a pi+/- as a parent and assembling them into Interactions.
  std::vector<AnaTrueParticleB*> pionTrajs;
  std::vector<Interaction*> interactions;

  for(UInt_t oai = 0; oai < trajOfInterest.size(); oai++){
    //If it's a pi+/- trajectory beginning in the volume of interest, save to pionTrajs.
    if (trajOfInterest[oai]->PDG == 211 || trajOfInterest[oai]->PDG == -211){
      if(InVOI(det, trajOfInterest[oai]->Position) )
        pionTrajs.push_back(trajOfInterest[oai]);
    }
    
    //If it lists a pion as a parent, assign it to an Interaction.
    if(trajOfInterest[oai]->ParentPDG == 211 || trajOfInterest[oai]->ParentPDG == -211){
      
      AnaTrueParticleB* parent = GetParent(trajOfInterest[oai],allTrajInTPCFGD,nTraj);
      Bool_t matchFound = false;
      for (UInt_t jsi = 0; jsi < interactions.size(); jsi++){        
        if (interactions[jsi]->IncludesTrajectory(trajOfInterest[oai])){          
          interactions[jsi]->AddTrajectory(trajOfInterest[oai],parent);
          matchFound = true;
          break;          
        }
      }
            
      //If an existing interaction does not match the trajectory, make a new interaction for
      //it to reside in, and add that to the interactions vector.
      if(!matchFound){        
        Interaction *newInteraction = new Interaction(trajOfInterest[oai], parent);
        interactions.push_back(newInteraction);
      }              
    } 
  }
      
      
  //So at this point we have all pion trajectories identified, and all non-primary trajectories assembled into
  //interactions.
  //Now, loop through the pion trajectories and the interactions, and assign them to a specific parent pion.
  //Then, use those trajectories to determine the secondary interaction type of the pion.
  //Loop over the pions.
  std::vector< std::pair<AnaTrueParticleB*,Int_t> > pionTrajAndInteraction;
      
  //Might as well count the interactions we're interested in now, so we can fill arrays later.
  Int_t nInteractions = 0;
  for (UInt_t p = 0; p < pionTrajs.size(); p++){
    
    AnaTrueParticleB* pionTraj = pionTrajs[p];
    TLorentzVector pos(pionTraj->Position[0],pionTraj->Position[1],pionTraj->Position[2],pionTraj->Position[3]);
    TLorentzVector posEnd(pionTraj->PositionEnd[0],pionTraj->PositionEnd[1],pionTraj->PositionEnd[2],pionTraj->PositionEnd[3]);
    //Only look for the interaction type if the pion trajectory ends inside of the volume of interest.    
#if (JMDEBUGVERSION == 0) || (JMDEBUGVERSION == 2) || (JMDEBUGVERSION == 3)
    if(InVOI(sel.GetDetectorFV(branch), pionTraj->PositionEnd)){
      
      //If there are two separate VOIs, since tracking ends at the end of the VOI, we
      //want to consider any pion trajectory that ends outside the VOI it started in
      //as having "escaped" (what this if statement elses to.)
#elif JMDEBUGVERSION == 1      
    if( (InVOI1( pionTraj->Position ) && InVOI1( pionTraj->PositionEnd )) || (InVOI2( pionTraj->Position ) && InVOI2( pionTraj->PositionEnd ))){
#endif
            
      //Define a vector to store the (relevant) interactions belonging to this pion.
      std::vector<Interaction*> thisPionInteractions;
      
      //Loop over the interactions vector and fill thisPionInteractions with the ones corresponding to 
      //this pion trajectory.
      for (UInt_t i = 0; i < interactions.size(); i++){
        
        //If the interaction lists the pion as a parent, add it to the list.
        //Also, we're only interested in interactions at or after the end of the pion trajectory.
        if (interactions[i]->parentID == pionTraj->ID
            && interactions[i]->position[3] >= pionTraj->PositionEnd[3]){
          thisPionInteractions.push_back(interactions[i]);
        }
      }
      
      //So, at this point should have all the Interactions considered a result of this pion.
      //Use this vector of interactions and the pion's position to determine its interaction type.
      //Go through the interactions:
      //A mu of correct charge originating from the trajectory end point => decay.
      //A mu originating in a later interactions indicates decay of a pion that need to
      //count as part of the original secondary interactions.
      //If it's pi0 decay within 10^-4 ns of the end point, add a pi0.
      //If it's outside of that time range, attribute it to a pion of the same type as the
      //incoming pion undergoing a charge exchange (so add another pi+/-)
      //Once have the full count of relevant particle types, use that to determine the
      //secondary interaction type.
      //Start by defining the counters.
      Int_t NPiPlus = 0;
      Int_t NPiZero = 0;
      Int_t NPiMinus = 0;
      Int_t NMuMinus = 0;
      Int_t NMuPlus = 0;
      Int_t NExotic = 0;
      
      for (UInt_t i = 0; i < thisPionInteractions.size(); i++){
        
        TLorentzVector thisPos(thisPionInteractions[i]->position[0],thisPionInteractions[i]->position[1],thisPionInteractions[i]->position[2],thisPionInteractions[i]->position[3]);
        //First check the interaction that shares the position with the end of the pion trajectory.
        //note that this condition only works if posEnd is a TLorentzVector!
        if(posEnd == thisPos){
          
          NPiPlus += thisPionInteractions[i]->NumberOfParticleType(211);
          NPiZero += thisPionInteractions[i]->NumberOfParticleType(111) 
            + thisPionInteractions[i]->NPiZeroFromDecayProducts();
          NPiMinus += thisPionInteractions[i]->NumberOfParticleType(-211);
          NMuPlus += thisPionInteractions[i]->NumberOfParticleType(-13);
          NMuMinus += thisPionInteractions[i]->NumberOfParticleType(13);
          NExotic += thisPionInteractions[i]->CountExoticParticles();     
        }
        
        //Otherwise it's time separated.
        //In the below case, assume the particles considered below originate from 
        //something coming out of the original pion trajectory ending vertex (as opposed
        //to some tertiary or beyond interaction of unsaved trajectories, which is unlikely.)
        else{
          
          //A time separated mu+ implies an unsaved pi+ trajectory.
          NPiPlus += thisPionInteractions[i]->NumberOfParticleType(-13);
          
          //A time separated mu- implies an unsaved pi- trajectory.
          NPiMinus += thisPionInteractions[i]->NumberOfParticleType(13);
          
          //If a pi0 decay was within 10^-4 ns of the pion trajectory end,
          //count it as being a product of the interaction the pion trajectory ended in.
          //If there was a saved pi0 trajectory in this time region, assume
          //that it originated from a time separated interaction. 
          Int_t TSPiZero = 0;
          if(thisPionInteractions[i]->position[3] - pionTraj->PositionEnd[3] < 10E-4){
            
            //Consider the decay products as originating from a pi0 from the
            //original pion's interaction.
            NPiZero += thisPionInteractions[i]->NPiZeroFromDecayProducts();
            
            //Consider any saved pi0 as originating from a time separated
            //interaction of an unsaved charged pion trajectory.
            TSPiZero += thisPionInteractions[i]->NumberOfParticleType(111); 
            
          }
          
          //Otherwise, count it or any saved pi0 as "Time Separated" pi0, which
          //we later assume originated from a time separated interaction of an
          //unsaved pi+/- trajectory.
          else{
            TSPiZero += thisPionInteractions[i]->NumberOfParticleType(111) 
              + thisPionInteractions[i]->NPiZeroFromDecayProducts();
          }
          
          //Any time separated pions or exotics are likely collectively due
          //to an unsaved pion trajectory of the same charge as the original pion
          //trajectory.  So count all pi+, time separated pi0, pi-, and exotics, and if the amount is
          //greater than 0, add a pion of the same charge as the parent of the
          //interaction.
          if (thisPionInteractions[i]->NumberOfParticleType(211)
              + TSPiZero
              + thisPionInteractions[i]->NumberOfParticleType(-211)
              + thisPionInteractions[i]->CountExoticParticles() > 0){
            
            if(pionTraj->PDG == 211){
              
              NPiPlus += 1;
            }
            
            else{
              NPiMinus += 1;
            }
            
            
          }
          
        }
        
      }
      
      //clean memory
      std::vector<Interaction*> interactions;
      for (std::vector<Interaction*>::iterator it = interactions.begin(); it != interactions.end(); it++){
        delete (*it); 
        
      }
      
      interactions.clear();
      
      
      
      //At this point we have all of the relevant particles in the pion's interaction counted.
      //Use these counts to determine the interaction type.
      Int_t interactionType;
      
      //First take care of decays.
      if((pionTraj->PDG == 211 && NMuPlus == 1)
         || (pionTraj->PDG == -211 && NMuMinus == 1)){
        
        interactionType = 3;
        
      }
      
      //Now slot anything with exotic particles in it into Other.
      else if(NExotic > 0){
        interactionType = 7;
      }
      
      //Multi-pi.
      //If it has more than 1 pion total.
      else if (NPiPlus + NPiZero + NPiMinus > 1){       
        interactionType = 5;
      }
      
      //CEX
      //Since we've excluded multi-pi events, only have to check the number
      //of the desired pion here.
      //This is CEX for both pi+ and pi-
      else if (NPiZero == 1){
        interactionType = 1;
      }
      
      //DCEX
      //Has a pi- if parent was pi+, or has a pi+ if parent was pi-
      else if ((pionTraj->PDG == 211 && NPiMinus == 1)
               || (pionTraj->PDG == -211 && NPiPlus == 1)){
        interactionType = 2;
      }
      
      //QE
      //If there is one pion of the same kind coming in going out.
      else if ((pionTraj->PDG == 211 && NPiPlus == 1)
               || (pionTraj->PDG == -211 && NPiMinus == 1)){
        interactionType = 4;
      }
      
      //ABS
      //If there's no pions or muons, it's absorption.
      else if (NPiPlus + NPiZero + NPiMinus
               + NMuMinus + NMuPlus == 0){
        interactionType = 0;
      }
      
      //As a catch all for weird stuff that shouldn't exist
      //(everything is expected to have been caught above.)
      else{
        interactionType = 7;
      } 
      
      
      //With the interaction type now determined, write this pion trajectory and its
      //interaction type out to the vector that stores them.
      pionTrajAndInteraction.push_back(std::make_pair(pionTraj,interactionType));
    } 
    
    //If it's not in the volume of interest, assign it the "Escaped" interaction type (8).
    else{
      pionTrajAndInteraction.push_back(std::make_pair(pionTraj,8));
      
    }
    
    //Now, check the recently just added pair to see whether it is one of the interaction
    //types we're interested in.  If it is, count as one of the nInteractions.
    Int_t thisIntType =  pionTrajAndInteraction.back().second;
    
    //If it's one of the interactions we're interested in, count it and
    //store its information for later processing.
    if(thisIntType == 0 || thisIntType == 1 || thisIntType == 4){
      
      nInteractions++;
    }
    
    
  } //Close loop over pion trajectories.
        
  //So, at this point we now have a vector of pion trajectories paired with the interaction that ends them
  //(or if this interaction occured outside the volume of interest, indication that they escaped.)
  //Now we need to step these pions through the detector.
  //Propagate from beginning to the end of the trajectory.  Elastic interactions aren't accounted for,
  //but the error due to that should be small.  Propagates out to length between start and end point,
  //either stopping when it hits that length or when it goes outside of the volume of interest.
        
  std::vector<Float_t> allPionSteps;
  
  Float_t AllPionsMCSumNSigmaStepLength = 0.0;
  Float_t AllPionsDataSumNSigmaStepLength = 0.0;
  
  //Set up as many variables as can.
  Int_t nPions = (Int_t)pionTrajAndInteraction.size();
  Bool_t* pionType = new Bool_t[(UInt_t)nPions];
  Int_t* nSteps = new Int_t[(UInt_t)nPions];
  Float_t* initMom = new Float_t[(UInt_t)nPions];
  Float_t* pInteraction =  new Float_t[(UInt_t)nInteractions];
  Int_t* typeInteraction = new Int_t[(UInt_t)nInteractions];
  Int_t interactionsSoFar = 0; // Since not every pion has an interaction to count.
  
  //Variables for validation purposes.
  std::vector<Float_t> testFinalMom;
  std::vector<TLorentzVector> testFinalPos;
  std::vector<Int_t> testInteractionType;
  Int_t* pionID = new Int_t[nPions];
  AnaTrueParticleB** TrueParticles = new AnaTrueParticleB*[nPions];
  AnaTrueParticleB** InteractionTrueParticles = new AnaTrueParticleB*[nInteractions];
  
  for(UInt_t ip = 0; ip < pionTrajAndInteraction.size(); ip++){
    
    AnaTrueParticleB* pionTraj = pionTrajAndInteraction[ip].first;
    Int_t intCode = pionTrajAndInteraction[ip].second;
    
    //Set whether it is a pi+ or a pi-.
    if(pionTraj->PDG == 211){
      
      pionType[ip] = true;
      
    }
    else{
      
      pionType[ip] = false;
      
    }
    
    //Store the initial momentum.
    //Remember: In highland the magnitude of the momentum
    //is stored as Momentum, and the direction elsewhere.
    initMom[ip] = pionTraj->Momentum;
    
    //Store the pion trajectory ID, for validation purposes.
    pionID[ip] = pionTraj->ID;

    TrueParticles[ip] = pionTraj;
    
    TVector3 dir(pionTraj->Direction[0], pionTraj->Direction[1], pionTraj->Direction[2]);
    TLorentzVector pos(pionTraj->Position[0],pionTraj->Position[1],pionTraj->Position[2],pionTraj->Position[3]);
    TLorentzVector posEnd(pionTraj->PositionEnd[0],pionTraj->PositionEnd[1],pionTraj->PositionEnd[2],pionTraj->PositionEnd[3]);
        
    SteppingResult sr = StepBetweenPoints(det, (pionTraj->PDG)/211, //Charge
                                          pos,
                                          (pionTraj->Momentum)*dir,
                                          posEnd,
                                          0.1, //Step length, 0.1 mm
                                          pionSIGeom); //The geometry.
        
    //Store the number of steps for this pion.
    nSteps[ip] = sr.stepLengths.size();
    //Fill the test variables.
    testFinalMom.push_back(sr.finalMom);
        
    testFinalPos.push_back(sr.finalPos);

    testInteractionType.push_back(intCode);
        
    //If it's one of the interactions we're interested in, store the information
    //we need to store to generate weights.
    if(intCode == 0 || intCode == 1 || intCode == 4){
      
          pInteraction[interactionsSoFar] = sr.finalMom;
          
          //Get the material this interaction occurred in.
          //(Use the true MC end as opposed to the propagated one, which is not exact.
          //The interaction necessarily occurred at the end of the trajectory
          TGeoNode *node = pionSIGeom->FindNode(pionTraj->PositionEnd[0],
                                            pionTraj->PositionEnd[1],
                                            pionTraj->PositionEnd[2]);
          
          TGeoVolume *volume = node->GetVolume();
          TGeoMaterial *material = volume->GetMaterial();
          std::string materialName = material->GetName();
          
          Int_t materialSeries = GetMaterialSeries(materialName);
          
          
          //Get the pionType Series to add to the interaction type.
          //0 for pi+, 100 for pi-
          //Int_t pionTypeSeries = (100)*(1-pionType[ip]);
          //pionType == true for pi+, false for pi-
          Int_t pionTypeSeries = 0;
          if(!pionType[ip]){
            pionTypeSeries = 100;
          }
          
          
          typeInteraction[interactionsSoFar] = intCode + materialSeries + pionTypeSeries;
          
      InteractionTrueParticles[interactionsSoFar] = pionTraj;
      
          //Increase the interaction counter to move on to
      //the next element.
          interactionsSoFar++;
    }
        
    
        
    //Loop over all the steps for this pion, and add them to allPionSteps
    for(UInt_t sIdx = 0; sIdx < sr.stepLengths.size() ; sIdx++){
      
          allPionSteps.push_back(sr.stepLengths[sIdx]);
      
    }
        
    //Add this pion's MC and Data sum N*sigma*StepLength to the total for
    //all pions in the event.
    AllPionsMCSumNSigmaStepLength += sr.MCSumNSigmaStepLength;
    AllPionsDataSumNSigmaStepLength += sr.DataSumNSigmaStepLength;
        
        
  }

  
  //Can't fill some output variables until know everything there is to know about the event.  Do that here.
  Int_t totalSteps = (Int_t)allPionSteps.size();
  
  //Can now create and fill the stepLengths array.
  Float_t* stepLengths = new Float_t[totalSteps];
  
  for(UInt_t i = 0; i < allPionSteps.size(); i++){
    
    stepLengths[i] = allPionSteps[i];
  }
  
  
  //Use the Step sums and the interactions arrays to compute the MC to data weights.
  //Start with the weight portion from the stepping.
  Float_t weightMCToData = TMath::Exp(-(AllPionsDataSumNSigmaStepLength - AllPionsMCSumNSigmaStepLength));
  
  //Now, go through the interactions and multiply the exponential by their contributions.
  for(Int_t ni = 0; ni < nInteractions; ni++){
    
    Int_t siType = typeInteraction[ni];
    Float_t siMom = pInteraction[ni];
    
    int momBin = int(siMom/10);
    if (momBin>800) momBin = 800;
    
    Double_t dataXSec = Interpolate(PiXSec::xsec_data[siType][momBin], PiXSec::xsec_data[siType][momBin+1],momBin,siMom);
    Double_t mcXSec = Interpolate(PiXSec::xsec_MC[siType][momBin], PiXSec::xsec_MC[siType][momBin+1],momBin,siMom);
    
    //Potential 2000 MeV/c cutoff: If want not to weight pions above 2000 MeV/c,
    //uncomment the next two commented lines.
    //if(pInteraction[ni] < 2000.0){
    if(mcXSec != 0.0){
      weightMCToData *= dataXSec/mcXSec;
    }
    
    //}
  }
    
  //Now all the variables that are needed have been filled.  They need to be written out to the microtree.
  //Pass them to our ad-hoc object.  May have a better solution later.
  //PionInteractionSystematic result;
  //This is declared at the beginning of this function now, so that the
  //geometry initialization code is called before try to use the geometry.
  result->nPions = nPions;
  result->pionType = pionType;
  result->totalSteps = totalSteps;
  result->nSteps = nSteps;
  result->initMom = initMom;
  result->stepLengths = stepLengths;
  result->nInteractions = nInteractions;
  result->pInteraction = pInteraction;
  result->typeInteraction = typeInteraction;
  result->testFinalMom = testFinalMom;
  result->testFinalPos = testFinalPos;
  result->testInteractionType = testInteractionType;
  result->weightMCToData = weightMCToData;
  result->MCSumNSigmaStepLength = AllPionsMCSumNSigmaStepLength;
  result->DataSumNSigmaStepLength = AllPionsDataSumNSigmaStepLength;
  result->MCProbNoInt = TMath::Exp(-AllPionsMCSumNSigmaStepLength);
  result->DataProbNoInt = TMath::Exp(-AllPionsDataSumNSigmaStepLength);
  result->pionID = pionID;
  result->TrueParticles = TrueParticles;
  result->InteractionTrueParticles = InteractionTrueParticles;
  
  return result;
    
}

---

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

• /hep/T2K/nd280rep/ANALYSIS/HIGHLAND2/ANAPARTICLE/psyche/psycheND280Utils/v3r12/src/PionInteractionSystematic.hxx
• /hep/T2K/nd280rep/ANALYSIS/HIGHLAND2/ANAPARTICLE/psyche/psycheND280Utils/v3r12/src/PionInteractionSystematic.cxx