From da0f9146a6f95b5d7c288ac7919436b38e2548c9 Mon Sep 17 00:00:00 2001 From: pbutti Date: Mon, 22 Apr 2024 19:11:32 +0200 Subject: [PATCH 1/3] Adding TrackBiasing tool --- utils/include/TrackBiasingTool.h | 43 ++++++++++++++++++++++++++++++++ utils/src/TrackBiasingTool.cxx | 31 +++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 utils/include/TrackBiasingTool.h create mode 100644 utils/src/TrackBiasingTool.cxx diff --git a/utils/include/TrackBiasingTool.h b/utils/include/TrackBiasingTool.h new file mode 100644 index 000000000..9d888c70d --- /dev/null +++ b/utils/include/TrackBiasingTool.h @@ -0,0 +1,43 @@ +#pragma once + +//------------------// +// C++ // +//------------------// +#include +#include +#include + +//------------------// +// hpstr // +//------------------// + +#include "Track.h" + +class TFile; +class TH1D; + +class TrackBiasingTool { + + public : + + TrackBiasingTool(const std::string& biasingfile, + const std::string& tracks = "KalmanFullTracks"); + + Track biasTrack(const Track& track); + +private: + + // General Normal distributions + + std::shared_ptr biasingfile_; + + //Biasing terms + + //This is per charge -1: electron +1: positron + TH1D* eop_h_top_; + TH1D* eop_h_bot_; + + // debug + bool debug_{false}; + +}; diff --git a/utils/src/TrackBiasingTool.cxx b/utils/src/TrackBiasingTool.cxx new file mode 100644 index 000000000..f886e382e --- /dev/null +++ b/utils/src/TrackBiasingTool.cxx @@ -0,0 +1,31 @@ +#include "TrackBiasingTool.h" +#include "TFile.h" +#include "TH1D.h" + +#include + +TrackBiasingTool::TrackBiasingTool(const std::string& biasingfile, + const std::string& tracks) { + + + + biasingfile_ = std::make_shared(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"); + +} + + +Track TrackBiasingTool::biasTrack(const Track& track) { + + Track trk; + + return trk; +} From cb4c7765f0ed558d6e2dcc7a61acf90349441c63 Mon Sep 17 00:00:00 2001 From: Pierfrancesco Butti Date: Wed, 24 Apr 2024 07:27:38 -0700 Subject: [PATCH 2/3] Added Track Biasing tool and interfaced to NewVertexAnaProcessor --- processors/include/NewVertexAnaProcessor.h | 3 ++ processors/src/NewVertexAnaProcessor.cxx | 33 ++++++++++++++++++-- utils/data/README.md | 24 +++++++++++++++ utils/data/track_pbias_eop_2019_10031.root | Bin 0 -> 4308 bytes utils/include/TrackBiasingTool.h | 5 ++- utils/src/TrackBiasingTool.cxx | 34 +++++++++++++++++++-- 6 files changed, 91 insertions(+), 8 deletions(-) create mode 100644 utils/data/track_pbias_eop_2019_10031.root diff --git a/processors/include/NewVertexAnaProcessor.h b/processors/include/NewVertexAnaProcessor.h index 5cb785025..1d8510c44 100644 --- a/processors/include/NewVertexAnaProcessor.h +++ b/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 smearingTool_; + std::string pBiasingFile_{""}; + std::shared_ptr biasingTool_; std::shared_ptr _vtx_histos; //!< description std::shared_ptr _mc_vtx_histos; //!< description diff --git a/processors/src/NewVertexAnaProcessor.cxx b/processors/src/NewVertexAnaProcessor.cxx index 969095607..27611c8d5 100644 --- a/processors/src/NewVertexAnaProcessor.cxx +++ b/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(pSmearingFile_,true); } + + if (not pBiasingFile_.empty()) { + biasingTool_ = std::make_shared(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(); diff --git a/utils/data/README.md b/utils/data/README.md index afaf4bc62..700f53a69 100644 --- a/utils/data/README.md +++ b/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: +_eop_vs_charge_ +Where can be for example the "KalmanFullTracks" and can be "top" or "bot" \ No newline at end of file diff --git a/utils/data/track_pbias_eop_2019_10031.root b/utils/data/track_pbias_eop_2019_10031.root new file mode 100644 index 0000000000000000000000000000000000000000..a25b8f5bc7fc71f49338db6b8f72425fb7e94f0c GIT binary patch literal 4308 zcmdT|XHZjJw+_V6OGN1)<)aED^rF;A??{W(Py$AJ2~9wXVCbOIq$3@rNEJ~+5kYzl zNRcWCNR7ZtksJ7PXYPGxzIWz+f9{!C=j=0kueJ8`oV91xbG$t~d;x&N5C8z+2mpY` z36@0|%Mpx-V2alWHzWW6Yyki$R{%s#)(XwUO4oP7e?BXx>iD!nnEgjr5a2S94qEiM z6ac-7OvnZRfK1d;u1IQol%3Dz0OEbQoB#0&09|?!Y>MD;PB2A+U3milw7LIoOGNIM znTat3|8WA4OHA{Z{?|-2#Z|EXK+||*pt|s1g4lX{nNDoh2KOb9%hj!l1NLybI z&wsGQLYIsX7efpPAo4*V9@!!M5D0Jk1y2Nm%pP$KMz&zMuixLR*Y6`V&`U0;tvd)w zG#G$lphAfTgPOwH1_{t!C?|}3*o20JuS6iZS6P}BcmoJL(BX`u>>-YwmJtmDd;kLP zL4b_OH`Tp?Ko>Q8H77BkAv`KPDm;9T1Af0RXb%ms{KKS|RsHNlg(PfBnOTH{gk%z^ zOni+E%i~TgkTlSt9l;P0ulf$mU_*7t;b}?o+zab&N3TLT60BJtg$)ZOs(&fqiLj@F zv&VN2uVAs9Jy;4VM{w9(qW(xQHUNpxSc=R%y5N?LwakuJ9*HldZ6uH136t$ z?k-42!o%m1$Sg_#nuk!f6BgF4Z07WecMreCxNY6>3Vqe+C~&ToT~wR%^0pwg!XqbtL z*4Vv^z=fu4zaw`M=ueGE6^Wn()tEqPf(Aixke{ z_u}va+vUAcKFzhAobNVOkA~a$0zmDLNqLCnAa$a zrd*BWZseU2Q~Y=Ozo_{+5eb>xYJ6oZ&2a_aVNXgQJR}Wx6uC}AB060O^E2&He9?~) zQ;0ekeh^q+<5=sPZP;|KEQ3zk*$CZ>HjTyK3r3qHW%o4DEEm_wivNNtoh_}6$_^{I!5L)UYmYj4M7?1g9sRau#G03SO$TP+lz^X#`|;`q4RmAnCsIflH0ZrT zs)}bzR9cv-U!{Bz(93yyqWqNm

