Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-14 02:38:16

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     int rT = 0;
0477     const TrajTrackAssociationCollection ttac = *(hTTAC.product());
0478     for (TrajTrackAssociationCollection::const_iterator it = ttac.begin(); it != ttac.end(); ++it) {
0479       ++rT;
0480       const edm::Ref<std::vector<Trajectory>> refTraj = it->key;
0481       auto track = it->val;
0482       trkIsHighPurity_ = track->quality(reco::TrackBase::highPurity);
0483       trkPt_ = track->pt();
0484       trkEta_ = track->eta();
0485       trkTheta_ = track->theta();
0486       trkPhi_ = track->phi();
0487 
0488       int iT = 0;
0489 #ifdef EDM_ML_DEBUG
0490       std::cout << " num of hits for track " << rT << " = " << track->recHitsSize() << std::endl;
0491 #endif
0492 
0493       std::vector<TrajectoryMeasurement> tmeasColl = refTraj->measurements();
0494       for (auto const& tmeasIt : tmeasColl) {
0495         if (!tmeasIt.updatedState().isValid())
0496           continue;
0497         if (!tmeasIt.recHit()->isValid())
0498           continue;
0499 
0500         const TrackingRecHit* hit = tmeasIt.recHit()->hit();
0501         const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(hit);
0502         if (pixhit == nullptr)
0503           continue;
0504         if (!pixhit->isValid())
0505           continue;
0506         ++iT;
0507         TrajectoryStateOnSurface tsos = tsoscomb(tmeasIt.forwardPredictedState(), tmeasIt.backwardPredictedState());
0508         const DetId& detId = hit->geographicalId();
0509         const GeomDet* geomDet(theGeometry->idToDet(detId));
0510 
0511         if (pixhit) {
0512           // get matched simhit
0513           matched.clear();
0514           matched = associate.associateHit(*pixhit);
0515 
0516           if (!matched.empty()) {
0517             float closest = 9999.9;
0518             LocalPoint lp = pixhit->localPosition();
0519             float rechit_x = lp.x();
0520             float rechit_y = lp.y();
0521 
0522             //loop over simhits and find closest
0523             //for (std::vector<PSimHit>::const_iterator m = matched.begin(); m<matched.end(); m++)
0524             for (auto const& m : matched) {
0525               float sim_x1 = m.entryPoint().x();
0526               float sim_x2 = m.exitPoint().x();
0527               float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0528               float sim_y1 = m.entryPoint().y();
0529               float sim_y2 = m.exitPoint().y();
0530               float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0531 
0532               float x_res = sim_xpos - rechit_x;
0533               float y_res = sim_ypos - rechit_y;
0534               float dist = sqrt(x_res * x_res + y_res * y_res);
0535               if (dist < closest) {
0536                 closest = dist;
0537                 closest_simhit = &m;
0538               }
0539             }  // end of simhit loop
0540           }    // end matched emtpy
0541 
0542           int num_simhit = matched.size();
0543 
0544           int layer_num = -99, ladder_num = -99, module_num = -99, disk_num = -99, blade_num = -99, panel_num = -99,
0545               side_num = -99;
0546 
0547           unsigned int subid = detId.subdetId();
0548           int detid_db = detId.rawId();
0549           if ((subid == PixelSubdetector::PixelBarrel) || (subid == PixelSubdetector::PixelEndcap)) {
0550             // 1 = PXB, 2 = PXF
0551             if (subid == PixelSubdetector::PixelBarrel) {
0552               layer_num = tTopo->pxbLayer(detId.rawId());
0553               ladder_num = tTopo->pxbLadder(detId.rawId());
0554               module_num = tTopo->pxbModule(detId.rawId());
0555 #ifdef EDM_ML_DEBUG
0556               std::cout << "\ndetId = " << subid << " : " << tTopo->pxbLayer(detId.rawId()) << " , "
0557                         << tTopo->pxbLadder(detId.rawId()) << " , " << tTopo->pxbModule(detId.rawId()) << std::endl;
0558 #endif
0559             } else if (subid == PixelSubdetector::PixelEndcap) {
0560               module_num = tTopo->pxfModule(detId());
0561               disk_num = tTopo->pxfDisk(detId());
0562               blade_num = tTopo->pxfBlade(detId());
0563               panel_num = tTopo->pxfPanel(detId());
0564               side_num = tTopo->pxfSide(detId());
0565             }
0566 
0567             recHit_.init();
0568             // fill on track rechits
0569             fillPRecHit(detid_db,
0570                         subid,
0571                         layer_num,
0572                         ladder_num,
0573                         module_num,
0574                         disk_num,
0575                         blade_num,
0576                         panel_num,
0577                         side_num,
0578                         hit,     // TransientTrackingRecHit *
0579                         pixhit,  // SiPixelRecHit *
0580                         num_simhit,
0581                         closest_simhit,
0582                         geomDet,
0583                         tsos);
0584             pixeltreeOnTrack_->Fill();
0585           }  // if ( (subid==1)||(subid==2) )
0586         }    // if cast is possible to SiPixelHit
0587       }      //end of loop on tracking rechits
0588     }        // end of loop on recotracks
0589   }          // else track collection is valid
0590 }  // end analyze function
0591 
0592 // Function for filling in all the rechits
0593 // I know it is lazy to pass everything, but I'm doing it anyway. -EB
0594 void Phase2PixelNtuple::fillPRecHit(const int detid_db,
0595                                     const int subid,
0596                                     const int layer_num,
0597                                     const int ladder_num,
0598                                     const int module_num,
0599                                     const int disk_num,
0600                                     const int blade_num,
0601                                     const int panel_num,
0602                                     const int side_num,
0603                                     SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0604                                     const int num_simhit,
0605                                     const PSimHit* closest_simhit,
0606                                     const GeomDet* PixGeom) {
0607   LocalPoint lp = pixeliter->localPosition();
0608   LocalError le = pixeliter->localPositionError();
0609 
0610   recHit_.x = lp.x();
0611   recHit_.y = lp.y();
0612   recHit_.xx = le.xx();
0613   recHit_.xy = le.xy();
0614   recHit_.yy = le.yy();
0615   GlobalPoint GP = PixGeom->surface().toGlobal(lp);
0616   recHit_.gx = GP.x();
0617   recHit_.gy = GP.y();
0618   recHit_.gz = GP.z();
0619   GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0, 0, 0));
0620   recHit_.theta = GP0.theta();
0621   recHit_.phi = GP0.phi();
0622 
0623   SiPixelRecHit::ClusterRef const& Cluster = pixeliter->cluster();
0624   recHit_.q = Cluster->charge();
0625   recHit_.spreadx = Cluster->sizeX();
0626   recHit_.spready = Cluster->sizeY();
0627 
0628   recHit_.subid = subid;
0629   recHit_.nsimhit = num_simhit;
0630 
0631   recHit_.layer = layer_num;
0632   recHit_.ladder = ladder_num;
0633   recHit_.module = module_num;
0634   recHit_.module = module_num;
0635   recHit_.disk = disk_num;
0636   recHit_.blade = blade_num;
0637   recHit_.panel = panel_num;
0638   recHit_.side = side_num;
0639 
0640   /*-- module topology --*/
0641   const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(PixGeom);
0642   const PixelTopology* topol = &(theGeomDet->specificTopology());
0643   recHit_.nRowsInDet = topol->nrows();
0644   recHit_.nColsInDet = topol->ncolumns();
0645   recHit_.pitchx = topol->pitch().first;
0646   recHit_.pitchy = topol->pitch().second;
0647   recHit_.thickness = theGeomDet->surface().bounds().thickness();
0648 
0649   MeasurementPoint mp = topol->measurementPosition(LocalPoint(recHit_.x, recHit_.y));
0650   recHit_.row = mp.x();
0651   recHit_.col = mp.y();
0652 
0653   if (Cluster.isNonnull()) {
0654     // compute local angles from det position
0655     std::pair<float, float> local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet);
0656     recHit_.cotAlphaFromDet = local_angles.first;
0657     recHit_.cotBetaFromDet = local_angles.second;
0658 
0659     // -- Get digis of this cluster
0660     const std::vector<SiPixelCluster::Pixel>& pixvector = Cluster->pixels();
0661 #ifdef EDM_ML_DEBUG
0662     std::cout << "  Found " << pixvector.size() << " pixels for this cluster " << std::endl;
0663 #endif
0664     for (unsigned int i = 0; i < pixvector.size(); ++i) {
0665       if (recHit_.fDgN > DGPERCLMAX - 1)
0666         break;
0667       SiPixelCluster::Pixel holdpix = pixvector[i];
0668 
0669       recHit_.fDgRow[recHit_.fDgN] = holdpix.x;
0670       recHit_.fDgCol[recHit_.fDgN] = holdpix.y;
0671 #ifdef EDM_ML_DEBUG
0672       std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl;
0673 #endif
0674       recHit_.fDgDetId[recHit_.fDgN] = detid_db;
0675       recHit_.fDgAdc[recHit_.fDgN] = -99.;
0676       recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc / 1000.;
0677       ++recHit_.fDgN;
0678     }
0679   }  // if ( Cluster.isNonnull() )
0680 
0681 #ifdef EDM_ML_DEBUG
0682   std::cout << "num_simhit = " << num_simhit << std::endl;
0683 #endif
0684   if (num_simhit > 0) {
0685     recHit_.pdgid = closest_simhit->particleType();
0686     recHit_.process = closest_simhit->processType();
0687 
0688     float sim_x1 = closest_simhit->entryPoint().x();
0689     float sim_x2 = closest_simhit->exitPoint().x();
0690     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0691     float sim_y1 = closest_simhit->entryPoint().y();
0692     float sim_y2 = closest_simhit->exitPoint().y();
0693     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0694 
0695     float time_to_detid_ns = GP0.mag() / (CLHEP::c_light * CLHEP::ns / CLHEP::cm);  // speed of light in ns
0696     recHit_.ht = closest_simhit->timeOfFlight() - time_to_detid_ns;
0697 
0698     recHit_.tx = closest_simhit->localDirection().x();
0699     recHit_.ty = closest_simhit->localDirection().y();
0700     recHit_.tz = closest_simhit->localDirection().z();
0701 
0702     MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy));
0703     recHit_.hrow = hmp.x();
0704     recHit_.hcol = hmp.y();
0705 
0706     // Leaving the comment below, useful for future reference
0707     // alpha: angle with respect to local x axis in local (x,z) plane
0708     // float cotalpha = sim_xdir/sim_zdir;
0709     // beta: angle with respect to local y axis in local (y,z) plane
0710     // float cotbeta = sim_ydir/sim_zdir;
0711 
0712 #ifdef EDM_ML_DEBUG
0713     std::cout << "num_simhit x, y = " << 0.5 * (sim_x1 + sim_x2) << " " << 0.5 * (sim_y1 + sim_y2) << std::endl;
0714 #endif
0715   }
0716 #ifdef EDM_ML_DEBUG
0717   std::cout << "Found RecHit in " << subid
0718             << " global x/y/z : " << PixGeom->surface().toGlobal(pixeliter->localPosition()).x() << " "
0719             << PixGeom->surface().toGlobal(pixeliter->localPosition()).y() << " "
0720             << PixGeom->surface().toGlobal(pixeliter->localPosition()).z() << std::endl;
0721 #endif
0722 }
0723 
0724 // Function for filling in on track rechits
0725 void Phase2PixelNtuple::fillPRecHit(const int detid_db,
0726                                     const int subid,
0727                                     const int layer_num,
0728                                     const int ladder_num,
0729                                     const int module_num,
0730                                     const int disk_num,
0731                                     const int blade_num,
0732                                     const int panel_num,
0733                                     const int side_num,
0734                                     const TrackingRecHit* recHit,
0735                                     const SiPixelRecHit* pixHit,
0736                                     const int num_simhit,
0737                                     const PSimHit* closest_simhit,
0738                                     const GeomDet* PixGeom,
0739                                     const TrajectoryStateOnSurface tsos) {
0740   LocalPoint lp = recHit->localPosition();
0741   LocalError le = recHit->localPositionError();
0742 
0743   recHit_.x = lp.x();
0744   recHit_.y = lp.y();
0745   recHit_.xx = le.xx();
0746   recHit_.xy = le.xy();
0747   recHit_.yy = le.yy();
0748 
0749   recHit_.probQ = pixHit->probabilityQ();
0750   recHit_.probXY = pixHit->probabilityXY();
0751   //std::cout << "printing pixHit_.xxloc " << recHit_.xxloc << std::endl;
0752 
0753   GlobalPoint GP = PixGeom->surface().toGlobal(recHit->localPosition());
0754   recHit_.gx = GP.x();
0755   recHit_.gy = GP.y();
0756   recHit_.gz = GP.z();
0757   GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0, 0, 0));
0758   recHit_.theta = GP0.theta();
0759   recHit_.phi = GP0.phi();
0760   recHit_.subid = subid;
0761 
0762   SiPixelRecHit::ClusterRef const& Cluster = pixHit->cluster();
0763   recHit_.q = Cluster->charge();
0764   recHit_.spreadx = Cluster->sizeX();
0765   recHit_.spready = Cluster->sizeY();
0766 
0767   recHit_.nsimhit = num_simhit;
0768 
0769   recHit_.layer = layer_num;
0770   recHit_.ladder = ladder_num;
0771   recHit_.module = module_num;
0772   recHit_.module = module_num;
0773   recHit_.disk = disk_num;
0774   recHit_.blade = blade_num;
0775   recHit_.panel = panel_num;
0776   recHit_.side = side_num;
0777 
0778   /*-- module topology --*/
0779   const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(PixGeom);
0780   const PixelTopology* topol = &(theGeomDet->specificTopology());
0781   recHit_.nRowsInDet = topol->nrows();
0782   recHit_.nColsInDet = topol->ncolumns();
0783   recHit_.pitchx = topol->pitch().first;
0784   recHit_.pitchy = topol->pitch().second;
0785   recHit_.thickness = theGeomDet->surface().bounds().thickness();
0786 
0787   if (Cluster.isNonnull()) {
0788     // compute local angles from det position
0789     std::pair<float, float> local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet);
0790     recHit_.cotAlphaFromDet = local_angles.first;
0791     recHit_.cotBetaFromDet = local_angles.second;
0792 
0793     // compute local angles from track trajectory
0794     recHit_.cotAlphaFromTrack = tsos.localParameters().dxdz();
0795     recHit_.cotBetaFromTrack = tsos.localParameters().dydz();
0796 
0797     // -- Get digis of this cluster
0798     const std::vector<SiPixelCluster::Pixel>& pixvector = Cluster->pixels();
0799 #ifdef EDM_ML_DEBUG
0800     std::cout << "  Found " << pixvector.size() << " pixels for this cluster " << std::endl;
0801 #endif
0802     for (unsigned int i = 0; i < pixvector.size(); ++i) {
0803       if (recHit_.fDgN > DGPERCLMAX - 1)
0804         break;
0805       SiPixelCluster::Pixel holdpix = pixvector[i];
0806 
0807       recHit_.fDgRow[recHit_.fDgN] = holdpix.x;
0808       recHit_.fDgCol[recHit_.fDgN] = holdpix.y;
0809 #ifdef EDM_ML_DEBUG
0810       std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl;
0811 #endif
0812       recHit_.fDgDetId[recHit_.fDgN] = detid_db;
0813       recHit_.fDgAdc[recHit_.fDgN] = -99.;
0814       recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc / 1000.;
0815       ++recHit_.fDgN;
0816     }
0817   }  // if ( Cluster.isNonnull() )
0818 
0819   if (num_simhit > 0) {
0820     recHit_.pdgid = closest_simhit->particleType();
0821     recHit_.process = closest_simhit->processType();
0822 
0823     float sim_x1 = closest_simhit->entryPoint().x();
0824     float sim_x2 = closest_simhit->exitPoint().x();
0825     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0826     float sim_y1 = closest_simhit->entryPoint().y();
0827     float sim_y2 = closest_simhit->exitPoint().y();
0828     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0829 
0830     float time_to_detid_ns = GP0.mag() / (CLHEP::c_light * CLHEP::ns / CLHEP::cm);  // speed of light in ns
0831     recHit_.ht = closest_simhit->timeOfFlight() - time_to_detid_ns;
0832 
0833     recHit_.tx = closest_simhit->localDirection().x();
0834     recHit_.ty = closest_simhit->localDirection().y();
0835     recHit_.tz = closest_simhit->localDirection().z();
0836 
0837     MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy));
0838     recHit_.hrow = hmp.x();
0839     recHit_.hcol = hmp.y();
0840 
0841     // Leaving the comment below, useful for future reference
0842     // alpha: angle with respect to local x axis in local (x,z) plane
0843     // float cotalpha = sim_xdir/sim_zdir;
0844     // beta: angle with respect to local y axis in local (y,z) plane
0845     // float cotbeta = sim_ydir/sim_zdir;
0846 
0847 #ifdef EDM_ML_DEBUG
0848     std::cout << "num_simhit x, y = " << 0.5 * (sim_x1 + sim_x2) << " " << 0.5 * (sim_y1 + sim_y2) << std::endl;
0849 #endif
0850   }
0851 }
0852 
0853 void Phase2PixelNtuple::fillEvt(const edm::Event& E) {
0854   evt_.run = E.id().run();
0855   evt_.evtnum = E.id().event();
0856 }
0857 
0858 void Phase2PixelNtuple::evt::init() {
0859   int dummy_int = 9999;
0860   run = dummy_int;
0861   evtnum = dummy_int;
0862 }
0863 
0864 void Phase2PixelNtuple::RecHit::init() {
0865   float dummy_float = 9999.0;
0866 
0867   pdgid = 0;
0868   process = 0;
0869   q = dummy_float;
0870   x = dummy_float;
0871   y = dummy_float;
0872   xx = dummy_float;
0873   xy = dummy_float;
0874   yy = dummy_float;
0875 
0876   row = dummy_float;
0877   col = dummy_float;
0878   gx = dummy_float;
0879   gy = dummy_float;
0880   gz = dummy_float;
0881   nsimhit = 0;
0882   subid = -99;
0883   module = -99;
0884   layer = -99;
0885   ladder = -99;
0886   disk = -99;
0887   blade = -99;
0888   panel = -99;
0889   side = -99;
0890   spreadx = 0;
0891   spready = 0;
0892   hx = dummy_float;
0893   hy = dummy_float;
0894   ht = dummy_float;
0895   tx = dummy_float;
0896   ty = dummy_float;
0897   tz = dummy_float;
0898   theta = dummy_float;
0899   phi = dummy_float;
0900 
0901   fDgN = DGPERCLMAX;
0902   for (int i = 0; i < fDgN; ++i) {
0903     fDgRow[i] = fDgCol[i] = -9999;
0904     fDgAdc[i] = fDgCharge[i] = -9999.;
0905     //    fDgRoc[i] = fDgRocR[i] = fDgRocC[i] = -9999;
0906   }
0907   fDgN = 0;
0908 }
0909 std::pair<float, float> Phase2PixelNtuple::computeAnglesFromDetPosition(const SiPixelCluster& cl,
0910                                                                         const PixelTopology& theTopol,
0911                                                                         const GeomDetUnit& theDet) const {
0912   // get cluster center of gravity (of charge)
0913   float xcenter = cl.x();
0914   float ycenter = cl.y();
0915 
0916   // get the cluster position in local coordinates (cm)
0917 
0918   // ggiurgiu@jhu.edu 12/09/2010 : This function is called without track info, therefore there are no track
0919   // angles to provide here. Call the default localPosition (without track info)
0920   LocalPoint lp = theTopol.localPosition(MeasurementPoint(xcenter, ycenter));
0921   const Local3DPoint origin = theDet.surface().toLocal(GlobalPoint(0, 0, 0));  // can be computed once...
0922 
0923   auto gvx = lp.x() - origin.x();
0924   auto gvy = lp.y() - origin.y();
0925   auto gvz = -1.f / origin.z();
0926   // normalization not required as only ratio used...
0927 
0928   // calculate angles
0929   float cotalpha_ = gvx * gvz;
0930   float cotbeta_ = gvy * gvz;
0931 
0932   return std::make_pair(cotalpha_, cotbeta_);
0933 }
0934 
0935 //define this as a plug-in
0936 DEFINE_FWK_MODULE(Phase2PixelNtuple);