Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-03 01:57:38

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