antiNumuCCMultiPiSelection.cxx

#include "antiNumuCCMultiPiSelection.hxx"
#include "antiNumuCCSelection.hxx"
#include "baseSelection.hxx"
#include "CutUtils.hxx"
#include "EventBoxUtils.hxx"
#include "SubDetId.hxx"
#include "SystId.hxx"
#include "VersioningUtils.hxx"
#include "SystematicUtils.hxx"
#include "SystematicTuning.hxx"


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

}

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

  // Copy all steps from the antiNumuCCSelection
  CopySteps(_antiNumuCCSelection);

  //Pawel - Pions sele
  //Additional actions for the multi-pi selection.
  AddStep(StepBase::kAction, "fill_summary antinu_pion", new FillSummaryAction_antinuCCMultiPi());
  AddStep(StepBase::kAction, "find_pions",           new FindPionsAction_antinuCCMultiPi());
  AddStep(StepBase::kAction, "find ECal photons",    new FindEcalPhotonsAction_antinuCCMultiPi());

  //Add a split to the trunk with 3 branches.
  AddSplit(3);

  //First branch is for CC-0pi
  AddStep(0, StepBase::kCut, "CC-0pi",        new NoPionCut_antinuCCMultiPi());
  //AddStep(0, StepBase::kCut, "ECal Pi0 veto", new EcalPi0VetoCut());

  //Second branch is for CC-1pi
  AddStep(1, StepBase::kCut, "CC-1pi",        new OnePionCut_antinuCCMultiPi());
  //AddStep(1, StepBase::kCut, "ECal Pi0 veto", new EcalPi0VetoCut());

  //Third branch is for CC-Other
  AddStep(2, StepBase::kCut, "CC-Other", new OthersCut_antinuCCMultiPi());

  // Set the branch aliases to the three branches
  SetBranchAlias(0,"CC-0pi",  0);
  SetBranchAlias(1,"CC-1pi",  1);
  SetBranchAlias(2,"CC-Other",2);

  // By default the preselection correspond to cuts 0-2
  SetPreSelectionAccumLevel(2);

  // Step and Cut numbers needed by CheckRedoSelection
  _MuonPIDCutIndex     = GetCutNumber("muon PID");
  _FindPionsStepIndex  = GetStepNumber("find_pions");

}

//********************************************************************
void antiNumuCCMultiPiSelection::DefineDetectorFV(){
  //********************************************************************

  // The detector in which the selection is applied
  SetDetectorFV(SubDetId::kFGD1);
}

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

  //taken from numuCCMultiPi
  //CC0pi
  if(allCutsPassed[0]){
    static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1AntiNuMuCC0Pi;
  }  
  //CC1pi
  else if (allCutsPassed[1]){
    static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1AntiNuMuCC1Pi;
  }
  //CCOther
  else if (allCutsPassed[2]){
    static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1AntiNuMuCCOther;
  }

  // otherwise kUnassigned is used from the EventSummary constructor
  return (static_cast<AnaEventSummaryB*>(event.Summary)->EventSample != SampleId::kUnassigned);
}

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

  // Cast the ToyBox to the appropriate type
  ToyBoxTracker& box = *static_cast<ToyBoxTracker*>(&boxB); 


  if(!box.HMPtrack) return 1;
  //from numuCCMultiPi
  static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1AntiNuMuCC0Pi] = box.HMPtrack;
  static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1AntiNuMuCC1Pi] = box.HMPtrack;
  static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1AntiNuMuCCOther] = box.HMPtrack;

  for(int i = 0; i < 4; ++i){
    static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1AntiNuMuCC0Pi][i] = box.HMPtrack->PositionStart[i];
    static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1AntiNuMuCC1Pi][i] = box.HMPtrack->PositionStart[i];
    static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1AntiNuMuCCOther][i] = box.HMPtrack->PositionStart[i];
  }
  if(box.HMPtrack->GetTrueParticle()){
    static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1AntiNuMuCC0Pi] = box.HMPtrack->GetTrueParticle()->TrueVertex;
    static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1AntiNuMuCC1Pi] = box.HMPtrack->GetTrueParticle()->TrueVertex;
    static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1AntiNuMuCCOther] = box.HMPtrack->GetTrueParticle()->TrueVertex;
  }
  return 1;
}

