/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEgamma/Examples/plugins/MCElectronAnalyzer.cc

Go to the documentation of this file.

#include <iostream>
//
#include "RecoEgamma/Examples/plugins/MCElectronAnalyzer.h"
#include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruthFinder.h"
#include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
//
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
//
#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
//
#include "DataFormats/Common/interface/Handle.h"
//
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
//
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"

//
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TTree.h"
#include "TVector3.h"
#include "TProfile.h"
//



using namespace std;


MCElectronAnalyzer::MCElectronAnalyzer( const edm::ParameterSet& pset )
   : fOutputFileName_( pset.getUntrackedParameter<string>("HistOutFile",std::string("TestConversions.root")) ),
     fOutputFile_(0)
{


}



MCElectronAnalyzer::~MCElectronAnalyzer() {


  delete theElectronMCTruthFinder_;

}


void MCElectronAnalyzer::beginJob()
{


  nEvt_=0;

  theElectronMCTruthFinder_ = new ElectronMCTruthFinder();


  fOutputFile_   = new TFile( fOutputFileName_.c_str(), "RECREATE" ) ;

  h_MCEleE_ = new TH1F("MCEleE","MC ele energy",100,0.,200.);
  h_MCElePhi_ = new TH1F("MCElePhi","MC ele phi",40,-3.14, 3.14);
  h_MCEleEta_ = new TH1F("MCEleEta","MC ele eta",40,-3., 3.);
  h_BremFrac_ = new TH1F("bremFrac","brem frac ", 100, 0., 1.);
  h_BremEnergy_ = new TH1F("BremE","Brem energy",100,0.,200.);

  p_BremVsR_ = new TProfile("BremVsR", " Mean Brem energy vs R ", 48, 0., 120.);
  p_BremVsEta_ = new TProfile("BremVsEta", " Mean Brem energy vs Eta ", 50, -2.5, 2.5);




  return ;
}


float MCElectronAnalyzer::etaTransformation(  float EtaParticle , float Zvertex)  {

//---Definitions
        const float PI    = 3.1415927;
        //const float TWOPI = 2.0*PI;

//---Definitions for ECAL
        const float R_ECAL           = 136.5;
        const float Z_Endcap         = 328.0;
        const float etaBarrelEndcap  = 1.479;

//---ETA correction

        float Theta = 0.0  ;
        float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;

        if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
        if(Theta<0.0) Theta = Theta+PI ;
        float ETA = - log(tan(0.5*Theta));

        if( std::abs(ETA) > etaBarrelEndcap )
          {
           float Zend = Z_Endcap ;
           if(EtaParticle<0.0 )  Zend = -Zend ;
           float Zlen = Zend - Zvertex ;
           float RR = Zlen/sinh(EtaParticle);
           Theta = atan(RR/Zend);
           if(Theta<0.0) Theta = Theta+PI ;
           ETA = - log(tan(0.5*Theta));
          }
//---Return the result
        return ETA;
//---end
}

float MCElectronAnalyzer::phiNormalization(float & phi)
{
//---Definitions
 const float PI    = 3.1415927;
 const float TWOPI = 2.0*PI;


 if(phi >  PI) {phi = phi - TWOPI;}
 if(phi < -PI) {phi = phi + TWOPI;}

 //  cout << " Float_t PHInormalization out " << PHI << endl;
 return phi;

}





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


  using namespace edm;
  //const float etaPhiDistance=0.01;
  // Fiducial region
  //const float TRK_BARL =0.9;
  //const float BARL = 1.4442; // DAQ TDR p.290
  //const float END_LO = 1.566;
  //const float END_HI = 2.5;
  // Electron mass
  //const Float_t mElec= 0.000511;


  nEvt_++;
  LogInfo("MCElectronAnalyzer") << "MCElectronAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
  //  LogDebug("MCElectronAnalyzer") << "MCElectronAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
  std::cout << "MCElectronAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";

  std::cout  << " MCElectronAnalyzer Looking for MC truth " << "\n";

  //get simtrack info
  std::vector<SimTrack> theSimTracks;
  std::vector<SimVertex> theSimVertices;

  edm::Handle<SimTrackContainer> SimTk;
  edm::Handle<SimVertexContainer> SimVtx;
  e.getByLabel("g4SimHits",SimTk);
  e.getByLabel("g4SimHits",SimVtx);

  theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
  theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());
  std::cout << " MCElectronAnalyzer This Event has " <<  theSimTracks.size() << " sim tracks " << std::endl;
  std::cout << " MCElectronAnalyzer This Event has " <<  theSimVertices.size() << " sim vertices " << std::endl;
  if (  ! theSimTracks.size() ) std::cout << " Event number " << e.id() << " has NO sim tracks " << std::endl;


  std::vector<ElectronMCTruth> MCElectronctrons=theElectronMCTruthFinder_->find (theSimTracks,  theSimVertices);
  std::cout << " MCElectronAnalyzer MCElectronctrons size " <<  MCElectronctrons.size() << std::endl;

  for ( std::vector<ElectronMCTruth>::const_iterator iEl=MCElectronctrons.begin(); iEl !=MCElectronctrons.end(); ++iEl ){

    h_MCEleE_->Fill  ( (*iEl).fourMomentum().e() );
    h_MCEleEta_->Fill  ( (*iEl).fourMomentum().pseudoRapidity() );
    h_MCElePhi_->Fill  ( (*iEl).fourMomentum().phi() );

    float totBrem=0 ;
    unsigned int iBrem ;
    for ( iBrem=0; iBrem < (*iEl).bremVertices().size(); ++iBrem ) {
      float rBrem= (*iEl).bremVertices()[iBrem].perp();
      float etaBrem= (*iEl).bremVertices()[iBrem].eta();
      if ( rBrem < 120 ) {
        totBrem +=  (*iEl).bremMomentum()[iBrem].e();
        p_BremVsR_ ->Fill ( rBrem, (*iEl).bremMomentum()[iBrem].e() );
        p_BremVsEta_ ->Fill ( etaBrem, (*iEl).bremMomentum()[iBrem].e() );

      }
    }


    h_BremFrac_->Fill( totBrem/(*iEl).fourMomentum().e() );
    h_BremEnergy_->Fill (  totBrem  );




  }




}




void MCElectronAnalyzer::endJob()
{




   fOutputFile_->Write() ;
   fOutputFile_->Close() ;

   edm::LogInfo("MCElectronAnalyzer") << "Analyzed " << nEvt_  << "\n";
   std::cout  << "MCElectronAnalyzer::endJob Analyzed " << nEvt_ << " events " << "\n";

   return ;
}