Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-22 01:53:50 

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   int rT = 0;
0365   for (View<reco::Track>::size_type i = 0; i < trackCollection->size(); ++i) {
0366     ++rT;
0367     RefToBase<reco::Track> track(trackCollection, i);
0368     //      std::cout << " num of hits for track " << rT << " = " << track->recHitsSize() << std::endl;
0369     for (trackingRecHit_iterator ih = track->recHitsBegin(); ih != track->recHitsEnd(); ++ih) {
0370       TrackingRecHit* hit = (*ih)->clone();
0371       const DetId& detId = hit->geographicalId();
0372       const GeomDet* geomDet(theGeometry->idToDet(detId));
0373 
0374       /////comment out begin
0375       /*
0376         unsigned int subdetId = detId.subdetId();
0377         int layerNumber=0;
0378         int ringNumber = 0;
0379         int stereo = 0;
0380         std::string detname;
0381         if ( subdetId == StripSubdetector::TIB) {
0382           detname = "TIB";
0383           
0384           layerNumber = tTopo->tibLayer(detId.rawId);
0385           stereo = tTopo->tibStereo(detId.rawId);
0386         } else if ( subdetId ==  StripSubdetector::TOB ) {
0387           detname = "TOB";
0388       
0389           layerNumber = tTopo->tobLayer(detId.rawId);
0390           stereo = tTopo->tobStereo(detId.rawId);
0391         } else if ( subdetId ==  StripSubdetector::TID) {
0392           detname = "TID";
0393           
0394           layerNumber = tTopo->tidWheel(detId.rawId);
0395           ringNumber = tTopo->tidRing(detId.rawId);
0396           stereo = tTopo->tidStereo(detId.rawId);
0397         } else if ( subdetId ==  StripSubdetector::TEC ) {
0398           detname = "TEC";
0399           
0400           layerNumber = tTopo->tecWheel(detId.rawId);
0401           ringNumber = tTopo->tecRing(detId.rawId);
0402           stereo = tTopo->tecStereo(detId.rawId);
0403         } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0404           detname = "PXB";
0405           
0406           layerNumber = tTopo->pxbLayer(detId.rawId);
0407           stereo = 1;
0408         } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0409           detname = "PXF";
0410           
0411           layerNumber = tTopo->pxfDisk(detId.rawId);
0412           stereo = 1;
0413         }
0414 */
0415       //        std::cout << "RecHit in " << detname << " from detid " << detId.rawId()
0416       //                  << " subdet = " << subdetId
0417       //                  << " layer = " << layerNumber
0418       //                  << " Stereo = " << stereo
0419       //                  << std::endl;
0420       if (hit->isValid()) {
0421         unsigned int subid = detId.subdetId();
0422         if ((subid == 1) || (subid == 2)) {
0423           // 1 = PXB, 2 = PXF
0424           fillPRecHit(subid, ih, geomDet);
0425           fillEvt(e);
0426           pixeltree2_->Fill();
0427           init();
0428           /*
0429             TrackingRecHit * hit = (*ih)->clone();
0430             LocalPoint lp = hit->localPosition();
0431             LocalError le = hit->localPositionError();
0432 //            std::cout << "   lp x,y = " << lp.x() << " " << lp.y() << " lpe xx,xy,yy = "
0433 //                  << le.xx() << " " << le.xy() << " " << le.yy() << std::endl;
0434             std::cout << "Found RecHit in " << detname << " from detid " << detId.rawId()
0435         << " subdet = " << subdetId
0436         << " layer = " << layerNumber
0437                 << "global x/y/z/r = "
0438                  << geomDet->surface().toGlobal(lp).x() << " " 
0439                  << geomDet->surface().toGlobal(lp).y() << " " 
0440                  << geomDet->surface().toGlobal(lp).z() << " " 
0441                  << geomDet->surface().toGlobal(lp).perp() 
0442                 << " err x/y = " << sqrt(le.xx()) << " " << sqrt(le.yy()) << std::endl;
0443 */
0444         }
0445       }
0446       delete hit;
0447     }  //end of loop on tracking rechits
0448   }    // end of loop on recotracks
0449 
0450   // now for strip rechits
0451   edm::Handle<SiStripRecHit2DCollection> rechitsrphi;
0452   edm::Handle<SiStripRecHit2DCollection> rechitsstereo;
0453   edm::Handle<SiStripMatchedRecHit2DCollection> rechitsmatched;
0454   e.getByLabel(rphiRecHits_, rechitsrphi);
0455   e.getByLabel(stereoRecHits_, rechitsstereo);
0456   e.getByLabel(matchedRecHits_, rechitsmatched);
0457 
0458   //std::cout << " Step A: Standard Strip RPHI RecHits found " << rechitsrphi.product()->dataSize() << std::endl;
0459   //std::cout << " Step A: Standard Strip Stereo RecHits found " << rechitsstereo.product()->dataSize() << std::endl;
0460   //std::cout << " Step A: Standard Strip Matched RecHits found " << rechitsmatched.product()->dataSize() << std::endl;
0461   if (rechitsrphi->size() > 0) {
0462     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0463     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorBegin = rechitsrphi->begin();
0464     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorEnd   = rechitsrphi->end();
0465     //    SiStripRecHit2DCollectionOld::const_iterator iterRecHit;
0466     SiStripRecHit2DCollection::const_iterator recHitIdIterator = (rechitsrphi.product())->begin();
0467     SiStripRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsrphi.product())->end();
0468 
0469     std::string detname;
0470 
0471     //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0472     //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0473     // Loop over Detector IDs
0474     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0475       SiStripRecHit2DCollection::DetSet detset = *recHitIdIterator;
0476 
0477       if (detset.empty())
0478         continue;
0479       DetId detId = DetId(detset.detId());  // Get the Detid object
0480 
0481       //      const DetId& detId =  iterRecHit->geographicalId();
0482       const GeomDet* geomDet(theGeometry->idToDet(detId));
0483 
0484       // Loop over rechits for this detid
0485       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0486       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0487       SiStripRecHit2DCollection::DetSet::const_iterator iterRecHit;
0488       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0489         /////comment out begin
0490         /*
0491     unsigned int subdetId = detId.subdetId();
0492     int layerNumber=0;
0493     int ringNumber = 0;
0494     int stereo = 0;
0495     if ( subdetId == StripSubdetector::TIB) {
0496       detname = "TIB";
0497       
0498       layerNumber = tTopo->tibLayer(detId.rawId);
0499       stereo = tTopo->tibStereo(detId.rawId);
0500     } else if ( subdetId ==  StripSubdetector::TOB ) {
0501       detname = "TOB";
0502       
0503       layerNumber = tTopo->tobLayer(detId.rawId);
0504       stereo = tTopo->tobStereo(detId.rawId);
0505     } else if ( subdetId ==  StripSubdetector::TID) {
0506       detname = "TID";
0507       
0508       layerNumber = tTopo->tidWheel(detId.rawId);
0509       ringNumber = tTopo->tidRing(detId.rawId);
0510       stereo = tTopo->tidStereo(detId.rawId);
0511     } else if ( subdetId ==  StripSubdetector::TEC ) {
0512       detname = "TEC";
0513       
0514       layerNumber = tTopo->tecWheel(detId.rawId);
0515       ringNumber = tTopo->tecRing(detId.rawId);
0516       stereo = tTopo->tecStereo(detId.rawId);
0517     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0518       detname = "PXB";
0519       
0520       layerNumber = tTopo->pxbLayer(detId.rawId);
0521       stereo = 1;
0522     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0523       detname = "PXF";
0524       
0525       layerNumber = tTopo->pxfDisk(detId.rawId);
0526       stereo = 1;
0527     }
0528 */
0529         //      std::cout << "Found SiStripRPhiRecHit in " << detname << " from detid " << detId.rawId()
0530         //                << " subdet = " << subdetId
0531         //                << " layer = " << layerNumber
0532         //                << " Stereo = " << stereo
0533         //                << std::endl;
0534         //      std::cout << "Rechit global x/y/z/r : "
0535         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0536         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0537         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0538         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0539         //comment out end
0540         unsigned int subid = detId.subdetId();
0541         fillSRecHit(subid, iterRecHit, geomDet);
0542         fillEvt(e);
0543         striptree_->Fill();
0544         init();
0545       }  // end of rechit loop
0546     }    // end of detid loop
0547   }      // end of loop test on rechit size
0548 
0549   // now stereo hits
0550   if (rechitsstereo.product()->dataSize() > 0) {
0551     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0552     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorBegin = rechitsstereo->begin();
0553     //    SiStripRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorEnd   = rechitsstereo->end();
0554     //    SiStripRecHit2DCollectionOld::const_iterator iterRecHit;
0555     SiStripRecHit2DCollection::const_iterator recHitIdIterator = (rechitsstereo.product())->begin();
0556     SiStripRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsstereo.product())->end();
0557 
0558     std::string detname;
0559 
0560     // Loop over Detector IDs
0561     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0562       SiStripRecHit2DCollection::DetSet detset = *recHitIdIterator;
0563       //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0564       //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0565 
0566       if (detset.empty())
0567         continue;
0568       DetId detId = DetId(detset.detId());  // Get the Detid object
0569 
0570       //      const DetId& detId =  iterRecHit->geographicalId();
0571       const GeomDet* geomDet(theGeometry->idToDet(detId));
0572 
0573       // Loop over rechits for this detid
0574       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0575       SiStripRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0576       SiStripRecHit2DCollection::DetSet::const_iterator iterRecHit;
0577       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0578         /////comment out begin
0579 
0580         /*  unsigned int subdetId = detId.subdetId();
0581     int layerNumber=0;
0582     int ringNumber = 0;
0583     int stereo = 0;
0584     if ( subdetId == StripSubdetector::TIB) {
0585       detname = "TIB";
0586       
0587       layerNumber = tTopo->tibLayer(detId.rawId);
0588       stereo = tTopo->tibStereo(detId.rawId);
0589     } else if ( subdetId ==  StripSubdetector::TOB ) {
0590       detname = "TOB";
0591       
0592       layerNumber = tTopo->tobLayer(detId.rawId);
0593       stereo = tTopo->tobStereo(detId.rawId);
0594     } else if ( subdetId ==  StripSubdetector::TID) {
0595       detname = "TID";
0596       
0597       layerNumber = tTopo->tidWheel(detId.rawId);
0598       ringNumber = tTopo->tidRing(detId.rawId);
0599       stereo = tTopo->tidStereo(detId.rawId);
0600     } else if ( subdetId ==  StripSubdetector::TEC ) {
0601       detname = "TEC";
0602       
0603       layerNumber = tTopo->tecWheel(detId.rawId);
0604       ringNumber = tTopo->tecRing(detId.rawId);
0605       stereo = tTopo->tecStereo(detId.rawId);
0606     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0607       detname = "PXB";
0608       
0609       layerNumber = tTopo->pxbLayer(detId.rawId);
0610       stereo = 1;
0611     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0612       detname = "PXF";
0613       
0614       layerNumber = tTopo->pxfDisk(detId.rawId);
0615       stereo = 1;
0616     }
0617 */
0618         //      std::cout << "Found SiStripStereoRecHit in " << detname << " from detid " << detId.rawId()
0619         //                << " subdet = " << subdetId
0620         //                << " layer = " << layerNumber
0621         //                << " Stereo = " << stereo
0622         //                << std::endl;
0623         //      std::cout << "Rechit global x/y/z/r : "
0624         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0625         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0626         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0627         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0628         //comment out end
0629         unsigned int subid = detId.subdetId();
0630         fillSRecHit(subid, iterRecHit, geomDet);
0631         fillEvt(e);
0632         striptree_->Fill();
0633         init();
0634       }  // end of rechit loop
0635     }    // end of detid loop
0636   }      // end of loop test on rechit size
0637 
0638   // now matched hits
0639   if (rechitsmatched.product()->dataSize() > 0) {
0640     //Loop over all rechits in RPHI collection (can also loop only over DetId)
0641     //SiStripMatchedRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorBegin = rechitsmatched->begin();
0642     //SiStripMatchedRecHit2DCollectionOld::const_iterator theRecHitRangeIteratorEnd   = rechitsmatched->end();
0643     //SiStripMatchedRecHit2DCollectionOld::const_iterator iterRecHit;
0644     SiStripMatchedRecHit2DCollection::const_iterator recHitIdIterator = (rechitsmatched.product())->begin();
0645     SiStripMatchedRecHit2DCollection::const_iterator recHitIdIteratorEnd = (rechitsmatched.product())->end();
0646 
0647     std::string detname;
0648 
0649     // Loop over Detector IDs
0650     for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0651       SiStripMatchedRecHit2DCollection::DetSet detset = *recHitIdIterator;
0652 
0653       if (detset.empty())
0654         continue;
0655       DetId detId = DetId(detset.detId());  // Get the Detid object
0656 
0657       //    for ( iterRecHit = theRecHitRangeIteratorBegin;
0658       //          iterRecHit != theRecHitRangeIteratorEnd; ++iterRecHit) {
0659 
0660       //      const DetId& detId =  iterRecHit->geographicalId();
0661       const GeomDet* geomDet(theGeometry->idToDet(detId));
0662 
0663       // Loop over rechits for this detid
0664       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin();
0665       SiStripMatchedRecHit2DCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0666       SiStripMatchedRecHit2DCollection::DetSet::const_iterator iterRecHit;
0667       for (iterRecHit = rechitRangeIteratorBegin; iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) {
0668         /////comment out begin
0669         /*
0670     unsigned int subdetId = detId.subdetId();
0671     int layerNumber=0;
0672     int ringNumber = 0;
0673     int stereo = 0;
0674     if ( subdetId == StripSubdetector::TIB) {
0675       detname = "TIB";
0676       
0677       layerNumber = tTopo->tibLayer(detId.rawId);
0678       stereo = tTopo->tibStereo(detId.rawId);
0679     } else if ( subdetId ==  StripSubdetector::TOB ) {
0680       detname = "TOB";
0681       
0682       layerNumber = tTopo->tobLayer(detId.rawId);
0683       stereo = tTopo->tobStereo(detId.rawId);
0684     } else if ( subdetId ==  StripSubdetector::TID) {
0685       detname = "TID";
0686       
0687       layerNumber = tTopo->tidWheel(detId.rawId);
0688       ringNumber = tTopo->tidRing(detId.rawId);
0689       stereo = tTopo->tidStereo(detId.rawId);
0690     } else if ( subdetId ==  StripSubdetector::TEC ) {
0691       detname = "TEC";
0692       
0693       layerNumber = tTopo->tecWheel(detId.rawId);
0694       ringNumber = tTopo->tecRing(detId.rawId);
0695       stereo = tTopo->tecStereo(detId.rawId);
0696     } else if ( subdetId ==  PixelSubdetector::PixelBarrel ) {
0697       detname = "PXB";
0698       
0699       layerNumber = tTopo->pxbLayer(detId.rawId);
0700       stereo = 1;
0701     } else if ( subdetId ==  PixelSubdetector::PixelEndcap ) {
0702       detname = "PXF";
0703       
0704       layerNumber = tTopo->pxfDisk(detId.rawId);
0705       stereo = 1;
0706     }
0707     
0708 */
0709         //      std::cout << "Found SiStripMatchedRecHit in " << detname << " from detid " << detId.rawId()
0710         //                << " subdet = " << subdetId
0711         //                << " layer = " << layerNumber
0712         //                << " Stereo = " << stereo
0713         //                << std::endl;
0714         //      std::cout << "Rechit global x/y/z/r : "
0715         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).x() << " "
0716         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).y() << " "
0717         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).z() << " "
0718         //                << geomDet->surface().toGlobal(iterRecHit->localPosition()).perp() << std::endl;
0719         //comment out end
0720         unsigned int subid = detId.subdetId();
0721         fillSRecHit(subid, iterRecHit, geomDet);
0722         fillEvt(e);
0723         striptree_->Fill();
0724         init();
0725       }  // end of rechit loop
0726     }    // end of detidt loop
0727   }      // end of loop test on rechit size
0728 
0729 }  // end analyze function
0730 
0731 void StdHitNtuplizer::fillSRecHit(const int subid,
0732                                   SiStripRecHit2DCollection::DetSet::const_iterator pixeliter,
0733                                   const GeomDet* theGeom) {
0734   LocalPoint lp = pixeliter->localPosition();
0735   LocalError le = pixeliter->localPositionError();
0736 
0737   striprecHit_.x = lp.x();
0738   striprecHit_.y = lp.y();
0739   striprecHit_.xx = le.xx();
0740   striprecHit_.xy = le.xy();
0741   striprecHit_.yy = le.yy();
0742   GlobalPoint GP = theGeom->surface().toGlobal(pixeliter->localPosition());
0743   striprecHit_.gx = GP.x();
0744   striprecHit_.gy = GP.y();
0745   striprecHit_.gz = GP.z();
0746   striprecHit_.subid = subid;
0747 }
0748 void StdHitNtuplizer::fillSRecHit(const int subid,
0749                                   SiStripMatchedRecHit2DCollection::DetSet::const_iterator pixeliter,
0750                                   const GeomDet* theGeom) {
0751   LocalPoint lp = pixeliter->localPosition();
0752   LocalError le = pixeliter->localPositionError();
0753 
0754   striprecHit_.x = lp.x();
0755   striprecHit_.y = lp.y();
0756   striprecHit_.xx = le.xx();
0757   striprecHit_.xy = le.xy();
0758   striprecHit_.yy = le.yy();
0759   GlobalPoint GP = theGeom->surface().toGlobal(pixeliter->localPosition());
0760   striprecHit_.gx = GP.x();
0761   striprecHit_.gy = GP.y();
0762   striprecHit_.gz = GP.z();
0763   striprecHit_.subid = subid;
0764 }
0765 void StdHitNtuplizer::fillSRecHit(const int subid, const FastTrackerRecHit& hit, const GeomDet* theGeom) {
0766   LocalPoint lp = hit.localPosition();
0767   LocalError le = hit.localPositionError();
0768 
0769   striprecHit_.x = lp.x();
0770   striprecHit_.y = lp.y();
0771   striprecHit_.xx = le.xx();
0772   striprecHit_.xy = le.xy();
0773   striprecHit_.yy = le.yy();
0774   //MeasurementPoint mp = topol->measurementPosition(LocalPoint(striprecHit_.x, striprecHit_.y));
0775   //striprecHit_.row = mp.x();
0776   //striprecHit_.col = mp.y();
0777   GlobalPoint GP = theGeom->surface().toGlobal(hit.localPosition());
0778   striprecHit_.gx = GP.x();
0779   striprecHit_.gy = GP.y();
0780   striprecHit_.gz = GP.z();
0781   striprecHit_.subid = subid;
0782 }
0783 void StdHitNtuplizer::fillPRecHit(const int subid,
0784                                   const int layer_num,
0785                                   SiPixelRecHitCollection::DetSet::const_iterator pixeliter,
0786                                   const int num_simhit,
0787                                   std::vector<PSimHit>::const_iterator closest_simhit,
0788                                   const GeomDet* PixGeom) {
0789   LocalPoint lp = pixeliter->localPosition();
0790   LocalError le = pixeliter->localPositionError();
0791 
0792   recHit_.x = lp.x();
0793   recHit_.y = lp.y();
0794   recHit_.xx = le.xx();
0795   recHit_.xy = le.xy();
0796   recHit_.yy = le.yy();
0797   //MeasurementPoint mp = topol->measurementPosition(LocalPoint(recHit_.x, recHit_.y));
0798   //recHit_.row = mp.x();
0799   //recHit_.col = mp.y();
0800   GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition());
0801   recHit_.gx = GP.x();
0802   recHit_.gy = GP.y();
0803   recHit_.gz = GP.z();
0804   recHit_.subid = subid;
0805   recHit_.layer = layer_num;
0806   recHit_.nsimhit = num_simhit;
0807   //std::cout << "num_simhit = " << num_simhit << std::endl;
0808   if (num_simhit > 0) {
0809     float sim_x1 = (*closest_simhit).entryPoint().x();
0810     float sim_x2 = (*closest_simhit).exitPoint().x();
0811     recHit_.hx = 0.5 * (sim_x1 + sim_x2);
0812     float sim_y1 = (*closest_simhit).entryPoint().y();
0813     float sim_y2 = (*closest_simhit).exitPoint().y();
0814     recHit_.hy = 0.5 * (sim_y1 + sim_y2);
0815     //std::cout << "num_simhit x, y = " << 0.5*(sim_x1+sim_x2) << " " << 0.5*(sim_y1+sim_y2) << std::endl;
0816   }
0817   /*
0818        std::cout << "Found RecHit in " << subid
0819                  << " global x/y/z : "
0820                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).x() << " " 
0821                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).y() << " " 
0822                  << PixGeom->surface().toGlobal(pixeliter->localPosition()).z() << std::endl;
0823 */
0824 }
0825 void StdHitNtuplizer::fillPRecHit(const int subid, trackingRecHit_iterator ih, const GeomDet* PixGeom) {
0826   TrackingRecHit* pixeliter = (*ih)->clone();
0827   LocalPoint lp = pixeliter->localPosition();
0828   LocalError le = pixeliter->localPositionError();
0829 
0830   recHit_.x = lp.x();
0831   recHit_.y = lp.y();
0832   recHit_.xx = le.xx();
0833   recHit_.xy = le.xy();
0834   recHit_.yy = le.yy();
0835   GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition());
0836   recHit_.gx = GP.x();
0837   recHit_.gy = GP.y();
0838   recHit_.gz = GP.z();
0839   delete pixeliter;
0840   recHit_.subid = subid;
0841 }
0842 
0843 void StdHitNtuplizer::fillEvt(const edm::Event& E) {
0844   evt_.run = E.id().run();
0845   evt_.evtnum = E.id().event();
0846 }
0847 
0848 void StdHitNtuplizer::init() {
0849   evt_.init();
0850   recHit_.init();
0851   striprecHit_.init();
0852 }
0853 
0854 void StdHitNtuplizer::evt::init() {
0855   int dummy_int = 9999;
0856   run = dummy_int;
0857   evtnum = dummy_int;
0858 }
0859 
0860 void StdHitNtuplizer::RecHit::init() {
0861   float dummy_float = 9999.0;
0862 
0863   x = dummy_float;
0864   y = dummy_float;
0865   xx = dummy_float;
0866   xy = dummy_float;
0867   yy = dummy_float;
0868   row = dummy_float;
0869   col = dummy_float;
0870   gx = dummy_float;
0871   gy = dummy_float;
0872   gz = dummy_float;
0873   layer = 0;
0874   nsimhit = 0;
0875   hx = dummy_float;
0876   hy = dummy_float;
0877   tx = dummy_float;
0878   ty = dummy_float;
0879   theta = dummy_float;
0880   phi = dummy_float;
0881 }
0882 
0883 //define this as a plug-in
0884 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!