//*********************************************************************
bool FindPionsAction_antinuCCMultiPi::Apply(AnaEventC& event, ToyBoxB& box) const{
  //*********************************************************************

  SubDetId::SubDetEnum det = static_cast<SubDetId::SubDetEnum>(box.DetectorFV);

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);
  EventBoxTracker* EventBox = static_cast<EventBoxTracker*>(event.EventBoxes[EventBoxId::kEventBoxTracker]);

  if (useTPCPions) antinumuCCMultiPiUtils::FindTPCPions(event, box, det, useOldSecondaryPID); 
  if (useME)       antinumuCCMultiPiUtils::FindMEPions(event,det, prod5Cut);
  if (useFGDPions) antinumuCCMultiPiUtils::FindIsoFGDPions(event, box, det);

  int nnegpions        =  ccmultipibox->nNegativePionTPCtracks; 
  int npospions        =  ccmultipibox->nPositivePionTPCtracks; 
  int nisofgdpions     =  ccmultipibox->nIsoFGDPiontracks;  
  int nmichelelectrons =  EventBox->nFGDMichelElectrons[det];
  int npi0 = ccmultipibox->nPosPi0TPCtracks + ccmultipibox->nElPi0TPCtracks;

  int pionFGD = 0;//nmichelelectrons;
  if (!nmichelelectrons && nisofgdpions>0) pionFGD = 1;

  /*
     ccmultipibox->nPosPions   = npospions+pionFGD;
     ccmultipibox->nOtherPions = nnegpions+npi0;
     */
  ccmultipibox->nPosPions   = npospions + nmichelelectrons;
  ccmultipibox->nNegPions = nnegpions + pionFGD;
  ccmultipibox->nOtherPions = ccmultipibox->nPosPions+npi0;

  return true;
}

//*********************************************************************
bool FindEcalPhotonsAction_antinuCCMultiPi::Apply( AnaEventC& event, ToyBoxB& box) const{
  //*********************************************************************
  // Find whether there is a pi0 track in the ECal or not 
  // In order to apply a pi0 veto cut later 
  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);
  SubDetId::SubDetEnum det = static_cast<SubDetId::SubDetEnum>(box.DetectorFV);

  // Get all Ecal objects
  EventBoxTracker::RecObjectGroupEnum groupID;
  groupID = EventBoxTracker::kTracksWithECal;
  EventBoxB* EventBox = event.EventBoxes[EventBoxId::kEventBoxTracker];

  AnaTrackB *highestObject = NULL;
  double higherEnergyObject = 0;

  //loop over all ECal objects
  for(int i = 0; i < EventBox->nRecObjectsInGroup[groupID]; i++ ) {
    AnaTrackB *allECALObjects = static_cast<AnaTrackB*>(EventBox->RecObjectsInGroup[groupID][i]);
    if(!(allECALObjects->nECALSegments)) continue;
    // Check for isolated ecal object
    if(!(anaUtils::TrackUsesOnlyDet(*allECALObjects,SubDetId::kDSECAL) ||
          anaUtils::TrackUsesOnlyDet(*allECALObjects,SubDetId::kTopTECAL) ||
          anaUtils::TrackUsesOnlyDet(*allECALObjects,SubDetId::kBottomTECAL) ||
          anaUtils::TrackUsesOnlyDet(*allECALObjects,SubDetId::kLeftTECAL) ||
          anaUtils::TrackUsesOnlyDet(*allECALObjects,SubDetId::kRightTECAL))
      )  continue;

    AnaECALParticleB* ecalBase = allECALObjects->ECALSegments[0];    
    AnaTrackB* object =  allECALObjects; 
    ccmultipibox->ISOEcal.push_back(object);   

    // Find the most energetic object in ECal 
    if(ecalBase->EMEnergy > higherEnergyObject){
      higherEnergyObject = ecalBase->EMEnergy;
      highestObject = object;
    } 
  }// end of loop on allECALTracks

  // Save the highest energetic ecal object in the box and the pi0 track
  if (highestObject){
    ccmultipibox->HObject  = highestObject;
    // Check if the highest energetic ecal object is shower-like
    if( antinumuCCMultiPiUtils::ECALPi0Selection(event, box, highestObject, MostUpstreamLayerHitCut, det) ) {
      ccmultipibox->Pi0Ecaltrack  = highestObject;
      ccmultipibox->nPi0Ecaltracks = ccmultipibox->nPi0Ecaltracks+1; 
    }
  }  
  return true;
}

//*********************************************************************
bool NoPionCut_antinuCCMultiPi::Apply(AnaEventC& event, ToyBoxB& box) const{
  //*********************************************************************

  (void)event;

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);

  // no negative[PP] and no other pions
  if( ccmultipibox->nNegPions + ccmultipibox->nOtherPions == 0 ) return true; 
  return false; 
}


