ATLAS Offline Software
Loading...
Searching...
No Matches
CaloHadDMCoeffCheck.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: CaloHadDMCoeffCheck.cxx,v 1.1 2009/03/03 17:30:23 pospelov Exp $
8//
9// Description: see CaloHadDMCoeffCheck.h
10//
11// Environment:
12// Software developed for the ATLAS Detector at CERN LHC
13//
14// Author List:
15// Gennady Pospelov
16//
17//-----------------------------------------------------------------------
18#include "CaloLocalHadCalib/CaloHadDMCoeffCheck.h"
19#include "CaloConditions/CaloLocalHadCoeff.h"
20#include "CaloConditions/CaloLocalHadDefs.h"
21#include "CaloGeoHelpers/CaloSampling.h"
22#include "CaloLocalHadCalib/CaloHadDMCoeffData.h"
23#include "CaloLocalHadCalib/CaloHadDMCoeffFit.h"
24#include "CaloLocalHadCalib/CaloLocalHadCoeffHelper.h"
25#include "CaloLocalHadCalib/GetLCSinglePionsPerf.h"
26#include <algorithm>
27#include <cmath>
28#include <fstream>
29#include <iostream>
30#include <sstream>
31
32#include <CLHEP/Vector/LorentzVector.h>
33#include <CLHEP/Units/SystemOfUnits.h>
34
35#include "TROOT.h"
36#include "TStyle.h"
37#include "TError.h"
38#include "TFile.h"
39#include "TChain.h"
40#include "TH1.h"
41#include "TH1F.h"
42#include "TF1.h"
43#include "TProfile.h"
44#include "TH2F.h"
45#include "TCanvas.h"
46#include "TPad.h"
47#include "TText.h"
48#include "TLatex.h"
49#include "TGraphErrors.h"
50
51
52using CLHEP::HepLorentzVector;
53using CLHEP::MeV;
54
55
56CaloHadDMCoeffCheck::CaloHadDMCoeffCheck() :
57 m_isTestbeam(false),
58 m_energyMin(200*MeV),
59 m_netabin(25), m_etamin(0.0), m_etamax(5.0),
60 m_deta(0),
61 m_nphibin(1), m_phimin(-M_PI), m_phimax(M_PI),
62 m_dphi(0),
63 m_nlogenerbin(22), m_logenermin(2.0), m_logenermax(6.4),
64 m_dlogener(0),
65 m_nrecobin(3), m_npdg(2), m_doTailRejection(true)
66{
67 m_data = nullptr;
68 m_HadDMCoeff = nullptr;
69 m_HadDMHelper = new CaloLocalHadCoeffHelper();
70
71}
72
73
78
79
80int CaloHadDMCoeffCheck::process(CaloHadDMCoeffData *myData, CaloLocalHadCoeff *myHadDMCoeff, bool isSingleParticle, bool tbflag)
81{
82 m_isTestbeam = tbflag;
83
84 std::cout << std::endl;
85 std::cout << std::endl;
86 std::cout << "--- CaloHadDMCoeffCheck::process() --- " << std::endl;
87
88 if(m_isTestbeam){
89 m_etamin = 2.5;
90 m_etamax = 5.0;
91 m_netabin = 25;
92 m_phimin = M_PI*3./4.-0.15;
93 m_phimax = M_PI*3./4.+0.15;
94 m_nphibin = 1;
95 m_logenermin = 2.1;
96 m_logenermax = 5.5;
97 m_nlogenerbin = 17;
98 }
99 m_deta = (m_etamax - m_etamin)/m_netabin;
100 m_dphi = (m_phimax - m_phimin)/m_nphibin;
101 m_dlogener = (m_logenermax - m_logenermin)/m_nlogenerbin;
102
103 m_data = myData;
104 m_HadDMCoeff = myHadDMCoeff;
105
106 // loading only necessary branches
107 m_data->fChain->SetBranchStatus("*",0);
108 m_data->fChain->SetBranchStatus("mc_ener",1);
109 m_data->fChain->SetBranchStatus("mc_eta",1);
110 m_data->fChain->SetBranchStatus("mc_phi",1);
111 m_data->fChain->SetBranchStatus("mc_pdg",1);
112 m_data->fChain->SetBranchStatus("ncls",1);
113 m_data->fChain->SetBranchStatus("cls_eta",1);
114 m_data->fChain->SetBranchStatus("cls_phi",1);
115 m_data->fChain->SetBranchStatus("cls_lambda",1);
116 m_data->fChain->SetBranchStatus("cls_calib_emfrac",1);
117 m_data->fChain->SetBranchStatus("cls_engcalib",1);
118 m_data->fChain->SetBranchStatus("engClusSumCalib",1);
119 m_data->fChain->SetBranchStatus("cls_ener_unw",1);
120 m_data->fChain->SetBranchStatus("narea",1);
121 m_data->fChain->SetBranchStatus("cls_eprep",1);
122 m_data->fChain->SetBranchStatus("cls_dmener",1);
123 m_data->fChain->SetBranchStatus("cls_smpener_unw",1);
124 m_data->fChain->SetBranchStatus("cls_oocener",1);
125 m_data->fChain->SetBranchStatus("cls_engcalibpres",1);
126
127 if( !m_data->GetEntries() ) {
128 std::cout << "CaloHadDMCoeffCheck::process -> Error! No entries in DeadMaterialTree." << std::endl;
129 return 0;
130 }
131
132 std::vector<std::vector<std::vector<CaloHadDMCoeffFit::PrepData *> > > ereco; // [sDM][m_netabin][m_nlogenerbin]
133 std::vector<std::vector<std::vector<CaloHadDMCoeffFit::PrepData *> > > etrue; // [sDM][m_netabin][m_nlogenerbin]
134 int nArea = m_HadDMCoeff->getSizeAreaSet()+1;
135 ereco.resize(nArea);
136 etrue.resize(nArea);
137 for(int i_area=0; i_area<nArea; i_area++){
138 ereco[i_area].resize(m_netabin);
139 etrue[i_area].resize(m_netabin);
140 for(int i_eta=0; i_eta<m_netabin; i_eta++){
141 ereco[i_area][i_eta].resize(m_nlogenerbin, nullptr);
142 etrue[i_area][i_eta].resize(m_nlogenerbin, nullptr);
143 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
144 ereco[i_area][i_eta][i_ener] = new CaloHadDMCoeffFit::PrepData();
145 etrue[i_area][i_eta][i_ener] = new CaloHadDMCoeffFit::PrepData();
146 } // i_ener
147 } // i_eta
148 } // i_area
149
150 /* ********************************************
151 first loop to calculate histogram limits
152 this code reproduces CaloUtils/src/CaloLCDeadMaterialTool.cxx
153 ********************************************* */
154 std::cout << "CaloHadDMCoeffCheck::process() -> Info. First loop to find histogram limits " << std::endl;
155 for(int i_ev=0; m_data->GetEntry(i_ev)>0;i_ev++) {
156 if(i_ev%10000==0) std::cout << " i_ev: " << i_ev << " '" << static_cast<TChain *>(m_data->fChain)->GetFile()->GetName() << "'" << std::endl;
157 if(m_data->m_mc_ener <= 0.0) {
158 std::cout << "CaloHadDMCoeffCheck::process() -> Warning! Unknown particle energy " << m_data->m_mc_ener << std::endl;
159 continue;
160 }
161
162 // checking event quality
163 if(isSingleParticle) {
164 bool GoodClusterFound(false);
165 if( m_data->m_ncls ) {
166 HepLorentzVector hlv_pion(1,0,0,1);
167 hlv_pion.setREtaPhi(1./cosh(m_data->m_mc_eta), m_data->m_mc_eta, m_data->m_mc_phi);
168 for(int i_cls=0; i_cls<m_data->m_ncls; i_cls++){ // loop over clusters
169 HepLorentzVector hlv_cls(1,0,0,1);
170 hlv_cls.setREtaPhi(1./cosh( (*m_data->m_cls_eta)[i_cls] ), (*m_data->m_cls_eta)[i_cls], (*m_data->m_cls_phi)[i_cls]);
171 double r = hlv_pion.angle(hlv_cls.vect());
172 if(r < 1.5*GetLCSinglePionsPerf::angle_mollier_factor(m_data->m_mc_eta)
173 && (*m_data->m_cls_engcalib)[i_cls] > 20.0*MeV
174 //&& (*m_data->m_cls_ener)[i_cls] > 0.01*m_data->m_mc_ener
175 ) {
176 GoodClusterFound = true;
177 break;
178 }
179 } // i_cls
180 }
181 if(!GoodClusterFound) continue;
182 }
183
184 int mc_enerbin = int( (log10(m_data->m_mc_ener) - m_logenermin)/m_dlogener);
185 int mc_etabin = int((fabs(m_data->m_mc_eta)-m_etamin)/m_deta);
186 int mc_phibin = int((fabs(m_data->m_mc_phi)-m_phimin)/m_dphi);
187
188 if(mc_etabin <0 || mc_etabin >= m_netabin || mc_enerbin <0 || mc_enerbin >= m_nlogenerbin || mc_phibin!=0) continue;
189 if(m_data->m_engClusSumCalib <=0.0) continue;
190
191 std::vector<std::vector<double > > engDmReco; // [ncls][narea]
192 getDmReco(engDmReco);
193 std::vector<double > engDmRecSumClus;
194 std::vector<double > engDmTrueSumClus;
195 engDmRecSumClus.resize(nArea, 0.0);
196 engDmTrueSumClus.resize(nArea, 0.0);
197 for(int i_cls=0; i_cls<m_data->m_ncls; i_cls++){ // loop over clusters
198 for(int i_area=0; i_area<m_HadDMCoeff->getSizeAreaSet(); i_area++){
199 engDmRecSumClus[i_area] += engDmReco[i_cls][i_area];
200 engDmTrueSumClus[i_area] += (*m_data->m_cls_dmener)[i_cls][i_area];
201 engDmRecSumClus[m_HadDMCoeff->getSizeAreaSet()] += engDmReco[i_cls][i_area];
202 engDmTrueSumClus[m_HadDMCoeff->getSizeAreaSet()] += (*m_data->m_cls_dmener)[i_cls][i_area];
203 }
204 }
205 for(int i_area=0; i_area<m_HadDMCoeff->getSizeAreaSet()+1; i_area++){
206 ereco[i_area][mc_etabin][mc_enerbin]->add(engDmRecSumClus[i_area]);
207 etrue[i_area][mc_etabin][mc_enerbin] -> add(engDmTrueSumClus[i_area]);
208 }
209 } // i_ev
210
211
212 /* ********************************************
213 histogram creation
214 ******************************************** */
215 char hname[256];
216 m_h2_etrue_vs_ereco.resize(nArea);
217 m_hp_etrue_vs_ereco.resize(nArea);
218 for(int i_area=0; i_area<nArea; i_area++){
219 m_h2_etrue_vs_ereco[i_area].resize(m_netabin);
220 m_hp_etrue_vs_ereco[i_area].resize(m_netabin);
221 for(int i_eta=0; i_eta<m_netabin; i_eta++){
222 m_h2_etrue_vs_ereco[i_area][i_eta].resize(m_nlogenerbin, nullptr);
223 m_hp_etrue_vs_ereco[i_area][i_eta].resize(m_nlogenerbin, nullptr);
224 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
225 float elimreco = ereco[i_area][i_eta][i_ener]->m_aver + 5.*sqrt(ereco[i_area][i_eta][i_ener]->m_rms);
226 float elimtrue = etrue[i_area][i_eta][i_ener]->m_aver + 5.*sqrt(etrue[i_area][i_eta][i_ener]->m_rms);
227 if(elimreco <=0.0) elimreco = 1.0;
228 if(elimtrue <=0.0) elimtrue = 1.0;
229
230 sprintf(hname,"h2_etrue_vs_ereco_dm%d_eta%d_ener%d",i_area, i_eta, i_ener);
231 m_h2_etrue_vs_ereco[i_area][i_eta][i_ener] = new TH2F(hname, hname, 30, 0.0, elimtrue, 30, 0.0, elimreco);
232 sprintf(hname,"hp_etrue_vs_ereco_dm%d_eta%d_ener%d",i_area, i_eta, i_ener);
233 m_hp_etrue_vs_ereco[i_area][i_eta][i_ener] = new TProfile(hname, hname, 30, 0.0, elimtrue);
234 m_hp_etrue_vs_ereco[i_area][i_eta][i_ener] -> SetMinimum(0.0);
235 m_hp_etrue_vs_ereco[i_area][i_eta][i_ener] -> SetMaximum(elimreco);
236 } // i_ener
237 } // i_eta
238 } // i_area
239
240 // reco over truth
241 m_hp_engRecoOverTruth_vs_eta.resize(m_npdg);
242 for(int i_pdg=0; i_pdg<m_npdg; i_pdg++){
243 m_hp_engRecoOverTruth_vs_eta[i_pdg].resize(nArea);
244 for(int i_area=0; i_area<nArea; i_area++){
245 m_hp_engRecoOverTruth_vs_eta[i_pdg][i_area].resize(m_nrecobin);
246 for(int i_reco=0; i_reco<m_nrecobin; i_reco++){
247 m_hp_engRecoOverTruth_vs_eta[i_pdg][i_area][i_reco].resize(m_nlogenerbin);
248 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
249 sprintf(hname,"m_hp_engRecoOverTruth_vs_eta_pdg%d_dm%d_reco%d_ener%d",i_pdg,i_area, i_reco, i_ener);
250 m_hp_engRecoOverTruth_vs_eta[i_pdg][i_area][i_reco][i_ener] = new TProfile(hname, hname, 50, m_etamin, m_etamax);
251 } // i_ener
252 } // m_nrecobin
253 } // i_area
254 } // i_pdg
255
256 // reco spect
257 m_engRecSpect.resize(m_npdg);
258 for(int i_pdg=0; i_pdg<m_npdg; i_pdg++){
259 m_engRecSpect[i_pdg].resize(nArea);
260 for(int i_area=0; i_area<nArea; i_area++){
261 m_engRecSpect[i_pdg][i_area].resize(m_nrecobin);
262 for(int i_reco=0; i_reco<m_nrecobin; i_reco++){
263 m_engRecSpect[i_pdg][i_area][i_reco].resize(m_nlogenerbin);
264 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
265 m_engRecSpect[i_pdg][i_area][i_reco][i_ener].resize(m_netabin, nullptr);
266 for(int i_eta=0; i_eta<m_netabin; i_eta++){
267 sprintf(hname,"m_engRecSpect_pdg%d_dm%d_reco%d_ener%d_eta%d",i_pdg,i_area, i_reco, i_ener, i_eta);
268 TH1F *h1 = new TH1F(hname, hname, 110, -0.2, 2.0);
269 m_engRecSpect[i_pdg][i_area][i_reco][i_ener][i_eta] = h1;
270 } // i_eta
271 } // i_ener
272 } // i_reco
273 } // i_area
274 } // i_pdg
275
276 /* ********************************************
277 histogram filling
278 ******************************************** */
279 std::cout << "CaloHadDMCoeffCheck::process() -> Info. Second loop to fill histogram " << std::endl;
280 for(int i_ev=0; m_data->GetEntry(i_ev)>0;i_ev++) {
281 if(i_ev%10000==0) std::cout << " i_ev: " << i_ev << " '" << static_cast<TChain *>(m_data->fChain)->GetFile()->GetName() << "'" << std::endl;
282
283 int mc_enerbin = int( (log10(m_data->m_mc_ener) - m_logenermin)/m_dlogener);
284 int mc_etabin = int((fabs(m_data->m_mc_eta)-m_etamin)/m_deta);
285 int mc_phibin = int((fabs(m_data->m_mc_phi)-m_phimin)/m_dphi);
286 if(mc_etabin <0 || mc_etabin >= m_netabin || mc_enerbin <0 || mc_enerbin >= m_nlogenerbin || mc_phibin!=0) continue;
287 if(m_data->m_engClusSumCalib <=0.0) continue;
288
289 int i_pdg=0;
290 if( abs(m_data->m_mc_pdg)==211) {
291 i_pdg=0;
292 }else{
293 i_pdg=1;
294 }
295
296 std::vector<std::vector<double > > engDmReco; // [ncls][narea]
297 getDmReco(engDmReco);
298 std::vector<double > engDmRecSumClus;
299 std::vector<double > engDmTrueSumClus;
300 engDmRecSumClus.resize(nArea, 0.0);
301 engDmTrueSumClus.resize(nArea, 0.0);
302 double engClusSumCalib = 0.0;
303 double engClusSumOOC = 0.0;
304 double engClusSumDM = 0.0;
305 double engClusSumCalibInPresampler = 0.0;
306
307 for(int i_cls=0; i_cls<m_data->m_ncls; i_cls++){ // loop over clusters
308 engClusSumCalib += m_data->m_cls_engcalib ? (*m_data->m_cls_engcalib)[i_cls] : 0.0;
309 engClusSumOOC += m_data->m_cls_oocener ? (*m_data->m_cls_oocener)[i_cls] : 0.0;
310 engClusSumCalibInPresampler += m_data->m_cls_engcalibpres ? (*m_data->m_cls_engcalibpres)[i_cls] : 0.0;
311 for(int i_area=0; i_area<m_HadDMCoeff->getSizeAreaSet(); i_area++){
312 engDmRecSumClus[i_area] += engDmReco[i_cls][i_area];
313 engDmTrueSumClus[i_area] += m_data->m_cls_dmener ? (*m_data->m_cls_dmener)[i_cls][i_area] : 0.0;
314 engDmRecSumClus[m_HadDMCoeff->getSizeAreaSet()] += engDmReco[i_cls][i_area];
315 engDmTrueSumClus[m_HadDMCoeff->getSizeAreaSet()] +=
316 m_data->m_cls_dmener ? (*m_data->m_cls_dmener)[i_cls][i_area] : 0.0;
317 engClusSumDM += m_data->m_cls_dmener ? (*m_data->m_cls_dmener)[i_cls][i_area] : 0.0;
318 }
319 }
320 for(int i_area=0; i_area<m_HadDMCoeff->getSizeAreaSet()+1; i_area++){
321 m_h2_etrue_vs_ereco[i_area][mc_etabin][mc_enerbin]->Fill(engDmTrueSumClus[i_area], engDmRecSumClus[i_area]);
322 m_hp_etrue_vs_ereco[i_area][mc_etabin][mc_enerbin]->Fill(engDmTrueSumClus[i_area], engDmRecSumClus[i_area]);
323 }
324
325 // engRecoOverTruth
326 double sum = engClusSumCalib + engClusSumOOC + engClusSumDM - engClusSumCalibInPresampler; // presampler is already included in the dead material
327 for(int i_area=0; i_area<m_HadDMCoeff->getSizeAreaSet()+1; i_area++){
328 for(int i_reco=0; i_reco<m_nrecobin; i_reco++){
329 double engReco = 0.0;
330 if(i_reco == kENG_RECO){ // red
331 engReco = sum - engDmTrueSumClus[i_area] + engDmRecSumClus[i_area];
332 }else if(i_reco == kENG_TRUTH){ // blue
333 engReco = sum;
334 }else if(i_reco == kENG_NORECO){ // green
335 engReco = sum - engDmTrueSumClus[i_area];
336 }else{
337 std::cout << "panic" << std::endl;
338 }
339 m_hp_engRecoOverTruth_vs_eta[i_pdg][i_area][i_reco][mc_enerbin]->Fill(fabs(m_data->m_mc_eta), engReco/m_data->m_mc_ener);
340 m_engRecSpect[i_pdg][i_area][i_reco][mc_enerbin][mc_etabin]->Fill(engReco/m_data->m_mc_ener);
341 } // i_reco
342 } // i_area
343
344 } // i_ev
345
346 /* ********************************************
347 histogram fitting
348 ******************************************** */
349 TF1 *fun_p1 = new TF1("fun_udeg3","[0]+[1]*x",1.,200000.);
350 for(int i_area=0; i_area<nArea; i_area++){
351 for(int i_eta=0; i_eta<m_netabin; i_eta++){
352 for(int i_ener=0; i_ener<m_nlogenerbin; i_ener++){
353 TProfile *hp = m_hp_etrue_vs_ereco[i_area][i_eta][i_ener];
354 if(hp->GetEntries() < 100) continue;
355 float fitlim1 = hp->GetBinCenter(2);
356 float fitlim2 = hp->GetBinCenter(30);
357 hp->Fit(fun_p1,"0Q","",fitlim1,fitlim2); // silent fit, but now drawing option for fit results is switched off
358 TF1 *f=hp->GetFunction(fun_p1->GetName());
359 if(f) {
360 f->ResetBit(TF1::kNotDraw);
361 f->SetLineWidth(1);
362 f->SetLineColor(4);
363 }
364 } // i_ener
365 } // i_eta
366 } // i_area
367
368
369
370 return 0;
371}
372
373
374
375/* ****************************************************************************
376Making nice postscript report
377**************************************************************************** */
378void CaloHadDMCoeffCheck::make_report(std::string &sreport)
379{
380 const char *a_title[sDM] = {
381 "EMB0_EMB1", // 0 before PreSamplerB, between PreSamplerB and EMB1
382 "EMB3_TILE0", // 1 between barrel and tile
383 "SCN", // 2 before TileGap3 (scintillator)
384 "EME0_EME12", // 3 before PreSamplerE, between PreSamplerE and EME1
385 "EME3_HEC0", // 4 between EME3 and HEC0
386 "FCAL", // 5 between HEC and FCAL, before FCAL
387 "LEAKAGE", // 6 leakage behind calorimeter
388 "UNCLASS", // 7 unclassified DM enegry
389 "DMTOT"
390 };
391
392 int a_enerbin_sel[] = {12, 15};
393 int a_etabin_sel[] = {2, 7, 10, 18};
394 //int a_area_sel[] = {sDM_EME3_HEC0, sDM_FCAL, sDM_LEAKAGE, sDM_UNCLASS, sDM_DMTOT};
395 int a_area_sel[] = {
396 sDM_EMB0_EMB1, // 0 before PreSamplerB, between PreSamplerB and EMB1
397 sDM_EMB3_TILE0, // 1 between barrel and tile
398 sDM_SCN, // 2 before TileGap3 (scintillator)
399 sDM_EME0_EME12, // 3 before PreSamplerE, between PreSamplerE and EME1
400 sDM_EME3_HEC0, // 4 between EME3 and HEC0
401 sDM_FCAL, // 5 between HEC and FCAL, before FCAL
402 sDM_LEAKAGE, // 6 leakage behind calorimeter
403 sDM_UNCLASS, // 7 unclassified DM enegry
404 sDM_DMTOT // 8 sum of DM energy over all zones
405 };
406
407 std::cout << "CaloHadDMCoeffCheck::make_report() -> Info. Making report..." << std::endl;
408 gStyle->SetCanvasColor(10);
409 gStyle->SetCanvasBorderMode(0);
410 gStyle->SetPadColor(10);
411 gStyle->SetPadBorderMode(0);
412 gStyle->SetPalette(1);
413 gStyle->SetTitleBorderSize(1);
414 gStyle->SetTitleFillColor(10);
415 int cc_xx = 768, cc_yy = 1024;
416 char str[1024];
417 gROOT->SetBatch(kTRUE);
418 gErrorIgnoreLevel=3; // root global variables to supress text output in ->Print() methods
419 std::string sfname = sreport;
420 TCanvas *ctmp = new TCanvas("ctmp","ctmp", cc_xx, cc_yy);
421 sfname += "[";
422 ctmp->Print(sfname.c_str());
423 sfname = sreport;
424
425 TLatex tex;
426 tex.SetNDC(1);
427 tex.SetTextSize(0.05);
428
429 TCanvas *c1ps = new TCanvas("c1ps","c1ps", cc_xx, cc_yy);
430 for(int i_ener=0; i_ener<int(sizeof(a_enerbin_sel)/sizeof(int)); i_ener++){
431 int i_ener_sel = a_enerbin_sel[i_ener];
432 for(int i_eta=0; i_eta<int(sizeof(a_etabin_sel)/sizeof(int)); i_eta++){
433 int i_eta_sel = a_etabin_sel[i_eta];
434 c1ps->Clear();
435 c1ps->Divide(3,3);
436 for(int i_area=0; i_area<int(sizeof(a_area_sel)/sizeof(int)); i_area++){
437 int i_area_sel = a_area_sel[i_area];
438 c1ps->cd(i_area+1);
439 gPad->SetLeftMargin(0.12);
440 gPad->SetTopMargin(0.20);
441 sprintf(str,"%s ener%d eta%d dm%d",a_title[i_area_sel], i_ener_sel, i_eta_sel, i_area);
442 TH1 *h = m_h2_etrue_vs_ereco[i_area_sel][i_eta_sel][i_ener_sel];
443 h->SetTitle(str);
444 h->GetXaxis()->SetTitle("E_{DM} True");
445 h->GetYaxis()->SetTitle("E_{DM} Reco");
446 h->GetYaxis()->SetTitleOffset(1.6);
447 h->GetYaxis()->SetNdivisions(508);
448 h->GetXaxis()->SetNdivisions(508);
449 h->Draw("box");
450 h = m_hp_etrue_vs_ereco[i_area_sel][i_eta_sel][i_ener_sel];
451 h->SetLineColor(kRed);
452 h->SetMarkerColor(kRed);
453 h->Draw("same");
454 sprintf(str,"E = %3.1f-%3.1f GeV",
455 pow(10,m_logenermin+m_dlogener*i_ener_sel)*1e-3,
456 pow(10,m_logenermin+m_dlogener*(i_ener_sel+1))*1e-3);
457 tex.SetTextSize(0.05);
458 tex.DrawLatex(0.2,0.85,str);
459 sprintf(str,"#eta = %3.1f-%3.1f",m_etamin+m_deta*i_eta_sel, m_etamin+m_deta*(i_eta_sel+1));
460 tex.DrawLatex(0.2,0.78,str);
461 } // i_area
462 c1ps->Print(sfname.c_str());
463 } // i_eta
464 } // i_ener
465
466 // engRecoOverTruth
467 int a_color[] = {kRed, kBlue, kGreen};
468 for(int i_pdg=0; i_pdg<m_npdg; i_pdg++){
469 for(int i_ener=0; i_ener<int(sizeof(a_enerbin_sel)/sizeof(int)); i_ener++){
470 int i_ener_sel = a_enerbin_sel[i_ener];
471
472 char cname[1024];
473 sprintf(cname,"c2ps_engRecoOverTruth_%d_%d",i_pdg, i_ener);
474 TCanvas *c1 = new TCanvas(cname, cname, cc_xx, cc_yy);
475 c1->cd();
476 sprintf(cname,"c2ps_engRecoOverTruth_%d_%d_pad1",i_pdg, i_ener);
477 TPad *c1_pad1 = new TPad(cname, cname, 0.0, 0.9, 1.0, 1.0); c1_pad1->Draw();
478 sprintf(cname,"c2ps_engRecoOverTruth_%d_%d_pad2",i_pdg, i_ener);
479 TPad *c1_pad2 = new TPad(cname, cname, 0.0, 0.0, 1.0, 0.9); c1_pad2->Draw();
480
481 c1_pad1->cd();
482 sprintf(str,"pdg:%d E = %3.1f-%3.1f GeV",i_pdg,
483 pow(10,m_logenermin+m_dlogener*i_ener_sel)*1e-3,
484 pow(10,m_logenermin+m_dlogener*(i_ener_sel+1))*1e-3);
485 tex.SetTextSize(0.12);
486 tex.DrawLatex(0.1, 0.5, str);
487 c1_pad2->cd();
488 c1_pad2->Divide(3,3);
489 for(int i_area=0; i_area<int(sizeof(a_area_sel)/sizeof(int)); i_area++){
490 int i_area_sel = a_area_sel[i_area];
491 c1_pad2->cd(i_area+1); gPad->SetGrid();
492 for(int i_reco=0; i_reco<m_nrecobin; i_reco++){
493 TProfile *hp = m_hp_engRecoOverTruth_vs_eta[i_pdg][i_area_sel][i_reco][i_ener_sel];
494 hp->SetLineColor(a_color[i_reco]);
495 hp->SetMinimum(0.5);
496 hp->SetMaximum(1.2);
497 hp->GetXaxis()->SetTitle("mc_eta");
498 hp->GetYaxis()->SetTitle("engReco / engTruth");
499 hp->SetTitle(a_title[i_area_sel]);
500 if(i_reco==0) {
501 hp->Draw();
502 }else{
503 hp->Draw("same");
504 }
505 } // i_reco
506 }
507 std::cout << "adding " << i_pdg << " " << i_ener << std::endl;
508 c1->Print(sfname.c_str());
509
510 sprintf(cname,"c2ps_engResolution_%d_%d",i_pdg, i_ener);
511 TCanvas *c2 = new TCanvas(cname, cname, cc_xx, cc_yy);
512 c2->cd();
513 sprintf(cname,"c2ps_engResolution_%d_%d_pad1",i_pdg, i_ener);
514 TPad *c2_pad1 = new TPad(cname, cname, 0.0, 0.9, 1.0, 1.0); c2_pad1->Draw();
515 sprintf(cname,"c2ps_engResolution_%d_%d_pad2",i_pdg, i_ener);
516 TPad *c2_pad2 = new TPad(cname, cname, 0.0, 0.0, 1.0, 0.9); c2_pad2->Draw();
517
518 c2_pad1->cd();
519 sprintf(str,"pdg:%d E = %3.1f-%3.1f GeV (resolution)",i_pdg,
520 pow(10,m_logenermin+m_dlogener*i_ener_sel)*1e-3,
521 pow(10,m_logenermin+m_dlogener*(i_ener_sel+1))*1e-3);
522 tex.SetTextSize(0.12);
523 tex.DrawLatex(0.1, 0.5, str);
524 c2_pad2->cd();
525 c2_pad2->Divide(3,3);
526 TH1F *h1ref = new TH1F("h1ref","h1ref",100, m_etamin, m_etamax);
527 h1ref->SetMaximum(0.25);
528 h1ref->SetMinimum(0.0);
529 h1ref->GetXaxis()->SetTitle("#eta");
530 h1ref->GetYaxis()->SetTitle("#sigma_{E}/E");
531 h1ref->SetStats(0);
532 for(int i_area=0; i_area<int(sizeof(a_area_sel)/sizeof(int)); i_area++){
533 int i_area_sel = a_area_sel[i_area];
534 c2_pad2->cd(i_area+1); gPad->SetGrid();
535 h1ref->SetTitle(a_title[i_area_sel]);
536 h1ref->DrawCopy();
537
538 for(int i_reco=0; i_reco<m_nrecobin; i_reco++){
539 TGraphErrors *gr = new TGraphErrors(m_netabin);
540 for(int i_eta=0; i_eta<m_netabin; i_eta++){
541 std::cout << " i_pdg:" << i_pdg
542 << " i_area_sel:" << i_area_sel
543 << " i_reco:" << i_reco
544 << " i_ener_sel:" << i_ener_sel
545 << " i_eta:" << i_eta
546 << std::endl;
547 TH1F *h1 = m_engRecSpect[i_pdg][i_area_sel][i_reco][i_ener_sel][i_eta];
548 if(h1 == nullptr) {
549 std::cout << " panic, h1==0" << std::endl;
550 continue;
551 }
552 if(h1->GetMean() == 0.0) {
553 std::cout << " Resolution() -> Warning! Skipping point, no energy i_ener:" << i_ener_sel
554 << " i_eta:" << i_eta << " " << (m_etamin + (i_eta+0.5)*m_deta)
555 << " mean2:" << h1->GetMean()
556 << " nent2:" << h1->GetEntries() << std::endl;
557 continue;
558 }
559
560 if(m_doTailRejection) {
561 float mean3(0), rms3(0);
562 GetRmsWithoutTails(h1, mean3, rms3);
563 gr->SetPoint(i_eta, m_etamin + (i_eta+0.5)*m_deta, rms3/mean3);
564 }else{
565 gr->SetPoint(i_eta, m_etamin + (i_eta+0.5)*m_deta, h1->GetRMS()/h1->GetMean());
566 }
567 } // i_eta
568 gr->SetLineColor(a_color[i_reco]);
569 gr->SetMarkerColor(a_color[i_reco]);
570 gr->SetMarkerSize(1.0);
571 gr->Draw("pl");
572 } // i_reco
573 } // i_area
574 c2->Print(sfname.c_str());
575
576 } // i_ener
577 } // i_pdg
578
579 sfname = sreport;
580 sfname += "]";
581 ctmp->Print(sfname.c_str());
582}
583
584
585/* ****************************************************************************
586Toy analog of CaloLCDeadMaterialTool.cxx
587**************************************************************************** */
588void CaloHadDMCoeffCheck::getDmReco(std::vector<std::vector<double > > &engDmReco) // [ncls][narea]
589{
590 engDmReco.clear();
591 engDmReco.resize(m_data->m_ncls);
592 for(int i_cls=0; i_cls<m_data->m_ncls; i_cls++){
593 engDmReco[i_cls].resize(m_HadDMCoeff->getSizeAreaSet(), 0.0);
594 }
595 for(int i_cls=0; i_cls<m_data->m_ncls; i_cls++){ // loop over clusters
596 if( (*m_data->m_cls_ener_unw)[i_cls] < m_energyMin) continue;
597 float clusEnerOrig = (*m_data->m_cls_ener_unw)[i_cls];
598 for(int i_dma=0; i_dma<m_data->m_narea; i_dma++){ // loop over DM areas
599 const CaloLocalHadCoeff::LocalHadArea *area = m_HadDMCoeff->getArea(i_dma);
600
601 float clusEner = (*m_data->m_cls_ener_unw)[i_cls];
602 float emax = pow(10,area->getDimension(CaloLocalHadCoeffHelper::DIM_ENER)->getXmax());
603 if( clusEner > emax) clusEner = emax-0.001;
604 (*m_data->m_cls_ener_unw)[i_cls] = clusEner;
605
606 std::vector<float> vars;
607 m_data->PackClusterVars(i_cls, vars);
608
609 (*m_data->m_cls_ener_unw)[i_cls] = clusEnerOrig;
610
611 const CaloLocalHadCoeff::LocalHadCoeff *pars = m_HadDMCoeff->getCoeff(i_dma, vars);
612 if(pars==nullptr) continue;
613 double eprep = (*m_data->m_cls_eprep)[i_cls][i_dma];
614 double engDmRec = 0.0;
615 if(area->getType() == CaloLocalHadDefs::AREA_DMFIT && area->getNpars() == 3) {
616 if(eprep > 0.0) engDmRec = (*pars)[0] + (*pars)[1]*pow(eprep, (*pars)[2]);
617 }else if(area->getType() == CaloLocalHadDefs::AREA_DMLOOKUP && area->getNpars() == 3){
618 if( (*pars)[1] > 40) {
619 //double isol = (*m_data->m_cls_isol)[i_cls];
620 double isol = 1.0;
621 engDmRec = isol*(*m_data->m_cls_ener_unw)[i_cls]*((*pars)[0] - 1.0 );
622 }
623 }else if(area->getType() == CaloLocalHadDefs::AREA_DMSMPW && area->getNpars() == CaloSampling::Unknown+1){
624 double ecalonew = 0.0;
625 double ecaloold = 0.0;
626 for(int i_smp=0; i_smp<CaloSampling::Unknown; i_smp++){
627 float smpener = m_data->m_cls_smpener_unw ? (*m_data->m_cls_smpener_unw)[i_cls][i_smp] : 0.;
628 ecaloold += smpener;
629 ecalonew += smpener * (*pars)[i_smp];
630 }
631 ecalonew += (*pars)[CaloSampling::Unknown]; // const
632 engDmRec = ecalonew - ecaloold;
633 }else{
634 std::cout << "Panic! Unknown DM area type" << std::endl;
635 }
636// double edmWrong = 0.0;
637// if(i_dma == sDM_EMB0_EMB1) {
638// edmWrong = (*m_data->m_cls_smpener_unw)[i_cls][CaloSampling::PreSamplerB];
639// }else if(i_dma == sDM_SCN) {
640// edmWrong = (*m_data->m_cls_smpener_unw)[i_cls][CaloSampling::TileGap3];
641// }else if(i_dma == sDM_EME0_EME12) {
642// edmWrong = (*m_data->m_cls_smpener_unw)[i_cls][CaloSampling::PreSamplerE];
643// }else{
644// edmWrong = 0.0;
645// }
646// engDmRec -= edmWrong;
647 if(engDmRec >0.0) engDmReco[i_cls][i_dma] = engDmRec;
648 } // i_dma
649 } // i_cls
650}
651
652
653
654/* **************************************************************************
655
656************************************************************************** */
657int CaloHadDMCoeffCheck::GetRmsWithoutTails(const TH1F *pH, float &mean, float &rms, float ths)
658{
659 mean = pH->GetMean();
660 rms = pH->GetRMS();
661 float lim1 = mean - ths*rms;
662 float lim2 = mean + ths*rms;
663
664 //double sum(0);
665 double d_aver(0);
666 double d_rms(0);
667 double d_sw(0);
668 for(int i=1; i<=(int)pH->GetNbinsX(); i++){
669 if( pH->GetBinCenter(i)>=lim1 && pH->GetBinCenter(i)<= lim2 ){
670 float xx = pH->GetBinCenter(i);
671 float w = pH->GetBinContent(i);
672 if(w == 0) continue;
673 //sum += xx;
674 d_rms = (d_sw/(d_sw+w))*(d_rms+(w/(d_sw+w))*(xx-d_aver)*(xx-d_aver));
675 d_aver = d_aver+(xx-d_aver)*w/(d_sw+w);
676 d_sw += w;
677 }
678 }
679 if(d_rms==0) {
680 d_rms = pH->GetRMS();
681 }else{
682 d_rms = sqrt(d_rms);
683 }
684
685 if(d_aver==0) {
686 d_aver = pH->GetMean();
687 }
688
689 mean = d_aver;
690 rms = d_rms;
691 return 0;
692}
M_PI
#define M_PI
Definition ActiveFraction.h:11
area
double area(double R)
Definition ConvertStaveServices.cxx:42
gr
#define gr
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
h
Header file for AthHistogramAlgorithm.
CaloHadDMCoeffCheck::m_HadDMHelper
CaloLocalHadCoeffHelper * m_HadDMHelper
Definition CaloHadDMCoeffCheck.h:67
CaloHadDMCoeffCheck::process
int process(CaloHadDMCoeffData *myData, CaloLocalHadCoeff *myHadDMCoeff, bool isSingleParticle=true, bool tbflag=false)
Definition CaloHadDMCoeffCheck.cxx:80
CaloHadDMCoeffCheck::m_h2_etrue_vs_ereco
std::vector< std::vector< std::vector< TH2F * > > > m_h2_etrue_vs_ereco
Definition CaloHadDMCoeffCheck.h:88
CaloHadDMCoeffCheck::getDmReco
void getDmReco(std::vector< std::vector< double > > &engDmReco)
Definition CaloHadDMCoeffCheck.cxx:588
CaloHadDMCoeffCheck::GetRmsWithoutTails
int GetRmsWithoutTails(const TH1F *h1, float &mean, float &rms, float ths=3.0)
Definition CaloHadDMCoeffCheck.cxx:657
CaloHadDMCoeffCheck::m_data
CaloHadDMCoeffData * m_data
Definition CaloHadDMCoeffCheck.h:66
CaloHadDMCoeffCheck::m_engRecSpect
std::vector< std::vector< std::vector< std::vector< std::vector< TH1F * > > > > > m_engRecSpect
Definition CaloHadDMCoeffCheck.h:91
CaloHadDMCoeffCheck::m_HadDMCoeff
CaloLocalHadCoeff * m_HadDMCoeff
Definition CaloHadDMCoeffCheck.h:68
CaloHadDMCoeffCheck::m_hp_engRecoOverTruth_vs_eta
std::vector< std::vector< std::vector< std::vector< TProfile * > > > > m_hp_engRecoOverTruth_vs_eta
Definition CaloHadDMCoeffCheck.h:90
CaloHadDMCoeffCheck::m_hp_etrue_vs_ereco
std::vector< std::vector< std::vector< TProfile * > > > m_hp_etrue_vs_ereco
Definition CaloHadDMCoeffCheck.h:89
CaloHadDMCoeffCheck::make_report
void make_report(std::string &sfname)
Definition CaloHadDMCoeffCheck.cxx:378
CaloHadDMCoeffData
Data to read from special DeadMaterialTree.
Definition CaloHadDMCoeffData.h:30
CaloLocalHadCoeff::LocalHadArea
Definition of correction area.
Definition CaloLocalHadCoeff.h:145
CaloLocalHadCoeff
Hold binned correction data for local hadronic calibration procedure.
Definition CaloLocalHadCoeff.h:41
CaloLocalHadCoeff::LocalHadCoeff
std::vector< float > LocalHadCoeff
Correction parameters for one general bin.
Definition CaloLocalHadCoeff.h:220
GetLCSinglePionsPerf::angle_mollier_factor
static double angle_mollier_factor(double x)
Definition GetLCSinglePionsPerf.cxx:1296
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition dependence.cxx:254
add
bool add(const std::string &hname, TKey *tobj)
Definition fastadd.cxx:55
r
int r
Definition globals.cxx:22
Generated on for ATLAS Offline Software by doxygen 1.14.0