HcalTestBeam/src/HcalTB06Histo.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     HcalTestBeam
0004 // Class  :     HcalTB06Histo
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author:
0010 //         Created:  Tue Oct 10 10:14:34 CEST 2006
0011 //
0012
0013 // system include files
0014 #include <iostream>
0015 #include <cmath>
0016
0017 // user include files
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06Histo.h"
0023
0024 //
0025 // constructors and destructor
0026 HcalTB06Histo::HcalTB06Histo(const edm::ParameterSet& ps) {
0027   verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
0028   double em1 = ps.getUntrackedParameter<double>("ETtotMax", 400.);
0029   double em2 = ps.getUntrackedParameter<double>("EHCalMax", 4.0);
0030   mkTree_ = ps.getUntrackedParameter<bool>("MakeTree", false);
0031   eBeam_ = 50.;
0032   mip_ = ps.getParameter<double>("MIP");
0033   edm::LogInfo("HcalTBSim") << "Verbose :" << verbose_ << " MakeTree: " << mkTree_ << " EMax: " << em1 << ":" << em2
0034                             << " MIP " << mip_;
0035
0036   // Book histograms
0037   edm::Service<TFileService> tfile;
0038
0039   if (!tfile.isAvailable())
0040     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0041                                       << "please add it to config file";
0042   iniE = tfile->make<TH1D>("iniE", "Incident Energy (GeV)", 4000, 0., em1);
0043   iEta = tfile->make<TH1D>("iEta", "Eta at incidence     ", 300, 0., 3.);
0044   iPhi = tfile->make<TH1D>("iPhi", "Phi at incidence     ", 300, -1., 1.);
0045   edepS = tfile->make<TH1D>("edepS", "Energy deposit == Total", 4000, 0., em1);
0046   edecS = tfile->make<TH1D>("edecS", "Energy deposit == ECal ", 300, -2., 28.);
0047   edhcS = tfile->make<TH1D>("edhcS", "Energy deposit == HCal ", 4000, 0., em2);
0048   edepN = tfile->make<TH1D>("edepN", "Etot/Ebeam   ", 200, -2.5, 2.5);
0049   edecN = tfile->make<TH1D>("edecN", "Eecal/Ebeam  ", 200, -2.5, 2.5);
0050   edhcN = tfile->make<TH1D>("edhcN", "Ehcal/Ebeam  ", 200, -2.5, 2.5);
0051   emhcN = tfile->make<TH1D>("emhcN", "Ehcal/Ebeam MIP in Ecal", 200, -2.5, 2.5);
0052   edehS = tfile->make<TH2D>("edehS", "Hcal vs Ecal", 100, 0., em1, 100, 0., em2);
0053
0054   if (mkTree_) {
0055     tree_ = tfile->make<TTree>("TB06Sim", "TB06Sim");
0056     tree_->Branch("eBeam_", &eBeam_, "eBeam_/D");
0057     tree_->Branch("etaBeam_", &etaBeam_, "etaBeam_/D");
0058     tree_->Branch("phiBeam_", &phiBeam_, "phiBeam_/D");
0059     tree_->Branch("edepEC_", &edepEC_, "edepEC_/D");
0060     tree_->Branch("edepHB_", &edepHB_, "edepHB_/D");
0061     tree_->Branch("edepHO_", &edepHO_, "edepHO_/D");
0062     tree_->Branch("noiseEC_", &noiseEC_, "noiseEC_/D");
0063     tree_->Branch("noiseHB_", &noiseHB_, "noiseHB_/D");
0064     tree_->Branch("noiseHO_", &noiseHO_, "noiseHO_/D");
0065     tree_->Branch("edepS1_", &edepS1_, "edepS1_/D");
0066     tree_->Branch("edepS2_", &edepS2_, "edepS2_/D");
0067     tree_->Branch("edepS3_", &edepS3_, "edepS3_/D");
0068     tree_->Branch("edepS4_", &edepS4_, "edepS4_/D");
0069     tree_->Branch("edepVC_", &edepVC_, "edepVC_/D");
0070     tree_->Branch("edepS7_", &edepS7_, "edepS7_/D");
0071     tree_->Branch("edepS8_", &edepS8_, "edepS8_/D");
0072   }
0073 }
0074
0075 HcalTB06Histo::~HcalTB06Histo() {}
0076
0077 //
0078 // member functions
0079 //
0080
0081 void HcalTB06Histo::fillPrimary(double energy, double eta, double phi) {
0082   if (verbose_)
0083     edm::LogInfo("HcalTBSim") << "HcalTB06Histo::fillPrimary: Energy " << energy << " Eta " << eta << " Phi " << phi;
0084   eBeam_ = energy;
0085   etaBeam_ = eta;
0086   phiBeam_ = phi;
0087   iniE->Fill(energy);
0088   iEta->Fill(eta);
0089   iPhi->Fill(phi);
0090 }
0091
0092 void HcalTB06Histo::fillEdep(double etots, double eecals, double ehcals) {
0093   if (verbose_)
0094     edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillEdep: Simulated Total " << etots << " ECal " << eecals << " HCal "
0095                               << ehcals;
0096   edepS->Fill(etots);
0097   edecS->Fill(eecals);
0098   edhcS->Fill(ehcals);
0099   edepN->Fill(etots / eBeam_);
0100   edecN->Fill(eecals / eBeam_);
0101   edhcN->Fill(ehcals / eBeam_);
0102   if (eecals <= mip_) {
0103     emhcN->Fill(etots / eBeam_);
0104   }
0105   edehS->Fill(eecals, ehcals);
0106 }
0107
0108 void HcalTB06Histo::fillTree(std::vector<double>& ecalo, std::vector<double>& etrig) {
0109   if (mkTree_) {
0110     edepEC_ = ecalo[0];
0111     noiseEC_ = ecalo[1];
0112     edepHB_ = ecalo[2];
0113     noiseHB_ = ecalo[3];
0114     edepHO_ = ecalo[4];
0115     noiseHO_ = ecalo[5];
0116     edepS1_ = etrig[0];
0117     edepS2_ = etrig[1];
0118     edepS3_ = etrig[2];
0119     edepS4_ = etrig[3];
0120     edepVC_ = etrig[4];
0121     edepS7_ = etrig[5];
0122     edepS8_ = etrig[6];
0123     tree_->Fill();
0124     if (verbose_)
0125       edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillTree: Energies " << edepEC_ << ":" << noiseEC_ << ":" << edepHB_
0126                                 << ":" << noiseHB_ << ":" << edepHO_ << ":" << noiseHO_ << " Trigger counters "
0127                                 << edepS1_ << ":" << edepS2_ << ":" << edepS3_ << ":" << edepS4_ << ":" << edepVC_
0128                                 << ":" << edepS7_ << ":" << edepS8_;
0129   }
0130 }