//*********************************************************************
bool OnePionCut_antinuCCMultiPi::Apply(AnaEventC& event, ToyBoxB& box) const{
  //*********************************************************************

  (void)event;

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);

  // one negative pion and no other pions
  if( ccmultipibox->nOtherPions != 0 ) return false; 
  if( ccmultipibox->nNegPions == 1 )   return true; //PP
  return false; 
}


//*********************************************************************
bool OthersCut_antinuCCMultiPi::Apply(AnaEventC& event, ToyBoxB& box) const{
  //*********************************************************************

  (void)event;

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);

  // More than one positive pion or more than zero other pions or event doesn't pass the pi0 veto
  if( ccmultipibox->nOtherPions != 0 )  return true; 
  if( ccmultipibox->nNegPions > 1 )     return true; 
  if( ccmultipibox->nPi0Ecaltracks > 0) return true;

  return false; 
}

//*********************************************************************
void antinumuCCMultiPiUtils::FindTPCPions(AnaEventC& event, ToyBoxB& box, SubDetId::SubDetEnum det, bool useOldSecondaryPID){
  //*********************************************************************

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);

  ccmultipibox->nPositivePionTPCtracks = 0;
  ccmultipibox->nPosPi0TPCtracks = 0;
  ccmultipibox->nNegativePionTPCtracks = 0;
  ccmultipibox->nElPi0TPCtracks = 0;

  EventBoxTracker::RecObjectGroupEnum groupID;
  if      (det==SubDetId::kFGD1) groupID = EventBoxTracker::kTracksWithGoodQualityTPCInFGD1FV;
  else if (det==SubDetId::kFGD2) groupID = EventBoxTracker::kTracksWithGoodQualityTPCInFGD2FV;
  else return;

  EventBoxB* EventBox = event.EventBoxes[EventBoxId::kEventBoxTracker];

  // Look for pions in positive tracks 
  for(int i = 0; i < EventBox->nRecObjectsInGroup[groupID]; i++ ) {
    AnaTrackB *ptrack = static_cast<AnaTrackB*>(EventBox->RecObjectsInGroup[groupID][i]);

    if (ptrack->Charge>0){
      if( ccmultipibox->HMPtrack == ptrack ) continue; // Same as the leading track. PP - for antineutrinos - leading positive

      if(useOldSecondaryPID){
        if ( antinumuCCMultiPiUtils::TPCpionSelection(ptrack) ) {
          ccmultipibox->PositivePionTPCtracks[ccmultipibox->nPositivePionTPCtracks] = ptrack;
          ccmultipibox->nPositivePionTPCtracks++;
        }
        else if ( antinumuCCMultiPiUtils::TPCElPi0Selection(ptrack) ) {
          ccmultipibox->PosPi0TPCtracks[ccmultipibox->nPosPi0TPCtracks] = ptrack; 
          ccmultipibox->nPosPi0TPCtracks++;
        }
      }
      else {
        Float_t PIDLikelihood[4];
        anaUtils::GetPIDLikelihood(*ptrack, PIDLikelihood);

        // For Positive tracks we distinguish pions, electrons and protons.
        double ElLklh = PIDLikelihood[1];  
        double ProtonLklh = PIDLikelihood[2];  
        double PionLklh = PIDLikelihood[3];  
        double norm = ElLklh+ProtonLklh+PionLklh;
        ProtonLklh /= norm; 
        ElLklh /= norm; 
        PionLklh /= norm; 

        if( ProtonLklh > ElLklh && ProtonLklh > PionLklh ) continue; // If the highest probability is a proton continue. 

        // Id associated to the largest of the two probabilities.
        if( PionLklh > ElLklh ){
          ccmultipibox->PositivePionTPCtracks[ccmultipibox->nPositivePionTPCtracks] = ptrack;
          ccmultipibox->nPositivePionTPCtracks++;
        }
        else {
          if( ptrack->Momentum > 900. ) continue; // This is mainly protons.
          ccmultipibox->PosPi0TPCtracks[ccmultipibox->nPosPi0TPCtracks] = ptrack; 
          ccmultipibox->nPosPi0TPCtracks++;
        }
      }
    }
    else{
      //if( ccmultipibox->HMNtrack == ptrack ) continue; // Same as the leading track. (pp deleted)

      if(useOldSecondaryPID) {
        if ( antinumuCCMultiPiUtils::TPCpionSelection(ptrack) ) {
          ccmultipibox->NegativePionTPCtracks[ccmultipibox->nNegativePionTPCtracks] = ptrack;
          ccmultipibox->nNegativePionTPCtracks++;
        } 
        else if (  antinumuCCMultiPiUtils::TPCElPi0Selection(ptrack) ) {
          ccmultipibox->ElPi0TPCtracks[ccmultipibox->nElPi0TPCtracks] = ptrack; 
          ccmultipibox->nElPi0TPCtracks++;
        }
      }
      else {
        // For Negative tracks we distinguish pions, electrons and protons.
        Float_t PIDLikelihood[4];
        anaUtils::GetPIDLikelihood(*ptrack, PIDLikelihood);

        // For Positive tracks we distinguish pions, electrons and protons.
        double ElLklh = PIDLikelihood[1];  
        double PionLklh = PIDLikelihood[3];  
        double norm = ElLklh+PionLklh;
        ElLklh /= norm; 
        PionLklh /= norm;

        if( PionLklh > 0.8 ){ // Id associated to the largest of the two probabilities.
          ccmultipibox->NegativePionTPCtracks[ccmultipibox->nNegativePionTPCtracks] = ptrack;
          ccmultipibox->nNegativePionTPCtracks++;
        }
        else{ 
          ccmultipibox->ElPi0TPCtracks[ccmultipibox->nElPi0TPCtracks] = ptrack; 
          ccmultipibox->nElPi0TPCtracks++;
        }
      }
    }
  }
}


