gammaSelection.cxx

#include "gammaSelection.hxx"

#include "baseSelection.hxx"
#include "CutUtils.hxx"
#include "EventBoxUtils.hxx"
#include "VersioningUtils.hxx"
#include "SystematicUtils.hxx"
#include "trackerSelectionUtils.hxx"
#include "SystematicTuning.hxx"

#include "nueCCUtils.hxx"
#include "nueCutUtils.hxx"

//********************************************************************
gammaSelection::gammaSelection(bool forceBreak): SelectionBase(forceBreak, EventBoxId::kEventBoxTracker) {
//********************************************************************

  _nueCCSelection.Initialize();

  //_elepullpri_min = (Float_t)ND::params().GetParameterD("gammaAnalysis.Cuts.PID.PullElecMin");
  //_elepullpri_max = (Float_t)ND::params().GetParameterD("gammaAnalysis.Cuts.PID.PullElecMax");
}

//********************************************************************
void gammaSelection::DefineSteps(){
//********************************************************************

  // last "true" means the step sequence is broken if cut is not passed (default is "false")
  AddStep(StepBase::kCut,    "Event quality",                       new EventQualityCut(), true); // 0
  AddStep(StepBase::kCut,    "> 0 tracks ",                         new MultiplicityCut(), true); // 1
  AddStep(StepBase::kAction, "Find leading tracks",                 new FindLeadingTracksAction_gamma() );
  AddStep(StepBase::kAction, "Find vertex",                         new FillVertex() );
  AddStep(StepBase::kAction, "Fill summary",                        new FillSummaryActionGamma() );

  AddStep(StepBase::kCut,    "# TPC nodes 1ary track",              new PrimaryTrackQualityCut_gamma() ); // 2

  AddStep(StepBase::kCut,    "Highest Momentum Track Momentum cut", new HighestMomentumCut() ); // 3
  AddStep(StepBase::kCut,    "Electron PID",                        new PIDCut_gamma() ); // 4
  AddStep(StepBase::kAction, "Find secondary tracks",               new FindPairsAction() );
  AddStep(StepBase::kCut,    "Secondary Track",                     new PairTrackCut_gamma() ); // 5
  AddStep(StepBase::kCut,    "Invariant Mass Cut",                  new MinvCut_gamma() ); // 6

  SetBranchAlias(0, "Gamma Selection");

  _FindLeadingTracksStepIndex = GetStepNumber("Find leading tracks");
  _TotalMultiplicityCutIndex  = GetCutNumber("> 0 tracks ");

  _ElecPIDCutIndex            = GetCutNumber("Electron PID");
  _ElecPIDStepIndex           = GetStepNumber("Electron PID");

  // By default the preselection correspond to cuts 0-3
  // It means that if any of the four first cuts (0,1,2,3) is not passed
  // the loop over toys will be broken ===> A single toy will be run
  SetPreSelectionAccumLevel(3);

}

//********************************************************************
void gammaSelection::DefineDetectorFV(){
//********************************************************************
  
  // The detector in which the selection is applied
  SetDetectorFV(SubDetId::kFGD1);

  // Set the detector field into the selection being used
  _nueCCSelection.SetDetectorFV(SubDetId::kFGD1);

}

//********************************************************************
bool FillSummaryActionGamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//********************************************************************

  ToyBoxNueCC& box = *static_cast<ToyBoxNueCC*>(&boxB);

  if(!box.MainTrack) return false;

  static_cast<AnaEventSummaryB*>(eventC.Summary)->LeptonCandidate[SampleId::kFGD1Gamma] = box.MainTrack;

  for(int i = 0; i < 4; ++i){
    static_cast<AnaEventSummaryB*>(eventC.Summary)->VertexPosition[SampleId::kFGD1Gamma][i] = box.MainTrack->PositionStart[i];
  }

  if(box.MainTrack->GetTrueParticle()){
    static_cast<AnaEventSummaryB*>(eventC.Summary)->TrueVertex[SampleId::kFGD1Gamma] = box.MainTrack->GetTrueParticle()->TrueVertex;
  }

  return true;

}

//********************************************************************
bool gammaSelection::FillEventSummary(AnaEventC& eventC, Int_t allCutsPassed[]){
//********************************************************************

  // The event sample corresponding to this selection
  if(allCutsPassed[0])
    static_cast<AnaEventSummaryB*>(eventC.Summary)->EventSample = SampleId::kFGD1Gamma;

  return (static_cast<AnaEventSummaryB*>(eventC.Summary)->EventSample != SampleId::kUnassigned);

}

