helperFunctions.cxx

Go to the documentation of this file.

#include "G4Rwgt/helperFunctions.h"


//-------------------------------------

// TODO I think maybe this whole thing needs to be restructured as a class


uint GetTrajPointIdx( TLorentzVector daughter_pos, 
                          std::vector<TLorentzVector> traj_points)
// Returns index of traj_points vector whose corresponding element
// is smallest distance from daughter_pos.This is a way of using daughter 
// process info to set primary processes.
//
// We set a threshold. If the separation between the points is less than the
// threshold, the association is made.

// This is hacky, but the reason it is necessary is that (for some reason)
// particle->Trajectory()->TrajectoryProcesses() is not filled in our
// simulation
{

  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;
}

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

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

  return length;

}

bool isPrimary(const sim::Particle* part)
// does this particle exit the nucleus and get handed to geant4
{
  std::string process = part->Process();

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

  return false;
}


G4ReweightTraj makeTraj(const sim::Particle* part, 
                        art::Event& evt,                       
                        art::ServiceHandle<cheat::BackTracker>  bt)
// Constructs trajectory object geant4reweight needs 
{


  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)
// Return true if any of the particles enter the detector,
// otherwise return false
{
  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)
// Return true if any of the particles have endpoint  in the detector,
// otherwise return false
{
  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)
// Return true if interaction vertex is inside detector
{
  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)
// Return true if all pions stay within detector,                          
// otherwise return false                                                  
{
  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;
}