SiPixelDigitizer/test/PixelMixedSimHitsTest.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PixelMixedSimHitsTest
0004 // Class:      PixelMixedSimHitsTest
0005 //
0006 /**\class PixelMixedSimHitsTest PixelMixedSimHitsTest.cc
0007
0008  Description: Test Barrel pixel simhits and fills root histograms.
0009  SimHits are loaded from CrossingFrame.
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  V.Chiochia
0016 //         Created:   CET 2006
0017 //
0018 //
0019 // system include files
0020 #include <iostream>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "DataFormats/Common/interface/Handle.h"
0029
0030 #include "FWCore/Utilities/interface/InputTag.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/ServiceRegistry/interface/Service.h"
0033 #include "FWCore/Framework/interface/ESHandle.h"
0034
0035 // my includes
0036 //#include "Geometry/CommonDetUnit/interface/GeomDet.h"
0037
0038 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0039 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0040 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0041 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"  //
0042
0043 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0044 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0045 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0046
0047 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0048 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0049 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0050 //#include "Geometry/Surface/interface/Surface.h"
0051
0052 // For ROOT
0053 #include <TROOT.h>
0054 #include <TChain.h>
0055 #include <TFile.h>
0056 #include <TF1.h>
0057 #include <TH2F.h>
0058 #include <TH1F.h>
0059 #include <TProfile.h>
0060
0061 using namespace std;
0062 using namespace edm;
0063
0064 //#define CHECK_GEOM
0065
0066 class PixelMixedSimHitsTest : public edm::one::EDAnalyzer {
0067 public:
0068   explicit PixelMixedSimHitsTest(const edm::ParameterSet &);
0069   ~PixelMixedSimHitsTest();
0070   virtual void beginJob();
0071   virtual void analyze(const edm::Event &, const edm::EventSetup &);
0072   virtual void endJob();
0073
0074 private:
0075   // ----------member data ---------------------------
0076   edm::ParameterSet conf_;
0077   const static bool PRINT = false;
0078   typedef std::vector<std::string> vstring;
0079   vstring trackerContainers;
0080   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken;
0081   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken;
0082   // Histograms
0083   TFile *hFile;
0084   TH1F *hBunchCrossing, *htof1m, *htof0, *htof1p, *htof2p;
0085   TH1F *heloss1, *heloss2, *heloss3, *hdetunit, *hpabs, *hpid, *htof, *htid;
0086   TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id;
0087   TH1F *hladder1idUp, *hladder2idUp, *hladder3idUp;
0088   TH1F *hthick1, *hthick2, *hthick3, *hlength1, *hlength2, *hlength3;
0089   TH1F *hwidth1, *hwidth2, *hwidth3;
0090   TH1F *hwidth1h, *hwidth2h, *hwidth3h;
0091   TH1F *hsimHitsPerDet1, *hsimHitsPerDet2, *hsimHitsPerDet3;
0092   TH1F *hsimHitsPerLay1, *hsimHitsPerLay2, *hsimHitsPerLay3;
0093   TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3;
0094   TH1F *heloss1e, *heloss2e, *heloss3e;
0095   TH1F *heloss1mu, *heloss2mu, *heloss3mu;
0096   TH1F *htheta1, *htheta2, *htheta3;
0097   TH1F *hphi1, *hphi1h, *hphi2, *hphi3;
0098   TH1F *hdetr, *hdetz, *hdetphi1, *hdetphi2, *hdetphi3;
0099   TH1F *hglobr1, *hglobr2, *hglobr3, *hglobz1, *hglobz2, *hglobz3;
0100   TH1F *hglobr1h;
0101   TH1F *hcolsB, *hrowsB, *hcolsF, *hrowsF;
0102   TH1F *hglox1;
0103   //
0104   TH2F *htest, *htest2, *htest3, *htest4, *htest5;
0105   //
0106   TProfile *hp1, *hp2, *hp3, *hp4, *hp5;
0107
0108 #ifdef CHECK_GEOM
0109   float modulePositionZ[3][44][8];
0110   float modulePositionR[3][44][8];
0111   float modulePositionPhi[3][44][8];
0112 #endif
0113 };
0114 //
0115 PixelMixedSimHitsTest::PixelMixedSimHitsTest(const edm::ParameterSet &iConfig) : conf_(iConfig) {
0116   trackerTopoToken = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0117   trackerGeomToken = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0118   trackerContainers.clear();
0119   trackerContainers = iConfig.getParameter<std::vector<std::string> >("ROUList");
0120   cout << " Construct PixelSimHitsTest " << endl;
0121 }
0122
0123 PixelMixedSimHitsTest::~PixelMixedSimHitsTest() {
0124   // do anything here that needs to be done at desctruction time
0125   // (e.g. close files, deallocate resources etc.)
0126   cout << " Destroy PixelSimHitsTest " << endl;
0127 }
0128
0129 // ------------ method called at the begining   ------------
0130 void PixelMixedSimHitsTest::beginJob() {
0131   using namespace edm;
0132   cout << "Initialize PixelSimHitsTest " << endl;
0133
0134   // put here whatever you want to do at the beginning of the job
0135   std::string outputFile = conf_.getParameter<std::string>("OutputFile");
0136   hFile = new TFile(outputFile.c_str(), "RECREATE");
0137
0138   const float max_charge = 200.;  // in ke
0139   hBunchCrossing = new TH1F("hBunchCrossing", "Bunch Crossing", 10, -5., 5.);
0140   heloss1 = new TH1F("heloss1", "Eloss l1", 100, 0., max_charge);
0141   heloss2 = new TH1F("heloss2", "Eloss l2", 100, 0., max_charge);
0142   heloss3 = new TH1F("heloss3", "Eloss l3", 100, 0., max_charge);
0143
0144   hdetunit = new TH1F("hdetunit", "Det unit", 1000, 302000000., 302300000.);
0145   hpabs = new TH1F("hpabs", "Pabs", 100, 0., 100.);
0146   htof = new TH1F("htof", "TOF", 50, -25., 25.);
0147   htof1m = new TH1F("htof1m", "TOF BC=-1", 250, -50., 75.);
0148   htof0 = new TH1F("htof0", "TOF BC=0", 250, -50., 75.);
0149   htof1p = new TH1F("htof1p", "TOF BC=+1", 250, -50., 75.);
0150   htof2p = new TH1F("htof2p", "TOF BC=+2", 250, -50., 75.);
0151
0152   hpid = new TH1F("hpid", "PID", 1000, 0., 1000.);
0153   htid = new TH1F("htid", "Track id", 100, 0., 100.);
0154
0155   hpixid = new TH1F("hpixid", "Pix det id", 10, 0., 10.);
0156   hpixsubid = new TH1F("hpixsubid", "Pix Barrel id", 10, 0., 10.);
0157   hlayerid = new TH1F("hlayerid", "Pix layer id", 10, 0., 10.);
0158   hladder1id = new TH1F("hladder1id", "Ladder L1 id", 100, -0.5, 49.5);
0159   hladder2id = new TH1F("hladder2id", "Ladder L2 id", 100, -0.5, 49.5);
0160   hladder3id = new TH1F("hladder3id", "Ladder L3 id", 100, -0.5, 49.5);
0161   hz1id = new TH1F("hz1id", "Z-index id L1", 10, 0., 10.);
0162   hz2id = new TH1F("hz2id", "Z-index id L2", 10, 0., 10.);
0163   hz3id = new TH1F("hz3id", "Z-index id L3", 10, 0., 10.);
0164
0165   hthick1 = new TH1F("hthick1", "Det 1 Thinckess", 400, 0., 0.04);
0166   hthick2 = new TH1F("hthick2", "Det 2 Thinckess", 400, 0., 0.04);
0167   hthick3 = new TH1F("hthick3", "Det 3 Thinckess", 400, 0., 0.04);
0168
0169   hlength1 = new TH1F("hlength1", "Det 1 Length", 700, -3.5, 3.5);
0170   hlength2 = new TH1F("hlength2", "Det 2 Length", 700, -3.5, 3.5);
0171   hlength3 = new TH1F("hlength3", "Det 3 Length", 700, -3.5, 3.5);
0172
0173   hwidth1 = new TH1F("hwidth1", "Det 1 Width", 200, -1., 1.);
0174   hwidth2 = new TH1F("hwidth2", "Det 2 Width", 200, -1., 1.);
0175   hwidth3 = new TH1F("hwidth3", "Det 3 Width", 200, -1., 1.);
0176
0177   hwidth1h = new TH1F("hwidth1h", "Det 1 Width half-m", 200, -1., 1.);
0178   hwidth2h = new TH1F("hwidth2h", "Det 2 Width half-m", 200, -1., 1.);
0179   hwidth3h = new TH1F("hwidth3h", "Det 3 Width half-m", 200, -1., 1.);
0180
0181   hsimHitsPerDet1 = new TH1F("hsimHitsPerDet1", "SimHits per det l1", 200, -0.5, 199.5);
0182   hsimHitsPerDet2 = new TH1F("hsimHitsPerDet2", "SimHits per det l2", 200, -0.5, 199.5);
0183   hsimHitsPerDet3 = new TH1F("hsimHitsPerDet3", "SimHits per det l3", 200, -0.5, 199.5);
0184   hsimHitsPerLay1 = new TH1F("hsimHitsPerLay1", "SimHits per layer l1", 2000, -0.5, 1999.5);
0185   hsimHitsPerLay2 = new TH1F("hsimHitsPerLay2", "SimHits per layer l2", 2000, -0.5, 1999.5);
0186   hsimHitsPerLay3 = new TH1F("hsimHitsPerLay3", "SimHits per layer l3", 2000, -0.5, 1999.5);
0187   hdetsPerLay1 = new TH1F("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5);
0188   hdetsPerLay3 = new TH1F("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5);
0189   hdetsPerLay2 = new TH1F("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5);
0190   heloss1e = new TH1F("heloss1e", "Eloss e l1", 100, 0., max_charge);
0191   heloss2e = new TH1F("heloss2e", "Eloss e l2", 100, 0., max_charge);
0192   heloss3e = new TH1F("heloss3e", "Eloss e l3", 100, 0., max_charge);
0193
0194   heloss1mu = new TH1F("heloss1mu", "Eloss mu l1", 100, 0., max_charge);
0195   heloss2mu = new TH1F("heloss2mu", "Eloss mu l2", 100, 0., max_charge);
0196   heloss3mu = new TH1F("heloss3mu", "Eloss mu l3", 100, 0., max_charge);
0197
0198   htheta1 = new TH1F("htheta1", "Theta det1", 350, 0.0, 3.5);
0199   htheta2 = new TH1F("htheta2", "Theta det2", 350, 0.0, 3.5);
0200   htheta3 = new TH1F("htheta3", "Theta det3", 350, 0.0, 3.5);
0201   hphi1 = new TH1F("hphi1", "phi l1", 1400, -3.5, 3.5);
0202   hphi2 = new TH1F("hphi2", "phi l2", 1400, -3.5, 3.5);
0203   hphi3 = new TH1F("hphi3", "phi l3", 1400, -3.5, 3.5);
0204   hphi1h = new TH1F("hphi1h", "phi l1", 1400, -3.5, 3.5);
0205
0206   hdetr = new TH1F("hdetr", "det r", 1500, 0., 15.);
0207   hdetz = new TH1F("hdetz", "det z", 5200, -26., 26.);
0208   hdetphi1 = new TH1F("hdetphi1", "det phi l1", 700, -3.5, 3.5);
0209   hdetphi2 = new TH1F("hdetphi2", "det phi l2", 700, -3.5, 3.5);
0210   hdetphi3 = new TH1F("hdetphi3", "det phi l3", 700, -3.5, 3.5);
0211
0212   hcolsB = new TH1F("hcolsB", "cols per bar det", 450, 0., 450.);
0213   hrowsB = new TH1F("hrowsB", "rows per bar det", 200, 0., 200.);
0214   hcolsF = new TH1F("hcolsF", "cols per for det", 300, 0., 300.);
0215   hrowsF = new TH1F("hrowsF", "rows per for det", 200, 0., 200.);
0216
0217   hladder1idUp = new TH1F("hladder1idUp", "Ladder L1 id", 100, -0.5, 49.5);
0218   hladder2idUp = new TH1F("hladder2idUp", "Ladder L2 id", 100, -0.5, 49.5);
0219   hladder3idUp = new TH1F("hladder3idUp", "Ladder L3 id", 100, -0.5, 49.5);
0220
0221   hglobr1 = new TH1F("hglobr1", "global r1", 150, 0., 15.);
0222   hglobz1 = new TH1F("hglobz1", "global z1", 540, -27., 27.);
0223   hglobr2 = new TH1F("hglobr2", "global r2", 150, 0., 15.);
0224   hglobz2 = new TH1F("hglobz2", "global z2", 540, -27., 27.);
0225   hglobr3 = new TH1F("hglobr3", "global r3", 150, 0., 15.);
0226   hglobz3 = new TH1F("hglobz3", "global z3", 540, -27., 27.);
0227
0228   hglox1 = new TH1F("hglox1", "global x in l1", 200, -1., 1.);
0229   hglobr1h = new TH1F("hglobr1h", "global r1", 700, 4.1, 4.8);
0230
0231   htest = new TH2F("htest", " ", 108, -27., 27., 35, -3.5, 3.5);
0232   htest2 = new TH2F("htest2", " ", 108, -27., 27., 60, 0., 600.);
0233   htest3 = new TH2F("htest3", " ", 240, -12., 12., 240, -12., 12.);
0234   htest4 = new TH2F("htest4", " ", 80, -4., 4., 100, -5., 5.);
0235
0236   hp1 = new TProfile("hp1", " ", 50, 0., 50.);               // default option
0237   hp2 = new TProfile("hp2", " ", 50, 0., 50., " ");          // option set to " "
0238   hp3 = new TProfile("hp3", " ", 50, 0., 50., -100., 100.);  // with y limits
0239   hp4 = new TProfile("hp4", " ", 50, 0., 50.);
0240   hp5 = new TProfile("hp5", " ", 50, 0., 50.);
0241
0242 #ifdef CHECK_GEOM
0243   // To get the module position
0244   for (int i = 0; i < 3; i++) {
0245     for (int n = 0; n < 44; n++) {
0246       for (int m = 0; m < 8; m++) {
0247         modulePositionR[i][n][m] = -1;
0248         modulePositionZ[i][n][m] = -1;
0249         modulePositionPhi[i][n][m] = -1;
0250       }
0251     }
0252   }
0253 #endif
0254 }
0255
0256 // ------------ method called to produce the data  ------------
0257 void PixelMixedSimHitsTest::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0258   //Retrieve tracker topology from geometry
0259   edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(trackerTopoToken);
0260
0261   const double PI = 3.142;
0262
0263   using namespace edm;
0264   if (PRINT)
0265     cout << " Analyze PixelSimHitsTest " << endl;
0266
0267   // Get event setup (to get global transformation)
0268   edm::ESHandle<TrackerGeometry> geom = iSetup.getHandle(trackerGeomToken);
0269   const TrackerGeometry &theTracker(*geom);
0270
0271   // Get input data
0272   int totalNumOfSimHits = 0;
0273   int totalNumOfSimHits1 = 0;
0274   int totalNumOfSimHits2 = 0;
0275   int totalNumOfSimHits3 = 0;
0276
0277   // To count simhits per det module
0278   //typedef std::map<unsigned int, std::vector<PSimHit>,
0279   //std::less<unsigned int>> simhit_map;
0280   //typedef simhit_map::iterator simhit_map_iterator;
0281   //simhit_map SimHitMap;
0282   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap1;
0283   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap2;
0284   map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap3;
0285
0286   //    Handle<PSimHitContainer> PixelBarrelHitsLowTof;
0287   //    Handle<PSimHitContainer> PixelBarrelHitsHighTof;
0288
0289   //    iEvent.getByLabel("SimG4Object","TrackerHitsPixelBarrelLowTof",PixelBarrelHitsLowTof);
0290   //    iEvent.getByLabel("SimG4Object","TrackerHitsPixelBarrelHighTof",PixelBarrelHitsHighTof);
0291
0292   // Step A: Get Inputs
0293   edm::Handle<CrossingFrame> cf;
0294   iEvent.getByType(cf);
0295
0296   std::unique_ptr<MixCollection<PSimHit> > allPixelTrackerHits(
0297       new MixCollection<PSimHit>(cf.product(), trackerContainers));
0298
0299   //   for(vector<PSimHit>::const_iterator isim = PixelBarrelHitsLowTof->begin();
0300   //       isim != PixelBarrelHitsLowTof->end(); ++isim){
0301
0302   MixCollection<PSimHit>::iterator isim;
0303   for (isim = allPixelTrackerHits->begin(); isim != allPixelTrackerHits->end(); isim++) {
0304     totalNumOfSimHits++;
0305     // Det id
0306     DetId detId = DetId((*isim).detUnitId());
0307     unsigned int detid = detId.det();       // for pixel=1
0308     unsigned int subid = detId.subdetId();  // barrel=1
0309
0310     if (detid != 1 && subid != 1)
0311       cout << " error in det id " << detid << " " << subid << endl;
0312     //if(PRINT) cout<<(*isim).detUnitId()<<" "<<detId.rawId()<<" "
0313     //     <<detId.null()<<" "<<detid<<" "<<subid<<endl;
0314
0315     //const GeomDetUnit * theGeomDet = theTracker.idToDet(detId);
0316     //const PixelGeomDetUnit * theGeomDet = theTracker.idToDet(detId);
0317     const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
0318     double detZ = theGeomDet->surface().position().z();
0319     double detR = theGeomDet->surface().position().perp();
0320     double detPhi = theGeomDet->surface().position().phi();
0321     hdetr->Fill(detR);
0322     hdetz->Fill(detZ);
0323
0324     //const BoundPlane plane = theGeomDet->surface(); // does not work
0325
0326     double detThick = theGeomDet->specificSurface().bounds().thickness();
0327     double detLength = theGeomDet->specificSurface().bounds().length();
0328     double detWidth = theGeomDet->specificSurface().bounds().width();
0329
0330     int cols = theGeomDet->specificTopology().ncolumns();
0331     int rows = theGeomDet->specificTopology().nrows();
0332
0333     hcolsB->Fill(float(cols));
0334     hrowsB->Fill(float(rows));
0335     if (PRINT)
0336       cout << "det z/r " << detZ << " " << detR << " " << detThick << " " << detLength << " " << detWidth << " " << cols
0337            << " " << rows << endl;
0338
0339     unsigned int layer = tTopo->pxbLayer(detId.rawId);
0340     unsigned int ladder = tTopo->pxbLadder(detId.rawId);
0341     unsigned int zindex = tTopo->pxbModule(detId.rawId);
0342     if (PRINT)
0343       cout << "det id " << detId.rawId() << " " << detid << " " << subid << " " << layer << " " << ladder << " "
0344            << zindex << endl;
0345 #ifdef CHECK_GEOM
0346     // To get the module position
0347     modulePositionR[layer - 1][ladder - 1][zindex - 1] = detR;
0348     modulePositionZ[layer - 1][ladder - 1][zindex - 1] = detZ;
0349     modulePositionPhi[layer - 1][ladder - 1][zindex - 1] = detPhi;
0350 #endif
0351
0352     // SimHit information
0353     float eloss = (*isim).energyLoss() * 1000000 / 3.7;  //convert GeV to ke
0354     float tof = (*isim).timeOfFlight();
0355     float p = (*isim).pabs();
0356     float theta = (*isim).thetaAtEntry();
0357     float phi = (*isim).phiAtEntry();
0358     int pid = abs((*isim).particleType());  // ignore sign
0359     int tid = (*isim).trackId();
0360     int procType = (*isim).processType();
0361
0362     float x = (*isim).entryPoint().x();  // width (row index, in col direction)
0363     float y = (*isim).entryPoint().y();  // length (col index, in row direction)
0364     float z = (*isim).entryPoint().z();  // thick
0365
0366     float x2 = (*isim).exitPoint().x();
0367     float y2 = (*isim).exitPoint().y();
0368     float z2 = (*isim).exitPoint().z();
0369     int bcrossing = (*isim).eventId().bunchCrossing();
0370
0371     float dz = abs(z - z2);
0372     bool moduleDirectionUp = (z < z2);  // for positive direction z2>z
0373
0374     float xpos = (x + x2) / 2.;
0375     float ypos = (y + y2) / 2.;
0376     float zpos = (z + z2) / 2.;
0377
0378     if (PRINT)
0379       cout << " simhit " << pid << " " << tid << " " << procType << " " << tof << " " << eloss << " " << p << " " << x
0380            << " " << y << " " << z << " " << dz << endl;
0381     if (PRINT)
0382       cout << "  pos " << xpos << " " << ypos << " " << zpos;
0383
0384     LocalPoint loc(xpos, ypos, zpos);
0385     //GlobalPoint glo = theGeomDet->surface().toGlobal(loc); // does not work!
0386     double gloX = theGeomDet->surface().toGlobal(loc).x();     //
0387     double gloY = theGeomDet->surface().toGlobal(loc).y();     //
0388     double gloR = theGeomDet->surface().toGlobal(loc).perp();  //
0389     double gloZ = theGeomDet->surface().toGlobal(loc).z();     //
0390     if (PRINT)
0391       cout << ", global " << gloX << " " << gloY << " " << gloR << " " << gloZ << endl;
0392
0393     htest3->Fill(gloX, gloY);
0394     hdetunit->Fill(float(detId.rawId()));
0395     hpabs->Fill(p);
0396     htof->Fill(tof);
0397     hpid->Fill(float(pid));
0398     htid->Fill(float(tid));
0399     hpixid->Fill(float(detid));
0400     hpixsubid->Fill(float(subid));
0401     hlayerid->Fill(float(layer));
0402     hBunchCrossing->Fill(float(bcrossing));
0403     if (bcrossing == -1)
0404       htof1m->Fill(tof);
0405     else if (bcrossing == 0)
0406       htof0->Fill(tof);
0407     else if (bcrossing == 1)
0408       htof1p->Fill(tof);
0409     else if (bcrossing == 2)
0410       htof2p->Fill(tof);
0411
0412     // Transform the theta from local module coordinates to global
0413     //if(theta<= PI/2.) theta = PI/2. - theta; // For +z global
0414     //else theta = (PI/2. + PI) - theta;
0415
0416     if (layer == 1) {
0417       //cout<<" layer "<<layer<<endl;
0418       totalNumOfSimHits1++;
0419       heloss1->Fill(eloss);
0420       if (pid == 11)
0421         heloss1e->Fill(eloss);
0422       else
0423         heloss1mu->Fill(eloss);
0424       hladder1id->Fill(float(ladder));
0425       hz1id->Fill(float(zindex));
0426       hthick1->Fill(dz);
0427       hlength1->Fill(y);
0428       if (ladder == 5 || ladder == 6 || ladder == 15 || ladder == 16) {
0429         // half-modules
0430         hwidth1h->Fill(x);
0431         if (pid == 13 && p > 1.) {  // select primary muons
0432           hphi1h->Fill(phi);
0433           hglox1->Fill(gloX);
0434           hglobr1h->Fill(gloR);
0435
0436           //       double gloX1 =
0437           //         theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).x(); //
0438           //       double gloY1 =
0439           //         theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).y(); //
0440           //       double gloR1 =
0441           //         theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).perp();
0442           //       cout<<" "<<ladder<<" "<<gloX1<<" "<<gloY1<<" "<<gloR1<<" "
0443           //           <<detR<<" "<<detPhi<<" "<<detZ<<" "<<gloX<<" "<<gloY<<" "
0444           //           <<xpos<<" "<<ypos<<" "<<zpos<<endl;
0445         }
0446       } else {
0447         hwidth1->Fill(x);
0448         if (pid == 13 && p > 1.)
0449           hphi1->Fill(phi);
0450       }
0451       SimHitMap1[detId.rawId()].push_back((*isim));
0452       htheta1->Fill(theta);
0453       hglobr1->Fill(gloR);
0454       hglobz1->Fill(gloZ);
0455
0456       // Check the coordinate system and counting
0457       htest->Fill(gloZ, ypos);
0458       if (pid != 11)
0459         htest2->Fill(gloZ, eloss);
0460
0461       if (pid != 11 && moduleDirectionUp)
0462         hladder1idUp->Fill(float(ladder));
0463
0464       if (ladder == 6)
0465         htest4->Fill(xpos, gloX);
0466       hp1->Fill(float(ladder), detR, 1);
0467       hp2->Fill(float(ladder), detPhi);
0468       hdetphi1->Fill(detPhi);
0469
0470     } else if (layer == 2) {
0471       //cout<<" layer "<<layer<<endl;
0472       totalNumOfSimHits2++;
0473       heloss2->Fill(eloss);
0474       if (pid == 11)
0475         heloss2e->Fill(eloss);
0476       else
0477         heloss2mu->Fill(eloss);
0478       hladder2id->Fill(float(ladder));
0479       hz2id->Fill(float(zindex));
0480       hthick2->Fill(dz);
0481       hlength2->Fill(y);
0482       if (ladder == 8 || ladder == 9 || ladder == 24 || ladder == 25) {
0483         hwidth2h->Fill(x);
0484       } else {
0485         hwidth2->Fill(x);
0486         if (pid == 13 && p > 1.)
0487           hphi2->Fill(phi);
0488       }
0489       SimHitMap2[detId.rawId()].push_back((*isim));
0490       hglobr2->Fill(gloR);
0491       hglobz2->Fill(gloZ);
0492       hdetphi2->Fill(detPhi);
0493       if (pid != 11 && moduleDirectionUp)
0494         hladder2idUp->Fill(float(ladder));
0495
0496     } else if (layer == 3) {
0497       //cout<<" layer "<<layer<<endl;
0498       totalNumOfSimHits3++;
0499       heloss3->Fill(eloss);
0500       if (pid == 11)
0501         heloss3e->Fill(eloss);
0502       else
0503         heloss3mu->Fill(eloss);
0504
0505       hladder3id->Fill(float(ladder));
0506       hz3id->Fill(float(zindex));
0507       hthick3->Fill(dz);
0508       hlength3->Fill(y);
0509       if (ladder == 11 || ladder == 12 || ladder == 33 || ladder == 34) {
0510         hwidth3h->Fill(x);
0511       } else {
0512         hwidth3->Fill(x);
0513         if (pid == 13 && p > 1.)
0514           hphi3->Fill(phi);
0515       }
0516       SimHitMap3[detId.rawId()].push_back((*isim));
0517       hglobr3->Fill(gloR);
0518       hglobz3->Fill(gloZ);
0519       hdetphi3->Fill(detPhi);
0520       if (pid != 11 && moduleDirectionUp)
0521         hladder3idUp->Fill(float(ladder));
0522     }
0523   }
0524
0525   hsimHitsPerLay1->Fill(float(totalNumOfSimHits1));
0526   hsimHitsPerLay2->Fill(float(totalNumOfSimHits2));
0527   hsimHitsPerLay3->Fill(float(totalNumOfSimHits3));
0528
0529   int numberOfDetUnits1 = SimHitMap1.size();
0530   int numberOfDetUnits2 = SimHitMap2.size();
0531   int numberOfDetUnits3 = SimHitMap3.size();
0532   int numberOfDetUnits = numberOfDetUnits1 + numberOfDetUnits2 + numberOfDetUnits3;
0533
0534   if (PRINT)
0535     cout << " Number of full det-units = " << numberOfDetUnits << " total simhits = " << totalNumOfSimHits << endl;
0536
0537   hdetsPerLay1->Fill(float(numberOfDetUnits1));
0538   hdetsPerLay2->Fill(float(numberOfDetUnits2));
0539   hdetsPerLay3->Fill(float(numberOfDetUnits3));
0540
0541   map<unsigned int, vector<PSimHit>, less<unsigned int> >::iterator simhit_map_iterator;
0542   for (simhit_map_iterator = SimHitMap1.begin(); simhit_map_iterator != SimHitMap1.end(); simhit_map_iterator++) {
0543     if (PRINT)
0544       cout << " Lay1 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
0545            << endl;
0546     hsimHitsPerDet1->Fill(float((simhit_map_iterator->second).size()));
0547   }
0548   for (simhit_map_iterator = SimHitMap2.begin(); simhit_map_iterator != SimHitMap2.end(); simhit_map_iterator++) {
0549     if (PRINT)
0550       cout << " Lay2 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
0551            << endl;
0552     hsimHitsPerDet2->Fill(float((simhit_map_iterator->second).size()));
0553   }
0554   for (simhit_map_iterator = SimHitMap3.begin(); simhit_map_iterator != SimHitMap3.end(); simhit_map_iterator++) {
0555     if (PRINT)
0556       cout << " Lay3 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
0557            << endl;
0558     hsimHitsPerDet3->Fill(float((simhit_map_iterator->second).size()));
0559   }
0560 }
0561 // ------------ method called to at the end of the job  ------------
0562 void PixelMixedSimHitsTest::endJob() {
0563   cout << " End PixelSimHitsTest " << endl;
0564   hFile->Write();
0565   hFile->Close();
0566
0567 #ifdef CHECK_GEOM
0568   // To get module positions
0569   cout << " Module position" << endl;
0570   cout << " Layer Ladder Zindex    R      Z      Phi " << endl;
0571   for (int i = 0; i < 3; i++) {
0572     int max_lad = 0;
0573     if (i == 0)
0574       max_lad = 20;
0575     else if (i == 1)
0576       max_lad = 32;
0577     else if (i == 2)
0578       max_lad = 44;
0579     for (int n = 0; n < max_lad; n++) {
0580       for (int m = 0; m < 8; m++) {
0581         cout << "   " << i + 1 << "      " << n + 1 << "      " << m + 1 << "    " << modulePositionR[i][n][m] << " "
0582              << modulePositionZ[i][n][m] << " " << modulePositionPhi[i][n][m] << endl;
0583       }
0584     }
0585   }
0586 #endif
0587 }
0588
0589 //define this as a plug-in
0590 DEFINE_FWK_MODULE(PixelMixedSimHitsTest);