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/PhotonsWithConversionsAnalyzer.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 "FWCore/Framework/interface/Event.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "FWCore/Framework/interface/ESHandle.h"
00020 #include "FWCore/Framework/interface/MakerMacros.h"
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 #include "FWCore/Utilities/interface/Exception.h"
00023 //
00024 #include "DataFormats/Common/interface/Handle.h"
00025 #include "DataFormats/TrackReco/interface/Track.h"
00026 #include "DataFormats/TrackReco/interface/TrackExtra.h"
00027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00028 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
00029 #include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
00030 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00031 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00032 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00033 //
00034 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00035 #include "MagneticField/Engine/interface/MagneticField.h"
00036 
00037 //
00038 #include "TFile.h"
00039 #include "TH1.h"
00040 #include "TH2.h"
00041 #include "TTree.h"
00042 #include "TVector3.h"
00043 #include "TProfile.h"
00044 //
00045 
00046 
00047 
00048 using namespace std;
00049 
00050 
00051 PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer( const edm::ParameterSet& pset )
00052   {
00053 
00054   photonCollectionProducer_ = pset.getParameter<std::string>("phoProducer");
00055   photonCollection_ = pset.getParameter<std::string>("photonCollection");
00056   //
00057 
00058 
00059 }
00060 
00061 
00062 
00063 PhotonsWithConversionsAnalyzer::~PhotonsWithConversionsAnalyzer() {
00064 
00065 
00066   delete thePhotonMCTruthFinder_;
00067 
00068 }
00069 
00070 
00071 void PhotonsWithConversionsAnalyzer::beginJob()
00072 {
00073 
00074 
00075   nEvt_=0;
00076   nMCPho_=0;
00077   nMatched_=0;
00078 
00079 
00080   thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
00081 
00082 
00083 
00084   edm::Service<TFileService> fs;
00085 
00086 
00088   h_ErecoEMC_  = fs->make<TH1F>("deltaE","    delta(reco-mc) energy",100,0.,2.);
00089   h_deltaPhi_ = fs->make<TH1F>("deltaPhi","  delta(reco-mc) phi",100,-0.1, 0.1);
00090   h_deltaEta_ = fs->make<TH1F>("deltaEta","  delta(reco-mc) eta",100,-0.05, 0.05);
00091 
00093   h_MCphoE_ = fs->make<TH1F>("MCphoE","MC photon energy",100,0.,100.);
00094   h_MCphoPhi_ = fs->make<TH1F>("MCphoPhi","MC photon phi",40,-3.14, 3.14);
00095   h_MCphoEta_ = fs->make<TH1F>("MCphoEta","MC photon eta",40,-3., 3.);
00096 
00097 
00099   h_MCConvE_ = fs->make<TH1F>("MCConvE","MC converted photon energy",100,0.,100.);
00100   h_MCConvPt_ = fs->make<TH1F>("MCConvPt","MC converted photon pt",100,0.,100.);
00101   h_MCConvEta_ = fs->make<TH1F>("MCConvEta","MC converted photon eta",50, 0., 2.5);
00102 
00104   h_scE_ = fs->make<TH1F>("scE","SC Energy ",100,0., 200.);
00105   h_scEt_ = fs->make<TH1F>("scEt","SC Et ",100,0., 200.);
00106   h_scEta_ = fs->make<TH1F>("scEta","SC Eta ",40,-3., 3.);
00107   h_scPhi_ = fs->make<TH1F>("scPhi","SC Phi ",40, -3.14, 3.14);
00108   //
00109   h_phoE_ = fs->make<TH1F>("phoE","Photon Energy ",100,0., 200.);
00110   h_phoEta_ = fs->make<TH1F>("phoEta","Photon Eta ",40,-3., 3.);
00111   h_phoPhi_ = fs->make<TH1F>("phoPhi","Photon  Phi ",40,  -3.14, 3.14);
00112 
00113   // Recontructed tracks from converted photon candidates
00114   h2_tk_nHitsVsR_ = fs->make<TH2F>("tknHitsVsR","Tracks Hits vs R  ", 12,0.,120.,20,0.5, 20.5);
00115   h2_tk_inPtVsR_  = fs->make<TH2F>("tkInPtvsR","Tracks inner Pt vs R  ", 12,0.,120.,100,0., 100.);
00116 
00117 
00118   return ;
00119 }
00120 
00121 
00122 float PhotonsWithConversionsAnalyzer::etaTransformation(  float EtaParticle , float Zvertex)  {
00123 
00124 //---Definitions
00125         const float PI    = 3.1415927;
00126         //UNUSED const float TWOPI = 2.0*PI;
00127 
00128 //---Definitions for ECAL
00129         const float R_ECAL           = 136.5;
00130         const float Z_Endcap         = 328.0;
00131         const float etaBarrelEndcap  = 1.479;
00132 
00133 //---ETA correction
00134 
00135         float Theta = 0.0  ;
00136         float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;
00137 
00138         if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
00139         if(Theta<0.0) Theta = Theta+PI ;
00140         float ETA = - log(tan(0.5*Theta));
00141 
00142         if( fabs(ETA) > etaBarrelEndcap )
00143           {
00144            float Zend = Z_Endcap ;
00145            if(EtaParticle<0.0 )  Zend = -Zend ;
00146            float Zlen = Zend - Zvertex ;
00147            float RR = Zlen/sinh(EtaParticle);
00148            Theta = atan(RR/Zend);
00149            if(Theta<0.0) Theta = Theta+PI ;
00150            ETA = - log(tan(0.5*Theta));
00151           }
00152 //---Return the result
00153         return ETA;
00154 //---end
00155 }
00156 
00157 
00158 
00159 
00160 
00161 
00162 void PhotonsWithConversionsAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& )
00163 {
00164 
00165 
00166   using namespace edm;
00167   const float etaPhiDistance=0.01;
00168   // Fiducial region
00169   //UNUSED const float TRK_BARL =0.9;
00170   //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
00171   //UNUSED const float END_LO = 1.566;
00172   //UNUSED const float END_HI = 2.5;
00173  // Electron mass
00174   //UNUSED const Float_t mElec= 0.000511;
00175 
00176 
00177   nEvt_++;
00178   LogInfo("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
00179   //  LogDebug("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
00180   std::cout << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
00181 
00182 
00184   Handle<reco::PhotonCollection> photonHandle;
00185   e.getByLabel(photonCollectionProducer_, photonCollection_ , photonHandle);
00186   const reco::PhotonCollection photonCollection = *(photonHandle.product());
00187   std::cout  << "ConvertedPhotonAnalyzer  Photons with conversions collection size " << photonCollection.size() << "\n";
00188 
00189 
00191   std::cout  << " ConvertedPhotonAnalyzer Looking for MC truth " << "\n";
00192 
00193   //get simtrack info
00194   std::vector<SimTrack> theSimTracks;
00195   std::vector<SimVertex> theSimVertices;
00196 
00197   edm::Handle<SimTrackContainer> SimTk;
00198   edm::Handle<SimVertexContainer> SimVtx;
00199   e.getByLabel("g4SimHits",SimTk);
00200   e.getByLabel("g4SimHits",SimVtx);
00201 
00202   theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
00203   theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());
00204 
00205   std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks,  theSimVertices);
00206   std::cout << " ConvertedPhotonAnalyzer mcPhotons size " <<  mcPhotons.size() << std::endl;
00207 
00208 
00209 
00210   // Loop over simulated photons
00211   //UNUSED int iDet=0;
00212   //UNUSED int iRadius=-1;
00213   //UNUSED int indPho=0;
00214 
00215   for ( std::vector<PhotonMCTruth>::const_iterator mcPho=mcPhotons.begin(); mcPho !=mcPhotons.end(); mcPho++) {
00216     float mcPhi= (*mcPho).fourMomentum().phi();
00217     float mcEta= (*mcPho).fourMomentum().pseudoRapidity();
00218     mcEta = etaTransformation(mcEta, (*mcPho).primaryVertex().z() );
00219 
00220     if ( (*mcPho).fourMomentum().et() < 20 ) continue;
00221     //    if ( ! (  fabs(mcEta) <= BARL || ( fabs(mcEta) >= END_LO && fabs(mcEta) <=END_HI ) ) ) {
00222     //     continue;
00223     //} // all ecal fiducial region
00224 
00225 
00226 
00227     h_MCphoE_->Fill(  (*mcPho).fourMomentum().e());
00228     h_MCphoEta_->Fill( (*mcPho).fourMomentum().eta());
00229     h_MCphoPhi_->Fill(  (*mcPho).fourMomentum().phi());
00230 
00231 
00232 
00233     // keep only visible conversions
00234     if (  (*mcPho).isAConversion() == 0 ) continue;
00235 
00236 
00237     nMCPho_++;
00238 
00239     h_MCConvEta_ ->Fill ( fabs( (*mcPho).fourMomentum().pseudoRapidity()) - 0.001);
00240 
00241 
00242     bool REJECTED;
00243 
00244 
00246     //std::cout   << " ConvertedPhotonAnalyzer  Starting loop over photon candidates " << "\n";
00247     for( reco::PhotonCollection::const_iterator  iPho = photonCollection.begin(); iPho != photonCollection.end(); iPho++) {
00248       REJECTED=false;
00249 
00250       //      std::cout  << " ConvertedPhotonAnalyzer Reco SC energy " << (*iPho).superCluster()->energy() <<  "\n";
00251 
00252       float phiClu=(*iPho).superCluster()->phi();
00253       float etaClu=(*iPho).superCluster()->eta();
00254       float deltaPhi = phiClu-mcPhi;
00255       float deltaEta = etaClu-mcEta;
00256 
00257 
00258       if ( deltaPhi > Geom::pi()  ) deltaPhi -= Geom::twoPi();
00259       if ( deltaPhi < -Geom::pi() ) deltaPhi += Geom::twoPi();
00260       deltaPhi=pow(deltaPhi,2);
00261       deltaEta=pow(deltaEta,2);
00262       float delta =  deltaPhi+deltaEta ;
00263       if (  delta >= etaPhiDistance  )  REJECTED=true;
00264 
00265 
00266       //      if ( ! (  fabs(etaClu) <= BARL || ( fabs(etaClu) >= END_LO && fabs(etaClu) <=END_HI ) ) ) REJECTED=true;
00267 
00268       if ( REJECTED ) continue;
00269       std::cout << " MATCHED " << std::endl;
00270       nMatched_++;
00271 
00272 
00273       //      std::cout << " ConvertedPhotonAnalyzer Matching candidate " << std::endl;
00274 
00275       // std::cout << " ConvertedPhotonAnalyzer Photons isAconversion " << (*mcPho).isAConversion() << " mcMatchingPhoton energy " <<  (*mcPho).fourMomentum().e()  << " ConvertedPhotonAnalyzer conversion vertex R " <<  (*mcPho).vertex().perp() << " Z " <<  (*mcPho).vertex().z() <<  std::endl;
00276 
00277 
00278       h_ErecoEMC_->Fill(   (*iPho).superCluster()->energy()/(*mcPho).fourMomentum().e());
00279       h_deltaPhi_-> Fill ( (*iPho).superCluster()->position().phi()- mcPhi);
00280       h_deltaEta_-> Fill ( (*iPho).superCluster()->position().eta()- mcEta);
00281 
00282       h_scE_->Fill( (*iPho).superCluster()->energy() );
00283       h_scEt_->Fill( (*iPho).superCluster()->energy()/cosh( (*iPho).superCluster()->position().eta()) );
00284       h_scEta_->Fill( (*iPho).superCluster()->position().eta() );
00285       h_scPhi_->Fill( (*iPho).superCluster()->position().phi() );
00286 
00287 
00288       h_phoE_->Fill( (*iPho).energy() );
00289       h_phoEta_->Fill( (*iPho).eta() );
00290       h_phoPhi_->Fill( (*iPho).phi() );
00291 
00292       if ( !(*iPho).hasConversionTracks() ) continue;
00293       //   std::cout << " This photons has " << (*iPho).conversions().size() << " conversions candidates " << std::endl;
00294       reco::ConversionRefVector conversions = (*iPho).conversions();
00295       //std::vector<reco::ConversionRef> conversions = (*iPho).conversions();
00296 
00297 
00298       for (unsigned int i=0; i<conversions.size(); i++) {
00299         //std::cout << " Conversion candidate Energy " << (*iPho).energy() << " number of tracks " << conversions[i]->nTracks() << std::endl;
00300         std::vector< edm::RefToBase<reco::Track> > tracks = conversions[i]->tracks();
00301 
00302         for (unsigned int i=0; i<tracks.size(); i++) {
00303 
00304           //      std::cout  << " ConvertedPhotonAnalyzer Reco Track charge " <<  tracks[i]->charge() << "  Num of RecHits " << tracks[i]->recHitsSize() << " inner momentum " <<  sqrt ( tracks[i]->innerMomentum().Mag2() )  <<  "\n";
00305 
00306 
00307           h2_tk_nHitsVsR_ -> Fill (  (*mcPho).vertex().perp(), tracks[i]->recHitsSize()   );
00308           h2_tk_inPtVsR_->Fill  (  (*mcPho).vertex().perp(),  sqrt( tracks[i]->innerMomentum().Mag2() ) );
00309         }
00310 
00311       }
00312 
00313 
00314     }
00315 
00316   }
00317 
00318 
00319 }
00320 
00321 
00322 
00323 
00324 void PhotonsWithConversionsAnalyzer::endJob()
00325 {
00326 
00327 
00328 
00329 
00330   //   fOutputFile_->Write() ;
00331   // fOutputFile_->Close() ;
00332 
00333    edm::LogInfo("ConvertedPhotonAnalyzer") << "Analyzed " << nEvt_  << "\n";
00334    // std::cout  << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
00335    std::cout  << "ConvertedPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events " << "\n";
00336 
00337    return ;
00338 }
00339 
00340 
00341
/data/refman/pasoursint/CMSSW_5_3_6/src/RecoEgamma/Examples/plugins/PhotonsWithConversionsAnalyzer.cc
