numuCCBwdSelection.cxx

#include "numuCCBwdSelection.hxx"

#include "baseSelection.hxx"
#include "EventBoxUtils.hxx"
#include "VersioningUtils.hxx"
#include "SystematicUtils.hxx"


//********************************************************************
numuCCBwdSelection::numuCCBwdSelection(bool forceBreak, SubDetId::SubDetEnum aFV): SelectionBase(forceBreak,EventBoxId::kEventBoxTracker) {
//********************************************************************

  _aFV=aFV;
}

//********************************************************************
void numuCCBwdSelection::DefineSteps(){
//********************************************************************
 
  // Muon Cand
  AddStep(StepBase::kCut,       "Event quality",                                new EventQualityCut(),  true ); 
 
  AddStep(StepBase::kCut,       "Find muon candidate",                          new numuCCBwd::FindMuonCandidateAction(), true ); 
  AddStep(StepBase::kAction,    "Find vertex",                                  new numuCCBwd::FindVertexAction() );
  AddStep(StepBase::kCut,       "Muon PID",                                     new numuCCBwd::MuonPIDCut() );
  AddStep(StepBase::kCut,       "ToF is available",                             new numuCCBwd::AvailableToF() );
  AddStep(StepBase::kCut,       "More than one node in second detector",        new numuCCBwd::NoShower() );
  AddStep(StepBase::kCut,       "Positive time of flight",                      new numuCCBwd::ToF() );
  AddStep(StepBase::kCut,       "Exclusive",                                    new numuCCBwd::Exclusive() ); 
  
  //Available ToF branch
  SetBranchAlias(0, "Bwd muon");
}

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

  // The detector in which the selection is applied
  SetDetectorFV(_aFV);
}

//**************************************************
bool numuCCBwd::Exclusive::Apply(AnaEventC& event, ToyBoxB& box ) const{
//**************************************************
  
  return !( static_cast<AnaEventSummaryB*>(event.Summary)->EventSample == SampleId::kFGD1NuMuCCBwd && box.DetectorFV == SubDetId::kFGD2 );
  
}


//********************************************************************
bool numuCCBwdSelection::FillEventSummary(AnaEventC& event, Int_t allCutsPassed[]){
//********************************************************************
  
  if(allCutsPassed[0]){
    static_cast<AnaEventSummaryB*>(event.Summary)->EventSample = GetSampleEnum();
  }
  return (static_cast<AnaEventSummaryB*>(event.Summary)->EventSample != SampleId::kUnassigned);
  
}

//**************************************************
void numuCCBwdSelection::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 );
  boxUtils::FillTrajsChargedInTPC( event );
}

//********************************************************************
SubDetId::SubDetEnum numuCCBwd::InFGDFV( const Float_t pos[4] ) {
//********************************************************************
  if( anaUtils::InFiducialVolume( SubDetId::kFGD1, pos, numuCCBwd::FVdefminFGD1, numuCCBwd::FVdefmaxFGD1 ) ){
    return SubDetId::kFGD1;
  }else if( anaUtils::InFiducialVolume( SubDetId::kFGD2, pos, numuCCBwd::FVdefminFGD2, numuCCBwd::FVdefmaxFGD2 ) ){
    return SubDetId::kFGD2;
  }
  return SubDetId::kInvalid;
}

