EcalSimHitsValidation.cc

Go to the documentation of this file.

/*
 * \file EcalSimHitsValidation.cc
 *
 * \author C.Rovelli
 *
*/

#include <DataFormats/EcalDetId/interface/EBDetId.h>
#include <DataFormats/EcalDetId/interface/EEDetId.h>
#include <DataFormats/EcalDetId/interface/ESDetId.h>
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "Validation/EcalHits/interface/EcalSimHitsValidation.h"
#include "DQMServices/Core/interface/DQMStore.h"

using namespace cms;
using namespace edm;
using namespace std;

EcalSimHitsValidation::EcalSimHitsValidation(const edm::ParameterSet& ps):
  g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
  HepMCToken(consumes<edm::HepMCProduct>(ps.getParameter<std::string>("moduleLabelMC")))
{
  EBHitsCollectionToken = consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("EBHitsCollection")));
  EEHitsCollectionToken = consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("EEHitsCollection")));
  ESHitsCollectionToken = consumes<edm::PCaloHitContainer>(edm::InputTag(g4InfoLabel, ps.getParameter<std::string>("ESHitsCollection")));

  // DQM ROOT output
  outputFile_ = ps.getUntrackedParameter<std::string>("outputFile", "");
 
  if ( !outputFile_.empty() ) {
    edm::LogInfo("OutputInfo") << " Ecal SimHits Task histograms will be saved to " << outputFile_.c_str();
  } else {
    edm::LogInfo("OutputInfo") << " Ecal SimHits Task histograms will NOT be saved";
  }
 
  // verbosity switch
  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
 
  // DQMServices                                                        
  dbe_ = nullptr;

  // get hold of back-end interface
  dbe_ = edm::Service<DQMStore>().operator->();           
  if ( dbe_ ) {
    if ( verbose_ ) { dbe_->setVerbose(1); } 
    else            { dbe_->setVerbose(0); }
  }
                                                                                                            
  if ( dbe_ ) {
    if ( verbose_ ) dbe_->showDirStructure();
  }
 
  meGunEnergy_ = nullptr;
  meGunEta_    = nullptr;   
  meGunPhi_    = nullptr;   
  meEBEnergyFraction_  = nullptr;
  meEEEnergyFraction_  = nullptr;
  meESEnergyFraction_  = nullptr;

  Char_t histo[200];
 
  
  if ( dbe_ ) {
    dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
  
    sprintf (histo, "EcalSimHitsValidation Gun Momentum" ) ;
    meGunEnergy_ = dbe_->book1D(histo, histo, 100, 0., 1000.);
  
    sprintf (histo, "EcalSimHitsValidation Gun Eta" ) ;
    meGunEta_ = dbe_->book1D(histo, histo, 700, -3.5, 3.5);
  
    sprintf (histo, "EcalSimHitsValidation Gun Phi" ) ;
    meGunPhi_ = dbe_->book1D(histo, histo, 360, 0., 360.);

    sprintf (histo, "EcalSimHitsValidation Barrel fraction of energy" ) ;
    meEBEnergyFraction_ = dbe_->book1D(histo, histo, 100 , 0. , 1.1);
  
    sprintf (histo, "EcalSimHitsValidation Endcap fraction of energy" ) ;
    meEEEnergyFraction_ = dbe_->book1D(histo, histo, 100 , 0. , 1.1);
  
    sprintf (histo, "EcalSimHitsValidation Preshower fraction of energy" ) ;
    meESEnergyFraction_ = dbe_->book1D(histo, histo, 60 , 0. , 0.001);
  }
 
}

EcalSimHitsValidation::~EcalSimHitsValidation(){
 
  if ( !outputFile_.empty() && dbe_ ) dbe_->save(outputFile_);

}

void EcalSimHitsValidation::beginJob(){

}

void EcalSimHitsValidation::endJob(){

}

void EcalSimHitsValidation::analyze(const edm::Event& e, const edm::EventSetup& c){

  edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
  
  std::vector<PCaloHit>  theEBCaloHits;
  std::vector<PCaloHit>  theEECaloHits;
  std::vector<PCaloHit>  theESCaloHits;

  edm::Handle<edm::HepMCProduct> MCEvt;
  edm::Handle<edm::PCaloHitContainer> EcalHitsEB;
  edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
  edm::Handle<edm::PCaloHitContainer> EcalHitsES;

  e.getByToken(HepMCToken, MCEvt);
  e.getByToken(EBHitsCollectionToken, EcalHitsEB);
  e.getByToken(EEHitsCollectionToken, EcalHitsEE);
  e.getByToken(ESHitsCollectionToken, EcalHitsES);

  for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
        p != MCEvt->GetEvent()->particles_end(); ++p ) {

    double htheta = (*p)->momentum().theta();
    double heta = -99999.;
    if( tan(htheta * 0.5) > 0 ) {
      heta = -log(tan(htheta * 0.5));
    }
    double hphi = (*p)->momentum().phi();
    hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
    hphi = hphi / M_PI * 180.;

    LogDebug("EventInfo") << "Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n" << "Energy = "<< (*p)->momentum().e() << " Eta = " << heta << " Phi = " << hphi;

    if (meGunEnergy_) meGunEnergy_->Fill((*p)->momentum().e());
    if (meGunEta_)    meGunEta_   ->Fill(heta);
    if (meGunPhi_)    meGunPhi_   ->Fill(hphi);

  }
  
  double EBEnergy_ = 0.;
  if ( EcalHitsEB.isValid() ) {
    theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
    for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin();
         isim != theEBCaloHits.end(); ++isim){
      EBEnergy_ += isim->energy();
    }
  }

  double EEEnergy_ = 0.;
  if ( EcalHitsEE.isValid() ) {
    theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
    for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin();
         isim != theEECaloHits.end(); ++isim){
      EEEnergy_ += isim->energy();
    }
  }
  
  double ESEnergy_ = 0.;
  if ( EcalHitsES.isValid() ) {
    theESCaloHits.insert(theESCaloHits.end(), EcalHitsES->begin(), EcalHitsES->end());
    for (std::vector<PCaloHit>::iterator isim = theESCaloHits.begin();
         isim != theESCaloHits.end(); ++isim){
      ESEnergy_ += isim->energy();
    }
  }
  
  double etot = EBEnergy_ + EEEnergy_ + ESEnergy_ ;
  double fracEB = 0.0;
  double fracEE = 0.0;
  double fracES = 0.0;
  
  if (etot>0.0) { 
    fracEB  = EBEnergy_/etot; 
    fracEE  = EEEnergy_/etot; 
    fracES  = ESEnergy_/etot; 
  }
  
  if (meEBEnergyFraction_) meEBEnergyFraction_->Fill(fracEB);
  
  if (meEEEnergyFraction_) meEEEnergyFraction_->Fill(fracEE);
  
  if (meESEnergyFraction_) meESEnergyFraction_->Fill(fracES);
}