CMS 3D CMS Logo
PhotonsWithConversionsAnalyzer.cc
Go to the documentation of this file.
1 #include <iostream>
2 //
3 #include "FWCore/ServiceRegistry/interface/Service.h"
4 #include "CommonTools/UtilAlgos/interface/TFileService.h"
5 //
6 #include "RecoEgamma/Examples/plugins/PhotonsWithConversionsAnalyzer.h"
7 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
8 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
9 #include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
10 //
11 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
12 #include "SimDataFormats/Track/interface/SimTrack.h"
13 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
14 #include "SimDataFormats/Vertex/interface/SimVertex.h"
15 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
16 //
17 #include "FWCore/Framework/interface/Event.h"
18 #include "FWCore/Framework/interface/EventSetup.h"
19 #include "FWCore/Framework/interface/ESHandle.h"
20 #include "FWCore/Framework/interface/MakerMacros.h"
21 #include "FWCore/MessageLogger/interface/MessageLogger.h"
22 #include "FWCore/Utilities/interface/Exception.h"
23 //
24 #include "DataFormats/Common/interface/Handle.h"
25 #include "DataFormats/TrackReco/interface/Track.h"
26 #include "DataFormats/TrackReco/interface/TrackExtra.h"
27 #include "DataFormats/TrackReco/interface/TrackFwd.h"
28 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
29 #include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
30 #include "DataFormats/EgammaCandidates/interface/Photon.h"
31 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
32 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
33 //
34 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
35 #include "MagneticField/Engine/interface/MagneticField.h"
36 
37 //
38 #include "TFile.h"
39 #include "TH1.h"
40 #include "TH2.h"
41 #include "TTree.h"
42 #include "TVector3.h"
43 #include "TProfile.h"
44 //
45 
46 using namespace std;
47 
48 PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer(const edm::ParameterSet& pset) {
49  photonCollectionProducer_ = pset.getParameter<std::string>("phoProducer");
50  photonCollection_ = pset.getParameter<std::string>("photonCollection");
51  //
52 }
53 
54 PhotonsWithConversionsAnalyzer::~PhotonsWithConversionsAnalyzer() { delete thePhotonMCTruthFinder_; }
55 
56 void PhotonsWithConversionsAnalyzer::beginJob() {
57  nEvt_ = 0;
58  nMCPho_ = 0;
59  nMatched_ = 0;
60 
61  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
62 
63  edm::Service<TFileService> fs;
64 
66  h_ErecoEMC_ = fs->make<TH1F>("deltaE", " delta(reco-mc) energy", 100, 0., 2.);
67  h_deltaPhi_ = fs->make<TH1F>("deltaPhi", " delta(reco-mc) phi", 100, -0.1, 0.1);
68  h_deltaEta_ = fs->make<TH1F>("deltaEta", " delta(reco-mc) eta", 100, -0.05, 0.05);
69 
71  h_MCphoE_ = fs->make<TH1F>("MCphoE", "MC photon energy", 100, 0., 100.);
72  h_MCphoPhi_ = fs->make<TH1F>("MCphoPhi", "MC photon phi", 40, -3.14, 3.14);
73  h_MCphoEta_ = fs->make<TH1F>("MCphoEta", "MC photon eta", 40, -3., 3.);
74 
76  h_MCConvE_ = fs->make<TH1F>("MCConvE", "MC converted photon energy", 100, 0., 100.);
77  h_MCConvPt_ = fs->make<TH1F>("MCConvPt", "MC converted photon pt", 100, 0., 100.);
78  h_MCConvEta_ = fs->make<TH1F>("MCConvEta", "MC converted photon eta", 50, 0., 2.5);
79 
81  h_scE_ = fs->make<TH1F>("scE", "SC Energy ", 100, 0., 200.);
82  h_scEt_ = fs->make<TH1F>("scEt", "SC Et ", 100, 0., 200.);
83  h_scEta_ = fs->make<TH1F>("scEta", "SC Eta ", 40, -3., 3.);
84  h_scPhi_ = fs->make<TH1F>("scPhi", "SC Phi ", 40, -3.14, 3.14);
85  //
86  h_phoE_ = fs->make<TH1F>("phoE", "Photon Energy ", 100, 0., 200.);
87  h_phoEta_ = fs->make<TH1F>("phoEta", "Photon Eta ", 40, -3., 3.);
88  h_phoPhi_ = fs->make<TH1F>("phoPhi", "Photon Phi ", 40, -3.14, 3.14);
89 
90  // Recontructed tracks from converted photon candidates
91  h2_tk_nHitsVsR_ = fs->make<TH2F>("tknHitsVsR", "Tracks Hits vs R ", 12, 0., 120., 20, 0.5, 20.5);
92  h2_tk_inPtVsR_ = fs->make<TH2F>("tkInPtvsR", "Tracks inner Pt vs R ", 12, 0., 120., 100, 0., 100.);
93 
94  return;
95 }
96 
97 float PhotonsWithConversionsAnalyzer::etaTransformation(float EtaParticle, float Zvertex) {
98  //---Definitions
99  const float PI = 3.1415927;
100  //UNUSED const float TWOPI = 2.0*PI;
101 
102  //---Definitions for ECAL
103  const float R_ECAL = 136.5;
104  const float Z_Endcap = 328.0;
105  const float etaBarrelEndcap = 1.479;
106 
107  //---ETA correction
108 
109  float Theta = 0.0;
110  float ZEcal = R_ECAL * sinh(EtaParticle) + Zvertex;
111 
112  if (ZEcal != 0.0)
113  Theta = atan(R_ECAL / ZEcal);
114  if (Theta < 0.0)
115  Theta = Theta + PI;
116  float ETA = -log(tan(0.5 * Theta));
117 
118  if (fabs(ETA) > etaBarrelEndcap) {
119  float Zend = Z_Endcap;
120  if (EtaParticle < 0.0)
121  Zend = -Zend;
122  float Zlen = Zend - Zvertex;
123  float RR = Zlen / sinh(EtaParticle);
124  Theta = atan(RR / Zend);
125  if (Theta < 0.0)
126  Theta = Theta + PI;
127  ETA = -log(tan(0.5 * Theta));
128  }
129  //---Return the result
130  return ETA;
131  //---end
132 }
133 
134 void PhotonsWithConversionsAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
135  using namespace edm;
136  const float etaPhiDistance = 0.01;
137  // Fiducial region
138  //UNUSED const float TRK_BARL =0.9;
139  //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
140  //UNUSED const float END_LO = 1.566;
141  //UNUSED const float END_HI = 2.5;
142  // Electron mass
143  //UNUSED const Float_t mElec= 0.000511;
144 
145  nEvt_++;
146  LogInfo("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id()
147  << " Global Counter " << nEvt_ << "\n";
148  // LogDebug("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
149  std::cout << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ << "\n";
150 
152  Handle<reco::PhotonCollection> photonHandle;
153  e.getByLabel(photonCollectionProducer_, photonCollection_, photonHandle);
154  const reco::PhotonCollection photonCollection = *(photonHandle.product());
155  std::cout << "ConvertedPhotonAnalyzer Photons with conversions collection size " << photonCollection.size() << "\n";
156 
158  std::cout << " ConvertedPhotonAnalyzer Looking for MC truth "
159  << "\n";
160 
161  //get simtrack info
162  std::vector<SimTrack> theSimTracks;
163  std::vector<SimVertex> theSimVertices;
164 
165  edm::Handle<SimTrackContainer> SimTk;
166  edm::Handle<SimVertexContainer> SimVtx;
167  e.getByLabel("g4SimHits", SimTk);
168  e.getByLabel("g4SimHits", SimVtx);
169 
170  theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
171  theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
172 
173  std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
174  std::cout << " ConvertedPhotonAnalyzer mcPhotons size " << mcPhotons.size() << std::endl;
175 
176  // Loop over simulated photons
177  //UNUSED int iDet=0;
178  //UNUSED int iRadius=-1;
179  //UNUSED int indPho=0;
180 
181  for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
182  float mcPhi = (*mcPho).fourMomentum().phi();
183  float mcEta = (*mcPho).fourMomentum().pseudoRapidity();
184  mcEta = etaTransformation(mcEta, (*mcPho).primaryVertex().z());
185 
186  if ((*mcPho).fourMomentum().et() < 20)
187  continue;
188  // if ( ! ( fabs(mcEta) <= BARL || ( fabs(mcEta) >= END_LO && fabs(mcEta) <=END_HI ) ) ) {
189  // continue;
190  //} // all ecal fiducial region
191 
192  h_MCphoE_->Fill((*mcPho).fourMomentum().e());
193  h_MCphoEta_->Fill((*mcPho).fourMomentum().eta());
194  h_MCphoPhi_->Fill((*mcPho).fourMomentum().phi());
195 
196  // keep only visible conversions
197  if ((*mcPho).isAConversion() == 0)
198  continue;
199 
200  nMCPho_++;
201 
202  h_MCConvEta_->Fill(fabs((*mcPho).fourMomentum().pseudoRapidity()) - 0.001);
203 
204  bool REJECTED;
205 
207  //std::cout << " ConvertedPhotonAnalyzer Starting loop over photon candidates " << "\n";
208  for (reco::PhotonCollection::const_iterator iPho = photonCollection.begin(); iPho != photonCollection.end();
209  iPho++) {
210  REJECTED = false;
211 
212  // std::cout << " ConvertedPhotonAnalyzer Reco SC energy " << (*iPho).superCluster()->energy() << "\n";
213 
214  float phiClu = (*iPho).superCluster()->phi();
215  float etaClu = (*iPho).superCluster()->eta();
216  float deltaPhi = phiClu - mcPhi;
217  float deltaEta = etaClu - mcEta;
218 
219  if (deltaPhi > Geom::pi())
220  deltaPhi -= Geom::twoPi();
221  if (deltaPhi < -Geom::pi())
222  deltaPhi += Geom::twoPi();
223  deltaPhi = pow(deltaPhi, 2);
224  deltaEta = pow(deltaEta, 2);
225  float delta = deltaPhi + deltaEta;
226  if (delta >= etaPhiDistance)
227  REJECTED = true;
228 
229  // if ( ! ( fabs(etaClu) <= BARL || ( fabs(etaClu) >= END_LO && fabs(etaClu) <=END_HI ) ) ) REJECTED=true;
230 
231  if (REJECTED)
232  continue;
233  std::cout << " MATCHED " << std::endl;
234  nMatched_++;
235 
236  // std::cout << " ConvertedPhotonAnalyzer Matching candidate " << std::endl;
237 
238  // 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;
239 
240  h_ErecoEMC_->Fill((*iPho).superCluster()->energy() / (*mcPho).fourMomentum().e());
241  h_deltaPhi_->Fill((*iPho).superCluster()->position().phi() - mcPhi);
242  h_deltaEta_->Fill((*iPho).superCluster()->position().eta() - mcEta);
243 
244  h_scE_->Fill((*iPho).superCluster()->energy());
245  h_scEt_->Fill((*iPho).superCluster()->energy() / cosh((*iPho).superCluster()->position().eta()));
246  h_scEta_->Fill((*iPho).superCluster()->position().eta());
247  h_scPhi_->Fill((*iPho).superCluster()->position().phi());
248 
249  h_phoE_->Fill((*iPho).energy());
250  h_phoEta_->Fill((*iPho).eta());
251  h_phoPhi_->Fill((*iPho).phi());
252 
253  if (!(*iPho).hasConversionTracks())
254  continue;
255  // std::cout << " This photons has " << (*iPho).conversions().size() << " conversions candidates " << std::endl;
256  reco::ConversionRefVector conversions = (*iPho).conversions();
257  //std::vector<reco::ConversionRef> conversions = (*iPho).conversions();
258 
259  for (unsigned int i = 0; i < conversions.size(); i++) {
260  //std::cout << " Conversion candidate Energy " << (*iPho).energy() << " number of tracks " << conversions[i]->nTracks() << std::endl;
261  std::vector<edm::RefToBase<reco::Track> > tracks = conversions[i]->tracks();
262 
263  for (unsigned int i = 0; i < tracks.size(); i++) {
264  // std::cout << " ConvertedPhotonAnalyzer Reco Track charge " << tracks[i]->charge() << " Num of RecHits " << tracks[i]->recHitsSize() << " inner momentum " << sqrt ( tracks[i]->innerMomentum().Mag2() ) << "\n";
265 
266  h2_tk_nHitsVsR_->Fill((*mcPho).vertex().perp(), tracks[i]->recHitsSize());
267  h2_tk_inPtVsR_->Fill((*mcPho).vertex().perp(), sqrt(tracks[i]->innerMomentum().Mag2()));
268  }
269  }
270 
271  }
272 
273  }
274 }
275 
276 void PhotonsWithConversionsAnalyzer::endJob() {
277  // fOutputFile_->Write() ;
278  // fOutputFile_->Close() ;
279 
280  edm::LogInfo("ConvertedPhotonAnalyzer") << "Analyzed " << nEvt_ << "\n";
281  // std::cout << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
282  std::cout << "ConvertedPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events "
283  << "\n";
284 
285  return;
286 }
PI
Definition: PayloadInspector.h:19
TrackExtra.h
Handle.h
PDWG_EXOHSCP_cff.tracks
tracks
Definition: PDWG_EXOHSCP_cff.py:28
mps_fire.i
i
Definition: mps_fire.py:355
MessageLogger.h
edm::Handle::product
T const * product() const
Definition: Handle.h:70
ESHandle.h
PI
#define PI
Definition: QcdUeDQM.h:37
PhotonsWithConversionsAnalyzer::endJob
void endJob() override
Definition: PhotonsWithConversionsAnalyzer.cc:276
edm
HLT enums.
Definition: AlignableModifier.h:19
gather_cfg.cout
cout
Definition: gather_cfg.py:144
edm::LogInfo
Definition: MessageLogger.h:254
PhotonFwd.h
edm::RefVector< ConversionCollection >
ConversionFwd.h
edm::Handle< reco::PhotonCollection >
MakerMacros.h
Photon.h
Track.h
TrackFwd.h
spr::deltaEta
static const double deltaEta
Definition: CaloConstants.h:8
ElectronMCTruth.h
PhotonsWithConversionsAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: PhotonsWithConversionsAnalyzer.cc:134
SiPixelRawToDigiRegional_cfi.deltaPhi
deltaPhi
Definition: SiPixelRawToDigiRegional_cfi.py:9
PhotonMCTruth.h
Service.h
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
SimVertex.h
IdealMagneticFieldRecord.h
Geom::pi
constexpr double pi()
Definition: Pi.h:31
Geom::twoPi
constexpr double twoPi()
Definition: Pi.h:32
PhotonMCTruthFinder
Definition: PhotonMCTruthFinder.h:20
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TFileService.h
PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer
PhotonsWithConversionsAnalyzer(const edm::ParameterSet &)
Definition: PhotonsWithConversionsAnalyzer.cc:48
PhotonMCTruthFinder.h
edm::ParameterSet
Definition: ParameterSet.h:36
Event.h
PhotonsWithConversionsAnalyzer::etaTransformation
float etaTransformation(float a, float b)
Definition: PhotonsWithConversionsAnalyzer.cc:97
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
edm::Service< TFileService >
MagneticField.h
edm::EventSetup
Definition: EventSetup.h:57
PhotonsWithConversionsAnalyzer::beginJob
void beginJob() override
Definition: PhotonsWithConversionsAnalyzer.cc:56
std
Definition: JetResolutionObject.h:76
PhotonsWithConversionsAnalyzer::~PhotonsWithConversionsAnalyzer
~PhotonsWithConversionsAnalyzer() override
Definition: PhotonsWithConversionsAnalyzer.cc:54
SuperCluster.h
ETA
#define ETA
Definition: GenericBenchmark.cc:28
EventSetup.h
Exception.h
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
reco::PhotonCollection
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:30
SimTrack.h
ExoticaDQM_cfi.photonCollection
photonCollection
Definition: ExoticaDQM_cfi.py:17
HepMCProduct.h
ZEcal
static constexpr float ZEcal
Definition: L1TkEmParticleProducer.cc:40
pwdgSkimBPark_cfi.conversions
conversions
Definition: pwdgSkimBPark_cfi.py:10
edm::Event
Definition: Event.h:73
PhotonsWithConversionsAnalyzer.h
SimTrackContainer.h
Z_Endcap
static constexpr float Z_Endcap
Definition: ECALPositionCalculator.cc:11
SimVertexContainer.h
R_ECAL
static constexpr float R_ECAL
Definition: ECALPositionCalculator.cc:10
TFileService::make
T * make(const Args &... args) const
make new ROOT object
Definition: TFileService.h:64
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
etaBarrelEndcap
static constexpr float etaBarrelEndcap
Definition: ECALPositionCalculator.cc:12
Conversion.h
Generated for CMSSW Reference Manual by doxygen 1.8.16