//********************************************************************
SubDetId::SubDetEnum numuCCBwd::EndsInFGDFV(const AnaTrackB& track) {
//********************************************************************
  return numuCCBwd::InFGDFV( track.PositionEnd );
}


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

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

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

        if(box.DetectorFV == SubDetId::kFGD1){
                if( track->ToF.Flag_P0D_FGD1 ){
                        return track->ToF.P0D_FGD1 > _P0D_FGD1;
                }else if(track->ToF.Flag_ECal_FGD1){
                        return track->ToF.ECal_FGD1 > _ECal_FGD1;                               
                }               
        }else if(box.DetectorFV == SubDetId::kFGD2){
                if( track->ToF.Flag_FGD1_FGD2 ){
                        return track->ToF.FGD1_FGD2 > _FGD1_FGD2;
                }else if(track->ToF.Flag_ECal_FGD2){
                        return track->ToF.ECal_FGD2 > _ECal_FGD2;                               
                }
        }

        return false;


}


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

        (void) event;

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

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

        if(box.DetectorFV == SubDetId::kFGD1){
                if( track->ToF.Flag_P0D_FGD1 && anaUtils::TrackUsesDet(*track, SubDetId::kP0D) ){
                        if( track->nP0DSegments>0 ){
                                return track->P0DSegments[0]->NNodes > 1 ;
                        }
                        std::cout<<"Flag but no segment in P0D"<<std::endl;
                }else if(track->ToF.Flag_ECal_FGD1 && anaUtils::TrackUsesDet(*track, SubDetId::kTECAL) ){
                        if( track->nECALSegments>0 ){
                                return track->ECALSegments[0]->NNodes > 1 ;
                        }
                        std::cout<<"Flag but no segment in ECal"<<std::endl;                            
                }               
        }else if(box.DetectorFV == SubDetId::kFGD2){
                if( track->ToF.Flag_FGD1_FGD2 && anaUtils::TrackUsesDet(*track, SubDetId::kFGD1) ){
                        for(int i = 0; i<track->nFGDSegments; i++){
                                if( SubDetId::GetDetectorUsed( track->FGDSegments[i]->Detector , SubDetId::kFGD1 ) ){return track->FGDSegments[i]->NNodes > 1 ;}
                        }
                        std::cout<<"Flag but no segment in FGD1"<<std::endl;
                }else if( track->ToF.Flag_ECal_FGD2 && anaUtils::TrackUsesDet(*track, SubDetId::kTECAL) ){
                        if( track->nECALSegments>0 ){
                                return track->ECALSegments[0]->NNodes > 1 ;
                        }
                        std::cout<<"Flag but no segment in ECal"<<std::endl;                            
                }
        }

        return false;

}

/// CUTS

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

        (void) event;

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

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

        bool available = false;
        if( box.DetectorFV == SubDetId::kFGD1 ){
                available = track->ToF.Flag_P0D_FGD1 || track->ToF.Flag_ECal_FGD1;
        }else if( box.DetectorFV == SubDetId::kFGD2 ){
                available = track->ToF.Flag_FGD1_FGD2 || track->ToF.Flag_ECal_FGD2;
        }

        return available;
}


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

  (void) event;

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


  AnaTrackB* track =  box.MainTrack;

  if(!track){ return false; }

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

  if (  PIDLikelihood[0] > _muon_pid ){
    return true; 
  }

  return false;
}

//Define vertex
//**************************************************
bool numuCCBwd::FindVertexAction::Apply(AnaEventC& event, ToyBoxB& boxB) const{
//**************************************************

  (void)event;

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

  // reset the vertex
  if (box.Vertex) delete box.Vertex;
  box.Vertex = NULL;

  if (! box.MainTrack) return false;

  box.Vertex = new AnaVertexB();
  anaUtils::CreateArray(box.Vertex->Particles, 1);

  box.Vertex->nParticles = 0;
  box.Vertex->Particles[box.Vertex->nParticles++] =  box.MainTrack;

  anaUtils::CopyArray( box.MainTrack->PositionEnd , box.Vertex->Position  , 4);

  if(   box.MainTrack->GetTrueParticle() ){ 
    box.TrueVertex = box.Vertex->TrueVertex =  box.MainTrack->GetTrueParticle()->TrueVertex;
  }

  return true;
}


//Modified track quality cut
//**************************************************
bool numuCCBwd::TrackQualityCut(const AnaTrackB& track, SubDetId::SubDetEnum tpc){
//**************************************************

    for(int i = 0; i < track.nTPCSegments; ++i){
        AnaTPCParticleB* tpc_track = track.TPCSegments[i];
        if( SubDetId::GetDetectorUsed( tpc_track->Detector , tpc ) ){
          if( tpc_track->NNodes > 18 ){return true;}
        }
    }

    return false;
}


