BFieldDistortionSystematics.cxx

#include "BFieldDistortionSystematics.hxx"
#include "ND280AnalysisUtils.hxx"
#include "EventBoxTracker.hxx"
#include "VersioningUtils.hxx"
#include <math.h>

//#define DEBUG

//********************************************************************
BFieldDistortionSystematics::BFieldDistortionSystematics():EventVariationBase(1),BinnedParams("BFieldDistortion",k1D_SYMMETRIC_NOMEAN){
//********************************************************************

  // Get the systematic source values 
  GetSigmaValueForBin(0, _mean_error);

  // This systematic uses a Uniform random throw
  _PDF = kUniform;
}

//********************************************************************
void BFieldDistortionSystematics::Apply(const ToyExperiment& toy, AnaEventC& event){
//********************************************************************

  // Get the SystBox for this event
  SystBoxB* box = GetSystBox(event);

#ifdef DEBUG 
  std::cout << "BFieldDistortionSystematics::Apply(): " << box->nRelevantRecObjects << std::endl;  
#endif

  // Loop over the relevant tracks for this systematic
  for (Int_t itrk=0;itrk<box->nRelevantRecObjects;itrk++){
#ifdef DEBUG 
    std::cout << "itrk = " << itrk << std::endl;
#endif
    AnaTrackB* track = static_cast<AnaTrackB*>(box->RelevantRecObjects[itrk]);    

    Float_t delta = 0;
    Float_t ntpcs = 0;
      // There is no more than 2 TPC segments
      // but we use the first one to correct.
      // The TPC segment
      //      AnaTPCParticleB* tpcTrack = track->TPCSegments[0];

      // Use the TPC segment with more nodes in closest TPC to the start position of the track
      AnaTPCParticleB* tpcTrack = static_cast<AnaTPCParticleB*>(anaUtils::GetSegmentWithMostNodesInClosestTpc(*track));      
      if (!tpcTrack) continue;
      Float_t gmom=track->Momentum;
      Float_t mom = tpcTrack->Momentum;
      Float_t refmom;
#ifdef DEBUG 
      std::cout<<" mom "<<mom<<std::endl;
#endif
      //the sign represents the charge in the prod6 files!
      if (versionUtils::prod6_systematics){
        refmom = fabs(tpcTrack->RefitMomentum);
        if (refmom == 99999 || !TMath::Finite(refmom)) continue;
      }else{
        refmom = tpcTrack->RefitMomentum;
        if (refmom < 0 || !TMath::Finite(refmom))      continue;
      }
      // Make sure we have a valid momentum
      if (mom < 0 || !TMath::Finite(mom))      continue;
      // the tpc mom should always be smaller
      // than the global momentum (gmom=tpcmom+fgdmom)
      if (mom > gmom || refmom > gmom) continue;
    
      delta += (refmom-mom)*_mean_error*toy.GetToyVariations(_index)->Variations[0];
    ntpcs += 1.0;      

#ifdef DEBUG 
    std::cout <<"  p0 = " << track->Momentum << std::endl;
#endif    
    
    // Apply the variation
    if(ntpcs>0) track->Momentum += delta/(ntpcs);
    
#ifdef DEBUG 
    std::cout << "p = " << track->Momentum << " delta "<<delta<<" ntpcs "<<ntpcs<<" var "<<toy.Variations[0]<< std::endl;
#endif
  }
}
  
//********************************************************************
bool BFieldDistortionSystematics::UndoSystematic(AnaEventC& event){
//********************************************************************

  // Get the SystBox for this event
  SystBoxB* box = GetSystBox(event);

  // Get the relevant tracks for this systematic
  AnaRecObjectC** tracks = box->RelevantRecObjects;

  for (Int_t itrk=0;itrk<box->nRelevantRecObjects;itrk++){
    AnaTrackB* track = static_cast<AnaTrackB*>(tracks[itrk]);
    // Go back to the corrected momentum
    track->Momentum = track->GetOriginalTrack()->Momentum;
  }

  // Don't reset the spill to corrected
  return false;
}

//********************************************************************
Int_t BFieldDistortionSystematics::GetRelevantRecObjectGroups(const SelectionBase& sel, Int_t* IDs) const{
//********************************************************************

  Int_t ngroups=0;
  for (UInt_t b=0; b<sel.GetNBranches(); b++){
    SubDetId_h det = sel.GetDetectorFV(b);
    if (det == SubDetId::kFGD1 || det == SubDetId::kFGD){
      IDs[ngroups++] = EventBoxTracker::kTracksWithTPCInFGD1FV;
    }
    if (det == SubDetId::kFGD2 || det == SubDetId::kFGD){
      IDs[ngroups++] = EventBoxTracker::kTracksWithTPCInFGD2FV;
    }
    if (det == SubDetId::kP0D){
      IDs[ngroups++] = EventBoxTracker::kTracksWithTPCInP0DFV;
    }
  }

  return ngroups;
}