Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:53 

0001 #ifndef SiStripHitEfficiencyHelpers_H
0002 #define SiStripHitEfficiencyHelpers_H
0003 
0004 // A bunch of helper functions to deal with menial tasks in the
0005 // hit efficiency computation for the PCL workflow
0006 
0007 // system includes
0008 #include <fmt/printf.h>
0009 #include <string>
0010 
0011 // user includes
0012 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0013 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0016 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0017 
0018 // ROOT includes
0019 #include "TEfficiency.h"
0020 #include "TProfile.h"
0021 #include "TString.h"
0022 
0023 namespace {
0024 
0025   enum bounds {
0026     k_LayersStart = 0,
0027     k_LayersAtTIBEnd = 4,
0028     k_LayersAtTOBEnd = 10,
0029     k_LayersAtTIDEnd = 13,
0030     k_LayersAtTECEnd = 22,
0031     k_END_OF_LAYERS = 23,
0032     k_END_OF_LAYS_AND_RINGS = 35
0033   };
0034 
0035   /*
0036    * for the trend plots of efficiency vs some variable
0037    */
0038   enum projections { k_vs_LUMI = 0, k_vs_PU = 1, k_vs_BX = 2, k_SIZE = 3 };
0039 
0040   const std::array<std::string, projections::k_SIZE> projFolder = {{"VsLumi", "VsPu", "VsBx"}};
0041   const std::array<std::string, projections::k_SIZE> projFoundHisto = {
0042       {"layerfound_vsLumi_layer_", "layerfound_vsPU_layer_", "foundVsBx_layer"}};
0043   const std::array<std::string, projections::k_SIZE> projTotalHisto = {
0044       {"layertotal_vsLumi_layer_", "layertotal_vsPU_layer_", "totalVsBx_layer"}};
0045   const std::array<std::string, projections::k_SIZE> projTitle = {{"Inst Lumi", "Pile-Up", "Bunch Crossing"}};
0046   const std::array<std::string, projections::k_SIZE> projXtitle = {
0047       {"instantaneous luminosity [Hz/cm^{2}]", "Pile-Up events", "Bunch Crossing number"}};
0048 
0049   inline void replaceInString(std::string& str, const std::string& from, const std::string& to) {
0050     if (from.empty())
0051       return;
0052     size_t start_pos = 0;
0053     while ((start_pos = str.find(from, start_pos)) != std::string::npos) {
0054       str.replace(start_pos, from.length(), to);
0055       start_pos += to.length();  // In case 'to' contains 'from', like replacing 'x' with 'yx'
0056     }
0057   }
0058 
0059   inline unsigned int checkLayer(unsigned int iidd, const TrackerTopology* tTopo) {
0060     switch (DetId(iidd).subdetId()) {
0061       case SiStripSubdetector::TIB:
0062         return tTopo->tibLayer(iidd);
0063       case SiStripSubdetector::TOB:
0064         return tTopo->tobLayer(iidd) + bounds::k_LayersAtTIBEnd;
0065       case SiStripSubdetector::TID:
0066         return tTopo->tidWheel(iidd) + bounds::k_LayersAtTOBEnd;
0067       case SiStripSubdetector::TEC:
0068         return tTopo->tecWheel(iidd) + bounds::k_LayersAtTIDEnd;
0069       default:
0070         return bounds::k_LayersStart;
0071     }
0072   }
0073 
0074   inline std::string layerName(unsigned int k, const bool showRings, const unsigned int nTEClayers) {
0075     const std::string ringlabel{showRings ? "R" : "D"};
0076     if (k > bounds::k_LayersStart && k <= bounds::k_LayersAtTIBEnd) {
0077       return fmt::format("TIB L{:d}", k);
0078     } else if (k > bounds::k_LayersAtTIBEnd && k <= bounds::k_LayersAtTOBEnd) {
0079       return fmt::format("TOB L{:d}", k - bounds::k_LayersAtTIBEnd);
0080     } else if (k > bounds::k_LayersAtTOBEnd && k <= bounds::k_LayersAtTIDEnd) {
0081       return fmt::format("TID {0}{1:d}", ringlabel, k - bounds::k_LayersAtTOBEnd);
0082     } else if (k > bounds::k_LayersAtTIDEnd && k <= bounds::k_LayersAtTIDEnd + nTEClayers) {
0083       return fmt::format("TEC {0}{1:d}", ringlabel, k - bounds::k_LayersAtTIDEnd);
0084     } else {
0085       return "should never be here!";
0086     }
0087   }
0088 
0089   inline std::string layerSideName(Long_t k, const bool showRings, const unsigned int nTEClayers) {
0090     const std::string ringlabel{showRings ? "R" : "D"};
0091     if (k > bounds::k_LayersStart && k <= bounds::k_LayersAtTIBEnd) {
0092       return fmt::format("TIB L{:d}", k);
0093     } else if (k > bounds::k_LayersAtTIBEnd && k <= bounds::k_LayersAtTOBEnd) {
0094       return fmt::format("TOB L{:d}", k - bounds::k_LayersAtTIBEnd);
0095     } else if (k > bounds::k_LayersAtTOBEnd && k < 14) {
0096       return fmt::format("TID- {0}{1:d}", ringlabel, k - bounds::k_LayersAtTOBEnd);
0097     } else if (k > 13 && k < 17) {
0098       return fmt::format("TID+ {0}{1:d}", ringlabel, k - 13);
0099     } else if (k > 16 && k < 17 + nTEClayers) {
0100       return fmt::format("TEC- {0}{1:d}", ringlabel, k - 16);
0101     } else if (k > 16 + nTEClayers) {
0102       return fmt::format("TEC+ {0}{1:d}", ringlabel, k - 16 - nTEClayers);
0103     } else {
0104       return "shoud never be here!";
0105     }
0106   }
0107 
0108   inline double checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters,
0109                                  double xx,
0110                                  double xerr) {
0111     double error = sqrt(parameters.second.xx() + xerr * xerr);
0112     double separation = abs(parameters.first.x() - xx);
0113     double consistency = separation / error;
0114     return consistency;
0115   }
0116 
0117   inline bool isDoubleSided(unsigned int iidd, const TrackerTopology* tTopo) {
0118     unsigned int layer;
0119     switch (DetId(iidd).subdetId()) {
0120       case SiStripSubdetector::TIB:
0121         layer = tTopo->tibLayer(iidd);
0122         return (layer == 1 || layer == 2);
0123       case SiStripSubdetector::TOB:
0124         layer = tTopo->tobLayer(iidd) + 4;
0125         return (layer == 5 || layer == 6);
0126       case SiStripSubdetector::TID:
0127         layer = tTopo->tidRing(iidd) + 10;
0128         return (layer == 11 || layer == 12);
0129       case SiStripSubdetector::TEC:
0130         layer = tTopo->tecRing(iidd) + 13;
0131         return (layer == 14 || layer == 15 || layer == 18);
0132       default:
0133         return false;
0134     }
0135   }
0136 
0137   inline bool check2DPartner(unsigned int iidd, const std::vector<TrajectoryMeasurement>& traj) {
0138     unsigned int partner_iidd = 0;
0139     bool found2DPartner = false;
0140     // first get the id of the other detector
0141     if ((iidd & 0x3) == 1)
0142       partner_iidd = iidd + 1;
0143     if ((iidd & 0x3) == 2)
0144       partner_iidd = iidd - 1;
0145     // next look in the trajectory measurements for a measurement from that detector
0146     // loop through trajectory measurements to find the partner_iidd
0147     for (const auto& tm : traj) {
0148       if (tm.recHit()->geographicalId().rawId() == partner_iidd) {
0149         found2DPartner = true;
0150       }
0151     }
0152     return found2DPartner;
0153   }
0154 
0155   inline bool isInBondingExclusionZone(
0156       unsigned int iidd, unsigned int TKlayers, double yloc, double yErr, const TrackerTopology* tTopo) {
0157     constexpr float exclusionWidth = 0.4;
0158     constexpr float TOBexclusion = 0.0;
0159     constexpr float TECexRing5 = -0.89;
0160     constexpr float TECexRing6 = -0.56;
0161     constexpr float TECexRing7 = 0.60;
0162 
0163     //Added by Chris Edelmaier to do TEC bonding exclusion
0164     const int subdetector = ((iidd >> 25) & 0x7);
0165     const int ringnumber = ((iidd >> 5) & 0x7);
0166 
0167     bool inZone = false;
0168     //New TOB and TEC bonding region exclusion zone
0169     if ((TKlayers >= 5 && TKlayers < 11) || ((subdetector == 6) && ((ringnumber >= 5) && (ringnumber <= 7)))) {
0170       //There are only 2 cases that we need to exclude for
0171       float highzone = 0.0;
0172       float lowzone = 0.0;
0173       float higherr = yloc + 5.0 * yErr;
0174       float lowerr = yloc - 5.0 * yErr;
0175       if (TKlayers >= 5 && TKlayers < 11) {
0176         //TOB zone
0177         highzone = TOBexclusion + exclusionWidth;
0178         lowzone = TOBexclusion - exclusionWidth;
0179       } else if (ringnumber == 5) {
0180         //TEC ring 5
0181         highzone = TECexRing5 + exclusionWidth;
0182         lowzone = TECexRing5 - exclusionWidth;
0183       } else if (ringnumber == 6) {
0184         //TEC ring 6
0185         highzone = TECexRing6 + exclusionWidth;
0186         lowzone = TECexRing6 - exclusionWidth;
0187       } else if (ringnumber == 7) {
0188         //TEC ring 7
0189         highzone = TECexRing7 + exclusionWidth;
0190         lowzone = TECexRing7 - exclusionWidth;
0191       }
0192       //Now that we have our exclusion region, we just have to properly identify it
0193       if ((highzone <= higherr) && (highzone >= lowerr))
0194         inZone = true;
0195       if ((lowzone >= lowerr) && (lowzone <= higherr))
0196         inZone = true;
0197       if ((higherr <= highzone) && (higherr >= lowzone))
0198         inZone = true;
0199       if ((lowerr >= lowzone) && (lowerr <= highzone))
0200         inZone = true;
0201     }
0202     return inZone;
0203   }
0204 
0205   struct ClusterInfo {
0206     float xResidual;
0207     float xResidualPull;
0208     float xLocal;
0209     ClusterInfo(float xRes, float xResPull, float xLoc) : xResidual(xRes), xResidualPull(xResPull), xLocal(xLoc) {}
0210   };
0211 
0212   inline float calcPhi(float x, float y) {
0213     float phi = 0;
0214     if ((x >= 0) && (y >= 0))
0215       phi = std::atan(y / x);
0216     else if ((x >= 0) && (y <= 0))
0217       phi = std::atan(y / x) + 2 * M_PI;
0218     else if ((x <= 0) && (y >= 0))
0219       phi = std::atan(y / x) + M_PI;
0220     else
0221       phi = std::atan(y / x) + M_PI;
0222     phi = phi * 180.0 / M_PI;
0223 
0224     return phi;
0225   }
0226 
0227   inline TProfile* computeEff(const TH1F* num, const TH1F* denum, const std::string nameHist) {
0228     std::string name = "eff_" + nameHist;
0229     std::string title = "SiStrip Hit Efficiency" + std::string(num->GetTitle());
0230     TProfile* efficHist = new TProfile(name.c_str(),
0231                                        title.c_str(),
0232                                        denum->GetXaxis()->GetNbins(),
0233                                        denum->GetXaxis()->GetXmin(),
0234                                        denum->GetXaxis()->GetXmax());
0235 
0236     for (int i = 1; i <= denum->GetNbinsX(); i++) {
0237       double nNum = num->GetBinContent(i);
0238       double nDenum = denum->GetBinContent(i);
0239       if (nDenum == 0 || nNum == 0) {
0240         continue;
0241       }
0242       if (nNum > nDenum) {
0243         edm::LogWarning("SiStripHitEfficiencyHelpers")
0244             << "Alert! specific bin's num is bigger than denum " << i << " " << nNum << " " << nDenum;
0245         nNum = nDenum;  // set the efficiency to 1
0246       }
0247       const double effVal = nNum / nDenum;
0248       efficHist->SetBinContent(i, effVal);
0249       efficHist->SetBinEntries(i, 1);
0250       const double errLo = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, false);
0251       const double errUp = TEfficiency::ClopperPearson((int)nDenum, (int)nNum, 0.683, true);
0252       const double errVal = (effVal - errLo > errUp - effVal) ? effVal - errLo : errUp - effVal;
0253       efficHist->SetBinError(i, sqrt(effVal * effVal + errVal * errVal));
0254 
0255       LogDebug("SiStripHitEfficiencyHelpers") << __PRETTY_FUNCTION__ << " " << nameHist << " bin:" << i
0256                                               << " err:" << sqrt(effVal * effVal + errVal * errVal);
0257     }
0258     return efficHist;
0259   }
0260 
0261   inline Measurement1D computeCPEfficiency(const double num, const double den) {
0262     if (den > 0) {
0263       const double effVal = num / den;
0264       const double errLo = TEfficiency::ClopperPearson((int)den, (int)num, 0.683, false);
0265       const double errUp = TEfficiency::ClopperPearson((int)den, (int)num, 0.683, true);
0266       const double errVal = (effVal - errLo > errUp - effVal) ? effVal - errLo : errUp - effVal;
0267       return Measurement1D(effVal, errVal);
0268     } else {
0269       return Measurement1D(0., 0.);
0270     }
0271   }
0272 }  // namespace
0273 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!