CMS 3D CMS Logo
MCPhotonAnalyzer.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/MCPhotonAnalyzer.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 "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
18 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
19 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
20 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
21 //
22 #include "DataFormats/Common/interface/Handle.h"
23 //
24 #include "FWCore/Framework/interface/Event.h"
25 #include "FWCore/Framework/interface/EventSetup.h"
26 #include "FWCore/Framework/interface/ESHandle.h"
27 #include "FWCore/Framework/interface/MakerMacros.h"
28 #include "FWCore/MessageLogger/interface/MessageLogger.h"
29 #include "FWCore/Utilities/interface/Exception.h"
30 //
31 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
32 #include "MagneticField/Engine/interface/MagneticField.h"
33 
34 //
35 #include "TFile.h"
36 #include "TH1.h"
37 #include "TH2.h"
38 #include "TTree.h"
39 #include "TVector3.h"
40 #include "TProfile.h"
41 //
42 
43 using namespace std;
44 
45 MCPhotonAnalyzer::MCPhotonAnalyzer(const edm::ParameterSet& pset) {}
46 
47 MCPhotonAnalyzer::~MCPhotonAnalyzer() { delete thePhotonMCTruthFinder_; }
48 
49 void MCPhotonAnalyzer::beginJob() {
50  nEvt_ = 0;
51 
52  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
53 
54  edm::Service<TFileService> fs;
55 
57  h_MCPhoE_ = fs->make<TH1F>("MCPhoE", "MC photon energy", 100, 0., 100.);
58  h_MCPhoPhi_ = fs->make<TH1F>("MCPhoPhi", "MC photon phi", 40, -3.14, 3.14);
59  h_MCPhoEta_ = fs->make<TH1F>("MCPhoEta", "MC photon eta", 25, 0., 2.5);
60  h_MCPhoEta1_ = fs->make<TH1F>("MCPhoEta1", "MC photon eta", 40, -3., 3.);
61  h_MCPhoEta2_ = fs->make<TH1F>("MCPhoEta2", "MC photon eta", 40, -3., 3.);
62  h_MCPhoEta3_ = fs->make<TH1F>("MCPhoEta3", "MC photon eta", 40, -3., 3.);
63  h_MCPhoEta4_ = fs->make<TH1F>("MCPhoEta4", "MC photon eta", 40, -3., 3.);
65  h_MCConvPhoE_ = fs->make<TH1F>("MCConvPhoE", "MC converted photon energy", 100, 0., 100.);
66  h_MCConvPhoPhi_ = fs->make<TH1F>("MCConvPhoPhi", "MC converted photon phi", 40, -3.14, 3.14);
67  h_MCConvPhoEta_ = fs->make<TH1F>("MCConvPhoEta", "MC converted photon eta", 25, 0., 2.5);
68  h_MCConvPhoR_ = fs->make<TH1F>("MCConvPhoR", "MC converted photon R", 120, 0., 120.);
69 
70  h_MCConvPhoREta1_ = fs->make<TH1F>("MCConvPhoREta1", "MC converted photon R", 120, 0., 120.);
71  h_MCConvPhoREta2_ = fs->make<TH1F>("MCConvPhoREta2", "MC converted photon R", 120, 0., 120.);
72  h_MCConvPhoREta3_ = fs->make<TH1F>("MCConvPhoREta3", "MC converted photon R", 120, 0., 120.);
73  h_MCConvPhoREta4_ = fs->make<TH1F>("MCConvPhoREta4", "MC converted photon R", 120, 0., 120.);
74 
75  h_convFracEta1_ = fs->make<TH1F>("convFracEta1", "Integrated(R) fraction of conversion |eta|=0.2", 120, 0., 120.);
76  h_convFracEta2_ = fs->make<TH1F>("convFracEta2", "Integrated(R) fraction of conversion |eta|=0.9", 120, 0., 120.);
77  h_convFracEta3_ = fs->make<TH1F>("convFracEta3", "Integrated(R) fraction of conversion |eta|=1.5", 120, 0., 120.);
78  h_convFracEta4_ = fs->make<TH1F>("convFracEta4", "Integrated(R) fraction of conversion |eta|=2.0", 120, 0., 120.);
80  h_MCConvPhoTwoTracksE_ =
81  fs->make<TH1F>("MCConvPhoTwoTracksE", "MC converted photon with 2 tracks energy", 100, 0., 100.);
82  h_MCConvPhoTwoTracksPhi_ =
83  fs->make<TH1F>("MCConvPhoTwoTracksPhi", "MC converted photon 2 tracks phi", 40, -3.14, 3.14);
84  h_MCConvPhoTwoTracksEta_ = fs->make<TH1F>("MCConvPhoTwoTracksEta", "MC converted photon 2 tracks eta", 40, -3., 3.);
85  h_MCConvPhoTwoTracksR_ = fs->make<TH1F>("MCConvPhoTwoTracksR", "MC converted photon 2 tracks eta", 48, 0., 120.);
86  // conversions with one track
87  h_MCConvPhoOneTrackE_ =
88  fs->make<TH1F>("MCConvPhoOneTrackE", "MC converted photon with 1 track energy", 100, 0., 100.);
89  h_MCConvPhoOneTrackPhi_ = fs->make<TH1F>("MCConvPhoOneTrackPhi", "MC converted photon 1 track phi", 40, -3.14, 3.14);
90  h_MCConvPhoOneTrackEta_ = fs->make<TH1F>("MCConvPhoOneTrackEta", "MC converted photon 1 track eta", 40, -3., 3.);
91  h_MCConvPhoOneTrackR_ = fs->make<TH1F>("MCConvPhoOneTrackR", "MC converted photon 1 track eta", 48, 0., 120.);
92 
94  h_MCEleE_ = fs->make<TH1F>("MCEleE", "MC ele energy", 100, 0., 200.);
95  h_MCElePhi_ = fs->make<TH1F>("MCElePhi", "MC ele phi", 40, -3.14, 3.14);
96  h_MCEleEta_ = fs->make<TH1F>("MCEleEta", "MC ele eta", 40, -3., 3.);
97  h_BremFrac_ = fs->make<TH1F>("bremFrac", "brem frac ", 100, 0., 1.);
98  h_BremEnergy_ = fs->make<TH1F>("BremE", "Brem energy", 100, 0., 200.);
99  h_EleEvsPhoE_ = fs->make<TH2F>("eleEvsPhoE", "eleEvsPhoE", 100, 0., 200., 100, 0., 200.);
100  h_bremEvsEleE_ = fs->make<TH2F>("bremEvsEleE", "bremEvsEleE", 100, 0., 200., 100, 0., 200.);
101 
102  p_BremVsR_ = fs->make<TProfile>("BremVsR", " Mean Brem Energy vs R ", 48, 0., 120.);
103  p_BremVsEta_ = fs->make<TProfile>("BremVsEta", " Mean Brem Energy vs Eta ", 50, -2.5, 2.5);
104 
105  p_BremVsConvR_ = fs->make<TProfile>("BremVsConvR", " Mean Brem Fraction vs conversion R ", 48, 0., 120.);
106  p_BremVsConvEta_ = fs->make<TProfile>("BremVsConvEta", " Mean Brem Fraction vs converion Eta ", 50, -2.5, 2.5);
107 
108  h_bremFracVsConvR_ = fs->make<TH2F>("bremFracVsConvR", "brem Fraction vs conversion R", 60, 0., 120., 100, 0., 1.);
109 
110  return;
111 }
112 
113 float MCPhotonAnalyzer::etaTransformation(float EtaParticle, float Zvertex) {
114  //---Definitions
115  const float PI = 3.1415927;
116  // const float TWOPI = 2.0*PI;
117 
118  //---Definitions for ECAL
119  const float R_ECAL = 136.5;
120  const float Z_Endcap = 328.0;
121  const float etaBarrelEndcap = 1.479;
122 
123  //---ETA correction
124 
125  float Theta = 0.0;
126  float ZEcal = R_ECAL * sinh(EtaParticle) + Zvertex;
127 
128  if (ZEcal != 0.0)
129  Theta = atan(R_ECAL / ZEcal);
130  if (Theta < 0.0)
131  Theta = Theta + PI;
132  float ETA = -log(tan(0.5 * Theta));
133 
134  if (fabs(ETA) > etaBarrelEndcap) {
135  float Zend = Z_Endcap;
136  if (EtaParticle < 0.0)
137  Zend = -Zend;
138  float Zlen = Zend - Zvertex;
139  float RR = Zlen / sinh(EtaParticle);
140  Theta = atan(RR / Zend);
141  if (Theta < 0.0)
142  Theta = Theta + PI;
143  ETA = -log(tan(0.5 * Theta));
144  }
145  //---Return the result
146  return ETA;
147  //---end
148 }
149 
150 float MCPhotonAnalyzer::phiNormalization(float& phi) {
151  //---Definitions
152  const float PI = 3.1415927;
153  const float TWOPI = 2.0 * PI;
154 
155  if (phi > PI) {
156  phi = phi - TWOPI;
157  }
158  if (phi < -PI) {
159  phi = phi + TWOPI;
160  }
161 
162  // cout << " Float_t PHInormalization out " << PHI << endl;
163  return phi;
164 }
165 
166 void MCPhotonAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
167  using namespace edm;
168  //UNUSED const float etaPhiDistance=0.01;
169  // Fiducial region
170  //UNUSED const float TRK_BARL =0.9;
171  //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
172  //UNUSED const float END_LO = 1.566;
173  //UNUSED const float END_HI = 2.5;
174  // Electron mass
175  //UNUSED const Float_t mElec= 0.000511;
176 
177  nEvt_++;
178  LogInfo("mcEleAnalyzer") << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_
179  << "\n";
180  // LogDebug("MCPhotonAnalyzer") << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
181  std::cout << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ << "\n";
182 
184  std::cout << " MCPhotonAnalyzer Looking for MC truth "
185  << "\n";
186 
187  //get simtrack info
188  std::vector<SimTrack> theSimTracks;
189  std::vector<SimVertex> theSimVertices;
190 
191  edm::Handle<SimTrackContainer> SimTk;
192  edm::Handle<SimVertexContainer> SimVtx;
193  e.getByLabel("g4SimHits", SimTk);
194  e.getByLabel("g4SimHits", SimVtx);
195 
196  theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
197  theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
198  std::cout << " MCPhotonAnalyzer This Event has " << theSimTracks.size() << " sim tracks " << std::endl;
199  std::cout << " MCPhotonAnalyzer This Event has " << theSimVertices.size() << " sim vertices " << std::endl;
200  if (theSimTracks.empty())
201  std::cout << " Event number " << e.id() << " has NO sim tracks " << std::endl;
202 
203  std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
204  std::cout << " MCPhotonAnalyzer mcPhotons size " << mcPhotons.size() << std::endl;
205 
206  for (std::vector<PhotonMCTruth>::const_iterator iPho = mcPhotons.begin(); iPho != mcPhotons.end(); ++iPho) {
207  if ((*iPho).fourMomentum().e() < 35)
208  continue;
209 
210  h_MCPhoE_->Fill((*iPho).fourMomentum().e());
211  // float correta = etaTransformation( (*iPho).fourMomentum().pseudoRapidity(), (*iPho).primaryVertex().z() );
212  float Theta = (*iPho).fourMomentum().theta();
213  float correta = -log(tan(0.5 * Theta));
214  correta = etaTransformation(correta, (*iPho).primaryVertex().z());
215  //h_MCPhoEta_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
216  h_MCPhoEta_->Fill(fabs(correta) - 0.001);
217  h_MCPhoPhi_->Fill((*iPho).fourMomentum().phi());
218 
219  /*
220  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15 )
221  h_MCPhoEta1_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
222  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85 )
223  h_MCPhoEta2_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
224  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55 )
225  h_MCPhoEta3_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
226  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95 )
227  h_MCPhoEta4_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
228  */
229 
230  if (fabs(correta) <= 0.3 && fabs(correta) > 0.2)
231  h_MCPhoEta1_->Fill(correta);
232  if (fabs(correta) <= 1.00 && fabs(correta) > 0.9)
233  h_MCPhoEta2_->Fill(correta);
234  if (fabs(correta) <= 1.6 && fabs(correta) > 1.5)
235  h_MCPhoEta3_->Fill(correta);
236  if (fabs(correta) <= 2. && fabs(correta) > 1.9)
237  h_MCPhoEta4_->Fill(correta);
238 
239  // if ( (*iPho).isAConversion() && (*iPho).vertex().perp()< 10 ) {
240  if ((*iPho).isAConversion()) {
241  h_MCConvPhoE_->Fill((*iPho).fourMomentum().e());
242  // h_MCConvPhoEta_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
243 
244  h_MCConvPhoEta_->Fill(fabs(correta) - 0.001);
245  h_MCConvPhoPhi_->Fill((*iPho).fourMomentum().phi());
246  h_MCConvPhoR_->Fill((*iPho).vertex().perp());
247 
248  /*
249  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15 )
250  h_MCConvPhoREta1_->Fill ( (*iPho).vertex().perp() );
251  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85 )
252  h_MCConvPhoREta2_->Fill ( (*iPho).vertex().perp() );
253  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55 )
254  h_MCConvPhoREta3_->Fill ( (*iPho).vertex().perp() );
255  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95 )
256  h_MCConvPhoREta4_->Fill ( (*iPho).vertex().perp() );
257  */
258 
259  if (fabs(correta) <= 0.3 && fabs(correta) > 0.2)
260  h_MCConvPhoREta1_->Fill((*iPho).vertex().perp());
261  if (fabs(correta) <= 1. && fabs(correta) > 0.9)
262  h_MCConvPhoREta2_->Fill((*iPho).vertex().perp());
263  if (fabs(correta) <= 1.6 && fabs(correta) > 1.5)
264  h_MCConvPhoREta3_->Fill((*iPho).vertex().perp());
265  if (fabs(correta) <= 2 && fabs(correta) > 1.9)
266  h_MCConvPhoREta4_->Fill((*iPho).vertex().perp());
267 
268  } // end conversions
269 
270  }
271 }
272 
273 void MCPhotonAnalyzer::endJob() {
274  Double_t s1 = 0;
275  Double_t s2 = 0;
276  Double_t s3 = 0;
277  Double_t s4 = 0;
278  int e1 = 0;
279  int e2 = 0;
280  int e3 = 0;
281  int e4 = 0;
282 
283  Double_t nTotEta1 = h_MCPhoEta1_->GetEntries();
284  Double_t nTotEta2 = h_MCPhoEta2_->GetEntries();
285  Double_t nTotEta3 = h_MCPhoEta3_->GetEntries();
286  Double_t nTotEta4 = h_MCPhoEta4_->GetEntries();
287 
288  for (int i = 1; i <= 120; ++i) {
289  e1 = (int)h_MCConvPhoREta1_->GetBinContent(i);
290  e2 = (int)h_MCConvPhoREta2_->GetBinContent(i);
291  e3 = (int)h_MCConvPhoREta3_->GetBinContent(i);
292  e4 = (int)h_MCConvPhoREta4_->GetBinContent(i);
293  s1 += e1;
294  s2 += e2;
295  s3 += e3;
296  s4 += e4;
297  h_convFracEta1_->SetBinContent(i, s1 * 100 / nTotEta1);
298  h_convFracEta2_->SetBinContent(i, s2 * 100 / nTotEta2);
299  h_convFracEta3_->SetBinContent(i, s3 * 100 / nTotEta3);
300  h_convFracEta4_->SetBinContent(i, s4 * 100 / nTotEta4);
301  }
302 
303  edm::LogInfo("MCPhotonAnalyzer") << "Analyzed " << nEvt_ << "\n";
304  std::cout << "MCPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events "
305  << "\n";
306 
307  return;
308 }
MCPhotonAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: MCPhotonAnalyzer.cc:166
PI
Definition: PayloadInspector.h:20
Handle.h
MCPhotonAnalyzer::phiNormalization
float phiNormalization(float &a)
Definition: MCPhotonAnalyzer.cc:150
mps_fire.i
i
Definition: mps_fire.py:428
MessageLogger.h
ESHandle.h
PI
#define PI
Definition: QcdUeDQM.h:37
edm
HLT enums.
Definition: AlignableModifier.h:19
CrossingFrame.h
gather_cfg.cout
cout
Definition: gather_cfg.py:144
MCPhotonAnalyzer::MCPhotonAnalyzer
MCPhotonAnalyzer(const edm::ParameterSet &)
Definition: MCPhotonAnalyzer.cc:45
indexGen.s2
s2
Definition: indexGen.py:107
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
edm::Handle
Definition: AssociativeIterator.h:50
MakerMacros.h
TrackingVertexContainer.h
ElectronMCTruth.h
MixCollection.h
PhotonMCTruth.h
Service.h
SimVertex.h
IdealMagneticFieldRecord.h
MCPhotonAnalyzer::endJob
void endJob() override
Definition: MCPhotonAnalyzer.cc:273
MCPhotonAnalyzer::~MCPhotonAnalyzer
~MCPhotonAnalyzer() override
Definition: MCPhotonAnalyzer.cc:47
PhotonMCTruthFinder
Definition: PhotonMCTruthFinder.h:20
TFileService.h
PhotonMCTruthFinder.h
edm::ParameterSet
Definition: ParameterSet.h:47
Event.h
TWOPI
#define TWOPI
Definition: DQMSourcePi0.cc:36
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
edm::Service< TFileService >
createfilelist.int
int
Definition: createfilelist.py:10
MagneticField.h
edm::EventSetup
Definition: EventSetup.h:57
SiPixelPhase1Clusters_cfi.e3
e3
Definition: SiPixelPhase1Clusters_cfi.py:9
MCPhotonAnalyzer::etaTransformation
float etaTransformation(float a, float b)
Definition: MCPhotonAnalyzer.cc:113
std
Definition: JetResolutionObject.h:76
MCPhotonAnalyzer::beginJob
void beginJob() override
Definition: MCPhotonAnalyzer.cc:49
TrackingParticleFwd.h
ETA
#define ETA
Definition: GenericBenchmark.cc:28
MCPhotonAnalyzer.h
EventSetup.h
Exception.h
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
SimTrack.h
HepMCProduct.h
ZEcal
static constexpr float ZEcal
Definition: L1TkEmParticleProducer.cc:40
edm::Event
Definition: Event.h:73
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
vertexPlots.e4
e4
Definition: vertexPlots.py:64
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
Generated for CMSSW Reference Manual by doxygen 1.8.16