//*********************************************************************
void antinumuCCMultiPiUtils::FindIsoFGDPions(AnaEventC& event, ToyBoxB& box, SubDetId::SubDetEnum det){
  //*********************************************************************

  ToyBoxAntiCCMultiPi * ccmultipibox = static_cast<ToyBoxAntiCCMultiPi*>(&box);
  ccmultipibox->nIsoFGDPiontracks = 0;
  ccmultipibox->nIsoFGDElPi0tracks = 0;

  EventBoxTracker::RecObjectGroupEnum groupID;
  if      (det==SubDetId::kFGD1) groupID = EventBoxTracker::kTracksWithFGD1AndNoTPC;
  else if (det==SubDetId::kFGD2) groupID = EventBoxTracker::kTracksWithFGD2AndNoTPC;
  else return;

  EventBoxB* EventBox = event.EventBoxes[EventBoxId::kEventBoxTracker];

  //loop over tracks
  for (int i=0;i<EventBox->nRecObjectsInGroup[groupID]; i++ ){
    AnaTrackB* track = static_cast<AnaTrackB*>(EventBox->RecObjectsInGroup[groupID][i]);
    // Apply the fiducial cut and PID
    if( antinumuCCMultiPiUtils::FGDpionSelection(track,det) ){ 
      ccmultipibox->IsoFGDPiontracks[ccmultipibox->nIsoFGDPiontracks] = track;
      ccmultipibox->nIsoFGDPiontracks++;
    }
    else if( antinumuCCMultiPiUtils::FGDElPi0Selection(event,box,track,det) ){
      ccmultipibox->IsoFGDElPi0tracks[ccmultipibox->nIsoFGDElPi0tracks] = track;
      ccmultipibox->nIsoFGDElPi0tracks++;
    }
  }
}

//*********************************************************************
void antinumuCCMultiPiUtils::FindMEPions(AnaEventC& eventBB, SubDetId::SubDetEnum det, bool prod5Cut){
  //*********************************************************************

  AnaEventB& event = *static_cast<AnaEventB*>(&eventBB); 

  EventBoxTracker* EventBox = static_cast<EventBoxTracker*>(event.EventBoxes[EventBoxId::kEventBoxTracker]);

  // Fill the box only the first time this is used
  if (EventBox->FGDMichelElectrons[det]) return;

  anaUtils::CreateArray(EventBox->FGDMichelElectrons[det], NMAXFGDTIMEBINS);
  EventBox->nFGDMichelElectrons[det] = anaUtils::GetFGDMichelElectrons(event, det, EventBox->FGDMichelElectrons[det], prod5Cut);
  anaUtils::ResizeArray(EventBox->FGDMichelElectrons[det], EventBox->nFGDMichelElectrons[det]);

}

