helperFunctions.h File Reference

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "ReweightBase/G4MultiReweighter.hh"
#include "ReweightBase/G4ReweightTraj.hh"
#include "ReweightBase/G4ReweightStep.hh"
#include "PropBase/G4ReweightParameterMaker.hh"
#include "Simulation/Particle.h"
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "TGeoMaterial.h"
#include "GeometryObjects/GeometryBase.h"

Go to the source code of this file.

Functions
uint	GetTrajPointIdx (TLorentzVector daughter_pos, std::vector< TLorentzVector > traj_points)
double	len (double x1, double x2, double x3)
bool	isPrimary (const sim::Particle *part)
G4ReweightTraj	makeTraj (const sim::Particle *part, art::Event &evt, art::ServiceHandle< cheat::BackTracker > bt)
bool	RockFilter_A (std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)
bool	RockFilter_B (std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)
bool	RockFilter_C (art::Ptr< simb::MCTruth > mctruth, art::ServiceHandle< geo::Geometry > geom)
bool	RockFilter_D (std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)

Function Documentation

uint GetTrajPointIdx	(	TLorentzVector	daughter_pos,
		std::vector< TLorentzVector >	traj_points
	)

Definition at line 9 of file helperFunctions.cxx.

References len(), and demo0::length.

Referenced by makeTraj().

{

  double thresh = 0.001; 

  TLorentzVector elem;
  double elem_X;
  double elem_Y;
  double elem_Z;

  double daughter_X = daughter_pos(0);
  double daughter_Y = daughter_pos(1);
  double daughter_Z = daughter_pos(2);  

  double deltaX;
  double deltaY;
  double deltaZ;

  double length;

  uint numPoints = traj_points.size();

  uint traj_idx = 12345; // Some ridiculous sentinel value
                            // Can't use -5 as it is unsigned int

  int min_len = 100.0;  // Some large initial value

  for(uint point_idx=0; point_idx < numPoints; point_idx++){
 
    elem = traj_points.at(point_idx);
    elem_X = elem(0);
    elem_Y = elem(1);
    elem_Z = elem(2);
    
    deltaX = elem_X - daughter_X;
    deltaY = elem_Y - daughter_Y;
    deltaZ = elem_Z - daughter_Z;

    length = len(deltaX, deltaY, deltaZ);

    if(length < min_len){ 
      min_len = length;
    }
    else{
      continue;
    }

    if(length < thresh){
      /*
      if(traj_idx != 12345){ // Has traj_idx already been set
        
        std::cout << "***********************************" << std::endl;
        std::cout << "WARNING - two points below threshold" << std::endl;
        std::cout << "Ambigous association in GetTrajPointIdx() " << std::endl;
        std::cout << "***********************************" << std::endl;
        
      }
      */
      traj_idx = point_idx;
      
    }
  }

  return traj_idx;
}

bool isPrimary

(

const sim::Particle *

part

)

Definition at line 100 of file helperFunctions.cxx.

References launch_batch_jobs::process, simb::MCParticle::Process(), and string.

{
  std::string process = part->Process();

  if(strcmp( process.c_str(), "Primary") == 0){
    return true;
  }

  return false;
}

double len	(	double	x1,
		double	x2,
		double	x3
	)

Definition at line 86 of file helperFunctions.cxx.

References demo0::length, cet::pow(), and std::sqrt().

Referenced by GetTrajPointIdx(), and makeTraj().

{

  // Returns length of vector
  double length = sqrt(
                      std::pow( x1, 2 )  +
                      std::pow( x2, 2 )  +
                      std::pow( x3, 2 )
                      );

  return length;

}

G4ReweightTraj makeTraj	(	const sim::Particle *	part,
		art::Event &	evt,
		art::ServiceHandle< cheat::BackTracker >	bt
	)

Definition at line 113 of file helperFunctions.cxx.

References simb::MCParticle::Daughter(), energy, art::EventID::event(), GetTrajPointIdx(), MECModelEnuComparisons::i, art::Event::id(), compare_h5_caf::idx, len(), demo0::length, make_pair(), simb::MCParticle::Mass(), make_training::momentum, simb::MCParticle::NumberDaughters(), simb::MCParticle::NumberTrajectoryPoints(), simb::MCParticle::PdgCode(), simb::MCParticle::Position(), cet::pow(), simb::MCParticle::Process(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), PandAna.Demos.demo1::range, std::sqrt(), string, simb::MCParticle::TrackId(), cheat::BackTracker::TrackIDToParticle(), X, Y, and Z.

Referenced by g4rwgt::MakeG4WeightTable::produce().

