EcalPreshowerSimHitsValidation.cc

Go to the documentation of this file.

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

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

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

EcalPreshowerSimHitsValidation::EcalPreshowerSimHitsValidation(const edm::ParameterSet& ps):

  HepMCLabel(ps.getParameter<std::string>("moduleLabelMC")),
  g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
  EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
  ESHitsCollection(ps.getParameter<std::string>("ESHitsCollection")){

  HepMCToken   = consumes <edm::HepMCProduct>      (HepMCLabel);
  EEHitsToken  = consumes <edm::PCaloHitContainer> (edm::InputTag(std::string(g4InfoLabel),std::string(EEHitsCollection)));
  ESHitsToken  = consumes <edm::PCaloHitContainer> (edm::InputTag(std::string(g4InfoLabel),std::string(ESHitsCollection)));
  // verbosity switch
  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
 
  // get hold of back-end interface
  dbe_ = nullptr;
  dbe_ = edm::Service<DQMStore>().operator->();
  if ( dbe_ ) {
    if ( verbose_ ) { dbe_->setVerbose(1); } 
    else            { dbe_->setVerbose(0); }
  }

  if ( dbe_ ) {
    if ( verbose_ ) dbe_->showDirStructure();
  }


  menESHits1zp_ = nullptr;     
  menESHits2zp_ = nullptr;     
  menESHits1zm_ = nullptr;     
  menESHits2zm_ = nullptr;     
                                    
  meESEnergyHits1zp_ = nullptr;
  meESEnergyHits2zp_ = nullptr;
  meESEnergyHits1zm_ = nullptr;
  meESEnergyHits2zm_ = nullptr;

  meEShitLog10Energy_       = nullptr;
  meEShitLog10EnergyNorm_   = nullptr;

  meE1alphaE2zp_ = nullptr;
  meE1alphaE2zm_ = nullptr;
  meEEoverESzp_  = nullptr;
  meEEoverESzm_  = nullptr;

  me2eszpOver1eszp_ = nullptr;
  me2eszmOver1eszm_ = nullptr;


  Char_t histo[200];
 
  if ( dbe_ ) {
    dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
  
    sprintf (histo, "ES hits layer 1 multiplicity z+" ) ;
    menESHits1zp_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;

    sprintf (histo, "ES hits layer 2 multiplicity z+" ) ;
    menESHits2zp_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;

    sprintf (histo, "ES hits layer 1 multiplicity z-" ) ;
    menESHits1zm_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;

    sprintf (histo, "ES hits layer 2 multiplicity z-" ) ;
    menESHits2zm_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;

    sprintf (histo, "ES hits energy layer 1 z+" ) ;
    meESEnergyHits1zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;

    sprintf (histo, "ES hits energy layer 2 z+" ) ;
    meESEnergyHits2zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;

    sprintf (histo, "ES hits energy layer 1 z-" ) ;
    meESEnergyHits1zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;

    sprintf (histo, "ES hits energy layer 2 z-" ) ;
    meESEnergyHits2zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;

    sprintf (histo, "ES hits log10energy spectrum" );
    meEShitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);

    sprintf (histo, "ES hits log10energy spectrum vs normalized energy" );
    meEShitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);

    sprintf (histo, "ES E1+07E2 z+" ) ;
    meE1alphaE2zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05);

    sprintf (histo, "ES E1+07E2 z-" ) ;
    meE1alphaE2zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05);

    sprintf (histo, "EE vs ES z+" ) ;
    meEEoverESzp_ = dbe_->bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);

    sprintf (histo, "EE vs ES z-" ) ;
    meEEoverESzm_ = dbe_->bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);

    sprintf (histo, "ES ene2oEne1 z+" ) ;
    me2eszpOver1eszp_ = dbe_->book1D(histo, histo, 50, 0., 10.);

    sprintf (histo, "ES ene2oEne1 z-" ) ;
    me2eszmOver1eszm_ = dbe_->book1D(histo, histo, 50, 0., 10.);
  }
}

EcalPreshowerSimHitsValidation::~EcalPreshowerSimHitsValidation(){

}

void EcalPreshowerSimHitsValidation::beginJob(){

}

