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.);
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.);
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.);
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.);
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.);
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);
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);
185 h_bremFracVsConvR_ = fs->
make<TH2F>(
"bremFracVsConvR",
"brem Fraction vs conversion R", 60, 0., 120., 100, 0., 1.);
192 const float PI = 3.1415927;
196 const float R_ECAL = 136.5;
203 float ZEcal =
R_ECAL * sinh(EtaParticle) + Zvertex;
213 if (EtaParticle < 0.0)
215 float Zlen = Zend - Zvertex;
216 float RR = Zlen / sinh(EtaParticle);
217 Theta = atan(RR / Zend);
229 const float PI = 3.1415927;
255 LogInfo(
"mcEleAnalyzer") <<
"MCPhotonAnalyzer Analyzing event number: " <<
e.id() <<
" Global Counter " << nEvt_
258 std::cout <<
"MCPhotonAnalyzer Analyzing event number: " <<
e.id() <<
" Global Counter " << nEvt_ <<
"\n";
261 std::cout <<
" MCPhotonAnalyzer Looking for MC truth "
265 std::vector<SimTrack> theSimTracks;
266 std::vector<SimVertex> theSimVertices;
270 e.getByLabel(
"g4SimHits", SimTk);
271 e.getByLabel(
"g4SimHits", SimVtx);
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;
280 std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
281 std::cout <<
" MCPhotonAnalyzer mcPhotons size " << mcPhotons.size() << std::endl;
283 for (std::vector<PhotonMCTruth>::const_iterator iPho = mcPhotons.begin(); iPho != mcPhotons.end(); ++iPho) {
284 if ((*iPho).fourMomentum().e() < 35)
287 h_MCPhoE_->Fill((*iPho).fourMomentum().e());
289 float Theta = (*iPho).fourMomentum().theta();
290 float correta = -
log(
tan(0.5 * Theta));
291 correta = etaTransformation(correta, (*iPho).primaryVertex().z());
293 h_MCPhoEta_->Fill(fabs(correta) - 0.001);
294 h_MCPhoPhi_->Fill((*iPho).fourMomentum().phi());
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);
317 if ((*iPho).isAConversion()) {
318 h_MCConvPhoE_->Fill((*iPho).fourMomentum().e());
321 h_MCConvPhoEta_->Fill(fabs(correta) - 0.001);
322 h_MCConvPhoPhi_->Fill((*iPho).fourMomentum().phi());
323 h_MCConvPhoR_->Fill((*iPho).vertex().perp());
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());
360 double nTotEta1 = h_MCPhoEta1_->GetEntries();
361 double nTotEta2 = h_MCPhoEta2_->GetEntries();
362 double nTotEta3 = h_MCPhoEta3_->GetEntries();
363 double nTotEta4 = h_MCPhoEta4_->GetEntries();
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);
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);
380 edm::LogInfo(
"MCPhotonAnalyzer") <<
"Analyzed " << nEvt_ <<
"\n";
381 std::cout <<
"MCPhotonAnalyzer::endJob Analyzed " << nEvt_ <<
" events "