{


  int event = evt.id().event();
 
  // ~~~ Assign values to variables ~~~

  int prim_pdg        = part->PdgCode(); 

  int prim_id    = part->TrackId();

  int prim_nDaughters = part->NumberDaughters();      

  // Don't know why these two are needed, but G4ReweightTraj needs them
  // The values I'm passing I assume do nothing
  int parent_id = 0;
  std::pair<int,int> range = std::make_pair(0,0);

  double mass;
  mass = part->Mass();
  mass *= 1.e3; // geant4reweight wants mass in MeV


  //Make a G4ReweightTraj -- This is the reweightable object
  G4ReweightTraj theTraj(prim_id, prim_pdg, parent_id, 
                        event, range);

  uint nTrajPoints = part->NumberTrajectoryPoints();
  std::vector<TLorentzVector> trajPoints;

  // Process string at each point in trajectory. Initialise values to "none"
  std::string processes[nTrajPoints];
  std::fill_n(processes, nTrajPoints, "none");

  std::vector<double> prim_X;
  std::vector<double> prim_Y;
  std::vector<double> prim_Z;

  std::vector<double> prim_Px;
  std::vector<double> prim_Py;
  std::vector<double> prim_Pz;

  //art::ServiceHandle<geo::Geometry> geom; //TODO probs bad to declare
                                          //     service handle in a
                                          //     function which will
                                          //     be called multiple times

  for( uint step_idx = 0; step_idx < nTrajPoints; 
       step_idx++){

    trajPoints.push_back( part->Position(step_idx) );

    double X = part->Position(step_idx).X();
    double Y = part->Position(step_idx).Y();
    double Z = part->Position(step_idx).Z();

    /*
    const TGeoMaterial* testmaterial1 = geom->Material( X,Y,Z );

    std::cout << "At point (" << X << ", " << Y << ", " << Z << ") " << testmaterial1->GetName()
              << " " << testmaterial1->GetZ() << " " << testmaterial1->GetA() 
              << " " << testmaterial1->GetDensity() << std::endl;    
    */
 
    prim_X.push_back(X);
    prim_Y.push_back(Y);
    prim_Z.push_back(Z);

    prim_Px.push_back(part->Px(step_idx));
    prim_Py.push_back(part->Py(step_idx));
    prim_Pz.push_back(part->Pz(step_idx));

  } // End loop over trajectory points

  const TLorentzVector prim_start = part->Position(0);
  const TLorentzVector prim_end = part->Position(nTrajPoints -1);

  for(int i=0; i < prim_nDaughters; i++){
    
    const int daughter_idx = part->Daughter(i);

    const sim::Particle* daughter = 
      bt->TrackIDToParticle(daughter_idx);

    const int daughter_pdg  = daughter->PdgCode(); 
    std::string daughter_proc = daughter->Process();

    theTraj.AddChild(new G4ReweightTraj(daughter_idx, daughter_pdg,
                                        prim_id, event, range) );


    const TLorentzVector daughter_start = daughter->Position(0);

    //****************
    uint trajPointIdx = GetTrajPointIdx(daughter_start, trajPoints);
    // TODO: There is redundancy here as e.g. at the end of the trajectory
    //       where there may be many daughters at the same point we are
    //       setting the process at that point once for each daughter.
    if(trajPointIdx != 12345){
      processes[trajPointIdx] = daughter->Process();
    }
    //****************
    

  }// end loop over pion daughters

  std::string endProc = processes[nTrajPoints -1];

  if( strcmp( endProc.c_str(), "hadElastic") == 0  ){
    // What is anyone supposed to do if they see this warning ??
    //std::cout << "WARNING: last step and elastic" << std::endl;
  }
  
    
  // Loop where we add steps to traj
  // This could possibly go in an existing loop for the sake of 
  // efficiency ?? 
  for(uint idx=1; idx < nTrajPoints; idx++){

    double deltaX = ( prim_X.at(idx) - prim_X.at(idx -1) );
    double deltaY = ( prim_Y.at(idx) - prim_Y.at(idx -1) );
    double deltaZ = ( prim_Z.at(idx) - prim_Z.at(idx -1) );

    double length = len(deltaX, deltaY, deltaZ);

    double preStepP[3] = {
      prim_Px.at(idx-1)*1.e3,  // geant4reweight wants Momentum in MeV
      prim_Py.at(idx-1)*1.e3,
      prim_Pz.at(idx-1)*1.e3
    };

    double postStepP[3] = {
      prim_Px.at(idx)*1.e3,
      prim_Py.at(idx)*1.e3,
      prim_Pz.at(idx)*1.e3
    };


    if(idx == 1){ 

      double momentum = sqrt( std::pow(preStepP[0], 2) +
                              std::pow(preStepP[1], 2) +
                              std::pow(preStepP[2], 2)
                            );


      // Set initial energy
      double energy = sqrt( std::pow(momentum, 2) +
                            std::pow(mass, 2)
                         );

      theTraj.SetEnergy(energy);
    }

    std::string proc  = processes[idx];


    
    // Events where a single geant4 step is huge
    // seem to spit out huge weights, so lets skip those
    //
    // TODO this could be an area of further investigation
    //      re: non-unitarity
    if(length > 250.0) continue;
    

    G4ReweightStep * theStep = new G4ReweightStep( prim_id, prim_pdg, 
                                                   parent_id,  event,
                                                   preStepP, postStepP,
                                                   length, proc );

    theStep->SetDeltaX(deltaX);
    theStep->SetDeltaY(deltaY);
    theStep->SetDeltaZ(deltaZ);

    theTraj.AddStep(theStep);

  } // end loop over trajPoints

  return theTraj;
}

