From c2fc610137c88964ec9f6ff7399b86f9132d682a Mon Sep 17 00:00:00 2001 From: Ernesto Migliore Date: Tue, 3 Mar 2015 00:28:28 +0100 Subject: [PATCH 1/2] add ntuplizer for pixel hits --- .../Geometry/plugins/BuildFile.xml | 14 + .../Geometry/plugins/StdPixelHitNtuplizer.cc | 923 ++++++++++++++++++ .../Geometry/test/StdPixelHitNtuplizer_cfg.py | 137 +++ 3 files changed, 1074 insertions(+) create mode 100644 SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml create mode 100644 SLHCUpgradeSimulations/Geometry/plugins/StdPixelHitNtuplizer.cc create mode 100644 SLHCUpgradeSimulations/Geometry/test/StdPixelHitNtuplizer_cfg.py diff --git a/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml b/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml new file mode 100644 index 0000000000000..a7e2456390993 --- /dev/null +++ b/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/SLHCUpgradeSimulations/Geometry/plugins/StdPixelHitNtuplizer.cc b/SLHCUpgradeSimulations/Geometry/plugins/StdPixelHitNtuplizer.cc new file mode 100644 index 0000000000000..710e076083864 --- /dev/null +++ b/SLHCUpgradeSimulations/Geometry/plugins/StdPixelHitNtuplizer.cc @@ -0,0 +1,923 @@ +/* + \class StdPixelHitNtuplizer +*/ + +// DataFormats +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSRecHit2DCollection.h" +#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" + +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h" + +#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h" +#include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h" +#include "DataFormats/TrackReco/interface/Track.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h" +#include "SimDataFormats/Track/interface/SimTrack.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" + +#include "FWCore/PluginManager/interface/ModuleDef.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +// Geometry +#include "MagneticField/Engine/interface/MagneticField.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" +#include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyBuilder.h" + +// For ROOT +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include + +// STD +#include +#include +#include + +using namespace std; +using namespace edm; +using namespace reco; + + +class StdPixelHitNtuplizer : public edm::EDAnalyzer +{ +public: + + explicit StdPixelHitNtuplizer(const edm::ParameterSet& conf); + virtual ~StdPixelHitNtuplizer(); + virtual void beginJob(); + virtual void endJob(); + virtual void analyze(const edm::Event& e, const edm::EventSetup& es); + +protected: + + void fillEvt(const edm::Event& ); + //void fillPRecHit(const int subid, SiPixelRecHitCollection::const_iterator pixeliter, + // const GeomDet* PixGeom); + void fillPRecHit(const int detid_db, const int subid, + const int layer_num,const int ladder_num,const int module_num, + const int disk_num,const int blade_num,const int panel_num,const int side_num, + SiPixelRecHitCollection::DetSet::const_iterator pixeliter, + const int num_simhit, + std::vector::const_iterator closest_simhit, + const GeomDet* PixGeom); + void fillPRecHit(const int detid_db, const int subid, + const int layer_num,const int ladder_num,const int module_num, + const int disk_num,const int blade_num,const int panel_num,const int side_num, + trackingRecHit_iterator pixeliter, + const int num_simhit, + std::vector::const_iterator closest_simhit, + const GeomDet* PixGeom); + std::pair computeAnglesFromDetPosition(const SiPixelCluster & cl, + const PixelTopology & top, + const GeomDetUnit & det ) const; + +private: + edm::ParameterSet conf_; + edm::InputTag src_; + bool verbose_; + bool picky_; + + static const int DGPERCLMAX = 100; + + //--- Structures for ntupling: + struct evt + { + int run; + int evtnum; + + void init(); + } evt_; + + struct RecHit + { + int pdgid; + int process; + float q; + float x; + float y; + float xx; + float xy; + float yy; + float row; + float col; + float gx; + float gy; + float gz; + int subid,module; + int layer,ladder; // BPix + int disk,blade,panel,side; // FPix + int nsimhit; + int spreadx,spready; + float hx, hy; + float hrow, hcol; + float tx, ty, tz; + float theta, phi; + + // detector topology + int nRowsInDet; + int nColsInDet; + float pitchx; + float pitchy; + float thickness; + + // local angles from det position + float cotAlphaFromDet, cotBetaFromDet; + + // digis + int fDgN; + int fDgRow[DGPERCLMAX], fDgCol[DGPERCLMAX]; + int fDgDetId[DGPERCLMAX]; + // int fDgRoc[DGPERCLMAX], fDgRocR[DGPERCLMAX], fDgRocC[DGPERCLMAX]; + float fDgAdc[DGPERCLMAX], fDgCharge[DGPERCLMAX]; + + void init(); + } recHit_; + + edm::Service tFileService; + TTree * pixeltree_; + TTree * pixeltreeOnTrack_; +}; + + +StdPixelHitNtuplizer::StdPixelHitNtuplizer(edm::ParameterSet const& conf) : + conf_(conf), + src_( conf.getParameter( "src" ) ), + verbose_( conf.getUntrackedParameter("verbose",false)), + picky_(conf.getUntrackedParameter("picky",false)), + pixeltree_(0), + pixeltreeOnTrack_(0) +{ +} + + +StdPixelHitNtuplizer::~StdPixelHitNtuplizer() { } + +void StdPixelHitNtuplizer::endJob() +{ + std::cout << " StdPixelHitNtuplizer::endJob" << std::endl; +} + + + +void StdPixelHitNtuplizer::beginJob() +{ + std::cout << " StdPixelHitNtuplizer::beginJob" << std::endl; + pixeltree_ = tFileService->make("PixelNtuple","Pixel hit analyzer ntuple"); + pixeltreeOnTrack_ = tFileService->make("PixelNtupleOnTrack","On-Track Pixel hit analyzer ntuple"); + + int bufsize = 64000; + //Common Branch + // pixeltree_->Branch("evt", &evt_, "run/I:evtnum/I", bufsize); + // pixeltree_->Branch("pixel_recHit", &recHit_, + // "pdgid/I:process:q/F:x:y:xx:xy:yy:row:col:gx:gy:gz:subid/I:module:layer:ladder:disk:blade:panel:side:nsimhit:spreadx:spready:hx/F:hy:tx:ty:tz:theta:phi:DgN/I:DgRow[DgN]/I:DgCol[DgN]/I:DgDetId[DgN]/I:DgAdc[DgN]/F:DgCharge[DgN]/F:DgClI[DgN]/I:ClN/I:ClDgN[ClN]/I:ClDgI[ClN][100]/I]", bufsize); + + //Common Branch + pixeltree_->Branch("evt", &evt_ , "run/I:evtnum/I", bufsize); + pixeltree_->Branch("pdgid", &recHit_.pdgid, "pdgid/I" ); + pixeltree_->Branch("process", &recHit_.process,"process/I"); + pixeltree_->Branch("q", &recHit_.q ,"q/F" ); + pixeltree_->Branch("x", &recHit_.x ,"x/F" ); + pixeltree_->Branch("y", &recHit_.y ,"y/F" ); + pixeltree_->Branch("xx", &recHit_.xx ,"xx/F" ); + pixeltree_->Branch("xy", &recHit_.xy ,"xy/F" ); + pixeltree_->Branch("yy", &recHit_.yy ,"yy/F" ); + pixeltree_->Branch("row", &recHit_.row ,"row/F" ); + pixeltree_->Branch("col", &recHit_.col ,"col/F" ); + pixeltree_->Branch("gx", &recHit_.gx ,"gx/F" ); + pixeltree_->Branch("gy", &recHit_.gy ,"gy/F" ); + pixeltree_->Branch("gz", &recHit_.gz ,"gz/F" ); + pixeltree_->Branch("subid", &recHit_.subid ,"subid/I" ); + pixeltree_->Branch("module", &recHit_.module ,"module/I" ); + pixeltree_->Branch("layer", &recHit_.layer ,"layer/I" ); + pixeltree_->Branch("ladder", &recHit_.ladder ,"ladder/I" ); + pixeltree_->Branch("disk", &recHit_.disk ,"disk/I" ); + pixeltree_->Branch("blade", &recHit_.blade ,"blade/I" ); + pixeltree_->Branch("panel", &recHit_.panel ,"panel/I" ); + pixeltree_->Branch("side", &recHit_.side ,"side/I" ); + pixeltree_->Branch("nsimhit", &recHit_.nsimhit ,"nsimhit/I"); + pixeltree_->Branch("spreadx", &recHit_.spreadx ,"spreadx/I"); + pixeltree_->Branch("spready", &recHit_.spready ,"spready/I"); + pixeltree_->Branch("hx", &recHit_.hx ,"hx/F"); + pixeltree_->Branch("hy", &recHit_.hy ,"hy/F"); + pixeltree_->Branch("hrow", &recHit_.hrow ,"hrow/F"); + pixeltree_->Branch("hcol", &recHit_.hcol ,"hcol/F"); + pixeltree_->Branch("tx", &recHit_.tx ,"tx/F"); + pixeltree_->Branch("ty", &recHit_.ty ,"ty/F"); + pixeltree_->Branch("tz", &recHit_.tz ,"tz/F"); + pixeltree_->Branch("theta", &recHit_.theta ,"theta/F" ); + pixeltree_->Branch("phi", &recHit_.phi ,"phi/F" ); + pixeltree_->Branch("nRowsInDet", &recHit_.nRowsInDet ,"nRowsInDet/I"); + pixeltree_->Branch("nColsInDet", &recHit_.nColsInDet ,"nColsInDet/I"); + pixeltree_->Branch("pitchx", &recHit_.pitchx ,"pitchx/F"); + pixeltree_->Branch("pitchy", &recHit_.pitchy ,"pitchy/F"); + pixeltree_->Branch("thickness", &recHit_.thickness ,"thickness/F"); + pixeltree_->Branch("cotAlphaFromDet", &recHit_.cotAlphaFromDet ,"cotAlphaFromDet/F"); + pixeltree_->Branch("cotBetaFromDet", &recHit_.cotBetaFromDet ,"cotBetaFraomDet/F"); + + pixeltree_->Branch("DgN", &recHit_.fDgN ,"DgN/I" ); + pixeltree_->Branch("DgRow", recHit_.fDgRow ,"DgRow[DgN]/I" ); + pixeltree_->Branch("DgCol", recHit_.fDgCol ,"DgCol[DgN]/I" ); + pixeltree_->Branch("DgDetId", recHit_.fDgDetId ,"DgDetId[DgN]/I" ); + pixeltree_->Branch("DgAdc", recHit_.fDgAdc ,"DgAdc[DgN]/F" ); + pixeltree_->Branch("DgCharge", recHit_.fDgCharge ,"DgCharge[DgN]/F"); + + //Common Branch (on-track) + pixeltreeOnTrack_->Branch("evt", &evt_ , "run/I:evtnum/I", bufsize); + pixeltreeOnTrack_->Branch("pdgid", &recHit_.pdgid, "pdgid/I" ); + pixeltreeOnTrack_->Branch("process", &recHit_.process,"process/I"); + pixeltreeOnTrack_->Branch("q", &recHit_.q ,"q/F" ); + pixeltreeOnTrack_->Branch("x", &recHit_.x ,"x/F" ); + pixeltreeOnTrack_->Branch("y", &recHit_.y ,"y/F" ); + pixeltreeOnTrack_->Branch("xx", &recHit_.xx ,"xx/F" ); + pixeltreeOnTrack_->Branch("xy", &recHit_.xy ,"xy/F" ); + pixeltreeOnTrack_->Branch("yy", &recHit_.yy ,"yy/F" ); + pixeltreeOnTrack_->Branch("row", &recHit_.row ,"row/F" ); + pixeltreeOnTrack_->Branch("col", &recHit_.col ,"col/F" ); + pixeltreeOnTrack_->Branch("gx", &recHit_.gx ,"gx/F" ); + pixeltreeOnTrack_->Branch("gy", &recHit_.gy ,"gy/F" ); + pixeltreeOnTrack_->Branch("gz", &recHit_.gz ,"gz/F" ); + pixeltreeOnTrack_->Branch("subid", &recHit_.subid ,"subid/I" ); + pixeltreeOnTrack_->Branch("module", &recHit_.module ,"module/I" ); + pixeltreeOnTrack_->Branch("layer", &recHit_.layer ,"layer/I" ); + pixeltreeOnTrack_->Branch("ladder", &recHit_.ladder ,"ladder/I" ); + pixeltreeOnTrack_->Branch("disk", &recHit_.disk ,"disk/I" ); + pixeltreeOnTrack_->Branch("blade", &recHit_.blade ,"blade/I" ); + pixeltreeOnTrack_->Branch("panel", &recHit_.panel ,"panel/I" ); + pixeltreeOnTrack_->Branch("side", &recHit_.side ,"side/I" ); + pixeltreeOnTrack_->Branch("nsimhit",&recHit_.nsimhit ,"nsimhit/I"); + pixeltreeOnTrack_->Branch("spreadx",&recHit_.spreadx ,"spreadx/I"); + pixeltreeOnTrack_->Branch("spready",&recHit_.spready ,"spready/I"); + pixeltreeOnTrack_->Branch("hx", &recHit_.hx ,"hx/F"); + pixeltreeOnTrack_->Branch("hy", &recHit_.hy ,"hy/F"); + pixeltreeOnTrack_->Branch("hrow", &recHit_.hrow ,"hrow/F"); + pixeltreeOnTrack_->Branch("hcol", &recHit_.hcol ,"hcol/F"); + pixeltreeOnTrack_->Branch("tx", &recHit_.tx ,"tx/F"); + pixeltreeOnTrack_->Branch("ty", &recHit_.ty ,"ty/F"); + pixeltreeOnTrack_->Branch("tz", &recHit_.tz ,"tz/F"); + pixeltreeOnTrack_->Branch("theta", &recHit_.theta ,"theta/F" ); + pixeltreeOnTrack_->Branch("phi", &recHit_.phi ,"phi/F" ); + pixeltreeOnTrack_->Branch("nRowsInDet", &recHit_.nRowsInDet ,"nRowsInDet/I"); + pixeltreeOnTrack_->Branch("nColsInDet", &recHit_.nColsInDet ,"nColsInDet/I"); + pixeltreeOnTrack_->Branch("pitchx", &recHit_.pitchx ,"pitchx/F"); + pixeltreeOnTrack_->Branch("pitchy", &recHit_.pitchy ,"pitchy/F"); + pixeltreeOnTrack_->Branch("thickness", &recHit_.thickness ,"thickness/F"); + pixeltreeOnTrack_->Branch("cotAlphaFromDet", &recHit_.cotAlphaFromDet ,"cotAlphaFromDet/F"); + pixeltreeOnTrack_->Branch("cotBetaFromDet", &recHit_.cotBetaFromDet ,"cotBetaFraomDet/F"); + + pixeltreeOnTrack_->Branch("DgN", &recHit_.fDgN ,"DgN/I" ); + pixeltreeOnTrack_->Branch("DgRow", recHit_.fDgRow ,"DgRow[DgN]/I" ); + pixeltreeOnTrack_->Branch("DgCol", recHit_.fDgCol ,"DgCol[DgN]/I" ); + pixeltreeOnTrack_->Branch("DgDetId", recHit_.fDgDetId ,"DgDetId[DgN]/I" ); + pixeltreeOnTrack_->Branch("DgAdc", recHit_.fDgAdc ,"DgAdc[DgN]/F" ); + pixeltreeOnTrack_->Branch("DgCharge", recHit_.fDgCharge ,"DgCharge[DgN]/F"); +} + +// Functions that gets called by framework every event +void StdPixelHitNtuplizer::analyze(const edm::Event& e, const edm::EventSetup& es) +{ + //Retrieve tracker topology from geometry + edm::ESHandle tTopoHandle; + es.get().get(tTopoHandle); + const TrackerTopology* const tTopo = tTopoHandle.product(); + + // geometry setup + edm::ESHandle geometry; + + es.get().get(geometry); + const TrackerGeometry* theGeometry = &(*geometry); + + // fastsim rechits + //edm::Handle theGSRecHits; + //edm::InputTag hitProducer; + //hitProducer = conf_.getParameter("HitProducer"); + //e.getByLabel(hitProducer, theGSRecHits); + + std::vector matched; + std::vector::const_iterator closest_simhit; + + edm::Handle recHitColl; + e.getByLabel( src_, recHitColl); + + // for finding matched simhit + TrackerHitAssociator associate( e, conf_ ); + + // std::cout << " Step A: Standard RecHits found " << (recHitColl.product())->dataSize() << std::endl; + if((recHitColl.product())->dataSize() > 0) { + SiPixelRecHitCollection::const_iterator recHitIdIterator = (recHitColl.product())->begin(); + SiPixelRecHitCollection::const_iterator recHitIdIteratorEnd = (recHitColl.product())->end(); + + std::string detname ; + + evt_.init(); + fillEvt(e); + + // Loop over Detector IDs + for ( ; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) { + SiPixelRecHitCollection::DetSet detset = *recHitIdIterator; + + if( detset.empty() ) continue; + DetId detId = DetId(detset.detId()); // Get the Detid object + + const GeomDet* geomDet( theGeometry->idToDet(detId) ); + + // Loop over rechits for this detid + SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorBegin = detset.begin(); + SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end(); + SiPixelRecHitCollection::DetSet::const_iterator iterRecHit; + for ( iterRecHit = rechitRangeIteratorBegin; + iterRecHit != rechitRangeIteratorEnd; ++iterRecHit) { + // get matched simhit + matched.clear(); + matched = associate.associateHit(*iterRecHit); + if ( !matched.empty() ) { + float closest = 9999.9; + std::vector::const_iterator closestit = matched.begin(); + LocalPoint lp = iterRecHit->localPosition(); + float rechit_x = lp.x(); + float rechit_y = lp.y(); + //loop over simhits and find closest + for (std::vector::const_iterator m = matched.begin(); mpxbLayer(detId.rawId()); + ladder_num = tTopo->pxbLadder(detId.rawId()); + module_num = tTopo->pxbModule(detId.rawId()); + // std::cout <<"\ndetId = "<pxbLayer(detId.rawId())<<" , "<pxbLadder(detId.rawId())<<" , "<< tTopo->pxbModule(detId.rawId()); + } else if ( subid == PixelSubdetector::PixelEndcap ) { + module_num = tTopo->pxfModule(detId()); + disk_num = tTopo->pxfDisk(detId()); + blade_num = tTopo->pxfBlade(detId()); + panel_num = tTopo->pxfPanel(detId()); + side_num = tTopo->pxfSide(detId()); + } + int num_simhit = matched.size(); + recHit_.init(); + fillPRecHit( detid_db, subid,layer_num,ladder_num,module_num,disk_num,blade_num,panel_num,side_num, + iterRecHit, num_simhit, closest_simhit, geomDet ); + pixeltree_->Fill(); + } + } // end of rechit loop + } // end of detid loop + } // end of loop test on recHitColl size + + // Now loop over recotracks + + edm::Handle > trackCollection; + edm::InputTag trackProducer; + trackProducer = conf_.getParameter("trackProducer"); + try {e.getByLabel(trackProducer, trackCollection);} + catch (cms::Exception& ex){ + if(verbose_){ + LogWarning("StdPixelHitNtuplizer") << "Track Collection" << trackProducer << "not found"; + } else {;} + } + + if (!trackCollection.isValid()){ + if(picky_){ + throw cms::Exception("ProductNotValid") << "TrackCollection product not valid"; + } else {;} + } else { + int rT = 0; + for(View::size_type i=0; isize(); ++i){ + ++rT; + RefToBase track(trackCollection, i); + int iT = 0; + //std::cout << " num of hits for track " << rT << " = " << track->recHitsSize() << std::endl; + for(trackingRecHit_iterator ih=track->recHitsBegin(); ih != track->recHitsEnd(); ++ih) { + ++iT; + //std::cout<<" analyzing hit n. "<clone(); + const DetId& detId = hit->geographicalId(); + const GeomDet* geomDet( theGeometry->idToDet(detId) ); + + const SiPixelRecHit *pixhit = dynamic_cast(hit); + + if(pixhit){ + if(pixhit->isValid() ) { + + // get matched simhit + matched.clear(); + matched = associate.associateHit(*pixhit); + + if ( !matched.empty() ) { + float closest = 9999.9; + std::vector::const_iterator closestit = matched.begin(); + LocalPoint lp = pixhit->localPosition(); + float rechit_x = lp.x(); + float rechit_y = lp.y(); + + //loop over simhits and find closest + for (std::vector::const_iterator m = matched.begin(); mpxbLayer(detId.rawId()); + ladder_num = tTopo->pxbLadder(detId.rawId()); + module_num = tTopo->pxbModule(detId.rawId()); + //std::cout <<"\ndetId = "<pxbLayer(detId.rawId())<<" , "<pxbLadder(detId.rawId())<<" , "<< tTopo->pxbModule(detId.rawId()) <pxfModule(detId()); + disk_num = tTopo->pxfDisk(detId()); + blade_num = tTopo->pxfBlade(detId()); + panel_num = tTopo->pxfPanel(detId()); + side_num = tTopo->pxfSide(detId()); + } + + recHit_.init(); + fillPRecHit(detid_db, subid, layer_num,ladder_num,module_num,disk_num,blade_num,panel_num,side_num, + ih, num_simhit, closest_simhit, geomDet ); + pixeltreeOnTrack_->Fill(); + /* + TrackingRecHit * hit = (*ih)->clone(); + LocalPoint lp = hit->localPosition(); + LocalError le = hit->localPositionError(); + // std::cout << " lp x,y = " << lp.x() << " " << lp.y() << " lpe xx,xy,yy = " + // << le.xx() << " " << le.xy() << " " << le.yy() << std::endl; + std::cout << "Found RecHit in " << detname << " from detid " << detId.rawId() + << " subdet = " << subdetId + << " layer = " << layerNumber + << "global x/y/z/r = " + << geomDet->surface().toGlobal(lp).x() << " " + << geomDet->surface().toGlobal(lp).y() << " " + << geomDet->surface().toGlobal(lp).z() << " " + << geomDet->surface().toGlobal(lp).perp() + << " err x/y = " << sqrt(le.xx()) << " " << sqrt(le.yy()) << std::endl; + */ + } // if ( (subid==1)||(subid==2) ) + } // if SiPixelHit is valid + } // if cast is possible to SiPixelHit + delete pixhit; + } //end of loop on tracking rechits + } // end of loop on recotracks + } // else track collection is valid +} // end analyze function + +// Function for filling in all the rechits +// I know it is lazy to pass everything, but I'm doing it anyway. -EB +void StdPixelHitNtuplizer::fillPRecHit(const int detid_db, const int subid, + const int layer_num,const int ladder_num,const int module_num, + const int disk_num,const int blade_num,const int panel_num,const int side_num, + SiPixelRecHitCollection::DetSet::const_iterator pixeliter, + const int num_simhit, + std::vector::const_iterator closest_simhit, + const GeomDet* PixGeom) +{ + LocalPoint lp = pixeliter->localPosition(); + LocalError le = pixeliter->localPositionError(); + + recHit_.x = lp.x(); + recHit_.y = lp.y(); + recHit_.xx = le.xx(); + recHit_.xy = le.xy(); + recHit_.yy = le.yy(); + GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition()); + recHit_.gx = GP.x(); + recHit_.gy = GP.y(); + recHit_.gz = GP.z(); + GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0,0,0)); + recHit_.theta = GP0.theta(); + recHit_.phi = GP0.phi(); + + SiPixelRecHit::ClusterRef const& Cluster = pixeliter->cluster(); + recHit_.q = Cluster->charge(); + recHit_.spreadx = Cluster->sizeX(); + recHit_.spready = Cluster->sizeY(); + + recHit_.subid = subid; + recHit_.nsimhit = num_simhit; + + recHit_.layer = layer_num; + recHit_.ladder= ladder_num; + recHit_.module= module_num; + recHit_.module= module_num; + recHit_.disk = disk_num; + recHit_.blade = blade_num; + recHit_.panel = panel_num; + recHit_.side = side_num; + + /*-- module topology --*/ + const PixelGeomDetUnit *theGeomDet = dynamic_cast (PixGeom ); + const PixelTopology * topol = &(theGeomDet->specificTopology()); + recHit_.nRowsInDet = topol->nrows(); + recHit_.nColsInDet = topol->ncolumns(); + recHit_.pitchx = topol->pitch().first; + recHit_.pitchy = topol->pitch().second; + recHit_.thickness = theGeomDet->surface().bounds().thickness(); + + MeasurementPoint mp = topol->measurementPosition(LocalPoint(recHit_.x, recHit_.y)); + recHit_.row = mp.x(); + recHit_.col = mp.y(); + + // if ( subid == 1 ) { + // std::cout << "all hits: BPIX Layer "<< layer_num << " RectangularPixelTopology " + // << " rechit row " << recHit_.row << + // << " , rechit col " << recHit_.col + // << " , nrow " << topol->nrows() << ", ncols " << topol->ncolumns() << ", pitchx " << topol->pitch().first << ", pitchy " << topol->pitch().second + // << std::endl; + // } + + if ( Cluster.isNonnull() ) { + + // compute local angles from det position + std::pair local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet); + recHit_.cotAlphaFromDet = local_angles.first; + recHit_.cotBetaFromDet = local_angles.second; + + + // -- Get digis of this cluster + const std::vector& pixvector = Cluster->pixels(); + // std::cout << " Found " << pixvector.size() << " pixels for this cluster " << std::endl; + for (unsigned int i = 0; i < pixvector.size(); ++i) { + if (recHit_.fDgN > DGPERCLMAX - 1) break; + SiPixelCluster::Pixel holdpix = pixvector[i]; + + recHit_.fDgRow[recHit_.fDgN] = holdpix.x; + recHit_.fDgCol[recHit_.fDgN] = holdpix.y; + // std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl; + recHit_.fDgDetId[recHit_.fDgN] = detid_db; + recHit_.fDgAdc[recHit_.fDgN] = -99.; + recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc/1000.; + ++recHit_.fDgN; + + } + } // if ( Cluster.isNonnull() ) + else + { + std::cout << "Pixel rechits with no associated cluster ?!?!?!!!!!!!!!!!!!!!!!!!!!!! " << std::endl; + // for (std::vector::const_iterator tmeasIt = tmeasColl.begin(); ... + } + /*-- --*/ + + //std::cout << "num_simhit = " << num_simhit << std::endl; + if(num_simhit > 0) { + + recHit_.pdgid = (*closest_simhit).particleType(); + recHit_.process= (*closest_simhit).processType(); + + float sim_x1 = (*closest_simhit).entryPoint().x(); + float sim_x2 = (*closest_simhit).exitPoint().x(); + recHit_.hx = 0.5*(sim_x1+sim_x2); + float sim_y1 = (*closest_simhit).entryPoint().y(); + float sim_y2 = (*closest_simhit).exitPoint().y(); + recHit_.hy = 0.5*(sim_y1+sim_y2); + + recHit_.tx = (*closest_simhit).localDirection().x(); + recHit_.ty = (*closest_simhit).localDirection().y(); + recHit_.tz = (*closest_simhit).localDirection().z(); + + MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy)); + recHit_.hrow = hmp.x(); + recHit_.hcol = hmp.y(); + + // alpha: angle with respect to local x axis in local (x,z) plane + // float cotalpha = sim_xdir/sim_zdir; + // beta: angle with respect to local y axis in local (y,z) plane + // float cotbeta = sim_ydir/sim_zdir; + + //std::cout << "num_simhit x, y = " << 0.5*(sim_x1+sim_x2) << " " << 0.5*(sim_y1+sim_y2) << std::endl; + } + /* + std::cout << "Found RecHit in " << subid + << " global x/y/z : " + << PixGeom->surface().toGlobal(pixeliter->localPosition()).x() << " " + << PixGeom->surface().toGlobal(pixeliter->localPosition()).y() << " " + << PixGeom->surface().toGlobal(pixeliter->localPosition()).z() << std::endl; + */ +} + +// Function for filling in on track rechits +void StdPixelHitNtuplizer::fillPRecHit(const int detid_db, const int subid, + const int layer_num,const int ladder_num,const int module_num, + const int disk_num,const int blade_num,const int panel_num,const int side_num, + trackingRecHit_iterator ih, + const int num_simhit, + std::vector::const_iterator closest_simhit, + const GeomDet* PixGeom) +{ + TrackingRecHit * pixeliter = (*ih)->clone(); + LocalPoint lp = pixeliter->localPosition(); + LocalError le = pixeliter->localPositionError(); + + recHit_.x = lp.x(); + recHit_.y = lp.y(); + recHit_.xx = le.xx(); + recHit_.xy = le.xy(); + recHit_.yy = le.yy(); + GlobalPoint GP = PixGeom->surface().toGlobal(pixeliter->localPosition()); + recHit_.gx = GP.x(); + recHit_.gy = GP.y(); + recHit_.gz = GP.z(); + GlobalPoint GP0 = PixGeom->surface().toGlobal(LocalPoint(0,0,0)); + recHit_.theta = GP0.theta(); + recHit_.phi = GP0.phi(); + recHit_.subid = subid; + + //std::cout<<"before getting the cluster"<(pixeliter)->cluster(); + recHit_.q = Cluster->charge(); + recHit_.spreadx = Cluster->sizeX(); + recHit_.spready = Cluster->sizeY(); + + recHit_.nsimhit = num_simhit; + + recHit_.layer = layer_num; + recHit_.ladder= ladder_num; + recHit_.module= module_num; + recHit_.module= module_num; + recHit_.disk = disk_num; + recHit_.blade = blade_num; + recHit_.panel = panel_num; + recHit_.side = side_num; + + /*-- module topology --*/ + const PixelGeomDetUnit *theGeomDet = dynamic_cast (PixGeom ); + const PixelTopology * topol = &(theGeomDet->specificTopology()); + recHit_.nRowsInDet = topol->nrows(); + recHit_.nColsInDet = topol->ncolumns(); + recHit_.pitchx = topol->pitch().first; + recHit_.pitchy = topol->pitch().second; + recHit_.thickness = theGeomDet->surface().bounds().thickness(); + + // if ( subid == 1 ) { + // std::cout << "on track only: BPIX Layer "<< layer_num << " RectangularPixelTopology " << "nrows " << topol->nrows() << ", ncols " << topol->ncolumns() << ", pitchx " << topol->pitch().first << ", pitchy " << topol->pitch().second << std::endl; + // } + + if ( Cluster.isNonnull() ) { + // compute local angles from det position + std::pair local_angles = computeAnglesFromDetPosition(*Cluster, *topol, *theGeomDet); + recHit_.cotAlphaFromDet = local_angles.first; + recHit_.cotBetaFromDet = local_angles.second; + + // -- Get digis of this cluster + const std::vector& pixvector = Cluster->pixels(); + // std::cout << " Found " << pixvector.size() << " pixels for this cluster " << std::endl; + for (unsigned int i = 0; i < pixvector.size(); ++i) { + if (recHit_.fDgN > DGPERCLMAX - 1) break; + SiPixelCluster::Pixel holdpix = pixvector[i]; + + recHit_.fDgRow[recHit_.fDgN] = holdpix.x; + recHit_.fDgCol[recHit_.fDgN] = holdpix.y; + // std::cout << "holdpix " << holdpix.x << " " << holdpix.y << std::endl; + recHit_.fDgDetId[recHit_.fDgN] = detid_db; + recHit_.fDgAdc[recHit_.fDgN] = -99.; + recHit_.fDgCharge[recHit_.fDgN] = holdpix.adc/1000.; + ++recHit_.fDgN; + } + } // if ( Cluster.isNonnull() ) + else + { + std::cout << "Pixel rechits with no associated cluster ?!?!?!!!!!!!!!!!!!!!!!!!!!!! " << std::endl; + // for (std::vector::const_iterator tmeasIt = tmeasColl.begin(); ... + } + + //std::cout << "num_simhit = " << num_simhit << std::endl; + if(num_simhit > 0) { + + recHit_.pdgid = (*closest_simhit).particleType(); + recHit_.process= (*closest_simhit).processType(); + + float sim_x1 = (*closest_simhit).entryPoint().x(); + float sim_x2 = (*closest_simhit).exitPoint().x(); + recHit_.hx = 0.5*(sim_x1+sim_x2); + float sim_y1 = (*closest_simhit).entryPoint().y(); + float sim_y2 = (*closest_simhit).exitPoint().y(); + recHit_.hy = 0.5*(sim_y1+sim_y2); + + recHit_.tx = (*closest_simhit).localDirection().x(); + recHit_.ty = (*closest_simhit).localDirection().y(); + recHit_.tz = (*closest_simhit).localDirection().z(); + + MeasurementPoint hmp = topol->measurementPosition(LocalPoint(recHit_.hx, recHit_.hy)); + recHit_.hrow = hmp.x(); + recHit_.hcol = hmp.y(); + + // alpha: angle with respect to local x axis in local (x,z) plane + // float cotalpha = sim_xdir/sim_zdir; + // beta: angle with respect to local y axis in local (y,z) plane + // float cotbeta = sim_ydir/sim_zdir; + + //std::cout << "num_simhit x, y = " << 0.5*(sim_x1+sim_x2) << " " << 0.5*(sim_y1+sim_y2) << std::endl; + } + + delete pixeliter; + +} + +void +StdPixelHitNtuplizer::fillEvt(const edm::Event& E) +{ + evt_.run = E.id().run(); + evt_.evtnum = E.id().event(); +} + + +void StdPixelHitNtuplizer::evt::init() +{ + int dummy_int = 9999; + run = dummy_int; + evtnum = dummy_int; +} + +void StdPixelHitNtuplizer::RecHit::init() +{ + float dummy_float = 9999.0; + + pdgid = 0; + process = 0; + q = dummy_float; + x = dummy_float; + y = dummy_float; + xx = dummy_float; + xy = dummy_float; + yy = dummy_float; + row = dummy_float; + col = dummy_float; + gx = dummy_float; + gy = dummy_float; + gz = dummy_float; + nsimhit = 0; + subid=-99; + module=-99; + layer=-99; + ladder=-99; + disk=-99; + blade=-99; + panel=-99; + side=-99; + spreadx=0; + spready=0; + hx = dummy_float; + hy = dummy_float; + tx = dummy_float; + ty = dummy_float; + tz = dummy_float; + theta = dummy_float; + phi = dummy_float; + + fDgN = DGPERCLMAX; + for (int i = 0; i < fDgN; ++i) { + fDgRow[i] = fDgCol[i] = -9999; + fDgAdc[i] = fDgCharge[i] = -9999.; + // fDgRoc[i] = fDgRocR[i] = fDgRocC[i] = -9999; + } + fDgN = 0; + +} + +std::pair StdPixelHitNtuplizer::computeAnglesFromDetPosition(const SiPixelCluster & cl, + const PixelTopology & theTopol, + const GeomDetUnit & theDet ) const + { + // EM check performed outside + // //--- This is a new det unit, so cache it + // PixelGeomDetUnit * theDet = dynamic_cast( &det ); + // if ( ! theDet ) + // { + // throw cms::Exception("StdPixelHitNtuplizer::computeAngleFromDetPosition") + // << " Wrong pointer to pixel detector !!!" << endl; + + // } + + // get cluster center of gravity (of charge) + float xcenter = cl.x(); + float ycenter = cl.y(); + + // get the cluster position in local coordinates (cm) + + // ggiurgiu@jhu.edu 12/09/2010 : This function is called without track info, therefore there are no track + // angles to provide here. Call the default localPosition (without track info) + LocalPoint lp = theTopol.localPosition( MeasurementPoint(xcenter, ycenter) ); + + // get the cluster position in global coordinates (cm) + GlobalPoint gp = theDet.surface().toGlobal( lp ); + float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() ); + + // normalize + float gpx = gp.x()/gp_mod; + float gpy = gp.y()/gp_mod; + float gpz = gp.z()/gp_mod; + + // make a global vector out of the global point; this vector will point from the + // origin of the detector to the cluster + GlobalVector gv(gpx, gpy, gpz); + + // make local unit vector along local X axis + const Local3DVector lvx(1.0, 0.0, 0.0); + + // get the unit X vector in global coordinates/ + GlobalVector gvx = theDet.surface().toGlobal( lvx ); + + // make local unit vector along local Y axis + const Local3DVector lvy(0.0, 1.0, 0.0); + + // get the unit Y vector in global coordinates + GlobalVector gvy = theDet.surface().toGlobal( lvy ); + + // make local unit vector along local Z axis + const Local3DVector lvz(0.0, 0.0, 1.0); + + // get the unit Z vector in global coordinates + GlobalVector gvz = theDet.surface().toGlobal( lvz ); + + // calculate the components of gv (the unit vector pointing to the cluster) + // in the local coordinate system given by the basis {gvx, gvy, gvz} + // note that both gv and the basis {gvx, gvy, gvz} are given in global coordinates + float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z(); + float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z(); + float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z(); + + + /* all the above is equivalent to + const Local3DPoint origin = theDet->surface().toLocal(GlobalPoint(0,0,0)); // can be computed once... + auto gvx = lp.x()-origin.x(); + auto gvy = lp.y()-origin.y(); + auto gvz = -origin.z(); + * normalization not required as only ratio used... + */ + + + // calculate angles + float cotalpha_ = gv_dot_gvx / gv_dot_gvz; + float cotbeta_ = gv_dot_gvy / gv_dot_gvz; + + return std::make_pair(cotalpha_, cotbeta_); + + } + + + //define this as a plug-in + DEFINE_FWK_MODULE(StdPixelHitNtuplizer); + diff --git a/SLHCUpgradeSimulations/Geometry/test/StdPixelHitNtuplizer_cfg.py b/SLHCUpgradeSimulations/Geometry/test/StdPixelHitNtuplizer_cfg.py new file mode 100644 index 0000000000000..4bb3e61696f3c --- /dev/null +++ b/SLHCUpgradeSimulations/Geometry/test/StdPixelHitNtuplizer_cfg.py @@ -0,0 +1,137 @@ +# Auto generated configuration file +# using: +# Revision: 1.172.2.5 +# Source: /cvs_server/repositories/CMSSW/CMSSW/Configuration/PyReleaseValidation/python/ConfigBuilder.py,v +# with command line options: step2 -s RECO -n 100 --conditions DESIGN_36_V10::All --datatier GEN-SIM-RECO --eventcontent RECOSIM --beamspot Gauss --fileout file:reco.root --filein file:raw.root --python_filename RecoMuon_Fullsim_cfg.py --no_exec +import FWCore.ParameterSet.Config as cms + +process = cms.Process('RECO') + +# import of standard configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load('Configuration.EventContent.EventContent_cff') + +process.load('SimGeneral.MixingModule.mixNoPU_cfi') + +process.load('Configuration.Geometry.GeometryExtended2023MuonReco_cff') +process.load('Configuration.Geometry.GeometryExtended2023Muon_cff') +process.load('Configuration.StandardSequences.MagneticField_38T_PostLS1_cff') + +process.load('Configuration.StandardSequences.Digi_cff') +process.load('Configuration.StandardSequences.SimL1Emulator_cff') +process.load('Configuration.StandardSequences.DigiToRaw_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + + +process.configurationMetadata = cms.untracked.PSet( + version = cms.untracked.string('$Revision: 1.14 $'), + annotation = cms.untracked.string('step2 nevts:100'), + name = cms.untracked.string('PyReleaseValidation') +) +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(False) + +) +# Input source +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + ' /store/relval/CMSSW_6_2_0_SLHC23_patch1/RelValMinBias_TuneZ2star_14TeV/GEN-SIM/PH2_1K_FB_V6_UPG2023Muon-v1/00000/22F2BC85-CAA7-E411-9C07-00248C0BE012.root' + ) +) + +# Output definition +process.output = cms.OutputModule("PoolOutputModule", + splitLevel = cms.untracked.int32(0), + outputCommands = process.RECOSIMEventContent.outputCommands, + fileName = cms.untracked.string('file:reco.root'), + dataset = cms.untracked.PSet( + dataTier = cms.untracked.string('GEN-SIM-RECO'), + filterName = cms.untracked.string('') + ) +) +#I'm only interested in the validation stuff +#process.output.outputCommands = cms.untracked.vstring('drop *','keep *_MEtoEDMConverter_*_*') + +#process.output = cms.OutputModule("PoolOutputModule", +# outputCommands = process.AODSIMEventContent.outputCommands, +# fileName = cms.untracked.string( +# 'file:/uscms_data/d2/brownson/slhc/quadMuon_RECO.root') +#) + + +# Additional output definition + +# Other statements +process.GlobalTag.globaltag = 'PH2_1K_FB_V6::All' + +### PhaseI Geometry and modifications ############################################### +# process.Timing = cms.Service("Timing") +## no playback when doing digis +#process.mix.playback = True +#process.MessageLogger.destinations = cms.untracked.vstring("detailedInfo_fullph1geom") + +### if pileup we need to set the number +#process.mix.input.nbPileupEvents = cms.PSet( +# averageNumber = cms.double(50.0) +#) +### if doing inefficiency at =50 +#process.simSiPixelDigis.AddPixelInefficiency = 20 + +process.ReadLocalMeasurement = cms.EDAnalyzer("StdPixelHitNtuplizer", + src = cms.InputTag("siPixelRecHits"), + ### if using simple (non-iterative) or old (as in 1_8_4) tracking + trackProducer = cms.InputTag("generalTracks"), + #verbose = cms.untracked.bool(True), + #picky = cms.untracked.bool(False), + ### for using track hit association + associatePixel = cms.bool(True), + associateStrip = cms.bool(False), + associateRecoTracks = cms.bool(False), + ROUList = cms.vstring('g4SimHitsTrackerHitsPixelBarrelLowTof', + 'g4SimHitsTrackerHitsPixelBarrelHighTof', + 'g4SimHitsTrackerHitsPixelEndcapLowTof', + 'g4SimHitsTrackerHitsPixelEndcapHighTof') +) + +# Path and EndPath definitions +process.digitisation_step = cms.Path(process.pdigi) +process.L1simulation_step = cms.Path(process.SimL1Emulator) + +process.digi2raw_step = cms.Path(process.DigiToRaw) +process.raw2digi_step = cms.Path(process.RawToDigi) +process.L1Reco_step = cms.Path(process.L1Reco) + +process.reconstruction_step = cms.Path(process.reconstruction) +process.user_step = cms.Path(process.ReadLocalMeasurement) +process.endjob_step = cms.Path(process.endOfProcess) +process.out_step = cms.EndPath(process.output) + +# Schedule definition +process.schedule = cms.Schedule(process.digitisation_step,process.L1simulation_step,process.digi2raw_step,process.raw2digi_step,process.L1Reco_step,process.reconstruction_step,process.user_step,process.endjob_step,process.out_step) + +# customisation of the process. +# Automatic addition of the customisation function from SLHCUpgradeSimulations.Configuration.combinedCustoms +from SLHCUpgradeSimulations.Configuration.combinedCustoms import cust_2023Muon +#call to customisation function cust_2023Muon imported from SLHCUpgradeSimulations.Configuration.combinedCustoms +process = cust_2023Muon(process) + +# Automatic addition of the customisation function from SimGeneral.MixingModule.fullMixCustomize_cff +from SimGeneral.MixingModule.fullMixCustomize_cff import setCrossingFrameOn +#call to customisation function setCrossingFrameOn imported from SimGeneral.MixingModule.fullMixCustomize_cff +process = setCrossingFrameOn(process) +# 620_SLHC17_patch1 Extended2023 MuonGEM added to digitizer +process.mix.mixObjects.mixSH.crossingFrames.extend(['MuonGEMHits', 'MuonME0Hits']) + +process.TFileService = cms.Service('TFileService', +fileName = cms.string("stdgrechitfullstdg_ntuple.root") +) +#print process.dumpPython() From 8fb963e9b38e6407dfe0b9e14c113b9459b2d7c7 Mon Sep 17 00:00:00 2001 From: Ernesto Migliore Date: Tue, 3 Mar 2015 14:36:06 +0100 Subject: [PATCH 2/2] fix buildfile in plugins to avoid multiple products --- SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml b/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml index a7e2456390993..d0ae229ee1617 100644 --- a/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml +++ b/SLHCUpgradeSimulations/Geometry/plugins/BuildFile.xml @@ -1,5 +1,5 @@ - +