SiPixelDigitizer/test/PixelSimHitsTest.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PixelSimHitsTest
0004 // Class:      PixelSimHitsTest
0005 //
0006 /**\class PixelSimHitsTest PixelSimHitsTest.cc
0007
0008  Description: Test pixel simhits. Barrel only. Uses root histos.
0009  Modifed for module() method in PXBDetId. 2/06
0010  Add global coordiantes. 21/2/06
0011  Update 28/2/12. Works with CMSSW_5_3_8, not in cvs, d.k.
0012  Needs updating for 6.2 (in cvs)
0013 New det-id.
0014  Implementation:
0015      <Notes on implementation>
0016 */
0017 //
0018 // Original Author:  d.k.
0019 //         Created:  Jan CET 2006
0020 // $Id: PixelSimHitsTest.cc,v 1.8 2009/11/13 14:14:23 fambrogl Exp $
0021 //
0022 //
0023 // system include files
0024 #include <iostream>
0025
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "DataFormats/Common/interface/Handle.h"
0033
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/ServiceRegistry/interface/Service.h"
0036 #include "FWCore/Framework/interface/ESHandle.h"
0037
0038 // my includes
0039 //#include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0043 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0044 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"  //
0045
0046 // for det id
0047 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0048 #include "DataFormats/TrackerCommon/interface/PixelBarrelName.h"
0049
0050 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0051
0052 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0053 //#include "Geometry/Surface/interface/Surface.h"
0054
0055 // To use root histos
0056 #include "FWCore/ServiceRegistry/interface/Service.h"
0057 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0058
0059 // For ROOT
0060 #include <TROOT.h>
0061 #include <TChain.h>
0062 #include <TFile.h>
0063 #include <TF1.h>
0064 #include <TH2F.h>
0065 #include <TH1F.h>
0066 #include <TProfile.h>
0067
0068 using namespace std;
0069 using namespace edm;
0070
0071 //#define CHECK_GEOM
0072
0073 class PixelSimHitsTest : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0074 public:
0075   explicit PixelSimHitsTest(const edm::ParameterSet &);
0076   ~PixelSimHitsTest() override;
0077   void beginJob() override;
0078   void analyze(const edm::Event &, const edm::EventSetup &) override;
0079   void endJob() override;
0080
0081 private:
0082   // ----------member data ---------------------------
0083   edm::ParameterSet conf_;
0084   bool PRINT;
0085   string mode_;  // select bpix/fpix
0086   edm::EDGetTokenT<PSimHitContainer> tPixelSimHits;
0087   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken;
0088   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken;
0089   int numEvents;
0090   double numSimHits, numSimHitsGood;
0091
0092   TH1F *heloss1, *heloss2, *heloss3, *hdetunit, *hpabs, *hpid, *htof, *htid;
0093   TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id;
0094   //TH1F *hladder1idUp, *hladder2idUp, *hladder3idUp;
0095   TH1F *hthick1, *hthick2, *hthick3, *hlength1, *hlength2, *hlength3;
0096   TH1F *hwidth1, *hwidth2, *hwidth3;
0097   TH1F *hsimHitsPerDet1, *hsimHitsPerDet2, *hsimHitsPerDet3;
0098   TH1F *hsimHitsPerLay1, *hsimHitsPerLay2, *hsimHitsPerLay3;
0099   TH1F *hsimHits, *hsimHitsGood;
0100   TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3;
0101   TH1F *heloss1e, *heloss2e, *heloss3e;
0102   TH1F *heloss1mu, *heloss2mu, *heloss3mu;
0103   TH1F *htheta1, *htheta2, *htheta3;
0104   TH1F *hphi1, *hphi2, *hphi3;
0105   TH1F *hdetr, *hdetz, *hdetphi1, *hdetphi2, *hdetphi3;
0106   TH1F *hglobr1, *hglobr2, *hglobr3, *hglobz1, *hglobz2, *hglobz3;
0107   TH1F *hcolsB, *hrowsB, *hcolsF, *hrowsF;
0108
0109   TH2F *htest, *htest2, *htest3, *htest4, *htest5;
0110   //TProfile *hp1, *hp2, *hp3, *hp4, *hp5;
0111
0112 #ifdef CHECK_GEOM
0113   float modulePositionZ[3][44][8];
0114   float modulePositionR[3][44][8];
0115   float modulePositionPhi[3][44][8];
0116 #endif
0117 };
0118 //
0119 PixelSimHitsTest::PixelSimHitsTest(const edm::ParameterSet &iConfig) : conf_(iConfig) {
0120   usesResource(TFileService::kSharedResource);
0121
0122   std::string src_, list_;
0123   src_ = iConfig.getParameter<std::string>("src");
0124   list_ = iConfig.getParameter<std::string>("list");
0125
0126   edm::InputTag tag(src_, list_);  // for the ByToken
0127   tPixelSimHits = consumes<PSimHitContainer>(tag);
0128
0129   trackerTopoToken = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0130   trackerGeomToken = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0131
0132   mode_ = iConfig.getUntrackedParameter<std::string>("mode", "bpix");  // select bpix or fpix
0133   PRINT = iConfig.getUntrackedParameter<bool>("Verbosity", false);     // printout
0134   cout << " Construct PixelSimHitsTest " << endl;
0135 }
0136
0137 PixelSimHitsTest::~PixelSimHitsTest() {
0138   // do anything here that needs to be done at desctruction time
0139   // (e.g. close files, deallocate resources etc.)
0140   cout << " Destroy PixelSimHitsTest " << endl;
0141 }
0142
0143 // ------------ method called at the begining   ------------
0144 void PixelSimHitsTest::beginJob() {
0145   using namespace edm;
0146   cout << "Initialize PixelSimHitsTest " << endl;
0147
0148   // put here whatever you want to do at the beginning of the job
0149   numEvents = 0;
0150   numSimHits = 0.0;
0151   numSimHitsGood = 0.0;
0152
0153   edm::Service<TFileService> fs;
0154
0155   const float max_charge = 200.;  // in ke
0156   heloss1 = fs->make<TH1F>("heloss1", "Eloss l1", 100, 0., max_charge);
0157   heloss2 = fs->make<TH1F>("heloss2", "Eloss l2", 100, 0., max_charge);
0158   heloss3 = fs->make<TH1F>("heloss3", "Eloss l3", 100, 0., max_charge);
0159
0160   hdetunit = fs->make<TH1F>("hdetunit", "Det unit", 1000, 302000000., 302300000.);
0161   hpabs = fs->make<TH1F>("hpabs", "Pabs", 100, 0., 110.);
0162   htof = fs->make<TH1F>("htof", "TOF", 50, -25., 25.);
0163   hpid = fs->make<TH1F>("hpid", "PID", 1000, 0., 1000.);
0164   htid = fs->make<TH1F>("htid", "Track id", 100, 0., 100.);
0165
0166   hpixid = fs->make<TH1F>("hpixid", "Pix det id", 10, 0., 10.);
0167   hpixsubid = fs->make<TH1F>("hpixsubid", "Pix Barrel id", 10, 0., 10.);
0168   hlayerid = fs->make<TH1F>("hlayerid", "Pix layer id", 10, 0., 10.);
0169   hladder1id = fs->make<TH1F>("hladder1id", "Ladder L1 id", 102, -25.5, 25.5);
0170   hladder2id = fs->make<TH1F>("hladder2id", "Ladder L2 id", 102, -25.5, 25.5);
0171   hladder3id = fs->make<TH1F>("hladder3id", "Ladder L3 id", 102, -25.5, 25.5);
0172   hz1id = fs->make<TH1F>("hz1id", "Z-index id L1", 10, -5., 5.);
0173   hz2id = fs->make<TH1F>("hz2id", "Z-index id L2", 10, -5., 5.);
0174   hz3id = fs->make<TH1F>("hz3id", "Z-index id L3", 10, -5., 5.);
0175
0176   hthick1 = fs->make<TH1F>("hthick1", "Det 1 Thinckess", 400, 0., 0.04);
0177   hthick2 = fs->make<TH1F>("hthick2", "Det 2 Thinckess", 400, 0., 0.04);
0178   hthick3 = fs->make<TH1F>("hthick3", "Det 3 Thinckess", 400, 0., 0.04);
0179
0180   hlength1 = fs->make<TH1F>("hlength1", "Det 1 Length", 700, -3.5, 3.5);
0181   hlength2 = fs->make<TH1F>("hlength2", "Det 2 Length", 700, -3.5, 3.5);
0182   hlength3 = fs->make<TH1F>("hlength3", "Det 3 Length", 700, -3.5, 3.5);
0183
0184   hwidth1 = fs->make<TH1F>("hwidth1", "Det 1 Width", 200, -1., 1.);
0185   hwidth2 = fs->make<TH1F>("hwidth2", "Det 2 Width", 200, -1., 1.);
0186   hwidth3 = fs->make<TH1F>("hwidth3", "Det 3 Width", 200, -1., 1.);
0187
0188   hsimHitsPerDet1 = fs->make<TH1F>("hsimHitsPerDet1", "SimHits per det l1", 200, -0.5, 199.5);
0189   hsimHitsPerDet2 = fs->make<TH1F>("hsimHitsPerDet2", "SimHits per det l2", 200, -0.5, 199.5);
0190   hsimHitsPerDet3 = fs->make<TH1F>("hsimHitsPerDet3", "SimHits per det l3", 200, -0.5, 199.5);
0191   hsimHitsPerLay1 = fs->make<TH1F>("hsimHitsPerLay1", "SimHits per layer l1", 2000, -0.5, 1999.5);
0192   hsimHitsPerLay2 = fs->make<TH1F>("hsimHitsPerLay2", "SimHits per layer l2", 2000, -0.5, 1999.5);
0193   hsimHitsPerLay3 = fs->make<TH1F>("hsimHitsPerLay3", "SimHits per layer l3", 2000, -0.5, 1999.5);
0194   hdetsPerLay1 = fs->make<TH1F>("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5);
0195   hdetsPerLay3 = fs->make<TH1F>("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5);
0196   hdetsPerLay2 = fs->make<TH1F>("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5);
0197   hsimHits = fs->make<TH1F>("hsimHits", "SimHits for bpix", 2000, -0.5, 1999.5);
0198   hsimHitsGood = fs->make<TH1F>("hsimHitsGood", "SimHits for bpix", 2000, -0.5, 1999.5);  // no deltas
0199
0200   heloss1e = fs->make<TH1F>("heloss1e", "Eloss e l1", 100, 0., max_charge);
0201   heloss2e = fs->make<TH1F>("heloss2e", "Eloss e l2", 100, 0., max_charge);
0202   heloss3e = fs->make<TH1F>("heloss3e", "Eloss e l3", 100, 0., max_charge);
0203
0204   heloss1mu = fs->make<TH1F>("heloss1mu", "Eloss mu l1", 100, 0., max_charge);
0205   heloss2mu = fs->make<TH1F>("heloss2mu", "Eloss mu l2", 100, 0., max_charge);
0206   heloss3mu = fs->make<TH1F>("heloss3mu", "Eloss mu l3", 100, 0., max_charge);
0207
0208   htheta1 = fs->make<TH1F>("htheta1", "Theta l1", 350, 0.0, 3.5);
0209   htheta2 = fs->make<TH1F>("htheta2", "Theta l2", 350, 0.0, 3.5);
0210   htheta3 = fs->make<TH1F>("htheta3", "Theta l3", 350, 0.0, 3.5);
0211   hphi1 = fs->make<TH1F>("hphi1", "phi l1", 1400, -3.5, 3.5);
0212   hphi2 = fs->make<TH1F>("hphi2", "phi l2", 1400, -3.5, 3.5);
0213   hphi3 = fs->make<TH1F>("hphi3", "phi l3", 1400, -3.5, 3.5);
0214
0215   hdetr = fs->make<TH1F>("hdetr", "det r", 1500, 0., 15.);
0216   hdetz = fs->make<TH1F>("hdetz", "det z", 5200, -26., 26.);
0217
0218   hdetphi1 = fs->make<TH1F>("hdetphi1", "det phi l1", 700, -3.5, 3.5);
0219   hdetphi2 = fs->make<TH1F>("hdetphi2", "det phi l2", 700, -3.5, 3.5);
0220   hdetphi3 = fs->make<TH1F>("hdetphi3", "det phi l3", 700, -3.5, 3.5);
0221
0222   hcolsB = fs->make<TH1F>("hcolsB", "cols per bar det", 450, 0., 450.);
0223   hrowsB = fs->make<TH1F>("hrowsB", "rows per bar det", 200, 0., 200.);
0224   hcolsF = fs->make<TH1F>("hcolsF", "cols per for det", 300, 0., 300.);
0225   hrowsF = fs->make<TH1F>("hrowsF", "rows per for det", 200, 0., 200.);
0226
0227   //hladder1idUp = fs->make<TH1F>( "hladder1idUp", "Ladder L1 id", 100, -0.5, 49.5);
0228   //hladder2idUp = fs->make<TH1F>( "hladder2idUp", "Ladder L2 id", 100, -0.5, 49.5);
0229   //hladder3idUp = fs->make<TH1F>( "hladder3idUp", "Ladder L3 id", 100, -0.5, 49.5);
0230
0231   hglobr1 = fs->make<TH1F>("hglobr1", "global r1", 150, 0., 15.);
0232   hglobz1 = fs->make<TH1F>("hglobz1", "global z1", 540, -27., 27.);
0233   hglobr2 = fs->make<TH1F>("hglobr2", "global r2", 150, 0., 15.);
0234   hglobz2 = fs->make<TH1F>("hglobz2", "global z2", 540, -27., 27.);
0235   hglobr3 = fs->make<TH1F>("hglobr3", "global r3", 150, 0., 15.);
0236   hglobz3 = fs->make<TH1F>("hglobz3", "global z3", 540, -27., 27.);
0237
0238   // layer 1 only
0239   htest = fs->make<TH2F>("htest", " ", 108, -27., 27., 35, -3.5, 3.5);     // global z versus local y
0240   htest2 = fs->make<TH2F>("htest2", " ", 108, -27., 27., 60, 0., 600.);    // global z versus eloss
0241   htest3 = fs->make<TH2F>("htest3", " ", 240, -12., 12., 240, -12., 12.);  // x-y plane
0242   //htest4 = fs->make<TH2F>("htest4"," ",80,-4.,4.,100,-5.,5.);
0243
0244   //hp1 = fs->make<TProfile>("hp1"," ",50,0.,50.);    // default option
0245   //hp2 = fs->make<TProfile>("hp2"," ",50,0.,50.," "); // option set to " "
0246   //hp3 = fs->make<TProfile>("hp3"," ",50,0.,50.,-100.,100.); //
0247
0248 #ifdef CHECK_GEOM
0249   // To get the module position
0250   for (int i = 0; i < 3; i++) {
0251     for (int n = 0; n < 44; n++) {
0252       for (int m = 0; m < 8; m++) {
0253         modulePositionR[i][n][m] = -1;
0254         modulePositionZ[i][n][m] = -1;
0255         modulePositionPhi[i][n][m] = -1;
0256       }
0257     }
0258   }
0259 #endif
0260 }
0261
0262 // ------------ method called to produce the data  ------------
0263 void PixelSimHitsTest::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0264   const bool DEBUG = false;
0265   //string mode = "bpix";
0266
0267   using namespace edm;
0268   if (PRINT)
0269     cout << " Analyze PixelSimHitsTest " << endl;
0270
0271   // Get event setup (to get global transformation)
0272   edm::ESHandle<TrackerGeometry> geom = iSetup.getHandle(trackerGeomToken);
0273   const TrackerGeometry &theTracker(*geom);
0274
0275   //Retrieve tracker topology from geometry (for det id)
0276   edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(trackerTopoToken);
0277
0278   // Get input data
0279   int totalNumOfSimHits = 0;
0280   int totalNumOfSimHits1 = 0;
0281   int totalNumOfSimHits2 = 0;
0282   int totalNumOfSimHits3 = 0;
0283   int goodHits = 0;  // above pt=0.1GeV
0284
0285   // To count simhits per det module
0286   //typedef std::map<unsigned int, std::vector<PSimHit>,
0287   //std::less<unsigned int>> simhit_map;
0288   //typedef simhit_map::iterator simhit_map_iterator;
0289   //simhit_map SimHitMap;
0290   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap1;
0291   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap2;
0292   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap3;
0293
0294   // Access hits containers
0295   if (DEBUG)
0296     cout << "Define simhit container" << endl;
0297   Handle<PSimHitContainer> PixelHits;
0298   //Handle<PSimHitContainer> PixelHitsLowTof;
0299   //Handle<PSimHitContainer> PixelHitsHighTof;
0300   iEvent.getByToken(tPixelSimHits, PixelHits);
0301
0302   //if(mode=="bpix") {
0303   //iEvent.getByLabel( src_ ,"TrackerHitsPixelBarrelLowTof" ,PixelHitsLowTof);
0304   //iEvent.getByLabel( src_ ,"TrackerHitsPixelBarrelHighTof",PixelHitsHighTof);
0305   //} else if(mode=="fpix") {
0306   //iEvent.getByLabel( src_ ,"TrackerHitsPixelEndcapLowTof",PixelHitsLowTof);
0307   //iEvent.getByLabel( src_ ,"TrackerHitsPixelEndcapHighTof",PixelHitsHighTof);
0308   //}
0309
0310   if (DEBUG)
0311     cout << "Loop over SimHits LowTof" << endl;
0312   //for(vector<PSimHit>::const_iterator isim = PixelHitsLowTof->begin();
0313   // isim != PixelHitsLowTof->end(); ++isim) {
0314   for (vector<PSimHit>::const_iterator isim = PixelHits->begin(); isim != PixelHits->end(); ++isim) {
0315     totalNumOfSimHits++;
0316     // Det id
0317     DetId detId = DetId((*isim).detUnitId());
0318     unsigned int dettype = detId.det();     // for pixel=1
0319     unsigned int subid = detId.subdetId();  // barrel=1
0320     unsigned int detid = detId.rawId();     // raw det id
0321
0322     if (dettype != 1 && subid != 1)
0323       cout << " error in det id " << dettype << " " << subid << endl;
0324     if (PRINT)
0325       cout << totalNumOfSimHits << " det id " << detid << " " << dettype << " " << subid << endl;
0326     if (DEBUG)
0327       cout << " det unit " << (*isim).detUnitId() << detId.null() << endl;
0328
0329     // Global variables
0330     const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
0331     double detZ = theGeomDet->surface().position().z();      // module z position
0332     double detR = theGeomDet->surface().position().perp();   //        r
0333     double detPhi = theGeomDet->surface().position().phi();  //        phi
0334     hdetr->Fill(detR);
0335     hdetz->Fill(detZ);
0336
0337     double detThick = theGeomDet->specificSurface().bounds().thickness();
0338     double detLength = theGeomDet->specificSurface().bounds().length();
0339     double detWidth = theGeomDet->specificSurface().bounds().width();
0340
0341     int cols = theGeomDet->specificTopology().ncolumns();
0342     int rows = theGeomDet->specificTopology().nrows();
0343
0344     if (DEBUG)
0345       cout << "det z/r " << detZ << "/" << detR << " thick/len/wid " << detThick << " " << detLength << " " << detWidth
0346            << " cols/rows " << cols << " " << rows << endl;
0347
0348     unsigned int layerC = 0;
0349     unsigned int ladderC = 0;
0350     unsigned int zindex = 0;
0351
0352     // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0353     int shell = 0;      // shell id
0354     int sector = 0;     // 1-8
0355     int ladder = 0;     // 1-22
0356     int layer = 0;      // 1-3
0357     int module = 0;     // 1-4
0358     bool half = false;  //
0359
0360     unsigned int disk = 0;   //1,2
0361     unsigned int blade = 0;  //1-24
0362     unsigned int side = 0;   //size=1 for -z, 2 for +z
0363     unsigned int panel = 0;  //panel=1
0364
0365     if (mode_ == "fpix") {
0366       disk = tTopo->pxfDisk(detid);      //1,2,3
0367       blade = tTopo->pxfBlade(detid);    //1-24
0368       zindex = tTopo->pxfModule(detid);  //
0369       side = tTopo->pxfSide(detid);      //size=1 for -z, 2 for +z
0370       panel = tTopo->pxfPanel(detid);    //panel=1
0371
0372       if (PRINT) {
0373         cout << "Forward det " << subid << ", disk " << disk << ", blade " << blade << ", module " << zindex
0374              << ", side " << side << ", panel " << panel << endl;
0375         //cout<<" col/row, pitch "<<cols<<" "<<rows<<" "
0376         //<<pitchX<<" "<<pitchY<<endl;
0377       }
0378
0379       hcolsF->Fill(float(cols));
0380       hrowsF->Fill(float(rows));
0381
0382     } else if (mode_ == "bpix") {  // Barrel
0383
0384       // Barell layer = 1,2,3
0385       layerC = tTopo->pxbLayer(detid);
0386       // Barrel ladder id 1-20,32,44.
0387       ladderC = tTopo->pxbLadder(detid);
0388       // Barrel Z-index=1,8
0389       zindex = tTopo->pxbModule(detid);
0390       // Convert to online
0391       PixelBarrelName pbn(detid);
0392       // Shell { mO = 1, mI = 2 , pO =3 , pI =4 };
0393       PixelBarrelName::Shell sh = pbn.shell();  //enum
0394       sector = pbn.sectorName();
0395       ladder = pbn.ladderName();
0396       layer = pbn.layerName();
0397       module = pbn.moduleName();
0398       half = pbn.isHalfModule();
0399       shell = int(sh);
0400       // change the module sign for z<0
0401       if (shell == 1 || shell == 2)
0402         module = -module;
0403       // change ladeer sign for Outer )x<0)
0404       if (shell == 1 || shell == 3)
0405         ladder = -ladder;
0406
0407       if (PRINT) {
0408         cout << " Barrel layer, ladder, module " << layerC << " " << ladderC << " " << zindex << " " << sh << "("
0409              << shell << ") " << sector << " " << layer << " " << ladder << " " << module << " " << half << endl;
0410         //cout<<" Barrel det, thick "<<detThick<<" "
0411         //  <<" layer, ladder, module "
0412         //  <<layer<<" "<<ladder<<" "<<zindex<<endl;
0413         //cout<<" col/row, pitch "<<cols<<" "<<rows<<" "
0414         //  <<pitchX<<" "<<pitchY<<endl;
0415       }
0416
0417       hcolsB->Fill(float(cols));
0418       hrowsB->Fill(float(rows));
0419
0420     }  // end fb-bar
0421
0422 #ifdef CHECK_GEOM
0423     // To get the module position
0424     modulePositionR[layer - 1][ladderC - 1][zindex - 1] = detR;
0425     modulePositionZ[layer - 1][ladderC - 1][zindex - 1] = detZ;
0426     modulePositionPhi[layer - 1][ladderC - 1][zindex - 1] = detPhi;
0427 #endif
0428
0429     // SimHit information
0430     float eloss = (*isim).energyLoss() * 1000000 / 3.7;  //convert GeV to ke
0431     float tof = (*isim).timeOfFlight();
0432     float p = (*isim).pabs();
0433     float pt = (*isim).momentumAtEntry().perp();
0434     float theta = (*isim).thetaAtEntry();
0435     float phi = (*isim).phiAtEntry();
0436     int pid = ((*isim).particleType());
0437     int tid = (*isim).trackId();
0438     int procType = (*isim).processType();
0439
0440     float x = (*isim).entryPoint().x();  // width (row index, in col direction)
0441     float y = (*isim).entryPoint().y();  // length (col index, in row direction)
0442     float z = (*isim).entryPoint().z();  // thick
0443
0444     float x2 = (*isim).exitPoint().x();
0445     float y2 = (*isim).exitPoint().y();
0446     float z2 = (*isim).exitPoint().z();
0447
0448     float dz = abs(z - z2);
0449     bool moduleDirectionUp = (z < z2);  // for positive direction z2>z
0450
0451     float xpos = (x + x2) / 2.;
0452     float ypos = (y + y2) / 2.;
0453     float zpos = (z + z2) / 2.;
0454
0455     if (pt > 0.1) {
0456       goodHits++;
0457     }
0458     if (PRINT) {
0459       if (pt > 0.1)
0460         cout << " simhit: ";
0461       else if (pid == 11)
0462         cout << " delta: ";
0463       else
0464         cout << " low pt (sec?): ";
0465       cout << " id " << pid << " " << tid << " proc " << procType << " tof " << tof << " de " << eloss << " pt " << pt
0466            << " entry " << x << "/" << y << "/" << z << " lenz " << dz << " " << moduleDirectionUp << endl;
0467     }
0468     if (DEBUG)
0469       cout << "  center pos " << xpos << " " << ypos << " " << zpos;
0470
0471     LocalPoint loc(xpos, ypos, zpos);
0472     double gloX = theGeomDet->surface().toGlobal(loc).x();     //
0473     double gloY = theGeomDet->surface().toGlobal(loc).y();     //
0474     double gloR = theGeomDet->surface().toGlobal(loc).perp();  //
0475     double gloZ = theGeomDet->surface().toGlobal(loc).z();     //
0476     if (DEBUG)
0477       cout << ", global pos " << gloX << " " << gloY << " " << gloR << " " << gloZ << endl;
0478
0479     htest3->Fill(gloX, gloY);
0480     hdetunit->Fill(float(detId.rawId()));
0481     hpabs->Fill(p);
0482     htof->Fill(tof);
0483     hpid->Fill(float(abs(pid)));
0484     htid->Fill(float(tid));
0485     hpixid->Fill(float(detid));
0486     hpixsubid->Fill(float(subid));
0487     hlayerid->Fill(float(layer));
0488
0489     // Transform the theta from local module coordinates to global
0490     //if(theta<= PI/2.) theta = PI/2. - theta; // For +z global
0491     //else theta = (PI/2. + PI) - theta;
0492
0493     if (mode_ == "fpix") {
0494       if (disk == 1) {
0495         //cout<<" disk "<<disk<<endl;
0496         totalNumOfSimHits1++;
0497         heloss1->Fill(eloss);
0498         if (pid == 11)
0499           heloss1e->Fill(eloss);
0500         else
0501           heloss1mu->Fill(eloss);
0502         hladder1id->Fill(float(blade));
0503         hz1id->Fill(float(zindex));
0504         hthick1->Fill(dz);
0505         hlength1->Fill(y);
0506         hwidth1->Fill(x);
0507
0508         //SimHitMap1[detId.rawId()].push_back((*isim));
0509         htheta1->Fill(theta);
0510         hglobr1->Fill(gloR);
0511         hglobz1->Fill(gloZ);
0512         hdetphi1->Fill(detPhi);
0513
0514       } else if (disk == 2) {
0515         //cout<<" disk "<<disk<<endl;
0516         totalNumOfSimHits2++;
0517         heloss2->Fill(eloss);
0518         if (pid == 11)
0519           heloss2e->Fill(eloss);
0520         else
0521           heloss2mu->Fill(eloss);
0522         hladder2id->Fill(float(blade));
0523         hz2id->Fill(float(zindex));
0524         hthick2->Fill(dz);
0525         hlength2->Fill(y);
0526         hwidth2->Fill(x);
0527
0528         //SimHitMap2[detId.rawId()].push_back((*isim));
0529         hglobr2->Fill(gloR);
0530         hglobz2->Fill(gloZ);
0531         hdetphi2->Fill(detPhi);
0532
0533       }  // end disks
0534
0535     } else if (mode_ == "bpix") {
0536       if (layer == 1) {
0537         //cout<<" layer "<<layer<<endl;
0538         totalNumOfSimHits1++;
0539         heloss1->Fill(eloss);
0540         if (abs(pid) == 11)
0541           heloss1e->Fill(eloss);
0542         else
0543           heloss1mu->Fill(eloss);
0544         hladder1id->Fill(float(ladder));
0545         hz1id->Fill(float(module));
0546         hthick1->Fill(dz);
0547         hlength1->Fill(y);
0548         hwidth1->Fill(x);
0549         if (abs(pid) == 13 && p > 1.)
0550           hphi1->Fill(phi);
0551
0552         // Test half modules
0553         //        if(ladder==5 || ladder==6 || ladder==15 || ladder==16 ) { // half-modules
0554         //   hwidth1h->Fill(x);
0555         //   if(pid==13 && p>1.) {  // select primary muons with mom above 1.
0556         //     hphi1h->Fill(phi);
0557         //     hglobr1h->Fill(gloR);
0558         //   }
0559         //        } else {
0560         //        }
0561
0562         SimHitMap1[detId.rawId()].push_back((*isim));
0563         htheta1->Fill(theta);
0564         hglobr1->Fill(gloR);
0565         hglobz1->Fill(gloZ);
0566
0567         // Check the coordinate system and counting
0568         htest->Fill(gloZ, ypos);
0569         if (abs(pid) != 11)
0570           htest2->Fill(gloZ, eloss);
0571
0572         //if(pid!=11 && moduleDirectionUp)  hladder1idUp->Fill(float(ladder));
0573
0574         //hp1->Fill(float(ladder),detR,1);
0575         //hp2->Fill(float(ladder),detPhi);
0576         hdetphi1->Fill(detPhi);
0577
0578       } else if (layer == 2) {
0579         //cout<<" layer "<<layer<<endl;
0580         totalNumOfSimHits2++;
0581         heloss2->Fill(eloss);
0582         if (abs(pid) == 11)
0583           heloss2e->Fill(eloss);
0584         else
0585           heloss2mu->Fill(eloss);
0586         hladder2id->Fill(float(ladder));
0587         hz2id->Fill(float(module));
0588         hthick2->Fill(dz);
0589         hlength2->Fill(y);
0590         hwidth2->Fill(x);
0591         if (abs(pid) == 13 && p > 1.)
0592           hphi2->Fill(phi);
0593
0594         // check half modules
0595         //        if(ladder==8 || ladder==9 || ladder==24 || ladder==25 ) {
0596         //   hwidth2h->Fill(x);
0597         //        } else {
0598         //        }
0599
0600         SimHitMap2[detId.rawId()].push_back((*isim));
0601         hglobr2->Fill(gloR);
0602         hglobz2->Fill(gloZ);
0603         hdetphi2->Fill(detPhi);
0604
0605         // check up/down modules
0606         //if(pid!=11 && moduleDirectionUp) hladder2idUp->Fill(float(ladder));
0607
0608       } else if (layer == 3) {
0609         //cout<<" layer "<<layer<<endl;
0610         totalNumOfSimHits3++;
0611         heloss3->Fill(eloss);
0612         if (abs(pid) == 11)
0613           heloss3e->Fill(eloss);
0614         else
0615           heloss3mu->Fill(eloss);
0616
0617         hladder3id->Fill(float(ladder));
0618         hz3id->Fill(float(module));
0619         hthick3->Fill(dz);
0620         hlength3->Fill(y);
0621         hwidth3->Fill(x);
0622         if (abs(pid) == 13 && p > 1.)
0623           hphi3->Fill(phi);
0624
0625         // check half modules
0626         //       if(ladder==11 || ladder==12 || ladder==33 || ladder==34 ) {
0627         //   hwidth3h->Fill(x);
0628         //        } else {
0629         //        }
0630
0631         SimHitMap3[detId.rawId()].push_back((*isim));
0632         hglobr3->Fill(gloR);
0633         hglobz3->Fill(gloZ);
0634         hdetphi3->Fill(detPhi);
0635         // check up/down modules
0636         //if(pid!=11 && moduleDirectionUp) hladder3idUp->Fill(float(ladder));
0637
0638       }  // layers
0639     }    // end fpix/bpix
0640   }
0641
0642   hsimHitsPerLay1->Fill(float(totalNumOfSimHits1));
0643   hsimHitsPerLay2->Fill(float(totalNumOfSimHits2));
0644   hsimHitsPerLay3->Fill(float(totalNumOfSimHits3));
0645   hsimHits->Fill(float(totalNumOfSimHits));
0646   hsimHitsGood->Fill(float(goodHits));
0647
0648   int numberOfDetUnits1 = SimHitMap1.size();
0649   int numberOfDetUnits2 = SimHitMap2.size();
0650   int numberOfDetUnits3 = SimHitMap3.size();
0651   int numberOfDetUnits = numberOfDetUnits1 + numberOfDetUnits2 + numberOfDetUnits3;
0652
0653   if (PRINT)
0654     cout << " Number of full det-units = " << numberOfDetUnits << " total simhits = " << totalNumOfSimHits
0655          << " good simhits (pt>0.1) " << goodHits << endl;
0656
0657   hdetsPerLay1->Fill(float(numberOfDetUnits1));
0658   hdetsPerLay2->Fill(float(numberOfDetUnits2));
0659   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0660
0661   numEvents++;
0662   numSimHits += totalNumOfSimHits;
0663   numSimHitsGood += goodHits;
0664
0665   if (mode_ == "bpix") {
0666     map<unsigned int, vector<PSimHit>, less<unsigned int> >::iterator simhit_map_iterator;
0667     for (simhit_map_iterator = SimHitMap1.begin(); simhit_map_iterator != SimHitMap1.end(); simhit_map_iterator++) {
0668       //if(PRINT) cout << " Lay1 det = "<<simhit_map_iterator->first <<" simHits = "
0669       //          << (simhit_map_iterator->second).size()<< endl;
0670       hsimHitsPerDet1->Fill(float((simhit_map_iterator->second).size()));
0671     }
0672     for (simhit_map_iterator = SimHitMap2.begin(); simhit_map_iterator != SimHitMap2.end(); simhit_map_iterator++) {
0673       //if(PRINT) cout << " Lay2 det = "<<simhit_map_iterator->first <<" simHits = "
0674       //          << (simhit_map_iterator->second).size()<< endl;
0675       hsimHitsPerDet2->Fill(float((simhit_map_iterator->second).size()));
0676     }
0677     for (simhit_map_iterator = SimHitMap3.begin(); simhit_map_iterator != SimHitMap3.end(); simhit_map_iterator++) {
0678       //if(PRINT) cout << " Lay3 det = "<<simhit_map_iterator->first <<" simHits = "
0679       //          << (simhit_map_iterator->second).size() << endl;
0680       hsimHitsPerDet3->Fill(float((simhit_map_iterator->second).size()));
0681     }
0682   }  // of bpix
0683 }
0684 // ------------ method called to at the end of the job  ------------
0685 void PixelSimHitsTest::endJob() {
0686   cout << " End PixelSimHitsTest " << endl;
0687
0688   numEvents++;
0689   numSimHits = numSimHits / float(numEvents);
0690   numSimHitsGood = numSimHitsGood / float(numEvents);
0691
0692   cout << " Events " << numEvents << " simhits " << numSimHits << " simhits > 0.1gev " << numSimHitsGood << endl;
0693
0694 #ifdef CHECK_GEOM
0695   // To get module positions
0696   cout << " Module position" << endl;
0697   cout << " Layer Ladder Zindex    R      Z      Phi " << endl;
0698   for (int i = 0; i < 3; i++) {
0699     int max_lad = 0;
0700     if (i == 0)
0701       max_lad = 20;
0702     else if (i == 1)
0703       max_lad = 32;
0704     else if (i == 2)
0705       max_lad = 44;
0706     for (int n = 0; n < max_lad; n++) {
0707       for (int m = 0; m < 8; m++) {
0708         cout << "   " << i + 1 << "      " << n + 1 << "      " << m + 1 << "    " << modulePositionR[i][n][m] << " "
0709              << modulePositionZ[i][n][m] << " " << modulePositionPhi[i][n][m] << endl;
0710       }
0711     }
0712   }
0713 #endif
0714 }
0715
0716 //define this as a plug-in
0717 DEFINE_FWK_MODULE(PixelSimHitsTest);