ATLAS Offline Software
Loading...
Searching...
No Matches
SigAnalysis.cxx
Go to the documentation of this file.
1
9
10
11#include "TrigInDetAnalysisExample/SigAnalysis.h"
12
13void SigAnalysis::initialise() {
14
15 std::cout << "SigAnalysis::initialise() " << name() << std::endl;
16
17 //+++ pT ranges
18 // double tmp_maxPt = 50000.;
19 double tmp_absResPt = 0.0005;
20
21 const int pTResBins = 100;
22
23 //+++ Eta ranges
24 double tmp_maxEta = 3.;
25 double tmp_absResEta = 0.04; // 0.0005;
26
27 //+++ Phi ranges
28 double tmp_maxPhi = 3.142;
29 double tmp_absResPhi = 0.02; // 0.0001;
30
31
32 const int etaBins = 60;
33 const int etaResBins = 300;
34
35 const int phiBins = 36;
36 const int phiResBins = 100;
37
38 const int zBins = 50;
39 const double zMax = 400;
40
41 const int zresBins = 100;
42 const double zresMax = 10;
43
44 const int d0Bins = 100;
45 const double d0Max = 10;
46
47 const int d0resBins = 100;
48 const double d0resMax = 10;
49
50 // beamspot corrected position
51
52 const int a0Bins = 100;
53 const double a0Max = 10;
54
55 const int a0resBins = 100;
56 const double a0resMax = 5;
57
58 //+++ Book histograms
59
60 // calculate a logarithmic binning in pt
61
62 const int ptnbins = 20;
63 // const int ptnbins = 30;
64 double ptbinlims[ptnbins+1];
65 for ( int i=0 ; i<=ptnbins ; i++ ) { ptbinlims[i] = std::pow(10, 2.0*i/ptnbins+2); }
66 // for ( int i=0 ; i<=ptnbins ; i++ ) { ptbinlims[i] = std::pow(10, 2.3*i/ptnbins+2); }
67
68
69
70 TDirectory* dir = gDirectory;
71
72 std::cout << "SigAnalysis::initialize() Directory " << gDirectory->GetName() << " " << name() << std::endl;
73
74 m_dir = new TIDDirectory(name());
75 m_dir->push();
76
77 std::cout << "SigAnalysis::initialize() Directory " << gDirectory->GetName() << " package directory, " << name() << std::endl;
78
79 Efficiency1D* heff[8];
80 Efficiency1D* hpurity[6];
81
82 // "reference" quantities
83 addHistogram( new TH1F( "pT", "pT", ptnbins, ptbinlims ) );
84 addHistogram( new TH1F( "eta", "eta", etaBins, -tmp_maxEta, tmp_maxEta ) );
85 addHistogram( new TH1F( "phi", "phi", phiBins, -tmp_maxPhi, tmp_maxPhi ) );
86 addHistogram( new TH1F( "z0", "z0", zBins, -zMax, zMax ) );
87 addHistogram( new TH1F( "d0", "d0", d0Bins, -d0Max, d0Max ) );
88 addHistogram( new TH1F( "a0", "a0", a0Bins, -a0Max, a0Max ) );
89
90 // efficienies and purities
91 heff[0] = new Efficiency1D( find("pT"), "pT_eff" );
92 heff[1] = new Efficiency1D( find("eta"), "eta_eff" );
93 heff[2] = new Efficiency1D( find("phi"), "phi_eff" );
94 heff[3] = new Efficiency1D( find("z0"), "z0_eff" );
95 heff[4] = new Efficiency1D( find("d0"), "d0_eff" );
96 heff[5] = new Efficiency1D( find("a0"), "a0_eff" );
97
98 heff[6] = new Efficiency1D( find("pT"), "pTm_eff" );
99 heff[7] = new Efficiency1D( find("pT"), "pTp_eff" );
100
101 m_eff_pt = heff[0];
102 m_eff_eta = heff[1];
103 m_eff_phi = heff[2];
104 m_eff_z0 = heff[3];
105 m_eff_d0 = heff[4];
106 m_eff_a0 = heff[5];
107
108 m_eff_ptm = heff[6];
109 m_eff_ptp = heff[7];
110
111 hpurity[0] = new Efficiency1D( find("pT"), "pT_pur" );
112 hpurity[1] = new Efficiency1D( find("eta"), "eta_pur" );
113 hpurity[2] = new Efficiency1D( find("phi"), "phi_pur" );
114 hpurity[3] = new Efficiency1D( find("z0"), "z0_pur" );
115 hpurity[4] = new Efficiency1D( find("d0"), "d0_pur" );
116 hpurity[5] = new Efficiency1D( find("a0"), "a0_pur" );
117
118 m_purity_pt = hpurity[0];
119 m_purity_eta = hpurity[1];
120 m_purity_phi = hpurity[2];
121 m_purity_z0 = hpurity[3];
122 m_purity_d0 = hpurity[4];
123 m_purity_a0 = hpurity[5];
124
125 // "test" quantities
126 addHistogram( new TH1F( "pT_rec", "pT_rec", ptnbins, ptbinlims ) );
127 addHistogram( new TH1F( "eta_rec", "eta_rec", etaBins, -tmp_maxEta, tmp_maxEta ) );
128 addHistogram( new TH1F( "phi_rec", "phi_rec", phiBins, -tmp_maxPhi, tmp_maxPhi ) );
129 addHistogram( new TH1F( "z0_rec", "z0_rec", zBins, -zMax, zMax ) );
130 addHistogram( new TH1F( "d0_rec", "d0_rec", d0Bins, -d0Max, d0Max ) );
131 addHistogram( new TH1F( "a0_rec", "a0_rec", a0Bins, -a0Max, a0Max ) );
132
133 // resolutions
134 addHistogram( new TH1F( "pT_res", "pT_res", pTResBins, -tmp_absResPt, tmp_absResPt ) );
135 addHistogram( new TH1F( "eta_res", "eta_res", etaResBins, -tmp_absResEta, tmp_absResEta ) );
136 addHistogram( new TH1F( "phi_res", "phi_res", phiResBins, -tmp_absResPhi, tmp_absResPhi ) );
137 addHistogram( new TH1F( "z0_res", "z0_res", zresBins, -zresMax, zresMax ) );
138 addHistogram( new TH1F( "d0_res", "d0_res", d0resBins, -0.5*d0resMax, 0.5*d0resMax ) );
139 addHistogram( new TH1F( "a0_res", "a0_res", a0resBins, -0.5*a0resMax, 0.5*a0resMax ) );
140
141 // hit occupancies
142
143 int NHits = 40;
144 int Ntracks = 50;
145
146 addHistogram( new TH1F( "nsct", "nsct", NHits, -0.5, float(NHits-0.5) ) );
147 addHistogram( new TH1F( "nsct_rec", "nsct_rec", NHits, -0.5, float(NHits-0.5) ) );
148
149 addHistogram( new TH1F( "npix", "npix", NHits, -0.5, float(NHits-0.5) ) );
150 addHistogram( new TH1F( "npix_rec", "npix_rec", NHits, -0.5, float(NHits-0.5) ) );
151
152 addHistogram( new TH1F( "ntrt", "ntrt", NHits, -0.5, float(NHits-0.5) ) );
153 addHistogram( new TH1F( "ntrt_rec", "ntrt_rec", NHits, -0.5, float(NHits-0.5) ) );
154
155 addHistogram( new TH1F( "nstraw", "nstraw", NHits*4, -0.5, float(4*NHits-0.5) ) );
156 addHistogram( new TH1F( "nstraw_rec", "nstraw_rec", NHits*4, -0.5, float(4*NHits-0.5) ) );
157
158 addHistogram( new TH1F( "ntracks", "ntracks", Ntracks, -0.5, float(Ntracks+0.5) ) );
159 addHistogram( new TH1F( "ntracks_rec", "ntracks_rec", Ntracks, -0.5, float(Ntracks+0.5) ) );
160
161
162 // beam offset fitting histos
163 m_h2 = new TH2F( "d0vphi", "d0vphi", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
164 m_h2r = new TH2F( "d0vphi_rec", "d0vphi_rec", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
165 m_h2m = new TH2F( "d0vphi_match", "d0vphi_match", phiBins, -3.142, 3.142, d0Bins, -d0Max, d0Max );
166
167 m_dir->pop();
168
169 dir->cd();
170
171 // std::cout << "initialize() Directory " << gDirectory->GetName() << " on leaving" << std::endl;
172
173}
174
175
176
177void SigAnalysis::finalise() {
178 std::cout << "SigAnalysis::finalise() " << name() << "\tNreco " << m_Nreco << " tracks" << std::endl;
179
180 // if ( m_Nreco==0 ) return;
181
182 m_dir->push();
183
184 std::cout << "SigAnalysis::finalise() " << gDirectory->GetName() << std::endl;
185
186 // std::map<std::string, TH1F*>::iterator hitr=m_histos.begin();
187 // std::map<std::string, TH1F*>::iterator hend=m_histos.end();
188 // for ( ; hitr!=hend ; hitr++ ) hitr->second->Write();
189
190 // std::cout << "DBG >" << m_eff_pt->Hist()->GetName() << "< DBG" << std::endl;
191
192 Efficiency1D* heff[8] = { m_eff_pt, m_eff_eta, m_eff_phi, m_eff_z0, m_eff_d0, m_eff_a0, m_eff_ptm, m_eff_ptp };
193 for ( int i=8 ; i-- ; ) { heff[i]->finalise(); } // heff[i]->Hist()->Write(); }
194
195 // std::cout << "DBG >" << m_purity_pt->Hist()->GetName() << "< DBG" << std::endl;
196
197 Efficiency1D* hpurity[6] = { m_purity_pt, m_purity_eta, m_purity_phi, m_purity_z0, m_purity_d0, m_purity_a0 };
198 for ( int i=6 ; i-- ; ) { hpurity[i]->finalise(); } // hpurity[i]->Hist()->Write(); }
199
200 m_dir->pop();
201
202}
203
204
205void SigAnalysis::execute(const std::vector<TIDA::Track*>& reftracks,
206 const std::vector<TIDA::Track*>& testtracks,
207 TrackAssociator* matcher )
208{
209 if ( m_print ) std::cout << "SigAnalysis::execute() \t " << name()
210 << "\tref " << reftracks.size()
211 << "\ttest " << testtracks.size() << std::endl;
212
213 // std::cout << "\tx " << m_xBeamReference << "\ty " << m_yBeamReference
214 // << "\tx " << m_xBeamTest << "\ty " << m_yBeamTest << std::endl;
215
216
217 // std::cout << "SigAnalysis (resolutions really) filling " << std::endl;
218
219 // should have these as a class variable
220 static const std::string varName[10] = { "pT", "eta", "phi", "z0", "d0", "a0", "nsct", "npix", "ntrt", "nstraw" };
221
222 // std::cout << "SigAnalysis ref size " << reftracks.size() << "\ttest size " << testtracks.size() << std::endl;
223
224 std::map<std::string, TH1F*>::iterator hmitr = m_histos.find("ntracks");
225 if ( hmitr!=m_histos.end() ) hmitr->second->Fill( reftracks.size() );
226
227 hmitr = m_histos.find("ntracks_rec");
228 if ( hmitr!=m_histos.end() ) hmitr->second->Fill( testtracks.size() );
229
230 bool dump = false;
231
232 for ( int i=reftracks.size() ; i-- ; ) {
233
234 double pTt = reftracks[i]->pT();
235 double z0t = reftracks[i]->z0();
236 double etat = reftracks[i]->eta();
237 double phit = reftracks[i]->phi();
238 double d0t = reftracks[i]->a0();
239 // this will be changed when we know the beam spot position
240 // double a0t = reftracks[i]->a0() + sin(phit)*m_xBeam - cos(phit)*m_yBeam;
241 double a0t = reftracks[i]->a0() + sin(phit)*m_xBeamReference - cos(phit)*m_yBeamReference;
242
243 double nsctt = reftracks[i]->sctHits();
244 double npixt = reftracks[i]->pixelHits();
245
246 double ntrtt = reftracks[i]->trHits();
247 double nstrawt = reftracks[i]->strawHits();
248
249 // std::cout << "Fill h2 " << " " << m_h2m << " " << *reftracks[i] << std::endl;
250
251 m_h2->Fill( phit, d0t );
252
253 const TIDA::Track* matchedreco = matcher->matched(reftracks[i]);
254
255 // std::cout << "\t\tSigAnalysis " << name() << "\t" << i << " " << *reftracks[i] << " -> ";
256
257 // raw reference track distributions
258 double vpart[10] = { std::fabs(pTt), etat, phit, z0t, d0t, a0t, nsctt, npixt, ntrtt, nstrawt };
259 for ( int it=0 ; it<10 ; it++ ) {
260 // std::string hname = varName[it];
261 // std::map<std::string, TH1F*>::iterator hmitr = m_histos.find(hname);
262 // if ( hmitr!=m_histos.end() ) hmitr->second->Fill( vpart[it] );
263
264 if ( TH1F* hptr = find( varName[it] ) ) hptr->Fill( vpart[it] );
265 else std::cerr << "hmmm histo " << varName[it] << " not found" << std::endl;
266
267 }
268
269
270 if ( matchedreco ) {
271
272 // efficiency histos
273 m_eff_pt->Fill(std::fabs(pTt));
274 m_eff_z0->Fill(z0t);
275 m_eff_eta->Fill(etat);
276 m_eff_phi->Fill(phit);
277 m_eff_d0->Fill(d0t);
278 m_eff_a0->Fill(a0t);
279
280 // signed pT
281 if ( pTt<0 ) m_eff_ptm->Fill(std::fabs(pTt));
282 else m_eff_ptp->Fill(std::fabs(pTt));
283
284 // residual histos
285 double pTr = matchedreco->pT();
286 double z0r = matchedreco->z0();
287 double etar = matchedreco->eta();
288 double phir = matchedreco->phi();
289 double d0r = matchedreco->a0();
290 double a0r = matchedreco->a0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
291
292 if ( m_h2m ) m_h2m->Fill( phit, d0t );
293
294 // if ( m_print ) std::cout << "SigAnalysis::execute() \t " << name() << "\t" << i << " "
295 if ( m_print ) std::cout << "SigAnalysis::execute() \t\t" << i << " "
296 << *reftracks[i] << " -> " << *matchedreco << "\t"
297 << pTr << " " << pTt << " " << d0r << " " << d0t << std::endl;
298
299 double vres[6] = { 1/pTt-1/pTr, etat-etar, phit-phir, z0t-z0r, d0t-d0r, a0t-a0r };
300 for ( int it=0 ; it<6 ; it++ ) {
301 if ( TH1F* hptr = find(varName[it]+"_res") ) hptr->Fill( vres[it] );
302 else std::cerr << "hmmm histo " << varName[it]+"_res" << " not found" << std::endl;
303 }
304
305 // in this loop over the reference tracks, could fill efficiency
306 // histograms
307
308 }
309 else {
310 // fill efficiencies with unmatched histos
311 // std::cout << "NULL" << std::endl;
312 m_eff_pt->FillDenom(std::fabs(pTt));
313 m_eff_z0->FillDenom(z0t);
314 m_eff_eta->FillDenom(etat);
315 m_eff_phi->FillDenom(phit);
316 m_eff_d0->FillDenom(d0t);
317 m_eff_a0->FillDenom(a0t);
318
319 // signed pT
320 if ( pTt<0 ) m_eff_ptm->FillDenom(std::fabs(pTt));
321 else m_eff_ptp->FillDenom(std::fabs(pTt));
322
323 if ( std::fabs(pTt)>4000 ) dump = true;
324 }
325
326 }
327
328
329 // for fake/purity histograms, loop over the test tracks
330 // and get the corresponding matched reference tracks from the
331 // reverse map in the TrackAscociator class - revmatched()
332
333 if ( m_print ) std::cout << "SigAnalysis::execute() \t " << name() << "\t " << m_icount << " events\t " << testtracks.size() << " tracks (" << m_Nreco << ")" << "\n---------------" << std::endl;
334
335 m_icount++;
336
337 m_Nreco += testtracks.size();
338
339 for ( int i=testtracks.size() ; i-- ; ) {
340
341 // std::cout << "\t\tSigAnalysis purity " << name() << "\t" << i << " " << *testtracks[i] << " -> ";
342
343 // double pTr = std::fabs(testtracks[i]->pT());
344 double pTr = testtracks[i]->pT();
345 double etar = testtracks[i]->eta();
346 double phir = testtracks[i]->phi();
347 double z0r = testtracks[i]->z0();
348 double d0r = testtracks[i]->a0();
349 double a0r = testtracks[i]->a0() + sin(phir)*m_xBeamTest - cos(phir)*m_yBeamTest; // this will be changed when we know the beam spot position
350 // double a0rp = testtracks[i]->a0() - sin(phir)*m_xBeam - cos(phir)*m_yBeam; // this will be changed when we know the beam spot position
351
352 // std::cout << "d0 " << d0r << "\tphi " << phir << "\tx " << m_xBeamTest << "\ty " << m_yBeamTest << std::endl;
353
354 double nsctr = testtracks[i]->sctHits();
355 double npixr = testtracks[i]->pixelHits();
356
357 double ntrtr = testtracks[i]->trHits();
358 double nstrawr = testtracks[i]->strawHits();
359
360 if ( m_h2r ) m_h2r->Fill( phir, d0r );
361
362 const TIDA::Track* matchedref = matcher->revmatched(testtracks[i]);
363
364 // if ( matchedref ) std::cout << *matchedref << std::endl;
365 // else std::cout << "NULL" << std::endl;
366
367 // raw test track distributions
368 double vpart[10] = { std::fabs(pTr), etar, phir, z0r, d0r, a0r, nsctr, npixr, ntrtr, nstrawr };
369 for ( int it=0 ; it<10 ; it++ ) {
370 // std::string hname = name()+"_"+varName[it]+"_rec";
371 // std::string hname = varName[it]+"_rec";
372 // std::map<std::string, TH1F*>::iterator hmitr = m_histos.find(hname);
373 // if ( hmitr!=m_histos.end() ) hmitr->second->Fill( vpar[it] );
374 // else std::cerr << "hmmm histo " << hname << " not found" << std::endl;
375 if ( TH1F* hptr = find(varName[it]+"_rec") ) hptr->Fill( vpart[it] );
376 else std::cerr << "hmmm histo " << varName[it]+"_rec" << " not found" << std::endl;
377 }
378
379 // purities
380 if ( matchedref ) {
381
382 // std::cout << *matchedref << std::endl;
383
384 m_purity_pt->Fill(std::fabs(pTr));
385 m_purity_z0->Fill(z0r);
386 m_purity_eta->Fill(etar);
387 m_purity_phi->Fill(phir);
388 m_purity_d0->Fill(d0r);
389 m_purity_a0->Fill(a0r);
390
391 }
392 else {
393 // std::cout << "NULL" << std::endl;
394 m_purity_pt->FillDenom(std::fabs(pTr));
395 m_purity_z0->FillDenom(z0r);
396 m_purity_eta->FillDenom(etar);
397 m_purity_phi->FillDenom(phir);
398 m_purity_d0->FillDenom(d0r);
399 m_purity_a0->FillDenom(a0r);
400 }
401
402 }
403
404 if ( dump && m_print ) {
405
406 std::cout << "SigAnalysis::execute() missed a high pT track - dumping tracks" << std::endl;
407
408 for ( int i=reftracks.size() ; i-- ; ) {
409
410 if ( std::fabs( reftracks[i]->pT() ) > 1000 ) {
411 std::cout << "\t dump " << *reftracks[i];
412 const TIDA::Track* matchedreco = matcher->matched(reftracks[i]);
413 if ( matchedreco ) std::cout << " <--> " << *matchedreco << std::endl;
414 else std::cout << std::endl;
415 }
416
417 }
418
419 for ( int i=testtracks.size() ; i-- ; ) {
420 const TIDA::Track* matchedref = matcher->revmatched(testtracks[i]);
421 if ( matchedref==0 ) std::cout << "\t\t\t\t\t " << *testtracks[i] << std::endl;
422 }
423
424 }
425
426 // std::cout << "SigAnalysis::execute() exiting" << std::endl;
427
428}
429
430
431
TrackAssociator
TIDA::Associator< TIDA::Track > TrackAssociator
Definition TrackAssociator.h:26
SigAnalysis::initialise
virtual void initialise()
standard operation interface
Definition SigAnalysis.cxx:13
SigAnalysis::m_dir
TIDDirectory * m_dir
Definition SigAnalysis.h:116
SigAnalysis::addHistogram
void addHistogram(TH1F *h)
Definition SigAnalysis.h:70
SigAnalysis::m_purity_d0
Efficiency1D * m_purity_d0
Definition SigAnalysis.h:99
SigAnalysis::finalise
virtual void finalise()
Definition SigAnalysis.cxx:177
SigAnalysis::m_purity_eta
Efficiency1D * m_purity_eta
Definition SigAnalysis.h:96
SigAnalysis::m_print
bool m_print
flag to print out the matched tracks etc
Definition SigAnalysis.h:111
SigAnalysis::m_eff_ptp
Efficiency1D * m_eff_ptp
Definition SigAnalysis.h:86
SigAnalysis::m_eff_a0
Efficiency1D * m_eff_a0
Definition SigAnalysis.h:93
SigAnalysis::m_purity_a0
Efficiency1D * m_purity_a0
Definition SigAnalysis.h:100
SigAnalysis::m_eff_d0
Efficiency1D * m_eff_d0
Definition SigAnalysis.h:92
SigAnalysis::m_purity_z0
Efficiency1D * m_purity_z0
Definition SigAnalysis.h:98
SigAnalysis::m_eff_z0
Efficiency1D * m_eff_z0
Definition SigAnalysis.h:91
SigAnalysis::m_histos
std::map< std::string, TH1F * > m_histos
Definition SigAnalysis.h:83
SigAnalysis::execute
virtual void execute(const std::vector< TIDA::Track * > &reftracks, const std::vector< TIDA::Track * > &testtracks, TrackAssociator *matcher)
Definition SigAnalysis.cxx:205
SigAnalysis::m_eff_phi
Efficiency1D * m_eff_phi
Definition SigAnalysis.h:90
SigAnalysis::m_icount
int m_icount
Event counter.
Definition SigAnalysis.h:114
SigAnalysis::m_eff_ptm
Efficiency1D * m_eff_ptm
Definition SigAnalysis.h:87
SigAnalysis::m_eff_eta
Efficiency1D * m_eff_eta
Definition SigAnalysis.h:89
SigAnalysis::m_purity_pt
Efficiency1D * m_purity_pt
Definition SigAnalysis.h:95
SigAnalysis::m_eff_pt
Efficiency1D * m_eff_pt
Definition SigAnalysis.h:85
SigAnalysis::find
TH1F * find(const std::string &n)
Definition SigAnalysis.h:75
SigAnalysis::m_Nreco
int m_Nreco
number of reconstructed tracks
Definition SigAnalysis.h:108
SigAnalysis::m_purity_phi
Efficiency1D * m_purity_phi
Definition SigAnalysis.h:97
TIDA::Associator::revmatched
virtual const T * revmatched(S *t)
Definition TIDAAssociator.h:52
TIDA::Associator::matched
virtual const S * matched(T *t)
Definition TIDAAssociator.h:45
T_Efficiency::finalise
void finalise(double scale=100)
actually calculate the efficiencies
Definition T_Efficiency.h:71
TrackAnalysis::m_yBeamReference
double m_yBeamReference
Definition TrackAnalysis.h:154
TrackAnalysis::m_xBeamReference
double m_xBeamReference
beamline positions reference sample
Definition TrackAnalysis.h:153
TrackAnalysis::name
const std::string & name() const
return identifier
Definition TrackAnalysis.h:52
TrackAnalysis::m_xBeamTest
double m_xBeamTest
test sample
Definition TrackAnalysis.h:158
dump
-event-from-file
Generated on for ATLAS Offline Software by doxygen 1.14.0