// -*- C++ -*-
//
// Package:    PixelMixedSimHitsTest
// Class:      PixelMixedSimHitsTest
//
/**\class PixelMixedSimHitsTest PixelMixedSimHitsTest.cc 

 Description: Test Barrel pixel simhits and fills root histograms. 
 SimHits are loaded from CrossingFrame.
  
 Implementation:
     <Notes on implementation>
*/
//
// Original Author:  V.Chiochia
//         Created:   CET 2006
//
//
// system include files
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DataFormats/Common/interface/Handle.h"

#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/ESHandle.h"

// my includes
//#include "Geometry/CommonDetUnit/interface/GeomDet.h"

#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
#include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"  //

#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"

#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
//#include "Geometry/Surface/interface/Surface.h"

// For ROOT
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TF1.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TProfile.h>

using namespace std;
using namespace edm;

//#define CHECK_GEOM

class PixelMixedSimHitsTest : public edm::one::EDAnalyzer {
public:
  explicit PixelMixedSimHitsTest(const edm::ParameterSet &);
  ~PixelMixedSimHitsTest();
  virtual void beginJob();
  virtual void analyze(const edm::Event &, const edm::EventSetup &);
  virtual void endJob();

private:
  // ----------member data ---------------------------
  edm::ParameterSet conf_;
  const static bool PRINT = false;
  typedef std::vector<std::string> vstring;
  vstring trackerContainers;
  edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken;
  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken;
  // Histograms
  TFile *hFile;
  TH1F *hBunchCrossing, *htof1m, *htof0, *htof1p, *htof2p;
  TH1F *heloss1, *heloss2, *heloss3, *hdetunit, *hpabs, *hpid, *htof, *htid;
  TH1F *hpixid, *hpixsubid, *hlayerid, *hladder1id, *hladder2id, *hladder3id, *hz1id, *hz2id, *hz3id;
  TH1F *hladder1idUp, *hladder2idUp, *hladder3idUp;
  TH1F *hthick1, *hthick2, *hthick3, *hlength1, *hlength2, *hlength3;
  TH1F *hwidth1, *hwidth2, *hwidth3;
  TH1F *hwidth1h, *hwidth2h, *hwidth3h;
  TH1F *hsimHitsPerDet1, *hsimHitsPerDet2, *hsimHitsPerDet3;
  TH1F *hsimHitsPerLay1, *hsimHitsPerLay2, *hsimHitsPerLay3;
  TH1F *hdetsPerLay1, *hdetsPerLay2, *hdetsPerLay3;
  TH1F *heloss1e, *heloss2e, *heloss3e;
  TH1F *heloss1mu, *heloss2mu, *heloss3mu;
  TH1F *htheta1, *htheta2, *htheta3;
  TH1F *hphi1, *hphi1h, *hphi2, *hphi3;
  TH1F *hdetr, *hdetz, *hdetphi1, *hdetphi2, *hdetphi3;
  TH1F *hglobr1, *hglobr2, *hglobr3, *hglobz1, *hglobz2, *hglobz3;
  TH1F *hglobr1h;
  TH1F *hcolsB, *hrowsB, *hcolsF, *hrowsF;
  TH1F *hglox1;
  //
  TH2F *htest, *htest2, *htest3, *htest4, *htest5;
  //
  TProfile *hp1, *hp2, *hp3, *hp4, *hp5;

#ifdef CHECK_GEOM
  float modulePositionZ[3][44][8];
  float modulePositionR[3][44][8];
  float modulePositionPhi[3][44][8];
#endif
};
//
PixelMixedSimHitsTest::PixelMixedSimHitsTest(const edm::ParameterSet &iConfig) : conf_(iConfig) {
  trackerTopoToken = esConsumes<TrackerTopology, TrackerTopologyRcd>();
  trackerGeomToken = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
  trackerContainers.clear();
  trackerContainers = iConfig.getParameter<std::vector<std::string> >("ROUList");
  cout << " Construct PixelSimHitsTest " << endl;
}

