Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:06 

0001 #include "FWCore/Framework/interface/ConsumesCollector.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 
0004 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0005 
0006 #include <vector>
0007 #include <numeric>
0008 
0009 namespace deDxTools {
0010   using namespace std;
0011   using namespace reco;
0012 
0013   bool shapeSelection(const SiStripCluster& clus) {
0014     // ----------------  Counting the number of maxima   --------------------------
0015     //----------------------------------------------------------------------------
0016     auto const& ampls = clus.amplitudes();
0017     Int_t NofMax = 0;
0018     Int_t recur255 = 1;
0019     Int_t recur254 = 1;
0020     bool MaxOnStart = false;
0021     bool MaxInMiddle = false, MaxOnEnd = false;
0022     Int_t MaxPos = 0;
0023 
0024     // Start with max
0025     if (ampls.size() != 1 &&
0026         ((ampls[0] > ampls[1]) ||
0027          (ampls.size() > 2 && ampls[0] == ampls[1] && ampls[1] > ampls[2] && ampls[0] != 254 && ampls[0] != 255) ||
0028          (ampls.size() == 2 && ampls[0] == ampls[1] && ampls[0] != 254 && ampls[0] != 255))) {
0029       NofMax = NofMax + 1;
0030       MaxOnStart = true;
0031     }
0032 
0033     // Maximum reached
0034     if (ampls.size() > 2) {
0035       for (unsigned int i = 1; i < ampls.size() - 1U; i++) {
0036         if ((ampls[i] > ampls[i - 1] && ampls[i] > ampls[i + 1]) ||
0037             (ampls.size() > 3 && i > 0 && i < ampls.size() - 2U && ampls[i] == ampls[i + 1] &&
0038              ampls[i] > ampls[i - 1] && ampls[i] > ampls[i + 2] && ampls[i] != 254 && ampls[i] != 255)) {
0039           NofMax = NofMax + 1;
0040           MaxInMiddle = true;
0041           MaxPos = i;
0042         }
0043         if (ampls[i] == 255 && ampls[i] == ampls[i - 1]) {
0044           recur255 = recur255 + 1;
0045           MaxPos = i - (recur255 / 2);
0046           if (ampls[i] > ampls[i + 1]) {
0047             NofMax = NofMax + 1;
0048             MaxInMiddle = true;
0049           }
0050         }
0051         if (ampls[i] == 254 && ampls[i] == ampls[i - 1]) {
0052           recur254 = recur254 + 1;
0053           MaxPos = i - (recur254 / 2);
0054           if (ampls[i] > ampls[i + 1]) {
0055             NofMax = NofMax + 1;
0056             MaxInMiddle = true;
0057           }
0058         }
0059       }
0060     }
0061 
0062     // End with a max
0063     if (ampls.size() > 1) {
0064       if (ampls[ampls.size() - 1] > ampls[ampls.size() - 2] ||
0065           (ampls.size() > 2 && ampls[ampls.size() - 1] == ampls[ampls.size() - 2] &&
0066            ampls[ampls.size() - 2] > ampls[ampls.size() - 3]) ||
0067           ampls[ampls.size() - 1] == 255) {
0068         NofMax = NofMax + 1;
0069         MaxOnEnd = true;
0070       }
0071     }
0072 
0073     // If only one strip touched
0074     if (ampls.size() == 1) {
0075       NofMax = 1;
0076     }
0077 
0078     // --------------  SHAPE-BASED SELECTION FOR UNIQUE MAXIMAS --------------
0079     //------------------------------------------------------------------------
0080     /*
0081                  ____
0082                 |    |____
0083             ____|    |    |
0084            |    |    |    |____
0085        ____|    |    |    |    |
0086       |    |    |    |    |    |____
0087     __|____|____|____|____|____|____|__
0088     C_Mnn C_Mn C_M  C_D  C_Dn C_Dnn
0089     */
0090     bool shapecdtn = false;
0091 
0092     Float_t C_M = 0.0;
0093     Float_t C_D = 0.0;
0094     Float_t C_Mn = 10000;
0095     Float_t C_Dn = 10000;
0096     Float_t C_Mnn = 10000;
0097     Float_t C_Dnn = 10000;
0098     Int_t CDPos;
0099     Float_t coeff1 = 1.7;
0100     Float_t coeff2 = 2.0;
0101     Float_t coeffn = 0.10;
0102     Float_t coeffnn = 0.02;
0103     Float_t noise = 4.0;
0104 
0105     if (NofMax == 1) {
0106       if (MaxOnStart == true) {
0107         C_M = (Float_t)ampls[0];
0108         C_D = (Float_t)ampls[1];
0109         if (ampls.size() < 3)
0110           shapecdtn = true;
0111         else if (ampls.size() == 3) {
0112           C_Dn = (Float_t)ampls[2];
0113           if (C_Dn <= coeff1 * coeffn * C_D + coeff2 * coeffnn * C_M + 2 * noise || C_D == 255)
0114             shapecdtn = true;
0115         } else if (ampls.size() > 3) {
0116           C_Dn = (Float_t)ampls[2];
0117           C_Dnn = (Float_t)ampls[3];
0118           if ((C_Dn <= coeff1 * coeffn * C_D + coeff2 * coeffnn * C_M + 2 * noise || C_D == 255) &&
0119               C_Dnn <= coeff1 * coeffn * C_Dn + coeff2 * coeffnn * C_D + 2 * noise) {
0120             shapecdtn = true;
0121           }
0122         }
0123       }
0124 
0125       if (MaxOnEnd == true) {
0126         C_M = (Float_t)ampls[ampls.size() - 1];
0127         C_D = (Float_t)ampls[ampls.size() - 2];
0128         if (ampls.size() < 3)
0129           shapecdtn = true;
0130         else if (ampls.size() == 3) {
0131           C_Dn = (Float_t)ampls[0];
0132           if (C_Dn <= coeff1 * coeffn * C_D + coeff2 * coeffnn * C_M + 2 * noise || C_D == 255)
0133             shapecdtn = true;
0134         } else if (ampls.size() > 3) {
0135           C_Dn = (Float_t)ampls[ampls.size() - 3];
0136           C_Dnn = (Float_t)ampls[ampls.size() - 4];
0137           if ((C_Dn <= coeff1 * coeffn * C_D + coeff2 * coeffnn * C_M + 2 * noise || C_D == 255) &&
0138               C_Dnn <= coeff1 * coeffn * C_Dn + coeff2 * coeffnn * C_D + 2 * noise) {
0139             shapecdtn = true;
0140           }
0141         }
0142       }
0143 
0144       if (MaxInMiddle == true) {
0145         C_M = (Float_t)ampls[MaxPos];
0146         int LeftOfMaxPos = MaxPos - 1;
0147         if (LeftOfMaxPos <= 0)
0148           LeftOfMaxPos = 0;
0149         int RightOfMaxPos = MaxPos + 1;
0150         if (RightOfMaxPos >= (int)ampls.size())
0151           RightOfMaxPos = ampls.size() - 1;
0152         if (ampls[LeftOfMaxPos] < ampls[RightOfMaxPos]) {
0153           C_D = (Float_t)ampls[RightOfMaxPos];
0154           C_Mn = (Float_t)ampls[LeftOfMaxPos];
0155           CDPos = RightOfMaxPos;
0156         } else {
0157           C_D = (Float_t)ampls[LeftOfMaxPos];
0158           C_Mn = (Float_t)ampls[RightOfMaxPos];
0159           CDPos = LeftOfMaxPos;
0160         }
0161         if (C_Mn < coeff1 * coeffn * C_M + coeff2 * coeffnn * C_D + 2 * noise || C_M == 255) {
0162           if (ampls.size() == 3)
0163             shapecdtn = true;
0164           else if (ampls.size() > 3) {
0165             if (CDPos > MaxPos) {
0166               if (ampls.size() - CDPos - 1 == 0) {
0167                 C_Dn = 0;
0168                 C_Dnn = 0;
0169               }
0170               if (ampls.size() - CDPos - 1 == 1) {
0171                 C_Dn = (Float_t)ampls[CDPos + 1];
0172                 C_Dnn = 0;
0173               }
0174               if (ampls.size() - CDPos - 1 > 1) {
0175                 C_Dn = (Float_t)ampls[CDPos + 1];
0176                 C_Dnn = (Float_t)ampls[CDPos + 2];
0177               }
0178               if (MaxPos >= 2) {
0179                 C_Mnn = (Float_t)ampls[MaxPos - 2];
0180               } else if (MaxPos < 2)
0181                 C_Mnn = 0;
0182             }
0183             if (CDPos < MaxPos) {
0184               if (CDPos == 0) {
0185                 C_Dn = 0;
0186                 C_Dnn = 0;
0187               }
0188               if (CDPos == 1) {
0189                 C_Dn = (Float_t)ampls[0];
0190                 C_Dnn = 0;
0191               }
0192               if (CDPos > 1) {
0193                 C_Dn = (Float_t)ampls[CDPos - 1];
0194                 C_Dnn = (Float_t)ampls[CDPos - 2];
0195               }
0196               if (ampls.size() - LeftOfMaxPos > 1 && MaxPos + 2 < (int)(ampls.size()) - 1) {
0197                 C_Mnn = (Float_t)ampls[MaxPos + 2];
0198               } else
0199                 C_Mnn = 0;
0200             }
0201             if ((C_Dn <= coeff1 * coeffn * C_D + coeff2 * coeffnn * C_M + 2 * noise || C_D == 255) &&
0202                 C_Mnn <= coeff1 * coeffn * C_Mn + coeff2 * coeffnn * C_M + 2 * noise &&
0203                 C_Dnn <= coeff1 * coeffn * C_Dn + coeff2 * coeffnn * C_D + 2 * noise) {
0204               shapecdtn = true;
0205             }
0206           }
0207         }
0208       }
0209     }
0210     if (ampls.size() == 1) {
0211       shapecdtn = true;
0212     }
0213 
0214     return shapecdtn;
0215   }
0216 
0217   int getCharge(const SiStripCluster* cluster,
0218                 int& nSatStrip,
0219                 const GeomDetUnit& detUnit,
0220                 const std::vector<std::vector<float> >& calibGains,
0221                 const unsigned int& offsetDU_) {
0222     const auto& Ampls = cluster->amplitudes();
0223 
0224     nSatStrip = 0;
0225     int charge = 0;
0226 
0227     if (calibGains.empty()) {
0228       for (unsigned int i = 0; i < Ampls.size(); i++) {
0229         int calibratedCharge = Ampls[i];
0230         charge += calibratedCharge;
0231         if (calibratedCharge >= 254)
0232           nSatStrip++;
0233       }
0234     } else {
0235       for (unsigned int i = 0; i < Ampls.size(); i++) {
0236         int calibratedCharge = Ampls[i];
0237 
0238         auto& gains = calibGains[detUnit.index() - offsetDU_];
0239         calibratedCharge = (int)(calibratedCharge / gains[(cluster->firstStrip() + i) / 128]);
0240         if (calibratedCharge >= 255) {
0241           if (calibratedCharge >= 1025)
0242             calibratedCharge = 255;
0243           else
0244             calibratedCharge = 254;
0245         }
0246 
0247         charge += calibratedCharge;
0248         if (calibratedCharge >= 254)
0249           nSatStrip++;
0250       }
0251     }
0252     return charge;
0253   }
0254 
0255   void makeCalibrationMap(const std::string& m_calibrationPath,
0256                           const TrackerGeometry& tkGeom,
0257                           std::vector<std::vector<float> >& calibGains,
0258                           const unsigned int& offsetDU_) {
0259     auto const& dus = tkGeom.detUnits();
0260     calibGains.resize(dus.size() - offsetDU_);
0261 
0262     TChain* t1 = new TChain("SiStripCalib/APVGain");
0263     t1->Add(m_calibrationPath.c_str());
0264 
0265     unsigned int tree_DetId;
0266     unsigned char tree_APVId;
0267     double tree_Gain;
0268     t1->SetBranchAddress("DetId", &tree_DetId);
0269     t1->SetBranchAddress("APVId", &tree_APVId);
0270     t1->SetBranchAddress("Gain", &tree_Gain);
0271 
0272     for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
0273       t1->GetEntry(ientry);
0274       auto& gains = calibGains[tkGeom.idToDetUnit(DetId(tree_DetId))->index() - offsetDU_];
0275       if (gains.size() < (size_t)(tree_APVId + 1)) {
0276         gains.resize(tree_APVId + 1);
0277       }
0278       gains[tree_APVId] = tree_Gain;
0279     }
0280     t1->Delete();
0281   }
0282 
0283   ESGetTokenH3DDVariant esConsumes(std::string const& Record, edm::ConsumesCollector& iCC) {
0284     if (Record == "SiStripDeDxMip_3D_Rcd") {
0285       return iCC.esConsumes<H3DD, SiStripDeDxMip_3D_Rcd, edm::Transition::BeginRun>();
0286     }
0287     if (Record == "SiStripDeDxPion_3D_Rcd") {
0288       return iCC.esConsumes<H3DD, SiStripDeDxPion_3D_Rcd, edm::Transition::BeginRun>();
0289     }
0290     if (Record == "SiStripDeDxKaon_3D_Rcd") {
0291       return iCC.esConsumes<H3DD, SiStripDeDxKaon_3D_Rcd, edm::Transition::BeginRun>();
0292     }
0293     if (Record == "SiStripDeDxProton_3D_Rcd") {
0294       return iCC.esConsumes<H3DD, SiStripDeDxProton_3D_Rcd, edm::Transition::BeginRun>();
0295     }
0296     if (Record == "SiStripDeDxElectron_3D_Rcd") {
0297       return iCC.esConsumes<H3DD, SiStripDeDxElectron_3D_Rcd, edm::Transition::BeginRun>();
0298     }
0299     throw cms::Exception("WrongRecord for dEdx") << "The record : " << Record << "is unknown\n";
0300   }
0301 
0302   PhysicsTools::Calibration::HistogramD3D const& getHistogramD3D(edm::EventSetup const& iES,
0303                                                                  ESGetTokenH3DDVariant const& iToken) {
0304     switch (iToken.index()) {
0305       case 0:
0306         return iES.getData(std::get<0>(iToken));
0307       case 1:
0308         return iES.getData(std::get<1>(iToken));
0309       case 2:
0310         return iES.getData(std::get<2>(iToken));
0311       case 3:
0312         return iES.getData(std::get<3>(iToken));
0313       case 4:
0314         return iES.getData(std::get<4>(iToken));
0315     }
0316     throw cms::Exception("HistogramD3DTokenUnset");
0317   }
0318 
0319   void buildDiscrimMap(PhysicsTools::Calibration::HistogramD3D const& deDxMap,
0320                        std::string const& ProbabilityMode,
0321                        TH3F*& Prob_ChargePath) {
0322     float xmin = deDxMap.rangeX().min;
0323     float xmax = deDxMap.rangeX().max;
0324     float ymin = deDxMap.rangeY().min;
0325     float ymax = deDxMap.rangeY().max;
0326     float zmin = deDxMap.rangeZ().min;
0327     float zmax = deDxMap.rangeZ().max;
0328 
0329     if (Prob_ChargePath)
0330       delete Prob_ChargePath;
0331     Prob_ChargePath = new TH3F("Prob_ChargePath",
0332                                "Prob_ChargePath",
0333                                deDxMap.numberOfBinsX(),
0334                                xmin,
0335                                xmax,
0336                                deDxMap.numberOfBinsY(),
0337                                ymin,
0338                                ymax,
0339                                deDxMap.numberOfBinsZ(),
0340                                zmin,
0341                                zmax);
0342 
0343     if (strcmp(ProbabilityMode.c_str(), "Accumulation") == 0) {
0344       for (int i = 0; i <= Prob_ChargePath->GetXaxis()->GetNbins() + 1; i++) {
0345         for (int j = 0; j <= Prob_ChargePath->GetYaxis()->GetNbins() + 1; j++) {
0346           float Ni = 0;
0347           for (int k = 0; k <= Prob_ChargePath->GetZaxis()->GetNbins() + 1; k++) {
0348             Ni += deDxMap.binContent(i, j, k);
0349           }
0350           for (int k = 0; k <= Prob_ChargePath->GetZaxis()->GetNbins() + 1; k++) {
0351             float tmp = 0;
0352             for (int l = 0; l <= k; l++) {
0353               tmp += deDxMap.binContent(i, j, l);
0354             }
0355             if (Ni > 0) {
0356               Prob_ChargePath->SetBinContent(i, j, k, tmp / Ni);
0357             } else {
0358               Prob_ChargePath->SetBinContent(i, j, k, 0);
0359             }
0360           }
0361         }
0362       }
0363     } else if (strcmp(ProbabilityMode.c_str(), "Integral") == 0) {
0364       for (int i = 0; i <= Prob_ChargePath->GetXaxis()->GetNbins() + 1; i++) {
0365         for (int j = 0; j <= Prob_ChargePath->GetYaxis()->GetNbins() + 1; j++) {
0366           float Ni = 0;
0367           for (int k = 0; k <= Prob_ChargePath->GetZaxis()->GetNbins() + 1; k++) {
0368             Ni += deDxMap.binContent(i, j, k);
0369           }
0370           for (int k = 0; k <= Prob_ChargePath->GetZaxis()->GetNbins() + 1; k++) {
0371             float tmp = deDxMap.binContent(i, j, k);
0372             if (Ni > 0) {
0373               Prob_ChargePath->SetBinContent(i, j, k, tmp / Ni);
0374             } else {
0375               Prob_ChargePath->SetBinContent(i, j, k, 0);
0376             }
0377           }
0378         }
0379       }
0380     } else {
0381       throw cms::Exception("WrongConfig for dEdx") << "The ProbabilityMode: " << ProbabilityMode << "is unknown\n";
0382     }
0383   }
0384 
0385   bool isSpanningOver2APV(unsigned int FirstStrip, unsigned int ClusterSize) {
0386     if (FirstStrip == 0)
0387       return true;
0388     if (FirstStrip == 128)
0389       return true;
0390     if (FirstStrip == 256)
0391       return true;
0392     if (FirstStrip == 384)
0393       return true;
0394     if (FirstStrip == 512)
0395       return true;
0396     if (FirstStrip == 640)
0397       return true;
0398 
0399     if (FirstStrip <= 127 && FirstStrip + ClusterSize > 127)
0400       return true;
0401     if (FirstStrip <= 255 && FirstStrip + ClusterSize > 255)
0402       return true;
0403     if (FirstStrip <= 383 && FirstStrip + ClusterSize > 383)
0404       return true;
0405     if (FirstStrip <= 511 && FirstStrip + ClusterSize > 511)
0406       return true;
0407     if (FirstStrip <= 639 && FirstStrip + ClusterSize > 639)
0408       return true;
0409 
0410     if (FirstStrip + ClusterSize == 127)
0411       return true;
0412     if (FirstStrip + ClusterSize == 255)
0413       return true;
0414     if (FirstStrip + ClusterSize == 383)
0415       return true;
0416     if (FirstStrip + ClusterSize == 511)
0417       return true;
0418     if (FirstStrip + ClusterSize == 639)
0419       return true;
0420     if (FirstStrip + ClusterSize == 767)
0421       return true;
0422 
0423     return false;
0424   }
0425 
0426   bool isFarFromBorder(const TrajectoryStateOnSurface& trajState, const GeomDetUnit* it) {
0427     if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
0428       edm::LogInfo("deDxTools") << "this detID doesn't seem to belong to the Tracker" << std::endl;
0429       return false;
0430     }
0431 
0432     LocalPoint HitLocalPos = trajState.localPosition();
0433     LocalError HitLocalError = trajState.localError().positionError();
0434 
0435     const BoundPlane plane = it->surface();
0436     const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0437     const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
0438 
0439     if (!trapezoidalBounds && !rectangularBounds)
0440       return false;
0441 
0442     double DistFromBorder = 1.0;
0443     double HalfLength =
0444         trapezoidalBounds ? (*trapezoidalBounds).parameters()[3] : it->surface().bounds().length() / 2.0;
0445 
0446     if (fabs(HitLocalPos.y()) + HitLocalError.yy() >= (HalfLength - DistFromBorder))
0447       return false;
0448 
0449     return true;
0450   }
0451 
0452 }  // namespace deDxTools
[ 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!