ATLAS Offline Software
Loading...
Searching...
No Matches
ConfAnalysis.cxx
Go to the documentation of this file.
1
9
10
11#include "ConfAnalysis.h"
12
13
14#include "ConfVtxAnalysis.h"
15
19
20#include <fstream>
21
22#include "BinConfig.h"
23
24#include "TF1.h"
25#include "TMath.h"
26
27#include "globals.h"
28
29bool PRINT_BRESIDUALS = false;
30
31std::ofstream dumpfile("dumpfile.log");
32
33void Normalise(TH1* h) {
34
35 for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
36
37 double del = h->GetBinLowEdge(i+1)-h->GetBinLowEdge(i);
38
39 h->SetBinContent( i, h->GetBinContent(i)/del );
40 h->SetBinError( i, h->GetBinError(i)/del );
41 }
42
43}
44
45
46
47
49
55
56
58
60 if ( !m_initialiseFirstEvent ) {
61 std::cout << "ConfAnalysis::initialise() " << std::endl;
63 }
64}
65
67
68 if ( m_initialised ) return;
69
70 m_initialised = true;
71
72 // std::cout << "ConfAnalysis::initialiseInternal() " << name() << std::endl;
73
74 BinConfig& binConfig = g_binConfig;
75
76 if ( name().find("_e")!=std::string::npos ) binConfig = electronBinConfig;
77 else if ( name().find("_mu")!=std::string::npos ) binConfig = muonBinConfig;
78 else if ( name().find("_tau")!=std::string::npos ) binConfig = tauBinConfig;
79 else if ( name().find("_b")!=std::string::npos ) binConfig = bjetBinConfig;
80 else if ( name().find("cosmic")!=std::string::npos ) binConfig = cosmicBinConfig;
81
82
83 //+++ pT ranges
84 // double tmp_maxPt = 50000.;
85 // double tmp_absResPt = 0.0005;
86 //double tmp_maxPt = 50.;
87 double tmp_absResPt = 0.5;
88
89 const int pTResBins = int(100*binConfig.ptres_NScale);
90
91
92 //+++ Eta ranges
95 double tmp_maxEta = 5.;
96 double tmp_absResEta = 0.04; // 0.0005;
97
98 //+++ Phi ranges
99 double tmp_maxPhi = 3.142;
100 double tmp_absResPhi = 0.02; // 0.0001;
101
102
103 // std::cout << "ConfAnalysis::initialise() " << name() << " config: " << binConfig << std::endl;
104
107 int etaBins = int(50*binConfig.eta_NScale);
108 const int etaResBins = int(600*binConfig.eta_NScale);
109
110 const int phiBins = int(30*binConfig.phi_NScale);
111 const int phiResBins = int(100*binConfig.phires_NScale);
112
113 const int zBins = int(150*binConfig.z0_NScale);
114 const double zMax = binConfig.z0Max;
115
116 const int zresBins = 100;
117 const double zresMax = 10;
118
119 const int d0Bins = int(100*binConfig.d0_NScale);
120 const double d0Max = binConfig.d0Max;
121
122 const int d0resBins = 100;
123 const double d0resMax = 5;
124
125 // beamspot corrected position
126 const int a0Bins = int(300*binConfig.a0_NScale);
127 const double a0Max = binConfig.a0Max;
128
129 const int a0resBins = 100;
130 const double a0resMax = 5;
131
132 //+++ Book histograms
133
134 // calculate a logarithmic binning in pt
135
136 int Npt = 0;
137 double pt_a = 1;
138 double pt_b = 1;
139
140 // Npt = int(40*binConfig.pt_NScale);
141 // pt_a = 3.5;
142 Npt = int(45*binConfig.pt_NScale);
143 pt_a = 4;
144 pt_b = 2;
145 // etaBins = 12;
146 // }
147 // else {
148 // Npt = 40;
149 // pt_a = 4;
150 // pt_b = 1;
151 // }
152
153
154
155 const int ptnbins = Npt;
156 std::vector<double> ptbinlimsv(ptnbins+1);
157 double* ptbinlims = &ptbinlimsv[0];
158 // for ( int i=0 ; i<=ptnbins ; i++ ) { ptbinlims[i] = std::pow(10, 2.0*i/ptnbins+2)/1000; }
159 // for ( int i=0 ; i<=ptnbins ; i++ ) { ptbinlims[i] = std::pow(10, 2.3*i/ptnbins+2); }
160 for ( int i=0 ; i<=ptnbins ; i++ ) ptbinlims[i] = std::pow(10, pt_a*i/ptnbins+pt_b)/1000;
161
162
163 // ADDED BY JK - FOR SIGNED PT PLOTS
164 //-----
165 const int ptnbins2 = (2*ptnbins);
166 // std::cout << "ptnbins2 = " << ptnbins2 << std::endl;
167 std::vector<double> ptbinlims2v(ptnbins2 + 1);
168 double* ptbinlims2 = &ptbinlims2v[0];
169 // std::cout << "ptbinlims2v.size() = " << ptbinlims2v.size() << std::endl;
170 int ptnbin_counter = 0;
171 for ( int i=ptnbins; i>0 ; i-- ) {
172 ptbinlims2[ptnbin_counter] = std::pow(10, pt_a*i/ptnbins+pt_b)/(-2000);
173 // std::cout << "ptbinlims[" << i << "] = " << ptbinlims[i] << " , so ptbinlims2[" << ptnbin_counter << "] = " << ptbinlims2[ptnbin_counter] << std::endl;
174 ptnbin_counter++;
175 }
176
177 for ( int i=0 ; i<ptnbins+1 ; i++ ) {
178 ptbinlims2[ptnbin_counter] = std::pow(10, pt_a*i/ptnbins+pt_b)/2000;
179 // std::cout << "ptbinlims[" << i << "] = " << ptbinlims[i] << " , so ptbinlims2[" << ptnbin_counter << "] = " << ptbinlims2[ptnbin_counter] << std::endl;
180 ptnbin_counter++;
181 }
182 //-----
183
184 const int iptnbins = 20;
185 const double minmaxipt=0.5;
186 std::vector<double> iptbinlimsv(iptnbins+1);
187 double* iptbinlims = &iptbinlimsv[0];
188 for ( int i=0 ; i<=iptnbins ; i++ ) {
189 iptbinlims[i] = -minmaxipt+i*minmaxipt*2./iptnbins;
190 }
191
192
193
194 TDirectory* dir = gDirectory;
195
196 // std::cout << "ConfAnalysis::initialize() Directory " << gDirectory->GetName() << " " << name() << std::endl;
197
198 m_dir = new TIDDirectory(name());
199 m_dir->push();
200
201 // std::cout << "ConfAnalysis::initialize() Directory " << gDirectory->GetName() << " " << name() << std::endl;
202
203 if ( !m_initialiseFirstEvent ) std::cout << "ConfAnalysis::initialize() Directory " << name() << std::endl;
204
205 if ( name() != gDirectory->GetName() ) std::cerr << "ConfAnalysis::initialize() Directory: problem with directory " << gDirectory->GetName() << " " << name() << std::endl;
206
207 // TIDDirectory d("histos");
208 // d.push();
209
210 // std::cout << "ConfAnalysis::initialize() Directory " << gDirectory->GetName() << " package directory, " << name() << std::endl;
211
212
213 m_res.push_back( m_rnpix_eta = new Resplot( "npix_eta", 2*etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
214 m_res.push_back( m_rnsct_eta = new Resplot( "nsct_eta", 2*etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
215 m_res.push_back( m_rntrt_eta = new Resplot( "ntrt_eta", 2*etaBins, -tmp_maxEta, tmp_maxEta, 100, -0.5, 99.5 ) );
216 m_res.push_back( m_rnsihit_eta= new Resplot( "nsihit_eta",etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
217
218 m_res.push_back( m_rnpix_lb = new Resplot( "npix_lb", 250, 0, 2500, 22, -0.5, 21.5 ) );
219
220 m_res.push_back( m_rnpix_phi = new Resplot( "npix_phi", etaBins, -M_PI, M_PI, 22, -0.5, 21.5 ) );
221 m_res.push_back( m_rnsct_phi = new Resplot( "nsct_phi", etaBins, -M_PI, M_PI, 22, -0.5, 21.5 ) );
222 m_res.push_back( m_rntrt_phi = new Resplot( "ntrt_phi", etaBins, -M_PI, M_PI, 100, -0.5, 99.5 ) );
223
224 m_res.push_back( m_rnpix_pt = new Resplot( "npix_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
225 m_res.push_back( m_rnsct_pt = new Resplot( "nsct_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
226 m_res.push_back( m_rntrt_pt = new Resplot( "ntrt_pt", ptnbins, ptbinlims, 100, -0.5, 99.5 ) );
227
228 m_res.push_back( m_rnpixh_pt = new Resplot( "npixh_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
229 m_res.push_back( m_rnscth_pt = new Resplot( "nscth_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
230
231
232 m_res.push_back( m_rnpix_pt_bad = new Resplot( "npix_pt_bad", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
233 m_res.push_back( m_rnsct_pt_bad = new Resplot( "nsct_pt_bad", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
234 m_res.push_back( m_rntrt_pt_bad = new Resplot( "ntrt_pt_bad", ptnbins, ptbinlims, 100, -0.5, 99.5 ) );
235
236
237 m_res.push_back( m_rnpix_eta_rec = new Resplot( "npix_eta_rec", 2*etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
238 m_res.push_back( m_rnsct_eta_rec = new Resplot( "nsct_eta_rec", 2*etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
239 m_res.push_back( m_rntrt_eta_rec = new Resplot( "ntrt_eta_rec", 2*etaBins, -tmp_maxEta, tmp_maxEta, 100, -0.5, 99.5 ) );
240 m_res.push_back( m_rnsihit_eta_rec= new Resplot( "nsihit_eta_rec", etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
241
242 m_res.push_back( m_rnpix_phi_rec = new Resplot( "npix_phi_rec", etaBins, -M_PI, M_PI, 22, -0.5, 21.5 ) );
243 m_res.push_back( m_rnsct_phi_rec = new Resplot( "nsct_phi_rec", etaBins, -M_PI, M_PI, 22, -0.5, 21.5 ) );
244 m_res.push_back( m_rntrt_phi_rec = new Resplot( "ntrt_phi_rec", etaBins, -M_PI, M_PI, 100, -0.5, 99.5 ) );
245
246 m_res.push_back( m_rnpix_pt_rec = new Resplot( "npix_pt_rec", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
247 m_res.push_back( m_rnsct_pt_rec = new Resplot( "nsct_pt_rec", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
248 m_res.push_back( m_rntrt_pt_rec = new Resplot( "ntrt_pt_rec", ptnbins, ptbinlims, 100, -0.5, 99.5 ) );
249
250 m_res.push_back( m_rnpixh_pt_rec = new Resplot( "npixh_pt_rec", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
251 m_res.push_back( m_rnscth_pt_rec = new Resplot( "nscth_pt_rec", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
252
253
254 m_res.push_back( m_rnsct_vs_npix = new Resplot( "nsct_vs_npix", 12, -0.5, 11.5, 22, -0.5, 21.5 ) );
255 m_res.push_back( m_rnsct_vs_npix_rec = new Resplot( "nsct_vs_npix_rec", 12, -0.5, 11.5, 22, -0.5, 21.5 ) );
256
257 m_res.push_back( m_rChi2prob = new Resplot( "Chi2prob", ptnbins, ptbinlims, 20, 0, 1 ) );
258 m_res.push_back( m_rChi2 = new Resplot( "Chi2", ptnbins, ptbinlims, 200, 0, 100 ) );
259 m_res.push_back( m_rChi2dof = new Resplot( "Chi2dof", ptnbins, ptbinlims, 100, 0, 10 ) );
260
261 m_res.push_back( m_rChi2prob_bad = new Resplot( "Chi2prob_bad", ptnbins, ptbinlims, 20, 0, 1 ) );
262 m_res.push_back( m_rChi2_bad = new Resplot( "Chi2_bad", ptnbins, ptbinlims, 200, 0, 100 ) );
263 m_res.push_back( m_rChi2dof_bad = new Resplot( "Chi2dof_bad", ptnbins, ptbinlims, 100, 0, 10 ) );
264
265 m_res.push_back( m_rChi2prob_rec = new Resplot( "Chi2prob_rec", ptnbins, ptbinlims, 20, 0, 1 ) );
266 m_res.push_back( m_rChi2_rec = new Resplot( "Chi2_rec", ptnbins, ptbinlims, 200, 0, 100 ) );
267 m_res.push_back( m_rChi2dof_rec = new Resplot( "Chi2dof_rec", ptnbins, ptbinlims, 100, 0, 10 ) );
268
269 m_res.push_back( m_rChi2d_vs_Chi2d = new Resplot( "Chi2d_vs_Chi2d", 200, 0, 100, 200, 0, 100 ) );
270
271 m_res.push_back( m_rDChi2dof = new Resplot( "DChi2dof", 200, 0, 100, 200, -100, 100 ) );
272
274
275 double d0bins[40] = { -5.0, -4.0, -3.0, -2.5,
276 -2.0, -1.8, -1.6, -1.4, -1.2,
277 -1.05, -0.95, -0.85, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
278 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, 1.05,
279 1.2, 1.4, 1.6, 1.8, 2.0,
280 2.5, 3.0, 4.0, 5.0 };
281
282
283
284 m_res.push_back( m_rnpix_d0 = new Resplot( "npix_d0", 39, d0bins, 22, -0.5, 21.5 ) );
285 m_res.push_back( m_rnsct_d0 = new Resplot( "nsct_d0", 39, d0bins, 22, -0.5, 21.5 ) );
286 m_res.push_back( m_rntrt_d0 = new Resplot( "ntrt_d0", 39, d0bins, 22, -0.5, 21.5 ) );
287
288 m_res.push_back( m_rnpixh_d0 = new Resplot( "npixh_d0", 39, d0bins, 22, -0.5, 21.5 ) );
289 m_res.push_back( m_rnscth_d0 = new Resplot( "nscth_d0", 39, d0bins, 22, -0.5, 21.5 ) );
290
291 m_res.push_back( m_rnsi_pt = new Resplot( "nsi_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
292 m_res.push_back( m_rnsih_pt = new Resplot( "nsih_pt", ptnbins, ptbinlims, 22, -0.5, 21.5 ) );
293
294 m_res.push_back( m_rnsi_d0 = new Resplot( "nsi_d0", 39, d0bins, 22, -0.5, 21.5 ) );
295 m_res.push_back( m_rnsih_d0 = new Resplot( "nsih_d0", 39, d0bins, 22, -0.5, 21.5 ) );
296
297 m_res.push_back( m_rnsi_eta = new Resplot( "nsi_eta", etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
298 m_res.push_back( m_rnsih_eta = new Resplot( "nsih_eta", etaBins, -tmp_maxEta, tmp_maxEta, 22, -0.5, 21.5 ) );
299
300 m_res.push_back( m_rnbl_d0 = new Resplot( "nbl_d0", 39, d0bins, 5, -0.5, 4.5 ) );
301 m_res.push_back( m_rnblh_d0 = new Resplot( "nblh_d0", 39, d0bins, 5, -0.5, 4.5 ) );
302
303
304 m_res.push_back( m_rnsct_lb = new Resplot( "nsct_lb", 250, 0, 2500, 22, -0.5, 21.5 ) );
305
306 m_res.push_back( m_rnpix_lb_rec = new Resplot( "npix_lb_rec", 250, 0, 2500, 22, -0.5, 21.5 ) );
307 m_res.push_back( m_rnsct_lb_rec = new Resplot( "nsct_lb_rec", 250, 0, 2500, 22, -0.5, 21.5 ) );
308
309 m_res.push_back( m_rnpix_d0_rec = new Resplot( "npix_d0_rec", 39, d0bins, 22, -0.5, 21.5 ) );
310 m_res.push_back( m_rnsct_d0_rec = new Resplot( "nsct_d0_rec", 39, d0bins, 22, -0.5, 21.5 ) );
311 m_res.push_back( m_rntrt_d0_rec = new Resplot( "ntrt_d0_rec", 39, d0bins, 22, -0.5, 21.5 ) );
312
313 // int Nptbins = 7;
314 // double _ptlims[8] = { 0, 500, 1000, 1500, 2000, 5000, 8000, 12000 };
315
316
317 // addHistogram( hchi2=new TH1F("chi2", "chi2", 100, 0, 20) );
318
319 // "reference" quantities
320 addHistogram( new TH1F( "pT", "pT", ptnbins, ptbinlims ) );
321 addHistogram( new TH1F( "eta", "eta", etaBins, -tmp_maxEta, tmp_maxEta ) );
322 addHistogram( new TH1F( "phi", "phi", phiBins, -tmp_maxPhi, tmp_maxPhi ) );
323 addHistogram( new TH1F( "z0", "z0", zBins, -zMax, zMax ) );
324 addHistogram( new TH1F( "d0", "d0", d0Bins, -d0Max, d0Max ) );
325 addHistogram( new TH1F( "a0", "a0", a0Bins, -a0Max, a0Max ) );
326
327 // error study histograms (reference)
328 addHistogram( new TH1F( "dpT", "dpT", 80, 0, 20 ) );
329 addHistogram( new TH1F( "deta", "deta", 50, 0, 1 ) );
330 addHistogram( new TH1F( "dphi", "dphi", 50, 0, 1 ) );
331 addHistogram( new TH1F( "dz0", "dz0", 100, 0, 2 ) );
332 addHistogram( new TH1F( "dd0", "dd0", 50, 0, 0.5 ) );
333 addHistogram( new TH1F( "da0", "da0", 50, 0, 0.5 ) );
334
335 addHistogram( new TH1F( "roi_deta", "roi_deta", 50, -1, 1 ) );
336 addHistogram( new TH1F( "roi_dphi", "roi_dphi", 50, -1, 1 ) );
337 addHistogram( new TH1F( "roi_dR", "roi_dR", 50, 0, 1 ) );
338
339 // tag and probe invariant mass histograms
340 if ( m_TnP_tool ) {
341 m_invmass = new TH1F( "invmass", "invariant mass;mass [GeV]", 320, 0, 200 );
342 m_invmassObj = new TH1F( "invmassObj", "invariant mass;mass [GeV]", 320, 0, 200 );
345 }
346
348 TH2F* heta_vs_pt = new TH2F( "eta_vs_pt", "p_{T} [GeV],#eta", ptnbins, ptbinlims, 10, -3, 3 );
349 TH2F* hd0_vs_pt = new TH2F( "d0_vs_pt", "p_{T} [GeV],d_{0} [mm]", ptnbins, ptbinlims, 30, -300, 300 );
350
351
352 m_eff_eta_vs_pt = new Efficiency2D( heta_vs_pt, "eff_eta_vs_pt" );
353 m_eff_d0_vs_pt = new Efficiency2D( hd0_vs_pt, "eff_d0_vs_pt" );
354
355 delete heta_vs_pt;
356 delete hd0_vs_pt;
357
358
359 // efficiencies and purities
360 m_eff_pt = new Efficiency1D( find("pT"), "pT_eff" );
361 m_eff_pt->Hist()->GetXaxis()->SetTitle("P_{T} [GeV]");
362 m_eff_pt->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
363
364 m_eff_eta = new Efficiency1D( find("eta"), "eta_eff" );
365 m_eff_eta->Hist()->GetXaxis()->SetTitle("#eta");
366 m_eff_eta->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
367
368 m_eff_phi = new Efficiency1D( find("phi"), "phi_eff" );
369 m_eff_phi->Hist()->GetXaxis()->SetTitle("#phi");
370 m_eff_phi->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
371
372 m_eff_z0 = new Efficiency1D( find("z0"), "z0_eff" );
373 m_eff_z0->Hist()->GetXaxis()->SetTitle("z0");
374 m_eff_z0->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
375
376 m_eff_d0 = new Efficiency1D( find("d0"), "d0_eff" );
377 m_eff_d0->Hist()->GetXaxis()->SetTitle("d0");
378 m_eff_d0->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
379
380 m_eff_a0 = new Efficiency1D( find("a0"), "a0_eff" );
381 m_eff_a0->Hist()->GetXaxis()->SetTitle("a0");
382 m_eff_a0->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
383
384 m_eff_ptm = new Efficiency1D( find("pT"), "pTm_eff" );
385 m_eff_ptm->Hist()->GetXaxis()->SetTitle("Negative P_{T} [GeV]");
386 m_eff_ptm->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
387
388 m_eff_ptp = new Efficiency1D( find("pT"), "pTp_eff" );
389 m_eff_ptp->Hist()->GetXaxis()->SetTitle("Positive P_{T} [GeV]");
390 m_eff_ptp->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
391
392 m_eff_roi_deta = new Efficiency1D( find("roi_deta"), "roi_deta_eff" );
393 m_eff_roi_deta->Hist()->GetXaxis()->SetTitle("RoI #Delta#eta");
394 m_eff_roi_deta->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
395
396 m_eff_roi_dphi = new Efficiency1D( find("roi_dphi"), "roi_dphi_eff" );
397 m_eff_roi_dphi->Hist()->GetXaxis()->SetTitle("RoI #Delta#phi");
398 m_eff_roi_dphi->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
399
400 m_eff_roi_dR = new Efficiency1D( find("roi_dR"), "roi_dR_eff" );
401 m_eff_roi_dR->Hist()->GetXaxis()->SetTitle("RoI #Delta R");
402 m_eff_roi_dR->Hist()->GetYaxis()->SetTitle("Efficiency [%]");
403
404 // addHistogram ( m_hDeltaR = new TH1F("DeltaR", "DeltaR", 100, 0, 0.1 ) );
405 addHistogram ( m_hDeltaR = new TH1F("DeltaR", "DeltaR", 100, 0, 0.2 ) );
406
407 m_purity_pt = new Efficiency1D( find("pT"), "pT_pur" );
408 m_purity_eta = new Efficiency1D( find("eta"), "eta_pur" );
409 m_purity_phi = new Efficiency1D( find("phi"), "phi_pur" );
410 m_purity_z0 = new Efficiency1D( find("z0"), "z0_pur" );
411 m_purity_d0 = new Efficiency1D( find("d0"), "d0_pur" );
412 m_purity_a0 = new Efficiency1D( find("a0"), "a0_pur" );
413
414 // "test" quantities
415 addHistogram( new TH1F( "pT_rec", "pT_rec", ptnbins, ptbinlims ) );
416 addHistogram( new TH1F( "eta_rec", "eta_rec", etaBins, -tmp_maxEta, tmp_maxEta ) );
417 addHistogram( new TH1F( "phi_rec", "phi_rec", phiBins, -tmp_maxPhi, tmp_maxPhi ) );
418 addHistogram( new TH1F( "z0_rec", "z0_rec", zBins, -zMax, zMax ) );
419 addHistogram( new TH1F( "d0_rec", "d0_rec", d0Bins, -d0Max, d0Max ) );
420 addHistogram( new TH1F( "a0_rec", "a0_rec", a0Bins, -a0Max, a0Max ) );
421
422 addHistogram2D( new TH2F( "eta_phi_rec", "eta_phi_rec", (tmp_maxEta+1)*30 , -tmp_maxEta-1, tmp_maxEta+1, (tmp_maxPhi+1)*30, -tmp_maxPhi-1, tmp_maxPhi+1 ) );
423 addHistogram2D( new TH2F( "phi_d0_rec", "phi_d0_rec", (2*tmp_maxPhi+2)*15, -tmp_maxPhi-1, tmp_maxPhi+1 ,d0Bins+20, -d0Max+7, d0Max-7 ));
424
425 // error study histograms (test)
426 addHistogram( new TH1F( "dpT_rec", "dpT_rec", 80, 0, 20.00 ) );
427 addHistogram( new TH1F( "deta_rec", "deta_rec", 50, 0, 0.02 ) );
428 addHistogram( new TH1F( "dphi_rec", "dphi_rec", 50, 0, 0.02 ) );
429 addHistogram( new TH1F( "dz0_rec", "dz0_rec", 100, 0, 1.5 ) );
430 addHistogram( new TH1F( "dd0_rec", "dd0_rec", 50, 0, 0.5 ) );
431 addHistogram( new TH1F( "da0_rec", "da0_rec", 50, 0, 0.5 ) );
432
433
434
435 // error study histograms (pull test-ref)
436 addHistogram( new TH1F("pT_pull", "pT_pull", 100, -10, 10) );
437 addHistogram( new TH1F("eta_pull", "eta_pull", 100, -10, 10) );
438 addHistogram( new TH1F("phi_pull", "phi_pull", 100, -10, 10) );
439 addHistogram( new TH1F("z0_pull", "z0_pull", 100, -10, 10) );
440 addHistogram( new TH1F("d0_pull", "d0_pull", 100, -10, 10) );
441 addHistogram( new TH1F("a0_pull", "a0_pull", 100, -10, 10) );
442
443 // error study histograms (pull test-ref) - SIMPLE VERSION
444 addHistogram( new TH1F("pT_pull_simple", "pT_pull_simple", 100, -10, 10) );
445 addHistogram( new TH1F("eta_pull_simple", "eta_pull_simple", 100, -10, 10) );
446 addHistogram( new TH1F("phi_pull_simple", "phi_pull_simple", 100, -10, 10) );
447 addHistogram( new TH1F("z0_pull_simple", "z0_pull_simple", 100, -10, 10) );
448 addHistogram( new TH1F("d0_pull_simple", "d0_pull_simple", 100, -10, 10) );
449 addHistogram( new TH1F("a0_pull_simple", "a0_pull_simple", 100, -10, 10) );
450
451
452 // resolutions
453 TH1F* pT_res = new TH1F( "pT_res", "pT_res", 4*pTResBins, -0.1, 0.1 );
454 pT_res->GetXaxis()->SetTitle("#Delta P_{T} [GeV]");
455 pT_res->GetYaxis()->SetTitle("Entries");
456 addHistogram( pT_res );
457
458
459 TH1F* spT_res = new TH1F( "spT_res", "spT_res", 4*pTResBins, -0.1, 0.1 );
460 spT_res->GetXaxis()->SetTitle("#Delta sP_{T} [GeV]");
461 spT_res->GetYaxis()->SetTitle("Entries");
462 addHistogram( spT_res );
463
464
465 TH1F* ipT_res = new TH1F( "ipT_res", "ipT_res", 4*pTResBins, -0.4, 0.4 );
466 ipT_res->GetXaxis()->SetTitle("#Delta 1/P_{T} [GeV^{-1}]");
467 ipT_res->GetYaxis()->SetTitle("Entries");
468 addHistogram( ipT_res );
469
470 TH1F* eta_res = new TH1F("eta_res", "eta_res", etaResBins, -2*tmp_absResEta, 2*tmp_absResEta );
471 eta_res->GetXaxis()->SetTitle("#Delta #eta");
472 eta_res->GetYaxis()->SetTitle("Entries");
473 addHistogram( eta_res );
474 addHistogram( new TH1F("etai_res", "etai_res", 1000, -0.04, 0.04 ) );
475
476 // TH1F* phi_res = new TH1F( "phi_res", "phi_res;#Delta #phi;Entries", 2*phiResBins, -2*tmp_absResPhi, 2*tmp_absResPhi );
477 // phi_res->GetXaxis()->SetTitle("#Delta #phi");
478 // phi_res->GetYaxis()->SetTitle("Entries");
479 // addHistogram( phi_res );
480
481 addHistogram( new TH1F( "phi_res", "phi_res;#Delta #phi;Entries", 2*phiResBins, -2*tmp_absResPhi, 2*tmp_absResPhi ) );
482 addHistogram( new TH1F( "z0_res", "z0_res;#Deltaz_{0};Entries", 16*zresBins, -8*zresMax, 8*zresMax ) );
483 addHistogram( new TH1F( "d0_res", "d0_res;#Deltad_{0};Entries", 4*d0resBins, -0.2*d0resMax, 0.2*d0resMax ) );
484 addHistogram( new TH1F( "a0_res", "a0_res;#Deltaa_{0};Entries", 4*a0resBins, -0.2*a0resMax, 0.2*a0resMax ) );
485
486
487 addHistogram( new TH1F( "dphi_res", "dphi_res;#Delta #phi;Entries", 2*phiResBins, -0.2*tmp_absResPhi, 0.2*tmp_absResPhi ) );
488 addHistogram( new TH1F( "dz0_res", "dz0_res;#Deltaz_{0};Entries", 8*zresBins, -0.8*zresMax, 0.8*zresMax ) );
489 addHistogram( new TH1F( "dd0_res", "dd0_res;#Deltad_{0};Entries", 4*d0resBins, -0.05*d0resMax, 0.05*d0resMax ) );
490 addHistogram( new TH1F( "da0_res", "da0_res;#Deltaa_{0};Entries", 4*a0resBins, -0.05*a0resMax, 0.05*a0resMax ) );
491
492
493 // std::cout << "booking resplots" << std::endl;
494
495 int factor = 10;
496
497 int wfactor = 2;
498 int zfactor = 16;
499
500
501 m_rDd0res.push_back( new Resplot("rDd0_vs_pt", ptnbins, ptbinlims, 1200, -0.1, 0.1 ) );
502 m_rDd0res.push_back( new Resplot("rDd0_vs_eta", 5*etaBins, -tmp_maxEta, tmp_maxEta, 1200, -0.1, 0.1 ) );
503 m_rDd0res.push_back( new Resplot("rDd0_vs_zed", zBins, -zMax, zMax, 1200, -0.1, 0.1 ) );
504 m_rDd0res.push_back( new Resplot("rDd0_vs_d0", 20, -1.5, 1.5, 1200, -0.1, 0.1 ) );
505 // m_rDd0res.push_back( new Resplot("rDd0_vs_phi", int(2*M_PI/0.1), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02), 1200, -0.01, 0.05 ) );
506 m_rDd0res.push_back( new Resplot("rDd0_vs_phi", 128, -M_PI, M_PI, 1200, -0.1, 0.1 ) );
507
508
509 m_rDa0res.push_back( new Resplot("rDa0_vs_pt", ptnbins, ptbinlims, 100, -0.1, 0.1 ) );
510 m_rDa0res.push_back( new Resplot("rDa0_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 100, -0.1, 0.1 ) );
511 m_rDa0res.push_back( new Resplot("rDa0_vs_zed", 0.2*zBins, -zMax, zMax, 100, -0.1, 0.1 ) );
512 m_rDa0res.push_back( new Resplot("rDa0_vs_da0", 20, -1.5, 1.5, 100, -0.1, 0.1 ) );
513 m_rDa0res.push_back( new Resplot("rDa0_vs_phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02), 100, -0.1, 0.1 ) );
514 m_rDa0res.push_back( new Resplot("rDa0_rec_vs_phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02), 100, -0.1, 0.1 ) );
515
516 m_rDa0res.push_back( new Resplot("rda0_vs_phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02), 100, -0.1, 0.1 ) );
517 m_rDa0res.push_back( new Resplot("rda0_rec_vs_phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02), 100, -0.1, 0.1 ) );
518
519 // m_rDd0res.push_back( new Resplot("rDd0_vs_ipt", iptnbins, iptbinlims, 100, -0.1, 0.1 ) );
520 // m_rDz0res.push_back( new Resplot("rDz0_vs_ipt", iptnbins, iptbinlims, 100, -1, 1 ) );
521
522
523 m_rDz0res.push_back( new Resplot("rDz0_vs_pt", ptnbins, ptbinlims, 501, -1, 1 ) );
524 m_rDz0res.push_back( new Resplot("rDz0_vs_eta", 5*etaBins, -tmp_maxEta, tmp_maxEta, 500, -1, 1 ) );
525 m_rDz0res.push_back( new Resplot("rDz0_vs_zed", zBins, -zMax, zMax, 500, -1, 1 ) );
526 m_rDz0res.push_back( new Resplot("rDz0_vs_phi", 128, -M_PI, M_PI, 500, -1, 1 ) );
527
528
529
532 m_retares.push_back( new Resplot("reta_vs_ipt", iptnbins, iptbinlims, 2*etaResBins, -wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
533 m_rphires.push_back( new Resplot("rphi_vs_ipt", iptnbins, iptbinlims, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
534 m_rzedres.push_back( new Resplot("rzed_vs_ipt", iptnbins, iptbinlims, 8*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
535 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_ipt", iptnbins, iptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
536 m_riptres.push_back( new Resplot("ript_vs_ipt", iptnbins, iptbinlims, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
537 m_rptres.push_back( new Resplot("rpt_vs_ipt", iptnbins, iptbinlims, 8*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
538 m_rd0res.push_back( new Resplot("rd0_vs_ipt", iptnbins, iptbinlims, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
539
540
541 m_retares.push_back( new Resplot("reta_vs_pt", ptnbins, ptbinlims, 8*etaResBins, -wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
542 m_rphires.push_back( new Resplot("rphi_vs_pt", ptnbins, ptbinlims, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
543 m_rzedres.push_back( new Resplot("rzed_vs_pt", ptnbins, ptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
544 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_pt", ptnbins, ptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
545 m_riptres.push_back( new Resplot("ript_vs_pt", ptnbins, ptbinlims, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
546 m_rptres.push_back( new Resplot("rpt_vs_pt", ptnbins, ptbinlims, 8*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
547 m_rd0res.push_back( new Resplot("rd0_vs_pt", ptnbins, ptbinlims, factor*24*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
548
549 m_retares.push_back( new Resplot("reta_vs_ET", ptnbins, ptbinlims, 8*etaResBins, -wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
550 m_rphires.push_back( new Resplot("rphi_vs_ET", ptnbins, ptbinlims, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
551 m_rzedres.push_back( new Resplot("rzed_vs_ET", ptnbins, ptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
552 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_ET", ptnbins, ptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
553 m_riptres.push_back( new Resplot("ript_vs_ET", ptnbins, ptbinlims, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
554 m_rptres.push_back( new Resplot("rpt_vs_ET", ptnbins, ptbinlims, 8*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
555 m_rd0res.push_back( new Resplot("rd0_vs_ET", ptnbins, ptbinlims, factor*24*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
556
557
558 // m_retares.push_back( new Resplot("reta_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*etaResBins, -tmp_absResEta, tmp_absResEta ) );
559 m_retares.push_back( new Resplot("reta_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*etaResBins, -wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
560 m_rphires.push_back( new Resplot("rphi_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
561 m_rzedres.push_back( new Resplot("rzed_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 12*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
562
563 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 24*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
564 //m_rzedres.push_back( new Resplot("rzed_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
565 m_riptres.push_back( new Resplot("ript_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 16*pTResBins, -tmp_absResPt, tmp_absResPt ) );
566 m_rptres.push_back( new Resplot("rpt_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 8*pTResBins, -tmp_absResPt, tmp_absResPt ) );
567 m_rd0res.push_back( new Resplot("rd0_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, factor*24*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
568
569
570 // rphivsDd0res = new Resplot( "rphi_vs_Dd0", 10, 0, 0.1, int(2*M_PI/0.02), -0.2*int(M_PI/0.02), 0.2*int(M_PI/0.02) );
571 m_hphivsDd0res[0] = new TH1F( "hphi_vs_Dd0_0", "phi for Dd0<0.1", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
572 m_hphivsDd0res[1] = new TH1F( "hphi_vs_Dd0_1", "phi for Dd0>0.1", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
573 m_hphivsDd0res[2] = new TH1F( "hphi_vs_Dd0_2", "all phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
574
575 m_hphivsDa0res[0] = new TH1F( "hphi_vs_Da0_0", "phi for Da0<0.1", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
576 m_hphivsDa0res[1] = new TH1F( "hphi_vs_Da0_1", "phi for Da0>0.1", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
577 m_hphivsDa0res[2] = new TH1F( "hphi_vs_Da0_2", "all phi", int(2*M_PI/0.02), -0.02*int(M_PI/0.02), 0.02*int(M_PI/0.02) );
578
579
580 for ( unsigned ih=0 ; ih<3 ; ih++ ) {
581 m_hphivsDd0res[ih]->SetDirectory(0);
582 m_hphivsDa0res[ih]->SetDirectory(0);
583 }
584
585 m_retares.push_back( new Resplot("reta_vs_zed", 0.2*zBins, -zMax, zMax, 2*etaResBins, -tmp_absResEta, tmp_absResEta ) );
586 m_rphires.push_back( new Resplot("rphi_vs_zed", 0.2*zBins, -zMax, zMax, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
587 m_rzedres.push_back( new Resplot("rzed_vs_zed", 0.2*zBins, -zMax, zMax, 8*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
588 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_zed", 0.2*zBins, -zMax, zMax, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
589 //m_rzedres.push_back( new Resplot("rzed_vs_zed", 0.2*zBins, -zMax, zMax, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
590 //m_rzedres.push_back( new Resplot("rzed_vs_zed", zBins, -zMax, zMax, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
591 m_riptres.push_back( new Resplot("ript_vs_zed", 0.2*zBins, -zMax, zMax, 2*pTResBins, -2*tmp_absResPt, 2*tmp_absResPt ) );
592 m_rptres.push_back( new Resplot("rpt_vs_zed", 0.2*zBins, -zMax, zMax, 8*pTResBins, -tmp_absResPt, tmp_absResPt ) );
593 m_rd0res.push_back( new Resplot("rd0_vs_zed", 0.2*zBins, -zMax, zMax, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
594
595
596 m_retares.push_back( new Resplot("reta_vs_nvtx", 24, 0, 72, 4*etaResBins, -tmp_absResEta, tmp_absResEta ) );
597 m_rphires.push_back( new Resplot("rphi_vs_nvtx", 24, 0, 72, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
598 m_rzedres.push_back( new Resplot("rzed_vs_nvtx", 24, 0, 72, 4*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
599 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_nvtx", 24, 0, 72, 24*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
600 m_riptres.push_back( new Resplot("ript_vs_nvtx", 24, 0, 72, 4*pTResBins, -tmp_absResPt, tmp_absResPt ) );
601 m_rptres.push_back( new Resplot("rpt_vs_nvtx", 24, 0, 72, 8*pTResBins, -tmp_absResPt, tmp_absResPt ) );
602 m_rd0res.push_back( new Resplot("rd0_vs_nvtx", 24, 0, 72, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
603
604
605 m_retares.push_back( new Resplot("reta_vs_ntracks", 60, 0, 600, 4*etaResBins, -tmp_absResEta, tmp_absResEta ) );
606 m_rphires.push_back( new Resplot("rphi_vs_ntracks", 60, 0, 600, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
607 m_rzedres.push_back( new Resplot("rzed_vs_ntracks", 60, 0, 600, 4*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
608 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_ntracks", 60, 0, 600, 24*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
609 //m_rzedres.push_back( new Resplot("rzed_vs_ntracks", 60, 0, 600, 4*zfactor*zresBins, -zfactor*0.5*zresMax, zfactor*0.5*zresMax ) );
610 m_riptres.push_back( new Resplot("ript_vs_ntracks", 60, 0, 600, 4*pTResBins, -tmp_absResPt, tmp_absResPt ) );
611 m_rptres.push_back( new Resplot("rpt_vs_ntracks", 60, 0, 600, 8*pTResBins, -tmp_absResPt, tmp_absResPt ) );
612 m_rd0res.push_back( new Resplot("rd0_vs_ntracks", 60, 0, 600, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
613
614
615 m_retares.push_back( new Resplot("reta_vs_phi", 128, -M_PI, M_PI, 2*etaResBins, -wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
616 m_rphires.push_back( new Resplot("rphi_vs_phi", 128, -M_PI, M_PI, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
617 m_rzedres.push_back( new Resplot("rzed_vs_phi", 128, -M_PI, M_PI, 8*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
618 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_phi", 128, -M_PI, M_PI, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
619 m_riptres.push_back( new Resplot("ript_vs_phi", 128, -M_PI, M_PI, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
620 m_rptres.push_back( new Resplot("rpt_vs_phi", 128, -M_PI, M_PI, 8*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
621 m_rd0res.push_back( new Resplot("rd0_vs_phi", 128, -M_PI, M_PI, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
622
623
624
625 m_retares.push_back( new Resplot("reta_vs_mu", 24, 0, 72, 4*etaResBins, -tmp_absResEta, tmp_absResEta ) );
626 m_rphires.push_back( new Resplot("rphi_vs_mu", 24, 0, 72, 8*phiResBins, -wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
627 m_rzedres.push_back( new Resplot("rzed_vs_mu", 24, 0, 72, 4*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
628 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_mu", 24, 0, 72, 24*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
629 m_riptres.push_back( new Resplot("ript_vs_mu", 24, 0, 72, 4*pTResBins, -tmp_absResPt, tmp_absResPt ) );
630 m_rptres.push_back( new Resplot("rpt_vs_mu", 24, 0, 72, 8*pTResBins, -tmp_absResPt, tmp_absResPt ) );
631 m_rd0res.push_back( new Resplot("rd0_vs_mu", 24, 0, 72, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
632
633
634 //ADDED BY JK
635 //-----
636 m_rzedres.push_back( new Resplot("rzed_vs_signed_pt", ptnbins2, ptbinlims2, 4*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
637 m_rzedres.push_back( new Resplot("rzed_vs_ABS_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
638 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_signed_pt", ptnbins2, ptbinlims2, 24*zfactor*zresBins, -zfactor*zresMax, zfactor*zresMax ) );
639 m_rzedthetares.push_back( new Resplot("rzedtheta_vs_ABS_pt", ptnbins, ptbinlims, 24*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax ) );
640 //m_rzedres.push_back( new Resplot("rzed_vs_signed_pt", ptnbins2, ptbinlims2, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
641 //m_rzedres.push_back( new Resplot("rzed_vs_ABS_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
642 m_rd0res.push_back( new Resplot("rd0_vs_signed_pt", ptnbins2, ptbinlims2, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
643 m_rd0res.push_back( new Resplot("rd0_vs_ABS_pt", ptnbins, ptbinlims, factor*8*a0resBins, -wfactor*a0resMax, wfactor*a0resMax ) );
644 //-----
645
646
647 // std::cout << "booked" << std::endl;
648
649 m_retaresPull.push_back( new Resplot("retaPull_vs_ipt", iptnbins, iptbinlims, 2*etaResBins, -5,5));//-wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
650 m_rphiresPull.push_back( new Resplot("rphiPull_vs_ipt", iptnbins, iptbinlims, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
651 m_rzedresPull.push_back( new Resplot("rzedPull_vs_ipt", iptnbins, iptbinlims, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax ) );
652 //m_riptresPull.push_back( new Resplot("ript_vs_ipt", iptnbins, iptbinlims, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
653 m_rptresPull.push_back( new Resplot("rptPull_vs_ipt", iptnbins, iptbinlims, 8*pTResBins, -5,5));//-wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
654 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_ipt", iptnbins, iptbinlims, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) ) ;
655
656 m_retaresPull.push_back( new Resplot("retaPull_vs_pt", ptnbins, ptbinlims, 2*etaResBins, -5,5));//-wfactor*tmp_absResEta, wfactor*tmp_absResEta ) );
657 m_rphiresPull.push_back( new Resplot("rphiPull_vs_pt", ptnbins, ptbinlims, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
658 m_rzedresPull.push_back( new Resplot("rzedPull_vs_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax ) );
659 //rzedres.push_back( new Resplot("rzed_vs_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
660 //m_riptresPull.push_back( new Resplot("ript_vs_pt", ptnbins, ptbinlims, 16*pTResBins, -wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
661 m_rptresPull.push_back( new Resplot("rptPull_vs_pt", ptnbins, ptbinlims, 8*pTResBins, -5,5));//-wfactor*tmp_absResPt, wfactor*tmp_absResPt ) );
662 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_pt", ptnbins, ptbinlims, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
663
664 m_retaresPull.push_back( new Resplot("retaPull_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*etaResBins, -5,5));//-tmp_absResEta, tmp_absResEta ) );
665 m_rphiresPull.push_back( new Resplot("rphiPull_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi) );
666 m_rzedresPull.push_back( new Resplot("rzedPull_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax) );
667 //rzedres.push_back( new Resplot("rzed_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax) );
668 //riptresPull.push_back( new Resplot("ript_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 16*pTResBins, -tmp_absResPt, tmp_absResPt ) );
669 m_rptresPull.push_back( new Resplot("rptPull_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, 8*pTResBins, -5,5));//-tmp_absResPt, tmp_absResPt ) );
670 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_eta", etaBins, -tmp_maxEta, tmp_maxEta, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ));
671 m_retaresPull.push_back( new Resplot("retaPull_vs_zed", 0.2*zBins, -zMax, zMax, 2*etaResBins, -5,5));//-tmp_absResEta, tmp_absResEta ) );
672 m_rphiresPull.push_back( new Resplot("rphiPull_vs_zed", 0.2*zBins, -zMax, zMax, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
673 m_rzedresPull.push_back( new Resplot("rzedPull_vs_zed", 0.2*zBins, -zMax, zMax, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax ) ) ;
674 //rzedres.push_back( new Resplot("rzed_vs_zed", 0.2*zBins, -zMax, zMax, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
675 //rzedres.push_back( new Resplot("rzed_vs_zed", zBins, -zMax, zMax, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
676 //riptresPull.push_back( new Resplot("ript_vs_zed", 0.2*zBins, -zMax, zMax, 2*pTResBins, -2*tmp_absResPt, 2*tmp_absResPt ) );
677 m_rptresPull.push_back( new Resplot("rptPull_vs_zed", 0.2*zBins, -zMax, zMax, 8*pTResBins, -5,5));//-tmp_absResPt, tmp_absResPt ) );
678 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_zed", 0.2*zBins, -zMax, zMax, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
679
680 m_retaresPull.push_back( new Resplot("retaPull_vs_nvtx", 12, 0, 36, 4*etaResBins, -5,5));//-tmp_absResEta, tmp_absResEta ) );
681 m_rphiresPull.push_back( new Resplot("rphiPull_vs_nvtx", 12, 0, 36, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
682 m_rzedresPull.push_back( new Resplot("rzedPull_vs_nvtx", 12, 0, 36, 4*zfactor*zresBins, -5,5));//-zfactor*zresMax, zfactor*zresMax ) );
683 //rzedres.push_back( new Resplot("rzed_vs_nvtx", 12, 0, 36, 4*zfactor*zresBins, -zfactor*0.5*zresMax, zfactor*0.5*zresMax ) );
684
685 //riptresPull.push_back( new Resplot("ript_vs_nvtx", 12, 0, 36, 4*pTResBins, -tmp_absResPt, tmp_absResPt ) );
686 m_rptresPull.push_back( new Resplot("rptPull_vs_nvtx", 12, 0, 36, 8*pTResBins, -5,5));//-tmp_absResPt, tmp_absResPt ) );
687 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_nvtx", 12, 0, 36, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
688
689 m_retaresPull.push_back( new Resplot("retaPull_vs_ntracks", 60, 0, 600, 4*etaResBins, -5,5));//-tmp_absResEta, tmp_absResEta ) );
690 m_rphiresPull.push_back( new Resplot("rphiPull_vs_ntracks", 60, 0, 600, 8*phiResBins, -5,5));//-wfactor*tmp_absResPhi, wfactor*tmp_absResPhi ) );
691 m_rzedresPull.push_back( new Resplot("rzedPull_vs_ntracks", 60, 0, 600, 4*zfactor*zresBins, -5,5));//-zfactor*zresMax, zfactor*zresMax ) );
692 //rzedres.push_back( new Resplot("rzed_vs_ntracks", 60, 0, 600, 4*zfactor*zresBins, -zfactor*0.5*zresMax, zfactor*0.5*zresMax ) );
693 //riptresPull.push_back( new Resplot("ript_vs_ntracks", 60, 0, 600, 4*pTResBins, -tmp_absResPt, tmp_absResPt ) );
694 m_rptresPull.push_back( new Resplot("rptPull_vs_ntracks", 60, 0, 600, 8*pTResBins, -5,5));//-tmp_absResPt, tmp_absResPt ) );
695 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_ntracks", 60, 0, 600, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
696
697 m_rzedresPull.push_back( new Resplot("rzedPull_vs_signed_pt", ptnbins2, ptbinlims2, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax ) );
698 m_rzedresPull.push_back( new Resplot("rzedPull_vs_ABS_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -5,5));//-2*zfactor*zresMax, 2*zfactor*zresMax ));
699 //rzedres.push_back( new Resplot("rzed_vs_signed_pt", ptnbins2, ptbinlims2, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ));
700 //rzedres.push_back( new Resplot("rzed_vs_ABS_pt", ptnbins, ptbinlims, 4*zfactor*zresBins, -2*zwidthfactor*zresMax, 2*zwidthfactor*zresMax ) );
701 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_signed_pt", ptnbins2, ptbinlims2, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
702 m_rd0resPull.push_back( new Resplot("rd0Pull_vs_ABS_pt", ptnbins, ptbinlims, factor*8*a0resBins, -5,5));//-wfactor*a0resMax, wfactor*a0resMax ) );
703
704 m_rzedreslb = new Resplot("rzed_vs_lb", 301, -0.5, 3009.5, 8*zfactor*zresBins, -2*zfactor*zresMax, 2*zfactor*zresMax );
705
706 m_rzedlb = new Resplot("zed_vs_lb", 301, -0.5, 3009.5, 100, -300, 300 );
707 m_rzedlb_rec = new Resplot("zed_vs_lb_rec", 301, -0.5, 3009.5, 100, -300, 300 );
708
709
710 m_rd0_vs_phi = new Resplot( "d0_vs_phi", 20, -M_PI, M_PI, 200, -10, 10 );
711 m_rd0_vs_phi_rec = new Resplot( "d0_vs_phi_rec", 20, -M_PI, M_PI, 200, -10, 10 );
712
713
715
716 m_rRoi_deta_vs_eta = new Resplot( "rRoi_deta_vs_eta", 20, -2.5, 2.5, 101, -0.5, 0.5 );
717 m_rRoi_dphi_vs_eta = new Resplot( "rRoi_dphi_vs_eta", 20, -2.5, 2.5, 101, -0.5, 0.5 );
718 m_rRoi_dzed_vs_eta = new Resplot( "rRoi_dzed_vs_eta", 20, -2.5, 2.5, 401, -30.0, 30.0 );
719
720 // std::cout << "ROI RESPLOTS " << m_rRoi_deta_vs_eta->Name() << " " << m_rRoi_dphi_vs_eta->Name() << std::endl;
721
722
723 // hit occupancies
724
725 int NHits = 40;
726 int Ntracks = 10000;
727
728 addHistogram( new TH1F( "nsct", "nsct", NHits, -0.5, float(NHits-0.5) ) );
729 addHistogram( new TH1F( "nsct_rec", "nsct_rec", NHits, -0.5, float(NHits-0.5) ) );
730
731 addHistogram( new TH1F( "npix", "npix", NHits, -0.5, float(NHits-0.5) ) );
732 addHistogram( new TH1F( "npix_rec", "npix_rec", NHits, -0.5, float(NHits-0.5) ) );
733
734 addHistogram( new TH1F( "nsi", "nsi", NHits, -0.5, float(NHits-0.5) ) );
735 addHistogram( new TH1F( "nsi_rec", "nsi_rec", NHits, -0.5, float(NHits-0.5) ) );
736 addHistogram( new TH1F( "nsi_matched", "nsi_matched", NHits, -0.5, float(NHits-0.5) ) );
737
738
739 addHistogram( new TH1F( "ntrt", "ntrt", NHits, -0.5, float(NHits-0.5) ) );
740 addHistogram( new TH1F( "ntrt_rec", "ntrt_rec", NHits, -0.5, float(NHits-0.5) ) );
741
742 addHistogram( new TH1F( "nstraw", "nstraw", NHits*4, -0.5, float(4*NHits-0.5) ) );
743 addHistogram( new TH1F( "nstraw_rec", "nstraw_rec", NHits*4, -0.5, float(4*NHits-0.5) ) );
744
745 addHistogram( new TH1F( "ntracks", "ntracks", Ntracks+1, -0.5, float(Ntracks+0.5) ) );
746 addHistogram( new TH1F( "ntracks_rec", "ntracks_rec", Ntracks+1, -0.5, float(Ntracks+0.5) ) );
747
748
749 // beam offset fitting histos
750 m_h2 = new Resplot( "d0vphi", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
751 m_h2r = new Resplot( "d0vphi_rec", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
752 m_h2m = new Resplot( "d0vphi_match", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
753
754 m_h2a0 = new Resplot( "a0vphi", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
755 m_h2a0r = new Resplot( "a0vphi_rec", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
756
758
759
760 TH1F heffvlb("eff vs lb", "eff vs lb", 301, -0.5, 3009.5 );
761
762 // 100,
763 // 1270515000, 1270560000
764 // 1272040000, 1272200000
765 // );
766 // 1270518944,
767 // 1270558762 );
768
769 // 1260470000,
770 // 1260680000 );
771 // 1260470000,
772 // 1260690000 );
773 // TH1F heffvlb("eff vs lb", "eff vs lb", 600, 0, 3000);
774
775 // TH1F heffvlb("eff vs lb", "eff vs lb", 600, 1260400000, 1260700000);
776
777
778 m_eff_vs_lb = new Efficiency1D( &heffvlb );
779
780 // m_z_vs_lb = new Resplot("z vs lb", 100, 1270515000, 1270560000, 100, -250, 250);
781 m_z_vs_lb = new Resplot("z vs lb", 301, -0.5, 3009.5, 100, -250, 250);
782
783 //m_rmap[142165] = 0;
784 //m_rmap[142166] = 500;
785 //m_rmap[142174] = 1000;
786 //m_rmap[142191] = 1500;
787 //m_rmap[142193] = 2000;
788 //m_rmap[142195] = 2500;
789
790
791 m_deltaR_v_eta = new Resplot("R v eta", 10, -2.5, 2.5, 100, 0, 0.1 );
792 m_deltaR_v_pt = new Resplot("R v pt", ptnbins, ptbinlims, 100, 0, 0.1 );
793
794
795 TH1F* eff_vs_mult = new TH1F( "eff_vs_mult", "eff_vs_mult", 25, 0, 25 );
796
797 m_eff_vs_mult = new Efficiency1D( eff_vs_mult, "eff_mult" );
798
799
800 m_n_vtx_tracks = new TH1F("nvtxtracks", "nvtxtracks", 150, 0, 600);
801 m_eff_vs_ntracks = new Efficiency1D( m_n_vtx_tracks, "eff_vs_ntracks");
802
803
804
805 //double oldnbins[21] = { 1, 5, 10, 15, 20,
806 // 25, 30, 35, 40, 45,
807 // 50, 60, 80, 100, 150,
808 // 200, 250, 300, 400, 500, 600 };
809
810
811 //double _nbins[23] = { 0, 29.4, 39.9179, 48.4358, 56.1813, 62.9524, 70.1377, 76.8907,
812 // 83.4667, 90.559, 97.4902, 105.737, 113.547, 121.281, 129.015, 139.824,
813 // 150.589, 164.093, 179.096, 206.867, 400, 500, 600 };
814
815 double nbins[23] = { 0, 29.5, 39.5, 48.5, 56.5, 63.5, 70.5, 77.5,
816 83.5, 91.5, 97.5, 106.5, 114.5, 121.5, 129.5, 140.5,
817 151.5, 164.5, 200.5, 250.5, 300.5, 400.5, 600 };
818
819 TH1F* n_vtx_tracks2 = new TH1F("nvtxtracks2", "nvtxtracks2", 22, nbins);
820 m_eff_vs_ntracks2 = new Efficiency1D( n_vtx_tracks2, "eff_vs_ntracks2");
821 delete n_vtx_tracks2;
822
823 m_n_vtx = new TH1F("nvtx", "nvtx", 81, -0.5, 80.5);
824 m_eff_vs_nvtx = new Efficiency1D( m_n_vtx, "eff_vs_nvtx");
825 //m_mu = new TH1F("mu", "mu", 3000, -0.5, 29.5);
826 m_mu = new TH1F("mu", "mu", 90, 0, 90);
827 m_eff_vs_mu = new Efficiency1D( m_mu, "eff_vs_mu");
828
829
831
832 double etovpt_bins[39] = {
833 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
834 1, 1.08571, 1.17877, 1.2798, 1.3895, 1.50859, 1.63789, 1.77828, 1.9307, 2.09618,
835 2.27585, 2.47091, 2.6827, 2.91263, 3.16228, 3.43332, 3.72759, 4.04709, 4.39397, 4.77058,
836 5.17947, 5.62341, 6.1054, 6.6287, 7.19686, 7.81371, 8.48343, 9.21055, 10
837 };
838
839 m_etovpt_raw = new TH1F("etovpt_raw", "ET / pT", 100, 0, 10 );
840
841 m_etovpt = new TH1F("etovpt", "ET / pT", 38, etovpt_bins );
842 m_eff_vs_etovpt = new Efficiency1D( m_etovpt, "eff_vs_etovpt");
843
844
845 m_et = new TH1F("ET", "ET; E_{T} [GeV]", ptnbins, ptbinlims );
846 m_eff_vs_et = new Efficiency1D( m_et, "eff_vs_ET" );
847
848
849 // std::cout << "initialize() Directory " << gDirectory->GetName() << " on leaving" << std::endl;
850
851 ConfVtxAnalysis* vtxanal = 0;
852 store().find( vtxanal, "rvtx" );
853 if ( vtxanal ) vtxanal->initialise();
854
855 // std::cout << "initialize() Directory " << gDirectory->GetName() << " on leaving" << std::endl;
856
857
858 m_dir->pop();
859
860 dir->cd();
861
862
863
864}
865
866
867
868
869
870// fit Gaussian to histogram
871
872TF1* FitFWGaussian(TH1D* s, double a, double b) {
873
874 // std::cout << "FitFWGaussian " << s->GetName() << std::endl;
875
876 // TF1* f1 = new TF1("gausr", "gaus");
877 TF1* f1 = new TF1("gausr", "[0]*exp(-(x-[1])*(x-[1])/([2]*[2]))");
878
879 f1->SetNpx(5000);
880 f1->SetLineWidth(1);
881
882 f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
883 f1->SetParameter(1,s->GetBinCenter(s->GetMaximumBin()));
884 // f1->SetParameter(2,s->GetRMS());
885 f1->SetParameter(2,1.5);
886 f1->FixParameter(2,1.5);
887
888 f1->SetParName(0, "Norm");
889 f1->SetParName(1, "Mean");
890 f1->SetParName(2, "Sigma");
891
892 int nbins=0;
893 for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)) nbins++;
894
895 if (nbins>2) {
896 if ( a==-999 || b==-999 ) s->Fit(f1,"Q");
897 else s->Fit(f1,"Q", "", a, b);
898 }
899 else for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j, 0);
900
901 // std::cout << "return" << std::endl;
902
903 return f1;
904}
905
906
907
908
909
910void fitSin( TH1D* h, const std::string& parent="" ) {
911
912 TF1* fsin = new TF1( "sinp", "sqrt([0]*[0])*sin([1]-x)+[2]" ); // , -M_PI, M_PI );
913
914 fsin->SetParameter(0,1);
915 fsin->SetParameter(1,0);
916 fsin->SetParameter(2,0);
917
918 fsin->SetLineWidth(1);
919
920 // h->SetTitle(";track #phi;d_{0} [mm]");
921
922 std::cout << __FUNCTION__ << " " << h->GetName() << std::endl;
923 h->Fit( fsin, "", "", -M_PI, M_PI );
924
925 double radius = fsin->GetParameter(0);
926 double phi = fsin->GetParameter(1);
927
928 double y = radius*sin(phi);
929 double x = radius*cos(phi);
930
931 std::cout << parent << "\t" << h->GetTitle() << "\tx = " << x << "\ty = " << y << std::endl;
932
933 delete fsin;
934}
935
936
937
938
940
942
943 // gDirectory->pwd();
944
945 if ( !m_initialised ) return;
946
947 std::cout << "ConfAnalysis::finalise() " << name();
948
949 if ( name().size()<19 ) std::cout << "\t";
950 if ( name().size()<30 ) std::cout << "\t";
951 if ( name().size()<41 ) std::cout << "\t";
952 if ( name().size()<52 ) std::cout << "\t";
953
954
955 std::cout << "\tNreco " << m_Nreco
956 << "\tNref " << m_Nref
957 << "\tNmatched " << m_Nmatched;
958
959 if (m_Nref) {
960 std::cout << " tracks approx " << (100*m_Nmatched)/m_Nref << "%" ;
961 }
962 std::cout << std::endl;
963
964 // if ( m_Nreco==0 ) return;
965
966 // TIDDirectory d( name() );
967 // d.push();
968
969 m_dir->push();
970
971 std::map<std::string, TH1F*>::iterator hitr=m_histos.begin();
972 std::map<std::string, TH1F*>::iterator hend=m_histos.end();
973 for ( ; hitr!=hend ; ++hitr ) hitr->second->Write();
974 // std::cout << "DBG >" << eff_pt->Hist()->GetName() << "< DBG" << std::endl;
975
976 // std::vector<Efficiency*> heff = { eff_pt,
977
978 const unsigned Neff = 11;
979 Efficiency1D* heff[Neff] = { m_eff_pt,
980 m_eff_eta,
981 m_eff_phi,
982 m_eff_z0,
983 m_eff_d0,
984 m_eff_a0,
985 m_eff_ptm,
986 m_eff_ptp,
989 m_eff_roi_dR };
990
991 for ( unsigned i=0 ; i<Neff ; i++ ) {
992 heff[i]->finalise();
993 heff[i]->Bayes()->Write( ( heff[i]->name()+"_tg" ).c_str() );
994 } // heff[i]->Hist()->Write(); }
995
996 // std::cout << "DBG >" << m_purity_pt->Hist()->GetName() << "< DBG" << std::endl;
997
998 m_eff_vs_mult->finalise();
999
1000 // Normalise(n_vtx_tracks);
1001
1002 m_eff_vs_ntracks->finalise();
1003 m_eff_vs_ntracks2->finalise();
1004
1005 m_eff_vs_nvtx->finalise();
1006 m_eff_vs_mu->finalise();
1007
1008 m_eff_vs_etovpt->finalise();
1009
1010 m_eff_vs_et->finalise();
1011
1012 m_eff_eta_vs_pt->finalise();
1013 m_eff_d0_vs_pt->finalise();
1014
1015
1016 const unsigned Npurity = 6;
1017 Efficiency1D* hpurity[Npurity] = {
1023 m_purity_a0 };
1024
1025 for ( unsigned i = 0 ; i<Npurity ; i++ ) hpurity[i]->finalise();
1026
1028 fitSin( m_rd0_vs_phi->Mean(), name()+"/rd0_vs_phi" );
1029 m_rd0_vs_phi->Write();
1030
1032 fitSin( m_rd0_vs_phi_rec->Mean(), name()+"/rd0_vs_phi_rec" );
1033 m_rd0_vs_phi_rec->Write();
1034
1035
1036 std::string spstr[5] = { "npix", "nsct", "nsi", "ntrt", "nbl" };
1037 for ( int i=m_res.size() ; i-- ; ) {
1038 TF1* (*resfit)(TH1D* s, double a, double b) = Resplot::FitNull95;
1039 for ( int ir=0 ; ir<5 ; ir++ ) if ( m_res[i]->Name().find(spstr[ir])!=std::string::npos ) { resfit = Resplot::FitNull; break; }
1040 m_res[i]->Finalise( resfit );
1041 m_res[i]->Write();
1042 delete m_res[i];
1043 }
1044
1045 m_deltaR_v_eta->Finalise(); m_deltaR_v_eta->Write();
1046 m_deltaR_v_pt->Finalise(); m_deltaR_v_pt->Write();
1047
1048 for ( unsigned i=m_rDd0res.size() ; i-- ; ) {
1049 m_rDd0res[i]->Finalise(Resplot::FitNull95);
1050 m_rDd0res[i]->Write();
1051 delete m_rDd0res[i];
1052 }
1053
1054 for ( unsigned i=m_rDa0res.size() ; i-- ; ) {
1055 m_rDa0res[i]->Finalise(Resplot::FitNull95);
1056 m_rDa0res[i]->Write();
1057 delete m_rDa0res[i];
1058 }
1059
1060 for ( unsigned i=m_rDz0res.size() ; i-- ; ) {
1061 m_rDz0res[i]->Finalise(Resplot::FitNull95);
1062 m_rDz0res[i]->Write();
1063 delete m_rDz0res[i];
1064 }
1065
1066
1067 for ( unsigned ih=0 ; ih<2 ; ih++ ) {
1068 m_hphivsDd0res[ih]->Divide( m_hphivsDd0res[2] );
1069 m_hphivsDa0res[ih]->Divide( m_hphivsDa0res[2] );
1070 }
1071
1072 for ( unsigned ih=0 ; ih<3 ; ih++ ) {
1073 m_hphivsDd0res[ih]->Write();
1074 m_hphivsDa0res[ih]->Write();
1075
1076 delete m_hphivsDd0res[ih];
1077 delete m_hphivsDa0res[ih];
1078 }
1079
1081
1085
1086 m_rRoi_deta_vs_eta->Write();
1087 m_rRoi_dphi_vs_eta->Write();
1088 m_rRoi_dzed_vs_eta->Write();
1089
1090 delete m_rRoi_deta_vs_eta;
1091 delete m_rRoi_dphi_vs_eta;
1092 delete m_rRoi_dzed_vs_eta;
1093
1095
1096
1097
1098 for ( unsigned i=m_retares.size() ; i-- ; ) {
1099
1100#if 1
1101
1102 m_retares[i]->Finalise(Resplot::FitNull95);
1103 m_rphires[i]->Finalise(Resplot::FitNull95);
1104 m_riptres[i]->Finalise(Resplot::FitNull95);
1105 m_rzedres[i]->Finalise(Resplot::FitNull95);
1106 m_rzedthetares[i]->Finalise(Resplot::FitNull95);
1107 // m_rptres[i]->Finalise(Resplot::FitBreit);
1108 //m_rptres[i]->Finalise(Resplot::FitNull95);
1109 m_rd0res[i]->Finalise(Resplot::FitNull95);
1110 // m_rd0res[i]->Finalise(Resplot::FitCentralGaussian);
1111 // m_rd0res_rms[i]->Finalise(Resplot::FitNull);
1112
1113#else
1114
1115 m_retares[i]->Finalise();
1116 m_rphires[i]->Finalise();
1117 m_rzedres[i]->Finalise();
1118 m_riptres[i]->Finalise();
1119 m_rptres[i]->Finalise();
1120 m_rd0res[i]->Finalise();
1121
1122#endif
1123
1124 m_retares[i]->Write();
1125 m_rphires[i]->Write();
1126 m_rzedres[i]->Write();
1127 m_rzedthetares[i]->Write();
1128 m_riptres[i]->Write();
1129 m_rptres[i]->Write();
1130 m_rd0res[i]->Write();
1131
1132 delete m_retares[i];
1133 delete m_rphires[i];
1134 delete m_rzedres[i];
1135 delete m_rzedthetares[i];
1136 delete m_riptres[i];
1137 delete m_rptres[i];
1138 delete m_rd0res[i];
1139
1140
1141 }
1142
1144
1145 m_rzedlb->Finalise(Resplot::FitNull95);
1147
1148 for ( unsigned i=m_retaresPull.size() ; i-- ; ) {
1149
1150 m_retaresPull[i]->Finalise(Resplot::FitNull);
1151 m_rphiresPull[i]->Finalise(Resplot::FitNull);
1152 m_rptresPull[i]->Finalise(Resplot::FitNull);
1153 m_rzedresPull[i]->Finalise(Resplot::FitNull);
1154 m_rd0resPull[i]->Finalise(Resplot::FitNull);
1155
1156 m_retaresPull[i]->Write();
1157 m_rphiresPull[i]->Write();
1158 m_rptresPull[i]->Write();
1159 m_rzedresPull[i]->Write();
1160 m_rd0resPull[i]->Write();
1161
1162 delete m_retaresPull[i];
1163 delete m_rphiresPull[i];
1164 delete m_rptresPull[i];
1165 delete m_rzedresPull[i];
1166 delete m_rd0resPull[i];
1167
1168 }
1169
1170
1171 m_rzedreslb->Write();
1172 delete m_rzedreslb;
1173
1174 m_rzedlb->Write();
1175 delete m_rzedlb;
1176
1177 m_rzedlb_rec->Write();
1178 delete m_rzedlb_rec;
1179
1180 // td->cd();
1181
1182 //ADDED BY JK
1183 //-----Only one more element in d0 and z0 vectors than eta now
1184 m_rzedres[m_rzedres.size()-1]->Finalise(Resplot::FitNull95);
1185 m_rzedres[m_rzedres.size()-1]->Write();
1186 m_rzedthetares.back()->Finalise(Resplot::FitNull95);
1187 m_rzedthetares.back()->Write();
1188 m_rd0res.back()->Finalise(Resplot::FitNull95);
1189 m_rd0res.back()->Write();
1190 //-----
1191
1192 m_eff_vs_lb->finalise();
1193
1194 m_z_vs_lb->Finalise(Resplot::FitNull95); m_z_vs_lb->Write();
1195 delete m_z_vs_lb;
1196
1197 // TH1F* hefflb = eff_vs_lb->Hist();
1198 // hefflb->Fit("pol0");
1199
1200 // ConfVtxAnalysis* vtxanal = 0;
1201 // store().find( vtxanal, "rvtx" );
1202 // if ( vtxanal ) vtxanal->finalise();
1203
1204 // d.pop();
1205
1206 ConfVtxAnalysis* vtxanal = 0;
1207 store().find( vtxanal, "rvtx" );
1208 if ( vtxanal ) vtxanal->finalise();
1209
1210 m_dir->pop();
1211
1212}
1213
1214extern int Nvtxtracks;
1215extern int NvtxCount;
1216
1217
1219
1220
1221double wrapphi( double phi ) {
1222 if ( phi<-M_PI ) phi += 2*M_PI;
1223 if ( phi> M_PI ) phi -= 2*M_PI;
1224 return phi;
1225}
1226
1227
1228void ConfAnalysis::execute( const std::vector<TIDA::Track*>& reftracks,
1229 const std::vector<TIDA::Track*>& testtracks,
1230 TrackAssociator* matcher,
1231 TrigObjectMatcher* objects ) {
1232
1233 // leave this commented code in for debug purposes ...
1234 // if ( objects ) std::cout << "TrigObjectMatcher: " << objects << std::endl;
1235
1237
1238 if ( m_print ) {
1239 std::cout << "ConfAnalysis::execute() \t " << name()
1240 << "\tref " << reftracks.size()
1241 << "\ttest " << testtracks.size();
1242 if ( groi ) std::cout << "\tgroi " << groi << " " << *groi;
1243 std::cout << std::endl;
1244 }
1245
1246 // std::cout << "ConfAnalysis (resolutions really) filling " << std::endl;
1247
1248 // should have these as a class variable
1249 static std::string varName[16] = { "pT", "eta", "phi", "z0", "d0", "a0",
1250 "nsct", "npix", "nsi", "ntrt", "nstraw",
1251 "dd0", "da0", "dz0", "deta", "dphi" };
1252
1253 std::map<std::string, TH1F*>::iterator hmitr = m_histos.find("ntracks");
1254 if ( hmitr!=m_histos.end() ) hmitr->second->Fill( reftracks.size() );
1255
1256 hmitr = m_histos.find("ntracks_rec");
1257 if ( hmitr!=m_histos.end() ) hmitr->second->Fill( testtracks.size() );
1258
1259 bool dump = false;
1260
1261 m_Nreco += testtracks.size();
1262 m_Nref += reftracks.size();
1263
1264
1265 // std::cout << "ConfAnalysis ref tracks " << std::endl;
1266
1267 // why don't we use this vertex position any more ???
1268 // m_Nref = 0;
1269 // for ( int i=reftracks.size() ; i-- ; ) {
1270 // double phit = reftracks[i]->phi();
1271 // // double a0t = reftracks[i]->a0() + sin(phit)*m_xBeamReference - cos(phit)*m_yBeamReference;
1272 // // if ( std::fabs(a0t)<a0 ) m_Nref++; ???
1273 // }
1274
1275 // if ( testtracks.size() ) std::cout << "NTRACKS " << testtracks.size() << std::endl;
1276
1277 for ( int i=reftracks.size() ; i-- ; ) {
1278
1280
1281 if ( groi!=0 ) {
1282 // std::cout << "ConfAnalysis::Fill() groi " << *groi << std::endl;
1283
1284 double deta = reftracks[i]->eta() - groi->eta();
1285 double dphi = wrapphi( reftracks[i]->phi() - groi->phi() );
1286 double dzed = reftracks[i]->z0() - groi->zed();
1287
1288 m_rRoi_deta_vs_eta->Fill( groi->eta(), deta );
1289 m_rRoi_dphi_vs_eta->Fill( groi->eta(), dphi );
1290 m_rRoi_dzed_vs_eta->Fill( groi->eta(), dzed );
1291 }
1292
1294 double ipTt = 1./(reftracks[i]->pT()/1000.);
1295 double pTt = reftracks[i]->pT()/1000;
1296
1297 double etat = reftracks[i]->eta();
1298 double phit = reftracks[i]->phi();
1299
1300 double thetat = 2*std::atan( exp( (-1)*etat ) );
1301
1306
1308
1309 double z0t = reftracks[i]->z0();
1310 double d0t = reftracks[i]->a0() - std::sin(phit)*m_xBeamReference + std::cos(phit)*m_yBeamReference;
1311 double a0t = reftracks[i]->a0();
1312
1313 if ( m_xBeamReference!=0 || m_yBeamReference!=0 ) {
1316 z0t = reftracks[i]->z0()+((std::cos(phit)*m_xBeamReference + std::sin(phit)*m_yBeamReference)/std::tan(thetat));
1317 d0t = reftracks[i]->a0();
1318 a0t = reftracks[i]->a0() + std::sin(phit)*m_xBeamReference - std::cos(phit)*m_yBeamReference;
1319 }
1320
1321 // if ( m_lfirst ) {
1322 // std::cout << "\na0t " << a0t << "(shifted " << d0t << ")" << std::endl;
1323 // std::cout << "\nz0t " << z0t << std::endl;
1324
1325 // std::cout << "\txBeamReference " << m_xBeamReference << "\tyBeamReference " << m_yBeamReference << std::endl;
1326 // }
1327
1328
1329
1330 // std::cout << "a0t " << a0t << std::endl;
1331
1332 m_rChi2prob->Fill( pTt, TMath::Prob(reftracks[i]->chi2(),reftracks[i]->dof()) );
1333 m_rChi2->Fill( pTt, reftracks[i]->chi2() );
1334 m_rChi2dof->Fill( pTt, reftracks[i]->chi2()/reftracks[i]->dof() );
1335
1336 // static double xbeamref = 0;
1337 // if ( m_lfirst || m_xBeamReference!=xbeamref) {
1338 // std::cout << __FUNCTION__ << "\tbeamline " << m_xBeamReference << " " << m_yBeamReference << " (ref)" << std::endl;
1339 // }
1340 // xbeamref = m_xBeamReference;
1341
1342
1344 double dpTt = reftracks[i]->dpT()/1000;
1345 double detat = reftracks[i]->deta();
1346 double dphit = reftracks[i]->dphi();
1347
1348 //RoI variables
1349 float droi_detat = 0.f;
1350 float droi_dphit = 0.f;
1351 float droi_dRt = 0.f;
1352 if (groi){ //checked for null on line 1281
1353 droi_detat = groi->eta() - reftracks[i]->eta();
1354 droi_dphit = groi->phi() - reftracks[i]->phi();
1355 if ( droi_dphit<-M_PI ) droi_dphit +=2*M_PI;
1356 if ( droi_dphit>M_PI ) droi_dphit -=2*M_PI;
1357 droi_dRt = std::sqrt(droi_dphit*droi_dphit + droi_detat*droi_detat);
1358 }
1359
1360 // double dz0t = reftracks[i]->dz0()+((std::cos(phit)*m_xBeamReference + std::sin(phit)*m_yBeamReference)/std::tan(thetat));
1361 // double dd0t = reftracks[i]->da0();
1362 double dz0t = reftracks[i]->dz0();
1363 double dd0t = reftracks[i]->da0() - std::sin(phit)*m_xBeamReference + std::cos(phit)*m_yBeamReference;
1364
1365 // this will be changed when we know the beam spot position
1366 // double a0t = reftracks[i]->a0() + sin(phit)*m_xBeam - cos(phit)*m_yBeam;
1367 double da0t = reftracks[i]->da0();
1368
1369#if 0
1370 std::cout << "etat = " << etat << " +/- " << detat << std::endl;
1371 std::cout << "phit = " << phit << " +/- " << dphit << std::endl;
1372 std::cout << "z0t = " << z0t << " +/- " << dz0t << std::endl;
1373 std::cout << "d0t = " << d0t << " +/- " << dd0t << std::endl;
1374 std::cout << "a0t = " << a0t << " +/- " << da0t << std::endl;
1375 std::cout << "pTt = " << pTt << " +/- " << dpTt << std::endl;
1376#endif
1377
1378 if ( std::fabs(a0t)>a0 ) continue;
1379
1380 // double chi2t = reftracks[i]->chi2();
1381 // hchi2->Fill( chi2t );
1382
1383 double nsctt = reftracks[i]->sctHits();
1384 double npixt = reftracks[i]->pixelHits()*0.5;
1385 double nsit = reftracks[i]->pixelHits()*0.5 + reftracks[i]->sctHits();
1386
1387 double nsctht = reftracks[i]->sctHoles();
1388 double npixht = reftracks[i]->pixelHoles();
1389 double nsiht = reftracks[i]->pixelHoles() + reftracks[i]->sctHoles();
1390
1391 double nbl = reftracks[i]->bLayerHits();
1392 double nblh = ( ( reftracks[i]->expectBL() && reftracks[i]->bLayerHits()<1 ) ? 1 : 0 );
1393
1394 // double ntrtt = reftracks[i]->trHits();
1395 double nstrawt = reftracks[i]->strawHits();
1396
1397 // double ts_scale = (ts-1260400000)*3000.0/(1260700000-1260400000);
1398
1399 // std::cout << "Fill h2 " << " " << m_h2m << " " << *reftracks[i] << std::endl;
1400
1401 m_h2->Fill( phit, d0t );
1402 m_h2a0->Fill( phit, a0t );
1403
1404 m_rd0_vs_phi->Fill( phit, a0t );
1405
1406 double mu_val = gevent->mu();
1407
1408 m_mu->Fill( mu_val );
1409
1410
1411 const TIDA::Track* matchedreco = matcher->matched(reftracks[i]);
1412
1413 // std::cout << "\t\tConfAnalysis " << name() << "\t" << i << " " << *reftracks[i] << " -> ";
1414
1415 // raw reference track distributions
1416 double vpart[16] = { std::fabs(pTt), etat, phit, z0t, d0t, a0t, nsctt, npixt, nsctt, nsit, nstrawt, dd0t, da0t, dz0t, detat, dphit };
1417
1420 for ( int it=0 ; it<11 ; it++ ) {
1421 // std::string hname = varName[it];
1422 // std::map<std::string, TH1F*>::iterator hmitr = m_histos.find(hname);
1423 // if ( hmitr!=m_histos.end() ) hmitr->second->Fill( vpart[it] );
1424
1425 if ( TH1F* hptr = find( varName[it] ) ) hptr->Fill( vpart[it] );
1426 else std::cerr << "hmmm histo " << varName[it] << " not found" << std::endl;
1427 }
1428
1429
1430 m_rnpix_eta->Fill( etat, npixt*1.0 );
1431 m_rnsct_eta->Fill( etat, nsctt*1.0 );
1432 m_rntrt_eta->Fill( etat, nstrawt*1.0 );
1433 m_rnsihit_eta->Fill( etat, npixt + nsctt*1.);
1434
1435 m_rnpix_phi->Fill( phit, npixt*1.0 );
1436 m_rnsct_phi->Fill( phit, nsctt*1.0 );
1437 m_rntrt_phi->Fill( phit, nstrawt*1.0 );
1438
1439 m_rnpix_pt->Fill( std::fabs(pTt), npixt*1.0 );
1440 m_rnsct_pt->Fill( std::fabs(pTt), nsctt*1.0 );
1441 m_rntrt_pt->Fill( std::fabs(pTt), nstrawt*1.0 );
1442
1443
1444 m_rnpix_d0->Fill( a0t, npixt*1.0 );
1445 m_rnsct_d0->Fill( a0t, nsctt*1.0 );
1446 m_rntrt_d0->Fill( a0t, nstrawt*1.0 );
1447
1448 m_rnpixh_d0->Fill( a0t, npixht );
1449 m_rnscth_d0->Fill( a0t, nsctht );
1450
1451 m_rnsi_pt->Fill( std::fabs(pTt), nsit );
1452 m_rnsih_pt->Fill( std::fabs(pTt), nsiht );
1453
1454 m_rnsi_d0->Fill( a0t, nsit );
1455 m_rnsih_d0->Fill( a0t, nsiht );
1456
1457 m_rnsi_eta->Fill( etat, nsit );
1458 m_rnsih_eta->Fill(etat, nsiht );
1459
1460 m_rnbl_d0->Fill( a0t, nbl );
1461 m_rnblh_d0->Fill( a0t, nblh );
1462
1463
1464 m_rnpixh_pt->Fill( std::fabs(pTt), npixht );
1465 m_rnscth_pt->Fill( std::fabs(pTt), nsctht );
1466
1467 m_rnpix_lb->Fill( gevent->lumi_block(), npixt*1.0 );
1468 m_rnsct_lb->Fill( gevent->lumi_block(), nsctt*1.0 );
1469
1470
1471 m_rnsct_vs_npix->Fill( npixt, nsctt );
1472
1473 double etovpt_val = 0;
1474 const TrackTrigObject* tobj = 0;
1475
1476 double ET=0;
1477
1478 if ( objects ) {
1479 tobj = objects->object( reftracks[i]->id() );
1480 if ( tobj ) {
1483 etovpt_val = std::fabs( tobj->pt()/reftracks[i]->pT() );
1484 m_etovpt->Fill( etovpt_val );
1485 m_etovpt_raw->Fill( etovpt_val );
1486 m_et->Fill( tobj->pt()*0.001 );
1487 ET = std::fabs(tobj->pt()*0.001);
1488 }
1489 }
1490
1491
1492 if ( matchedreco ) {
1493
1494 m_Nmatched++;
1495
1496 // efficiency histos
1497 m_eff_pt->Fill(std::fabs(pTt));
1498 m_eff_z0->Fill(z0t);
1499 m_eff_eta->Fill(etat);
1500 m_eff_phi->Fill(phit);
1501 m_eff_d0->Fill(d0t);
1502 m_eff_a0->Fill(a0t);
1503 m_eff_roi_deta->Fill(droi_detat);
1504 m_eff_roi_dphi->Fill(droi_dphit);
1505 m_eff_roi_dR->Fill(droi_dRt);
1506
1507 // signed pT
1508 if ( pTt<0 ) m_eff_ptm->Fill(std::fabs(pTt));
1509 else m_eff_ptp->Fill(std::fabs(pTt));
1510
1511 m_eff_vs_mult->Fill( m_Nref );
1512
1513
1514 m_eff_eta_vs_pt->Fill( std::fabs(pTt), etat );
1515 m_eff_d0_vs_pt->Fill( std::fabs(pTt), a0t );
1516
1517
1518
1519 // m_eff_vs_lb->Fill( m_rmap[r]+lb );
1520 // m_eff_vs_lb->Fill( ts_scale );
1521 // m_eff_vs_lb->Fill( ts );
1522
1523 m_eff_vs_lb->Fill( gevent->lumi_block() );
1524
1525 if ( tobj ) {
1526 m_eff_vs_etovpt->Fill(etovpt_val);
1527 m_eff_vs_et->Fill( std::fabs(tobj->pt()*0.001) );
1528 }
1529
1530
1532
1534 double pTr = matchedreco->pT()/1000;
1535 double etar = matchedreco->eta();
1536 double phir = matchedreco->phi();
1537 //double z0r = matchedreco->z0() + std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest; ;
1538 double thetar = 2*std::atan( exp( (-1)*etar) );
1539
1540 // double z0r = matchedreco->z0()+((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
1541 // double d0r = matchedreco->a0();
1542 // double a0r = matchedreco->a0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
1543
1544
1545
1546
1547
1548 // static bool tfirst = true;
1549 // static double xbeamtest = 0;
1550 // if ( m_lfirst || xbeamtest!=m_xBeamTest) {
1551 // std::cout << __FUNCTION__ << "\tbeamline " << m_xBeamTest << " " << m_yBeamTest << " (test)" << std::endl;
1552 // }
1553 // xbeamtest = m_xBeamTest;
1554
1555 double d0r = 0;
1556 double a0r = 0;
1557 double z0r = 0;
1558
1559 d0r = matchedreco->a0() - sin(phir)*m_xBeamTest + cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
1560 a0r = matchedreco->a0();
1561 z0r = matchedreco->z0();
1562
1563 if ( m_xBeamTest!=0 || m_yBeamTest!=0 ) {
1564 d0r = matchedreco->a0();
1565 a0r = matchedreco->a0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
1566 z0r = matchedreco->z0()+((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
1567 }
1568
1569
1570 // if ( m_lfirst ) {
1571 // std::cout << "\na0r " << a0r << "(shifted " << d0r << ")" << std::endl;
1572 // std::cout << "\nz0r " << z0r << std::endl;
1573 // std::cout << "\txBeamReference " << m_xBeamReference << "\tyBeamReference " << m_yBeamReference << std::endl;
1574 // }
1575
1576
1577 // static int it=0;
1578
1579 // std::cout << "it " << it++ << std::endl;
1580
1581 // m_lfirst = false;
1582
1583 double nsctr = matchedreco->sctHits();
1584 double npixr = matchedreco->pixelHits()*0.5;
1585 double nsir = matchedreco->pixelHits()*0.5 + matchedreco->sctHits();
1586
1587 double nscthr = matchedreco->sctHoles();
1588 double npixhr = matchedreco->pixelHoles();
1589
1590
1591 m_rnpix_lb_rec->Fill( gevent->lumi_block(), npixr*1.0 );
1592 m_rnsct_lb_rec->Fill( gevent->lumi_block(), nsctr*1.0 );
1593
1594
1595 //double ntrtr = matchedreco->trHits();
1596 double nstrawr = matchedreco->strawHits();
1597
1599
1600 double dpTr = matchedreco->dpT()/1000;
1601 double detar = matchedreco->deta();
1602 double dphir = matchedreco->dphi();
1603
1604 // double dz0r = matchedreco->dz0()+((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
1605 // double dd0r = matchedreco->da0();
1606 // double da0r = matchedreco->da0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest;
1607
1608 double dz0r = matchedreco->dz0();
1609 double dd0r = matchedreco->da0() - sin(phir)*m_xBeamTest + cos(phir)*m_yBeamTest;
1610 double da0r = matchedreco->da0();
1611
1612#if 0
1613 std::cout << "etar = " << etar << " +/- " << detar << std::endl;
1614 std::cout << "phir = " << phir << " +/- " << dphir << std::endl;
1615 std::cout << "pTr = " << pTr << " +/- " << dpTr << std::endl;
1616 std::cout << "a0r = " << a0r << " +/- " << da0r << std::endl;
1617 std::cout << "d0r = " << d0r << " +/- " << dd0r << std::endl;
1618 std::cout << "z0r = " << z0r << " +/- " << dz0r << std::endl;
1619#endif
1620
1621 if ( m_h2m ) m_h2m->Fill( phit, d0t );
1622
1623 m_rd0_vs_phi_rec->Fill( phir, a0r );
1624
1625
1627 double resfiller[9] = { std::fabs(ipTt), std::fabs(pTt), ET, etat, z0t, double(NvtxCount), double(Nvtxtracks), phit, mu_val };
1628
1629 double Delta_ipt = 1.0/pTr - 1.0/pTt;
1630 double Delta_pt = pTr - pTt;
1631
1632 if ( pTt<0 ) {
1633 Delta_ipt *= -1;
1634 Delta_pt *= -1;
1635 }
1636
1637 for ( int irfill=0 ; irfill<9 ; irfill++ ) {
1638 m_retares[irfill]->Fill( resfiller[irfill], etar-etat );
1639 m_rphires[irfill]->Fill( resfiller[irfill], phir-phit );
1640 m_rzedres[irfill]->Fill( resfiller[irfill], z0r-z0t );
1641 m_rzedthetares[irfill]->Fill( resfiller[irfill], z0r*std::sin(thetar)-z0t*std::sin(thetat) );
1642 m_riptres[irfill]->Fill( resfiller[irfill], Delta_ipt );
1643 m_rptres[irfill]->Fill( resfiller[irfill], Delta_pt );
1644 m_rd0res[irfill]->Fill( resfiller[irfill], a0r-a0t );
1645 }
1646
1647 double lb = gevent->lumi_block();
1648
1649 m_rzedreslb->Fill( lb, z0r-z0t );
1650 m_rzedlb->Fill( lb, z0t );
1651 m_rzedlb_rec->Fill( lb, z0r );
1652
1653 for ( int irfill=0 ; irfill<6 ; irfill++ ) {
1654 m_rphiresPull[irfill]->Fill( resfiller[irfill], (phir - phit) / sqrt( (dphit*dphit) + (dphir*dphir) ) );
1655 m_retaresPull[irfill]->Fill( resfiller[irfill], (etar - etat) / sqrt( (detat*detat) + (detar*detar) ) );
1656 m_rptresPull[irfill]->Fill( resfiller[irfill], Delta_pt / sqrt( (dpTt*dpTt) + (dpTr*dpTr) ) );
1657 m_rzedresPull[irfill]->Fill( resfiller[irfill], (z0r - z0t) / sqrt( (dz0t*dz0t) + (dz0r*dz0r) ) );
1658 m_rd0resPull[irfill]->Fill( resfiller[irfill], (a0r - a0t) / sqrt( (da0t*da0t) + (da0r*da0r) ) );
1659 }
1660
1661 m_rDz0res[0]->Fill( std::fabs(pTt), dz0r-dz0t );
1662 m_rDz0res[1]->Fill( etat, dz0r-dz0t );
1663 m_rDz0res[2]->Fill( z0t, dz0r-dz0t );
1664 m_rDz0res[3]->Fill( phit, dz0r-dz0t );
1665
1666 if ( dumpflag ) {
1667 std::ostream& dumpstream = dumpfile;
1668 if ( dz0t>0 && std::fabs( dz0r-dz0t )>0.04 ) {
1669 dump = true;
1670 dumpstream << "POOR sigma(z0) agreement \n\trefrack: " << *reftracks[i] << "\n\ttestrack: " << *matchedreco << std::endl;
1671 // std::cout << "dz0r dz0t" << dz0r << "\t" << dz0t << std::endl;
1672 }
1673 }
1674
1676 m_rDd0res[0]->Fill( std::fabs(pTt), dd0r-dd0t );
1677 m_rDd0res[1]->Fill( etat, dd0r-dd0t );
1678 m_rDd0res[2]->Fill( z0t, dd0r-dd0t );
1679 m_rDd0res[3]->Fill( d0t, dd0r-dd0t );
1680 m_rDd0res[4]->Fill( phit, dd0r-dd0t );
1681
1682 m_rDa0res[0]->Fill( std::fabs(pTt), da0r-da0t );
1683 m_rDa0res[1]->Fill( etat, da0r-da0t );
1684 m_rDa0res[2]->Fill( z0t, da0r-da0t );
1685 m_rDa0res[3]->Fill( da0t, da0r-da0t );
1686
1687 const double Deltaphi = 2*M_PI/NMod;
1688
1689 double phistart = 11.0819;
1690 if ( NMod==22 ) phistart = 7.05803;
1691
1692 double phit_wrap = phit - phistart*M_PI/180;
1693
1694 if ( phit_wrap<-M_PI ) phit_wrap += 2*M_PI;
1695
1696 double iphi = phit_wrap - (Deltaphi*int((phit_wrap+M_PI)/Deltaphi) - M_PI);
1697
1698 // double iphi = phit - M_PI*int(7*(phit+M_PI)/M_PI)/7 - 11.0819*M_PI/180 + M_PI;
1699 // double iphi = phit - M_PI*int(7*phit/M_PI)/7.0 - 11.0819*M_PI/180;
1700
1701 m_rDa0res[4]->Fill( phit, da0r-da0t );
1702 m_rDa0res[5]->Fill( iphi, da0r-da0t );
1703
1704 m_rDa0res[6]->Fill( iphi, da0t );
1705 m_rDa0res[7]->Fill( iphi, da0r );
1706
1707 //ADDED BY JK
1708 //-----
1711 m_rzedres[m_rphires.size()-1]->Fill( resfiller[1], z0r-z0t );
1712 m_rzedres[m_rphires.size()+1]->Fill( resfiller[1], z0r-z0t );
1713
1714 m_rzedthetares[m_rphires.size()-1]->Fill( resfiller[1], z0r*std::sin(thetar)-z0t*std::sin(thetat) );
1715 m_rzedthetares[m_rphires.size()]->Fill( resfiller[1], z0r*std::sin(thetar)-z0t*std::sin(thetat) );
1716
1717 m_rd0res[m_rphires.size()]->Fill( resfiller[1], a0r-a0t );
1718 m_rd0res[m_rphires.size()+1]->Fill( fabs(resfiller[1]), a0r-a0t ); //
1719
1720 m_rzedresPull[m_rphires.size()-3]->Fill( resfiller[1], (z0r - z0t) / std::sqrt( (dz0t*dz0t) + (dz0r*dz0r) ) );
1721 m_rzedresPull[m_rphires.size()-2]->Fill( fabs(resfiller[1]), (z0r - z0t) / std::sqrt( (dz0t*dz0t) + (dz0r*dz0r) ) );
1722 m_rd0resPull[m_rphires.size()-3]->Fill( resfiller[1], (a0r - a0t) / std::sqrt( (da0t*da0t) + (da0r*da0r) ) );
1723 m_rd0resPull[m_rphires.size()-2]->Fill( fabs(resfiller[1]), (a0r - a0t) / std::sqrt( (da0t*da0t) + (da0r*da0r) ) );
1724
1725 //-----
1726
1727 m_rnpix_eta_rec->Fill( etat, npixr*1.0 );
1728 m_rnsct_eta_rec->Fill( etat, nsctr*1.0 );
1729 m_rntrt_eta_rec->Fill( etat, nstrawr*1.0 );
1730 m_rnsihit_eta_rec->Fill( etat, npixr*0.5 + nsctr*1.0);
1731
1732 m_rnpix_phi_rec->Fill( phit, npixr*1.0 );
1733 m_rnsct_phi_rec->Fill( phit, nsctr*1.0 );
1734 m_rntrt_phi_rec->Fill( phit, nstrawr*1.0 );
1735
1736 m_rnpix_pt_rec->Fill( std::fabs(pTt), npixr*1.0 );
1737 m_rnsct_pt_rec->Fill( std::fabs(pTt), nsctr*1.0 );
1738 m_rntrt_pt_rec->Fill( std::fabs(pTt), nstrawr*1.0 );
1739
1740 m_rnpixh_pt_rec->Fill( std::fabs(pTt), npixhr*0.5 );
1741 m_rnscth_pt_rec->Fill( std::fabs(pTt), nscthr*1.0 );
1742
1743 m_rnpix_d0_rec->Fill( a0t, npixr*1.0 );
1744 m_rnsct_d0_rec->Fill( a0t, nsctr*1.0 );
1745 m_rntrt_d0_rec->Fill( a0t, nstrawr*1.0 );
1746
1749 m_n_vtx_tracks->Fill( Nvtxtracks );
1750
1751 m_eff_vs_nvtx->Fill( NvtxCount );
1752 m_n_vtx->Fill( NvtxCount );
1753
1754 m_eff_vs_mu->Fill( mu_val );
1755
1756 double vres[6] = { Delta_ipt, etar-etat, phir-phit, z0r-z0t, d0r-d0t, a0r-a0t };
1757 for ( int it=0 ; it<6 ; it++ ) {
1758 if ( it==0 ) {
1759 find("ipT_res")->Fill( vres[0] );
1760 find("spT_res")->Fill( 1.0/pTr-1.0/pTt );
1761 }
1762 if ( TH1F* hptr = find(varName[it]+"_res") ) hptr->Fill( vres[it] );
1763 else std::cerr << "hmmm histo " << varName[it]+"_res" << " not found" << std::endl;
1764 }
1765 //2D plot
1766 if ( TH2F* hptr = find2D("eta_phi_rec") ) {
1767 hptr->Fill( etar,phir );
1768 hptr->GetXaxis()->SetTitle("#eta");
1769 hptr->GetYaxis()->SetTitle("#phi");
1770 //hptr->SetFillStyle("COLZ");
1771 }
1772 if ( TH2F* hptr = find2D("phi_d0_Rec") ) {
1773 hptr->Fill( phir,d0r );
1774 hptr->GetXaxis()->SetTitle("#phi");
1775 hptr->GetYaxis()->SetTitle("d_{0} [mm]");
1776 //hptr->SetFillStyle("COLZ");
1777 }
1778
1779 // raw matched test track errors
1780 if ( TH1F* hptr = find("dpT_rec") ) hptr->Fill(dpTr);
1781 if ( TH1F* hptr =find("deta_rec")) hptr->Fill(detar);
1782 if ( TH1F* hptr =find("dphi_rec")) hptr->Fill(dphir);
1783 if ( TH1F* hptr =find("dz0_rec")) hptr->Fill(dz0r);
1784 if ( TH1F* hptr =find("dd0_rec")) hptr->Fill(dd0r);
1785 if ( TH1F* hptr =find("da0_rec")) hptr->Fill(da0r);
1786
1787 // raw matched reference track errors
1788 if ( TH1F* hptr = find("dpT") ) hptr->Fill(dpTt);
1789 if ( TH1F* hptr = find("deta")) hptr->Fill(detat);
1790 if ( TH1F* hptr = find("dphi")) hptr->Fill(dphit);
1791 if ( TH1F* hptr = find("dz0")) hptr->Fill(dz0t);
1792 if ( TH1F* hptr = find("dd0")) hptr->Fill(dd0t);
1793 if ( TH1F* hptr = find("da0")) hptr->Fill(da0t);
1794
1795
1796 if ( TH1F* hptr = find("dd0_res")) hptr->Fill(dd0r-dd0t);
1797 if ( TH1F* hptr = find("da0_res")) hptr->Fill(da0r-da0t);
1798 if ( TH1F* hptr = find("dz0_res")) hptr->Fill(dz0r-dz0t);
1799
1800 double Dd0 = dd0r-dd0t;
1801 double Da0 = da0r-da0t;
1802
1803 double Ddphi = dphir - dphit;
1804
1805 m_hphivsDd0res[2]->Fill( phit );
1806
1807 // if ( matchedreco->bLayerHits()<=3 ) std::cout << "\nov2\t" << Dd0 << " " << *reftracks[i] << std::endl;
1808 // else std::cout << "\nov4\t" << Dd0 << " " << *reftracks[i] << std::endl;
1809
1810
1811
1812 if ( PRINT_BRESIDUALS ) {
1813 if ( matchedreco->bLayerHits()<=3 ) std::cout << "\nov2\t" << Dd0 << " " << *matchedreco << std::endl;
1814 else std::cout << "\nov4\t" << Dd0 << " " << *matchedreco << std::endl;
1815 }
1816
1817 if ( std::fabs(Dd0)<0.01 ) {
1818 m_hphivsDd0res[0]->Fill( phit );
1819
1820 if ( PRINT_BRESIDUALS ) {
1821 std::cout << "close residual " << Dd0 << " " << Ddphi
1822 << " "<< reftracks[i]->bLayerHits()-matchedreco->bLayerHits()
1823 << " "<< reftracks[i]->pixelHits()-matchedreco->pixelHits();
1824 std::cout << "\nccr\t" << Dd0 << " " << Ddphi << " " << *reftracks[i];
1825 std::cout << "\ncct\t" << Dd0 << " " << Ddphi << " " << *matchedreco << std::endl;
1826 }
1827 }
1828 else {
1829 m_hphivsDd0res[1]->Fill( phit );
1830 if ( PRINT_BRESIDUALS ) {
1831 std::cout << "far residual " << Dd0 << " " << Ddphi
1832 << " "<< reftracks[i]->bLayerHits()-matchedreco->bLayerHits()
1833 << " "<< reftracks[i]->pixelHits()-matchedreco->pixelHits();
1834 std::cout << "\nffr\t" << Dd0 << " " << Ddphi << " " << *reftracks[i];
1835 std::cout << "\nfft\t" << Dd0 << " " << Ddphi << " " << *matchedreco << std::endl;
1836 }
1837 }
1838
1839
1840 m_hphivsDa0res[2]->Fill( iphi );
1841 if ( std::fabs(Da0)<0.01 ) m_hphivsDa0res[0]->Fill( iphi );
1842 else m_hphivsDa0res[1]->Fill( iphi );
1843
1844
1845 // pull stats
1846 double pull_pt = Delta_pt / std::sqrt( (dpTt*dpTt) + (dpTr*dpTr) );
1847 double pull_eta = (etar - etat) / std::sqrt( (detat*detat) + (detar*detar) );
1848 double pull_phi = (phir - phit) / std::sqrt( (dphit*dphit) + (dphir*dphir) );
1849 double pull_z0 = (z0r - z0t) / std::sqrt( (dz0t*dz0t) + (dz0r*dz0r) );
1850 double pull_d0 = (d0r - d0t) / std::sqrt( (dd0t*dd0t) + (dd0r*dd0r) );
1851 double pull_a0 = (a0r - a0t) / std::sqrt( (da0t*da0t) + (da0r*da0r) );
1852
1853 if ( TH1F* hptr = find("pT_pull") ) hptr->Fill(pull_pt);
1854 if ( TH1F* hptr = find("eta_pull")) hptr->Fill(pull_eta);
1855 if ( TH1F* hptr = find("phi_pull")) hptr->Fill(pull_phi);
1856 if ( TH1F* hptr = find("z0_pull")) hptr->Fill(pull_z0);
1857 if ( TH1F* hptr = find("d0_pull")) hptr->Fill(pull_d0);
1858 if ( TH1F* hptr = find("a0_pull")) hptr->Fill(pull_a0);
1859
1860 // pull stats - SIMPLE VERSION
1861 double pull_pt_simp = Delta_pt / sqrt( dpTr*dpTr );
1862 double pull_eta_simp = (etar - etat) / sqrt( detar*detar );
1863 double pull_phi_simp = (phir - phit) / sqrt( dphir*dphir );
1864 double pull_z0_simp = (z0r - z0t) / sqrt( dz0r*dz0r );
1865 double pull_d0_simp = (d0r - d0t) / sqrt( dd0r*dd0r );
1866 double pull_a0_simp = (a0r - a0t) / sqrt( da0r*da0r );
1867
1868 if ( TH1F* hptr = find("pT_pull_simple") ) hptr->Fill(pull_pt_simp);
1869 if ( TH1F* hptr = find("eta_pull_simple")) hptr->Fill(pull_eta_simp);
1870 if ( TH1F* hptr = find("phi_pull_simple")) hptr->Fill(pull_phi_simp);
1871 if ( TH1F* hptr = find("z0_pull_simple")) hptr->Fill(pull_z0_simp);
1872 if ( TH1F* hptr = find("d0_pull_simple")) hptr->Fill(pull_d0_simp);
1873 if ( TH1F* hptr = find("a0_pull_simple")) hptr->Fill(pull_a0_simp);
1874
1875
1876 if ( TH1F* hptr = find("etai_res") ) hptr->Fill( etat-etar );
1877
1878
1879 double Delphi = phit-phir;
1880 double Deleta = etat-etar;
1881
1882 if ( Delphi<-M_PI ) Delphi+=2*M_PI;
1883 if ( Delphi>M_PI ) Delphi -=2*M_PI;
1884
1885 double DeltaR = std::sqrt(Delphi*Delphi+Deleta*Deleta);
1886
1887 m_hDeltaR->Fill(DeltaR);
1888
1889 m_deltaR_v_eta->Fill(etat, DeltaR);
1890 m_deltaR_v_pt->Fill(std::fabs(pTt), DeltaR);
1891
1892 // in this loop over the reference tracks, could fill efficiency
1893 // histograms
1894
1895 // m_eff_vs_lb->Fill( m_rmap[r]+lb );
1896
1897 if ( TH1F* hptr = find("nsi_matched")) hptr->Fill(nsir);
1898
1900
1901
1902 m_rChi2prob_rec->Fill( std::fabs(pTr), TMath::Prob(matchedreco->chi2(),matchedreco->dof()) );
1903 m_rChi2_rec->Fill( std::fabs(pTr), matchedreco->chi2() );
1904 m_rChi2dof_rec->Fill( std::fabs(pTr), matchedreco->chi2()/matchedreco->dof() );
1905
1906 m_rChi2d_vs_Chi2d->Fill( reftracks[i]->chi2()/reftracks[i]->dof(),
1907 matchedreco->chi2()/matchedreco->dof() );
1908
1909 m_rDChi2dof->Fill( reftracks[i]->chi2()/reftracks[i]->dof(),
1910 (matchedreco->chi2()/matchedreco->dof())-(reftracks[i]->chi2()/reftracks[i]->dof()) );
1911
1912 }
1913 else {
1914
1916
1917 m_rnpix_pt_bad->Fill( std::fabs(pTt), npixt*0.5 );
1918 m_rnsct_pt_bad->Fill( std::fabs(pTt), nsctt*1.0 );
1919 m_rntrt_pt_bad->Fill( std::fabs(pTt), nstrawt*1.0 );
1920
1921 m_rChi2prob_bad->Fill( std::fabs(pTt), TMath::Prob(reftracks[i]->chi2(),reftracks[i]->dof()) );
1922 m_rChi2_bad->Fill( std::fabs(pTt), reftracks[i]->chi2() );
1923 m_rChi2dof_bad->Fill( std::fabs(pTt), reftracks[i]->chi2()/reftracks[i]->dof() );
1924
1925
1926 // fill efficiencies with unmatched histos
1927 // std::cout << "NULL" << std::endl;
1928 m_eff_pt->FillDenom(std::fabs(pTt));
1929 m_eff_z0->FillDenom(z0t);
1930 m_eff_eta->FillDenom(etat);
1931 m_eff_phi->FillDenom(phit);
1932 m_eff_d0->FillDenom(d0t);
1933 m_eff_a0->FillDenom(a0t);
1934
1935 // signed pT
1936 if ( pTt<0 ) m_eff_ptm->FillDenom(std::fabs(pTt));
1937 else m_eff_ptp->FillDenom(std::fabs(pTt));
1938
1939 m_eff_roi_deta->FillDenom(droi_detat);
1940 m_eff_roi_dphi->FillDenom(droi_dphit);
1941 m_eff_roi_dR->FillDenom(droi_dRt);
1942
1943 m_eff_vs_mult->FillDenom( m_Nref );
1944
1945 m_eff_eta_vs_pt->FillDenom( std::fabs(pTt), etat );
1946 m_eff_d0_vs_pt->FillDenom( std::fabs(pTt), a0t );
1947
1948 dump = false;
1949
1950 m_eff_vs_ntracks->FillDenom( Nvtxtracks );
1951 m_eff_vs_ntracks2->FillDenom( Nvtxtracks );
1952 m_n_vtx_tracks->Fill( Nvtxtracks );
1953
1954
1955 m_eff_vs_nvtx->FillDenom( NvtxCount );
1956 m_n_vtx->Fill( NvtxCount );
1957
1958 double mu_val = gevent->mu();
1959
1960 m_eff_vs_mu->FillDenom(mu_val);
1961
1962 if ( tobj ) {
1963 m_eff_vs_etovpt->FillDenom(etovpt_val);
1964 m_eff_vs_et->FillDenom( std::fabs(tobj->pt()*0.001) );
1965 }
1966
1967 if ( dumpflag ) {
1968 std::ostream& dumpstream = dumpfile;
1969
1970 if ( std::fabs(pTt)>1 ) {
1971 dump = true;
1972
1973 hipt = true;
1974 dumpstream << m_name << "\tMISSING TRACK run " << r << "\tevent " << ev
1975 << "\tlb " << lb << "\tN vertices " << NvtxCount << std::endl;
1976 dumpstream << m_name << "\tMISSING TRACK RoI " << *groi << std::endl;
1977 dumpstream << m_name << "\tMISSING TRACK Track " << *reftracks[i];
1978 if ( std::fabs(pTt)>=30 ) dumpstream << "\tvery high pt";
1979 if ( std::fabs(pTt)>4 &&
1980 std::fabs(pTt)<30 ) dumpstream << "\t high pt";
1981 dumpstream << std::endl;
1982
1983 if ( std::fabs(pTt)>=20 ){
1984 dumpstream << "Test tracks " << std::endl;
1985 for (unsigned int ii=0; ii<testtracks.size(); ii++){
1986 dumpstream << *testtracks[ii] << std::endl;
1987 }
1988 }
1989 }
1990 }
1991
1992
1993 // m_eff_vs_lb->FillDenom( ts );
1994 m_eff_vs_lb->FillDenom( gevent->lumi_block() );
1995 }
1996
1997 }
1998
1999 // return;
2000
2001
2002 // for fake/purity histograms, loop over the test tracks
2003 // and get the corresponding matched reference tracks from the
2004 // reverse map in the TrackAssociator class - revmatched()
2005
2006 static int icount = 0;
2007
2008 // if ( icount%1000 ) std::cout << "chain " << name() << "\t " << m_Nreco << " tracks" << std::endl;
2009 // if ( icount%1000 )
2010 if ( m_print ) std::cout << "ConfAnalysis::execute() \t " << name() << "\t " << icount << " events\t " << testtracks.size() << " tracks (" << m_Nreco << ")" << "\n---------------" << std::endl;
2011
2012 icount++;
2013
2014 for ( int i=testtracks.size() ; i-- ; ) {
2015
2016 // std::cout << "\t\tConfAnalysis purity " << name() << "\t" << i << " " << *testtracks[i] << " -> ";
2017
2018 // double pTr = std::fabs(testtracks[i]->pT());
2019 double pTr = testtracks[i]->pT()/1000;
2020 double etar = testtracks[i]->eta();
2021 double phir = testtracks[i]->phi();
2022 double thetar = 2*std::atan( exp( (-1)*etar) );
2023
2024 double z0r = testtracks[i]->z0(); // + ((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
2025 double d0r = testtracks[i]->a0() - sin(phir)*m_xBeamTest + cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
2026 double a0r = testtracks[i]->a0();
2027 // double a0rp = testtracks[i]->a0() - sin(phir)*m_xBeam - cos(phir)*m_yBeam; // this will be changed when we know the beam spot position
2028
2029 if ( m_xBeamTest!=0 && m_yBeamTest!=0 ) {
2030 d0r = testtracks[i]->a0();
2031 a0r = testtracks[i]->a0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
2032 z0r = testtracks[i]->z0()+((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
2033 }
2034
2035
2036 // std::cout << "d0 " << d0r << "\tphi " << phir << "\tx " << m_xBeamTest << "\ty " << m_yBeamTest << std::endl;
2037
2038 double nsctr = testtracks[i]->sctHits();
2039 double npixr = testtracks[i]->pixelHits()*0.5;
2040 double nsir = testtracks[i]->pixelHits()*0.5 + testtracks[i]->sctHits();
2041
2042 double ntrtr = testtracks[i]->trHits();
2043 double nstrawr = testtracks[i]->strawHits();
2044
2045
2046 m_rnsct_vs_npix_rec->Fill( npixr, nsctr );
2047
2048
2049#if 0
2050 double dpTr_b = testtracks[i]->dpT()/1000;
2051 double detar_b = testtracks[i]->deta();
2052 double dphir_b = testtracks[i]->dphi();
2053 double dz0r_b = testtracks[i]->dz0(); // + ((std::cos(phir)*m_xBeamTest + std::sin(phir)*m_yBeamTest)/std::tan(thetar));
2054 double dd0r_b = testtracks[i]->da0() - sin(phir)*m_xBeamTest + cos(phir)*m_yBeamTest;
2055 double da0r_b = testtracks[i]->da0();
2056
2057
2058 std::cout << "pTr_b = " << pTr << " +/- " << dpTr_b << std::endl;
2059 std::cout << "phir_b = " << phir << " +/- " << dphir_b << std::endl;
2060 std::cout << "z0r_b = " << z0r << " +/- " << dz0r_b << std::endl;
2061 std::cout << "d0r_b = " << d0r << " +/- " << dd0r_b << std::endl;
2062 std::cout << "a0r_b = " << a0r << " +/- " << da0r_b << std::endl;
2063 std::cout << "etar_b = " << etar << " +/- " << detar_b << std::endl;
2064#endif
2065
2066 // double ts_scale = (ts-1260400000)*3000.0/(1260700000-1260400000);
2067
2068 // m_z_vs_lb->Fill( m_rmap[r]+lb, z0r );
2069 // m_z_vs_lb->Fill( ts, z0r );
2070 m_z_vs_lb->Fill( gevent->lumi_block(), z0r );
2071
2072 // hnpix_v_sct_rec->Fill( nsctr*0.5, npixr*1.0 );
2073
2074 if ( m_h2r ) m_h2r->Fill( phir, d0r );
2075 if ( m_h2a0r ) m_h2a0r->Fill( phir, a0r );
2076
2077 const TIDA::Track* matchedref = matcher->revmatched(testtracks[i]);
2078
2079 // if ( matchedref ) std::cout << *matchedref << std::endl;
2080 // else std::cout << "NULL" << std::endl;
2081
2082#if 1
2083 // raw test track distributions
2084 double vpart[11] = { std::fabs(pTr), etar, phir, z0r, d0r, a0r, nsctr, npixr, nsir, ntrtr, nstrawr };
2085 for ( int it=0 ; it<11 ; it++ ) {
2086 // std::string hname = name()+"_"+varName[it]+"_rec";
2087 // std::string hname = varName[it]+"_rec";
2088 // std::map<std::string, TH1F*>::iterator hmitr = m_histos.find(hname);
2089 // if ( hmitr!=m_histos.end() ) hmitr->second->Fill( vpar[it] );
2090 // else std::cerr << "hmmm histo " << hname << " not found" << std::endl;
2091 if ( TH1F* hptr = find(varName[it]+"_rec") ) hptr->Fill( vpart[it] );
2092 else std::cerr << "hmmm histo " << varName[it]+"_rec" << " not found" << std::endl;
2093 }
2094 //2D plot
2095 if ( TH2F* hptr = find2D("eta_phi_rec") ) {
2096 hptr->Fill( etar,phir );
2097 hptr->GetXaxis()->SetTitle("#eta");
2098 hptr->GetYaxis()->SetTitle("#phi");
2099 //hptr->SetFillStyle("COLZ");
2100 }
2101 if ( TH2F* hptr = find2D("phi_d0_rec") ) {
2102 hptr->Fill( phir,d0r );
2103 hptr->GetXaxis()->SetTitle("#phi");
2104 hptr->GetYaxis()->SetTitle("d_{0} [mm]");
2105 //hptr->SetFillStyle("COLZ");
2106 }
2107#endif
2108
2109
2110 // purities
2111 if ( matchedref ) {
2112
2113 // std::cout << *matchedref << std::endl;
2114
2115 m_purity_pt->Fill(std::fabs(pTr));
2116 m_purity_z0->Fill(z0r);
2117 m_purity_eta->Fill(etar);
2118 m_purity_phi->Fill(phir);
2119 m_purity_d0->Fill(d0r);
2120 m_purity_a0->Fill(a0r);
2121
2122 // hnpix_v_sct_match->Fill( nsctr*0.5, npixr*0.5 );
2123
2124 }
2125 else {
2126 // std::cout << "NULL" << std::endl;
2127 m_purity_pt->FillDenom(std::fabs(pTr));
2128 m_purity_z0->FillDenom(z0r);
2129 m_purity_eta->FillDenom(etar);
2130 m_purity_phi->FillDenom(phir);
2131 m_purity_d0->FillDenom(d0r);
2132 m_purity_a0->FillDenom(a0r);
2133 }
2134
2135 }
2136
2137 if ( dump && m_print ) {
2138
2139 std::cout << "ConfAnalysis::execute() missed a high pT track - dumping tracks" << std::endl;
2140
2141 for ( int i=reftracks.size() ; i-- ; ) {
2142
2143 if ( std::fabs( reftracks[i]->pT() ) > 1000 ) {
2144 std::cout << "\t dump " << *reftracks[i];
2145 const TIDA::Track* matchedreco = matcher->matched(reftracks[i]);
2146 if ( matchedreco ) std::cout << " <--> " << *matchedreco << std::endl;
2147 else std::cout << std::endl;
2148 }
2149
2150 }
2151
2152 for ( int i=testtracks.size() ; i-- ; ) {
2153 const TIDA::Track* matchedref = matcher->revmatched(testtracks[i]);
2154 if ( matchedref==0 ) std::cout << "\t\t\t\t\t " << *testtracks[i] << std::endl;
2155 }
2156
2157 }
2158
2159 if ( m_print ) std::cout << "ConfAnalysis::execute() exiting" << std::endl;
2160
2161}
2162
2163
2164
2165
2166
#define M_PI
Scalar phi() const
phi method
scales to modify the binning for histograms
int Nvtxtracks
Definition globals.cxx:17
void fitSin(TH1D *h, const std::string &parent="")
BinConfig cosmicBinConfig("cosmic")
TF1 * FitFWGaussian(TH1D *s, double a, double b)
BinConfig electronBinConfig("electron")
double wrapphi(double phi)
fill all the histograms - matched histograms, efficiencies etc
int NvtxCount
Definition globals.cxx:18
BinConfig muonBinConfig("muon")
bool PRINT_BRESIDUALS
stack trace headers
BinConfig g_binConfig("standard")
BinConfig bjetBinConfig("bjet")
std::ofstream dumpfile("dumpfile.log")
BinConfig tauBinConfig("tau")
void Normalise(TH1 *h)
TIDA::Event * gevent
Definition globals.cxx:13
JetDumper::Name Name
Definition JetDumper.cxx:19
static Double_t a
Basic event class to contain a vector of chains for trigger analysis.
TIDA::Associator< TIDA::Track > TrackAssociator
base class for a single track selection filter allowing parameter setting for complex track selection
#define y
#define x
Header file for AthHistogramAlgorithm.
Efficiency1D * m_eff_roi_deta
TH1F * m_n_vtx_tracks
Resplot * m_rnsih_pt
Resplot * m_rChi2prob_rec
Resplot * m_h2r
Efficiency1D * m_eff_ptp
Efficiency1D * m_eff_phi
Resplot * m_deltaR_v_eta
Efficiency1D * m_purity_phi
Resplot * m_rnpix_d0
Resplot * m_rnscth_pt
Resplot * m_rRoi_dzed_vs_eta
Resplot * m_rntrt_d0
Resplot * m_rnpix_eta_rec
Resplot * m_rnsih_eta
Resplot * m_rnsct_pt_rec
Efficiency1D * m_eff_vs_ntracks2
Resplot * m_rnsct_d0_rec
Resplot * m_rntrt_eta_rec
Efficiency1D * m_eff_vs_mu
Resplot * m_rnscth_pt_rec
Efficiency1D * m_eff_a0
Efficiency1D * m_purity_a0
Resplot * m_rnpixh_pt_rec
TH1F * m_etovpt_raw
electron specific ET/PT related stuff
Efficiency1D * m_eff_vs_lb
Resplot * m_rnsihit_eta
Resplot * m_rzedreslb
Resplot * m_h2a0
Efficiency1D * m_purity_pt
Resplot * m_rnsct_phi_rec
Efficiency1D * m_eff_vs_nvtx
Efficiency2D * m_eff_eta_vs_pt
std::vector< Resplot * > m_retaresPull
Resplot * m_h2m
int m_Nreco
number of reconstructed tracks
std::vector< Resplot * > m_rd0resPull
Resplot * m_deltaR_v_pt
TH1F * m_hphivsDd0res[3]
Resplot * m_rChi2_rec
void addHistogram(TH1F *h)
Resplot * m_rnpix_lb
std::vector< Resplot * > m_rptres
Resplot * m_rnpix_pt_bad
Resplot * m_rnpix_d0_rec
TH1F * find(const std::string &n)
Resplot * m_rntrt_pt_bad
Efficiency1D * m_eff_vs_ntracks
bool m_initialiseFirstEvent
Efficiency1D * m_eff_vs_etovpt
virtual void execute(const std::vector< TIDA::Track * > &reftracks, const std::vector< TIDA::Track * > &testtracks, TrackAssociator *matcher, TrigObjectMatcher *objects)
Resplot * m_rnblh_d0
Resplot * m_rnpix_eta
Efficiency1D * m_eff_eta
Resplot * m_rChi2dof
TagNProbe * m_TnP_tool
Resplot * m_rnsct_eta
Resplot * m_rChi2
std::vector< Resplot * > m_rDz0res
Resplot * m_rnbl_d0
TH1F * m_invmassObj
std::vector< Resplot * > m_rzedthetares
Resplot * m_rnpix_phi_rec
Efficiency1D * m_eff_vs_mult
virtual void initialise()
book all the histograms
Resplot * m_rChi2prob
Resplot * m_rnsi_eta
Resplot * m_rnsihit_eta_rec
Resplot * m_rnsct_lb_rec
Efficiency1D * m_eff_roi_dphi
Resplot * m_rnsct_phi
Efficiency1D * m_eff_roi_dR
Resplot * m_rntrt_phi_rec
std::vector< Resplot * > m_rptresPull
std::vector< Resplot * > m_res
Resplot * m_rnpix_lb_rec
Resplot * m_rd0_vs_phi
beam spot dependent
Resplot * m_rDChi2dof
Resplot * m_rntrt_eta
Resplot * m_rnsct_vs_npix_rec
void addHistogram2D(TH2F *h)
Resplot * m_rChi2dof_rec
TH1F * m_hphivsDa0res[3]
Efficiency1D * m_purity_d0
Efficiency1D * m_eff_pt
std::vector< Resplot * > m_rd0res
Resplot * m_rnsct_pt_bad
Efficiency1D * m_eff_vs_et
Resplot * m_rnsct_eta_rec
Resplot * m_rChi2prob_bad
std::vector< Resplot * > m_retares
Resplot * m_rnsct_lb
Resplot * m_rnpixh_d0
Resplot * m_rnsi_pt
Efficiency1D * m_eff_z0
Resplot * m_rnsih_d0
Resplot * m_rChi2d_vs_Chi2d
Resplot * m_rntrt_pt
Efficiency2D * m_eff_d0_vs_pt
Resplot * m_rd0_vs_phi_rec
virtual void finalise()
calculate the efficiencies and write them out with all the histograms
std::vector< Resplot * > m_riptres
Resplot * m_z_vs_lb
Resplot * m_rntrt_d0_rec
Resplot * m_rnscth_d0
Resplot * m_rnpixh_pt
Efficiency1D * m_purity_eta
std::vector< Resplot * > m_rzedresPull
Resplot * m_rnpix_phi
Resplot * m_rChi2_bad
Resplot * m_rnpix_pt_rec
Resplot * m_rnsct_d0
std::vector< Resplot * > m_rphires
Resplot * m_rChi2dof_bad
std::vector< Resplot * > m_rphiresPull
Resplot * m_h2a0r
virtual void initialiseInternal()
Resplot * m_h2
Resplot * m_rzedlb_rec
Resplot * m_rntrt_phi
Resplot * m_rnsi_d0
std::map< std::string, TH1F * > m_histos
bool m_print
flag to print out the matched tracks etc
Resplot * m_rzedlb
Resplot * m_rnpix_pt
Resplot * m_rntrt_pt_rec
TH2F * find2D(const std::string &n)
Resplot * m_rRoi_dphi_vs_eta
Resplot * m_rRoi_deta_vs_eta
Residuals.
TIDDirectory * m_dir
std::vector< Resplot * > m_rzedres
Efficiency1D * m_eff_d0
Resplot * m_rnsct_pt
std::vector< Resplot * > m_rDa0res
Efficiency1D * m_purity_z0
std::vector< Resplot * > m_rDd0res
Efficiency1D * m_eff_ptm
Resplot * m_rnsct_vs_npix
TGraphAsymmErrors * Bayes(double scale=100)
evaluate the uncertainties correctly ...
static TF1 * FitNull(TH1D *s, double a=-999, double b=-999)
Definition Resplot.cxx:1267
static TF1 * FitNull95(TH1D *s, double a=0, double b=0)
Definition Resplot.cxx:1652
virtual const T * revmatched(S *t)
virtual const S * matched(T *t)
int find(T *&t, const std::string &key)
void finalise(double scale=100)
actually calculate the efficiencies
TIDA::FeatureStore & store()
double m_yBeamReference
double m_xBeamReference
beamline positions reference sample
std::string m_name
identifier of the of the analysis - also used for the root directory into which the histograms are pu...
const std::string & name() const
return identifier
double m_xBeamTest
test sample
double pt() const
double chi2(TH1 *h0, TH1 *h1)
int ir
counter of the current depth
Definition fastadd.cxx:49
bool dumpflag
Definition globals.cxx:30
TIDARoiDescriptor * groi
all these externals initialised in globals.cxx
Definition globals.cxx:14
int NMod
Definition globals.cxx:20
int lb
Definition globals.cxx:23
bool hipt
Definition globals.cxx:29
int r
Definition globals.cxx:22
double a0
Definition globals.cxx:27
int ev
Definition globals.cxx:25
-event-from-file
double a0_NScale
Definition BinConfig.h:73
double phires_NScale
Definition BinConfig.h:79
double eta_NScale
Definition BinConfig.h:70
double a0Max
Definition BinConfig.h:85
double phi_NScale
Definition BinConfig.h:71
double d0Max
Definition BinConfig.h:84
double z0_NScale
Definition BinConfig.h:74
double z0Max
Definition BinConfig.h:86
double d0_NScale
Definition BinConfig.h:72
double ptres_NScale
scales for the residuals
Definition BinConfig.h:77
double pt_NScale
Definition BinConfig.h:69