Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 // File: ReadPixClusters.cc
0002 // Description: To test the pixel clusters.
0003 // Author: Danek Kotlinski
0004 // Creation Date:  Initial version. 3/06
0005 // Modify to work with CMSSW620, 8/13, CMSSW700, 10/13 d.k.
0006 // Add ByToken data access.
0007 //--------------------------------------------
0008 #include <memory>
0009 #include <string>
0010 #include <iostream>
0011 
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0018 
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 
0026 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0027 #include "DataFormats/Common/interface/DetSetVector.h"
0028 #include "DataFormats/Common/interface/Ref.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 
0031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0032 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0033 
0034 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0035 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0036 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0037 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0038 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0039 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0041 
0042 // For L1
0043 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0044 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0045 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0046 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0047 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0048 
0049 // For HLT
0050 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0051 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0052 #include "DataFormats/Common/interface/TriggerResults.h"
0053 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0054 #include "FWCore/Common/interface/TriggerNames.h"
0055 
0056 // To use root histos
0057 #include "FWCore/ServiceRegistry/interface/Service.h"
0058 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0059 
0060 // For ROOT
0061 #include <TROOT.h>
0062 #include <TChain.h>
0063 #include <TFile.h>
0064 #include <TF1.h>
0065 #include <TH2F.h>
0066 #include <TH1F.h>
0067 
0068 #define NEW_ID
0069 #ifdef NEW_ID
0070 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0071 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0072 #else
0073 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0074 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0075 #endif
0076 
0077 #define HISTOS
0078 
0079 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0080 
0081 using namespace std;
0082 
0083 //=============================================================================
0084 
0085 class ReadPixClusters : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0086 public:
0087   explicit ReadPixClusters(const edm::ParameterSet &conf);
0088   virtual ~ReadPixClusters();
0089   virtual void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0090   virtual void beginJob() override;
0091   virtual void endJob() override;
0092 
0093 private:
0094   edm::InputTag src_;
0095   edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> tPixelCluster;
0096   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0097   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0098   bool printLocal;
0099   int countEvents, countAllEvents;
0100   double sumClusters;
0101 
0102   //TFile* hFile;
0103   TH1F *hdetunit;
0104   TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id;
0105 
0106   TH1F *hcharge1, *hcharge2, *hcharge3, *hcharge4, *hcharge5;
0107   TH1F *hpixcharge1, *hpixcharge2, *hpixcharge3, *hpixcharge4, *hpixcharge5;
0108   TH1F *hcols1, *hcols2, *hcols3, *hrows1, *hrows2, *hrows3;
0109   TH1F *hpcols1, *hpcols2, *hpcols3, *hprows1, *hprows2, *hprows3;
0110   TH1F *hsize1, *hsize2, *hsize3, *hsize4, *hsize5, *hsizex1, *hsizex2, *hsizex3, *hsizex4, *hsizex5, *hsizey1,
0111       *hsizey2, *hsizey3, *hsizey4, *hsizey5;
0112 
0113   TH1F *hclusPerDet1, *hclusPerDet2, *hclusPerDet3;
0114   TH1F *hpixPerDet1, *hpixPerDet2, *hpixPerDet3;
0115   TH1F *hpixPerLink1, *hpixPerLink2, *hpixPerLink3;
0116   TH1F *hclusPerLay1, *hclusPerLay2, *hclusPerLay3;
0117   TH1F *hpixPerLay1, *hpixPerLay2, *hpixPerLay3;
0118   TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3;
0119   TH1F *hclus, *hclusBPix, *hclusFPix, *hdigis, *hdigisB, *hdigisF;
0120 
0121   TH1F *hdetr, *hdetz;
0122   //   TH1F *hcolsB,  *hrowsB,  *hcolsF,  *hrowsF;
0123   //   TH2F *htest1, *htest2;
0124   TH2F *hDetMap1, *hDetMap2, *hDetMap3;
0125   TH2F *hpDetMap1, *hpDetMap2, *hpDetMap3;
0126   TH2F *hpixDetMap1, *hpixDetMap2, *hpixDetMap3;
0127   TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3;
0128 
0129   TH1F *hevent, *horbit, *hlumi, *hlumi0, *hlumi1;
0130   TH1F *hbx, *hbx0, *hbx1;
0131   TH1F *hdets, *hmbits1, *hmbits2, *hmbits3, *hmaxPixPerDet;
0132 
0133   TH1F *hclusPerDisk1, *hclusPerDisk2, *hclusPerDisk3, *hclusPerDisk4;
0134   TH1F *htest;
0135 };
0136 
0137 /////////////////////////////////////////////////////////////////
0138 // Contructor, empty.
0139 ReadPixClusters::ReadPixClusters(edm::ParameterSet const &conf) : src_(conf.getParameter<edm::InputTag>("src")) {
0140   usesResource(TFileService::kSharedResource);
0141   printLocal = conf.getUntrackedParameter<bool>("Verbosity", false);
0142   //src_ =  conf.getParameter<edm::InputTag>( "src" );
0143   edm::LogPrint("ReadPixClusters") << " Construct " << printLocal;
0144   tPixelCluster = consumes<edmNew::DetSetVector<SiPixelCluster>>(src_);
0145   trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0146   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0147 }
0148 // Virtual destructor needed.
0149 ReadPixClusters::~ReadPixClusters() = default;
0150 
0151 // ------------ method called at the begining   ------------
0152 void ReadPixClusters::beginJob() {
0153   edm::LogPrint("ReadPixClusters") << "Initialize PixelClusterTest " << printLocal;
0154 
0155 #ifdef HISTOS
0156 
0157   // NEW way to use root (from 2.0.0?)
0158   edm::Service<TFileService> fs;
0159 
0160   //=====================================================================
0161 
0162   hladder1id = fs->make<TH1F>("hladder1id", "Ladder L1 id", 23, -11.5, 11.5);
0163   hladder2id = fs->make<TH1F>("hladder2id", "Ladder L2 id", 35, -17.5, 17.5);
0164   hladder3id = fs->make<TH1F>("hladder3id", "Ladder L3 id", 47, -23.5, 23.5);
0165   hz1id = fs->make<TH1F>("hz1id", "Z-index id L1", 11, -5.5, 5.5);
0166   hz2id = fs->make<TH1F>("hz2id", "Z-index id L2", 11, -5.5, 5.5);
0167   hz3id = fs->make<TH1F>("hz3id", "Z-index id L3", 11, -5.5, 5.5);
0168 
0169   int sizeH = 200;
0170   float lowH = -0.5;
0171   float highH = 199.5;
0172 
0173   hclusPerDet1 = fs->make<TH1F>("hclusPerDet1", "Clus per det l1", sizeH, lowH, highH);
0174   hclusPerDet2 = fs->make<TH1F>("hclusPerDet2", "Clus per det l2", sizeH, lowH, highH);
0175   hclusPerDet3 = fs->make<TH1F>("hclusPerDet3", "Clus per det l3", sizeH, lowH, highH);
0176 
0177   sizeH = 1000;
0178   highH = 1999.5;
0179   hpixPerDet1 = fs->make<TH1F>("hpixPerDet1", "Pix per det l1", sizeH, lowH, highH);
0180   hpixPerDet2 = fs->make<TH1F>("hpixPerDet2", "Pix per det l2", sizeH, lowH, highH);
0181   hpixPerDet3 = fs->make<TH1F>("hpixPerDet3", "Pix per det l3", sizeH, lowH, highH);
0182 
0183   sizeH = 1000;
0184   highH = 999.5;
0185   hpixPerLink1 = fs->make<TH1F>("hpixPerLink1", "Pix per link l1", sizeH, lowH, highH);
0186   hpixPerLink2 = fs->make<TH1F>("hpixPerLink2", "Pix per link l2", sizeH, lowH, highH);
0187   hpixPerLink3 = fs->make<TH1F>("hpixPerLink3", "Pix per link l3", sizeH, lowH, highH);
0188 
0189   sizeH = 2000;
0190   highH = 1999.5;
0191   hclusPerLay1 = fs->make<TH1F>("hclusPerLay1", "Clus per layer l1", sizeH, lowH, highH);
0192   hclusPerLay2 = fs->make<TH1F>("hclusPerLay2", "Clus per layer l2", sizeH, lowH, highH);
0193   hclusPerLay3 = fs->make<TH1F>("hclusPerLay3", "Clus per layer l3", sizeH, lowH, highH);
0194 
0195   hclusPerDisk1 = fs->make<TH1F>("hclusPerDisk1", "Clus per disk1", sizeH, lowH, highH);
0196   hclusPerDisk2 = fs->make<TH1F>("hclusPerDisk2", "Clus per disk2", sizeH, lowH, highH);
0197   hclusPerDisk3 = fs->make<TH1F>("hclusPerDisk3", "Clus per disk3", sizeH, lowH, highH);
0198   hclusPerDisk4 = fs->make<TH1F>("hclusPerDisk4", "Clus per disk4", sizeH, lowH, highH);
0199 
0200   sizeH = 2000;
0201   highH = 9999.5;
0202   hpixPerLay1 = fs->make<TH1F>("hpixPerLay1", "Pix per layer l1", sizeH, lowH, highH);
0203   hpixPerLay2 = fs->make<TH1F>("hpixPerLay2", "Pix per layer l2", sizeH, lowH, highH);
0204   hpixPerLay3 = fs->make<TH1F>("hpixPerLay3", "Pix per layer l3", sizeH, lowH, highH);
0205 
0206   hclus = fs->make<TH1F>("hclus", "Clus per event", sizeH, lowH, highH);
0207   hdigis = fs->make<TH1F>("hdigis", "All Digis in clus per event", 2000, lowH, 19999.5);
0208   hdigisB = fs->make<TH1F>("hdigisB", "BPix Digis in clus per event", 8000, lowH, 19999.5);
0209   hdigisF = fs->make<TH1F>("hdigisF", "FPix Digis in clus per event", 2000, lowH, 7999.5);
0210   hclusBPix = fs->make<TH1F>("hclusBPix", "Bpix Clus per event", 2000, 0., 2000.);
0211   hclusFPix = fs->make<TH1F>("hclusFPix", "Fpix Clus per event", 2000, 0., 2000.);
0212   hdets = fs->make<TH1F>("hdets", "Dets per event", 2000, -0.5, 1999.5);
0213 
0214   hmaxPixPerDet = fs->make<TH1F>("hmaxPixPerDet", "Max pixels per det", 1000, -0.5, 999.5);
0215 
0216   sizeH = 1000;
0217   highH = 1999.5;
0218 
0219   hdetsPerLay1 = fs->make<TH1F>("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5);
0220   hdetsPerLay3 = fs->make<TH1F>("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5);
0221   hdetsPerLay2 = fs->make<TH1F>("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5);
0222 
0223   sizeH = 120;
0224   lowH = 0.;
0225   highH = 121.0;                                                             // charge limit in kelec
0226   hcharge1 = fs->make<TH1F>("hcharge1", "Clu charge l1", sizeH, 0., highH);  //in ke
0227   hcharge2 = fs->make<TH1F>("hcharge2", "Clu charge l2", sizeH, 0., highH);
0228   hcharge3 = fs->make<TH1F>("hcharge3", "Clu charge l3", sizeH, 0., highH);
0229   hcharge4 = fs->make<TH1F>("hcharge4", "Clu charge d1", sizeH, 0., highH);
0230   hcharge5 = fs->make<TH1F>("hcharge5", "Clu charge d2", sizeH, 0., highH);
0231 
0232   sizeH = 90;
0233   highH = 61.0;                                                                    // charge limit in kelec
0234   hpixcharge1 = fs->make<TH1F>("hpixcharge1", "Pix charge l1", sizeH, 0., highH);  //in ke
0235   hpixcharge2 = fs->make<TH1F>("hpixcharge2", "Pix charge l2", sizeH, 0., highH);
0236   hpixcharge3 = fs->make<TH1F>("hpixcharge3", "Pix charge l3", sizeH, 0., highH);
0237   hpixcharge4 = fs->make<TH1F>("hpixcharge4", "Pix charge d1", sizeH, 0., highH);
0238   hpixcharge5 = fs->make<TH1F>("hpixcharge5", "Pix charge d2", sizeH, 0., highH);
0239 
0240   hcols1 = fs->make<TH1F>("hcols1", "Layer 1 cols", 500, -0.5, 499.5);
0241   hcols2 = fs->make<TH1F>("hcols2", "Layer 2 cols", 500, -0.5, 499.5);
0242   hcols3 = fs->make<TH1F>("hcols3", "Layer 3 cols", 500, -0.5, 499.5);
0243 
0244   hrows1 = fs->make<TH1F>("hrows1", "Layer 1 rows", 200, -0.5, 199.5);
0245   hrows2 = fs->make<TH1F>("hrows2", "Layer 2 rows", 200, -0.5, 199.5);
0246   hrows3 = fs->make<TH1F>("hrows3", "layer 3 rows", 200, -0.5, 199.5);
0247 
0248   hpcols1 = fs->make<TH1F>("hpcols1", "Layer 1 pix cols", 500, -0.5, 499.5);
0249   hpcols2 = fs->make<TH1F>("hpcols2", "Layer 2 pix cols", 500, -0.5, 499.5);
0250   hpcols3 = fs->make<TH1F>("hpcols3", "Layer 3 pix cols", 500, -0.5, 499.5);
0251 
0252   hprows1 = fs->make<TH1F>("hprows1", "Layer 1 pix rows", 200, -0.5, 199.5);
0253   hprows2 = fs->make<TH1F>("hprows2", "Layer 2 pix rows", 200, -0.5, 199.5);
0254   hprows3 = fs->make<TH1F>("hprows3", "layer 3 pix rows", 200, -0.5, 199.5);
0255 
0256   sizeH = 1000;
0257   highH = 999.5;  // charge limit in kelec
0258   hsize1 = fs->make<TH1F>("hsize1", "layer 1 clu size", sizeH, -0.5, highH);
0259   hsize2 = fs->make<TH1F>("hsize2", "layer 2 clu size", sizeH, -0.5, highH);
0260   hsize3 = fs->make<TH1F>("hsize3", "layer 3 clu size", sizeH, -0.5, highH);
0261   hsize4 = fs->make<TH1F>("hsize4", "disk 1 clu size", sizeH, -0.5, highH);
0262   hsize5 = fs->make<TH1F>("hsize5", "disk 2 clu size", sizeH, -0.5, highH);
0263 
0264   hsizex1 = fs->make<TH1F>("hsizex1", "lay1 clu size in x", 10, -0.5, 9.5);
0265   hsizex2 = fs->make<TH1F>("hsizex2", "lay2 clu size in x", 10, -0.5, 9.5);
0266   hsizex3 = fs->make<TH1F>("hsizex3", "lay3 clu size in x", 10, -0.5, 9.5);
0267   hsizex4 = fs->make<TH1F>("hsizex4", "d1 clu size in x", 10, -0.5, 9.5);
0268   hsizex5 = fs->make<TH1F>("hsizex5", "d2 clu size in x", 10, -0.5, 9.5);
0269 
0270   hsizey1 = fs->make<TH1F>("hsizey1", "lay1 clu size in y", 20, -0.5, 19.5);
0271   hsizey2 = fs->make<TH1F>("hsizey2", "lay2 clu size in y", 20, -0.5, 19.5);
0272   hsizey3 = fs->make<TH1F>("hsizey3", "lay3 clu size in y", 20, -0.5, 19.5);
0273   hsizey4 = fs->make<TH1F>("hsizey4", "d1 clu size in y", 20, -0.5, 19.5);
0274   hsizey5 = fs->make<TH1F>("hsizey5", "d2 clu size in y", 20, -0.5, 19.5);
0275 
0276   hevent = fs->make<TH1F>("hevent", "event", 1000, 0, 10000000.);
0277   horbit = fs->make<TH1F>("horbit", "orbit", 100, 0, 100000000.);
0278 
0279   hlumi1 = fs->make<TH1F>("hlumi1", "lumi", 2000, 0, 2000.);
0280   hlumi0 = fs->make<TH1F>("hlumi0", "lumi", 2000, 0, 2000.);
0281   hlumi = fs->make<TH1F>("hlumi", "lumi", 2000, 0, 2000.);
0282   hbx1 = fs->make<TH1F>("hbx1", "bx", 4000, 0, 4000.);
0283   hbx0 = fs->make<TH1F>("hbx0", "bx", 4000, 0, 4000.);
0284   hbx = fs->make<TH1F>("hbx", "bx", 4000, 0, 4000.);
0285 
0286   hDetMap1 = fs->make<TH2F>("hDetMap1", " ", 9, -4.5, 4.5, 21, -10.5, 10.5);
0287   hDetMap1->SetOption("colz");
0288   hDetMap2 = fs->make<TH2F>("hDetMap2", " ", 9, -4.5, 4.5, 33, -16.5, 16.5);
0289   hDetMap2->SetOption("colz");
0290   hDetMap3 = fs->make<TH2F>("hDetMap3", " ", 9, -4.5, 4.5, 45, -22.5, 22.5);
0291   hDetMap3->SetOption("colz");
0292 
0293   hpDetMap1 = fs->make<TH2F>("hpDetMap1", " ", 9, -4.5, 4.5, 21, -10.5, 10.5);
0294   hpDetMap1->SetOption("colz");
0295   hpDetMap2 = fs->make<TH2F>("hpDetMap2", " ", 9, -4.5, 4.5, 33, -16.5, 16.5);
0296   hpDetMap2->SetOption("colz");
0297   hpDetMap3 = fs->make<TH2F>("hpDetMap3", " ", 9, -4.5, 4.5, 45, -22.5, 22.5);
0298   hpDetMap3->SetOption("colz");
0299 
0300   hpixDetMap1 = fs->make<TH2F>("hpixDetMap1", "pix det layer 1", 416, 0., 416., 160, 0., 160.);
0301   hpixDetMap2 = fs->make<TH2F>("hpixDetMap2", "pix det layer 2", 416, 0., 416., 160, 0., 160.);
0302   hpixDetMap3 = fs->make<TH2F>("hpixDetMap3", "pix det layer 3", 416, 0., 416., 160, 0., 160.);
0303 
0304   hcluDetMap1 = fs->make<TH2F>("hcluDetMap1", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0305   hcluDetMap2 = fs->make<TH2F>("hcluDetMap2", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0306   hcluDetMap3 = fs->make<TH2F>("hcluDetMap3", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0307   htest = fs->make<TH1F>("htest", "FPix R", 300, -15., 15.);
0308 
0309 #endif
0310 
0311   countEvents = 0;
0312   countAllEvents = 0;
0313   sumClusters = 0.;
0314 }
0315 // ------------ method called to at the end of the job  ------------
0316 void ReadPixClusters::endJob() {
0317   sumClusters = sumClusters / float(countEvents);
0318   edm::LogPrint("ReadPixClusters") << " End PixelClusTest, events all/with hits=  " << countAllEvents << "/"
0319                                    << countEvents << " " << sumClusters << " " << printLocal;
0320 }
0321 //////////////////////////////////////////////////////////////////
0322 // Functions that gets called by framework every event
0323 void ReadPixClusters::analyze(const edm::Event &e, const edm::EventSetup &es) {
0324   using namespace edm;
0325 
0326   // Get event setup
0327   edm::ESHandle<TrackerGeometry> geom = es.getHandle(trackerGeomToken_);
0328   const TrackerGeometry &theTracker(*geom);
0329 
0330 #ifdef NEW_ID
0331   //Retrieve tracker topology from geometry
0332   edm::ESHandle<TrackerTopology> tTopo = es.getHandle(trackerTopoToken_);
0333 #endif
0334 
0335   countAllEvents++;
0336   RunNumber_t const run = e.id().run();
0337   EventNumber_t const event = e.id().event();
0338   LuminosityBlockNumber_t const lumiBlock = e.luminosityBlock();
0339 
0340   int bx = e.bunchCrossing();
0341   int orbit = e.orbitNumber();
0342 
0343   // Get Cluster Collection from InputTag
0344   edm::Handle<edmNew::DetSetVector<SiPixelCluster>> clusters;
0345   e.getByToken(tPixelCluster, clusters);
0346 
0347   const edmNew::DetSetVector<SiPixelCluster> &input = *clusters;
0348   int numOf = input.size();
0349 
0350   hbx0->Fill(float(bx));
0351   hlumi0->Fill(float(lumiBlock));
0352 
0353   if (printLocal)
0354     edm::LogPrint("ReadPixClusters") << "run " << run << " event " << event << " bx " << bx << " lumi " << lumiBlock
0355                                      << " orbit " << orbit << " " << numOf;
0356 
0357   hdets->Fill(float(numOf));  // number of modules with pix
0358 
0359   // Select events with pixels
0360   //if(numOf<1) return; // skip events with  pixel dets
0361   //if(numOf<4) return; // skip events with few pixel dets
0362 
0363   hevent->Fill(float(event));
0364   hlumi->Fill(float(lumiBlock));
0365   hbx->Fill(float(bx));
0366   horbit->Fill(float(orbit));
0367 
0368   countEvents++;
0369   int numberOfDetUnits = 0;
0370   int numberOfClusters = 0;
0371   int numberOfPixels = 0;
0372   int numberOfDetUnits1 = 0;
0373   int numOfClustersPerDet1 = 0;
0374   int numOfClustersPerLay1 = 0;
0375   int numberOfDetUnits2 = 0;
0376   int numOfClustersPerDet2 = 0;
0377   int numOfClustersPerLay2 = 0;
0378   int numberOfDetUnits3 = 0;
0379   int numOfClustersPerDet3 = 0;
0380   int numOfClustersPerLay3 = 0;
0381 
0382   int numOfPixPerLay1 = 0;
0383   int numOfPixPerLay2 = 0;
0384   int numOfPixPerLay3 = 0;
0385 
0386   int numOfPixPerDet1 = 0;
0387   int numOfPixPerDet2 = 0;
0388   int numOfPixPerDet3 = 0;
0389 
0390   int numOfPixPerLink11 = 0;
0391   int numOfPixPerLink12 = 0;
0392   int numOfPixPerLink21 = 0;
0393   int numOfPixPerLink22 = 0;
0394   //SK:unused  int numOfPixPerLink3=0;
0395 
0396   int maxClusPerDet = 0;
0397   int maxPixPerDet = 0;
0398   unsigned int maxPixPerClu = 0;
0399 
0400   int numOfClustersPerDisk1 = 0;
0401   int numOfClustersPerDisk2 = 0;
0402   int numOfClustersPerDisk3 = 0;
0403   int numOfClustersPerDisk4 = 0;
0404   int numOfPixPerDisk1 = 0;
0405   int numOfPixPerDisk2 = 0;
0406   int numOfPixPerDisk3 = 0;
0407   int numOfPixPerDisk4 = 0;
0408 
0409   float aveCharge1 = 0., aveCharge2 = 0., aveCharge3 = 0., aveCharge4 = 0., aveCharge5 = 0.;
0410 
0411   static int module1[416][160] = {{0}};
0412   static int module2[416][160] = {{0}};
0413   static int module3[416][160] = {{0}};
0414 
0415   // get vector of detunit ids
0416   //--- Loop over detunits.
0417   edmNew::DetSetVector<SiPixelCluster>::const_iterator DSViter = input.begin();
0418   for (; DSViter != input.end(); DSViter++) {
0419     bool valid = false;
0420     unsigned int detid = DSViter->detId();
0421     // Det id
0422     DetId detId = DetId(detid);             // Get the Detid object
0423     unsigned int detType = detId.det();     // det type, pixel=1
0424     unsigned int subid = detId.subdetId();  //subdetector type, barrel=1
0425 
0426     if (printLocal)
0427       edm::LogPrint("ReadPixClusters") << "Det: " << detId.rawId() << " " << detId.null() << " " << detType << " "
0428                                        << subid;
0429 
0430 #ifdef HISTOS
0431       //hdetunit->Fill(float(detid));
0432       //hpixid->Fill(float(detType));
0433       //hpixsubid->Fill(float(subid));
0434 #endif  // HISTOS
0435 
0436     if (detType != 1)
0437       continue;  // look only at pixels
0438     ++numberOfDetUnits;
0439 
0440     //const GeomDetUnit * genericDet = geom->idToDet(detId);
0441     //const PixelGeomDetUnit * pixDet =
0442     //dynamic_cast<const PixelGeomDetUnit*>(genericDet);
0443 
0444     // Get the geom-detector
0445     const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
0446     double detZ = theGeomDet->surface().position().z();
0447     double detR = theGeomDet->surface().position().perp();
0448 
0449     //const BoundPlane& plane = theGeomDet->surface(); //for transf.
0450 
0451     //double detThick = theGeomDet->specificSurface().bounds().thickness();
0452     //int cols = theGeomDet->specificTopology().ncolumns();
0453     //int rows = theGeomDet->specificTopology().nrows();
0454 
0455     const PixelTopology *topol = &(theGeomDet->specificTopology());
0456 
0457     // barrel ids
0458     unsigned int layerC = 0;
0459     unsigned int ladderC = 0;
0460     unsigned int zindex = 0;
0461     int shell = 0;      // shell id // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0462     int sector = 0;     // 1-8
0463     int ladder = 0;     // 1-22
0464     int layer = 0;      // 1-3
0465     int module = 0;     // 1-4
0466     bool half = false;  //
0467 
0468     // Endcap ids
0469     unsigned int disk = 0;     //1,2,3
0470     unsigned int blade = 0;    //1-24
0471     unsigned int moduleF = 0;  // plaquette 1,2,3,4
0472     unsigned int side = 0;     //size=1 for -z, 2 for +z
0473     unsigned int panel = 0;    //panel=1
0474 
0475     edmNew::DetSet<SiPixelCluster>::const_iterator clustIt;
0476 
0477     // Subdet id, pix barrel=1, forward=2
0478     if (subid == 2) {  // forward
0479 
0480 #ifdef NEW_ID
0481       disk = tTopo->pxfDisk(detid);      //1,2,3
0482       blade = tTopo->pxfBlade(detid);    //1-24
0483       zindex = tTopo->pxfModule(detid);  //
0484       side = tTopo->pxfSide(detid);      //size=1 for -z, 2 for +z
0485       panel = tTopo->pxfPanel(detid);    //panel=1
0486 #else
0487       PXFDetId pdetId = PXFDetId(detid);
0488       disk = pdetId.disk();       //1,2,3
0489       blade = pdetId.blade();     //1-24
0490       moduleF = pdetId.module();  // plaquette
0491       side = pdetId.side();       //size=1 for -z, 2 for +z
0492       panel = pdetId.panel();     //panel=1
0493 #endif
0494 
0495       if (printLocal)
0496         edm::LogPrint("ReadPixClusters") << " forward det, disk " << disk << ", blade " << blade << ", module "
0497                                          << moduleF << ", side " << side << ", panel " << panel << " pos = " << detZ
0498                                          << " " << detR;
0499 
0500       bool fpixInner = (((panel == 1) && (moduleF <= 2)) || ((panel == 2) && (moduleF <= 1)));  // make split at 10cm
0501       if (fpixInner)
0502         htest->Fill(detR);
0503 
0504     } else if (subid == 1) {  // barrel
0505 
0506 #ifdef HISTOS
0507       //hdetr->Fill(detR);
0508       //hdetz->Fill(detZ);
0509 #endif  // HISTOS
0510 
0511       //hcolsB->Fill(float(cols));
0512       //hrowsB->Fill(float(rows));
0513 
0514 #ifdef NEW_ID
0515       layerC = tTopo->pxbLayer(detid);
0516       ladderC = tTopo->pxbLadder(detid);
0517       zindex = tTopo->pxbModule(detid);
0518       PixelBarrelName pbn(detid);
0519 #else
0520       PXBDetId pdetId = PXBDetId(detid);
0521       //unsigned int detTypeP=pdetId.det();
0522       //unsigned int subidP=pdetId.subdetId();
0523       // Barell layer = 1,2,3
0524       layerC = pdetId.layer();
0525       // Barrel ladder id 1-20,32,44.
0526       ladderC = pdetId.ladder();
0527       // Barrel Z-index=1,8
0528       zindex = pdetId.module();
0529       // Convert to online
0530       PixelBarrelName pbn(pdetId);
0531 #endif
0532 
0533       // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0534       PixelBarrelName::Shell sh = pbn.shell();  //enum
0535       sector = pbn.sectorName();
0536       ladder = pbn.ladderName();
0537       layer = pbn.layerName();
0538       module = pbn.moduleName();
0539       half = pbn.isHalfModule();
0540       shell = int(sh);
0541       // change the module sign for z<0
0542       if (shell == 1 || shell == 2)
0543         module = -module;
0544       // change ladeer sign for Outer )x<0)
0545       if (shell == 1 || shell == 3)
0546         ladder = -ladder;
0547 
0548       if (printLocal) {
0549         edm::LogPrint("ReadPixClusters") << " Barrel layer, ladder, module " << layerC << " " << ladderC << " "
0550                                          << zindex << " " << sh << "(" << shell << ") " << sector << " " << layer << " "
0551                                          << ladder << " " << module << " " << half;
0552         //edm::LogPrint("ReadPixClusters")<<" Barrel det, thick "<<detThick<<" "
0553         //  <<" layer, ladder, module "
0554         //  <<layer<<" "<<ladder<<" "<<zindex;
0555         //edm::LogPrint("ReadPixClusters")<<" col/row, pitch "<<cols<<" "<<rows<<" "
0556         //  <<pitchX<<" "<<pitchY;
0557       }
0558 
0559     }  // if subid
0560 
0561     if (printLocal) {
0562       edm::LogPrint("ReadPixClusters") << "List clusters : ";
0563       edm::LogPrint("ReadPixClusters") << "Num Charge Size SizeX SizeY X Y Xmin Xmax Ymin Ymax Edge";
0564     }
0565 
0566     // Loop over clusters
0567     for (clustIt = DSViter->begin(); clustIt != DSViter->end(); clustIt++) {
0568       sumClusters++;
0569       numberOfClusters++;
0570       float ch = float(clustIt->charge()) / 1000.;  // convert ke to electrons
0571       int size = clustIt->size();
0572       int sizeX = clustIt->sizeX();  //x=row=rfi, starting from CMSSW6 this is limited to 64
0573       int sizeY = clustIt->sizeY();  //y=col=z_global, starting from CMSSW6 this is limited to 64
0574       float x = clustIt->x();        // cluster position as float (int+0.5)
0575       float y = clustIt->y();        // analog average
0576 
0577       // Returns int index of the cluster min/max
0578       // This parameters are affected by the 64*64 limit on the cluster grid, so max-min<=64.
0579       int minPixelRow = clustIt->minPixelRow();  //x
0580       int maxPixelRow = clustIt->maxPixelRow();
0581       int minPixelCol = clustIt->minPixelCol();  //y
0582       int maxPixelCol = clustIt->maxPixelCol();
0583 
0584       //unsigned int geoId = clustIt->geographicalId(); // always 0?!
0585 
0586       // edge method moved to topologu class
0587       bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
0588       bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
0589 
0590       if (printLocal)
0591         edm::LogPrint("ReadPixClusters") << numberOfClusters << " " << ch << " " << size << " " << sizeX << " " << sizeY
0592                                          << " " << x << " " << y << " " << minPixelRow << " " << maxPixelRow << " "
0593                                          << minPixelCol << " " << maxPixelCol << " " << edgeHitX << " " << edgeHitY;
0594 
0595       // Get the pixels in the Cluster
0596       const vector<SiPixelCluster::Pixel> &pixelsVec = clustIt->pixels();
0597       if (printLocal)
0598         edm::LogPrint("ReadPixClusters") << " Pixels in this cluster ";
0599 
0600       bool bigInX = false, bigInY = false;
0601       bool edgeHitX2 = false;  // edge method moved
0602       bool edgeHitY2 = false;  // to topologu class
0603       //bool edgeInX = false; // edge method moved
0604       //bool edgeInY = false; // to topology class
0605       //SK:unused      bool cluBigInX = false; // does this clu include a big pixel
0606       //SK:unused      bool cluBigInY = false; // does this clu include a big pixel
0607       //int noisy = 0;
0608 
0609       if (pixelsVec.size() > maxPixPerClu)
0610         maxPixPerClu = pixelsVec.size();
0611 
0612       // Look at pixels in this cluster. ADC is calibrated, in electrons
0613       for (unsigned int i = 0; i < pixelsVec.size(); ++i) {  // loop over pixels
0614         numberOfPixels++;
0615         float pixx = pixelsVec[i].x;  // index as float=iteger, row index, 0-159
0616         float pixy = pixelsVec[i].y;  // same, col index, 0-415
0617         float adc = (float(pixelsVec[i].adc) / 1000.);
0618 
0619         //int chan = PixelChannelIdentifier::pixelToChannel(int(pixx),int(pixy));
0620         //bool binInX = (RectangularPixelTopology::isItBigPixelInX(int(pixx)));
0621         //bool bigInY = (RectangularPixelTopology::isItBigPixelInY(int(pixy)));
0622         //int roc = rocId(int(pixy),int(pixx));  // column, row
0623 
0624 #ifdef HISTOS
0625         // Pixel histos
0626         if (subid == 1) {  // barrel
0627           if (layer == 1) {
0628             numOfPixPerDet1++;
0629             numOfPixPerLay1++;
0630             valid = valid || true;
0631             hpixcharge1->Fill(adc);
0632             hpixDetMap1->Fill(pixy, pixx);
0633             hpDetMap1->Fill(float(module), float(ladder));
0634             module1[int(pixx)][int(pixy)]++;
0635 
0636             hpcols1->Fill(pixy);
0637             hprows1->Fill(pixx);
0638 
0639             if (pixx < 80.)
0640               numOfPixPerLink11++;
0641             else
0642               numOfPixPerLink12++;
0643 
0644           } else if (layer == 2) {
0645             numOfPixPerDet2++;
0646             numOfPixPerLay2++;
0647 
0648             hpcols2->Fill(pixy);
0649             hprows2->Fill(pixx);
0650 
0651             if (pixx < 80.)
0652               numOfPixPerLink21++;
0653             else
0654               numOfPixPerLink22++;
0655 
0656             hpixcharge2->Fill(adc);
0657             hpixDetMap2->Fill(pixy, pixx);
0658             hpDetMap2->Fill(float(module), float(ladder));
0659             module2[int(pixx)][int(pixy)]++;
0660 
0661           } else if (layer == 3) {
0662             numOfPixPerDet3++;
0663             numOfPixPerLay3++;
0664             valid = valid || true;
0665             hpixcharge3->Fill(adc);
0666             hpixDetMap3->Fill(pixy, pixx);
0667             hpDetMap3->Fill(float(module), float(ladder));
0668             module3[int(pixx)][int(pixy)]++;
0669 
0670             hpcols3->Fill(pixy);
0671             hprows3->Fill(pixx);
0672 
0673           }  // if layer
0674 
0675         } else if (subid == 2) {  // endcap
0676           // pixels
0677 
0678           if (disk == 1) {  // disk1 -+z
0679             if (side == 1)
0680               numOfPixPerDisk2++;  // d1,-z
0681             else if (side == 2)
0682               numOfPixPerDisk3++;  // d1, +z
0683             else
0684               edm::LogPrint("ReadPixClusters") << " unknown side " << side;
0685 
0686             hpixcharge4->Fill(adc);
0687 
0688           } else if (disk == 2) {  // disk2 -+z
0689 
0690             if (side == 1)
0691               numOfPixPerDisk1++;  // d2, -z
0692             else if (side == 2)
0693               numOfPixPerDisk4++;  // d2, +z
0694             else
0695               edm::LogPrint("ReadPixClusters") << " unknown side " << side;
0696 
0697             hpixcharge5->Fill(adc);
0698 
0699           } else
0700             edm::LogPrint("ReadPixClusters") << " unknown disk " << disk;
0701 
0702         }  // end if subdet (pixel loop)
0703 
0704 #endif  // HISTOS
0705 
0706         // This looks if there is an edge pixel in the cluster
0707         bool edgeInX = topol->isItEdgePixelInX(int(pixx));
0708         bool edgeInY = topol->isItEdgePixelInY(int(pixy));
0709         if (edgeInX)
0710           edgeHitX2 = true;
0711         if (edgeInY)
0712           edgeHitY2 = true;
0713 
0714         if (printLocal)
0715           edm::LogPrint("ReadPixClusters") << i << " " << pixx << " " << pixy << " " << adc << " " << bigInX << " "
0716                                            << bigInY << " " << edgeInX << " " << edgeInY;
0717 
0718         //SK:unused if(bigInX) cluBigInX=true;
0719         //SK:unused if(bigInY) cluBigInY=true;
0720 
0721       }  // pixel loop
0722 
0723       // if(edgeHitX ||  edgeHitX2 || edgeHitY ||  edgeHitY2)
0724       //    edm::LogPrint("ReadPixClusters")<<" egde "<<edgeHitX<<" "<<edgeHitX2<<" "<<edgeHitY<<" "<<edgeHitY2<<" "
0725       //        <<size<<" "<<sizeX<<" "<<sizeY<<" "
0726       //        <<x<<" "<<y<<" "<<minPixelRow<<" "<<maxPixelRow<<" "<<minPixelCol<<" "
0727       //        <<maxPixelCol<<" "<<edgeHitX<<" "<<edgeHitY;
0728 
0729       // This will happen  for clusters wider than 64
0730       if ((edgeHitX != edgeHitX2) && sizeX < 64)
0731         edm::LogPrint("ReadPixClusters") << " wrong egdeX " << edgeHitX << " "
0732                                          << edgeHitX2
0733                                          //<<" "<<event<<" "<<detid<<" "<<numberOfClusters
0734                                          << " size = " << size << " " << sizeX << " " << sizeY << " " << x << " " << y
0735                                          << " " << minPixelRow << " " << maxPixelRow << " " << minPixelCol << " "
0736                                          << maxPixelCol << " " << edgeHitX << " " << edgeHitY;
0737 
0738       if ((edgeHitY != edgeHitY2) && sizeY < 64)
0739         edm::LogPrint("ReadPixClusters") << " wrong egdeY " << edgeHitY << " "
0740                                          << edgeHitY2
0741                                          //<<" "<<event<<" "<<detid<<" "<<numberOfClusters
0742                                          << " size = " << size << " " << sizeX << " " << sizeY << " " << x << " " << y
0743                                          << " " << minPixelRow << " " << maxPixelRow << " " << minPixelCol << " "
0744                                          << maxPixelCol << " " << edgeHitX << " " << edgeHitY;
0745 
0746 #ifdef HISTOS
0747 
0748       // Cluster histos
0749       if (subid == 1) {  // barrel
0750         //if (subid==1) {  // barrel
0751         if (layer == 1) {  // layer
0752 
0753           hDetMap1->Fill(float(module), float(ladder));
0754           hcluDetMap1->Fill(y, x);
0755           hcharge1->Fill(ch);
0756           hcols1->Fill(y);
0757           hrows1->Fill(x);
0758           hsize1->Fill(float(size));
0759           hsizex1->Fill(float(sizeX));
0760           hsizey1->Fill(float(sizeY));
0761           numOfClustersPerDet1++;
0762           numOfClustersPerLay1++;
0763 
0764           aveCharge1 += ch;
0765 
0766         } else if (layer == 2) {
0767           hDetMap2->Fill(float(module), float(ladder));
0768           hcluDetMap2->Fill(y, x);
0769           hcharge2->Fill(ch);
0770           hcols2->Fill(y);
0771           hrows2->Fill(x);
0772           hsize2->Fill(float(size));
0773           hsizex2->Fill(float(sizeX));
0774           hsizey2->Fill(float(sizeY));
0775           numOfClustersPerDet2++;
0776           numOfClustersPerLay2++;
0777 
0778           aveCharge2 += ch;
0779 
0780         } else if (layer == 3) {
0781           hDetMap3->Fill(float(module), float(ladder));
0782           hcluDetMap3->Fill(y, x);
0783           hcharge3->Fill(ch);
0784           hcols3->Fill(y);
0785           hrows3->Fill(x);
0786           hsize3->Fill(float(size));
0787           hsizex3->Fill(float(sizeX));
0788           hsizey3->Fill(float(sizeY));
0789           numOfClustersPerDet3++;
0790           numOfClustersPerLay3++;
0791 
0792           aveCharge3 += ch;
0793 
0794         }  // end if bpix layer
0795 
0796       } else if (subid == 2) {  // endcap
0797 
0798         //edm::LogPrint("ReadPixClusters")<<disk<<" "<<side;
0799         if (disk == 1) {  // disk1 -+z
0800           if (side == 1)
0801             numOfClustersPerDisk2++;  // d1,-z
0802           else if (side == 2)
0803             numOfClustersPerDisk3++;  // d1, +z
0804           else
0805             edm::LogPrint("ReadPixClusters") << " unknown side " << side;
0806 
0807           hcharge4->Fill(ch);
0808           aveCharge4 += ch;
0809           hsize4->Fill(float(size));
0810           hsizex4->Fill(float(sizeX));
0811           hsizey4->Fill(float(sizeY));
0812 
0813         } else if (disk == 2) {  // disk2 -+z
0814 
0815           if (side == 1)
0816             numOfClustersPerDisk1++;  // d2, -z
0817           else if (side == 2)
0818             numOfClustersPerDisk4++;  // d2, +z
0819           else
0820             edm::LogPrint("ReadPixClusters") << " unknown side " << side;
0821 
0822           hcharge5->Fill(ch);
0823           aveCharge5 += ch;
0824           hsize5->Fill(float(size));
0825           hsizex5->Fill(float(sizeX));
0826           hsizey5->Fill(float(sizeY));
0827 
0828         } else
0829           edm::LogPrint("ReadPixClusters") << " unknown disk " << disk;  // end fpix disk
0830 
0831       }  // end if barrel/forward
0832 
0833 #endif  // HISTOS
0834 
0835     }  // clusters
0836 
0837     if (numOfClustersPerDet1 > maxClusPerDet)
0838       maxClusPerDet = numOfClustersPerDet1;
0839     if (numOfClustersPerDet2 > maxClusPerDet)
0840       maxClusPerDet = numOfClustersPerDet2;
0841     if (numOfClustersPerDet3 > maxClusPerDet)
0842       maxClusPerDet = numOfClustersPerDet3;
0843 
0844     if (printLocal) {
0845       if (layer == 1)
0846         edm::LogPrint("ReadPixClusters") << "Lay1: number of clusters per det = " << numOfClustersPerDet1;
0847       else if (layer == 2)
0848         edm::LogPrint("ReadPixClusters") << "Lay2: number of clusters per det = " << numOfClustersPerDet1;
0849       else if (layer == 3)
0850         edm::LogPrint("ReadPixClusters") << "Lay3: number of clusters per det = " << numOfClustersPerDet1;
0851     }  // end if printLocal
0852 
0853 #ifdef HISTOS
0854     if (subid == 1) {  // barrel
0855 
0856       // Det histos
0857       if (layer == 1) {
0858         hladder1id->Fill(float(ladder));
0859         hz1id->Fill(float(module));
0860         ++numberOfDetUnits1;
0861         hclusPerDet1->Fill(float(numOfClustersPerDet1));
0862         hpixPerDet1->Fill(float(numOfPixPerDet1));
0863         if (numOfPixPerDet1 > maxPixPerDet)
0864           maxPixPerDet = numOfPixPerDet1;
0865         numOfClustersPerDet1 = 0;
0866         numOfPixPerDet1 = 0;
0867         //if(numOfPixPerLink11>798 || numOfPixPerLink12>798) select=true;
0868         hpixPerLink1->Fill(float(numOfPixPerLink11));
0869         hpixPerLink1->Fill(float(numOfPixPerLink12));
0870         numOfPixPerLink11 = 0;
0871         numOfPixPerLink12 = 0;
0872 
0873       } else if (layer == 2) {
0874         hladder2id->Fill(float(ladder));
0875         hz2id->Fill(float(module));
0876         ++numberOfDetUnits2;
0877         hclusPerDet2->Fill(float(numOfClustersPerDet2));
0878         hpixPerDet2->Fill(float(numOfPixPerDet2));
0879         if (numOfPixPerDet2 > maxPixPerDet)
0880           maxPixPerDet = numOfPixPerDet2;
0881         numOfClustersPerDet2 = 0;
0882         numOfPixPerDet2 = 0;
0883         hpixPerLink2->Fill(float(numOfPixPerLink21));
0884         hpixPerLink2->Fill(float(numOfPixPerLink22));
0885         numOfPixPerLink21 = 0;
0886         numOfPixPerLink22 = 0;
0887 
0888       } else if (layer == 3) {
0889         hladder3id->Fill(float(ladder));
0890         hz3id->Fill(float(module));
0891         ++numberOfDetUnits3;
0892         hclusPerDet3->Fill(float(numOfClustersPerDet3));
0893         hpixPerDet3->Fill(float(numOfPixPerDet3));
0894         if (numOfPixPerDet3 > maxPixPerDet)
0895           maxPixPerDet = numOfPixPerDet3;
0896         numOfClustersPerDet3 = 0;
0897         numOfPixPerDet3 = 0;
0898         //SK:unused numOfPixPerLink3=0;
0899 
0900       }  // layer
0901 
0902     }  // end barrel/forward
0903 
0904 #endif  // HISTOS
0905 
0906   }  // detunits loop
0907 
0908   if (0) {
0909     edm::LogPrint("ReadPixClusters") << "run " << run << " event " << event << " bx " << bx << " lumi " << lumiBlock
0910                                      << " orbit " << orbit << " num " << countEvents;
0911     edm::LogPrint("ReadPixClusters") << "Num of pix " << numberOfPixels << " num of clus " << numberOfClusters
0912                                      << " num of dets " << numOf << " max clus per det " << maxClusPerDet
0913                                      << " max pix per clu " << maxPixPerClu << " count " << countEvents;
0914     edm::LogPrint("ReadPixClusters") << "Number of clusters per      Lay1,2,3: " << numOfClustersPerLay1 << " "
0915                                      << numOfClustersPerLay2 << " " << numOfClustersPerLay3;
0916     edm::LogPrint("ReadPixClusters") << "Number of pixels per        Lay1,2,3: " << numOfPixPerLay1 << " "
0917                                      << numOfPixPerLay2 << " " << numOfPixPerLay3;
0918     edm::LogPrint("ReadPixClusters") << "Number of dets with clus in Lay1,2,3: " << numberOfDetUnits1 << " "
0919                                      << numberOfDetUnits2 << " " << numberOfDetUnits3;
0920     edm::LogPrint("ReadPixClusters") << "Number of clus in disks 1,2,3,4     : " << numOfClustersPerDisk1 << " "
0921                                      << numOfClustersPerDisk2 << " " << numOfClustersPerDisk3 << " "
0922                                      << numOfClustersPerDisk4;
0923     aveCharge1 /= float(numOfClustersPerLay1);
0924     aveCharge2 /= float(numOfClustersPerLay2);
0925     aveCharge3 /= float(numOfClustersPerLay3);
0926     aveCharge4 /= float(numOfClustersPerDisk2 + numOfClustersPerDisk3);
0927     aveCharge5 /= float(numOfClustersPerDisk1 + numOfClustersPerDisk4);
0928     edm::LogPrint("ReadPixClusters") << " Average charge " << aveCharge1 << " " << aveCharge2 << " " << aveCharge3
0929                                      << " " << aveCharge4 << " " << aveCharge5;
0930   }
0931 
0932 #ifdef HISTOS
0933 
0934   hdigis->Fill(float(numberOfPixels));  // all pix
0935   int tmp1 = numOfPixPerLay1 + numOfPixPerLay2 + numOfPixPerLay3;
0936   hdigisB->Fill(float(tmp1));  // pix in bpix
0937   tmp1 = numOfPixPerDisk1 + numOfPixPerDisk2 + numOfPixPerDisk3 + numOfPixPerDisk4;
0938   hdigisF->Fill(float(tmp1));  // pix in fpix
0939 
0940   hclus->Fill(float(numberOfClusters));  // clusters fpix+bpix
0941   tmp1 = numOfClustersPerLay1 + numOfClustersPerLay2 + numOfClustersPerLay3;
0942   hclusBPix->Fill(float(tmp1));  // clusters in bpix
0943   int tmp2 = numOfClustersPerDisk1 + numOfClustersPerDisk2 + numOfClustersPerDisk3 + numOfClustersPerDisk4;
0944   hclusFPix->Fill(float(tmp2));  // clusters in fpix
0945 
0946   hclusPerLay1->Fill(float(numOfClustersPerLay1));
0947   hclusPerLay2->Fill(float(numOfClustersPerLay2));
0948   hclusPerLay3->Fill(float(numOfClustersPerLay3));
0949   hpixPerLay1->Fill(float(numOfPixPerLay1));
0950   hpixPerLay2->Fill(float(numOfPixPerLay2));
0951   hpixPerLay3->Fill(float(numOfPixPerLay3));
0952   if (numOfClustersPerLay1 > 0)
0953     hdetsPerLay1->Fill(float(numberOfDetUnits1));
0954   if (numOfClustersPerLay2 > 0)
0955     hdetsPerLay2->Fill(float(numberOfDetUnits2));
0956   if (numOfClustersPerLay3 > 0)
0957     hdetsPerLay3->Fill(float(numberOfDetUnits3));
0958 
0959   hclusPerDisk1->Fill(float(numOfClustersPerDisk1));
0960   hclusPerDisk2->Fill(float(numOfClustersPerDisk2));
0961   hclusPerDisk3->Fill(float(numOfClustersPerDisk3));
0962   hclusPerDisk4->Fill(float(numOfClustersPerDisk4));
0963 
0964   hmaxPixPerDet->Fill(float(maxPixPerDet));
0965 
0966   // Select only events with more tha 3 clusters
0967   if (numberOfClusters > 3) {
0968     hbx1->Fill(float(bx));
0969     hlumi1->Fill(float(lumiBlock));
0970   }  // if num of clusters
0971 
0972 #endif  // HISTOS
0973 
0974 }  // end
0975 
0976 //define this as a plug-in
0977 DEFINE_FWK_MODULE(ReadPixClusters);