Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /** \class StdHitNtuplizer
0002  * File: StdHitNtuplizer.cc
0003  * Authors: H. Cheung
0004  ************************************************************/
0005 
0006 // For ROOT
0007 #include <TROOT.h>
0008 #include <TTree.h>
0009 #include <TFile.h>
0010 #include <TF1.h>
0011 #include <TH2F.h>
0012 #include <TH1F.h>
0013 
0014 // STD
0015 #include <memory>
0016 #include <string>
0017 #include <iostream>
0018 
0019 // USER INCLUDES
0020 #include "DataFormats/Common/interface/DetSetVector.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/Common/interface/Ref.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0025 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0026 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0027 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0030 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0031 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0032 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0033 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/Framework/interface/EventSetup.h"
0036 #include "FWCore/Framework/interface/MakerMacros.h"
0037 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0038 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "FWCore/PluginManager/interface/ModuleDef.h"
0041 #include "FWCore/Utilities/interface/InputTag.h"
0042 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0043 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0044 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0045 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0046 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0047 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h"
0048 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0049 #include "MagneticField/Engine/interface/MagneticField.h"
0050 #include "SimDataFormats/Track/interface/SimTrack.h"
0051 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0052 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0053 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0054 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0055 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0056 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0057 
0058 class StdHitNtuplizer : public edm::one::EDAnalyzer<> {
0059 public:
0060   explicit StdHitNtuplizer(const edm::ParameterSet& conf);
0061   virtual ~StdHitNtuplizer();
0062   virtual void beginJob();
0063   virtual void endJob();
0064   virtual void analyze(const edm::Event& e, const edm::EventSetup& es);
0065 
0066 protected:
0067   void fillEvt(const edm::Event&);
0068   void fillSRecHit(const int subid,
0069                    SiStripRecHit2DCollection::DetSet::const_iterator pixeliter,
0070                    const GeomDet* theGeom);
0071   void fillSRecHit(const int subid,
0072                    SiStripMatchedRecHit2DCollection::DetSet::const_iterator pixeliter,
0073                    const GeomDet* theGeom);
0074   void fillSRecHit(const int subid, const FastTrackerRecHit& hit, const GeomDet* theGeom);
0075   //void fillPRecHit(const int subid, SiPixelRecHitCollection::const_iterator pixeliter,
0076   //                 const GeomDet* PixGeom);
0077   void fillPRecHit(const int subid,
0078                    const int layer_num,
0079                    SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0080                    const int num_simhit,
0081                    std::vector<PSimHit>::const_iterator closest_simhit,
0082                    const GeomDet* PixGeom);
0083   void fillPRecHit(const int subid, trackingRecHit_iterator pixeliter, const GeomDet* PixGeom);
0084 
0085 private:
0086   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geom_esToken;
0087   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topo_esToken;
0088   edm::ParameterSet conf_;
0089   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0090   edm::InputTag src_;
0091   edm::InputTag rphiRecHits_;
0092   edm::InputTag stereoRecHits_;
0093   edm::InputTag matchedRecHits_;
0094 
0095   void init();
0096 
0097   //--- Structures for ntupling:
0098   struct evt {
0099     int run;
0100     int evtnum;
0101 
0102     void init();
0103   } evt_, stripevt_;
0104 
0105   struct RecHit {
0106     float x;
0107     float y;
0108     float xx;
0109     float xy;
0110     float yy;
0111     float row;
0112     float col;
0113     float gx;
0114     float gy;
0115     float gz;
0116     int subid;
0117     int layer;
0118     int nsimhit;
0119     float hx, hy;
0120     float tx, ty;
0121     float theta, phi;
0122 
0123     void init();
0124   } recHit_, striprecHit_;
0125 
0126   TFile* tfile_;
0127   TTree* pixeltree_;
0128   TTree* striptree_;
0129   TTree* pixeltree2_;
0130 };
0131 
0132 using namespace std;
0133 using namespace edm;
0134 using namespace reco;
0135 
0136 StdHitNtuplizer::StdHitNtuplizer(edm::ParameterSet const& conf)
0137     : geom_esToken(esConsumes()),
0138       topo_esToken(esConsumes()),
0139       conf_(conf),
0140       trackerHitAssociatorConfig_(conf, consumesCollector()),
0141       src_(conf.getParameter<edm::InputTag>("src")),
0142       rphiRecHits_(conf.getParameter<edm::InputTag>("rphiRecHits")),
0143       stereoRecHits_(conf.getParameter<edm::InputTag>("stereoRecHits")),
0144       matchedRecHits_(conf.getParameter<edm::InputTag>("matchedRecHits")),
0145       tfile_(0),
0146       pixeltree_(0),
0147       striptree_(0),
0148       pixeltree2_(0) {}
0149 
0150 StdHitNtuplizer::~StdHitNtuplizer() = default;
0151 
0152 void StdHitNtuplizer::endJob() {
0153   std::cout << " StdHitNtuplizer::endJob" << std::endl;
0154   tfile_->Write();
0155   tfile_->Close();
0156 }
0157 
0158 void StdHitNtuplizer::beginJob() {
0159   std::cout << " StdHitNtuplizer::beginJob" << std::endl;
0160   std::string outputFile = conf_.getParameter<std::string>("OutputFile");
0161 
0162   tfile_ = new TFile(outputFile.c_str(), "RECREATE");
0163   pixeltree_ = new TTree("PixelNtuple", "Pixel hit analyzer ntuple");
0164   striptree_ = new TTree("StripNtuple", "Strip hit analyzer ntuple");
0165   pixeltree2_ = new TTree("Pixel2Ntuple", "Track Pixel hit analyzer ntuple");
0166 
0167   int bufsize = 64000;
0168 
0169   //Common Branch
0170   pixeltree_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize);
0171   pixeltree_->Branch("pixel_recHit",
0172                      &recHit_,
0173                      "x/F:y:xx:xy:yy:row:col:gx:gy:gz:subid/I:layer:nsimhit:hx/F:hy:tx:ty:theta:phi",
0174                      bufsize);
0175   pixeltree2_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize);
0176   pixeltree2_->Branch("pixel_recHit",
0177                       &recHit_,
0178                       "x/F:y:xx:xy:yy:row:col:gx:gy:gz:subid/I:layer:nsimhit:hx/F:hy:tx:ty:theta:phi",
0179                       bufsize);
0180 
0181   // Strip Branches
0182   striptree_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize);
0183   striptree_->Branch("strip_recHit",
0184                      &striprecHit_,
0185                      "x/F:y:xx:xy:yy:row:col:gx:gy:gz:subid/I:layer:nsimhit:hx/F:hy:tx:ty:theta:phi",
0186                      bufsize);
0187 }
0188 
0189 // Functions that gets called by framework every event
0190 void StdHitNtuplizer::analyze(const edm::Event& e, const edm::EventSetup& es) {
0191   //Retrieve tracker topology from geometry
0192   const TrackerTopology* const tTopo = &es.getData(topo_esToken);
0193 
0194   // geometry setup
0195   const TrackerGeometry* theGeometry = &es.getData(geom_esToken);
0196 
0197   // fastsim rechits
0198   //edm::Handle<SiTrackerGSRecHit2DCollection> theGSRecHits;
0199   //edm::InputTag hitProducer;
0200   //hitProducer = conf_.getParameter<edm::InputTag>("HitProducer");
0201   //e.getByLabel(hitProducer, theGSRecHits);
0202 
0203   edm::Handle<SiPixelRecHitCollection> recHitColl;
0204   e.getByLabel(src_, recHitColl);
0205 
0206   // for finding matched simhit
0207   TrackerHitAssociator associate(e, trackerHitAssociatorConfig_);
0208 
0209   //  std::cout << " Step A: Standard RecHits found " << (recHitColl.product())->dataSize() << std::endl;
0210   if ((recHitColl.product())->dataSize() > 0) {
0211     SiPixelRecHitCollection::const_iterator recHitIdIterator = (recHitColl.product())->begin();
0212     SiPixelRecHitCollection::const_iterator recHitIdIteratorEnd = (recHitColl.product())->end();
0213 
0214     std::string detname;
0215     std::vector<PSimHit> matched;
0216     std::vector<PSimHit>::const_iterator closest_simhit;
0217 
0218     // Loop over Detector IDs
0219     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0220       SiPixelRecHitCollection::DetSet detset = *recHitIdIterator;
0221 
0222       if (detset.empty())
0223         continue;
0224       DetId detId = DetId(detset.detId());  // Get the Detid object
0225       //unsigned int detType=detId.det();    // det type, tracker=1
0226 
0227       const GeomDet* geomDet(theGeometry->idToDet(detId));
0228 
0229       // Loop over rechits for this detid
0230       SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0231       SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0232       SiPixelRecHitCollection::DetSet::const_iterator iterRecHit;
0233       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0234         // get matched simhit
0235         matched.clear();
0236         matched = associate.associateHit(*iterRecHit);
0237         if (!matched.empty()) {
0238           float closest = 9999.9;
0239           std::vector<PSimHit>::const_iterator closestit = matched.begin();
0240           LocalPoint lp = iterRecHit->localPosition();
0241           float rechit_x = lp.x();
0242           float rechit_y = lp.y();
0243           //loop over simhits and find closest
0244           for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++) {
0245             float sim_x1 = (*m).entryPoint().x();
0246             float sim_x2 = (*m).exitPoint().x();
0247             float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0248             float sim_y1 = (*m).entryPoint().y();
0249             float sim_y2 = (*m).exitPoint().y();
0250             float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0251 
0252             float x_res = fabs(sim_xpos - rechit_x);
0253             float y_res = fabs(sim_ypos - rechit_y);
0254             float dist = sqrt(x_res * x_res + y_res * y_res);
0255             if (dist < closest) {
0256               closest = dist;
0257               closestit = m;
0258             }
0259           }  // end of simhit loop
0260           closest_simhit = closestit;
0261         }  // end matched emtpy
0262         /////comment out begin
0263         /*
0264     unsigned int subdetId = detId.subdetId();
0265     int layerNumber=0;
0266     int ringNumber = 0;
0267     int stereo = 0;
0268     if ( subdetId == StripSubdetector::TIB) {
0269       detname = "TIB";
0270       
0271       layerNumber = tTopo->tibLayer(detId.rawId);
0272       stereo = tTopo->tibStereo(detId.rawId);
0273     } else if ( subdetId ==  StripSubdetector::TOB ) {
0274       detname = "TOB";
0275       
0276       layerNumber = tTopo->tobLayer(detId.rawId);
0277       stereo = tTopo->tobStereo(detId.rawId);
0278     } else if ( subdetId ==  StripSubdetector::TID) {
0279       detname = "TID";
0280       
0281       layerNumber = tTopo->tidWheel(detId.rawId);
0282       ringNumber = tTopo->tidRing(detId.rawId);
0283       stereo = tTopo->tidStereo(detId.rawId);
0284     } else if ( subdetId ==  StripSubdetector::TEC ) {
0285       detname = "TEC";
0286       
0287       layerNumber = tTopo->tecWheel(detId.rawId);
0288       ringNumber = tTopo->tecRing(detId.rawId);
0289       stereo = tTopo->tecStereo(detId.rawId);
0290     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0291       detname = "PXB";
0292       
0293       layerNumber = tTopo->pxbLayer(detId.rawId);
0294       stereo = 1;
0295     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0296       detname = "PXF";
0297       
0298       layerNumber = tTopo->pxfDisk(detId.rawId);
0299       stereo = 1;
0300     }
0301     
0302         std::cout << "Found SiPixelRecHit in " << detname << " from detid " << detId.rawId()
0303                   << " subdet = " << subdetId
0304                   << " layer = " << layerNumber
0305                   << " Stereo = " << stereo
0306                   << std::endl;
0307         std::cout << "Rechit global x/y/z/r : "
0308                   << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " " 
0309                   << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " " 
0310                   << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " " 
0311                   << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0312 */
0313         //comment out end
0314         unsigned int subid = detId.subdetId();
0315         int layer_num = 0;
0316         if ((subid == 1) || (subid == 2)) {
0317           // 1 = PXB, 2 = PXF
0318           if (subid == PixelSubdetector::PixelBarrel) {
0319             layer_num = tTopo->pxbLayer(detId.rawId());
0320           } else if (subid == PixelSubdetector::PixelEndcap) {
0321             layer_num = tTopo->pxfDisk(detId.rawId());
0322           }
0323           int num_simhit = matched.size();
0324           fillPRecHit(subid, layer_num, iterRecHit, num_simhit, closest_simhit, geomDet);
0325           fillEvt(e);
0326           pixeltree_->Fill();
0327           init();
0328           // more info
0329           /*
0330           LocalPoint lp = iterRecHit->localPosition();
0331           LocalError le = iterRecHit->localPositionError();
0332           std::cout << "Filled SiPixelRecHit in " << detname << " from detid " << detId.rawId()
0333                     << " subdet = " << subdetId
0334                     << " layer = " << layerNumber
0335                     << "global x/y/z/r = "
0336                     << geomDet->surface().toGlobal(lp).x() << " " 
0337                     << geomDet->surface().toGlobal(lp).y() << " " 
0338                     << geomDet->surface().toGlobal(lp).z() << " " 
0339                     << geomDet->surface().toGlobal(lp).perp() 
0340                     << " err x/y = " << sqrt(le.xx()) << " " << sqrt(le.yy()) 
0341                     << " and num matched simhits = " << num_simhit << std::endl;
0342 */
0343         }
0344       }  // end of rechit loop
0345     }  // end of detid loop
0346   }  // end of loop test on recHitColl size
0347 
0348   // Now loop over recotracks
0349 
0350   edm::Handle<View<reco::Track> > trackCollection;
0351   edm::InputTag trackProducer;
0352   trackProducer = conf_.getParameter<edm::InputTag>("trackProducer");
0353   e.getByLabel(trackProducer, trackCollection);
0354 
0355   /*
0356   std::cout << " num of reco::Tracks with "
0357             << trackProducer.process()<<":"
0358             << trackProducer.label()<<":"
0359             << trackProducer.instance()
0360             << ": " << trackCollection->size() << "\n";
0361 */
0362   for (View<reco::Track>::size_type i = 0; i < trackCollection->size(); ++i) {
0363     RefToBase<reco::Track> track(trackCollection, i);
0364     for (trackingRecHit_iterator ih = track->recHitsBegin(); ih != track->recHitsEnd(); ++ih) {
0365       TrackingRecHit* hit = (*ih)->clone();
0366       const DetId& detId = hit->geographicalId();
0367       const GeomDet* geomDet(theGeometry->idToDet(detId));
0368 
0369       /////comment out begin
0370       /*
0371         unsigned int subdetId = detId.subdetId();
0372         int layerNumber=0;
0373         int ringNumber = 0;
0374         int stereo = 0;
0375         std::string detname;
0376         if ( subdetId == StripSubdetector::TIB) {
0377           detname = "TIB";
0378           
0379           layerNumber = tTopo->tibLayer(detId.rawId);
0380           stereo = tTopo->tibStereo(detId.rawId);
0381         } else if ( subdetId ==  StripSubdetector::TOB ) {
0382           detname = "TOB";
0383       
0384           layerNumber = tTopo->tobLayer(detId.rawId);
0385           stereo = tTopo->tobStereo(detId.rawId);
0386         } else if ( subdetId ==  StripSubdetector::TID) {
0387           detname = "TID";
0388           
0389           layerNumber = tTopo->tidWheel(detId.rawId);
0390           ringNumber = tTopo->tidRing(detId.rawId);
0391           stereo = tTopo->tidStereo(detId.rawId);
0392         } else if ( subdetId ==  StripSubdetector::TEC ) {
0393           detname = "TEC";
0394           
0395           layerNumber = tTopo->tecWheel(detId.rawId);
0396           ringNumber = tTopo->tecRing(detId.rawId);
0397           stereo = tTopo->tecStereo(detId.rawId);
0398         } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0399           detname = "PXB";
0400           
0401           layerNumber = tTopo->pxbLayer(detId.rawId);
0402           stereo = 1;
0403         } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0404           detname = "PXF";
0405           
0406           layerNumber = tTopo->pxfDisk(detId.rawId);
0407           stereo = 1;
0408         }
0409 */
0410       //        std::cout << "RecHit in " << detname << " from detid " << detId.rawId()
0411       //                  << " subdet = " << subdetId
0412       //                  << " layer = " << layerNumber
0413       //                  << " Stereo = " << stereo
0414       //                  << std::endl;
0415       if (hit->isValid()) {
0416         unsigned int subid = detId.subdetId();
0417         if ((subid == 1) || (subid == 2)) {
0418           // 1 = PXB, 2 = PXF
0419           fillPRecHit(subid, ih, geomDet);
0420           fillEvt(e);
0421           pixeltree2_->Fill();
0422           init();
0423           /*
0424             TrackingRecHit * hit = (*ih)->clone();
0425             LocalPoint lp = hit->localPosition();
0426             LocalError le = hit->localPositionError();
0427 //            std::cout << "   lp x,y = " << lp.x() << " " << lp.y() << " lpe xx,xy,yy = "
0428 //                  << le.xx() << " " << le.xy() << " " << le.yy() << std::endl;
0429             std::cout << "Found RecHit in " << detname << " from detid " << detId.rawId()
0430         << " subdet = " << subdetId
0431         << " layer = " << layerNumber
0432                 << "global x/y/z/r = "
0433                  << geomDet->surface().toGlobal(lp).x() << " " 
0434                  << geomDet->surface().toGlobal(lp).y() << " " 
0435                  << geomDet->surface().toGlobal(lp).z() << " " 
0436                  << geomDet->surface().toGlobal(lp).perp() 
0437                 << " err x/y = " << sqrt(le.xx()) << " " << sqrt(le.yy()) << std::endl;
0438 */
0439         }
0440       }
0441       delete hit;
0442     }  //end of loop on tracking rechits
0443   }  // end of loop on recotracks
0444 
0445   // now for strip rechits
0446   edm::Handle<SiStripRecHit2DCollection> rechitsrphi;
0447   edm::Handle<SiStripRecHit2DCollection> rechitsstereo;
0448   edm::Handle<SiStripMatchedRecHit2DCollection> rechitsmatched;
0449   e.getByLabel(rphiRecHits_, rechitsrphi);
0450   e.getByLabel(stereoRecHits_, rechitsstereo);
0451   e.getByLabel(matchedRecHits_, rechitsmatched);
0452 
0453   //std::cout << " Step A: Standard Strip RPHI RecHits found " << rechitsrphi.product()->dataSize() << std::endl;
0454   //std::cout << " Step A: Standard Strip Stereo RecHits found " << rechitsstereo.product()->dataSize() << std::endl;
0455   //std::cout << " Step A: Standard Strip Matched RecHits found " << rechitsmatched.product()->dataSize() << std::endl;
0456   if (rechitsrphi->size() > 0) {
0457     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0458     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorBegin = rechitsrphi->begin();
0459     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorEnd   = rechitsrphi->end();
0460     //    SiStripRecHit2DCollectionOld::const_iterator iterRecHit;
0461     SiStripRecHit2DCollection::const_iterator recHitIdIterator = (rechitsrphi.product())->begin();
0462     SiStripRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsrphi.product())->end();
0463 
0464     std::string detname;
0465 
0466     //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0467     //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0468     // Loop over Detector IDs
0469     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0470       SiStripRecHit2DCollection::DetSet detset = *recHitIdIterator;
0471 
0472       if (detset.empty())
0473         continue;
0474       DetId detId = DetId(detset.detId());  // Get the Detid object
0475 
0476       //      const DetId& detId =  iterRecHit->geographicalId();
0477       const GeomDet* geomDet(theGeometry->idToDet(detId));
0478 
0479       // Loop over rechits for this detid
0480       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0481       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0482       SiStripRecHit2DCollection::DetSet::const_iterator iterRecHit;
0483       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0484         /////comment out begin
0485         /*
0486     unsigned int subdetId = detId.subdetId();
0487     int layerNumber=0;
0488     int ringNumber = 0;
0489     int stereo = 0;
0490     if ( subdetId == StripSubdetector::TIB) {
0491       detname = "TIB";
0492       
0493       layerNumber = tTopo->tibLayer(detId.rawId);
0494       stereo = tTopo->tibStereo(detId.rawId);
0495     } else if ( subdetId ==  StripSubdetector::TOB ) {
0496       detname = "TOB";
0497       
0498       layerNumber = tTopo->tobLayer(detId.rawId);
0499       stereo = tTopo->tobStereo(detId.rawId);
0500     } else if ( subdetId ==  StripSubdetector::TID) {
0501       detname = "TID";
0502       
0503       layerNumber = tTopo->tidWheel(detId.rawId);
0504       ringNumber = tTopo->tidRing(detId.rawId);
0505       stereo = tTopo->tidStereo(detId.rawId);
0506     } else if ( subdetId ==  StripSubdetector::TEC ) {
0507       detname = "TEC";
0508       
0509       layerNumber = tTopo->tecWheel(detId.rawId);
0510       ringNumber = tTopo->tecRing(detId.rawId);
0511       stereo = tTopo->tecStereo(detId.rawId);
0512     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0513       detname = "PXB";
0514       
0515       layerNumber = tTopo->pxbLayer(detId.rawId);
0516       stereo = 1;
0517     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0518       detname = "PXF";
0519       
0520       layerNumber = tTopo->pxfDisk(detId.rawId);
0521       stereo = 1;
0522     }
0523 */
0524         //      std::cout << "Found SiStripRPhiRecHit in " << detname << " from detid " << detId.rawId()
0525         //                << " subdet = " << subdetId
0526         //                << " layer = " << layerNumber
0527         //                << " Stereo = " << stereo
0528         //                << std::endl;
0529         //      std::cout << "Rechit global x/y/z/r : "
0530         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0531         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0532         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0533         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0534         //comment out end
0535         unsigned int subid = detId.subdetId();
0536         fillSRecHit(subid, iterRecHit, geomDet);
0537         fillEvt(e);
0538         striptree_->Fill();
0539         init();
0540       }  // end of rechit loop
0541     }  // end of detid loop
0542   }  // end of loop test on rechit size
0543 
0544   // now stereo hits
0545   if (rechitsstereo.product()->dataSize() > 0) {
0546     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0547     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorBegin = rechitsstereo->begin();
0548     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorEnd   = rechitsstereo->end();
0549     //    SiStripRecHit2DCollectionOld::const_iterator iterRecHit;
0550     SiStripRecHit2DCollection::const_iterator recHitIdIterator = (rechitsstereo.product())->begin();
0551     SiStripRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsstereo.product())->end();
0552 
0553     std::string detname;
0554 
0555     // Loop over Detector IDs
0556     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0557       SiStripRecHit2DCollection::DetSet detset = *recHitIdIterator;
0558       //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0559       //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0560 
0561       if (detset.empty())
0562         continue;
0563       DetId detId = DetId(detset.detId());  // Get the Detid object
0564 
0565       //      const DetId& detId =  iterRecHit->geographicalId();
0566       const GeomDet* geomDet(theGeometry->idToDet(detId));
0567 
0568       // Loop over rechits for this detid
0569       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0570       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0571       SiStripRecHit2DCollection::DetSet::const_iterator iterRecHit;
0572       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0573         /////comment out begin
0574 
0575         /*  unsigned int subdetId = detId.subdetId();
0576     int layerNumber=0;
0577     int ringNumber = 0;
0578     int stereo = 0;
0579     if ( subdetId == StripSubdetector::TIB) {
0580       detname = "TIB";
0581       
0582       layerNumber = tTopo->tibLayer(detId.rawId);
0583       stereo = tTopo->tibStereo(detId.rawId);
0584     } else if ( subdetId ==  StripSubdetector::TOB ) {
0585       detname = "TOB";
0586       
0587       layerNumber = tTopo->tobLayer(detId.rawId);
0588       stereo = tTopo->tobStereo(detId.rawId);
0589     } else if ( subdetId ==  StripSubdetector::TID) {
0590       detname = "TID";
0591       
0592       layerNumber = tTopo->tidWheel(detId.rawId);
0593       ringNumber = tTopo->tidRing(detId.rawId);
0594       stereo = tTopo->tidStereo(detId.rawId);
0595     } else if ( subdetId ==  StripSubdetector::TEC ) {
0596       detname = "TEC";
0597       
0598       layerNumber = tTopo->tecWheel(detId.rawId);
0599       ringNumber = tTopo->tecRing(detId.rawId);
0600       stereo = tTopo->tecStereo(detId.rawId);
0601     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0602       detname = "PXB";
0603       
0604       layerNumber = tTopo->pxbLayer(detId.rawId);
0605       stereo = 1;
0606     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0607       detname = "PXF";
0608       
0609       layerNumber = tTopo->pxfDisk(detId.rawId);
0610       stereo = 1;
0611     }
0612 */
0613         //      std::cout << "Found SiStripStereoRecHit in " << detname << " from detid " << detId.rawId()
0614         //                << " subdet = " << subdetId
0615         //                << " layer = " << layerNumber
0616         //                << " Stereo = " << stereo
0617         //                << std::endl;
0618         //      std::cout << "Rechit global x/y/z/r : "
0619         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0620         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0621         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0622         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0623         //comment out end
0624         unsigned int subid = detId.subdetId();
0625         fillSRecHit(subid, iterRecHit, geomDet);
0626         fillEvt(e);
0627         striptree_->Fill();
0628         init();
0629       }  // end of rechit loop
0630     }  // end of detid loop
0631   }  // end of loop test on rechit size
0632 
0633   // now matched hits
0634   if (rechitsmatched.product()->dataSize() > 0) {
0635     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0636     SiStripMatchedRecHit2DCollection::const_iterator recHitIdIterator = (rechitsmatched.product())->begin();
0637     SiStripMatchedRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsmatched.product())->end();
0638 
0639     std::string detname;
0640 
0641     // Loop over Detector IDs
0642     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0643       SiStripMatchedRecHit2DCollection::DetSet detset = *recHitIdIterator;
0644 
0645       if (detset.empty())
0646         continue;
0647       DetId detId = DetId(detset.detId());  // Get the Detid object
0648 
0649       //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0650       //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0651 
0652       //      const DetId& detId =  iterRecHit->geographicalId();
0653       const GeomDet* geomDet(theGeometry->idToDet(detId));
0654 
0655       // Loop over rechits for this detid
0656       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0657       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0658       SiStripMatchedRecHit2DCollection::DetSet::const_iterator iterRecHit;
0659       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0660         /////comment out begin
0661         /*
0662     unsigned int subdetId = detId.subdetId();
0663     int layerNumber=0;
0664     int ringNumber = 0;
0665     int stereo = 0;
0666     if ( subdetId == StripSubdetector::TIB) {
0667       detname = "TIB";
0668       
0669       layerNumber = tTopo->tibLayer(detId.rawId);
0670       stereo = tTopo->tibStereo(detId.rawId);
0671     } else if ( subdetId ==  StripSubdetector::TOB ) {
0672       detname = "TOB";
0673       
0674       layerNumber = tTopo->tobLayer(detId.rawId);
0675       stereo = tTopo->tobStereo(detId.rawId);
0676     } else if ( subdetId ==  StripSubdetector::TID) {
0677       detname = "TID";
0678       
0679       layerNumber = tTopo->tidWheel(detId.rawId);
0680       ringNumber = tTopo->tidRing(detId.rawId);
0681       stereo = tTopo->tidStereo(detId.rawId);
0682     } else if ( subdetId ==  StripSubdetector::TEC ) {
0683       detname = "TEC";
0684       
0685       layerNumber = tTopo->tecWheel(detId.rawId);
0686       ringNumber = tTopo->tecRing(detId.rawId);
0687       stereo = tTopo->tecStereo(detId.rawId);
0688     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0689       detname = "PXB";
0690       
0691       layerNumber = tTopo->pxbLayer(detId.rawId);
0692       stereo = 1;
0693     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0694       detname = "PXF";
0695       
0696       layerNumber = tTopo->pxfDisk(detId.rawId);
0697       stereo = 1;
0698     }
0699     
0700 */
0701         //      std::cout << "Found SiStripMatchedRecHit in " << detname << " from detid " << detId.rawId()
0702         //                << " subdet = " << subdetId
0703         //                << " layer = " << layerNumber
0704         //                << " Stereo = " << stereo
0705         //                << std::endl;
0706         //      std::cout << "Rechit global x/y/z/r : "
0707         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0708         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0709         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0710         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0711         //comment out end
0712         unsigned int subid = detId.subdetId();
0713         fillSRecHit(subid, iterRecHit, geomDet);
0714         fillEvt(e);
0715         striptree_->Fill();
0716         init();
0717       }  // end of rechit loop
0718     }  // end of detidt loop
0719   }  // end of loop test on rechit size
0720 
0721 }  // end analyze function
0722 
0723 void StdHitNtuplizer::fillSRecHit(const int subid,
0724                                   SiStripRecHit2DCollection::DetSet::const_iterator pixeliter,
0725                                   const GeomDet* theGeom) {
0726   LocalPoint lp = pixeliter->localPosition();
0727   LocalError le = pixeliter->localPositionError();
0728 
0729   striprecHit_.x = lp.x();
0730   striprecHit_.y = lp.y();
0731   striprecHit_.xx = le.xx();
0732   striprecHit_.xy = le.xy();
0733   striprecHit_.yy = le.yy();
0734   GlobalPoint GP = theGeom->surface().toGlobal(pixeliter->localPosition());
0735   striprecHit_.gx = GP.x();
0736   striprecHit_.gy = GP.y();
0737   striprecHit_.gz = GP.z();
0738   striprecHit_.subid = subid;
0739 }
0740 void StdHitNtuplizer::fillSRecHit(const int subid,
0741                                   SiStripMatchedRecHit2DCollection::DetSet::const_iterator pixeliter,
0742                                   const GeomDet* theGeom) {
0743   LocalPoint lp = pixeliter->localPosition();
0744   LocalError le = pixeliter->localPositionError();
0745 
0746   striprecHit_.x = lp.x();
0747   striprecHit_.y = lp.y();
0748   striprecHit_.xx = le.xx();
0749   striprecHit_.xy = le.xy();
0750   striprecHit_.yy = le.yy();
0751   GlobalPoint GP = theGeom->surface().toGlobal(pixeliter->localPosition());
0752   striprecHit_.gx = GP.x();
0753   striprecHit_.gy = GP.y();
0754   striprecHit_.gz = GP.z();
0755   striprecHit_.subid = subid;
0756 }
0757 void StdHitNtuplizer::fillSRecHit(const int subid, const FastTrackerRecHit& hit, const GeomDet* theGeom) {
0758   LocalPoint lp = hit.localPosition();
0759   LocalError le = hit.localPositionError();
0760 
0761   striprecHit_.x = lp.x();
0762   striprecHit_.y = lp.y();
0763   striprecHit_.xx = le.xx();
0764   striprecHit_.xy = le.xy();
0765   striprecHit_.yy = le.yy();
0766   //MeasurementPoint mp = topol->measurementPosition(LocalPoint(striprecHit_.x, striprecHit_.y));
0767   //striprecHit_.row = mp.x();
0768   //striprecHit_.col = mp.y();
0769   GlobalPoint GP = theGeom->surface().toGlobal(hit.localPosition());
0770   striprecHit_.gx = GP.x();
0771   striprecHit_.gy = GP.y();
0772   striprecHit_.gz = GP.z();
0773   striprecHit_.subid = subid;
0774 }
0775 void StdHitNtuplizer::fillPRecHit(const int subid,
0776                                   const int layer_num,
0777                                   SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0778                                   const int num_simhit,
0779                                   std::vector<PSimHit>::const_iterator closest_simhit,
0780                                   const GeomDet* PixGeom) {
0781   LocalPoint lp = pixeliter->localPosition();
0782   LocalError le = pixeliter->localPositionError();
0783 
0784   recHit_.x = lp.x();
0785   recHit_.y = lp.y();
0786   recHit_.xx = le.xx();
0787   recHit_.xy = le.xy();
0788   recHit_.yy = le.yy();
0789   //MeasurementPoint mp = topol->measurementPosition(LocalPoint(recHit_.x, recHit_.y));
0790   //recHit_.row = mp.x();
0791   //recHit_.col = mp.y();
0792   GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition());
0793   recHit_.gx = GP.x();
0794   recHit_.gy = GP.y();
0795   recHit_.gz = GP.z();
0796   recHit_.subid = subid;
0797   recHit_.layer = layer_num;
0798   recHit_.nsimhit = num_simhit;
0799   //std::cout << "num_simhit = " << num_simhit << std::endl;
0800   if (num_simhit > 0) {
0801     float sim_x1 = (*closest_simhit).entryPoint().x();
0802     float sim_x2 = (*closest_simhit).exitPoint().x();
0803     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0804     float sim_y1 = (*closest_simhit).entryPoint().y();
0805     float sim_y2 = (*closest_simhit).exitPoint().y();
0806     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0807     //std::cout << "num_simhit x, y = " << 0.5*(sim_x1+sim_x2) << " " << 0.5*(sim_y1+sim_y2) << std::endl;
0808   }
0809   /*
0810        std::cout << "Found RecHit in " << subid
0811                  << " global x/y/z : "
0812                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).x() << " " 
0813                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).y() << " " 
0814                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).z() << std::endl;
0815 */
0816 }
0817 void StdHitNtuplizer::fillPRecHit(const int subid, trackingRecHit_iterator ih, const GeomDet* PixGeom) {
0818   TrackingRecHit* pixeliter = (*ih)->clone();
0819   LocalPoint lp = pixeliter->localPosition();
0820   LocalError le = pixeliter->localPositionError();
0821 
0822   recHit_.x = lp.x();
0823   recHit_.y = lp.y();
0824   recHit_.xx = le.xx();
0825   recHit_.xy = le.xy();
0826   recHit_.yy = le.yy();
0827   GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition());
0828   recHit_.gx = GP.x();
0829   recHit_.gy = GP.y();
0830   recHit_.gz = GP.z();
0831   delete pixeliter;
0832   recHit_.subid = subid;
0833 }
0834 
0835 void StdHitNtuplizer::fillEvt(const edm::Event& E) {
0836   evt_.run = E.id().run();
0837   evt_.evtnum = E.id().event();
0838 }
0839 
0840 void StdHitNtuplizer::init() {
0841   evt_.init();
0842   recHit_.init();
0843   striprecHit_.init();
0844 }
0845 
0846 void StdHitNtuplizer::evt::init() {
0847   int dummy_int = 9999;
0848   run = dummy_int;
0849   evtnum = dummy_int;
0850 }
0851 
0852 void StdHitNtuplizer::RecHit::init() {
0853   float dummy_float = 9999.0;
0854 
0855   x = dummy_float;
0856   y = dummy_float;
0857   xx = dummy_float;
0858   xy = dummy_float;
0859   yy = dummy_float;
0860   row = dummy_float;
0861   col = dummy_float;
0862   gx = dummy_float;
0863   gy = dummy_float;
0864   gz = dummy_float;
0865   layer = 0;
0866   nsimhit = 0;
0867   hx = dummy_float;
0868   hy = dummy_float;
0869   tx = dummy_float;
0870   ty = dummy_float;
0871   theta = dummy_float;
0872   phi = dummy_float;
0873 }
0874 
0875 //define this as a plug-in
0876 DEFINE_FWK_MODULE(StdHitNtuplizer);
[ 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!