/
githubmirror
/
cmssw
Обзор
Документация
Войти
/
githubmirror
/
cmssw
Код
Запросы
0
Пакеты
0
Релизы
0
Аналитика
Безопасность
master
RecoLocalTracker/SiPixelRecHits/test/ReadPixelRecHit.cc
687 строк
30 KB
Shahzad Malik Muzaffar
[RECONSTRUCTION] Drop Geometry/CommonDetUnit package
28 май 2026, 12:57
28 май 2026, 12:57
685f036
Код
Авторство
О чём код?
//-------------------------------------------- // File: ReadPixelRecHit.cc // Description: see ReadPixelRecHit.h // Author: J.Sheav (JHU) // 11/8/06: New loop over rechits and InputTag, V.Chiochia // Creation Date: OGU Aug. 1 2005 Initial version. // Add occupancy histograms. D.Kotlinski. 10/06 //-------------------------------------------- // system includes #include <memory> #include <string> #include <iostream> // user includes #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "Geometry/CommonTopologies/interface/GeomDetType.h" #include "Geometry/CommonTopologies/interface/GeomDet.h" #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h" #include "Geometry/CommonTopologies/interface/PixelGeomDetType.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "Geometry/CommonTopologies/interface/PixelTopology.h" #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h" #include "DataFormats/DetId/interface/DetId.h" #include "DataFormats/SiPixelDetId/interface/PXBDetId.h" #include "DataFormats/SiPixelDetId/interface/PXFDetId.h" #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" #include "FWCore/Framework/interface/one/EDAnalyzer.h" #include "FWCore/Framework/interface/Event.h" #include "DataFormats/Common/interface/Handle.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #define DO_HISTO #ifdef DO_HISTO // For ROOT #include <TROOT.h> #include <TChain.h> #include <TFile.h> #include <TF1.h> #include <TH2F.h> #include <TH1F.h> #include <TProfile.h> #endif class ReadPixelRecHit : public edm::one::EDAnalyzer<edm::one::SharedResources> { public: explicit ReadPixelRecHit(const edm::ParameterSet &conf); virtual ~ReadPixelRecHit(); virtual void analyze(const edm::Event &e, const edm::EventSetup &c); virtual void beginJob(); virtual void endJob(); private: const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> esTokenGeom_; edm::InputTag src_; const edm::EDGetTokenT<SiPixelRecHitCollection> siPixelRecHitsToken_; bool print; #ifdef DO_HISTO TFile *hFile; TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id; TH1F *hcharge1, *hcharge2, *hcharge3; TH1F *hadcCharge1, *hadcCharge2, *hadcCharge3, *hadcCharge1big; TH1F *hxpos1, *hxpos2, *hxpos3, *hypos1, *hypos2, *hypos3; TH1F *hsize1, *hsize2, *hsize3, *hsizex1, *hsizex2, *hsizex3, *hsizey1, *hsizey2, *hsizey3; TH1F *hrecHitsPerDet1, *hrecHitsPerDet2, *hrecHitsPerDet3; TH1F *hrecHitsPerLay1, *hrecHitsPerLay2, *hrecHitsPerLay3; TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3; TH1F *hdetr, *hdetz; // Forward endcaps TH1F *hdetrF, *hdetzF; TH1F *hdisk, *hblade, *hmodule, *hpanel, *hside; TH1F *hcharge1F, *hcharge2F; TH1F *hadcCharge1F, *hadcCharge2F; TH1F *hxpos1F, *hxpos2F, *hypos1F, *hypos2F; TH1F *hsize1F, *hsize2F, *hsizex1F, *hsizex2F, *hsizey1F, *hsizey2F; TH1F *hrecHitsPerDet1F, *hrecHitsPerDet2F; TH1F *hrecHitsPerLay1F, *hrecHitsPerLay2F; TH1F *hdetsPerLay1F, *hdetsPerLay2F; TH1F *hAlignErrorX1, *hAlignErrorX2, *hAlignErrorX3; TH1F *hAlignErrorX4, *hAlignErrorX5, *hAlignErrorX6, *hAlignErrorX7; TH1F *hAlignErrorY1, *hAlignErrorY2, *hAlignErrorY3; TH1F *hAlignErrorY4, *hAlignErrorY5, *hAlignErrorY6, *hAlignErrorY7; TH1F *hErrorX1, *hErrorX2, *hErrorX3, *hErrorX4, *hErrorX5, *hErrorX6, *hErrorX7; TH1F *hErrorY1, *hErrorY2, *hErrorY3, *hErrorY4, *hErrorY5, *hErrorY6, *hErrorY7; TProfile *hErrorXB, *hErrorXF, *hErrorYB, *hErrorYF; TProfile *hAErrorXB, *hAErrorXF, *hAErrorYB, *hAErrorYF; #endif }; using namespace std; //---------------------------------------------------------------------- ReadPixelRecHit::ReadPixelRecHit(edm::ParameterSet const &conf) : esTokenGeom_(esConsumes()), src_(conf.getParameter<edm::InputTag>("src")), siPixelRecHitsToken_(consumes<SiPixelRecHitCollection>(src_)), print(conf.getUntrackedParameter<bool>("Verbosity", false)) { usesResource(TFileService::kSharedResource); edm::LogPrint("ReadPixelRecHit") << " Verbosity " << print; } //---------------------------------------------------------------------- // Virtual destructor needed. ReadPixelRecHit::~ReadPixelRecHit() = default; //--------------------------------------------------------------------- // ------------ method called at the begining ------------ void ReadPixelRecHit::beginJob() { edm::LogPrint("ReadPixelRecHit") << "Initialize PixelRecHitTest "; #ifdef DO_HISTO // put here whatever you want to do at the beginning of the job // hFile = new TFile ( "histo.root", "RECREATE" ); // NEW way to use root (from 2.0.0?) edm::Service<TFileService> fs; // Histos go to a subdirectory "PixRecHits") //TFileDirectory subDir = fs->mkdir( "mySubDirectory" ); //TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" ); //h_pt = fs->make<TH1F>( "pt" , "p_{t}", 100, 0., ptMax_ ); hpixid = fs->make<TH1F>("hpixid", "Pix det id", 10, 0., 10.); hpixsubid = fs->make<TH1F>("hpixsubid", "Pix Barrel id", 10, 0., 10.); hlayerid = fs->make<TH1F>("hlayerid", "Pix layer id", 10, 0., 10.); hladder1id = fs->make<TH1F>("hladder1id", "Ladder L1 id", 50, 0., 50.); hladder2id = fs->make<TH1F>("hladder2id", "Ladder L2 id", 50, 0., 50.); hladder3id = fs->make<TH1F>("hladder3id", "Ladder L3 id", 50, 0., 50.); hz1id = fs->make<TH1F>("hz1id", "Z-index id L1", 10, 0., 10.); hz2id = fs->make<TH1F>("hz2id", "Z-index id L2", 10, 0., 10.); hz3id = fs->make<TH1F>("hz3id", "Z-index id L3", 10, 0., 10.); hrecHitsPerDet1 = fs->make<TH1F>("hrecHitsPerDet1", "RecHits per det l1", 200, -0.5, 199.5); hrecHitsPerDet2 = fs->make<TH1F>("hrecHitsPerDet2", "RecHits per det l2", 200, -0.5, 199.5); hrecHitsPerDet3 = fs->make<TH1F>("hrecHitsPerDet3", "RecHits per det l3", 200, -0.5, 199.5); hrecHitsPerLay1 = fs->make<TH1F>("hrecHitsPerLay1", "RecHits per layer l1", 2000, -0.5, 1999.5); hrecHitsPerLay2 = fs->make<TH1F>("hrecHitsPerLay2", "RecHits per layer l2", 2000, -0.5, 1999.5); hrecHitsPerLay3 = fs->make<TH1F>("hrecHitsPerLay3", "RecHits per layer l3", 2000, -0.5, 1999.5); hdetsPerLay1 = fs->make<TH1F>("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5); hdetsPerLay3 = fs->make<TH1F>("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5); hdetsPerLay2 = fs->make<TH1F>("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5); hcharge1 = fs->make<TH1F>("hcharge1", "Clu charge l1", 200, 0., 200.); //in ke hcharge2 = fs->make<TH1F>("hcharge2", "Clu charge l2", 200, 0., 200.); hcharge3 = fs->make<TH1F>("hcharge3", "Clu charge l3", 200, 0., 200.); hadcCharge1 = fs->make<TH1F>("hadcCharge1", "pix charge l1", 50, 0., 50.); //in ke hadcCharge2 = fs->make<TH1F>("hadcCharge2", "pix charge l2", 50, 0., 50.); //in ke hadcCharge3 = fs->make<TH1F>("hadcCharge3", "pix charge l3", 50, 0., 50.); //in ke hadcCharge1big = fs->make<TH1F>("hadcCharge1big", "big pix charge l1", 50, 0., 50.); //in ke hxpos1 = fs->make<TH1F>("hxpos1", "Layer 1 cols", 700, -3.5, 3.5); hxpos2 = fs->make<TH1F>("hxpos2", "Layer 2 cols", 700, -3.5, 3.5); hxpos3 = fs->make<TH1F>("hxpos3", "Layer 3 cols", 700, -3.5, 3.5); hypos1 = fs->make<TH1F>("hypos1", "Layer 1 rows", 200, -1., 1.); hypos2 = fs->make<TH1F>("hypos2", "Layer 2 rows", 200, -1., 1.); hypos3 = fs->make<TH1F>("hypos3", "layer 3 rows", 200, -1., 1.); hsize1 = fs->make<TH1F>("hsize1", "layer 1 clu size", 100, -0.5, 99.5); hsize2 = fs->make<TH1F>("hsize2", "layer 2 clu size", 100, -0.5, 99.5); hsize3 = fs->make<TH1F>("hsize3", "layer 3 clu size", 100, -0.5, 99.5); hAlignErrorX1 = fs->make<TH1F>("hAlignErrorX1", "Align error Layer 1 X", 100, 0.0, 100.); hAlignErrorY1 = fs->make<TH1F>("hAlignErrorY1", "Align error Layer 1 Y", 100, 0.0, 100.); hAlignErrorX2 = fs->make<TH1F>("hAlignErrorX2", "Align error Layer 2 X", 100, 0.0, 100.); hAlignErrorY2 = fs->make<TH1F>("hAlignErrorY2", "Align error Layer 2 Y", 100, 0.0, 100.); hAlignErrorX3 = fs->make<TH1F>("hAlignErrorX3", "Align error Layer 3 X", 100, 0.0, 100.); hAlignErrorY3 = fs->make<TH1F>("hAlignErrorY3", "Align error Layer 3 Y", 100, 0.0, 100.); hAlignErrorX4 = fs->make<TH1F>("hAlignErrorX4", "Align error Disk -1 X", 100, 0.0, 100.); hAlignErrorY4 = fs->make<TH1F>("hAlignErrorY4", "Align error Disk -1 Y", 100, 0.0, 100.); hAlignErrorX5 = fs->make<TH1F>("hAlignErrorX5", "Align error Disk -2 X", 100, 0.0, 100.); hAlignErrorY5 = fs->make<TH1F>("hAlignErrorY5", "Align error Disk -2 Y", 100, 0.0, 100.); hAlignErrorX6 = fs->make<TH1F>("hAlignErrorX6", "Align error Disk +1 X", 100, 0.0, 100.); hAlignErrorY6 = fs->make<TH1F>("hAlignErrorY6", "Align error Disk +1 Y", 100, 0.0, 100.); hAlignErrorX7 = fs->make<TH1F>("hAlignErrorX7", "Align error Disk +2 X", 100, 0.0, 100.); hAlignErrorY7 = fs->make<TH1F>("hAlignErrorY7", "Align error Disk +2 Y", 100, 0.0, 100.); hErrorX1 = fs->make<TH1F>("hErrorX1", "Error Layer 1 X", 100, 0.0, 100.); hErrorY1 = fs->make<TH1F>("hErrorY1", "Error Layer 1 Y", 100, 0.0, 100.); hErrorX2 = fs->make<TH1F>("hErrorX2", "Error Layer 2 X", 100, 0.0, 100.); hErrorY2 = fs->make<TH1F>("hErrorY2", "Error Layer 2 Y", 100, 0.0, 100.); hErrorX3 = fs->make<TH1F>("hErrorX3", "Error Layer 3 X", 100, 0.0, 100.); hErrorY3 = fs->make<TH1F>("hErrorY3", "Error Layer 3 Y", 100, 0.0, 100.); hErrorX4 = fs->make<TH1F>("hErrorX4", "Error Disk -1 X", 100, 0.0, 100.); hErrorY4 = fs->make<TH1F>("hErrorY4", "Error Disk -1 Y", 100, 0.0, 100.); hErrorX5 = fs->make<TH1F>("hErrorX5", "Error Disk -2 X", 100, 0.0, 100.); hErrorY5 = fs->make<TH1F>("hErrorY5", "Error Disk -2 Y", 100, 0.0, 100.); hErrorX6 = fs->make<TH1F>("hErrorX6", "Error Disk +1 X", 100, 0.0, 100.); hErrorY6 = fs->make<TH1F>("hErrorY6", "Error Disk +1 Y", 100, 0.0, 100.); hErrorX7 = fs->make<TH1F>("hErrorX7", "Error Disk +2 X", 100, 0.0, 100.); hErrorY7 = fs->make<TH1F>("hErrorY7", "Error Disk +2 Y", 100, 0.0, 100.); hsizex1 = fs->make<TH1F>("hsizex1", "lay1 clu size in x", 10, -0.5, 9.5); hsizex2 = fs->make<TH1F>("hsizex2", "lay2 clu size in x", 10, -0.5, 9.5); hsizex3 = fs->make<TH1F>("hsizex3", "lay3 clu size in x", 10, -0.5, 9.5); hsizey1 = fs->make<TH1F>("hsizey1", "lay1 clu size in y", 20, -0.5, 19.5); hsizey2 = fs->make<TH1F>("hsizey2", "lay2 clu size in y", 20, -0.5, 19.5); hsizey3 = fs->make<TH1F>("hsizey3", "lay3 clu size in y", 20, -0.5, 19.5); hdetr = fs->make<TH1F>("hdetr", "det r", 150, 0., 15.); hdetz = fs->make<TH1F>("hdetz", "det z", 520, -26., 26.); // Forward edcaps hdetrF = fs->make<TH1F>("hdetrF", "Fdet r", 150, 5., 20.); hdetzF = fs->make<TH1F>("hdetzF", "Fdet z", 600, -60., 60.); hdisk = fs->make<TH1F>("hdisk", "FPix disk id", 10, 0., 10.); hblade = fs->make<TH1F>("hblade", "FPix blade id", 30, 0., 30.); hmodule = fs->make<TH1F>("hmodule", "FPix plaq. id", 10, 0., 10.); hpanel = fs->make<TH1F>("hpanel", "FPix panel id", 10, 0., 10.); hside = fs->make<TH1F>("hside", "FPix size id", 10, 0., 10.); hcharge1F = fs->make<TH1F>("hcharge1F", "Clu charge 21", 200, 0., 200.); //in ke hcharge2F = fs->make<TH1F>("hcharge2F", "Clu charge 22", 200, 0., 200.); hxpos1F = fs->make<TH1F>("hxpos1F", "Disk 1 cols", 700, -3.5, 3.5); hxpos2F = fs->make<TH1F>("hxpos2F", "Disk 2 cols", 700, -3.5, 3.5); hypos1F = fs->make<TH1F>("hypos1F", "Disk 1 rows", 200, -1., 1.); hypos2F = fs->make<TH1F>("hypos2F", "Disk 2 rows", 200, -1., 1.); hsize1F = fs->make<TH1F>("hsize1F", "Disk 1 clu size", 100, -0.5, 99.5); hsize2F = fs->make<TH1F>("hsize2F", "Disk 2 clu size", 100, -0.5, 99.5); hsizex1F = fs->make<TH1F>("hsizex1F", "d1 clu size in x", 10, -0.5, 9.5); hsizex2F = fs->make<TH1F>("hsizex2F", "d2 clu size in x", 10, -0.5, 9.5); hsizey1F = fs->make<TH1F>("hsizey1F", "d1 clu size in y", 20, -0.5, 19.5); hsizey2F = fs->make<TH1F>("hsizey2F", "d2 clu size in y", 20, -0.5, 19.5); hadcCharge1F = fs->make<TH1F>("hadcCharge1F", "pix charge d1", 50, 0., 50.); //in ke hadcCharge2F = fs->make<TH1F>("hadcCharge2F", "pix charge d2", 50, 0., 50.); //in ke hrecHitsPerDet1F = fs->make<TH1F>("hrecHitsPerDet1F", "RecHits per det l1", 200, -0.5, 199.5); hrecHitsPerDet2F = fs->make<TH1F>("hrecHitsPerDet2F", "RecHits per det l2", 200, -0.5, 199.5); hrecHitsPerLay1F = fs->make<TH1F>("hrecHitsPerLay1F", "RecHits per layer l1", 2000, -0.5, 1999.5); hrecHitsPerLay2F = fs->make<TH1F>("hrecHitsPerLay2F", "RecHits per layer l2", 2000, -0.5, 1999.5); hdetsPerLay1F = fs->make<TH1F>("hdetsPerLay1F", "Full dets per layer l1", 161, -0.5, 160.5); hdetsPerLay2F = fs->make<TH1F>("hdetsPerLay2F", "Full dets per layer l2", 257, -0.5, 256.5); hErrorXB = fs->make<TProfile>("hErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.); hErrorXF = fs->make<TProfile>("hErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.); hErrorYB = fs->make<TProfile>("hErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.); hErrorYF = fs->make<TProfile>("hErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.); hAErrorXB = fs->make<TProfile>("hAErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.); hAErrorXF = fs->make<TProfile>("hAErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.); hAErrorYB = fs->make<TProfile>("hAErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.); hAErrorYF = fs->make<TProfile>("hAErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.); edm::LogPrint("ReadPixelRecHit") << " book histos "; #endif } //----------------------------------------------------------------------- void ReadPixelRecHit::endJob() { edm::LogPrint("ReadPixelRecHit") << " End PixelRecHitTest "; } //--------------------------------------------------------------------- // Functions that gets called by framework every event void ReadPixelRecHit::analyze(const edm::Event &e, const edm::EventSetup &es) { using namespace edm; const bool localPrint = false; //const bool localPrint = true; // Get event setup (to get global transformation) const TrackerGeometry &theTracker = es.getData(esTokenGeom_); edm::Handle<SiPixelRecHitCollection> recHitColl; e.getByToken(siPixelRecHitsToken_, recHitColl); if (print) edm::LogPrint("ReadPixelRecHit") << " FOUND " << (recHitColl.product())->dataSize() << " Pixel Hits"; SiPixelRecHitCollection::const_iterator recHitIdIterator = (recHitColl.product())->begin(); SiPixelRecHitCollection::const_iterator recHitIdIteratorEnd = (recHitColl.product())->end(); int numOfRecHits = 0; int numberOfDetUnits1 = 0; int numOfRecHitsPerDet1 = 0; int numOfRecHitsPerLay1 = 0; int numberOfDetUnits2 = 0; int numOfRecHitsPerDet2 = 0; int numOfRecHitsPerLay2 = 0; int numberOfDetUnits3 = 0; int numOfRecHitsPerDet3 = 0; int numOfRecHitsPerLay3 = 0; int numberOfDetUnits1F = 0; int numOfRecHitsPerDet1F = 0; int numOfRecHitsPerLay1F = 0; int numberOfDetUnits2F = 0; int numOfRecHitsPerDet2F = 0; int numOfRecHitsPerLay2F = 0; // Loop over Detector IDs for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) { SiPixelRecHitCollection::DetSet detset = *recHitIdIterator; DetId detId = DetId(detset.detId()); // Get the Detid object unsigned int detType = detId.det(); // det type, tracker=1 unsigned int subid = detId.subdetId(); //subdetector type, barrel=1, fpix=2 if (detType != 1) edm::LogPrint("ReadPixelRecHit") << " wrong det id " << detType; //look only at tracker if (!((subid == 1) || (subid == 2))) edm::LogPrint("ReadPixelRecHit") << " wrong sub det id " << subid; // look only at bpix&fpix if (print) edm::LogPrint("ReadPixelRecHit") << " Det ID " << detId.rawId(); // Get rechits if (detset.empty()) continue; // Get the geometrical information for this det const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId)); double detZ = theGeomDet->surface().position().z(); double detR = theGeomDet->surface().position().perp(); //const BoundPlane& plane = theGeomDet->surface(); //for transf. unused //double detThick = theGeomDet->specificSurface().bounds().thickness(); unused const PixelTopology *topol = &(theGeomDet->specificTopology()); //int cols = theGeomDet->specificTopology().ncolumns(); UNUSED //int rows = theGeomDet->specificTopology().nrows(); float alignErrorX = 0., alignErrorY = 0.; LocalError lape = theGeomDet->localAlignmentError(); //edm::LogPrint("ReadPixelRecHit")<<lape.valid(); if (lape.valid()) { if (lape.xx() > 0.) alignErrorX = sqrt(lape.xx()) * 1E4; //float tmp12= sqrt(lape.xy())*1E4; if (lape.yy() > 0.) alignErrorY = sqrt(lape.yy()) * 1E4; if (print) edm::LogPrint("ReadPixelRecHit") << " Alignment errors " << alignErrorX << " " << alignErrorY; //edm::LogPrint("ReadPixelRecHit")<<" Alignment errors "<<alignErrorX<<" "<<alignErrorY; } unsigned int layer = 0, disk = 0, ladder = 0, zindex = 0, blade = 0, panel = 0, side = 0; if (subid == 1) { // Subdet it, pix barrel=1 PXBDetId pdetId = PXBDetId(detId); //unsigned int detTypeP=pdetId.det(); unused //unsigned int subidP=pdetId.subdetId(); unused // Barell layer = 1,2,3 layer = pdetId.layer(); // Barrel ladder id 1-20,32,44. ladder = pdetId.ladder(); // Barrel Z-index=1,8 zindex = pdetId.module(); if (zindex < 5) side = 1; else side = 2; if (localPrint) edm::LogPrint("ReadPixelRecHit") << " Layer " << layer << " ladder " << ladder << " z " << zindex; //<<pdetId.rawId()<<" "<<pdetId.null()<<detTypeP<<" "<<subidP<<" "; #ifdef DO_HISTO hdetr->Fill(detR); hdetz->Fill(detZ); //hcolsB->Fill(float(cols)); //hyposB->Fill(float(rows)); hlayerid->Fill(float(layer)); #endif } else { // FPIX ------------------------------------- // test cols & rows //edm::LogPrint("ReadPixelRecHit")<<" det z/r "<<detZ<<" "<<detR<<" "<<detThick<<" " // <<cols<<" "<<rows; PXFDetId pdetId = PXFDetId(detId.rawId()); disk = pdetId.disk(); //1,2,3 blade = pdetId.blade(); //1-24 side = pdetId.side(); //size=1 for -z, 2 for +z panel = pdetId.panel(); //panel=1,2 unsigned int module = pdetId.module(); // plaquette if (print) edm::LogPrint("ReadPixelRecHit") << " forward hit " << side << " " << disk << " " << blade << " " << panel << " " << module; #ifdef DO_HISTO hdetrF->Fill(detR); hdetzF->Fill(detZ); hdisk->Fill(float(disk)); hblade->Fill(float(blade)); hmodule->Fill(float(module)); hpanel->Fill(float(panel)); hside->Fill(float(side)); #endif } // end BPix FPix if #ifdef DO_HISTO if (layer == 1) { hladder1id->Fill(float(ladder)); hz1id->Fill(float(zindex)); hAlignErrorX1->Fill(alignErrorX); hAlignErrorY1->Fill(alignErrorY); hAErrorXB->Fill(float(ladder + (110 * (side - 1))), alignErrorX); hAErrorYB->Fill(float(ladder + (110 * (side - 1))), alignErrorY); ++numberOfDetUnits1; numOfRecHitsPerDet1 = 0; } else if (layer == 2) { hladder2id->Fill(float(ladder)); hz2id->Fill(float(zindex)); hAlignErrorX2->Fill(alignErrorX); hAlignErrorY2->Fill(alignErrorY); hAErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorX); hAErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorY); ++numberOfDetUnits2; numOfRecHitsPerDet2 = 0; } else if (layer == 3) { hladder3id->Fill(float(ladder)); hz3id->Fill(float(zindex)); hAlignErrorX3->Fill(alignErrorX); hAlignErrorY3->Fill(alignErrorY); hAErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorX); hAErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorY); ++numberOfDetUnits3; numOfRecHitsPerDet3 = 0; } else if (disk == 1) { ++numberOfDetUnits1F; numOfRecHitsPerDet1F = 0; if (side == 1) { hAlignErrorX4->Fill(alignErrorX); hAlignErrorY4->Fill(alignErrorY); hAErrorXF->Fill(float(blade + 25), alignErrorX); hAErrorYF->Fill(float(blade + 25), alignErrorY); } else { hAlignErrorX6->Fill(alignErrorX); hAlignErrorY6->Fill(alignErrorY); hAErrorXF->Fill(float(blade + 50), alignErrorX); hAErrorYF->Fill(float(blade + 50), alignErrorY); } } else if (disk == 2) { ++numberOfDetUnits2F; numOfRecHitsPerDet2F = 0; if (side == 1) { hAlignErrorX5->Fill(alignErrorX); hAlignErrorY5->Fill(alignErrorY); hAErrorXF->Fill(float(blade), alignErrorX); hAErrorYF->Fill(float(blade), alignErrorY); } else { hAlignErrorX7->Fill(alignErrorX); hAlignErrorY7->Fill(alignErrorY); hAErrorXF->Fill(float(blade + 75), alignErrorX); hAErrorYF->Fill(float(blade + 75), alignErrorY); } } #endif //----Loop over rechits for this detId SiPixelRecHitCollection::DetSet::const_iterator pixeliter = detset.begin(); SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end(); for (; pixeliter != rechitRangeIteratorEnd; ++pixeliter) { //loop on the rechit numOfRecHits++; // RecHit local position is now transient, // one needs to run tracking to get position OR rerun localreco LocalPoint lp = pixeliter->localPosition(); LocalError le = pixeliter->localPositionError(); float xRecHit = lp.x(); float yRecHit = lp.y(); float xerror = sqrt(le.xx()) * 1E4; float yerror = sqrt(le.yy()) * 1E4; if (print) edm::LogPrint("ReadPixelRecHit") << " RecHit: " << numOfRecHits << " x/y " << xRecHit << " " << yRecHit << " errors x/y " << xerror << " " << yerror; //MeasurementPoint mp = topol->measurementPosition(xRecHit,yRecHit); //GlobalPoint GP = PixGeom->surface().toGlobal(Local3DPoint(lp)); // Get cluster edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = pixeliter->cluster(); float ch = (clust->charge()) / 1000.; // convert ke to electrons int size = clust->size(); int sizeX = clust->sizeX(); int sizeY = clust->sizeY(); float xClu = clust->x(); float yClu = clust->y(); int maxPixelCol = clust->maxPixelCol(); int maxPixelRow = clust->maxPixelRow(); int minPixelCol = clust->minPixelCol(); int minPixelRow = clust->minPixelRow(); // unsigned int geoId = clust->geographicalId(); // alsways 0 // edge method moved to topologu class bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow)); bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol)); if (print) edm::LogPrint("ReadPixelRecHit") << "Clu: charge " << ch << " size " << size << " size x/y " << sizeX << " " << sizeY << " meas. " << xClu << " " << yClu << " edge " << edgeHitX << " " << edgeHitY; if (print) edm::LogPrint("ReadPixelRecHit") << " pixels:"; // Get the pixels in the Cluster const vector<SiPixelCluster::Pixel> &pixelsVec = clust->pixels(); //if(localPrint) edm::LogPrint("ReadPixelRecHit")<<" Pixels in this cluster "; map<unsigned int, float, less<unsigned int> > chanMap; // Channel map // Look at pixels in this cluster. ADC is calibrated, in electrons for (unsigned int i = 0; i < pixelsVec.size(); ++i) { float pixx = pixelsVec[i].x; // index as a float so = i+0.5 float pixy = pixelsVec[i].y; float adc = ((pixelsVec[i].adc) / 1000); // in kelec. // OLD way //int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy)); bool bigInX = topol->isItBigPixelInX(int(pixx)); bool bigInY = topol->isItBigPixelInY(int(pixy)); bool edgeInX = topol->isItEdgePixelInX(int(pixx)); bool edgeInY = topol->isItEdgePixelInY(int(pixy)); if (print) edm::LogPrint("ReadPixelRecHit") << i << " index " << pixx << " " << pixy << " adc " << adc << " edge " << edgeInX << " " << edgeInY << " big " << bigInX << " " << bigInY; //if(print && sizeX==1 && bigInX) //edm::LogPrint("ReadPixelRecHit")<<" single big x "<<xClu<<" "<<pixx<<" "; //if(print && sizeY==1 && bigInY) //edm::LogPrint("ReadPixelRecHit")<<" single big y "<<yClu<<" "<<pixy<<" "; #ifdef DO_HISTO if (layer == 1) { hadcCharge1->Fill(adc); //if(bigInX || bigInY) hadcCharge1big->Fill(adc); } else if (layer == 2) { hadcCharge2->Fill(adc); } else if (layer == 3) { hadcCharge3->Fill(adc); } else if (disk == 1) { hadcCharge1F->Fill(adc); } else if (disk == 2) { hadcCharge2F->Fill(adc); } #endif } // End pixel loop #ifdef DO_HISTO if (layer == 1) { hcharge1->Fill(ch); hxpos1->Fill(yRecHit); hypos1->Fill(xRecHit); hsize1->Fill(float(size)); hsizex1->Fill(float(sizeX)); hsizey1->Fill(float(sizeY)); numOfRecHitsPerDet1++; numOfRecHitsPerLay1++; hErrorX1->Fill(xerror); hErrorY1->Fill(yerror); hErrorXB->Fill(float(ladder + (110 * (side - 1))), xerror); hErrorYB->Fill(float(ladder + (110 * (side - 1))), yerror); } else if (layer == 2) { // layer 2 hcharge2->Fill(ch); hxpos2->Fill(yRecHit); hypos2->Fill(xRecHit); hsize2->Fill(float(size)); hsizex2->Fill(float(sizeX)); hsizey2->Fill(float(sizeY)); numOfRecHitsPerDet2++; numOfRecHitsPerLay2++; hErrorX2->Fill(xerror); hErrorY2->Fill(yerror); hErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), xerror); hErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), yerror); } else if (layer == 3) { // Layer 3 hcharge3->Fill(ch); hxpos3->Fill(yRecHit); hypos3->Fill(xRecHit); hsize3->Fill(float(size)); hsizex3->Fill(float(sizeX)); hsizey3->Fill(float(sizeY)); numOfRecHitsPerDet3++; numOfRecHitsPerLay3++; hErrorX3->Fill(xerror); hErrorY3->Fill(yerror); hErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), xerror); hErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), yerror); } else if (disk == 1) { hcharge1F->Fill(ch); hxpos1F->Fill(yRecHit); hypos1F->Fill(xRecHit); hsize1F->Fill(float(size)); hsizex1F->Fill(float(sizeX)); hsizey1F->Fill(float(sizeY)); if (side == 1) { // -z hErrorX4->Fill(xerror); hErrorY4->Fill(yerror); hErrorXF->Fill(float(blade + 25), xerror); hErrorYF->Fill(float(blade + 25), yerror); } else { // +z hErrorX6->Fill(xerror); hErrorY6->Fill(yerror); hErrorXF->Fill(float(blade + 50), xerror); hErrorYF->Fill(float(blade + 50), yerror); } numOfRecHitsPerDet1F++; numOfRecHitsPerLay1F++; } else if (disk == 2) { // disk 2 hcharge2F->Fill(ch); hxpos2F->Fill(yRecHit); hypos2F->Fill(xRecHit); hsize2F->Fill(float(size)); hsizex2F->Fill(float(sizeX)); hsizey2F->Fill(float(sizeY)); numOfRecHitsPerDet2F++; numOfRecHitsPerLay2F++; if (side == 1) { // -z hErrorX5->Fill(xerror); hErrorY5->Fill(yerror); hErrorXF->Fill(float(blade), xerror); hErrorYF->Fill(float(blade), yerror); } else { // +z hErrorX7->Fill(xerror); hErrorY7->Fill(yerror); hErrorXF->Fill(float(blade + 75), xerror); hErrorYF->Fill(float(blade + 75), yerror); } } #endif } // End RecHit loop #ifdef DO_HISTO if (layer == 1) hrecHitsPerDet1->Fill(float(numOfRecHitsPerDet1)); else if (layer == 2) hrecHitsPerDet2->Fill(float(numOfRecHitsPerDet2)); else if (layer == 3) hrecHitsPerDet3->Fill(float(numOfRecHitsPerDet3)); else if (disk == 1) hrecHitsPerDet1F->Fill(float(numOfRecHitsPerDet1F)); else if (disk == 2) hrecHitsPerDet2F->Fill(float(numOfRecHitsPerDet2F)); #endif } // End Det loop #ifdef DO_HISTO if (print) edm::LogPrint("ReadPixelRecHit") << " Rechits per layer " << numOfRecHitsPerLay1 << " " << numOfRecHitsPerLay2 << " " << numOfRecHitsPerLay3 << " " << numOfRecHitsPerLay1F << " " << numOfRecHitsPerLay2F; hrecHitsPerLay1->Fill(float(numOfRecHitsPerLay1)); hrecHitsPerLay2->Fill(float(numOfRecHitsPerLay2)); hrecHitsPerLay3->Fill(float(numOfRecHitsPerLay3)); hdetsPerLay1->Fill(float(numberOfDetUnits1)); hdetsPerLay2->Fill(float(numberOfDetUnits2)); hdetsPerLay3->Fill(float(numberOfDetUnits3)); hrecHitsPerLay1F->Fill(float(numOfRecHitsPerLay1F)); hrecHitsPerLay2F->Fill(float(numOfRecHitsPerLay2F)); hdetsPerLay1F->Fill(float(numberOfDetUnits1F)); hdetsPerLay2F->Fill(float(numberOfDetUnits2F)); hdetsPerLay3->Fill(float(numberOfDetUnits3)); #endif } //define this as a plug-in DEFINE_FWK_MODULE(ReadPixelRecHit);