Reference Manual

Go to the documentation of this file.
00001 #include <iostream>
00002 //
00003 #include "FWCore/ServiceRegistry/interface/Service.h"
00004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00005 //
00006 #include "RecoEgamma/Examples/plugins/MCPhotonAnalyzer.h"
00007 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
00008 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
00009 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
00010 //
00011 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00012 #include "SimDataFormats/Track/interface/SimTrack.h"
00013 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00014 #include "SimDataFormats/Vertex/interface/SimVertex.h"
00015 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00016 //
00017 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
00018 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
00019 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
00020 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
00021 //
00022 #include "DataFormats/Common/interface/Handle.h"
00023 //
00024 #include "FWCore/Framework/interface/Event.h"
00025 #include "FWCore/Framework/interface/EventSetup.h"
00026 #include "FWCore/Framework/interface/ESHandle.h"
00027 #include "FWCore/Framework/interface/MakerMacros.h"
00028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00029 #include "FWCore/Utilities/interface/Exception.h"
00030 //
00031 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00032 #include "MagneticField/Engine/interface/MagneticField.h"
00033 
00034 //
00035 #include "TFile.h"
00036 #include "TH1.h"
00037 #include "TH2.h"
00038 #include "TTree.h"
00039 #include "TVector3.h"
00040 #include "TProfile.h"
00041 //
00042 
00043 
00044 
00045 using namespace std;
00046 
00047 
00048 MCPhotonAnalyzer::MCPhotonAnalyzer( const edm::ParameterSet& pset ){}
00049 
00050 
00051 
00052 MCPhotonAnalyzer::~MCPhotonAnalyzer() {
00053 
00054 
00055   delete thePhotonMCTruthFinder_;
00056 
00057 }
00058 
00059 
00060 void MCPhotonAnalyzer::beginJob()
00061 {
00062 
00063 
00064   nEvt_=0;
00065 
00066   thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
00067 
00068   edm::Service<TFileService> fs;
00069 
00070 
00072   h_MCPhoE_ = fs->make<TH1F>("MCPhoE","MC photon energy",100,0.,100.);
00073   h_MCPhoPhi_ = fs->make<TH1F>("MCPhoPhi","MC photon phi",40,-3.14, 3.14);
00074   h_MCPhoEta_ = fs->make<TH1F>("MCPhoEta","MC photon eta",25,0., 2.5);
00075   h_MCPhoEta1_ = fs->make<TH1F>("MCPhoEta1","MC photon eta",40,-3., 3.);
00076   h_MCPhoEta2_ = fs->make<TH1F>("MCPhoEta2","MC photon eta",40,-3., 3.);
00077   h_MCPhoEta3_ = fs->make<TH1F>("MCPhoEta3","MC photon eta",40,-3., 3.);
00078   h_MCPhoEta4_ = fs->make<TH1F>("MCPhoEta4","MC photon eta",40,-3., 3.);
00080   h_MCConvPhoE_ = fs->make<TH1F>("MCConvPhoE","MC converted photon energy",100,0.,100.);
00081   h_MCConvPhoPhi_ = fs->make<TH1F>("MCConvPhoPhi","MC converted photon phi",40,-3.14, 3.14);
00082   h_MCConvPhoEta_ = fs->make<TH1F>("MCConvPhoEta","MC converted photon eta",25,0., 2.5);
00083   h_MCConvPhoR_ = fs->make<TH1F>("MCConvPhoR","MC converted photon R",120,0.,120.);
00084 
00085   h_MCConvPhoREta1_ = fs->make<TH1F>("MCConvPhoREta1","MC converted photon R",120,0.,120.);
00086   h_MCConvPhoREta2_ = fs->make<TH1F>("MCConvPhoREta2","MC converted photon R",120,0.,120.);
00087   h_MCConvPhoREta3_ = fs->make<TH1F>("MCConvPhoREta3","MC converted photon R",120,0.,120.);
00088   h_MCConvPhoREta4_ = fs->make<TH1F>("MCConvPhoREta4","MC converted photon R",120,0.,120.);
00089 
00090   h_convFracEta1_ = fs->make<TH1F>("convFracEta1","Integrated(R) fraction of conversion |eta|=0.2",120,0.,120.);
00091   h_convFracEta2_ = fs->make<TH1F>("convFracEta2","Integrated(R) fraction of conversion |eta|=0.9",120,0.,120.);
00092   h_convFracEta3_ = fs->make<TH1F>("convFracEta3","Integrated(R) fraction of conversion |eta|=1.5",120,0.,120.);
00093   h_convFracEta4_ = fs->make<TH1F>("convFracEta4","Integrated(R) fraction of conversion |eta|=2.0",120,0.,120.);
00095   h_MCConvPhoTwoTracksE_ = fs->make<TH1F>("MCConvPhoTwoTracksE","MC converted photon with 2 tracks  energy",100,0.,100.);
00096   h_MCConvPhoTwoTracksPhi_ = fs->make<TH1F>("MCConvPhoTwoTracksPhi","MC converted photon 2 tracks  phi",40,-3.14, 3.14);
00097   h_MCConvPhoTwoTracksEta_ = fs->make<TH1F>("MCConvPhoTwoTracksEta","MC converted photon 2 tracks eta",40,-3., 3.);
00098   h_MCConvPhoTwoTracksR_ = fs->make<TH1F>("MCConvPhoTwoTracksR","MC converted photon 2 tracks eta",48,0.,120.);
00099   // conversions with one track
00100   h_MCConvPhoOneTrackE_ = fs->make<TH1F>("MCConvPhoOneTrackE","MC converted photon with 1 track  energy",100,0.,100.);
00101   h_MCConvPhoOneTrackPhi_ = fs->make<TH1F>("MCConvPhoOneTrackPhi","MC converted photon 1 track  phi",40,-3.14, 3.14);
00102   h_MCConvPhoOneTrackEta_ = fs->make<TH1F>("MCConvPhoOneTrackEta","MC converted photon 1 track eta",40,-3., 3.);
00103   h_MCConvPhoOneTrackR_ = fs->make<TH1F>("MCConvPhoOneTrackR","MC converted photon 1 track eta",48,0.,120.);
00104 
00106   h_MCEleE_ = fs->make<TH1F>("MCEleE","MC ele energy",100,0.,200.);
00107   h_MCElePhi_ = fs->make<TH1F>("MCElePhi","MC ele phi",40,-3.14, 3.14);
00108   h_MCEleEta_ = fs->make<TH1F>("MCEleEta","MC ele eta",40,-3., 3.);
00109   h_BremFrac_ = fs->make<TH1F>("bremFrac","brem frac ", 100, 0., 1.);
00110   h_BremEnergy_ = fs->make<TH1F>("BremE","Brem energy",100,0.,200.);
00111   h_EleEvsPhoE_ =  fs->make<TH2F>("eleEvsPhoE","eleEvsPhoE",100,0.,200.,100,0.,200.);
00112   h_bremEvsEleE_ = fs->make<TH2F>("bremEvsEleE","bremEvsEleE",100,0.,200.,100,0.,200.);
00113 
00114   p_BremVsR_ = fs->make<TProfile>("BremVsR", " Mean Brem Energy vs R ", 48, 0., 120.);
00115   p_BremVsEta_ = fs->make<TProfile>("BremVsEta", " Mean Brem Energy vs Eta ", 50, -2.5, 2.5);
00116 
00117   p_BremVsConvR_ = fs->make<TProfile>("BremVsConvR", " Mean Brem Fraction vs conversion R ", 48, 0., 120.);
00118   p_BremVsConvEta_ = fs->make<TProfile>("BremVsConvEta", " Mean Brem Fraction vs converion Eta ", 50, -2.5, 2.5);
00119 
00120   h_bremFracVsConvR_ = fs->make<TH2F>("bremFracVsConvR","brem Fraction vs conversion R",60,0.,120.,100,0.,1.);
00121 
00122   return ;
00123 }
00124 
00125 
00126 float MCPhotonAnalyzer::etaTransformation(  float EtaParticle , float Zvertex)  {
00127 
00128 //---Definitions
00129         const float PI    = 3.1415927;
00130 //      const float TWOPI = 2.0*PI;
00131 
00132 //---Definitions for ECAL
00133         const float R_ECAL           = 136.5;
00134         const float Z_Endcap         = 328.0;
00135         const float etaBarrelEndcap  = 1.479;
00136 
00137 //---ETA correction
00138 
00139         float Theta = 0.0  ;
00140         float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;
00141 
00142         if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
00143         if(Theta<0.0) Theta = Theta+PI ;
00144         float ETA = - log(tan(0.5*Theta));
00145 
00146         if( fabs(ETA) > etaBarrelEndcap )
00147           {
00148            float Zend = Z_Endcap ;
00149            if(EtaParticle<0.0 )  Zend = -Zend ;
00150            float Zlen = Zend - Zvertex ;
00151            float RR = Zlen/sinh(EtaParticle);
00152            Theta = atan(RR/Zend);
00153            if(Theta<0.0) Theta = Theta+PI ;
00154            ETA = - log(tan(0.5*Theta));
00155           }
00156 //---Return the result
00157         return ETA;
00158 //---end
00159 }
00160 
00161 float MCPhotonAnalyzer::phiNormalization(float & phi)
00162 {
00163 //---Definitions
00164  const float PI    = 3.1415927;
00165  const float TWOPI = 2.0*PI;
00166 
00167 
00168  if(phi >  PI) {phi = phi - TWOPI;}
00169  if(phi < -PI) {phi = phi + TWOPI;}
00170 
00171  //  cout << " Float_t PHInormalization out " << PHI << endl;
00172  return phi;
00173 
00174 }
00175 
00176 
00177 
00178 
00179 
00180 void MCPhotonAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& )
00181 {
00182 
00183 
00184   using namespace edm;
00185   //UNUSED const float etaPhiDistance=0.01;
00186   // Fiducial region
00187   //UNUSED const float TRK_BARL =0.9;
00188   //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
00189   //UNUSED const float END_LO = 1.566;
00190   //UNUSED const float END_HI = 2.5;
00191  // Electron mass
00192   //UNUSED const Float_t mElec= 0.000511;
00193 
00194 
00195   nEvt_++;
00196   LogInfo("mcEleAnalyzer") << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
00197   //  LogDebug("MCPhotonAnalyzer") << "MCPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
00198   std::cout << "MCPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
00199 
00200 
00201 
00203   std::cout  << " MCPhotonAnalyzer Looking for MC truth " << "\n";
00204 
00205   //get simtrack info
00206   std::vector<SimTrack> theSimTracks;
00207   std::vector<SimVertex> theSimVertices;
00208 
00209   edm::Handle<SimTrackContainer> SimTk;
00210   edm::Handle<SimVertexContainer> SimVtx;
00211   e.getByLabel("g4SimHits",SimTk);
00212   e.getByLabel("g4SimHits",SimVtx);
00213 
00214   theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
00215   theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());
00216   std::cout << " MCPhotonAnalyzer This Event has " <<  theSimTracks.size() << " sim tracks " << std::endl;
00217   std::cout << " MCPhotonAnalyzer This Event has " <<  theSimVertices.size() << " sim vertices " << std::endl;
00218   if (  ! theSimTracks.size() ) std::cout << " Event number " << e.id() << " has NO sim tracks " << std::endl;
00219 
00220 
00221   std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks,  theSimVertices);
00222   std::cout << " MCPhotonAnalyzer mcPhotons size " <<  mcPhotons.size() << std::endl;
00223 
00224 
00225 
00226   for ( std::vector<PhotonMCTruth>::const_iterator iPho=mcPhotons.begin(); iPho !=mcPhotons.end(); ++iPho ){
00227 
00228     if ( (*iPho).fourMomentum().e() < 35 ) continue;
00229 
00230     h_MCPhoE_->Fill  ( (*iPho).fourMomentum().e() );
00231     //    float correta = etaTransformation( (*iPho).fourMomentum().pseudoRapidity(),  (*iPho).primaryVertex().z() );
00232     float Theta= (*iPho).fourMomentum().theta();
00233     float correta = - log(tan(0.5*Theta));
00234     correta = etaTransformation( correta,  (*iPho).primaryVertex().z() );
00235          //h_MCPhoEta_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00236     h_MCPhoEta_->Fill  ( fabs(correta)-0.001 );
00237     h_MCPhoPhi_->Fill  ( (*iPho).fourMomentum().phi() );
00238 
00239     /*
00240     if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15  )
00241       h_MCPhoEta1_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00242     if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85  )
00243       h_MCPhoEta2_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00244     if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55  )
00245       h_MCPhoEta3_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00246     if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95  )
00247       h_MCPhoEta4_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00248     */
00249 
00250 
00251     if ( fabs(correta ) <= 0.3 &&  fabs(correta ) >0.2  )
00252       h_MCPhoEta1_->Fill  ( correta );
00253     if ( fabs(correta ) <= 1.00  &&  fabs( correta ) >0.9  )
00254       h_MCPhoEta2_->Fill  ( correta );
00255     if ( fabs( correta ) <= 1.6 &&  fabs(correta ) >1.5  )
00256       h_MCPhoEta3_->Fill  ( correta );
00257     if ( fabs(correta ) <= 2.  &&  fabs(correta ) >1.9  )
00258       h_MCPhoEta4_->Fill  ( correta );
00259 
00260 
00261 
00262     //    if ( (*iPho).isAConversion()  && (*iPho).vertex().perp()< 10 ) {
00263         if ( (*iPho).isAConversion() ) {
00264 
00265 
00266 
00267       h_MCConvPhoE_->Fill  ( (*iPho).fourMomentum().e() );
00268       //      h_MCConvPhoEta_->Fill  ( (*iPho).fourMomentum().pseudoRapidity() );
00269 
00270       h_MCConvPhoEta_->Fill  ( fabs(correta)-0.001 );
00271       h_MCConvPhoPhi_->Fill  ( (*iPho).fourMomentum().phi() );
00272       h_MCConvPhoR_->Fill  ( (*iPho).vertex().perp() );
00273 
00274       /*
00275       if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15  )
00276         h_MCConvPhoREta1_->Fill  ( (*iPho).vertex().perp() );
00277       if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85  )
00278         h_MCConvPhoREta2_->Fill  ( (*iPho).vertex().perp() );
00279       if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55  )
00280         h_MCConvPhoREta3_->Fill  ( (*iPho).vertex().perp() );
00281       if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05  &&  fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95  )
00282         h_MCConvPhoREta4_->Fill  ( (*iPho).vertex().perp() );
00283       */
00284 
00285 
00286       if ( fabs(correta ) <= 0.3 &&  fabs(correta ) >0.2  )
00287         h_MCConvPhoREta1_->Fill  ( (*iPho).vertex().perp() );
00288       if ( fabs(correta ) <= 1.  &&  fabs( correta ) >0.9  )
00289         h_MCConvPhoREta2_->Fill  ( (*iPho).vertex().perp() );
00290       if ( fabs( correta ) <= 1.6  &&  fabs(correta ) >1.5  )
00291         h_MCConvPhoREta3_->Fill  ( (*iPho).vertex().perp() );
00292       if ( fabs(correta ) <= 2  &&  fabs(correta ) >1.9  )
00293         h_MCConvPhoREta4_->Fill  ( (*iPho).vertex().perp() );
00294 
00295 
00296 
00297 
00298 
00299     } // end conversions
00300 
00301 
00302 
00303 
00304 
00305   }   
00306 
00307 
00308 
00309 
00310 }
00311 
00312 
00313 
00314 
00315 void MCPhotonAnalyzer::endJob()
00316 {
00317 
00318 
00319   Double_t s1=0;
00320   Double_t s2=0;
00321   Double_t s3=0;
00322   Double_t s4=0;
00323   int e1=0;
00324   int e2=0;
00325   int e3=0;
00326   int e4=0;
00327 
00328   Double_t nTotEta1 = h_MCPhoEta1_->GetEntries();
00329   Double_t nTotEta2 = h_MCPhoEta2_->GetEntries();
00330   Double_t nTotEta3 = h_MCPhoEta3_->GetEntries();
00331   Double_t nTotEta4 = h_MCPhoEta4_->GetEntries();
00332 
00333   for ( int i=1; i<=120; ++i) {
00334   e1 = (int)h_MCConvPhoREta1_->GetBinContent(i);
00335   e2 = (int)h_MCConvPhoREta2_->GetBinContent(i);
00336   e3 = (int)h_MCConvPhoREta3_->GetBinContent(i);
00337   e4 = (int)h_MCConvPhoREta4_->GetBinContent(i);
00338   s1+=e1;
00339   s2+=e2;
00340   s3+=e3;
00341   s4+=e4;
00342   h_convFracEta1_->SetBinContent(i,s1*100/nTotEta1);
00343   h_convFracEta2_->SetBinContent(i,s2*100/nTotEta2);
00344   h_convFracEta3_->SetBinContent(i,s3*100/nTotEta3);
00345   h_convFracEta4_->SetBinContent(i,s4*100/nTotEta4);
00346 
00347 
00348 
00349 
00350 
00351 }
00352 
00353 
00354 
00355    edm::LogInfo("MCPhotonAnalyzer") << "Analyzed " << nEvt_  << "\n";
00356    std::cout  << "MCPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events " << "\n";
00357 
00358    return ;
00359 }
00360
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoEgamma/Examples/plugins/MCPhotonAnalyzer.cc
