TruthUtils.cxx

#include "TruthUtils.hxx"
#include "ConstituentsUtils.hxx"


//**************************************************
Int_t anaUtils::GetNMichelElectrons(const AnaTrueVertexB& trueVertex, SubDetId::SubDetEnum det){
  //**************************************************

  // out of FV
  if ( ! InFiducialVolume(det, trueVertex.Position)) return -999;

  /// particles coming out the vertex
  Int_t NMich=0;

  for (int tit = 0; tit < trueVertex.nTrueParticles; tit++) {
    AnaTrueParticleB* true_track = trueVertex.TrueParticles[tit];

    //if not a pi+/-, a e+/-, a mu+/-
    if( abs(true_track->GParentPDG)!=211) continue;
    if( abs(true_track->ParentPDG)!=13  ) continue;
    if( abs(true_track->PDG)!=11  )       continue;

    NMich++;    
  }

  //  std::cout<<" nb of michel electrons "<<NMich << " " << trueVertex.nTrueParticles <<std::endl;

  return  NMich;
}

//**************************************************
Float_t anaUtils::GetTrueLinearLengthInTPC(const AnaTrueParticleB& trueTrack, Float_t& distz){
  //**************************************************

  Float_t dist = -9999999;
  distz = -999;
  for (Int_t idet = 0; idet < trueTrack.nDetCrossings; idet++) {
    //i.e crossing the active part of the tpc
    if (!trueTrack.DetCrossings[idet]->Detector)
      continue;
    if (SubDetId::GetDetectorUsed(trueTrack.DetCrossings[idet]->Detector, SubDetId::kTPC) && trueTrack.DetCrossings[idet]->InActive) {
      Float_t dist_temp = anaUtils::GetSeparationSquared(trueTrack.DetCrossings[idet]->ExitPosition, trueTrack.DetCrossings[idet]->EntrancePosition);
      Float_t distz_temp = fabs(trueTrack.DetCrossings[idet]->ExitPosition[2] - trueTrack.DetCrossings[idet]->EntrancePosition[2]);
      if (dist_temp > dist)
        dist = dist_temp;
      if (distz_temp > distz)
        distz = distz_temp;

    }
  }

  return sqrt(dist);
}

//********************************************************************
int anaUtils::GetAllChargedTrajsInFgdECalInBunch(const AnaEventB& event, AnaTrueParticleB* outTraj[], SubDetId::SubDetEnum det){
  //********************************************************************
  int count = 0;

  if (!SubDetId::IsFGDDetector(det))
    return false;
  
  AnaTrueParticleB* tracks[NMAXTRUEPARTICLES];
  int ntracks = anaUtils::GetAllChargedTrajInBunch(event, tracks);
  if((UInt_t)ntracks>NMAXTRUEPARTICLES) ntracks = NMAXTRUEPARTICLES;

  for (int i=0; i<ntracks; i++){
    AnaTrueParticleB* track = tracks[i];
    if (!track) continue;

    //should satisfy some minimum criteria in FGD and ECal
    if (anaUtils::TrueParticleEntersDet(track, det) && anaUtils::TrueParticleEntersDet(track, SubDetId::kTECAL) && !anaUtils::TrueParticleCrossesTPC(track))
      outTraj[count++] = track;

  }

  return count;
}


//********************************************************************
int anaUtils::GetAllTrajsInECalInBunch(const AnaEventB& event, AnaTrueParticleB* outTraj[]){
  //********************************************************************
  int count = 0;


  for(Int_t i = 0; i < event.nTrueParticles;i++){
   
    AnaTrueParticleB* track = event.TrueParticles[i];
    
    if (!track) continue;
    
    if(!track->TrueVertex) continue;
    if(track->TrueVertex->Bunch != event.Bunch) continue;
  
    SubDetId::SubDetEnum ecal_det[20];
    
    int necal_det =  anaUtils::GetECalDetCrossed(track, ecal_det);

    if (necal_det > 0)
      outTraj[count++] = track;
  }

  return count;
}

//********************************************************************
bool anaUtils::TrueParticleEntersDet(const AnaTrueParticleB* track, SubDetId::SubDetEnum det){
  //********************************************************************

  if(!track)
    return false;
  
  
  //this is a check to account for no splitting of truth into sub-detectors for SMRD and ECa
  if (det > SubDetId::kDSECAL && det < SubDetId::kTPC){
    std::cout<< "anaUtils::TrueParticleEntersDet() WARNING: " << det << " is not available from true info " << std::endl;
  } 

  //check that a trajectory crosses volume
  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    if(!SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, det)) continue;

    return true;

  }

  return false;

}

