numuCC4piMultiPiSelection.cxx

#include "numuCC4piMultiPiSelection.hxx"
#include "baseSelection.hxx"
#include "SystematicTuning.hxx"
#include "CutUtils.hxx"
#include "EventBoxUtils.hxx"
#include "Parameters.hxx"
#include "SubDetId.hxx"
#include "SystId.hxx"
#include "SystematicUtils.hxx"

//********************************************************************
numuCC4piMultiPiSelection::numuCC4piMultiPiSelection(bool forceBreak): SelectionBase(forceBreak,EventBoxId::kEventBoxTracker) {
    //********************************************************************
    
    // Initialize systematic tuning parameters
    systTuning::Initialize();
}

//********************************************************************
void numuCC4piMultiPiSelection::DefineSteps(){
    //********************************************************************
    
    AddStep(StepBase::kCut,       "event quality",          new EventQualityCut(),      true); //if passed accum_level=1
    AddStep(StepBase::kCut,       "> 0 tracks ",            new numuCC4pi::TotalMultiplicityCut(), true); //if passed accum_level=2
    AddStep(StepBase::kAction,    "Sort TPC tracks",        new numuCC4pi::SortTracksAction());
    AddStep(StepBase::kCut,       "quality+fiducial",       new numuCC4pi::TrackGQandFVCut(),      true); //if passed accum_level=3
    AddStep(StepBase::kAction,    "veto Action",            new numuCC4pi::VetoAction());
    AddStep(StepBase::kAction,    "muon PID Action",        new numuCC4pi::PIDAction());
    AddStep(StepBase::kAction,    "find vertex",            new numuCC4pi::FindVertexAction());
    AddStep(StepBase::kAction,    "fill summary",           new FillSummaryAction_numuCC4piMultiPi());
    AddStep(StepBase::kAction,    "find_pions",             new numuCC4piMultiPi::FindPionsAction());
    AddStep(StepBase::kAction,    "find ECal photons",      new numuCC4piMultiPi::FindEcalPhotonsAction());
        
    AddSplit(4);
    
    AddStep(0, StepBase::kCut,    "Fwd Quality Cut",        new numuCC4pi::Fwd_Quality());         //if passed accum_level=4
    AddStep(0, StepBase::kCut,    "Fwd Veto Cut",           new numuCC4pi::Fwd_Veto());            //if passed accum_level=5       
    AddStep(0, StepBase::kCut,    "Fwd PID Cut",            new numuCC4pi::Fwd_PID());             //if passed accum_level=6
    AddStep(0, StepBase::kCut,    "Fwd 4piMultiPi Cut",     new numuCC4pi::Fwd());             //if passed accum_level=7        
        
        AddSplit(3,0);
        AddStep(0,0, StepBase::kCut,  "Fwd CC-0pi",             new numuCC4piMultiPi::NoPionCut());
    AddStep(0,0, StepBase::kCut,  "Fwd ECal Pi0 veto",      new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(0,1, StepBase::kCut,  "Fwd CC-1pi",             new numuCC4piMultiPi::OnePionCut());
    AddStep(0,1, StepBase::kCut,  "Fwd ECal Pi0 veto",      new numuCC4piMultiPi::EcalPi0VetoCut());
   
    AddStep(0,2, StepBase::kCut,  "Fwd CC-Other",           new numuCC4piMultiPi::OthersCut());
        
    AddStep(1, StepBase::kCut,    "Bwd Quality Cut",        new numuCC4pi::Bwd_Quality());
    AddStep(1, StepBase::kCut,    "Bwd Veto Cut",           new numuCC4pi::Bwd_Veto());
    AddStep(1, StepBase::kCut,    "Bwd PID Cut",            new numuCC4pi::Bwd_PID());
    AddStep(1, StepBase::kCut,    "Bwd 4piMultiPi Cut",     new numuCC4pi::Bwd());
        
        AddSplit(3,1);
        AddStep(1,0, StepBase::kCut,  "Bwd CC-0pi",             new numuCC4piMultiPi::NoPionCut());
    AddStep(1,0, StepBase::kCut,  "Bwd ECal Pi0 veto",      new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(1,1, StepBase::kCut,  "Bwd CC-1pi",             new numuCC4piMultiPi::OnePionCut());
    AddStep(1,1, StepBase::kCut,  "Bwd ECal Pi0 veto",      new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(1,2, StepBase::kCut,  "Bwd CC-Other",           new numuCC4piMultiPi::OthersCut());
           
    AddStep(2, StepBase::kCut,    "HAFwd Quality Cut",      new numuCC4pi::HAFwd_Quality());
    AddStep(2, StepBase::kCut,    "HAFwd Veto Cut",         new numuCC4pi::HAFwd_Veto());
    AddStep(2, StepBase::kCut,    "HAFwd PID Cut",          new numuCC4pi::HAFwd_PID());
    AddStep(2, StepBase::kCut,    "HAFwd 4piMultiPi Cut",   new numuCC4pi::HAFwd());
        
        AddSplit(3,2);
        AddStep(2,0, StepBase::kCut,  "HAFwd CC-0pi",           new numuCC4piMultiPi::NoPionCut());
    AddStep(2,0, StepBase::kCut,  "HAFwd ECal Pi0 veto",    new numuCC4piMultiPi::EcalPi0VetoCut());
   
    AddStep(2,1, StepBase::kCut,  "HAFwd CC-1pi",           new numuCC4piMultiPi::OnePionCut());
    AddStep(2,1, StepBase::kCut,  "HAFwd ECal Pi0 veto",    new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(2,2, StepBase::kCut,  "HAFwd CC-Other",         new numuCC4piMultiPi::OthersCut());
           
    AddStep(3, StepBase::kCut,    "HABwd Quality Cut",      new numuCC4pi::HABwd_Quality());
    AddStep(3, StepBase::kCut,    "HABwd Veto Cut",         new numuCC4pi::HABwd_Veto());
    AddStep(3, StepBase::kCut,    "HABwd PID Cut",          new numuCC4pi::HABwd_PID());
    AddStep(3, StepBase::kCut,    "HABwd 4piMultiPi Cut",   new numuCC4pi::HABwd());
        
        AddSplit(3,3);
        AddStep(3,0, StepBase::kCut,  "HABwd CC-0pi",           new numuCC4piMultiPi::NoPionCut());
    AddStep(3,0, StepBase::kCut,  "HABwd ECal Pi0 veto",    new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(3,1, StepBase::kCut,  "HABwd CC-1pi",           new numuCC4piMultiPi::OnePionCut());
    AddStep(3,1, StepBase::kCut,  "HABwd ECal Pi0 veto",    new numuCC4piMultiPi::EcalPi0VetoCut());
    
        AddStep(3,2, StepBase::kCut, "HABwd CC-Other",          new numuCC4piMultiPi::OthersCut());
        
    SetBranchAlias(0, "Fwd CC0pi",    0, 0);
    SetBranchAlias(1, "Fwd CC1pi",    0, 1);
    SetBranchAlias(2, "Fwd CCOther",  0, 2);
    
        SetBranchAlias(3, "Bwd CC0pi",     1,0);
    SetBranchAlias(4, "Bwd CC1pi",     1,1);
    SetBranchAlias(5, "Bwd CCOther",   1,2);
        
        SetBranchAlias(6, "HAFwd CC0pi",   2,0);
    SetBranchAlias(7, "HAFwd CC1pi",   2,1);
    SetBranchAlias(8, "HAFwd CCOther", 2,2);

        SetBranchAlias(9, "HABwd CC0pi",   3,0);
    SetBranchAlias(10, "HABwd CC1pi",  3,1);
    SetBranchAlias(11, "HABwd CCOther",3,2);

    SetPreSelectionAccumLevel(2);
    
}

//********************************************************************
void numuCC4piMultiPiSelection::DefineDetectorFV(){
    //********************************************************************
    
    // Change FV definition to take all thickness
    // Note! this will affect all downstream stuff
    FVDef::FVdefminFGD1[2] = 0;
    FVDef::FVdefmaxFGD1[2] = 0;
    
    // The detector in which the selection is applied
    SetDetectorFV(SubDetId::kFGD1);
    
}


//**************************************************
void numuCC4piMultiPiSelection::InitializeEvent(AnaEventC& eventBB){
//**************************************************

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

        // Create the appropriate EventBox if it does not exist yet
        if (!event.EventBoxes[EventBoxId::kEventBoxTracker]) event.EventBoxes[EventBoxId::kEventBoxTracker] = new EventBoxTracker();

        boxUtils::FillTracksWithTPC(event,  static_cast<SubDetId::SubDetEnum>(GetDetectorFV()));
        boxUtils::FillTracksWithFGD(event,  static_cast<SubDetId::SubDetEnum>(GetDetectorFV()));
        boxUtils::FillTracksWithECal(event);

        boxUtils::FillTrajsChargedInTPC(event);
        boxUtils::FillTrajsChargedInFGDAndNoTPC(event,  static_cast<SubDetId::SubDetEnum>(GetDetectorFV()));
        boxUtils::FillTrajsChargedHATracker(event,  static_cast<SubDetId::SubDetEnum>(GetDetectorFV()));
        boxUtils::FillTrajsInECal(event);

}

//********************************************************************
bool numuCC4piMultiPiSelection::FillEventSummary(AnaEventC& event, Int_t allCutsPassed[]){
    //********************************************************************
    
    if(allCutsPassed[0]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC0Pi;
    if(allCutsPassed[1]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC1Pi;
    if(allCutsPassed[2]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCCOther;

    if(allCutsPassed[3]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC0Pi;
    if(allCutsPassed[4]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC1Pi;
    if(allCutsPassed[5]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCCOther;

    if(allCutsPassed[6]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC0Pi;
    if(allCutsPassed[7]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC1Pi;
    if(allCutsPassed[8]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCCOther;

    if(allCutsPassed[9]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC0Pi;
    if(allCutsPassed[10]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCC1Pi;
    if(allCutsPassed[11]) static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = SampleId::kFGD1NuMuCCOther;
        
    return (static_cast<AnaEventSummaryB*>(event.Summary)->EventSample != SampleId::kUnassigned);
}

//********************************************************************
bool FillSummaryAction_numuCC4piMultiPi::Apply(AnaEventC& event, ToyBoxB& boxB) const{
    //********************************************************************
    
    // Cast the ToyBox to the appropriate type
    ToyBoxCC4piMultiPi* cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&boxB);
    
    if(!cc4piMultiPibox->MainTrack) return 1;
    
    static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1NuMuCC0Pi] = cc4piMultiPibox->MainTrack;
    static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1NuMuCC1Pi] = cc4piMultiPibox->MainTrack;
    static_cast<AnaEventSummaryB*>(event.Summary)->LeptonCandidate[SampleId::kFGD1NuMuCCOther] = cc4piMultiPibox->MainTrack;
    
    anaUtils::CopyArray( cc4piMultiPibox->MainTrack->PositionStart, static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1NuMuCC0Pi], 4);
    anaUtils::CopyArray( cc4piMultiPibox->MainTrack->PositionStart, static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1NuMuCC1Pi], 4);
    anaUtils::CopyArray( cc4piMultiPibox->MainTrack->PositionStart, static_cast<AnaEventSummaryB*>(event.Summary)->VertexPosition[SampleId::kFGD1NuMuCCOther], 4);
    
    if(cc4piMultiPibox->MainTrack->GetTrueParticle()) static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1NuMuCC0Pi] = cc4piMultiPibox->MainTrack->GetTrueParticle()->TrueVertex;
    if(cc4piMultiPibox->MainTrack->GetTrueParticle()) static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1NuMuCC1Pi] = cc4piMultiPibox->MainTrack->GetTrueParticle()->TrueVertex;
    if(cc4piMultiPibox->MainTrack->GetTrueParticle()) static_cast<AnaEventSummaryB*>(event.Summary)->TrueVertex[SampleId::kFGD1NuMuCCOther] = cc4piMultiPibox->MainTrack->GetTrueParticle()->TrueVertex;
    
    return true;
}

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

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

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

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

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

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

    cc4piMultiPibox->nPosPions   = npospions+pionFGD;
    cc4piMultiPibox->nOtherPions = nnegpions+npi0;
    return true;
}

//*********************************************************************
bool numuCC4piMultiPi::FindEcalPhotonsAction::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 
    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&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; 
        cc4piMultiPibox->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){
        cc4piMultiPibox->HObject  = highestObject;
        // Check if the highest energetic ecal object is shower-like
        if( numuCC4piMultiPiUtils::ECALPi0Selection(event, box, highestObject, MostUpstreamLayerHitCut, det) ) {
            cc4piMultiPibox->Pi0Ecaltrack  = highestObject;
            cc4piMultiPibox->nPi0Ecaltracks = cc4piMultiPibox->nPi0Ecaltracks+1; 
        }
    }  
    return true;
}

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

    (void)event;

    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);

    // no positive and no other pions
    if( cc4piMultiPibox->nPosPions + cc4piMultiPibox->nOtherPions == 0 ) {
                return true; 
        }
    return false; 
}


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

    (void)event;

    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);
    
    // one Positive pion and no other pions
    if( cc4piMultiPibox->nOtherPions != 0 ) return false; 
    if( cc4piMultiPibox->nPosPions == 1 ){
        return true; 
    }  
    return false; 
}


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

    (void)event;

    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);

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

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

    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);

    cc4piMultiPibox->nPositivePionTPCtracks = 0;
    cc4piMultiPibox->nPosPi0TPCtracks = 0;
    cc4piMultiPibox->nNegativePionTPCtracks = 0;
    cc4piMultiPibox->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(useOldSecondaryPID){
            if ( numuCC4piMultiPiUtils::TPCpionSelection(ptrack) ) {
                cc4piMultiPibox->PositivePionTPCtracks[cc4piMultiPibox->nPositivePionTPCtracks] = ptrack;
                cc4piMultiPibox->nPositivePionTPCtracks++;
            }
            else if ( numuCC4piMultiPiUtils::TPCElPi0Selection(ptrack) ) {
                cc4piMultiPibox->PosPi0TPCtracks[cc4piMultiPibox->nPosPi0TPCtracks] = ptrack; 
                cc4piMultiPibox->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 ){
                cc4piMultiPibox->PositivePionTPCtracks[cc4piMultiPibox->nPositivePionTPCtracks] = ptrack;
                cc4piMultiPibox->nPositivePionTPCtracks++;
            }
            else {
                if( ptrack->Momentum > 900. ) continue; // This is mainly protons.
                cc4piMultiPibox->PosPi0TPCtracks[cc4piMultiPibox->nPosPi0TPCtracks] = ptrack; 
                cc4piMultiPibox->nPosPi0TPCtracks++;
            }
        }
      }
      else{
        if( cc4piMultiPibox->MainTrack == ptrack ) continue; // Same as the leading track.

        if(useOldSecondaryPID) {
            if ( numuCC4piMultiPiUtils::TPCpionSelection(ptrack) ) {
                cc4piMultiPibox->NegativePionTPCtracks[cc4piMultiPibox->nNegativePionTPCtracks] = ptrack;
                cc4piMultiPibox->nNegativePionTPCtracks++;
            } 
            else if (  numuCC4piMultiPiUtils::TPCElPi0Selection(ptrack) ) {
                cc4piMultiPibox->ElPi0TPCtracks[cc4piMultiPibox->nElPi0TPCtracks] = ptrack; 
                cc4piMultiPibox->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.
                cc4piMultiPibox->NegativePionTPCtracks[cc4piMultiPibox->nNegativePionTPCtracks] = ptrack;
                cc4piMultiPibox->nNegativePionTPCtracks++;
            }
            else{ 
                cc4piMultiPibox->ElPi0TPCtracks[cc4piMultiPibox->nElPi0TPCtracks] = ptrack; 
                cc4piMultiPibox->nElPi0TPCtracks++;
            }
        }
      }
    }
}

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

    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);
    cc4piMultiPibox->nIsoFGDPiontracks = 0;
    cc4piMultiPibox->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( numuCC4piMultiPiUtils::FGDpionSelection(track,det) ){ 
            cc4piMultiPibox->IsoFGDPiontracks[cc4piMultiPibox->nIsoFGDPiontracks] = track;
            cc4piMultiPibox->nIsoFGDPiontracks++;
        }
        else if( numuCC4piMultiPiUtils::FGDElPi0Selection(event,box,track,det) ){
            cc4piMultiPibox->IsoFGDElPi0tracks[cc4piMultiPibox->nIsoFGDElPi0tracks] = track;
            cc4piMultiPibox->nIsoFGDElPi0tracks++;
        }
    }
}