bool RockFilter_A	(	std::vector< const sim::Particle * >	parts,
		art::ServiceHandle< geo::Geometry >	geom
	)

Definition at line 299 of file helperFunctions.cxx.

References MECModelEnuComparisons::i, geo::GeometryBase::isInsideDetectorBigBox(), calib::j, simb::MCParticle::NumberTrajectoryPoints(), part, and simb::MCParticle::Position().

{
  for(uint i=0; i<parts.size(); i++){
    const sim::Particle *part = parts.at(i);
    uint nTrajPoints = part->NumberTrajectoryPoints();
    for(uint j=0; j<nTrajPoints; j++){
      if(geom->isInsideDetectorBigBox( part->Position(j).X(),
                                       part->Position(j).Y(),
                                       part->Position(j).Z() ) 
         ){
        return true;
      }

    }// end for(trajPoints)

  }// end for(parts)
  
  return false;
}

bool RockFilter_B	(	std::vector< const sim::Particle * >	parts,
		art::ServiceHandle< geo::Geometry >	geom
	)

Definition at line 323 of file helperFunctions.cxx.

References MECModelEnuComparisons::i, geo::GeometryBase::isInsideDetectorBigBox(), simb::MCParticle::NumberTrajectoryPoints(), part, and simb::MCParticle::Position().

{
  for(uint i=0; i<parts.size(); i++){
    const sim::Particle *part = parts.at(i);
    uint nTrajPoints = part->NumberTrajectoryPoints();

    if(geom->isInsideDetectorBigBox( part->Position(nTrajPoints -1).X(),
                                     part->Position(nTrajPoints -1).Y(),
                                     part->Position(nTrajPoints -1).Z() ) 
       ){
      return true;
    }



  }// end for(parts)
  
  return false;
}

bool RockFilter_C	(	art::Ptr< simb::MCTruth >	mctruth,
		art::ServiceHandle< geo::Geometry >	geom
	)

Definition at line 346 of file helperFunctions.cxx.

References simb::MCTruth::GetNeutrino(), geo::GeometryBase::isInsideDetectorBigBox(), simb::MCNeutrino::Nu(), and simb::MCParticle::Position().

Referenced by g4rwgt::MakeG4WeightTable::produce().

{
  const TVector3 vtx = mctruth->GetNeutrino().Nu().Position().Vect();
  
  if(geom->isInsideDetectorBigBox(vtx.X(), vtx.Y(), vtx.Z())){
    return true;
  }
  return false;

}

bool RockFilter_D	(	std::vector< const sim::Particle * >	parts,
		art::ServiceHandle< geo::Geometry >	geom
	)

Definition at line 359 of file helperFunctions.cxx.

References abs(), MECModelEnuComparisons::i, geo::GeometryBase::isInsideDetectorBigBox(), calib::j, simb::MCParticle::NumberTrajectoryPoints(), part, simb::MCParticle::PdgCode(), and simb::MCParticle::Position().

{
  for(uint i=0; i<parts.size(); i++){
    const sim::Particle *part = parts.at(i);

    if(abs(part->PdgCode()) != 211) continue;

    uint nTrajPoints = part->NumberTrajectoryPoints();
    for(uint j=0; j<nTrajPoints; j++){
      if(!geom->isInsideDetectorBigBox( part->Position(j).X(),
                                        part->Position(j).Y(),
                                        part->Position(j).Z() )
         ){
        return false;
      }

    }// end for(trajPoints)                                                

  }// end for(parts)                                                       

  return true;
}