Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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