CbmL1Performance.cxx 94.3 KB
Newer Older
1
2
/* Copyright (C) 2006-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
3
   Authors: Ivan Kisel,  Sergey Gorbunov [committer], Igor Kulakov, Valentina Akishina, Grigory Kozlov */
4

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction
 *
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de
 *
 *====================================================================
 *
 *  L1 Fit performance
 *
 *====================================================================
 */
20
21
#include "CbmKF.h"
#include "CbmKFMath.h"
Administrator's avatar
Administrator committed
22
#include "CbmKFTrack.h"  // for vertex pulls
23
#include "CbmL1.h"
24
#include "CbmL1Constants.h"
25
26
#include "CbmL1Counters.h"
#include "CbmMatch.h"
27
28
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
Sergei Zharko's avatar
updates    
Sergei Zharko committed
29
#include "CbmQaTable.h"
30
31
32
#include "CbmStsDigi.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
33
34
#include "CbmStsSetup.h"
#include "CbmStsStation.h"
35
36
37
38
#include "CbmTofHit.h"
#include "CbmTofPoint.h"
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
39

40
#include "FairTrackParam.h"  // for vertex pulls
41

42
43
44
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
45
46
47
48
49
#include <TFile.h>

#include <iostream>
#include <list>
#include <map>
Administrator's avatar
Administrator committed
50
#include <vector>
51

52
53
54
55
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Extrapolation.h"  // for vertex pulls
#include "L1Algo/L1Fit.h"            // for vertex pulls

56
57
58
59
using std::cout;
using std::endl;
using std::ios;
using std::map;
Administrator's avatar
Administrator committed
60
using std::pair;
61
using std::setw;
Administrator's avatar
Administrator committed
62
using std::vector;
63

64
65
void CbmL1::TrackMatch()
{
Administrator's avatar
Administrator committed
66
67
  map<int, CbmL1MCTrack*> pMCTrackMap;
  pMCTrackMap.clear();
68

Administrator's avatar
Administrator committed
69
  // fill pMCTrackMap
70
  for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) {
Administrator's avatar
Administrator committed
71
    CbmL1MCTrack& MC = *i;
72

73
74
    if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()) { pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC)); }
    else {
75
76
77
      cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
    }
  }
Administrator's avatar
Administrator committed
78
  // -- prepare information about reconstructed tracks
79
  const int nRTracks = vRTracks.size();
Administrator's avatar
Administrator committed
80
  for (int iR = 0; iR < nRTracks; iR++) {
81
82
    CbmL1Track* prtra = &(vRTracks[iR]);

Administrator's avatar
Administrator committed
83
84
85
    //  cout<<iR<<" iR"<<endl;

    int hitsum = prtra->StsHits.size();  // number of hits in track
86

Administrator's avatar
Administrator committed
87
    // count how many hits from each mcTrack belong to current recoTrack
88
89
    map<int, int>& hitmap = prtra->hitMap;  // how many hits from each mcTrack belong to current recoTrack
    for (vector<int>::iterator ih = (prtra->StsHits).begin(); ih != (prtra->StsHits).end(); ++ih) {
90
91

      const int nMCPoints = vStsHits[*ih].mcPointIds.size();
Administrator's avatar
Administrator committed
92
      for (int iP = 0; iP < nMCPoints; iP++) {
93
        int iMC = vStsHits[*ih].mcPointIds[iP];
Administrator's avatar
Administrator committed
94
95

        //     cout<<iMC<<" iMC"<<endl;
96
97
        int ID = -1;
        if (iMC >= 0) ID = vMCPoints[iMC].ID;
98
        if (hitmap.find(ID) == hitmap.end()) hitmap[ID] = 1;
Administrator's avatar
Administrator committed
99
        else {
100
101
          hitmap[ID] += 1;
        }
Administrator's avatar
Administrator committed
102
103
      }  // for iPoint
    }    // for iHit
104

Administrator's avatar
Administrator committed
105
    // RTrack <-> MCTrack identification
106
    double max_percent = 0.0;  // [%]. maximum persent of hits, which belong to one mcTrack.
Administrator's avatar
Administrator committed
107
108
    for (map<int, int>::iterator posIt = hitmap.begin(); posIt != hitmap.end();
         posIt++) {  // loop over all touched MCTracks
109

Administrator's avatar
Administrator committed
110
      if (posIt->first < 0) continue;  // not a MC track - based on fake hits
111

Administrator's avatar
Administrator committed
112
      // count max-purity
113
      if (double(posIt->second) > max_percent * double(hitsum)) max_percent = double(posIt->second) / double(hitsum);
114

Administrator's avatar
Administrator committed
115
      // set relation to the mcTrack
116
117
      if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
      CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
118

119
      if (double(posIt->second) >= CbmL1Constants::MinPurity * double(hitsum)) {  // found correspondent MCTrack
120
121
        pmtra->AddRecoTrack(prtra);
        prtra->AddMCTrack(pmtra);
122
123
      }
      else {
124
125
        pmtra->AddTouchTrack(prtra);
      }
Administrator's avatar
Administrator committed
126
    }  // for hitmap
127
    prtra->SetMaxPurity(max_percent);
Administrator's avatar
Administrator committed
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
  }  // for reco tracks
}  //