//*********************************************************************
void numuCC4piMultiPiUtils::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 numuCC4piMultiPiUtils::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 numuCC4piMultiPiUtils::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 numuCC4piMultiPiUtils::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 numuCC4piMultiPiUtils::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 

    //    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&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 numuCC4piMultiPi::EcalPi0VetoCut::Apply( AnaEventC& event, ToyBoxB& box) const{
//*********************************************************************

    (void)event;

    // Select event without pi0 track (see FindEcalPhotonsAction)
    ToyBoxCC4piMultiPi * cc4piMultiPibox = static_cast<ToyBoxCC4piMultiPi*>(&box);
    if(cc4piMultiPibox->nPi0Ecaltracks > 0) return false;
    else return true;
}

//***********************************************************************
bool numuCC4piMultiPiUtils::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 numuCC4piMultiPiSelection::IsRelevantSystematic(const AnaEventC& event, const ToyBoxB& box, SystId_h systId, Int_t branch) const{
    //**************************************************
    
    (void)event;
    (void)branch;
    (void)box;
    (void)systId;
    
    return true;
    
}

//**************************************************
bool numuCC4piMultiPiSelection::IsRelevantRecObjectForSystematic(const AnaEventC& event, AnaRecObjectC* recObj, SystId_h systId, Int_t branch) const{
    //**************************************************
    
    (void)event;
    (void)branch;
    (void)systId;
    
    AnaTrackB* track = static_cast<AnaTrackB*>(recObj);
    
    if (!track) return false;
    if (!track->TrueObject) return false;
    
    return true;
    
}