bool HMFirst(AnaTrackB* a, AnaTrackB* b){//Is "a" before "b" ?
  
    AnaTPCParticleB* aTPC = NULL;
    if( numuCCBwd::EndsInFGDFV(*a) == SubDetId::kFGD1 ){
      aTPC = static_cast<AnaTPCParticleB*>( anaUtils::GetSegmentInDet( *a, SubDetId::kTPC1) );
    }else if( numuCCBwd::EndsInFGDFV(*a) == SubDetId::kFGD2 ){
      aTPC = static_cast<AnaTPCParticleB*>( anaUtils::GetSegmentInDet( *a, SubDetId::kTPC2) );
    }
    
    AnaTPCParticleB* bTPC = NULL;
    if( numuCCBwd::EndsInFGDFV(*b) == SubDetId::kFGD1 ){
      bTPC = static_cast<AnaTPCParticleB*>( anaUtils::GetSegmentInDet( *b, SubDetId::kTPC1) );
    }else if( numuCCBwd::EndsInFGDFV(*b) == SubDetId::kFGD2 ){
      bTPC = static_cast<AnaTPCParticleB*>( anaUtils::GetSegmentInDet( *b, SubDetId::kTPC2) );
    }
    
    if(!aTPC) return false;
    if(!bTPC) return true;
    
    return aTPC->Momentum > bTPC->Momentum;
    
}

//Find muon candidate
//**************************************************
bool numuCCBwd::FindMuonCandidateAction::Apply(AnaEventC& event, ToyBoxB& boxB) const{
//**************************************************

  /*
  Can't use the usual find leading track action since it requires the upstream and not the downstream end to be inside the FGD FV
  Muon candidate is taken as the highest mom "positive" track whose downstream end is inside the FGD FV, >18 TPC hits needed
  <=> HMPT FGD FV end <=> HMNT FGD FV start     
  */
  
  // Cast the ToyBox to the appropriate type
  ToyBoxTracker& box = *static_cast<ToyBoxTracker*>(&boxB); 

  box.MainTrack = NULL;

  EventBoxB* EventBox = event.EventBoxes[EventBoxId::kEventBoxTracker];
  AnaRecObjectC** selTracksNS = EventBox->RecObjectsInGroup[EventBoxTracker::kTracksWithTPC];
  int nTPC = EventBox->nRecObjectsInGroup[EventBoxTracker::kTracksWithTPC];
  
  //Sort according to TPC mom TODO : put it in a correction
  std::vector<AnaTrackB*> selTracksV;
  for (Int_t i=0;i<nTPC; ++i){
    selTracksV.push_back( static_cast<AnaTrackB*>(selTracksNS[i]) );
  }
  std::sort(selTracksV.begin(), selTracksV.end(), HMFirst);
  
  
  //Copy vector in an array
  AnaTrackB** selTracks = new AnaTrackB*[nTPC];
  std::copy(selTracksV.begin(), selTracksV.end(), selTracks);
  
  
  
  //loop over tpc tracks
  for (Int_t i=0;i<nTPC; ++i){
    AnaTrackB* track = selTracks[i];    

    if( numuCCBwd::EndsInFGDFV(*track) != box.DetectorFV ) continue;

    //apply the quality cut (>18 hits in TPC)
    if ( ( box.DetectorFV == SubDetId::kFGD1 && !numuCCBwd::TrackQualityCut(*track, SubDetId::kTPC1 ) ) 
      || ( box.DetectorFV == SubDetId::kFGD2 && !numuCCBwd::TrackQualityCut(*track, SubDetId::kTPC2 ) ) ){
      continue;
    }


    //Require a "positive" charge (implies a negative one in the backward hypothesis)
    if( track->Charge == 1 ){
      //Take the first (HM) track which satifsies the previous requirements
      box.MainTrack = track;
      
      delete[] selTracks;
      return true;
    }
  }
  
  delete[] selTracks;
  return false;
  
}

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

  (void) event;
  (void) track;
  (void) systId;
  (void) branch;
  
  return true;  
}

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

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

  return true;
}

//Redo selection ?
//********************************************************************
bool numuCCBwdSelection::CheckRedoSelection(const AnaEventC& event, const ToyBoxB& PreviousToyBox, Int_t& redoFromStep){
//********************************************************************

  (void)event;
  (void)PreviousToyBox;
  (void)redoFromStep;

  return true;
  
}