Added Track Biasing tool and interfaced to NewVertexAnaProcessor

processors/include/NewVertexAnaProcessor.h

@@ -18,6 +18,7 @@
 #include "utilities.h"
 
 #include "TrackSmearingTool.h"
+#include "TrackBiasingTool.h"
 #include "FlatTupleMaker.h"
 #include "AnaHelpers.h"
 
@@ -117,6 +118,8 @@ class NewVertexAnaProcessor : public Processor {
 
         std::string pSmearingFile_{""};
         std::shared_ptr<TrackSmearingTool> smearingTool_;
+        std::string pBiasingFile_{""};
+        std::shared_ptr<TrackBiasingTool> biasingTool_;
 
         std::shared_ptr<TrackHistos> _vtx_histos; //!< description
         std::shared_ptr<MCAnaHistos> _mc_vtx_histos; //!< description

processors/src/NewVertexAnaProcessor.cxx

@@ -39,6 +39,7 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) {
         analysis_        = parameters.getString("analysis");
 
         pSmearingFile_ = parameters.getString("pSmearingFile",pSmearingFile_);
+	pBiasingFile_  = parameters.getString("pBiasingFile",pBiasingFile_);
 
         //region definitions
         regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_);
@@ -251,6 +252,11 @@ void NewVertexAnaProcessor::initialize(TTree* tree) {
       // just using the same seed=42 for now
       smearingTool_    =   std::make_shared<TrackSmearingTool>(pSmearingFile_,true);
     }
+
+    if (not pBiasingFile_.empty()) {
+      biasingTool_ = std::make_shared<TrackBiasingTool>(pBiasingFile_);
+    }
+
 }
 
 bool NewVertexAnaProcessor::process(IEvent* ievent) {
@@ -328,7 +334,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
     // Loop over vertices in event and make selections
     for ( int i_vtx = 0; i_vtx <  vtxs_->size(); i_vtx++ ) {
         vtxSelector->getCutFlowHisto()->Fill(0.,weight);
-
+	
         Vertex* vtx = vtxs_->at(i_vtx);
         Particle* ele = nullptr;
         Particle* pos = nullptr;
@@ -357,6 +363,13 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
         ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
         pos_trk.applyCorrection("track_time", posTrackTimeBias_);
 
+	//Correct for the momentum bias
+
+	if (biasingTool_) {
+	  biasingTool_->updateWithBiasP(ele_trk);
+          biasingTool_->updateWithBiasP(pos_trk);
+	}
+
         double invm_smear = 1.;
         if (smearingTool_) {
           double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
@@ -365,6 +378,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
           double smeared_prod = ele_trk.getP()*pos_trk.getP();
           invm_smear = sqrt(smeared_prod/unsmeared_prod);
         }
+
 
         //Add the momenta to the tracks - do not do that
         //ele_trk.setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]);
@@ -657,7 +671,12 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
             //Track Time Corrections
             ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
             pos_trk.applyCorrection("track_time", posTrackTimeBias_);
-
+
+	    if (biasingTool_) {
+	      biasingTool_->updateWithBiasP(ele_trk);
+	      biasingTool_->updateWithBiasP(pos_trk);
+	    }
+
             double invm_smear = 1.;
             if (smearingTool_) {
               double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
@@ -1052,6 +1071,16 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
             //Track Time Corrections
             ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
             pos_trk.applyCorrection("track_time", posTrackTimeBias_);
+
+
+	    // Track Momentum bias
+
+	    if (biasingTool_) {
+	      biasingTool_->updateWithBiasP(ele_trk);
+	      biasingTool_->updateWithBiasP(pos_trk);
+	    }
+
+
             double invm_smear = 1.;
             if (smearingTool_) {
               double unsmeared_prod = ele_trk.getP()*pos_trk.getP();

utils/data/README.md

@@ -22,3 +22,27 @@ The histogram storing the smearing terms as function of the number of hits is na
 ```
 python3.9 smearingPlots.py -i fee_2pt3_recon.root -m fee_2pt3_recon_mc.root
 ```
+
+
+
+
+
+## TrackBiasingTool
+
+The track biasing tool can be used to bias some of the track fit bias parameters or to evaluating tracking systematics
+on the analysis.
+
+### Workflow
+
+Currently the EoP average biases are stored in two histograms for top and bottom volumes separated by charge (e-/e+).
+They can be loaded in the trackBiasing tool that is designed to take in input a modifiable track and will perform the
+bias of (for the moment) the momentum.
+
+The TrackBiasingTool can be loaded in the anaProcessors and used to bias tracks that will then be used to form the flat ntuples
+or final plots.
+
+### Histograms
+
+The biasing tool for the moment needs two histograms:
+<Tracks>_eop_vs_charge_<volume>
+Where <Tracks> can be for example the "KalmanFullTracks" and <volume> can be "top" or "bot"

utils/include/TrackBiasingTool.h

@@ -23,12 +23,11 @@ class TrackBiasingTool {
   TrackBiasingTool(const std::string& biasingfile,
 		   const std::string& tracks = "KalmanFullTracks");
 
-  Track biasTrack(const Track& track);
+  double biasTrackP(const Track& track);
+  void updateWithBiasP(Track& trk);
 
 private:
 
-  // General Normal distributions
-
   std::shared_ptr<TFile> biasingfile_;
 
   //Biasing terms

utils/src/TrackBiasingTool.cxx

@@ -22,10 +22,38 @@ TrackBiasingTool::TrackBiasingTool(const std::string& biasingfile,
 
 }
 
+double TrackBiasingTool::biasTrackP(const Track& trk) {
+
+  double p = trk.getP();
+  double isTop = trk.getTanLambda() > 0. ? true : false;
+  double q = trk.getCharge();
+
+  TH1D* bias_histo_ = isTop ?  eop_h_top_ : eop_h_bot_;
+
+  int binN = bias_histo_->GetXaxis()->FindBin(q);
+  if (debug_)
+    std::cout<<"Track charge="<<q<<" bin="<<binN<<std::endl;
+
+  double eop_bias = bias_histo_->GetBinContent(binN);
+
+  double pcorr = p * eop_bias;
+
+  if (debug_)
+    std::cout<<"Original p = " << p <<" corrected p="<< pcorr <<std::endl;
+
+  return pcorr;
+}
+
+void TrackBiasingTool::updateWithBiasP(Track& trk) {
 
-Track TrackBiasingTool::biasTrack(const Track& track) {
+  double biased_p = biasTrackP(trk);
 
-  Track trk;
+  std::vector<double> momentum = trk.getMomentum();
+  double unbiasedp = trk.getP();
+
+  for (double& coordinate : momentum)
+    coordinate *= (biased_p / unbiasedp);
+
+  trk.setMomentum(momentum);
 
-  return trk;
 }