//*********************************************************************
bool antinumuCCMultiPiUtils::TPCpionSelection(AnaTrackB *track){
  //*********************************************************************

  Float_t PIDLikelihood[4];
  anaUtils::GetPIDLikelihood(*track, PIDLikelihood);

  if ( PIDLikelihood[3] < 0.3  ) return false; 

  double cut1 = (PIDLikelihood[0]+PIDLikelihood[3])/(1.-PIDLikelihood[2]); 

  if( track->Momentum < 500.  && cut1 < 0.8 ) return false;

  return true;


}
//*********************************************************************
bool antinumuCCMultiPiUtils::TPCElPi0Selection(AnaTrackB *track){
  //*********************************************************************

  if( track->nTPCSegments == 0 ) return false; // There is no TPC segment

  bool seltrack = false; 

  //    if (!cutUtils::TrackQualityCut(*track)) return seltrack;

  if( track->Momentum < 50. ) return seltrack; 

  //  if( track->Charge > 0. && track->Momentum > 800.) return seltrack; 

  Float_t pulls[4];

  for(int i = 0; i < track->nTPCSegments ; i++ ){

    AnaTPCParticleB *tpcTrack = track->TPCSegments[i];

    if( !tpcTrack  ) continue; // Extra protection.

    // Pulls are: Muon, Electron, Proton, Pion
    anaUtils::ComputeTPCPull(*tpcTrack,pulls);

    if (TMath::Abs(pulls[0]) > 1.e+6 // muon pull
        || TMath::Abs(pulls[1]) > 1.e+6 // electron pull
        || TMath::Abs(pulls[2]) > 1.e+6 // proton pull
        || TMath::Abs(pulls[3]) > 1.e+6) continue; // pion pull

    if( pulls[1] < -2.0 || pulls[1] > 2.0  ) break;

    if( track->Charge > 0. &&  ( pulls[2] > -4.0 && pulls[2] < 8.0 )  ) break;

    seltrack = true; 

  }
  return seltrack;
}

//*********************************************************************
bool antinumuCCMultiPiUtils::FGDpionSelection(AnaTrackB *track, SubDetId::SubDetEnum det){
  //*********************************************************************
  if( track->nFGDSegments != 1 ) return false; // There is no FGD segment or it has more than one FGD 

  AnaFGDParticleB *fgdTrack = track->FGDSegments[0];

  if( !fgdTrack  ) return false; // Extra protection. 

  if( TMath::Abs(fgdTrack->Pullp) > 1.e+6 || TMath::Abs(fgdTrack->Pullmu) > 1.e+6 || TMath::Abs(fgdTrack->Pullpi) > 1.e+6  ) return false; 

  // Check the pullnotset variable. 
  // This tells us if the pull was calculated by the algorithm. 
  // If it was not calculated this variable will be 1, so for calculated pulls fgdpulls!=1. 
  if(  fgdTrack->Pullno == 1 ) return false;

  // This tell us something about the geometry. 
  // If the algorithm says the particle started and stopped in 1st/ 2nd fgd this is equal 1 or 2. 
  //  If the particle stopped in 1st/2nd fgd, this is equal -1/-2. Other cases: 0.
  if( fgdTrack->Containment == 0 ) return false;
  if     ( det == SubDetId::kFGD1 && fgdTrack->Containment != 1 ) return false;
  else if( det == SubDetId::kFGD2 && fgdTrack->Containment != 2 ) return false; 
  else if( det == SubDetId::kFGD  && fgdTrack->Containment != 1 && fgdTrack->Containment != 2) return false; 

  double cosFGDpion = fgdTrack->DirectionStart[2];

  if(cosFGDpion > -0.3 && cosFGDpion < 0.3) return false;/// only for systematics issues. Must be remove in the future!!

  if( fgdTrack->Pullpi < 2.5 && fgdTrack->Pullpi > -2.0 ) return true; 

  return false;


}

