ATLAS Offline Software
Loading...
Searching...
No Matches
GetLCSinglePionsPerf.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5//-----------------------------------------------------------------------
6// File and Version Information:
7// $Id: CheckSinglePionsReco.cxx,v 1.1.1.1 2008/11/04 08:56:11 menke Exp $
8//
9// Description: see CheckSinglePionsReco.h
10//
11// Environment:
12// Software developed for the ATLAS Detector at CERN LHC
13//
14// Author List:
15// Gennady Pospelov
16//
17//-----------------------------------------------------------------------
18
19//-----------------------
20// This Class's Header --
21//-----------------------
23
24#include <iterator>
25#include <cmath>
26
28#include "GaudiKernel/ISvcLocator.h"
29#include "GaudiKernel/StatusCode.h"
30
31#include "CaloEvent/CaloCell.h"
37#include <CLHEP/Vector/LorentzVector.h>
40
43
44#include "TProfile.h"
45#include "TProfile2D.h"
46#include "TFile.h"
47#include "TH1F.h"
48#include "TVector3.h"
49#include "TTree.h"
50
51
52using CLHEP::HepLorentzVector;
53using CLHEP::MeV;
54
55
56/* ****************************************************************************
57
58**************************************************************************** */
59GetLCSinglePionsPerf::GetLCSinglePionsPerf(const std::string& name, ISvcLocator* pSvcLocator)
60 : AthAlgorithm(name, pSvcLocator),
61 m_id_helper(nullptr),
62 m_calo_id(nullptr),
63 m_caloDmDescrManager(nullptr),
64 m_outputFileName("CheckSinglePionsReco.root"),
65 m_outputFile(nullptr),
66 m_distance_cut(1.5),
67 m_netabin(25), m_etamin(0.0), m_etamax(5.0),
68 m_deta(0),
70 m_dphi(0),
73 m_dlogener(0),
74 m_nnormtype(3),
75 m_ncluscoll(4),
76 m_ntagcases(3),
77 m_nmoments(0),
78 m_nmomsums(6),
81 m_doEngTag(true),
82 m_doEngRecSpect(true),
83 m_doEngNoiseClus(true),
84 m_doClusMoments(true),
87 m_useGoodClus(false),
89 m_mc_ener(0),
90 m_mc_eta(0),
91 m_mc_phi(0),
92 m_mc_etabin(0),
93 m_mc_enerbin(0),
95{
96 // collection name where cluster always have moments
97 declareProperty("ClusterBasicCollName",m_clusterBasicCollName);
98
99 // collections names after different correction steps
100 declareProperty("ClusterCollectionNames",m_clusterCollNames);
101
102 // distance cut for clusters to calcuate engReco
103 declareProperty("DistanceCut",m_distance_cut);
104
105 // Name of output file to save histograms in
106 declareProperty("OutputFileName",m_outputFileName);
107
108 // Names of CalibrationHitContainers to use
109 declareProperty("CalibrationHitContainerNames",m_CalibrationHitContainerNames);
110 declareProperty("DMCalibrationHitContainerNames",m_DMCalibrationHitContainerNames);
111
112 // to create different types of validation histograms
113 declareProperty("doEngRecOverTruth",m_doEngRecOverTruth);
114 declareProperty("doEngTag",m_doEngTag);
115 declareProperty("doEngRecSpect",m_doEngRecSpect);
116 declareProperty("doEngNoiseClus",m_doEngNoiseClus);
117 declareProperty("doClusMoments",m_doClusMoments);
118 declareProperty("doRecoEfficiency",m_doRecoEfficiency);
119 declareProperty("useRecoEfficiency",m_useRecoEfficiency);
120 declareProperty("useGoodClus",m_useGoodClus);
121
122 declareProperty("etamin",m_etamin);
123 declareProperty("etamax",m_etamax);
124 declareProperty("netabin",m_netabin);
125
126 declareProperty("phimin",m_phimin);
127 declareProperty("phimax",m_phimax);
128 declareProperty("nphibin",m_nphibin);
129
130 declareProperty("logenermin",m_logenermin);
131 declareProperty("logenermax",m_logenermax);
132 declareProperty("nlogenerbin",m_nlogenerbin);
133
134 declareProperty("doCalibHitsValidation",m_doCalibHitsValidation);
135// Only use clusters that contain at least m_truthPiEngFraction % of the true pion energy
136 declareProperty("usePionClustersOnly",m_usePionClusters);
137
138 // dead material identifier description manager
140
141 m_validMoments.resize(0);
142 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_TOT"),xAOD::CaloCluster::ENG_CALIB_TOT));
143 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_OUT_L"),xAOD::CaloCluster::ENG_CALIB_OUT_L));
144 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TOT"),xAOD::CaloCluster::ENG_CALIB_DEAD_TOT));
145 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_EMB0"),xAOD::CaloCluster::ENG_CALIB_DEAD_EMB0));
146 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TILE0"),xAOD::CaloCluster::ENG_CALIB_DEAD_TILE0));
147 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TILEG3"),xAOD::CaloCluster::ENG_CALIB_DEAD_TILEG3));
148 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_EME0"),xAOD::CaloCluster::ENG_CALIB_DEAD_EME0));
149 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_HEC0"),xAOD::CaloCluster::ENG_CALIB_DEAD_HEC0));
150 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_FCAL"),xAOD::CaloCluster::ENG_CALIB_DEAD_FCAL));
151 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_LEAKAGE"),xAOD::CaloCluster::ENG_CALIB_DEAD_LEAKAGE));
152 m_validMoments.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_UNCLASS"),xAOD::CaloCluster::ENG_CALIB_DEAD_UNCLASS));
153
154 m_nmoments = m_validMoments.size();
155
156}
157
158
159
160/* ****************************************************************************
161
162**************************************************************************** */
167
168
169
170/* ****************************************************************************
171
172**************************************************************************** */
174{
175 ATH_CHECK( detStore()->retrieve(m_calo_id,"CaloCell_ID") );
176 ATH_CHECK( detStore()->retrieve(m_id_helper) );
177
178
182
183 m_outputFile = new TFile(m_outputFileName.c_str(),"RECREATE");
184 m_outputFile->cd();
185
186 // defining histograms
187 char hname[256];
188 float h_eta_min = -5.0;
189 float h_eta_max=5.0;
190 int h_nch_eta = 50;
191 if(m_etamin != 0.0) {
192 h_eta_min = m_etamin;
193 h_eta_max = m_etamax;
194 }
195 const int n_logener_bins = m_nlogenerbin*2;
196 double *xbins = new double[n_logener_bins+1];
197 const double inv_n_logener_bins = 1. / static_cast<double> (n_logener_bins);
198 for(int i_bin=0; i_bin<=n_logener_bins; i_bin++) {
199 xbins[i_bin] = pow(10,m_logenermin+i_bin*(m_logenermax-m_logenermin)*inv_n_logener_bins);
200 }
202
203 /* **************************************************************************
204 histograms to check classification performance
205 *************************************************************************** */
206 if(m_doEngTag) {
207 // as a function of pion eta for several energy bins
209 for(int i_tag=0; i_tag<m_ntagcases; i_tag++) {
210 m_engTag_vs_eta[i_tag].resize(m_nlogenerbin, nullptr);
211 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
212 sprintf(hname, "hp_engTag_vs_eta_tag%d_ener%d",i_tag, i_ener);
213 TProfile *hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
214 hp->GetXaxis()->SetTitle("#eta");
215 hp->GetYaxis()->SetTitle("E_{class}/E_{tot}");
216 m_engTag_vs_eta[i_tag][i_ener] = hp;
217 } // i_ener
218 } // i_tag
219 // as a function of pion energy for several eta bins
221 for(int i_tag=0; i_tag<m_ntagcases; i_tag++){
222 m_engTag_vs_ebeam[i_tag].resize(m_netabin, nullptr);
223 for(int i_eta=0; i_eta<m_netabin; i_eta++) {
224 sprintf(hname, "hp_engTag_vs_ebeam_tag%d_eta%d",i_tag, i_eta);
225 TProfile *hp = new TProfile(hname, hname, n_logener_bins, xbins);
226 hp->GetXaxis()->SetTitle("E_{#pi}, MeV");
227 hp->GetYaxis()->SetTitle("E_{class}/E_{tot}");
228 m_engTag_vs_ebeam[i_tag][i_eta] = hp;
229 }
230 }
231 } // m_doEngTag
232
233
234 /* **************************************************************************
235 histograms to check reconstructed energy after each step of local hadronic
236 calibration, reco energy is defined as a sum of clusters energies within certain
237 angle around particle direction
238 *************************************************************************** */
240 // as a function of pion eta for several energy bins
242 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
244 for(int i_norm=0; i_norm<m_nnormtype; i_norm++){
245 m_engRecOverTruth_vs_eta[i_coll][i_norm].resize(m_nlogenerbin, nullptr);
246 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
247 double ylow(0), yup(0);
248 if(i_norm != kLEVEL_PARTICLE) {
249 ylow = 0.0;
250 yup = 100.0;
251 }
252 sprintf(hname, "hp_engRecOverTruth_vs_eta_coll%d_norm%d_ener%d",i_coll, i_norm, i_ener);
253 TProfile *hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max, ylow, yup);
254 hp->GetXaxis()->SetTitle("#eta");
255 if(i_norm==kLEVEL_PARTICLE) {
256 hp->GetYaxis()->SetTitle("E_{reco}/E_{#pi}");
257 }else{
258 hp->GetYaxis()->SetTitle("E_{reco}/E_{calib}");
259 }
260 m_engRecOverTruth_vs_eta[i_coll][i_norm][i_ener] = hp;
261 } // i_ener
262 } // i_norm
263 } // i_coll
264 // as a function of pion energy for several eta bins
266 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
268 for(int i_norm=0; i_norm<m_nnormtype; i_norm++){
269 m_engRecOverTruth_vs_ebeam[i_coll][i_norm].resize(m_netabin, nullptr);
270 for(int i_eta=0; i_eta<m_netabin; i_eta++) {
271 double ylow(0), yup(0);
272 if(i_norm != kLEVEL_PARTICLE) {
273 ylow = 0.0;
274 yup = 100.0;
275 }
276 sprintf(hname, "hp_m_engRecOverTruth_vs_ebeam_coll%d_norm%d_eta%d",i_coll, i_norm, i_eta);
277 TProfile *hp = new TProfile(hname, hname, n_logener_bins, xbins, ylow, yup);
278 hp->GetXaxis()->SetTitle("E_{#pi}, MeV");
279 if(i_norm==kLEVEL_PARTICLE) {
280 hp->GetYaxis()->SetTitle("E_{reco}/E_{#pi}");
281 }else{
282 hp->GetYaxis()->SetTitle("E_{reco}/E_{calib}");
283 }
284 m_engRecOverTruth_vs_ebeam[i_coll][i_norm][i_eta] = hp;
285 } // i_eta
286 } // i_norm
287 } // i_coll
288 } // m_doEngRecOverTruth
289
290
291 /* **************************************************************************
292 histograms to save reco and noise clusters energy spectras
293 *************************************************************************** */
294 if(m_doEngRecSpect) {
296 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
297 m_engRecSpect[i_coll].resize(m_nlogenerbin);
298 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
299 m_engRecSpect[i_coll][i_ener].resize(m_netabin, nullptr);
300 for(int i_eta=0; i_eta<m_netabin; i_eta++){
301 sprintf(hname, "hp_engRecSpect_coll%d_ener%d_eta%d",i_coll, i_ener, i_eta);
302 TH1F *h1 = new TH1F(hname, hname, 110, -0.2, 2.0);
303 h1->GetXaxis()->SetTitle("E_{reco}/E_{#pi}");
304 m_engRecSpect[i_coll][i_ener][i_eta] = h1;
305 } // i_eta
306 } // i_ener
307 } // i_coll
308 // energy spectra of incident pions
310 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
311 m_engPionSpect[i_ener].resize(m_netabin, nullptr);
312 for(int i_eta=0; i_eta<m_netabin; i_eta++){
313 sprintf(hname, "hp_engPionSpect_ener%d_eta%d",i_ener, i_eta);
314 double x1 = pow(10,m_logenermin+i_ener*m_dlogener);
315 double x2 = pow(10,m_logenermin+(i_ener+1)*m_dlogener);
316 TH1F *h1 = new TH1F(hname, hname, 20, x1, x2);
317 h1->GetXaxis()->SetTitle("E_{#pi}, MeV");
318 m_engPionSpect[i_ener][i_eta] = h1;
319 } // i_eta
320 } // i_ener
321 } // m_doEngRecSpect
322
323 /* **************************************************************************
324 histograms to save noise clusters spectras
325 *************************************************************************** */
326 if(m_doEngNoiseClus) {
327 // noise clusters spectras
329 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
330 m_engNoiseClusSpect[i_coll].resize(m_netabin, nullptr);
331 for(int i_eta=0; i_eta<m_netabin; i_eta++){
332 sprintf(hname, "hp_engNoiseClusSpect_coll%d_eta%d",i_coll, i_eta);
333 TH1F *h1 = new TH1F(hname, hname, 100, -2000., 2000.);
334 h1->GetXaxis()->SetTitle("E_{clus}, MeV");
335 m_engNoiseClusSpect[i_coll][i_eta] = h1;
336 } // i_eta
337 } // i_coll
338 // average
339 m_engNoiseClus_vs_eta.resize(m_ncluscoll, nullptr);
340 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
341 sprintf(hname, "hp_engNoiseClus_vs_eta_coll%d",i_coll);
342 TProfile *hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
343 hp->GetXaxis()->SetTitle("#eta");
344 hp->GetYaxis()->SetTitle("E_{clus}, MeV");
345 m_engNoiseClus_vs_eta[i_coll] = hp;
346 //hp->SetErrorOption("s");
347 }
348 } // m_doEngNoiseClus
349
350 /* **************************************************************************
351 histograms to save cluster moments
352 *************************************************************************** */
353 if(m_doClusMoments) {
355 for(int i_mom=0; i_mom<m_nmoments; i_mom++){
356 m_clusMoment_vs_eta[i_mom].resize(m_nmomsums);
357 for(int i_sum=0; i_sum<m_nmomsums; i_sum++){
358 m_clusMoment_vs_eta[i_mom][i_sum].resize(m_nlogenerbin, nullptr);
359 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
360 sprintf(hname, "hp_clusMoment_vs_eta_mom%d_sum%d_ener%d",i_mom, i_sum, i_ener);
361 TProfile *hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
362 hp->GetXaxis()->SetTitle("#eta");
363 hp->GetYaxis()->SetTitle("Mom_{cls}/Mom_{sum}");
364 hp->SetTitle(m_validMoments[i_mom].first.c_str());
365 m_clusMoment_vs_eta[i_mom][i_sum][i_ener] = hp;
366 } // i_ener
367 } // i_sum
368 } // i_mom
369
371 for(int i_mom=0; i_mom<m_nmoments; i_mom++){
372 m_clusMoment_vs_ebeam[i_mom].resize(m_nmomsums);
373 for(int i_sum=0; i_sum<m_nmomsums; i_sum++){
374 m_clusMoment_vs_ebeam[i_mom][i_sum].resize(m_netabin, nullptr);
375 for(int i_eta=0; i_eta<m_netabin; i_eta++){
376 sprintf(hname, "hp_clusMoment_vs_ebeam_mom%d_sum%d_eta%d",i_mom, i_sum, i_eta);
377 TProfile *hp = new TProfile(hname, hname, n_logener_bins, xbins);
378 hp->GetXaxis()->SetTitle("E_{#pi}");
379 hp->GetYaxis()->SetTitle("Mom_{cls}/Mom_{sum}");
380 hp->SetTitle(m_validMoments[i_mom].first.c_str());
381 m_clusMoment_vs_ebeam[i_mom][i_sum][i_eta] = hp;
382 } // i_eta
383 } // i_sum
384 } // i_mom
385
386 } // m_doClusMoments
387
388
389 /* **************************************************************************
390 histograms to save reconstruction efficiency
391 *************************************************************************** */
393 // as a function of eta for several energy bins
394 for(int i_eff=0; i_eff<2; i_eff++){ // 0-reconstructed, 1-all events
395 m_RecoEfficiency_vs_eta[i_eff].resize(m_nlogenerbin, nullptr);
396 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
397 sprintf(hname, "hp_m_RecoEfficiency_eff%d_ener%d", i_eff, i_ener);
398 TH1F *h1 = new TH1F(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
399 h1->GetXaxis()->SetTitle("#eta");
400 h1->GetYaxis()->SetTitle("Efficiency");
401 m_RecoEfficiency_vs_eta[i_eff][i_ener] = h1;
402 } // i_eta
403 } // i_eff
404 // as a function of energy for several eta bins
405 for(int i_eff=0; i_eff<2; i_eff++){ // 0-reconstructed, 1-all events
406 m_RecoEfficiency_vs_ebeam[i_eff].resize(m_netabin, nullptr);
407 for(int i_eta=0; i_eta<m_netabin; i_eta++){
408 sprintf(hname, "hp_m_RecoEfficiency_eff%d_eta%d", i_eff, i_eta);
409 TH1F *h1 = new TH1F(hname, hname, n_logener_bins, xbins);
410 h1->GetXaxis()->SetTitle("E_{#pi}, MeV");
411 h1->GetYaxis()->SetTitle("Efficiency");
412 m_RecoEfficiency_vs_ebeam[i_eff][i_eta] = h1;
413 } // i_eta
414 } // i_eff
415 }
416
417 /* **************************************************************************
418 histograms to check calibration hits machinery
419 *************************************************************************** */
425 for(int i_ener=0; i_ener<m_nlogenerbin2; i_ener++){
426 sprintf(hname, "hp_SumCalibHitOverEbeam_vs_eta_ener%d", i_ener);
427 TProfile *hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
428 hp->GetXaxis()->SetTitle("#eta");
429 hp->GetYaxis()->SetTitle("#sum E_{calib}/E_{#pi}");
431 sprintf(hname, "hp_DefaultCalculatorOverEbeam_vs_eta_ener%d", i_ener);
432 hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
433 hp->GetXaxis()->SetTitle("#eta");
434 hp->GetYaxis()->SetTitle("#sum E_{DefCalc}/E_{#pi}");
436 sprintf(hname, "hp_SumCalibHitAssignedOverEbeam_vs_eta_ener%d", i_ener);
437 hp = new TProfile(hname, hname, h_nch_eta, h_eta_min, h_eta_max);
438 hp->GetXaxis()->SetTitle("#eta");
439 hp->GetYaxis()->SetTitle("#sum E_{assigned}/E_{#pi}");
441 sprintf(hname, "hp_SumCalibHitAssignedOverEbeam_vs_etaphi_ener%d", i_ener);
442 TProfile2D *hp2 = new TProfile2D(hname, hname, 25, h_eta_min, h_eta_max, 25, -M_PI, M_PI);
443 hp2->GetXaxis()->SetTitle("#eta");
444 hp2->GetYaxis()->SetTitle("#phi");
445 hp2->GetZaxis()->SetTitle("#sum E_{assigned}/E_{#pi}");
447 } // i_ener
448
449// m_SumCalibHitOverEbeam_vs_ebeam.resize(m_netabin, 0);
450// // m_DefaultCalculatorOverEbeam_vs_ebeam.resize(m_netabin, 0);
451// for(int i_eta=0; i_eta<m_netabin; i_eta++){
452// sprintf(hname, "hp_SumCalibHitOverEbeam_vs_ebeam_eta%d", i_eta);
453// TProfile *hp = new TProfile(hname, hname, n_logener_bins, xbins);
454// hp->GetXaxis()->SetTitle("#eta");
455// hp->GetYaxis()->SetTitle("#sum E_{calib}/E_{#pi}");
456// m_SumCalibHitOverEbeam_vs_ebeam[i_eta] = hp;
457// // sprintf(hname, "hp_DefaultCalculatorOverEbeam_vs_ebeam_eta%d", i_eta);
458// // hp = new TProfile(hname, hname, n_logener_bins, xbins);
459// // hp->GetXaxis()->SetTitle("#eta");
460// // hp->GetYaxis()->SetTitle("#sum E_{DefCalc}/E_{#pi}");
461// // m_DefaultCalculatorOverEbeam_vs_ebeam[i_eta] = hp;
462// } // i_eta
463 } // m_doCalibHitsValidation
464
465
466 delete [] xbins;
467
468 ATH_CHECK( m_clusterBasicCollName.initialize() );
469 ATH_CHECK( m_clusterCollNames.initialize() );
472
473 return StatusCode::SUCCESS;
474}
475
476
477
478/* ****************************************************************************
479writing histograms into the file
480**************************************************************************** */
482{
483 ATH_MSG_INFO( "Writing out histograms" );
484
485 // ereco over truth after each step of calibration
487 m_outputFile->cd();
488 m_outputFile->mkdir("engRecOverTruth");
489 m_outputFile->cd("engRecOverTruth");
490 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
491 for(int i_norm=0; i_norm<m_nnormtype; i_norm++){
492 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
493 if(m_engRecOverTruth_vs_eta[i_coll][i_norm][i_ener]) m_engRecOverTruth_vs_eta[i_coll][i_norm][i_ener]->Write();
494 } // i_ener
495 } // i_norm
496 } // i_coll
497 //
498 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
499 for(int i_norm=0; i_norm<m_nnormtype; i_norm++){
500 for(int i_eta=0; i_eta<m_netabin; i_eta++) {
501 if(m_engRecOverTruth_vs_ebeam[i_coll][i_norm][i_eta]) m_engRecOverTruth_vs_ebeam[i_coll][i_norm][i_eta]->Write();
502 } // i_eta
503 } // i_norm
504 } // i_coll
505 } // m_doEngRecOverTruth
506
507 // classification performance
508 if(m_doEngTag) {
509 m_outputFile->cd();
510 m_outputFile->mkdir("engTag");
511 m_outputFile->cd("engTag");
512 for(int i_tag=0; i_tag<m_ntagcases; i_tag++) {
513 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
514 if(m_engTag_vs_eta[i_tag][i_ener]) m_engTag_vs_eta[i_tag][i_ener]->Write();
515 } // i_ener
516 } // i_tag
517 // as a function of pion energy for several eta bins
518 for(int i_tag=0; i_tag<m_ntagcases; i_tag++){
519 for(int i_eta=0; i_eta<m_netabin; i_eta++) {
520 if(m_engTag_vs_ebeam[i_tag][i_eta]) m_engTag_vs_ebeam[i_tag][i_eta]->Write();
521 }
522 }
523 } // m_doEngTag
524
525 // reco energy spactras
526 if(m_doEngRecSpect) {
527 m_outputFile->cd();
528 m_outputFile->mkdir("engRecSpect");
529 m_outputFile->cd("engRecSpect");
530 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
531 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
532 for(int i_eta=0; i_eta<m_netabin; i_eta++){
533 if(m_engRecSpect[i_coll][i_ener][i_eta]) m_engRecSpect[i_coll][i_ener][i_eta]->Write();
534 } // i_eta
535 } // i_ener
536 } // i_coll
537 // pione energy spectra
538 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
539 for(int i_eta=0; i_eta<m_netabin; i_eta++){
540 if(m_engPionSpect[i_ener][i_eta]) m_engPionSpect[i_ener][i_eta]->Write();
541 } // i_eta
542 } // i_ener
543 } // m_doEngRecSpect
544
545 // noise clusters spectra
546 if(m_doEngNoiseClus) {
547 m_outputFile->cd();
548 m_outputFile->mkdir("engNoiseClus");
549 m_outputFile->cd("engNoiseClus");
550 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
551 for(int i_eta=0; i_eta<m_netabin; i_eta++){
552 if(m_engNoiseClusSpect[i_coll][i_eta]) m_engNoiseClusSpect[i_coll][i_eta]->Write();
553 } // i_eta
554 } // i_coll
555 // average noise
556 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
557 if(m_engNoiseClus_vs_eta[i_coll]) m_engNoiseClus_vs_eta[i_coll] -> Write();
558 }
559 } // m_doEngRecSpect
560
561 if(m_doClusMoments) {
562 m_outputFile->cd();
563 m_outputFile->mkdir("clusMoment");
564 m_outputFile->cd("clusMoment");
565 for(int i_mom=0; i_mom<m_nmoments; i_mom++){
566 for(int i_sum=0; i_sum<m_nmomsums; i_sum++){
567 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
568 if(m_clusMoment_vs_eta[i_mom][i_sum][i_ener]) m_clusMoment_vs_eta[i_mom][i_sum][i_ener]->Write();
569 } // i_ener
570 } // i_sum
571 } // i_mom
572 for(int i_mom=0; i_mom<m_nmoments; i_mom++){
573 for(int i_sum=0; i_sum<m_nmomsums; i_sum++){
574 for(int i_eta=0; i_eta<m_netabin; i_eta++){
575 if(m_clusMoment_vs_ebeam[i_mom][i_sum][i_eta]) m_clusMoment_vs_ebeam[i_mom][i_sum][i_eta]->Write();
576 } // i_eta
577 } // i_sum
578 } // i_mom
579
580 } // m_doClusMoments
581
583 m_outputFile->cd();
584 m_outputFile->mkdir("RecoEff");
585 m_outputFile->cd("RecoEff");
586 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
587 if(m_RecoEfficiency_vs_eta[0][i_ener] && m_RecoEfficiency_vs_eta[1][i_ener]) {
588 m_RecoEfficiency_vs_eta[0][i_ener]->Write();
589 m_RecoEfficiency_vs_eta[1][i_ener]->Write();
590 }
591 } // i_eta
592 for(int i_eta=0; i_eta<m_netabin; i_eta++){
593 if(m_RecoEfficiency_vs_ebeam[0][i_eta] && m_RecoEfficiency_vs_ebeam[1][i_eta]) {
594 m_RecoEfficiency_vs_ebeam[0][i_eta]->Write();
595 m_RecoEfficiency_vs_ebeam[1][i_eta]->Write();
596 }
597 } // i_eta
598 } // m_doRecoEfficiency
599
600 /* **************************************************************************
601 histograms to check calibration hits machinery
602 *************************************************************************** */
604 m_outputFile->cd();
605 m_outputFile->mkdir("CalibHitsValidation");
606 m_outputFile->cd("CalibHitsValidation");
607 for(int i_ener=0; i_ener<m_nlogenerbin2; i_ener++){
612 } // i_ener
613// for(int i_eta=0; i_eta<m_netabin; i_eta++){
614// if(m_SumCalibHitOverEbeam_vs_ebeam[i_eta]) m_SumCalibHitOverEbeam_vs_ebeam[i_eta]->Write();
615// // if(m_DefaultCalculatorOverEbeam_vs_ebeam[i_eta]) m_DefaultCalculatorOverEbeam_vs_ebeam[i_eta]->Write();
616// } // i_eta
617
618 } // m_doCalibHitsValidation
619
620 /* ********************************************
621 writing tree with histogram parameters
622 ******************************************** */
623 m_outputFile->cd();
624 m_outputFile->mkdir("Parameters");
625 m_outputFile->cd("Parameters");
626 TTree *tree = new TTree("ParamTree","ParamTree");
627 tree->Branch("m_netabin",&m_netabin,"m_netabin/I");
628 tree->Branch("m_etamin",&m_etamin,"m_etamin/F");
629 tree->Branch("m_etamax",&m_etamax,"m_etamax/F");
630 tree->Branch("m_nphibin",&m_nphibin,"m_nphibin/I");
631 tree->Branch("m_phimin",&m_phimin,"m_phimin/F");
632 tree->Branch("m_phimax",&m_phimax,"m_phimax/F");
633 tree->Branch("m_nlogenerbin",&m_nlogenerbin,"m_nlogenerbin/I");
634 tree->Branch("m_logenermin",&m_logenermin,"m_logenermin/F");
635 tree->Branch("m_logenermax",&m_logenermax,"m_logenermax/F");
636 tree->Branch("m_nnormtype",&m_nnormtype,"m_nnormtype/I");
637 tree->Branch("m_ncluscoll",&m_ncluscoll,"m_ncluscoll/I");
638 tree->Branch("m_ntagcases",&m_ntagcases,"m_ntagcases/I");
639 tree->Branch("m_nmoments",&m_nmoments,"m_nmoments/I");
640 tree->Branch("m_nmomsums",&m_nmomsums,"m_nmomsums/I");
641 tree->Fill();
642 tree->Write();
643
644 m_outputFile->Close();
645
646 return StatusCode::SUCCESS;
647}
648
649
650
651/* ****************************************************************************
652
653**************************************************************************** */
655{
656 const EventContext& ctx = getContext();
657
658 /* ********************************************
659 reading particles
660 ******************************************** */
661 const McEventCollection* truthEvent=nullptr;
662 ATH_CHECK( evtStore()->retrieve(truthEvent, "TruthEvent") );
663
664 if( truthEvent->at(0)->particles_empty() ){
665 ATH_MSG_ERROR( "No particles in McEventCollection" );
666 return StatusCode::FAILURE;
667 }
668 HepMC::ConstGenParticlePtr gen{nullptr};
669 for (auto p: *truthEvent->at(0)) {
670 if (MC::isGenStable(p)) {
671 if (std::abs(p->pdg_id()) == MC::PIPLUS ||
672 std::abs(p->pdg_id()) == MC::PI0) {
673 gen = std::move(p);
674 break;
675 }
676 }
677 }
678 if(!gen){
679 ATH_MSG_ERROR( "No final stable pion in McEventCollection" );
680 return StatusCode::FAILURE;
681 }
682 m_mc_eta = gen->momentum().pseudoRapidity();
683 m_mc_phi = gen->momentum().phi();
684 m_mc_ener = gen->momentum().e();
685
686 // energy and eta bins for particle
688 m_mc_etabin = int((fabs(m_mc_eta)-m_etamin)/m_deta);
690
692
694 if(clusColl->empty()) {
695 ATH_MSG_WARNING( "Empty ClusterContainer " << m_clusterBasicCollName );
696 return StatusCode::SUCCESS;
697 }
698
699 // check that we have at least one meaningfull cluster
701 bool PionWasReconstructed = false;
702 float engClusCalibThs = 1.0*MeV;
703 HepLorentzVector hlv_pion(1,0,0,1);
704 hlv_pion.setREtaPhi(1./cosh(m_mc_eta),m_mc_eta,m_mc_phi);
705 for (const xAOD::CaloCluster* theCluster : *clusColl) {
706 double mx_calib_tot;
707 if( !theCluster->retrieveMoment( xAOD::CaloCluster::ENG_CALIB_TOT, mx_calib_tot) ) {
708 ATH_MSG_ERROR( "Moment ENG_CALIB_TOT is absent" );
709 }
710 double mx_dens, mx_center_lambda;
711 if (!theCluster->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA, mx_center_lambda))
712 mx_center_lambda = 0;
713 if (!theCluster->retrieveMoment(xAOD::CaloCluster::FIRST_ENG_DENS, mx_dens))
714 mx_dens = 0;
715
716 HepLorentzVector hlv_cls(1,0,0,1);
717 hlv_cls.setREtaPhi(1./cosh(theCluster->eta()), theCluster->eta(), theCluster->phi());
718 double r = hlv_pion.angle(hlv_cls.vect());
720 && mx_calib_tot > engClusCalibThs
721 && theCluster->size() > 1
722 && mx_dens > 0
723 && mx_center_lambda != 0.0 ) {
724 PionWasReconstructed = true;
725 break;
726 }
727 }
728 if(PionWasReconstructed){
731 }
734 if(m_useRecoEfficiency && !PionWasReconstructed) return StatusCode::SUCCESS;
735 }
736
737 if(m_doEngRecSpect) {
739 }
740
742 fill_reco (*clusColl, ctx);
743 }
744
745 if(m_doClusMoments) {
746 fill_moments (*clusColl, ctx);
747 }
748
750 fill_calibhits (*clusColl, ctx);
751 }
752
753 return StatusCode::SUCCESS;
754}
755
756
757
758/* ****************************************************************************
759+ average reconstructed energy after different steps of calibration, energy spectras
760+ classification performance
761+ noise clusters
762**************************************************************************** */
764 const EventContext& ctx)
765{
766 /* ********************************************
767 reading cluster collections for each calibration level
768 ******************************************** */
769 std::vector<const xAOD::CaloClusterContainer *> clusCollVector;
772 clusCollVector.push_back(pColl.cptr());
773 }
774
775 HepLorentzVector hlv_pion(1,0,0,1);
776 hlv_pion.setREtaPhi(1./cosh(m_mc_eta),m_mc_eta,m_mc_phi);
777
778 std::vector<double > engClusSum; // for 4 collections
779 engClusSum.resize(m_ncluscoll, 0.0);
780 std::vector<double > engClusSumPresOnly; // for 4 collections
781 engClusSumPresOnly.resize(m_ncluscoll, 0.0);
782
783 double engClusSumCalib(0);
784 double engClusSumCalibPresOnly(0);
785 double engClusSumTrueOOC(0);
786 double engClusSumTrueDM(0);
787 std::vector<double> engClusSumCalibTagged;
788 engClusSumCalibTagged.resize(m_ntagcases, 0.0);
789
790 int nGoodClus = 0;
791 unsigned int iClus = -1;
792 for (const xAOD::CaloCluster* theCluster : clusColl) {
793 ++iClus;
794 double mx_calib_tot;
795 if( !theCluster->retrieveMoment( xAOD::CaloCluster::ENG_CALIB_TOT, mx_calib_tot) ) {
796 ATH_MSG_ERROR( "Moment ENG_CALIB_TOT is absent" );
797 }
798 double mx_calib_emb0 = 0;
799 double mx_calib_eme0 = 0;
800 double mx_calib_tileg3 = 0;
801 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EMB0, mx_calib_emb0)
802 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EME0, mx_calib_eme0)
803 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TILEG3, mx_calib_tileg3)){
804 ATH_MSG_ERROR( "One of the moment ENG_CALIB_EMB0, ENG_CALIB_EME0, ENG_CALIB_TILEG3 is absent" );
805 }
806 double mx_calib_out = 0;
807 double mx_calib_dead_tot = 0;
808 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_OUT_L, mx_calib_out)
809 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_DEAD_TOT, mx_calib_dead_tot) ){
810 ATH_MSG_ERROR( "One of the moment ENG_CALIB_OUT_L, ENG_CALIB_DEAD_TOT is absent" );
811 }
812
813 double mx_dens, mx_center_lambda;
814 if (!theCluster->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA, mx_center_lambda))
815 mx_center_lambda = 0;
816 if (!theCluster->retrieveMoment(xAOD::CaloCluster::FIRST_ENG_DENS, mx_dens))
817 mx_dens = 0;
818
819 bool clusIsGood(true);//Fabian
820 if(m_useGoodClus){
821 if(theCluster->size() <= 1 || mx_dens <=0 || mx_center_lambda ==0 || ( m_usePionClusters && mx_calib_tot <= m_truthPiEngFraction * m_mc_ener )) {
822 clusIsGood = false;
823 }
824 }
825
826 HepLorentzVector hlv_cls(1,0,0,1);
827 hlv_cls.setREtaPhi(1./cosh(theCluster->eta()), theCluster->eta(), theCluster->phi());
828 double r = hlv_pion.angle(hlv_cls.vect());
829
831 ( ( m_usePionClusters && r < m_distance_cut*angle_mollier_factor(m_mc_eta)) && clusIsGood ) ) {
832
833 theCluster = clusCollVector[kTOPO]->at(iClus); // this collection has RecoStatus
834 CaloRecoStatus recoStatus = theCluster->recoStatus();
835
836 if(recoStatus.checkStatus(CaloRecoStatus::TAGGEDEM) ) {
837 engClusSumCalibTagged[kTAGEM] += mx_calib_tot;
838 }else if(recoStatus.checkStatus(CaloRecoStatus::TAGGEDHAD) ) {
839 engClusSumCalibTagged[kTAGHAD] += mx_calib_tot;
840 }else if(recoStatus.checkStatus(CaloRecoStatus::TAGGEDUNKNOWN) ) {
841 engClusSumCalibTagged[kTAGUNK] += mx_calib_tot;
842 }else{
843 ATH_MSG_ERROR("CheckSinglePionsReco::execute() -> Error! Unkown classification status " << recoStatus.getStatusWord());
844 }
845
846 if(clusIsGood){
847 nGoodClus++;
848 engClusSumCalib += mx_calib_tot;
849 engClusSumCalibPresOnly += (mx_calib_emb0+mx_calib_eme0+mx_calib_tileg3);
850 engClusSumTrueOOC += mx_calib_out;
851 engClusSumTrueDM += mx_calib_dead_tot;
852 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
853 theCluster = clusCollVector[i_coll]->at(iClus);
854 engClusSum[i_coll] += theCluster->e();
855 engClusSumPresOnly[i_coll] += (theCluster->eSample(CaloSampling::PreSamplerB)+theCluster->eSample(CaloSampling::PreSamplerE)+theCluster->eSample(CaloSampling::TileGap3));
856 }
857 }
858
859 }else{
861 // to save noise clusters
862 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
863 theCluster = clusCollVector[i_coll]->at(iClus);
864 m_engNoiseClusSpect[i_coll][m_mc_etabin]->Fill(theCluster->e());
865 m_engNoiseClus_vs_eta[i_coll]->Fill(theCluster->eta(), theCluster->e());
866 }
867 } // m_doEngNoiseClus
868 } // r
869 // xxx
870 } // iclus
871
872 if(nGoodClus == 0) return 0;
873
874 /* ********************************************
875 filling histograms for ereco over truth after each step of calibration
876 ******************************************** */
878 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++) {
879 // with respect to the initial pion energy
880 int i_norm = 0;
881 double enom = engClusSum[i_coll];
882 if(i_coll != kTOPO) {
883 enom = enom - engClusSumPresOnly[i_coll];
884 }
885 if(m_mc_ener != 0.0) {
886 m_engRecOverTruth_vs_eta[i_coll][i_norm][m_mc_enerbin]->Fill(m_mc_eta, enom/m_mc_ener);
887 m_engRecOverTruth_vs_ebeam[i_coll][i_norm][m_mc_etabin]->Fill(m_mc_ener, enom/m_mc_ener);
888 }
889 // with respect to the calibration level
890 i_norm = 1;
891 double edenom = engClusSumCalib - engClusSumCalibPresOnly; // for em or had
892 if(i_coll == kTOPO_OOC){ // ooc
893 edenom += engClusSumTrueOOC;
894 }else if(i_coll == kTOPO){ // dm
895 edenom += (engClusSumTrueOOC + engClusSumTrueDM);
896 }
897 if(edenom != 0.0) {
898 m_engRecOverTruth_vs_eta[i_coll][i_norm][m_mc_enerbin]->Fill(m_mc_eta, enom/edenom);
899 m_engRecOverTruth_vs_ebeam[i_coll][i_norm][m_mc_etabin]->Fill(m_mc_ener, enom/edenom);
900 }
901 // with respect to the calibration level with "ideal" previous step
902 i_norm = 2;
903 enom = engClusSum[i_coll] - engClusSumPresOnly[i_coll];
904 if(i_coll == kTOPO_OOC){ // ooc
905 double engClusSumRecoOOC = engClusSum[kTOPO_OOC] - engClusSum[kTOPO_W];
906 enom = engClusSumCalib - engClusSumCalibPresOnly + engClusSumRecoOOC;
907 }else if(i_coll == kTOPO){ // dm
908 double engClusSumRecoDM = engClusSum[kTOPO] - (engClusSum[kTOPO_OOC]-engClusSumPresOnly[kTOPO_OOC]);
909 enom = (engClusSumCalib - engClusSumCalibPresOnly) + engClusSumTrueOOC + engClusSumRecoDM;
910 }
911 if(edenom != 0.0) {
912 m_engRecOverTruth_vs_eta[i_coll][i_norm][m_mc_enerbin]->Fill(m_mc_eta, enom/edenom);
913 m_engRecOverTruth_vs_ebeam[i_coll][i_norm][m_mc_etabin]->Fill(m_mc_ener, enom/edenom);
914 }
915 } // i_coll
916 } // m_doEngRecOverTruth
917
918 /* ********************************************
919 filling histograms with classification results
920 ******************************************** */
921 if(m_doEngTag) {
922 if(engClusSumCalib!=0.0) {
923 const double inv_engClusSumCalib = 1. / engClusSumCalib;
924 for(int i_tag=0; i_tag<m_ntagcases; i_tag++) {
925 m_engTag_vs_eta[i_tag][m_mc_enerbin]->Fill(m_mc_eta, engClusSumCalibTagged[i_tag]*inv_engClusSumCalib);
926 m_engTag_vs_ebeam[i_tag][m_mc_etabin]->Fill(m_mc_ener, engClusSumCalibTagged[i_tag]*inv_engClusSumCalib);
927 }
928 }
929 } // m_doEngTag
930
931 /* ********************************************
932 filling histograms for reco energy spectras
933 ******************************************** */
934 if(m_doEngRecSpect && m_mc_ener > 0){
935 for(int i_coll=0; i_coll<m_ncluscoll; i_coll++){
936 double enom = engClusSum[i_coll];
937 if(i_coll != kTOPO) {
938 enom = enom - engClusSumPresOnly[i_coll];
939 }
940 m_engRecSpect[i_coll][m_mc_enerbin][m_mc_etabin]->Fill(enom/m_mc_ener);
941 }
942 } // m_doEngRecSpect
943
944 return 0;
945}
946
947
948
949/* ****************************************************************************
950to check moments assignment
951**************************************************************************** */
953 const EventContext& ctx)
954{
955 ATH_MSG_DEBUG("GetLCSinglePionsPerf::fill_moments got cont. size: "<<clusColl.size());
956 /* ********************************************
957 reading calibration containers
958 ******************************************** */
959 std::vector<const CaloCalibrationHitContainer *> v_cchc;
962 v_cchc.push_back(cchc.cptr());
963 }
964
965 std::vector<const CaloCalibrationHitContainer *> v_dmcchc;
968 v_dmcchc.push_back(cchc.cptr());
969 }
970
971 // calculate total calibration energy int active+inactive hits
972 double engCalibTot = 0.0;
973 std::vector<double > engCalibTotSmp;
974 engCalibTotSmp.resize(CaloSampling::Unknown, 0.0);
975 for (const CaloCalibrationHitContainer* cchc : v_cchc) {
976 //loop over cells in calibration container
977 for (const CaloCalibrationHit* hit : *cchc) {
978 Identifier myId = hit->cellID();
979 if(!myId.is_valid()) {
980 ATH_MSG_WARNING("Error! Bad identifier " << myId << " in container '" << cchc->Name() << "',"
981 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
982 continue;
983 }
984 if(!m_id_helper->is_lar(myId) && !m_id_helper->is_tile(myId)) {
985 ATH_MSG_WARNING("Error! Bad identifier (nor lar or tile) " << myId << " in container '" << cchc->Name() << "',"
986 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
987 continue;
988 }
989
991 engCalibTot += hit->energyTotal();
992 engCalibTotSmp[nsmp] += hit->energyTotal();
993 }
994 }
995
996
997 double engCalibDeadTot = 0.0;
998 std::vector<double > engCalibDeadTotInArea;
999 engCalibDeadTotInArea.resize(CaloDmDescrArea::DMA_MAX,0.0);
1000 for (const CaloCalibrationHitContainer* dmcchc : v_dmcchc) {
1001 //loop over cells in calibration container
1002 for (const CaloCalibrationHit* hit : *dmcchc) {
1003 Identifier myId = hit->cellID();
1004 if(!myId.is_valid()) {
1005 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Bad dead material identifier " << myId << " in container '" << dmcchc->Name() << "',"
1006 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1007 continue;
1008 }
1009 if(!m_id_helper->is_lar_dm(myId) && !m_id_helper->is_tile_dm(myId)) {
1010 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Bad dead material identifier (nor lar_dm or tile_dm) " << myId << " in container '" << dmcchc->Name() << "',"
1011 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1012 continue;
1013 }
1014 engCalibDeadTot += hit->energyTotal();
1015 //
1016 CaloDmDescrElement* myCDDE(nullptr);
1017 myCDDE = m_caloDmDescrManager->get_element(myId);
1018 if( !myCDDE ) {
1019 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Absent DM description element!");
1020 continue;
1021 }
1022 int nDmArea = m_caloDmDescrManager->get_dm_area(myId);
1023 engCalibDeadTotInArea[nDmArea] += hit->energyTotal();
1024 engCalibDeadTotInArea[CaloDmDescrArea::DMA_ALL] += hit->energyTotal();
1025 }
1026 }
1027// double engCalibDeadTotWithPres = engCalibDeadTot + engCalibTotSmp[CaloSampling::PreSamplerB]
1028// + engCalibTotSmp[CaloSampling::PreSamplerE]
1029// + engCalibTotSmp[CaloSampling::TileGap3];
1030
1031 /* ********************************************
1032 reading clusters moments
1033 ******************************************** */
1034 double engClusSumCalib(0);
1035 double engClusSumCalibPresOnly(0);
1036 std::vector<std::vector<double > > clsMoments; // [iClus][m_nmoments]
1037 clsMoments.resize(clusColl.size());
1038 std::vector<double > clsMomentsSum; // [m_nmoments]
1039 clsMomentsSum.resize(m_nmoments,0.0);
1040
1041 // retrieving cluster moments
1042 unsigned int iClus = -1;
1043 for (const xAOD::CaloCluster* theCluster : clusColl) {
1044 ++iClus;
1045 clsMoments[iClus].resize(m_validMoments.size(),0.0);
1046 int iMoment=0;
1047 for(moment_name_vector::const_iterator im=m_validMoments.begin(); im!=m_validMoments.end(); ++im, ++iMoment){
1048 double mx;
1049 if( theCluster->retrieveMoment((*im).second, mx) ) {
1050 clsMoments[iClus][iMoment] = mx;
1051 }else{
1052 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Can't retrieve moment " << (*im).first << " " << (*im).second);
1053 }
1054 clsMomentsSum[iMoment] += clsMoments[iClus][iMoment];
1055 }
1056 double mx;
1057 if(theCluster->retrieveMoment( xAOD::CaloCluster::ENG_CALIB_TOT, mx ) ){
1058 engClusSumCalib += mx;
1059 }else{
1060 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Can't retrieve moment xAOD::CaloCluster::ENG_CALIB_TOT ");
1061 }
1062 double mx_calib_emb0, mx_calib_eme0, mx_calib_tileg3;
1063 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EMB0, mx_calib_emb0)
1064 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EME0, mx_calib_eme0)
1065 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TILEG3, mx_calib_tileg3)){
1066 ATH_MSG_WARNING("One of the moment ENG_CALIB_EMB0, ENG_CALIB_EME0, ENG_CALIB_TILEG3 is absent");
1067 }else{
1068 engClusSumCalibPresOnly += (mx_calib_emb0+mx_calib_eme0+mx_calib_tileg3);
1069 }
1070 } // iClus
1071
1072 double engCalibDeadTotWithPres = engCalibDeadTot + engClusSumCalibPresOnly;
1073
1074 // DM energy before barrel presampler, inside it, and between presampler and strips
1075 double eng_calib_dead_emb0 = engCalibDeadTotInArea[CaloDmDescrArea::DMA_EMB0]
1076 + engCalibDeadTotInArea[CaloDmDescrArea::DMA_EMB1]
1077 + engCalibTotSmp[CaloSampling::PreSamplerB];
1078 // DM energy between barrel and tile
1079 double eng_calib_dead_tile0 = engCalibDeadTotInArea[CaloDmDescrArea::DMA_EMB3_TILE0];
1080 // DM energy before scintillator and inside scintillator
1081 double eng_calib_dead_tileg3 = engCalibDeadTotInArea[CaloDmDescrArea::DMA_SCN]
1082 + engCalibTotSmp[CaloSampling::TileGap3];
1083 // DM energy beforee endcap presampler, inside it and between presampler and strips
1084 double eng_calib_dead_eme0 = engCalibDeadTotInArea[CaloDmDescrArea::DMA_EME0]
1085 + engCalibDeadTotInArea[CaloDmDescrArea::DMA_EME12]
1086 + engCalibTotSmp[CaloSampling::PreSamplerE];
1087 // DM energy between emec and hec
1088 double eng_calib_dead_hec0 = engCalibDeadTotInArea[CaloDmDescrArea::DMA_EME3_HEC0];
1089 // DM energy before FCAL and between HEC and FCAL
1090 double eng_calib_dead_fcal = engCalibDeadTotInArea[CaloDmDescrArea::DMA_FCAL0]
1091 + engCalibDeadTotInArea[CaloDmDescrArea::DMA_HEC_FCAL];
1092 // DM leakage behind the calorimeter
1093 double eng_calib_dead_leakage = engCalibDeadTotInArea[CaloDmDescrArea::DMA_LEAK];
1094 // the rest of DM energy which remains unclassified
1095 double eng_calib_dead_unclass = engCalibDeadTotWithPres - eng_calib_dead_emb0 - eng_calib_dead_tile0
1096 - eng_calib_dead_tileg3 - eng_calib_dead_eme0 - eng_calib_dead_hec0 - eng_calib_dead_fcal
1097 - eng_calib_dead_leakage;
1098
1099 // filling histograms
1100 int iMoment=0;
1101 for(moment_name_vector::const_iterator im=m_validMoments.begin(); im!=m_validMoments.end(); ++im, ++iMoment){
1102 double xnorm = 0.0;
1103 switch ( (*im).second ) {
1105 xnorm = clsMomentsSum[iMoment];
1106 break;
1108 xnorm = engCalibTot - engClusSumCalib;
1109 break;
1111 xnorm = engCalibTot - engClusSumCalib;
1112 break;
1114 xnorm = engCalibTot - engClusSumCalib;
1115 break;
1117 xnorm = engCalibDeadTotWithPres;
1118 break;
1120 xnorm = eng_calib_dead_emb0;
1121 break;
1123 xnorm = eng_calib_dead_tile0;
1124 break;
1126 xnorm = eng_calib_dead_tileg3;
1127 break;
1129 xnorm = eng_calib_dead_eme0;
1130 break;
1132 xnorm = eng_calib_dead_hec0;
1133 break;
1135 xnorm = eng_calib_dead_fcal;
1136 break;
1138 xnorm = eng_calib_dead_leakage;
1139 break;
1141 xnorm = eng_calib_dead_unclass;
1142 break;
1143 default:
1144 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Not implemented for " << (*im).first << " " << (*im).second);
1145 break;
1146 }
1147 if(m_doClusMoments && xnorm > m_mc_ener*0.0001 && xnorm > 0) {
1148 // moments assigned to first 3 maximum clusters
1149 const double inv_xnorm = 1. / xnorm;
1150 for(unsigned int i_cls=0; i_cls<clusColl.size(); i_cls++){
1151 m_clusMoment_vs_eta[iMoment][i_cls][m_mc_enerbin]->Fill(m_mc_eta, clsMoments[i_cls][iMoment]*inv_xnorm);
1152 m_clusMoment_vs_ebeam[iMoment][i_cls][m_mc_etabin]->Fill(m_mc_ener, clsMoments[i_cls][iMoment]*inv_xnorm);
1153 if(i_cls>=2) break;
1154 }
1155 // sum of moments assigned to clusters wrt total sum of moment
1156 m_clusMoment_vs_eta[iMoment][kMOM_SUMCLS][m_mc_enerbin]->Fill(m_mc_eta, clsMomentsSum[iMoment]/xnorm);
1157 m_clusMoment_vs_ebeam[iMoment][kMOM_SUMCLS][m_mc_etabin]->Fill(m_mc_ener, clsMomentsSum[iMoment]/xnorm);
1158 // unassigned energy wrt pion energy
1159 m_clusMoment_vs_eta[iMoment][kMOM_MISS][m_mc_enerbin]->Fill(m_mc_eta, (xnorm-clsMomentsSum[iMoment])/m_mc_ener);
1160 m_clusMoment_vs_ebeam[iMoment][kMOM_MISS][m_mc_etabin]->Fill(m_mc_ener, (xnorm-clsMomentsSum[iMoment])/m_mc_ener);
1161 // unassigned energy wrt pion energy
1164 }
1165 } // i_mom
1166
1167 return 0;
1168}
1169
1170
1171
1172/* ****************************************************************************
1173calibration hits validation
1174**************************************************************************** */
1176 const EventContext& ctx)
1177{
1178 /* ********************************************
1179 reading calibration containers
1180 ******************************************** */
1181 std::vector<const CaloCalibrationHitContainer *> v_cchc;
1184 v_cchc.push_back(cchc.cptr());
1185 }
1186
1187 std::vector<const CaloCalibrationHitContainer *> v_dmcchc;
1190 v_dmcchc.push_back(cchc.cptr());
1191 }
1192
1193 // calculate total calibration energy int active+inactive hits
1194 double engCalibTot = 0.0;
1195 for (const CaloCalibrationHitContainer* cchc : v_cchc) {
1196 //loop over cells in calibration container
1197 for (const CaloCalibrationHit* hit : *cchc) {
1198 Identifier myId = hit->cellID();
1199 if(!myId.is_valid()) {
1200 ATH_MSG_ERROR("Error! Bad identifier " << myId << " in container '" << cchc->Name() << "',"
1201 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1202 continue;
1203 }
1204 if(!m_id_helper->is_lar(myId) && !m_id_helper->is_tile(myId)) {
1205 ATH_MSG_ERROR("Error! Bad identifier (nor lar or tile) " << myId << " in container '" << cchc->Name() << "',"
1206 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1207 continue;
1208 }
1209 engCalibTot += hit->energyTotal();
1210 }
1211 }
1212
1213 // calibration energy in dead material
1214 double engCalibDeadTot = 0.0;
1215 double engDefaultCalculator = 0.0;
1216 for (const CaloCalibrationHitContainer* dmcchc : v_dmcchc) {
1217 //loop over cells in calibration container
1218 for (const CaloCalibrationHit* hit : *dmcchc) {
1219 Identifier myId = hit->cellID();
1220 if(!myId.is_valid()) {
1221 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Bad dead material identifier " << myId << " in container '" << dmcchc->Name() << "',"
1222 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1223 continue;
1224 }
1225 if(!m_id_helper->is_lar_dm(myId) && !m_id_helper->is_tile_dm(myId)) {
1226 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Bad dead material identifier (nor lar_dm or tile_dm) " << myId << " in container '" << dmcchc->Name() << "',"
1227 << " AtlasDetectorID says '" << m_id_helper->show_to_string(myId) << "'");
1228 continue;
1229 }
1230 engCalibDeadTot += hit->energyTotal();
1231 //
1232 CaloDmDescrElement* myCDDE(nullptr);
1233 myCDDE = m_caloDmDescrManager->get_element(myId);
1234 if( !myCDDE ) {
1235 ATH_MSG_ERROR("GetLCSinglePionsPerf::fill_moments() -> Error! Absent DM description element!");
1236 continue;
1237 }
1238 int nDmArea = m_caloDmDescrManager->get_dm_area(myId);
1239 if(nDmArea == CaloDmDescrArea::DMA_DEFCALC) engDefaultCalculator += hit->energyTotal();
1240 }
1241 }
1242
1243 // calibration energy assigned to one or another cluster
1244 // retrieving cluster moments
1245 double engCalibAssigned = 0.0;
1246 double engClusSumCalibPresOnly = 0.0;
1247 for (const xAOD::CaloCluster* theCluster : clusColl) {
1248 double mx_calib_tot, mx_calib_ooc, mx_calib_dm;
1249// if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT, mx_calib_tot)
1250// || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_OUT_L, mx_calib_ooc)
1251// || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_DEAD_TOT mx_calib_dm)){
1252// ATH_MSG_ERROR("One of the moment ENG_CALIB_TOT, ENG_CALIB_OUT_L, ENG_CALIB_DEAD_TOT is absent");
1253// }
1254 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT, mx_calib_tot) ){
1255 ATH_MSG_WARNING("Moment ENG_CALIB_TOT is absent");
1256 }
1257 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_OUT_L, mx_calib_ooc) ){
1258 ATH_MSG_WARNING("Moment ENG_CALIB_OUT_L is absent" );
1259 }
1260 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_DEAD_TOT, mx_calib_dm) ){
1261 ATH_MSG_WARNING("Moment ENG_CALIB_DEAD_TOT is absent");
1262 }
1263 engCalibAssigned += (mx_calib_tot + mx_calib_ooc + mx_calib_dm);
1264
1265 double mx_calib_emb0, mx_calib_eme0, mx_calib_tileg3;
1266 if( !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EMB0, mx_calib_emb0)
1267 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_EME0, mx_calib_eme0)
1268 || !theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TILEG3, mx_calib_tileg3)){
1269 ATH_MSG_WARNING("One of the moment ENG_CALIB_EMB0, ENG_CALIB_EME0, ENG_CALIB_TILEG3 is absent");
1270 }else{
1271 engClusSumCalibPresOnly += (mx_calib_emb0+mx_calib_eme0+mx_calib_tileg3);
1272 }
1273
1274 } // iClus
1275 // energy in presamplers is included already in the dead material
1276 engCalibAssigned -= engClusSumCalibPresOnly;
1277
1278 int n_enerbin = int( (log10(m_mc_ener) - m_logenermin)/ ((m_logenermax - m_logenermin)/m_nlogenerbin2));
1279 // sum of moments assigned to clusters wrt total sum of moment
1280 m_SumCalibHitOverEbeam_vs_eta[n_enerbin]->Fill(m_mc_eta, (engCalibTot+engCalibDeadTot)/m_mc_ener);
1281// m_SumCalibHitOverEbeam_vs_ebeam[m_mc_etabin]->Fill(m_mc_ener, (engCalibTot+engCalibDeadTot)/m_mc_ener);
1282 m_SumCalibHitAssignedOverEbeam_vs_eta[n_enerbin]->Fill(m_mc_eta, engCalibAssigned/m_mc_ener);
1283 m_SumCalibHitAssignedOverEbeam_vs_etaphi[n_enerbin]->Fill(m_mc_eta, m_mc_phi, engCalibAssigned/m_mc_ener);
1284
1285 m_DefaultCalculatorOverEbeam_vs_eta[n_enerbin]->Fill(m_mc_eta, engDefaultCalculator/m_mc_ener);
1286// m_DefaultCalculatorOverEbeam_vs_ebeam[m_mc_etabin]->Fill(m_mc_ener, engDefaultCalculator/m_mc_ener);
1287
1288 return 0;
1289}
1290
1291
1292
1293/* ****************************************************************************
1294
1295**************************************************************************** */
1297{
1298 double eta = fabs(x);
1299 double ff;
1300 if(eta<1.6){
1301 ff = atan(5.0*1.7/(200.0*cosh(eta)));
1302 }else if(eta<3.2){
1303 ff = atan(5.0*1.6/(420./tanh(eta)));
1304 }else{
1305 ff = atan(5.0*0.95/(505./tanh(eta)));
1306 }
1307 return ff*(1./atan(5.0*1.7/200.0));
1308}
1309
1310
1311
#define M_PI
Scalar eta() const
pseudorapidity method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
Helpers for checking error return status codes and reporting errors.
ATLAS-specific HepMC functions.
#define x
constexpr int pow(int base, int exp) noexcept
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
const ServiceHandle< StoreGateSvc > & detStore() const
Class to store calorimeter calibration hit.
static const CaloDmDescrManager * instance()
reconstruction status indicator
virtual const store_type & getStatusWord() const
retrieve the entire status word
virtual bool checkStatus(const StatusIndicator &statusIndicator) const
Check status.
const T * at(size_type n) const
Access an element, as an rvalue.
size_type size() const noexcept
Returns the number of elements in the collection.
std::vector< TProfile2D * > m_SumCalibHitAssignedOverEbeam_vs_etaphi
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_DMCalibrationHitContainerNames
SG::ReadHandleKeyArray< xAOD::CaloClusterContainer > m_clusterCollNames
virtual StatusCode finalize()
int fill_moments(const xAOD::CaloClusterContainer &clusColl, const EventContext &ctx)
std::vector< TProfile * > m_DefaultCalculatorOverEbeam_vs_eta
virtual StatusCode initialize()
const CaloCell_ID * m_calo_id
moment_name_vector m_validMoments
std::vector< TH1F * > m_RecoEfficiency_vs_eta[2]
std::vector< TH1F * > m_RecoEfficiency_vs_ebeam[2]
const AtlasDetectorID * m_id_helper
std::vector< std::vector< std::vector< TProfile * > > > m_engRecOverTruth_vs_ebeam
std::vector< std::vector< TH1F * > > m_engPionSpect
std::vector< std::vector< std::vector< TProfile * > > > m_engRecOverTruth_vs_eta
static double angle_mollier_factor(double x)
std::vector< TProfile * > m_SumCalibHitAssignedOverEbeam_vs_eta
SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterBasicCollName
std::vector< std::vector< TH1F * > > m_engNoiseClusSpect
std::pair< std::string, xAOD::CaloCluster::MomentType > moment_name_pair
std::vector< std::vector< TProfile * > > m_engTag_vs_eta
std::vector< std::vector< std::vector< TProfile * > > > m_clusMoment_vs_eta
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_CalibrationHitContainerNames
std::vector< TProfile * > m_SumCalibHitOverEbeam_vs_eta
std::vector< std::vector< std::vector< TH1F * > > > m_engRecSpect
GetLCSinglePionsPerf(const std::string &name, ISvcLocator *pSvcLocator)
int fill_reco(const xAOD::CaloClusterContainer &clusColl, const EventContext &ctx)
const CaloDmDescrManager * m_caloDmDescrManager
virtual StatusCode execute()
std::vector< TProfile * > m_engNoiseClus_vs_eta
std::vector< std::vector< std::vector< TProfile * > > > m_clusMoment_vs_ebeam
int fill_calibhits(const xAOD::CaloClusterContainer &clusColl, const EventContext &ctx)
std::vector< std::vector< TProfile * > > m_engTag_vs_ebeam
bool is_valid() const
Check if id is in a valid state.
This defines the McEventCollection, which is really just an ObjectVector of McEvent objectsFile: Gene...
Property holding a SG store/key/clid from which a ReadHandle is made.
const_pointer_type cptr()
Dereference the pointer.
@ ENG_CALIB_OUT_M
Attached Calibration Hit energy outside clusters but inside the calorimeter with medium matching (Ang...
@ ENG_CALIB_OUT_L
Attached Calibration Hit energy outside clusters but inside the calorimeter with loose matching (Angl...
@ ENG_CALIB_DEAD_UNCLASS
Attached Calibration Hit energy in dead material in unclassified areas of the detector.
@ ENG_CALIB_DEAD_HEC0
Attached Calibration Hit energy in dead material between EME3 and HEC0.
@ FIRST_ENG_DENS
First Moment in E/V.
@ ENG_CALIB_DEAD_TILEG3
Attached Calibration Hit energy in dead material before scintillator.
@ CENTER_LAMBDA
Shower depth at Cluster Centroid.
@ ENG_CALIB_DEAD_EME0
Attached Calibration Hit energy in dead material before EME0, between EME0 and EME1.
@ ENG_CALIB_DEAD_TILE0
Attached Calibration Hit energy in dead material between EMB3 and TILE0.
@ ENG_CALIB_DEAD_FCAL
Attached Calibration Hit energy in dead material before FCAL, between FCAL and HEC.
@ ENG_CALIB_TOT
Calibration Hit energy inside the cluster.
@ ENG_CALIB_OUT_T
Attached Calibration Hit energy outside clusters but inside the calorimeter with tight matching (Angl...
@ ENG_CALIB_DEAD_LEAKAGE
Attached Calibration Hit energy in dead material behind calorimeters.
@ ENG_CALIB_EMB0
Calibration Hit energy inside the cluster barrel presampler.
@ ENG_CALIB_EME0
Calibration Hit energy inside the cluster endcap presampler.
@ ENG_CALIB_DEAD_EMB0
Attached Calibration Hit energy in dead material before EMB0, between EMB0 and EMB1.
@ ENG_CALIB_DEAD_TOT
Attached Calibration Hit energy in dead material.
@ ENG_CALIB_TILEG3
Calibration Hit energy inside the cluster scintillator.
int r
Definition globals.cxx:22
const GenParticle * ConstGenParticlePtr
Definition GenParticle.h:38
static const int PI0
static const int PIPLUS
bool isGenStable(const T &p)
Determine if the particle is stable at the generator (not det-sim) level,.
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
CaloClusterContainer_v1 CaloClusterContainer
Define the latest version of the calorimeter cluster container.
TChain * tree