CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
MCPhotonAnalyzer.cc
Go to the documentation of this file.
20 
21 #include "TFile.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "TTree.h"
25 #include "TVector3.h"
26 #include "TProfile.h"
27 
28 #include <iostream>
29 #include <map>
30 #include <vector>
31 
33 public:
34  //
35  explicit MCPhotonAnalyzer(const edm::ParameterSet&);
36  ~MCPhotonAnalyzer() override;
37 
38  void analyze(const edm::Event&, const edm::EventSetup&) override;
39  void beginJob() override;
40  void endJob() override;
41 
42 private:
43  float etaTransformation(float a, float b);
44  float phiNormalization(float& a);
45 
46  //
48 
50 
52  TFile* fOutputFile_;
53 
54  int nEvt_;
55  int nMatched_;
56 
58 
59  double mcPhi_;
60  double mcEta_;
61 
66 
67  // all photons
68  TH1F* h_MCPhoE_;
69  TH1F* h_MCPhoEta_;
70  TH1F* h_MCPhoEta1_;
71  TH1F* h_MCPhoEta2_;
72  TH1F* h_MCPhoEta3_;
73  TH1F* h_MCPhoEta4_;
74  TH1F* h_MCPhoPhi_;
75  // Conversion
88 
99 
100  TH1F* h_MCEleE_;
101  TH1F* h_MCEleEta_;
102  TH1F* h_MCElePhi_;
103  TH1F* h_BremFrac_;
107 
108  TProfile* p_BremVsR_;
109  TProfile* p_BremVsEta_;
110 
111  TProfile* p_BremVsConvR_;
112  TProfile* p_BremVsConvEta_;
113 
115 };
116 
119 
120 using namespace std;
121 
123 
124 MCPhotonAnalyzer::~MCPhotonAnalyzer() { delete thePhotonMCTruthFinder_; }
125 
127  nEvt_ = 0;
128 
129  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
130 
132 
134  h_MCPhoE_ = fs->make<TH1F>("MCPhoE", "MC photon energy", 100, 0., 100.);
135  h_MCPhoPhi_ = fs->make<TH1F>("MCPhoPhi", "MC photon phi", 40, -3.14, 3.14);
136  h_MCPhoEta_ = fs->make<TH1F>("MCPhoEta", "MC photon eta", 25, 0., 2.5);
137  h_MCPhoEta1_ = fs->make<TH1F>("MCPhoEta1", "MC photon eta", 40, -3., 3.);
138  h_MCPhoEta2_ = fs->make<TH1F>("MCPhoEta2", "MC photon eta", 40, -3., 3.);
139  h_MCPhoEta3_ = fs->make<TH1F>("MCPhoEta3", "MC photon eta", 40, -3., 3.);
140  h_MCPhoEta4_ = fs->make<TH1F>("MCPhoEta4", "MC photon eta", 40, -3., 3.);
142  h_MCConvPhoE_ = fs->make<TH1F>("MCConvPhoE", "MC converted photon energy", 100, 0., 100.);
143  h_MCConvPhoPhi_ = fs->make<TH1F>("MCConvPhoPhi", "MC converted photon phi", 40, -3.14, 3.14);
144  h_MCConvPhoEta_ = fs->make<TH1F>("MCConvPhoEta", "MC converted photon eta", 25, 0., 2.5);
145  h_MCConvPhoR_ = fs->make<TH1F>("MCConvPhoR", "MC converted photon R", 120, 0., 120.);
146 
147  h_MCConvPhoREta1_ = fs->make<TH1F>("MCConvPhoREta1", "MC converted photon R", 120, 0., 120.);
148  h_MCConvPhoREta2_ = fs->make<TH1F>("MCConvPhoREta2", "MC converted photon R", 120, 0., 120.);
149  h_MCConvPhoREta3_ = fs->make<TH1F>("MCConvPhoREta3", "MC converted photon R", 120, 0., 120.);
150  h_MCConvPhoREta4_ = fs->make<TH1F>("MCConvPhoREta4", "MC converted photon R", 120, 0., 120.);
151 
152  h_convFracEta1_ = fs->make<TH1F>("convFracEta1", "Integrated(R) fraction of conversion |eta|=0.2", 120, 0., 120.);
153  h_convFracEta2_ = fs->make<TH1F>("convFracEta2", "Integrated(R) fraction of conversion |eta|=0.9", 120, 0., 120.);
154  h_convFracEta3_ = fs->make<TH1F>("convFracEta3", "Integrated(R) fraction of conversion |eta|=1.5", 120, 0., 120.);
155  h_convFracEta4_ = fs->make<TH1F>("convFracEta4", "Integrated(R) fraction of conversion |eta|=2.0", 120, 0., 120.);
157  h_MCConvPhoTwoTracksE_ =
158  fs->make<TH1F>("MCConvPhoTwoTracksE", "MC converted photon with 2 tracks energy", 100, 0., 100.);
159  h_MCConvPhoTwoTracksPhi_ =
160  fs->make<TH1F>("MCConvPhoTwoTracksPhi", "MC converted photon 2 tracks phi", 40, -3.14, 3.14);
161  h_MCConvPhoTwoTracksEta_ = fs->make<TH1F>("MCConvPhoTwoTracksEta", "MC converted photon 2 tracks eta", 40, -3., 3.);
162  h_MCConvPhoTwoTracksR_ = fs->make<TH1F>("MCConvPhoTwoTracksR", "MC converted photon 2 tracks eta", 48, 0., 120.);
163  // conversions with one track
164  h_MCConvPhoOneTrackE_ =
165  fs->make<TH1F>("MCConvPhoOneTrackE", "MC converted photon with 1 track energy", 100, 0., 100.);
166  h_MCConvPhoOneTrackPhi_ = fs->make<TH1F>("MCConvPhoOneTrackPhi", "MC converted photon 1 track phi", 40, -3.14, 3.14);
167  h_MCConvPhoOneTrackEta_ = fs->make<TH1F>("MCConvPhoOneTrackEta", "MC converted photon 1 track eta", 40, -3., 3.);
168  h_MCConvPhoOneTrackR_ = fs->make<TH1F>("MCConvPhoOneTrackR", "MC converted photon 1 track eta", 48, 0., 120.);
169 
171  h_MCEleE_ = fs->make<TH1F>("MCEleE", "MC ele energy", 100, 0., 200.);
172  h_MCElePhi_ = fs->make<TH1F>("MCElePhi", "MC ele phi", 40, -3.14, 3.14);
173  h_MCEleEta_ = fs->make<TH1F>("MCEleEta", "MC ele eta", 40, -3., 3.);
174  h_BremFrac_ = fs->make<TH1F>("bremFrac", "brem frac ", 100, 0., 1.);
175  h_BremEnergy_ = fs->make<TH1F>("BremE", "Brem energy", 100, 0., 200.);
176  h_EleEvsPhoE_ = fs->make<TH2F>("eleEvsPhoE", "eleEvsPhoE", 100, 0., 200., 100, 0., 200.);
177  h_bremEvsEleE_ = fs->make<TH2F>("bremEvsEleE", "bremEvsEleE", 100, 0., 200., 100, 0., 200.);
178 
179  p_BremVsR_ = fs->make<TProfile>("BremVsR", " Mean Brem Energy vs R ", 48, 0., 120.);
180  p_BremVsEta_ = fs->make<TProfile>("BremVsEta", " Mean Brem Energy vs Eta ", 50, -2.5, 2.5);
181 
182  p_BremVsConvR_ = fs->make<TProfile>("BremVsConvR", " Mean Brem Fraction vs conversion R ", 48, 0., 120.);
183  p_BremVsConvEta_ = fs->make<TProfile>("BremVsConvEta", " Mean Brem Fraction vs converion Eta ", 50, -2.5, 2.5);
184 
185  h_bremFracVsConvR_ = fs->make<TH2F>("bremFracVsConvR", "brem Fraction vs conversion R", 60, 0., 120., 100, 0., 1.);
186 
187  return;
188 }
189 
190 float MCPhotonAnalyzer::etaTransformation(float EtaParticle, float Zvertex) {
191  //---Definitions
192  const float PI = 3.1415927;
193  // const float TWOPI = 2.0*PI;
194 
195  //---Definitions for ECAL
196  const float R_ECAL = 136.5;
197  const float Z_Endcap = 328.0;
198  const float etaBarrelEndcap = 1.479;
199 
200  //---ETA correction
201 
202  float Theta = 0.0;
203  float ZEcal = R_ECAL * sinh(EtaParticle) + Zvertex;
204 
205  if (ZEcal != 0.0)
206  Theta = atan(R_ECAL / ZEcal);
207  if (Theta < 0.0)
208  Theta = Theta + PI;
209  float ETA = -log(tan(0.5 * Theta));
210 
211  if (fabs(ETA) > etaBarrelEndcap) {
212  float Zend = Z_Endcap;
213  if (EtaParticle < 0.0)
214  Zend = -Zend;
215  float Zlen = Zend - Zvertex;
216  float RR = Zlen / sinh(EtaParticle);
217  Theta = atan(RR / Zend);
218  if (Theta < 0.0)
219  Theta = Theta + PI;
220  ETA = -log(tan(0.5 * Theta));
221  }
222  //---Return the result
223  return ETA;
224  //---end
225 }
226 
228  //---Definitions
229  const float PI = 3.1415927;
230  const float TWOPI = 2.0 * PI;
231 
232  if (phi > PI) {
233  phi = phi - TWOPI;
234  }
235  if (phi < -PI) {
236  phi = phi + TWOPI;
237  }
238 
239  // cout << " Float_t PHInormalization out " << PHI << endl;
240  return phi;
241 }
242 
244  using namespace edm;
245  //UNUSED const float etaPhiDistance=0.01;
246  // Fiducial region
247  //UNUSED const float TRK_BARL =0.9;
248  //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
249  //UNUSED const float END_LO = 1.566;
250  //UNUSED const float END_HI = 2.5;
251  // Electron mass
252  //UNUSED const Float_t mElec= 0.000511;
253 
254  nEvt_++;
255  LogInfo("mcEleAnalyzer") << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_
256  << "\n";
257  // LogDebug("MCPhotonAnalyzer") << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
258  std::cout << "MCPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ << "\n";
259 
261  std::cout << " MCPhotonAnalyzer Looking for MC truth "
262  << "\n";
263 
264  //get simtrack info
265  std::vector<SimTrack> theSimTracks;
266  std::vector<SimVertex> theSimVertices;
267 
270  e.getByLabel("g4SimHits", SimTk);
271  e.getByLabel("g4SimHits", SimVtx);
272 
273  theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
274  theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
275  std::cout << " MCPhotonAnalyzer This Event has " << theSimTracks.size() << " sim tracks " << std::endl;
276  std::cout << " MCPhotonAnalyzer This Event has " << theSimVertices.size() << " sim vertices " << std::endl;
277  if (theSimTracks.empty())
278  std::cout << " Event number " << e.id() << " has NO sim tracks " << std::endl;
279 
280  std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
281  std::cout << " MCPhotonAnalyzer mcPhotons size " << mcPhotons.size() << std::endl;
282 
283  for (std::vector<PhotonMCTruth>::const_iterator iPho = mcPhotons.begin(); iPho != mcPhotons.end(); ++iPho) {
284  if ((*iPho).fourMomentum().e() < 35)
285  continue;
286 
287  h_MCPhoE_->Fill((*iPho).fourMomentum().e());
288  // float correta = etaTransformation( (*iPho).fourMomentum().pseudoRapidity(), (*iPho).primaryVertex().z() );
289  float Theta = (*iPho).fourMomentum().theta();
290  float correta = -log(tan(0.5 * Theta));
291  correta = etaTransformation(correta, (*iPho).primaryVertex().z());
292  //h_MCPhoEta_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
293  h_MCPhoEta_->Fill(fabs(correta) - 0.001);
294  h_MCPhoPhi_->Fill((*iPho).fourMomentum().phi());
295 
296  /*
297  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15 )
298  h_MCPhoEta1_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
299  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85 )
300  h_MCPhoEta2_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
301  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55 )
302  h_MCPhoEta3_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
303  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95 )
304  h_MCPhoEta4_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
305  */
306 
307  if (fabs(correta) <= 0.3 && fabs(correta) > 0.2)
308  h_MCPhoEta1_->Fill(correta);
309  if (fabs(correta) <= 1.00 && fabs(correta) > 0.9)
310  h_MCPhoEta2_->Fill(correta);
311  if (fabs(correta) <= 1.6 && fabs(correta) > 1.5)
312  h_MCPhoEta3_->Fill(correta);
313  if (fabs(correta) <= 2. && fabs(correta) > 1.9)
314  h_MCPhoEta4_->Fill(correta);
315 
316  // if ( (*iPho).isAConversion() && (*iPho).vertex().perp()< 10 ) {
317  if ((*iPho).isAConversion()) {
318  h_MCConvPhoE_->Fill((*iPho).fourMomentum().e());
319  // h_MCConvPhoEta_->Fill ( (*iPho).fourMomentum().pseudoRapidity() );
320 
321  h_MCConvPhoEta_->Fill(fabs(correta) - 0.001);
322  h_MCConvPhoPhi_->Fill((*iPho).fourMomentum().phi());
323  h_MCConvPhoR_->Fill((*iPho).vertex().perp());
324 
325  /*
326  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.25 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.15 )
327  h_MCConvPhoREta1_->Fill ( (*iPho).vertex().perp() );
328  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 0.95 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=0.85 )
329  h_MCConvPhoREta2_->Fill ( (*iPho).vertex().perp() );
330  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 1.65 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.55 )
331  h_MCConvPhoREta3_->Fill ( (*iPho).vertex().perp() );
332  if ( fabs((*iPho).fourMomentum().pseudoRapidity() ) <= 2.05 && fabs((*iPho).fourMomentum().pseudoRapidity() ) >=1.95 )
333  h_MCConvPhoREta4_->Fill ( (*iPho).vertex().perp() );
334  */
335 
336  if (fabs(correta) <= 0.3 && fabs(correta) > 0.2)
337  h_MCConvPhoREta1_->Fill((*iPho).vertex().perp());
338  if (fabs(correta) <= 1. && fabs(correta) > 0.9)
339  h_MCConvPhoREta2_->Fill((*iPho).vertex().perp());
340  if (fabs(correta) <= 1.6 && fabs(correta) > 1.5)
341  h_MCConvPhoREta3_->Fill((*iPho).vertex().perp());
342  if (fabs(correta) <= 2 && fabs(correta) > 1.9)
343  h_MCConvPhoREta4_->Fill((*iPho).vertex().perp());
344 
345  } // end conversions
346 
347  }
348 }
349 
351  double s1 = 0;
352  double s2 = 0;
353  double s3 = 0;
354  double s4 = 0;
355  int e1 = 0;
356  int e2 = 0;
357  int e3 = 0;
358  int e4 = 0;
359 
360  double nTotEta1 = h_MCPhoEta1_->GetEntries();
361  double nTotEta2 = h_MCPhoEta2_->GetEntries();
362  double nTotEta3 = h_MCPhoEta3_->GetEntries();
363  double nTotEta4 = h_MCPhoEta4_->GetEntries();
364 
365  for (int i = 1; i <= 120; ++i) {
366  e1 = (int)h_MCConvPhoREta1_->GetBinContent(i);
367  e2 = (int)h_MCConvPhoREta2_->GetBinContent(i);
368  e3 = (int)h_MCConvPhoREta3_->GetBinContent(i);
369  e4 = (int)h_MCConvPhoREta4_->GetBinContent(i);
370  s1 += e1;
371  s2 += e2;
372  s3 += e3;
373  s4 += e4;
374  h_convFracEta1_->SetBinContent(i, s1 * 100 / nTotEta1);
375  h_convFracEta2_->SetBinContent(i, s2 * 100 / nTotEta2);
376  h_convFracEta3_->SetBinContent(i, s3 * 100 / nTotEta3);
377  h_convFracEta4_->SetBinContent(i, s4 * 100 / nTotEta4);
378  }
379 
380  edm::LogInfo("MCPhotonAnalyzer") << "Analyzed " << nEvt_ << "\n";
381  std::cout << "MCPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events "
382  << "\n";
383 
384  return;
385 }
float etaTransformation(float a, float b)
TH1F * h_MCConvPhoOneTrackE_
Conversions with one track.
double mcPhi_
global variable for the MC photon
static std::vector< std::string > checklist log
~MCPhotonAnalyzer() override
TH1F * h_MCConvPhoTwoTracksPhi_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const TrackerGeometry * trackerGeom
void endJob() override
TH1F * h_MCConvPhoOneTrackPhi_
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
static constexpr float R_ECAL
TProfile * p_BremVsConvR_
std::string fOutputFileName_
#define ETA
TH1F * h_MCConvPhoTwoTracksE_
Conversions with two tracks.
TH1F * h_MCConvPhoOneTrackEta_
TProfile * p_BremVsConvEta_
std::string SimTkLabel
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
#define TWOPI
static constexpr float etaBarrelEndcap
MCPhotonAnalyzer(const edm::ParameterSet &)
#define PI
Definition: QcdUeDQM.h:37
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
void beginJob() override
std::string HepMCLabel
std::string SimVtxLabel
Log< level::Info, false > LogInfo
void analyze(const edm::Event &, const edm::EventSetup &) override
double b
Definition: hdecay.h:118
std::string SimHitLabel
TH1F * h_MCConvPhoTwoTracksEta_
edm::EventID id() const
Definition: EventBase.h:59
double a
Definition: hdecay.h:119
static constexpr float ZEcal
PhotonMCTruthFinder * thePhotonMCTruthFinder_
tuple cout
Definition: gather_cfg.py:144
float phiNormalization(float &a)
static constexpr float Z_Endcap