PixelMixedSimHitsTest::~PixelMixedSimHitsTest() {
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  cout << " Destroy PixelSimHitsTest " << endl;
}

// ------------ method called at the begining   ------------
void PixelMixedSimHitsTest::beginJob() {
  using namespace edm;
  cout << "Initialize PixelSimHitsTest " << endl;

  // put here whatever you want to do at the beginning of the job
  std::string outputFile = conf_.getParameter<std::string>("OutputFile");
  hFile = new TFile(outputFile.c_str(), "RECREATE");

  const float max_charge = 200.;  // in ke
  hBunchCrossing = new TH1F("hBunchCrossing", "Bunch Crossing", 10, -5., 5.);
  heloss1 = new TH1F("heloss1", "Eloss l1", 100, 0., max_charge);
  heloss2 = new TH1F("heloss2", "Eloss l2", 100, 0., max_charge);
  heloss3 = new TH1F("heloss3", "Eloss l3", 100, 0., max_charge);

  hdetunit = new TH1F("hdetunit", "Det unit", 1000, 302000000., 302300000.);
  hpabs = new TH1F("hpabs", "Pabs", 100, 0., 100.);
  htof = new TH1F("htof", "TOF", 50, -25., 25.);
  htof1m = new TH1F("htof1m", "TOF BC=-1", 250, -50., 75.);
  htof0 = new TH1F("htof0", "TOF BC=0", 250, -50., 75.);
  htof1p = new TH1F("htof1p", "TOF BC=+1", 250, -50., 75.);
  htof2p = new TH1F("htof2p", "TOF BC=+2", 250, -50., 75.);

  hpid = new TH1F("hpid", "PID", 1000, 0., 1000.);
  htid = new TH1F("htid", "Track id", 100, 0., 100.);

  hpixid = new TH1F("hpixid", "Pix det id", 10, 0., 10.);
  hpixsubid = new TH1F("hpixsubid", "Pix Barrel id", 10, 0., 10.);
  hlayerid = new TH1F("hlayerid", "Pix layer id", 10, 0., 10.);
  hladder1id = new TH1F("hladder1id", "Ladder L1 id", 100, -0.5, 49.5);
  hladder2id = new TH1F("hladder2id", "Ladder L2 id", 100, -0.5, 49.5);
  hladder3id = new TH1F("hladder3id", "Ladder L3 id", 100, -0.5, 49.5);
  hz1id = new TH1F("hz1id", "Z-index id L1", 10, 0., 10.);
  hz2id = new TH1F("hz2id", "Z-index id L2", 10, 0., 10.);
  hz3id = new TH1F("hz3id", "Z-index id L3", 10, 0., 10.);

  hthick1 = new TH1F("hthick1", "Det 1 Thinckess", 400, 0., 0.04);
  hthick2 = new TH1F("hthick2", "Det 2 Thinckess", 400, 0., 0.04);
  hthick3 = new TH1F("hthick3", "Det 3 Thinckess", 400, 0., 0.04);

  hlength1 = new TH1F("hlength1", "Det 1 Length", 700, -3.5, 3.5);
  hlength2 = new TH1F("hlength2", "Det 2 Length", 700, -3.5, 3.5);
  hlength3 = new TH1F("hlength3", "Det 3 Length", 700, -3.5, 3.5);

  hwidth1 = new TH1F("hwidth1", "Det 1 Width", 200, -1., 1.);
  hwidth2 = new TH1F("hwidth2", "Det 2 Width", 200, -1., 1.);
  hwidth3 = new TH1F("hwidth3", "Det 3 Width", 200, -1., 1.);

  hwidth1h = new TH1F("hwidth1h", "Det 1 Width half-m", 200, -1., 1.);
  hwidth2h = new TH1F("hwidth2h", "Det 2 Width half-m", 200, -1., 1.);
  hwidth3h = new TH1F("hwidth3h", "Det 3 Width half-m", 200, -1., 1.);

  hsimHitsPerDet1 = new TH1F("hsimHitsPerDet1", "SimHits per det l1", 200, -0.5, 199.5);
  hsimHitsPerDet2 = new TH1F("hsimHitsPerDet2", "SimHits per det l2", 200, -0.5, 199.5);
  hsimHitsPerDet3 = new TH1F("hsimHitsPerDet3", "SimHits per det l3", 200, -0.5, 199.5);
  hsimHitsPerLay1 = new TH1F("hsimHitsPerLay1", "SimHits per layer l1", 2000, -0.5, 1999.5);
  hsimHitsPerLay2 = new TH1F("hsimHitsPerLay2", "SimHits per layer l2", 2000, -0.5, 1999.5);
  hsimHitsPerLay3 = new TH1F("hsimHitsPerLay3", "SimHits per layer l3", 2000, -0.5, 1999.5);
  hdetsPerLay1 = new TH1F("hdetsPerLay1", "Full dets per layer l1", 161, -0.5, 160.5);
  hdetsPerLay3 = new TH1F("hdetsPerLay3", "Full dets per layer l3", 353, -0.5, 352.5);
  hdetsPerLay2 = new TH1F("hdetsPerLay2", "Full dets per layer l2", 257, -0.5, 256.5);
  heloss1e = new TH1F("heloss1e", "Eloss e l1", 100, 0., max_charge);
  heloss2e = new TH1F("heloss2e", "Eloss e l2", 100, 0., max_charge);
  heloss3e = new TH1F("heloss3e", "Eloss e l3", 100, 0., max_charge);

  heloss1mu = new TH1F("heloss1mu", "Eloss mu l1", 100, 0., max_charge);
  heloss2mu = new TH1F("heloss2mu", "Eloss mu l2", 100, 0., max_charge);
  heloss3mu = new TH1F("heloss3mu", "Eloss mu l3", 100, 0., max_charge);

  htheta1 = new TH1F("htheta1", "Theta det1", 350, 0.0, 3.5);
  htheta2 = new TH1F("htheta2", "Theta det2", 350, 0.0, 3.5);
  htheta3 = new TH1F("htheta3", "Theta det3", 350, 0.0, 3.5);
  hphi1 = new TH1F("hphi1", "phi l1", 1400, -3.5, 3.5);
  hphi2 = new TH1F("hphi2", "phi l2", 1400, -3.5, 3.5);
  hphi3 = new TH1F("hphi3", "phi l3", 1400, -3.5, 3.5);
  hphi1h = new TH1F("hphi1h", "phi l1", 1400, -3.5, 3.5);

  hdetr = new TH1F("hdetr", "det r", 1500, 0., 15.);
  hdetz = new TH1F("hdetz", "det z", 5200, -26., 26.);
  hdetphi1 = new TH1F("hdetphi1", "det phi l1", 700, -3.5, 3.5);
  hdetphi2 = new TH1F("hdetphi2", "det phi l2", 700, -3.5, 3.5);
  hdetphi3 = new TH1F("hdetphi3", "det phi l3", 700, -3.5, 3.5);

  hcolsB = new TH1F("hcolsB", "cols per bar det", 450, 0., 450.);
  hrowsB = new TH1F("hrowsB", "rows per bar det", 200, 0., 200.);
  hcolsF = new TH1F("hcolsF", "cols per for det", 300, 0., 300.);
  hrowsF = new TH1F("hrowsF", "rows per for det", 200, 0., 200.);

  hladder1idUp = new TH1F("hladder1idUp", "Ladder L1 id", 100, -0.5, 49.5);
  hladder2idUp = new TH1F("hladder2idUp", "Ladder L2 id", 100, -0.5, 49.5);
  hladder3idUp = new TH1F("hladder3idUp", "Ladder L3 id", 100, -0.5, 49.5);

  hglobr1 = new TH1F("hglobr1", "global r1", 150, 0., 15.);
  hglobz1 = new TH1F("hglobz1", "global z1", 540, -27., 27.);
  hglobr2 = new TH1F("hglobr2", "global r2", 150, 0., 15.);
  hglobz2 = new TH1F("hglobz2", "global z2", 540, -27., 27.);
  hglobr3 = new TH1F("hglobr3", "global r3", 150, 0., 15.);
  hglobz3 = new TH1F("hglobz3", "global z3", 540, -27., 27.);

  hglox1 = new TH1F("hglox1", "global x in l1", 200, -1., 1.);
  hglobr1h = new TH1F("hglobr1h", "global r1", 700, 4.1, 4.8);

  htest = new TH2F("htest", " ", 108, -27., 27., 35, -3.5, 3.5);
  htest2 = new TH2F("htest2", " ", 108, -27., 27., 60, 0., 600.);
  htest3 = new TH2F("htest3", " ", 240, -12., 12., 240, -12., 12.);
  htest4 = new TH2F("htest4", " ", 80, -4., 4., 100, -5., 5.);

  hp1 = new TProfile("hp1", " ", 50, 0., 50.);               // default option
  hp2 = new TProfile("hp2", " ", 50, 0., 50., " ");          // option set to " "
  hp3 = new TProfile("hp3", " ", 50, 0., 50., -100., 100.);  // with y limits
  hp4 = new TProfile("hp4", " ", 50, 0., 50.);
  hp5 = new TProfile("hp5", " ", 50, 0., 50.);

#ifdef CHECK_GEOM
  // To get the module position
  for (int i = 0; i < 3; i++) {
    for (int n = 0; n < 44; n++) {
      for (int m = 0; m < 8; m++) {
        modulePositionR[i][n][m] = -1;
        modulePositionZ[i][n][m] = -1;
        modulePositionPhi[i][n][m] = -1;
      }
    }
  }
#endif
}