void EcalPreshowerSimHitsValidation::endJob(){

}

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

  edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
  
  edm::Handle<edm::HepMCProduct> MCEvt;
  e.getByToken(HepMCToken, MCEvt);

  edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
  e.getByToken(EEHitsToken,EcalHitsEE);

  edm::Handle<edm::PCaloHitContainer> EcalHitsES;
  e.getByToken(ESHitsToken,EcalHitsES);

  std::vector<PCaloHit> theEECaloHits;  
  if( EcalHitsEE.isValid() ) {
    theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
  }

  std::vector<PCaloHit> theESCaloHits;
  if( EcalHitsES.isValid() ) {
    theESCaloHits.insert(theESCaloHits.end(), EcalHitsES->begin(), EcalHitsES->end());
  }

  double ESEnergy_ = 0.;
  //std::map<unsigned int, std::vector<PCaloHit>,std::less<unsigned int> > CaloHitMap;

  // endcap
  double EEetzp_ = 0.;
  double EEetzm_ = 0.;
  for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin(); isim != theEECaloHits.end(); ++isim){
    EEDetId eeid (isim->id()) ;
    if (eeid.zside() > 0 ) EEetzp_ += isim->energy();
    if (eeid.zside() < 0 ) EEetzm_ += isim->energy();
  }
  
  
  uint32_t nESHits1zp = 0;
  uint32_t nESHits1zm = 0;
  uint32_t nESHits2zp = 0;
  uint32_t nESHits2zm = 0;
  double ESet1zp_ = 0.;
  double ESet2zp_ = 0.;
  double ESet1zm_ = 0.;
  double ESet2zm_ = 0.;
  std::vector<double> econtr(140, 0. );
  
  for (std::vector<PCaloHit>::iterator isim = theESCaloHits.begin();
       isim != theESCaloHits.end(); ++isim){
    //CaloHitMap[(*isim).id()].push_back((*isim));
    
    ESDetId esid (isim->id()) ;
    
    LogDebug("HitInfo")
      << " CaloHit " << isim->getName() << "\n" 
      << " DetID = "<<isim->id()<< " ESDetId: z side " << esid.zside() << "  plane " << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
      << " Time = " << isim->time() << "\n"
      << " Track Id = " << isim->geantTrackId() << "\n"
      << " Energy = " << isim->energy();
    
    ESEnergy_ += isim->energy();
    if( isim->energy() > 0 ) {
      meEShitLog10Energy_->Fill(log10(isim->energy()));
      int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
      if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
    }

    if (esid.plane() == 1 ) { 
      if (esid.zside() > 0 ) {
    nESHits1zp++ ;
    ESet1zp_ += isim->energy();
    if (meESEnergyHits1zp_) meESEnergyHits1zp_->Fill(isim->energy()) ; 
      }
      else if (esid.zside() < 0 ) {
    nESHits1zm++ ; 
    ESet1zm_ += isim->energy();
    if (meESEnergyHits1zm_) meESEnergyHits1zm_->Fill(isim->energy()) ; 
      }
    }
    else if (esid.plane() == 2 ) {
      if (esid.zside() > 0 ) {
    nESHits2zp++ ; 
    ESet2zp_ += isim->energy();
    if (meESEnergyHits2zp_) meESEnergyHits2zp_->Fill(isim->energy()) ; 
      }
      else if (esid.zside() < 0 ) {
    nESHits2zm++ ; 
    ESet2zm_ += isim->energy();
    if (meESEnergyHits2zm_) meESEnergyHits2zm_->Fill(isim->energy()) ; 
      }
    }
    
  }
  
  if (menESHits1zp_) menESHits1zp_->Fill(nESHits1zp);
  if (menESHits1zm_) menESHits1zm_->Fill(nESHits1zm);
  
  if (menESHits2zp_) menESHits2zp_->Fill(nESHits2zp);
  if (menESHits2zm_) menESHits2zm_->Fill(nESHits2zm);
  
  if( meEShitLog10EnergyNorm_ && ESEnergy_ != 0 ) {
    for( int i=0; i<140; i++ ) {
      meEShitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/ESEnergy_ );
    }
  }


  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));
    }

    if ( heta > 1.653 && heta < 2.6 ) {

      if (meE1alphaE2zp_) meE1alphaE2zp_->Fill(ESet1zp_+0.7*ESet2zp_);
      if (meEEoverESzp_)  meEEoverESzp_ ->Fill((ESet1zp_+0.7*ESet2zp_)/0.00009, EEetzp_);
      if ((me2eszpOver1eszp_) && (ESet1zp_ != 0.)) me2eszpOver1eszp_->Fill(ESet2zp_/ESet1zp_);
    }
    if ( heta < -1.653 && heta > -2.6 ) {
      if (meE1alphaE2zm_) meE1alphaE2zm_->Fill(ESet1zm_+0.7*ESet2zm_);
      if (meEEoverESzm_)  meEEoverESzm_ ->Fill((ESet1zm_+0.7*ESet2zm_)/0.00009, EEetzm_);
      if ((me2eszmOver1eszm_) && (ESet1zm_ != 0.)) me2eszmOver1eszm_->Fill(ESet2zm_/ESet1zm_);
    }
  }
}


//  LocalWords:  EcalSimHitsValidation