diff --git a/CMakeLists.txt b/CMakeLists.txt index 4d701c5cb..f5b2d5b8d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -65,7 +65,7 @@ include(MacroModule) # import macro for declaring external dependencies include(MacroExtDeps) -set(MODULES event analysis processing processors) +set(MODULES event utils analysis processing processors) # build each module in the list foreach(module ${MODULES}) diff --git a/analysis/plotconfigs/fee_smearing/feeSmearing.json b/analysis/plotconfigs/fee_smearing/feeSmearing.json new file mode 100644 index 000000000..57bf7657d --- /dev/null +++ b/analysis/plotconfigs/fee_smearing/feeSmearing.json @@ -0,0 +1,262 @@ +{ + "n_tracks_h" : { + "bins" : 10, + "minX" : 0, + "maxX" : 10, + "xtitle" : "N_{tracks}", + "ytitle" : "Events" + }, + "nHits_2d_h" : { + "bins" : 15, + "minX" : 0, + "maxX" : 15, + "xtitle" : "N_{2Dhits}", + "ytitle" : "Tracks" + }, + "nShared_h" : { + "bins" : 8, + "minX" : -0.5, + "maxX" : 7.5, + "xtitle" : "N_{shared}^{hits}", + "ytitle" : "Tracks" + }, + "d0_h" : { + "bins" : 200, + "minX" : -10, + "maxX" : 10, + "xtitle" : "d_{0} [mm]", + "ytitle" : "Tracks" + }, + "Phi_h" : { + "bins" : 100, + "minX" : -0.5, + "maxX" : 0.5, + "xtitle" : "phi_{0}", + "ytitle" : "Tracks" + }, + "Omega_h" : { + "bins" : 100, + "minX" : -0.002, + "maxX" : 0.002, + "xtitle" : "#omega", + "ytitle" : "Tracks" + }, + "TanLambda_h" : { + "bins" : 200, + "minX" : -0.3, + "maxX" : 0.3, + "xtitle" : "tan(#lambda)", + "ytitle" : "Tracks" + }, + "Z0_h" : { + "bins" : 200, + "minX" : -5, + "maxX" : 5, + "xtitle" : "z_0 [mm]", + "ytitle" : "Tracks" + }, + "time_h" : { + "bins" : 200, + "minX" : -40, + "maxX" : 40, + "xtitle" : "track time [ns]", + "ytitle" : "Tracks" + }, + "chi2_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 100, + "xtitle" : "track #chi^{2}", + "ytitle" : "Tracks" + }, + "p_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 4, + "xtitle" : "p_{e^{-}} [GeV]", + "ytitle" : "Tracks" + }, + "p_vs_nHits_hh" : { + "binsX" : 5, + "minX" : 8, + "maxX" : 13, + "binsY" : 200, + "minY" : 0, + "maxY" : 4, + "xtitle" : "Nhits", + "ytitle" : "p_{e^{-}} [GeV]" + }, + "p_vs_nHits_top_hh" : { + "binsX" : 5, + "minX" : 8, + "maxX" : 13, + "binsY" : 200, + "minY" : 0, + "maxY" : 4, + "xtitle" : "Nhits", + "ytitle" : "top p_{e^{-}} [GeV]" + }, + "p_vs_nHits_bot_hh" : { + "binsX" : 5, + "minX" : 8, + "maxX" : 13, + "binsY" : 200, + "minY" : 0, + "maxY" : 4, + "xtitle" : "Nhits", + "ytitle" : "bot p_{e^{-}} [GeV]" + }, + "p_vs_TanLambda_hh" : { + "binsX" : 100, + "minX" : -0.1, + "maxX" : 0.1, + "binsY" : 200, + "minY" : 0, + "maxY" : 4, + "xtitle" : "tan(#lambda)", + "ytitle" : "p_{e^{-}} [GeV]" + }, + "p_vs_Phi_hh" : { + "binsX" : 100, + "minX" : -0.2, + "maxX" : 0.2, + "binsY" : 200, + "minY" : 0, + "maxY" : 4, + "xtitle" : "#phi_{0}", + "ytitle" : "p_{e^{-}} [GeV]" + }, + "p_vs_TanLambda_Phi_hhh" : { + "binsX" : 100, + "minX" : -0.2, + "maxX" : 0.2, + "binsY" : 100, + "minY" : -0.1, + "maxY" : 0.1, + "binsZ" : 200, + "minZ" : 0, + "maxZ" : 4, + "xtitle" : "#phi_{0}", + "ytitle" : "tan(#lambda)", + "ztitle" : "p_{e^{-}} [GeV]" + }, + "p_vs_TanLambda_nHits_hhh" : { + "binsX" : 100, + "minX" : -0.2, + "maxX" : 0.2, + "binsY" : 5, + "minY" : 8, + "maxY" : 13, + "binsZ" : 200, + "minZ" : 0, + "maxZ" : 4, + "xtitle" : "#phi_{0}", + "ytitle" : "tan(#lambda)", + "ztitle" : "p_{e^{-}} [GeV]" + }, + "chi2ndf_h" : { + "bins" : 200, + "minX" : 0, + "maxX" : 30, + "xtitle" : "track #chi^{2} / ndf", + "ytitle" : "Tracks" + }, + "d0_vs_p_hh" : { + "binsX" : 100, + "minX" : 0, + "maxX" : 4, + "binsY" : 200, + "minY" : -10, + "maxY" : 10, + "xtitle" : "trk p [GeV]", + "ytitle" : "trk d0 [mm]" + }, + "d0_vs_phi0_hh" : { + "binsX" : 100, + "minX" : -0.4, + "maxX" : 0.4, + "binsY" : 200, + "minY" : -10, + "maxY" : 10, + "xtitle" : "trk #phi_{0}", + "ytitle" : "trk d0 [mm]" + }, + + "d0_vs_tanlambda_hh" : { + "binsX" : 200, + "minX" : -0.1, + "maxX" : 0.1, + "binsY" : 200, + "minY" : -10, + "maxY" : 10, + "xtitle" : "trk(#lamdba)", + "ytitle" : "trk d0 [mm]" + }, + + "z0_vs_p_hh" : { + "binsX" : 100, + "minX" : 0, + "maxX" : 4, + "binsY" : 200, + "minY" : -5, + "maxY" : 5, + "xtitle" : "trk p [GeV]", + "ytitle" : "trk z0 [mm]" + }, + + "z0_vs_phi0_hh" : { + "binsX" : 100, + "minX" : -0.4, + "maxX" : 0.4, + "binsY" : 200, + "minY" : -5, + "maxY" : 5, + "xtitle" : "trk #phi_{0}", + "ytitle" : "trk z0 [mm]" + }, + + "z0_vs_tanlambda_hh" : { + "binsX" : 200, + "minX" : -0.1, + "maxX" : 0.1, + "binsY" : 200, + "minY" : -5, + "maxY" : 5, + "xtitle" : "trk(#lamdba)", + "ytitle" : "trk z0 [mm]" + }, + "tanlambda_vs_phi0_hh" : { + "binsX" : 100, + "minX" : -0.4, + "maxX" : 0.4, + "binsY" : 100, + "minY" : -0.2, + "maxY" : 0.2, + "xtitle" : "trk #phi_{0}", + "ytitle" : "trk tan(#lambda)" + }, + "xpos_at_ecal_h" : { + "bins" : 400, + "minX" : -100, + "maxX" : 100, + "xtitle" : "track x pos [mm]", + "ytitle" : "Tracks" + }, + "ypos_at_ecal_h" : { + "bins" : 400, + "minX" : -100, + "maxX" : 100, + "xtitle" : "track y pos [mm]", + "ytitle" : "Tracks" + }, + "xypos_at_ecal_hh" : { + "binsX" : 100, + "minX" : -200, + "maxX" : 200, + "binsY" : 100, + "minY" : -200, + "maxY" : 200, + "xtitle" : "trk ecal x [mm]", + "ytitle" : "trk ecal y [mm]" + } +} diff --git a/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json b/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json index b107e45b8..47744a905 100644 --- a/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json +++ b/analysis/plotconfigs/tracking/simps/vtxAnalysis_2016_simp_reach.json @@ -1012,6 +1012,13 @@ "xtitle" : "vtx E_{sum} [GeV]", "ytitle" : "Vertices" }, + "vtx_smear_InvM_h" : { + "bins" : 200, + "minX" : 0.0, + "maxX" : 0.2, + "xtitle" : "M_{vtx} [GeV]", + "ytitle" : "Vertices" + }, "vtx_InvM_vtx_z_hh" : { "binsX" : 200, "minX" : 0.0, diff --git a/event/include/Track.h b/event/include/Track.h index 685bfebf2..af7ab1e03 100644 --- a/event/include/Track.h +++ b/event/include/Track.h @@ -288,9 +288,9 @@ class Track : public TObject { * @return momentum magnitude */ - double getP(){return sqrt(px_*px_ + py_*py_ + pz_*pz_);}; + double getP() const {return sqrt(px_*px_ + py_*py_ + pz_*pz_);}; - double getPt() {return sqrt(px_*px_ + pz_*pz_);} + double getPt() const {return sqrt(px_*px_ + py_*py_);} /** * Set the lambda kink of the given layer. diff --git a/processors/CMakeLists.txt b/processors/CMakeLists.txt index 96c41806d..1df85e56c 100644 --- a/processors/CMakeLists.txt +++ b/processors/CMakeLists.txt @@ -2,6 +2,6 @@ # Declare processors module module( NAME processors - DEPENDENCIES event processing analysis + DEPENDENCIES event processing analysis utils EXTERNAL_DEPENDENCIES ROOT LCIO ) diff --git a/processors/config/anaFeeSmearing_cfg.py b/processors/config/anaFeeSmearing_cfg.py new file mode 100644 index 000000000..e045112f4 --- /dev/null +++ b/processors/config/anaFeeSmearing_cfg.py @@ -0,0 +1,47 @@ +import HpstrConf +import os +import sys +import baseConfig as base + +options = base.parser.parse_args() + +# Use the input file to set the output file name +inFilename = options.inFilename +outFilename = options.outFilename + +# The seed is the file number. +smearingSeed = options.seed + +print('Input file: %s' % inFilename) +print('Output file: %s' % outFilename) + +p = HpstrConf.Process() + +p.run_mode = 1 +p.skip_events = options.skip_events +p.max_events = options.nevents + + +p.add_library("libprocessors") + +#anaTrks = HpstrConf.Processor('anaTrks', 'TrackHitAnaProcessor') +anaTrks = HpstrConf.Processor('anaTrks', 'TrackingAnaProcessor') +anaTrks.parameters["debug"] = 0 +anaTrks.parameters["seed"] = smearingSeed +anaTrks.parameters["trkCollName"] = 'KalmanFullTracks' +anaTrks.parameters["histCfg"] = os.environ['HPSTR_BASE']+'/analysis/plotconfigs/fee_smearing/feeSmearing.json' +anaTrks.parameters["selectionjson"] = os.environ['HPSTR_BASE']+'/analysis/selections/trackHit/trackHitAna.json' +anaTrks.parameters["isData"] = options.isData + +#SmearingClosureTest +anaTrks.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+"/utils/data/smearingFile_2016_all_12112023.root" + +RegionPath = os.environ['HPSTR_BASE']+"/analysis/selections/feeSmearing/" +anaTrks.parameters["regionDefinitions"] = [] + +p.sequence = [anaTrks] + +p.input_files = inFilename +p.output_files = [outFilename] + +p.printProcess() diff --git a/processors/config/anaSimps_2016_cfg.py b/processors/config/anaSimps_2016_cfg.py index 898ebb05d..911b66881 100644 --- a/processors/config/anaSimps_2016_cfg.py +++ b/processors/config/anaSimps_2016_cfg.py @@ -51,6 +51,7 @@ vtxana.parameters["isRadPDG"] = options.isRadPDG vtxana.parameters["makeFlatTuple"] = options.makeFlatTuple vtxana.parameters["beamPosCfg"] = "" +vtxana.parameters["pSmearingFile"] = os.environ['HPSTR_BASE']+"/utils/data/smearingFile_2016_all_12112023.root" if options.isData: vtxana.parameters["v0ProjectionFitsCfg"] = os.environ['HPSTR_BASE']+'/analysis/data/v0_projection_2016_config.json' else: diff --git a/processors/config/baseConfig.py b/processors/config/baseConfig.py index 23a344a58..98a5f3258 100644 --- a/processors/config/baseConfig.py +++ b/processors/config/baseConfig.py @@ -22,6 +22,7 @@ help="Which analysis is being run ", metavar="analysis", default="vertex") parser.add_argument('--infile', '-i', type=str, dest="inFilename", metavar='infiles', nargs="+", help="Input files, specify on or more.") +parser.add_argument('--seed', '-s', type=int, dest='seed',metavar='seed',default=42) #options = parser.parse_args() diff --git a/processors/include/NewVertexAnaProcessor.h b/processors/include/NewVertexAnaProcessor.h index 111a0d1ba..5cb785025 100644 --- a/processors/include/NewVertexAnaProcessor.h +++ b/processors/include/NewVertexAnaProcessor.h @@ -17,6 +17,7 @@ #include "MCAnaHistos.h" #include "utilities.h" +#include "TrackSmearingTool.h" #include "FlatTupleMaker.h" #include "AnaHelpers.h" @@ -114,6 +115,9 @@ class NewVertexAnaProcessor : public Processor { int makeFlatTuple_{0}; //!< make true in config to save flat tuple TTree* tree_{nullptr}; //!< description + std::string pSmearingFile_{""}; + std::shared_ptr smearingTool_; + std::shared_ptr _vtx_histos; //!< description std::shared_ptr _mc_vtx_histos; //!< description diff --git a/processors/include/TrackingAnaProcessor.h b/processors/include/TrackingAnaProcessor.h index ec0745a6c..43d9d5cd8 100644 --- a/processors/include/TrackingAnaProcessor.h +++ b/processors/include/TrackingAnaProcessor.h @@ -11,6 +11,8 @@ // ROOT // //----------// #include "TClonesArray.h" +#include "TH1D.h" +#include "TH2D.h" //-----------// // hpstr // @@ -22,6 +24,7 @@ #include "CalCluster.h" #include "EventHeader.h" #include "TrackHistos.h" +#include "TrackSmearingTool.h" // Forward declarations class TTree; @@ -111,8 +114,22 @@ class TrackingAnaProcessor : public Processor { bool doTruth_{false}; //!< description int isData_{1}; //! is data int debug_{0}; //!< debug level + int seed_{0}; //!< seed float time_offset_{0}; //! time offset - + + //Momentum smearing closure test + std::shared_ptr smearingTool_; + std::shared_ptr smearingToolRel_; + std::string pSmearingFile_{""}; + TH1D* psmear_h_; + TH2D* psmear_vs_nHits_hh_; + TH2D* psmear_vs_nHits_top_hh_; + TH2D* psmear_vs_nHits_bot_hh_; + + TH1D* psmear_rel_h_; + TH2D* psmear_vs_nHits_rel_hh_; + TH2D* psmear_vs_nHits_top_rel_hh_; + TH2D* psmear_vs_nHits_bot_rel_hh_; }; // TrackingAnaProcessor diff --git a/processors/src/NewVertexAnaProcessor.cxx b/processors/src/NewVertexAnaProcessor.cxx index 5f10ddb9d..969095607 100644 --- a/processors/src/NewVertexAnaProcessor.cxx +++ b/processors/src/NewVertexAnaProcessor.cxx @@ -38,6 +38,8 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) { isData_ = parameters.getInteger("isData",isData_); analysis_ = parameters.getString("analysis"); + pSmearingFile_ = parameters.getString("pSmearingFile",pSmearingFile_); + //region definitions regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_); @@ -244,6 +246,11 @@ void NewVertexAnaProcessor::initialize(TTree* tree) { //tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); if (brMap_.find(hitColl_.c_str()) != brMap_.end()) tree_->SetBranchAddress(hitColl_.c_str(), &hits_ , &bhits_); if(!isData_ && !mcColl_.empty()) tree_->SetBranchAddress(mcColl_.c_str() , &mcParts_, &bmcParts_); + + if (not pSmearingFile_.empty()) { + // just using the same seed=42 for now + smearingTool_ = std::make_shared(pSmearingFile_,true); + } } bool NewVertexAnaProcessor::process(IEvent* ievent) { @@ -350,6 +357,15 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { ele_trk.applyCorrection("track_time",eleTrackTimeBias_); pos_trk.applyCorrection("track_time", posTrackTimeBias_); + double invm_smear = 1.; + if (smearingTool_) { + double unsmeared_prod = ele_trk.getP()*pos_trk.getP(); + smearingTool_->updateWithSmearP(ele_trk); + smearingTool_->updateWithSmearP(pos_trk); + 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]); //pos_trk.setMomentum(pos->getMomentum()[0],pos->getMomentum()[1],pos->getMomentum()[2]); @@ -530,6 +546,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { _vtx_histos->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight); _vtx_histos->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight); _vtx_histos->Fill1DHisto("vtx_Esum_h", ele_E + pos_E, weight); + _vtx_histos->Fill1DHisto("vtx_smear_InvM_h", invm_smear*(vtx->getInvMass()), weight); _vtx_histos->Fill1DHisto("ele_pos_clusTimeDiff_h", (corr_eleClusterTime - corr_posClusterTime), weight); _vtx_histos->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk.getIsolation(0), ele_trk.getIsolation(1)), vtx->getZ(), weight); _vtx_histos->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk.getIsolation(0), pos_trk.getIsolation(1)), vtx->getZ(), weight); @@ -640,6 +657,15 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //Track Time Corrections ele_trk.applyCorrection("track_time",eleTrackTimeBias_); pos_trk.applyCorrection("track_time", posTrackTimeBias_); + + double invm_smear = 1.; + if (smearingTool_) { + double unsmeared_prod = ele_trk.getP()*pos_trk.getP(); + smearingTool_->updateWithSmearP(ele_trk); + smearingTool_->updateWithSmearP(pos_trk); + double smeared_prod = ele_trk.getP()*pos_trk.getP(); + invm_smear = sqrt(smeared_prod/unsmeared_prod); + } //Add the momenta to the tracks //ele_trk.setMomentum(ele->getMomentum()[0],ele->getMomentum()[1],ele->getMomentum()[2]); @@ -1026,6 +1052,14 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { //Track Time Corrections ele_trk.applyCorrection("track_time",eleTrackTimeBias_); pos_trk.applyCorrection("track_time", posTrackTimeBias_); + double invm_smear = 1.; + if (smearingTool_) { + double unsmeared_prod = ele_trk.getP()*pos_trk.getP(); + smearingTool_->updateWithSmearP(ele_trk); + smearingTool_->updateWithSmearP(pos_trk); + double smeared_prod = ele_trk.getP()*pos_trk.getP(); + invm_smear = sqrt(smeared_prod/unsmeared_prod); + } //Get the layers hit on each track std::vector ele_hit_layers = ele_trk.getHitLayers(); @@ -1147,6 +1181,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) { _reg_vtx_histos[region]->Fill1DHisto("pos_track_n2dhits_h", pos2dHits, weight); _reg_vtx_histos[region]->Fill1DHisto("vtx_Psum_h", p_ele.P()+p_pos.P(), weight); _reg_vtx_histos[region]->Fill1DHisto("vtx_Esum_h", eleClus.getEnergy()+posClus.getEnergy(), weight); + _reg_vtx_histos[region]->Fill1DHisto("vtx_smear_InvM_h", invm_smear*(vtx->getInvMass()), weight); _reg_vtx_histos[region]->Fill2DHisto("ele_vtxZ_iso_hh", TMath::Min(ele_trk.getIsolation(0), ele_trk.getIsolation(1)), vtx->getZ(), weight); _reg_vtx_histos[region]->Fill2DHisto("pos_vtxZ_iso_hh", TMath::Min(pos_trk.getIsolation(0), pos_trk.getIsolation(1)), vtx->getZ(), weight); _reg_vtx_histos[region]->Fill2DTrack(&ele_trk,weight,"ele_"); diff --git a/processors/src/TrackingAnaProcessor.cxx b/processors/src/TrackingAnaProcessor.cxx index d641a706b..e8a651084 100644 --- a/processors/src/TrackingAnaProcessor.cxx +++ b/processors/src/TrackingAnaProcessor.cxx @@ -16,6 +16,7 @@ void TrackingAnaProcessor::configure(const ParameterSet& parameters) { try { debug_ = parameters.getInteger("debug",debug_); + seed_ = parameters.getInteger("seed",seed_); trkCollName_ = parameters.getString("trkCollName",trkCollName_); histCfgFilename_ = parameters.getString("histCfg",histCfgFilename_); doTruth_ = (bool) parameters.getInteger("doTruth",doTruth_); @@ -24,6 +25,10 @@ void TrackingAnaProcessor::configure(const ParameterSet& parameters) { isData_ = parameters.getInteger("isData",isData_); ecalCollName_ = parameters.getString("ecalCollName",ecalCollName_); regionSelections_ = parameters.getVString("regionDefinitions",regionSelections_); + + //Momentum smearing closure test + pSmearingFile_ = parameters.getString("pSmearingFile",pSmearingFile_); + } catch (std::runtime_error& error) { @@ -91,6 +96,50 @@ void TrackingAnaProcessor::initialize(TTree* tree) { if (!ecalCollName_.empty()) tree->SetBranchAddress(ecalCollName_.c_str(),&ecal_, &becal_); + + //Momentum smearing closure test + if (!pSmearingFile_.empty()) { + + std::cout<<"Smearing Tool Seed "<(pSmearingFile_,false,seed_); + smearingToolRel_ = std::make_shared(pSmearingFile_,true, seed_); + + psmear_h_ = new TH1D("psmear_h", + "psmear_h",200,0,4); + + psmear_vs_nHits_hh_ = new TH2D("psmear_vs_nHits_hh", + "psmear_vs_nHits_hh", + 5,8,13, + 200,0,4); + psmear_vs_nHits_top_hh_ = new TH2D("psmear_vs_nHits_top_hh", + "psmear_vs_nHits_top_hh", + 5,8,13, + 200,0,4); + psmear_vs_nHits_bot_hh_ = new TH2D("psmear_vs_nHits_bot_hh", + "psmear_vs_nHits_bot_hh", + 5,8,13, + 200,0,4); + + + psmear_rel_h_ = new TH1D("psmear_rel_h", + "psmear_rel_h",200,0,4); + + psmear_vs_nHits_rel_hh_ = new TH2D("psmear_vs_nHits_rel_hh", + "psmear_vs_nHits_rel_hh", + 5,8,13, + 200,0,4); + psmear_vs_nHits_top_rel_hh_ = new TH2D("psmear_vs_nHits_top_rel_hh", + "psmear_vs_nHits_top_rel_hh", + 5,8,13, + 200,0,4); + psmear_vs_nHits_bot_rel_hh_ = new TH2D("psmear_vs_nHits_bot_rel_hh", + "psmear_vs_nHits_bot_rel_hh", + 5,8,13, + 200,0,4); + + } + + } bool TrackingAnaProcessor::process(IEvent* ievent) { @@ -116,8 +165,8 @@ bool TrackingAnaProcessor::process(IEvent* ievent) { maxTime = 50; } - if (ecal_->size() <= 2) - return true; + //if (ecal_->size() <= 2) + // return true; bool foundFeeCluster = false; @@ -223,15 +272,62 @@ bool TrackingAnaProcessor::process(IEvent* ievent) { trkHistos_->Fill2DHisto("p_vs_nHits_hh", track->getTrackerHitCount(), track->getP()); - + + if (track->getTanLambda() > 0 ) + trkHistos_->Fill2DHisto("p_vs_nHits_top_hh", + track->getTrackerHitCount(), + track->getP()); + else + trkHistos_->Fill2DHisto("p_vs_nHits_bot_hh", + track->getTrackerHitCount(), + track->getP()); + trkHistos_->Fill3DHisto("p_vs_TanLambda_nHits_hhh", track->getTanLambda(), track->getTrackerHitCount(), track->getP()); - }//Loop on tracks + //pSmearing closure Test + if (!isData_ && !pSmearingFile_.empty()) { + + + //Check that I get a gaussian as expected + double p_base = 1.; + + + double psmear = smearingTool_->smearTrackP(*track); + double psmear_rel = smearingToolRel_->smearTrackP(*track); + + double nhits = track->getTrackerHitCount(); + double isTop = track->getTanLambda() > 0 ? true : false; + + psmear_h_->Fill(psmear); + psmear_vs_nHits_hh_->Fill(nhits,psmear); + + if (isTop) { + psmear_vs_nHits_top_hh_->Fill(nhits,psmear); + } + else { + psmear_vs_nHits_bot_hh_->Fill(nhits,psmear); + } + + + psmear_rel_h_->Fill(psmear_rel); + psmear_vs_nHits_rel_hh_->Fill(nhits,psmear_rel); + + if (isTop) { + psmear_vs_nHits_top_rel_hh_->Fill(nhits,psmear_rel); + } + else { + psmear_vs_nHits_bot_rel_hh_->Fill(nhits,psmear_rel); + } + + } // closer test + + }//Loop on tracks + trkHistos_->Fill1DHisto("n_tracks_h",n_sel_tracks); return true; @@ -258,6 +354,29 @@ void TrackingAnaProcessor::finalize() { reg_selectors_[it->first]->getCutFlowHisto()->Write(); } + if (!pSmearingFile_.empty()) { + outF_->cd(trkCollName_.c_str()); + psmear_h_->Write(); + + psmear_vs_nHits_hh_->Write(); + psmear_vs_nHits_top_hh_->Write(); + psmear_vs_nHits_bot_hh_->Write(); + delete psmear_h_; + delete psmear_vs_nHits_hh_; + delete psmear_vs_nHits_top_hh_; + delete psmear_vs_nHits_bot_hh_; + + psmear_rel_h_->Write(); + psmear_vs_nHits_rel_hh_->Write(); + psmear_vs_nHits_top_rel_hh_->Write(); + psmear_vs_nHits_bot_rel_hh_->Write(); + delete psmear_rel_h_; + delete psmear_vs_nHits_rel_hh_; + delete psmear_vs_nHits_top_rel_hh_; + delete psmear_vs_nHits_bot_rel_hh_; + + } + //trkHistos_->Clear(); } diff --git a/scripts/run_jobPool.py b/scripts/run_jobPool.py index 9b561c831..506a0e7ce 100644 --- a/scripts/run_jobPool.py +++ b/scripts/run_jobPool.py @@ -68,6 +68,7 @@ def launchTestsArgs(options, infilename, fileN, jobN): "-i", infilename, "-o", outfilename, "-t", str(options.isData), + "-s", str(jobN), #Each job gets a different seed options.extraFlags, ] print(cmd) diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt new file mode 100644 index 000000000..d54069b67 --- /dev/null +++ b/utils/CMakeLists.txt @@ -0,0 +1,7 @@ + +# Declare processors module +module( + NAME utils + DEPENDENCIES event + EXTERNAL_DEPENDENCIES ROOT +) diff --git a/utils/data/README.md b/utils/data/README.md new file mode 100644 index 000000000..afaf4bc62 --- /dev/null +++ b/utils/data/README.md @@ -0,0 +1,24 @@ +# Utils data area + +## TrackSmearingTool + +The smearing terms from the fee analysis are stored in .root files. +The basic smearing terms are the stored in 1D histograms as function of the number of hits on tracks, separated for top and bottom SVT halves. + +### Files naming convention: +`feeSmearing___.root` +- _year_: is for which dataset to be used +- _iov_: the interval of validity of the smearing terms within a certain dataset. "all" for the full dataset. +- _productionDate_: when the file has been done. + +### Histograms naming convention: +The histogram storing the smearing terms as function of the number of hits is named: +``` +_p_vs_nHits_hh_smearing +``` +- _Tracks_: name of the tracks used for evaluating the smearing terms. "KalmanFullTracks" is the default. + +### Generate smeering terms +``` +python3.9 smearingPlots.py -i fee_2pt3_recon.root -m fee_2pt3_recon_mc.root +``` diff --git a/utils/data/smearingFile_2016_all_12112023.root b/utils/data/smearingFile_2016_all_12112023.root new file mode 100644 index 000000000..63031e483 Binary files /dev/null and b/utils/data/smearingFile_2016_all_12112023.root differ diff --git a/utils/include/TrackSmearingTool.h b/utils/include/TrackSmearingTool.h new file mode 100644 index 000000000..819be035e --- /dev/null +++ b/utils/include/TrackSmearingTool.h @@ -0,0 +1,51 @@ +#pragma once + +//------------------// +// C++ // +//------------------// +#include +#include +#include + +//------------------// +// hpstr // +//------------------// + +#include "Track.h" + +class TFile; +class TH1D; + +class TrackSmearingTool { + + public : + + // The seed needs to be set accordingly for each instance / job of the smearing tool + TrackSmearingTool(const std::string& smearingfile, + const bool relSmearing = true, + const int seed = 42, + const std::string& tracks = "KalmanFullTracks"); + + double smearTrackP(const Track& trk); + void updateWithSmearP(Track& trk); + + private: + + //Random engine + std::shared_ptr generator_; + + // General Normal distributions + + std::shared_ptr> normal_; + + std::shared_ptr smearingfile_; + + //Smearing terms + TH1D* smearing_histo_top_; + TH1D* smearing_histo_bot_; + + // debug + bool debug_{false}; + bool relSmearing_{false}; + +}; diff --git a/utils/src/TrackSmearingTool.cxx b/utils/src/TrackSmearingTool.cxx new file mode 100644 index 000000000..769e1dedb --- /dev/null +++ b/utils/src/TrackSmearingTool.cxx @@ -0,0 +1,89 @@ +#include "TrackSmearingTool.h" +#include "TFile.h" +#include "TH1D.h" + +#include + +TrackSmearingTool::TrackSmearingTool(const std::string& smearingfile, + const bool relSmearing, + const int seed, + const std::string& tracks){ + + + relSmearing_ = relSmearing; + std::string hsuffix = relSmearing_ ? "_rel" : ""; + smearingfile_ = std::make_shared(smearingfile.c_str()); + + if (!smearingfile_) + throw std::invalid_argument( "Provided input smearing file does not exists"); + + //cache the smearing histograms + smearing_histo_top_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix).c_str()); + smearing_histo_bot_ = (TH1D*) smearingfile_->Get((tracks+"_p_vs_nHits_hh_smearing"+hsuffix).c_str()); + + if (!smearing_histo_top_ || !smearing_histo_bot_) + throw std::invalid_argument("Top and Bottom smearing histograms not found in smearing file"); + + //setup random engine + if (debug_) + std::cout<<"Setting up random engine with seed "<(seed); + + normal_ = std::make_shared>(0.,1.); + +} + +double TrackSmearingTool::smearTrackP(const Track& track) { + + double p = track.getP(); + double nhits = track.getTrackerHitCount(); + bool isTop = track.getTanLambda() > 0. ? true : false; + int binN = smearing_histo_top_->FindBin(nhits); + + if (debug_) + std::cout<<"Track nhits="< smearing_histo_top_->GetXaxis()->GetNbins()) { + throw std::invalid_argument("Bin not found in smearing histogram"); + } + + double rel_smear = (*normal_)(*generator_); + double sp = 0.; + + if (isTop) + sp = rel_smear * smearing_histo_top_->GetBinContent(binN); + else + sp = rel_smear * smearing_histo_bot_->GetBinContent(binN); + + double psmear = 0.; + + if (relSmearing_) + psmear = p + sp*p; + else + psmear = p + sp; + + + if (debug_) { + std::cout<<"Track isTop: "< unsmeared_sampler, + int trials = 100, + int min_nhits = 8, + int max_nhits = 12 +) { + TrackSmearingTool tst{smearing_file}; + TrackerHit empty_hit; + Track track; + + track.setTrackParameters( + 0. /* d0 */, + 0. /* phi0 */, + 1. /* omega */, + 1. /* tan_lambda - choosing positive (top) */, + 0. /* z0 */ + ); + for (int i{0}; i < min_nhits-1; i++) track.addHit(&empty_hit); + + bool top = true; + int nhits = min_nhits; + double unsmeared_momentum = 1.; + double smeared_momentum = 1.; + + output_d->cd(); + TTree data{name.c_str(), name.c_str()}; + data.Branch("top", &top); + data.Branch("nhits", &nhits); + data.Branch("unsmeared_momentum", &unsmeared_momentum); + data.Branch("smeared_momentum", &smeared_momentum); + + for (nhits = min_nhits; nhits < max_nhits+1; nhits++) { + track.addHit(&empty_hit); + for (int _i{0}; _i < trials; _i++) { + unsmeared_momentum = unsmeared_sampler(); + track.setMomentum(0. /*px*/, 0. /*py*/, unsmeared_momentum /*pz*/); + tst.updateWithSmearP(track); + smeared_momentum = track.getP(); + data.Fill(); + } + } + data.Write(); +} + +double fixed_momentum_sampler() { + return 1.0; +} + +double normal_momentum() { + static std::normal_distribution dist(1.0, 0.05); + static std::mt19937 gen; + return dist(gen); +} + +int main(int nargs, char* argv[]) try { + /* + std::cout << nargs << " : "; + for (int i{0}; i < nargs; i++) std::cout << argv[i] << " "; + std::cout << std::endl; + */ + if (nargs < 2) { + usage(); + return 1; + } else if (nargs < 3) { + usage(); + std::cout << "\n ERROR: Need two files for testing." << std::endl; + } + std::string smearing_file{argv[1]}, output_file{argv[2]}; + TFile output{argv[2], "RECREATE"}; + + smearing_test( + &output, + smearing_file, + "fixed", + fixed_momentum_sampler + ); + + smearing_test( + &output, + smearing_file, + "gaussian", + normal_momentum, + 10000 /* num trials */ + ); + + output.Write(); + output.Close(); +} catch (const std::runtime_error& rte) { + std::cerr << rte.what() << std::endl; + return 127; +}