// ------------ method called to produce the data  ------------
void PixelMixedSimHitsTest::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(trackerTopoToken);

  const double PI = 3.142;

  using namespace edm;
  if (PRINT)
    cout << " Analyze PixelSimHitsTest " << endl;

  // Get event setup (to get global transformation)
  edm::ESHandle<TrackerGeometry> geom = iSetup.getHandle(trackerGeomToken);
  const TrackerGeometry &theTracker(*geom);

  // Get input data
  int totalNumOfSimHits = 0;
  int totalNumOfSimHits1 = 0;
  int totalNumOfSimHits2 = 0;
  int totalNumOfSimHits3 = 0;

  // To count simhits per det module
  //typedef std::map<unsigned int, std::vector<PSimHit>,
  //std::less<unsigned int>> simhit_map;
  //typedef simhit_map::iterator simhit_map_iterator;
  //simhit_map SimHitMap;
  map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap1;
  map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap2;
  map<unsigned int, vector<PSimHit>, less<unsigned int> > SimHitMap3;

  //    Handle<PSimHitContainer> PixelBarrelHitsLowTof;
  //    Handle<PSimHitContainer> PixelBarrelHitsHighTof;

  //    iEvent.getByLabel("SimG4Object","TrackerHitsPixelBarrelLowTof",PixelBarrelHitsLowTof);
  //    iEvent.getByLabel("SimG4Object","TrackerHitsPixelBarrelHighTof",PixelBarrelHitsHighTof);

  // Step A: Get Inputs
  edm::Handle<CrossingFrame> cf;
  iEvent.getByType(cf);

  std::unique_ptr<MixCollection<PSimHit> > allPixelTrackerHits(
      new MixCollection<PSimHit>(cf.product(), trackerContainers));

  //   for(vector<PSimHit>::const_iterator isim = PixelBarrelHitsLowTof->begin();
  //       isim != PixelBarrelHitsLowTof->end(); ++isim){

  MixCollection<PSimHit>::iterator isim;
  for (isim = allPixelTrackerHits->begin(); isim != allPixelTrackerHits->end(); isim++) {
    totalNumOfSimHits++;
    // Det id
    DetId detId = DetId((*isim).detUnitId());
    unsigned int detid = detId.det();       // for pixel=1
    unsigned int subid = detId.subdetId();  // barrel=1

    if (detid != 1 && subid != 1)
      cout << " error in det id " << detid << " " << subid << endl;
    //if(PRINT) cout<<(*isim).detUnitId()<<" "<<detId.rawId()<<" "
    //	   <<detId.null()<<" "<<detid<<" "<<subid<<endl;

    //const GeomDetUnit * theGeomDet = theTracker.idToDet(detId);
    //const PixelGeomDetUnit * theGeomDet = theTracker.idToDet(detId);
    const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(detId));
    double detZ = theGeomDet->surface().position().z();
    double detR = theGeomDet->surface().position().perp();
    double detPhi = theGeomDet->surface().position().phi();
    hdetr->Fill(detR);
    hdetz->Fill(detZ);

    //const BoundPlane plane = theGeomDet->surface(); // does not work

    double detThick = theGeomDet->specificSurface().bounds().thickness();
    double detLength = theGeomDet->specificSurface().bounds().length();
    double detWidth = theGeomDet->specificSurface().bounds().width();

    int cols = theGeomDet->specificTopology().ncolumns();
    int rows = theGeomDet->specificTopology().nrows();

    hcolsB->Fill(float(cols));
    hrowsB->Fill(float(rows));
    if (PRINT)
      cout << "det z/r " << detZ << " " << detR << " " << detThick << " " << detLength << " " << detWidth << " " << cols
           << " " << rows << endl;

    unsigned int layer = tTopo->pxbLayer(detId.rawId);
    unsigned int ladder = tTopo->pxbLadder(detId.rawId);
    unsigned int zindex = tTopo->pxbModule(detId.rawId);
    if (PRINT)
      cout << "det id " << detId.rawId() << " " << detid << " " << subid << " " << layer << " " << ladder << " "
           << zindex << endl;