YgkkNPHmLWv!v+B^^YExR37m<}LnY~PZVr3B% za`=+Df5z{jS`r5f2P-S<@71N6H{*rO>&3hh7#REJWW#Mm7eiyxmPru~$X`U8p4RQv z$}wjyLXxYAWPXq6c}Dt@9!9p%$hA!lLxd$jgl96l|G&@bVSZInttJ|n3JlAg(#<7i2Z~I zmJ>6+G{l66jNJKrKyRbAtJ}`nxAYOy0Do&rTV~USY?0DQ>v^}Sk+N)h$Dx*8^Gh7p z@*jtL?t&Acr>z_&&i!;&Kgcd5t=Afv6Lxsub^r2bZZ}L=x5r6r|u6R&{q<&hs#s-Quch2>t;D8eg#LGoq<|mO6erq zwY37~{I~QvgoMp)ghgQ{oEW4r5jCn5ulnk_&EwEaq~cjHo%xA(WMFPrEoeyW z8hE~s^ULWPh~sf3y(p1Ka`Zt3#-2}4IZwg(KwB4fJYFlYLCyZkQw7}Db$YzaZvXO2 z@%u$(R;+*OwQWK{pt=54CrA6uII#(wIr5$}%ruyeYKMB4*W#@&OJd~`I9ykJOxU& zji)d~xF8R+FU_yjsr^cy+3d#@jTy>$^Sq_w0nY$ui-2qp3#kfODI7?lx&KyQ6i#h` zbC6=YA+@Dd@{A=YwM&yc^7nA$;xO?1!vk8DFfk6buiI9mo=D1(%+%(wPb+dlQZ9Z$ zZG|8Ba6#(&d4(dg?aBj`IkC@*B=jK_6gawBmx!#GseO2mQFH0K>K)(3G-T}rEM$WA zD}Z%_ESdM2VaCDsY^{yVS!y^@yLC(0Y5gY-y%JDJh_rL7 za|n+1gy%?P&p;lr7)LzZSdO|( zTzw7ph*hI&9Q)Hyoy5p!&nHUK>NC0}9=c-L+N#)DO$ugMI1gUM(#i9KnQSMd~o8eG467`P%^z!g{9g<#la?p{ybxg;D1{? zPYTJcG>K#tRI;U039?+*m?*wYVK?Tejio!Rbc=7L~Pc^%!q9GKb$c#A!*eAlG%NlMDJAKBlLoUr&_68ivi^N-NStz3ETH z?6c@Ra-^(IdP@6X@f--^U;f@ zHO@D4!S=sgX_0DP%gIjpp#psAq7;59LWDenG|gtrh^Chh4+EnhUm)FVJNEwz&ihi?Vaf6zG|0G4GUN zDD9%OXQVm)GhpdIa%3G@R>eIeHR?%ECNf1YZL`g2;nZ9NR*-X3pS%%~=t0!%WzsaM zzv)n#m}J1Fgx2QBlMPu~`6)K|md~p`>({sZMYsIYgN>mQ)Ou7!Tp1gc_n+SE7H@L0 zR`dmE+B{&MDxD>;vNCzwI45WX3KX8jg;ks&6W9CL5f$>xU8|&yb?f_yy6^G(y=}S2xHVtl_;XL36+cWuZPV5s7hhi%M@y zitqcG^y^-JG#$1}lc`vNM~7>UU;L&{wu*`i#^Cvo2e52nJ9+OeP)x@QeG%T=2M^8S ziV3GJh%K-GG8k~5hqX0VcAI%+OxkNVMX$e5mBrAO-u$`X`8y%Ox$L@D?RrT%e|x6p za`Qd0__u{2jg9rN!@{0Shv@`Y%UAR?OzlX7zWM&&G`7FHZgXjkmKp z#v)O8t9(s)8~ljP}NZ*wyZjX9kQ3BIr&Ee}dA0gHbe*b!maR;M+P^u7(_BoA7FIGOkiks7g&W*#6T12HE9E?R)#G0&XM*Ov>hb?ZITI#`|4xMc n$0hIICHaplJpvlJEASNm4gTuq1VP~Z9e9HN0gwAjy#)Lh-Xdw= literal 0 HcmV?d00001 diff --git a/utils/include/TrackBiasingTool.h b/utils/include/TrackBiasingTool.h index 9d888c70d..c232210a8 100644 --- a/utils/include/TrackBiasingTool.h +++ b/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 biasingfile_; //Biasing terms diff --git a/utils/src/TrackBiasingTool.cxx b/utils/src/TrackBiasingTool.cxx index f886e382e..e9254d30f 100644 --- a/utils/src/TrackBiasingTool.cxx +++ b/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="<getTrack(); Track pos_trk = pos->getTrack(); + + if (debug_) { + std::cout<<"Check Ele/Pos Track momenta"< 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"< 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. @@ -671,8 +696,15 @@ 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); } @@ -1076,6 +1108,17 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { // 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); } diff --git a/utils/include/TrackBiasingTool.h b/utils/include/TrackBiasingTool.h index c232210a8..b931da491 100644 --- a/utils/include/TrackBiasingTool.h +++ b/utils/include/TrackBiasingTool.h @@ -24,6 +24,11 @@ class TrackBiasingTool { 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: diff --git a/utils/src/TrackBiasingTool.cxx b/utils/src/TrackBiasingTool.cxx index e9254d30f..6bdb9d187 100644 --- a/utils/src/TrackBiasingTool.cxx +++ b/utils/src/TrackBiasingTool.cxx @@ -57,3 +57,15 @@ void TrackBiasingTool::updateWithBiasP(Track& trk) { 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 momentum = trk.getMomentum(); + + for (double& coordinate : momentum) + coordinate *= scaleFactor ; + + trk.setMomentum(momentum); + +}