CbmRoot
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 "CbmL1.h"
17 
18 #include "CbmKFTrack.h" // for vertex pulls
19 #include "CbmL1Constants.h"
20 #include "FairTrackParam.h" // for vertex pulls
21 #include "L1Algo/L1AddMaterial.h" // for vertex pulls
22 #include "L1Algo/L1Algo.h"
23 #include "L1Algo/L1Extrapolation.h" // for vertex pulls
24 #include "L1Algo/L1StsHit.h"
25 
26 #include "CbmMuchPixelHit.h"
27 #include "CbmMuchPoint.h"
28 
29 #include "CbmTrdHit.h"
30 #include "CbmTrdPoint.h"
31 
32 #include "CbmTofHit.h"
33 #include "CbmTofPoint.h"
34 
35 #include "CbmKF.h"
36 #include "CbmKFMath.h"
37 
38 #include "CbmStsSetup.h"
39 #include "CbmStsStation.h"
40 
41 #include "CbmMatch.h"
42 
43 #include "CbmL1Counters.h"
44 
45 #include <TFile.h>
46 
47 #include <iostream>
48 #include <list>
49 #include <map>
50 #include <vector>
51 
52 using std::cout;
53 using std::endl;
54 using std::ios;
55 using std::map;
56 using std::pair;
57 using std::setw;
58 using std::vector;
59 
60 void CbmL1::TrackMatch() {
61  map<int, CbmL1MCTrack*> pMCTrackMap;
62  pMCTrackMap.clear();
63 
64  // fill pMCTrackMap
65  for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin();
66  i != vMCTracks.end();
67  ++i) {
68  CbmL1MCTrack& MC = *i;
69 
70  if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()) {
71  pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC));
72  } else {
73  cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
74  }
75  }
76  // -- prepare information about reconstructed tracks
77  const int nRTracks = vRTracks.size();
78  for (int iR = 0; iR < nRTracks; iR++) {
79  CbmL1Track* prtra = &(vRTracks[iR]);
80 
81  // cout<<iR<<" iR"<<endl;
82 
83  int hitsum = prtra->StsHits.size(); // number of hits in track
84 
85  // count how many hits from each mcTrack belong to current recoTrack
86  map<int, int>& hitmap =
87  prtra
88  ->hitMap; // how many hits from each mcTrack belong to current recoTrack
89  for (vector<int>::iterator ih = (prtra->StsHits).begin();
90  ih != (prtra->StsHits).end();
91  ++ih) {
92 
93  const int nMCPoints = vStsHits[*ih].mcPointIds.size();
94  for (int iP = 0; iP < nMCPoints; iP++) {
95  int iMC = vStsHits[*ih].mcPointIds[iP];
96 
97  // cout<<iMC<<" iMC"<<endl;
98  int ID = -1;
99  if (iMC >= 0) ID = vMCPoints[iMC].ID;
100  if (hitmap.find(ID) == hitmap.end())
101  hitmap[ID] = 1;
102  else {
103  hitmap[ID] += 1;
104  }
105  } // for iPoint
106  } // for iHit
107 
108  // RTrack <-> MCTrack identification
109  double max_percent =
110  0.0; // [%]. maximum persent of hits, which belong to one mcTrack.
111  for (map<int, int>::iterator posIt = hitmap.begin(); posIt != hitmap.end();
112  posIt++) { // loop over all touched MCTracks
113 
114  if (posIt->first < 0) continue; // not a MC track - based on fake hits
115 
116  // count max-purity
117  if (double(posIt->second) > max_percent * double(hitsum))
118  max_percent = double(posIt->second) / double(hitsum);
119 
120  // set relation to the mcTrack
121  if (double(posIt->second)
122  >= CbmL1Constants::MinPurity
123  * double(hitsum)) { // found correspondent MCTrack
124  if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
125  CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
126 
127  pmtra->AddRecoTrack(prtra);
128  prtra->AddMCTrack(pmtra);
129  } else {
130  if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
131  CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
132 
133  pmtra->AddTouchTrack(prtra);
134  }
135  } // for hitmap
136  prtra->SetMaxPurity(max_percent);
137  } // for reco tracks
138 } //
139 
140 
141 struct TL1PerfEfficiencies : public TL1Efficiencies {
142  TL1PerfEfficiencies()
143  : TL1Efficiencies()
144  , ratio_killed()
145  , ratio_clone()
146  , ratio_length()
147  , ratio_fakes()
148  , killed()
149  , clone()
150  , reco_length()
151  , reco_fakes()
152  , mc_length()
153  , mc_length_hits() {
154  // add total efficiency
155  AddCounter("long_fast_prim", "LongRPrim efficiency");
156  AddCounter("fast_prim", "RefPrim efficiency");
157  AddCounter("fast_sec", "RefSec efficiency");
158  AddCounter("fast", "Refset efficiency");
159  AddCounter("total", "Allset efficiency");
160  AddCounter("slow_prim", "ExtraPrim efficiency");
161  AddCounter("slow_sec", "ExtraSec efficiency");
162  AddCounter("slow", "Extra efficiency");
163  AddCounter("d0", "D0 efficiency");
164  AddCounter("short", "Short123s efficiency");
165  AddCounter("shortPion", "Short123s pion eff");
166  AddCounter("shortProton", "Short123s proton eff");
167  AddCounter("shortKaon", "Short123s kaon eff");
168  AddCounter("shortE", "Short123s e eff");
169  AddCounter("shortRest", "Short123s rest eff");
170 
171  AddCounter("fast_sec_e", "RefSec E efficiency");
172  AddCounter("fast_e", "Refset E efficiency");
173  AddCounter("total_e", "Allset E efficiency");
174  AddCounter("slow_sec_e", "ExtraSecE efficiency");
175  AddCounter("slow_e", "Extra E efficiency");
176  }
177 
178  virtual ~TL1PerfEfficiencies() {};
179 
180  virtual void AddCounter(TString shortname, TString name) {
181  TL1Efficiencies::AddCounter(shortname, name);
182  ratio_killed.AddCounter();
183  ratio_clone.AddCounter();
184  ratio_length.AddCounter();
185  ratio_fakes.AddCounter();
186  killed.AddCounter();
187  clone.AddCounter();
188  reco_length.AddCounter();
189  reco_fakes.AddCounter();
190  mc_length.AddCounter();
191  mc_length_hits.AddCounter();
192  };
193 
194  TL1PerfEfficiencies& operator+=(TL1PerfEfficiencies& a) {
195  TL1Efficiencies::operator+=(a);
196  killed += a.killed;
197  clone += a.clone;
198  reco_length += a.reco_length;
199  reco_fakes += a.reco_fakes;
200  mc_length += a.mc_length;
201  mc_length_hits += a.mc_length_hits;
202  return *this;
203  };
204 
205  void CalcEff() {
206  TL1Efficiencies::CalcEff();
207  ratio_killed = killed / mc;
208  ratio_clone = clone / mc;
209  TL1TracksCatCounters<int> allReco = reco + clone;
210  ratio_length = reco_length / allReco;
211  ratio_fakes = reco_fakes / allReco;
212  };
213 
214  void Inc(bool isReco,
215  bool isKilled,
216  double _ratio_length,
217  double _ratio_fakes,
218  int _nclones,
219  int _mc_length,
220  int _mc_length_hits,
221  TString name) {
222  TL1Efficiencies::Inc(isReco, name);
223 
224  const int index = indices[name];
225 
226  if (isKilled) killed.counters[index]++;
227  reco_length.counters[index] += _ratio_length;
228  reco_fakes.counters[index] += _ratio_fakes;
229  clone.counters[index] += _nclones;
230  mc_length.counters[index] += _mc_length;
231  mc_length_hits.counters[index] += _mc_length_hits;
232  };
233 
234  void PrintEff() {
235  L1_assert(nEvents != 0);
236 
237  cout.setf(ios::fixed);
238  cout.setf(ios::showpoint);
239  cout.precision(3);
240  cout.setf(ios::right);
241  cout << "Track category : "
242  << " Eff "
243  << " / "
244  << "Killed"
245  << " / "
246  << "Length"
247  << " / "
248  << "Fakes "
249  << " / "
250  << "Clones"
251  << " / "
252  << "All Reco"
253  << " | "
254  << " All MC "
255  << " / "
256  << "MCl(hits)"
257  << " / "
258  << "MCl(MCps)" << endl;
259 
260  int NCounters = mc.NCounters;
261  for (int iC = 0; iC < NCounters; iC++) {
262  if ((names[iC] != "D0 efficiency") || (mc.counters[iC] != 0))
263  cout
264  << names[iC] << " : " << ratio_reco.counters[iC] << " / "
265  << ratio_killed.counters
266  [iC] // tracks with aren't reco because other tracks takes their hit(-s)
267  << " / " << ratio_length.counters[iC] // nRecoMCHits/nMCHits
268  << " / " << ratio_fakes.counters[iC] // nFakeHits/nRecoAllHits
269  << " / " << ratio_clone.counters[iC] // nCloneTracks/nMCTracks
270  << " / " << setw(8) << reco.counters[iC] / double(nEvents) << " | "
271  << setw(8) << mc.counters[iC] / double(nEvents) << " / "
272  << mc_length_hits.counters[iC] / double(mc.counters[iC]) << " / "
273  << mc_length.counters[iC] / double(mc.counters[iC]) << endl;
274  }
275  cout << "Ghost probability : " << ratio_ghosts << " | " << ghosts
276  << endl;
277  };
278 
279  TL1TracksCatCounters<double> ratio_killed;
280  TL1TracksCatCounters<double> ratio_clone;
281  TL1TracksCatCounters<double> ratio_length;
282  TL1TracksCatCounters<double> ratio_fakes;
283 
284  TL1TracksCatCounters<int> killed;
285  TL1TracksCatCounters<int> clone;
286  TL1TracksCatCounters<double> reco_length;
287  TL1TracksCatCounters<double> reco_fakes;
288  TL1TracksCatCounters<int> mc_length;
289  TL1TracksCatCounters<int> mc_length_hits;
290 };
291 
292 
293 void CbmL1::EfficienciesPerformance() {
294  static TL1PerfEfficiencies L1_NTRA; // average efficiencies
295 
296  static int L1_NEVENTS = 0;
297  static double L1_CATIME = 0.0;
298 
299 
300  TL1PerfEfficiencies ntra; // efficiencies for current event
301 
302  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin();
303  rtraIt != vRTracks.end();
304  ++rtraIt) {
305  ntra.ghosts += rtraIt->IsGhost();
306  // if(rtraIt->IsGhost()){ // Debug.
307  // cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
308  // for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++ ){
309  // cout << " (" << posIt->second << " " << posIt->first << ") ";
310  // }
311  // cout << endl;
312  // }
313  }
314 
315  int sta_nhits[algo->NStations], sta_nfakes[algo->NStations];
316 
317  for (int i = 0; i < algo->NStations; i++) {
318  sta_nhits[i] = 0;
319  sta_nfakes[i] = 0;
320  }
321 
322  for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin();
323  mtraIt != vMCTracks.end();
324  mtraIt++) {
325  CbmL1MCTrack& mtra = *(mtraIt);
326 
327 
328  // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
329  if (!mtra.IsReconstructable() && !mtra.IsAdditional()) continue;
330 
331  // -- find used constans --
332  // is track reconstructed
333  const bool reco = (mtra.IsReconstructed());
334  // is track killed. At least one hit of it belong to some recoTrack
335  const bool killed = !reco && mtra.IsDisturbed();
336  // ration length for current mcTrack
337  vector<CbmL1Track*>& rTracks =
338  mtra.GetRecoTracks(); // for length calculations
339  double ratio_length = 0;
340  double ratio_fakes = 0;
341  double mc_length_hits = mtra.NStations();
342 
343 
344  int mc_length = mtra.NMCStations();
345  if (reco) {
346  for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
347  ratio_length += static_cast<double>(rTracks[irt]->GetNOfHits())
348  * rTracks[irt]->GetMaxPurity() / mc_length_hits;
349  ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
350  }
351  }
352 
353  // number of clones
354  int nclones = 0;
355  if (reco) nclones = mtra.GetNClones();
356  // if (nclones){ // Debug. Look at clones
357  // for (int irt = 0; irt < rTracks.size(); irt++){
358  // const int ista = vHitStore[rTracks[irt]->StsHits[0]].iStation;
359  // cout << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
360  // }
361  // cout << mtra.NStations() << endl;
362  // }
363 
364  if (mtra.IsAdditional()) { // short
365  ntra.Inc(reco,
366  killed,
367  ratio_length,
368  ratio_fakes,
369  nclones,
370  mc_length,
371  mc_length_hits,
372  "short");
373  switch (mtra.pdg) {
374  case 11:
375  case -11:
376  ntra.Inc(reco,
377  killed,
378  ratio_length,
379  ratio_fakes,
380  nclones,
381  mc_length,
382  mc_length_hits,
383  "shortE");
384  break;
385  case 211:
386  case -211:
387  ntra.Inc(reco,
388  killed,
389  ratio_length,
390  ratio_fakes,
391  nclones,
392  mc_length,
393  mc_length_hits,
394  "shortPion");
395  break;
396  case 321:
397  case -321:
398  ntra.Inc(reco,
399  killed,
400  ratio_length,
401  ratio_fakes,
402  nclones,
403  mc_length,
404  mc_length_hits,
405  "shortKaon");
406  break;
407  case 2212:
408  case -2212:
409  ntra.Inc(reco,
410  killed,
411  ratio_length,
412  ratio_fakes,
413  nclones,
414  mc_length,
415  mc_length_hits,
416  "shortProton");
417  break;
418  default:
419  ntra.Inc(reco,
420  killed,
421  ratio_length,
422  ratio_fakes,
423  nclones,
424  mc_length,
425  mc_length_hits,
426  "shortRest");
427  }
428  } else { // separate all efficiecies from short eff
429 
430 
431  ntra.Inc(reco,
432  killed,
433  ratio_length,
434  ratio_fakes,
435  nclones,
436  mc_length,
437  mc_length_hits,
438  "total");
439 
440  if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
441  ntra.Inc(reco,
442  killed,
443  ratio_length,
444  ratio_fakes,
445  nclones,
446  mc_length,
447  mc_length_hits,
448  "d0");
449  }
450 
451  if (mtra.p > CbmL1Constants::MinRefMom) { // reference tracks
452  ntra.Inc(reco,
453  killed,
454  ratio_length,
455  ratio_fakes,
456  nclones,
457  mc_length,
458  mc_length_hits,
459  "fast");
460 
461  if (mtra.IsPrimary()) { // reference primary
462  if (mtra.NStations() == NStation) { // long reference primary
463  ntra.Inc(reco,
464  killed,
465  ratio_length,
466  ratio_fakes,
467  nclones,
468  mc_length,
469  mc_length_hits,
470  "long_fast_prim");
471  }
472  ntra.Inc(reco,
473  killed,
474  ratio_length,
475  ratio_fakes,
476  nclones,
477  mc_length,
478  mc_length_hits,
479  "fast_prim");
480  } else { // reference secondary
481  ntra.Inc(reco,
482  killed,
483  ratio_length,
484  ratio_fakes,
485  nclones,
486  mc_length,
487  mc_length_hits,
488  "fast_sec");
489  }
490  } else { // extra set of tracks
491  ntra.Inc(reco,
492  killed,
493  ratio_length,
494  ratio_fakes,
495  nclones,
496  mc_length,
497  mc_length_hits,
498  "slow");
499 
500  if (mtra.IsPrimary()) { // extra primary
501  ntra.Inc(reco,
502  killed,
503  ratio_length,
504  ratio_fakes,
505  nclones,
506  mc_length,
507  mc_length_hits,
508  "slow_prim");
509  } else {
510  ntra.Inc(reco,
511  killed,
512  ratio_length,
513  ratio_fakes,
514  nclones,
515  mc_length,
516  mc_length_hits,
517  "slow_sec");
518  }
519  } // if extra
520 
521  if (mtra.pdg == 11 || mtra.pdg == -11) {
522  ntra.Inc(reco,
523  killed,
524  ratio_length,
525  ratio_fakes,
526  nclones,
527  mc_length,
528  mc_length_hits,
529  "total_e");
530 
531  if (mtra.p > CbmL1Constants::MinRefMom) { // reference tracks
532  ntra.Inc(reco,
533  killed,
534  ratio_length,
535  ratio_fakes,
536  nclones,
537  mc_length,
538  mc_length_hits,
539  "fast_e");
540 
541  if (mtra.IsPrimary()) { // reference primary
542  } else { // reference secondary
543  ntra.Inc(reco,
544  killed,
545  ratio_length,
546  ratio_fakes,
547  nclones,
548  mc_length,
549  mc_length_hits,
550  "fast_sec_e");
551  }
552  } else { // extra set of tracks
553  ntra.Inc(reco,
554  killed,
555  ratio_length,
556  ratio_fakes,
557  nclones,
558  mc_length,
559  mc_length_hits,
560  "slow_e");
561 
562  if (mtra.IsPrimary()) { // extra primary
563  } else {
564  ntra.Inc(reco,
565  killed,
566  ratio_length,
567  ratio_fakes,
568  nclones,
569  mc_length,
570  mc_length_hits,
571  "slow_sec_e");
572  }
573  } // if extra
574  }
575  }
576 
577  } // for mcTracks
578 
579  L1_CATIME += algo->CATime;
580  L1_NEVENTS++;
581  ntra.IncNEvents();
582  L1_NTRA += ntra;
583 
584  ntra.CalcEff();
585  L1_NTRA.CalcEff();
586 
587  // cout.precision(3);
588  if (fVerbose) {
589  if (fVerbose > 1) {
590  ntra.PrintEff();
591  cout << "Number of true and fake hits in stations: " << endl;
592  for (int i = 0; i < algo->NStations; i++) {
593  cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
594  }
595  cout << endl;
596  } // fVerbose > 1
597  cout << endl;
598  cout << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS " << endl
599  << endl;
600  L1_NTRA.PrintEff();
601  cout << "MC tracks/event found : "
602  << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]])
603  / double(L1_NEVENTS))
604  << endl;
605  cout << endl;
606  cout << "CA Track Finder: " << L1_CATIME / L1_NEVENTS << " s/ev" << endl
607  << endl;
608  }
609 } // void CbmL1::Performance()
610 
611 void CbmL1::
612  HistoPerformance() // TODO: check if works correctly. Change vHitRef on match data in CbmL1**Track classes
613 {
614 
615  //CbmKF &KF = *CbmKF::Instance();
616 
617  static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom,
618  *p_eff_d0_vs_mom, *p_eff_prim_vs_theta, *p_eff_all_vs_pt, *p_eff_prim_vs_pt,
619  *p_eff_all_vs_nhits, *p_eff_prim_vs_nhits, *p_eff_sec_vs_nhits;
620 
621  static TH1F *h_reg_MCmom, *h_acc_MCmom, *h_reco_MCmom, *h_ghost_Rmom;
622  static TH1F *h_reg_prim_MCmom, *h_acc_prim_MCmom, *h_reco_prim_MCmom;
623  static TH1F *h_reg_sec_MCmom, *h_acc_sec_MCmom, *h_reco_sec_MCmom;
624 
625  static TH1F* h_acc_mom_short123s;
626 
627  static TH1F *h_reg_mom_prim, *h_reg_mom_sec, *h_reg_nhits_prim,
628  *h_reg_nhits_sec;
629  static TH1F *h_acc_mom_prim, *h_acc_mom_sec, *h_acc_nhits_prim,
630  *h_acc_nhits_sec;
631  static TH1F *h_reco_mom_prim, *h_reco_mom_sec, *h_reco_nhits_prim,
632  *h_reco_nhits_sec;
633  static TH1F *h_rest_mom_prim, *h_rest_mom_sec, *h_rest_nhits_prim,
634  *h_rest_nhits_sec;
635 
636  //static TH1F *h_hit_density[10];
637 
638  static TH1F *h_ghost_mom, *h_ghost_nhits, *h_ghost_fstation, *h_ghost_chi2,
639  *h_ghost_prob, *h_ghost_tx, *h_ghost_ty;
640  static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station,
641  *h_reco_chi2, *h_reco_prob, *h_rest_prob, *h_reco_clean, *h_reco_time;
642  static TProfile *h_reco_timeNtr, *h_reco_timeNhit;
643  static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit;
644  static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;
645 
646  static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station,
647  *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
648 
649  static TH1F *h_mcp, *h_nmchits;
650  // static TH1F *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
651 
652  // static TH2F *P_vs_P ;
653 
654  static TH2F *h2_vertex, *h2_prim_vertex, *h2_sec_vertex;
655  //static TH3F *h3_vertex, *h3_prim_vertex, *h3_sec_vertex;
656 
657  static TH2F *h2_reg_nhits_vs_mom_prim, *h2_reg_nhits_vs_mom_sec,
658  *h2_reg_fstation_vs_mom_prim, *h2_reg_fstation_vs_mom_sec,
659  *h2_reg_lstation_vs_fstation_prim, *h2_reg_lstation_vs_fstation_sec;
660  static TH2F *h2_acc_nhits_vs_mom_prim, *h2_acc_nhits_vs_mom_sec,
661  *h2_acc_fstation_vs_mom_prim, *h2_acc_fstation_vs_mom_sec,
662  *h2_acc_lstation_vs_fstation_prim, *h2_acc_lstation_vs_fstation_sec;
663  static TH2F *h2_reco_nhits_vs_mom_prim, *h2_reco_nhits_vs_mom_sec,
664  *h2_reco_fstation_vs_mom_prim, *h2_reco_fstation_vs_mom_sec,
665  *h2_reco_lstation_vs_fstation_prim, *h2_reco_lstation_vs_fstation_sec;
666  static TH2F *h2_ghost_nhits_vs_mom, *h2_ghost_fstation_vs_mom,
667  *h2_ghost_lstation_vs_fstation;
668  static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec,
669  *h2_rest_fstation_vs_mom_prim, *h2_rest_fstation_vs_mom_sec,
670  *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
671 
672  static bool first_call = 1;
673 
674  if (first_call) {
675  first_call = 0;
676 
677  TDirectory* curdir = gDirectory;
678  gDirectory = histodir;
679  gDirectory->cd("L1");
680 
681  p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom",
682  "AllSet Efficiency vs Momentum",
683  100,
684  0.0,
685  5.0,
686  0.0,
687  100.0);
688  p_eff_prim_vs_mom = new TProfile("p_eff_prim_vs_mom",
689  "Primary Set Efficiency vs Momentum",
690  100,
691  0.0,
692  5.0,
693  0.0,
694  100.0);
695  p_eff_sec_vs_mom = new TProfile("p_eff_sec_vs_mom",
696  "Secondary Set Efficiency vs Momentum",
697  100,
698  0.0,
699  5.0,
700  0.0,
701  100.0);
702  p_eff_d0_vs_mom = new TProfile("p_eff_d0_vs_mom",
703  "D0 Secondary Tracks Efficiency vs Momentum",
704  150,
705  0.0,
706  15.0,
707  0.0,
708  100.0);
709  p_eff_prim_vs_theta = new TProfile("p_eff_prim_vs_theta",
710  "All Primary Set Efficiency vs Theta",
711  100,
712  0.0,
713  30.0,
714  0.0,
715  100.0);
716  p_eff_all_vs_pt = new TProfile(
717  "p_eff_all_vs_pt", "AllSet Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
718  p_eff_prim_vs_pt = new TProfile("p_eff_prim_vs_pt",
719  "Primary Set Efficiency vs Pt",
720  100,
721  0.0,
722  5.0,
723  0.0,
724  100.0);
725 
726  p_eff_all_vs_nhits = new TProfile("p_eff_all_vs_nhits",
727  "AllSet Efficiency vs NMCHits",
728  8,
729  3.0,
730  11.0,
731  0.0,
732  100.0);
733  p_eff_prim_vs_nhits = new TProfile("p_eff_prim_vs_nhits",
734  "PrimSet Efficiency vs NMCHits",
735  8,
736  3.0,
737  11.0,
738  0.0,
739  100.0);
740  p_eff_sec_vs_nhits = new TProfile("p_eff_sec_vs_nhits",
741  "SecSet Efficiency vs NMCHits",
742  8,
743  3.0,
744  11.0,
745  0.0,
746  100.0);
747 
748  h_reg_MCmom =
749  new TH1F("h_reg_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
750  h_acc_MCmom =
751  new TH1F("h_acc_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
752  h_reco_MCmom =
753  new TH1F("h_reco_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
754  h_ghost_Rmom = new TH1F("h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0);
755  h_reg_prim_MCmom = new TH1F(
756  "h_reg_prim_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
757  h_acc_prim_MCmom =
758  new TH1F("h_acc_prim_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
759  h_reco_prim_MCmom =
760  new TH1F("h_reco_prim_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
761  h_reg_sec_MCmom = new TH1F(
762  "h_reg_sec_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
763  h_acc_sec_MCmom =
764  new TH1F("h_acc_sec_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
765  h_reco_sec_MCmom =
766  new TH1F("h_reco_sec_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
767 
768  h_acc_mom_short123s =
769  new TH1F("h_acc_mom_short123s",
770  "Momentum of accepted tracks with 3 hits on first stations",
771  500,
772  0.0,
773  5.0);
774 
775  h_reg_mom_prim = new TH1F(
776  "h_reg_mom_prim", "Momentum of registered primary tracks", 500, 0.0, 5.0);
777  h_reg_mom_sec = new TH1F("h_reg_mom_sec",
778  "Momentum of registered secondary tracks",
779  500,
780  0.0,
781  5.0);
782  h_acc_mom_prim = new TH1F(
783  "h_acc_mom_prim", "Momentum of accepted primary tracks", 500, 0.0, 5.0);
784  h_acc_mom_sec = new TH1F(
785  "h_acc_mom_sec", "Momentum of accepted secondary tracks", 500, 0.0, 5.0);
786  h_reco_mom_prim = new TH1F("h_reco_mom_prim",
787  "Momentum of reconstructed primary tracks",
788  500,
789  0.0,
790  5.0);
791  h_reco_mom_sec = new TH1F("h_reco_mom_sec",
792  "Momentum of reconstructed secondary tracks",
793  500,
794  0.0,
795  5.0);
796  h_rest_mom_prim = new TH1F(
797  "h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
798  h_rest_mom_sec = new TH1F("h_rest_mom_sec",
799  "Momentum of not found secondary tracks",
800  500,
801  0.0,
802  5.0);
803 
804  h_reg_nhits_prim = new TH1F("h_reg_nhits_prim",
805  "Number of hits in registered primary track",
806  51,
807  -0.1,
808  10.1);
809  h_reg_nhits_sec = new TH1F("h_reg_nhits_sec",
810  "Number of hits in registered secondary track",
811  51,
812  -0.1,
813  10.1);
814  h_acc_nhits_prim = new TH1F("h_acc_nhits_prim",
815  "Number of hits in accepted primary track",
816  51,
817  -0.1,
818  10.1);
819  h_acc_nhits_sec = new TH1F("h_acc_nhits_sec",
820  "Number of hits in accepted secondary track",
821  51,
822  -0.1,
823  10.1);
824  h_reco_nhits_prim =
825  new TH1F("h_reco_nhits_prim",
826  "Number of hits in reconstructed primary track",
827  51,
828  -0.1,
829  10.1);
830  h_reco_nhits_sec =
831  new TH1F("h_reco_nhits_sec",
832  "Number of hits in reconstructed secondary track",
833  51,
834  -0.1,
835  10.1);
836  h_rest_nhits_prim = new TH1F("h_rest_nhits_prim",
837  "Number of hits in not found primary track",
838  51,
839  -0.1,
840  10.1);
841  h_rest_nhits_sec = new TH1F("h_rest_nhits_sec",
842  "Number of hits in not found secondary track",
843  51,
844  -0.1,
845  10.1);
846 
847  h_ghost_mom =
848  new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
849  h_ghost_nhits = new TH1F(
850  "h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
851  h_ghost_fstation = new TH1F(
852  "h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
853  h_ghost_chi2 =
854  new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
855  h_ghost_prob =
856  new TH1F("h_ghost_prob", "Prob of ghost track", 505, 0., 1.01);
857  h_ghost_r =
858  new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0);
859  h_ghost_tx = new TH1F(
860  "h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
861  h_ghost_ty = new TH1F(
862  "h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);
863 
864  h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
865  h_reco_nhits =
866  new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
867  h_reco_station =
868  new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
869  h_reco_chi2 =
870  new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
871  h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 505, 0., 1.01);
872  h_rest_prob =
873  new TH1F("h_rest_prob", "Prob of reco rest track", 505, 0., 1.01);
874  h_reco_clean =
875  new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5);
876  h_reco_time =
877  new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
878  h_reco_timeNtr = new TProfile("h_reco_timeNtr",
879  "CA Track Finder Time (s/ev) vs N Tracks",
880  200,
881  0.0,
882  1000.0);
883  h_reco_timeNhit = new TProfile("h_reco_timeNhit",
884  "CA Track Finder Time (s/ev) vs N Hits",
885  200,
886  0.0,
887  30000.0);
888  h_reco_fakeNtr = new TProfile(
889  "h_reco_fakeNtr", "N Fake hits vs N Tracks", 200, 0.0, 1000.0);
890  h_reco_fakeNhit = new TProfile(
891  "h_reco_fakeNhit", "N Fake hits vs N Hits", 200, 0.0, 30000.0);
892 
893  h_reco_d0_mom =
894  new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);
895 
896  // h_hit_density[0] = new TH1F("h_hit_density0", "Hit density station 1", 50, 0.0, 5.00);
897  // h_hit_density[1] = new TH1F("h_hit_density1", "Hit density station 2", 100, 0.0, 10.00);
898  // h_hit_density[2] = new TH1F("h_hit_density2", "Hit density station 3", 140, 0.0, 13.99);
899  // h_hit_density[3] = new TH1F("h_hit_density3", "Hit density station 4", 180, 0.0, 18.65);
900  // h_hit_density[4] = new TH1F("h_hit_density4", "Hit density station 5", 230, 0.0, 23.32);
901  // h_hit_density[5] = new TH1F("h_hit_density5", "Hit density station 6", 270, 0.0, 27.98);
902  // h_hit_density[6] = new TH1F("h_hit_density6", "Hit density station 7", 340, 0.0, 34.97);
903  // h_hit_density[7] = new TH1F("h_hit_density7", "Hit density station 8", 460, 0.0, 46.63);
904  // h_hit_density[8] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
905  // h_hit_density[9] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
906  // h_hit_density[10] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
907 
908  h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
909  h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);
910 
911  h_sec_r =
912  new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);
913 
914  h_notfound_mom =
915  new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
916  h_notfound_nhits = new TH1F(
917  "h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
918  h_notfound_station = new TH1F(
919  "h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
920  h_notfound_r = new TH1F(
921  "h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
922  h_notfound_tx = new TH1F(
923  "h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
924  h_notfound_ty = new TH1F(
925  "h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);
926 
927  /*
928  h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
929  h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
930  VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
931  VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
932 */
933  h_mcp = new TH1F("h_mcp", "MC P ", 500, 0.0, 5.0);
934  /*
935  MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
936  MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
937  MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
938 */
939  h_nmchits = new TH1F("h_nmchits", "N STS hits on MC track", 50, 0.0, 10.0);
940 
941  // P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);
942 
943  h2_vertex = new TH2F(
944  "h2_vertex", "2D vertex distribution", 105, -5., 100., 100, -50., 50.);
945  h2_prim_vertex = new TH2F("h2_primvertex",
946  "2D primary vertex distribution",
947  105,
948  -5.,
949  100.,
950  100,
951  -50.,
952  50.);
953  h2_sec_vertex = new TH2F("h2_sec_vertex",
954  "2D secondary vertex distribution",
955  105,
956  -5.,
957  100.,
958  100,
959  -50.,
960  50.);
961 
962  //h3_vertex = new TH3F("h3_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
963  //h3_prim_vertex = new TH3F("h3_primvertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
964  //h3_sec_vertex = new TH3F("h3_sec_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
965 
966  h2_reg_nhits_vs_mom_prim =
967  new TH2F("h2_reg_nhits_vs_mom_prim",
968  "NHits vs. Momentum (Reg. Primary Tracks)",
969  51,
970  -0.05,
971  5.05,
972  11,
973  -0.5,
974  10.5);
975  h2_reg_nhits_vs_mom_sec =
976  new TH2F("h2_reg_nhits_vs_mom_sec",
977  "NHits vs. Momentum (Reg. Secondary Tracks)",
978  51,
979  -0.05,
980  5.05,
981  11,
982  -0.5,
983  10.5);
984  h2_acc_nhits_vs_mom_prim =
985  new TH2F("h2_acc_nhits_vs_mom_prim",
986  "NHits vs. Momentum (Acc. Primary Tracks)",
987  51,
988  -0.05,
989  5.05,
990  11,
991  -0.5,
992  10.5);
993  h2_acc_nhits_vs_mom_sec =
994  new TH2F("h2_acc_nhits_vs_mom_sec",
995  "NHits vs. Momentum (Acc. Secondary Tracks)",
996  51,
997  -0.05,
998  5.05,
999  11,
1000  -0.5,
1001  10.5);
1002  h2_reco_nhits_vs_mom_prim =
1003  new TH2F("h2_reco_nhits_vs_mom_prim",
1004  "NHits vs. Momentum (Reco Primary Tracks)",
1005  51,
1006  -0.05,
1007  5.05,
1008  11,
1009  -0.5,
1010  10.5);
1011  h2_reco_nhits_vs_mom_sec =
1012  new TH2F("h2_reco_nhits_vs_mom_sec",
1013  "NHits vs. Momentum (Reco Secondary Tracks)",
1014  51,
1015  -0.05,
1016  5.05,
1017  11,
1018  -0.5,
1019  10.5);
1020  h2_ghost_nhits_vs_mom = new TH2F("h2_ghost_nhits_vs_mom",
1021  "NHits vs. Momentum (Ghost Tracks)",
1022  51,
1023  -0.05,
1024  5.05,
1025  11,
1026  -0.5,
1027  10.5);
1028  h2_rest_nhits_vs_mom_prim =
1029  new TH2F("h2_rest_nhits_vs_mom_prim",
1030  "NHits vs. Momentum (!Found Primary Tracks)",
1031  51,
1032  -0.05,
1033  5.05,
1034  11,
1035  -0.5,
1036  10.5);
1037  h2_rest_nhits_vs_mom_sec =
1038  new TH2F("h2_rest_nhits_vs_mom_sec",
1039  "NHits vs. Momentum (!Found Secondary Tracks)",
1040  51,
1041  -0.05,
1042  5.05,
1043  11,
1044  -0.5,
1045  10.5);
1046 
1047  h2_reg_fstation_vs_mom_prim =
1048  new TH2F("h2_reg_fstation_vs_mom_prim",
1049  "First Station vs. Momentum (Reg. Primary Tracks)",
1050  51,
1051  -0.05,
1052  5.05,
1053  11,
1054  -0.5,
1055  10.5);
1056  h2_reg_fstation_vs_mom_sec =
1057  new TH2F("h2_reg_fstation_vs_mom_sec",
1058  "First Station vs. Momentum (Reg. Secondary Tracks)",
1059  51,
1060  -0.05,
1061  5.05,
1062  11,
1063  -0.5,
1064  10.5);
1065  h2_acc_fstation_vs_mom_prim =
1066  new TH2F("h2_acc_fstation_vs_mom_prim",
1067  "First Station vs. Momentum (Acc. Primary Tracks)",
1068  51,
1069  -0.05,
1070  5.05,
1071  11,
1072  -0.5,
1073  10.5);
1074  h2_acc_fstation_vs_mom_sec =
1075  new TH2F("h2_acc_fstation_vs_mom_sec",
1076  "First Station vs. Momentum (Acc. Secondary Tracks)",
1077  51,
1078  -0.05,
1079  5.05,
1080  11,
1081  -0.5,
1082  10.5);
1083  h2_reco_fstation_vs_mom_prim =
1084  new TH2F("h2_reco_fstation_vs_mom_prim",
1085  "First Station vs. Momentum (Reco Primary Tracks)",
1086  51,
1087  -0.05,
1088  5.05,
1089  11,
1090  -0.5,
1091  10.5);
1092  h2_reco_fstation_vs_mom_sec =
1093  new TH2F("h2_reco_fstation_vs_mom_sec",
1094  "First Station vs. Momentum (Reco Secondary Tracks)",
1095  51,
1096  -0.05,
1097  5.05,
1098  11,
1099  -0.5,
1100  10.5);
1101  h2_ghost_fstation_vs_mom =
1102  new TH2F("h2_ghost_fstation_vs_mom",
1103  "First Station vs. Momentum (Ghost Tracks)",
1104  51,
1105  -0.05,
1106  5.05,
1107  11,
1108  -0.5,
1109  10.5);
1110  h2_rest_fstation_vs_mom_prim =
1111  new TH2F("h2_rest_fstation_vs_mom_prim",
1112  "First Station vs. Momentum (!Found Primary Tracks)",
1113  51,
1114  -0.05,
1115  5.05,
1116  11,
1117  -0.5,
1118  10.5);
1119  h2_rest_fstation_vs_mom_sec =
1120  new TH2F("h2_rest_fstation_vs_mom_sec",
1121  "First Station vs. Momentum (!Found Secondary Tracks)",
1122  51,
1123  -0.05,
1124  5.05,
1125  11,
1126  -0.5,
1127  10.5);
1128 
1129  h2_reg_lstation_vs_fstation_prim =
1130  new TH2F("h2_reg_lstation_vs_fstation_prim",
1131  "Last vs. First Station (Reg. Primary Tracks)",
1132  11,
1133  -0.5,
1134  10.5,
1135  11,
1136  -0.5,
1137  10.5);
1138  h2_reg_lstation_vs_fstation_sec =
1139  new TH2F("h2_reg_lstation_vs_fstation_sec",
1140  "Last vs. First Station (Reg. Secondary Tracks)",
1141  11,
1142  -0.5,
1143  10.5,
1144  11,
1145  -0.5,
1146  10.5);
1147  h2_acc_lstation_vs_fstation_prim =
1148  new TH2F("h2_acc_lstation_vs_fstation_prim",
1149  "Last vs. First Station (Acc. Primary Tracks)",
1150  11,
1151  -0.5,
1152  10.5,
1153  11,
1154  -0.5,
1155  10.5);
1156  h2_acc_lstation_vs_fstation_sec =
1157  new TH2F("h2_acc_lstation_vs_fstation_sec",
1158  "Last vs. First Station (Acc. Secondary Tracks)",
1159  11,
1160  -0.5,
1161  10.5,
1162  11,
1163  -0.5,
1164  10.5);
1165  h2_reco_lstation_vs_fstation_prim =
1166  new TH2F("h2_reco_lstation_vs_fstation_prim",
1167  "Last vs. First Station (Reco Primary Tracks)",
1168  11,
1169  -0.5,
1170  10.5,
1171  11,
1172  -0.5,
1173  10.5);
1174  h2_reco_lstation_vs_fstation_sec =
1175  new TH2F("h2_reco_lstation_vs_fstation_sec",
1176  "Last vs. First Station (Reco Secondary Tracks)",
1177  11,
1178  -0.5,
1179  10.5,
1180  11,
1181  -0.5,
1182  10.5);
1183  h2_ghost_lstation_vs_fstation =
1184  new TH2F("h2_ghost_lstation_vs_fstation",
1185  "Last vs. First Station (Ghost Tracks)",
1186  11,
1187  -0.5,
1188  10.5,
1189  11,
1190  -0.5,
1191  10.5);
1192  h2_rest_lstation_vs_fstation_prim =
1193  new TH2F("h2_rest_lstation_vs_fstation_prim",
1194  "Last vs. First Station (!Found Primary Tracks)",
1195  11,
1196  -0.5,
1197  10.5,
1198  11,
1199  -0.5,
1200  10.5);
1201  h2_rest_lstation_vs_fstation_sec =
1202  new TH2F("h2_rest_lstation_vs_fstation_sec",
1203  "Last vs. First Station (!Found Secondary Tracks)",
1204  11,
1205  -0.5,
1206  10.5,
1207  11,
1208  -0.5,
1209  10.5);
1210 
1211  //maindir->cd();
1212 
1213  // ----- Create list of all histogram pointers
1214 
1215  gDirectory = curdir;
1216 
1217  } // first_call
1218 
1219 
1220  static int NEvents = 0;
1221  if (NEvents > 0) {
1222  h_reg_MCmom->Scale(NEvents);
1223  h_acc_MCmom->Scale(NEvents);
1224  h_reco_MCmom->Scale(NEvents);
1225  h_ghost_Rmom->Scale(NEvents);
1226  h_reg_prim_MCmom->Scale(NEvents);
1227  h_acc_prim_MCmom->Scale(NEvents);
1228  h_reco_prim_MCmom->Scale(NEvents);
1229  h_reg_sec_MCmom->Scale(NEvents);
1230  h_acc_sec_MCmom->Scale(NEvents);
1231  h_reco_sec_MCmom->Scale(NEvents);
1232  }
1233  NEvents++;
1234 
1235  // //hit density calculation: h_hit_density[10]
1236  //
1237  // for (vector<CbmL1HitStore>::iterator hIt = vHitStore.begin(); hIt != vHitStore.end(); ++hIt){
1238  // float x = hIt->x;
1239  // float y = hIt->y;
1240  // float r = sqrt(x*x+y*y);
1241  // h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
1242  // }
1243 
1244  //
1245  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin();
1246  rtraIt != vRTracks.end();
1247  ++rtraIt) {
1248  CbmL1Track* prtra = &(*rtraIt);
1249  if ((prtra->StsHits).size() < 1) continue;
1250  { // fill histos
1251  if (fabs(prtra->T[4]) > 1.e-10) h_reco_mom->Fill(fabs(1.0 / prtra->T[4]));
1252  h_reco_nhits->Fill((prtra->StsHits).size());
1253  CbmL1HitStore& mh = vHitStore[prtra->StsHits[0]];
1254  h_reco_station->Fill(mh.iStation);
1255  }
1256 
1257  h_reco_clean->Fill(prtra->GetMaxPurity());
1258 
1259  if (prtra->NDF > 0) {
1260  if (prtra->IsGhost()) {
1261  h_ghost_chi2->Fill(prtra->chi2 / prtra->NDF);
1262  h_ghost_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1263  } else {
1264  if (prtra->GetMCTrack()[0].IsReconstructable()) {
1265  h_reco_chi2->Fill(prtra->chi2 / prtra->NDF);
1266  h_reco_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1267  } else {
1268  // h_rest_chi2->Fill(prtra->chi2/prtra->NDF);
1269  h_rest_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
1270  }
1271  }
1272  }
1273 
1274 
1275  // fill ghost histos
1276  if (prtra->IsGhost()) {
1277  if (fabs(prtra->T[4]) > 1.e-10) {
1278  h_ghost_mom->Fill(fabs(1.0 / prtra->T[4]));
1279  h_ghost_Rmom->Fill(fabs(1.0 / prtra->T[4]));
1280  }
1281  h_ghost_nhits->Fill((prtra->StsHits).size());
1282  CbmL1HitStore& h1 = vHitStore[prtra->StsHits[0]];
1283  CbmL1HitStore& h2 = vHitStore[prtra->StsHits[1]];
1284  h_ghost_fstation->Fill(h1.iStation);
1285  h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
1286  double z1 = algo->vStations[h1.iStation].z[0];
1287  double z2 = algo->vStations[h2.iStation].z[0];
1288  if (fabs(z2 - z1) > 1.e-4) {
1289  h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
1290  h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
1291  }
1292 
1293  if (fabs(prtra->T[4]) > 1.e-10)
1294  h2_ghost_nhits_vs_mom->Fill(fabs(1.0 / prtra->T[4]),
1295  (prtra->StsHits).size());
1296  CbmL1HitStore& hf = vHitStore[prtra->StsHits[0]];
1297  CbmL1HitStore& hl =
1298  vHitStore[prtra->StsHits[(prtra->StsHits).size() - 1]];
1299  if (fabs(prtra->T[4]) > 1.e-10)
1300  h2_ghost_fstation_vs_mom->Fill(fabs(1.0 / prtra->T[4]),
1301  hf.iStation + 1);
1302  if (hl.iStation >= hf.iStation)
1303  h2_ghost_lstation_vs_fstation->Fill(hf.iStation + 1, hl.iStation + 1);
1304  }
1305 
1306  } // for reco tracks
1307 
1308  int mc_total = 0; // total amount of reconstructable mcTracks
1309  for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin();
1310  mtraIt != vMCTracks.end();
1311  mtraIt++) {
1312  CbmL1MCTrack& mtra = *(mtraIt);
1313  // if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
1314 
1315  // No Sts hits?
1316  int nmchits = mtra.StsHits.size();
1317  if (nmchits == 0) continue;
1318 
1319  double momentum = mtra.p;
1320  double pt = sqrt(mtra.px * mtra.px + mtra.py * mtra.py);
1321  double theta = acos(mtra.pz / mtra.p) * 180 / 3.1415;
1322 
1323  h_mcp->Fill(momentum);
1324  h_nmchits->Fill(nmchits);
1325 
1326  int nSta = mtra.NStations();
1327 
1328  CbmL1HitStore& fh = vHitStore[*(mtra.StsHits.begin())];
1329  CbmL1HitStore& lh = vHitStore[*(mtra.StsHits.rbegin())];
1330 
1331  h_reg_MCmom->Fill(momentum);
1332  if (mtra.IsPrimary()) {
1333  h_reg_mom_prim->Fill(momentum);
1334  h_reg_prim_MCmom->Fill(momentum);
1335  h_reg_nhits_prim->Fill(nSta);
1336  h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
1337  h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1338  if (lh.iStation >= fh.iStation)
1339  h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation + 1,
1340  lh.iStation + 1);
1341  } else {
1342  h_reg_mom_sec->Fill(momentum);
1343  h_reg_sec_MCmom->Fill(momentum);
1344  h_reg_nhits_sec->Fill(nSta);
1345  if (momentum > 0.01) {
1346  h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
1347  h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1348  if (lh.iStation >= fh.iStation)
1349  h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation + 1,
1350  lh.iStation + 1);
1351  }
1352  }
1353 
1354  if (mtra.IsAdditional()) { h_acc_mom_short123s->Fill(momentum); }
1355 
1356  if (!mtra.IsReconstructable()) continue;
1357  mc_total++;
1358 
1359  h_acc_MCmom->Fill(momentum);
1360  if (mtra.IsPrimary()) {
1361  h_acc_mom_prim->Fill(momentum);
1362  h_acc_prim_MCmom->Fill(momentum);
1363  h_acc_nhits_prim->Fill(nSta);
1364  h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
1365  h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1366  if (lh.iStation >= fh.iStation)
1367  h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation + 1,
1368  lh.iStation + 1);
1369  } else {
1370  h_acc_mom_sec->Fill(momentum);
1371  h_acc_sec_MCmom->Fill(momentum);
1372  h_acc_nhits_sec->Fill(nSta);
1373  if (momentum > 0.01) {
1374  h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
1375  h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1376  if (lh.iStation >= fh.iStation)
1377  h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation + 1,
1378  lh.iStation + 1);
1379  }
1380  }
1381 
1382 
1383  // vertex distribution of primary and secondary tracks
1384  h2_vertex->Fill(mtra.z, mtra.y);
1385  //h3_vertex->Fill(mtra.z, mtra.x, mtra.y);
1386  if (mtra.mother_ID < 0) { // primary
1387  h2_prim_vertex->Fill(mtra.z, mtra.y);
1388  //h3_prim_vertex->Fill(mtra.z, mtra.x, mtra.y);
1389  } else {
1390  h2_sec_vertex->Fill(mtra.z, mtra.y);
1391  //h3_sec_vertex->Fill(mtra.z, mtra.x, mtra.y);
1392  }
1393 
1394 
1395  int iph = mtra.StsHits[0];
1396  CbmL1HitStore& ph = vHitStore[iph];
1397 
1398  h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1399  if (fabs(mtra.pz) > 1.e-8) {
1400  h_tx->Fill(mtra.px / mtra.pz);
1401  h_ty->Fill(mtra.py / mtra.pz);
1402  }
1403 
1404  // reconstructed tracks
1405  bool reco = (mtra.IsReconstructed());
1406  int nMCHits = mtra.NStations();
1407 
1408  if (reco) {
1409  p_eff_all_vs_mom->Fill(momentum, 100.0);
1410  p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
1411  p_eff_all_vs_pt->Fill(pt, 100.0);
1412  h_reco_MCmom->Fill(momentum);
1413  if (mtra.mother_ID < 0) { // primary
1414  p_eff_prim_vs_mom->Fill(momentum, 100.0);
1415  p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
1416  p_eff_prim_vs_pt->Fill(pt, 100.0);
1417  p_eff_prim_vs_theta->Fill(theta, 100.0);
1418  } else {
1419  p_eff_sec_vs_mom->Fill(momentum, 100.0);
1420  p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
1421  }
1422  if (mtra.mother_ID < 0) { // primary
1423  h_reco_mom_prim->Fill(momentum);
1424  h_reco_prim_MCmom->Fill(momentum);
1425  h_reco_nhits_prim->Fill(nSta);
1426  h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
1427  h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1428  if (lh.iStation >= fh.iStation)
1429  h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation + 1,
1430  lh.iStation + 1);
1431  } else {
1432  h_reco_mom_sec->Fill(momentum);
1433  h_reco_sec_MCmom->Fill(momentum);
1434  h_reco_nhits_sec->Fill(nSta);
1435  if (momentum > 0.01) {
1436  h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
1437  h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1438  if (lh.iStation >= fh.iStation)
1439  h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation + 1,
1440  lh.iStation + 1);
1441  }
1442  }
1443  } else {
1444  h_notfound_mom->Fill(momentum);
1445  p_eff_all_vs_mom->Fill(momentum, 0.0);
1446  p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
1447  p_eff_all_vs_pt->Fill(pt, 0.0);
1448  if (mtra.mother_ID < 0) { // primary
1449  p_eff_prim_vs_mom->Fill(momentum, 0.0);
1450  p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
1451  p_eff_prim_vs_pt->Fill(pt, 0.0);
1452  p_eff_prim_vs_theta->Fill(theta, 0.0);
1453  } else {
1454  p_eff_sec_vs_mom->Fill(momentum, 0.0);
1455  p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
1456  }
1457  if (mtra.mother_ID < 0) { // primary
1458  h_rest_mom_prim->Fill(momentum);
1459  h_rest_nhits_prim->Fill(nSta);
1460  h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
1461  h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
1462  if (lh.iStation >= fh.iStation)
1463  h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation + 1,
1464  lh.iStation + 1);
1465  } else {
1466  h_rest_mom_sec->Fill(momentum);
1467  h_rest_nhits_sec->Fill(nSta);
1468  if (momentum > 0.01) {
1469  h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
1470  h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1471  if (lh.iStation >= fh.iStation)
1472  h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation + 1,
1473  lh.iStation + 1);
1474  }
1475  }
1476  }
1477 
1478  // killed tracks. At least one hit of they belong to some recoTrack
1479  // bool killed = 0;
1480  if (!reco) {
1481  h_notfound_nhits->Fill(nmchits);
1482  h_notfound_station->Fill(ph.iStation);
1483  h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1484 
1485  if (mtra.Points.size() > 0) {
1486  CbmL1MCPoint& pMC = vMCPoints[mtra.Points[0]];
1487  h_notfound_tx->Fill(pMC.px / pMC.pz);
1488  h_notfound_ty->Fill(pMC.py / pMC.pz);
1489  }
1490 
1491  // CbmL1HitStore &ph21 = vHitStore[mtra.StsHits[0]];
1492  // CbmL1HitStore &ph22 = vHitStore[mtra.StsHits[1]];
1493 
1494  // double z21 = algo->vStations[ph21.iStation].z[0];
1495  // double z22 = algo->vStations[ph22.iStation].z[0];
1496  // if( fabs(z22-z21)>1.e-4 ){
1497  // h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
1498  // h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
1499  // }
1500 
1501  // if( mtra.IsDisturbed() ) killed = 1;
1502  }
1503 
1504 
1505  if ((mtra.IsPrimary()) && (mtra.z > 0)) { // D0
1506  h_reco_d0_mom->Fill(momentum);
1507  if (reco)
1508  p_eff_d0_vs_mom->Fill(momentum, 100.0);
1509  else
1510  p_eff_d0_vs_mom->Fill(momentum, 0.0);
1511  }
1512 
1513  } // for mcTracks
1514 
1515  int NFakes = 0;
1516  for (unsigned int ih = 0; ih < algo->vStsHits->size(); ih++) {
1517  int iMC = vHitMCRef[ih]; // TODO2: adapt to linking
1518  if (iMC < 0) NFakes++;
1519  }
1520 
1521  h_reco_time->Fill(algo->CATime);
1522  h_reco_timeNtr->Fill(mc_total, algo->CATime);
1523  h_reco_timeNhit->Fill(algo->vStsHits->size(), algo->CATime);
1524 
1525  h_reco_fakeNtr->Fill(mc_total, NFakes);
1526  h_reco_fakeNhit->Fill(algo->vStsHits->size() - NFakes, NFakes);
1527 
1528 
1529  h_reg_MCmom->Scale(1.f / NEvents);
1530  h_acc_MCmom->Scale(1.f / NEvents);
1531  h_reco_MCmom->Scale(1.f / NEvents);
1532  h_ghost_Rmom->Scale(1.f / NEvents);
1533  h_reg_prim_MCmom->Scale(1.f / NEvents);
1534  h_acc_prim_MCmom->Scale(1.f / NEvents);
1535  h_reco_prim_MCmom->Scale(1.f / NEvents);
1536  h_reg_sec_MCmom->Scale(1.f / NEvents);
1537  h_acc_sec_MCmom->Scale(1.f / NEvents);
1538  h_reco_sec_MCmom->Scale(1.f / NEvents);
1539 
1540 } // void CbmL1::HistoPerformance()
1541 
1542 
1543 void CbmL1::TrackFitPerformance() {
1544  const int Nh_fit = 14;
1545  static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit],
1546  *h_fitPV[Nh_fit], *h_fit_chi2; //, *h_smoothed[12][Nh_fit];
1547 
1548  static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
1549 
1550  static bool first_call = 1;
1551 
1552  if (first_call) {
1553  first_call = 0;
1554 
1555  TDirectory* currentDir = gDirectory;
1556  gDirectory = histodir;
1557  gDirectory->cd("L1");
1558  gDirectory->mkdir("Fit");
1559  gDirectory->cd("Fit");
1560  {
1561  PRes2D =
1562  new TH2F("PRes2D", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1563  PRes2DPrim = new TH2F(
1564  "PRes2DPrim", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1565  PRes2DSec = new TH2F(
1566  "PRes2DSec", "Resolution P vs P", 100, 0., 20., 100, -15., 15.);
1567 
1568  struct {
1569  const char* name;
1570  const char* title;
1571  Int_t n;
1572  Double_t l, r;
1573  } Table[Nh_fit] = {
1574  {"x", "Residual X [#mum]", 140, -40., 40.},
1575  {"y", "Residual Y [#mum]", 100, -450., 450.},
1576  //{"y", "Residual Y [#mum]", 100, -45., 45.},
1577  {"tx", "Residual Tx [mrad]", 100, -4., 4.},
1578  //{"tx", "Residual Tx [mrad]", 100, -7., 7.},
1579  //{"tx", "Residual Tx [mrad]", 100, -2.5, 2.5},
1580  {"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
1581  //{"ty", "Residual Ty [mrad]", 100, -2.5, 2.5},
1582  {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1},
1583  //{"P", "Resolution P/Q [100%]", 100, -0.2, 0.2 },
1584  {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1585  {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1586  {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1587  {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1588  {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1589  {"QPreco", "Reco Q/P ", 100, -10., 10.},
1590  {"QPmc", "MC Q/P ", 100, -10., 10.},
1591  {"t", "Residual T [ns]", 100, -6., 6.},
1592  {"pt", "Pull T [residual/estimated_error]", 100, -6., 6.}};
1593 
1594  struct Tab {
1595  const char* name;
1596  const char* title;
1597  Int_t n;
1598  Double_t l, r;
1599  };
1600  Tab TableVertex[Nh_fit] = {
1601  //{"x", "Residual X [cm]", 200, -0.01, 0.01},
1602  {"x", "Residual X [cm]", 2000, -1., 1.},
1603  //{"y", "Residual Y [cm]", 200, -0.01, 0.01},
1604  {"y", "Residual Y [cm]", 2000, -1., 1.},
1605  //{"tx", "Residual Tx [mrad]", 100, -3., 3.},
1606  {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1607  //{"ty", "Residual Ty [mrad]", 100, -3., 3.},
1608  {"ty", "Residual Ty [mrad]", 100, -2., 2.},
1609  {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1},
1610  {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1611  {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1612  {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1613  {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1614  {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1615  {"QPreco", "Reco Q/P ", 100, -10., 10.},
1616  {"QPmc", "MC Q/P ", 100, -10., 10.},
1617  {"t", "Residual T [ns]", 100, -10., 10.},
1618  {"pt", "Pull T [residual/estimated_error]", 100, -6., 6.}};
1619 
1620  for (int i = 0; i < Nh_fit; i++) {
1621  char n[225], t[255];
1622  sprintf(n, "fst_%s", Table[i].name);
1623  sprintf(t, "First point %s", Table[i].title);
1624  h_fit[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1625  sprintf(n, "lst_%s", Table[i].name);
1626  sprintf(t, "Last point %s", Table[i].title);
1627  h_fitL[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1628  sprintf(n, "svrt_%s", TableVertex[i].name);
1629  sprintf(t, "Secondary vertex point %s", TableVertex[i].title);
1630  h_fitSV[i] =
1631  new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1632  sprintf(n, "pvrt_%s", TableVertex[i].name);
1633  sprintf(t, "Primary vertex point %s", TableVertex[i].title);
1634  h_fitPV[i] =
1635  new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1636 
1637  for (int j = 0; j < 12; j++) {
1638  sprintf(n, "smth_%s_sta_%i", Table[i].name, j);
1639  sprintf(t, "Station %i %s", j, Table[i].title);
1640  // h_smoothed[j][i] = new TH1F(n,t, Table[i].n, Table[i].l, Table[i].r);
1641  }
1642  }
1643  h_fit_chi2 = new TH1F("h_fit_chi2", "Chi2/NDF", 50, -0.5, 10.0);
1644  }
1645  gDirectory = currentDir;
1646  } // if first call
1647 
1648  for (vector<CbmL1Track>::iterator it = vRTracks.begin(); it != vRTracks.end();
1649  ++it) {
1650 
1651  if (it->IsGhost()) continue;
1652 
1653  { // first hit
1654 #define L1FSTPARAMEXTRAPOLATE
1655 #ifdef L1FSTPARAMEXTRAPOLATE
1656 
1657  const int last_station = vHitStore[it->StsHits.back()].iStation;
1658 
1659  CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
1660  L1TrackPar trPar(it->T, it->C);
1661  L1FieldRegion fld _fvecalignment;
1662  L1FieldValue B[3], targB _fvecalignment;
1663  float z[3] = {0.f, 0.f, 0.f};
1664  int ih = 0;
1665  for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
1666  const int iMCP = mc.Points[iMCPoint];
1667  CbmL1MCPoint& mcP = vMCPoints[iMCP];
1668  L1Station& st = algo->vStations[mcP.iStation];
1669  z[ih] = st.z[0];
1670  if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
1671  st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
1672  ih++;
1673  if (ih == 3) break;
1674  }
1675  if (ih < 3) continue;
1676  CbmL1MCPoint& mcP = vMCPoints[mc.Points[0]];
1677  targB = algo->vtxFieldValue;
1678  fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1679  L1Extrapolate(trPar, mcP.zIn, trPar.qp, fld);
1680 
1681  h_fit[0]->Fill((trPar.x[0] - mcP.xIn) * 1.e4);
1682  h_fit[1]->Fill((trPar.y[0] - mcP.yIn) * 1.e4);
1683  h_fit[2]->Fill((trPar.tx[0] - mcP.pxIn / mcP.pzIn) * 1.e3);
1684  h_fit[3]->Fill((trPar.ty[0] - mcP.pyIn / mcP.pzIn) * 1.e3);
1685  h_fit[4]->Fill(fabs(1. / trPar.qp[0]) / mcP.p - 1);
1686 
1687  PRes2D->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
1688 
1689  CbmL1MCTrack mcTrack = *(it->GetMCTracks()[0]);
1690  if (mcTrack.IsPrimary())
1691  PRes2DPrim->Fill(mcP.p,
1692  (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
1693  else
1694  PRes2DSec->Fill(mcP.p, (1. / fabs(trPar.qp[0]) - mcP.p) / mcP.p * 100.);
1695 
1696  if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1697  h_fit[5]->Fill((trPar.x[0] - mcP.xIn) / sqrt(trPar.C00[0]));
1698  if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1699  h_fit[6]->Fill((trPar.y[0] - mcP.yIn) / sqrt(trPar.C11[0]));
1700  if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1701  h_fit[7]->Fill((trPar.tx[0] - mcP.pxIn / mcP.pzIn)
1702  / sqrt(trPar.C22[0]));
1703  if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1704  h_fit[8]->Fill((trPar.ty[0] - mcP.pyIn / mcP.pzIn)
1705  / sqrt(trPar.C33[0]));
1706  if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1707  h_fit[9]->Fill((trPar.qp[0] - mcP.q / mcP.p) / sqrt(trPar.C44[0]));
1708  h_fit[10]->Fill(trPar.qp[0]);
1709  h_fit[11]->Fill(mcP.q / mcP.p);
1710  if (last_station > NMvdStations) h_fit[12]->Fill(trPar.t[0] - mcP.time);
1711  if (last_station > NMvdStations)
1712  if (finite(trPar.C55[0]) && trPar.C55[0] > 0)
1713  h_fit[13]->Fill((trPar.t[0] - mcP.time) / sqrt(trPar.C55[0]));
1714 
1715 #else
1716  int iMC = vHitMCRef[it->StsHits.front()]; // TODO2: adapt to linking
1717  if (iMC < 0) continue;
1718  CbmL1MCPoint& mc = vMCPoints[iMC];
1719  // if( !( mc.pdg == -11 && mc.mother_ID == -1 ) ) continue; // electrons only
1720  h_fit[0]->Fill((mc.x - it->T[0]) * 1.e4);
1721  h_fit[1]->Fill((mc.y - it->T[1]) * 1.e4);
1722  h_fit[2]->Fill((mc.px / mc.pz - it->T[2]) * 1.e3);
1723  h_fit[3]->Fill((mc.py / mc.pz - it->T[3]) * 1.e3);
1724  h_fit[4]->Fill(it->T[4] / mc.q * mc.p - 1);
1725 
1726  PRes2D->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1727 
1728  CbmL1MCTrack mcTrack = *(it->GetMCTracks()[0]);
1729  if (mcTrack.IsPrimary())
1730  PRes2DPrim->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1731  else
1732  PRes2DSec->Fill(mc.p, (1. / fabs(it->T[4]) - mc.p) / mc.p * 100.);
1733 
1734  if (finite(it->C[0]) && it->C[0] > 0)
1735  h_fit[5]->Fill((mc.x - it->T[0]) / sqrt(it->C[0]));
1736  if (finite(it->C[2]) && it->C[2] > 0)
1737  h_fit[6]->Fill((mc.y - it->T[1]) / sqrt(it->C[2]));
1738  if (finite(it->C[5]) && it->C[5] > 0)
1739  h_fit[7]->Fill((mc.px / mc.pz - it->T[2]) / sqrt(it->C[5]));
1740  if (finite(it->C[9]) && it->C[9] > 0)
1741  h_fit[8]->Fill((mc.py / mc.pz - it->T[3]) / sqrt(it->C[9]));
1742  if (finite(it->C[14]) && it->C[14] > 0)
1743  h_fit[9]->Fill((mc.q / mc.p - it->T[4]) / sqrt(it->C[14]));
1744  h_fit[10]->Fill(it->T[4]);
1745  h_fit[11]->Fill(mc.q / mc.p);
1746  h_fit[12]->Fill(mc.time - it->T[6]);
1747  if (finite(it->C[20]) && it->C[20] > 0)
1748  h_fit[13]->Fill((mc.time - it->T[6]) / sqrt(it->C[20]));
1749 #endif
1750  }
1751 
1752 
1753  { // last hit
1754  int iMC = vHitMCRef[it->StsHits.back()]; // TODO2: adapt to linking
1755  if (iMC < 0) continue;
1756 
1757 #define L1FSTPARAMEXTRAPOLATE
1758 #ifdef L1FSTPARAMEXTRAPOLATE
1759 
1760  const int last_station = vHitStore[it->StsHits.back()].iStation;
1761 
1762  CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
1763  L1TrackPar trPar(it->TLast, it->CLast);
1764  L1FieldRegion fld _fvecalignment;
1765  L1FieldValue B[3], targB _fvecalignment;
1766  float z[3] = {0.f, 0.f, 0.f};
1767  int ih = 0;
1768  for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
1769  const int iMCP = mc.Points[iMCPoint];
1770  CbmL1MCPoint& mcP = vMCPoints[iMCP];
1771  L1Station& st = algo->vStations[mcP.iStation];
1772  z[ih] = st.z[0];
1773  if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
1774  st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
1775  ih++;
1776  if (ih == 3) break;
1777  }
1778  if (ih < 3) continue;
1779  CbmL1MCPoint& mcP = vMCPoints[iMC];
1780  targB = algo->vtxFieldValue;
1781  fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1782  L1Extrapolate(trPar, mcP.zOut, trPar.qp, fld);
1783 
1784  h_fitL[0]->Fill((trPar.x[0] - mcP.xOut) * 1.e4);
1785  h_fitL[1]->Fill((trPar.y[0] - mcP.yOut) * 1.e4);
1786  h_fitL[2]->Fill((trPar.tx[0] - mcP.pxOut / mcP.pzOut) * 1.e3);
1787  h_fitL[3]->Fill((trPar.ty[0] - mcP.pyOut / mcP.pzOut) * 1.e3);
1788  h_fitL[4]->Fill(fabs(1. / trPar.qp[0]) / mcP.p - 1);
1789  if (last_station > NMvdStations) h_fitL[12]->Fill(trPar.t[0] - mcP.time);
1790 
1791 
1792  if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1793  h_fitL[5]->Fill((trPar.x[0] - mcP.xOut) / sqrt(trPar.C00[0]));
1794  if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1795  h_fitL[6]->Fill((trPar.y[0] - mcP.yOut) / sqrt(trPar.C11[0]));
1796  if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1797  h_fitL[7]->Fill((trPar.tx[0] - mcP.pxOut / mcP.pzOut)
1798  / sqrt(trPar.C22[0]));
1799  if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1800  h_fitL[8]->Fill((trPar.ty[0] - mcP.pyOut / mcP.pzOut)
1801  / sqrt(trPar.C33[0]));
1802  if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1803  h_fitL[9]->Fill((trPar.qp[0] - mcP.q / mcP.p) / sqrt(trPar.C44[0]));
1804  h_fitL[10]->Fill(trPar.qp[0]);
1805  h_fitL[11]->Fill(mcP.q / mcP.p);
1806  if (last_station > NMvdStations)
1807  if (finite(trPar.C55[0]) && trPar.C55[0] > 0)
1808  h_fitL[13]->Fill((trPar.t[0] - mcP.time) / sqrt(trPar.C55[0]));
1809 #else
1810  CbmL1MCPoint& mc = vMCPoints[iMC];
1811 
1812  h_fitL[0]->Fill((it->TLast[0] - mc.x) * 1.e4);
1813  h_fitL[1]->Fill((it->TLast[1] - mc.y) * 1.e4);
1814  h_fitL[2]->Fill((it->TLast[2] - mc.px / mc.pz) * 1.e3);
1815  h_fitL[3]->Fill((it->TLast[3] - mc.py / mc.pz) * 1.e3);
1816  h_fitL[4]->Fill(fabs(1. / it->TLast[4]) / mc.p - 1);
1817  if (finite(it->CLast[0]) && it->CLast[0] > 0)
1818  h_fitL[5]->Fill((it->TLast[0] - mc.x) / sqrt(it->CLast[0]));
1819  if (finite(it->CLast[2]) && it->CLast[2] > 0)
1820  h_fitL[6]->Fill((it->TLast[1] - mc.y) / sqrt(it->CLast[2]));
1821  if (finite(it->CLast[5]) && it->CLast[5] > 0)
1822  h_fitL[7]->Fill((it->TLast[2] - mc.px / mc.pz) / sqrt(it->CLast[5]));
1823  if (finite(it->CLast[9]) && it->CLast[9] > 0)
1824  h_fitL[8]->Fill((it->TLast[3] - mc.py / mc.pz) / sqrt(it->CLast[9]));
1825  if (finite(it->CLast[14]) && it->CLast[14] > 0)
1826  h_fitL[9]->Fill((it->TLast[4] - mc.q / mc.p) / sqrt(it->CLast[14]));
1827  h_fitL[10]->Fill(it->TLast[4]);
1828  h_fitL[11]->Fill(mc.q / mc.p);
1829  h_fitL[12]->Fill(it->TLast[6] - mc.time);
1830  if (finite(it->CLast[20]) && it->CLast[20] > 0)
1831  h_fitL[13]->Fill((it->TLast[6] - mc.time) / sqrt(it->CLast[20]));
1832 
1833 #endif
1834  }
1835 
1836  { // vertex
1837  CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
1838  L1TrackPar trPar(it->T, it->C);
1839 
1840  if (mc.mother_ID != -1) { // secondary
1841 
1842  { // extrapolate to vertex
1843  L1FieldRegion fld _fvecalignment;
1844  L1FieldValue B[3] _fvecalignment;
1845  float z[3];
1846  for (unsigned int ih = 0; ih < 3; ih++) {
1847  if (ih >= mc.Points.size()) continue; //If nofMCPoints in track < 3
1848  const int iMCP = mc.Points[ih];
1849  CbmL1MCPoint& mcP = vMCPoints[iMCP];
1850  L1Station& st = algo->vStations[mcP.iStation];
1851  z[ih] = st.z[0];
1852  st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
1853  };
1854  fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1855 
1856  L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1857  // add material
1858  const int fSta = vHitStore[it->StsHits[0]].iStation;
1859  const int dir = int((mc.z - algo->vStations[fSta].z[0])
1860  / fabs(mc.z - algo->vStations[fSta].z[0]));
1861  // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
1862  for (int iSta = fSta /*+dir*/;
1863  (iSta >= 0) && (iSta < NStation)
1864  && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0);
1865  iSta += dir) {
1866  // cout << iSta << " " << dir << endl;
1867  L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp);
1868  if (iSta + dir == NMvdStations - 1)
1869  L1AddPipeMaterial(trPar, trPar.qp);
1870  }
1871  }
1872  if (mc.z != trPar.z[0]) continue;
1873  // static int good = 0;
1874  // static int bad = 0;
1875  // if (mc.z != trPar.z[0]){
1876  // bad++;
1877  // continue;
1878  // }
1879  // else good++;
1880  // cout << "bad\\good" << bad << " " << good << endl;
1881 
1882  // calculate pulls
1883  //h_fitSV[0]->Fill( (mc.x-trPar.x[0]) *1.e4);
1884  //h_fitSV[1]->Fill( (mc.y-trPar.y[0]) *1.e4);
1885  h_fitSV[0]->Fill((trPar.x[0] - mc.x));
1886  h_fitSV[1]->Fill((trPar.y[0] - mc.y));
1887  h_fitSV[2]->Fill((trPar.tx[0] - mc.px / mc.pz) * 1.e3);
1888  h_fitSV[3]->Fill((trPar.ty[0] - mc.py / mc.pz) * 1.e3);
1889  h_fitSV[4]->Fill(fabs(1. / trPar.qp[0]) / mc.p - 1);
1890  if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1891  h_fitSV[5]->Fill((trPar.x[0] - mc.x) / sqrt(trPar.C00[0]));
1892  if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1893  h_fitSV[6]->Fill((trPar.y[0] - mc.y) / sqrt(trPar.C11[0]));
1894  if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1895  h_fitSV[7]->Fill((trPar.tx[0] - mc.px / mc.pz) / sqrt(trPar.C22[0]));
1896  if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1897  h_fitSV[8]->Fill((trPar.ty[0] - mc.py / mc.pz) / sqrt(trPar.C33[0]));
1898  if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1899  h_fitSV[9]->Fill((trPar.qp[0] - mc.q / mc.p) / sqrt(trPar.C44[0]));
1900  h_fitSV[10]->Fill(trPar.qp[0]);
1901  h_fitSV[11]->Fill(mc.q / mc.p);
1902  h_fitSV[12]->Fill(trPar.t[0] - mc.time);
1903  if (finite(trPar.C55[0]) && trPar.C55[0] > 0)
1904  h_fitSV[13]->Fill((trPar.t[0] - mc.time) / sqrt(trPar.C55[0]));
1905  } else { // primary
1906 
1907 #define L1EXTRAPOLATE
1908 #ifdef L1EXTRAPOLATE
1909  { // extrapolate to vertex
1910  L1FieldRegion fld _fvecalignment;
1911  L1FieldValue B[3], targB _fvecalignment;
1912  float z[3];
1913 
1914  targB = algo->vtxFieldValue;
1915 
1916  int ih = 1;
1917  for (unsigned int iHit = 0; iHit < it->StsHits.size(); iHit++) {
1918  const int iStation = vHitStore[it->StsHits[iHit]].iStation;
1919  L1Station& st = algo->vStations[iStation];
1920  z[ih] = st.z[0];
1921  st.fieldSlice.GetFieldValue(vHitStore[it->StsHits[iHit]].x,
1922  vHitStore[it->StsHits[iHit]].y,
1923  B[ih]);
1924  ih++;
1925  if (ih == 3) break;
1926  }
1927  if (ih < 3) continue;
1928 
1929  // add material
1930  const int fSta = vHitStore[it->StsHits[0]].iStation;
1931 
1932  const int dir = (mc.z - algo->vStations[fSta].z[0])
1933  / abs(mc.z - algo->vStations[fSta].z[0]);
1934  // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
1935 
1936  for (int iSta = fSta + dir;
1937  (iSta >= 0) && (iSta < NStation)
1938  && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0);
1939  iSta += dir) {
1940 
1941  z[0] = algo->vStations[iSta].z[0];
1942  float dz = z[1] - z[0];
1943  algo->vStations[iSta].fieldSlice.GetFieldValue(
1944  trPar.x - trPar.tx * dz, trPar.y - trPar.ty * dz, B[0]);
1945  fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1946 
1947  fvec mass2 = 0.1395679f * 0.1395679f;
1948  L1Extrapolate(trPar, algo->vStations[iSta].z[0], trPar.qp, fld);
1949  L1AddMaterial(trPar,
1950  algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y),
1951  trPar.qp);
1952  EnergyLossCorrection(
1953  trPar,
1954  mass2,
1955  algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y),
1956  trPar.qp,
1957  fvec(1.f));
1958  if (iSta + dir == NMvdStations - 1) {
1959  L1AddPipeMaterial(trPar, trPar.qp);
1960  EnergyLossCorrection(
1961  trPar, mass2, PipeRadThick, trPar.qp, fvec(1.f));
1962  }
1963  B[2] = B[1];
1964  z[2] = z[1];
1965  B[1] = B[0];
1966  z[1] = z[0];
1967  }
1968 
1969  z[0] = 0;
1970  B[0] = targB;
1971  fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
1972  L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1973  }
1974  if (mc.z != trPar.z[0]) continue;
1975 
1976 
1977  // calculate pulls
1978  h_fitPV[0]->Fill((mc.x - trPar.x[0]));
1979  h_fitPV[1]->Fill((mc.y - trPar.y[0]));
1980  h_fitPV[2]->Fill((mc.px / mc.pz - trPar.tx[0]) * 1.e3);
1981  h_fitPV[3]->Fill((mc.py / mc.pz - trPar.ty[0]) * 1.e3);
1982  h_fitPV[4]->Fill(fabs(1 / trPar.qp[0]) / mc.p - 1);
1983  if (finite(trPar.C00[0]) && trPar.C00[0] > 0)
1984  h_fitPV[5]->Fill((mc.x - trPar.x[0]) / sqrt(trPar.C00[0]));
1985  if (finite(trPar.C11[0]) && trPar.C11[0] > 0)
1986  h_fitPV[6]->Fill((mc.y - trPar.y[0]) / sqrt(trPar.C11[0]));
1987  if (finite(trPar.C22[0]) && trPar.C22[0] > 0)
1988  h_fitPV[7]->Fill((mc.px / mc.pz - trPar.tx[0]) / sqrt(trPar.C22[0]));
1989  if (finite(trPar.C33[0]) && trPar.C33[0] > 0)
1990  h_fitPV[8]->Fill((mc.py / mc.pz - trPar.ty[0]) / sqrt(trPar.C33[0]));
1991  if (finite(trPar.C44[0]) && trPar.C44[0] > 0)
1992  h_fitPV[9]->Fill((mc.q / mc.p - trPar.qp[0]) / sqrt(trPar.C44[0]));
1993  h_fitPV[10]->Fill(trPar.qp[0]);
1994  h_fitPV[11]->Fill(mc.q / mc.p);
1995  h_fitPV[12]->Fill(mc.time - trPar.t[0]);
1996  if (finite(trPar.C55[0]) && trPar.C55[0] > 0)
1997  h_fitPV[13]->Fill((mc.time - trPar.t[0]) / sqrt(trPar.C55[0]));
1998 #else
1999  FairTrackParam fTP;
2000 
2001  CbmKFMath::CopyTC2TrackParam(&fTP, it->T, it->C);
2002 
2003  CbmKFTrack kfTr;
2004  kfTr.SetTrackParam(fTP);
2005 
2006  kfTr.Extrapolate(mc.z);
2007  CbmL1Track it2;
2008  for (int ipar = 0; ipar < 6; ipar++)
2009  it2.T[ipar] = kfTr.GetTrack()[ipar];
2010  for (int ipar = 0; ipar < 15; ipar++)
2011  it2.C[ipar] = kfTr.GetCovMatrix()[ipar];
2012 
2013 
2014  // calculate pulls
2015  // h_fitPV[0]->Fill( (mc.x-it2.T[0]) *1.e4);
2016  // h_fitPV[1]->Fill( (mc.y-it2.T[1]) *1.e4);
2017  h_fitPV[0]->Fill((mc.x - it2.T[0]));
2018  h_fitPV[1]->Fill((mc.y - it2.T[1]));
2019  h_fitPV[2]->Fill((mc.px / mc.pz - it2.T[2]) * 1.e3);
2020  h_fitPV[3]->Fill((mc.py / mc.pz - it2.T[3]) * 1.e3);
2021  h_fitPV[4]->Fill(it2.T[4] / mc.q * mc.p - 1);
2022  if (finite(it2.C[0]) && it2.C[0] > 0)
2023  h_fitPV[5]->Fill((mc.x - it2.T[0]) / sqrt(it2.C[0]));
2024  if (finite(it2.C[2]) && it2.C[2] > 0)
2025  h_fitPV[6]->Fill((mc.y - it2.T[1]) / sqrt(it2.C[2]));
2026  if (finite(it2.C[5]) && it2.C[5] > 0)
2027  h_fitPV[7]->Fill((mc.px / mc.pz - it2.T[2]) / sqrt(it2.C[5]));
2028  if (finite(it2.C[9]) && it2.C[9] > 0)
2029  h_fitPV[8]->Fill((mc.py / mc.pz - it2.T[3]) / sqrt(it2.C[9]));
2030  if (finite(it2.C[14]) && it2.C[14] > 0)
2031  h_fitPV[9]->Fill((mc.q / mc.p - it2.T[4]) / sqrt(it2.C[14]));
2032  h_fitPV[10]->Fill(it2.T[4]);
2033  h_fitPV[11]->Fill(mc.q / mc.p);
2034  h_fitPV[12]->Fill(mc.time - it2.T[6]);
2035  if (finite(it2.C[20]) && it2.C[20] > 0)
2036  h_fitPV[13]->Fill((mc.time - it2.T[6]) / sqrt(it2.C[20]));
2037 #endif // L1EXTRAPOLATE
2038  }
2039  }
2040 
2041  h_fit_chi2->Fill(it->chi2 / it->NDF);
2042  }
2043 
2044 } // void CbmL1::TrackFitPerformance()
2045 
2046 
2047 void CbmL1::FieldApproxCheck() {
2048  TDirectory* curr = gDirectory;
2049  TFile* currentFile = gFile;
2050  TFile* fout = new TFile("FieldApprox.root", "RECREATE");
2051  fout->cd();
2052 
2053  FairField* MF = CbmKF::Instance()->GetMagneticField();
2054  for (int ist = 0; ist < NStation; ist++) {
2055  double z = 0;
2056  double Xmax = -100, Ymax = -100;
2057  if (ist < NMvdStations) {
2058  CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
2059  z = t.z;
2060  Xmax = Ymax = t.R;
2061  } else {
2062  CbmStsStation* station =
2063  CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
2064  z = station->GetZ();
2065 
2066  Xmax = station->GetXmax();
2067  Ymax = station->GetYmax();
2068  } // if mvd
2069 
2070 
2071  // float step = 1.;
2072 
2073  int NbinsX = 100; //static_cast<int>(2*Xmax/step);
2074  int NbinsY = 100; //static_cast<int>(2*Ymax/step);
2075  float ddx = 2 * Xmax / NbinsX;
2076  float ddy = 2 * Ymax / NbinsY;
2077 
2078  TH2F* stB = new TH2F(Form("station %i, dB", ist + 1),
2079  Form("station %i, dB, z = %0.f cm", ist + 1, z),
2080  static_cast<int>(NbinsX + 1),
2081  -(Xmax + ddx / 2.),
2082  (Xmax + ddx / 2.),
2083  static_cast<int>(NbinsY + 1),
2084  -(Ymax + ddy / 2.),
2085  (Ymax + ddy / 2.));
2086  TH2F* stBx = new TH2F(Form("station %i, dBx", ist + 1),
2087  Form("station %i, dBx, z = %0.f cm", ist + 1, z),
2088  static_cast<int>(NbinsX + 1),
2089  -(Xmax + ddx / 2.),
2090  (Xmax + ddx / 2.),
2091  static_cast<int>(NbinsY + 1),
2092  -(Ymax + ddy / 2.),
2093  (Ymax + ddy / 2.));
2094  TH2F* stBy = new TH2F(Form("station %i, dBy", ist + 1),
2095  Form("station %i, dBy, z = %0.f cm", ist + 1, z),
2096  static_cast<int>(NbinsX + 1),
2097  -(Xmax + ddx / 2.),
2098  (Xmax + ddx / 2.),
2099  static_cast<int>(NbinsY + 1),
2100  -(Ymax + ddy / 2.),
2101  (Ymax + ddy / 2.));
2102  TH2F* stBz = new TH2F(Form("station %i, dBz", ist + 1),
2103  Form("station %i, dBz, z = %0.f cm", ist + 1, z),
2104  static_cast<int>(NbinsX + 1),
2105  -(Xmax + ddx / 2.),
2106  (Xmax + ddx / 2.),
2107  static_cast<int>(NbinsY + 1),
2108  -(Ymax + ddy / 2.),
2109  (Ymax + ddy / 2.));
2110 
2111  Double_t r[3], B[3];
2112  L1FieldSlice FSl;
2113  L1FieldValue B_L1;
2114  Double_t bbb, bbb_L1;
2115 
2116  const int M = 5; // polinom order
2117  const int N = (M + 1) * (M + 2) / 2;
2118  L1Station& st = algo->vStations[ist];
2119  for (int i = 0; i < N; i++) {
2120  FSl.cx[i] = st.fieldSlice.cx[i][0];
2121  FSl.cy[i] = st.fieldSlice.cy[i][0];
2122  FSl.cz[i] = st.fieldSlice.cz[i][0];
2123  }
2124 
2125  Int_t i = 1, j = 1;
2126 
2127  double x, y;
2128  for (int ii = 1; ii <= NbinsX + 1; ii++) {
2129  j = 1;
2130  x = -Xmax + (ii - 1) * ddx;
2131  for (int jj = 1; jj <= NbinsY + 1; jj++) {
2132  y = -Ymax + (jj - 1) * ddy;
2133  double rrr = sqrt(fabs(x * x / Xmax / Xmax + y / Ymax * y / Ymax));
2134  if (rrr > 1.) {
2135  j++;
2136  continue;
2137  }
2138  r[2] = z;
2139  r[0] = x;
2140  r[1] = y;
2141  MF->GetFieldValue(r, B);
2142  bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
2143 
2144  FSl.GetFieldValue(x, y, B_L1);
2145  bbb_L1 = sqrt(B_L1.x[0] * B_L1.x[0] + B_L1.y[0] * B_L1.y[0]
2146  + B_L1.z[0] * B_L1.z[0]);
2147 
2148  stB->SetBinContent(ii, jj, (bbb - bbb_L1));
2149  stBx->SetBinContent(ii, jj, (B[0] - B_L1.x[0]));
2150  stBy->SetBinContent(ii, jj, (B[1] - B_L1.y[0]));
2151  stBz->SetBinContent(ii, jj, (B[2] - B_L1.z[0]));
2152  j++;
2153  }
2154  i++;
2155  }
2156 
2157  stB->GetXaxis()->SetTitle("X, cm");
2158  stB->GetYaxis()->SetTitle("Y, cm");
2159  stB->GetXaxis()->SetTitleOffset(1);
2160  stB->GetYaxis()->SetTitleOffset(1);
2161  stB->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
2162  stB->GetZaxis()->SetTitleOffset(1.3);
2163 
2164  stBx->GetXaxis()->SetTitle("X, cm");
2165  stBx->GetYaxis()->SetTitle("Y, cm");
2166  stBx->GetXaxis()->SetTitleOffset(1);
2167  stBx->GetYaxis()->SetTitleOffset(1);
2168  stBx->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
2169  stBx->GetZaxis()->SetTitleOffset(1.3);
2170 
2171  stBy->GetXaxis()->SetTitle("X, cm");
2172  stBy->GetYaxis()->SetTitle("Y, cm");
2173  stBy->GetXaxis()->SetTitleOffset(1);
2174  stBy->GetYaxis()->SetTitleOffset(1);
2175  stBy->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
2176  stBy->GetZaxis()->SetTitleOffset(1.3);
2177 
2178  stBz->GetXaxis()->SetTitle("X, cm");
2179  stBz->GetYaxis()->SetTitle("Y, cm");
2180  stBz->GetXaxis()->SetTitleOffset(1);
2181  stBz->GetYaxis()->SetTitleOffset(1);
2182  stBz->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
2183  stBz->GetZaxis()->SetTitleOffset(1.3);
2184 
2185  stB->Write();
2186  stBx->Write();
2187  stBy->Write();
2188  stBz->Write();
2189 
2190  } // for ista
2191 
2192  fout->Close();
2193  fout->Delete();
2194  gFile = currentFile;
2195  gDirectory = curr;
2196 } // void CbmL1::FieldApproxCheck()
2197 
2198 #include "TMath.h"
2199 void CbmL1::FieldIntegralCheck() {
2200  TDirectory* curr = gDirectory;
2201  TFile* currentFile = gFile;
2202  TFile* fout = new TFile("FieldApprox.root", "RECREATE");
2203  fout->cd();
2204 
2205  FairField* MF = CbmKF::Instance()->GetMagneticField();
2206 
2207  int nPointsZ = 1000;
2208  int nPointsPhi = 100;
2209  int nPointsTheta = 100;
2210  double startZ = 0, endZ = 100.;
2211  double startPhi = 0, endPhi = 2 * TMath::Pi();
2212  double startTheta = -30. / 180. * TMath::Pi(),
2213  endTheta = 30. / 180. * TMath::Pi();
2214 
2215  double DZ = endZ - startZ;
2216  double DP = endPhi - startPhi;
2217  double DT = endTheta - startTheta;
2218 
2219  float ddp = endPhi / nPointsPhi;
2220  float ddt = 2 * endTheta / nPointsTheta;
2221 
2222  TH2F* hSb = new TH2F("Field Integral",
2223  "Field Integral",
2224  static_cast<int>(nPointsPhi),
2225  -(startPhi + ddp / 2.),
2226  (endPhi + ddp / 2.),
2227  static_cast<int>(nPointsTheta),
2228  (startTheta - ddt / 2.),
2229  (endTheta + ddt / 2.));
2230 
2231  for (int iP = 0; iP < nPointsPhi; iP++) {
2232  double phi = startPhi + iP * DP / nPointsPhi;
2233  for (int iT = 0; iT < nPointsTheta; iT++) {
2234  double theta = startTheta + iT * DT / nPointsTheta;
2235 
2236  double Sb = 0;
2237  for (int iZ = 1; iZ < nPointsZ; iZ++) {
2238  double z = startZ + DZ * iZ / nPointsZ;
2239  double x = z * TMath::Tan(theta) * TMath::Cos(phi);
2240  double y = z * TMath::Tan(theta) * TMath::Sin(phi);
2241  double r[3] = {x, y, z};
2242  double b[3];
2243  MF->GetFieldValue(r, b);
2244  double B = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
2245  Sb += B * DZ / nPointsZ / 100. / 10.;
2246  }
2247  hSb->SetBinContent(iP + 1, iT + 1, Sb);
2248  }
2249  }
2250 
2251  hSb->GetXaxis()->SetTitle("#phi [rad]");
2252  hSb->GetYaxis()->SetTitle("#theta [rad]");
2253  hSb->GetXaxis()->SetTitleOffset(1);
2254  hSb->GetYaxis()->SetTitleOffset(1);
2255  hSb->GetZaxis()->SetTitle("Field Integral [T#dotm]");
2256  hSb->GetZaxis()->SetTitleOffset(1.3);
2257 
2258  hSb->Write();
2259 
2260 
2261  fout->Close();
2262  fout->Delete();
2263  gFile = currentFile;
2264  gDirectory = curr;
2265 } // void CbmL1::FieldApproxCheck()
2266 
2267 void CbmL1::InputPerformance() {
2268  // static TH1I *nStripFHits, *nStripBHits, *nStripFMC, *nStripBMC;
2269 
2270  static TH1F *resXsts, *resYsts, *resTsts, *resXmvd, *resYmvd /*, *pullZ*/;
2271  static TH1F *pullXsts, *pullYsts, *pullTsts, *pullXmvd,
2272  *pullYmvd /*, *pullZ*/;
2273  static TH1F *pullXmuch, *pullYmuch, *pullTmuch, *resXmuch, *resYmuch,
2274  *resTmuch /*, *pullZ*/;
2275  static TH1F *pullXtrd, *pullYtrd, *pullTtrd, *resXtrd, *resYtrd,
2276  *resTtrd /*, *pullZ*/;
2277  static TH1F *pullXtof, *pullYtof, *pullTtof, *resXtof, *resYtof,
2278  *resTtof /*, *pullZ*/;
2279 
2280 
2281  static bool first_call = 1;
2282  if (first_call) {
2283  first_call = 0;
2284 
2285  TDirectory* currentDir = gDirectory;
2286  gDirectory = histodir;
2287  gDirectory->cd("L1");
2288  gDirectory->mkdir("Input");
2289  gDirectory->cd("Input");
2290  gDirectory->mkdir("STS");
2291  gDirectory->cd("STS");
2292 
2293  // nStripFHits = new TH1I("nHits_f", "NHits On Front Strip", 10, 0, 10);
2294  // nStripBHits = new TH1I("nHits_b", "NHits On Back Strip", 10, 0, 10);
2295  // nStripFMC = new TH1I("nMC_f", "N MC Points On Front Strip", 10, 0, 10);
2296  // nStripBMC = new TH1I("nMC_b", "N MC Points On Back Strip", 10, 0, 10);
2297 
2298  pullXsts = new TH1F("Px_sts", "STS Pull x", 100, -5, 5);
2299  pullYsts = new TH1F("Py_sts", "STS Pull y", 100, -5, 5);
2300  pullTsts = new TH1F("Pt_sts", "STS Pull t", 100, -5, 5);
2301 
2302  resXsts = new TH1F("x_sts", "STS Residual x", 100, -50, 50);
2303  resYsts = new TH1F("y_sts", "STS Residual y", 100, -500, 500);
2304  resTsts = new TH1F("t_sts", "STS Residual t", 100, -20, 20);
2305 
2306  gDirectory->cd("..");
2307  gDirectory->mkdir("MVD");
2308  gDirectory->cd("MVD");
2309 
2310  pullXmvd = new TH1F("Px_mvd", "MVD Pull x", 100, -5, 5);
2311  pullYmvd = new TH1F("Py_mvd", "MVD Pull y", 100, -5, 5);
2312 
2313  resXmvd = new TH1F("x_mvd", "MVD Residual x", 100, -50, 50);
2314  resYmvd = new TH1F("y_mvd", "MVD Residual y", 100, -50, 50);
2315 
2316  TH1* histo;
2317  histo = resXsts;
2318  histo->GetXaxis()->SetTitle("Residual x, um");
2319  histo = resYsts;
2320  histo->GetXaxis()->SetTitle("Residual y, um");
2321  histo = resTsts;
2322  histo->GetXaxis()->SetTitle("Residual t, ns");
2323  histo = resXmvd;
2324  histo->GetXaxis()->SetTitle("Residual x, um");
2325  histo = resYmvd;
2326  histo->GetXaxis()->SetTitle("Residual y, um");
2327  histo = pullXsts;
2328  histo->GetXaxis()->SetTitle("Pull x");
2329  histo = pullYsts;
2330  histo->GetXaxis()->SetTitle("Pull y");
2331  histo = pullTsts;
2332  histo->GetXaxis()->SetTitle("Pull t");
2333  histo = pullXmvd;
2334  histo->GetXaxis()->SetTitle("Pull x");
2335  histo = pullYmvd;
2336  histo->GetXaxis()->SetTitle("Pull y");
2337 
2338 
2339  gDirectory->cd("..");
2340  gDirectory->mkdir("MUCH");
2341  gDirectory->cd("MUCH");
2342 
2343  pullXmuch = new TH1F("Px_much", "MUCH Pull x", 100, -10, 10);
2344  pullYmuch = new TH1F("Py_much", "MUCH Pull y", 100, -10, 10);
2345  pullTmuch = new TH1F("Pt_much", "MUCH Pull t", 100, -10, 10);
2346 
2347  resXmuch = new TH1F("x_much", "MUCH Residual x", 100, -100000, 100000);
2348  resYmuch = new TH1F("y_much", "MUCH Residual y", 100, -100000, 100000);
2349  resTmuch = new TH1F("t_much", "MUCH Residual t", 100, -40, 40);
2350 
2351 
2352  histo = resXmuch;
2353  histo->GetXaxis()->SetTitle("Residual x, um");
2354  histo = resYmuch;
2355  histo->GetXaxis()->SetTitle("Residual y, um");
2356  histo = resTmuch;
2357  histo->GetXaxis()->SetTitle("Residual t, ns");
2358  histo = pullXmuch;
2359  histo->GetXaxis()->SetTitle("Pull x");
2360  histo = pullYmuch;
2361  histo->GetXaxis()->SetTitle("Pull y");
2362  histo = pullTmuch;
2363  histo->GetXaxis()->SetTitle("Pull t");
2364 
2365  gDirectory->cd("..");
2366  gDirectory->mkdir("TRD");
2367  gDirectory->cd("TRD");
2368 
2369  pullXtrd = new TH1F("Px_trd", "TRD Pull x", 100, -5, 5);
2370  pullYtrd = new TH1F("Py_trd", "TRD Pull y", 100, -5, 5);
2371  pullTtrd = new TH1F("Pt_trd", "TRD Pull t", 100, -5, 5);
2372 
2373  resXtrd = new TH1F("x_trd", "TRD Residual x", 100, -200000, 200000);
2374  resYtrd = new TH1F("y_trd", "TRD Residual y", 100, -200000, 200000);
2375  resTtrd = new TH1F("t_trd", "TRD Residual t", 100, -40, 40);
2376 
2377 
2378  histo = resXtrd;
2379  histo->GetXaxis()->SetTitle("Residual x, um");
2380  histo = resYtrd;
2381  histo->GetXaxis()->SetTitle("Residual y, um");
2382  histo = resTtrd;
2383  histo->GetXaxis()->SetTitle("Residual t, ns");
2384  histo = pullXtrd;
2385  histo->GetXaxis()->SetTitle("Pull x");
2386  histo = pullYtrd;
2387  histo->GetXaxis()->SetTitle("Pull y");
2388  histo = pullTtrd;
2389  histo->GetXaxis()->SetTitle("Pull t");
2390 
2391  gDirectory->cd("..");
2392  gDirectory->mkdir("TOF");
2393  gDirectory->cd("TOF");
2394 
2395  pullXtof = new TH1F("Px_tof", "TOF Pull x", 100, -5, 5);
2396  pullYtof = new TH1F("Py_tof", "TOF Pull y", 100, -5, 5);
2397  pullTtof = new TH1F("Pt_tof", "TOF Pull t", 100, -5, 5);
2398 
2399  resXtof = new TH1F("x_tof", "TOF Residual x", 100, -50000, 50000);
2400  resYtof = new TH1F("y_tof", "TOF Residual y", 100, -50000, 50000);
2401  resTtof = new TH1F("t_tof", "TOF Residual t", 100, -0.4, 0.4);
2402 
2403 
2404  histo = resXtof;
2405  histo->GetXaxis()->SetTitle("Residual x, um");
2406  histo = resYtof;
2407  histo->GetXaxis()->SetTitle("Residual y, um");
2408  histo = resTtof;
2409  histo->GetXaxis()->SetTitle("Residual t, ns");
2410  histo = pullXtof;
2411  histo->GetXaxis()->SetTitle("Pull x");
2412  histo = pullYtof;
2413  histo->GetXaxis()->SetTitle("Pull y");
2414  histo = pullTtof;
2415  histo->GetXaxis()->SetTitle("Pull t");
2416 
2417  gDirectory = currentDir;
2418  } // first_call
2419 
2420  // std::map<unsigned int, unsigned int> stripFToNHitMap,stripBToNHitMap;
2421  // std::map<unsigned int, unsigned int> stripFToNMCMap,stripBToNMCMap;
2422 
2423  map<unsigned int, unsigned int>::iterator it;
2424  Int_t nMC = -1;
2425  if (listStsPts) { nMC = listStsPts->GetEntries(); }
2426 
2427  if (listStsHits && listStsHitMatch) {
2428  for (unsigned int iH = 0; iH < vStsHits.size(); iH++) {
2429  const CbmL1StsHit& h = vStsHits[iH];
2430 
2431  if (h.Det != 1) continue; // mvd hit
2432  const CbmStsHit* sh =
2433  L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(h.extIndex));
2434 
2435 
2436  // int iMCPoint = -1;
2437  CbmLink link;
2438  CbmStsPoint* pt = 0;
2439 
2440  if (listStsClusterMatch) {
2441  const CbmMatch* frontClusterMatch = static_cast<const CbmMatch*>(
2442  listStsClusterMatch->At(sh->GetFrontClusterId()));
2443  const CbmMatch* backClusterMatch = static_cast<const CbmMatch*>(
2444  listStsClusterMatch->At(sh->GetBackClusterId()));
2445  CbmMatch stsHitMatch;
2446 
2447  for (Int_t iFrontLink = 0;
2448  iFrontLink < frontClusterMatch->GetNofLinks();
2449  iFrontLink++) {
2450  const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);
2451 
2452  for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks();
2453  iBackLink++) {
2454  const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
2455  if (frontLink == backLink) {
2456  stsHitMatch.AddLink(frontLink);
2457  stsHitMatch.AddLink(backLink);
2458  }
2459  }
2460  }
2461 
2462 
2463  if (stsHitMatch.GetNofLinks() > 0) {
2464  Float_t bestWeight = 0.f;
2465  for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
2466  if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
2467  bestWeight = stsHitMatch.GetLink(iLink).GetWeight();
2468  Int_t iFile = stsHitMatch.GetLink(iLink).GetFile();
2469  Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry();
2470 
2471  link = stsHitMatch.GetLink(iLink);
2472 
2473 
2474  pt = (CbmStsPoint*) fStsPoints->Get(
2475  iFile, iEvent, stsHitMatch.GetLink(iLink).GetIndex());
2476  }
2477  }
2478  }
2479 
2480  if (pt == 0) continue;
2481 
2482  double mcTime = pt->GetTime();
2483 
2484  if (fTimesliceMode)
2485  mcTime += fEventList->GetEventTime(link.GetEntry(), link.GetFile());
2486 
2487  // hit pulls and residuals
2488 
2489  TVector3 hitPos, mcPos, hitErr;
2490  sh->Position(hitPos);
2491  sh->PositionError(hitErr);
2492 
2493  // pt->Position(mcPos); // this is wrong!
2494  mcPos.SetX(pt->GetX(hitPos.Z()));
2495  mcPos.SetY(pt->GetY(hitPos.Z()));
2496  mcPos.SetZ(hitPos.Z());
2497 
2498 #if 0 // standard errors
2499  if (hitErr.X() != 0) pullXsts->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2500  if (hitErr.Y() != 0) pullYsts->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2501 #elif 1 // qa errors
2502  if (hitErr.X() != 0)
2503  pullXsts->Fill((hitPos.X() - mcPos.X()) / sh->GetDx()); // qa errors
2504  if (hitErr.Y() != 0)
2505  pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
2506 
2507  pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
2508 #else // errors used in TF
2509  if (hitErr.X() != 0)
2510  pullXsts->Fill((hitPos.X() - mcPos.X())
2511  / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
2512  if (hitErr.Y() != 0)
2513  pullYsts->Fill((hitPos.Y() - mcPos.Y())
2514  / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
2515 #endif
2516 
2517  resXsts->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
2518  resYsts->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
2519  resTsts->Fill((sh->GetTime() - mcTime));
2520  }
2521  }
2522  } // sts
2523 
2524 
2525  if (listMvdHits && listMvdHitMatches) {
2526  Int_t nEnt = listMvdHits->GetEntries();
2527  for (int j = 0; j < nEnt; j++) {
2528 
2529  CbmMvdHit* sh = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(j));
2530  CbmMatch* hm = L1_DYNAMIC_CAST<CbmMatch*>(listMvdHitMatches->At(j));
2531 
2532  int iMC = -1;
2533  // float mcWeight = -1.f;
2534  // for(int iDigiLink=0; iDigiLink<hm->GetNofLinks(); iDigiLink++)
2535  // {
2536  // if( hm->GetLink(iDigiLink).GetWeight() > mcWeight)
2537  // {
2538  // mcWeight = hm->GetLink(iDigiLink).GetWeight();
2539  // iMC = hm->GetLink(iDigiLink).GetIndex();
2540  // }
2541  // }
2542  if (hm->GetNofLinks() > 0) iMC = hm->GetLink(0).GetIndex();
2543 
2544 
2545  if (iMC < 0) continue;
2546  // hit pulls and residuals
2547 
2548 
2549  TVector3 hitPos, mcPos, hitErr;
2550  sh->Position(hitPos);
2551  sh->PositionError(hitErr);
2552 
2553  CbmMvdPoint* pt = 0;
2554  nMC = listMvdPts->GetEntriesFast();
2555 
2556  if (iMC >= 0 && iMC < nMC)
2557  pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(listMvdPts->At(iMC));
2558 
2559  if (!pt) {
2560  // cout << " No MC points! " << "iMC=" << iMC << endl;
2561  continue;
2562  }
2563 
2564  mcPos.SetX((pt->GetX() + pt->GetXOut()) / 2.);
2565  mcPos.SetY((pt->GetY() + pt->GetYOut()) / 2.);
2566  mcPos.SetZ(hitPos.Z());
2567 
2568  // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2569  // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2570  // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
2571  // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
2572  if (hitErr.X() != 0)
2573  pullXmvd->Fill(
2574  (hitPos.X() - mcPos.X())
2575  / sqrt(algo->vStations[0].XYInfo.C00[0])); // errors used in TF
2576  if (hitErr.Y() != 0)
2577  pullYmvd->Fill((hitPos.Y() - mcPos.Y())
2578  / sqrt(algo->vStations[0].XYInfo.C11[0]));
2579 
2580  resXmvd->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
2581  resYmvd->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
2582  }
2583  } // mvd
2584 
2585 
2586  if (fMuchPixelHits && listMuchHitMatches) {
2587  for (unsigned int iH = 0; iH < vStsHits.size(); iH++) {
2588  const CbmL1StsHit& h = vStsHits[iH];
2589 
2590  if (h.Det != 2) continue; // mvd hit
2591 
2592  const CbmMuchPixelHit* sh =
2593  L1_DYNAMIC_CAST<CbmMuchPixelHit*>(fMuchPixelHits->At(h.extIndex));
2594  CbmMatch* hm =
2595  L1_DYNAMIC_CAST<CbmMatch*>(listMuchHitMatches->At(h.extIndex));
2596 
2597 
2598  if (hm->GetNofLinks() == 0) continue;
2599  Float_t bestWeight = 0.f;
2600  Float_t totalWeight = 0.f;
2601  int iMCPoint = -1;
2602  CbmLink link;
2603 
2604  for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
2605  totalWeight += hm->GetLink(iLink).GetWeight();
2606  if (hm->GetLink(iLink).GetWeight() > bestWeight) {
2607  bestWeight = hm->GetLink(iLink).GetWeight();
2608  iMCPoint = hm->GetLink(iLink).GetIndex();
2609  link = hm->GetLink(iLink);
2610  }
2611  }
2612  if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
2613 
2614  CbmMuchPoint* pt = (CbmMuchPoint*) fMuchPoints->Get(
2615  link.GetFile(), link.GetEntry(), link.GetIndex());
2616  double mcTime = pt->GetTime();
2617 
2618  if (fTimesliceMode)
2619  mcTime += fEventList->GetEventTime(link.GetEntry(), link.GetFile());
2620  // mcTime+=20;
2621 
2622  // hit pulls and residuals
2623 
2624 
2625  TVector3 hitPos, mcPos, hitErr;
2626  sh->Position(hitPos);
2627  sh->PositionError(hitErr);
2628 
2629  // pt->Position(mcPos); // this is wrong!
2630  // mcPos.SetX( pt->GetX( hitPos.Z() ) );
2631  // mcPos.SetY( pt->GetY( hitPos.Z() ) );
2632  // mcPos.SetZ( hitPos.Z() );
2633 
2634  mcPos.SetX(0.5 * (pt->GetXIn() + pt->GetXOut()));
2635  mcPos.SetY(0.5 * (pt->GetYIn() + pt->GetYOut()));
2636  mcPos.SetZ(hitPos.Z());
2637 
2638 #if 0 // standard errors
2639  if (hitErr.X() != 0) pullXmuch->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2640  if (hitErr.Y() != 0) pullYmuch->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2641 #elif 1 // qa errors
2642  if (hitErr.X() != 0)
2643  pullXmuch->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
2644  if (hitErr.Y() != 0) pullYmuch->Fill((h.y - mcPos.Y()) / sh->GetDy());
2645 
2646  pullTmuch->Fill((h.t - mcTime) / sh->GetTimeError());
2647 #else // errors used in TF
2648  if (hitErr.X() != 0)
2649  pullXmuch->Fill((hitPos.X() - mcPos.X())
2650  / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
2651  if (hitErr.Y() != 0)
2652  pullYmuch->Fill((hitPos.Y() - mcPos.Y())
2653  / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
2654 #endif
2655 
2656  resXmuch->Fill((h.x - mcPos.X()) * 10 * 1000);
2657  resYmuch->Fill((h.y - mcPos.Y()) * 10 * 1000);
2658  resTmuch->Fill((h.t - mcTime));
2659  }
2660  } // much
2661 
2662 
2663  if (listTrdHits && fTrdHitMatches) {
2664  for (unsigned int iH = 0; iH < vStsHits.size(); iH++) {
2665  const CbmL1StsHit& h = vStsHits[iH];
2666 
2667  if (h.Det != 3) continue; // mvd hit
2668  const CbmTrdHit* sh =
2669  L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(h.extIndex));
2670  CbmMatch* hm = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(h.extIndex));
2671 
2672 
2673  if (hm->GetNofLinks() == 0) continue;
2674  Float_t bestWeight = 0.f;
2675  Float_t totalWeight = 0.f;
2676  int iMCPoint = -1;
2677  CbmLink link;
2678 
2679  for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
2680  totalWeight += hm->GetLink(iLink).GetWeight();
2681  if (hm->GetLink(iLink).GetWeight() > bestWeight) {
2682  bestWeight = hm->GetLink(iLink).GetWeight();
2683  iMCPoint = hm->GetLink(iLink).GetIndex();
2684  link = hm->GetLink(iLink);
2685  }
2686  }
2687  if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
2688 
2689  CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get(
2690  link.GetFile(), link.GetEntry(), link.GetIndex());
2691  double mcTime = pt->GetTime();
2692 
2693  if (fTimesliceMode)
2694  mcTime += fEventList->GetEventTime(link.GetEntry(), link.GetFile());
2695 
2696  // hit pulls and residuals
2697  // if ((sh->GetPlaneId()) == 0) continue;
2698  // if ((sh->GetPlaneId()) == 2) continue;
2699  // if ((sh->GetPlaneId()) == 4) continue;
2700 
2701  TVector3 hitPos, mcPos, hitErr;
2702  sh->Position(hitPos);
2703  sh->PositionError(hitErr);
2704 
2705  // pt->Position(mcPos); // this is wrong!
2706  mcPos.SetX((pt->GetXIn() + pt->GetXOut()) / 2.);
2707  mcPos.SetY((pt->GetYIn() + pt->GetYOut()) / 2.);
2708  mcPos.SetZ(hitPos.Z());
2709 
2710 #if 0 // standard errors
2711  if (hitErr.X() != 0) pullXtrd->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2712  if (hitErr.Y() != 0) pullYtrd->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2713 #elif 1 // qa errors
2714  if (hitErr.X() != 0)
2715  pullXtrd->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
2716  if (hitErr.Y() != 0) pullYtrd->Fill((h.y - mcPos.Y()) / sh->GetDy());
2717 
2718  pullTtrd->Fill((h.t - mcTime) / sh->GetTimeError());
2719 #else // errors used in TF
2720  if (hitErr.X() != 0)
2721  pullXtrd->Fill((hitPos.X() - mcPos.X())
2722  / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
2723  if (hitErr.Y() != 0)
2724  pullYtrd->Fill((hitPos.Y() - mcPos.Y())
2725  / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
2726 #endif
2727 
2728  resXtrd->Fill((h.x - mcPos.X()) * 10 * 1000);
2729  resYtrd->Fill((h.y - mcPos.Y()) * 10 * 1000);
2730  resTtrd->Fill((h.t - mcTime));
2731  }
2732  } // much
2733 
2734 
2735  if (fTofHits && fTofHitDigiMatches) {
2736  for (unsigned int iH = 0; iH < vStsHits.size(); iH++) {
2737  const CbmL1StsHit& h = vStsHits[iH];
2738 
2739  if (h.Det != 4) continue; // mvd hit
2740 
2741  CbmTofHit* sh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(h.extIndex));
2742  CbmMatch* hm =
2743  L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(h.extIndex));
2744 
2745 
2746  if (hm->GetNofLinks() == 0) continue;
2747  Float_t bestWeight = 0.f;
2748  Float_t totalWeight = 0.f;
2749  int iMCPoint = -1;
2750  CbmLink link;
2751 
2752  for (int iLink = 0; iLink < hm->GetNofLinks(); iLink++) {
2753  totalWeight += hm->GetLink(iLink).GetWeight();
2754  if (hm->GetLink(iLink).GetWeight() > bestWeight) {
2755  bestWeight = hm->GetLink(iLink).GetWeight();
2756  iMCPoint = hm->GetLink(iLink).GetIndex();
2757  link = hm->GetLink(iLink);
2758  }
2759  }
2760 
2761  // if (bestWeight / totalWeight < 0.7 || iMCPoint < 0) continue;
2762 
2763  if (iMCPoint < 0) continue;
2764 
2765  CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(
2766  link.GetFile(), link.GetEntry(), link.GetIndex());
2767  double mcTime = pt->GetTime();
2768 
2769  if (fTimesliceMode)
2770  mcTime += fEventList->GetEventTime(link.GetEntry(), link.GetFile());
2771 
2772  // hit pulls and residuals
2773 
2774 
2775  TVector3 hitPos, mcPos, hitErr;
2776  sh->Position(hitPos);
2777  sh->PositionError(hitErr);
2778 
2779  // pt->Position(mcPos); // this is wrong!
2780  mcPos.SetX((pt->GetX()));
2781  mcPos.SetY((pt->GetY()));
2782  mcPos.SetZ(hitPos.Z());
2783 
2784 #if 0 // standard errors
2785  if (hitErr.X() != 0) pullXmuch->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2786  if (hitErr.Y() != 0) pullYmuch->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2787 #elif 1 // qa errors
2788  if (hitErr.X() != 0)
2789  pullXtof->Fill((h.x - mcPos.X()) / sh->GetDx()); // qa errors
2790  if (hitErr.Y() != 0) pullYtof->Fill((h.y - mcPos.Y()) / sh->GetDy());
2791 
2792  pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
2793 #else // errors used in TF
2794  if (hitErr.X() != 0)
2795  pullXtof->Fill((hitPos.X() - mcPos.X())
2796  / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
2797  if (hitErr.Y() != 0)
2798  pullYtof->Fill((hitPos.Y() - mcPos.Y())
2799  / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
2800 #endif
2801 
2802  resXtof->Fill((h.x - mcPos.X()) * 10 * 1000);
2803  resYtof->Fill((h.y - mcPos.Y()) * 10 * 1000);
2804  resTtof->Fill((sh->GetTime() - mcTime));
2805  }
2806  } // much
2807 
2808  // for (it = stripFToNHitMap.begin(); it != stripFToNHitMap.end(); it++){
2809  // nStripFHits->Fill(it->second);
2810  // }
2811  // for (it = stripBToNHitMap.begin(); it != stripBToNHitMap.end(); it++){
2812  // nStripBHits->Fill(it->second);
2813  // }
2814  // for (it = stripFToNMCMap.begin(); it != stripFToNMCMap.end(); it++){
2815  // nStripFMC->Fill(it->second);
2816  // }
2817  // for (it = stripBToNMCMap.begin(); it != stripBToNMCMap.end(); it++){
2818  // nStripBMC->Fill(it->second);
2819  // }
2820 
2821  // strips Not ended
2822  // if( listCbmStsDigi ){
2823  // Int_t nEnt = listCbmStsDigi->GetEntries();
2824  // for (int j=0; j < nEnt; j++ ){
2825  // CbmStsDigi *digi = (CbmStsDigi*) listCbmStsDigi->At(j);
2826  // // = sh->GetNLinks(0);
2827  // // find position of mc point
2828  // FairLink link = digi->GetLink(0);
2829  // int iMCPoint = link.GetIndex();
2830  // CbmStsPoint *mcPoint = (CbmStsPoint*) listStsPts->At(iMCPoint);
2831  // TVector3 mcPos;
2832  // if (digi->GetSide() == 0)
2833  // mcPoint->PositionIn(mcPos);
2834  // else
2835  // mcPoint->PositionOut(mcPos);
2836  //
2837  // CbmStsStation_old *sta = StsDigi.GetStation(digi->GetStationNr());
2838  // CbmStsSector* sector = sta->GetSector(digi->GetSectorNr());
2839  // digi->GetChannelNr();
2840  //
2841  // }
2842  // } // listCbmStsDigi
2843 
2844 
2845 } // void CbmL1::InputPerformance()
CbmL1::fTofHits
TClonesArray * fTofHits
Definition: CbmL1.h:314
CbmMuchPoint::GetXOut
Double_t GetXOut() const
Definition: CbmMuchPoint.h:73
CbmStsStation::GetZ
Double_t GetZ() const
Definition: CbmStsStation.h:127
L1Algo::vStsHits
const vector< L1StsHit > * vStsHits
Definition: L1Algo.h:344
CbmL1MCPoint::q
double q
Definition: CbmL1MCPoint.h:59
CbmL1::listTrdHits
TClonesArray * listTrdHits
Definition: CbmL1.h:308
CbmKFTube::z
Double_t z
Definition: CbmKFMaterial.h:92
CbmMatch
Definition: CbmMatch.h:22
L1TrackPar::qp
fvec qp
Definition: L1TrackPar.h:9
CbmL1TrackPar::NDF
int NDF
Definition: CbmL1TrackPar.h:18
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
CbmTrdPoint::GetYOut
Double_t GetYOut() const
Definition: CbmTrdPoint.h:67
L1TrackPar::t
fvec t
Definition: L1TrackPar.h:9
TL1TracksCatCounters::counters
vector< T > counters
Definition: CbmL1Counters.h:114
CbmL1Track::StsHits
vector< int > StsHits
Definition: CbmL1Track.h:67
CbmL1::InputPerformance
void InputPerformance()
Definition: CbmL1Performance.cxx:2267
CbmMuchPoint
Definition: CbmMuchPoint.h:21
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
L1Algo.h
CbmPixelHit::Position
void Position(TVector3 &pos) const
Copies hit position to pos.
Definition: CbmPixelHit.cxx:64
_fvecalignment
#define _fvecalignment
Definition: L1/vectors/P4_F32vec4.h:254
TL1PerfEfficiencies::killed
TL1TracksCatCounters< int > killed
Definition: CbmL1Performance.cxx:284
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
TL1PerfEfficiencies::ratio_killed
TL1TracksCatCounters< double > ratio_killed
Definition: CbmL1Performance.cxx:277
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
TL1PerfEfficiencies::ratio_clone
TL1TracksCatCounters< double > ratio_clone
Definition: CbmL1Performance.cxx:280
CbmL1::fEventList
CbmMCEventList * fEventList
Definition: CbmL1.h:263
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
CbmStsStation::GetYmax
Double_t GetYmax() const
Definition: CbmStsStation.h:119
CbmKF.h
CbmL1MCPoint::x
double x
Definition: CbmL1MCPoint.h:56
CbmL1::vStsHits
vector< CbmL1StsHit > vStsHits
Used data = Repacked input data.
Definition: CbmL1.h:335
CbmL1MCTrack::q
double q
Definition: CbmL1MCTrack.h:31
CbmL1::algo
L1Algo * algo
Definition: CbmL1.h:119
CbmL1MCPoint
Definition: CbmL1MCPoint.h:23
TL1PerfEfficiencies::mc_length
TL1TracksCatCounters< int > mc_length
Definition: CbmL1Performance.cxx:288
CbmL1MCPoint::time
double time
Definition: CbmL1MCPoint.h:59
CbmKF::GetMagneticField
FairField * GetMagneticField()
Definition: CbmKF.h:88
CbmL1HitStore::iStation
int iStation
Definition: CbmL1.h:97
CbmStsSetup.h
L1Station
Definition: L1Station.h:9
L1Station::fieldSlice
L1FieldSlice fieldSlice
Definition: L1Station.h:30
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmL1::fTrdPoints
CbmMCDataArray * fTrdPoints
Definition: CbmL1.h:307
CbmL1Constants::MinPurity
const double MinPurity
Definition: CbmL1Constants.h:10
L1FieldValue::z
fvec z
Definition: L1Field.h:17
CbmL1MCTrack::NMCStations
int NMCStations() const
Definition: CbmL1MCTrack.h:108
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, const L1FieldRegion &F, fvec *w=0)
Definition: L1Extrapolation.h:314
TL1PerfEfficiencies::reco_length
TL1TracksCatCounters< double > reco_length
Definition: CbmL1Performance.cxx:286
CbmL1Track::GetMCTrack
CbmL1MCTrack * GetMCTrack()
Definition: CbmL1Track.h:48
CbmL1MCPoint::y
double y
Definition: CbmL1MCPoint.h:56
CbmKFTrack::GetTrack
Double_t * GetTrack()
Is it electron.
Definition: CbmKFTrack.h:58
CbmL1::listStsHitMatch
TClonesArray * listStsHitMatch
Definition: CbmL1.h:277
CbmL1Constants::MinRefMom
const double MinRefMom
Definition: CbmL1Constants.h:8
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmL1HitStore::x
double x
Definition: CbmL1.h:98
CbmTrdHit
data class for a reconstructed Energy-4D measurement in the TRD
Definition: CbmTrdHit.h:35
TL1PerfEfficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition: CbmL1Performance.cxx:180
CbmL1::listMvdHits
TClonesArray * listMvdHits
Definition: CbmL1.h:288
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition: CbmStsSetup.cxx:293
CbmL1Track::hitMap
map< int, int > hitMap
Definition: CbmL1Track.h:74
CbmStsStation.h
CbmL1MCTrack::y
double y
Definition: CbmL1MCTrack.h:31
CbmL1MCTrack::x
double x
Definition: CbmL1MCTrack.h:31
acos
friend F32vec4 acos(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:138
CbmL1MCTrack::IsDisturbed
bool IsDisturbed() const
Definition: CbmL1MCTrack.h:120
CbmTrdPoint::GetXOut
Double_t GetXOut() const
Definition: CbmTrdPoint.h:66
CbmL1::listMvdHitMatches
TClonesArray * listMvdHitMatches
Definition: CbmL1.h:290
TL1Efficiencies::names
vector< TString > names
Definition: CbmL1Counters.h:145
L1TrackPar::z
fvec z
Definition: L1TrackPar.h:9
CbmL1HitStore
Definition: CbmL1.h:94
CbmStsHit::GetFrontClusterId
Int_t GetFrontClusterId() const
Definition: CbmStsHit.h:93
TL1PerfEfficiencies::clone
TL1TracksCatCounters< int > clone
Definition: CbmL1Performance.cxx:285
CbmL1.h
TL1PerfEfficiencies::ratio_fakes
TL1TracksCatCounters< double > ratio_fakes
Definition: CbmL1Performance.cxx:282
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:214
CbmKF::vMvdMaterial
std::vector< CbmKFTube > vMvdMaterial
Definition: CbmKF.h:73
L1FieldSlice::cz
fvec cz[21]
Definition: L1Field.h:34
CbmL1MCTrack::px
double px
Definition: CbmL1MCTrack.h:31
CbmStsPoint
Definition: CbmStsPoint.h:27
CbmL1::listStsPts
TClonesArray * listStsPts
Definition: CbmL1.h:272
CbmMatch.h
L1FieldSlice::cy
fvec cy[21]
Definition: L1Field.h:34
L1TrackPar::ty
fvec ty
Definition: L1TrackPar.h:9
TL1PerfEfficiencies::~TL1PerfEfficiencies
virtual ~TL1PerfEfficiencies()
Definition: CbmL1Performance.cxx:178
CbmKFTube::R
Double_t R
Definition: CbmKFMaterial.h:93
TL1Efficiencies::Inc
void Inc(bool isReco, TString name)
Definition: CbmL1Counters.h:191
L1AddMaterial
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.10565f *0.10565f)
Definition: L1AddMaterial.h:559
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
TL1Efficiencies::reco
TL1TracksCatCounters< int > reco
Definition: CbmL1Counters.h:154
CbmL1Constants.h
CbmHit::GetTimeError
Double_t GetTimeError() const
Definition: CbmHit.h:76
CbmL1::fTofHitDigiMatches
TClonesArray * fTofHitDigiMatches
Definition: CbmL1.h:313
CbmL1MCTrack::pdg
int pdg
Definition: CbmL1MCTrack.h:32
CbmL1::NMvdStations
int NMvdStations
Definition: CbmL1.h:244
CbmL1::listMvdPts
TClonesArray * listMvdPts
Definition: CbmL1.h:287
CbmMvdHit
Definition: CbmMvdHit.h:29
TL1PerfEfficiencies::operator+=
TL1PerfEfficiencies & operator+=(TL1PerfEfficiencies &a)
Definition: CbmL1Performance.cxx:194
CbmL1TrackPar::T
double T[7]
Definition: CbmL1TrackPar.h:17
CbmMuchPoint::GetYOut
Double_t GetYOut() const
Definition: CbmMuchPoint.h:74
L1TrackPar::y
fvec y
Definition: L1TrackPar.h:9
CbmL1MCPoint::pz
double pz
Definition: CbmL1MCPoint.h:56
TL1Efficiencies::operator+=
TL1Efficiencies & operator+=(TL1Efficiencies &a)
Definition: CbmL1Counters.h:181
CbmL1MCTrack::py
double py
Definition: CbmL1MCTrack.h:31
TL1TracksCatCounters< int >
CbmKFMath.h
L1_assert
#define L1_assert(v)
Definition: CbmL1Def.h:56
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition: CbmStsHit.h:31
L1FieldValue::y
fvec y
Definition: L1Field.h:17
L1FieldSlice::cx
fvec cx[21]
Definition: L1Field.h:34
TL1TracksCatCounters::AddCounter
void AddCounter()
Definition: CbmL1Counters.h:24
L1Algo::NStations
int NStations
Definition: L1Algo.h:333
CbmKF::Instance
static CbmKF * Instance()
Definition: CbmKF.h:39
finite
T finite(T x)
Definition: CbmL1Def.h:21
L1TrackPar::C44
fvec C44
Definition: L1TrackPar.h:10
L1FieldRegion
Definition: L1Field.h:85
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
CbmMvdPoint
Definition: CbmMvdPoint.h:28
EnergyLossCorrection
void EnergyLossCorrection(L1TrackPar &T, const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1, fvec=0)
Definition: L1AddMaterial.h:235
L1Station::z
fvec z
Definition: L1Station.h:28
TL1PerfEfficiencies::reco_fakes
TL1TracksCatCounters< double > reco_fakes
Definition: CbmL1Performance.cxx:287
L1AddPipeMaterial
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.10565f *0.10565f)
Definition: L1AddMaterial.h:737
TL1PerfEfficiencies::CalcEff
void CalcEff()
Definition: CbmL1Performance.cxx:205
CbmL1::TrackMatch
void TrackMatch()
Reconstruction Performance.
Definition: CbmL1Performance.cxx:60
CbmMuchPoint.h
CbmL1::fMuchPoints
CbmMCDataArray * fMuchPoints
Definition: CbmL1.h:294
TL1Efficiencies::nEvents
int nEvents
Definition: CbmL1Counters.h:157
CbmL1MCPoint::yIn
double yIn
Definition: CbmL1MCPoint.h:57
CbmL1MCTrack::p
double p
Definition: CbmL1MCTrack.h:31
CbmL1::fTimesliceMode
int fTimesliceMode
Definition: CbmL1.h:365
TL1Efficiencies::ghosts
int ghosts
Definition: CbmL1Counters.h:155
TL1PerfEfficiencies
Definition: CbmL1Performance.cxx:141
TL1PerfEfficiencies::ratio_length
TL1TracksCatCounters< double > ratio_length
Definition: CbmL1Performance.cxx:281
TL1Efficiencies::CalcEff
void CalcEff()
Definition: CbmL1Counters.h:169
CbmL1MCTrack::Points
vector< int > Points
Definition: CbmL1MCTrack.h:33
CbmL1MCTrack::mother_ID
int mother_ID
Definition: CbmL1MCTrack.h:32
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmL1MCTrack::z
double z
Definition: CbmL1MCTrack.h:31
L1TrackPar::tx
fvec tx
Definition: L1TrackPar.h:9
CbmL1::HistoPerformance
void HistoPerformance()
Definition: CbmL1Performance.cxx:612
CbmL1MCTrack::IsPrimary
bool IsPrimary() const
Definition: CbmL1MCTrack.h:104
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition: CbmKFTrack.h:59
CbmMCEventList::GetEventTime
Double_t GetEventTime(UInt_t event, UInt_t file)
Event start time.
Definition: CbmMCEventList.cxx:86
CbmL1MCTrack::ID
int ID
Definition: CbmL1MCTrack.h:32
CbmMatch::AddLink
void AddLink(const CbmLink &newLink)
Definition: CbmMatch.cxx:42
CbmL1MCPoint::iStation
int iStation
Definition: CbmL1MCPoint.h:61
CbmL1MCPoint::xOut
double xOut
Definition: CbmL1MCPoint.h:58
TL1PerfEfficiencies::PrintEff
void PrintEff()
Definition: CbmL1Performance.cxx:234
CbmL1MCPoint::p
double p
Definition: CbmL1MCPoint.h:59
CbmL1MCTrack::AddTouchTrack
void AddTouchTrack(CbmL1Track *tTr)
Definition: CbmL1MCTrack.h:119
L1FieldSlice
Definition: L1Field.h:31
CbmL1MCTrack::IsAdditional
bool IsAdditional() const
Definition: CbmL1MCTrack.h:106
TL1Efficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition: CbmL1Counters.h:160
CbmL1MCPoint::pyIn
double pyIn
Definition: CbmL1MCPoint.h:57
CbmTrdHit.h
Class for hits in TRD detector.
L1TrackPar::x
fvec x
Definition: L1TrackPar.h:9
CbmKFTrack.h
CbmL1MCTrack::IsReconstructed
bool IsReconstructed() const
Definition: CbmL1MCTrack.h:117
L1TrackPar::C00
fvec C00
Definition: L1TrackPar.h:9
CbmL1::fStsPoints
CbmMCDataArray * fStsPoints
Definition: CbmL1.h:266
CbmL1MCTrack::GetRecoTracks
vector< CbmL1Track * > & GetRecoTracks()
Definition: CbmL1MCTrack.h:115
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmL1Track::GetMaxPurity
double GetMaxPurity()
Definition: CbmL1Track.h:53
TL1Efficiencies::ratio_reco
TL1TracksCatCounters< double > ratio_reco
Definition: CbmL1Counters.h:149
L1FieldSlice::GetFieldValue
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition: L1Field.h:41
L1Extrapolation.h
CbmL1Counters.h
TL1Efficiencies::indices
map< TString, int > indices
Definition: CbmL1Counters.h:147
CbmL1Track
Definition: CbmL1Track.h:33
CbmL1Track::AddMCTrack
void AddMCTrack(CbmL1MCTrack *mcTr)
Definition: CbmL1Track.h:46
CbmL1MCPoint::zIn
double zIn
Definition: CbmL1MCPoint.h:57
CbmTrdPoint
Definition: CbmTrdPoint.h:23
CbmL1::NStation
int NStation
Definition: CbmL1.h:244
CbmStsSetup::GetStation
CbmStsStation * GetStation(Int_t stationId) const
Definition: CbmStsSetup.h:97
CbmL1MCPoint::pxOut
double pxOut
Definition: CbmL1MCPoint.h:58
CbmL1::TrackFitPerformance
void TrackFitPerformance()
Definition: CbmL1Performance.cxx:1543
L1TrackPar::C33
fvec C33
Definition: L1TrackPar.h:9
CbmL1::EfficienciesPerformance
void EfficienciesPerformance()
Definition: CbmL1Performance.cxx:293
L1Algo::CATime
double CATime
Definition: L1Algo.h:353
CbmL1MCTrack::AddRecoTrack
void AddRecoTrack(CbmL1Track *rTr)
Definition: CbmL1MCTrack.h:114
CbmStsHit::GetBackClusterId
Int_t GetBackClusterId() const
Definition: CbmStsHit.h:69
CbmTrdPoint::GetXIn
Double_t GetXIn() const
Definition: CbmTrdPoint.h:63
CbmL1StsHit
Definition: CbmL1StsHit.h:8
TL1PerfEfficiencies::mc_length_hits
TL1TracksCatCounters< int > mc_length_hits
Definition: CbmL1Performance.cxx:289
CbmKFTube
Definition: CbmKFMaterial.h:77
L1AddMaterial.h
CbmL1MCTrack::time
double time
Definition: CbmL1MCTrack.h:31
CbmL1::vRTracks
vector< CbmL1Track > vRTracks
Definition: CbmL1.h:125
CbmL1MCPoint::pyOut
double pyOut
Definition: CbmL1MCPoint.h:58
CbmMuchPoint::GetYIn
Double_t GetYIn() const
Definition: CbmMuchPoint.h:71
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmMvdPoint::GetXOut
Double_t GetXOut() const
Definition: CbmMvdPoint.h:70
CbmStsPoint::GetX
Double_t GetX(Double_t z) const
Definition: CbmStsPoint.cxx:96
L1TrackPar
Definition: L1TrackPar.h:6
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
CbmTofPoint.h
CbmKFMath::CopyTC2TrackParam
static void CopyTC2TrackParam(FairTrackParam *par, Double_t T[], Double_t C[])
Definition: CbmKFMath.cxx:809
L1StsHit.h
TL1PerfEfficiencies::TL1PerfEfficiencies
TL1PerfEfficiencies()
Definition: CbmL1Performance.cxx:142
L1TrackPar::C55
fvec C55
Definition: L1TrackPar.h:10
CbmL1HitStore::y
double y
Definition: CbmL1.h:98
CbmMuchPixelHit.h
Class for pixel hits in MUCH detector.
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmL1MCTrack::pz
double pz
Definition: CbmL1MCTrack.h:31
CbmTrdPoint.h
TL1Efficiencies
Definition: CbmL1Counters.h:117
CbmMvdPoint::GetYOut
Double_t GetYOut() const
Definition: CbmMvdPoint.h:71
L1TrackPar::C11
fvec C11
Definition: L1TrackPar.h:9
CbmMuchPoint::GetXIn
Double_t GetXIn() const
Definition: CbmMuchPoint.h:70
CbmL1::FieldIntegralCheck
void FieldIntegralCheck()
Definition: CbmL1Performance.cxx:2199
CbmL1::vHitMCRef
vector< int > vHitMCRef
Definition: CbmL1.h:339
ID
int ID
Definition: CbmMvdSensorDigiToHitTask.cxx:71
TL1Efficiencies::ratio_ghosts
double ratio_ghosts
Definition: CbmL1Counters.h:150
CbmL1MCTrack::StsHits
vector< int > StsHits
Definition: CbmL1MCTrack.h:34
CbmTofPoint
Geometric intersection of a MC track with a TOFb detector.
Definition: CbmTofPoint.h:40
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
CbmL1MCPoint::pzOut
double pzOut
Definition: CbmL1MCPoint.h:58
CbmL1::fMuchPixelHits
TClonesArray * fMuchPixelHits
Definition: CbmL1.h:299
CbmL1MCPoint::px
double px
Definition: CbmL1MCPoint.h:56
CbmMuchPixelHit
Definition: CbmMuchPixelHit.h:17
CbmL1MCPoint::yOut
double yOut
Definition: CbmL1MCPoint.h:58
CbmL1MCTrack
Definition: CbmL1MCTrack.h:29
CbmL1TrackPar::C
double C[21]
Definition: CbmL1TrackPar.h:17
CbmL1::vMCPoints
vector< CbmL1MCPoint > vMCPoints
Definition: CbmL1.h:197
TL1TracksCatCounters::NCounters
int NCounters
Definition: CbmL1Counters.h:110
CbmStsStation::GetXmax
Double_t GetXmax() const
Definition: CbmStsStation.h:117
CbmL1::FieldApproxCheck
void FieldApproxCheck()
Definition: CbmL1Performance.cxx:2047
CbmL1MCPoint::pxIn
double pxIn
Definition: CbmL1MCPoint.h:57
z1
Double_t z1[nSects1]
Definition: pipe_v16a_mvdsts100.h:6
CbmL1::fTrdHitMatches
TClonesArray * fTrdHitMatches
Definition: CbmL1.h:309
CbmKFTrack
Definition: CbmKFTrack.h:21
PipeRadThick
const fvec PipeRadThick
Definition: L1AddMaterial.h:11
L1Algo::fRadThick
vector< L1Material > fRadThick
Definition: L1Algo.h:337
CbmL1::listStsClusterMatch
TClonesArray * listStsClusterMatch
Definition: CbmL1.h:278
CbmL1MCPoint::xIn
double xIn
Definition: CbmL1MCPoint.h:57
PairAnalysisStyler::Fill
static Int_t Fill[]
Definition: PairAnalysisStyleDefs.h:82
CbmTrdPoint::GetYIn
Double_t GetYIn() const
Definition: CbmTrdPoint.h:64
TL1Efficiencies::mc
TL1TracksCatCounters< int > mc
Definition: CbmL1Counters.h:153
CbmL1MCPoint::py
double py
Definition: CbmL1MCPoint.h:56
CbmL1Track::IsGhost
bool IsGhost()
Definition: CbmL1Track.h:50
CbmL1::vHitStore
vector< CbmL1HitStore > vHitStore
Definition: CbmL1.h:180
L1FieldValue
Definition: L1Field.h:11
CbmL1::vMCTracks
vector< CbmL1MCTrack > vMCTracks
Definition: CbmL1.h:337
CbmL1::listStsHits
TClonesArray * listStsHits
Definition: CbmL1.h:276
CbmPixelHit::PositionError
void PositionError(TVector3 &dpos) const
Copies hit position error to pos.
Definition: CbmPixelHit.cxx:68
CbmL1MCTrack::NStations
int NStations() const
Definition: CbmL1MCTrack.h:107
CbmKFTrack::SetTrackParam
void SetTrackParam(const FairTrackParam &track)
Definition: CbmKFTrack.cxx:34
CbmStsPoint::GetY
Double_t GetY(Double_t z) const
Definition: CbmStsPoint.cxx:106
CbmL1MCPoint::zOut
double zOut
Definition: CbmL1MCPoint.h:58
CbmStsStation
Class representing a station of the StsSystem.
Definition: CbmStsStation.h:29
CbmL1MCTrack::IsReconstructable
bool IsReconstructable() const
Definition: CbmL1MCTrack.h:105
CbmL1::fTofPoints
CbmMCDataArray * fTofPoints
Definition: CbmL1.h:312
TL1Efficiencies::IncNEvents
void IncNEvents()
Definition: CbmL1Counters.h:140
L1FieldValue::x
fvec x
Definition: L1Field.h:15
L1TrackPar::C22
fvec C22
Definition: L1TrackPar.h:9
CbmL1MCPoint::pzIn
double pzIn
Definition: CbmL1MCPoint.h:57
CbmL1MCTrack::GetNClones
int GetNClones() const
Definition: CbmL1MCTrack.h:116
CbmL1::listMuchHitMatches
TClonesArray * listMuchHitMatches
Definition: CbmL1.h:295
CbmL1Track::SetMaxPurity
void SetMaxPurity(double maxPurity_)
Definition: CbmL1Track.h:52
CbmL1TrackPar::chi2
double chi2
Definition: CbmL1TrackPar.h:17
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0)
Access to i-th hit.
Definition: CbmKFTrackInterface.cxx:39
CbmL1::histodir
TDirectory * histodir
Definition: CbmL1.h:352