Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:38:11 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!