Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-05-10 03:53:01

0001 #include "Alignment/CommonAlignmentMonitor/plugins/AlignmentStats.h"
0002 
0003 #include "DataFormats/Common/interface/View.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0007 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0008 #include "Alignment/CommonAlignment/interface/Alignable.h"
0009 #include "Alignment/CommonAlignment/interface/Utilities.h"
0010 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0011 
0012 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0015 #include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h"
0016 #include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h"
0017 
0018 #include "DataFormats/Alignment/interface/AlignmentClusterFlag.h"
0019 
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0023 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0024 #include "Utilities/General/interface/ClassName.h"
0025 
0026 using namespace std;
0027 
0028 AlignmentStats::AlignmentStats(const edm::ParameterSet &iConfig)
0029     : esTokenTTopo_(esConsumes()),
0030       esTokenTkGeo_(esConsumes()),
0031       src_(iConfig.getParameter<edm::InputTag>("src")),
0032       overlapAM_(iConfig.getParameter<edm::InputTag>("OverlapAssoMap")),
0033       keepTrackStats_(iConfig.getParameter<bool>("keepTrackStats")),
0034       keepHitPopulation_(iConfig.getParameter<bool>("keepHitStats")),
0035       statsTreeName_(iConfig.getParameter<string>("TrkStatsFileName")),
0036       hitsTreeName_(iConfig.getParameter<string>("HitStatsFileName")),
0037       prescale_(iConfig.getParameter<uint32_t>("TrkStatsPrescale")),
0038       trackToken_(consumes<reco::TrackCollection>(src_)),
0039       mapToken_(consumes<AliClusterValueMap>(overlapAM_)) {
0040   //sanity checks
0041 
0042   //init
0043   outtree_ = nullptr;
0044 
0045 }  //end constructor
0046 
0047 void AlignmentStats::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0048   edm::ParameterSetDescription desc;
0049   desc.add<edm::InputTag>("src", edm::InputTag("generalTracks"));
0050   desc.add<edm::InputTag>("OverlapAssoMap", edm::InputTag("OverlapAssoMap"));
0051   desc.add<bool>("keepTrackStats", false);
0052   desc.add<bool>("keepHitStats", false);
0053   desc.add<std::string>("TrkStatsFileName", "tracks_statistics.root");
0054   desc.add<std::string>("HitStatsFileName", "HitMaps.root");
0055   desc.add<unsigned int>("TrkStatsPrescale", 1);
0056   descriptions.add("AlignmentStats", desc);
0057 }
0058 
0059 void AlignmentStats::beginJob() {  // const edm::EventSetup &iSetup
0060 
0061   //book track stats tree
0062   treefile_ = new TFile(statsTreeName_.c_str(), "RECREATE");
0063   treefile_->cd();
0064   outtree_ = new TTree("AlignmentTrackStats", "Statistics of Tracks used for Alignment");
0065   // int nHitsinPXB[MAXTRKS_], nHitsinPXE[MAXTRKS_], nHitsinTEC[MAXTRKS_], nHitsinTIB[MAXTRKS_],nHitsinTOB[MAXTRKS_],nHitsinTID[MAXTRKS_];
0066 
0067   outtree_->Branch("Ntracks", &ntracks, "Ntracks/i");
0068   outtree_->Branch("Run_", &run_, "Run_Nr/I");
0069   outtree_->Branch("Event", &event_, "EventNr/I");
0070   outtree_->Branch("Eta", &Eta, "Eta[Ntracks]/F");
0071   outtree_->Branch("Phi", &Phi, "Phi[Ntracks]/F");
0072   outtree_->Branch("P", &P, "P[Ntracks]/F");
0073   outtree_->Branch("Pt", &Pt, "Pt[Ntracks]/F");
0074   outtree_->Branch("Chi2n", &Chi2n, "Chi2n[Ntracks]/F");
0075   outtree_->Branch("Nhits", &Nhits, "Nhits[Ntracks][7]/I");
0076   /*
0077     outtree_->Branch("NhitsPXB"       , ,);
0078     outtree_->Branch("NhitsPXE"       , ,);
0079     outtree_->Branch("NhitsTIB"       , ,);
0080     outtree_->Branch("NhitsTID"       , ,);
0081     outtree_->Branch("NhitsTOB"       , ,);
0082     outtree_->Branch("NhitsTOB"       , ,);
0083   */
0084 
0085   tmpPresc_ = prescale_;
0086 }  //end beginJob
0087 
0088 void AlignmentStats::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0089   //load list of detunits needed then in endJob
0090   if (!trackerGeometry_) {
0091     trackerGeometry_ = std::make_unique<TrackerGeometry>(iSetup.getData(esTokenTkGeo_));
0092   }
0093 
0094   if (!trackerTopology_) {
0095     trackerTopology_ = std::make_unique<TrackerTopology>(iSetup.getData(esTokenTTopo_));
0096   }
0097 
0098   //take trajectories and tracks to loop on
0099   // edm::Handle<TrajTrackAssociationCollection> TrackAssoMap;
0100   const edm::Handle<reco::TrackCollection> &Tracks = iEvent.getHandle(trackToken_);
0101 
0102   //take overlap HitAssomap
0103   const edm::Handle<AliClusterValueMap> &hMap = iEvent.getHandle(mapToken_);
0104   const AliClusterValueMap &OverlapMap = *hMap;
0105 
0106   // Initialise
0107   run_ = 1;
0108   event_ = 1;
0109   ntracks = 0;
0110   run_ = iEvent.id().run();
0111   event_ = iEvent.id().event();
0112   ntracks = Tracks->size();
0113   if (ntracks > 1)
0114     edm::LogVerbatim("AlignmenStats") << "~~~~~~~~~~~~\n For this event processing " << ntracks << " tracks";
0115 
0116   unsigned int trk_cnt = 0;
0117 
0118   for (int j = 0; j < MAXTRKS_; j++) {
0119     Eta[j] = -9999.0;
0120     Phi[j] = -8888.0;
0121     P[j] = -7777.0;
0122     Pt[j] = -6666.0;
0123     Chi2n[j] = -2222.0;
0124     for (int k = 0; k < 7; k++) {
0125       Nhits[j][k] = 0;
0126     }
0127   }
0128 
0129   // int npxbhits=0;
0130 
0131   //loop on tracks
0132   for (std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk;
0133        ++ittrk) {
0134     Eta[trk_cnt] = ittrk->eta();
0135     Phi[trk_cnt] = ittrk->phi();
0136     Chi2n[trk_cnt] = ittrk->normalizedChi2();
0137     P[trk_cnt] = ittrk->p();
0138     Pt[trk_cnt] = ittrk->pt();
0139     Nhits[trk_cnt][0] = ittrk->numberOfValidHits();
0140 
0141     if (ntracks > 1)
0142       edm::LogVerbatim("AlignmenStats") << "Track #" << trk_cnt + 1 << " params:    Eta=" << Eta[trk_cnt]
0143                                         << "  Phi=" << Phi[trk_cnt] << "  P=" << P[trk_cnt]
0144                                         << "   Nhits=" << Nhits[trk_cnt][0];
0145 
0146     //loop on tracking rechits
0147     //edm::LogVerbatim("AlignmenStats") << "   loop on hits of track #" << (itt - tracks->begin());
0148     for (auto const &hit : ittrk->recHits()) {
0149       if (!hit->isValid())
0150         continue;
0151       DetId detid = hit->geographicalId();
0152       int subDet = detid.subdetId();
0153       uint32_t rawId = hit->geographicalId().rawId();
0154 
0155       //  if(subDet==1)npxbhits++;
0156 
0157       //look if you find this detid in the map
0158       DetHitMap::iterator mapiter;
0159       mapiter = hitmap_.find(rawId);
0160       if (mapiter != hitmap_.end()) {  //present, increase its value by one
0161         //  hitmap_[rawId]=hitmap_[rawId]+1;
0162         ++(hitmap_[rawId]);
0163       } else {  //not present, let's add this key to the map with value=1
0164         hitmap_.insert(pair<uint32_t, uint32_t>(rawId, 1));
0165       }
0166 
0167       AlignmentClusterFlag inval;
0168 
0169       bool hitInPixel = (subDet == PixelSubdetector::PixelBarrel) || (subDet == PixelSubdetector::PixelEndcap);
0170       bool hitInStrip = (subDet == SiStripDetId::TIB) || (subDet == SiStripDetId::TID) ||
0171                         (subDet == SiStripDetId::TOB) || (subDet == SiStripDetId::TEC);
0172 
0173       if (!(hitInPixel || hitInStrip)) {
0174         //skip only this hit, don't stop everything throwing an exception
0175         edm::LogError("AlignmentStats") << "Hit not belonging neither to pixels nor to strips! Skipping it. SubDet= "
0176                                         << subDet;
0177         continue;
0178       }
0179 
0180       //check also if the hit is an overlap. If yes fill a dedicated hitmap
0181       if (hitInStrip) {
0182         const std::type_info &type = typeid(*hit);
0183 
0184         if (type == typeid(SiStripRecHit1D)) {
0185           //Notice the difference respect to when one loops on Trajectories: the recHit is a TrackingRecHit and not a TransientTrackingRecHit
0186           const SiStripRecHit1D *striphit = dynamic_cast<const SiStripRecHit1D *>(hit);
0187           if (striphit != nullptr) {
0188             SiStripRecHit1D::ClusterRef stripclust(striphit->cluster());
0189             inval = OverlapMap[stripclust];
0190           } else {
0191             //  edm::LogError("AlignmentStats")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit1D failed!   TypeId of the RecHit: "<<className(*hit);
0192             throw cms::Exception("NullPointerError")
0193                 << "ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit1D failed!   TypeId of the RecHit: "
0194                 << className(*hit);
0195           }
0196         }  //end if sistriprechit1D
0197         if (type == typeid(SiStripRecHit2D)) {
0198           const SiStripRecHit2D *striphit = dynamic_cast<const SiStripRecHit2D *>(hit);
0199           if (striphit != nullptr) {
0200             SiStripRecHit2D::ClusterRef stripclust(striphit->cluster());
0201             inval = OverlapMap[stripclust];
0202             //edm::LogVerbatim("AlignmenStats")<<"Taken the Strip Cluster with ProdId "<<stripclust.id() <<"; the Value in the map is "<<inval<<"  (DetId is "<<hit->geographicalId().rawId()<<")";
0203           } else {
0204             throw cms::Exception("NullPointerError")
0205                 << "ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit2D failed!   TypeId of the RecHit: "
0206                 << className(*hit);
0207             //    edm::LogError("AlignmentStats")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit2D failed!   TypeId of the RecHit: "<<className(*hit);
0208           }
0209         }  //end if sistriprechit2D
0210 
0211       }  //end if hit in Strips
0212       else {
0213         const SiPixelRecHit *pixelhit = dynamic_cast<const SiPixelRecHit *>(hit);
0214         if (pixelhit != nullptr) {
0215           SiPixelClusterRefNew pixclust(pixelhit->cluster());
0216           inval = OverlapMap[pixclust];
0217           //edm::LogVerbatim("AlignmenStats")<<"Taken the Pixel Cluster with ProdId "<<pixclust.id() <<"; the Value in the map is "<<inval<<"  (DetId is "<<hit->geographicalId().rawId()<<")";
0218         } else {
0219           edm::LogError("AlignmentStats")
0220               << "ERROR in <AlignmentStats::analyze>: Dynamic cast of Pixel RecHit failed!   TypeId of the RecHit: "
0221               << className(*hit);
0222         }
0223       }  //end else hit is in Pixel
0224 
0225       bool isOverlapHit(inval.isOverlap());
0226 
0227       if (isOverlapHit) {
0228         edm::LogVerbatim("AlignmenStats") << "This hit is an overlap !";
0229         DetHitMap::iterator overlapiter;
0230         overlapiter = overlapmap_.find(rawId);
0231 
0232         if (overlapiter !=
0233             overlapmap_.end()) {  //the det already collected at least an overlap, increase its value by one
0234           overlapmap_[rawId] = overlapmap_[rawId] + 1;
0235         } else {  //first overlap on det unit, let's add it to the map
0236           overlapmap_.insert(pair<uint32_t, uint32_t>(rawId, 1));
0237         }
0238       }  //end if the hit is an overlap
0239 
0240       int subdethit = static_cast<int>(hit->geographicalId().subdetId());
0241       if (ntracks > 1)
0242         edm::LogVerbatim("AlignmenStats") << "Hit in SubDet=" << subdethit;
0243       Nhits[trk_cnt][subdethit] = Nhits[trk_cnt][subdethit] + 1;
0244     }  //end loop on trackingrechits
0245     trk_cnt++;
0246 
0247   }  //end loop on tracks
0248 
0249   //edm::LogVerbatim("AlignmenStats")<<"Total number of pixel hits is " << npxbhits;
0250 
0251   tmpPresc_--;
0252   if (tmpPresc_ < 1) {
0253     outtree_->Fill();
0254     tmpPresc_ = prescale_;
0255   }
0256   if (trk_cnt != ntracks)
0257     edm::LogError("AlignmentStats") << "\nERROR! trk_cnt=" << trk_cnt << "   ntracks=" << ntracks;
0258 
0259   return;
0260 }
0261 
0262 void AlignmentStats::endJob() {
0263   treefile_->cd();
0264   edm::LogInfo("AlignmentStats") << "Writing out the TrackStatistics in " << gDirectory->GetPath();
0265   outtree_->Write();
0266   delete outtree_;
0267 
0268   //create tree with hit maps (hitstree)
0269   //book track stats tree
0270   TFile *hitsfile = new TFile(hitsTreeName_.c_str(), "RECREATE");
0271   hitsfile->cd();
0272   TTree *hitstree = new TTree("AlignmentHitMap", "Maps of Hits used for Alignment");
0273 
0274   unsigned int id = 0, nhits = 0, noverlaps = 0;
0275   float posX(-99999.0), posY(-77777.0), posZ(-88888.0);
0276   float posEta(-6666.0), posPhi(-5555.0), posR(-4444.0);
0277   int subdet = 0;
0278   unsigned int layer = 0;
0279   bool is2D = false, isStereo = false;
0280   hitstree->Branch("DetId", &id, "DetId/i");
0281   hitstree->Branch("Nhits", &nhits, "Nhits/i");
0282   hitstree->Branch("Noverlaps", &noverlaps, "Noverlaps/i");
0283   hitstree->Branch("SubDet", &subdet, "SubDet/I");
0284   hitstree->Branch("Layer", &layer, "Layer/i");
0285   hitstree->Branch("is2D", &is2D, "is2D/B");
0286   hitstree->Branch("isStereo", &isStereo, "isStereo/B");
0287   hitstree->Branch("posX", &posX, "posX/F");
0288   hitstree->Branch("posY", &posY, "posY/F");
0289   hitstree->Branch("posZ", &posZ, "posZ/F");
0290   hitstree->Branch("posR", &posR, "posR/F");
0291   hitstree->Branch("posEta", &posEta, "posEta/F");
0292   hitstree->Branch("posPhi", &posPhi, "posPhi/F");
0293 
0294   /*
0295   TTree *overlapstree=new TTree("OverlapHitMap","Maps of Overlaps used for Alignment");
0296   hitstree->Branch("DetId",   &id ,     "DetId/i");
0297   hitstree->Branch("NOverlaps",   &nhits ,  "Nhits/i");
0298   hitstree->Branch("SubDet",  &subdet,  "SubDet/I");
0299   hitstree->Branch("Layer",   &layer,   "Layer/i");
0300   hitstree->Branch("is2D" ,   &is2D,    "is2D/B");
0301   hitstree->Branch("isStereo",&isStereo,"isStereo/B");
0302   hitstree->Branch("posX",    &posX,    "posX/F");
0303   hitstree->Branch("posY",    &posY,    "posY/F");
0304   hitstree->Branch("posZ",    &posZ,    "posZ/F");
0305   hitstree->Branch("posR",    &posR,    "posR/F");
0306   hitstree->Branch("posEta",  &posEta,  "posEta/F");
0307   hitstree->Branch("posPhi",  &posPhi,  "posPhi/F");
0308   */
0309 
0310   std::unique_ptr<AlignableTracker> theAliTracker =
0311       std::make_unique<AlignableTracker>(trackerGeometry_.get(), trackerTopology_.get());
0312   const auto &Detunitslist = theAliTracker->deepComponents();
0313   int ndetunits = Detunitslist.size();
0314   edm::LogInfo("AlignmentStats") << "Number of DetUnits in the AlignableTracker: " << ndetunits;
0315 
0316   for (int det_cnt = 0; det_cnt < ndetunits; ++det_cnt) {
0317     //re-initialize for safety
0318     id = 0;
0319     nhits = 0;
0320     noverlaps = 0;
0321     posX = -99999.0;
0322     posY = -77777.0;
0323     posZ = -88888.0;
0324     posEta = -6666.0;
0325     posPhi = -5555.0;
0326     posR = -4444.0;
0327     subdet = 0;
0328     layer = 0;
0329     is2D = false;
0330     isStereo = false;
0331 
0332     //if detunit in vector is found also in the map, look for how many hits were collected
0333     //and save in the tree this number
0334     id = static_cast<uint32_t>(Detunitslist[det_cnt]->id());
0335     if (hitmap_.find(id) != hitmap_.end()) {
0336       nhits = hitmap_[id];
0337     }
0338     //if not, save nhits=0
0339     else {
0340       nhits = 0;
0341       hitmap_.insert(pair<uint32_t, uint32_t>(id, 0));
0342     }
0343 
0344     if (overlapmap_.find(id) != overlapmap_.end()) {
0345       noverlaps = overlapmap_[id];
0346     }
0347     //if not, save nhits=0
0348     else {
0349       noverlaps = 0;
0350       overlapmap_.insert(pair<uint32_t, uint32_t>(id, 0));
0351     }
0352 
0353     //take other geometrical infos from the det
0354     posX = Detunitslist[det_cnt]->globalPosition().x();
0355     posY = Detunitslist[det_cnt]->globalPosition().y();
0356     posZ = Detunitslist[det_cnt]->globalPosition().z();
0357 
0358     align::GlobalVector vec(posX, posY, posZ);
0359     posR = vec.perp();
0360     posPhi = vec.phi();
0361     posEta = vec.eta();
0362     //   posPhi = atan2(posY,posX);
0363 
0364     DetId detid(id);
0365     subdet = detid.subdetId();
0366 
0367     //get layers, petals, etc...
0368     if (subdet == PixelSubdetector::PixelBarrel) {  //PXB
0369 
0370       layer = trackerTopology_->pxbLayer(id);
0371       is2D = true;
0372       isStereo = false;
0373     } else if (subdet == PixelSubdetector::PixelEndcap) {
0374       layer = trackerTopology_->pxfDisk(id);
0375       is2D = true;
0376       isStereo = false;
0377     } else if (subdet == SiStripDetId::TIB) {
0378       layer = trackerTopology_->tibLayer(id);
0379       is2D = trackerTopology_->tibIsDoubleSide(id);
0380       isStereo = trackerTopology_->tibIsStereo(id);
0381     } else if (subdet == SiStripDetId::TID) {
0382       layer = trackerTopology_->tidWheel(id);
0383       is2D = trackerTopology_->tidIsDoubleSide(id);
0384       isStereo = trackerTopology_->tidIsStereo(id);
0385     } else if (subdet == SiStripDetId::TOB) {
0386       layer = trackerTopology_->tobLayer(id);
0387       is2D = trackerTopology_->tobIsDoubleSide(id);
0388       isStereo = trackerTopology_->tobIsStereo(id);
0389     } else if (subdet == SiStripDetId::TEC) {
0390       layer = trackerTopology_->tecWheel(id);
0391       is2D = trackerTopology_->tecIsDoubleSide(id);
0392       isStereo = trackerTopology_->tecIsStereo(id);
0393     } else {
0394       edm::LogError("AlignmentStats") << "Detector not belonging neither to pixels nor to strips! Skipping it. SubDet= "
0395                                       << subdet;
0396     }
0397 
0398     //write in the hitstree
0399     hitstree->Fill();
0400   }  //end loop over detunits
0401 
0402   //save hitstree
0403   hitstree->Write();
0404   delete hitstree;
0405   //delete Detunitslist;
0406   hitmap_.clear();
0407   overlapmap_.clear();
0408   delete hitsfile;
0409 }
0410 // ========= MODULE DEF ==============
0411 #include "FWCore/PluginManager/interface/ModuleDef.h"
0412 #include "FWCore/Framework/interface/MakerMacros.h"
0413 DEFINE_FWK_MODULE(AlignmentStats);