Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-17 23:31:50

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 numberOfDetUnits = 0;
0283   int numOfRecHits = 0;
0284 
0285   int numberOfDetUnits1 = 0;
0286   int numOfRecHitsPerDet1 = 0;
0287   int numOfRecHitsPerLay1 = 0;
0288   int numberOfDetUnits2 = 0;
0289   int numOfRecHitsPerDet2 = 0;
0290   int numOfRecHitsPerLay2 = 0;
0291   int numberOfDetUnits3 = 0;
0292   int numOfRecHitsPerDet3 = 0;
0293   int numOfRecHitsPerLay3 = 0;
0294   int numberOfDetUnits1F = 0;
0295   int numOfRecHitsPerDet1F = 0;
0296   int numOfRecHitsPerLay1F = 0;
0297   int numberOfDetUnits2F = 0;
0298   int numOfRecHitsPerDet2F = 0;
0299   int numOfRecHitsPerLay2F = 0;
0300 
0301   // Loop over Detector IDs
0302   for (; recHitIdIterator != recHitIdIteratorEnd; recHitIdIterator++) {
0303     SiPixelRecHitCollection::DetSet detset = *recHitIdIterator;
0304 
0305     DetId detId = DetId(detset.detId());    // Get the Detid object
0306     unsigned int detType = detId.det();     // det type, tracker=1
0307     unsigned int subid = detId.subdetId();  //subdetector type, barrel=1, fpix=2
0308 
0309     if (detType != 1)
0310       edm::LogPrint("ReadPixelRecHit") << " wrong det id " << detType;  //look only at tracker
0311     if (!((subid == 1) || (subid == 2)))
0312       edm::LogPrint("ReadPixelRecHit") << " wrong sub det id " << subid;  // look only at bpix&fpix
0313 
0314     if (print)
0315       edm::LogPrint("ReadPixelRecHit") << "     Det ID " << detId.rawId();
0316 
0317     //  Get rechits
0318     if (detset.empty())
0319       continue;
0320 
0321     // Get the geometrical information for this det
0322     const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
0323     double detZ = theGeomDet->surface().position().z();
0324     double detR = theGeomDet->surface().position().perp();
0325 
0326     //const BoundPlane& plane = theGeomDet->surface(); //for transf.  unused
0327     //double detThick = theGeomDet->specificSurface().bounds().thickness(); unused
0328 
0329     //const RectangularPixelTopology * topol =
0330     //dynamic_cast<const RectangularPixelTopology*>(&(theGeomDet->specificTopology()));
0331 
0332     const PixelTopology *topol = &(theGeomDet->specificTopology());
0333 
0334     //int cols = theGeomDet->specificTopology().ncolumns(); UNUSED
0335     //int rows = theGeomDet->specificTopology().nrows();
0336 
0337     float alignErrorX = 0., alignErrorY = 0.;
0338     LocalError lape = theGeomDet->localAlignmentError();
0339     //edm::LogPrint("ReadPixelRecHit")<<lape.valid();
0340     if (lape.valid()) {
0341       if (lape.xx() > 0.)
0342         alignErrorX = sqrt(lape.xx()) * 1E4;
0343       //float tmp12= sqrt(lape.xy())*1E4;
0344       if (lape.yy() > 0.)
0345         alignErrorY = sqrt(lape.yy()) * 1E4;
0346       if (print)
0347         edm::LogPrint("ReadPixelRecHit") << " Alignment errors " << alignErrorX << " " << alignErrorY;
0348       //edm::LogPrint("ReadPixelRecHit")<<" Alignment errors "<<alignErrorX<<" "<<alignErrorY;
0349     }
0350 
0351     unsigned int layer = 0, disk = 0, ladder = 0, zindex = 0, blade = 0, panel = 0, side = 0;
0352     if (subid == 1) {  // Subdet it, pix barrel=1
0353       ++numberOfDetUnits;
0354 
0355       PXBDetId pdetId = PXBDetId(detId);
0356       //unsigned int detTypeP=pdetId.det();   unused
0357       //unsigned int subidP=pdetId.subdetId(); unused
0358       // Barell layer = 1,2,3
0359       layer = pdetId.layer();
0360       // Barrel ladder id 1-20,32,44.
0361       ladder = pdetId.ladder();
0362       // Barrel Z-index=1,8
0363       zindex = pdetId.module();
0364       if (zindex < 5)
0365         side = 1;
0366       else
0367         side = 2;
0368 
0369       if (localPrint)
0370         edm::LogPrint("ReadPixelRecHit") << " Layer " << layer << " ladder " << ladder << " z " << zindex;
0371         //<<pdetId.rawId()<<" "<<pdetId.null()<<detTypeP<<" "<<subidP<<" ";
0372 
0373 #ifdef DO_HISTO
0374       hdetr->Fill(detR);
0375       hdetz->Fill(detZ);
0376       //hcolsB->Fill(float(cols));
0377       //hyposB->Fill(float(rows));
0378       hlayerid->Fill(float(layer));
0379 #endif
0380 
0381     } else {  // FPIX -------------------------------------
0382 
0383       // test cols & rows
0384       //edm::LogPrint("ReadPixelRecHit")<<" det z/r "<<detZ<<" "<<detR<<" "<<detThick<<" "
0385       //  <<cols<<" "<<rows;
0386 
0387       PXFDetId pdetId = PXFDetId(detId.rawId());
0388       disk = pdetId.disk();                   //1,2,3
0389       blade = pdetId.blade();                 //1-24
0390       side = pdetId.side();                   //size=1 for -z, 2 for +z
0391       panel = pdetId.panel();                 //panel=1,2
0392       unsigned int module = pdetId.module();  // plaquette
0393 
0394       if (print)
0395         edm::LogPrint("ReadPixelRecHit") << " forward hit " << side << " " << disk << " " << blade << " " << panel
0396                                          << " " << module;
0397 
0398 #ifdef DO_HISTO
0399       hdetrF->Fill(detR);
0400       hdetzF->Fill(detZ);
0401       hdisk->Fill(float(disk));
0402       hblade->Fill(float(blade));
0403       hmodule->Fill(float(module));
0404       hpanel->Fill(float(panel));
0405       hside->Fill(float(side));
0406 #endif
0407 
0408     }  // end BPix FPix if
0409 
0410 #ifdef DO_HISTO
0411     if (layer == 1) {
0412       hladder1id->Fill(float(ladder));
0413       hz1id->Fill(float(zindex));
0414       hAlignErrorX1->Fill(alignErrorX);
0415       hAlignErrorY1->Fill(alignErrorY);
0416       hAErrorXB->Fill(float(ladder + (110 * (side - 1))), alignErrorX);
0417       hAErrorYB->Fill(float(ladder + (110 * (side - 1))), alignErrorY);
0418       ++numberOfDetUnits1;
0419       numOfRecHitsPerDet1 = 0;
0420 
0421     } else if (layer == 2) {
0422       hladder2id->Fill(float(ladder));
0423       hz2id->Fill(float(zindex));
0424       hAlignErrorX2->Fill(alignErrorX);
0425       hAlignErrorY2->Fill(alignErrorY);
0426       hAErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorX);
0427       hAErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), alignErrorY);
0428       ++numberOfDetUnits2;
0429       numOfRecHitsPerDet2 = 0;
0430 
0431     } else if (layer == 3) {
0432       hladder3id->Fill(float(ladder));
0433       hz3id->Fill(float(zindex));
0434       hAlignErrorX3->Fill(alignErrorX);
0435       hAlignErrorY3->Fill(alignErrorY);
0436       hAErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorX);
0437       hAErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), alignErrorY);
0438       ++numberOfDetUnits3;
0439       numOfRecHitsPerDet3 = 0;
0440 
0441     } else if (disk == 1) {
0442       ++numberOfDetUnits1F;
0443       numOfRecHitsPerDet1F = 0;
0444       if (side == 1) {
0445         hAlignErrorX4->Fill(alignErrorX);
0446         hAlignErrorY4->Fill(alignErrorY);
0447         hAErrorXF->Fill(float(blade + 25), alignErrorX);
0448         hAErrorYF->Fill(float(blade + 25), alignErrorY);
0449       } else {
0450         hAlignErrorX6->Fill(alignErrorX);
0451         hAlignErrorY6->Fill(alignErrorY);
0452         hAErrorXF->Fill(float(blade + 50), alignErrorX);
0453         hAErrorYF->Fill(float(blade + 50), alignErrorY);
0454       }
0455 
0456     } else if (disk == 2) {
0457       ++numberOfDetUnits2F;
0458       numOfRecHitsPerDet2F = 0;
0459       if (side == 1) {
0460         hAlignErrorX5->Fill(alignErrorX);
0461         hAlignErrorY5->Fill(alignErrorY);
0462         hAErrorXF->Fill(float(blade), alignErrorX);
0463         hAErrorYF->Fill(float(blade), alignErrorY);
0464       } else {
0465         hAlignErrorX7->Fill(alignErrorX);
0466         hAlignErrorY7->Fill(alignErrorY);
0467         hAErrorXF->Fill(float(blade + 75), alignErrorX);
0468         hAErrorYF->Fill(float(blade + 75), alignErrorY);
0469       }
0470     }
0471 #endif
0472 
0473     //----Loop over rechits for this detId
0474     SiPixelRecHitCollection::DetSet::const_iterator pixeliter = detset.begin();
0475     SiPixelRecHitCollection::DetSet::const_iterator rechitRangeIteratorEnd = detset.end();
0476     for (; pixeliter != rechitRangeIteratorEnd; ++pixeliter) {  //loop on the rechit
0477 
0478       numOfRecHits++;
0479 
0480       // RecHit local position is now transient,
0481       // one needs to run tracking to get position OR rerun localreco
0482       LocalPoint lp = pixeliter->localPosition();
0483       LocalError le = pixeliter->localPositionError();
0484       float xRecHit = lp.x();
0485       float yRecHit = lp.y();
0486       float xerror = sqrt(le.xx()) * 1E4;
0487       float yerror = sqrt(le.yy()) * 1E4;
0488       if (print)
0489         edm::LogPrint("ReadPixelRecHit") << " RecHit: " << numOfRecHits << " x/y " << xRecHit << " " << yRecHit
0490                                          << " errors x/y " << xerror << " " << yerror;
0491 
0492       //MeasurementPoint mp = topol->measurementPosition(xRecHit,yRecHit);
0493       //GlobalPoint GP = PixGeom->surface().toGlobal(Local3DPoint(lp));
0494 
0495       // Get cluster
0496       edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = pixeliter->cluster();
0497 
0498       float ch = (clust->charge()) / 1000.;  // convert ke to electrons
0499       int size = clust->size();
0500       int sizeX = clust->sizeX();
0501       int sizeY = clust->sizeY();
0502       float xClu = clust->x();
0503       float yClu = clust->y();
0504       int maxPixelCol = clust->maxPixelCol();
0505       int maxPixelRow = clust->maxPixelRow();
0506       int minPixelCol = clust->minPixelCol();
0507       int minPixelRow = clust->minPixelRow();
0508 
0509       // unsigned int geoId = clust->geographicalId(); // alsways 0
0510 
0511       // edge method moved to topologu class
0512       bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
0513       bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
0514 
0515       if (print)
0516         edm::LogPrint("ReadPixelRecHit") << "Clu: charge " << ch << " size " << size << " size x/y " << sizeX << " "
0517                                          << sizeY << " meas. " << xClu << " " << yClu << " edge " << edgeHitX << " "
0518                                          << edgeHitY;
0519 
0520       if (print)
0521         edm::LogPrint("ReadPixelRecHit") << " pixels:";
0522 
0523       // Get the pixels in the Cluster
0524       const vector<SiPixelCluster::Pixel> &pixelsVec = clust->pixels();
0525       //if(localPrint) edm::LogPrint("ReadPixelRecHit")<<" Pixels in this cluster ";
0526       map<unsigned int, float, less<unsigned int> > chanMap;  // Channel map
0527       // Look at pixels in this cluster. ADC is calibrated, in electrons
0528       for (unsigned int i = 0; i < pixelsVec.size(); ++i) {
0529         float pixx = pixelsVec[i].x;  // index as a float so = i+0.5
0530         float pixy = pixelsVec[i].y;
0531         float adc = ((pixelsVec[i].adc) / 1000);  // in kelec.
0532 
0533         // OLD way
0534         //int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy));
0535         bool bigInX = topol->isItBigPixelInX(int(pixx));
0536         bool bigInY = topol->isItBigPixelInY(int(pixy));
0537 
0538         bool edgeInX = topol->isItEdgePixelInX(int(pixx));
0539         bool edgeInY = topol->isItEdgePixelInY(int(pixy));
0540 
0541         if (print)
0542           edm::LogPrint("ReadPixelRecHit") << i << " index " << pixx << " " << pixy << " adc " << adc << " edge "
0543                                            << edgeInX << " " << edgeInY << " big " << bigInX << " " << bigInY;
0544 
0545           //if(print && sizeX==1 && bigInX)
0546           //edm::LogPrint("ReadPixelRecHit")<<" single big x "<<xClu<<" "<<pixx<<" ";
0547           //if(print && sizeY==1 && bigInY)
0548           //edm::LogPrint("ReadPixelRecHit")<<" single big y "<<yClu<<" "<<pixy<<" ";
0549 #ifdef DO_HISTO
0550         if (layer == 1) {
0551           hadcCharge1->Fill(adc);
0552           //if(bigInX || bigInY) hadcCharge1big->Fill(adc);
0553         } else if (layer == 2) {
0554           hadcCharge2->Fill(adc);
0555         } else if (layer == 3) {
0556           hadcCharge3->Fill(adc);
0557         } else if (disk == 1) {
0558           hadcCharge1F->Fill(adc);
0559         } else if (disk == 2) {
0560           hadcCharge2F->Fill(adc);
0561         }
0562 #endif
0563       }  // End pixel loop
0564 
0565 #ifdef DO_HISTO
0566       if (layer == 1) {
0567         hcharge1->Fill(ch);
0568         hxpos1->Fill(yRecHit);
0569         hypos1->Fill(xRecHit);
0570         hsize1->Fill(float(size));
0571         hsizex1->Fill(float(sizeX));
0572         hsizey1->Fill(float(sizeY));
0573         numOfRecHitsPerDet1++;
0574         numOfRecHitsPerLay1++;
0575         hErrorX1->Fill(xerror);
0576         hErrorY1->Fill(yerror);
0577         hErrorXB->Fill(float(ladder + (110 * (side - 1))), xerror);
0578         hErrorYB->Fill(float(ladder + (110 * (side - 1))), yerror);
0579 
0580       } else if (layer == 2) {  // layer 2
0581 
0582         hcharge2->Fill(ch);
0583         hxpos2->Fill(yRecHit);
0584         hypos2->Fill(xRecHit);
0585         hsize2->Fill(float(size));
0586         hsizex2->Fill(float(sizeX));
0587         hsizey2->Fill(float(sizeY));
0588         numOfRecHitsPerDet2++;
0589         numOfRecHitsPerLay2++;
0590         hErrorX2->Fill(xerror);
0591         hErrorY2->Fill(yerror);
0592         hErrorXB->Fill(float(ladder + 25 + (110 * (side - 1))), xerror);
0593         hErrorYB->Fill(float(ladder + 25 + (110 * (side - 1))), yerror);
0594 
0595       } else if (layer == 3) {  // Layer 3
0596 
0597         hcharge3->Fill(ch);
0598         hxpos3->Fill(yRecHit);
0599         hypos3->Fill(xRecHit);
0600         hsize3->Fill(float(size));
0601         hsizex3->Fill(float(sizeX));
0602         hsizey3->Fill(float(sizeY));
0603         numOfRecHitsPerDet3++;
0604         numOfRecHitsPerLay3++;
0605         hErrorX3->Fill(xerror);
0606         hErrorY3->Fill(yerror);
0607         hErrorXB->Fill(float(ladder + 60 + (110 * (side - 1))), xerror);
0608         hErrorYB->Fill(float(ladder + 60 + (110 * (side - 1))), yerror);
0609 
0610       } else if (disk == 1) {
0611         hcharge1F->Fill(ch);
0612         hxpos1F->Fill(yRecHit);
0613         hypos1F->Fill(xRecHit);
0614         hsize1F->Fill(float(size));
0615         hsizex1F->Fill(float(sizeX));
0616         hsizey1F->Fill(float(sizeY));
0617         if (side == 1) {  // -z
0618           hErrorX4->Fill(xerror);
0619           hErrorY4->Fill(yerror);
0620           hErrorXF->Fill(float(blade + 25), xerror);
0621           hErrorYF->Fill(float(blade + 25), yerror);
0622         } else {  // +z
0623           hErrorX6->Fill(xerror);
0624           hErrorY6->Fill(yerror);
0625           hErrorXF->Fill(float(blade + 50), xerror);
0626           hErrorYF->Fill(float(blade + 50), yerror);
0627         }
0628         numOfRecHitsPerDet1F++;
0629         numOfRecHitsPerLay1F++;
0630       } else if (disk == 2) {  // disk 2
0631 
0632         hcharge2F->Fill(ch);
0633         hxpos2F->Fill(yRecHit);
0634         hypos2F->Fill(xRecHit);
0635         hsize2F->Fill(float(size));
0636         hsizex2F->Fill(float(sizeX));
0637         hsizey2F->Fill(float(sizeY));
0638         numOfRecHitsPerDet2F++;
0639         numOfRecHitsPerLay2F++;
0640         if (side == 1) {  // -z
0641           hErrorX5->Fill(xerror);
0642           hErrorY5->Fill(yerror);
0643           hErrorXF->Fill(float(blade), xerror);
0644           hErrorYF->Fill(float(blade), yerror);
0645         } else {  // +z
0646           hErrorX7->Fill(xerror);
0647           hErrorY7->Fill(yerror);
0648           hErrorXF->Fill(float(blade + 75), xerror);
0649           hErrorYF->Fill(float(blade + 75), yerror);
0650         }
0651       }
0652 #endif
0653 
0654     }  // End RecHit loop
0655 
0656 #ifdef DO_HISTO
0657     if (layer == 1)
0658       hrecHitsPerDet1->Fill(float(numOfRecHitsPerDet1));
0659     else if (layer == 2)
0660       hrecHitsPerDet2->Fill(float(numOfRecHitsPerDet2));
0661     else if (layer == 3)
0662       hrecHitsPerDet3->Fill(float(numOfRecHitsPerDet3));
0663     else if (disk == 1)
0664       hrecHitsPerDet1F->Fill(float(numOfRecHitsPerDet1F));
0665     else if (disk == 2)
0666       hrecHitsPerDet2F->Fill(float(numOfRecHitsPerDet2F));
0667 #endif
0668 
0669   }  // End Det loop
0670 
0671 #ifdef DO_HISTO
0672   if (print)
0673     edm::LogPrint("ReadPixelRecHit") << " Rechits per layer " << numOfRecHitsPerLay1 << " " << numOfRecHitsPerLay2
0674                                      << " " << numOfRecHitsPerLay3 << " " << numOfRecHitsPerLay1F << " "
0675                                      << numOfRecHitsPerLay2F;
0676   hrecHitsPerLay1->Fill(float(numOfRecHitsPerLay1));
0677   hrecHitsPerLay2->Fill(float(numOfRecHitsPerLay2));
0678   hrecHitsPerLay3->Fill(float(numOfRecHitsPerLay3));
0679   hdetsPerLay1->Fill(float(numberOfDetUnits1));
0680   hdetsPerLay2->Fill(float(numberOfDetUnits2));
0681   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0682 
0683   hrecHitsPerLay1F->Fill(float(numOfRecHitsPerLay1F));
0684   hrecHitsPerLay2F->Fill(float(numOfRecHitsPerLay2F));
0685   hdetsPerLay1F->Fill(float(numberOfDetUnits1F));
0686   hdetsPerLay2F->Fill(float(numberOfDetUnits2F));
0687   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0688 #endif
0689 }
0690 
0691 //define this as a plug-in
0692 DEFINE_FWK_MODULE(ReadPixelRecHit);