Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:00:41 

0001 //--------------------------------------------
0002 // File: ReadPixelRecHit.cc
0003 // Description:  see ReadPixelRecHit.h
0004 // Author:  J.Sheav (JHU)
0005 //          11/8/06: New loop over rechits and InputTag, V.Chiochia
0006 // Creation Date:  OGU Aug. 1 2005 Initial version.
0007 // Add occupancy histograms. D.Kotlinski. 10/06
0008 //--------------------------------------------
0009 
0010 // system includes
0011 #include <memory>
0012 #include <string>
0013 #include <iostream>
0014 
0015 // user includes
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0022 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0023 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0024 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0027 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0030 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0032 #include "FWCore/ServiceRegistry/interface/Service.h"
0033 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0034 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "DataFormats/Common/interface/Handle.h"
0037 #include "FWCore/Framework/interface/EventSetup.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 
0040 #define DO_HISTO
0041 #ifdef DO_HISTO
0042 // For ROOT
0043 #include <TROOT.h>
0044 #include <TChain.h>
0045 #include <TFile.h>
0046 #include <TF1.h>
0047 #include <TH2F.h>
0048 #include <TH1F.h>
0049 #include <TProfile.h>
0050 #endif
0051 
0052 class ReadPixelRecHit : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0053 public:
0054   explicit ReadPixelRecHit(const edm::ParameterSet &conf);
0055   virtual ~ReadPixelRecHit();
0056 
0057   virtual void analyze(const edm::Event &e, const edm::EventSetup &c);
0058   virtual void beginJob();
0059   virtual void endJob();
0060 
0061 private:
0062   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> esTokenGeom_;
0063   edm::InputTag src_;
0064   const edm::EDGetTokenT<SiPixelRecHitCollection> siPixelRecHitsToken_;
0065   bool print;
0066 
0067 #ifdef DO_HISTO
0068   TFile *hFile;
0069   TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id;
0070   TH1F *hcharge1, *hcharge2, *hcharge3;
0071   TH1F *hadcCharge1, *hadcCharge2, *hadcCharge3, *hadcCharge1big;
0072   TH1F *hxpos1, *hxpos2, *hxpos3, *hypos1, *hypos2, *hypos3;
0073   TH1F *hsize1, *hsize2, *hsize3, *hsizex1, *hsizex2, *hsizex3, *hsizey1, *hsizey2, *hsizey3;
0074 
0075   TH1F *hrecHitsPerDet1, *hrecHitsPerDet2, *hrecHitsPerDet3;
0076   TH1F *hrecHitsPerLay1, *hrecHitsPerLay2, *hrecHitsPerLay3;
0077   TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3;
0078 
0079   TH1F *hdetr, *hdetz;
0080 
0081   // Forward endcaps
0082   TH1F *hdetrF, *hdetzF;
0083   TH1F *hdisk, *hblade, *hmodule, *hpanel, *hside;
0084   TH1F *hcharge1F, *hcharge2F;
0085   TH1F *hadcCharge1F, *hadcCharge2F;
0086   TH1F *hxpos1F, *hxpos2F, *hypos1F, *hypos2F;
0087   TH1F *hsize1F, *hsize2F, *hsizex1F, *hsizex2F, *hsizey1F, *hsizey2F;
0088   TH1F *hrecHitsPerDet1F, *hrecHitsPerDet2F;
0089   TH1F *hrecHitsPerLay1F, *hrecHitsPerLay2F;
0090   TH1F *hdetsPerLay1F, *hdetsPerLay2F;
0091 
0092   TH1F *hAlignErrorX1, *hAlignErrorX2, *hAlignErrorX3;
0093   TH1F *hAlignErrorX4, *hAlignErrorX5, *hAlignErrorX6, *hAlignErrorX7;
0094   TH1F *hAlignErrorY1, *hAlignErrorY2, *hAlignErrorY3;
0095   TH1F *hAlignErrorY4, *hAlignErrorY5, *hAlignErrorY6, *hAlignErrorY7;
0096   TH1F *hErrorX1, *hErrorX2, *hErrorX3, *hErrorX4, *hErrorX5, *hErrorX6, *hErrorX7;
0097   TH1F *hErrorY1, *hErrorY2, *hErrorY3, *hErrorY4, *hErrorY5, *hErrorY6, *hErrorY7;
0098   TProfile *hErrorXB, *hErrorXF, *hErrorYB, *hErrorYF;
0099   TProfile *hAErrorXB, *hAErrorXF, *hAErrorYB, *hAErrorYF;
0100 #endif
0101 };
0102 
0103 using namespace std;
0104 
0105 //----------------------------------------------------------------------
0106 ReadPixelRecHit::ReadPixelRecHit(edm::ParameterSet const &conf)
0107     : esTokenGeom_(esConsumes()),
0108       src_(conf.getParameter<edm::InputTag>("src")),
0109       siPixelRecHitsToken_(consumes<SiPixelRecHitCollection>(src_)),
0110       print(conf.getUntrackedParameter<bool>("Verbosity", false)) {
0111   usesResource(TFileService::kSharedResource);
0112   edm::LogPrint("ReadPixelRecHit") << " Verbosity " << print;
0113 }
0114 //----------------------------------------------------------------------
0115 // Virtual destructor needed.
0116 ReadPixelRecHit::~ReadPixelRecHit() = default;
0117 //---------------------------------------------------------------------
0118 // ------------ method called at the begining   ------------
0119 void ReadPixelRecHit::beginJob() {
0120   edm::LogPrint("ReadPixelRecHit") << "Initialize PixelRecHitTest ";
0121 
0122 #ifdef DO_HISTO
0123   // put here whatever you want to do at the beginning of the job
0124   // hFile = new TFile ( "histo.root", "RECREATE" );
0125 
0126   // NEW way to use root (from 2.0.0?)
0127   edm::Service<TFileService> fs;
0128 
0129   // Histos go to a subdirectory "PixRecHits")
0130   //TFileDirectory subDir = fs->mkdir( "mySubDirectory" );
0131   //TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
0132   //h_pt    = fs->make<TH1F>( "pt"  , "p_{t}", 100,  0., ptMax_ );
0133 
0134   hpixid = fs->make<TH1F>("hpixid", "Pix det id", 10, 0., 10.);
0135   hpixsubid = fs->make<TH1F>("hpixsubid", "Pix Barrel id", 10, 0., 10.);
0136   hlayerid = fs->make<TH1F>("hlayerid", "Pix layer id", 10, 0., 10.);
0137   hladder1id = fs->make<TH1F>("hladder1id", "Ladder L1 id", 50, 0., 50.);
0138   hladder2id = fs->make<TH1F>("hladder2id", "Ladder L2 id", 50, 0., 50.);
0139   hladder3id = fs->make<TH1F>("hladder3id", "Ladder L3 id", 50, 0., 50.);
0140   hz1id = fs->make<TH1F>("hz1id", "Z-index id L1", 10, 0., 10.);
0141   hz2id = fs->make<TH1F>("hz2id", "Z-index id L2", 10, 0., 10.);
0142   hz3id = fs->make<TH1F>("hz3id", "Z-index id L3", 10, 0., 10.);
0143 
0144   hrecHitsPerDet1 = fs->make<TH1F>("hrecHitsPerDet1", "RecHits per det l1", 200, -0.5, 199.5);
0145   hrecHitsPerDet2 = fs->make<TH1F>("hrecHitsPerDet2", "RecHits per det l2", 200, -0.5, 199.5);
0146   hrecHitsPerDet3 = fs->make<TH1F>("hrecHitsPerDet3", "RecHits per det l3", 200, -0.5, 199.5);
0147   hrecHitsPerLay1 = fs->make<TH1F>("hrecHitsPerLay1", "RecHits per layer l1", 2000, -0.5, 1999.5);
0148   hrecHitsPerLay2 = fs->make<TH1F>("hrecHitsPerLay2", "RecHits per layer l2", 2000, -0.5, 1999.5);
0149 
0150   hrecHitsPerLay3 = fs->make<TH1F>("hrecHitsPerLay3", "RecHits per layer l3", 2000, -0.5, 1999.5);
0151   hdetsPerLay1 = fs->make<TH1F>("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5);
0152   hdetsPerLay3 = fs->make<TH1F>("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5);
0153   hdetsPerLay2 = fs->make<TH1F>("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5);
0154 
0155   hcharge1 = fs->make<TH1F>("hcharge1", "Clu charge l1", 200, 0., 200.);  //in ke
0156   hcharge2 = fs->make<TH1F>("hcharge2", "Clu charge l2", 200, 0., 200.);
0157   hcharge3 = fs->make<TH1F>("hcharge3", "Clu charge l3", 200, 0., 200.);
0158   hadcCharge1 = fs->make<TH1F>("hadcCharge1", "pix charge l1", 50, 0., 50.);            //in ke
0159   hadcCharge2 = fs->make<TH1F>("hadcCharge2", "pix charge l2", 50, 0., 50.);            //in ke
0160   hadcCharge3 = fs->make<TH1F>("hadcCharge3", "pix charge l3", 50, 0., 50.);            //in ke
0161   hadcCharge1big = fs->make<TH1F>("hadcCharge1big", "big pix charge l1", 50, 0., 50.);  //in ke
0162 
0163   hxpos1 = fs->make<TH1F>("hxpos1", "Layer 1 cols", 700, -3.5, 3.5);
0164   hxpos2 = fs->make<TH1F>("hxpos2", "Layer 2 cols", 700, -3.5, 3.5);
0165   hxpos3 = fs->make<TH1F>("hxpos3", "Layer 3 cols", 700, -3.5, 3.5);
0166 
0167   hypos1 = fs->make<TH1F>("hypos1", "Layer 1 rows", 200, -1., 1.);
0168   hypos2 = fs->make<TH1F>("hypos2", "Layer 2 rows", 200, -1., 1.);
0169   hypos3 = fs->make<TH1F>("hypos3", "layer 3 rows", 200, -1., 1.);
0170 
0171   hsize1 = fs->make<TH1F>("hsize1", "layer 1 clu size", 100, -0.5, 99.5);
0172   hsize2 = fs->make<TH1F>("hsize2", "layer 2 clu size", 100, -0.5, 99.5);
0173   hsize3 = fs->make<TH1F>("hsize3", "layer 3 clu size", 100, -0.5, 99.5);
0174 
0175   hAlignErrorX1 = fs->make<TH1F>("hAlignErrorX1", "Align error Layer 1 X", 100, 0.0, 100.);
0176   hAlignErrorY1 = fs->make<TH1F>("hAlignErrorY1", "Align error Layer 1 Y", 100, 0.0, 100.);
0177   hAlignErrorX2 = fs->make<TH1F>("hAlignErrorX2", "Align error Layer 2 X", 100, 0.0, 100.);
0178   hAlignErrorY2 = fs->make<TH1F>("hAlignErrorY2", "Align error Layer 2 Y", 100, 0.0, 100.);
0179   hAlignErrorX3 = fs->make<TH1F>("hAlignErrorX3", "Align error Layer 3 X", 100, 0.0, 100.);
0180   hAlignErrorY3 = fs->make<TH1F>("hAlignErrorY3", "Align error Layer 3 Y", 100, 0.0, 100.);
0181   hAlignErrorX4 = fs->make<TH1F>("hAlignErrorX4", "Align error Disk -1 X", 100, 0.0, 100.);
0182   hAlignErrorY4 = fs->make<TH1F>("hAlignErrorY4", "Align error Disk -1 Y", 100, 0.0, 100.);
0183   hAlignErrorX5 = fs->make<TH1F>("hAlignErrorX5", "Align error Disk -2 X", 100, 0.0, 100.);
0184   hAlignErrorY5 = fs->make<TH1F>("hAlignErrorY5", "Align error Disk -2 Y", 100, 0.0, 100.);
0185   hAlignErrorX6 = fs->make<TH1F>("hAlignErrorX6", "Align error Disk +1 X", 100, 0.0, 100.);
0186   hAlignErrorY6 = fs->make<TH1F>("hAlignErrorY6", "Align error Disk +1 Y", 100, 0.0, 100.);
0187   hAlignErrorX7 = fs->make<TH1F>("hAlignErrorX7", "Align error Disk +2 X", 100, 0.0, 100.);
0188   hAlignErrorY7 = fs->make<TH1F>("hAlignErrorY7", "Align error Disk +2 Y", 100, 0.0, 100.);
0189 
0190   hErrorX1 = fs->make<TH1F>("hErrorX1", "Error Layer 1 X", 100, 0.0, 100.);
0191   hErrorY1 = fs->make<TH1F>("hErrorY1", "Error Layer 1 Y", 100, 0.0, 100.);
0192   hErrorX2 = fs->make<TH1F>("hErrorX2", "Error Layer 2 X", 100, 0.0, 100.);
0193   hErrorY2 = fs->make<TH1F>("hErrorY2", "Error Layer 2 Y", 100, 0.0, 100.);
0194   hErrorX3 = fs->make<TH1F>("hErrorX3", "Error Layer 3 X", 100, 0.0, 100.);
0195   hErrorY3 = fs->make<TH1F>("hErrorY3", "Error Layer 3 Y", 100, 0.0, 100.);
0196   hErrorX4 = fs->make<TH1F>("hErrorX4", "Error Disk -1 X", 100, 0.0, 100.);
0197   hErrorY4 = fs->make<TH1F>("hErrorY4", "Error Disk -1 Y", 100, 0.0, 100.);
0198   hErrorX5 = fs->make<TH1F>("hErrorX5", "Error Disk -2 X", 100, 0.0, 100.);
0199   hErrorY5 = fs->make<TH1F>("hErrorY5", "Error Disk -2 Y", 100, 0.0, 100.);
0200   hErrorX6 = fs->make<TH1F>("hErrorX6", "Error Disk +1 X", 100, 0.0, 100.);
0201   hErrorY6 = fs->make<TH1F>("hErrorY6", "Error Disk +1 Y", 100, 0.0, 100.);
0202   hErrorX7 = fs->make<TH1F>("hErrorX7", "Error Disk +2 X", 100, 0.0, 100.);
0203   hErrorY7 = fs->make<TH1F>("hErrorY7", "Error Disk +2 Y", 100, 0.0, 100.);
0204 
0205   hsizex1 = fs->make<TH1F>("hsizex1", "lay1 clu size in x", 10, -0.5, 9.5);
0206   hsizex2 = fs->make<TH1F>("hsizex2", "lay2 clu size in x", 10, -0.5, 9.5);
0207   hsizex3 = fs->make<TH1F>("hsizex3", "lay3 clu size in x", 10, -0.5, 9.5);
0208   hsizey1 = fs->make<TH1F>("hsizey1", "lay1 clu size in y", 20, -0.5, 19.5);
0209   hsizey2 = fs->make<TH1F>("hsizey2", "lay2 clu size in y", 20, -0.5, 19.5);
0210   hsizey3 = fs->make<TH1F>("hsizey3", "lay3 clu size in y", 20, -0.5, 19.5);
0211 
0212   hdetr = fs->make<TH1F>("hdetr", "det r", 150, 0., 15.);
0213   hdetz = fs->make<TH1F>("hdetz", "det z", 520, -26., 26.);
0214 
0215   // Forward edcaps
0216   hdetrF = fs->make<TH1F>("hdetrF", "Fdet r", 150, 5., 20.);
0217   hdetzF = fs->make<TH1F>("hdetzF", "Fdet z", 600, -60., 60.);
0218 
0219   hdisk = fs->make<TH1F>("hdisk", "FPix disk id", 10, 0., 10.);
0220   hblade = fs->make<TH1F>("hblade", "FPix blade id", 30, 0., 30.);
0221   hmodule = fs->make<TH1F>("hmodule", "FPix plaq. id", 10, 0., 10.);
0222   hpanel = fs->make<TH1F>("hpanel", "FPix panel id", 10, 0., 10.);
0223   hside = fs->make<TH1F>("hside", "FPix size id", 10, 0., 10.);
0224 
0225   hcharge1F = fs->make<TH1F>("hcharge1F", "Clu charge 21", 200, 0., 200.);  //in ke
0226   hcharge2F = fs->make<TH1F>("hcharge2F", "Clu charge 22", 200, 0., 200.);
0227   hxpos1F = fs->make<TH1F>("hxpos1F", "Disk 1 cols", 700, -3.5, 3.5);
0228   hxpos2F = fs->make<TH1F>("hxpos2F", "Disk 2 cols", 700, -3.5, 3.5);
0229   hypos1F = fs->make<TH1F>("hypos1F", "Disk 1 rows", 200, -1., 1.);
0230   hypos2F = fs->make<TH1F>("hypos2F", "Disk 2 rows", 200, -1., 1.);
0231   hsize1F = fs->make<TH1F>("hsize1F", "Disk 1 clu size", 100, -0.5, 99.5);
0232   hsize2F = fs->make<TH1F>("hsize2F", "Disk 2 clu size", 100, -0.5, 99.5);
0233   hsizex1F = fs->make<TH1F>("hsizex1F", "d1 clu size in x", 10, -0.5, 9.5);
0234   hsizex2F = fs->make<TH1F>("hsizex2F", "d2 clu size in x", 10, -0.5, 9.5);
0235   hsizey1F = fs->make<TH1F>("hsizey1F", "d1 clu size in y", 20, -0.5, 19.5);
0236   hsizey2F = fs->make<TH1F>("hsizey2F", "d2 clu size in y", 20, -0.5, 19.5);
0237   hadcCharge1F = fs->make<TH1F>("hadcCharge1F", "pix charge d1", 50, 0., 50.);  //in ke
0238   hadcCharge2F = fs->make<TH1F>("hadcCharge2F", "pix charge d2", 50, 0., 50.);  //in ke
0239 
0240   hrecHitsPerDet1F = fs->make<TH1F>("hrecHitsPerDet1F", "RecHits per det l1", 200, -0.5, 199.5);
0241   hrecHitsPerDet2F = fs->make<TH1F>("hrecHitsPerDet2F", "RecHits per det l2", 200, -0.5, 199.5);
0242   hrecHitsPerLay1F = fs->make<TH1F>("hrecHitsPerLay1F", "RecHits per layer l1", 2000, -0.5, 1999.5);
0243   hrecHitsPerLay2F = fs->make<TH1F>("hrecHitsPerLay2F", "RecHits per layer l2", 2000, -0.5, 1999.5);
0244   hdetsPerLay1F = fs->make<TH1F>("hdetsPerLay1F", "Full dets per layer l1", 161, -0.5, 160.5);
0245   hdetsPerLay2F = fs->make<TH1F>("hdetsPerLay2F", "Full dets per layer l2", 257, -0.5, 256.5);
0246 
0247   hErrorXB = fs->make<TProfile>("hErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.);
0248   hErrorXF = fs->make<TProfile>("hErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.);
0249   hErrorYB = fs->make<TProfile>("hErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.);
0250   hErrorYF = fs->make<TProfile>("hErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.);
0251 
0252   hAErrorXB = fs->make<TProfile>("hAErrorXB", "bpix x errors per ladder", 220, 0., 220., 0.0, 1000.);
0253   hAErrorXF = fs->make<TProfile>("hAErrorXF", "fpix x errors per ladder", 100, 0., 100., 0.0, 1000.);
0254   hAErrorYB = fs->make<TProfile>("hAErrorYB", "bpix y errors per ladder", 220, 0., 220., 0.0, 1000.);
0255   hAErrorYF = fs->make<TProfile>("hAErrorYF", "fpix y errors per ladder", 100, 0., 100., 0.0, 1000.);
0256 
0257   edm::LogPrint("ReadPixelRecHit") << " book histos ";
0258 
0259 #endif
0260 }
0261 //-----------------------------------------------------------------------
0262 void ReadPixelRecHit::endJob() { edm::LogPrint("ReadPixelRecHit") << " End PixelRecHitTest "; }
0263 //---------------------------------------------------------------------
0264 // Functions that gets called by framework every event
0265 void ReadPixelRecHit::analyze(const edm::Event &e, const edm::EventSetup &es) {
0266   using namespace edm;
0267   const bool localPrint = false;
0268   //const bool localPrint = true;
0269 
0270   // Get event setup (to get global transformation)
0271   const TrackerGeometry &theTracker = es.getData(esTokenGeom_);
0272 
0273   edm::Handle<SiPixelRecHitCollection> recHitColl;
0274   e.getByToken(siPixelRecHitsToken_, recHitColl);
0275 
0276   if (print)
0277     edm::LogPrint("ReadPixelRecHit") << " FOUND " << (recHitColl.product())->dataSize() << " Pixel Hits";
0278 
0279   SiPixelRecHitCollection::const_iterator recHitIdIterator = (recHitColl.product())->begin();
0280   SiPixelRecHitCollection::const_iterator recHitIdIteratorEnd = (recHitColl.product())->end();
0281 
0282   int numOfRecHits = 0;
0283 
0284   int numberOfDetUnits1 = 0;
0285   int numOfRecHitsPerDet1 = 0;
0286   int numOfRecHitsPerLay1 = 0;
0287   int numberOfDetUnits2 = 0;
0288   int numOfRecHitsPerDet2 = 0;
0289   int numOfRecHitsPerLay2 = 0;
0290   int numberOfDetUnits3 = 0;
0291   int numOfRecHitsPerDet3 = 0;
0292   int numOfRecHitsPerLay3 = 0;
0293   int numberOfDetUnits1F = 0;
0294   int numOfRecHitsPerDet1F = 0;
0295   int numOfRecHitsPerLay1F = 0;
0296   int numberOfDetUnits2F = 0;
0297   int numOfRecHitsPerDet2F = 0;
0298   int numOfRecHitsPerLay2F = 0;
0299 
0300   // Loop over Detector IDs
0301   for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0302     SiPixelRecHitCollection::DetSet detset = *recHitIdIterator;
0303 
0304     DetId detId = DetId(detset.detId());    // Get the Detid object
0305     unsigned int detType = detId.det();     // det type, tracker=1
0306     unsigned int subid = detId.subdetId();  //subdetector type, barrel=1, fpix=2
0307 
0308     if (detType != 1)
0309       edm::LogPrint("ReadPixelRecHit") << " wrong det id " << detType;  //look only at tracker
0310     if (!((subid == 1) || (subid == 2)))
0311       edm::LogPrint("ReadPixelRecHit") << " wrong sub det id " << subid;  // look only at bpix&fpix
0312 
0313     if (print)
0314       edm::LogPrint("ReadPixelRecHit") << "     Det ID " << detId.rawId();
0315 
0316     //  Get rechits
0317     if (detset.empty())
0318       continue;
0319 
0320     // Get the geometrical information for this det
0321     const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
0322     double detZ = theGeomDet->surface().position().z();
0323     double detR = theGeomDet->surface().position().perp();
0324 
0325     //const BoundPlane& plane = theGeomDet->surface(); //for transf.  unused
0326     //double detThick = theGeomDet->specificSurface().bounds().thickness(); unused
0327 
0328     //const RectangularPixelTopology * topol =
0329     //dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
0330 
0331     const PixelTopology *topol = &(theGeomDet->specificTopology());
0332 
0333     //int cols = theGeomDet->specificTopology().ncolumns(); UNUSED
0334     //int rows = theGeomDet->specificTopology().nrows();
0335 
0336     float alignErrorX = 0., alignErrorY = 0.;
0337     LocalError lape = theGeomDet->localAlignmentError();
0338     //edm::LogPrint("ReadPixelRecHit")<<lape.valid();
0339     if (lape.valid()) {
0340       if (lape.xx() > 0.)
0341         alignErrorX = sqrt(lape.xx()) * 1E4;
0342       //float tmp12= sqrt(lape.xy())*1E4;
0343       if (lape.yy() > 0.)
0344         alignErrorY = sqrt(lape.yy()) * 1E4;
0345       if (print)
0346         edm::LogPrint("ReadPixelRecHit") << " Alignment errors " << alignErrorX << " " << alignErrorY;
0347       //edm::LogPrint("ReadPixelRecHit")<<" Alignment errors "<<alignErrorX<<" "<<alignErrorY;
0348     }
0349 
0350     unsigned int layer = 0, disk = 0, ladder = 0, zindex = 0, blade = 0, panel = 0, side = 0;
0351     if (subid == 1) {  // Subdet it, pix barrel=1
0352 
0353       PXBDetId pdetId = PXBDetId(detId);
0354       //unsigned int detTypeP=pdetId.det();   unused
0355       //unsigned int subidP=pdetId.subdetId(); unused
0356       // Barell layer = 1,2,3
0357       layer = pdetId.layer();
0358       // Barrel ladder id 1-20,32,44.
0359       ladder = pdetId.ladder();
0360       // Barrel Z-index=1,8
0361       zindex = pdetId.module();
0362       if (zindex < 5)
0363         side = 1;
0364       else
0365         side = 2;
0366 
0367       if (localPrint)
0368         edm::LogPrint("ReadPixelRecHit") << " Layer " << layer << " ladder " << ladder << " z " << zindex;
0369         //<<pdetId.rawId()<<" "<<pdetId.null()<<detTypeP<<" "<<subidP<<" ";
0370 
0371 #ifdef DO_HISTO
0372       hdetr->Fill(detR);
0373       hdetz->Fill(detZ);
0374       //hcolsB->Fill(float(cols));
0375       //hyposB->Fill(float(rows));
0376       hlayerid->Fill(float(layer));
0377 #endif
0378 
0379     } else {  // FPIX -------------------------------------
0380 
0381       // test cols & rows
0382       //edm::LogPrint("ReadPixelRecHit")<<" det z/r "<<detZ<<" "<<detR<<" "<<detThick<<" "
0383       //  <<cols<<" "<<rows;
0384 
0385       PXFDetId pdetId = PXFDetId(detId.rawId());
0386       disk = pdetId.disk();                   //1,2,3
0387       blade = pdetId.blade();                 //1-24
0388       side = pdetId.side();                   //size=1 for -z, 2 for +z
0389       panel = pdetId.panel();                 //panel=1,2
0390       unsigned int module = pdetId.module();  // plaquette
0391 
0392       if (print)
0393         edm::LogPrint("ReadPixelRecHit") << " forward hit " << side << " " << disk << " " << blade << " " << panel
0394                                          << " " << module;
0395 
0396 #ifdef DO_HISTO
0397       hdetrF->Fill(detR);
0398       hdetzF->Fill(detZ);
0399       hdisk->Fill(float(disk));
0400       hblade->Fill(float(blade));
0401       hmodule->Fill(float(module));
0402       hpanel->Fill(float(panel));
0403       hside->Fill(float(side));
0404 #endif
0405 
0406     }  // end BPix FPix if
0407 
0408 #ifdef DO_HISTO
0409     if (layer == 1) {
0410       hladder1id->Fill(float(ladder));
0411       hz1id->Fill(float(zindex));
0412       hAlignErrorX1->Fill(alignErrorX);
0413       hAlignErrorY1->Fill(alignErrorY);
0414       hAErrorXB->Fill(float(ladder + (110 * (side - 1))), alignErrorX);
0415       hAErrorYB->Fill(float(ladder + (110 * (side - 1))), alignErrorY);
0416       ++numberOfDetUnits1;
0417       numOfRecHitsPerDet1 = 0;
0418 
0419     } else if (layer == 2) {
0420       hladder2id->Fill(float(ladder));
0421       hz2id->Fill(float(zindex));
0422       hAlignErrorX2->Fill(alignErrorX);
0423       hAlignErrorY2->Fill(alignErrorY);
0424       hAErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorX);
0425       hAErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorY);
0426       ++numberOfDetUnits2;
0427       numOfRecHitsPerDet2 = 0;
0428 
0429     } else if (layer == 3) {
0430       hladder3id->Fill(float(ladder));
0431       hz3id->Fill(float(zindex));
0432       hAlignErrorX3->Fill(alignErrorX);
0433       hAlignErrorY3->Fill(alignErrorY);
0434       hAErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorX);
0435       hAErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorY);
0436       ++numberOfDetUnits3;
0437       numOfRecHitsPerDet3 = 0;
0438 
0439     } else if (disk == 1) {
0440       ++numberOfDetUnits1F;
0441       numOfRecHitsPerDet1F = 0;
0442       if (side == 1) {
0443         hAlignErrorX4->Fill(alignErrorX);
0444         hAlignErrorY4->Fill(alignErrorY);
0445         hAErrorXF->Fill(float(blade + 25), alignErrorX);
0446         hAErrorYF->Fill(float(blade + 25), alignErrorY);
0447       } else {
0448         hAlignErrorX6->Fill(alignErrorX);
0449         hAlignErrorY6->Fill(alignErrorY);
0450         hAErrorXF->Fill(float(blade + 50), alignErrorX);
0451         hAErrorYF->Fill(float(blade + 50), alignErrorY);
0452       }
0453 
0454     } else if (disk == 2) {
0455       ++numberOfDetUnits2F;
0456       numOfRecHitsPerDet2F = 0;
0457       if (side == 1) {
0458         hAlignErrorX5->Fill(alignErrorX);
0459         hAlignErrorY5->Fill(alignErrorY);
0460         hAErrorXF->Fill(float(blade), alignErrorX);
0461         hAErrorYF->Fill(float(blade), alignErrorY);
0462       } else {
0463         hAlignErrorX7->Fill(alignErrorX);
0464         hAlignErrorY7->Fill(alignErrorY);
0465         hAErrorXF->Fill(float(blade + 75), alignErrorX);
0466         hAErrorYF->Fill(float(blade + 75), alignErrorY);
0467       }
0468     }
0469 #endif
0470 
0471     //----Loop over rechits for this detId
0472     SiPixelRecHitCollection::DetSet::const_iterator pixeliter = detset.begin();
0473     SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0474     for (; pixeliter != rechitRangeIteratorEnd; ++pixeliter) {  //loop on the rechit
0475 
0476       numOfRecHits++;
0477 
0478       // RecHit local position is now transient,
0479       // one needs to run tracking to get position OR rerun localreco
0480       LocalPoint lp = pixeliter->localPosition();
0481       LocalError le = pixeliter->localPositionError();
0482       float xRecHit = lp.x();
0483       float yRecHit = lp.y();
0484       float xerror = sqrt(le.xx()) * 1E4;
0485       float yerror = sqrt(le.yy()) * 1E4;
0486       if (print)
0487         edm::LogPrint("ReadPixelRecHit") << " RecHit: " << numOfRecHits << " x/y " << xRecHit << " " << yRecHit
0488                                          << " errors x/y " << xerror << " " << yerror;
0489 
0490       //MeasurementPoint mp = topol->measurementPosition(xRecHit,yRecHit);
0491       //GlobalPoint GP = PixGeom->surface().toGlobal(Local3DPoint(lp));
0492 
0493       // Get cluster
0494       edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = pixeliter->cluster();
0495 
0496       float ch = (clust->charge()) / 1000.;  // convert ke to electrons
0497       int size = clust->size();
0498       int sizeX = clust->sizeX();
0499       int sizeY = clust->sizeY();
0500       float xClu = clust->x();
0501       float yClu = clust->y();
0502       int maxPixelCol = clust->maxPixelCol();
0503       int maxPixelRow = clust->maxPixelRow();
0504       int minPixelCol = clust->minPixelCol();
0505       int minPixelRow = clust->minPixelRow();
0506 
0507       // unsigned int geoId = clust->geographicalId(); // alsways 0
0508 
0509       // edge method moved to topologu class
0510       bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
0511       bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
0512 
0513       if (print)
0514         edm::LogPrint("ReadPixelRecHit") << "Clu: charge " << ch << " size " << size << " size x/y " << sizeX << " "
0515                                          << sizeY << " meas. " << xClu << " " << yClu << " edge " << edgeHitX << " "
0516                                          << edgeHitY;
0517 
0518       if (print)
0519         edm::LogPrint("ReadPixelRecHit") << " pixels:";
0520 
0521       // Get the pixels in the Cluster
0522       const vector<SiPixelCluster::Pixel> &pixelsVec = clust->pixels();
0523       //if(localPrint) edm::LogPrint("ReadPixelRecHit")<<" Pixels in this cluster ";
0524       map<unsigned int, float, less<unsigned int> > chanMap;  // Channel map
0525       // Look at pixels in this cluster. ADC is calibrated, in electrons
0526       for (unsigned int i = 0; i < pixelsVec.size(); ++i) {
0527         float pixx = pixelsVec[i].x;  // index as a float so = i+0.5
0528         float pixy = pixelsVec[i].y;
0529         float adc = ((pixelsVec[i].adc) / 1000);  // in kelec.
0530 
0531         // OLD way
0532         //int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy));
0533         bool bigInX = topol->isItBigPixelInX(int(pixx));
0534         bool bigInY = topol->isItBigPixelInY(int(pixy));
0535 
0536         bool edgeInX = topol->isItEdgePixelInX(int(pixx));
0537         bool edgeInY = topol->isItEdgePixelInY(int(pixy));
0538 
0539         if (print)
0540           edm::LogPrint("ReadPixelRecHit") << i << " index " << pixx << " " << pixy << " adc " << adc << " edge "
0541                                            << edgeInX << " " << edgeInY << " big " << bigInX << " " << bigInY;
0542 
0543           //if(print && sizeX==1 && bigInX)
0544           //edm::LogPrint("ReadPixelRecHit")<<" single big x "<<xClu<<" "<<pixx<<" ";
0545           //if(print && sizeY==1 && bigInY)
0546           //edm::LogPrint("ReadPixelRecHit")<<" single big y "<<yClu<<" "<<pixy<<" ";
0547 #ifdef DO_HISTO
0548         if (layer == 1) {
0549           hadcCharge1->Fill(adc);
0550           //if(bigInX || bigInY) hadcCharge1big->Fill(adc);
0551         } else if (layer == 2) {
0552           hadcCharge2->Fill(adc);
0553         } else if (layer == 3) {
0554           hadcCharge3->Fill(adc);
0555         } else if (disk == 1) {
0556           hadcCharge1F->Fill(adc);
0557         } else if (disk == 2) {
0558           hadcCharge2F->Fill(adc);
0559         }
0560 #endif
0561       }  // End pixel loop
0562 
0563 #ifdef DO_HISTO
0564       if (layer == 1) {
0565         hcharge1->Fill(ch);
0566         hxpos1->Fill(yRecHit);
0567         hypos1->Fill(xRecHit);
0568         hsize1->Fill(float(size));
0569         hsizex1->Fill(float(sizeX));
0570         hsizey1->Fill(float(sizeY));
0571         numOfRecHitsPerDet1++;
0572         numOfRecHitsPerLay1++;
0573         hErrorX1->Fill(xerror);
0574         hErrorY1->Fill(yerror);
0575         hErrorXB->Fill(float(ladder + (110 * (side - 1))), xerror);
0576         hErrorYB->Fill(float(ladder + (110 * (side - 1))), yerror);
0577 
0578       } else if (layer == 2) {  // layer 2
0579 
0580         hcharge2->Fill(ch);
0581         hxpos2->Fill(yRecHit);
0582         hypos2->Fill(xRecHit);
0583         hsize2->Fill(float(size));
0584         hsizex2->Fill(float(sizeX));
0585         hsizey2->Fill(float(sizeY));
0586         numOfRecHitsPerDet2++;
0587         numOfRecHitsPerLay2++;
0588         hErrorX2->Fill(xerror);
0589         hErrorY2->Fill(yerror);
0590         hErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), xerror);
0591         hErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), yerror);
0592 
0593       } else if (layer == 3) {  // Layer 3
0594 
0595         hcharge3->Fill(ch);
0596         hxpos3->Fill(yRecHit);
0597         hypos3->Fill(xRecHit);
0598         hsize3->Fill(float(size));
0599         hsizex3->Fill(float(sizeX));
0600         hsizey3->Fill(float(sizeY));
0601         numOfRecHitsPerDet3++;
0602         numOfRecHitsPerLay3++;
0603         hErrorX3->Fill(xerror);
0604         hErrorY3->Fill(yerror);
0605         hErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), xerror);
0606         hErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), yerror);
0607 
0608       } else if (disk == 1) {
0609         hcharge1F->Fill(ch);
0610         hxpos1F->Fill(yRecHit);
0611         hypos1F->Fill(xRecHit);
0612         hsize1F->Fill(float(size));
0613         hsizex1F->Fill(float(sizeX));
0614         hsizey1F->Fill(float(sizeY));
0615         if (side == 1) {  // -z
0616           hErrorX4->Fill(xerror);
0617           hErrorY4->Fill(yerror);
0618           hErrorXF->Fill(float(blade + 25), xerror);
0619           hErrorYF->Fill(float(blade + 25), yerror);
0620         } else {  // +z
0621           hErrorX6->Fill(xerror);
0622           hErrorY6->Fill(yerror);
0623           hErrorXF->Fill(float(blade + 50), xerror);
0624           hErrorYF->Fill(float(blade + 50), yerror);
0625         }
0626         numOfRecHitsPerDet1F++;
0627         numOfRecHitsPerLay1F++;
0628       } else if (disk == 2) {  // disk 2
0629 
0630         hcharge2F->Fill(ch);
0631         hxpos2F->Fill(yRecHit);
0632         hypos2F->Fill(xRecHit);
0633         hsize2F->Fill(float(size));
0634         hsizex2F->Fill(float(sizeX));
0635         hsizey2F->Fill(float(sizeY));
0636         numOfRecHitsPerDet2F++;
0637         numOfRecHitsPerLay2F++;
0638         if (side == 1) {  // -z
0639           hErrorX5->Fill(xerror);
0640           hErrorY5->Fill(yerror);
0641           hErrorXF->Fill(float(blade), xerror);
0642           hErrorYF->Fill(float(blade), yerror);
0643         } else {  // +z
0644           hErrorX7->Fill(xerror);
0645           hErrorY7->Fill(yerror);
0646           hErrorXF->Fill(float(blade + 75), xerror);
0647           hErrorYF->Fill(float(blade + 75), yerror);
0648         }
0649       }
0650 #endif
0651 
0652     }  // End RecHit loop
0653 
0654 #ifdef DO_HISTO
0655     if (layer == 1)
0656       hrecHitsPerDet1->Fill(float(numOfRecHitsPerDet1));
0657     else if (layer == 2)
0658       hrecHitsPerDet2->Fill(float(numOfRecHitsPerDet2));
0659     else if (layer == 3)
0660       hrecHitsPerDet3->Fill(float(numOfRecHitsPerDet3));
0661     else if (disk == 1)
0662       hrecHitsPerDet1F->Fill(float(numOfRecHitsPerDet1F));
0663     else if (disk == 2)
0664       hrecHitsPerDet2F->Fill(float(numOfRecHitsPerDet2F));
0665 #endif
0666 
0667   }  // End Det loop
0668 
0669 #ifdef DO_HISTO
0670   if (print)
0671     edm::LogPrint("ReadPixelRecHit") << " Rechits per layer " << numOfRecHitsPerLay1 << " " << numOfRecHitsPerLay2
0672                                      << " " << numOfRecHitsPerLay3 << " " << numOfRecHitsPerLay1F << " "
0673                                      << numOfRecHitsPerLay2F;
0674   hrecHitsPerLay1->Fill(float(numOfRecHitsPerLay1));
0675   hrecHitsPerLay2->Fill(float(numOfRecHitsPerLay2));
0676   hrecHitsPerLay3->Fill(float(numOfRecHitsPerLay3));
0677   hdetsPerLay1->Fill(float(numberOfDetUnits1));
0678   hdetsPerLay2->Fill(float(numberOfDetUnits2));
0679   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0680 
0681   hrecHitsPerLay1F->Fill(float(numOfRecHitsPerLay1F));
0682   hrecHitsPerLay2F->Fill(float(numOfRecHitsPerLay2F));
0683   hdetsPerLay1F->Fill(float(numberOfDetUnits1F));
0684   hdetsPerLay2F->Fill(float(numberOfDetUnits2F));
0685   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0686 #endif
0687 }
0688 
0689 //define this as a plug-in
0690 DEFINE_FWK_MODULE(ReadPixelRecHit);
[ 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!