//********************************************************************
bool anaUtils::TrueParticleCrossesTECal(const AnaTrueParticleB* track, Float_t min_sep){
  //********************************************************************

  if(!track)
    return false;

  Float_t surf_pos[4] = {
    DetDef::tecalTLmin[1],
    DetDef::tecalBLmax[1],
    DetDef::tecalLmin[0],
    DetDef::tecalRmax[0] 
  };

  //check that a trajectory crosses BarrelECal active volume
  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    AnaDetCrossingB* cross = track->DetCrossings[idet];
    if(!cross)
      continue;

    if(!(cross->Detector & (1<<SubDetId::kTECAL)) || !cross->InActive) continue;

    //since truth info does not allow splitting the ecal sub-volumes --> use hard-coded values 

    // top 
    Float_t sep = 0.;
    if (cross->ExitPosition[1]>1265.){
      sep = cross->ExitPosition[1]-surf_pos[0];
    }
    //bottom
    else if (cross->ExitPosition[1]<-1265.){
      sep = surf_pos[1]-cross->ExitPosition[1];
    }
    //left
    else if (cross->ExitPosition[0]>1300.){
      sep = cross->ExitPosition[0]-surf_pos[2];
    }
    //right
    else if (cross->ExitPosition[0]<-1300.){
      sep = surf_pos[3]-cross->ExitPosition[0];
    }

    if(sep>min_sep)
      return true;
  }

  return false;

}


//********************************************************************
bool anaUtils::TrueParticleCrossesSMRD(const AnaTrueParticleB* track, Float_t min_sep){
  //********************************************************************

  if(!track)
    return false;

  Float_t surf_pos[4] = {
    DetDef::smrdTLmin[1],
    DetDef::smrdBLmax[1],
    DetDef::smrd6Lmin[0],
    DetDef::smrd6Rmax[0] 
  };

  //loop through det crossings
  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    AnaDetCrossingB* cross = track->DetCrossings[idet];
    if(!cross)
      continue;

    //InActive flag is not set for SMRD in oaAnalysis
    if(!(cross->Detector & (1<<SubDetId::kSMRD)))
      continue;

    // top 
    Float_t sep = 0.;
    if (cross->ExitPosition[1]>2010.){
      sep = cross->ExitPosition[1]-surf_pos[0];
    }
    //bottom
    else if (cross->ExitPosition[1]<-2010.){
      sep = surf_pos[1]-cross->ExitPosition[1];
    }
    //left
    else if (cross->ExitPosition[0]>1830.){
      sep = cross->ExitPosition[0]-surf_pos[2];
    }
    //right
    else if (cross->ExitPosition[0]<-1830.){
      sep = surf_pos[3]-cross->ExitPosition[0];
    }

    if(sep>min_sep)
      return true;
  }

  return false;

}

//********************************************************************
bool anaUtils::TrueParticleCrossesFGD(const AnaTrueParticleB* track, SubDetId::SubDetEnum det){
  //********************************************************************

  if(!track)
    return false;

  //for the moment consider the minimum separation in the direction perpendicular to the beam as:
  // 48 mm (iron) mm (we are intersted for tracks that come from inside,  so should cross the first iron block)
  // as it comes from oaAnalysis only one point is saved for SMRD,  as I understand the code it should be the exit point
  // so make a cut in X/Y distance between the point and particular surface inner border  

  if  (det!=SubDetId::kFGD && !SubDetId::IsFGDDetector(det))
    return false;


  //loop through det crossings
  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    AnaDetCrossingB* cross = track->DetCrossings[idet];
    if(!cross)
      continue;

    // i.e crossing the active part of the FGD
    if (!SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, det) || !track->DetCrossings[idet]->InActive)
      continue;

    //the separation should be done using the z position, since the fgd is separated by layer in z,
    //making the z position the reconstructed quantity to be cut on
    Float_t sep = fabs(track->DetCrossings[idet]->EntrancePosition[2] - track->DetCrossings[idet]->ExitPosition[2]);

    //should be at least two layers
    if(sep>DetDef::fgdXYModuleWidth)
      return true;

  }

  return false;

}

//********************************************************************
bool anaUtils::TrueParticleCrossesTPC(const AnaTrueParticleB* track, SubDetId::SubDetEnum det){
  //********************************************************************

  if(!track)
    return false;

  if  (det!=SubDetId::kTPC && !SubDetId::IsTPCDetector(det))
    return false;

  Float_t dist=-999999.;

  //loop through det crossings
  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    AnaDetCrossingB* cross = track->DetCrossings[idet];
    if(!cross)
      continue;

    // i.e crossing the active part of the FGD
    if (!SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, det) || !track->DetCrossings[idet]->InActive)
      continue;
    
    //i.e crossing the active part of the tpc
    if(SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, det) && track->DetCrossings[idet]->InActive){
        Float_t sep = anaUtils::GetSeparationSquared(track->DetCrossings[idet]->EntrancePosition, track->DetCrossings[idet]->ExitPosition);
        if(sep>dist) dist=sep;
      }
  }
  
  if(dist>62500)//bigger than the ~1/4 of the width of the TPC
      return true;
      
  return false;

}

