Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-22 04:03:20 

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