BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmL1Performance.cxx
Go to the documentation of this file.
1/*
2 *====================================================================
3 *
4 * CBM Level 1 Reconstruction
5 *
6 * Authors: I.Kisel, S.Gorbunov
7 *
8 * e-mail : ikisel@kip.uni-heidelberg.de
9 *
10 *====================================================================
11 *
12 * L1 Fit performance
13 *
14 *====================================================================
15 */
16#include "CbmKF.h"
17#include "CbmKFMath.h"
18#include "CbmKFTrack.h" // for vertex pulls
19#include "CbmL1.h"
20#include "CbmL1Constants.h"
21#include "CbmL1Counters.h"
22#include "CbmL1PartEfficiencies.h"
23#include "CbmL1ParticlesFinder.h"
24#include "CbmMCTrack.h"
25#include "CbmMvdHit.h"
26#include "CbmMvdHitMatch.h"
27#include "CbmMvdPoint.h"
28#include "CbmStsDigi.h"
29#include "CbmStsHit.h"
30#include "CbmStsPoint.h"
31#include "CbmStsSector.h" // for field FieldCheck.
32#include "CbmStsSensor.h" // for field FieldCheck.
33#include "CbmStsStation.h" // for field FieldCheck.
34#include "FairTrackParam.h" // for vertex pulls
35#include "L1Algo/L1AddMaterial.h" // for vertex pulls
36#include "L1Algo/L1Algo.h"
37#include "L1Algo/L1Extrapolation.h" // for vertex pulls
38#include "L1Algo/L1StsHit.h"
39#include "TFile.h"
40#include "TH1.h"
41#include "TH2.h"
42#include "TProfile.h"
43
44#include <iostream>
45#include <list>
46#include <map>
47#include <vector>
48
49using std::cout;
50using std::endl;
51using std::ios;
52using std::map;
53using std::pair;
54using std::vector;
55
56void CbmL1::TrackMatch()
57{
58 map<int, CbmL1MCTrack*> pMCTrackMap;
59 pMCTrackMap.clear();
60
61 // fill pMCTrackMap
62 for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) {
63 CbmL1MCTrack& MC = *i;
64
65 if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()) {
66 pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC));
67 } else {
68 cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
69 }
70 }
71 // -- prepare information about reconstructed tracks
72 const int nRTracks = vRTracks.size();
73 for (int iR = 0; iR < nRTracks; iR++) {
74 CbmL1Track* prtra = &(vRTracks[iR]);
75
76 int hitsum = prtra->StsHits.size(); // number of hits in track
77
78 // count how many hits from each mcTrack belong to current recoTrack
79 map<int, int>& hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack
80 for (vector<int>::iterator ih = (prtra->StsHits).begin(); ih != (prtra->StsHits).end(); ++ih) {
81
82 const int nMCPoints = vStsHits[*ih].mcPointIds.size();
83 for (int iP = 0; iP < nMCPoints; iP++) {
84 int iMC = vStsHits[*ih].mcPointIds[iP];
85 int ID = -1;
86 if (iMC >= 0)
87 ID = vMCPoints[iMC].ID;
88 if (hitmap.find(ID) == hitmap.end())
89 hitmap[ID] = 1;
90 else {
91 hitmap[ID] += 1;
92 }
93 } // for iPoint
94 } // for iHit
95
96 // RTrack <-> MCTrack identification
97 double max_percent = 0.0; // [%]. maximum persent of hits, which belong to one mcTrack.
98 for (map<int, int>::iterator posIt = hitmap.begin(); posIt != hitmap.end(); posIt++)
99 { // loop over all touched MCTracks
100
101 if (posIt->first < 0)
102 continue; // not a MC track - based on fake hits
103
104 // count max-purity
105 if (double(posIt->second) > max_percent * double(hitsum))
106 max_percent = double(posIt->second) / double(hitsum);
107
108 // set relation to the mcTrack
109 if (double(posIt->second) > CbmL1Constants::MinPurity * double(hitsum)) { // found correspondent MCTrack
110 if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end())
111 continue;
112 CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
113
114 pmtra->AddRecoTrack(prtra);
115 prtra->AddMCTrack(pmtra);
116 } else {
117 if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end())
118 continue;
119 CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
120
121 pmtra->AddTouchTrack(prtra);
122 }
123 } // for hitmap
124 prtra->SetMaxPurity(max_percent);
125 } // for reco tracks
126} //
127
128struct TL1PerfEfficiencies : public TL1Efficiencies
129{
130 TL1PerfEfficiencies()
131 : TL1Efficiencies()
132 , ratio_killed()
133 , ratio_clone()
134 , ratio_length()
135 , ratio_fakes()
136 , killed()
137 , clone()
138 , reco_length()
139 , reco_fakes()
140 , mc_length()
141 , mc_length_hits()
142 {
143 // add total efficiency
144 AddCounter("long_fast_prim", "LongRPrim efficiency");
145 AddCounter("fast_prim", "RefPrim efficiency");
146 AddCounter("fast_sec", "RefSec efficiency");
147 AddCounter("fast", "Refset efficiency");
148 AddCounter("total", "Allset efficiency");
149 AddCounter("slow_prim", "ExtraPrim efficiency");
150 AddCounter("slow_sec", "ExtraSec efficiency");
151 AddCounter("slow", "Extra efficiency");
152 AddCounter("d0", "D0 efficiency");
153 AddCounter("short", "Short123s efficiency");
154 AddCounter("shortPion", "Short123s pion eff");
155 AddCounter("shortProton", "Short123s proton eff");
156 AddCounter("shortKaon", "Short123s kaon eff");
157 AddCounter("shortE", "Short123s e eff");
158 AddCounter("shortRest", "Short123s rest eff");
159
160 AddCounter("fast_sec_e", "RefSec E efficiency");
161 AddCounter("fast_e", "Refset E efficiency");
162 AddCounter("total_e", "Allset E efficiency");
163 AddCounter("slow_sec_e", "ExtraSecE efficiency");
164 AddCounter("slow_e", "Extra E efficiency");
165 }
166
167 virtual ~TL1PerfEfficiencies() {};
168
169 virtual void AddCounter(TString shortname, TString name)
170 {
171 TL1Efficiencies::AddCounter(shortname, name);
172 ratio_killed.AddCounter();
173 ratio_clone.AddCounter();
174 ratio_length.AddCounter();
175 ratio_fakes.AddCounter();
176 killed.AddCounter();
177 clone.AddCounter();
178 reco_length.AddCounter();
179 reco_fakes.AddCounter();
180 mc_length.AddCounter();
181 mc_length_hits.AddCounter();
182 };
183
184 TL1PerfEfficiencies& operator+=(TL1PerfEfficiencies& a)
185 {
186 TL1Efficiencies::operator+=(a);
187 killed += a.killed;
188 clone += a.clone;
189 reco_length += a.reco_length;
190 reco_fakes += a.reco_fakes;
191 mc_length += a.mc_length;
192 mc_length_hits += a.mc_length_hits;
193 return *this;
194 };
195
196 void CalcEff()
197 {
198 TL1Efficiencies::CalcEff();
199 ratio_killed = killed / mc;
200 ratio_clone = clone / mc;
201 TL1TracksCatCounters<int> allReco = reco + clone;
202 ratio_length = reco_length / allReco;
203 ratio_fakes = reco_fakes / allReco;
204 };
205
206 void Inc(bool isReco,
207 bool isKilled,
208 double _ratio_length,
209 double _ratio_fakes,
210 int _nclones,
211 int _mc_length,
212 int _mc_length_hits,
213 TString name)
214 {
215 TL1Efficiencies::Inc(isReco, name);
216
217 const int index = indices[name];
218
219 if (isKilled)
220 killed.counters[index]++;
221 reco_length.counters[index] += _ratio_length;
222 reco_fakes.counters[index] += _ratio_fakes;
223 clone.counters[index] += _nclones;
224 mc_length.counters[index] += _mc_length;
225 mc_length_hits.counters[index] += _mc_length_hits;
226 };
227
228 void PrintEff()
229 {
230 L1_assert(nEvents != 0);
231
232 cout.setf(ios::fixed);
233 cout.setf(ios::showpoint);
234 cout.precision(3);
235 cout.setf(ios::right);
236 cout << "Track category : " << " Eff " << " / " << "Killed" << " / " << "Length" << " / " << "Fakes "
237 << " / " << "Clones" << " / " << "All Reco" << " | " << " All MC " << " / " << "MCl(hits)" << " / "
238 << "MCl(MCps)" << endl;
239
240 int NCounters = mc.NCounters;
241 for (int iC = 0; iC < NCounters; iC++) {
242 if ((names[iC] != "D0 efficiency") || (mc.counters[iC] != 0))
243 cout << names[iC] << " : " << ratio_reco.counters[iC] << " / "
244 << ratio_killed.counters[iC] // tracks with aren't reco because other tracks takes their hit(-s)
245 << " / " << ratio_length.counters[iC] // nRecoMCHits/nMCHits
246 << " / " << ratio_fakes.counters[iC] // nFakeHits/nRecoAllHits
247 << " / " << ratio_clone.counters[iC] // nCloneTracks/nMCTracks
248 << " / " << setw(8) << reco.counters[iC] / double(nEvents) << " | " << setw(8)
249 << mc.counters[iC] / double(nEvents) << " / "
250 << mc_length_hits.counters[iC] / double(mc.counters[iC]) << " / "
251 << mc_length.counters[iC] / double(mc.counters[iC]) << endl;
252 }
253 cout << "Ghost probability : " << ratio_ghosts << " | " << ghosts << endl;
254 };
255
256 TL1TracksCatCounters<double> ratio_killed;
257 TL1TracksCatCounters<double> ratio_clone;
258 TL1TracksCatCounters<double> ratio_length;
259 TL1TracksCatCounters<double> ratio_fakes;
260
261 TL1TracksCatCounters<int> killed;
262 TL1TracksCatCounters<int> clone;
263 TL1TracksCatCounters<double> reco_length;
264 TL1TracksCatCounters<double> reco_fakes;
265 TL1TracksCatCounters<int> mc_length;
266 TL1TracksCatCounters<int> mc_length_hits;
267};
268
269void CbmL1::EfficienciesPerformance()
270{
271 static TL1PerfEfficiencies L1_NTRA; // average efficiencies
272
273 static int L1_NEVENTS = 0;
274 static double L1_CATIME = 0.0;
275
276 TL1PerfEfficiencies ntra; // efficiencies for current event
277
278 for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt) {
279 ntra.ghosts += rtraIt->IsGhost();
280 // if(rtraIt->IsGhost()){ // Debug.
281 // cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " |
282 // "; for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++
283 // ){
284 // cout << " (" << posIt->second << " " << posIt->first << ") ";
285 // }
286 // cout << endl;
287 // }
288 }
289
290 int sta_nhits[algo->NStations], sta_nfakes[algo->NStations];
291
292 for (int i = 0; i < algo->NStations; i++) {
293 sta_nhits[i] = 0;
294 sta_nfakes[i] = 0;
295 }
296
297 for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++) {
298 CbmL1MCTrack& mtra = *(mtraIt);
299
300 // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
301 if (!mtra.IsReconstructable() && !mtra.IsAdditional())
302 continue;
303
304 // -- find used constans --
305 // is track reconstructed
306 const bool reco = (mtra.IsReconstructed());
307 // is track killed. At least one hit of it belong to some recoTrack
308 const bool killed = !reco && mtra.IsDisturbed();
309 // ration length for current mcTrack
310 vector<CbmL1Track*>& rTracks = mtra.GetRecoTracks(); // for length calculations
311 double ratio_length = 0;
312 double ratio_fakes = 0;
313 double mc_length_hits = mtra.NStations();
314 int mc_length = mtra.NMCStations();
315 if (reco) {
316 for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
317 ratio_length +=
318 static_cast<double>(rTracks[irt]->GetNOfHits()) * rTracks[irt]->GetMaxPurity() / mc_length_hits;
319 ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
320 }
321 }
322 // number of clones
323 int nclones = 0;
324 if (reco)
325 nclones = mtra.GetNClones();
326 // if (nclones){ // Debug. Look at clones
327 // for (int irt = 0; irt < rTracks.size(); irt++){
328 // const int ista = vHitStore[rTracks[irt]->StsHits[0]].iStation;
329 // cout << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
330 // }
331 // cout << mtra.NStations() << endl;
332 // }
333
334 if (mtra.IsAdditional()) { // short
335 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
336 switch (mtra.pdg) {
337 case 11:
338 case -11:
339 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
340 break;
341 case 211:
342 case -211:
343 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortPion");
344 break;
345 case 321:
346 case -321:
347 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortKaon");
348 break;
349 case 2212:
350 case -2212:
351 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits,
352 "shortProton");
353 break;
354 default:
355 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortRest");
356 }
357 } else { // separate all efficiecies from short eff
358
359 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
360
361 if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
362 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
363 }
364
365 if (mtra.p > CbmL1Constants::MinRefMom) { // reference tracks
366 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast");
367
368 if (mtra.IsPrimary()) { // reference primary
369 if (mtra.NStations() == NStation) { // long reference primary
370 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits,
371 "long_fast_prim");
372 }
373 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_prim");
374 } else { // reference secondary
375 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec");
376 }
377 } else { // extra set of tracks
378 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow");
379
380 if (mtra.IsPrimary()) { // extra primary
381 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_prim");
382 } else {
383 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec");
384 }
385 } // if extra
386
387 if (mtra.pdg == 11 || mtra.pdg == -11) {
388 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");
389
390 if (mtra.p > CbmL1Constants::MinRefMom) { // reference tracks
391 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
392
393 if (mtra.IsPrimary()) { // reference primary
394 } else { // reference secondary
395 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits,
396 "fast_sec_e");
397 }
398 } else { // extra set of tracks
399 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
400
401 if (mtra.IsPrimary()) { // extra primary
402 } else {
403 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits,
404 "slow_sec_e");
405 }
406 } // if extra
407 }
408 }
409
410 } // for mcTracks
411
412 L1_CATIME += algo->CATime;
413 L1_NEVENTS++;
414 ntra.IncNEvents();
415 L1_NTRA += ntra;
416
417 ntra.CalcEff();
418 L1_NTRA.CalcEff();
419
420 // cout.precision(3);
421 if (fVerbose) {
422 if (fVerbose > 1) {
423 ntra.PrintEff();
424 cout << "Number of true and fake hits in stations: " << endl;
425 for (int i = 0; i < algo->NStations; i++) {
426 cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
427 }
428 cout << endl;
429 } // fVerbose > 1
430 cout << endl;
431 cout << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS " << endl << endl;
432 L1_NTRA.PrintEff();
433 cout << "MC tracks/event found : "
434 << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS)) << endl;
435 cout << endl;
436 cout << "CA Track Finder: " << L1_CATIME / L1_NEVENTS << " s/ev" << endl << endl;
437 }
438} // void CbmL1::Performance()
439
440void CbmL1::GetMCParticles()
441{
442 vMCParticles.clear();
443 // convert
444 for (int i = 0; i < listMCTracks->GetEntriesFast(); i++) {
445 CbmMCTrack& mtra = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(i)));
446 KFMCParticle part;
447 part.SetMCTrackID(i);
448 part.SetMotherId(mtra.GetMotherId());
449 part.SetPDG(mtra.GetPdgCode());
450 vMCParticles.push_back(part);
451 }
452
453 // find relations
454 const unsigned int nParticles = vMCParticles.size();
455 for (unsigned int iP = 0; iP < nParticles; iP++) {
456 KFMCParticle& part = vMCParticles[iP];
457 for (unsigned int iP2 = 0; iP2 < nParticles; iP2++) {
458 KFMCParticle& part2 = vMCParticles[iP2];
459
460 if (part.GetMotherId() == part2.GetMCTrackID()) {
461 part2.AddDaughter(iP);
462 }
463 }
464 }
465 // // correct idices according to the new array
466 // for ( unsigned int iP = 0; iP < nParticles; iP++ ) {
467 // KFMCParticle &part = vMCParticles[iP];
468 // for(unsigned int iP2 = 0; iP2 < nParticles; iP2++) {
469 // KFMCParticle &part2 = vMCParticles[iP2];
470
471 // if(part.GetMotherId() == part2.GetMCTrackID()) {
472 // part.SetMotherId (iP2);
473 // }
474 // }
475 // }
476}
477
478void CbmL1::FindReconstructableMCParticles()
479{
480 const unsigned int nParticles = vMCParticles.size();
481
482 for (unsigned int iP = 0; iP < nParticles; iP++) {
483 KFMCParticle& part = vMCParticles[iP];
484 CheckMCParticleIsReconstructable(part);
485 }
486}
487
488void CbmL1::CheckMCParticleIsReconstructable(KFMCParticle& part)
489{
490
491 if (part.IsReconstructable())
492 return;
493
494 const int nParticles = 5;
495 int partPDG[nParticles] = {310, 3122, 3312, -3312, 3334};
496 vector<int> partDaughterPdg[nParticles];
497
498 partDaughterPdg[0].push_back(211);
499 partDaughterPdg[0].push_back(-211);
500
501 partDaughterPdg[1].push_back(2212);
502 partDaughterPdg[1].push_back(-211);
503
504 partDaughterPdg[2].push_back(3122);
505 partDaughterPdg[2].push_back(-211);
506
507 partDaughterPdg[3].push_back(3122);
508 partDaughterPdg[3].push_back(211);
509
510 partDaughterPdg[4].push_back(3122);
511 partDaughterPdg[4].push_back(-321);
512
513 // tracks
514 if (/*part.NDaughters() == 0*/ part.GetPDG() == -211 || part.GetPDG() == 211 || part.GetPDG() == 2212
515 || part.GetPDG() == 321 || part.GetPDG() == -321)
516 { // TODO other particles
517
518 switch (fFindParticlesMode) {
519 case 1:
520 part.SetAsReconstructable();
521 break;
522 case 2:
523 unsigned int i;
524 for (i = 0; i < vMCTracks.size(); i++) // TODO opt
525 if (vMCTracks[i].ID == part.GetMCTrackID())
526 break;
527 if (i != vMCTracks.size())
528 if (vMCTracks[i].IsReconstructable())
529 part.SetAsReconstructable();
530 break;
531 case 3:
532 // int i;
533 for (i = 0; i < vMCTracks.size(); i++)
534 if (vMCTracks[i].ID == part.GetMCTrackID())
535 break;
536 if (i != vMCTracks.size())
537 if (vMCTracks[i].IsReconstructed())
538 part.SetAsReconstructable();
539 break;
540 default:
541 L1_assert(0);
542 };
543 }
544
545 // mother particles
546 else
547 {
548 for (int iPart = 0; iPart < nParticles; iPart++) {
549 if (part.GetPDG() == partPDG[iPart]) {
550 const unsigned int nDaughters = partDaughterPdg[iPart].size();
551 if (part.GetDaughterIds().size() != nDaughters)
552 return;
553 int pdg[nDaughters];
554
555 for (unsigned int iD = 0; iD < nDaughters; iD++)
556 pdg[iD] = vMCParticles[part.GetDaughterIds()[iD]].GetPDG();
557
558 bool isDaughterFound[nDaughters];
559 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
560 isDaughterFound[iDMC] = 0;
561
562 bool isReco = 1;
563 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
564 for (unsigned int iD = 0; iD < nDaughters; iD++)
565 if (pdg[iD] == partDaughterPdg[iPart][iDMC])
566 isDaughterFound[iDMC] = 1;
567
568 for (unsigned int iDMC = 0; iDMC < nDaughters; iDMC++)
569 isReco = isReco && isDaughterFound[iDMC];
570
571 if (!isReco)
572 return;
573 }
574 }
575
576 const vector<int>& dIds = part.GetDaughterIds();
577 const unsigned int nD = dIds.size();
578 bool reco = 1;
579 for (unsigned int iD = 0; iD < nD && reco; iD++) {
580 KFMCParticle& dp = vMCParticles[dIds[iD]];
581 CheckMCParticleIsReconstructable(dp);
582 reco &= dp.IsReconstructable();
583 }
584 if (reco)
585 part.SetAsReconstructable();
586 }
587}
588
589/// Procedure for match Reconstructed and MC Particles. Should be called before Performances
590void CbmL1::MatchParticles()
591{
592 // get all reco particles ( temp )
593 MCtoRParticleId.clear();
594 RtoMCParticleId.clear();
595 MCtoRParticleId.resize(vMCParticles.size());
596 RtoMCParticleId.resize(vRParticles.size());
597
598 // match tracks ( particles which are direct copy of tracks )
599 for (unsigned int iRP = 0; iRP < vRParticles.size(); iRP++) {
600 CbmKFParticle& rPart = vRParticles[iRP];
601 L1_assert(rPart.NDaughters() > 0);
602 L1_ASSERT(static_cast<unsigned int>(rPart.Id()) == iRP, rPart.Id() << " != " << iRP);
603 if (rPart.NDaughters() != 1)
604 continue;
605
606 const int rTrackId = rPart.DaughterIds()[0];
607
608 if (vRTracks[rTrackId].IsGhost())
609 continue;
610 const int mcTrackId = vRTracks[rTrackId].GetMCTrack()->ID;
611
612 for (unsigned int iMP = 0; iMP < vMCParticles.size(); iMP++) {
613 KFMCParticle& mPart = vMCParticles[iMP];
614 if (mPart.GetMCTrackID() == mcTrackId) { // match is found
615 if (mPart.GetPDG() == rPart.GetPDG()) {
616 MCtoRParticleId[iMP].ids.push_back(iRP);
617 RtoMCParticleId[iRP].ids.push_back(iMP);
618 } else {
619 MCtoRParticleId[iMP].idsMI.push_back(iRP);
620 RtoMCParticleId[iRP].idsMI.push_back(iMP);
621 }
622 }
623 }
624 }
625
626 // match created mother particles
627 for (unsigned int iRP = 0; iRP < vRParticles.size(); iRP++) {
628 CbmKFParticle& rPart = vRParticles[iRP];
629 const unsigned int NRDaughters = rPart.NDaughters();
630 if (NRDaughters < 2)
631 continue;
632
633 unsigned int iD = 0;
634 int mmId = -2; // mother MC id
635 {
636 const int rdId = rPart.DaughterIds()[iD];
637 if (!RtoMCParticleId[rdId].IsMatched())
638 continue;
639 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
640 mmId = vMCParticles[mdId].GetMotherId();
641 }
642 iD++;
643 for (; iD < NRDaughters; iD++) {
644 const int rdId = rPart.DaughterIds()[iD];
645 if (!RtoMCParticleId[rdId].IsMatched())
646 break;
647 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
648 if (vMCParticles[mdId].GetMotherId() != mmId)
649 break;
650 }
651 if (iD == NRDaughters && mmId != -1) { // match is found and it is not primary vertex
652 KFMCParticle& mmPart = vMCParticles[mmId];
653 if (mmPart.GetPDG() == rPart.GetPDG() && mmPart.NDaughters() == rPart.NDaughters()) {
654 MCtoRParticleId[mmId].ids.push_back(iRP);
655 RtoMCParticleId[iRP].ids.push_back(mmId);
656 } else {
657 MCtoRParticleId[mmId].idsMI.push_back(iRP);
658 RtoMCParticleId[iRP].idsMI.push_back(mmId);
659 }
660 }
661 }
662}
663
664void CbmL1::PartEffPerformance()
665{
666 static int NEVENTS = 0;
667
668 static CbmL1PartEfficiencies PARTEFF; // average efficiencies
669
670 CbmKFPartEfficiencies partEff; // efficiencies for current event
671
672 const int NRP = vRParticles.size();
673 for (int iP = 0; iP < NRP; ++iP) {
674 const CbmKFParticle& part = vRParticles[iP];
675 const int pdg = part.GetPDG();
676
677 const bool isBG = RtoMCParticleId[iP].idsMI.size() != 0;
678 const bool isGhost = !RtoMCParticleId[iP].IsMatched();
679
680 for (int iPart = 0; iPart < PARTEFF.nParticles; iPart++)
681 if (pdg == PARTEFF.partPDG[iPart])
682 partEff.IncReco(isGhost, isBG, PARTEFF.partName[iPart].Data());
683 }
684
685 const int NMP = vMCParticles.size();
686 for (int iP = 0; iP < NMP; ++iP) {
687 const KFMCParticle& part = vMCParticles[iP];
688 if (!part.IsReconstructable())
689 continue;
690 const int pdg = part.GetPDG();
691 const int mId = part.GetMotherId();
692
693 const bool isReco = MCtoRParticleId[iP].ids.size() != 0;
694 const int nClones = MCtoRParticleId[iP].ids.size() - 1;
695
696 for (int iPart = 0; iPart < PARTEFF.nParticles; iPart++) {
697 if (pdg == PARTEFF.partPDG[iPart]) {
698 partEff.Inc(isReco, nClones, PARTEFF.partName[iPart].Data());
699 if (mId == -1)
700 partEff.Inc(isReco, nClones, (PARTEFF.partName[iPart] + "_prim").Data());
701 else
702 partEff.Inc(isReco, nClones, (PARTEFF.partName[iPart] + "_sec").Data());
703 }
704 }
705 }
706
707 NEVENTS++;
708
709 PARTEFF += partEff;
710
711 partEff.CalcEff();
712 PARTEFF.CalcEff();
713
714 // cout.precision(3);
715 if (fVerbose) {
716 if (NEVENTS % 500 == 0) {
717 cout << " ---- KF Particle finder --- " << endl;
718 // cout << "L1 STAT : " << fNEvents << " EVENT " << endl << endl;
719 // partEff.PrintEff();
720 // cout << endl;
721 cout << "ACCUMULATED STAT : " << NEVENTS << " EVENTS " << endl << endl;
722 PARTEFF.PrintEff();
723
724 cout << endl;
725 // cout<<"CA Track Finder: " << L1_CATIME/L1_fNEvents << " s/ev" << endl << endl;
726 }
727 }
728}
729
730void CbmL1::PartHistoPerformance()
731{
732 static const int NParticles = 2; // Ks, Lambda
733
734 static const int nFitQA = 16;
735 static TH1F* hFitDaughtersQA[NParticles][nFitQA];
736 static TH1F* hFitQA[NParticles][nFitQA];
737
738 static const int nHistoPartParam = 5;
739 static TH1F* hPartParam[NParticles][nHistoPartParam]; // mass, decay length, c*tau, chi/ndf, prob ...
740 static TH1F* hPartParamBG[NParticles][nHistoPartParam];
741 static TH1F* hPartParamSignal[NParticles][nHistoPartParam];
742 static const int nHistoPartParamQA = 3;
743 // static TH1F *hPartParamQA[NParticles][nHistoPartParamQA*2]; // residuals and pulls of these parameters
744
745 static const int nHistosPV = 6;
746 static TH1F* hPVFitQa[nHistosPV];
747
748 TString patName[NParticles] = {"K0s", "Lambda"};
749
750 static bool first_call = 1;
751
752 if (first_call) {
753 first_call = 0;
754
755 TDirectory* currentDir = gDirectory;
756 gDirectory = histodir;
757 gDirectory->cd("L1");
758 gDirectory->mkdir("Particles");
759 gDirectory->cd("Particles");
760 {
761 for (int iPart = 0; iPart < NParticles; ++iPart) {
762 gDirectory->mkdir(patName[iPart].Data());
763 gDirectory->cd(patName[iPart].Data());
764 {
765 TString res = "res";
766 TString pull = "pull";
767
768 gDirectory->mkdir("DaughtersQA");
769 gDirectory->cd("DaughtersQA");
770 {
771 TString parName[nFitQA / 2] = {"X", "Y", "Z", "Px", "Py", "Pz", "E", "M"};
772 int nBins = 100;
773 float xMax[nFitQA / 2] = {0.15, 0.15, 0.03, 0.01, 0.01, 0.06, 0.06, 0.01};
774
775 for (int iH = 0; iH < nFitQA / 2; iH++) {
776 hFitDaughtersQA[iPart][iH] = new TH1F(
777 (res + parName[iH]).Data(), (res + parName[iH]).Data(), nBins, -xMax[iH], xMax[iH]);
778 hFitDaughtersQA[iPart][iH + 8] =
779 new TH1F((pull + parName[iH]).Data(), (pull + parName[iH]).Data(), nBins, -6, 6);
780 }
781 }
782 gDirectory->cd(".."); // particle directory
783
784 gDirectory->mkdir("FitQA");
785 gDirectory->cd("FitQA");
786 {
787 TString parName[nFitQA / 2] = {"X", "Y", "Z", "Px", "Py", "Pz", "E", "M"};
788 int nBins = 50;
789 float xMax[nFitQA / 2] = {0.15, 0.15, 1.2, 0.02, 0.02, 0.15, 0.15, 0.006};
790
791 for (int iH = 0; iH < nFitQA / 2; iH++) {
792 hFitQA[iPart][iH] = new TH1F((res + parName[iH]).Data(), (res + parName[iH]).Data(), nBins,
793 -xMax[iH], xMax[iH]);
794 hFitQA[iPart][iH + 8] =
795 new TH1F((pull + parName[iH]).Data(), (pull + parName[iH]).Data(), nBins, -6, 6);
796 }
797 }
798 gDirectory->cd(".."); // particle directory
799
800 gDirectory->mkdir("Parameters");
801 gDirectory->cd("Parameters");
802 {
803 TString parName[nHistoPartParam] = {"M", "DL", "c#tau", "chi2ndf", "prob"};
804 int nBins[nHistoPartParam] = {1000, 100, 100, 100, 100};
805 float xMin[nHistoPartParam] = {0.3f, 0., 0., 0., 0.};
806 float xMax[nHistoPartParam] = {1.3f, 30., 30., 20., 1.};
807 if (iPart == 1) {
808 xMin[0] = 1.0;
809 xMax[0] = 2.0;
810 }
811
812 for (int iH = 0; iH < nHistoPartParam; iH++)
813 hPartParam[iPart][iH] =
814 new TH1F(parName[iH].Data(), parName[iH].Data(), nBins[iH], xMin[iH], xMax[iH]);
815
816 gDirectory->mkdir("Signal");
817 gDirectory->cd("Signal");
818 {
819 for (int iH = 0; iH < nHistoPartParam; iH++)
820 hPartParamSignal[iPart][iH] =
821 new TH1F((parName[iH]).Data(), (parName[iH]).Data(), nBins[iH], xMin[iH], xMax[iH]);
822
823 gDirectory->mkdir("QA");
824 gDirectory->cd("QA");
825 {
826 // int nBinsQA = 50;
827 // float xMaxQA[nHistoPartParamQA] = {0.01,0.001,0.001};
828 for (int iH = 0; iH < nHistoPartParamQA; iH++) {
829 // hPartParamQA[iPart][iH] =
830 // new TH1F((res+parName[iH]).Data(), (res+parName[iH]).Data(), nBinsQA,
831 // -xMaxQA[iH],xMaxQA[iH]);
832 // hPartParamQA[iPart][iH+nHistoPartParamQA] =
833 // new TH1F((pull+parName[iH]).Data(), (pull+parName[iH]).Data(), nBinsQA, -6,6);
834 }
835 }
836 gDirectory->cd(".."); // particle directory / Parameters / Signal
837 }
838 gDirectory->cd(".."); // particle directory / Parameters
839 gDirectory->mkdir("Background");
840 gDirectory->cd("Background");
841 {
842 for (int iH = 0; iH < nHistoPartParam; iH++)
843 hPartParamBG[iPart][iH] =
844 new TH1F((parName[iH]).Data(), (parName[iH]).Data(), nBins[iH], xMin[iH], xMax[iH]);
845 }
846 gDirectory->cd(".."); // particle directory
847 }
848 gDirectory->cd(".."); // particle directory
849 }
850 gDirectory->cd(".."); // L1
851 }
852 }
853 gDirectory->cd(".."); // L1
854 gDirectory->mkdir("PrimaryVertexQA");
855 gDirectory->cd("PrimaryVertexQA");
856 {
857 struct
858 {
859 string name;
860 string title;
861 Int_t n;
862 Double_t l, r;
863 } Table[nHistosPV] = {{"PVResX", "x_{rec}-x_{mc}, cm", 100, -0.006f, 0.006f},
864 {"PVResY", "y_{rec}-y_{mc}, cm", 100, -0.006f, 0.006f},
865 {"PVResZ", "z_{rec}-z_{mc}, cm", 100, -0.06f, 0.06f},
866 {"PVPullX", "Pull X", 100, -6.f, 6.f},
867 {"PVPullY", "Pull Y", 100, -6.f, 6.f},
868 {"PVPullZ", "Pull Z", 100, -6.f, 6.f}};
869 for (int iHPV = 0; iHPV < nHistosPV; ++iHPV) {
870 hPVFitQa[iHPV] = new TH1F(Table[iHPV].name.data(), Table[iHPV].title.data(), Table[iHPV].n,
871 Table[iHPV].l, Table[iHPV].r);
872 }
873 }
874 gDirectory = currentDir;
875 }
876
877 for (unsigned int iP = 0; iP < vRParticles.size(); iP++) {
878 int iParticle = -1;
879 switch (vRParticles[iP].GetPDG()) {
880 case 310: {
881 iParticle = 0;
882 } break;
883 case 3122: {
884 iParticle = 1;
885 }
886 }
887 if (iParticle < 0)
888 continue;
889
890 Double_t M, ErrM;
891 Double_t dL, ErrdL; // decay length
892 Double_t cT, ErrcT; // c*tau
893 CbmKFParticle TempPart = vRParticles[iP];
894 TempPart.GetMass(M, ErrM);
895 TempPart.GetDecayLength(dL, ErrdL);
896 TempPart.GetLifeTime(cT, ErrcT);
897 Double_t chi2 = TempPart.GetChi2();
898 Int_t ndf = TempPart.GetNDF();
899 Double_t prob = TMath::Prob(chi2, ndf);
900
901 hPartParam[iParticle][0]->Fill(M);
902 hPartParam[iParticle][1]->Fill(dL);
903 hPartParam[iParticle][2]->Fill(cT);
904 hPartParam[iParticle][3]->Fill(chi2 / ndf);
905 hPartParam[iParticle][4]->Fill(prob);
906
907 if (!RtoMCParticleId[iP].IsMatchedWithPdg()) // background
908 {
909 hPartParamBG[iParticle][0]->Fill(M);
910 hPartParamBG[iParticle][1]->Fill(dL);
911 hPartParamBG[iParticle][2]->Fill(cT);
912 hPartParamBG[iParticle][3]->Fill(chi2 / ndf);
913 hPartParamBG[iParticle][4]->Fill(prob);
914 continue;
915 }
916 hPartParamSignal[iParticle][0]->Fill(M);
917 hPartParamSignal[iParticle][1]->Fill(dL);
918 hPartParamSignal[iParticle][2]->Fill(cT);
919 hPartParamSignal[iParticle][3]->Fill(chi2 / ndf);
920 hPartParamSignal[iParticle][4]->Fill(prob);
921
922 int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
923 KFMCParticle& mcPart = vMCParticles[iMCPart];
924 // Fit quality of the mother particle
925 {
926 int iMCTrack = mcPart.GetMCTrackID();
927 CbmMCTrack& mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(iMCTrack)));
928 int mcDaughterId = mcPart.GetDaughterIds()[0];
929 CbmMCTrack& mcDaughter = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
930
931 Double_t decayVtx[3] = {mcDaughter.GetStartX(), mcDaughter.GetStartY(), mcDaughter.GetStartZ()};
932 Double_t recParam[8] = {0};
933 Double_t errParam[8] = {0};
934
935 for (int iPar = 0; iPar < 3; iPar++) {
936 recParam[iPar] = TempPart.GetParameter(iPar);
937 Double_t error = TempPart.GetCovariance(iPar, iPar);
938 if (error < 0.) {
939 error = 1.e20;
940 }
941 errParam[iPar] = TMath::Sqrt(error);
942 }
943 TempPart.Extrapolate(TempPart.r, TempPart.GetDStoPoint(decayVtx));
944 for (int iPar = 3; iPar < 7; iPar++) {
945 recParam[iPar] = TempPart.GetParameter(iPar);
946 Double_t error = TempPart.GetCovariance(iPar, iPar);
947 if (error < 0.) {
948 error = 1.e20;
949 }
950 errParam[iPar] = TMath::Sqrt(error);
951 }
952
953 Double_t Emc = sqrt(mcTrack.GetP() * mcTrack.GetP() + mcTrack.GetMass() * mcTrack.GetMass());
954 Double_t res[8] = {0}, pull[8] = {0},
955 mcParam[8] = {decayVtx[0], decayVtx[1], decayVtx[2], mcTrack.GetPx(),
956 mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass()};
957 for (int iPar = 0; iPar < 7; iPar++) {
958 res[iPar] = recParam[iPar] - mcParam[iPar];
959 if (fabs(errParam[iPar]) > 1.e-20)
960 pull[iPar] = res[iPar] / errParam[iPar];
961 }
962 res[7] = M - mcParam[7];
963 if (fabs(ErrM) > 1.e-20)
964 pull[7] = res[7] / ErrM;
965
966 for (int iPar = 0; iPar < 8; iPar++) {
967 hFitQA[iParticle][iPar]->Fill(res[iPar]);
968 hFitQA[iParticle][iPar + 8]->Fill(pull[iPar]);
969 }
970
971 // Double_t mcT = mcDaughter.GetStartT() - mcTrack.GetStartT();
972 }
973 // Fit quality of daughters
974 for (int iD = 0; iD < mcPart.NDaughters(); ++iD) {
975 int mcDaughterId = mcPart.GetDaughterIds()[iD];
976 if (!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg())
977 continue;
978 CbmMCTrack& mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
979 int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
980 CbmKFParticle Daughter = vRParticles[recDaughterId];
981 Daughter.GetMass(M, ErrM);
982
983 Double_t decayVtx[3] = {mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ()};
984 Daughter.Extrapolate(Daughter.r, Daughter.GetDStoPoint(decayVtx));
985
986 Double_t Emc = sqrt(mcTrack.GetP() * mcTrack.GetP() + mcTrack.GetMass() * mcTrack.GetMass());
987 Double_t res[8] = {0}, pull[8] = {0}, mcParam[8] = {mcTrack.GetStartX(),
988 mcTrack.GetStartY(),
989 mcTrack.GetStartZ(),
990 mcTrack.GetPx(),
991 mcTrack.GetPy(),
992 mcTrack.GetPz(),
993 Emc,
994 mcTrack.GetMass()};
995 for (int iPar = 0; iPar < 7; iPar++) {
996 Double_t error = Daughter.GetCovariance(iPar, iPar);
997 if (error < 0.) {
998 error = 1.e20;
999 }
1000 error = TMath::Sqrt(error);
1001 res[iPar] = Daughter.GetParameter(iPar) - mcParam[iPar];
1002 if (fabs(error) > 1.e-20)
1003 pull[iPar] = res[iPar] / error;
1004 }
1005 res[7] = M - mcParam[7];
1006 if (fabs(ErrM) > 1.e-20)
1007 pull[7] = res[7] / ErrM;
1008
1009 for (int iPar = 0; iPar < 8; iPar++) {
1010 hFitDaughtersQA[iParticle][iPar]->Fill(res[iPar]);
1011 hFitDaughtersQA[iParticle][iPar + 8]->Fill(pull[iPar]);
1012 }
1013 }
1014 // for(int iD=0; iD<mcPart.NDaughters(); ++iD)
1015 // {
1016 // int mcDaughterId = mcPart.GetDaughterIds()[iD];
1017 // if(!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg()) continue;
1018 // CbmMCTrack &mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
1019 // int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
1020 //
1021 // CbmKFParticleInterface Daughter1;
1022 // CbmKFParticle_simd Daughter = vRParticles[recDaughterId];
1023 // fvec M,ErrM;
1024 // Daughter.GetMass(M,ErrM);
1025 //
1026 // fvec decayVtx[3] = {mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ()};
1027 // Daughter1.Extrapolate(&Daughter, Daughter.r , Daughter1.GetDStoPoint(&Daughter,decayVtx));
1028 //
1029 // fvec Emc = sqrt(mcTrack.GetP()*mcTrack.GetP() + mcTrack.GetMass()*mcTrack.GetMass());
1030 // fvec res[8] = {0},
1031 // pull[8] = {0},
1032 // mcParam[8] = { mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ(),
1033 // mcTrack.GetPx(), mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass() };
1034 // for(int iPar=0; iPar < 7; iPar++ )
1035 // {
1036 // fvec error = Daughter.C[Daughter1.IJ(iPar, iPar)];
1037 // if(error[0] < 0.) { error = 1.e20;}
1038 // error = TMath::Sqrt(error[0]);
1039 // res[iPar] = Daughter.r[iPar] - mcParam[iPar];
1040 // if(fabs(error)[0] > 1.e-20) pull[iPar] = res[iPar]/error;
1041 // }
1042 // res[7] = M - mcParam[7];
1043 // if(fabs(ErrM)[0] > 1.e-20) pull[7] = res[7]/ErrM;
1044 //
1045 // for(int iPar=0; iPar < 8; iPar++ )
1046 // {
1047 // hFitDaughtersQA[iParticle][iPar]->Fill(res[iPar][0]);
1048 // hFitDaughtersQA[iParticle][iPar+8]->Fill(pull[iPar][0]);
1049 // }
1050 // }
1051 }
1052 // Find MC parameters of the primary vertex
1053 float mcPVx[3] = {0.f};
1054 for (unsigned iMC = 0; iMC < vMCTracks.size(); ++iMC) {
1055 if (vMCTracks[iMC].IsPrimary()) {
1056 mcPVx[0] = vMCTracks[iMC].x;
1057 mcPVx[1] = vMCTracks[iMC].y;
1058 mcPVx[2] = vMCTracks[iMC].z;
1059 break;
1060 }
1061 }
1062 float dRPVr[3] = {(float)(PF->GetPV()->GetRefX() - mcPVx[0]), (float)(PF->GetPV()->GetRefY() - mcPVx[1]),
1063 (float)(PF->GetPV()->GetRefZ() - mcPVx[2])};
1064 float dRPVp[3] = {(float)(dRPVr[0] / sqrt(PF->GetPV()->GetCovMatrix()[0])),
1065 (float)(dRPVr[1] / sqrt(PF->GetPV()->GetCovMatrix()[2])),
1066 (float)(dRPVr[2] / sqrt(PF->GetPV()->GetCovMatrix()[5]))};
1067 for (unsigned int iHPV = 0; iHPV < 3; ++iHPV)
1068 hPVFitQa[iHPV]->Fill(dRPVr[iHPV]);
1069 for (unsigned int iHPV = 3; iHPV < 6; ++iHPV)
1070 hPVFitQa[iHPV]->Fill(dRPVp[iHPV - 3]);
1071
1072} // void CbmL1::ParticlesEfficienciesPerformance()
1073
1074void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRef on match data in CbmL1**Track classes
1075{
1076
1077 // CbmKF &KF = *CbmKF::Instance();
1078
1079 static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom, *p_eff_prim_vs_theta,
1080 *p_eff_all_vs_pt, *p_eff_prim_vs_pt, *p_eff_all_vs_nhits, *p_eff_prim_vs_nhits, *p_eff_sec_vs_nhits;
1081
1082 static TH1F *h_reg_MCmom, *h_acc_MCmom, *h_reco_MCmom, *h_ghost_Rmom;
1083 static TH1F *h_reg_prim_MCmom, *h_acc_prim_MCmom, *h_reco_prim_MCmom;
1084 static TH1F *h_reg_sec_MCmom, *h_acc_sec_MCmom, *h_reco_sec_MCmom;
1085
1086 static TH1F* h_acc_mom_short123s;
1087
1088 static TH1F *h_reg_mom_prim, *h_reg_mom_sec, *h_reg_nhits_prim, *h_reg_nhits_sec;
1089 static TH1F *h_acc_mom_prim, *h_acc_mom_sec, *h_acc_nhits_prim, *h_acc_nhits_sec;
1090 static TH1F *h_reco_mom_prim, *h_reco_mom_sec, *h_reco_nhits_prim, *h_reco_nhits_sec;
1091 static TH1F *h_rest_mom_prim, *h_rest_mom_sec, *h_rest_nhits_prim, *h_rest_nhits_sec;
1092
1093 // static TH1F *h_hit_density[10];
1094
1095 static TH1F *h_ghost_mom, *h_ghost_nhits, *h_ghost_fstation, *h_ghost_chi2, *h_ghost_prob, *h_ghost_tx, *h_ghost_ty;
1096 static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station, *h_reco_chi2, *h_reco_prob, *h_rest_prob,
1097 *h_reco_clean, *h_reco_time;
1098 static TProfile *h_reco_timeNtr, *h_reco_timeNhit;
1099 static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit;
1100 static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;
1101
1102 static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
1103
1104 static TH1F *h_mcp, *h_nmchits; //, *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
1105
1106 // static TH2F *P_vs_P ;
1107
1108 static TH2F *h2_vertex, *h2_prim_vertex, *h2_sec_vertex;
1109 // static TH3F *h3_vertex, *h3_prim_vertex, *h3_sec_vertex;
1110
1111 static TH2F *h2_reg_nhits_vs_mom_prim, *h2_reg_nhits_vs_mom_sec, *h2_reg_fstation_vs_mom_prim,
1112 *h2_reg_fstation_vs_mom_sec, *h2_reg_lstation_vs_fstation_prim, *h2_reg_lstation_vs_fstation_sec;
1113 static TH2F *h2_acc_nhits_vs_mom_prim, *h2_acc_nhits_vs_mom_sec, *h2_acc_fstation_vs_mom_prim,
1114 *h2_acc_fstation_vs_mom_sec, *h2_acc_lstation_vs_fstation_prim, *h2_acc_lstation_vs_fstation_sec;
1115 static TH2F *h2_reco_nhits_vs_mom_prim, *h2_reco_nhits_vs_mom_sec, *h2_reco_fstation_vs_mom_prim,
1116 *h2_reco_fstation_vs_mom_sec, *h2_reco_lstation_vs_fstation_prim, *h2_reco_lstation_vs_fstation_sec;
1117 static TH2F *h2_ghost_nhits_vs_mom, *h2_ghost_fstation_vs_mom, *h2_ghost_lstation_vs_fstation;
1118 static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec, *h2_rest_fstation_vs_mom_prim,
1119 *h2_rest_fstation_vs_mom_sec, *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
1120
1121 static bool first_call = 1;
1122
1123 if (first_call) {
1124 first_call = 0;
1125
1126 TDirectory* curdir = gDirectory;
1127 gDirectory = histodir;
1128 gDirectory->cd("L1");
1129
1130 p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
1131 p_eff_prim_vs_mom =
1132 new TProfile("p_eff_prim_vs_mom", "Primary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
1133 p_eff_sec_vs_mom =
1134 new TProfile("p_eff_sec_vs_mom", "Secondary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
1135 p_eff_d0_vs_mom =
1136 new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum", 150, 0.0, 15.0, 0.0, 100.0);
1137 p_eff_prim_vs_theta =
1138 new TProfile("p_eff_prim_vs_theta", "All Primary Set Efficiency vs Theta", 100, 0.0, 30.0, 0.0, 100.0);
1139 p_eff_all_vs_pt = new TProfile("p_eff_all_vs_pt", "AllSet Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
1140 p_eff_prim_vs_pt = new TProfile("p_eff_prim_vs_pt", "Primary Set Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
1141
1142 p_eff_all_vs_nhits =
1143 new TProfile("p_eff_all_vs_nhits", "AllSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
1144 p_eff_prim_vs_nhits =
1145 new TProfile("p_eff_prim_vs_nhits", "PrimSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
1146 p_eff_sec_vs_nhits =
1147 new TProfile("p_eff_sec_vs_nhits", "SecSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
1148
1149 h_reg_MCmom = new TH1F("h_reg_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1150 h_acc_MCmom = new TH1F("h_acc_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1151 h_reco_MCmom = new TH1F("h_reco_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1152 h_ghost_Rmom = new TH1F("h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0);
1153 h_reg_prim_MCmom = new TH1F("h_reg_prim_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1154 h_acc_prim_MCmom = new TH1F("h_acc_prim_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1155 h_reco_prim_MCmom = new TH1F("h_reco_prim_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1156 h_reg_sec_MCmom = new TH1F("h_reg_sec_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1157 h_acc_sec_MCmom = new TH1F("h_acc_sec_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1158 h_reco_sec_MCmom = new TH1F("h_reco_sec_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1159
1160 h_acc_mom_short123s =
1161 new TH1F("h_acc_mom_short123s", "Momentum of accepted tracks with 3 hits on first stations", 500, 0.0, 5.0);
1162
1163 h_reg_mom_prim = new TH1F("h_reg_mom_prim", "Momentum of registered primary tracks", 500, 0.0, 5.0);
1164 h_reg_mom_sec = new TH1F("h_reg_mom_sec", "Momentum of registered secondary tracks", 500, 0.0, 5.0);
1165 h_acc_mom_prim = new TH1F("h_acc_mom_prim", "Momentum of accepted primary tracks", 500, 0.0, 5.0);
1166 h_acc_mom_sec = new TH1F("h_acc_mom_sec", "Momentum of accepted secondary tracks", 500, 0.0, 5.0);
1167 h_reco_mom_prim = new TH1F("h_reco_mom_prim", "Momentum of reconstructed primary tracks", 500, 0.0, 5.0);
1168 h_reco_mom_sec = new TH1F("h_reco_mom_sec", "Momentum of reconstructed secondary tracks", 500, 0.0, 5.0);
1169 h_rest_mom_prim = new TH1F("h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
1170 h_rest_mom_sec = new TH1F("h_rest_mom_sec", "Momentum of not found secondary tracks", 500, 0.0, 5.0);
1171
1172 h_reg_nhits_prim = new TH1F("h_reg_nhits_prim", "Number of hits in registered primary track", 51, -0.1, 10.1);
1173 h_reg_nhits_sec = new TH1F("h_reg_nhits_sec", "Number of hits in registered secondary track", 51, -0.1, 10.1);
1174 h_acc_nhits_prim = new TH1F("h_acc_nhits_prim", "Number of hits in accepted primary track", 51, -0.1, 10.1);
1175 h_acc_nhits_sec = new TH1F("h_acc_nhits_sec", "Number of hits in accepted secondary track", 51, -0.1, 10.1);
1176 h_reco_nhits_prim =
1177 new TH1F("h_reco_nhits_prim", "Number of hits in reconstructed primary track", 51, -0.1, 10.1);
1178 h_reco_nhits_sec =
1179 new TH1F("h_reco_nhits_sec", "Number of hits in reconstructed secondary track", 51, -0.1, 10.1);
1180 h_rest_nhits_prim = new TH1F("h_rest_nhits_prim", "Number of hits in not found primary track", 51, -0.1, 10.1);
1181 h_rest_nhits_sec = new TH1F("h_rest_nhits_sec", "Number of hits in not found secondary track", 51, -0.1, 10.1);
1182
1183 h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
1184 h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
1185 h_ghost_fstation = new TH1F("h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
1186 h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
1187 h_ghost_prob = new TH1F("h_ghost_prob", "Prob of ghost track", 505, 0., 1.01);
1188 h_ghost_r = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0);
1189 h_ghost_tx = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
1190 h_ghost_ty = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);
1191
1192 h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
1193 h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
1194 h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
1195 h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
1196 h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 505, 0., 1.01);
1197 h_rest_prob = new TH1F("h_rest_prob", "Prob of reco rest track", 505, 0., 1.01);
1198 h_reco_clean = new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5);
1199 h_reco_time = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
1200 h_reco_timeNtr = new TProfile("h_reco_timeNtr", "CA Track Finder Time (s/ev) vs N Tracks", 200, 0.0, 1000.0);
1201 h_reco_timeNhit = new TProfile("h_reco_timeNhit", "CA Track Finder Time (s/ev) vs N Hits", 200, 0.0, 30000.0);
1202 h_reco_fakeNtr = new TProfile("h_reco_fakeNtr", "N Fake hits vs N Tracks", 200, 0.0, 1000.0);
1203 h_reco_fakeNhit = new TProfile("h_reco_fakeNhit", "N Fake hits vs N Hits", 200, 0.0, 30000.0);
1204
1205 h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);
1206
1207 // h_hit_density[0] = new TH1F("h_hit_density0", "Hit density station 1", 50, 0.0, 5.00);
1208 // h_hit_density[1] = new TH1F("h_hit_density1", "Hit density station 2", 100, 0.0, 10.00);
1209 // h_hit_density[2] = new TH1F("h_hit_density2", "Hit density station 3", 140, 0.0, 13.99);
1210 // h_hit_density[3] = new TH1F("h_hit_density3", "Hit density station 4", 180, 0.0, 18.65);
1211 // h_hit_density[4] = new TH1F("h_hit_density4", "Hit density station 5", 230, 0.0, 23.32);
1212 // h_hit_density[5] = new TH1F("h_hit_density5", "Hit density station 6", 270, 0.0, 27.98);
1213 // h_hit_density[6] = new TH1F("h_hit_density6", "Hit density station 7", 340, 0.0, 34.97);
1214 // h_hit_density[7] = new TH1F("h_hit_density7", "Hit density station 8", 460, 0.0, 46.63);
1215 // h_hit_density[8] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1216 // h_hit_density[9] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1217 // h_hit_density[10] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1218
1219 h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
1220 h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);
1221
1222 h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);
1223
1224 h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
1225 h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
1226 h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
1227 h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
1228 h_notfound_tx = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
1229 h_notfound_ty = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);
1230
1231 // h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
1232 // h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
1233 // VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
1234 // VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
1235 h_mcp = new TH1F("h_mcp", "MC P ", 500, 0.0, 5.0);
1236 // MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
1237 // MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
1238 // MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
1239 h_nmchits = new TH1F("h_nmchits", "N STS hits on MC track", 50, 0.0, 10.0);
1240
1241 // P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);
1242
1243 h2_vertex = new TH2F("h2_vertex", "2D vertex distribution", 105, -5., 100., 100, -50., 50.);
1244 h2_prim_vertex = new TH2F("h2_primvertex", "2D primary vertex distribution", 105, -5., 100., 100, -50., 50.);
1245 h2_sec_vertex = new TH2F("h2_sec_vertex", "2D secondary vertex distribution", 105, -5., 100., 100, -50., 50.);
1246
1247 // h3_vertex = new TH3F("h3_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
1248 // h3_prim_vertex = new TH3F("h3_primvertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100,
1249 // -50., 50.); h3_sec_vertex = new TH3F("h3_sec_vertex", "3D vertex distribution", 20, -5., 100., 100,
1250 // -50., 50., 100, -50., 50.);
1251
1252 h2_reg_nhits_vs_mom_prim = new TH2F("h2_reg_nhits_vs_mom_prim", "NHits vs. Momentum (Reg. Primary Tracks)", 51,
1253 -0.05, 5.05, 11, -0.5, 10.5);
1254 h2_reg_nhits_vs_mom_sec = new TH2F("h2_reg_nhits_vs_mom_sec", "NHits vs. Momentum (Reg. Secondary Tracks)", 51,
1255 -0.05, 5.05, 11, -0.5, 10.5);
1256 h2_acc_nhits_vs_mom_prim = new TH2F("h2_acc_nhits_vs_mom_prim", "NHits vs. Momentum (Acc. Primary Tracks)", 51,
1257 -0.05, 5.05, 11, -0.5, 10.5);
1258 h2_acc_nhits_vs_mom_sec = new TH2F("h2_acc_nhits_vs_mom_sec", "NHits vs. Momentum (Acc. Secondary Tracks)", 51,
1259 -0.05, 5.05, 11, -0.5, 10.5);
1260 h2_reco_nhits_vs_mom_prim = new TH2F("h2_reco_nhits_vs_mom_prim", "NHits vs. Momentum (Reco Primary Tracks)",
1261 51, -0.05, 5.05, 11, -0.5, 10.5);
1262 h2_reco_nhits_vs_mom_sec = new TH2F("h2_reco_nhits_vs_mom_sec", "NHits vs. Momentum (Reco Secondary Tracks)",
1263 51, -0.05, 5.05, 11, -0.5, 10.5);
1264 h2_ghost_nhits_vs_mom =
1265 new TH2F("h2_ghost_nhits_vs_mom", "NHits vs. Momentum (Ghost Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1266 h2_rest_nhits_vs_mom_prim = new TH2F("h2_rest_nhits_vs_mom_prim", "NHits vs. Momentum (!Found Primary Tracks)",
1267 51, -0.05, 5.05, 11, -0.5, 10.5);
1268 h2_rest_nhits_vs_mom_sec = new TH2F("h2_rest_nhits_vs_mom_sec", "NHits vs. Momentum (!Found Secondary Tracks)",
1269 51, -0.05, 5.05, 11, -0.5, 10.5);
1270
1271 h2_reg_fstation_vs_mom_prim =
1272 new TH2F("h2_reg_fstation_vs_mom_prim", "First Station vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05,
1273 11, -0.5, 10.5);
1274 h2_reg_fstation_vs_mom_sec =
1275 new TH2F("h2_reg_fstation_vs_mom_sec", "First Station vs. Momentum (Reg. Secondary Tracks)", 51, -0.05,
1276 5.05, 11, -0.5, 10.5);
1277 h2_acc_fstation_vs_mom_prim =
1278 new TH2F("h2_acc_fstation_vs_mom_prim", "First Station vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05,
1279 11, -0.5, 10.5);
1280 h2_acc_fstation_vs_mom_sec =
1281 new TH2F("h2_acc_fstation_vs_mom_sec", "First Station vs. Momentum (Acc. Secondary Tracks)", 51, -0.05,
1282 5.05, 11, -0.5, 10.5);
1283 h2_reco_fstation_vs_mom_prim =
1284 new TH2F("h2_reco_fstation_vs_mom_prim", "First Station vs. Momentum (Reco Primary Tracks)", 51, -0.05,
1285 5.05, 11, -0.5, 10.5);
1286 h2_reco_fstation_vs_mom_sec =
1287 new TH2F("h2_reco_fstation_vs_mom_sec", "First Station vs. Momentum (Reco Secondary Tracks)", 51, -0.05,
1288 5.05, 11, -0.5, 10.5);
1289 h2_ghost_fstation_vs_mom = new TH2F("h2_ghost_fstation_vs_mom", "First Station vs. Momentum (Ghost Tracks)", 51,
1290 -0.05, 5.05, 11, -0.5, 10.5);
1291 h2_rest_fstation_vs_mom_prim =
1292 new TH2F("h2_rest_fstation_vs_mom_prim", "First Station vs. Momentum (!Found Primary Tracks)", 51, -0.05,
1293 5.05, 11, -0.5, 10.5);
1294 h2_rest_fstation_vs_mom_sec =
1295 new TH2F("h2_rest_fstation_vs_mom_sec", "First Station vs. Momentum (!Found Secondary Tracks)", 51, -0.05,
1296 5.05, 11, -0.5, 10.5);
1297
1298 h2_reg_lstation_vs_fstation_prim =
1299 new TH2F("h2_reg_lstation_vs_fstation_prim", "Last vs. First Station (Reg. Primary Tracks)", 11, -0.5, 10.5,
1300 11, -0.5, 10.5);
1301 h2_reg_lstation_vs_fstation_sec =
1302 new TH2F("h2_reg_lstation_vs_fstation_sec", "Last vs. First Station (Reg. Secondary Tracks)", 11, -0.5,
1303 10.5, 11, -0.5, 10.5);
1304 h2_acc_lstation_vs_fstation_prim =
1305 new TH2F("h2_acc_lstation_vs_fstation_prim", "Last vs. First Station (Acc. Primary Tracks)", 11, -0.5, 10.5,
1306 11, -0.5, 10.5);
1307 h2_acc_lstation_vs_fstation_sec =
1308 new TH2F("h2_acc_lstation_vs_fstation_sec", "Last vs. First Station (Acc. Secondary Tracks)", 11, -0.5,
1309 10.5, 11, -0.5, 10.5);
1310 h2_reco_lstation_vs_fstation_prim =
1311 new TH2F("h2_reco_lstation_vs_fstation_prim", "Last vs. First Station (Reco Primary Tracks)", 11, -0.5,
1312 10.5, 11, -0.5, 10.5);
1313 h2_reco_lstation_vs_fstation_sec =
1314 new TH2F("h2_reco_lstation_vs_fstation_sec", "Last vs. First Station (Reco Secondary Tracks)", 11, -0.5,
1315 10.5, 11, -0.5, 10.5);
1316 h2_ghost_lstation_vs_fstation = new TH2F(
1317 "h2_ghost_lstation_vs_fstation", "Last vs. First Station (Ghost Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1318 h2_rest_lstation_vs_fstation_prim =
1319 new TH2F("h2_rest_lstation_vs_fstation_prim", "Last vs. First Station (!Found Primary Tracks)", 11, -0.5,
1320 10.5, 11, -0.5, 10.5);
1321 h2_rest_lstation_vs_fstation_sec =
1322 new TH2F("h2_rest_lstation_vs_fstation_sec", "Last vs. First Station (!Found Secondary Tracks)", 11, -0.5,
1323 10.5, 11, -0.5, 10.5);
1324
1325 // maindir->cd();
1326
1327 // ----- Create list of all histogram pointers
1328
1329 gDirectory = curdir;
1330
1331 } // first_call
1332
1333 static int NEvents = 0;
1334 if (NEvents > 0) {
1335 h_reg_MCmom->Scale(NEvents);
1336 h_acc_MCmom->Scale(NEvents);
1337 h_reco_MCmom->Scale(NEvents);
1338 h_ghost_Rmom->Scale(NEvents);
1339 h_reg_prim_MCmom->Scale(NEvents);
1340 h_acc_prim_MCmom->Scale(NEvents);
1341 h_reco_prim_MCmom->Scale(NEvents);
1342 h_reg_sec_MCmom->Scale(NEvents);
1343 h_acc_sec_MCmom->Scale(NEvents);
1344 h_reco_sec_MCmom->Scale(NEvents);
1345 }
1346 NEvents++;
1347
1348 // //hit density calculation: h_hit_density[10]
1349 //
1350 // for (vector<CbmL1HitStore>::iterator hIt = vHitStore.begin(); hIt != vHitStore.end(); ++hIt){
1351 // float x = hIt->x;
1352 // float y = hIt->y;
1353 // float r = sqrt(x*x+y*y);
1354 // h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
1355 // }
1356
1357 //
1358 for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt) {
1359 CbmL1Track* prtra = &(*rtraIt);
1360 if ((prtra->StsHits).size() < 1)
1361 continue;
1362 { // fill histos
1363 if (fabs(prtra->T[4]) > 1.e-10)
1364 h_reco_mom->Fill(fabs(1.0 / prtra->T[4]));
1365 h_reco_nhits->Fill((prtra->StsHits).size());
1366 CbmL1HitStore& mh = vHitStore[prtra->StsHits[0]];
1367 h_reco_station->Fill(mh.iStation);
1368 }
1369
1370 h_reco_clean->Fill(prtra->GetMaxPurity());
1371
1372 if (prtra->NDF > 0) {
1373 if (prtra->IsGhost()) {
1374 h_ghost_chi2->Fill(prtra->chi2 / prtra->NDF);
1375 h_ghost_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1376 } else {
1377 if (prtra->GetMCTrack()[0].IsReconstructable()) {
1378 h_reco_chi2->Fill(prtra->chi2 / prtra->NDF);
1379 h_reco_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1380 } else {
1381 // h_rest_chi2->Fill(prtra->chi2/prtra->NDF);
1382 h_rest_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1383 }
1384 }
1385 }
1386
1387 // fill ghost histos
1388 if (prtra->IsGhost()) {
1389 if (fabs(prtra->T[4]) > 1.e-10) {
1390 h_ghost_mom->Fill(fabs(1.0 / prtra->T[4]));
1391 h_ghost_Rmom->Fill(fabs(1.0 / prtra->T[4]));
1392 }
1393 h_ghost_nhits->Fill((prtra->StsHits).size());
1394 CbmL1HitStore& h1 = vHitStore[prtra->StsHits[0]];
1395 CbmL1HitStore& h2 = vHitStore[prtra->StsHits[1]];
1396 h_ghost_fstation->Fill(h1.iStation);
1397 h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
1398 double z1 = algo->vStations[h1.iStation].z[0];
1399 double z2 = algo->vStations[h2.iStation].z[0];
1400 if (fabs(z2 - z1) > 1.e-4) {
1401 h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
1402 h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
1403 }
1404
1405 if (fabs(prtra->T[4]) > 1.e-10)
1406 h2_ghost_nhits_vs_mom->Fill(fabs(1.0 / prtra->T[4]), (prtra->StsHits).size());
1407 CbmL1HitStore& hf = vHitStore[prtra->StsHits[0]];
1408 CbmL1HitStore& hl = vHitStore[prtra->StsHits[(prtra->StsHits).size() - 1]];
1409 if (fabs(prtra->T[4]) > 1.e-10)
1410 h2_ghost_fstation_vs_mom->Fill(fabs(1.0 / prtra->T[4]), hf.iStation + 1);
1411 if (hl.iStation >= hf.iStation)
1412 h2_ghost_lstation_vs_fstation->Fill(hf.iStation + 1, hl.iStation + 1);
1413 }
1414
1415 } // for reco tracks
1416
1417 int mc_total = 0; // total amount of reconstructable mcTracks
1418 for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++) {
1419 CbmL1MCTrack& mtra = *(mtraIt);
1420 // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
1421
1422 // No Sts hits?
1423 int nmchits = mtra.StsHits.size();
1424 if (nmchits == 0)
1425 continue;
1426
1427 double momentum = mtra.p;
1428 double pt = sqrt(mtra.px * mtra.px + mtra.py * mtra.py);
1429 double theta = acos(mtra.pz / mtra.p) * 180 / 3.1415;
1430
1431 h_mcp->Fill(momentum);
1432 h_nmchits->Fill(nmchits);
1433
1434 int nSta = mtra.NStations();
1435
1436 CbmL1HitStore& fh = vHitStore[*(mtra.StsHits.begin())];
1437 CbmL1HitStore& lh = vHitStore[*(mtra.StsHits.rbegin())];
1438
1439 h_reg_MCmom->Fill(momentum);
1440 if (mtra.IsPrimary()) {
1441 h_reg_mom_prim->Fill(momentum);
1442 h_reg_prim_MCmom->Fill(momentum);
1443 h_reg_nhits_prim->Fill(nSta);
1444 h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
1445 h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1446 if (lh.iStation >= fh.iStation)
1447 h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
1448 } else {
1449 h_reg_mom_sec->Fill(momentum);
1450 h_reg_sec_MCmom->Fill(momentum);
1451 h_reg_nhits_sec->Fill(nSta);
1452 if (momentum > 0.01) {
1453 h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
1454 h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1455 if (lh.iStation >= fh.iStation)
1456 h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1457 }
1458 }
1459
1460 if (mtra.IsAdditional()) {
1461 h_acc_mom_short123s->Fill(momentum);
1462 }
1463
1464 if (!mtra.IsReconstructable())
1465 continue;
1466 mc_total++;
1467
1468 h_acc_MCmom->Fill(momentum);
1469 if (mtra.IsPrimary()) {
1470 h_acc_mom_prim->Fill(momentum);
1471 h_acc_prim_MCmom->Fill(momentum);
1472 h_acc_nhits_prim->Fill(nSta);
1473 h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
1474 h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1475 if (lh.iStation >= fh.iStation)
1476 h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
1477 } else {
1478 h_acc_mom_sec->Fill(momentum);
1479 h_acc_sec_MCmom->Fill(momentum);
1480 h_acc_nhits_sec->Fill(nSta);
1481 if (momentum > 0.01) {
1482 h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
1483 h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1484 if (lh.iStation >= fh.iStation)
1485 h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1486 }
1487 }
1488
1489 // vertex distribution of primary and secondary tracks
1490 h2_vertex->Fill(mtra.z, mtra.y);
1491 // h3_vertex->Fill(mtra.z, mtra.x, mtra.y);
1492 if (mtra.mother_ID < 0) { // primary
1493 h2_prim_vertex->Fill(mtra.z, mtra.y);
1494 // h3_prim_vertex->Fill(mtra.z, mtra.x, mtra.y);
1495 } else {
1496 h2_sec_vertex->Fill(mtra.z, mtra.y);
1497 // h3_sec_vertex->Fill(mtra.z, mtra.x, mtra.y);
1498 }
1499
1500 int iph = mtra.StsHits[0];
1501 CbmL1HitStore& ph = vHitStore[iph];
1502
1503 h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1504 if (fabs(mtra.pz) > 1.e-8) {
1505 h_tx->Fill(mtra.px / mtra.pz);
1506 h_ty->Fill(mtra.py / mtra.pz);
1507 }
1508
1509 // reconstructed tracks
1510 bool reco = (mtra.IsReconstructed());
1511 int nMCHits = mtra.NStations();
1512
1513 if (reco) {
1514 p_eff_all_vs_mom->Fill(momentum, 100.0);
1515 p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
1516 p_eff_all_vs_pt->Fill(pt, 100.0);
1517 h_reco_MCmom->Fill(momentum);
1518 if (mtra.mother_ID < 0) { // primary
1519 p_eff_prim_vs_mom->Fill(momentum, 100.0);
1520 p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
1521 p_eff_prim_vs_pt->Fill(pt, 100.0);
1522 p_eff_prim_vs_theta->Fill(theta, 100.0);
1523 } else {
1524 p_eff_sec_vs_mom->Fill(momentum, 100.0);
1525 p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
1526 }
1527 if (mtra.mother_ID < 0) { // primary
1528 h_reco_mom_prim->Fill(momentum);
1529 h_reco_prim_MCmom->Fill(momentum);
1530 h_reco_nhits_prim->Fill(nSta);
1531 h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
1532 h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1533 if (lh.iStation >= fh.iStation)
1534 h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
1535 } else {
1536 h_reco_mom_sec->Fill(momentum);
1537 h_reco_sec_MCmom->Fill(momentum);
1538 h_reco_nhits_sec->Fill(nSta);
1539 if (momentum > 0.01) {
1540 h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
1541 h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1542 if (lh.iStation >= fh.iStation)
1543 h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1544 }
1545 }
1546 } else {
1547 h_notfound_mom->Fill(momentum);
1548 p_eff_all_vs_mom->Fill(momentum, 0.0);
1549 p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
1550 p_eff_all_vs_pt->Fill(pt, 0.0);
1551 if (mtra.mother_ID < 0) { // primary
1552 p_eff_prim_vs_mom->Fill(momentum, 0.0);
1553 p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
1554 p_eff_prim_vs_pt->Fill(pt, 0.0);
1555 p_eff_prim_vs_theta->Fill(theta, 0.0);
1556 } else {
1557 p_eff_sec_vs_mom->Fill(momentum, 0.0);
1558 p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
1559 }
1560 if (mtra.mother_ID < 0) { // primary
1561 h_rest_mom_prim->Fill(momentum);
1562 h_rest_nhits_prim->Fill(nSta);
1563 h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
1564 h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1565 if (lh.iStation >= fh.iStation)
1566 h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
1567 } else {
1568 h_rest_mom_sec->Fill(momentum);
1569 h_rest_nhits_sec->Fill(nSta);
1570 if (momentum > 0.01) {
1571 h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
1572 h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1573 if (lh.iStation >= fh.iStation)
1574 h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1575 }
1576 }
1577 }
1578
1579 // killed tracks. At least one hit of they belong to some recoTrack
1580 // bool killed = 0;
1581 if (!reco) {
1582 h_notfound_nhits->Fill(nmchits);
1583 h_notfound_station->Fill(ph.iStation);
1584 h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1585
1586 if (mtra.Points.size() > 0) {
1587 CbmL1MCPoint& pMC = vMCPoints[mtra.Points[0]];
1588 h_notfound_tx->Fill(pMC.px / pMC.pz);
1589 h_notfound_ty->Fill(pMC.py / pMC.pz);
1590 }
1591
1592 // CbmL1HitStore &ph21 = vHitStore[mtra.StsHits[0]];
1593 // CbmL1HitStore &ph22 = vHitStore[mtra.StsHits[1]];
1594
1595 // double z21 = algo->vStations[ph21.iStation].z[0];
1596 // double z22 = algo->vStations[ph22.iStation].z[0];
1597 // if( fabs(z22-z21)>1.e-4 ){
1598 // h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
1599 // h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
1600 // }
1601
1602 // if( mtra.IsDisturbed() ) killed = 1;
1603 }
1604
1605 if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
1606 h_reco_d0_mom->Fill(momentum);
1607 if (reco)
1608 p_eff_d0_vs_mom->Fill(momentum, 100.0);
1609 else
1610 p_eff_d0_vs_mom->Fill(momentum, 0.0);
1611 }
1612
1613 } // for mcTracks
1614
1615 int NFakes = 0;
1616 for (unsigned int ih = 0; ih < algo->vStsHits.size(); ih++) {
1617 int iMC = vHitMCRef[ih]; // TODO2: adapt to linking
1618 if (iMC < 0)
1619 NFakes++;
1620 }
1621
1622 h_reco_time->Fill(algo->CATime);
1623 h_reco_timeNtr->Fill(mc_total, algo->CATime);
1624 h_reco_timeNhit->Fill(algo->vStsHits.size(), algo->CATime);
1625
1626 h_reco_fakeNtr->Fill(mc_total, NFakes);
1627 h_reco_fakeNhit->Fill(algo->vStsHits.size() - NFakes, NFakes);
1628
1629 h_reg_MCmom->Scale(1.f / NEvents);
1630 h_acc_MCmom->Scale(1.f / NEvents);
1631 h_reco_MCmom->Scale(1.f / NEvents);
1632 h_ghost_Rmom->Scale(1.f / NEvents);
1633 h_reg_prim_MCmom->Scale(1.f / NEvents);
1634 h_acc_prim_MCmom->Scale(1.f / NEvents);
1635 h_reco_prim_MCmom->Scale(1.f / NEvents);
1636 h_reg_sec_MCmom->Scale(1.f / NEvents);
1637 h_acc_sec_MCmom->Scale(1.f / NEvents);
1638 h_reco_sec_MCmom->Scale(1.f / NEvents);
1639
1640} // void CbmL1::HistoPerformance()
1641
1642void CbmL1::TrackFitPerformance()
1643{
1644 const int Nh_fit = 12;
1645 static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit], *h_fit_chi2;
1646
1647 static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
1648
1649 static bool first_call = 1;
1650
1651 if (first_call) {
1652 first_call = 0;
1653
1654 TDirectory* currentDir = gDirectory;
1655 gDirectory = histodir;
1656 gDirectory->cd("L1");
1657 gDirectory->mkdir("Fit");
1658 gDirectory->cd("Fit");
1659 {
1660 PRes2D = new TH2F("PRes2D", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1661 PRes2DPrim = new TH2F("PRes2DPrim", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1662 PRes2DSec = new TH2F("PRes2DSec", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1663
1664 struct
1665 {
1666 const char* name;
1667 const char* title;
1668 Int_t n;
1669 Double_t l, r;
1670 } Table[Nh_fit] = {{"x", "Residual X [#mum]", 100, -45., 45.},
1671 {"y", "Residual Y [#mum]", 100, -450., 450.},
1672 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1673 {"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
1674 {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1},
1675 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1676 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1677 {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1678 {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1679 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1680 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1681 {"QPmc", "MC Q/P ", 100, -10., 10.}};
1682
1683 struct Tab
1684 {
1685 const char* name;
1686 const char* title;
1687 Int_t n;
1688 Double_t l, r;
1689 };
1690 Tab TableVertex[Nh_fit] = {{"x", "Residual X [cm]", 2000, -1., 1.},
1691 {"y", "Residual Y [cm]", 2000, -1., 1.},
1692 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1693 {"ty", "Residual Ty [mrad]", 100, -2., 2.},
1694 {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1},
1695 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1696 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1697 {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1698 {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1699 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1700 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1701 {"QPmc", "MC Q/P ", 100, -10., 10.}};
1702
1703 for (int i = 0; i < Nh_fit; i++) {
1704 char n[225], t[255];
1705 sprintf(n, "fst_%s", Table[i].name);
1706 sprintf(t, "First point %s", Table[i].title);
1707 h_fit[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1708 sprintf(n, "lst_%s", Table[i].name);
1709 sprintf(t, "Last point %s", Table[i].title);
1710 h_fitL[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1711 sprintf(n, "svrt_%s", TableVertex[i].name);
1712 sprintf(t, "Secondary vertex point %s", TableVertex[i].title);
1713 h_fitSV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1714 sprintf(n, "pvrt_%s", TableVertex[i].name);
1715 sprintf(t, "Primary vertex point %s", TableVertex[i].title);
1716 h_fitPV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1717 }
1718 h_fit_chi2 = new TH1F("h_fit_chi2", "Chi2/NDF", 50, -0.5, 10.0);
1719 }
1720 gDirectory = currentDir;
1721 } // if first call
1722
1723 for (vector<CbmL1Track>::iterator it = vRTracks.begin(); it != vRTracks.end(); ++it) {
1724
1725 if (it->IsGhost())
1726 continue;
1727 { // first hit
1728 int iMC = vHitMCRef[it->StsHits.front()]; // TODO2: adapt to linking
1729 if (iMC < 0)
1730 continue;
1731 CbmL1MCPoint& mc = vMCPoints[iMC];
1732 // if( !( mc.pdg == -11 && mc.mother_ID == -1 ) ) continue; // electrons only
1733 h_fit[0]->Fill((mc.x - it->T[0]) * 1.e4);
1734 h_fit[1]->Fill((mc.y - it->T[1]) * 1.e4);
1735 h_fit[2]->Fill((mc.px / mc.pz - it->T[2]) * 1.e3);
1736 h_fit[3]->Fill((mc.py / mc.pz - it->T[3]) * 1.e3);
1737 h_fit[4]->Fill(it->T[4] / mc.q * mc.p - 1);
1738
1739 PRes2D->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1740
1741 CbmL1MCTrack mcTrack = *(it->GetMCTracks()[0]);
1742 if (mcTrack.IsPrimary())
1743 PRes2DPrim->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1744 else
1745 PRes2DSec->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1746
1747 if (finite(it->C[0]) && it->C[0] > 0)
1748 h_fit[5]->Fill((mc.x - it->T[0]) / sqrt(it->C[0]));
1749 if (finite(it->C[2]) && it->C[2] > 0)
1750 h_fit[6]->Fill((mc.y - it->T[1]) / sqrt(it->C[2]));
1751 if (finite(it->C[5]) && it->C[5] > 0)
1752 h_fit[7]->Fill((mc.px / mc.pz - it->T[2]) / sqrt(it->C[5]));
1753 if (finite(it->C[9]) && it->C[9] > 0)
1754 h_fit[8]->Fill((mc.py / mc.pz - it->T[3]) / sqrt(it->C[9]));
1755 if (finite(it->C[14]) && it->C[14] > 0)
1756 h_fit[9]->Fill((mc.q / mc.p - it->T[4]) / sqrt(it->C[14]));
1757 h_fit[10]->Fill(it->T[4]);
1758 h_fit[11]->Fill(mc.q / mc.p);
1759 }
1760
1761 { // last hit
1762 int iMC = vHitMCRef[it->StsHits.back()]; // TODO2: adapt to linking
1763 if (iMC < 0)
1764 continue;
1765 CbmL1MCPoint& mc = vMCPoints[iMC];
1766 h_fitL[0]->Fill((mc.x - it->TLast[0]) * 1.e4);
1767 h_fitL[1]->Fill((mc.y - it->TLast[1]) * 1.e4);
1768 h_fitL[2]->Fill((mc.px / mc.pz - it->TLast[2]) * 1.e3);
1769 h_fitL[3]->Fill((mc.py / mc.pz - it->TLast[3]) * 1.e3);
1770 h_fitL[4]->Fill(it->T[4] / mc.q * mc.p - 1);
1771 if (finite(it->CLast[0]) && it->CLast[0] > 0)
1772 h_fitL[5]->Fill((mc.x - it->TLast[0]) / sqrt(it->CLast[0]));
1773 if (finite(it->CLast[2]) && it->CLast[2] > 0)
1774 h_fitL[6]->Fill((mc.y - it->TLast[1]) / sqrt(it->CLast[2]));
1775 if (finite(it->CLast[5]) && it->CLast[5] > 0)
1776 h_fitL[7]->Fill((mc.px / mc.pz - it->TLast[2]) / sqrt(it->CLast[5]));
1777 if (finite(it->CLast[9]) && it->CLast[9] > 0)
1778 h_fitL[8]->Fill((mc.py / mc.pz - it->TLast[3]) / sqrt(it->CLast[9]));
1779 if (finite(it->CLast[14]) && it->CLast[14] > 0)
1780 h_fitL[9]->Fill((mc.q / mc.p - it->TLast[4]) / sqrt(it->CLast[14]));
1781 h_fitL[10]->Fill(it->TLast[4]);
1782 h_fitL[11]->Fill(mc.q / mc.p);
1783 }
1784
1785 { // vertex
1786 CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
1787 L1TrackPar trPar(it->T, it->C);
1788
1789 if (mc.mother_ID != -1) { // secondary
1790
1791 { // extrapolate to vertex
1792 L1FieldRegion fld _fvecalignment;
1793 L1FieldValue B[3] _fvecalignment;
1794 float z[3];
1795 for (int ih = 0; ih < 3; ih++) {
1796 const int iMCP = mc.Points[ih];
1797 CbmL1MCPoint& mcP = vMCPoints[iMCP];
1798 L1Station& st = algo->vStations[mcP.iStation];
1799 z[ih] = st.z[0];
1800 st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
1801 };
1802 fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1803
1804 L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1805 // add material
1806 const int fSta = vHitStore[it->StsHits[0]].iStation;
1807 const int dir = int((mc.z - algo->vStations[fSta].z[0]) / fabs(mc.z - algo->vStations[fSta].z[0]));
1808 // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large
1809 // distance
1810 for (int iSta = fSta /*+dir*/;
1811 (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0);
1812 iSta += dir)
1813 {
1814 // cout << iSta << " " << dir << endl;
1815 L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp);
1816 if (iSta + dir == NMvdStations - 1)
1817 L1AddPipeMaterial(trPar, trPar.qp);
1818 }
1819 }
1820 if (mc.z != trPar.z[0])
1821 continue;
1822 // static int good = 0;
1823 // static int bad = 0;
1824 // if (mc.z != trPar.z[0]){
1825 // bad++;
1826 // continue;
1827 // }
1828 // else good++;
1829 // cout << "bad\\good" << bad << " " << good << endl;
1830
1831 // calculate pulls
1832 // h_fitSV[0]->Fill( (mc.x-trPar.x[0]) *1.e4);
1833 // h_fitSV[1]->Fill( (mc.y-trPar.y[0]) *1.e4);
1834 h_fitSV[0]->Fill((mc.x - trPar.x[0]));
1835 h_fitSV[1]->Fill((mc.y - trPar.y[0]));
1836 h_fitSV[2]->Fill((mc.px / mc.pz - trPar.tx[0]) * 1.e3);
1837 h_fitSV[3]->Fill((mc.py / mc.pz - trPar.ty[0]) * 1.e3);
1838 h_fitSV[4]->Fill(trPar.qp[0] / mc.q * mc.p - 1);
1839 if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1840 h_fitSV[5]->Fill((mc.x - trPar.x[0]) / sqrt(trPar.C00[0]));
1841 if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1842 h_fitSV[6]->Fill((mc.y - trPar.y[0]) / sqrt(trPar.C11[0]));
1843 if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1844 h_fitSV[7]->Fill((mc.px / mc.pz - trPar.tx[0]) / sqrt(trPar.C22[0]));
1845 if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1846 h_fitSV[8]->Fill((mc.py / mc.pz - trPar.ty[0]) / sqrt(trPar.C33[0]));
1847 if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1848 h_fitSV[9]->Fill((mc.q / mc.p - trPar.qp[0]) / sqrt(trPar.C44[0]));
1849 h_fitSV[10]->Fill(trPar.qp[0]);
1850 h_fitSV[11]->Fill(mc.q / mc.p);
1851 } else { // primary
1852 if (vHitStore[it->StsHits[0]].iStation != 0)
1853 continue; // only mvd tracks
1854
1855// #define L1EXTRAPOLATE
1856#ifdef L1EXTRAPOLATE
1857 { // extrapolate to vertex
1858 L1FieldRegion fld _fvecalignment;
1859 L1FieldValue B[3], targB _fvecalignment;
1860 float z[3];
1861 for (int ih = 0; ih < 3; ih++) {
1862 const int iMCP = mc.Points[ih];
1863 CbmL1MCPoint& mcP = vMCPoints[iMCP];
1864 L1Station& st = algo->vStations[mcP.iStation];
1865 z[ih] = st.z[0];
1866 // cout << mcP.z << endl;
1867 st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
1868 };
1869 // cout << mc.x << " " << mc.y << " " << mc.z << endl;
1870
1871 // targetFieldSlice->GetFieldValue( mc.x, mc.y, targB );
1872 targB = algo->vtxFieldValue;
1873 fld.Set(targB, 0., B[0], z[0], B[1], z[1]);
1874
1875 L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1876 // L1Extrapolate0(trPar, mc.z, fld);
1877 // add material
1878 const int fSta = vHitStore[it->StsHits[0]].iStation;
1879
1880 const int dir = (mc.z - algo->vStations[fSta].z[0]) / abs(mc.z - algo->vStations[fSta].z[0]);
1881 // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large
1882 // distance
1883 for (int iSta = fSta /*+dir*/;
1884 (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0);
1885 iSta += dir)
1886 {
1887 // cout << iSta << " " << dir << endl;
1888 L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp);
1889 if (iSta + dir == NMvdStations - 1)
1890 L1AddPipeMaterial(trPar, trPar.qp);
1891 }
1892 }
1893 if (mc.z != trPar.z[0])
1894 continue;
1895 // static int good = 0;
1896 // static int bad = 0;
1897 // if (mc.z != trPar.z[0]){
1898 // bad++;
1899 // continue;
1900 // }
1901 // else good++;
1902 // cout << "bad\\good" << bad << " " << good << endl;
1903
1904 // calculate pulls
1905 // h_fitPV[0]->Fill( (mc.x-trPar.x[0]) *1.e4);
1906 // h_fitPV[1]->Fill( (mc.y-trPar.y[0]) *1.e4);
1907 h_fitPV[0]->Fill((mc.x - trPar.x[0]));
1908 h_fitPV[1]->Fill((mc.y - trPar.y[0]));
1909 h_fitPV[2]->Fill((mc.px / mc.pz - trPar.tx[0]) * 1.e3);
1910 h_fitPV[3]->Fill((mc.py / mc.pz - trPar.ty[0]) * 1.e3);
1911 h_fitPV[4]->Fill(trPar.qp[0] / mc.q * mc.p - 1);
1912 if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1913 h_fitPV[5]->Fill((mc.x - trPar.x[0]) / sqrt(trPar.C00[0]));
1914 if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1915 h_fitPV[6]->Fill((mc.y - trPar.y[0]) / sqrt(trPar.C11[0]));
1916 if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1917 h_fitPV[7]->Fill((mc.px / mc.pz - trPar.tx[0]) / sqrt(trPar.C22[0]));
1918 if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1919 h_fitPV[8]->Fill((mc.py / mc.pz - trPar.ty[0]) / sqrt(trPar.C33[0]));
1920 if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1921 h_fitPV[9]->Fill((mc.q / mc.p - trPar.qp[0]) / sqrt(trPar.C44[0]));
1922 h_fitPV[10]->Fill(trPar.qp[0]);
1923 h_fitPV[11]->Fill(mc.q / mc.p);
1924#else
1925 FairTrackParam fTP;
1926
1927 CbmKFMath::CopyTC2TrackParam(&fTP, it->T, it->C);
1928
1929 CbmKFTrack kfTr;
1930 kfTr.SetTrackParam(fTP);
1931
1932 kfTr.Extrapolate(mc.z);
1933 CbmL1Track it2;
1934 for (int ipar = 0; ipar < 6; ipar++)
1935 it2.T[ipar] = kfTr.GetTrack()[ipar];
1936 for (int ipar = 0; ipar < 15; ipar++)
1937 it2.C[ipar] = kfTr.GetCovMatrix()[ipar];
1938
1939 // calculate pulls
1940 // h_fitPV[0]->Fill( (mc.x-it2.T[0]) *1.e4);
1941 // h_fitPV[1]->Fill( (mc.y-it2.T[1]) *1.e4);
1942 h_fitPV[0]->Fill((mc.x - it2.T[0]));
1943 h_fitPV[1]->Fill((mc.y - it2.T[1]));
1944 h_fitPV[2]->Fill((mc.px / mc.pz - it2.T[2]) * 1.e3);
1945 h_fitPV[3]->Fill((mc.py / mc.pz - it2.T[3]) * 1.e3);
1946 h_fitPV[4]->Fill(it2.T[4] / mc.q * mc.p - 1);
1947 if (finite(it2.C[0]) && it2.C[0] > 0)
1948 h_fitPV[5]->Fill((mc.x - it2.T[0]) / sqrt(it2.C[0]));
1949 if (finite(it2.C[2]) && it2.C[2] > 0)
1950 h_fitPV[6]->Fill((mc.y - it2.T[1]) / sqrt(it2.C[2]));
1951 if (finite(it2.C[5]) && it2.C[5] > 0)
1952 h_fitPV[7]->Fill((mc.px / mc.pz - it2.T[2]) / sqrt(it2.C[5]));
1953 if (finite(it2.C[9]) && it2.C[9] > 0)
1954 h_fitPV[8]->Fill((mc.py / mc.pz - it2.T[3]) / sqrt(it2.C[9]));
1955 if (finite(it2.C[14]) && it2.C[14] > 0)
1956 h_fitPV[9]->Fill((mc.q / mc.p - it2.T[4]) / sqrt(it2.C[14]));
1957 h_fitPV[10]->Fill(it2.T[4]);
1958 h_fitPV[11]->Fill(mc.q / mc.p);
1959#endif // L1EXTRAPOLATE
1960 }
1961 }
1962
1963 h_fit_chi2->Fill(it->chi2 / it->NDF);
1964 }
1965
1966} // void CbmL1::TrackFitPerformance()
1967
1968void CbmL1::FieldApproxCheck()
1969{
1970 TDirectory* curr = gDirectory;
1971 TFile* currentFile = gFile;
1972 TFile* fout = new TFile("FieldApprox.root", "RECREATE");
1973 fout->cd();
1974
1975 FairField* MF = CbmKF::Instance()->GetMagneticField();
1976 for (int ist = 0; ist < NStation; ist++) {
1977 double z = 0;
1978 double Xmax = -100, Ymax = -100;
1979 if (ist < NMvdStations) {
1980 CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
1981 z = t.z;
1982 Xmax = Ymax = t.R;
1983 } else {
1984 // AZ CbmStsStation *st = StsDigi.GetStation(ist - NMvdStations);
1985 CbmStsStation* st = StsDigi->GetStation(ist - NMvdStations);
1986 z = st->GetZ();
1987
1988 double x, y;
1989 for (int isec = 0; isec < st->GetNSectors(); isec++) {
1990 CbmStsSector* sect = L1_DYNAMIC_CAST<CbmStsSector*>(st->GetSector(isec));
1991 for (int isen = 0; isen < sect->GetNSensors(); isen++) {
1992 x = sect->GetSensor(isen)->GetX0() + sect->GetSensor(isen)->GetLx() / 2.;
1993 y = sect->GetSensor(isen)->GetY0() + sect->GetSensor(isen)->GetLy() / 2.;
1994 if (x > Xmax)
1995 Xmax = x;
1996 if (y > Ymax)
1997 Ymax = y;
1998 }
1999 }
2000 cout << "Station " << ist << ", Xmax " << Xmax << ", Ymax" << Ymax << endl;
2001
2002 } // if mvd
2003
2004 // float step = 1.;
2005
2006 int NbinsX = 100; // static_cast<int>(2*Xmax/step);
2007 int NbinsY = 100; // static_cast<int>(2*Ymax/step);
2008 float ddx = 2 * Xmax / NbinsX;
2009 float ddy = 2 * Ymax / NbinsY;
2010
2011 TH2F* stB = new TH2F(Form("station %i, dB", ist + 1), Form("station %i, dB, z = %0.f cm", ist + 1, z),
2012 static_cast<int>(NbinsX + 1), -(Xmax + ddx / 2.), (Xmax + ddx / 2.),
2013 static_cast<int>(NbinsY + 1), -(Ymax + ddy / 2.), (Ymax + ddy / 2.));
2014 TH2F* stBx = new TH2F(Form("station %i, dBx", ist + 1), Form("station %i, dBx, z = %0.f cm", ist + 1, z),
2015 static_cast<int>(NbinsX + 1), -(Xmax + ddx / 2.), (Xmax + ddx / 2.),
2016 static_cast<int>(NbinsY + 1), -(Ymax + ddy / 2.), (Ymax + ddy / 2.));
2017 TH2F* stBy = new TH2F(Form("station %i, dBy", ist + 1), Form("station %i, dBy, z = %0.f cm", ist + 1, z),
2018 static_cast<int>(NbinsX + 1), -(Xmax + ddx / 2.), (Xmax + ddx / 2.),
2019 static_cast<int>(NbinsY + 1), -(Ymax + ddy / 2.), (Ymax + ddy / 2.));
2020 TH2F* stBz = new TH2F(Form("station %i, dBz", ist + 1), Form("station %i, dBz, z = %0.f cm", ist + 1, z),
2021 static_cast<int>(NbinsX + 1), -(Xmax + ddx / 2.), (Xmax + ddx / 2.),
2022 static_cast<int>(NbinsY + 1), -(Ymax + ddy / 2.), (Ymax + ddy / 2.));
2023
2024 Double_t r[3], B[3];
2025 L1FieldSlice FSl;
2026 L1FieldValue B_L1;
2027 Double_t bbb, bbb_L1;
2028
2029 const int M = 5; // polinom order
2030 const int N = (M + 1) * (M + 2) / 2;
2031 L1Station& st = algo->vStations[ist];
2032 for (int i = 0; i < N; i++) {
2033 FSl.cx[i] = st.fieldSlice.cx[i][0];
2034 FSl.cy[i] = st.fieldSlice.cy[i][0];
2035 FSl.cz[i] = st.fieldSlice.cz[i][0];
2036 }
2037
2038 Int_t i = 1, j = 1;
2039
2040 double x, y;
2041 for (int ii = 1; ii <= NbinsX + 1; ii++) {
2042 j = 1;
2043 x = -Xmax + (ii - 1) * ddx;
2044 for (int jj = 1; jj <= NbinsY + 1; jj++) {
2045 y = -Ymax + (jj - 1) * ddy;
2046 double rrr = sqrt(fabs(x * x / Xmax / Xmax + y / Ymax * y / Ymax));
2047 if (rrr > 1.) {
2048 j++;
2049 continue;
2050 }
2051 r[2] = z;
2052 r[0] = x;
2053 r[1] = y;
2054 MF->GetFieldValue(r, B);
2055 bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
2056
2057 FSl.GetFieldValue(x, y, B_L1);
2058 bbb_L1 = sqrt(B_L1.x[0] * B_L1.x[0] + B_L1.y[0] * B_L1.y[0] + B_L1.z[0] * B_L1.z[0]);
2059
2060 stB->SetBinContent(ii, jj, (bbb - bbb_L1));
2061 stBx->SetBinContent(ii, jj, (B[0] - B_L1.x[0]));
2062 stBy->SetBinContent(ii, jj, (B[1] - B_L1.y[0]));
2063 stBz->SetBinContent(ii, jj, (B[2] - B_L1.z[0]));
2064 j++;
2065 }
2066 i++;
2067 }
2068
2069 stB->GetXaxis()->SetTitle("X, cm");
2070 stB->GetYaxis()->SetTitle("Y, cm");
2071 stB->GetXaxis()->SetTitleOffset(1);
2072 stB->GetYaxis()->SetTitleOffset(1);
2073 stB->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
2074 stB->GetZaxis()->SetTitleOffset(1.3);
2075
2076 stBx->GetXaxis()->SetTitle("X, cm");
2077 stBx->GetYaxis()->SetTitle("Y, cm");
2078 stBx->GetXaxis()->SetTitleOffset(1);
2079 stBx->GetYaxis()->SetTitleOffset(1);
2080 stBx->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
2081 stBx->GetZaxis()->SetTitleOffset(1.3);
2082
2083 stBy->GetXaxis()->SetTitle("X, cm");
2084 stBy->GetYaxis()->SetTitle("Y, cm");
2085 stBy->GetXaxis()->SetTitleOffset(1);
2086 stBy->GetYaxis()->SetTitleOffset(1);
2087 stBy->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
2088 stBy->GetZaxis()->SetTitleOffset(1.3);
2089
2090 stBz->GetXaxis()->SetTitle("X, cm");
2091 stBz->GetYaxis()->SetTitle("Y, cm");
2092 stBz->GetXaxis()->SetTitleOffset(1);
2093 stBz->GetYaxis()->SetTitleOffset(1);
2094 stBz->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
2095 stBz->GetZaxis()->SetTitleOffset(1.3);
2096
2097 stB->Write();
2098 stBx->Write();
2099 stBy->Write();
2100 stBz->Write();
2101
2102 } // for ista
2103
2104 fout->Close();
2105 fout->Delete();
2106 gFile = currentFile;
2107 gDirectory = curr;
2108} // void CbmL1::FieldApproxCheck()
2109
2110#include "TMath.h"
2111void CbmL1::FieldIntegralCheck()
2112{
2113 TDirectory* curr = gDirectory;
2114 TFile* currentFile = gFile;
2115 TFile* fout = new TFile("FieldApprox.root", "RECREATE");
2116 fout->cd();
2117
2118 FairField* MF = CbmKF::Instance()->GetMagneticField();
2119
2120 int nPointsZ = 1000;
2121 int nPointsPhi = 100;
2122 int nPointsTheta = 100;
2123 double startZ = 0, endZ = 100.;
2124 double startPhi = 0, endPhi = 2 * TMath::Pi();
2125 double startTheta = -30. / 180. * TMath::Pi(), endTheta = 30. / 180. * TMath::Pi();
2126
2127 double DZ = endZ - startZ;
2128 double DP = endPhi - startPhi;
2129 double DT = endTheta - startTheta;
2130
2131 float ddp = endPhi / nPointsPhi;
2132 float ddt = 2 * endTheta / nPointsTheta;
2133
2134 TH2F* hSb =
2135 new TH2F("Field Integral", "Field Integral", static_cast<int>(nPointsPhi), -(startPhi + ddp / 2.),
2136 (endPhi + ddp / 2.), static_cast<int>(nPointsTheta), (startTheta - ddt / 2.), (endTheta + ddt / 2.));
2137
2138 for (int iP = 0; iP < nPointsPhi; iP++) {
2139 double phi = startPhi + iP * DP / nPointsPhi;
2140 for (int iT = 0; iT < nPointsTheta; iT++) {
2141 double theta = startTheta + iT * DT / nPointsTheta;
2142
2143 double Sb = 0;
2144 for (int iZ = 1; iZ < nPointsZ; iZ++) {
2145 double z = startZ + DZ * iZ / nPointsZ;
2146 double x = z * TMath::Tan(theta) * TMath::Cos(phi);
2147 double y = z * TMath::Tan(theta) * TMath::Sin(phi);
2148 double r[3] = {x, y, z};
2149 double b[3];
2150 MF->GetFieldValue(r, b);
2151 double B = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
2152 Sb += B * DZ / nPointsZ / 100. / 10.;
2153 }
2154 hSb->SetBinContent(iP + 1, iT + 1, Sb);
2155 }
2156 }
2157
2158 hSb->GetXaxis()->SetTitle("#phi [rad]");
2159 hSb->GetYaxis()->SetTitle("#theta [rad]");
2160 hSb->GetXaxis()->SetTitleOffset(1);
2161 hSb->GetYaxis()->SetTitleOffset(1);
2162 hSb->GetZaxis()->SetTitle("Field Integral [T#dotm]");
2163 hSb->GetZaxis()->SetTitleOffset(1.3);
2164
2165 hSb->Write();
2166
2167 fout->Close();
2168 fout->Delete();
2169 gFile = currentFile;
2170 gDirectory = curr;
2171} // void CbmL1::FieldApproxCheck()
2172
2173void CbmL1::InputPerformance()
2174{
2175 static TH1I *nStripFHits, *nStripBHits, *nStripFMC, *nStripBMC;
2176
2177 static TH1F *resX, *resY /*, *pullZ*/;
2178 static TH1F *pullX, *pullY /*, *pullZ*/;
2179
2180 static bool first_call = 1;
2181 if (first_call) {
2182 first_call = 0;
2183
2184 TDirectory* currentDir = gDirectory;
2185 gDirectory = histodir;
2186 gDirectory->cd("L1");
2187 gDirectory->mkdir("Input");
2188 gDirectory->cd("Input");
2189
2190 nStripFHits = new TH1I("nHits_f", "NHits On Front Strip", 10, 0, 10);
2191 nStripBHits = new TH1I("nHits_b", "NHits On Back Strip", 10, 0, 10);
2192 nStripFMC = new TH1I("nMC_f", "N MC Points On Front Strip", 10, 0, 10);
2193 nStripBMC = new TH1I("nMC_b", "N MC Points On Back Strip", 10, 0, 10);
2194
2195 pullX = new TH1F("Px", "Pull x", 50, -5, 5);
2196 pullY = new TH1F("Py", "Pull y", 50, -5, 5);
2197
2198 resX = new TH1F("x", "Residual x", 50, -50, 50);
2199 resY = new TH1F("y", "Residual y", 50, -500, 500);
2200 TH1* histo;
2201 histo = resX;
2202 histo->GetXaxis()->SetTitle("Residual x, um");
2203 histo = resY;
2204 histo->GetXaxis()->SetTitle("Residual y, um");
2205
2206 gDirectory = currentDir;
2207 } // first_call
2208
2209 std::map<unsigned int, unsigned int> stripFToNHitMap, stripBToNHitMap;
2210 std::map<unsigned int, unsigned int> stripFToNMCMap, stripBToNMCMap;
2211
2212 map<unsigned int, unsigned int>::iterator it;
2213 Int_t nMC = -1;
2214 if (listStsPts) {
2215 nMC = listStsPts->GetEntries();
2216 }
2217 if (listStsHits) {
2218 for (unsigned int iH = 0; iH < vStsHits.size(); iH++) {
2219 const CbmL1StsHit& h = vStsHits[iH];
2220
2221 if (h.extIndex < 0)
2222 continue; // mvd hit
2223 const CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(h.extIndex));
2224 // strip - MC correspondence
2225 int iStripF = sh->GetDigi(0);
2226 int iStripB = sh->GetDigi(1);
2227 if (iStripF >= 0)
2228 stripFToNHitMap[iStripF]++;
2229 if (iStripB >= 0)
2230 stripBToNHitMap[iStripB]++;
2231
2232 if (h.mcPointIds.size() == 0)
2233 continue; // fake
2234 int iMC = vMCPoints[h.mcPointIds[0]].pointId; // TODO: adapt to linking
2235 if (!(iMC >= 0 && iMC < nMC))
2236 continue;
2237
2238 if (iStripF >= 0)
2239 stripFToNMCMap[iStripF]++;
2240 if (iStripB >= 0)
2241 stripBToNMCMap[iStripB]++;
2242
2243 // hit pulls and residuals
2244
2245 TVector3 hitPos, mcPos, hitErr;
2246 sh->Position(hitPos);
2247 sh->PositionError(hitErr);
2248 CbmStsPoint* pt = 0;
2249 pt = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(iMC));
2250 if (!pt) {
2251 // cout << " No MC points! " << "iMC=" << iMC << endl;
2252 continue;
2253 }
2254 // pt->Position(mcPos); // this is wrong!
2255 mcPos.SetX(pt->GetX(hitPos.Z()));
2256 mcPos.SetY(pt->GetY(hitPos.Z()));
2257 mcPos.SetZ(hitPos.Z());
2258
2259#if 0 // standard errors
2260 if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2261 if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2262#elif 1 // qa errors
2263 if (hitErr.X() != 0)
2264 pullX->Fill((hitPos.X() - mcPos.X()) / sh->GetDx()); // qa errors
2265 if (hitErr.Y() != 0)
2266 pullY->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
2267#else // errors used in TF
2268 if (hitErr.X() != 0)
2269 pullX->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
2270 if (hitErr.Y() != 0)
2271 pullY->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
2272#endif
2273
2274 resX->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
2275 resY->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
2276 }
2277 } // sts
2278 if (listMvdHits) {
2279 Int_t nEnt = listMvdHits->GetEntries();
2280 // Int_t nMC = listStsPts->GetEntries();
2281 for (int j = 0; j < nEnt; j++) {
2282 CbmMvdHit* sh = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(j));
2283 CbmMvdHitMatch* hm = L1_DYNAMIC_CAST<CbmMvdHitMatch*>(listMvdHitMatches->At(j));
2284
2285 // int iMC = sh->GetRefIndex();
2286 int iMC = hm->GetPointId();
2287
2288 // hit pulls and residuals
2289
2290 TVector3 hitPos, mcPos, hitErr;
2291 sh->Position(hitPos);
2292 sh->PositionError(hitErr);
2293 CbmMvdPoint* pt = 0;
2294 if (iMC >= 0 && iMC < nMC)
2295 pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(listMvdPts->At(iMC));
2296 if (!pt) {
2297 // cout << " No MC points! " << "iMC=" << iMC << endl;
2298 continue;
2299 }
2300 // pt->Position(mcPos); // this is wrong!
2301 mcPos.SetX((pt->GetX() + pt->GetXOut()) / 2.);
2302 mcPos.SetY((pt->GetY() + pt->GetYOut()) / 2);
2303 mcPos.SetZ(hitPos.Z());
2304
2305 // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2306 // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2307 // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
2308 // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
2309 if (hitErr.X() != 0)
2310 pullX->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[0].XYInfo.C00[0])); // errors used in TF
2311 if (hitErr.Y() != 0)
2312 pullY->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[0].XYInfo.C11[0]));
2313
2314 resX->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
2315 resY->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
2316 }
2317 } // mvd
2318
2319 for (it = stripFToNHitMap.begin(); it != stripFToNHitMap.end(); it++) {
2320 nStripFHits->Fill(it->second);
2321 }
2322 for (it = stripBToNHitMap.begin(); it != stripBToNHitMap.end(); it++) {
2323 nStripBHits->Fill(it->second);
2324 }
2325 for (it = stripFToNMCMap.begin(); it != stripFToNMCMap.end(); it++) {
2326 nStripFMC->Fill(it->second);
2327 }
2328 for (it = stripBToNMCMap.begin(); it != stripBToNMCMap.end(); it++) {
2329 nStripBMC->Fill(it->second);
2330 }
2331
2332 // strips Not ended
2333 // if( listCbmStsDigi ){
2334 // Int_t nEnt = listCbmStsDigi->GetEntries();
2335 // for (int j=0; j < nEnt; j++ ){
2336 // CbmStsDigi *digi = (CbmStsDigi*) listCbmStsDigi->At(j);
2337 // // = sh->GetNLinks(0);
2338 // // find position of mc point
2339 // FairLink link = digi->GetLink(0);
2340 // int iMCPoint = link.GetIndex();
2341 // CbmStsPoint *mcPoint = (CbmStsPoint*) listStsPts->At(iMCPoint);
2342 // TVector3 mcPos;
2343 // if (digi->GetSide() == 0)
2344 // mcPoint->PositionIn(mcPos);
2345 // else
2346 // mcPoint->PositionOut(mcPos);
2347 //
2348 // CbmStsStation *sta = StsDigi.GetStation(digi->GetStationNr());
2349 // CbmStsSector* sector = sta->GetSector(digi->GetSectorNr());
2350 // digi->GetChannelNr();
2351 //
2352 // }
2353 // } // listCbmStsDigi
2354
2355} // void CbmL1::InputPerformance()
L1_DYNAMIC_CAST
#define L1_DYNAMIC_CAST
Definition CbmL1Def.h:26
L1_ASSERT
#define L1_ASSERT(v, msg)
Definition CbmL1Def.h:42
L1_assert
#define L1_assert(v)
Definition CbmL1Def.h:49
L1AddPipeMaterial
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:91
L1AddMaterial
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:11
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, L1FieldRegion &F, fvec *w=0)
Definition L1Extrapolation.h:11
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
_fvecalignment
#define _fvecalignment
Definition P4_F32vec4.h:236
i
int i
Definition P4_F32vec4.h:22
Int_t
int Int_t
Double_t
double Double_t
gDirectory
#define gDirectory
gFile
#define gFile
Form
char * Form(const char *fmt,...)
iZ
@ iZ
Definition BmnAlignDefines.h:87
CbmKFMath::CopyTC2TrackParam
static void CopyTC2TrackParam(FairTrackParam *par, Double_t T[], Double_t C[])
Definition CbmKFMath.cxx:773
CbmKFPartEfficiencies::Inc
void Inc(bool isReco, int nClones, TString name)
Definition CbmKFPartEfficiencies.h:444
CbmKFPartEfficiencies::IncReco
void IncReco(bool isGhost, bool isBg, TString name)
Definition CbmKFPartEfficiencies.h:455
CbmKFPartEfficiencies::partName
TString partName[nParticles]
Definition CbmKFPartEfficiencies.h:531
CbmKFParticle::NDaughters
int NDaughters() const
Definition CbmKFParticle.h:24
CbmKFParticle::Extrapolate
void Extrapolate(Double_t r0[], double T)
Definition CbmKFParticle.cxx:1191
CbmKFParticle::GetMass
void GetMass(Double_t &M, Double_t &Error)
Definition CbmKFParticle.cxx:1446
CbmKFParticle::r
Double_t r[8]
Definition CbmKFParticle.h:40
CbmKFParticle::GetNDF
Int_t GetNDF() const
Definition CbmKFParticle.h:91
CbmKFParticle::GetLifeTime
void GetLifeTime(Double_t &T, Double_t &Error)
Definition CbmKFParticle.cxx:1476
CbmKFParticle::GetDecayLength
void GetDecayLength(Double_t &L, Double_t &Error)
Definition CbmKFParticle.cxx:1459
CbmKFParticle::GetParameter
Double_t GetParameter(Int_t i) const
Definition CbmKFParticle.h:97
CbmKFParticle::GetChi2
Double_t GetChi2() const
Definition CbmKFParticle.h:90
CbmKFParticle::GetPDG
int GetPDG() const
Definition CbmKFParticle.h:139
CbmKFParticle::GetCovariance
Double_t GetCovariance(Int_t i) const
Definition CbmKFParticle.h:98
CbmKFParticle::Id
int Id() const
Definition CbmKFParticle.h:23
CbmKFParticle::DaughterIds
const vector< int > & DaughterIds() const
Definition CbmKFParticle.h:25
CbmKFParticle::GetDStoPoint
Double_t GetDStoPoint(const Double_t xyz[]) const
Definition CbmKFParticle.cxx:1502
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0, Bool_t line=false)
Access to i-th hit.
Definition CbmKFTrackInterface.cxx:48
CbmKFTrack::SetTrackParam
void SetTrackParam(FairTrackParam &track)
Definition CbmKFTrack.cxx:35
CbmKFTrack::GetTrack
Double_t * GetTrack()
Is it electron.
Definition CbmKFTrack.h:74
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition CbmKFTrack.h:75
CbmKFTube::R
Double_t R
Definition CbmKFMaterial.h:90
CbmKFTube::z
Double_t z
Definition CbmKFMaterial.h:89
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:34
CbmKF::GetMagneticField
BmnNewFieldMap * GetMagneticField()
Definition CbmKF.h:72
CbmKF::vMvdMaterial
std::vector< CbmKFTube > vMvdMaterial
Definition CbmKF.h:59
CbmL1HitStore::y
double y
Definition CbmL1.h:44
CbmL1HitStore::x
double x
Definition CbmL1.h:44
CbmL1HitStore::iStation
int iStation
Definition CbmL1.h:43
CbmL1MCTrack::GetNClones
int GetNClones() const
Definition CbmL1MCTrack.h:105
CbmL1MCTrack::AddRecoTrack
void AddRecoTrack(CbmL1Track *rTr)
Definition CbmL1MCTrack.h:103
CbmL1MCTrack::IsReconstructed
bool IsReconstructed() const
Definition CbmL1MCTrack.h:106
CbmL1MCTrack::IsReconstructable
bool IsReconstructable() const
Definition CbmL1MCTrack.h:95
CbmL1MCTrack::NStations
int NStations() const
Definition CbmL1MCTrack.h:97
CbmL1MCTrack::StsHits
vector< int > StsHits
Definition CbmL1MCTrack.h:36
CbmL1MCTrack::AddTouchTrack
void AddTouchTrack(CbmL1Track *tTr)
Definition CbmL1MCTrack.h:108
CbmL1MCTrack::GetRecoTracks
vector< CbmL1Track * > & GetRecoTracks()
Definition CbmL1MCTrack.h:104
CbmL1MCTrack::IsAdditional
bool IsAdditional() const
Definition CbmL1MCTrack.h:96
CbmL1MCTrack::IsDisturbed
bool IsDisturbed() const
Definition CbmL1MCTrack.h:109
CbmL1MCTrack::NMCStations
int NMCStations() const
Definition CbmL1MCTrack.h:98
CbmL1MCTrack::IsPrimary
bool IsPrimary() const
Definition CbmL1MCTrack.h:94
CbmL1MCTrack::Points
vector< int > Points
Definition CbmL1MCTrack.h:35
CbmL1StsHit::mcPointIds
vector< int > mcPointIds
Definition CbmL1StsHit.h:22
CbmL1Track::GetMaxPurity
double GetMaxPurity()
Definition CbmL1Track.h:52
CbmL1Track::GetMCTrack
CbmL1MCTrack * GetMCTrack()
Definition CbmL1Track.h:47
CbmL1Track::hitMap
map< int, int > hitMap
Definition CbmL1Track.h:61
CbmL1Track::AddMCTrack
void AddMCTrack(CbmL1MCTrack *mcTr)
Definition CbmL1Track.h:45
CbmL1Track::SetMaxPurity
void SetMaxPurity(double maxPurity_)
Definition CbmL1Track.h:51
CbmL1Track::StsHits
vector< int > StsHits
Definition CbmL1Track.h:58
CbmL1Track::IsGhost
bool IsGhost()
Definition CbmL1Track.h:49
CbmL1::algo
L1Algo * algo
Definition CbmL1.h:54
CbmL1::vHitStore
vector< CbmL1HitStore > vHitStore
Definition CbmL1.h:95
CbmL1::NMvdStations
int NMvdStations
Definition CbmL1.h:143
CbmL1::vMCParticles
vector< KFMCParticle > vMCParticles
reco particles
Definition CbmL1.h:179
CbmL1::FindReconstructableMCParticles
void FindReconstructableMCParticles()
Definition CbmL1Performance.cxx:478
CbmL1::CheckMCParticleIsReconstructable
void CheckMCParticleIsReconstructable(KFMCParticle &part)
Definition CbmL1Performance.cxx:488
CbmL1::TrackMatch
void TrackMatch()
Reconstruction Performance.
Definition CbmL1Performance.cxx:56
CbmL1::listMCTracks
TClonesArray * listMCTracks
Input data.
Definition CbmL1.h:161
CbmL1::vMCTracks
vector< CbmL1MCTrack > vMCTracks
Definition CbmL1.h:174
CbmL1::vMCPoints
vector< CbmL1MCPoint > vMCPoints
hits with hit-mcpoint match information
Definition CbmL1.h:173
CbmL1::listMvdHitMatches
TClonesArray * listMvdHitMatches
Definition CbmL1.h:169
CbmL1::vStsHits
vector< CbmL1StsHit > vStsHits
Used data = Repacked input data.
Definition CbmL1.h:172
CbmL1::PartHistoPerformance
void PartHistoPerformance()
Definition CbmL1Performance.cxx:730
CbmL1::MCtoRParticleId
vector< KFMatchParticles > MCtoRParticleId
MC particles.
Definition CbmL1.h:180
CbmL1::listMvdPts
TClonesArray * listMvdPts
Definition CbmL1.h:167
CbmL1::vRTracks
vector< CbmL1Track > vRTracks
Definition CbmL1.h:57
CbmL1::InputPerformance
void InputPerformance()
Definition CbmL1Performance.cxx:2173
CbmL1::vRParticles
vector< CbmKFParticle > vRParticles
Definition CbmL1.h:178
CbmL1::FieldIntegralCheck
void FieldIntegralCheck()
Definition CbmL1Performance.cxx:2111
CbmL1::listStsPts
TClonesArray * listStsPts
Definition CbmL1.h:162
CbmL1::fFindParticlesMode
int fFindParticlesMode
Definition CbmL1.h:191
CbmL1::listMvdHits
TClonesArray * listMvdHits
Mvd MC points.
Definition CbmL1.h:168
CbmL1::vHitMCRef
vector< int > vHitMCRef
Definition CbmL1.h:175
CbmL1::TrackFitPerformance
void TrackFitPerformance()
Definition CbmL1Performance.cxx:1642
CbmL1::histodir
TDirectory * histodir
Definition CbmL1.h:183
CbmL1::PF
CbmL1ParticlesFinder * PF
for access to L1 Algorithm from L1::Instance
Definition CbmL1.h:55
CbmL1::MatchParticles
void MatchParticles()
Procedure for match Reconstructed and MC Particles. Should be called before Performances.
Definition CbmL1Performance.cxx:590
CbmL1::listStsHits
TClonesArray * listStsHits
Definition CbmL1.h:165
CbmL1::EfficienciesPerformance
void EfficienciesPerformance()
Definition CbmL1Performance.cxx:269
CbmL1::CbmL1MCTrack
friend class CbmL1MCTrack
Definition CbmL1.h:103
CbmL1::RtoMCParticleId
vector< KFMatchParticles > RtoMCParticleId
array for match
Definition CbmL1.h:181
CbmL1::HistoPerformance
void HistoPerformance()
Definition CbmL1Performance.cxx:1074
CbmL1::FieldApproxCheck
void FieldApproxCheck()
Definition CbmL1Performance.cxx:1968
CbmL1::PartEffPerformance
void PartEffPerformance()
Definition CbmL1Performance.cxx:664
CbmL1::NStation
int NStation
Definition CbmL1.h:143
CbmL1::GetMCParticles
void GetMCParticles()
Definition CbmL1Performance.cxx:440
CbmL1::StsDigi
CbmStsDigiScheme * StsDigi
Definition CbmL1.h:156
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetMass
Double_t GetMass() const
Definition CbmMCTrack.cxx:127
CbmMCTrack::GetPy
Double_t GetPy() const
Definition CbmMCTrack.h:59
CbmMCTrack::GetPz
Double_t GetPz() const
Definition CbmMCTrack.h:60
CbmMCTrack::GetPx
Double_t GetPx() const
Definition CbmMCTrack.h:58
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition CbmMCTrack.h:63
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition CbmMCTrack.h:62
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition CbmMCTrack.h:61
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Accessors.
Definition CbmMCTrack.h:56
CbmMvdHitMatch
CbmMvdHitMatch header file.
Definition CbmMvdHitMatch.h:20
CbmMvdHitMatch::GetPointId
Int_t GetPointId() const
Definition CbmMvdHitMatch.h:36
CbmMvdHit
CbmMvdHit.h.
Definition CbmMvdHit.h:26
CbmMvdPoint
CbmMvdPoint.h.
Definition CbmMvdPoint.h:27
CbmMvdPoint::GetYOut
Double_t GetYOut() const
Definition CbmMvdPoint.h:64
CbmMvdPoint::GetXOut
Double_t GetXOut() const
Accessors.
Definition CbmMvdPoint.h:63
CbmStsHit::GetDigi
Int_t GetDigi(Int_t iSide) const
Definition CbmStsHit.cxx:120
CbmStsPoint::GetY
Double_t GetY(Double_t z) const
Definition CbmStsPoint.cxx:103
CbmStsPoint::GetX
Double_t GetX(Double_t z) const
Point coordinates at given z from linear extrapolation.
Definition CbmStsPoint.cxx:92
CbmStsSector::GetSensor
CbmStsSensor * GetSensor(Int_t iSensor)
Definition CbmStsSector.h:74
CbmStsSector::GetNSensors
Int_t GetNSensors() const
Definition CbmStsSector.h:71
CbmStsSensor::GetX0
Double_t GetX0() const
Definition CbmStsSensor.h:108
CbmStsSensor::GetLy
Double_t GetLy() const
Definition CbmStsSensor.h:113
CbmStsSensor::GetY0
Double_t GetY0() const
Definition CbmStsSensor.h:109
CbmStsSensor::GetLx
Double_t GetLx() const
Definition CbmStsSensor.h:112
CbmStsStation::GetNSectors
Int_t GetNSectors() const
Definition CbmStsStation.h:62
CbmStsStation::GetZ
Double_t GetZ(Int_t it=0)
Definition CbmStsStation.cxx:151
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
KFMCParticle::IsReconstructable
bool IsReconstructable() const
Definition KFMCParticle.h:42
KFMCParticle::GetPDG
int GetPDG() const
Definition KFMCParticle.h:41
KFMCParticle::SetMCTrackID
void SetMCTrackID(int id)
Definition KFMCParticle.h:36
KFMCParticle::GetDaughterIds
const std::vector< int > & GetDaughterIds() const
Definition KFMCParticle.h:32
KFMCParticle::NDaughters
int NDaughters() const
Definition KFMCParticle.h:31
KFMCParticle::GetMotherId
int GetMotherId() const
Definition KFMCParticle.h:40
KFMCParticle::SetAsReconstructable
void SetAsReconstructable()
Definition KFMCParticle.h:43
KFMCParticle::SetPDG
void SetPDG(int pdg)
Definition KFMCParticle.h:35
KFMCParticle::SetMotherId
void SetMotherId(int id)
Definition KFMCParticle.h:37
KFMCParticle::GetMCTrackID
int GetMCTrackID() const
Definition KFMCParticle.h:39
KFMCParticle::AddDaughter
void AddDaughter(int i)
Definition KFMCParticle.cxx:29
L1FieldSlice::GetFieldValue
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition L1Field.h:47
L1FieldSlice::cy
fvec cy[21]
Definition L1Field.h:39
L1FieldSlice::cx
fvec cx[21]
Definition L1Field.h:39
L1FieldSlice::cz
fvec cz[21]
Definition L1Field.h:39
L1Station::fieldSlice
L1FieldSlice fieldSlice
Definition L1Station.h:29
TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
TDirectoryFile::cd
Bool_t cd() override
TFile::Delete
void Delete(const char *namecycle="") override
TFile::Close
void Close(Option_t *option="") override
TH1::GetZaxis
TAxis * GetZaxis()
TH1::GetXaxis
TAxis * GetXaxis()
TH1::GetYaxis
TAxis * GetYaxis()
TH1::Fill
virtual Int_t Fill(const char *name, Double_t w)
TH1::Scale
virtual void Scale(Double_t c1=1, Option_t *option="")
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TH2::SetBinContent
void SetBinContent(Int_t bin, Double_t content) override
TNamed::SetTitle
virtual void SetTitle(const char *title="")
TObject::Write
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
TProfile::Fill
Int_t Fill(const char *namex, Double_t y) override
TString::Data
const char * Data() const
TVector3::Z
Double_t Z() const
TVector3::SetY
void SetY(Double_t)
TVector3::Y
Double_t Y() const
TVector3::X
Double_t X() const
TVector3::SetX
void SetX(Double_t)
TVector3::SetZ
void SetZ(Double_t)
abs
RVec< PromoteType< T > > abs(const RVec< T > &v)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
acos
RVec< PromoteType< T > > acos(const RVec< T > &v)
n
*brief user defined string literal for JSON values *sa std::size_t n
Definition json.hpp:20362
CbmL1Constants::MinPurity
const double MinPurity
Definition CbmL1Constants.h:9
CbmL1Constants::MinRefMom
const double MinRefMom
Definition CbmL1Constants.h:8
TMath::Prob
Double_t Prob(Double_t chi2, Int_t ndf)
TMath::Sqrt
Double_t Sqrt(Double_t x)
TMath::Cos
Double_t Cos(Double_t)
TMath::Pi
constexpr Double_t Pi()
TMath::Sin
Double_t Sin(Double_t)
TMath::Tan
Double_t Tan(Double_t)
CbmL1TrackPar::C
double C[15]
Definition CbmL1TrackPar.h:22
TL1Efficiencies::reco
TL1TracksCatCounters< int > reco
Definition CbmL1Counters.h:165
TL1Efficiencies::mc
TL1TracksCatCounters< int > mc
Definition CbmL1Counters.h:164
TL1Efficiencies::names
vector< TString > names
Definition CbmL1Counters.h:157
TL1Efficiencies::indices
map< TString, int > indices
Definition CbmL1Counters.h:158
TL1Efficiencies::operator+=
TL1Efficiencies & operator+=(TL1Efficiencies &a)
Definition CbmL1Counters.h:194
TL1Efficiencies::ratio_reco
TL1TracksCatCounters< double > ratio_reco
Definition CbmL1Counters.h:160
TL1Efficiencies::Inc
void Inc(bool isReco, TString name)
Definition CbmL1Counters.h:205
TL1Efficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition CbmL1Counters.h:171
TL1PerfEfficiencies::ratio_clone
TL1TracksCatCounters< double > ratio_clone
Definition CbmL1Performance.cxx:257
TL1PerfEfficiencies::mc_length_hits
TL1TracksCatCounters< int > mc_length_hits
Definition CbmL1Performance.cxx:266
TL1PerfEfficiencies::ratio_length
TL1TracksCatCounters< double > ratio_length
Definition CbmL1Performance.cxx:258
TL1PerfEfficiencies::ratio_fakes
TL1TracksCatCounters< double > ratio_fakes
Definition CbmL1Performance.cxx:259
TL1PerfEfficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition CbmL1Performance.cxx:169
TL1PerfEfficiencies::reco_length
TL1TracksCatCounters< double > reco_length
Definition CbmL1Performance.cxx:263
TL1PerfEfficiencies::mc_length
TL1TracksCatCounters< int > mc_length
Definition CbmL1Performance.cxx:265
TL1PerfEfficiencies::reco_fakes
TL1TracksCatCounters< double > reco_fakes
Definition CbmL1Performance.cxx:264
TL1PerfEfficiencies::Inc
void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length, int _mc_length_hits, TString name)
Definition CbmL1Performance.cxx:206
TL1PerfEfficiencies::killed
TL1TracksCatCounters< int > killed
Definition CbmL1Performance.cxx:261
TL1PerfEfficiencies::ratio_killed
TL1TracksCatCounters< double > ratio_killed
Definition CbmL1Performance.cxx:256
TL1PerfEfficiencies::clone
TL1TracksCatCounters< int > clone
Definition CbmL1Performance.cxx:262
TL1PerfEfficiencies::operator+=
TL1PerfEfficiencies & operator+=(TL1PerfEfficiencies &a)
Definition CbmL1Performance.cxx:184
TL1TracksCatCounters
counters used for efficiency calculation
Definition CbmL1Counters.h:19