gammaFGD2Selection.cxx

#include "gammaFGD2Selection.hxx"
#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"

//********************************************************************
gammaFGD2Selection::gammaFGD2Selection(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 gammaFGD2Selection::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

  AddStep(StepBase::kAction, "find best FGD1 e+e- pair",            new FindFGD1PairsAction_gamma() );
  AddStep(StepBase::kCut,    "Pair Veto",                           new FGD1PairCut_gamma() ); // 7  
  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);

}

//**************************************************
bool FindFGD1PairsAction_gamma::Apply(AnaEventC& eventC, ToyBoxB& boxB) const {
//**************************************************
  
  AnaEventB& event = *static_cast<AnaEventB*>(&eventC);
  ToyBoxNueCC& nuebox = *static_cast<ToyBoxNueCC*>(&boxB);

  AnaTrackB* MainTrackFGD1 = nueCCUtils::FindLeadingTrackInDet(event, SubDetId::kFGD1);

  if(!MainTrackFGD1){
    nuebox.VetoTrack = NULL;
    return true;
  }

  if(nuebox.VetoTrack)
    std::cout << "nuebox.VetoTrack already used!" << std::endl;

  nuebox.VetoTrack = nueCCUtils::FindPairTrack(event, *MainTrackFGD1, _delta_pos_max, _pullel_min,  _pullel_max);

  return true;

}

//**************************************************
bool FGD1PairCut_gamma::Apply(AnaEventC& event, ToyBoxB& boxB) const {
//**************************************************

  (void) event;
  
  ToyBoxNueCC& nuebox = *static_cast<ToyBoxNueCC*>(&boxB);
  
  if(nuebox.VetoTrack) return false;

  return true;

}

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

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

//********************************************************************
bool gammaFGD2Selection::FillEventSummary(AnaEventC& event, Int_t allCutsPassed[]){
//********************************************************************

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

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

}

//*********************************************************************
void gammaFGD2Selection::InitializeEvent(AnaEventC& eventC){
//*********************************************************************

  _nueCCSelection.InitializeEvent(eventC);

}

//*********************************************************************
bool gammaFGD2Selection::IsRelevantRecObjectForSystematic(const AnaEventC& event, AnaRecObjectC* recObj, SystId_h systId, Int_t branch) const{
//*********************************************************************

  (void)event;
  (void)branch;

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

  // 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 gammaFGD2Selection::IsRelevantTrueObjectForSystematic(const AnaEventC& eventC, AnaTrueObjectC* trueObj, SystId_h systId, Int_t branch) const{
//*********************************************************************

  return _nueCCSelection.IsRelevantTrueObjectForSystematic(eventC, trueObj, systId, branch);

}

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

  (void)event;
  (void)branch;
  
  const ToyBoxNueCC &nuebox = *static_cast<const ToyBoxNueCC*>(&boxB);

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

  // cast the box to the appropriate type

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

  return true;

}

//**************************************************
bool gammaFGD2Selection::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->TrueObject){
      // 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;      

      // Otherwise don't consider this TrueObject
      return false;
    }
    
    // Now check the paired track
    if(box.PairTrack->TrueObject){
      // At first order the inclusive selection only depends on the tracking efficiency of the electron candidate.
      if(truePart->ID == box.PairTrack->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.PairTrack->GetTrueParticle()->PDG) != 11) return true;      

      // Otherwise don't consider this TrueObject
      return false;
    }
  } 
  else if(systId == SystId::kSIPion){
    if(box.MainTrack->TrueObject && 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->TrueObject && 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->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 gammaFGD2Selection::IsRelevantSystematic(const AnaEventC& event, const ToyBoxB& boxB, SystId_h systId, Int_t branch) const{
//*********************************************************************

  (void)event;
  (void)branch;
  (void)boxB;

  if     (systId == SystId::kFgdTrackEff)   // No FGD track eff
    return false;
  else if(systId == SystId::kFgdHybridTrackEff)  // No FGD hybrid track eff
    return false;
  else if(systId == SystId::kMichelEleEff)  // No michel electron Systematic
    return false;
  else if(systId == SystId::kSIPion){ 
    // check whether should be applied what the main track is a true pion or its product
    if (!systTuning::APPLY_SYST_FINE_TUNING) return false;
  }
  else if(systId == SystId::kSIProton){  
    //check whether should be still applied when the main track is a true proton or its product: TESTING PURPOSES
    if(!systTuning::APPLY_SYST_FINE_TUNING) return false;
  }
  else if(systId == SystId::kECalTrackEff)  // No ecal related syst
    return false;
  else if(systId == SystId::kTpcECalMatchEff)  // No ecal related syst
    return false;
  else if(systId == SystId::kECalEmHipPID)  // No ecal related syst
    return false;
  else if(systId == SystId::kFGD2Shower)  // No fgd2 shower syst
    return false;
  else if(systId == SystId::kECalPID)  // No ecal related syst
    return false;
  else if(systId == SystId::kNuETPCPileUp)  // No pile-up syst 
    return true; ///Changed that because I added the TPC2 veto (PL 21/07)
  else if(systId == SystId::kNuEP0DPileUp)  // No pile-up syst 
    return false;
  else if(systId == SystId::kNuEECalPileUp)  // No pile-up syst
    return false;
  else if(systId == SystId::kPileUp)  // No pile-up syst
    return false;
  else if(systId == SystId::kOOFV) // Different oofv syst for gamma
    return false;
  
  return true;

}

//*********************************************************************
bool gammaFGD2Selection::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;
    //}
  }

  return false;
}