Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:58 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!