//*********************************************************************
bool FindLeadingTracksAction_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//*********************************************************************

  ToyBoxNueCC& box = *static_cast<ToyBoxNueCC*>(&boxB);

  trackerSelUtils::FindLeadingTracks(eventC, boxB);

  if( box.HMtrack) box.MainTrack = box.HMtrack;
  //if(!box.MainTrack) return false;
  
  return true;
  
}

//*********************************************************************
bool PairTrackCut_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//*********************************************************************

  (void) eventC;

  ToyBoxNueCC& nuebox   = *static_cast<ToyBoxNueCC*>(&boxB);
  bool accept = (nuebox.PairTrack);

  return accept;

}

//*********************************************************************
bool PrimaryTrackQualityCut_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//*********************************************************************

  (void) eventC;

  ToyBoxNueCC& box = *static_cast<ToyBoxNueCC*>(&boxB);

  // No vertex - reject
  if(!box.Vertex) return false;

  AnaTrackB* track = box.MainTrack;
  if(!track) return false;

  return nueCutUtils::TrackQuality(track, _num_tpc_nodes, _num_tpc_nodes);

}


//*********************************************************************
bool MinvCut_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//*********************************************************************

  (void) eventC;

  ToyBoxNueCC& nuebox   = *static_cast<ToyBoxNueCC*>(&boxB);

  if (!nuebox.MainTrack)
    return false;

  if (!nuebox.PairTrack)
    return false;

  return (nueCCUtils::GetInvMass(*nuebox.MainTrack, *nuebox.PairTrack) < _inv_mass_min);
  
}

//*********************************************************************
bool PIDCut_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const{
//*********************************************************************

  (void) eventC;
  ToyBoxNueCC& box = *static_cast<ToyBoxNueCC*>(&boxB);

  AnaTrackB* track_pri = box.MainTrack;
  //  AnaTrackB* track_sec = box.PairTrack;

  bool accept = nueCutUtils::TPCElectronPull(track_pri, _elepullpri_min, _elepullpri_max,
                                             _elepullpri_min, _elepullpri_max, 0.0);
    
  return accept;

}



//*********************************************************************
void gammaSelection::InitializeEvent(AnaEventC& eventC){
//*********************************************************************
  
  _nueCCSelection.InitializeEvent(eventC);
}

//*********************************************************************
bool gammaSelection::CheckRedoSelection(const AnaEventC& eventC, const ToyBoxB& PreviousToyBoxB, Int_t& redoFromStep){
//*********************************************************************

  const AnaEventB& event = *static_cast<const AnaEventB*>(&eventC); 

  // cast the box to the appropriate type
  const ToyBoxNueCC& PreviousToyBox = *static_cast<const ToyBoxNueCC*>(&PreviousToyBoxB);

  if(!PreviousToyBox.HMtrack) return false;

  // TODO: The code below implies calling cutUtils::FindLeadingTracks(event, box) twice. Can we avoid that ?
  if(PreviousToyBox.MaxAccumLevel > _TotalMultiplicityCutIndex){
    ToyBoxNueCC box;
    trackerSelUtils::FindLeadingTracks(event, box);
    
    // Redo the selection if any of the leading tracks changes identity
    if(PreviousToyBox.SHMNtrack != box.SHMNtrack  || 
       PreviousToyBox.SHMPtrack != box.SHMPtrack  ||
       PreviousToyBox.HMNtrack  != box.HMNtrack   ||
       PreviousToyBox.HMPtrack  != box.HMPtrack   ||
       PreviousToyBox.SHMtrack  != box.SHMtrack   ||
       PreviousToyBox.HMtrack   != box.HMtrack){
      
      redoFromStep = _FindLeadingTracksStepIndex;
      return true;
    }
  }

  if(PreviousToyBox.MaxAccumLevel >= _ElecPIDCutIndex){
    //bool pidCut = nueCutUtils::TPCElectronPull(PreviousToyBox.HMtrack, _elepullpri_min, _elepullpri_max, _elepullpri_min, _elepullpri_max, 0.0);
    //if(( pidCut && (PreviousToyBox.AccumLevel[0] == _ElecPIDCutIndex)) ||
    // (!pidCut && (PreviousToyBox.AccumLevel[0]  > _ElecPIDCutIndex))){
    // We should redo the selection starting from the PIDCut step
    redoFromStep = _ElecPIDStepIndex;
    return true;
    //}
  }
    
    /*
    AnaTrackB* selTracks[NMAXPARTICLES];
    int nTPC = anaUtils::GetAllTracksUsingTPC(event, selTracks);
    double distance_min, distance;
    for (int i = 0; i < nTPC; i++){
      AnaTrackB* secondary = selTracks[i];

      if(secondary == box.SHMNtrack ||
         secondary == box.SHMPtrack ||
         secondary == box.HMNtrack  ||
         secondary == box.HMPtrack  ||
         secondary == box.SHMtrack  ||
         secondary == box.HMtrack)
        continue;
    
      if(!cutUtils::FiducialCut(*secondary, static_cast<SubDetId::SubDetEnum>(box.DetectorFV)))
        continue;
    
      if(box.MainTrack->Charge * secondary->Charge > 0)
        continue;

      TVector3 PrimaryPos = anaUtils::ArrayToTVector3(box.MainTrack->PositionStart);
      TVector3 SecondaryPos = anaUtils::ArrayToTVector3(secondary->PositionStart);

      //Calculate distance between tracks
      distance = (PrimaryPos - SecondaryPos).Mag();
      if(distance < distance_min) {
        nuebox.PairTrack = secondary; //Fill secondary track box with closer track
        distance_min = distance;
      }
    }
    if(PreviousNuebox.PairTrack != nuebox.PairTrack){
      redoFromStep = _FindPairTracksStepIndex;
      return true;
    }
  }
    */

  return false;

}