struct TL1PerfEfficiencies : public TL1Efficiencies {
  TL1PerfEfficiencies()
    : TL1Efficiencies()
    , ratio_killed()
    , ratio_clone()
    , ratio_length()
    , ratio_fakes()
    , killed()
    , clone()
    , reco_length()
    , reco_fakes()
    , mc_length()
144
145
    , mc_length_hits()
  {
Administrator's avatar
Administrator committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    // add total efficiency
    AddCounter("long_fast_prim", "LongRPrim efficiency");
    AddCounter("fast_prim", "RefPrim   efficiency");
    AddCounter("fast_sec", "RefSec    efficiency");
    AddCounter("fast", "Refset    efficiency");
    AddCounter("total", "Allset    efficiency");
    AddCounter("slow_prim", "ExtraPrim efficiency");
    AddCounter("slow_sec", "ExtraSec  efficiency");
    AddCounter("slow", "Extra     efficiency");
    AddCounter("d0", "D0        efficiency");
    AddCounter("short", "Short123s efficiency");
    AddCounter("shortPion", "Short123s pion   eff");
    AddCounter("shortProton", "Short123s proton eff");
    AddCounter("shortKaon", "Short123s kaon   eff");
    AddCounter("shortE", "Short123s e      eff");
    AddCounter("shortRest", "Short123s rest   eff");

    AddCounter("fast_sec_e", "RefSec  E efficiency");
    AddCounter("fast_e", "Refset  E efficiency");
    AddCounter("total_e", "Allset  E efficiency");
    AddCounter("slow_sec_e", "ExtraSecE efficiency");
    AddCounter("slow_e", "Extra   E efficiency");
168
169
  }

Administrator's avatar
Administrator committed
170
  virtual ~TL1PerfEfficiencies() {};
171

172
173
  virtual void AddCounter(TString shortname, TString name)
  {
174
175
176
177
178
179
180
181
182
183
184
185
186
    TL1Efficiencies::AddCounter(shortname, name);
    ratio_killed.AddCounter();
    ratio_clone.AddCounter();
    ratio_length.AddCounter();
    ratio_fakes.AddCounter();
    killed.AddCounter();
    clone.AddCounter();
    reco_length.AddCounter();
    reco_fakes.AddCounter();
    mc_length.AddCounter();
    mc_length_hits.AddCounter();
  };

187
188
  TL1PerfEfficiencies& operator+=(TL1PerfEfficiencies& a)
  {
189
    TL1Efficiencies::operator+=(a);
Administrator's avatar
Administrator committed
190
191
    killed += a.killed;
    clone += a.clone;
192
    reco_length += a.reco_length;
Administrator's avatar
Administrator committed
193
194
    reco_fakes += a.reco_fakes;
    mc_length += a.mc_length;
195
196
197
198
    mc_length_hits += a.mc_length_hits;
    return *this;
  };

199
200
  void CalcEff()
  {
201
    TL1Efficiencies::CalcEff();
Administrator's avatar
Administrator committed
202
203
    ratio_killed                      = killed / mc;
    ratio_clone                       = clone / mc;
204
    TL1TracksCatCounters<int> allReco = reco + clone;
Administrator's avatar
Administrator committed
205
206
    ratio_length                      = reco_length / allReco;
    ratio_fakes                       = reco_fakes / allReco;
207
208
  };

209
210
211
  void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length,
           int _mc_length_hits, TString name)
  {
212
213
214
215
216
217
218
219
220
221
222
223
    TL1Efficiencies::Inc(isReco, name);

    const int index = indices[name];

    if (isKilled) killed.counters[index]++;
    reco_length.counters[index] += _ratio_length;
    reco_fakes.counters[index] += _ratio_fakes;
    clone.counters[index] += _nclones;
    mc_length.counters[index] += _mc_length;
    mc_length_hits.counters[index] += _mc_length_hits;
  };

224
  void PrintEff(bool ifPrintTableToLog = false, const std::string& nameOfTable = "efficiency_table")
225
  {
226
227
    L1_assert(nEvents != 0);

Sergei Zharko's avatar
update    
Sergei Zharko committed
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
    // cout.setf(ios::fixed);
    // cout.setf(ios::showpoint);
    // cout.precision(3);
    // cout.setf(ios::right);
    // cout << "Track category         : "
    //      << " Eff  "
    //      << " / "
    //      << "Killed"
    //      << " / "
    //      << "Length"
    //      << " / "
    //      << "Fakes "
    //      << " / "
    //      << "Clones"
    //      << " / "
    //      << "All Reco"
    //      << " | "
    //      << "  All MC "
    //      << " / "
    //      << "MCl(hits)"
    //      << " / "
    //      << "MCl(MCps)" << endl;
250
    int NCounters = mc.GetNcounters();
Sergei Zharko's avatar
update    
Sergei Zharko committed
251
252
253
254
    std::vector<std::string> rowNames(20);
    for (int iC = 0; iC < NCounters; ++iC) {
      rowNames[iC] = std::string(names[iC].Data());
    }
Sergei Zharko's avatar
updates    
Sergei Zharko committed
255
256
    std::vector<std::string> colNames = {"Eff.",     "Killed", "Length",    "Fakes",    "Clones",
                                         "All Reco", "All MC", "MCl(hits)", "MCl(MCps)"};
Sergei Zharko's avatar
update    
Sergei Zharko committed
257
258
259
260

    TDirectory* curdir = gDirectory;
    gDirectory         = fOutDir;

Sergei Zharko's avatar
updates    
Sergei Zharko committed
261
    CbmQaTable* aTable = new CbmQaTable(nameOfTable.c_str(), "Track Efficiency", 20, 9);
Sergei Zharko's avatar
update    
Sergei Zharko committed
262
263
    aTable->SetNamesOfRows(rowNames);
    aTable->SetNamesOfCols(colNames);
Administrator's avatar
Administrator committed
264
    for (int iC = 0; iC < NCounters; iC++) {
Sergei Zharko's avatar
update    
Sergei Zharko committed
265
266
267
268
269
270
271
272
273
      aTable->SetCell(iC, 0, ratio_reco.counters[iC]);
      aTable->SetCell(iC, 1, ratio_killed.counters[iC]);
      aTable->SetCell(iC, 2, ratio_length.counters[iC]);
      aTable->SetCell(iC, 3, ratio_fakes.counters[iC]);
      aTable->SetCell(iC, 4, ratio_clone.counters[iC]);
      aTable->SetCell(iC, 5, reco.counters[iC] / double(nEvents));
      aTable->SetCell(iC, 6, mc.counters[iC] / double(nEvents));
      aTable->SetCell(iC, 7, mc_length_hits.counters[iC] / double(mc.counters[iC]));
      aTable->SetCell(iC, 8, mc_length.counters[iC] / double(mc.counters[iC]));
274
    }
275
276
277
    if (ifPrintTableToLog) {
      cout << *aTable;  // print a table to log
    }
278
    cout << "Ghost     probability  : " << ratio_ghosts << "  | " << ghosts << endl;
Sergei Zharko's avatar
updates    
Sergei Zharko committed
279

Sergei Zharko's avatar
update    
Sergei Zharko committed
280
    gDirectory = curdir;
281
282
  };

Sergei Zharko's avatar
update    
Sergei Zharko committed
283

284
285
286
287
288
289
290
291
292
293
294
  TL1TracksCatCounters<double> ratio_killed;
  TL1TracksCatCounters<double> ratio_clone;
  TL1TracksCatCounters<double> ratio_length;
  TL1TracksCatCounters<double> ratio_fakes;

  TL1TracksCatCounters<int> killed;
  TL1TracksCatCounters<int> clone;
  TL1TracksCatCounters<double> reco_length;
  TL1TracksCatCounters<double> reco_fakes;
  TL1TracksCatCounters<int> mc_length;
  TL1TracksCatCounters<int> mc_length_hits;