//*********************************************************************
bool antinumuCCMultiPiUtils::FGDElPi0Selection(AnaEventC& event, ToyBoxB& box, AnaTrackB *track, SubDetId::SubDetEnum det){
  //*********************************************************************

  (void)box;

  if( track->nFGDSegments != 1 ) return false; // There is no FGD segment or it has more than one FGD 

  //    ToyBoxCCMultiPi * ccmultipibox = static_cast<ToyBoxCCMultiPi*>(&box);

  AnaFGDParticleB *fgdTrack = track->FGDSegments[0];

  if( !fgdTrack  ) return false; // Extra protection. 

  if( TMath::Abs(fgdTrack->Pullp) > 1.e+6 || TMath::Abs(fgdTrack->Pullmu) > 1.e+6 || TMath::Abs(fgdTrack->Pullpi) > 1.e+6  ) return false; 

  // Check the pullnotset variable. 
  // This tells us if the pull was calculated by the algorithm. 
  // If it was not calculated this variable will be 1, so for calculated pulls fgdpulls!=1. 
  if(  fgdTrack->Pullno == 1 ) return false;

  // This tell us something about the geometry. 
  // If the algorithm says the particle started and stopped in 1st/ 2nd fgd this is equal 1 or 2. 
  //  If the particle stopped in 1st/2nd fgd, this is equal -1/-2. Other cases: 0.
  if( fgdTrack->Containment == 0 ) return false; 
  if     ( det == SubDetId::kFGD1 && fgdTrack->Containment != 1 ) return false; 
  else if( det == SubDetId::kFGD2 && fgdTrack->Containment != 2 ) return false; 
  else if( det == SubDetId::kFGD  && fgdTrack->Containment != 1 && fgdTrack->Containment != 2) return false; 

  EventBoxTracker* EventBox = static_cast<EventBoxTracker*>(event.EventBoxes[EventBoxId::kEventBoxTracker]);

  if(EventBox->nFGDMichelElectrons[det] > 0  && fgdTrack->Pullpi < -3.0 ) return true;
  if(EventBox->nFGDMichelElectrons[det] == 0 && fgdTrack->Pullpi < -1.5 ) return true;


  return false;
}


//***********************************************************************
bool antinumuCCMultiPiUtils::ECALPi0Selection(AnaEventC& event, ToyBoxB& box, AnaTrackB *HO, int MostUpstreamLayerHitCut, SubDetId::SubDetEnum det){
  //***********************************************************************

  (void)event;
  (void)box;
  (void)det;

  // Identify if the highest ECal object comes from a TPC positive track
  // or the muon candidate.
  // If not, then the ecal object might comes from a pi0 (if the object is shower like)

  bool Ecalpi0 = false;

  double EMENERGY = -9999.;
  double MIPEM = -9999.;
  double EM_UPSTREAM = -9999;
  AnaECALParticleB* EcalSegment = HO->ECALSegments[0];
  if( !EcalSegment ) return false; 

  // Check that the Highest energetic ecal obj is actually shower-like (MIPEM>0)
  // Cut on EMEnergy because of big spike around 26 MeV (feature of the reconstruction)
  MIPEM    = EcalSegment->PIDMipEm;
  EMENERGY = EcalSegment->EMEnergy;
  if (EMENERGY > 30. && MIPEM > 0.) Ecalpi0 = true;

  // This is the cut on MostUpstreamLayerHit variable.
  // The cut threshold has to be set in the parameter file
  // For more information on this variable, see slide 10: http://www.t2k.org/nd280/physics/convener/2015/November2/ECal_pi0_veto
  EM_UPSTREAM = EcalSegment->MostUpStreamLayerHit;
  if(EM_UPSTREAM > MostUpstreamLayerHitCut ) Ecalpi0 = false;

  return Ecalpi0;
}



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

  return _antiNumuCCSelection.IsRelevantRecObjectForSystematic(event,track,systId,branch);
}

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

  return _antiNumuCCSelection.IsRelevantTrueObjectForSystematic(event,trueTrack,systId,branch);
}

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

  (void)event;
  (void)branch;
  (void)box;  

  switch(systId){
    case SystId::kFgdTrackEff: return false;
                               //kSIProton: check whether should be still applied when the main track is a true proton or its product: TESTING PURPOSES
    case SystId::kSIProton: if (!systTuning::APPLY_SYST_FINE_TUNING) return false;//new highland does not have SIProton_APPLY...
    default: return true; //the rest is included (including KFgdHybridTrackEff, kMichelEleEff)    
  }
}

//**************************************************
void antiNumuCCMultiPiSelection::InitializeEvent(AnaEventC& event){
  //**************************************************

  _antiNumuCCSelection.InitializeEvent(event);   
}

//********************************************************************
bool antiNumuCCMultiPiSelection::CheckRedoSelection(const AnaEventC& event, const ToyBoxB& PreviousToyBox, Int_t& redoFromStep){
  //********************************************************************

  // Must redo selection if antiNumuCCSelection decides so
  if( _antiNumuCCSelection.CheckRedoSelection(event,PreviousToyBox,redoFromStep)) return true;

  // Otherwise selection should not be redone since the number of tracks with TPC and FGD will not be changed by systematics
  return false;
}