//*********************************************************************
bool gammaSelection::IsRelevantRecObjectForSystematic(const AnaEventC& event, AnaRecObjectC* trackRecC, SystId_h systId, Int_t branch) const{
//*********************************************************************
  
  (void)event;
  (void)branch;

  AnaTrackB* track = static_cast<AnaTrackB*>(trackRecC);

  // True track should always exist
  if(!track->TrueObject) return false;

  if(systId == SystId::kTpcClusterEff){
    if(track->nTPCSegments == 0) return false;
    //    AnaTPCParticleB* tpcTrack = static_cast<AnaTPCParticleB*>(track->TPCSegments[0]);
    AnaTPCParticleB* tpcTrack = static_cast<AnaTPCParticleB*>(anaUtils::GetSegmentWithMostNodesInClosestTpc(*track));      
    if(!tpcTrack) return false;

    if(tpcTrack->NNodes == 17 || tpcTrack->NNodes == 18)
      return true;

    return false;
  }

  return true;
  
}


//*********************************************************************
bool gammaSelection::IsRelevantTrueObjectForSystematic(const AnaEventC& event, AnaTrueObjectC* trueTrack, SystId_h systId, Int_t branch) const{
//*********************************************************************

  return _nueCCSelection.IsRelevantTrueObjectForSystematic(event, trueTrack, systId, branch);

}


//*********************************************************************
bool gammaSelection::IsRelevantRecObjectForSystematicInToy(const AnaEventC& event, const ToyBoxB& boxB, AnaRecObjectC* recObj, SystId_h systId, Int_t branch) const{
//*********************************************************************

  (void)event;
  (void)branch;
  
  // cast the box to the appropriate type
  const ToyBoxNueCC &nuebox = *static_cast<const ToyBoxNueCC*>(&boxB);

  AnaTrackB* track = static_cast<AnaTrackB*>(recObj);

  if(systId == SystId::kChargeIDEff){
    if (track != nuebox.MainTrack && track != nuebox.PairTrack) return false;
  }

  return true;

}