Sergei Zharko's avatar
updates    
Sergei Zharko committed
295
296

  TDirectory* fOutDir {nullptr};  // Specified for saving tables
297
298
299
};


300
301
void CbmL1::EfficienciesPerformance()
{
Administrator's avatar
Administrator committed
302
  static TL1PerfEfficiencies L1_NTRA;  // average efficiencies
303

Administrator's avatar
Administrator committed
304
305
  static int L1_NEVENTS   = 0;
  static double L1_CATIME = 0.0;
306

307

Administrator's avatar
Administrator committed
308
  TL1PerfEfficiencies ntra;  // efficiencies for current event
Sergei Zharko's avatar
updates    
Sergei Zharko committed
309
310
  ntra.fOutDir    = fHistoDir;  // Setup a pointer for output directory
  L1_NTRA.fOutDir = fHistoDir;  // Save average efficiencies
311

312
  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt) {
313
    ntra.ghosts += rtraIt->IsGhost();
Administrator's avatar
Administrator committed
314
315
316
317
318
319
320
    //     if(rtraIt->IsGhost()){ // Debug.
    //       cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
    //       for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++ ){
    //         cout << " (" << posIt->second << " " << posIt->first << ") ";
    //       }
    //       cout << endl;
    //     }
321
322
  }

323
  int sta_nhits[algo->GetNstations()], sta_nfakes[algo->GetNstations()];
324

325
  for (int i = 0; i < algo->GetNstations(); i++) {
326
327
328
329
    sta_nhits[i]  = 0;
    sta_nfakes[i] = 0;
  }

330
  for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++) {
Administrator's avatar
Administrator committed
331
    CbmL1MCTrack& mtra = *(mtraIt);
332

333

Administrator's avatar
Administrator committed
334
335
    //    if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
    if (!mtra.IsReconstructable() && !mtra.IsAdditional()) continue;
336

Administrator's avatar
Administrator committed
337
338
    // -- find used constans --
    // is track reconstructed
339
    const bool reco = (mtra.IsReconstructed());
Administrator's avatar
Administrator committed
340
    // is track killed. At least one hit of it belong to some recoTrack
341
    const bool killed = !reco && mtra.IsDisturbed();
Administrator's avatar
Administrator committed
342
    // ration length for current mcTrack
343
344
345
346
    L1Vector<CbmL1Track*>& rTracks = mtra.GetRecoTracks();  // for length calculations
    double ratio_length            = 0;
    double ratio_fakes             = 0;
    double mc_length_hits          = mtra.NStations();
347

Administrator's avatar
Administrator committed
348

349
    int mc_length = mtra.NMCStations();
Administrator's avatar
Administrator committed
350
    if (reco) {
351
      for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
352
        ratio_length += static_cast<double>(rTracks[irt]->GetNOfHits()) * rTracks[irt]->GetMaxPurity() / mc_length_hits;
353
354
355
        ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
      }
    }
Administrator's avatar
Administrator committed
356
357

    // number of clones
358
359
    int nclones = 0;
    if (reco) nclones = mtra.GetNClones();
Administrator's avatar
Administrator committed
360
361
362
363
364
365
366
367
368
    //     if (nclones){ // Debug. Look at clones
    //       for (int irt = 0; irt < rTracks.size(); irt++){
    //         const int ista = vHitStore[rTracks[irt]->StsHits[0]].iStation;
    //         cout << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
    //       }
    //       cout << mtra.NStations() << endl;
    //     }

    if (mtra.IsAdditional()) {  // short
369
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
Administrator's avatar
Administrator committed
370
      switch (mtra.pdg) {
371
372
        case 11:
        case -11:
373
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
374
375
376
          break;
        case 211:
        case -211:
377
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortPion");
378
379
380
          break;
        case 321:
        case -321:
381
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortKaon");
382
383
384
          break;
        case 2212:
        case -2212:
385
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortProton");
386
          break;
387
        default: ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortRest");
388
      }
389
390
    }
    else {  // separate all efficiecies from short eff
Administrator's avatar
Administrator committed
391
392


393
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
Administrator's avatar
Administrator committed
394

395
      if (mtra.isSignal) {  // D0
396
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
397
398
      }

Administrator's avatar
Administrator committed
399
      if (mtra.p > CbmL1Constants::MinRefMom) {  // reference tracks
400
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast");
Administrator's avatar
Administrator committed
401
402
403

        if (mtra.IsPrimary()) {                // reference primary
          if (mtra.NStations() == NStation) {  // long reference primary
404
            ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "long_fast_prim");
405
          }
406
407
408
409
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_prim");
        }
        else {  // reference secondary
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec");
410
        }
411
412
413
      }
      else {  // extra set of tracks
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow");
Administrator's avatar
Administrator committed
414
415

        if (mtra.IsPrimary()) {  // extra primary
416
417
418
419
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_prim");
        }
        else {
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec");
420
        }
Administrator's avatar
Administrator committed
421
422
423
      }  // if extra

      if (mtra.pdg == 11 || mtra.pdg == -11) {
424
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");
Administrator's avatar
Administrator committed
425
426

        if (mtra.p > CbmL1Constants::MinRefMom) {  // reference tracks
427
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
Administrator's avatar
Administrator committed
428
429

          if (mtra.IsPrimary()) {  // reference primary
430
          }
431
432
433
434
435
436
          else {  // reference secondary
            ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec_e");
          }
        }
        else {  // extra set of tracks
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
Administrator's avatar
Administrator committed
437
438

          if (mtra.IsPrimary()) {  // extra primary
439
440
441
          }
          else {
            ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec_e");
442
          }
Administrator's avatar
Administrator committed
443
        }  // if extra
444
445
446
      }
    }

Administrator's avatar
Administrator committed
447
  }  // for mcTracks
448

449
  L1_CATIME += fTrackingTime;
450
451
452
453
454
455
456
457
  L1_NEVENTS++;
  ntra.IncNEvents();
  L1_NTRA += ntra;

  ntra.CalcEff();
  L1_NTRA.CalcEff();

  //   cout.precision(3);
Administrator's avatar
Administrator committed
458
459
  if (fVerbose) {
    if (fVerbose > 1) {
460
461
      ntra.PrintEff();
      cout << "Number of true and fake hits in stations: " << endl;
462
      for (int i = 0; i < algo->GetNstations(); i++) {
Administrator's avatar
Administrator committed
463
        cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << "   ";
464
465
      }
      cout << endl;
Administrator's avatar
Administrator committed
466
    }  // fVerbose > 1
467
    cout << endl;
468
    cout << "L1 ACCUMULATED STAT    : " << L1_NEVENTS << " EVENTS " << endl << endl;
469
    L1_NTRA.PrintEff(/*ifPrintTableToLog = */ true);
Administrator's avatar
Administrator committed
470
    cout << "MC tracks/event found  : "
471
         << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS)) << endl;
