Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // -*- C++ -*-
0002 //
0003 // Package:    DeDxDiscriminatorLearner
0004 // Class:      DeDxDiscriminatorLearner
0005 //
0006 /**\class DeDxDiscriminatorLearner DeDxDiscriminatorLearner.cc RecoTracker/DeDxDiscriminatorLearner/src/DeDxDiscriminatorLearner.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Loic Quertenmont(querten)
0015 //         Created:  Mon October 20 10:09:02 CEST 2008
0016 //
0017 
0018 // system include files
0019 #include <memory>
0020 
0021 #include "CalibTracker/SiStripChannelGain/plugins/DeDxDiscriminatorLearner.h"
0022 
0023 #include "DataFormats/TrackReco/interface/Track.h"
0024 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0025 
0026 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0029 
0030 using namespace reco;
0031 using namespace std;
0032 using namespace edm;
0033 
0034 DeDxDiscriminatorLearner::DeDxDiscriminatorLearner(const edm::ParameterSet& iConfig)
0035     : ConditionDBWriter<PhysicsTools::Calibration::HistogramD3D>(iConfig) {
0036   m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0037   m_trajTrackAssociationTag =
0038       consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation"));
0039   m_tkGeomToken = esConsumes<edm::Transition::BeginRun>();
0040 
0041   P_Min = iConfig.getParameter<double>("P_Min");
0042   P_Max = iConfig.getParameter<double>("P_Max");
0043   P_NBins = iConfig.getParameter<int>("P_NBins");
0044   Path_Min = iConfig.getParameter<double>("Path_Min");
0045   Path_Max = iConfig.getParameter<double>("Path_Max");
0046   Path_NBins = iConfig.getParameter<int>("Path_NBins");
0047   Charge_Min = iConfig.getParameter<double>("Charge_Min");
0048   Charge_Max = iConfig.getParameter<double>("Charge_Max");
0049   Charge_NBins = iConfig.getParameter<int>("Charge_NBins");
0050 
0051   MinTrackMomentum = iConfig.getUntrackedParameter<double>("minTrackMomentum", 5.0);
0052   MaxTrackMomentum = iConfig.getUntrackedParameter<double>("maxTrackMomentum", 99999.0);
0053   MinTrackEta = iConfig.getUntrackedParameter<double>("minTrackEta", -5.0);
0054   MaxTrackEta = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0);
0055   MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips", 255);
0056   MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits", 4);
0057 
0058   algoMode = iConfig.getUntrackedParameter<string>("AlgoMode", "SingleJob");
0059   HistoFile = iConfig.getUntrackedParameter<string>("HistoFile", "out.root");
0060   VInputFiles = iConfig.getUntrackedParameter<vector<string> >("InputFiles");
0061 
0062   shapetest = iConfig.getParameter<bool>("ShapeTest");
0063   useCalibration = iConfig.getUntrackedParameter<bool>("UseCalibration");
0064   m_calibrationPath = iConfig.getUntrackedParameter<string>("calibrationPath");
0065 }
0066 
0067 DeDxDiscriminatorLearner::~DeDxDiscriminatorLearner() {}
0068 
0069 // ------------ method called once each job just before starting event loop  ------------
0070 
0071 void DeDxDiscriminatorLearner::algoBeginJob(const edm::EventSetup& iSetup) {
0072   Charge_Vs_Path = new TH3F("Charge_Vs_Path",
0073                             "Charge_Vs_Path",
0074                             P_NBins,
0075                             P_Min,
0076                             P_Max,
0077                             Path_NBins,
0078                             Path_Min,
0079                             Path_Max,
0080                             Charge_NBins,
0081                             Charge_Min,
0082                             Charge_Max);
0083 
0084   if (useCalibration && calibGains.empty()) {
0085     const auto& tkGeom = iSetup.getData(m_tkGeomToken);
0086 
0087     m_off = tkGeom.offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel
0088 
0089     deDxTools::makeCalibrationMap(m_calibrationPath, tkGeom, calibGains, m_off);
0090   }
0091 
0092   //Read the calibTree if in calibTree mode
0093   if (strcmp(algoMode.c_str(), "CalibTree") == 0)
0094     algoAnalyzeTheTree(iSetup);
0095 }
0096 
0097 // ------------ method called once each job just after ending the event loop  ------------
0098 
0099 void DeDxDiscriminatorLearner::algoEndJob() {
0100   if (strcmp(algoMode.c_str(), "MultiJob") == 0) {
0101     TFile* Output = new TFile(HistoFile.c_str(), "RECREATE");
0102     Charge_Vs_Path->Write();
0103     Output->Write();
0104     Output->Close();
0105   } else if (strcmp(algoMode.c_str(), "WriteOnDB") == 0) {
0106     TFile* Input = new TFile(HistoFile.c_str());
0107     Charge_Vs_Path = (TH3F*)(Input->FindObjectAny("Charge_Vs_Path"))->Clone();
0108     Input->Close();
0109   } else if (strcmp(algoMode.c_str(), "CalibTree") == 0) {
0110     TFile* Output = new TFile(HistoFile.c_str(), "RECREATE");
0111     Charge_Vs_Path->Write();
0112     Output->Write();
0113     Output->Close();
0114     TFile* Input = new TFile(HistoFile.c_str());
0115     Charge_Vs_Path = (TH3F*)(Input->FindObjectAny("Charge_Vs_Path"))->Clone();
0116     Input->Close();
0117   }
0118 }
0119 
0120 void DeDxDiscriminatorLearner::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0121   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
0122   iEvent.getByToken(m_trajTrackAssociationTag, trajTrackAssociationHandle);
0123 
0124   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0125   iEvent.getByToken(m_tracksTag, trackCollectionHandle);
0126 
0127   unsigned track_index = 0;
0128   for (TrajTrackAssociationCollection::const_iterator it = trajTrackAssociationHandle->begin();
0129        it != trajTrackAssociationHandle->end();
0130        ++it, track_index++) {
0131     const Track& track = *it->val;
0132     const Trajectory& traj = *it->key;
0133 
0134     if (track.eta() < MinTrackEta || track.eta() > MaxTrackEta) {
0135       continue;
0136     }
0137     if (track.pt() < MinTrackMomentum || track.pt() > MaxTrackMomentum) {
0138       continue;
0139     }
0140     if (track.found() < MinTrackHits) {
0141       continue;
0142     }
0143 
0144     const vector<TrajectoryMeasurement>& measurements = traj.measurements();
0145     for (vector<TrajectoryMeasurement>::const_iterator it = measurements.begin(); it != measurements.end(); it++) {
0146       TrajectoryStateOnSurface trajState = it->updatedState();
0147       if (!trajState.isValid())
0148         continue;
0149 
0150       const TrackingRecHit* recHit = (*it->recHit()).hit();
0151       if (!recHit || !recHit->isValid())
0152         continue;
0153       LocalVector trackDirection = trajState.localDirection();
0154       float cosine = trackDirection.z() / trackDirection.mag();
0155 
0156       processHit(recHit, trajState.localMomentum().mag(), cosine, trajState);
0157     }
0158   }
0159 }
0160 
0161 void DeDxDiscriminatorLearner::processHit(const TrackingRecHit* recHit,
0162                                           float trackMomentum,
0163                                           float& cosine,
0164                                           const TrajectoryStateOnSurface& trajState) {
0165   auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
0166   if (!thit.isValid())
0167     return;
0168 
0169   auto const& clus = thit.firstClusterRef();
0170   if (!clus.isValid())
0171     return;
0172 
0173   int NSaturating = 0;
0174   if (clus.isPixel()) {
0175     return;
0176   } else if (clus.isStrip() && !thit.isMatched()) {
0177     auto& detUnit = *(recHit->detUnit());
0178     auto& cluster = clus.stripCluster();
0179     if (cluster.amplitudes().size() > MaxNrStrips) {
0180       return;
0181     }
0182     if (deDxTools::isSpanningOver2APV(cluster.firstStrip(), cluster.amplitudes().size())) {
0183       return;
0184     }
0185     if (!deDxTools::isFarFromBorder(trajState, &detUnit)) {
0186       return;
0187     }
0188     float pathLen = 10.0 * detUnit.surface().bounds().thickness() / fabs(cosine);
0189     float chargeAbs = deDxTools::getCharge(&cluster, NSaturating, detUnit, calibGains, m_off);
0190     float charge = chargeAbs / pathLen;
0191     if (!shapetest || (shapetest && deDxTools::shapeSelection(cluster))) {
0192       Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
0193     }
0194   } else if (clus.isStrip() && thit.isMatched()) {
0195     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0196     if (!matchedHit)
0197       return;
0198     const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(matchedHit->det());
0199 
0200     auto& detUnitM = *(gdet->monoDet());
0201     auto& clusterM = matchedHit->monoCluster();
0202     if (clusterM.amplitudes().size() > MaxNrStrips) {
0203       return;
0204     }
0205     if (deDxTools::isSpanningOver2APV(clusterM.firstStrip(), clusterM.amplitudes().size())) {
0206       return;
0207     }
0208     if (!deDxTools::isFarFromBorder(trajState, &detUnitM)) {
0209       return;
0210     }
0211     float pathLen = 10.0 * detUnitM.surface().bounds().thickness() / fabs(cosine);
0212     float chargeAbs = deDxTools::getCharge(&clusterM, NSaturating, detUnitM, calibGains, m_off);
0213     float charge = chargeAbs / pathLen;
0214     if (!shapetest || (shapetest && deDxTools::shapeSelection(clusterM))) {
0215       Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
0216     }
0217 
0218     auto& detUnitS = *(gdet->stereoDet());
0219     auto& clusterS = matchedHit->stereoCluster();
0220     if (clusterS.amplitudes().size() > MaxNrStrips) {
0221       return;
0222     }
0223     if (deDxTools::isSpanningOver2APV(clusterS.firstStrip(), clusterS.amplitudes().size())) {
0224       return;
0225     }
0226     if (!deDxTools::isFarFromBorder(trajState, &detUnitS)) {
0227       return;
0228     }
0229     pathLen = 10.0 * detUnitS.surface().bounds().thickness() / fabs(cosine);
0230     chargeAbs = deDxTools::getCharge(&clusterS, NSaturating, detUnitS, calibGains, m_off);
0231     charge = chargeAbs / pathLen;
0232     if (!shapetest || (shapetest && deDxTools::shapeSelection(clusterS))) {
0233       Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
0234     }
0235   }
0236 }
0237 
0238 //this function is only used when we run over a calibTree instead of running over EDM files
0239 void DeDxDiscriminatorLearner::algoAnalyzeTheTree(const edm::EventSetup& iSetup) {
0240   const auto& tkGeom = iSetup.getData(m_tkGeomToken);
0241 
0242   unsigned int NEvent = 0;
0243   for (unsigned int i = 0; i < VInputFiles.size(); i++) {
0244     printf("Openning file %3i/%3i --> %s\n", i + 1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str()));
0245     fflush(stdout);
0246     TChain* tree = new TChain("gainCalibrationTree/tree");
0247     tree->Add(VInputFiles[i].c_str());
0248 
0249     TString EventPrefix("");
0250     TString EventSuffix("");
0251 
0252     TString TrackPrefix("track");
0253     TString TrackSuffix("");
0254 
0255     TString CalibPrefix("GainCalibration");
0256     TString CalibSuffix("");
0257 
0258     unsigned int eventnumber = 0;
0259     tree->SetBranchAddress(EventPrefix + "event" + EventSuffix, &eventnumber, nullptr);
0260     unsigned int runnumber = 0;
0261     tree->SetBranchAddress(EventPrefix + "run" + EventSuffix, &runnumber, nullptr);
0262     std::vector<bool>* TrigTech = nullptr;
0263     tree->SetBranchAddress(EventPrefix + "TrigTech" + EventSuffix, &TrigTech, nullptr);
0264 
0265     std::vector<double>* trackchi2ndof = nullptr;
0266     tree->SetBranchAddress(TrackPrefix + "chi2ndof" + TrackSuffix, &trackchi2ndof, nullptr);
0267     std::vector<float>* trackp = nullptr;
0268     tree->SetBranchAddress(TrackPrefix + "momentum" + TrackSuffix, &trackp, nullptr);
0269     std::vector<float>* trackpt = nullptr;
0270     tree->SetBranchAddress(TrackPrefix + "pt" + TrackSuffix, &trackpt, nullptr);
0271     std::vector<double>* tracketa = nullptr;
0272     tree->SetBranchAddress(TrackPrefix + "eta" + TrackSuffix, &tracketa, nullptr);
0273     std::vector<double>* trackphi = nullptr;
0274     tree->SetBranchAddress(TrackPrefix + "phi" + TrackSuffix, &trackphi, nullptr);
0275     std::vector<unsigned int>* trackhitsvalid = nullptr;
0276     tree->SetBranchAddress(TrackPrefix + "hitsvalid" + TrackSuffix, &trackhitsvalid, nullptr);
0277 
0278     std::vector<int>* trackindex = nullptr;
0279     tree->SetBranchAddress(CalibPrefix + "trackindex" + CalibSuffix, &trackindex, nullptr);
0280     std::vector<unsigned int>* rawid = nullptr;
0281     tree->SetBranchAddress(CalibPrefix + "rawid" + CalibSuffix, &rawid, nullptr);
0282     std::vector<unsigned short>* firststrip = nullptr;
0283     tree->SetBranchAddress(CalibPrefix + "firststrip" + CalibSuffix, &firststrip, nullptr);
0284     std::vector<unsigned short>* nstrips = nullptr;
0285     tree->SetBranchAddress(CalibPrefix + "nstrips" + CalibSuffix, &nstrips, nullptr);
0286     std::vector<unsigned int>* charge = nullptr;
0287     tree->SetBranchAddress(CalibPrefix + "charge" + CalibSuffix, &charge, nullptr);
0288     std::vector<float>* path = nullptr;
0289     tree->SetBranchAddress(CalibPrefix + "path" + CalibSuffix, &path, nullptr);
0290     std::vector<unsigned char>* amplitude = nullptr;
0291     tree->SetBranchAddress(CalibPrefix + "amplitude" + CalibSuffix, &amplitude, nullptr);
0292     std::vector<double>* gainused = nullptr;
0293     tree->SetBranchAddress(CalibPrefix + "gainused" + CalibSuffix, &gainused, nullptr);
0294 
0295     printf("Number of Events = %i + %i = %i\n",
0296            NEvent,
0297            (unsigned int)tree->GetEntries(),
0298            (unsigned int)(NEvent + tree->GetEntries()));
0299     NEvent += tree->GetEntries();
0300     printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
0301     printf("Looping on the Tree          :");
0302     int TreeStep = tree->GetEntries() / 50;
0303     if (TreeStep <= 1)
0304       TreeStep = 1;
0305     for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
0306       if (ientry % TreeStep == 0) {
0307         printf(".");
0308         fflush(stdout);
0309       }
0310       tree->GetEntry(ientry);
0311 
0312       int FirstAmplitude = 0;
0313       for (unsigned int c = 0; c < (*path).size(); c++) {
0314         FirstAmplitude += (*nstrips)[c];
0315         int t = (*trackindex)[c];
0316         if ((*trackpt)[t] < 5)
0317           continue;
0318         if ((*trackhitsvalid)[t] < 5)
0319           continue;
0320 
0321         int Charge = 0;
0322         if (useCalibration) {
0323           auto& gains = calibGains[tkGeom.idToDetUnit(DetId((*rawid)[c]))->index() - m_off];
0324           auto& gain = gains[(*firststrip)[c] / 128];
0325           for (unsigned int s = 0; s < (*nstrips)[c]; s++) {
0326             int StripCharge = (*amplitude)[FirstAmplitude - (*nstrips)[c] + s];
0327             if (StripCharge < 254) {
0328               StripCharge = (int)(StripCharge / gain);
0329               if (StripCharge >= 1024) {
0330                 StripCharge = 255;
0331               } else if (StripCharge >= 254) {
0332                 StripCharge = 254;
0333               }
0334             }
0335             Charge += StripCharge;
0336           }
0337         } else {
0338           Charge = (*charge)[c];
0339         }
0340 
0341         //          printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[c],Charge,Gains[(*rawid)[c]]);
0342         double ClusterChargeOverPath = ((double)Charge) / (*path)[c];
0343         Charge_Vs_Path->Fill((*trackp)[t], (*path)[c], ClusterChargeOverPath);
0344       }
0345     }
0346     printf("\n");
0347   }
0348 }
0349 
0350 std::unique_ptr<PhysicsTools::Calibration::HistogramD3D> DeDxDiscriminatorLearner::getNewObject() {
0351   auto obj = std::make_unique<PhysicsTools::Calibration::HistogramD3D>(Charge_Vs_Path->GetNbinsX(),
0352                                                                        Charge_Vs_Path->GetXaxis()->GetXmin(),
0353                                                                        Charge_Vs_Path->GetXaxis()->GetXmax(),
0354                                                                        Charge_Vs_Path->GetNbinsY(),
0355                                                                        Charge_Vs_Path->GetYaxis()->GetXmin(),
0356                                                                        Charge_Vs_Path->GetYaxis()->GetXmax(),
0357                                                                        Charge_Vs_Path->GetNbinsZ(),
0358                                                                        Charge_Vs_Path->GetZaxis()->GetXmin(),
0359                                                                        Charge_Vs_Path->GetZaxis()->GetXmax());
0360 
0361   for (int ix = 0; ix <= Charge_Vs_Path->GetNbinsX() + 1; ix++) {
0362     for (int iy = 0; iy <= Charge_Vs_Path->GetNbinsY() + 1; iy++) {
0363       for (int iz = 0; iz <= Charge_Vs_Path->GetNbinsZ() + 1; iz++) {
0364         obj->setBinContent(ix, iy, iz, Charge_Vs_Path->GetBinContent(ix, iy, iz));
0365         //          if(Charge_Vs_Path->GetBinContent(ix,iy)!=0)printf("%i %i %i --> %f\n",ix,iy, iz, Charge_Vs_Path->GetBinContent(ix,iy,iz));
0366       }
0367     }
0368   }
0369 
0370   return obj;
0371 }
0372 
0373 //define this as a plug-in
0374 DEFINE_FWK_MODULE(DeDxDiscriminatorLearner);
[ 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!