//*********************************************************************
bool gammaSelection::IsRelevantTrueObjectForSystematicInToy(const AnaEventC& event, const ToyBoxB& boxB, AnaTrueObjectC* trueObj, SystId_h systId, Int_t branch) const{
//*********************************************************************

  (void)event;
  (void)branch;

  // cast the box to the appropriate type
  const ToyBoxNueCC &box = *static_cast<const ToyBoxNueCC*>(&boxB);

  AnaTrueParticleB* truePart = static_cast<AnaTrueParticleB*>(trueObj);

  if(systId == SystId::kTpcTrackEff){   
    if(box.MainTrack->GetTrueParticle()){
      // At first order the inclusive selection only depends on the tracking efficiency of the electron candidate.
      if(truePart->ID  == box.MainTrack->GetTrueParticle()->ID) return true; 
      // Consider also the case in which the electron candidate is not a true electron but this track it is
      if(abs(truePart->PDG) == 11 && abs(box.MainTrack->GetTrueParticle()->PDG) != 11) return true; 

      // Now check the paired track
      if(box.PairTrack && box.PairTrack->GetTrueParticle()){
        if(truePart->ID == box.PairTrack->GetTrueParticle()->ID) return true;
        
        if(abs(truePart->PDG) == 11 && abs(box.PairTrack->GetTrueParticle()->PDG) != 11) return true;
      }
    }
    // Apply for electrons(?)
    else if (abs(truePart->PDG) == 11)  return true;

    return false;
  } 
  else if(systId == SystId::kSIPion){
    if(box.MainTrack->GetTrueParticle() && truePart->ID == box.MainTrack->GetTrueParticle()->ID){
      if(abs(box.MainTrack->GetTrueParticle()->PDG)        == 211) return true;
      if(abs(box.MainTrack->GetTrueParticle()->ParentPDG)  == 211) return true;
      if(abs(box.MainTrack->GetTrueParticle()->GParentPDG) == 211) return true;
      return false;
    }
    // If this truePart is NOT associated to the MainTrack, consider the posibility of this truePart to become the MainTrack and be identified as electron
    // For the moment assume a negative pion may become the MainTrack if its momentum its above 90% of the current MainTrack momentum. 
    // Ideally we should check that this true pion is not associated to any reconstructed track, but this is not possible now without looping.
    // Multiply by the charge of the MainTrack such that it can be use for antiNumu selection
    if (truePart->PDG == 211*((Int_t)box.MainTrack->Charge) && truePart->Momentum > 0.9*box.MainTrack->Momentum) return true;
 
    // Now check the paired track - it doesn't have to be the most energetic track
    if(box.PairTrack && box.PairTrack->GetTrueParticle() && truePart->ID == box.PairTrack->GetTrueParticle()->ID){
      if(abs(box.PairTrack->GetTrueParticle()->PDG)        == 211) return true;
      if(abs(box.PairTrack->GetTrueParticle()->ParentPDG)  == 211) return true;
      if(abs(box.PairTrack->GetTrueParticle()->GParentPDG) == 211) return true;
      
      return false;
    }
  }
  else if(systId == SystId::kSIProton){
    // If this truePart is associated to the MainTrack
    if(box.MainTrack->GetTrueParticle() && truePart->ID == box.MainTrack->GetTrueParticle()->ID){
      if(box.MainTrack->GetTrueParticle()->PDG        == 2212) return true;
      if(box.MainTrack->GetTrueParticle()->ParentPDG  == 2212) return true;
      if(box.MainTrack->GetTrueParticle()->GParentPDG == 2212) return true;
      return false;
    }
    // If this truePart is NOT associated to the MainTrack, consider the posibility of this truePart to become the MainTrack
    //if(truePart->PDG == 2212 && truePart->Momentum > 0.9*box.MainTrack->Momentum) return true;
 
    // Now check the paired track - it doesn't have to be the most energetic track
    if(box.PairTrack && box.PairTrack->GetTrueParticle() && truePart->ID == box.PairTrack->GetTrueParticle()->ID){
      if(box.PairTrack->GetTrueParticle()->PDG        == 2212) return true;
      if(box.PairTrack->GetTrueParticle()->ParentPDG  == 2212) return true;
      if(box.PairTrack->GetTrueParticle()->GParentPDG == 2212) return true;
      return false;
    }
    
    return false;
  }

  return true;

}


//*********************************************************************
bool gammaSelection::IsRelevantSystematic(const AnaEventC& event, const ToyBoxB& box, SystId_h systId, Int_t branch) const{
//*********************************************************************

  if     (systId == SystId::kTpcECalMatchEff)  // No ecal related syst
    return false;
  else if(systId == SystId::kECalPID)  // No ecal related syst
    return false;
  else if(systId == SystId::kECalEMResol)  // No ecal related syst
    return false;
  else if(systId == SystId::kECalEMScale)  // No ecal related syst
    return false;
  else if(systId == SystId::kNuETPCPileUp)  // No pile-up syst 
    return false;
  else if(systId == SystId::kNuEP0DPileUp)  // No pile-up syst 
    return false;
  else if(systId == SystId::kNuEECalPileUp)  // No pile-up syst
    return false;

  // The rest are disabled here
  return _nueCCSelection.IsRelevantSystematic(event, box, systId, branch);
}