Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:05:17

0001 /*
0002    class Phase2PixelNtuple
0003    */
0004 
0005 // DataFormats
0006 #include "DataFormats/Common/interface/DetSetVector.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0010 
0011 #include "DataFormats/DetId/interface/DetId.h"
0012 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0013 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0014 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0015 
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 
0018 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0019 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0020 
0021 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0022 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0023 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0024 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0025 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0026 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0027 
0028 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0029 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0031 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0032 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0033 
0034 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0035 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0036 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0037 #include "SimDataFormats/Track/interface/SimTrack.h"
0038 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0039 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0040 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0041 
0042 #include "FWCore/Framework/interface/MakerMacros.h"
0043 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0044 #include "FWCore/Framework/interface/Event.h"
0045 #include "FWCore/Framework/interface/EventSetup.h"
0046 #include "FWCore/Utilities/interface/InputTag.h"
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0049 
0050 // Geometry
0051 #include "MagneticField/Engine/interface/MagneticField.h"
0052 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0053 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0054 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0055 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0056 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0057 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0058 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0059 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0060 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0061 
0062 // For ROOT
0063 #include "FWCore/ServiceRegistry/interface/Service.h"
0064 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0065 #include <TTree.h>
0066 
0067 // STD
0068 #include <memory>
0069 #include <string>
0070 #include <iostream>
0071 
0072 // CLHEP (for speed of light)
0073 #include "CLHEP/Units/PhysicalConstants.h"
0074 #include "CLHEP/Units/SystemOfUnits.h"
0075 
0076 //#define EDM_ML_DEBUG
0077 
0078 using namespace std;
0079 using namespace edm;
0080 using namespace reco;
0081 
0082 class Phase2PixelNtuple : public edm::one::EDAnalyzer<> {
0083 public:
0084   explicit Phase2PixelNtuple(const edm::ParameterSet& conf);
0085   virtual ~Phase2PixelNtuple();
0086   virtual void beginJob();
0087   virtual void endJob();
0088   virtual void analyze(const edm::Event& e, const edm::EventSetup& es);
0089 
0090 protected:
0091   void fillEvt(const edm::Event&);
0092   void fillPRecHit(const int detid_db,
0093                    const int subid,
0094                    const int layer_num,
0095                    const int ladder_num,
0096                    const int module_num,
0097                    const int disk_num,
0098                    const int blade_num,
0099                    const int panel_num,
0100                    const int side_num,
0101                    SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0102                    const int num_simhit,
0103                    const PSimHit* closest_simhit,
0104                    const GeomDet* PixGeom);
0105   void fillPRecHit(const int detid_db,
0106                    const int subid,
0107                    const int layer_num,
0108                    const int ladder_num,
0109                    const int module_num,
0110                    const int disk_num,
0111                    const int blade_num,
0112                    const int panel_num,
0113                    const int side_num,
0114                    const TrackingRecHit* transRecHit,
0115                    const SiPixelRecHit* pixHit,
0116                    const int num_simhit,
0117                    const PSimHit* closest_simhit,
0118                    const GeomDet* PixGeom,
0119                    const TrajectoryStateOnSurface tsos);
0120 
0121   std::pair<float, float> computeAnglesFromDetPosition(const SiPixelCluster& cl,
0122                                                        const PixelTopology& top,
0123                                                        const GeomDetUnit& det) const;
0124 
0125 private:
0126   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geom_esToken;
0127   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topo_esToken;
0128   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> ttkb_esToken;
0129   edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> tthb_esToken;
0130 
0131   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0132   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelRecHit>> pixelRecHits_token_;
0133   edm::EDGetTokenT<edm::View<reco::Track>> recoTracks_token_;
0134   edm::EDGetTokenT<TrajTrackAssociationCollection> tta_token_;
0135 
0136   bool verbose_;
0137   bool picky_;
0138   static const int DGPERCLMAX = 100;
0139 
0140   float trkPt_, trkEta_, trkTheta_, trkPhi_;
0141   int trkIsHighPurity_;
0142 
0143   //--- Structures for ntupling:
0144   struct evt {
0145     int run;
0146     int evtnum;
0147 
0148     void init();
0149   } evt_;
0150 
0151   struct RecHit {
0152     int pdgid;
0153     int process;
0154     float q;
0155     float x;
0156     float y;
0157     float xx;
0158     float xy;
0159     float yy;
0160 
0161     float probQ;
0162     float probXY;
0163 
0164     float row;
0165     float col;
0166     float gx;
0167     float gy;
0168     float gz;
0169     int subid, module;
0170     int layer, ladder;             // BPix
0171     int disk, blade, panel, side;  // FPix
0172     int nsimhit;
0173     int spreadx, spready;
0174     float hx, hy, ht;
0175     float hrow, hcol;
0176     float tx, ty, tz;
0177     float theta, phi;
0178 
0179     // detector topology
0180     int nRowsInDet;
0181     int nColsInDet;
0182     float pitchx;
0183     float pitchy;
0184     float thickness;
0185 
0186     // local angles from det position
0187     float cotAlphaFromDet, cotBetaFromDet;
0188     // local angles from track trajectory
0189     float cotAlphaFromTrack, cotBetaFromTrack;
0190 
0191     // digis
0192     int fDgN;
0193     int fDgRow[DGPERCLMAX], fDgCol[DGPERCLMAX];
0194     int fDgDetId[DGPERCLMAX];
0195     float fDgAdc[DGPERCLMAX], fDgCharge[DGPERCLMAX];
0196 
0197     void init();
0198   } recHit_;
0199 
0200   edm::Service<TFileService> tFileService;
0201   TTree* pixeltree_;
0202   TTree* pixeltreeOnTrack_;
0203 };
0204 
0205 Phase2PixelNtuple::Phase2PixelNtuple(edm::ParameterSet const& conf)
0206     : geom_esToken(esConsumes()),
0207       topo_esToken(esConsumes()),
0208       ttkb_esToken(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0209       tthb_esToken(esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("ttrhBuilder")))),
0210       trackerHitAssociatorConfig_(conf, consumesCollector()),
0211       pixelRecHits_token_(consumes<edmNew::DetSetVector<SiPixelRecHit>>(edm::InputTag("siPixelRecHits"))),
0212       recoTracks_token_(consumes<edm::View<reco::Track>>(conf.getParameter<edm::InputTag>("trackProducer"))),
0213       tta_token_(consumes<TrajTrackAssociationCollection>(conf.getParameter<InputTag>("trajectoryInput"))),
0214       verbose_(conf.getUntrackedParameter<bool>("verbose", false)),
0215       picky_(conf.getUntrackedParameter<bool>("picky", false)),
0216       pixeltree_(0),
0217       pixeltreeOnTrack_(0) {}
0218 
0219 Phase2PixelNtuple::~Phase2PixelNtuple() {}
0220 
0221 void Phase2PixelNtuple::endJob() {}
0222 
0223 void Phase2PixelNtuple::beginJob() {
0224   pixeltree_ = tFileService->make<TTree>("PixelNtuple", "Pixel hit analyzer ntuple");
0225   pixeltreeOnTrack_ = tFileService->make<TTree>("PixelNtupleOnTrack", "On-Track Pixel hit analyzer ntuple");
0226 
0227   int bufsize = 64000;
0228   //Common Branch
0229   pixeltree_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize);
0230   pixeltree_->Branch("pdgid", &recHit_.pdgid, "pdgid/I");
0231   pixeltree_->Branch("process", &recHit_.process, "process/I");
0232   pixeltree_->Branch("q", &recHit_.q, "q/F");
0233   pixeltree_->Branch("x", &recHit_.x, "x/F");
0234   pixeltree_->Branch("y", &recHit_.y, "y/F");
0235   pixeltree_->Branch("xx", &recHit_.xx, "xx/F");
0236   pixeltree_->Branch("xy", &recHit_.xy, "xy/F");
0237   pixeltree_->Branch("yy", &recHit_.yy, "yy/F");
0238   pixeltree_->Branch("row", &recHit_.row, "row/F");
0239   pixeltree_->Branch("col", &recHit_.col, "col/F");
0240   pixeltree_->Branch("gx", &recHit_.gx, "gx/F");
0241   pixeltree_->Branch("gy", &recHit_.gy, "gy/F");
0242   pixeltree_->Branch("gz", &recHit_.gz, "gz/F");
0243   pixeltree_->Branch("subid", &recHit_.subid, "subid/I");
0244   pixeltree_->Branch("module", &recHit_.module, "module/I");
0245   pixeltree_->Branch("layer", &recHit_.layer, "layer/I");
0246   pixeltree_->Branch("ladder", &recHit_.ladder, "ladder/I");
0247   pixeltree_->Branch("disk", &recHit_.disk, "disk/I");
0248   pixeltree_->Branch("blade", &recHit_.blade, "blade/I");
0249   pixeltree_->Branch("panel", &recHit_.panel, "panel/I");
0250   pixeltree_->Branch("side", &recHit_.side, "side/I");
0251   pixeltree_->Branch("nsimhit", &recHit_.nsimhit, "nsimhit/I");
0252   pixeltree_->Branch("spreadx", &recHit_.spreadx, "spreadx/I");
0253   pixeltree_->Branch("spready", &recHit_.spready, "spready/I");
0254   pixeltree_->Branch("hx", &recHit_.hx, "hx/F");
0255   pixeltree_->Branch("hy", &recHit_.hy, "hy/F");
0256   pixeltree_->Branch("ht", &recHit_.ht, "ht/F");
0257   pixeltree_->Branch("hrow", &recHit_.hrow, "hrow/F");
0258   pixeltree_->Branch("hcol", &recHit_.hcol, "hcol/F");
0259   pixeltree_->Branch("tx", &recHit_.tx, "tx/F");
0260   pixeltree_->Branch("ty", &recHit_.ty, "ty/F");
0261   pixeltree_->Branch("tz", &recHit_.tz, "tz/F");
0262   pixeltree_->Branch("theta", &recHit_.theta, "theta/F");
0263   pixeltree_->Branch("phi", &recHit_.phi, "phi/F");
0264   pixeltree_->Branch("nRowsInDet", &recHit_.nRowsInDet, "nRowsInDet/I");
0265   pixeltree_->Branch("nColsInDet", &recHit_.nColsInDet, "nColsInDet/I");
0266   pixeltree_->Branch("pitchx", &recHit_.pitchx, "pitchx/F");
0267   pixeltree_->Branch("pitchy", &recHit_.pitchy, "pitchy/F");
0268   pixeltree_->Branch("thickness", &recHit_.thickness, "thickness/F");
0269   pixeltree_->Branch("cotAlphaFromDet", &recHit_.cotAlphaFromDet, "cotAlphaFromDet/F");
0270   pixeltree_->Branch("cotBetaFromDet", &recHit_.cotBetaFromDet, "cotBetaFromDet/F");
0271 
0272   pixeltree_->Branch("DgN", &recHit_.fDgN, "DgN/I");
0273   pixeltree_->Branch("DgRow", recHit_.fDgRow, "DgRow[DgN]/I");
0274   pixeltree_->Branch("DgCol", recHit_.fDgCol, "DgCol[DgN]/I");
0275   pixeltree_->Branch("DgDetId", recHit_.fDgDetId, "DgDetId[DgN]/I");
0276   pixeltree_->Branch("DgAdc", recHit_.fDgAdc, "DgAdc[DgN]/F");
0277   pixeltree_->Branch("DgCharge", recHit_.fDgCharge, "DgCharge[DgN]/F");
0278 
0279   //Common Branch (on-track)
0280   pixeltreeOnTrack_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize);
0281   pixeltreeOnTrack_->Branch("pdgid", &recHit_.pdgid, "pdgid/I");
0282   pixeltreeOnTrack_->Branch("process", &recHit_.process, "process/I");
0283   pixeltreeOnTrack_->Branch("q", &recHit_.q, "q/F");
0284   pixeltreeOnTrack_->Branch("x", &recHit_.x, "x/F");
0285   pixeltreeOnTrack_->Branch("y", &recHit_.y, "y/F");
0286   pixeltreeOnTrack_->Branch("xx", &recHit_.xx, "xx/F");
0287   pixeltreeOnTrack_->Branch("xy", &recHit_.xy, "xy/F");
0288   pixeltreeOnTrack_->Branch("yy", &recHit_.yy, "yy/F");
0289   pixeltreeOnTrack_->Branch("probQ", &recHit_.probQ);
0290   pixeltreeOnTrack_->Branch("probXY", &recHit_.probXY);
0291 
0292   pixeltreeOnTrack_->Branch("row", &recHit_.row, "row/F");
0293   pixeltreeOnTrack_->Branch("col", &recHit_.col, "col/F");
0294   pixeltreeOnTrack_->Branch("gx", &recHit_.gx, "gx/F");
0295   pixeltreeOnTrack_->Branch("gy", &recHit_.gy, "gy/F");
0296   pixeltreeOnTrack_->Branch("gz", &recHit_.gz, "gz/F");
0297   pixeltreeOnTrack_->Branch("subid", &recHit_.subid, "subid/I");
0298   pixeltreeOnTrack_->Branch("module", &recHit_.module, "module/I");
0299   pixeltreeOnTrack_->Branch("layer", &recHit_.layer, "layer/I");
0300   pixeltreeOnTrack_->Branch("ladder", &recHit_.ladder, "ladder/I");
0301   pixeltreeOnTrack_->Branch("disk", &recHit_.disk, "disk/I");
0302   pixeltreeOnTrack_->Branch("blade", &recHit_.blade, "blade/I");
0303   pixeltreeOnTrack_->Branch("panel", &recHit_.panel, "panel/I");
0304   pixeltreeOnTrack_->Branch("side", &recHit_.side, "side/I");
0305   pixeltreeOnTrack_->Branch("nsimhit", &recHit_.nsimhit, "nsimhit/I");
0306   pixeltreeOnTrack_->Branch("spreadx", &recHit_.spreadx, "spreadx/I");
0307   pixeltreeOnTrack_->Branch("spready", &recHit_.spready, "spready/I");
0308   pixeltreeOnTrack_->Branch("hx", &recHit_.hx, "hx/F");
0309   pixeltreeOnTrack_->Branch("hy", &recHit_.hy, "hy/F");
0310   pixeltreeOnTrack_->Branch("ht", &recHit_.ht, "ht/F");
0311   pixeltreeOnTrack_->Branch("hrow", &recHit_.hrow, "hrow/F");
0312   pixeltreeOnTrack_->Branch("hcol", &recHit_.hcol, "hcol/F");
0313   pixeltreeOnTrack_->Branch("tx", &recHit_.tx, "tx/F");
0314   pixeltreeOnTrack_->Branch("ty", &recHit_.ty, "ty/F");
0315   pixeltreeOnTrack_->Branch("tz", &recHit_.tz, "tz/F");
0316   pixeltreeOnTrack_->Branch("theta", &recHit_.theta, "theta/F");
0317   pixeltreeOnTrack_->Branch("phi", &recHit_.phi, "phi/F");
0318   pixeltreeOnTrack_->Branch("nRowsInDet", &recHit_.nRowsInDet, "nRowsInDet/I");
0319   pixeltreeOnTrack_->Branch("nColsInDet", &recHit_.nColsInDet, "nColsInDet/I");
0320   pixeltreeOnTrack_->Branch("pitchx", &recHit_.pitchx, "pitchx/F");
0321   pixeltreeOnTrack_->Branch("pitchy", &recHit_.pitchy, "pitchy/F");
0322   pixeltreeOnTrack_->Branch("thickness", &recHit_.thickness, "thickness/F");
0323   pixeltreeOnTrack_->Branch("cotAlphaFromDet", &recHit_.cotAlphaFromDet, "cotAlphaFromDet/F");
0324   pixeltreeOnTrack_->Branch("cotBetaFromDet", &recHit_.cotBetaFromDet, "cotBetaFromDet/F");
0325 
0326   pixeltreeOnTrack_->Branch("trkPt", &trkPt_);
0327   pixeltreeOnTrack_->Branch("trkEta", &trkEta_);
0328   pixeltreeOnTrack_->Branch("trkTheta", &trkTheta_);
0329   pixeltreeOnTrack_->Branch("trkPhi", &trkPhi_);
0330   pixeltreeOnTrack_->Branch("trkIsHighPurity", &trkIsHighPurity_);
0331   pixeltreeOnTrack_->Branch("cotAlphaFromTrack", &recHit_.cotAlphaFromTrack, "cotAlphaFromTrack/F");
0332   pixeltreeOnTrack_->Branch("cotBetaFromTrack", &recHit_.cotBetaFromTrack, "cotBetaFromTrack/F");
0333 
0334   pixeltreeOnTrack_->Branch("DgN", &recHit_.fDgN, "DgN/I");
0335   pixeltreeOnTrack_->Branch("DgRow", recHit_.fDgRow, "DgRow[DgN]/I");
0336   pixeltreeOnTrack_->Branch("DgCol", recHit_.fDgCol, "DgCol[DgN]/I");
0337   pixeltreeOnTrack_->Branch("DgDetId", recHit_.fDgDetId, "DgDetId[DgN]/I");
0338   pixeltreeOnTrack_->Branch("DgAdc", recHit_.fDgAdc, "DgAdc[DgN]/F");
0339   pixeltreeOnTrack_->Branch("DgCharge", recHit_.fDgCharge, "DgCharge[DgN]/F");
0340 }
0341 
0342 // Functions that gets called by framework every event
0343 void Phase2PixelNtuple::analyze(const edm::Event& e, const edm::EventSetup& es) {
0344   //Retrieve tracker topology from geometry
0345   const TrackerTopology* const tTopo = &es.getData(topo_esToken);
0346 
0347   // geometry setup
0348   const TrackerGeometry* theGeometry = &es.getData(geom_esToken);
0349 
0350   std::vector<PSimHit> matched;
0351   const PSimHit* closest_simhit = nullptr;
0352 
0353   edm::Handle<SiPixelRecHitCollection> recHitColl;
0354   e.getByToken(pixelRecHits_token_, recHitColl);
0355 
0356   // for finding matched simhit
0357   TrackerHitAssociator associate(e, trackerHitAssociatorConfig_);
0358 
0359   //Transient Rechit Builders
0360   const TkTransientTrackingRecHitBuilder* builder =
0361       static_cast<TkTransientTrackingRecHitBuilder const*>(&es.getData(tthb_esToken));
0362   auto hitCloner = builder->cloner();
0363 
0364   if ((recHitColl.product())->dataSize() > 0) {
0365     std::string detname;
0366 
0367     evt_.init();
0368     fillEvt(e);
0369 
0370     // Loop over Detector IDs
0371     for (auto recHitIdIterator : *(recHitColl.product())) {
0372       SiPixelRecHitCollection::DetSet detset = recHitIdIterator;
0373 
0374       if (detset.empty())
0375         continue;
0376       DetId detId = DetId(detset.detId());  // Get the Detid object
0377 
0378       const GeomDet* geomDet(theGeometry->idToDet(detId));
0379 
0380       // Loop over rechits for this detid
0381       for (auto iterRecHit : detset) {
0382         // get matched simhit
0383         matched.clear();
0384         matched = associate.associateHit(iterRecHit);
0385         if (!matched.empty()) {
0386           float closest = 9999.9;
0387           LocalPoint lp = iterRecHit.localPosition();
0388           float rechit_x = lp.x();
0389           float rechit_y = lp.y();
0390           //loop over simhits and find closest
0391           for (auto const& m : matched) {
0392             float sim_x1 = m.entryPoint().x();
0393             float sim_x2 = m.exitPoint().x();
0394             float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0395             float sim_y1 = m.entryPoint().y();
0396             float sim_y2 = m.exitPoint().y();
0397             float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0398 
0399             float x_res = sim_xpos - rechit_x;
0400             float y_res = sim_ypos - rechit_y;
0401             float dist = sqrt(x_res * x_res + y_res * y_res);
0402             if (dist < closest) {
0403               closest = dist;
0404               closest_simhit = &m;
0405             }
0406           }  // end of simhit loop
0407         }    // end matched emtpy
0408         unsigned int subid = detId.subdetId();
0409         int detid_db = detId.rawId();
0410         int layer_num = -99, ladder_num = -99, module_num = -99, disk_num = -99, blade_num = -99, panel_num = -99,
0411             side_num = -99;
0412         if ((subid == PixelSubdetector::PixelBarrel) || (subid == PixelSubdetector::PixelEndcap)) {
0413           // 1 = PXB, 2 = PXF
0414           if (subid == PixelSubdetector::PixelBarrel) {
0415             layer_num = tTopo->pxbLayer(detId.rawId());
0416             ladder_num = tTopo->pxbLadder(detId.rawId());
0417             module_num = tTopo->pxbModule(detId.rawId());
0418 #ifdef EDM_ML_DEBUG
0419             std::cout << "\ndetId = " << subid << " : " << tTopo->pxbLayer(detId.rawId()) << " , "
0420                       << tTopo->pxbLadder(detId.rawId()) << " , " << tTopo->pxbModule(detId.rawId());
0421 #endif
0422           } else if (subid == PixelSubdetector::PixelEndcap) {
0423             module_num = tTopo->pxfModule(detId());
0424             disk_num = tTopo->pxfDisk(detId());
0425             blade_num = tTopo->pxfBlade(detId());
0426             panel_num = tTopo->pxfPanel(detId());
0427             side_num = tTopo->pxfSide(detId());
0428           }
0429           int num_simhit = matched.size();
0430           recHit_.init();
0431           // filling in on ALL track rechits
0432           fillPRecHit(detid_db,
0433                       subid,
0434                       layer_num,
0435                       ladder_num,
0436                       module_num,
0437                       disk_num,
0438                       blade_num,
0439                       panel_num,
0440                       side_num,
0441                       &iterRecHit,
0442                       num_simhit,
0443                       closest_simhit,
0444                       geomDet);
0445           pixeltree_->Fill();
0446         }
0447       }  // end of rechit loop
0448     }    // end of detid loop
0449   }      // end of loop test on recHitColl size
0450 
0451   // Now loop over recotracks
0452   edm::Handle<View<reco::Track>> trackCollection;
0453   e.getByToken(recoTracks_token_, trackCollection);
0454 
0455   // -- Track trajectory association map
0456   edm::Handle<TrajTrackAssociationCollection> hTTAC;
0457   e.getByToken(tta_token_, hTTAC);
0458   TrajectoryStateCombiner tsoscomb;
0459 
0460   if (!trackCollection.isValid()) {
0461     if (picky_) {
0462       throw cms::Exception("ProductNotValid") << "TrackCollection product not valid";
0463     } else {
0464       std::cout << "TrackCollection product not valid" << endl;
0465       ;
0466     }
0467 
0468   } else if (!hTTAC.isValid()) {
0469     if (picky_) {
0470       throw cms::Exception("ProductNotValid") << "TrajectoryAssociationCollection product not valid";
0471     } else {
0472       std::cout << "TrajectoryAssociationCollection product not valid" << endl;
0473     }
0474 
0475   } else {
0476 #ifdef EDM_ML_DEBUG
0477     int rT = 0;
0478 #endif
0479     const TrajTrackAssociationCollection ttac = *(hTTAC.product());
0480     for (TrajTrackAssociationCollection::const_iterator it = ttac.begin(); it != ttac.end(); ++it) {
0481 #ifdef EDM_ML_DEBUG
0482       ++rT;
0483 #endif
0484       const edm::Ref<std::vector<Trajectory>> refTraj = it->key;
0485       auto track = it->val;
0486       trkIsHighPurity_ = track->quality(reco::TrackBase::highPurity);
0487       trkPt_ = track->pt();
0488       trkEta_ = track->eta();
0489       trkTheta_ = track->theta();
0490       trkPhi_ = track->phi();
0491 
0492 #ifdef EDM_ML_DEBUG
0493       std::cout << " num of hits for track " << rT << " = " << track->recHitsSize() << std::endl;
0494 #endif
0495 
0496       std::vector<TrajectoryMeasurement> tmeasColl = refTraj->measurements();
0497       for (auto const& tmeasIt : tmeasColl) {
0498         if (!tmeasIt.updatedState().isValid())
0499           continue;
0500         if (!tmeasIt.recHit()->isValid())
0501           continue;
0502 
0503         const TrackingRecHit* hit = tmeasIt.recHit()->hit();
0504         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
0505         if (pixhit == nullptr)
0506           continue;
0507         if (!pixhit->isValid())
0508           continue;
0509         TrajectoryStateOnSurface tsos = tsoscomb(tmeasIt.forwardPredictedState(), tmeasIt.backwardPredictedState());
0510         const DetId& detId = hit->geographicalId();
0511         const GeomDet* geomDet(theGeometry->idToDet(detId));
0512 
0513         if (pixhit) {
0514           // get matched simhit
0515           matched.clear();
0516           matched = associate.associateHit(*pixhit);
0517 
0518           if (!matched.empty()) {
0519             float closest = 9999.9;
0520             LocalPoint lp = pixhit->localPosition();
0521             float rechit_x = lp.x();
0522             float rechit_y = lp.y();
0523 
0524             //loop over simhits and find closest
0525             //for (std::vector<PSimHit>::const_iterator m = matched.begin(); m<matched.end(); m++)
0526             for (auto const& m : matched) {
0527               float sim_x1 = m.entryPoint().x();
0528               float sim_x2 = m.exitPoint().x();
0529               float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0530               float sim_y1 = m.entryPoint().y();
0531               float sim_y2 = m.exitPoint().y();
0532               float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0533 
0534               float x_res = sim_xpos - rechit_x;
0535               float y_res = sim_ypos - rechit_y;
0536               float dist = sqrt(x_res * x_res + y_res * y_res);
0537               if (dist < closest) {
0538                 closest = dist;
0539                 closest_simhit = &m;
0540               }
0541             }  // end of simhit loop
0542           }    // end matched emtpy
0543 
0544           int num_simhit = matched.size();
0545 
0546           int layer_num = -99, ladder_num = -99, module_num = -99, disk_num = -99, blade_num = -99, panel_num = -99,
0547               side_num = -99;
0548 
0549           unsigned int subid = detId.subdetId();
0550           int detid_db = detId.rawId();
0551           if ((subid == PixelSubdetector::PixelBarrel) || (subid == PixelSubdetector::PixelEndcap)) {
0552             // 1 = PXB, 2 = PXF
0553             if (subid == PixelSubdetector::PixelBarrel) {
0554               layer_num = tTopo->pxbLayer(detId.rawId());
0555               ladder_num = tTopo->pxbLadder(detId.rawId());
0556               module_num = tTopo->pxbModule(detId.rawId());
0557 #ifdef EDM_ML_DEBUG
0558               std::cout << "\ndetId = " << subid << " : " << tTopo->pxbLayer(detId.rawId()) << " , "
0559                         << tTopo->pxbLadder(detId.rawId()) << " , " << tTopo->pxbModule(detId.rawId()) << std::endl;
0560 #endif
0561             } else if (subid == PixelSubdetector::PixelEndcap) {
0562               module_num = tTopo->pxfModule(detId());
0563               disk_num = tTopo->pxfDisk(detId());
0564               blade_num = tTopo->pxfBlade(detId());
0565               panel_num = tTopo->pxfPanel(detId());
0566               side_num = tTopo->pxfSide(detId());
0567             }
0568 
0569             recHit_.init();
0570             // fill on track rechits
0571             fillPRecHit(detid_db,
0572                         subid,
0573                         layer_num,
0574                         ladder_num,
0575                         module_num,
0576                         disk_num,
0577                         blade_num,
0578                         panel_num,
0579                         side_num,
0580                         hit,     // TransientTrackingRecHit *
0581                         pixhit,  // SiPixelRecHit *
0582                         num_simhit,
0583                         closest_simhit,
0584                         geomDet,
0585                         tsos);
0586             pixeltreeOnTrack_->Fill();
0587           }  // if ( (subid==1)||(subid==2) )
0588         }    // if cast is possible to SiPixelHit
0589       }      //end of loop on tracking rechits
0590     }        // end of loop on recotracks
0591   }          // else track collection is valid
0592 }  // end analyze function
0593 
0594 // Function for filling in all the rechits
0595 // I know it is lazy to pass everything, but I'm doing it anyway. -EB
0596 void Phase2PixelNtuple::fillPRecHit(const int detid_db,
0597                                     const int subid,
0598                                     const int layer_num,
0599                                     const int ladder_num,
0600                                     const int module_num,
0601                                     const int disk_num,
0602                                     const int blade_num,
0603                                     const int panel_num,
0604                                     const int side_num,
0605                                     SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0606                                     const int num_simhit,
0607                                     const PSimHit* closest_simhit,
0608                                     const GeomDet* PixGeom) {
0609   LocalPoint lp = pixeliter->localPosition();
0610   LocalError le = pixeliter->localPositionError();
0611 
0612   recHit_.x = lp.x();
0613   recHit_.y = lp.y();
0614   recHit_.xx = le.xx();
0615   recHit_.xy = le.xy();
0616   recHit_.yy = le.yy();
0617   GlobalPoint GP = PixGeom->surface().toGlobal(lp);
0618   recHit_.gx = GP.x();
0619   recHit_.gy = GP.y();
0620   recHit_.gz = GP.z();
0621   GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0, 0, 0));
0622   recHit_.theta = GP0.theta();
0623   recHit_.phi = GP0.phi();
0624 
0625   SiPixelRecHit::ClusterRef const& Cluster = pixeliter->cluster();
0626   recHit_.q = Cluster->charge();
0627   recHit_.spreadx = Cluster->sizeX();
0628   recHit_.spready = Cluster->sizeY();
0629 
0630   recHit_.subid = subid;
0631   recHit_.nsimhit = num_simhit;
0632 
0633   recHit_.layer = layer_num;
0634   recHit_.ladder = ladder_num;
0635   recHit_.module = module_num;
0636   recHit_.module = module_num;
0637   recHit_.disk = disk_num;
0638   recHit_.blade = blade_num;
0639   recHit_.panel = panel_num;
0640   recHit_.side = side_num;
0641 
0642   /*-- module topology --*/
0643   const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(PixGeom);
0644   const PixelTopology* topol = &(theGeomDet->specificTopology());
0645   recHit_.nRowsInDet = topol->nrows();
0646   recHit_.nColsInDet = topol->ncolumns();
0647   recHit_.pitchx = topol->pitch().first;
0648   recHit_.pitchy = topol->pitch().second;
0649   recHit_.thickness = theGeomDet->surface().bounds().thickness();
0650 
0651   MeasurementPoint mp = topol->measurementPosition(LocalPoint(recHit_.x, recHit_.y));
0652   recHit_.row = mp.x();
0653   recHit_.col = mp.y();
0654 
0655   if (Cluster.isNonnull()) {
0656     // compute local angles from det position
0657     std::pair<float, float> local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet);
0658     recHit_.cotAlphaFromDet = local_angles.first;
0659     recHit_.cotBetaFromDet = local_angles.second;
0660 
0661     // -- Get digis of this cluster
0662     const std::vector<SiPixelCluster::Pixel>& pixvector = Cluster->pixels();
0663 #ifdef EDM_ML_DEBUG
0664     std::cout << "  Found " << pixvector.size() << " pixels for this cluster " << std::endl;
0665 #endif
0666     for (unsigned int i = 0; i < pixvector.size(); ++i) {
0667       if (recHit_.fDgN > DGPERCLMAX - 1)
0668         break;
0669       SiPixelCluster::Pixel holdpix = pixvector[i];
0670 
0671       recHit_.fDgRow[recHit_.fDgN] = holdpix.x;
0672       recHit_.fDgCol[recHit_.fDgN] = holdpix.y;
0673 #ifdef EDM_ML_DEBUG
0674       std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl;
0675 #endif
0676       recHit_.fDgDetId[recHit_.fDgN] = detid_db;
0677       recHit_.fDgAdc[recHit_.fDgN] = -99.;
0678       recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc / 1000.;
0679       ++recHit_.fDgN;
0680     }
0681   }  // if ( Cluster.isNonnull() )
0682 
0683 #ifdef EDM_ML_DEBUG
0684   std::cout << "num_simhit = " << num_simhit << std::endl;
0685 #endif
0686   if (num_simhit > 0) {
0687     recHit_.pdgid = closest_simhit->particleType();
0688     recHit_.process = closest_simhit->processType();
0689 
0690     float sim_x1 = closest_simhit->entryPoint().x();
0691     float sim_x2 = closest_simhit->exitPoint().x();
0692     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0693     float sim_y1 = closest_simhit->entryPoint().y();
0694     float sim_y2 = closest_simhit->exitPoint().y();
0695     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0696 
0697     float time_to_detid_ns = GP0.mag() / (CLHEP::c_light * CLHEP::ns / CLHEP::cm);  // speed of light in ns
0698     recHit_.ht = closest_simhit->timeOfFlight() - time_to_detid_ns;
0699 
0700     recHit_.tx = closest_simhit->localDirection().x();
0701     recHit_.ty = closest_simhit->localDirection().y();
0702     recHit_.tz = closest_simhit->localDirection().z();
0703 
0704     MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy));
0705     recHit_.hrow = hmp.x();
0706     recHit_.hcol = hmp.y();
0707 
0708     // Leaving the comment below, useful for future reference
0709     // alpha: angle with respect to local x axis in local (x,z) plane
0710     // float cotalpha = sim_xdir/sim_zdir;
0711     // beta: angle with respect to local y axis in local (y,z) plane
0712     // float cotbeta = sim_ydir/sim_zdir;
0713 
0714 #ifdef EDM_ML_DEBUG
0715     std::cout << "num_simhit x, y = " << 0.5 * (sim_x1 + sim_x2) << " " << 0.5 * (sim_y1 + sim_y2) << std::endl;
0716 #endif
0717   }
0718 #ifdef EDM_ML_DEBUG
0719   std::cout << "Found RecHit in " << subid
0720             << " global x/y/z : " << PixGeom->surface().toGlobal(pixeliter->localPosition()).x() << " "
0721             << PixGeom->surface().toGlobal(pixeliter->localPosition()).y() << " "
0722             << PixGeom->surface().toGlobal(pixeliter->localPosition()).z() << std::endl;
0723 #endif
0724 }
0725 
0726 // Function for filling in on track rechits
0727 void Phase2PixelNtuple::fillPRecHit(const int detid_db,
0728                                     const int subid,
0729                                     const int layer_num,
0730                                     const int ladder_num,
0731                                     const int module_num,
0732                                     const int disk_num,
0733                                     const int blade_num,
0734                                     const int panel_num,
0735                                     const int side_num,
0736                                     const TrackingRecHit* recHit,
0737                                     const SiPixelRecHit* pixHit,
0738                                     const int num_simhit,
0739                                     const PSimHit* closest_simhit,
0740                                     const GeomDet* PixGeom,
0741                                     const TrajectoryStateOnSurface tsos) {
0742   LocalPoint lp = recHit->localPosition();
0743   LocalError le = recHit->localPositionError();
0744 
0745   recHit_.x = lp.x();
0746   recHit_.y = lp.y();
0747   recHit_.xx = le.xx();
0748   recHit_.xy = le.xy();
0749   recHit_.yy = le.yy();
0750 
0751   recHit_.probQ = pixHit->probabilityQ();
0752   recHit_.probXY = pixHit->probabilityXY();
0753   //std::cout << "printing pixHit_.xxloc " << recHit_.xxloc << std::endl;
0754 
0755   GlobalPoint GP = PixGeom->surface().toGlobal(recHit->localPosition());
0756   recHit_.gx = GP.x();
0757   recHit_.gy = GP.y();
0758   recHit_.gz = GP.z();
0759   GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0, 0, 0));
0760   recHit_.theta = GP0.theta();
0761   recHit_.phi = GP0.phi();
0762   recHit_.subid = subid;
0763 
0764   SiPixelRecHit::ClusterRef const& Cluster = pixHit->cluster();
0765   recHit_.q = Cluster->charge();
0766   recHit_.spreadx = Cluster->sizeX();
0767   recHit_.spready = Cluster->sizeY();
0768 
0769   recHit_.nsimhit = num_simhit;
0770 
0771   recHit_.layer = layer_num;
0772   recHit_.ladder = ladder_num;
0773   recHit_.module = module_num;
0774   recHit_.module = module_num;
0775   recHit_.disk = disk_num;
0776   recHit_.blade = blade_num;
0777   recHit_.panel = panel_num;
0778   recHit_.side = side_num;
0779 
0780   /*-- module topology --*/
0781   const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(PixGeom);
0782   const PixelTopology* topol = &(theGeomDet->specificTopology());
0783   recHit_.nRowsInDet = topol->nrows();
0784   recHit_.nColsInDet = topol->ncolumns();
0785   recHit_.pitchx = topol->pitch().first;
0786   recHit_.pitchy = topol->pitch().second;
0787   recHit_.thickness = theGeomDet->surface().bounds().thickness();
0788 
0789   if (Cluster.isNonnull()) {
0790     // compute local angles from det position
0791     std::pair<float, float> local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet);
0792     recHit_.cotAlphaFromDet = local_angles.first;
0793     recHit_.cotBetaFromDet = local_angles.second;
0794 
0795     // compute local angles from track trajectory
0796     recHit_.cotAlphaFromTrack = tsos.localParameters().dxdz();
0797     recHit_.cotBetaFromTrack = tsos.localParameters().dydz();
0798 
0799     // -- Get digis of this cluster
0800     const std::vector<SiPixelCluster::Pixel>& pixvector = Cluster->pixels();
0801 #ifdef EDM_ML_DEBUG
0802     std::cout << "  Found " << pixvector.size() << " pixels for this cluster " << std::endl;
0803 #endif
0804     for (unsigned int i = 0; i < pixvector.size(); ++i) {
0805       if (recHit_.fDgN > DGPERCLMAX - 1)
0806         break;
0807       SiPixelCluster::Pixel holdpix = pixvector[i];
0808 
0809       recHit_.fDgRow[recHit_.fDgN] = holdpix.x;
0810       recHit_.fDgCol[recHit_.fDgN] = holdpix.y;
0811 #ifdef EDM_ML_DEBUG
0812       std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl;
0813 #endif
0814       recHit_.fDgDetId[recHit_.fDgN] = detid_db;
0815       recHit_.fDgAdc[recHit_.fDgN] = -99.;
0816       recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc / 1000.;
0817       ++recHit_.fDgN;
0818     }
0819   }  // if ( Cluster.isNonnull() )
0820 
0821   if (num_simhit > 0) {
0822     recHit_.pdgid = closest_simhit->particleType();
0823     recHit_.process = closest_simhit->processType();
0824 
0825     float sim_x1 = closest_simhit->entryPoint().x();
0826     float sim_x2 = closest_simhit->exitPoint().x();
0827     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0828     float sim_y1 = closest_simhit->entryPoint().y();
0829     float sim_y2 = closest_simhit->exitPoint().y();
0830     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0831 
0832     float time_to_detid_ns = GP0.mag() / (CLHEP::c_light * CLHEP::ns / CLHEP::cm);  // speed of light in ns
0833     recHit_.ht = closest_simhit->timeOfFlight() - time_to_detid_ns;
0834 
0835     recHit_.tx = closest_simhit->localDirection().x();
0836     recHit_.ty = closest_simhit->localDirection().y();
0837     recHit_.tz = closest_simhit->localDirection().z();
0838 
0839     MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy));
0840     recHit_.hrow = hmp.x();
0841     recHit_.hcol = hmp.y();
0842 
0843     // Leaving the comment below, useful for future reference
0844     // alpha: angle with respect to local x axis in local (x,z) plane
0845     // float cotalpha = sim_xdir/sim_zdir;
0846     // beta: angle with respect to local y axis in local (y,z) plane
0847     // float cotbeta = sim_ydir/sim_zdir;
0848 
0849 #ifdef EDM_ML_DEBUG
0850     std::cout << "num_simhit x, y = " << 0.5 * (sim_x1 + sim_x2) << " " << 0.5 * (sim_y1 + sim_y2) << std::endl;
0851 #endif
0852   }
0853 }
0854 
0855 void Phase2PixelNtuple::fillEvt(const edm::Event& E) {
0856   evt_.run = E.id().run();
0857   evt_.evtnum = E.id().event();
0858 }
0859 
0860 void Phase2PixelNtuple::evt::init() {
0861   int dummy_int = 9999;
0862   run = dummy_int;
0863   evtnum = dummy_int;
0864 }
0865 
0866 void Phase2PixelNtuple::RecHit::init() {
0867   float dummy_float = 9999.0;
0868 
0869   pdgid = 0;
0870   process = 0;
0871   q = dummy_float;
0872   x = dummy_float;
0873   y = dummy_float;
0874   xx = dummy_float;
0875   xy = dummy_float;
0876   yy = dummy_float;
0877 
0878   row = dummy_float;
0879   col = dummy_float;
0880   gx = dummy_float;
0881   gy = dummy_float;
0882   gz = dummy_float;
0883   nsimhit = 0;
0884   subid = -99;
0885   module = -99;
0886   layer = -99;
0887   ladder = -99;
0888   disk = -99;
0889   blade = -99;
0890   panel = -99;
0891   side = -99;
0892   spreadx = 0;
0893   spready = 0;
0894   hx = dummy_float;
0895   hy = dummy_float;
0896   ht = dummy_float;
0897   tx = dummy_float;
0898   ty = dummy_float;
0899   tz = dummy_float;
0900   theta = dummy_float;
0901   phi = dummy_float;
0902 
0903   fDgN = DGPERCLMAX;
0904   for (int i = 0; i < fDgN; ++i) {
0905     fDgRow[i] = fDgCol[i] = -9999;
0906     fDgAdc[i] = fDgCharge[i] = -9999.;
0907     //    fDgRoc[i] = fDgRocR[i] = fDgRocC[i] = -9999;
0908   }
0909   fDgN = 0;
0910 }
0911 std::pair<float, float> Phase2PixelNtuple::computeAnglesFromDetPosition(const SiPixelCluster& cl,
0912                                                                         const PixelTopology& theTopol,
0913                                                                         const GeomDetUnit& theDet) const {
0914   // get cluster center of gravity (of charge)
0915   float xcenter = cl.x();
0916   float ycenter = cl.y();
0917 
0918   // get the cluster position in local coordinates (cm)
0919 
0920   // ggiurgiu@jhu.edu 12/09/2010 : This function is called without track info, therefore there are no track
0921   // angles to provide here. Call the default localPosition (without track info)
0922   LocalPoint lp = theTopol.localPosition(MeasurementPoint(xcenter, ycenter));
0923   const Local3DPoint origin = theDet.surface().toLocal(GlobalPoint(0, 0, 0));  // can be computed once...
0924 
0925   auto gvx = lp.x() - origin.x();
0926   auto gvy = lp.y() - origin.y();
0927   auto gvz = -1.f / origin.z();
0928   // normalization not required as only ratio used...
0929 
0930   // calculate angles
0931   float cotalpha_ = gvx * gvz;
0932   float cotbeta_ = gvy * gvz;
0933 
0934   return std::make_pair(cotalpha_, cotbeta_);
0935 }
0936 
0937 //define this as a plug-in
0938 DEFINE_FWK_MODULE(Phase2PixelNtuple);