CommonAlignmentMonitor/plugins/AlignmentStats.cc

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