Skip to content

Commit

Permalink
Merge pull request #195 from JeffersonLab/biasingTool
Browse files Browse the repository at this point in the history
Biasing tool
  • Loading branch information
pbutti authored Apr 25, 2024
2 parents 8f94ce9 + ca4d7f7 commit db1af6c
Show file tree
Hide file tree
Showing 6 changed files with 224 additions and 5 deletions.
5 changes: 5 additions & 0 deletions processors/include/NewVertexAnaProcessor.h
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "utilities.h"

#include "TrackSmearingTool.h"
#include "TrackBiasingTool.h"
#include "FlatTupleMaker.h"
#include "AnaHelpers.h"

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -149,6 +152,8 @@ class NewVertexAnaProcessor : public Processor {
json v0proj_fits_;//!< json object v0proj
double eleTrackTimeBias_ = 0.0;
double posTrackTimeBias_ = 0.0;

double bFieldScaleFactor_ = -1;
int current_run_number_{-999}; //!< track current run number
};

Expand Down
82 changes: 77 additions & 5 deletions processors/src/NewVertexAnaProcessor.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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_);
Expand All @@ -52,6 +53,9 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) {
eleTrackTimeBias_ = parameters.getDouble("eleTrackTimeBias",eleTrackTimeBias_);
posTrackTimeBias_ = parameters.getDouble("posTrackTimeBias",posTrackTimeBias_);

//bField scale factor (ONLY to correct for mistakes in ntuplizing). If < 0 is not used
bFieldScaleFactor_ = parameters.getDouble("bFieldScaleFactor",-1);

}
catch (std::runtime_error& error)
{
Expand Down Expand Up @@ -251,6 +255,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) {
Expand Down Expand Up @@ -328,7 +337,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;
Expand All @@ -349,14 +358,42 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
if (debug_) std::cout << "got parts" << std::endl;
Track ele_trk = ele->getTrack();
Track pos_trk = pos->getTrack();

//Beam Position Corrections

if (debug_) {
std::cout<<"Check Ele/Pos Track momenta"<<std::endl;
std::cout<<ele_trk.getP()<<" "<<pos_trk.getP()<<std::endl;
std::cout<<ele_trk.getOmega()<<" "<<pos_trk.getOmega()<<std::endl;
}

// Beam Position Corrections
ele_trk.applyCorrection("z0", beamPosCorrections_.at(1));
pos_trk.applyCorrection("z0", beamPosCorrections_.at(1));
//Track Time Corrections
// Track Time Corrections
ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
pos_trk.applyCorrection("track_time", posTrackTimeBias_);

// Correct for the momentum bias

if (biasingTool_) {

// Correct for wrong track momentum - Bug Fix
// In case there was mis-configuration during reco/hpstr-ntuple step, correct
// the momentum magnitude here using the right bField for the data taking year

if (bFieldScaleFactor_ > 0) {
biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
}

biasingTool_->updateWithBiasP(ele_trk);
biasingTool_->updateWithBiasP(pos_trk);
}

if (debug_) {
std::cout<<"Corrected Ele/Pos Track momenta"<<std::endl;
std::cout<<ele_trk.getP()<<" "<<pos_trk.getP()<<std::endl;
}

double invm_smear = 1.;
if (smearingTool_) {
double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
Expand All @@ -365,6 +402,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]);
Expand Down Expand Up @@ -606,6 +644,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
//TODO Clean this up => Cuts should be implemented in each region?
//TODO Bring the preselection out of this stupid loop



if (debug_) std::cout << "start regions" << std::endl;
//TODO add yields. => Quite terrible way to loop.
Expand Down Expand Up @@ -657,7 +696,19 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
//Track Time Corrections
ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
pos_trk.applyCorrection("track_time", posTrackTimeBias_);


if (biasingTool_) {

//Correct the wrong Bfield first
if (bFieldScaleFactor_ > 0) {
biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
}

biasingTool_->updateWithBiasP(ele_trk);
biasingTool_->updateWithBiasP(pos_trk);
}

double invm_smear = 1.;
if (smearingTool_) {
double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
Expand Down Expand Up @@ -1052,6 +1103,27 @@ 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_) {

// Correct for wrong track momentum - Bug Fix
// In case there was mis-configuration during reco/hpstr-ntuple step, correct
// the momentum magnitude here using the right bField for the data taking year

if (bFieldScaleFactor_ > 0) {
biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
}


biasingTool_->updateWithBiasP(ele_trk);
biasingTool_->updateWithBiasP(pos_trk);
}


double invm_smear = 1.;
if (smearingTool_) {
double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
Expand Down
24 changes: 24 additions & 0 deletions utils/data/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Binary file added utils/data/track_pbias_eop_2019_10031.root
Binary file not shown.
47 changes: 47 additions & 0 deletions utils/include/TrackBiasingTool.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
#pragma once

//------------------//
// C++ //
//------------------//
#include <iostream>
#include <random>
#include <memory>

//------------------//
// hpstr //
//------------------//

#include "Track.h"

class TFile;
class TH1D;

class TrackBiasingTool {

public :

TrackBiasingTool(const std::string& biasingfile,
const std::string& tracks = "KalmanFullTracks");

double biasTrackP(const Track& track);

//Update the track P with a specific scale Factor
void updateWithBiasP(Track& trk, double scaleFactor);

//Update the track P with scale Factors according to the internal calibration plots
void updateWithBiasP(Track& trk);

private:

std::shared_ptr<TFile> biasingfile_;

//Biasing terms

//This is per charge -1: electron +1: positron
TH1D* eop_h_top_;
TH1D* eop_h_bot_;

// debug
bool debug_{false};

};
71 changes: 71 additions & 0 deletions utils/src/TrackBiasingTool.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
#include "TrackBiasingTool.h"
#include "TFile.h"
#include "TH1D.h"

#include <stdexcept>

TrackBiasingTool::TrackBiasingTool(const std::string& biasingfile,
const std::string& tracks) {



biasingfile_ = std::make_shared<TFile>(biasingfile.c_str());

if (!biasingfile_)
throw std::invalid_argument("Provided input biasing file doesn't exists");

eop_h_top_ = (TH1D*) biasingfile_->Get((tracks+"_eop_vs_charge_top").c_str());
eop_h_bot_ = (TH1D*) biasingfile_->Get((tracks+"_eop_vs_charge_bot").c_str());

if (!eop_h_top_ || !eop_h_bot_)
throw std::invalid_argument("Top and Bottom biasing histograms not found in smearing file");

}

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) {

double biased_p = biasTrackP(trk);

std::vector<double> momentum = trk.getMomentum();
double unbiasedp = trk.getP();

for (double& coordinate : momentum)
coordinate *= (biased_p / unbiasedp);

trk.setMomentum(momentum);

}

// This will recompute the track momentum from curvature and store it in the track
void TrackBiasingTool::updateWithBiasP(Track& trk, double scaleFactor) {

std::vector<double> momentum = trk.getMomentum();

for (double& coordinate : momentum)
coordinate *= scaleFactor ;

trk.setMomentum(momentum);

}

0 comments on commit db1af6c

Please sign in to comment.