//**************************************************
bool numuCC4piMultiPiSelection::IsRelevantTrueObjectForSystematic(const AnaEventC& event, AnaTrueObjectC* trueObj, SystId_h systId, Int_t branch) const{
    //**************************************************
    
    (void)event;
    (void)branch;
    (void)systId;
    
    AnaTrueParticleB* trueTrack = static_cast<AnaTrueParticleB*>(trueObj);
    
    if (!trueTrack) return false;
    
    return true;
    
}

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

        (void)event;
        (void)branch;

        const ToyBoxCC4piMultiPi& cc4piMultiPibox = *static_cast<const ToyBoxCC4piMultiPi*>(&boxB); 

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

        if(systId == SystId::kTpcClusterEff){
                if (track->nTPCSegments > 0) {
                        AnaTPCParticleB* tpcTrack = static_cast<AnaTPCParticleB*>(anaUtils::GetSegmentWithMostNodesInClosestTpc(*track));
                        if (tpcTrack) {
                                if (tpcTrack->NNodes > 16 && tpcTrack->NNodes < 19) return true;
                        }
                }
                return false;
        }
        if(systId == SystId::kChargeIDEff){
                if (track == cc4piMultiPibox.MainTrack) return true;
                return false;
        }
        if(systId == SystId::kTpcFgdMatchEff){
                if (track == cc4piMultiPibox.MainTrack) return true;
                return false;
        }
        if(systId == SystId::kECalPID){
                if (track == cc4piMultiPibox.MainTrack) {
                        if (branch==2 || branch==3) return true;
                        else if (branch==0 || branch==4 || branch==5) {
                                if ( anaUtils::TrackUsesDet(*track,SubDetId::kTECAL) ) return true;
                                if ( anaUtils::InFiducialVolume(SubDetId::kDSECAL, track->PositionEnd, numuCC4pi_utils::_FVdefminDsECal, numuCC4pi_utils::_FVdefmaxDsECal) ) return true;
                        }
                }
                return false;
        }

        return true;

}

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

        (void)event;
        (void)branch;

        const ToyBoxCC4piMultiPi& cc4piMultiPibox = *static_cast<const ToyBoxCC4piMultiPi*>(&boxB);

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

        if ( !cc4piMultiPibox.MainTrack->GetTrueParticle() ) return false;

        if(systId == SystId::kTpcTrackEff){
                if (trueTrack->PDG == 13 && cc4piMultiPibox.MainTrack->GetTrueParticle()->PDG!=13) return true;
                if (branch!=2 && branch!=3) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kECalTrackEff){
                if (trueTrack->PDG == 13 && cc4piMultiPibox.MainTrack->GetTrueParticle()->PDG!=13) return true;
                if (branch==2 || branch==3) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kTpcP0dMatchEff){
                if (branch==1) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kTpcECalMatchEff){
                if (branch>-1 && branch<6) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kFgdECalMatchEff){
                if (branch==2 || branch==3) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kFgdECalSmrdMatchEff){
                if (branch==2 || branch==3) { if (trueTrack->ID  == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID) return true; }
                return false;
        }
        if(systId == SystId::kSIPion){
                // If this trueTrack is associated to the MainTrack
                if (trueTrack->ID == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID){
                        if (abs(cc4piMultiPibox.MainTrack->GetTrueParticle()->PDG)        == 211) return true;
                        if (abs(cc4piMultiPibox.MainTrack->GetTrueParticle()->ParentPDG)  == 211) return true;
                        if (abs(cc4piMultiPibox.MainTrack->GetTrueParticle()->GParentPDG) == 211) return true;
                        return false;
                }
                // If this trueTrack is NOT associated to the MainTrack, consider the posibility of this trueTrack to become the MainTrack and be identified as muon
                // 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 (trueTrack->PDG == 211*((Int_t)cc4piMultiPibox.MainTrack->Charge) && trueTrack->Momentum > 0.9*cc4piMultiPibox.MainTrack->Momentum) return true;
                return false;
        }
        if(systId == SystId::kSIProton){
                // If this trueTrack is associated to the MainTrack
                if (trueTrack->ID == cc4piMultiPibox.MainTrack->GetTrueParticle()->ID){
                        if (cc4piMultiPibox.MainTrack->GetTrueParticle()->PDG        == 2212) return true;
                        if (cc4piMultiPibox.MainTrack->GetTrueParticle()->ParentPDG  == 2212) return true;
                        if (cc4piMultiPibox.MainTrack->GetTrueParticle()->GParentPDG == 2212) return true;
                        return false;
                }
                // If this trueTrack is NOT associated to the MainTrack, consider the posibility of this trueTrack to become the MainTrack
                if (trueTrack->PDG == 2212 && trueTrack->Momentum > 0.9*cc4piMultiPibox.MainTrack->Momentum) return true;
                return false;
        }

        return true;

}

//********************************************************************
bool numuCC4piMultiPiSelection::CheckRedoSelection(const AnaEventC& eventC, const ToyBoxB& PreviousToyBoxB, Int_t& redoFromStep){
    //********************************************************************
    
    (void)eventC;
    (void)PreviousToyBoxB;
    (void)redoFromStep;
    
    return true;
    
}