Administrator's avatar
Administrator committed
472
    cout << endl;
473
    cout << "CA Track Finder: " << L1_CATIME / L1_NEVENTS << (fLegacyEventMode ? " s/ev" : " s/time slice") << endl
474
         << endl;
475
  }
Administrator's avatar
Administrator committed
476
}  // void CbmL1::Performance()
477

478
void CbmL1::HistoPerformance()  // TODO: check if works correctly. Change vHitRef on match data in CbmL1**Track classes
479
480
481
482
{

  //CbmKF &KF = *CbmKF::Instance();

483
484
  static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom, *p_eff_prim_vs_theta,
    *p_eff_all_vs_pt, *p_eff_prim_vs_pt, *p_eff_all_vs_nhits, *p_eff_prim_vs_nhits, *p_eff_sec_vs_nhits;
485
486
487
488
489

  static TH1F *h_reg_MCmom, *h_acc_MCmom, *h_reco_MCmom, *h_ghost_Rmom;
  static TH1F *h_reg_prim_MCmom, *h_acc_prim_MCmom, *h_reco_prim_MCmom;
  static TH1F *h_reg_sec_MCmom, *h_acc_sec_MCmom, *h_reco_sec_MCmom;

Administrator's avatar
Administrator committed
490
  static TH1F* h_acc_mom_short123s;
491

492
493
494
495
  static TH1F *h_reg_mom_prim, *h_reg_mom_sec, *h_reg_nhits_prim, *h_reg_nhits_sec;
  static TH1F *h_acc_mom_prim, *h_acc_mom_sec, *h_acc_nhits_prim, *h_acc_nhits_sec;
  static TH1F *h_reco_mom_prim, *h_reco_mom_sec, *h_reco_nhits_prim, *h_reco_nhits_sec;
  static TH1F *h_rest_mom_prim, *h_rest_mom_sec, *h_rest_nhits_prim, *h_rest_nhits_sec;
496
497
498

  //static TH1F *h_hit_density[10];

499
500
501
  static TH1F *h_ghost_mom, *h_ghost_nhits, *h_ghost_fstation, *h_ghost_chi2, *h_ghost_prob, *h_ghost_tx, *h_ghost_ty;
  static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station, *h_reco_chi2, *h_reco_prob, *h_rest_prob,
    *h_reco_clean, *h_reco_time;
Administrator's avatar
Administrator committed
502
503
  static TProfile *h_reco_timeNtr, *h_reco_timeNhit;
  static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit;
504
505
  static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;

506
  static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
507
508

  static TH1F *h_mcp, *h_nmchits;
Administrator's avatar
Administrator committed
509
  //  static TH1F *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
510

Administrator's avatar
Administrator committed
511
  //  static TH2F *P_vs_P ;
512
513
514
515

  static TH2F *h2_vertex, *h2_prim_vertex, *h2_sec_vertex;
  //static TH3F *h3_vertex, *h3_prim_vertex, *h3_sec_vertex;

516
517
518
519
520
521
522
523
524
  static TH2F *h2_reg_nhits_vs_mom_prim, *h2_reg_nhits_vs_mom_sec, *h2_reg_fstation_vs_mom_prim,
    *h2_reg_fstation_vs_mom_sec, *h2_reg_lstation_vs_fstation_prim, *h2_reg_lstation_vs_fstation_sec;
  static TH2F *h2_acc_nhits_vs_mom_prim, *h2_acc_nhits_vs_mom_sec, *h2_acc_fstation_vs_mom_prim,
    *h2_acc_fstation_vs_mom_sec, *h2_acc_lstation_vs_fstation_prim, *h2_acc_lstation_vs_fstation_sec;
  static TH2F *h2_reco_nhits_vs_mom_prim, *h2_reco_nhits_vs_mom_sec, *h2_reco_fstation_vs_mom_prim,
    *h2_reco_fstation_vs_mom_sec, *h2_reco_lstation_vs_fstation_prim, *h2_reco_lstation_vs_fstation_sec;
  static TH2F *h2_ghost_nhits_vs_mom, *h2_ghost_fstation_vs_mom, *h2_ghost_lstation_vs_fstation;
  static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec, *h2_rest_fstation_vs_mom_prim,
    *h2_rest_fstation_vs_mom_sec, *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
525
526
527

  static bool first_call = 1;

Administrator's avatar
Administrator committed
528
  if (first_call) {
529
530
    first_call = 0;

Administrator's avatar
Administrator committed
531
    TDirectory* curdir = gDirectory;
Sergey Gorbunov's avatar
cleanup    
Sergey Gorbunov committed
532
    gDirectory         = fHistoDir;
533

534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
    p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
    p_eff_prim_vs_mom =
      new TProfile("p_eff_prim_vs_mom", "Primary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
    p_eff_sec_vs_mom =
      new TProfile("p_eff_sec_vs_mom", "Secondary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
    p_eff_d0_vs_mom =
      new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum", 150, 0.0, 15.0, 0.0, 100.0);
    p_eff_prim_vs_theta =
      new TProfile("p_eff_prim_vs_theta", "All Primary Set Efficiency vs Theta", 100, 0.0, 30.0, 0.0, 100.0);
    p_eff_all_vs_pt  = new TProfile("p_eff_all_vs_pt", "AllSet Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
    p_eff_prim_vs_pt = new TProfile("p_eff_prim_vs_pt", "Primary Set Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);

    p_eff_all_vs_nhits = new TProfile("p_eff_all_vs_nhits", "AllSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
    p_eff_prim_vs_nhits =
      new TProfile("p_eff_prim_vs_nhits", "PrimSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
    p_eff_sec_vs_nhits = new TProfile("p_eff_sec_vs_nhits", "SecSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);

    h_reg_MCmom       = new TH1F("h_reg_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
    h_acc_MCmom       = new TH1F("h_acc_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
    h_reco_MCmom      = new TH1F("h_reco_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
    h_ghost_Rmom      = new TH1F("h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0);
    h_reg_prim_MCmom  = new TH1F("h_reg_prim_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
    h_acc_prim_MCmom  = new TH1F("h_acc_prim_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
    h_reco_prim_MCmom = new TH1F("h_reco_prim_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
    h_reg_sec_MCmom   = new TH1F("h_reg_sec_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
    h_acc_sec_MCmom   = new TH1F("h_acc_sec_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
    h_reco_sec_MCmom  = new TH1F("h_reco_sec_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
Administrator's avatar
Administrator committed
561
562

    h_acc_mom_short123s =
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
      new TH1F("h_acc_mom_short123s", "Momentum of accepted tracks with 3 hits on first stations", 500, 0.0, 5.0);

    h_reg_mom_prim  = new TH1F("h_reg_mom_prim", "Momentum of registered primary tracks", 500, 0.0, 5.0);
    h_reg_mom_sec   = new TH1F("h_reg_mom_sec", "Momentum of registered secondary tracks", 500, 0.0, 5.0);
    h_acc_mom_prim  = new TH1F("h_acc_mom_prim", "Momentum of accepted primary tracks", 500, 0.0, 5.0);
    h_acc_mom_sec   = new TH1F("h_acc_mom_sec", "Momentum of accepted secondary tracks", 500, 0.0, 5.0);
    h_reco_mom_prim = new TH1F("h_reco_mom_prim", "Momentum of reconstructed primary tracks", 500, 0.0, 5.0);
    h_reco_mom_sec  = new TH1F("h_reco_mom_sec", "Momentum of reconstructed secondary tracks", 500, 0.0, 5.0);
    h_rest_mom_prim = new TH1F("h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
    h_rest_mom_sec  = new TH1F("h_rest_mom_sec", "Momentum of not found secondary tracks", 500, 0.0, 5.0);

    h_reg_nhits_prim  = new TH1F("h_reg_nhits_prim", "Number of hits in registered primary track", 51, -0.1, 10.1);
    h_reg_nhits_sec   = new TH1F("h_reg_nhits_sec", "Number of hits in registered secondary track", 51, -0.1, 10.1);
    h_acc_nhits_prim  = new TH1F("h_acc_nhits_prim", "Number of hits in accepted primary track", 51, -0.1, 10.1);
    h_acc_nhits_sec   = new TH1F("h_acc_nhits_sec", "Number of hits in accepted secondary track", 51, -0.1, 10.1);
    h_reco_nhits_prim = new TH1F("h_reco_nhits_prim", "Number of hits in reconstructed primary track", 51, -0.1, 10.1);
    h_reco_nhits_sec  = new TH1F("h_reco_nhits_sec", "Number of hits in reconstructed secondary track", 51, -0.1, 10.1);
    h_rest_nhits_prim = new TH1F("h_rest_nhits_prim", "Number of hits in not found primary track", 51, -0.1, 10.1);
    h_rest_nhits_sec  = new TH1F("h_rest_nhits_sec", "Number of hits in not found secondary track", 51, -0.1, 10.1);

    h_ghost_mom      = new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
    h_ghost_nhits    = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
    h_ghost_fstation = new TH1F("h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
    h_ghost_chi2     = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
    h_ghost_prob     = new TH1F("h_ghost_prob", "Prob of ghost track", 505, 0., 1.01);
    h_ghost_r        = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0);
    h_ghost_tx       = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
    h_ghost_ty       = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);

    h_reco_mom      = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
    h_reco_nhits    = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
    h_reco_station  = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
    h_reco_chi2     = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
    h_reco_prob     = new TH1F("h_reco_prob", "Prob of reco track", 505, 0., 1.01);
    h_rest_prob     = new TH1F("h_rest_prob", "Prob of reco rest track", 505, 0., 1.01);
    h_reco_clean    = new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5);
    h_reco_time     = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
    h_reco_timeNtr  = new TProfile("h_reco_timeNtr", "CA Track Finder Time (s/ev) vs N Tracks", 200, 0.0, 1000.0);
    h_reco_timeNhit = new TProfile("h_reco_timeNhit", "CA Track Finder Time (s/ev) vs N Hits", 200, 0.0, 30000.0);
    h_reco_fakeNtr  = new TProfile("h_reco_fakeNtr", "N Fake hits vs N Tracks", 200, 0.0, 1000.0);
    h_reco_fakeNhit = new TProfile("h_reco_fakeNhit", "N Fake hits vs N Hits", 200, 0.0, 30000.0);

    h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);
Administrator's avatar
Administrator committed
606
607
608
609
610
611
612
613
614
615
616
617

    //     h_hit_density[0] = new TH1F("h_hit_density0", "Hit density station 1", 50, 0.0,  5.00);
    //     h_hit_density[1] = new TH1F("h_hit_density1", "Hit density station 2", 100, 0.0, 10.00);
    //     h_hit_density[2] = new TH1F("h_hit_density2", "Hit density station 3", 140, 0.0, 13.99);
    //     h_hit_density[3] = new TH1F("h_hit_density3", "Hit density station 4", 180, 0.0, 18.65);
    //     h_hit_density[4] = new TH1F("h_hit_density4", "Hit density station 5", 230, 0.0, 23.32);
    //     h_hit_density[5] = new TH1F("h_hit_density5", "Hit density station 6", 270, 0.0, 27.98);
    //     h_hit_density[6] = new TH1F("h_hit_density6", "Hit density station 7", 340, 0.0, 34.97);
    //     h_hit_density[7] = new TH1F("h_hit_density7", "Hit density station 8", 460, 0.0, 46.63);
    //     h_hit_density[8] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
    //     h_hit_density[9] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
    //     h_hit_density[10] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
618
619
620
621

    h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
    h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);

622
623
624
625
626
627
628
629
    h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);

    h_notfound_mom     = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
    h_notfound_nhits   = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
    h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
    h_notfound_r       = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
    h_notfound_tx      = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
    h_notfound_ty      = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);
Administrator's avatar
Administrator committed
630
631

    /*
632
633
634
635
636
637
    h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
    h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
    VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
    VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
*/
    h_mcp = new TH1F("h_mcp", "MC P ", 500, 0.0, 5.0);
Administrator's avatar
Administrator committed
638
    /*
639
640
641
642
643
644
    MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
    MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
    MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
*/
    h_nmchits = new TH1F("h_nmchits", "N STS hits on MC track", 50, 0.0, 10.0);

Administrator's avatar
Administrator committed
645
646
    //    P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);

647
648
649
    h2_vertex      = new TH2F("h2_vertex", "2D vertex distribution", 105, -5., 100., 100, -50., 50.);
    h2_prim_vertex = new TH2F("h2_primvertex", "2D primary vertex distribution", 105, -5., 100., 100, -50., 50.);
    h2_sec_vertex  = new TH2F("h2_sec_vertex", "2D secondary vertex distribution", 105, -5., 100., 100, -50., 50.);
Administrator's avatar
Administrator committed
650
651
652
653
654
655

    //h3_vertex = new TH3F("h3_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
    //h3_prim_vertex = new TH3F("h3_primvertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
    //h3_sec_vertex = new TH3F("h3_sec_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);

    h2_reg_nhits_vs_mom_prim =
656
657
658
      new TH2F("h2_reg_nhits_vs_mom_prim", "NHits vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
    h2_reg_nhits_vs_mom_sec = new TH2F("h2_reg_nhits_vs_mom_sec", "NHits vs. Momentum (Reg. Secondary Tracks)", 51,
                                       -0.05, 5.05, 11, -0.5, 10.5);
Administrator's avatar
Administrator committed
659
    h2_acc_nhits_vs_mom_prim =
660
661
662
663
664
665
666
667
668
669
670
671
672
      new TH2F("h2_acc_nhits_vs_mom_prim", "NHits vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
    h2_acc_nhits_vs_mom_sec   = new TH2F("h2_acc_nhits_vs_mom_sec", "NHits vs. Momentum (Acc. Secondary Tracks)", 51,
                                       -0.05, 5.05, 11, -0.5, 10.5);
    h2_reco_nhits_vs_mom_prim = new TH2F("h2_reco_nhits_vs_mom_prim", "NHits vs. Momentum (Reco Primary Tracks)", 51,
                                         -0.05, 5.05, 11, -0.5, 10.5);
    h2_reco_nhits_vs_mom_sec  = new TH2F("h2_reco_nhits_vs_mom_sec", "NHits vs. Momentum (Reco Secondary Tracks)", 51,
                                        -0.05, 5.05, 11, -0.5, 10.5);
    h2_ghost_nhits_vs_mom =
      new TH2F("h2_ghost_nhits_vs_mom", "NHits vs. Momentum (Ghost Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
    h2_rest_nhits_vs_mom_prim = new TH2F("h2_rest_nhits_vs_mom_prim", "NHits vs. Momentum (!Found Primary Tracks)", 51,
                                         -0.05, 5.05, 11, -0.5, 10.5);
    h2_rest_nhits_vs_mom_sec  = new TH2F("h2_rest_nhits_vs_mom_sec", "NHits vs. Momentum (!Found Secondary Tracks)", 51,
                                        -0.05, 5.05, 11, -0.5, 10.5);
Administrator's avatar
Administrator committed
673
674

    h2_reg_fstation_vs_mom_prim =
675
676
      new TH2F("h2_reg_fstation_vs_mom_prim", "First Station vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
677
    h2_reg_fstation_vs_mom_sec =
678
679
      new TH2F("h2_reg_fstation_vs_mom_sec", "First Station vs. Momentum (Reg. Secondary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
680
    h2_acc_fstation_vs_mom_prim =
681
682
      new TH2F("h2_acc_fstation_vs_mom_prim", "First Station vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
683
    h2_acc_fstation_vs_mom_sec =
684
685
      new TH2F("h2_acc_fstation_vs_mom_sec", "First Station vs. Momentum (Acc. Secondary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
686
    h2_reco_fstation_vs_mom_prim =
687
688
      new TH2F("h2_reco_fstation_vs_mom_prim", "First Station vs. Momentum (Reco Primary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
689
    h2_reco_fstation_vs_mom_sec =
690
691
692
693
      new TH2F("h2_reco_fstation_vs_mom_sec", "First Station vs. Momentum (Reco Secondary Tracks)", 51, -0.05, 5.05, 11,
               -0.5, 10.5);
    h2_ghost_fstation_vs_mom = new TH2F("h2_ghost_fstation_vs_mom", "First Station vs. Momentum (Ghost Tracks)", 51,
                                        -0.05, 5.05, 11, -0.5, 10.5);
Administrator's avatar
Administrator committed
694
    h2_rest_fstation_vs_mom_prim =
695
696
      new TH2F("h2_rest_fstation_vs_mom_prim", "First Station vs. Momentum (!Found Primary Tracks)", 51, -0.05, 5.05,
               11, -0.5, 10.5);
Administrator's avatar
Administrator committed
697
    h2_rest_fstation_vs_mom_sec =
698
699
      new TH2F("h2_rest_fstation_vs_mom_sec", "First Station vs. Momentum (!Found Secondary Tracks)", 51, -0.05, 5.05,
               11, -0.5, 10.5);
Administrator's avatar
Administrator committed
700
701

    h2_reg_lstation_vs_fstation_prim =
702
703
      new TH2F("h2_reg_lstation_vs_fstation_prim", "Last vs. First Station (Reg. Primary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
704
    h2_reg_lstation_vs_fstation_sec =
705
706
      new TH2F("h2_reg_lstation_vs_fstation_sec", "Last vs. First Station (Reg. Secondary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
707
    h2_acc_lstation_vs_fstation_prim =
708
709
      new TH2F("h2_acc_lstation_vs_fstation_prim", "Last vs. First Station (Acc. Primary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
710
    h2_acc_lstation_vs_fstation_sec =
711
712
      new TH2F("h2_acc_lstation_vs_fstation_sec", "Last vs. First Station (Acc. Secondary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
713
    h2_reco_lstation_vs_fstation_prim =
714
715
      new TH2F("h2_reco_lstation_vs_fstation_prim", "Last vs. First Station (Reco Primary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
Administrator's avatar
Administrator committed
716
    h2_reco_lstation_vs_fstation_sec =
717
718
719
720
      new TH2F("h2_reco_lstation_vs_fstation_sec", "Last vs. First Station (Reco Secondary Tracks)", 11, -0.5, 10.5, 11,
               -0.5, 10.5);
    h2_ghost_lstation_vs_fstation = new TH2F("h2_ghost_lstation_vs_fstation", "Last vs. First Station (Ghost Tracks)",
                                             11, -0.5, 10.5, 11, -0.5, 10.5);
Administrator's avatar
Administrator committed
721
    h2_rest_lstation_vs_fstation_prim =
722
723
      new TH2F("h2_rest_lstation_vs_fstation_prim", "Last vs. First Station (!Found Primary Tracks)", 11, -0.5, 10.5,
               11, -0.5, 10.5);
Administrator's avatar
Administrator committed
724
    h2_rest_lstation_vs_fstation_sec =
725
726
      new TH2F("h2_rest_lstation_vs_fstation_sec", "Last vs. First Station (!Found Secondary Tracks)", 11, -0.5, 10.5,
               11, -0.5, 10.5);
Administrator's avatar
Administrator committed
727
728
729
730

    //maindir->cd();

    // ----- Create list of all histogram pointers
731
732
733

    gDirectory = curdir;

Administrator's avatar
Administrator committed
734
  }  // first_call
735
736
737


  static int NEvents = 0;
Administrator's avatar
Administrator committed
738
  if (NEvents > 0) {
739
740
741
742
743
744
745
746
747
748
749
750
751
    h_reg_MCmom->Scale(NEvents);
    h_acc_MCmom->Scale(NEvents);
    h_reco_MCmom->Scale(NEvents);
    h_ghost_Rmom->Scale(NEvents);
    h_reg_prim_MCmom->Scale(NEvents);
    h_acc_prim_MCmom->Scale(NEvents);
    h_reco_prim_MCmom->Scale(NEvents);
    h_reg_sec_MCmom->Scale(NEvents);
    h_acc_sec_MCmom->Scale(NEvents);
    h_reco_sec_MCmom->Scale(NEvents);
  }
  NEvents++;

Administrator's avatar
Administrator committed
752
  //   //hit density calculation: h_hit_density[10]
753
  //
Administrator's avatar
Administrator committed
754
755
756
757
758
759
  //   for (vector<CbmL1HitStore>::iterator hIt = vHitStore.begin(); hIt != vHitStore.end(); ++hIt){
  //     float x = hIt->x;
  //     float y = hIt->y;
  //     float r = sqrt(x*x+y*y);
  //     h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
  //   }
760
761

  //
762
  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt) {
763
    CbmL1Track* prtra = &(*rtraIt);
Administrator's avatar
Administrator committed
764
    if ((prtra->StsHits).size() < 1) continue;
765
    {  // fill histos
Administrator's avatar
Administrator committed
766
      if (fabs(prtra->T[4]) > 1.e-10) h_reco_mom->Fill(fabs(1.0 / prtra->T[4]));
767
      h_reco_nhits->Fill((prtra->StsHits).size());
Administrator's avatar
Administrator committed
768
      CbmL1HitStore& mh = vHitStore[prtra->StsHits[0]];
769
770
771
      h_reco_station->Fill(mh.iStation);
    }

Administrator's avatar
Administrator committed
772
    h_reco_clean->Fill(prtra->GetMaxPurity());
773
774

    if (prtra->NDF > 0) {
Administrator's avatar
Administrator committed
775
776
777
      if (prtra->IsGhost()) {
        h_ghost_chi2->Fill(prtra->chi2 / prtra->NDF);
        h_ghost_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
778
779
      }
      else {
Administrator's avatar
Administrator committed
780
781
782
        if (prtra->GetMCTrack()[0].IsReconstructable()) {
          h_reco_chi2->Fill(prtra->chi2 / prtra->NDF);
          h_reco_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
783
784
        }
        else {
Administrator's avatar
Administrator committed
785
786
          //          h_rest_chi2->Fill(prtra->chi2/prtra->NDF);
          h_rest_prob->Fill(TMath::Prob(prtra->chi2, prtra->NDF));
787
788
789
790
791
        }
      }
    }


Administrator's avatar
Administrator committed
792
793
794
795
796
    // fill ghost histos
    if (prtra->IsGhost()) {
      if (fabs(prtra->T[4]) > 1.e-10) {
        h_ghost_mom->Fill(fabs(1.0 / prtra->T[4]));
        h_ghost_Rmom->Fill(fabs(1.0 / prtra->T[4]));
797
798
      }
      h_ghost_nhits->Fill((prtra->StsHits).size());
Administrator's avatar
Administrator committed
799
800
      CbmL1HitStore& h1 = vHitStore[prtra->StsHits[0]];
      CbmL1HitStore& h2 = vHitStore[prtra->StsHits[1]];
801
      h_ghost_fstation->Fill(h1.iStation);
Administrator's avatar
Administrator committed
802
      h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
803
804
      double z1 = algo->GetStations()[h1.iStation].z[0];
      double z2 = algo->GetStations()[h2.iStation].z[0];
Administrator's avatar
Administrator committed
805
806
807
      if (fabs(z2 - z1) > 1.e-4) {
        h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
        h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
808
809
      }

810
      if (fabs(prtra->T[4]) > 1.e-10) h2_ghost_nhits_vs_mom->Fill(fabs(1.0 / prtra->T[4]), (prtra->StsHits).size());
Administrator's avatar
Administrator committed
811
      CbmL1HitStore& hf = vHitStore[prtra->StsHits[0]];
812
813
814
      CbmL1HitStore& hl = vHitStore[prtra->StsHits[(prtra->StsHits).size() - 1]];
      if (fabs(prtra->T[4]) > 1.e-10) h2_ghost_fstation_vs_mom->Fill(fabs(1.0 / prtra->T[4]), hf.iStation + 1);
      if (hl.iStation >= hf.iStation) h2_ghost_lstation_vs_fstation->Fill(hf.iStation + 1, hl.iStation + 1);
815
816
    }

Administrator's avatar
Administrator committed
817
  }  // for reco tracks
818

Administrator's avatar
Administrator committed
819
  int mc_total = 0;  // total amount of reconstructable mcTracks
820
  for (vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++) {
Administrator's avatar
Administrator committed
821
822
    CbmL1MCTrack& mtra = *(mtraIt);
    //    if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
823
824
825
826
827
828

    // No Sts hits?
    int nmchits = mtra.StsHits.size();
    if (nmchits == 0) continue;

    double momentum = mtra.p;
Administrator's avatar
Administrator committed
829
830
    double pt       = sqrt(mtra.px * mtra.px + mtra.py * mtra.py);
    double theta    = acos(mtra.pz / mtra.p) * 180 / 3.1415;
831
832
833
834
835
836

    h_mcp->Fill(momentum);
    h_nmchits->Fill(nmchits);

    int nSta = mtra.NStations();

Administrator's avatar
Administrator committed
837
838
    CbmL1HitStore& fh = vHitStore[*(mtra.StsHits.begin())];
    CbmL1HitStore& lh = vHitStore[*(mtra.StsHits.rbegin())];
839
840

    h_reg_MCmom->Fill(momentum);
Administrator's avatar
Administrator committed
841
    if (mtra.IsPrimary()) {
842
843
844
845
      h_reg_mom_prim->Fill(momentum);
      h_reg_prim_MCmom->Fill(momentum);
      h_reg_nhits_prim->Fill(nSta);
      h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
846
      h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
847
848
849
      if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
    }
    else {
850
851
852
      h_reg_mom_sec->Fill(momentum);
      h_reg_sec_MCmom->Fill(momentum);
      h_reg_nhits_sec->Fill(nSta);
Administrator's avatar
Administrator committed
853
      if (momentum > 0.01) {
854
        h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
855
        h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
856
        if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
857
858
859
      }
    }

Administrator's avatar
Administrator committed
860
    if (mtra.IsAdditional()) { h_acc_mom_short123s->Fill(momentum); }
861

Administrator's avatar
Administrator committed
862
    if (!mtra.IsReconstructable()) continue;
863
864
865
    mc_total++;

    h_acc_MCmom->Fill(momentum);
Administrator's avatar
Administrator committed
866
    if (mtra.IsPrimary()) {
867
868
869
870
      h_acc_mom_prim->Fill(momentum);
      h_acc_prim_MCmom->Fill(momentum);
      h_acc_nhits_prim->Fill(nSta);
      h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
871
      h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
872
873
874
      if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
    }
    else {
875
876
877
      h_acc_mom_sec->Fill(momentum);
      h_acc_sec_MCmom->Fill(momentum);
      h_acc_nhits_sec->Fill(nSta);
Administrator's avatar
Administrator committed
878
      if (momentum > 0.01) {
879
        h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
880
        h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
881
        if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
882
883
884
885
      }
    }


Administrator's avatar
Administrator committed
886
    // vertex distribution of primary and secondary tracks
887
888
    h2_vertex->Fill(mtra.z, mtra.y);
    //h3_vertex->Fill(mtra.z, mtra.x, mtra.y);
Administrator's avatar
Administrator committed
889
    if (mtra.mother_ID < 0) {  // primary
890
891
      h2_prim_vertex->Fill(mtra.z, mtra.y);
      //h3_prim_vertex->Fill(mtra.z, mtra.x, mtra.y);
892
893
    }
    else {
894
895
896
897
898
      h2_sec_vertex->Fill(mtra.z, mtra.y);
      //h3_sec_vertex->Fill(mtra.z, mtra.x, mtra.y);
    }


Administrator's avatar
Administrator committed
899
900
    int iph           = mtra.StsHits[0];
    CbmL1HitStore& ph = vHitStore[iph];
901

Administrator's avatar
Administrator committed
902
903
904
905
    h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
    if (fabs(mtra.pz) > 1.e-8) {
      h_tx->Fill(mtra.px / mtra.pz);
      h_ty->Fill(mtra.py / mtra.pz);
906
907
    }

Administrator's avatar
Administrator committed
908
909
    // reconstructed tracks
    bool reco   = (mtra.IsReconstructed());
910
911
    int nMCHits = mtra.NStations();

Administrator's avatar
Administrator committed
912
    if (reco) {
913
914
915
916
      p_eff_all_vs_mom->Fill(momentum, 100.0);
      p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
      p_eff_all_vs_pt->Fill(pt, 100.0);
      h_reco_MCmom->Fill(momentum);
Administrator's avatar
Administrator committed
917
      if (mtra.mother_ID < 0) {  // primary
918
919
920
921
        p_eff_prim_vs_mom->Fill(momentum, 100.0);
        p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
        p_eff_prim_vs_pt->Fill(pt, 100.0);
        p_eff_prim_vs_theta->Fill(theta, 100.0);
922
923
      }
      else {
924
925
926
        p_eff_sec_vs_mom->Fill(momentum, 100.0);
        p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
      }
Administrator's avatar
Administrator committed
927
      if (mtra.mother_ID < 0) {  // primary
928
929
930
931
        h_reco_mom_prim->Fill(momentum);
        h_reco_prim_MCmom->Fill(momentum);
        h_reco_nhits_prim->Fill(nSta);
        h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
932
        h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
933
934
935
        if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
      }
      else {
936
937
938
        h_reco_mom_sec->Fill(momentum);
        h_reco_sec_MCmom->Fill(momentum);
        h_reco_nhits_sec->Fill(nSta);
Administrator's avatar
Administrator committed
939
        if (momentum > 0.01) {
940
          h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
941
          h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
942
          if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
943
944
        }
      }
945
946
    }
    else {
947
948
949
950
      h_notfound_mom->Fill(momentum);
      p_eff_all_vs_mom->Fill(momentum, 0.0);
      p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
      p_eff_all_vs_pt->Fill(pt, 0.0);
Administrator's avatar
Administrator committed
951
      if (mtra.mother_ID < 0) {  // primary
952
953
954
955
        p_eff_prim_vs_mom->Fill(momentum, 0.0);
        p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
        p_eff_prim_vs_pt->Fill(pt, 0.0);
        p_eff_prim_vs_theta->Fill(theta, 0.0);
956
957
      }
      else {
958
959
960
        p_eff_sec_vs_mom->Fill(momentum, 0.0);
        p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
      }
Administrator's avatar
Administrator committed
961
      if (mtra.mother_ID < 0) {  // primary
962
963
964
        h_rest_mom_prim->Fill(momentum);
        h_rest_nhits_prim->Fill(nSta);
        h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
965
        h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
966
967
968
        if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
      }
      else {
969
970
        h_rest_mom_sec->Fill(momentum);
        h_rest_nhits_sec->Fill(nSta);
Administrator's avatar
Administrator committed
971
        if (momentum > 0.01) {
972
          h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
Administrator's avatar
Administrator committed
973
          h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
974
          if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
975
976
977
978
        }
      }
    }

Administrator's avatar
Administrator committed
979
980
981
    // killed tracks. At least one hit of they belong to some recoTrack
    //     bool killed = 0;
    if (!reco) {
982
983
      h_notfound_nhits->Fill(nmchits);
      h_notfound_station->Fill(ph.iStation);
Administrator's avatar
Administrator committed
984
      h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
985

Administrator's avatar
Administrator committed
986
987
988
989
      if (mtra.Points.size() > 0) {
        CbmL1MCPoint& pMC = vMCPoints[mtra.Points[0]];
        h_notfound_tx->Fill(pMC.px / pMC.pz);
        h_notfound_ty->Fill(pMC.py / pMC.pz);
990
991
      }

Administrator's avatar
Administrator committed
992
993
      //      CbmL1HitStore &ph21 = vHitStore[mtra.StsHits[0]];
      //      CbmL1HitStore &ph22 = vHitStore[mtra.StsHits[1]];
994

995
996
      //      double z21 = algo->GetStations()[ph21.iStation].z[0];
      //      double z22 = algo->GetStations()[ph22.iStation].z[0];
Administrator's avatar
Administrator committed
997
998
999
1000
      //      if( fabs(z22-z21)>1.e-4 ){
      //        h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
      //        h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
      //      }