HcalTB06Histo.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     HcalTestBeam
// Class  :     HcalTB06Histo
//
// Implementation:
//     <Notes on implementation>
//
// Original Author:
//         Created:  Tue Oct 10 10:14:34 CEST 2006
//

// system include files
#include <iostream>
#include <cmath>

// user include files
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "SimG4CMS/HcalTestBeam/interface/HcalTB06Histo.h"

//
// constructors and destructor
HcalTB06Histo::HcalTB06Histo(const edm::ParameterSet& ps) {
  verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
  double em1 = ps.getUntrackedParameter<double>("ETtotMax", 400.);
  double em2 = ps.getUntrackedParameter<double>("EHCalMax", 4.0);
  mkTree_ = ps.getUntrackedParameter<bool>("MakeTree", false);
  eBeam_ = 50.;
  mip_ = ps.getParameter<double>("MIP");
  edm::LogInfo("HcalTBSim") << "Verbose :" << verbose_ << " MakeTree: " << mkTree_ << " EMax: " << em1 << ":" << em2
                            << " MIP " << mip_;

  // Book histograms
  edm::Service<TFileService> tfile;

  if (!tfile.isAvailable())
    throw cms::Exception("BadConfig") << "TFileService unavailable: "
                                      << "please add it to config file";
  iniE = tfile->make<TH1D>("iniE", "Incident Energy (GeV)", 4000, 0., em1);
  iEta = tfile->make<TH1D>("iEta", "Eta at incidence     ", 300, 0., 3.);
  iPhi = tfile->make<TH1D>("iPhi", "Phi at incidence     ", 300, -1., 1.);
  edepS = tfile->make<TH1D>("edepS", "Energy deposit == Total", 4000, 0., em1);
  edecS = tfile->make<TH1D>("edecS", "Energy deposit == ECal ", 300, -2., 28.);
  edhcS = tfile->make<TH1D>("edhcS", "Energy deposit == HCal ", 4000, 0., em2);
  edepN = tfile->make<TH1D>("edepN", "Etot/Ebeam   ", 200, -2.5, 2.5);
  edecN = tfile->make<TH1D>("edecN", "Eecal/Ebeam  ", 200, -2.5, 2.5);
  edhcN = tfile->make<TH1D>("edhcN", "Ehcal/Ebeam  ", 200, -2.5, 2.5);
  emhcN = tfile->make<TH1D>("emhcN", "Ehcal/Ebeam MIP in Ecal", 200, -2.5, 2.5);
  edehS = tfile->make<TH2D>("edehS", "Hcal vs Ecal", 100, 0., em1, 100, 0., em2);

  if (mkTree_) {
    tree_ = tfile->make<TTree>("TB06Sim", "TB06Sim");
    tree_->Branch("eBeam_", &eBeam_, "eBeam_/D");
    tree_->Branch("etaBeam_", &etaBeam_, "etaBeam_/D");
    tree_->Branch("phiBeam_", &phiBeam_, "phiBeam_/D");
    tree_->Branch("edepEC_", &edepEC_, "edepEC_/D");
    tree_->Branch("edepHB_", &edepHB_, "edepHB_/D");
    tree_->Branch("edepHO_", &edepHO_, "edepHO_/D");
    tree_->Branch("noiseEC_", &noiseEC_, "noiseEC_/D");
    tree_->Branch("noiseHB_", &noiseHB_, "noiseHB_/D");
    tree_->Branch("noiseHO_", &noiseHO_, "noiseHO_/D");
    tree_->Branch("edepS1_", &edepS1_, "edepS1_/D");
    tree_->Branch("edepS2_", &edepS2_, "edepS2_/D");
    tree_->Branch("edepS3_", &edepS3_, "edepS3_/D");
    tree_->Branch("edepS4_", &edepS4_, "edepS4_/D");
    tree_->Branch("edepVC_", &edepVC_, "edepVC_/D");
    tree_->Branch("edepS7_", &edepS7_, "edepS7_/D");
    tree_->Branch("edepS8_", &edepS8_, "edepS8_/D");
  }
}

HcalTB06Histo::~HcalTB06Histo() {}

//
// member functions
//

void HcalTB06Histo::fillPrimary(double energy, double eta, double phi) {
  if (verbose_)
    edm::LogInfo("HcalTBSim") << "HcalTB06Histo::fillPrimary: Energy " << energy << " Eta " << eta << " Phi " << phi;
  eBeam_ = energy;
  etaBeam_ = eta;
  phiBeam_ = phi;
  iniE->Fill(energy);
  iEta->Fill(eta);
  iPhi->Fill(phi);
}

void HcalTB06Histo::fillEdep(double etots, double eecals, double ehcals) {
  if (verbose_)
    edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillEdep: Simulated Total " << etots << " ECal " << eecals << " HCal "
                              << ehcals;
  edepS->Fill(etots);
  edecS->Fill(eecals);
  edhcS->Fill(ehcals);
  edepN->Fill(etots / eBeam_);
  edecN->Fill(eecals / eBeam_);
  edhcN->Fill(ehcals / eBeam_);
  if (eecals <= mip_) {
    emhcN->Fill(etots / eBeam_);
  }
  edehS->Fill(eecals, ehcals);
}

void HcalTB06Histo::fillTree(std::vector<double>& ecalo, std::vector<double>& etrig) {
  if (mkTree_) {
    edepEC_ = ecalo[0];
    noiseEC_ = ecalo[1];
    edepHB_ = ecalo[2];
    noiseHB_ = ecalo[3];
    edepHO_ = ecalo[4];
    noiseHO_ = ecalo[5];
    edepS1_ = etrig[0];
    edepS2_ = etrig[1];
    edepS3_ = etrig[2];
    edepS4_ = etrig[3];
    edepVC_ = etrig[4];
    edepS7_ = etrig[5];
    edepS8_ = etrig[6];
    tree_->Fill();
    if (verbose_)
      edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillTree: Energies " << edepEC_ << ":" << noiseEC_ << ":" << edepHB_
                                << ":" << noiseHB_ << ":" << edepHO_ << ":" << noiseHO_ << " Trigger counters "
                                << edepS1_ << ":" << edepS2_ << ":" << edepS3_ << ":" << edepS4_ << ":" << edepVC_
                                << ":" << edepS7_ << ":" << edepS8_;
  }
}