Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-08 23:51:02 

0001 // system include files
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <vector>
0006 #include <map>
0007 #include <string>
0008 
0009 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0011 
0012 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023 #include "FWCore/Utilities/interface/transform.h"
0024 
0025 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0026 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0027 
0028 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0029 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0031 
0032 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0033 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0034 
0035 class HGCalSimHitValidation : public DQMEDAnalyzer {
0036 public:
0037   struct energysum {
0038     energysum() {
0039       etotal = 0;
0040       for (int i = 0; i < 6; ++i)
0041         eTime[i] = 0.;
0042     }
0043     double eTime[6], etotal;
0044   };
0045 
0046   struct hitsinfo {
0047     hitsinfo() {
0048       x = y = z = phi = eta = 0.0;
0049       cell = cell2 = sector = sector2 = type = layer = 0;
0050     }
0051     double x, y, z, phi, eta;
0052     int cell, cell2, sector, sector2, type, layer;
0053   };
0054 
0055   explicit HGCalSimHitValidation(const edm::ParameterSet&);
0056   ~HGCalSimHitValidation() override = default;
0057 
0058   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0059 
0060 protected:
0061   void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0062   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0063   void analyze(const edm::Event&, const edm::EventSetup&) override;
0064 
0065 private:
0066   void analyzeHits(std::vector<PCaloHit>& hits);
0067   void fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer);
0068   void fillHitsInfo(std::pair<hitsinfo, energysum> hit_, unsigned int itimeslice, double esum);
0069 
0070   TH1F* createHisto(std::string histname, const int nbins, float minIndexX, float maxIndexX, bool isLogX = true);
0071   void histoSetting(TH1F*& histo,
0072                     const char* xTitle,
0073                     const char* yTitle = "",
0074                     Color_t lineColor = kBlack,
0075                     Color_t markerColor = kBlack,
0076                     int linewidth = 1);
0077   void histoSetting(TH2F*& histo,
0078                     const char* xTitle,
0079                     const char* yTitle = "",
0080                     Color_t lineColor = kBlack,
0081                     Color_t markerColor = kBlack,
0082                     int linewidth = 1);
0083   void fillMuonTomoHistos(int partialType, std::pair<hitsinfo, energysum> hit_);
0084 
0085   // ----------member data ---------------------------
0086   const std::string nameDetector_, caloHitSource_;
0087   const HGCalDDDConstants* hgcons_;
0088   const std::vector<double> times_;
0089   const int verbosity_;
0090   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tok_hgcal_;
0091   const edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
0092   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0093   unsigned int layers_;
0094   int firstLayer_;
0095 
0096   std::vector<MonitorElement*> HitOccupancy_Plus_, HitOccupancy_Minus_;
0097   std::vector<MonitorElement*> EtaPhi_Plus_, EtaPhi_Minus_;
0098   MonitorElement *MeanHitOccupancy_Plus_, *MeanHitOccupancy_Minus_;
0099   static const unsigned int maxTime_ = 6;
0100   std::vector<MonitorElement*> energy_[maxTime_];
0101   std::vector<MonitorElement*> energyFWF_, energyFWCN_, energyFWCK_;
0102   std::vector<MonitorElement*> energyPWF_, energyPWCN_, energyPWCK_;
0103   std::vector<MonitorElement*> hitXYFWF_, hitXYFWCN_, hitXYFWCK_, hitXYB_;
0104   unsigned int nTimes_;
0105 };
0106 
0107 HGCalSimHitValidation::HGCalSimHitValidation(const edm::ParameterSet& iConfig)
0108     : nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
0109       caloHitSource_(iConfig.getParameter<std::string>("CaloHitSource")),
0110       times_(iConfig.getParameter<std::vector<double> >("TimeSlices")),
0111       verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0112       tok_hgcal_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0113           edm::ESInputTag{"", nameDetector_})),
0114       tok_hepMC_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0115       tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", caloHitSource_))),
0116       firstLayer_(1) {
0117   nTimes_ = (times_.size() > maxTime_) ? maxTime_ : times_.size();
0118 }
0119 
0120 void HGCalSimHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0121   edm::ParameterSetDescription desc;
0122   std::vector<double> times = {25.0, 1000.0};
0123   desc.add<std::string>("DetectorName", "HGCalEESensitive");
0124   desc.add<std::string>("CaloHitSource", "HGCHitsEE");
0125   desc.add<std::vector<double> >("TimeSlices", times);
0126   desc.addUntracked<int>("Verbosity", 0);
0127   desc.addUntracked<bool>("TestNumber", true);
0128   descriptions.add("hgcalSimHitValidationEE", desc);
0129 }
0130 
0131 void HGCalSimHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0132   //Generator input
0133   if (verbosity_ > 0) {
0134     const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_hepMC_);
0135     if (!evtMC.isValid()) {
0136       edm::LogVerbatim("HGCalValidation") << "no HepMCProduct found";
0137     } else {
0138       const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0139       unsigned int k(0);
0140       for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0141            ++p, ++k) {
0142         edm::LogVerbatim("HGCalValidation") << "Particle[" << k << "] with pt " << (*p)->momentum().perp() << " eta "
0143                                             << (*p)->momentum().eta() << " phi " << (*p)->momentum().phi();
0144       }
0145     }
0146   }
0147 
0148   //Now the hits
0149   const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hits_);
0150   if (theCaloHitContainers.isValid()) {
0151     if (verbosity_ > 0)
0152       edm::LogVerbatim("HGCalValidation") << " PcalohitItr = " << theCaloHitContainers->size();
0153     std::vector<PCaloHit> caloHits;
0154     caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0155     analyzeHits(caloHits);
0156   } else if (verbosity_ > 0) {
0157     edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not exist!";
0158   }
0159 }
0160 
0161 void HGCalSimHitValidation::analyzeHits(std::vector<PCaloHit>& hits) {
0162   std::map<int, int> OccupancyMap_plus, OccupancyMap_minus;
0163   OccupancyMap_plus.clear();
0164   OccupancyMap_minus.clear();
0165 
0166   std::map<uint32_t, std::pair<hitsinfo, energysum> > map_hits;
0167   map_hits.clear();
0168 
0169   if (verbosity_ > 0)
0170     edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << hits.size() << " PcaloHit elements";
0171   unsigned int nused(0);
0172   for (unsigned int i = 0; i < hits.size(); i++) {
0173     double energy = hits[i].energy();
0174     double time = hits[i].time();
0175     uint32_t id_ = hits[i].id();
0176     int cell, sector, subsector(0), layer, zside;
0177     int cell2(0), type(0);
0178     if (hgcons_->waferHexagon8()) {
0179       HGCSiliconDetId detId = HGCSiliconDetId(id_);
0180       cell = detId.cellU();
0181       cell2 = detId.cellV();
0182       sector = detId.waferU();
0183       subsector = detId.waferV();
0184       type = detId.type();
0185       layer = detId.layer();
0186       zside = detId.zside();
0187     } else if (hgcons_->tileTrapezoid()) {
0188       HGCScintillatorDetId detId = HGCScintillatorDetId(id_);
0189       sector = detId.ietaAbs();
0190       cell = detId.iphi();
0191       subsector = 1;
0192       type = detId.type();
0193       layer = detId.layer();
0194       zside = detId.zside();
0195     } else {
0196       edm::LogError("HGCalValidation") << "Wrong geometry mode " << hgcons_->geomMode();
0197       continue;
0198     }
0199     nused++;
0200     if (verbosity_ > 1)
0201       edm::LogVerbatim("HGCalValidation")
0202           << "Detector " << nameDetector_ << " zside = " << zside << " sector|wafer = " << sector << ":" << subsector
0203           << " type = " << type << " layer = " << layer << " cell = " << cell << ":" << cell2 << " energy = " << energy
0204           << " energyem = " << hits[i].energyEM() << " energyhad = " << hits[i].energyHad() << " time = " << time;
0205 
0206     HepGeom::Point3D<float> gcoord;
0207     std::pair<float, float> xy;
0208     if (hgcons_->waferHexagon8()) {
0209       xy = hgcons_->locateCell(zside, layer, sector, subsector, cell, cell2, false, true, false, false);
0210     } else {
0211       xy = hgcons_->locateCellTrap(zside, layer, sector, cell, false, false);
0212     }
0213     double zp = hgcons_->waferZ(layer, false);
0214     if (zside < 0)
0215       zp = -zp;
0216     float xp = (zp < 0) ? -xy.first : xy.first;
0217     gcoord = HepGeom::Point3D<float>(xp, xy.second, zp);
0218     double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0219     if (verbosity_ > 1)
0220       edm::LogVerbatim("HGCalValidation")
0221           << std::hex << id_ << std::dec << " global coordinate " << gcoord << " time " << time << ":" << tof;
0222     time -= tof;
0223 
0224     energysum esum;
0225     hitsinfo hinfo;
0226     if (map_hits.count(id_) != 0) {
0227       hinfo = map_hits[id_].first;
0228       esum = map_hits[id_].second;
0229     } else {
0230       hinfo.x = gcoord.x();
0231       hinfo.y = gcoord.y();
0232       hinfo.z = gcoord.z();
0233       hinfo.sector = sector;
0234       hinfo.sector2 = subsector;
0235       hinfo.cell = cell;
0236       hinfo.cell2 = cell2;
0237       hinfo.type = type;
0238       hinfo.layer = layer - firstLayer_;
0239       hinfo.phi = gcoord.getPhi();
0240       hinfo.eta = gcoord.getEta();
0241     }
0242     esum.etotal += energy;
0243     for (unsigned int k = 0; k < nTimes_; ++k) {
0244       if (time > 0 && time < times_[k])
0245         esum.eTime[k] += energy;
0246     }
0247 
0248     if (verbosity_ > 1)
0249       edm::LogVerbatim("HGCalValidation") << " -----------------------   gx = " << hinfo.x << " gy = " << hinfo.y
0250                                           << " gz = " << hinfo.z << " phi = " << hinfo.phi << " eta = " << hinfo.eta;
0251     map_hits[id_] = std::pair<hitsinfo, energysum>(hinfo, esum);
0252   }
0253   if (verbosity_ > 0)
0254     edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << map_hits.size()
0255                                         << " detector elements being hit";
0256 
0257   std::map<uint32_t, std::pair<hitsinfo, energysum> >::iterator itr;
0258   for (itr = map_hits.begin(); itr != map_hits.end(); ++itr) {
0259     hitsinfo hinfo = (*itr).second.first;
0260     energysum esum = (*itr).second.second;
0261     int layer = hinfo.layer;
0262     double eta = hinfo.eta;
0263     int type, part, orient;
0264     int partialType = -1;
0265     if ((nameDetector_ == "HGCalEESensitive") || (nameDetector_ == "HGCalHESiliconSensitive")) {
0266       HGCSiliconDetId detId = HGCSiliconDetId((*itr).first);
0267       std::tie(type, part, orient) = hgcons_->waferType(detId, false);
0268       partialType = part;
0269     }
0270 
0271     for (unsigned int itimeslice = 0; itimeslice < nTimes_; itimeslice++) {
0272       fillHitsInfo((*itr).second, itimeslice, esum.eTime[itimeslice]);
0273     }
0274 
0275     if (eta > 0.0)
0276       fillOccupancyMap(OccupancyMap_plus, layer);
0277     else
0278       fillOccupancyMap(OccupancyMap_minus, layer);
0279 
0280     fillMuonTomoHistos(partialType, (*itr).second);
0281   }
0282   if (verbosity_ > 0)
0283     edm::LogVerbatim("HGCalValidation") << "With map:used:total " << hits.size() << "|" << nused << "|"
0284                                         << map_hits.size() << " hits";
0285 
0286   for (auto const& itr : OccupancyMap_plus) {
0287     int layer = itr.first;
0288     int occupancy = itr.second;
0289     HitOccupancy_Plus_.at(layer)->Fill(occupancy);
0290   }
0291   for (auto const& itr : OccupancyMap_minus) {
0292     int layer = itr.first;
0293     int occupancy = itr.second;
0294     HitOccupancy_Minus_.at(layer)->Fill(occupancy);
0295   }
0296 }
0297 
0298 void HGCalSimHitValidation::fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer) {
0299   if (OccupancyMap.find(layer) != OccupancyMap.end()) {
0300     ++OccupancyMap[layer];
0301   } else {
0302     OccupancyMap[layer] = 1;
0303   }
0304 }
0305 
0306 void HGCalSimHitValidation::fillHitsInfo(std::pair<hitsinfo, energysum> hits, unsigned int itimeslice, double esum) {
0307   unsigned int ilayer = hits.first.layer;
0308   if (ilayer < layers_) {
0309     energy_[itimeslice].at(ilayer)->Fill(esum);
0310     if (itimeslice == 0) {
0311       EtaPhi_Plus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
0312       EtaPhi_Minus_.at(ilayer)->Fill(hits.first.eta, hits.first.phi);
0313     }
0314   } else {
0315     if (verbosity_ > 0)
0316       edm::LogVerbatim("HGCalValidation")
0317           << "Problematic Hit for " << nameDetector_ << " at sector " << hits.first.sector << ":" << hits.first.sector2
0318           << " layer " << hits.first.layer << " cell " << hits.first.cell << ":" << hits.first.cell2 << " energy "
0319           << hits.second.etotal;
0320   }
0321 }
0322 
0323 void HGCalSimHitValidation::fillMuonTomoHistos(int partialType, std::pair<hitsinfo, energysum> hits) {
0324   hitsinfo hinfo = hits.first;
0325   energysum esum = hits.second;
0326   double edep =
0327       esum.eTime[0] * CLHEP::GeV /
0328       CLHEP::keV;  //index 0 and 1 corresponds to 25 ns and 1000 ns, respectively. In addititon, chaging energy loss unit to keV.
0329 
0330   unsigned int ilayer = hinfo.layer;
0331   double x = hinfo.x * CLHEP::mm / CLHEP::cm;  // chaging length unit to cm.
0332   double y = hinfo.y * CLHEP::mm / CLHEP::cm;
0333   if (ilayer < layers_) {
0334     if (nameDetector_ == "HGCalEESensitive" or nameDetector_ == "HGCalHESiliconSensitive") {
0335       // Fill the energy loss histograms for MIP
0336       if (!TMath::AreEqualAbs(edep, 0.0, 1.e-5)) {  //to avoid peak at zero due Eloss less than 10 mili eV.
0337         if (hinfo.type == HGCSiliconDetId::HGCalFine) {
0338           if (partialType == 0)
0339             energyFWF_.at(ilayer)->Fill(edep);
0340           if (partialType > 0)
0341             energyPWF_.at(ilayer)->Fill(edep);
0342         }
0343         if (hinfo.type == HGCSiliconDetId::HGCalCoarseThin) {
0344           if (partialType == 0)
0345             energyFWCN_.at(ilayer)->Fill(edep);
0346           if (partialType > 0)
0347             energyPWCN_.at(ilayer)->Fill(edep);
0348         }
0349         if (hinfo.type == HGCSiliconDetId::HGCalCoarseThick) {
0350           if (partialType == 0)
0351             energyFWCK_.at(ilayer)->Fill(edep);
0352           if (partialType > 0)
0353             energyPWCK_.at(ilayer)->Fill(edep);
0354         }
0355       }
0356 
0357       // Fill the XY distribution of detector hits
0358       if (hinfo.type == HGCSiliconDetId::HGCalFine)
0359         hitXYFWF_.at(ilayer)->Fill(x, y);
0360 
0361       if (hinfo.type == HGCSiliconDetId::HGCalCoarseThin)
0362         hitXYFWCN_.at(ilayer)->Fill(x, y);
0363 
0364       if (hinfo.type == HGCSiliconDetId::HGCalCoarseThick)
0365         hitXYFWCK_.at(ilayer)->Fill(x, y);
0366 
0367     }  //is Silicon
0368     if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0369       hitXYB_.at(ilayer)->Fill(x, y);
0370     }  //is Scintillator
0371   }    //layer condition
0372 }
0373 
0374 // ------------ method called when starting to processes a run  ------------
0375 void HGCalSimHitValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0376   hgcons_ = &iSetup.getData(tok_hgcal_);
0377   layers_ = hgcons_->layers(false);
0378   firstLayer_ = hgcons_->firstLayer();
0379   if (verbosity_ > 0)
0380     edm::LogVerbatim("HGCalValidation") << nameDetector_ << " defined with " << layers_ << " Layers with first at "
0381                                         << firstLayer_;
0382 }
0383 
0384 void HGCalSimHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0385   iB.setCurrentFolder("HGCAL/HGCalSimHitsV/" + nameDetector_);
0386 
0387   std::ostringstream histoname;
0388   for (unsigned int il = 0; il < layers_; ++il) {
0389     int ilayer = firstLayer_ + static_cast<int>(il);
0390     auto istr1 = std::to_string(ilayer);
0391     while (istr1.size() < 2) {
0392       istr1.insert(0, "0");
0393     }
0394     histoname.str("");
0395     histoname << "HitOccupancy_Plus_layer_" << istr1;
0396     HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Plus", 501, -0.5, 500.5));
0397     histoname.str("");
0398     histoname << "HitOccupancy_Minus_layer_" << istr1;
0399     HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "HitOccupancy_Minus", 501, -0.5, 500.5));
0400 
0401     histoname.str("");
0402     histoname << "EtaPhi_Plus_"
0403               << "layer_" << istr1;
0404     EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
0405     histoname.str("");
0406     histoname << "EtaPhi_Minus_"
0407               << "layer_" << istr1;
0408     EtaPhi_Minus_.push_back(
0409         iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
0410 
0411     for (unsigned int itimeslice = 0; itimeslice < nTimes_; itimeslice++) {
0412       histoname.str("");
0413       histoname << "energy_time_" << itimeslice << "_layer_" << istr1;
0414       energy_[itimeslice].push_back(iB.book1D(histoname.str().c_str(), "energy_", 100, 0, 0.1));
0415     }
0416 
0417     ///////////// Histograms for Energy loss in full wafers////////////
0418     if ((nameDetector_ == "HGCalEESensitive") || (nameDetector_ == "HGCalHESiliconSensitive")) {
0419       histoname.str("");
0420       histoname << "energy_FullWafer_Fine_layer_" << istr1;
0421       TH1F* hEdepFWF = createHisto(histoname.str(), 100, 0., 400., false);
0422       histoSetting(hEdepFWF, "Eloss (keV)", "", kRed, kRed, 2);
0423       energyFWF_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWF));
0424       hEdepFWF->Delete();
0425 
0426       histoname.str("");
0427       histoname << "energy_FullWafer_CoarseThin_layer_" << istr1;
0428       TH1F* hEdepFWCN = createHisto(histoname.str(), 100, 0., 400., false);
0429       histoSetting(hEdepFWCN, "Eloss (keV)", "", kGreen + 1, kGreen + 1, 2);
0430       energyFWCN_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWCN));
0431       hEdepFWCN->Delete();
0432 
0433       histoname.str("");
0434       histoname << "energy_FullWafer_CoarseThick_layer_" << istr1;
0435       TH1F* hEdepFWCK = createHisto(histoname.str(), 100, 0., 400., false);
0436       histoSetting(hEdepFWCK, "Eloss (keV)", "", kMagenta, kMagenta, 2);
0437       energyFWCK_.push_back(iB.book1D(histoname.str().c_str(), hEdepFWCK));
0438       hEdepFWCK->Delete();
0439     }
0440 
0441     ///////////////////////////////////////////////////////////////////
0442 
0443     ///////////// Histograms for Energy loss in partial wafers////////////
0444     if ((nameDetector_ == "HGCalEESensitive") || (nameDetector_ == "HGCalHESiliconSensitive")) {
0445       histoname.str("");
0446       histoname << "energy_PartialWafer_Fine_layer_" << istr1;
0447       TH1F* hEdepPWF = createHisto(histoname.str(), 100, 0., 400., false);
0448       histoSetting(hEdepPWF, "Eloss (keV)", "", kRed, kRed, 2);
0449       energyPWF_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWF));
0450       hEdepPWF->Delete();
0451 
0452       histoname.str("");
0453       histoname << "energy_PartialWafer_CoarseThin_layer_" << istr1;
0454       TH1F* hEdepPWCN = createHisto(histoname.str(), 100, 0., 400., false);
0455       histoSetting(hEdepPWCN, "Eloss (keV)", "", kGreen + 1, kGreen + 1, 2);
0456       energyPWCN_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWCN));
0457       hEdepPWCN->Delete();
0458 
0459       histoname.str("");
0460       histoname << "energy_PartialWafer_CoarseThick_layer_" << istr1;
0461       TH1F* hEdepPWCK = createHisto(histoname.str(), 100, 0., 400., false);
0462       histoSetting(hEdepPWCK, "Eloss (keV)", "", kMagenta, kMagenta, 2);
0463       energyPWCK_.push_back(iB.book1D(histoname.str().c_str(), hEdepPWCK));
0464       hEdepPWCK->Delete();
0465     }
0466     ///////////////////////////////////////////////////////////////////
0467 
0468     // ///////////// Histograms for the XY distribution of fired cells/scintillator tiles ///////////////
0469     if ((nameDetector_ == "HGCalEESensitive") || (nameDetector_ == "HGCalHESiliconSensitive")) {
0470       histoname.str("");
0471       histoname << "hitXY_FullWafer_Fine_layer_" << istr1;
0472       TH2F* hitXYFWF = new TH2F(
0473           Form("hitXYFWF_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0474       histoSetting(hitXYFWF, "x (cm)", "y (cm)", kRed, kRed);
0475       hitXYFWF_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWF));
0476       hitXYFWF->Delete();
0477 
0478       histoname.str("");
0479       histoname << "hitXY_FullWafer_CoarseThin_layer_" << istr1;
0480       TH2F* hitXYFWCN = new TH2F(
0481           Form("hitXYFWCN_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0482       histoSetting(hitXYFWCN, "x (cm)", "y (cm)", kGreen + 1, kGreen + 1);
0483       hitXYFWCN_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWCN));
0484       hitXYFWCN->Delete();
0485 
0486       histoname.str("");
0487       histoname << "hitXY_FullWafer_CoarseThick_layer_" << istr1;
0488       TH2F* hitXYFWCK = new TH2F(
0489           Form("hitXYFWCK_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0490       histoSetting(hitXYFWCK, "x (cm)", "y (cm)", kMagenta, kMagenta);
0491       hitXYFWCK_.push_back(iB.book2D(histoname.str().c_str(), hitXYFWCK));
0492       hitXYFWCK->Delete();
0493     }
0494 
0495     if (nameDetector_ == "HGCalHEScintillatorSensitive") {
0496       histoname.str("");
0497       histoname << "hitXY_Scintillator_layer_" << istr1;
0498       TH2F* hitXYB = new TH2F(
0499           Form("hitXYB_%s", histoname.str().c_str()), histoname.str().c_str(), 100, -300., 300., 100, -300., 300.);
0500       histoSetting(hitXYB, "x (cm)", "y (cm)", kBlue, kBlue);
0501       hitXYB_.push_back(iB.book2D(histoname.str().c_str(), hitXYB));
0502       hitXYB->Delete();
0503     }
0504     //////////////////////////////////////////////////////////////////////////////////////////////////
0505   }
0506   MeanHitOccupancy_Plus_ = iB.book1D("MeanHitOccupancy_Plus", "MeanHitOccupancy_Plus", layers_, 0.5, layers_ + 0.5);
0507   MeanHitOccupancy_Minus_ = iB.book1D("MeanHitOccupancy_Minus", "MeanHitOccupancy_Minus", layers_, 0.5, layers_ + 0.5);
0508 }
0509 
0510 TH1F* HGCalSimHitValidation::createHisto(
0511     std::string histname, const int nbins, float minIndexX, float maxIndexX, bool isLogX) {
0512   TH1F* hist = nullptr;
0513   if (isLogX) {
0514     Double_t xbins[nbins + 1];
0515     double dx = (maxIndexX - minIndexX) / nbins;
0516     for (int i = 0; i <= nbins; i++) {
0517       xbins[i] = TMath::Power(10, (minIndexX + i * dx));
0518     }
0519     hist = new TH1F(Form("hEdep_%s", histname.c_str()), histname.c_str(), nbins, xbins);
0520   } else {
0521     hist = new TH1F(Form("hEdep_%s", histname.c_str()), histname.c_str(), nbins, minIndexX, maxIndexX);
0522   }
0523   return hist;
0524 }
0525 
0526 void HGCalSimHitValidation::histoSetting(
0527     TH1F*& histo, const char* xTitle, const char* yTitle, Color_t lineColor, Color_t markerColor, int lineWidth) {
0528   histo->SetStats();
0529   histo->SetLineColor(lineColor);
0530   histo->SetLineWidth(lineWidth);
0531   histo->SetMarkerColor(markerColor);
0532   histo->GetXaxis()->SetTitle(xTitle);
0533   histo->GetYaxis()->SetTitle(yTitle);
0534 }
0535 
0536 void HGCalSimHitValidation::histoSetting(
0537     TH2F*& histo, const char* xTitle, const char* yTitle, Color_t lineColor, Color_t markerColor, int lineWidth) {
0538   histo->SetStats();
0539   histo->SetLineColor(lineColor);
0540   histo->SetLineWidth(lineWidth);
0541   histo->SetMarkerColor(markerColor);
0542   histo->SetMarkerStyle(kFullCircle);
0543   histo->SetMarkerSize(0.6);
0544   histo->GetXaxis()->SetTitle(xTitle);
0545   histo->GetYaxis()->SetTitle(yTitle);
0546 }
0547 #include "FWCore/Framework/interface/MakerMacros.h"
0548 //define this as a plug-in
0549 DEFINE_FWK_MODULE(HGCalSimHitValidation);
[ 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!