//**************************************************
int anaUtils::GetECalDetCrossed(const AnaTrueParticleB* track, SubDetId::SubDetEnum det[]) {
  //**************************************************

  int count = 0;

  if (!track)
    return count;

  for (Int_t idet = 0; idet < track->nDetCrossings; idet++){

    if((track->DetCrossings[idet]->Detector & (1<<SubDetId::kTECAL)) && track->DetCrossings[idet]->InActive) {

      //use geometry info to split between sub-detectors
      if (track->DetCrossings[idet]->EntrancePosition[1] > 1265.        || track->DetCrossings[idet]->ExitPosition[1] > 1265)
        det[count++] = SubDetId::kTopTECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[1] < -1265.  || track->DetCrossings[idet]->ExitPosition[1] < -1265.)
        det[count++] = SubDetId::kBottomTECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[0] > 1300.   || track->DetCrossings[idet]->ExitPosition[0] > 1300.)
        det[count++] = SubDetId::kLeftTECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[0] < -1300.  || track->DetCrossings[idet]->ExitPosition[0] < -1300.) 
        det[count++] = SubDetId::kRightTECAL;

    }
    
    else if((track->DetCrossings[idet]->Detector & (1<<SubDetId::kPECAL)) && track->DetCrossings[idet]->InActive) {

      //use geometry info to split between sub-detectors
      if (track->DetCrossings[idet]->EntrancePosition[1] > 1265.        || track->DetCrossings[idet]->ExitPosition[1] > 1265)
        det[count++] = SubDetId::kTopPECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[1 ] < -1265. || track->DetCrossings[idet]->ExitPosition[1] < -1265.)
        det[count++] = SubDetId::kBottomPECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[0] > 1300.   || track->DetCrossings[idet]->ExitPosition[0] > 1300.)
        det[count++] = SubDetId::kLeftPECAL;

      else if (track->DetCrossings[idet]->EntrancePosition[0] < -1300.  || track->DetCrossings[idet]->ExitPosition[0] < -1300.) 
        det[count++] = SubDetId::kRightPECAL;
    }
    else if((track->DetCrossings[idet]->Detector & (1<<SubDetId::kDSECAL)) && track->DetCrossings[idet]->InActive) {
        det[count++] = SubDetId::kDSECAL;
    }
  }
  return count;
}



//**************************************************
int anaUtils::GetSMRDDetCrossed(const AnaTrueParticleB* track, SubDetId::SubDetEnum det[]) {
  //**************************************************
  int count = 0;

  if (!track)
    return count;

  for(Int_t idet=0;idet<track->nDetCrossings;idet++){
    //InActive flag is not set for SMRD in oaAnalysis

    if(!(track->DetCrossings[idet]->Detector & (1<<SubDetId::kSMRD))) 
      continue;

    //use geometry info to split between sub-detectors
    if (track->DetCrossings[idet]->EntrancePosition[1]>2010.)
      det[count++] = SubDetId::kTopSMRD;

    if (track->DetCrossings[idet]->EntrancePosition[1]<-2010.)
      det[count++] = SubDetId::kBottomSMRD;

    if (track->DetCrossings[idet]->EntrancePosition[0]>1830.)
      det[count++] = SubDetId::kLeftSMRD;

    if (track->DetCrossings[idet]->EntrancePosition[0]<-1830.) 
      det[count++] = SubDetId::kRightSMRD;
  }

  return count;
}

//**************************************************
int anaUtils::GetFGDDetCrossed(const AnaTrueParticleB* track, SubDetId::SubDetEnum det[]) {
  //**************************************************
  int count = 0;

  if (!track)
    return count;

  for(Int_t idet=0;idet<track->nDetCrossings;idet++){

    // i.e crossing the active part of the FGD
    if (SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, SubDetId::kFGD1) 
        && track->DetCrossings[idet]->InActive)
      det[count++] = SubDetId::kFGD1;

    if (SubDetId::GetDetectorUsed(track->DetCrossings[idet]->Detector, SubDetId::kFGD2) 
        && track->DetCrossings[idet]->InActive)
      det[count++] = SubDetId::kFGD2;

  }

  return count;
}

//**************************************************
AnaDetCrossingB* anaUtils::GetAnaDetCrossing(const AnaTrueParticleB* track, SubDetId::SubDetEnum det) {
  //**************************************************
 
  AnaDetCrossingB* cross = NULL; 
  
  
  if (!track) return cross; 
  
  for (int i = 0; i < track->nDetCrossings; i++){
    
    if (!track->DetCrossings[i]) continue;
    
    SubDetId::SubDetEnum det_bit = SubDetId::GetSubdetectorEnum(!track->DetCrossings[i]->Detector);
    
    if (det_bit == SubDetId::kInvalid) continue;
    
    if (det_bit != det) continue;
    
    if (!track->DetCrossings[i]->InActive) continue;
    
    cross = track->DetCrossings[i]; 

  }
  
  return cross;

}

//**************************************************
Float_t anaUtils::GetEntranceMomentum(const AnaDetCrossingB& cross){
  //************************************************** 
  return sqrt(pow(cross.EntranceMomentum[0],2) +
              pow(cross.EntranceMomentum[1],2) + 
              pow(cross.EntranceMomentum[2],2));  
}

//**************************************************
Float_t anaUtils::GetExitMomentum(const AnaDetCrossingB& cross){
  //**************************************************
 
  return sqrt(pow(cross.ExitMomentum[0],2) +
              pow(cross.ExitMomentum[1],2) + 
              pow(cross.ExitMomentum[2],2));  
}