#ifdef CHECK_GEOM
    // To get the module position
    modulePositionR[layer - 1][ladder - 1][zindex - 1] = detR;
    modulePositionZ[layer - 1][ladder - 1][zindex - 1] = detZ;
    modulePositionPhi[layer - 1][ladder - 1][zindex - 1] = detPhi;
#endif

    // SimHit information
    float eloss = (*isim).energyLoss() * 1000000 / 3.7;  //convert GeV to ke
    float tof = (*isim).timeOfFlight();
    float p = (*isim).pabs();
    float theta = (*isim).thetaAtEntry();
    float phi = (*isim).phiAtEntry();
    int pid = abs((*isim).particleType());  // ignore sign
    int tid = (*isim).trackId();
    int procType = (*isim).processType();

    float x = (*isim).entryPoint().x();  // width (row index, in col direction)
    float y = (*isim).entryPoint().y();  // length (col index, in row direction)
    float z = (*isim).entryPoint().z();  // thick

    float x2 = (*isim).exitPoint().x();
    float y2 = (*isim).exitPoint().y();
    float z2 = (*isim).exitPoint().z();
    int bcrossing = (*isim).eventId().bunchCrossing();

    float dz = abs(z - z2);
    bool moduleDirectionUp = (z < z2);  // for positive direction z2>z

    float xpos = (x + x2) / 2.;
    float ypos = (y + y2) / 2.;
    float zpos = (z + z2) / 2.;

    if (PRINT)
      cout << " simhit " << pid << " " << tid << " " << procType << " " << tof << " " << eloss << " " << p << " " << x
           << " " << y << " " << z << " " << dz << endl;
    if (PRINT)
      cout << "  pos " << xpos << " " << ypos << " " << zpos;

    LocalPoint loc(xpos, ypos, zpos);
    //GlobalPoint glo = theGeomDet->surface().toGlobal(loc); // does not work!
    double gloX = theGeomDet->surface().toGlobal(loc).x();     //
    double gloY = theGeomDet->surface().toGlobal(loc).y();     //
    double gloR = theGeomDet->surface().toGlobal(loc).perp();  //
    double gloZ = theGeomDet->surface().toGlobal(loc).z();     //
    if (PRINT)
      cout << ", global " << gloX << " " << gloY << " " << gloR << " " << gloZ << endl;

    htest3->Fill(gloX, gloY);
    hdetunit->Fill(float(detId.rawId()));
    hpabs->Fill(p);
    htof->Fill(tof);
    hpid->Fill(float(pid));
    htid->Fill(float(tid));
    hpixid->Fill(float(detid));
    hpixsubid->Fill(float(subid));
    hlayerid->Fill(float(layer));
    hBunchCrossing->Fill(float(bcrossing));
    if (bcrossing == -1)
      htof1m->Fill(tof);
    else if (bcrossing == 0)
      htof0->Fill(tof);
    else if (bcrossing == 1)
      htof1p->Fill(tof);
    else if (bcrossing == 2)
      htof2p->Fill(tof);

    // Transform the theta from local module coordinates to global
    //if(theta<= PI/2.) theta = PI/2. - theta; // For +z global
    //else theta = (PI/2. + PI) - theta;

    if (layer == 1) {
      //cout<<" layer "<<layer<<endl;
      totalNumOfSimHits1++;
      heloss1->Fill(eloss);
      if (pid == 11)
        heloss1e->Fill(eloss);
      else
        heloss1mu->Fill(eloss);
      hladder1id->Fill(float(ladder));
      hz1id->Fill(float(zindex));
      hthick1->Fill(dz);
      hlength1->Fill(y);
      if (ladder == 5 || ladder == 6 || ladder == 15 || ladder == 16) {
        // half-modules
        hwidth1h->Fill(x);
        if (pid == 13 && p > 1.) {  // select primary muons
          hphi1h->Fill(phi);
          hglox1->Fill(gloX);
          hglobr1h->Fill(gloR);

          // 	   double gloX1 =
          // 	     theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).x(); //
          // 	   double gloY1 =
          // 	     theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).y(); //
          // 	   double gloR1 =
          // 	     theGeomDet->surface().toGlobal(LocalPoint(0,0,0)).perp();
          // 	   cout<<" "<<ladder<<" "<<gloX1<<" "<<gloY1<<" "<<gloR1<<" "
          // 	       <<detR<<" "<<detPhi<<" "<<detZ<<" "<<gloX<<" "<<gloY<<" "
          // 	       <<xpos<<" "<<ypos<<" "<<zpos<<endl;
        }
      } else {
        hwidth1->Fill(x);
        if (pid == 13 && p > 1.)
          hphi1->Fill(phi);
      }
      SimHitMap1[detId.rawId()].push_back((*isim));
      htheta1->Fill(theta);
      hglobr1->Fill(gloR);
      hglobz1->Fill(gloZ);

      // Check the coordinate system and counting
      htest->Fill(gloZ, ypos);
      if (pid != 11)
        htest2->Fill(gloZ, eloss);

      if (pid != 11 && moduleDirectionUp)
        hladder1idUp->Fill(float(ladder));

      if (ladder == 6)
        htest4->Fill(xpos, gloX);
      hp1->Fill(float(ladder), detR, 1);
      hp2->Fill(float(ladder), detPhi);
      hdetphi1->Fill(detPhi);

    } else if (layer == 2) {
      //cout<<" layer "<<layer<<endl;
      totalNumOfSimHits2++;
      heloss2->Fill(eloss);
      if (pid == 11)
        heloss2e->Fill(eloss);
      else
        heloss2mu->Fill(eloss);
      hladder2id->Fill(float(ladder));
      hz2id->Fill(float(zindex));
      hthick2->Fill(dz);
      hlength2->Fill(y);
      if (ladder == 8 || ladder == 9 || ladder == 24 || ladder == 25) {
        hwidth2h->Fill(x);
      } else {
        hwidth2->Fill(x);
        if (pid == 13 && p > 1.)
          hphi2->Fill(phi);
      }
      SimHitMap2[detId.rawId()].push_back((*isim));
      hglobr2->Fill(gloR);
      hglobz2->Fill(gloZ);
      hdetphi2->Fill(detPhi);
      if (pid != 11 && moduleDirectionUp)
        hladder2idUp->Fill(float(ladder));

    } else if (layer == 3) {
      //cout<<" layer "<<layer<<endl;
      totalNumOfSimHits3++;
      heloss3->Fill(eloss);
      if (pid == 11)
        heloss3e->Fill(eloss);
      else
        heloss3mu->Fill(eloss);

      hladder3id->Fill(float(ladder));
      hz3id->Fill(float(zindex));
      hthick3->Fill(dz);
      hlength3->Fill(y);
      if (ladder == 11 || ladder == 12 || ladder == 33 || ladder == 34) {
        hwidth3h->Fill(x);
      } else {
        hwidth3->Fill(x);
        if (pid == 13 && p > 1.)
          hphi3->Fill(phi);
      }
      SimHitMap3[detId.rawId()].push_back((*isim));
      hglobr3->Fill(gloR);
      hglobz3->Fill(gloZ);
      hdetphi3->Fill(detPhi);
      if (pid != 11 && moduleDirectionUp)
        hladder3idUp->Fill(float(ladder));
    }
  }

  hsimHitsPerLay1->Fill(float(totalNumOfSimHits1));
  hsimHitsPerLay2->Fill(float(totalNumOfSimHits2));
  hsimHitsPerLay3->Fill(float(totalNumOfSimHits3));

  int numberOfDetUnits1 = SimHitMap1.size();
  int numberOfDetUnits2 = SimHitMap2.size();
  int numberOfDetUnits3 = SimHitMap3.size();
  int numberOfDetUnits = numberOfDetUnits1 + numberOfDetUnits2 + numberOfDetUnits3;

  if (PRINT)
    cout << " Number of full det-units = " << numberOfDetUnits << " total simhits = " << totalNumOfSimHits << endl;

  hdetsPerLay1->Fill(float(numberOfDetUnits1));
  hdetsPerLay2->Fill(float(numberOfDetUnits2));
  hdetsPerLay3->Fill(float(numberOfDetUnits3));

  map<unsigned int, vector<PSimHit>, less<unsigned int> >::iterator simhit_map_iterator;
  for (simhit_map_iterator = SimHitMap1.begin(); simhit_map_iterator != SimHitMap1.end(); simhit_map_iterator++) {
    if (PRINT)
      cout << " Lay1 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
           << endl;
    hsimHitsPerDet1->Fill(float((simhit_map_iterator->second).size()));
  }
  for (simhit_map_iterator = SimHitMap2.begin(); simhit_map_iterator != SimHitMap2.end(); simhit_map_iterator++) {
    if (PRINT)
      cout << " Lay2 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
           << endl;
    hsimHitsPerDet2->Fill(float((simhit_map_iterator->second).size()));
  }
  for (simhit_map_iterator = SimHitMap3.begin(); simhit_map_iterator != SimHitMap3.end(); simhit_map_iterator++) {
    if (PRINT)
      cout << " Lay3 det = " << simhit_map_iterator->first << " simHits = " << (simhit_map_iterator->second).size()
           << endl;
    hsimHitsPerDet3->Fill(float((simhit_map_iterator->second).size()));
  }
}
// ------------ method called to at the end of the job  ------------
void PixelMixedSimHitsTest::endJob() {
  cout << " End PixelSimHitsTest " << endl;
  hFile->Write();
  hFile->Close();

#ifdef CHECK_GEOM
  // To get module positions
  cout << " Module position" << endl;
  cout << " Layer Ladder Zindex    R      Z      Phi " << endl;
  for (int i = 0; i < 3; i++) {
    int max_lad = 0;
    if (i == 0)
      max_lad = 20;
    else if (i == 1)
      max_lad = 32;
    else if (i == 2)
      max_lad = 44;
    for (int n = 0; n < max_lad; n++) {
      for (int m = 0; m < 8; m++) {
        cout << "   " << i + 1 << "      " << n + 1 << "      " << m + 1 << "    " << modulePositionR[i][n][m] << " "
             << modulePositionZ[i][n][m] << " " << modulePositionPhi[i][n][m] << endl;
      }
    }
  }
#endif
}

//define this as a plug-in
DEFINE_FWK_MODULE(PixelMixedSimHitsTest);
