TrackFindingTMTT/src/Histos.cc

0001 #include "L1Trigger/TrackFindingTMTT/interface/Histos.h"
0002 #include "L1Trigger/TrackFindingTMTT/interface/InputData.h"
0003 #include "L1Trigger/TrackFindingTMTT/interface/Sector.h"
0004 #include "L1Trigger/TrackFindingTMTT/interface/HTrphi.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/Make3Dtracks.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/TrkRZfilter.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
0008 #include "L1Trigger/TrackFindingTMTT/interface/Utility.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0010
0011 #include "DataFormats/Math/interface/deltaPhi.h"
0012 #include "DataFormats/Math/interface/deltaR.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014
0015 #include <TH1F.h>
0016 #include <TH2F.h>
0017 #include <TH2Poly.h>
0018 #include <TF1.h>
0019 #include <TPad.h>
0020 #include <TProfile.h>
0021 #include <TGraphAsymmErrors.h>
0022 #include <TGraph.h>
0023 #include <TEfficiency.h>
0024
0025 #include <algorithm>
0026 #include <array>
0027 #include <unordered_set>
0028 #include <list>
0029 #include <sstream>
0030 #include <memory>
0031 #include <mutex>
0032
0033 using namespace std;
0034
0035 namespace tmtt {
0036
0037   //=== Store cfg parameters.
0038
0039   Histos::Histos(const Settings* settings) : settings_(settings), oldSumW2opt_(false), bApproxMistake_(false) {
0040     genMinStubLayers_ = settings->genMinStubLayers();
0041     numPhiSectors_ = settings->numPhiSectors();
0042     numEtaRegions_ = settings->numEtaRegions();
0043     houghMinPt_ = settings->houghMinPt();
0044     houghNbinsPt_ = settings->houghNbinsPt();
0045     houghNbinsPhi_ = settings->houghNbinsPhi();
0046     chosenRofZ_ = settings->chosenRofZ();
0047     trackFitters_ = settings->trackFitters();
0048     useRZfilter_ = settings->useRZfilter();
0049     ranRZfilter_ = (not useRZfilter_.empty());  // Was any r-z track filter run?
0050     resPlotOpt_ = settings->resPlotOpt();       // Only use signal events for helix resolution plots?
0051   }
0052
0053   //=== Book all histograms
0054
0055   void Histos::book() {
0056     // Don't bother booking histograms if user didn't request them via TFileService in their cfg.
0057     if (not this->enabled())
0058       return;
0059
0060     oldSumW2opt_ = TH1::GetDefaultSumw2();
0061     TH1::SetDefaultSumw2(true);
0062
0063     // Book histograms about input data.
0064     this->bookInputData();
0065     // Book histograms checking if (eta,phi) sector definition choices are good.
0066     this->bookEtaPhiSectors();
0067     // Book histograms checking filling of r-phi HT array.
0068     this->bookRphiHT();
0069     // Book histograms about r-z track filters.
0070     if (ranRZfilter_)
0071       this->bookRZfilters();
0072     // Book histograms studying 3D track candidates found after HT.
0073     this->bookTrackCands("HT");
0074     // Book histograms studying 3D track candidates found after r-z track filter.
0075     if (ranRZfilter_)
0076       this->bookTrackCands("RZ");
0077     // Book histograms studying track fitting performance
0078     this->bookTrackFitting();
0079   }
0080
0081   //=== Fill all histograms
0082
0083   void Histos::fill(const InputData& inputData,
0084                     const Array2D<unique_ptr<Sector>>& mSectors,
0085                     const Array2D<unique_ptr<HTrphi>>& mHtRphis,
0086                     const Array2D<unique_ptr<Make3Dtracks>>& mMake3Dtrks,
0087                     const std::map<std::string, std::list<const L1fittedTrack*>>& mapFinalTracks) {
0088     // Each function here protected by a mytex lock, so only one thread can run it at a time.
0089
0090     // Don't bother filling histograms if user didn't request them via TFileService in their cfg.
0091     if (not this->enabled())
0092       return;
0093
0094     // Fill histograms about input data.
0095     this->fillInputData(inputData);
0096
0097     // Fill histograms checking if (eta,phi) sector definition choices are good.
0098     this->fillEtaPhiSectors(inputData, mSectors);
0099
0100     // Fill histograms checking filling of r-phi HT array.
0101     this->fillRphiHT(mHtRphis);
0102
0103     // Fill histograms about r-z track filters.
0104     if (ranRZfilter_)
0105       this->fillRZfilters(mMake3Dtrks);
0106
0107     // Fill histograms studying 3D track candidates found after HT.
0108     this->fillTrackCands(inputData, mMake3Dtrks, "HT");
0109
0110     // Fill histograms studying 3D track candidates found after r-z track filter.
0111     if (ranRZfilter_) {
0112       this->fillTrackCands(inputData, mMake3Dtrks, "RZ");
0113     }
0114
0115     // Fill histograms studying track fitting performance
0116     this->fillTrackFitting(inputData, mapFinalTracks);
0117   }
0118
0119   //=== Book histograms using input stubs and tracking particles.
0120
0121   TFileDirectory Histos::bookInputData() {
0122     TFileDirectory inputDir = fs_->mkdir("InputData");
0123
0124     hisStubsVsEta_ = inputDir.make<TH1F>("StubsVsEta", "; #eta; No. stubs in tracker", 30, -3.0, 3.0);
0125     hisStubsVsR_ = inputDir.make<TH1F>("StubsVsR", "; radius (cm); No. stubs in tracker", 1200, 0., 120.);
0126
0127     hisNumLayersPerTP_ =
0128         inputDir.make<TH1F>("NumLayersPerTP", "; Number of layers per TP for alg. eff.", 20, -0.5, 19.5);
0129     hisNumPSLayersPerTP_ =
0130         inputDir.make<TH1F>("NumPSLayersPerTP", "; Number of PS layers per TP for alg. eff.", 20, -0.5, 19.5);
0131
0132     // Study efficiency of tightened front end-electronics cuts.
0133
0134     hisStubKillFE_ = inputDir.make<TProfile>(
0135         "StubKillFE", "; barrelLayer or 10+endcapRing; Stub fraction rejected by FE chip", 30, -0.5, 29.5);
0136     hisStubIneffiVsInvPt_ =
0137         inputDir.make<TProfile>("StubIneffiVsPt", "; 1/Pt; Inefficiency of FE chip for good stubs", 25, 0.0, 0.5);
0138     hisStubIneffiVsEta_ =
0139         inputDir.make<TProfile>("StubIneffiVsEta", "; |#eta|; Inefficiency of FE chip for good stubs", 15, 0.0, 3.0);
0140
0141     // Study stub resolution.
0142
0143     hisBendStub_ = inputDir.make<TH1F>("BendStub", "; Stub bend in units of strips", 59, -7.375, 7.375);
0144     hisBendResStub_ = inputDir.make<TH1F>("BendResStub", "; Stub bend minus TP bend in units of strips", 100, -5., 5.);
0145
0146     // Histos for denominator of tracking efficiency
0147     hisTPinvptForEff_ = inputDir.make<TH1F>("TPinvptForEff", "; TP 1/Pt (for effi.);", 50, 0., 0.5);
0148     hisTPetaForEff_ = inputDir.make<TH1F>("TPetaForEff", "; TP #eta (for effi.);", 20, -3., 3.);
0149     hisTPphiForEff_ = inputDir.make<TH1F>("TPphiForEff", "; TP #phi (for effi.);", 20, -M_PI, M_PI);
0150     hisTPd0ForEff_ = inputDir.make<TH1F>("TPd0ForEff", "; TP d0 (for effi.);", 40, 0., 4.);
0151     hisTPz0ForEff_ = inputDir.make<TH1F>("TPz0ForEff", "; TP z0 (for effi.);", 50, 0., 25.);
0152     //
0153     hisTPinvptForAlgEff_ = inputDir.make<TH1F>("TPinvptForAlgEff", "; TP 1/Pt (for alg. effi.);", 50, 0., 0.5);
0154     hisTPetaForAlgEff_ = inputDir.make<TH1F>("TPetaForAlgEff", "; TP #eta (for alg. effi.);", 20, -3., 3.);
0155     hisTPphiForAlgEff_ = inputDir.make<TH1F>("TPphiForAlgEff", "; TP #phi (for alg. effi.);", 20, -M_PI, M_PI);
0156     hisTPd0ForAlgEff_ = inputDir.make<TH1F>("TPd0ForAlgEff", "; TP d0 (for alg. effi.);", 40, 0., 4.);
0157     hisTPz0ForAlgEff_ = inputDir.make<TH1F>("TPz0ForAlgEff", "; TP z0 (for alg. effi.);", 50, 0., 25.);
0158
0159     return inputDir;
0160   }
0161
0162   //=== Fill histograms using input stubs and tracking particles.
0163
0164   void Histos::fillInputData(const InputData& inputData) {
0165     // Allow only one thread to run this function at a time
0166     static std::mutex myMutex;
0167     std::lock_guard<std::mutex> myGuard(myMutex);
0168
0169     const list<const Stub*>& vStubs = inputData.stubsConst();
0170     const list<TP>& vTPs = inputData.getTPs();
0171
0172     for (const Stub* stub : vStubs) {
0173       hisStubsVsEta_->Fill(stub->eta());
0174       hisStubsVsR_->Fill(stub->r());
0175     }
0176
0177     // Study efficiency of stubs to pass front-end electronics cuts.
0178
0179     const list<Stub>& vAllStubs = inputData.allStubs();  // Get all stubs prior to FE cuts to do this.
0180     for (const Stub& s : vAllStubs) {
0181       unsigned int layerOrTenPlusRing = s.barrel() ? s.layerId() : 10 + s.trackerModule()->endcapRing();
0182       // Fraction of all stubs (good and bad) failing tightened front-end electronics cuts.
0183       hisStubKillFE_->Fill(layerOrTenPlusRing, (!s.frontendPass()));
0184     }
0185
0186     // Study efficiency for good stubs of tightened front end-electronics cuts.
0187     for (const TP& tp : vTPs) {
0188       if (tp.useForAlgEff()) {  // Only bother for stubs that are on TP that we have a chance of reconstructing.
0189         const vector<const Stub*>& stubs = tp.assocStubs();
0190         for (const Stub* s : stubs) {
0191           hisStubIneffiVsInvPt_->Fill(1. / tp.pt(), (!s->frontendPass()));
0192           hisStubIneffiVsEta_->Fill(std::abs(tp.eta()), (!s->frontendPass()));
0193         }
0194       }
0195     }
0196
0197     // Plot stub bend-derived information.
0198     for (const Stub* stub : vStubs) {
0199       hisBendStub_->Fill(stub->bend());
0200     }
0201
0202     // Look at stub resolution.
0203     for (const TP& tp : vTPs) {
0204       if (tp.useForAlgEff()) {
0205         const vector<const Stub*>& assStubs = tp.assocStubs();
0206
0207         for (const Stub* stub : assStubs) {
0208           hisBendResStub_->Fill(stub->bend() - tp.dphi(stub->r()) / stub->dphiOverBend());
0209         }
0210
0211         if (std::abs(tp.eta()) < 0.5) {
0212           double nLayersOnTP = Utility::countLayers(settings_, assStubs, true, false);
0213           double nPSLayersOnTP = Utility::countLayers(settings_, assStubs, true, true);
0214           hisNumLayersPerTP_->Fill(nLayersOnTP);
0215           hisNumPSLayersPerTP_->Fill(nPSLayersOnTP);
0216         }
0217       }
0218     }
0219
0220     // Determine r (z) range of each barrel layer (endcap wheel).
0221
0222     for (const Stub* stub : vStubs) {
0223       unsigned int layer = stub->layerId();
0224       if (stub->barrel()) {
0225         // Get range in r of each barrel layer.
0226         float r = stub->r();
0227         if (mapBarrelLayerMinR_.find(layer) == mapBarrelLayerMinR_.end()) {
0228           mapBarrelLayerMinR_[layer] = r;
0229           mapBarrelLayerMaxR_[layer] = r;
0230         } else {
0231           if (mapBarrelLayerMinR_[layer] > r)
0232             mapBarrelLayerMinR_[layer] = r;
0233           if (mapBarrelLayerMaxR_[layer] < r)
0234             mapBarrelLayerMaxR_[layer] = r;
0235         }
0236       } else {
0237         layer = layer % 10;
0238         // Range in |z| of each endcap wheel.
0239         float z = std::abs(stub->z());
0240         if (mapEndcapWheelMinZ_.find(layer) == mapEndcapWheelMinZ_.end()) {
0241           mapEndcapWheelMinZ_[layer] = z;
0242           mapEndcapWheelMaxZ_[layer] = z;
0243         } else {
0244           if (mapEndcapWheelMinZ_[layer] > z)
0245             mapEndcapWheelMinZ_[layer] = z;
0246           if (mapEndcapWheelMaxZ_[layer] < z)
0247             mapEndcapWheelMaxZ_[layer] = z;
0248         }
0249       }
0250     }
0251
0252     // Determine Range in (r,|z|) of each module type.
0253
0254     for (const Stub* stub : vStubs) {
0255       float r = stub->r();
0256       float z = std::abs(stub->z());
0257       unsigned int modType = stub->trackerModule()->moduleTypeID();
0258       // Do something ugly, as modules in 1-2nd & 3-4th endcap wheels are different to those in wheel 5 ...
0259       // And boundary between flat & tilted modules in barrel layers 1-3 varies in z.
0260       if (stub->barrel() && stub->layerId() == 1) {  // barrel layer 1
0261         if (mapExtraAModuleTypeMinR_.find(modType) == mapExtraAModuleTypeMinR_.end()) {
0262           mapExtraAModuleTypeMinR_[modType] = r;
0263           mapExtraAModuleTypeMaxR_[modType] = r;
0264           mapExtraAModuleTypeMinZ_[modType] = z;
0265           mapExtraAModuleTypeMaxZ_[modType] = z;
0266         } else {
0267           if (mapExtraAModuleTypeMinR_[modType] > r)
0268             mapExtraAModuleTypeMinR_[modType] = r;
0269           if (mapExtraAModuleTypeMaxR_[modType] < r)
0270             mapExtraAModuleTypeMaxR_[modType] = r;
0271           if (mapExtraAModuleTypeMinZ_[modType] > z)
0272             mapExtraAModuleTypeMinZ_[modType] = z;
0273           if (mapExtraAModuleTypeMaxZ_[modType] < z)
0274             mapExtraAModuleTypeMaxZ_[modType] = z;
0275         }
0276       } else if (stub->barrel() && stub->layerId() == 2) {  // barrel layer 2
0277         if (mapExtraBModuleTypeMinR_.find(modType) == mapExtraBModuleTypeMinR_.end()) {
0278           mapExtraBModuleTypeMinR_[modType] = r;
0279           mapExtraBModuleTypeMaxR_[modType] = r;
0280           mapExtraBModuleTypeMinZ_[modType] = z;
0281           mapExtraBModuleTypeMaxZ_[modType] = z;
0282         } else {
0283           if (mapExtraBModuleTypeMinR_[modType] > r)
0284             mapExtraBModuleTypeMinR_[modType] = r;
0285           if (mapExtraBModuleTypeMaxR_[modType] < r)
0286             mapExtraBModuleTypeMaxR_[modType] = r;
0287           if (mapExtraBModuleTypeMinZ_[modType] > z)
0288             mapExtraBModuleTypeMinZ_[modType] = z;
0289           if (mapExtraBModuleTypeMaxZ_[modType] < z)
0290             mapExtraBModuleTypeMaxZ_[modType] = z;
0291         }
0292       } else if (!stub->barrel() && (stub->layerId() % 10 == 1 || stub->layerId() % 10 == 2)) {  // endcap wheel 1-2
0293         if (mapExtraCModuleTypeMinR_.find(modType) == mapExtraCModuleTypeMinR_.end()) {
0294           mapExtraCModuleTypeMinR_[modType] = r;
0295           mapExtraCModuleTypeMaxR_[modType] = r;
0296           mapExtraCModuleTypeMinZ_[modType] = z;
0297           mapExtraCModuleTypeMaxZ_[modType] = z;
0298         } else {
0299           if (mapExtraCModuleTypeMinR_[modType] > r)
0300             mapExtraCModuleTypeMinR_[modType] = r;
0301           if (mapExtraCModuleTypeMaxR_[modType] < r)
0302             mapExtraCModuleTypeMaxR_[modType] = r;
0303           if (mapExtraCModuleTypeMinZ_[modType] > z)
0304             mapExtraCModuleTypeMinZ_[modType] = z;
0305           if (mapExtraCModuleTypeMaxZ_[modType] < z)
0306             mapExtraCModuleTypeMaxZ_[modType] = z;
0307         }
0308       } else if (!stub->barrel() && (stub->layerId() % 10 == 3 || stub->layerId() % 10 == 4)) {  // endcap wheel 3-4
0309         if (mapExtraDModuleTypeMinR_.find(modType) == mapExtraDModuleTypeMinR_.end()) {
0310           mapExtraDModuleTypeMinR_[modType] = r;
0311           mapExtraDModuleTypeMaxR_[modType] = r;
0312           mapExtraDModuleTypeMinZ_[modType] = z;
0313           mapExtraDModuleTypeMaxZ_[modType] = z;
0314         } else {
0315           if (mapExtraDModuleTypeMinR_[modType] > r)
0316             mapExtraDModuleTypeMinR_[modType] = r;
0317           if (mapExtraDModuleTypeMaxR_[modType] < r)
0318             mapExtraDModuleTypeMaxR_[modType] = r;
0319           if (mapExtraDModuleTypeMinZ_[modType] > z)
0320             mapExtraDModuleTypeMinZ_[modType] = z;
0321           if (mapExtraDModuleTypeMaxZ_[modType] < z)
0322             mapExtraDModuleTypeMaxZ_[modType] = z;
0323         }
0324       } else {  // barrel layer 3-6 or endcap wheel 5.
0325         if (mapModuleTypeMinR_.find(modType) == mapModuleTypeMinR_.end()) {
0326           mapModuleTypeMinR_[modType] = r;
0327           mapModuleTypeMaxR_[modType] = r;
0328           mapModuleTypeMinZ_[modType] = z;
0329           mapModuleTypeMaxZ_[modType] = z;
0330         } else {
0331           if (mapModuleTypeMinR_[modType] > r)
0332             mapModuleTypeMinR_[modType] = r;
0333           if (mapModuleTypeMaxR_[modType] < r)
0334             mapModuleTypeMaxR_[modType] = r;
0335           if (mapModuleTypeMinZ_[modType] > z)
0336             mapModuleTypeMinZ_[modType] = z;
0337           if (mapModuleTypeMaxZ_[modType] < z)
0338             mapModuleTypeMaxZ_[modType] = z;
0339         }
0340       }
0341     }
0342
0343     //=== Make denominator of tracking efficiency plots
0344
0345     for (const TP& tp : vTPs) {
0346       if (tp.useForEff()) {  // Check TP is good for efficiency measurement.
0347         // Plot kinematics of all good TP.
0348         hisTPinvptForEff_->Fill(1. / tp.pt());
0349         hisTPetaForEff_->Fill(tp.eta());
0350         hisTPphiForEff_->Fill(tp.phi0());
0351         // Plot also production point of all good TP.
0352         hisTPd0ForEff_->Fill(std::abs(tp.d0()));
0353         hisTPz0ForEff_->Fill(std::abs(tp.z0()));
0354
0355         if (tp.useForAlgEff()) {  // Check TP is good for algorithmic efficiency measurement.
0356           hisTPinvptForAlgEff_->Fill(1. / tp.pt());
0357           hisTPetaForAlgEff_->Fill(tp.eta());
0358           hisTPphiForAlgEff_->Fill(tp.phi0());
0359           // Plot also production point of all good TP.
0360           hisTPd0ForAlgEff_->Fill(std::abs(tp.d0()));
0361           hisTPz0ForAlgEff_->Fill(std::abs(tp.z0()));
0362         }
0363       }
0364     }
0365   }
0366
0367   //=== Book histograms checking if (eta,phi) sector defis(nition choices are good.
0368
0369   TFileDirectory Histos::bookEtaPhiSectors() {
0370     TFileDirectory inputDir = fs_->mkdir("CheckSectors");
0371
0372     // Check if stubs excessively duplicated between overlapping sectors.
0373     hisNumEtaSecsPerStub_ =
0374         inputDir.make<TH1F>("NumEtaSecPerStub", "; No. of #eta sectors each stub in", 20, -0.5, 19.5);
0375     hisNumPhiSecsPerStub_ =
0376         inputDir.make<TH1F>("NumPhiSecPerStub", "; No. of #phi sectors each stub in", 20, -0.5, 19.5);
0377
0378     // Count stubs per (eta,phi) sector.
0379     hisNumStubsPerSec_ = inputDir.make<TH1F>("NumStubsPerSec", "; No. of stubs per sector", 150, -0.5, 299.5);
0380
0381     return inputDir;
0382   }
0383
0384   //=== Fill histograms checking if (eta,phi) sector definition choices are good.
0385
0386   void Histos::fillEtaPhiSectors(const InputData& inputData, const Array2D<unique_ptr<Sector>>& mSectors) {
0387     // Allow only one thread to run this function at a time
0388     static std::mutex myMutex;
0389     std::lock_guard<std::mutex> myGuard(myMutex);
0390
0391     const list<const Stub*>& vStubs = inputData.stubsConst();
0392     //const list<TP>& vTPs = inputData.getTPs();
0393
0394     //=== Loop over all stubs, counting how many sectors each one appears in.
0395
0396     for (const Stub* stub : vStubs) {
0397       // Number of (eta,phi), phi & eta sectors containing this stub.
0398       unsigned int nEtaSecs = 0;
0399       unsigned int nPhiSecs = 0;
0400
0401       // Loop over (eta, phi) sectors.
0402       for (unsigned int iPhiSec = 0; iPhiSec < numPhiSectors_; iPhiSec++) {
0403         for (unsigned int iEtaReg = 0; iEtaReg < numEtaRegions_; iEtaReg++) {
0404           const Sector* sector = mSectors(iPhiSec, iEtaReg).get();
0405
0406           // Check if sector contains stub stub, and if so count it.
0407           // Take care to just use one eta (phi) typical region when counting phi (eta) sectors.
0408           if (iPhiSec == 0 && sector->insideEta(stub))
0409             nEtaSecs++;
0410           if (iEtaReg == 0 && sector->insidePhi(stub))
0411             nPhiSecs++;
0412         }
0413       }
0414
0415       // Plot number of sectors each stub appears in.
0416       hisNumEtaSecsPerStub_->Fill(nEtaSecs);
0417       hisNumPhiSecsPerStub_->Fill(nPhiSecs);
0418     }
0419
0420     //=== Loop over all sectors, counting the stubs in each one.
0421     for (unsigned int iEtaReg = 0; iEtaReg < numEtaRegions_; iEtaReg++) {
0422       for (unsigned int iPhiSec = 0; iPhiSec < numPhiSectors_; iPhiSec++) {
0423         const Sector* sector = mSectors(iPhiSec, iEtaReg).get();
0424
0425         unsigned int nStubs = 0;
0426         for (const Stub* stub : vStubs) {
0427           if (sector->inside(stub))
0428             nStubs++;
0429         }
0430         hisNumStubsPerSec_->Fill(nStubs);
0431       }
0432     }
0433   }
0434
0435   //=== Book histograms checking filling of r-phi HT array.
0436
0437   TFileDirectory Histos::bookRphiHT() {
0438     TFileDirectory inputDir = fs_->mkdir("HTrphi");
0439
0440     return inputDir;
0441   }
0442
0443   //=== Fill histograms checking filling of r-phi HT array.
0444
0445   void Histos::fillRphiHT(const Array2D<unique_ptr<HTrphi>>& mHtRphis) {
0446     //--- Loop over (eta,phi) sectors, counting the number of stubs in the HT array of each.
0447
0448     // Allow only one thread to run this function at a time (UNCOMMENT IF YOU ADD HISTOS HERE)
0449     //static std::mutex myMutex;
0450     //std::lock_guard<std::mutex> myGuard(myMutex);
0451   }
0452
0453   //=== Book histograms about r-z track filters (or other filters applied after r-phi HT array).
0454
0455   TFileDirectory Histos::bookRZfilters() {
0456     TFileDirectory inputDir = fs_->mkdir("RZfilters");
0457
0458     return inputDir;
0459   }
0460
0461   //=== Fill histograms about r-z track filters.
0462
0463   void Histos::fillRZfilters(const Array2D<unique_ptr<Make3Dtracks>>& mMake3Dtrks) {
0464     // Allow only one thread to run this function at a time (UNCOMMENT IF YOU ADD HISTOS HERE)
0465     //static std::mutex myMutex;
0466     //std::lock_guard<std::mutex> myGuard(myMutex);
0467   }
0468
0469   //=== Book histograms studying track candidates found by Hough Transform.
0470
0471   TFileDirectory Histos::bookTrackCands(const string& tName) {
0472     // Now book histograms for studying tracking in general.
0473
0474     auto addn = [tName](const string& s) { return TString::Format("%s_%s", s.c_str(), tName.c_str()); };
0475
0476     TFileDirectory inputDir = fs_->mkdir(addn("TrackCands").Data());
0477
0478     bool TMTT = (tName == "HT" || tName == "RZ");
0479
0480     // Count tracks in various ways (including/excluding duplicates, excluding fakes ...)
0481     profNumTrackCands_[tName] =
0482         inputDir.make<TProfile>(addn("NumTrackCands"), "; class; N. of tracks in tracker", 7, 0.5, 7.5);
0483     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(7, "TP for eff recoed");
0484     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(6, "TP recoed");
0485     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(5, "TP recoed x #eta sector dups");
0486     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(4, "TP recoed x sector dups");
0487     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(2, "TP recoed x track dups");
0488     profNumTrackCands_[tName]->GetXaxis()->SetBinLabel(1, "reco tracks including fakes");
0489     profNumTrackCands_[tName]->LabelsOption("d");
0490
0491     hisNumTrksPerNon_[tName] = inputDir.make<TH1F>(addn("NumTrksPerNon"), "; No. tracks per nonant;", 100, -0.5, 399.5);
0492
0493     unsigned int nEta = numEtaRegions_;
0494     hisNumTracksVsQoverPt_[tName] =
0495         inputDir.make<TH1F>(addn("NumTracksVsQoverPt"), "; Q/Pt; No. of tracks in tracker", 100, -0.5, 0.5);
0496     if (TMTT) {
0497       profNumTracksVsEta_[tName] = inputDir.make<TProfile>(
0498           addn("NumTracksVsEta"), "; #eta region; No. of tracks in tracker", nEta, -0.5, nEta - 0.5);
0499     }
0500
0501     // Count stubs per event assigned to tracks (determines HT data output rate)
0502
0503     profStubsOnTracks_[tName] =
0504         inputDir.make<TProfile>(addn("StubsOnTracks"), "; ; No. of stubs on tracks per event", 1, 0.5, 1.5);
0505     hisStubsOnTracksPerNon_[tName] =
0506         inputDir.make<TH1F>(addn("StubsOnTracksPerNon"), "; No. of stubs on tracks per nonant", 100, -0.5, 4999.5);
0507
0508     hisStubsPerTrack_[tName] = inputDir.make<TH1F>(addn("StubsPerTrack"), ";No. of stubs per track;", 50, -0.5, 49.5);
0509     hisLayersPerTrack_[tName] =
0510         inputDir.make<TH1F>(addn("LayersPerTrack"), ";No. of layers with stubs per track;", 20, -0.5, 19.5);
0511
0512     if (TMTT) {
0513       hisNumStubsPerLink_[tName] =
0514           inputDir.make<TH1F>(addn("NumStubsPerLink"), "; Mean #stubs per MHT output opto-link;", 50, -0.5, 249.5);
0515       profMeanStubsPerLink_[tName] =
0516           inputDir.make<TProfile>(addn("MeanStubsPerLink"), "; Mean #stubs per MHT output opto-link;", 36, -0.5, 35.5);
0517     }
0518
0519     hisFracMatchStubsOnTracks_[tName] = inputDir.make<TH1F>(
0520         addn("FracMatchStubsOnTracks"), "; Frac. of stubs per trk matching best TP;", 101, -0.005, 1.005);
0521
0522     if (TMTT) {
0523       // Study duplication of tracks within an individual HT array.
0524       profDupTracksVsEta_[tName] =
0525           inputDir.make<TProfile>(addn("DupTracksVsTPeta"), "; #eta; No. of dup. trks per TP;", 15, 0.0, 3.0);
0526       profDupTracksVsInvPt_[tName] =
0527           inputDir.make<TProfile>(addn("DupTracksVsInvPt"), "; 1/Pt; No. of dup. trks per TP", 25, 0., 0.5);
0528     }
0529
0530     // Histos of track params.
0531     hisQoverPt_[tName] = inputDir.make<TH1F>(addn("QoverPt"), "; track q/Pt", 100, -0.5, 0.5);
0532     hisPhi0_[tName] = inputDir.make<TH1F>(addn("Phi0"), "; track #phi0", 70, -3.5, 3.5);
0533     hisEta_[tName] = inputDir.make<TH1F>(addn("Eta"), "; track #eta", 60, -3.0, 3.0);
0534     hisZ0_[tName] = inputDir.make<TH1F>(addn("Z0"), "; track z0", 100, -25.0, 25.0);
0535
0536     // Histos of track parameter resolution
0537     hisQoverPtRes_[tName] = inputDir.make<TH1F>(addn("QoverPtRes"), "; track resolution in q/Pt", 100, -0.06, 0.06);
0538     hisPhi0Res_[tName] = inputDir.make<TH1F>(addn("Phi0Res"), "; track resolution in #phi0", 100, -0.04, 0.04);
0539     hisEtaRes_[tName] = inputDir.make<TH1F>(addn("EtaRes"), "; track resolution in #eta", 100, -1.0, 1.0);
0540     hisZ0Res_[tName] = inputDir.make<TH1F>(addn("Z0Res"), "; track resolution in z0", 100, -10.0, 10.0);
0541
0542     // Histos for tracking efficiency vs. TP kinematics
0543     hisRecoTPinvptForEff_[tName] =
0544         inputDir.make<TH1F>(addn("RecoTPinvptForEff"), "; TP 1/Pt of recoed tracks (for effi.);", 50, 0., 0.5);
0545     hisRecoTPetaForEff_[tName] =
0546         inputDir.make<TH1F>(addn("RecoTPetaForEff"), "; TP #eta of recoed tracks (for effi.);", 20, -3., 3.);
0547     hisRecoTPphiForEff_[tName] =
0548         inputDir.make<TH1F>(addn("RecoTPphiForEff"), "; TP #phi of recoed tracks (for effi.);", 20, -M_PI, M_PI);
0549
0550     // Histo for efficiency to reconstruct track perfectly (no incorrect hits).
0551     hisPerfRecoTPinvptForEff_[tName] = inputDir.make<TH1F>(
0552         addn("PerfRecoTPinvptForEff"), "; TP 1/Pt of recoed tracks (for perf. effi.);", 50, 0., 0.5);
0553     hisPerfRecoTPetaForEff_[tName] =
0554         inputDir.make<TH1F>(addn("PerfRecoTPetaForEff"), "; TP #eta of recoed tracks (for perf. effi.);", 20, -3., 3.);
0555
0556     // Histos for  tracking efficiency vs. TP production point
0557     hisRecoTPd0ForEff_[tName] =
0558         inputDir.make<TH1F>(addn("RecoTPd0ForEff"), "; TP d0 of recoed tracks (for effi.);", 40, 0., 4.);
0559     hisRecoTPz0ForEff_[tName] =
0560         inputDir.make<TH1F>(addn("RecoTPz0ForEff"), "; TP z0 of recoed tracks (for effi.);", 50, 0., 25.);
0561
0562     // Histos for algorithmic tracking efficiency vs. TP kinematics
0563     hisRecoTPinvptForAlgEff_[tName] =
0564         inputDir.make<TH1F>(addn("RecoTPinvptForAlgEff"), "; TP 1/Pt of recoed tracks (for alg. effi.);", 50, 0., 0.5);
0565     hisRecoTPetaForAlgEff_[tName] =
0566         inputDir.make<TH1F>(addn("RecoTPetaForAlgEff"), "; TP #eta of recoed tracks (for alg. effi.);", 20, -3., 3.);
0567     hisRecoTPphiForAlgEff_[tName] = inputDir.make<TH1F>(
0568         addn("RecoTPphiForAlgEff"), "; TP #phi of recoed tracks (for alg. effi.);", 20, -M_PI, M_PI);
0569
0570     // Histo for efficiency to reconstruct track perfectly (no incorrect hits).
0571     hisPerfRecoTPinvptForAlgEff_[tName] = inputDir.make<TH1F>(
0572         addn("PerfRecoTPinvptForAlgEff"), "; TP 1/Pt of recoed tracks (for perf. alg. effi.);", 50, 0., 0.5);
0573     hisPerfRecoTPetaForAlgEff_[tName] =
0574         inputDir.make<TH1F>(addn("PerfRecoTPetaForAlgEff"), "; TP #eta (for perf. alg. effi.);", 20, -3., 3.);
0575
0576     // Histos for algorithmic tracking efficiency vs. TP production point
0577     hisRecoTPd0ForAlgEff_[tName] =
0578         inputDir.make<TH1F>(addn("RecoTPd0ForAlgEff"), "; TP d0 of recoed tracks (for alg. effi.);", 40, 0., 4.);
0579     hisRecoTPz0ForAlgEff_[tName] =
0580         inputDir.make<TH1F>(addn("RecoTPz0ForAlgEff"), "; TP z0 of recoed tracks (for alg. effi.);", 50, 0., 25.);
0581
0582     return inputDir;
0583   }
0584
0585   //=== Fill histograms studying track candidates found before track fit is run.
0586
0587   void Histos::fillTrackCands(const InputData& inputData,
0588                               const Array2D<std::unique_ptr<Make3Dtracks>>& mMake3Dtrks,
0589                               const string& tName) {
0590     // Allow only one thread to run this function at a time
0591     static std::mutex myMutex;
0592     std::lock_guard<std::mutex> myGuard(myMutex);
0593
0594     vector<L1track3D> tracks;
0595     bool withRZfilter = (tName == "RZ") ? true : false;
0596     for (unsigned int iEtaReg = 0; iEtaReg < numEtaRegions_; iEtaReg++) {
0597       for (unsigned int iPhiSec = 0; iPhiSec < numPhiSectors_; iPhiSec++) {
0598         const Make3Dtracks* make3Dtrk = mMake3Dtrks(iPhiSec, iEtaReg).get();
0599         const std::list<L1track3D>& tracksSec = make3Dtrk->trackCands3D(withRZfilter);
0600         tracks.insert(tracks.end(), tracksSec.begin(), tracksSec.end());
0601       }
0602     }
0603     this->fillTrackCands(inputData, tracks, tName);
0604   }
0605
0606   //=== Fill histograms studying track candidates found before track fit is run.
0607
0608   void Histos::fillTrackCands(const InputData& inputData, const vector<L1track3D>& tracks, const string& tName) {
0609     bool withRZfilter = (tName == "RZ");
0610
0611     bool algoTMTT = (tName == "HT" || tName == "RZ");  // Check if running TMTT or Hybrid L1 tracking.
0612
0613     const list<TP>& vTPs = inputData.getTPs();
0614
0615     //=== Count track candidates found in the tracker.
0616
0617     const unsigned int numPhiNonants = settings_->numPhiNonants();
0618     vector<unsigned int> nTrksPerEtaReg(numEtaRegions_, 0);
0619     vector<unsigned int> nTrksPerNonant(numPhiNonants, 0);
0620     for (const L1track3D& t : tracks) {
0621       unsigned int iNonant = floor((t.iPhiSec()) * numPhiNonants / (numPhiSectors_));  // phi nonant number
0622       nTrksPerEtaReg[t.iEtaReg()]++;
0623       nTrksPerNonant[iNonant]++;
0624     }
0625
0626     profNumTrackCands_[tName]->Fill(1.0, tracks.size());  // Plot mean number of tracks/event.
0627     if (algoTMTT) {
0628       for (unsigned int iEtaReg = 0; iEtaReg < numEtaRegions_; iEtaReg++) {
0629         profNumTracksVsEta_[tName]->Fill(iEtaReg, nTrksPerEtaReg[iEtaReg]);
0630       }
0631     }
0632     for (unsigned int iNonant = 0; iNonant < numPhiNonants; iNonant++) {
0633       hisNumTrksPerNon_[tName]->Fill(nTrksPerNonant[iNonant]);
0634     }
0635
0636     //=== Count stubs per event assigned to track candidates in the Tracker
0637
0638     unsigned int nStubsOnTracks = 0;
0639     vector<unsigned int> nStubsOnTracksInNonant(numPhiNonants, 0);
0640
0641     for (const L1track3D& t : tracks) {
0642       const vector<const Stub*>& stubs = t.stubsConst();
0643       unsigned int nStubs = stubs.size();
0644       unsigned int iNonant = floor((t.iPhiSec()) * numPhiNonants / (numPhiSectors_));  // phi nonant number
0645       // Count stubs on all tracks in this sector & nonant.
0646       nStubsOnTracks += nStubs;
0647       nStubsOnTracksInNonant[iNonant] += nStubs;
0648     }
0649
0650     profStubsOnTracks_[tName]->Fill(1.0, nStubsOnTracks);
0651
0652     for (unsigned int iNonant = 0; iNonant < numPhiNonants; iNonant++) {
0653       hisStubsOnTracksPerNon_[tName]->Fill(nStubsOnTracksInNonant[iNonant]);
0654     }
0655
0656     // Plot number of tracks & number of stubs per output HT opto-link.
0657
0658     if (algoTMTT && not withRZfilter) {
0659       //const unsigned int numPhiSecPerNon = numPhiSectors_ / numPhiNonants;
0660       // Hard-wired bodge
0661       const unsigned int nLinks = houghNbinsPt_ / 2;  // Hard-wired to number of course HT bins. Check.
0662
0663       for (unsigned int iPhiNon = 0; iPhiNon < numPhiNonants; iPhiNon++) {
0664         // Each nonant has a separate set of links.
0665         vector<unsigned int> stubsToLinkCount(nLinks, 0);  // Must use vectors to count links with zero entries.
0666         for (const L1track3D& trk : tracks) {
0667           unsigned int iNonantTrk = floor((trk.iPhiSec()) * numPhiNonants / (numPhiSectors_));  // phi nonant number
0668           if (iPhiNon == iNonantTrk) {
0669             unsigned int link = trk.optoLinkID();
0670             if (link < nLinks) {
0671               stubsToLinkCount[link] += trk.numStubs();
0672             } else {
0673               std::stringstream text;
0674               text << "\n ===== HISTOS MESS UP: Increase size of nLinks ===== " << link << "\n";
0675               static std::once_flag printOnce;
0676               std::call_once(
0677                   printOnce, [](string t) { edm::LogWarning("L1track") << t; }, text.str());
0678             }
0679           }
0680         }
0681
0682         for (unsigned int link = 0; link < nLinks; link++) {
0683           unsigned int nstbs = stubsToLinkCount[link];
0684           hisNumStubsPerLink_[tName]->Fill(nstbs);
0685           profMeanStubsPerLink_[tName]->Fill(link, nstbs);
0686         }
0687       }
0688     }
0689
0690     // Plot q/pt spectrum of track candidates, and number of stubs/tracks
0691     for (const L1track3D& trk : tracks) {
0692       hisNumTracksVsQoverPt_[tName]->Fill(trk.qOverPt());  // Plot reconstructed q/Pt of track cands.
0693       hisStubsPerTrack_[tName]->Fill(trk.numStubs());      // Stubs per track.
0694       hisLayersPerTrack_[tName]->Fill(trk.numLayers());
0695     }
0696
0697     // Count fraction of stubs on each track matched to a TP that are from same TP.
0698
0699     for (const L1track3D& trk : tracks) {
0700       // Only consider tracks that match a tracking particle for the alg. efficiency measurement.
0701       const TP* tp = trk.matchedTP();
0702       if (tp != nullptr) {
0703         if (tp->useForAlgEff()) {
0704           hisFracMatchStubsOnTracks_[tName]->Fill(trk.purity());
0705         }
0706       }
0707     }
0708
0709     // Count total number of tracking particles in the event that were reconstructed,
0710     // counting also how many of them were reconstructed multiple times (duplicate tracks).
0711
0712     unsigned int nRecoedTPsForEff = 0;     // Total no. of TPs for effi measurement recoed as >= 1 track.
0713     unsigned int nRecoedTPs = 0;           // Total no. of TPs recoed as >= 1 one track.
0714     unsigned int nEtaSecsMatchingTPs = 0;  // Total no. of eta sectors that all TPs were reconstructed in
0715     unsigned int nSecsMatchingTPs = 0;     // Total no. of eta x phi sectors that all TPs were reconstructed in
0716     unsigned int nTrksMatchingTPs = 0;     // Total no. of tracks that all TPs were reconstructed as
0717
0718     for (const TP& tp : vTPs) {
0719       vector<const L1track3D*> matchedTrks;
0720       for (const L1track3D& trk : tracks) {
0721         const TP* tpAssoc = trk.matchedTP();
0722         if (tpAssoc != nullptr) {
0723           if (tpAssoc->index() == tp.index())
0724             matchedTrks.push_back(&trk);
0725         }
0726       }
0727       unsigned int nTrk = matchedTrks.size();
0728
0729       bool tpRecoed = false;
0730
0731       if (nTrk > 0) {
0732         tpRecoed = true;           // This TP was reconstructed at least once in tracker.
0733         nTrksMatchingTPs += nTrk;  // Increment sum by no. of tracks this TP was reconstructed as
0734
0735         set<unsigned int> iEtaRegRecoed;
0736         for (const L1track3D* trk : matchedTrks)
0737           iEtaRegRecoed.insert(trk->iEtaReg());
0738         nEtaSecsMatchingTPs = iEtaRegRecoed.size();
0739
0740         set<pair<unsigned int, unsigned int>> iSecRecoed;
0741         for (const L1track3D* trk : matchedTrks)
0742           iSecRecoed.insert({trk->iPhiSec(), trk->iEtaReg()});
0743         nSecsMatchingTPs = iSecRecoed.size();
0744
0745         if (algoTMTT) {
0746           for (const auto& p : iSecRecoed) {
0747             unsigned int nTrkInSec = 0;
0748             for (const L1track3D* trk : matchedTrks) {
0749               if (trk->iPhiSec() == p.first && trk->iEtaReg() == p.second)
0750                 nTrkInSec++;
0751             }
0752             if (nTrkInSec > 0) {
0753               profDupTracksVsEta_[tName]->Fill(
0754                   std::abs(tp.eta()), nTrkInSec - 1);  // Study duplication of tracks within an individual HT array.
0755               profDupTracksVsInvPt_[tName]->Fill(
0756                   std::abs(tp.qOverPt()), nTrkInSec - 1);  // Study duplication of tracks within an individual HT array.
0757             }
0758           }
0759         }
0760       }
0761
0762       if (tpRecoed) {
0763         // Increment sum each time a TP is reconstructed at least once inside Tracker
0764         if (tp.useForEff())
0765           nRecoedTPsForEff++;
0766         nRecoedTPs++;
0767       }
0768     }
0769
0770     //--- Plot mean number of tracks/event, counting number due to different kinds of duplicates
0771
0772     // Plot number of TPs for the efficiency measurement that are reconstructed.
0773     profNumTrackCands_[tName]->Fill(7.0, nRecoedTPsForEff);
0774     // Plot number of TPs that are reconstructed.
0775     profNumTrackCands_[tName]->Fill(6.0, nRecoedTPs);
0776     // Plot number of TPs that are reconstructed. Count +1 for each eta sector they are reconstructed in.
0777     profNumTrackCands_[tName]->Fill(5.0, nEtaSecsMatchingTPs);
0778     // Plot number of TPs that are reconstructed. Count +1 for each (eta,phi) sector they are reconstructed in.
0779     profNumTrackCands_[tName]->Fill(4.0, nSecsMatchingTPs);
0780     // Plot number of TP that are reconstructed. Count +1 for each track they are reconstructed as.
0781     profNumTrackCands_[tName]->Fill(2.0, nTrksMatchingTPs);
0782
0783     // Histos of track helix params.
0784     for (const L1track3D& trk : tracks) {
0785       hisQoverPt_[tName]->Fill(trk.qOverPt());
0786       hisPhi0_[tName]->Fill(trk.phi0());
0787       hisEta_[tName]->Fill(trk.eta());
0788       hisZ0_[tName]->Fill(trk.z0());
0789     }
0790
0791     // Histos of track parameter resolution
0792
0793     for (const TP& tp : vTPs) {
0794       if ((resPlotOpt_ && tp.useForAlgEff()) ||
0795           (not resPlotOpt_)) {  // Check TP is good for efficiency measurement (& also comes from signal event if requested)
0796
0797         // For each tracking particle, find the corresponding reconstructed track(s).
0798         for (const L1track3D& trk : tracks) {
0799           const TP* tpAssoc = trk.matchedTP();
0800           if (tpAssoc != nullptr) {
0801             if (tpAssoc->index() == tp.index()) {
0802               hisQoverPtRes_[tName]->Fill(trk.qOverPt() - tp.qOverPt());
0803               hisPhi0Res_[tName]->Fill(reco::deltaPhi(trk.phi0(), tp.phi0()));
0804               hisEtaRes_[tName]->Fill(trk.eta() - tp.eta());
0805               hisZ0Res_[tName]->Fill(trk.z0() - tp.z0());
0806             }
0807           }
0808         }
0809       }
0810     }
0811
0812     //=== Study tracking efficiency by looping over tracking particles.
0813
0814     for (const TP& tp : vTPs) {
0815       if (tp.useForEff()) {  // Check TP is good for efficiency measurement.
0816
0817         // Check if this TP was reconstructed anywhere in the tracker..
0818         bool tpRecoed = false;
0819         bool tpRecoedPerfect = false;
0820         for (const L1track3D& trk : tracks) {
0821           const TP* tpAssoc = trk.matchedTP();
0822           if (tpAssoc != nullptr) {
0823             if (tpAssoc->index() == tp.index()) {
0824               tpRecoed = true;
0825               if (trk.purity() == 1.)
0826                 tpRecoedPerfect = true;
0827             }
0828           }
0829         }
0830
0831         // If TP was reconstucted by HT, then plot its kinematics.
0832         if (tpRecoed) {
0833           hisRecoTPinvptForEff_[tName]->Fill(1. / tp.pt());
0834           hisRecoTPetaForEff_[tName]->Fill(tp.eta());
0835           hisRecoTPphiForEff_[tName]->Fill(tp.phi0());
0836           // Plot also production point of all good reconstructed TP.
0837           hisRecoTPd0ForEff_[tName]->Fill(std::abs(tp.d0()));
0838           hisRecoTPz0ForEff_[tName]->Fill(std::abs(tp.z0()));
0839           // Also plot efficiency to perfectly reconstruct the track (no fake hits)
0840           if (tpRecoedPerfect) {
0841             hisPerfRecoTPinvptForEff_[tName]->Fill(1. / tp.pt());
0842             hisPerfRecoTPetaForEff_[tName]->Fill(tp.eta());
0843           }
0844           if (tp.useForAlgEff()) {  // Check TP is good for algorithmic efficiency measurement.
0845             hisRecoTPinvptForAlgEff_[tName]->Fill(1. / tp.pt());
0846             hisRecoTPetaForAlgEff_[tName]->Fill(tp.eta());
0847             hisRecoTPphiForAlgEff_[tName]->Fill(tp.phi0());
0848             // Plot also production point of all good reconstructed TP.
0849             hisRecoTPd0ForAlgEff_[tName]->Fill(std::abs(tp.d0()));
0850             hisRecoTPz0ForAlgEff_[tName]->Fill(std::abs(tp.z0()));
0851
0852             // Also plot efficiency to perfectly reconstruct the track (no fake hits)
0853             if (tpRecoedPerfect) {
0854               hisPerfRecoTPinvptForAlgEff_[tName]->Fill(1. / tp.pt());
0855               hisPerfRecoTPetaForAlgEff_[tName]->Fill(tp.eta());
0856             }
0857           }
0858         }
0859       }
0860     }
0861   }
0862
0863   //=== Book histograms for studying track fitting.
0864
0865   map<string, TFileDirectory> Histos::bookTrackFitting() {
0866     map<string, TFileDirectory> inputDirMap;
0867
0868     for (const string& fitName : trackFitters_) {
0869       // Define lambda function to facilitate adding "fitName" histogram names.
0870       auto addn = [fitName](const string& s) { return TString::Format("%s_%s", s.c_str(), fitName.c_str()); };
0871
0872       TFileDirectory inputDir = fs_->mkdir(fitName);
0873       inputDirMap[fitName] = inputDir;
0874
0875       profNumFitTracks_[fitName] =
0876           inputDir.make<TProfile>(addn("NumFitTracks"), "; class; # of fitted tracks", 11, 0.5, 11.5, -0.5, 9.9e6);
0877       profNumFitTracks_[fitName]->GetXaxis()->SetBinLabel(7, "TP for eff fitted");
0878       profNumFitTracks_[fitName]->GetXaxis()->SetBinLabel(6, "TP fitted");
0879       profNumFitTracks_[fitName]->GetXaxis()->SetBinLabel(2, "Fit tracks that are genuine");
0880       profNumFitTracks_[fitName]->GetXaxis()->SetBinLabel(1, "Fit tracks including fakes");
0881       profNumFitTracks_[fitName]->LabelsOption("d");
0882
0883       hisNumFitTrks_[fitName] =
0884           inputDir.make<TH1F>(addn("NumFitTrks"), "; No. fitted tracks in tracker;", 200, -0.5, 399.5);
0885       hisNumFitTrksPerNon_[fitName] =
0886           inputDir.make<TH1F>(addn("NumFitTrksPerNon"), "; No. fitted tracks per nonant;", 200, -0.5, 199.5);
0887
0888       hisStubsPerFitTrack_[fitName] =
0889           inputDir.make<TH1F>(addn("StubsPerFitTrack"), "; No. of stubs per fitted track", 20, -0.5, 19.5);
0890       profStubsOnFitTracks_[fitName] = inputDir.make<TProfile>(
0891           addn("StubsOnFitTracks"), "; ; No. of stubs on all fitted tracks per event", 1, 0.5, 1.5);
0892
0893       hisFitQinvPtMatched_[fitName] =
0894           inputDir.make<TH1F>(addn("FitQinvPtMatched"), "Fitted q/p_{T} for matched tracks", 120, -0.6, 0.6);
0895       hisFitPhi0Matched_[fitName] =
0896           inputDir.make<TH1F>(addn("FitPhi0Matched"), "Fitted #phi_{0} for matched tracks", 70, -3.5, 3.5);
0897       hisFitD0Matched_[fitName] =
0898           inputDir.make<TH1F>(addn("FitD0Matched"), "Fitted d_{0} for matched tracks", 100, -2., 2.);
0899       hisFitZ0Matched_[fitName] =
0900           inputDir.make<TH1F>(addn("FitZ0Matched"), "Fitted z_{0} for matched tracks", 100, -25., 25.);
0901       hisFitEtaMatched_[fitName] =
0902           inputDir.make<TH1F>(addn("FitEtaMatched"), "Fitted #eta for matched tracks", 70, -3.5, 3.5);
0903
0904       hisFitQinvPtUnmatched_[fitName] =
0905           inputDir.make<TH1F>(addn("FitQinvPtUnmatched"), "Fitted q/p_{T} for unmatched tracks", 120, -0.6, 0.6);
0906       hisFitPhi0Unmatched_[fitName] =
0907           inputDir.make<TH1F>(addn("FitPhi0Unmatched"), "Fitted #phi_{0} for unmatched tracks", 70, -3.5, 3.5);
0908       hisFitD0Unmatched_[fitName] =
0909           inputDir.make<TH1F>(addn("FitD0Unmatched"), "Fitted d_{0} for unmatched tracks", 100, -2., 2.);
0910       hisFitZ0Unmatched_[fitName] =
0911           inputDir.make<TH1F>(addn("FitZ0Unmatched"), "Fitted z_{0} for unmatched tracks", 100, -25., 25.);
0912       hisFitEtaUnmatched_[fitName] =
0913           inputDir.make<TH1F>(addn("FitEtaUnmatched"), "Fitted #eta for unmatched tracks", 70, -3.5, 3.5);
0914
0915       const unsigned int nBinsChi2 = 39;
0916       const float chi2dofBins[nBinsChi2 + 1] = {0.0,  0.2,  0.4,  0.6,   0.8,   1.0,   1.2,   1.4,   1.6,   1.8,
0917                                                 2.0,  2.4,  2.8,  3.2,   3.6,   4.0,   4.5,   5.0,   6.0,   7.0,
0918                                                 8.0,  9.0,  10.0, 12.0,  14.0,  16.0,  18.0,  20.0,  25.0,  30.0,
0919                                                 40.0, 50.0, 75.0, 100.0, 150.0, 200.0, 250.0, 350.0, 500.0, 1000.0};
0920
0921       hisFitChi2DofRphiMatched_[fitName] =
0922           inputDir.make<TH1F>(addn("FitChi2DofRphiMatched"), ";#chi^{2}rphi;", nBinsChi2, chi2dofBins);
0923       hisFitChi2DofRzMatched_[fitName] =
0924           inputDir.make<TH1F>(addn("FitChi2DofRzMatched"), ";#chi^{2}rz/DOF;", nBinsChi2, chi2dofBins);
0925       profFitChi2DofRphiVsInvPtMatched_[fitName] =
0926           inputDir.make<TProfile>(addn("FitChi2DofRphiVsInvPtMatched"), "; 1/p_{T}; Fit #chi^{2}rphi/dof", 25, 0., 0.5);
0927
0928       hisFitChi2DofRphiUnmatched_[fitName] =
0929           inputDir.make<TH1F>(addn("FitChi2DofRphiUnmatched"), ";#chi^{2}rphi/DOF;", nBinsChi2, chi2dofBins);
0930       hisFitChi2DofRzUnmatched_[fitName] =
0931           inputDir.make<TH1F>(addn("FitChi2DofRzUnmatched"), ";#chi^{2}rz/DOF;", nBinsChi2, chi2dofBins);
0932       profFitChi2DofRphiVsInvPtUnmatched_[fitName] = inputDir.make<TProfile>(
0933           addn("FitChi2DofRphiVsInvPtUnmatched"), "; 1/p_{T}; Fit #chi^{2}rphi/dof", 25, 0., 0.5);
0934
0935       // Monitoring specific track fit algorithms.
0936       if (fitName.find("KF") != string::npos) {
0937         hisKalmanNumUpdateCalls_[fitName] =
0938             inputDir.make<TH1F>(addn("KalmanNumUpdateCalls"), "; Calls to KF updator;", 100, -0.5, 99.5);
0939
0940         hisKalmanChi2DofSkipLay0Matched_[fitName] = inputDir.make<TH1F>(
0941             addn("KalmanChi2DofSkipLay0Matched"), ";#chi^{2} for nSkippedLayers = 0;", nBinsChi2, chi2dofBins);
0942         hisKalmanChi2DofSkipLay1Matched_[fitName] = inputDir.make<TH1F>(
0943             addn("KalmanChi2DofSkipLay1Matched"), ";#chi^{2} for nSkippedLayers = 1;", nBinsChi2, chi2dofBins);
0944         hisKalmanChi2DofSkipLay2Matched_[fitName] = inputDir.make<TH1F>(
0945             addn("KalmanChi2DofSkipLay2Matched"), ";#chi^{2} for nSkippedLayers = 2;", nBinsChi2, chi2dofBins);
0946         hisKalmanChi2DofSkipLay0Unmatched_[fitName] = inputDir.make<TH1F>(
0947             addn("KalmanChi2DofSkipLay0Unmatched"), ";#chi^{2} for nSkippedLayers = 0;", nBinsChi2, chi2dofBins);
0948         hisKalmanChi2DofSkipLay1Unmatched_[fitName] = inputDir.make<TH1F>(
0949             addn("KalmanChi2DofSkipLay1Unmatched"), ";#chi^{2} for nSkippedLayers = 1;", nBinsChi2, chi2dofBins);
0950         hisKalmanChi2DofSkipLay2Unmatched_[fitName] = inputDir.make<TH1F>(
0951             addn("KalmanChi2DofSkipLay2Unmatched"), ";#chi^{2} for nSkippedLayers = 2;", nBinsChi2, chi2dofBins);
0952       }
0953
0954       // Plots of helix param resolution.
0955
0956       hisQoverPtResVsTrueEta_[fitName] = inputDir.make<TProfile>(
0957           addn("QoverPtResVsTrueEta"), "q/p_{T} resolution; |#eta|; q/p_{T} resolution", 30, 0.0, 3.0);
0958       hisPhi0ResVsTrueEta_[fitName] = inputDir.make<TProfile>(
0959           addn("PhiResVsTrueEta"), "#phi_{0} resolution; |#eta|; #phi_{0} resolution", 30, 0.0, 3.0);
0960       hisEtaResVsTrueEta_[fitName] =
0961           inputDir.make<TProfile>(addn("EtaResVsTrueEta"), "#eta resolution; |#eta|; #eta resolution", 30, 0.0, 3.0);
0962       hisZ0ResVsTrueEta_[fitName] =
0963           inputDir.make<TProfile>(addn("Z0ResVsTrueEta"), "z_{0} resolution; |#eta|; z_{0} resolution", 30, 0.0, 3.0);
0964       hisD0ResVsTrueEta_[fitName] =
0965           inputDir.make<TProfile>(addn("D0ResVsTrueEta"), "d_{0} resolution; |#eta|; d_{0} resolution", 30, 0.0, 3.0);
0966
0967       hisQoverPtResVsTrueInvPt_[fitName] = inputDir.make<TProfile>(
0968           addn("QoverPtResVsTrueInvPt"), "q/p_{T} resolution; 1/p_{T}; q/p_{T} resolution", 25, 0.0, 0.5);
0969       hisPhi0ResVsTrueInvPt_[fitName] = inputDir.make<TProfile>(
0970           addn("PhiResVsTrueInvPt"), "#phi_{0} resolution; 1/p_{T}; #phi_{0} resolution", 25, 0.0, 0.5);
0971       hisEtaResVsTrueInvPt_[fitName] =
0972           inputDir.make<TProfile>(addn("EtaResVsTrueInvPt"), "#eta resolution; 1/p_{T}; #eta resolution", 25, 0.0, 0.5);
0973       hisZ0ResVsTrueInvPt_[fitName] = inputDir.make<TProfile>(
0974           addn("Z0ResVsTrueInvPt"), "z_{0} resolution; 1/p_{T}; z_{0} resolution", 25, 0.0, 0.5);
0975       hisD0ResVsTrueInvPt_[fitName] = inputDir.make<TProfile>(
0976           addn("D0ResVsTrueInvPt"), "d_{0} resolution; 1/p_{T}; d_{0} resolution", 25, 0.0, 0.5);
0977
0978       // Duplicate track histos.
0979       profDupFitTrksVsEta_[fitName] =
0980           inputDir.make<TProfile>(addn("DupFitTrksVsEta"), "; #eta; No. of duplicate tracks per TP", 12, 0., 3.);
0981       profDupFitTrksVsInvPt_[fitName] =
0982           inputDir.make<TProfile>(addn("DupFitTrksVsInvPt"), "; 1/Pt; No. of duplicate tracks per TP", 25, 0., 0.5);
0983
0984       // Histos for tracking efficiency vs. TP kinematics. (Binning must match similar histos in bookTrackCands()).
0985       hisFitTPinvptForEff_[fitName] =
0986           inputDir.make<TH1F>(addn("FitTPinvptForEff"), "; TP 1/Pt of fitted tracks (for effi.);", 50, 0., 0.5);
0987       hisFitTPetaForEff_[fitName] =
0988           inputDir.make<TH1F>(addn("FitTPetaForEff"), "; TP #eta of fitted tracks (for effi.);", 20, -3., 3.);
0989       hisFitTPphiForEff_[fitName] =
0990           inputDir.make<TH1F>(addn("FitTPphiForEff"), "; TP #phi of fitted tracks (for effi.);", 20, -M_PI, M_PI);
0991
0992       // Histo for efficiency to reconstruct track perfectly (no incorrect hits). (Binning must match similar histos in bookTrackCands()).
0993       hisPerfFitTPinvptForEff_[fitName] = inputDir.make<TH1F>(
0994           addn("PerfFitTPinvptForEff"), "; TP 1/Pt of fitted tracks (for perf. effi.);", 50, 0., 0.5);
0995       hisPerfFitTPetaForEff_[fitName] = inputDir.make<TH1F>(
0996           addn("PerfFitTPetaForEff"), "; TP #eta of fitted tracks (for perfect effi.);", 20, -3., 3.);
0997
0998       // Histos for tracking efficiency vs. TP production point. (Binning must match similar histos in bookTrackCands()).
0999       hisFitTPd0ForEff_[fitName] =
1000           inputDir.make<TH1F>(addn("FitTPd0ForEff"), "; TP d0 of fitted tracks (for effi.);", 40, 0., 4.);
1001       hisFitTPz0ForEff_[fitName] =
1002           inputDir.make<TH1F>(addn("FitTPz0ForEff"), "; TP z0 of fitted tracks (for effi.);", 50, 0., 25.);
1003
1004       // Histos for algorithmic tracking efficiency vs. TP kinematics. (Binning must match similar histos in bookTrackCands()).
1005       hisFitTPinvptForAlgEff_[fitName] =
1006           inputDir.make<TH1F>(addn("FitTPinvptForAlgEff"), "; TP 1/Pt of fitted tracks (for alg. effi.);", 50, 0., 0.5);
1007       hisFitTPetaForAlgEff_[fitName] =
1008           inputDir.make<TH1F>(addn("FitTPetaForAlgEff"), "; TP #eta of fitted tracks (for alg. effi.);", 20, -3., 3.);
1009       hisFitTPphiForAlgEff_[fitName] = inputDir.make<TH1F>(
1010           addn("FitTPphiForAlgEff"), "; TP #phi of fitted tracks (for alg. effi.);", 20, -M_PI, M_PI);
1011
1012       // Histo for efficiency to reconstruct track perfectly (no incorrect hits). (Binning must match similar histos in bookTrackCands()).
1013       hisPerfFitTPinvptForAlgEff_[fitName] = inputDir.make<TH1F>(
1014           addn("PerfFitTPinvptForAlgEff"), "; TP 1/Pt of fitted tracks (for perf. alg. effi.);", 50, 0., 0.5);
1015       hisPerfFitTPetaForAlgEff_[fitName] =
1016           inputDir.make<TH1F>(addn("PerfFitTPetaForAlgEff"), "; TP #eta (for perf. alg. effi.);", 20, -3., 3.);
1017
1018       // Histos for algorithmic tracking efficiency vs. TP production point. (Binning must match similar histos in bookTrackCands()).
1019       hisFitTPd0ForAlgEff_[fitName] =
1020           inputDir.make<TH1F>(addn("FitTPd0ForAlgEff"), "; TP d0 of fitted tracks (for alg. effi.);", 40, 0., 4.);
1021       hisFitTPz0ForAlgEff_[fitName] =
1022           inputDir.make<TH1F>(addn("FitTPz0ForAlgEff"), "; TP z0 of fitted tracks (for alg. effi.);", 50, 0., 25.);
1023     }
1024     return inputDirMap;
1025   }
1026
1027   //=== Fill histograms for studying track fitting.
1028
1029   void Histos::fillTrackFitting(const InputData& inputData,
1030                                 const map<string, list<const L1fittedTrack*>>& mapFinalTracks) {
1031     // Allow only one thread to run this function at a time
1032     static std::mutex myMutex;
1033     std::lock_guard<std::mutex> myGuard(myMutex);
1034
1035     const list<TP>& vTPs = inputData.getTPs();
1036
1037     // Loop over all the fitting algorithms we are trying.
1038     for (const string& fitName : trackFitters_) {
1039       const list<const L1fittedTrack*>& fittedTracks = mapFinalTracks.at(fitName);  // Get fitted tracks.
1040
1041       // Count tracks
1042       unsigned int nFitTracks = 0;
1043       unsigned int nFitsMatchingTP = 0;
1044
1045       const unsigned int numPhiNonants = settings_->numPhiNonants();
1046       vector<unsigned int> nFitTracksPerNonant(numPhiNonants, 0);
1047
1048       for (const L1fittedTrack* fitTrk : fittedTracks) {
1049         nFitTracks++;
1050
1051         // Get matched truth particle, if any.
1052         const TP* tp = fitTrk->matchedTP();
1053         if (tp != nullptr)
1054           nFitsMatchingTP++;
1055         // Count fitted tracks per nonant.
1056         unsigned int iNonant = (numPhiSectors_ > 0) ? floor(fitTrk->iPhiSec() * numPhiNonants / (numPhiSectors_))
1057                                                     : 0;  // phi nonant number
1058         nFitTracksPerNonant[iNonant]++;
1059       }
1060
1061       profNumFitTracks_[fitName]->Fill(1, nFitTracks);
1062       profNumFitTracks_[fitName]->Fill(2, nFitsMatchingTP);
1063
1064       hisNumFitTrks_[fitName]->Fill(nFitTracks);
1065       for (const unsigned int& num : nFitTracksPerNonant) {
1066         hisNumFitTrksPerNon_[fitName]->Fill(num);
1067       }
1068
1069       // Count stubs assigned to fitted tracks.
1070       unsigned int nTotStubs = 0;
1071       for (const L1fittedTrack* fitTrk : fittedTracks) {
1072         unsigned int nStubs = fitTrk->numStubs();
1073         hisStubsPerFitTrack_[fitName]->Fill(nStubs);
1074         nTotStubs += nStubs;
1075       }
1076       profStubsOnFitTracks_[fitName]->Fill(1., nTotStubs);
1077
1078       // Note truth particles that are successfully fitted. And which give rise to duplicate tracks.
1079
1080       map<const TP*, bool> tpRecoedMap;  // Note which truth particles were successfully fitted.
1081       map<const TP*, bool>
1082           tpPerfRecoedMap;  // Note which truth particles were successfully fitted with no incorrect hits.
1083       map<const TP*, unsigned int> tpRecoedDup;  // Note that this TP gave rise to duplicate tracks.
1084       for (const TP& tp : vTPs) {
1085         tpRecoedMap[&tp] = false;
1086         tpPerfRecoedMap[&tp] = false;
1087         unsigned int nMatch = 0;
1088         for (const L1fittedTrack* fitTrk : fittedTracks) {
1089           const TP* assocTP = fitTrk->matchedTP();  // Get the TP the fitted track matches to, if any.
1090           if (assocTP == &tp) {
1091             tpRecoedMap[&tp] = true;
1092             if (fitTrk->purity() == 1.)
1093               tpPerfRecoedMap[&tp] = true;
1094             nMatch++;
1095           }
1096         }
1097         tpRecoedDup[&tp] = nMatch;
1098       }
1099
1100       // Count truth particles that are successfully fitted.
1101
1102       unsigned int nFittedTPs = 0;
1103       unsigned int nFittedTPsForEff = 0;
1104       for (const TP& tp : vTPs) {
1105         if (tpRecoedMap[&tp]) {  // Was this truth particle successfully fitted?
1106           nFittedTPs++;
1107           if (tp.useForEff())
1108             nFittedTPsForEff++;
1109         }
1110       }
1111
1112       profNumFitTracks_[fitName]->Fill(6, nFittedTPs);
1113       profNumFitTracks_[fitName]->Fill(7, nFittedTPsForEff);
1114
1115       // Loop over fitted tracks again.
1116
1117       for (const L1fittedTrack* fitTrk : fittedTracks) {
1118         // Info for specific track fit algorithms.
1119         unsigned int nSkippedLayers = 0;
1120         unsigned int numUpdateCalls = 0;
1121         if (fitName.find("KF") != string::npos) {
1122           fitTrk->infoKF(nSkippedLayers, numUpdateCalls);
1123           hisKalmanNumUpdateCalls_[fitName]->Fill(numUpdateCalls);
1124         }
1125
1126         //--- Compare fitted tracks that match truth particles to those that don't.
1127
1128         // Get matched truth particle, if any.
1129         const TP* tp = fitTrk->matchedTP();
1130
1131         if (tp != nullptr) {
1132           hisFitQinvPtMatched_[fitName]->Fill(fitTrk->qOverPt());
1133           hisFitPhi0Matched_[fitName]->Fill(fitTrk->phi0());
1134           hisFitD0Matched_[fitName]->Fill(fitTrk->d0());
1135           hisFitZ0Matched_[fitName]->Fill(fitTrk->z0());
1136           hisFitEtaMatched_[fitName]->Fill(fitTrk->eta());
1137
1138           // Only plot matched chi2 for tracks with no incorrect stubs.
1139           if (fitTrk->purity() == 1.) {
1140             hisFitChi2DofRphiMatched_[fitName]->Fill(fitTrk->chi2rphi() / fitTrk->numDOFrphi());
1141             hisFitChi2DofRzMatched_[fitName]->Fill(fitTrk->chi2rz() / fitTrk->numDOFrz());
1142             profFitChi2DofRphiVsInvPtMatched_[fitName]->Fill(std::abs(fitTrk->qOverPt()),
1143                                                              (fitTrk->chi2rphi() / fitTrk->numDOFrphi()));
1144
1145             if (fitName.find("KF") != string::npos) {
1146               // No. of skipped layers on track during Kalman track fit.
1147               if (nSkippedLayers == 0) {
1148                 hisKalmanChi2DofSkipLay0Matched_[fitName]->Fill(fitTrk->chi2dof());
1149               } else if (nSkippedLayers == 1) {
1150                 hisKalmanChi2DofSkipLay1Matched_[fitName]->Fill(fitTrk->chi2dof());
1151               } else if (nSkippedLayers >= 2) {
1152                 hisKalmanChi2DofSkipLay2Matched_[fitName]->Fill(fitTrk->chi2dof());
1153               }
1154             }
1155           }
1156
1157         } else {
1158           hisFitQinvPtUnmatched_[fitName]->Fill(fitTrk->qOverPt());
1159           hisFitPhi0Unmatched_[fitName]->Fill(fitTrk->phi0());
1160           hisFitD0Unmatched_[fitName]->Fill(fitTrk->d0());
1161           hisFitZ0Unmatched_[fitName]->Fill(fitTrk->z0());
1162           hisFitEtaUnmatched_[fitName]->Fill(fitTrk->eta());
1163
1164           hisFitChi2DofRphiUnmatched_[fitName]->Fill(fitTrk->chi2rphi() / fitTrk->numDOFrphi());
1165           hisFitChi2DofRzUnmatched_[fitName]->Fill(fitTrk->chi2rz() / fitTrk->numDOFrz());
1166           profFitChi2DofRphiVsInvPtUnmatched_[fitName]->Fill(std::abs(fitTrk->qOverPt()),
1167                                                              (fitTrk->chi2rphi() / fitTrk->numDOFrphi()));
1168
1169           if (fitName.find("KF") != string::npos) {
1170             // No. of skipped layers on track during Kalman track fit.
1171             if (nSkippedLayers == 0) {
1172               hisKalmanChi2DofSkipLay0Unmatched_[fitName]->Fill(fitTrk->chi2dof());
1173             } else if (nSkippedLayers == 1) {
1174               hisKalmanChi2DofSkipLay1Unmatched_[fitName]->Fill(fitTrk->chi2dof());
1175             } else if (nSkippedLayers >= 2) {
1176               hisKalmanChi2DofSkipLay2Unmatched_[fitName]->Fill(fitTrk->chi2dof());
1177             }
1178           }
1179         }
1180       }
1181
1182       // Study helix param resolution.
1183
1184       for (const L1fittedTrack* fitTrk : fittedTracks) {
1185         const TP* tp = fitTrk->matchedTP();
1186         if (tp != nullptr) {
1187           // IRT
1188           if ((resPlotOpt_ && tp->useForAlgEff()) ||
1189               (not resPlotOpt_)) {  // Check TP is good for efficiency measurement (& also comes from signal event if requested)
1190
1191             // Plot helix parameter resolution against eta or Pt.
1192             hisQoverPtResVsTrueEta_[fitName]->Fill(std::abs(tp->eta()), std::abs(fitTrk->qOverPt() - tp->qOverPt()));
1193             hisPhi0ResVsTrueEta_[fitName]->Fill(std::abs(tp->eta()),
1194                                                 std::abs(reco::deltaPhi(fitTrk->phi0(), tp->phi0())));
1195             hisEtaResVsTrueEta_[fitName]->Fill(std::abs(tp->eta()), std::abs(fitTrk->eta() - tp->eta()));
1196             hisZ0ResVsTrueEta_[fitName]->Fill(std::abs(tp->eta()), std::abs(fitTrk->z0() - tp->z0()));
1197             hisD0ResVsTrueEta_[fitName]->Fill(std::abs(tp->eta()), std::abs(fitTrk->d0() - tp->d0()));
1198
1199             hisQoverPtResVsTrueInvPt_[fitName]->Fill(std::abs(tp->qOverPt()),
1200                                                      std::abs(fitTrk->qOverPt() - tp->qOverPt()));
1201             hisPhi0ResVsTrueInvPt_[fitName]->Fill(std::abs(tp->qOverPt()),
1202                                                   std::abs(reco::deltaPhi(fitTrk->phi0(), tp->phi0())));
1203             hisEtaResVsTrueInvPt_[fitName]->Fill(std::abs(tp->qOverPt()), std::abs(fitTrk->eta() - tp->eta()));
1204             hisZ0ResVsTrueInvPt_[fitName]->Fill(std::abs(tp->qOverPt()), std::abs(fitTrk->z0() - tp->z0()));
1205             hisD0ResVsTrueInvPt_[fitName]->Fill(std::abs(tp->qOverPt()), std::abs(fitTrk->d0() - tp->d0()));
1206           }
1207         }
1208       }
1209
1210       //=== Study duplicate tracks.
1211
1212       for (const TP& tp : vTPs) {
1213         if (tpRecoedMap[&tp]) {  // Was this truth particle successfully fitted?
1214           profDupFitTrksVsEta_[fitName]->Fill(std::abs(tp.eta()), tpRecoedDup[&tp] - 1);
1215           profDupFitTrksVsInvPt_[fitName]->Fill(std::abs(tp.qOverPt()), tpRecoedDup[&tp] - 1);
1216         }
1217       }
1218
1219       //=== Study tracking efficiency by looping over tracking particles.
1220
1221       for (const TP& tp : vTPs) {
1222         if (tp.useForEff()) {  // Check TP is good for efficiency measurement.
1223
1224           // If TP was reconstucted by HT, then plot its kinematics.
1225           if (tpRecoedMap[&tp]) {  // This truth particle was successfully fitted.
1226             hisFitTPinvptForEff_[fitName]->Fill(1. / tp.pt());
1227             hisFitTPetaForEff_[fitName]->Fill(tp.eta());
1228             hisFitTPphiForEff_[fitName]->Fill(tp.phi0());
1229             // Plot also production point of all good reconstructed TP.
1230             hisFitTPd0ForEff_[fitName]->Fill(std::abs(tp.d0()));
1231             hisFitTPz0ForEff_[fitName]->Fill(std::abs(tp.z0()));
1232             // Also plot efficiency to perfectly reconstruct the track (no fake hits)
1233             if (tpPerfRecoedMap[&tp]) {  // This truth particle was successfully fitted with no incorrect hits.
1234               hisPerfFitTPinvptForEff_[fitName]->Fill(1. / tp.pt());
1235               hisPerfFitTPetaForEff_[fitName]->Fill(tp.eta());
1236             }
1237             if (tp.useForAlgEff()) {  // Check TP is good for algorithmic efficiency measurement.
1238               hisFitTPinvptForAlgEff_[fitName]->Fill(1. / tp.pt());
1239               hisFitTPetaForAlgEff_[fitName]->Fill(tp.eta());
1240               hisFitTPphiForAlgEff_[fitName]->Fill(tp.phi0());
1241               // Plot also production point of all good reconstructed TP.
1242               hisFitTPd0ForAlgEff_[fitName]->Fill(std::abs(tp.d0()));
1243               hisFitTPz0ForAlgEff_[fitName]->Fill(std::abs(tp.z0()));
1244               // Also plot efficiency to perfectly reconstruct the track (no fake hits)
1245               if (tpPerfRecoedMap[&tp]) {
1246                 hisPerfFitTPinvptForAlgEff_[fitName]->Fill(1. / tp.pt());
1247                 hisPerfFitTPetaForAlgEff_[fitName]->Fill(tp.eta());
1248               }
1249             }
1250           }
1251         }
1252       }
1253     }
1254   }
1255
1256   //=== Produce plots of tracking efficiency after HT or after r-z track filter (run at end of job).
1257
1258   TFileDirectory Histos::plotTrackEfficiency(const string& tName) {
1259     // Define lambda function to facilitate adding "tName" to directory & histogram names.
1260     auto addn = [tName](const string& s) { return TString::Format("%s_%s", s.c_str(), tName.c_str()); };
1261
1262     TFileDirectory inputDir = fs_->mkdir(addn("Effi").Data());
1263     // Plot tracking efficiency
1264     makeEfficiencyPlot(inputDir,
1265                        teffEffVsInvPt_[tName],
1266                        hisRecoTPinvptForEff_[tName],
1267                        hisTPinvptForEff_,
1268                        addn("EffVsInvPt"),
1269                        "; 1/Pt; Tracking efficiency");
1270     makeEfficiencyPlot(inputDir,
1271                        teffEffVsEta_[tName],
1272                        hisRecoTPetaForEff_[tName],
1273                        hisTPetaForEff_,
1274                        addn("EffVsEta"),
1275                        "; #eta; Tracking efficiency");
1276
1277     makeEfficiencyPlot(inputDir,
1278                        teffEffVsPhi_[tName],
1279                        hisRecoTPphiForEff_[tName],
1280                        hisTPphiForEff_,
1281                        addn("EffVsPhi"),
1282                        "; #phi; Tracking efficiency");
1283
1284     makeEfficiencyPlot(inputDir,
1285                        teffEffVsD0_[tName],
1286                        hisRecoTPd0ForEff_[tName],
1287                        hisTPd0ForEff_,
1288                        addn("EffVsD0"),
1289                        "; d0 (cm); Tracking efficiency");
1290     makeEfficiencyPlot(inputDir,
1291                        teffEffVsZ0_[tName],
1292                        hisRecoTPz0ForEff_[tName],
1293                        hisTPz0ForEff_,
1294                        addn("EffVsZ0"),
1295                        "; z0 (cm); Tracking efficiency");
1296
1297     // Also plot efficiency to reconstruct track perfectly.
1298     makeEfficiencyPlot(inputDir,
1299                        teffPerfEffVsInvPt_[tName],
1300                        hisPerfRecoTPinvptForEff_[tName],
1301                        hisTPinvptForEff_,
1302                        addn("PerfEffVsInvPt"),
1303                        "; 1/Pt; Tracking perfect efficiency");
1304     makeEfficiencyPlot(inputDir,
1305                        teffPerfEffVsEta_[tName],
1306                        hisPerfRecoTPetaForEff_[tName],
1307                        hisTPetaForEff_,
1308                        addn("PerfEffVsEta"),
1309                        "; #eta; Tracking perfect efficiency");
1310
1311     // Plot algorithmic tracking efficiency
1312     makeEfficiencyPlot(inputDir,
1313                        teffAlgEffVsInvPt_[tName],
1314                        hisRecoTPinvptForAlgEff_[tName],
1315                        hisTPinvptForAlgEff_,
1316                        addn("AlgEffVsInvPt"),
1317                        "; 1/Pt; Tracking efficiency");
1318     makeEfficiencyPlot(inputDir,
1319                        teffAlgEffVsEta_[tName],
1320                        hisRecoTPetaForAlgEff_[tName],
1321                        hisTPetaForAlgEff_,
1322                        addn("AlgEffVsEta"),
1323                        "; #eta; Tracking efficiency");
1324     makeEfficiencyPlot(inputDir,
1325                        teffAlgEffVsPhi_[tName],
1326                        hisRecoTPphiForAlgEff_[tName],
1327                        hisTPphiForAlgEff_,
1328                        addn("AlgEffVsPhi"),
1329                        "; #phi; Tracking efficiency");
1330
1331     makeEfficiencyPlot(inputDir,
1332                        teffAlgEffVsD0_[tName],
1333                        hisRecoTPd0ForAlgEff_[tName],
1334                        hisTPd0ForAlgEff_,
1335                        addn("AlgEffVsD0"),
1336                        "; d0 (cm); Tracking efficiency");
1337     makeEfficiencyPlot(inputDir,
1338                        teffAlgEffVsZ0_[tName],
1339                        hisRecoTPz0ForAlgEff_[tName],
1340                        hisTPz0ForAlgEff_,
1341                        addn("AlgEffVsZ0"),
1342                        "; z0 (cm); Tracking efficiency");
1343
1344     // Also plot algorithmic efficiency to reconstruct track perfectly.
1345     makeEfficiencyPlot(inputDir,
1346                        teffPerfAlgEffVsInvPt_[tName],
1347                        hisPerfRecoTPinvptForAlgEff_[tName],
1348                        hisTPinvptForAlgEff_,
1349                        addn("PerfAlgEffVsInvPt"),
1350                        "; 1/Pt; Tracking perfect efficiency");
1351     makeEfficiencyPlot(inputDir,
1352                        teffPerfAlgEffVsEta_[tName],
1353                        hisPerfRecoTPetaForAlgEff_[tName],
1354                        hisTPetaForAlgEff_,
1355                        addn("PerfAlgEffVsEta"),
1356                        "; #eta; Tracking perfect efficiency");
1357
1358     return inputDir;
1359   }
1360
1361   //=== Produce plots of tracking efficiency after track fit (run at end of job).
1362
1363   TFileDirectory Histos::plotTrackEffAfterFit(const string& fitName) {
1364     // Define lambda function to facilitate adding "fitName" to directory & histogram names.
1365     auto addn = [fitName](const string& s) { return TString::Format("%s_%s", s.c_str(), fitName.c_str()); };
1366
1367     TFileDirectory inputDir = fs_->mkdir(addn("Effi").Data());
1368     // Plot tracking efficiency
1369     makeEfficiencyPlot(inputDir,
1370                        teffEffFitVsInvPt_[fitName],
1371                        hisFitTPinvptForEff_[fitName],
1372                        hisTPinvptForEff_,
1373                        addn("EffFitVsInvPt"),
1374                        "; 1/Pt; Tracking efficiency");
1375     makeEfficiencyPlot(inputDir,
1376                        teffEffFitVsEta_[fitName],
1377                        hisFitTPetaForEff_[fitName],
1378                        hisTPetaForEff_,
1379                        addn("EffFitVsEta"),
1380                        "; #eta; Tracking efficiency");
1381     makeEfficiencyPlot(inputDir,
1382                        teffEffFitVsPhi_[fitName],
1383                        hisFitTPphiForEff_[fitName],
1384                        hisTPphiForEff_,
1385                        addn("EffFitVsPhi"),
1386                        "; #phi; Tracking efficiency");
1387
1388     makeEfficiencyPlot(inputDir,
1389                        teffEffFitVsD0_[fitName],
1390                        hisFitTPd0ForEff_[fitName],
1391                        hisTPd0ForEff_,
1392                        addn("EffFitVsD0"),
1393                        "; d0 (cm); Tracking efficiency");
1394     makeEfficiencyPlot(inputDir,
1395                        teffEffFitVsZ0_[fitName],
1396                        hisFitTPz0ForEff_[fitName],
1397                        hisTPz0ForEff_,
1398                        addn("EffFitVsZ0"),
1399                        "; z0 (cm); Tracking efficiency");
1400
1401     // Also plot efficiency to reconstruct track perfectly.
1402     makeEfficiencyPlot(inputDir,
1403                        teffPerfEffFitVsInvPt_[fitName],
1404                        hisPerfFitTPinvptForEff_[fitName],
1405                        hisTPinvptForEff_,
1406                        addn("PerfEffFitVsInvPt"),
1407                        "; 1/Pt; Tracking perfect efficiency");
1408     makeEfficiencyPlot(inputDir,
1409                        teffPerfEffFitVsEta_[fitName],
1410                        hisPerfFitTPetaForEff_[fitName],
1411                        hisTPetaForEff_,
1412                        addn("PerfEffFitVsEta"),
1413                        "; #eta; Tracking perfect efficiency");
1414
1415     // Plot algorithmic tracking efficiency
1416     makeEfficiencyPlot(inputDir,
1417                        teffAlgEffFitVsInvPt_[fitName],
1418                        hisFitTPinvptForAlgEff_[fitName],
1419                        hisTPinvptForAlgEff_,
1420                        addn("AlgEffFitVsInvPt"),
1421                        "; 1/Pt; Tracking efficiency");
1422     makeEfficiencyPlot(inputDir,
1423                        teffAlgEffFitVsEta_[fitName],
1424                        hisFitTPetaForAlgEff_[fitName],
1425                        hisTPetaForAlgEff_,
1426                        addn("AlgEffFitVsEta"),
1427                        "; #eta; Tracking efficiency");
1428     makeEfficiencyPlot(inputDir,
1429                        teffAlgEffFitVsPhi_[fitName],
1430                        hisFitTPphiForAlgEff_[fitName],
1431                        hisTPphiForAlgEff_,
1432                        addn("AlgEffFitVsPhi"),
1433                        "; #phi; Tracking efficiency");
1434
1435     makeEfficiencyPlot(inputDir,
1436                        teffAlgEffFitVsD0_[fitName],
1437                        hisFitTPd0ForAlgEff_[fitName],
1438                        hisTPd0ForAlgEff_,
1439                        addn("AlgEffFitVsD0"),
1440                        "; d0 (cm); Tracking efficiency");
1441     makeEfficiencyPlot(inputDir,
1442                        teffAlgEffFitVsZ0_[fitName],
1443                        hisFitTPz0ForAlgEff_[fitName],
1444                        hisTPz0ForAlgEff_,
1445                        addn("AlgEffFitVsZ0"),
1446                        "; z0 (cm); Tracking efficiency");
1447
1448     // Also plot algorithmic efficiency to reconstruct track perfectly.
1449     makeEfficiencyPlot(inputDir,
1450                        teffPerfAlgEffFitVsInvPt_[fitName],
1451                        hisPerfFitTPinvptForAlgEff_[fitName],
1452                        hisTPinvptForAlgEff_,
1453                        addn("PerfAlgEffFitVsInvPt"),
1454                        "; 1/Pt; Tracking perfect efficiency");
1455     makeEfficiencyPlot(inputDir,
1456                        teffPerfAlgEffFitVsEta_[fitName],
1457                        hisPerfFitTPetaForAlgEff_[fitName],
1458                        hisTPetaForAlgEff_,
1459                        addn("Perf AlgEffFitVsEta"),
1460                        "; #eta; Tracking perfect efficiency");
1461     return inputDir;
1462   }
1463
1464   void Histos::makeEfficiencyPlot(
1465       TFileDirectory& inputDir, TEfficiency* outputEfficiency, TH1F* pass, TH1F* all, TString name, TString title) {
1466     outputEfficiency = inputDir.make<TEfficiency>(*pass, *all);
1467     outputEfficiency->SetName(name);
1468     outputEfficiency->SetTitle(title);
1469   }
1470
1471   //=== Print summary of track-finding performance after track pattern reco.
1472
1473   void Histos::printTrackPerformance(const string& tName) {
1474     float numTrackCands = profNumTrackCands_[tName]->GetBinContent(1);   // No. of track cands
1475     float numTrackCandsErr = profNumTrackCands_[tName]->GetBinError(1);  // No. of track cands uncertainty
1476     float numMatchedTrackCandsIncDups =
1477         profNumTrackCands_[tName]->GetBinContent(2);  // Ditto, counting only those matched to TP
1478     float numMatchedTrackCandsExcDups = profNumTrackCands_[tName]->GetBinContent(6);  // Ditto, but excluding duplicates
1479     float numFakeTracks = numTrackCands - numMatchedTrackCandsIncDups;
1480     float numExtraDupTracks = numMatchedTrackCandsIncDups - numMatchedTrackCandsExcDups;
1481     float fracFake = numFakeTracks / (numTrackCands + 1.0e-6);
1482     float fracDup = numExtraDupTracks / (numTrackCands + 1.0e-6);
1483
1484     float numStubsOnTracks = profStubsOnTracks_[tName]->GetBinContent(1);
1485     float meanStubsPerTrack =
1486         numStubsOnTracks / (numTrackCands + 1.0e-6);  //protection against demoninator equals zero.
1487     unsigned int numRecoTPforAlg = hisRecoTPinvptForAlgEff_[tName]->GetEntries();
1488     // Histograms of input truth particles (e.g. hisTPinvptForAlgEff_), used for denominator of efficiencies, are identical,
1489     // irrespective of whether made after HT or after r-z track filter, so always use the former.
1490     unsigned int numTPforAlg = hisTPinvptForAlgEff_->GetEntries();
1491     unsigned int numPerfRecoTPforAlg = hisPerfRecoTPinvptForAlgEff_[tName]->GetEntries();
1492     float algEff = float(numRecoTPforAlg) / (numTPforAlg + 1.0e-6);  //protection against demoninator equals zero.
1493     float algEffErr = sqrt(algEff * (1 - algEff) / (numTPforAlg + 1.0e-6));  // uncertainty
1494     float algPerfEff =
1495         float(numPerfRecoTPforAlg) / (numTPforAlg + 1.0e-6);  //protection against demoninator equals zero.
1496     float algPerfEffErr = sqrt(algPerfEff * (1 - algPerfEff) / (numTPforAlg + 1.0e-6));  // uncertainty
1497
1498     PrintL1trk() << "=========================================================================";
1499     if (tName == "HT") {
1500       PrintL1trk() << "               TRACK-FINDING SUMMARY AFTER HOUGH TRANSFORM             ";
1501     } else if (tName == "RZ") {
1502       PrintL1trk() << "               TRACK-FINDING SUMMARY AFTER R-Z TRACK FILTER            ";
1503     } else if (tName == "TRACKLET") {
1504       PrintL1trk() << "               TRACK-FINDING SUMMARY AFTER TRACKLET PATTERN RECO       ";
1505     }
1506     PrintL1trk() << "Number of track candidates found per event = " << numTrackCands << " +- " << numTrackCandsErr;
1507     PrintL1trk() << "                     with mean stubs/track = " << meanStubsPerTrack;
1508     PrintL1trk() << "Fraction of track cands that are fake = " << fracFake;
1509     PrintL1trk() << "Fraction of track cands that are genuine, but extra duplicates = " << fracDup;
1510
1511     PrintL1trk() << std::fixed << std::setprecision(4) << "Algorithmic tracking efficiency = " << numRecoTPforAlg << "/"
1512                  << numTPforAlg << " = " << algEff << " +- " << algEffErr;
1513     PrintL1trk() << "Perfect algorithmic tracking efficiency = " << numPerfRecoTPforAlg << "/" << numTPforAlg << " = "
1514                  << algPerfEff << " +- " << algPerfEffErr << " (no incorrect hits)";
1515   }
1516
1517   //=== Print summary of track-finding performance after helix fit for given track fitter.
1518
1519   void Histos::printFitTrackPerformance(const string& fitName) {
1520     float numFitTracks = profNumFitTracks_[fitName]->GetBinContent(1);   // No. of track cands
1521     float numFitTracksErr = profNumFitTracks_[fitName]->GetBinError(1);  // No. of track cands uncertainty
1522     float numMatchedFitTracksIncDups =
1523         profNumFitTracks_[fitName]->GetBinContent(2);  // Ditto, counting only those matched to TP
1524     float numMatchedFitTracksExcDups = profNumFitTracks_[fitName]->GetBinContent(6);  // Ditto, but excluding duplicates
1525     float numFakeFitTracks = numFitTracks - numMatchedFitTracksIncDups;
1526     float numExtraDupFitTracks = numMatchedFitTracksIncDups - numMatchedFitTracksExcDups;
1527     float fracFakeFit = numFakeFitTracks / (numFitTracks + 1.0e-6);
1528     float fracDupFit = numExtraDupFitTracks / (numFitTracks + 1.0e-6);
1529
1530     float numStubsOnFitTracks = profStubsOnFitTracks_[fitName]->GetBinContent(1);
1531     float meanStubsPerFitTrack =
1532         numStubsOnFitTracks / (numFitTracks + 1.0e-6);  //protection against demoninator equals zero.
1533     unsigned int numFitTPforAlg = hisFitTPinvptForAlgEff_[fitName]->GetEntries();
1534     // Histograms of input truth particles (e.g. hisTPinvptForAlgEff_), used for denominator of efficiencies, are identical,
1535     // irrespective of whether made after HT or after r-z track filter, so always use the former.
1536     unsigned int numTPforAlg = hisTPinvptForAlgEff_->GetEntries();
1537     unsigned int numPerfFitTPforAlg = hisPerfFitTPinvptForAlgEff_[fitName]->GetEntries();
1538     float fitEff = float(numFitTPforAlg) / (numTPforAlg + 1.0e-6);  //protection against demoninator equals zero.
1539     float fitEffErr = sqrt(fitEff * (1 - fitEff) / (numTPforAlg + 1.0e-6));  // uncertainty
1540     float fitPerfEff =
1541         float(numPerfFitTPforAlg) / (numTPforAlg + 1.0e-6);  //protection against demoninator equals zero.
1542     float fitPerfEffErr = sqrt(fitPerfEff * (1 - fitPerfEff) / (numTPforAlg + 1.0e-6));  // uncertainty
1543
1544     // Does this fitter require r-z track filter to be run before it?
1545     bool useRZfilt = (std::count(useRZfilter_.begin(), useRZfilter_.end(), fitName) > 0);
1546
1547     PrintL1trk() << "=========================================================================";
1548     PrintL1trk() << "                    TRACK FIT SUMMARY FOR: " << fitName;
1549     PrintL1trk() << "Number of fitted track candidates found per event = " << numFitTracks << " +- " << numFitTracksErr;
1550     PrintL1trk() << "                     with mean stubs/track = " << meanStubsPerFitTrack;
1551     PrintL1trk() << "Fraction of fitted tracks that are fake = " << fracFakeFit;
1552     PrintL1trk() << "Fraction of fitted tracks that are genuine, but extra duplicates = " << fracDupFit;
1553     PrintL1trk() << "Algorithmic fitting efficiency = " << numFitTPforAlg << "/" << numTPforAlg << " = " << fitEff
1554                  << " +- " << fitEffErr;
1555     PrintL1trk() << "Perfect algorithmic fitting efficiency = " << numPerfFitTPforAlg << "/" << numTPforAlg << " = "
1556                  << fitPerfEff << " +- " << fitPerfEffErr << " (no incorrect hits)";
1557     if (useRZfilt)
1558       PrintL1trk() << "(The above fitter used the '" << settings_->rzFilterName() << "' r-z track filter.)";
1559   }
1560
1561   //=== Print tracking performance summary & make tracking efficiency histograms.
1562
1563   void Histos::endJobAnalysis(const HTrphi::ErrorMonitor* htRphiErrMon) {
1564     // Don't bother producing summary if user didn't request histograms via TFileService in their cfg.
1565     if (not this->enabled())
1566       return;
1567
1568     // Protection when running in wierd mixed hybrid-TMTT modes.
1569     bool wierdMixedMode = (hisRecoTPinvptForEff_.find("TRACKLET") == hisRecoTPinvptForEff_.end());
1570
1571     if (settings_->hybrid() && not wierdMixedMode) {
1572       // Produce plots of tracking efficieny after tracklet pattern reco.
1573       this->plotTrackletSeedEfficiency();
1574       this->plotTrackEfficiency("TRACKLET");
1575       this->plotHybridDupRemovalEfficiency();
1576
1577     } else {
1578       // Produce plots of tracking efficiency using track candidates found after HT.
1579       this->plotTrackEfficiency("HT");
1580
1581       // Optionally produce plots of tracking efficiency using track candidates found after r-z track filter.
1582       if (ranRZfilter_)
1583         this->plotTrackEfficiency("RZ");
1584     }
1585
1586     // Produce more plots of tracking efficiency using track candidates after track fit.
1587     for (auto& fitName : trackFitters_) {
1588       this->plotTrackEffAfterFit(fitName);
1589     }
1590
1591     PrintL1trk() << "=========================================================================";
1592
1593     // Print r (z) range in which each barrel layer (endcap wheel) appears.
1594     // (Needed by firmware).
1595     PrintL1trk();
1596     PrintL1trk() << "--- r range in which stubs in each barrel layer appear ---";
1597     for (const auto& p : mapBarrelLayerMinR_) {
1598       unsigned int layer = p.first;
1599       PrintL1trk() << "   layer = " << layer << " : " << mapBarrelLayerMinR_[layer] << " < r < "
1600                    << mapBarrelLayerMaxR_[layer];
1601     }
1602     PrintL1trk() << "--- |z| range in which stubs in each endcap wheel appear ---";
1603     for (const auto& p : mapEndcapWheelMinZ_) {
1604       unsigned int layer = p.first;
1605       PrintL1trk() << "   wheel = " << layer << " : " << mapEndcapWheelMinZ_[layer] << " < |z| < "
1606                    << mapEndcapWheelMaxZ_[layer];
1607     }
1608
1609     // Print (r,|z|) range in which each module type (defined in DigitalStub) appears.
1610     // (Needed by firmware).
1611     PrintL1trk();
1612     PrintL1trk() << "--- (r,|z|) range in which each module type (defined in DigitalStub) appears ---";
1613     for (const auto& p : mapModuleTypeMinR_) {
1614       unsigned int modType = p.first;
1615       PrintL1trk() << "   Module type = " << modType << setprecision(1) << " : r range = ("
1616                    << mapModuleTypeMinR_[modType] << "," << mapModuleTypeMaxR_[modType] << "); z range = ("
1617                    << mapModuleTypeMinZ_[modType] << "," << mapModuleTypeMaxZ_[modType] << ")";
1618     }
1619     // Ugly bodge to allow for modules in barrel layers 1-2 & endcap wheels 3-5 being different.
1620     PrintL1trk() << "and in addition";
1621     for (const auto& p : mapExtraAModuleTypeMinR_) {
1622       unsigned int modType = p.first;
1623       PrintL1trk() << "   Module type = " << modType << setprecision(1) << " : r range = ("
1624                    << mapExtraAModuleTypeMinR_[modType] << "," << mapExtraAModuleTypeMaxR_[modType] << "); z range = ("
1625                    << mapExtraAModuleTypeMinZ_[modType] << "," << mapExtraAModuleTypeMaxZ_[modType] << ")";
1626     }
1627     PrintL1trk() << "and in addition";
1628     for (const auto& p : mapExtraBModuleTypeMinR_) {
1629       unsigned int modType = p.first;
1630       PrintL1trk() << "   Module type = " << modType << setprecision(1) << " : r range = ("
1631                    << mapExtraBModuleTypeMinR_[modType] << "," << mapExtraBModuleTypeMaxR_[modType] << "); z range = ("
1632                    << mapExtraBModuleTypeMinZ_[modType] << "," << mapExtraBModuleTypeMaxZ_[modType] << ")";
1633     }
1634     PrintL1trk() << "and in addition";
1635     for (const auto& p : mapExtraCModuleTypeMinR_) {
1636       unsigned int modType = p.first;
1637       PrintL1trk() << "   Module type = " << modType << setprecision(1) << " : r range = ("
1638                    << mapExtraCModuleTypeMinR_[modType] << "," << mapExtraCModuleTypeMaxR_[modType] << "); z range = ("
1639                    << mapExtraCModuleTypeMinZ_[modType] << "," << mapExtraCModuleTypeMaxZ_[modType] << ")";
1640     }
1641     PrintL1trk() << "and in addition";
1642     for (const auto& p : mapExtraDModuleTypeMinR_) {
1643       unsigned int modType = p.first;
1644       PrintL1trk() << "   Module type = " << modType << setprecision(1) << " : r range = ("
1645                    << mapExtraDModuleTypeMinR_[modType] << "," << mapExtraDModuleTypeMaxR_[modType] << "); z range = ("
1646                    << mapExtraDModuleTypeMinZ_[modType] << "," << mapExtraDModuleTypeMaxZ_[modType] << ")";
1647     }
1648     PrintL1trk();
1649
1650     if (settings_->hybrid() && not wierdMixedMode) {
1651       //--- Print summary of tracklet pattern reco
1652       this->printTrackletSeedFindingPerformance();
1653       this->printTrackPerformance("TRACKLET");
1654       this->printHybridDupRemovalPerformance();
1655     } else {
1656       //--- Print summary of track-finding performance after HT
1657       this->printTrackPerformance("HT");
1658       //--- Optionally print summary of track-finding performance after r-z track filter.
1659       if (ranRZfilter_)
1660         this->printTrackPerformance("RZ");
1661     }
1662
1663     //--- Print summary of track-finding performance after helix fit, for each track fitting algorithm used.
1664     for (const string& fitName : trackFitters_) {
1665       this->printFitTrackPerformance(fitName);
1666     }
1667     PrintL1trk() << "=========================================================================";
1668
1669     if (htRphiErrMon != nullptr && not settings_->hybrid()) {
1670       // Check that stub filling was consistent with known limitations of HT firmware design.
1671
1672       PrintL1trk() << "Max. |gradients| of stub lines in HT array is: r-phi = " << htRphiErrMon->maxLineGradient;
1673
1674       if (htRphiErrMon->maxLineGradient > 1.) {
1675         PrintL1trk()
1676             << "WARNING: Line |gradient| exceeds 1, which firmware will not be able to cope with! Please adjust HT "
1677                "array size to avoid this.";
1678
1679       } else if (htRphiErrMon->numErrorsTypeA > 0.) {
1680         float frac = float(htRphiErrMon->numErrorsTypeA) / float(htRphiErrMon->numErrorsNorm);
1681         PrintL1trk()
1682             << "WARNING: Despite line gradients being less than one, some fraction of HT columns have filled cells "
1683                "with no filled neighbours in W, SW or NW direction. Firmware will object to this! ";
1684         PrintL1trk() << "This fraction = " << frac << " for r-phi HT";
1685
1686       } else if (htRphiErrMon->numErrorsTypeB > 0.) {
1687         float frac = float(htRphiErrMon->numErrorsTypeB) / float(htRphiErrMon->numErrorsNorm);
1688         PrintL1trk()
1689             << "WARNING: Despite line gradients being less than one, some fraction of HT columns recorded individual "
1690                "stubs being added to more than two cells! Thomas firmware will object to this! ";
1691         PrintL1trk() << "This fraction = " << frac << " for r-phi HT";
1692       }
1693     }
1694
1695     // Check if GP B approximation cfg params are inconsistent.
1696     if (bApproxMistake_)
1697       PrintL1trk() << "\n WARNING: BApprox cfg params are inconsistent - see printout above.";
1698
1699     // Restore original ROOT default cfg.
1700     TH1::SetDefaultSumw2(oldSumW2opt_);
1701   }
1702
1703   //=== Determine "B" parameter, used in GP firmware to allow for tilted modules.
1704
1705   void Histos::trackerGeometryAnalysis(const list<TrackerModule>& listTrackerModule) {
1706     // Don't bother producing summary if user didn't request histograms via TFileService in their cfg.
1707     if (not this->enabled())
1708       return;
1709
1710     map<float, float> dataForGraph;
1711     for (const TrackerModule& trackerModule : listTrackerModule) {
1712       if (trackerModule.tiltedBarrel()) {
1713         float paramB = trackerModule.paramB();
1714         float zOverR = std::abs(trackerModule.minZ()) / trackerModule.minR();
1715         dataForGraph[paramB] = zOverR;  // Only store each unique paramB once, to get better histo weights.
1716       }
1717     }
1718
1719     const int nEntries = dataForGraph.size();
1720     vector<float> vecParamB(nEntries);
1721     vector<float> vecZoverR(nEntries);
1722     unsigned int i = 0;
1723     for (const auto& p : dataForGraph) {
1724       vecParamB[i] = p.first;
1725       vecZoverR[i] = p.second;
1726       i++;
1727     }
1728
1729     PrintL1trk() << "=========================================================================";
1730     PrintL1trk() << "--- Fit to cfg params for FPGA-friendly approximation to B parameter in GP & KF ---";
1731     PrintL1trk() << "--- (used to allowed for tilted barrel modules)                                 ---";
1732
1733     TFileDirectory inputDir = fs_->mkdir("InputDataB");
1734     graphBVsZoverR_ = inputDir.make<TGraph>(nEntries, &vecZoverR[0], &vecParamB[0]);
1735     graphBVsZoverR_->SetNameTitle("B vs module Z/R", "; Module Z/R; B");
1736     graphBVsZoverR_->Fit("pol1", "q");
1737     TF1* fittedFunction = graphBVsZoverR_->GetFunction("pol1");
1738     double gradient = fittedFunction->GetParameter(1);
1739     double intercept = fittedFunction->GetParameter(0);
1740     double gradient_err = fittedFunction->GetParError(1);
1741     double intercept_err = fittedFunction->GetParError(0);
1742     PrintL1trk() << "         BApprox_gradient (fitted)  = " << gradient << " +- " << gradient_err;
1743     PrintL1trk() << "         BApprox_intercept (fitted) = " << intercept << " +- " << intercept_err;
1744     PrintL1trk() << "=========================================================================";
1745
1746     // Check fitted params consistent with those assumed in cfg file.
1747     if (settings_->useApproxB()) {
1748       double gradientDiff = std::abs(gradient - settings_->bApprox_gradient());
1749       double interceptDiff = std::abs(intercept - settings_->bApprox_intercept());
1750       constexpr unsigned int nSigma = 5;
1751       if (gradientDiff > nSigma * gradient_err ||
1752           interceptDiff > nSigma * intercept_err) {  // Uncertainty independent of number of events
1753         PrintL1trk() << "\n WARNING: fitted parameters inconsistent with those specified in cfg file:";
1754         PrintL1trk() << "         BApprox_gradient  (cfg) = " << settings_->bApprox_gradient();
1755         PrintL1trk() << "         BApprox_intercept (cfg) = " << settings_->bApprox_intercept();
1756         bApproxMistake_ = true;  // Note that problem has occurred.
1757       }
1758     }
1759   }
1760
1761 }  // namespace tmtt