CbmL1.cxx 86 KB
Newer Older
1
2
/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
Valentina Akishina's avatar
Valentina Akishina committed
3
   Authors: Ivan Kisel, Sergey Gorbunov [committer], Valentina Akishina, Igor Kulakov, Maksym Zyzak */
4

5
6
7
/*
 *====================================================================
 *
8
9
 *  CBM Level 1 Reconstruction
 *
10
11
 *  Authors: I.Kisel,  S.Gorbunov
 *
12
 *  e-mail : ikisel@kip.uni-heidelberg.de
13
14
15
16
17
18
19
 *
 *====================================================================
 *
 *  L1 main program
 *
 *====================================================================
 */
Sergei Zharko's avatar
update    
Sergei Zharko committed
20

21
22
#include "CbmL1.h"

23
#include "CbmKF.h"
24
#include "CbmKFVertex.h"
Administrator's avatar
Administrator committed
25
#include "CbmL1PFFitter.h"
26
#include "CbmMCDataManager.h"
27
28
29
30
31
32
#include "CbmMuchGeoScheme.h"
#include "CbmMuchModuleGem.h"
#include "CbmMuchPad.h"
#include "CbmMuchStation.h"
#include "CbmMvdDetector.h"
#include "CbmMvdStationPar.h"
33
#include "CbmSetup.h"
Sergei Zharko's avatar
Sergei Zharko committed
34

Sergei Zharko's avatar
updates    
Sergei Zharko committed
35
#include <boost/filesystem.hpp>
36
37
// TODO: include of CbmSetup.h creates problems on Mac
// #include "CbmSetup.h"
38
#include "CbmMCDataObject.h"
Administrator's avatar
Administrator committed
39
#include "CbmStsFindTracks.h"
40
#include "CbmStsHit.h"
Administrator's avatar
Administrator committed
41
#include "CbmStsParSetModule.h"
42
43
#include "CbmStsParSetSensor.h"
#include "CbmStsParSetSensorCond.h"
Administrator's avatar
Administrator committed
44
45
#include "CbmStsSetup.h"
#include "CbmStsStation.h"
46
47
48
#include "CbmTofCell.h"
#include "CbmTofDigiBdfPar.h"
#include "CbmTofDigiPar.h"     // in tof/TofParam
Administrator's avatar
Administrator committed
49
50
#include "CbmTrdParModDigi.h"  // for CbmTrdModule
#include "CbmTrdParSetDigi.h"  // for CbmTrdParSetDigi
51

52
53
#include "FairEventHeader.h"
#include "FairRunAna.h"
54
55
56

#include "TGeoArb8.h"
#include "TGeoBoolNode.h"
Administrator's avatar
Administrator committed
57
#include "TGeoCompositeShape.h"
58
59
60
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMatrixD.h"
61
#include "TProfile2D.h"
62
63
#include "TROOT.h"
#include "TRandom3.h"
Administrator's avatar
Administrator committed
64
65
#include "TVector3.h"
#include "TVectorD.h"
66
67
68
69
70
#include <TFile.h>

#include <fstream>
#include <iomanip>
#include <iostream>
71
#include <list>
72

73
#include "L1Algo/L1Algo.h"
Sergei Zharko's avatar
Sergei Zharko committed
74
#include "L1Algo/L1Assert.h"
75
76
#include "L1Algo/L1Branch.h"
#include "L1Algo/L1Field.h"
Sergey Gorbunov's avatar
Sergey Gorbunov committed
77
#include "L1Algo/L1Hit.h"
78
79
80
#include "L1AlgoInputData.h"
#include "L1Event.h"

81
82
83
84
85
86
87
88
using std::cout;
using std::endl;
using std::ios;

#include "KFTopoPerformance.h"

ClassImp(CbmL1)

Sergei Zharko's avatar
Sergei Zharko committed
89
  static L1Algo algo_static _fvecalignment;  // TODO: gAlgo
90
91
92

//L1AlgoInputData* fData_static _fvecalignment;

Administrator's avatar
Administrator committed
93
94
95
CbmL1* CbmL1::fInstance = 0;


Sergey Gorbunov's avatar
Sergey Gorbunov committed
96
CbmL1::CbmL1() : FairTask("L1")
97
{
Administrator's avatar
Administrator committed
98
  if (!fInstance) fInstance = this;
Sergei Zharko's avatar
Sergei Zharko committed
99
  if (!fpInitManager) { fpInitManager = algo_static.GetInitManager(); }
100
101
}

102
CbmL1::CbmL1(const char* name, Int_t iVerbose, Int_t _fPerformance, int fSTAPDataMode_, TString fSTAPDataDir_,
Administrator's avatar
Administrator committed
103
104
             int findParticleMode_)
  : FairTask(name, iVerbose)
Sergey Gorbunov's avatar
Sergey Gorbunov committed
105
  , fPerformance(_fPerformance)
Administrator's avatar
Administrator committed
106
107
108
  , fSTAPDataMode(fSTAPDataMode_)
  , fSTAPDataDir(fSTAPDataDir_)
  , fFindParticlesMode(findParticleMode_)
109
{
Administrator's avatar
Administrator committed
110
  if (!fInstance) fInstance = this;
Sergei Zharko's avatar
Sergei Zharko committed
111
  if (!fpInitManager) { fpInitManager = algo_static.GetInitManager(); }
112
113
}

114
115
CbmL1::~CbmL1()
{
Sergey Gorbunov's avatar
Sergey Gorbunov committed
116
  if (fInstance == this) fInstance = nullptr;
117
118
}

119
120
void CbmL1::CheckDetectorPresence()
{
121
122
123
124
125
126
127
  Bool_t IsMuch = CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch);
  Bool_t IsTrd  = CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd);
  Bool_t IsTof  = CbmSetup::Instance()->IsActive(ECbmModuleId::kTof);
  //Bool_t IsSts        = CbmSetup::Instance()->IsActive(ECbmModuleId::kSts);
  Bool_t IsMvd = CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);

  /*
128
129
130
131
132
133
134
135
  TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();

  for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
    TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));

    if (TString(topNode->GetName()).Contains("much")) IsMuch = 1;
    if (TString(topNode->GetName()).Contains("trd")) IsTrd = 1;
    if (TString(topNode->GetName()).Contains("tof")) IsTof = 1;
136
    //if (TString(topNode->GetName()).Contains("sts")) IsSts = 1;
137
138
    if (TString(topNode->GetName()).Contains("mvd")) IsMvd = 1;
  }
139
  */
140
141
142
143
144
145
146

  fUseMUCH = (fUseMUCH && IsMuch);
  fUseTRD  = fUseTRD && IsTrd;
  fUseTOF  = fUseTOF && IsTof;
  fUseMVD  = fUseMVD && IsMvd;
}

147

148
149
void CbmL1::SetParContainers()
{
Administrator's avatar
Administrator committed
150
151
152
153
154
  FairRunAna* ana     = FairRunAna::Instance();
  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
  fDigiPar            = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
  //  fTrdDigiPar = (CbmTrdDigiPar*)(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdDigiPar"));
  fTrdDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
155
156
157
158

  // Trigger loading of STS parameters
  // If L1 runs standalone the parameters are not required by any STS task
  // but they are needed by CbmStsSetup
159
160
161
162
163
164
  fStsParSetModule     = dynamic_cast<CbmStsParSetModule*>(rtdb->getContainer("CbmStsParSetModule"));
  fStsParSetSensor     = dynamic_cast<CbmStsParSetSensor*>(rtdb->getContainer("CbmStsParSetSensor"));
  fStsParSetSensorCond = dynamic_cast<CbmStsParSetSensorCond*>(rtdb->getContainer("CbmStsParSetSensorCond"));

  fTofDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar"));
  rtdb->initContainers(ana->GetRunId());  // needed to get tracking stations for ToF hits
165
166
167
}


168
169
InitStatus CbmL1::ReInit()
{
170
171
172
173
  SetParContainers();
  return Init();
}

174
175
InitStatus CbmL1::Init()
{
176
177
  fData = new L1AlgoInputData();

Administrator's avatar
Administrator committed
178
179
180
181
182
183
184
185
  if (fVerbose > 1) {
    char y[20] = " [0;33;44m";         // yellow
    char Y[20] = " [1;33;44m";         // yellow bold
    char W[20] = " [1;37;44m";         // white bold
    char o[20] = " [0m\n";             // color off
    Y[0] = y[0] = W[0] = o[0] = 0x1B;  // escape character

    cout << endl << endl;
186
187
188
189
190
191
192
193
194
195
196
197
198
    cout << "  " << W << "                                                                 " << o;
    cout << "  " << W << "  ===////======================================================  " << o;
    cout << "  " << W << "  =                                                           =  " << o;
    cout << "  " << W << "  =                   " << Y << "L1 on-line finder" << W << "                       =  " << o;
    cout << "  " << W << "  =                                                           =  " << o;
    cout << "  " << W << "  =  " << W << "Cellular Automaton 3.1 Vector" << y << " with " << W << "KF Quadro" << y
         << " technology" << W << "  =  " << o;
    cout << "  " << W << "  =                                                           =  " << o;
    cout << "  " << W << "  =  " << y << "Designed for CBM collaboration" << W << "                           =  " << o;
    cout << "  " << W << "  =  " << y << "All rights reserved" << W << "                                      =  " << o;
    cout << "  " << W << "  =                                                           =  " << o;
    cout << "  " << W << "  ========================================================////=  " << o;
    cout << "  " << W << "                                                                 " << o;
Administrator's avatar
Administrator committed
199
    cout << endl << endl;
200
201
  }

Sergey Gorbunov's avatar
Sergey Gorbunov committed
202
203
204
205
206
207
208
209
  if (fVerbose > 1) {
#ifdef _OPENMP
    LOG(info) << "L1 is running with OpenMP" << endl;
#else
    LOG(info) << "L1 is running without OpenMP" << endl;
#endif
  }

Administrator's avatar
Administrator committed
210
  FairRootManager* fManger = FairRootManager::Instance();
211

Administrator's avatar
Administrator committed
212
  FairRunAna* Run = FairRunAna::Instance();
213
  {
214
215
    fUseMVD                    = 1;
    CbmStsFindTracks* FindTask = L1_DYNAMIC_CAST<CbmStsFindTracks*>(Run->GetTask("STSFindTracks"));
Administrator's avatar
Administrator committed
216
    if (FindTask) fUseMVD = FindTask->MvdUsage();
217
218
219
220
    // TODO: include of CbmSetup.h creates problems on Mac
    // if (!CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd)) { fUseMVD = false; }
    // N Mvd stations is read from the KF material
    if (CbmKF::Instance()->vMvdMaterial.size() == 0) { fUseMVD = false; }
221
222
  }

Sergey Gorbunov's avatar
cleanup    
Sergey Gorbunov committed
223
  fHistoDir = gROOT->mkdir("L1");
224

Administrator's avatar
Administrator committed
225

226
  // turn on reconstruction in sub-detectors
Administrator's avatar
Administrator committed
227

228
  fUseMUCH = 0;
Administrator's avatar
Administrator committed
229
230
231
  fUseTRD  = 0;
  fUseTOF  = 0;

232
  if (fTrackingMode == L1Algo::TrackingMode::kMcbm) {
233
    fUseMUCH = 1;
234
235
236
    fUseTRD  = 1;
    fUseTOF  = 1;
  }
Administrator's avatar
Administrator committed
237

238
  if (fTrackingMode == L1Algo::TrackingMode::kGlobal) {
239
    fUseMUCH = 0;
240
    fUseTRD  = 1;
241
    fUseTOF  = 0;
242
  }
243

244
245
  CheckDetectorPresence();

Administrator's avatar
Administrator committed
246
247
248
249
250
251
  fStsPoints  = 0;
  fMvdPoints  = 0;
  fMuchPoints = 0;
  fTrdPoints  = 0;
  fTofPoints  = 0;
  fMCTracks   = 0;
252

253
254
255
256
257
  listMvdHitMatches  = 0;
  fTrdHitMatches     = 0;
  listMuchHitMatches = 0;
  fTofHitDigiMatches = 0;

258
  listStsClusters = 0;
259

260
261
  vFileEvent.clear();

262
  if (!fLegacyEventMode) {  //  Time-slice mode selected
263
264
265
    LOG(info) << GetName() << ": running in time-slice mode.";
    fTimeSlice = NULL;
    fTimeSlice = (CbmTimeSlice*) fManger->GetObject("TimeSlice.");
266
    if (fTimeSlice == NULL) LOG(fatal) << GetName() << ": No time slice branch in the tree!";
267

Administrator's avatar
Administrator committed
268
  }  //? time-slice mode
269

Administrator's avatar
Administrator committed
270
  else  // event mode
271
272
273
    LOG(info) << GetName() << ": running in event mode.";


274
275
276
  listStsClusters     = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsCluster"));
  listStsHitMatch     = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHitMatch"));
  listStsClusterMatch = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsClusterMatch"));
Administrator's avatar
Administrator committed
277
278
279
280
281
282
283
284
285
286
287
288

  listStsHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHit"));

  if (!fUseMUCH) {
    fMuchPixelHits = 0;

    fDigisMuch       = 0;
    fDigiMatchesMuch = 0;
    fClustersMuch    = 0;

    fMuchPoints        = 0;
    listMuchHitMatches = 0;
289
290
  }
  else {
Administrator's avatar
Administrator committed
291
292
293


    fMuchPixelHits = (TClonesArray*) fManger->GetObject("MuchPixelHit");
294
  }
Administrator's avatar
Administrator committed
295
296
297

  if (!fUseTRD) {
    fTrdPoints     = 0;
298
    fTrdHitMatches = 0;
Administrator's avatar
Administrator committed
299
300
    fTrdPoints     = 0;
    listTrdHits    = 0;
301
302
  }
  else {
Administrator's avatar
Administrator committed
303
304

    listTrdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("TrdHit"));
305
  }
Administrator's avatar
Administrator committed
306
307
308
309
310

  if (!fUseTOF) {
    fTofPoints         = 0;
    fTofHitDigiMatches = 0;
    fTofHits           = 0;
311
312
  }
  else {
Administrator's avatar
Administrator committed
313
    fTofHits = (TClonesArray*) fManger->GetObject("TofHit");
314
315
  }

Administrator's avatar
Administrator committed
316
  if (fPerformance) {
317
    CbmMCDataManager* mcManager = (CbmMCDataManager*) fManger->GetObject("MCDataManager");
Administrator's avatar
Administrator committed
318
319
    if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";

320
    fStsPoints = mcManager->InitBranch("StsPoint");
321

322
323
    fMcEventHeader = mcManager->GetObject("MCEventHeader.");

324
    fMCTracks = mcManager->InitBranch("MCTrack");
325

Administrator's avatar
Administrator committed
326
327
    if (NULL == fStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
    if (NULL == fMCTracks) LOG(fatal) << GetName() << ": No MCTrack data!";
328
    if (NULL == fMcEventHeader) LOG(fatal) << GetName() << ": No MC event header data!";
Administrator's avatar
Administrator committed
329

330
    if (!fLegacyEventMode) {
Administrator's avatar
Administrator committed
331
      fEventList = (CbmMCEventList*) fManger->GetObject("MCEventList.");
332
      if (NULL == fEventList) LOG(fatal) << GetName() << ": No MCEventList data!";
Administrator's avatar
Administrator committed
333
334
    }

335
336
    if (fUseMVD) {
      fMvdPoints         = mcManager->InitBranch("MvdPoint");
337
338
      listMvdDigiMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdDigiMatch"));
      listMvdHitMatches  = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHitMatch"));
339
      if (!listMvdHitMatches) { LOG(error) << "No listMvdHitMatches provided, performance is not done correctly"; }
340
    }
Administrator's avatar
Administrator committed
341
342
343
344
345

    if (!fUseTRD) {
      fTrdPoints     = 0;
      fTrdHitMatches = 0;
      fTrdPoints     = 0;
346
347
    }
    else {
348
      fTrdHitMatches = (TClonesArray*) fManger->GetObject("TrdHitMatch");
Administrator's avatar
Administrator committed
349
      fTrdPoints     = mcManager->InitBranch("TrdPoint");
350
    }
Administrator's avatar
Administrator committed
351
352
353

    if (!fUseMUCH) {
      fMuchPoints        = 0;
354
      listMuchHitMatches = 0;
355
356
    }
    else {
357

358
359
360
361
362
      fDigisMuch         = (TClonesArray*) fManger->GetObject("MuchDigi");
      fDigiMatchesMuch   = (TClonesArray*) fManger->GetObject("MuchDigiMatch");
      fClustersMuch      = (TClonesArray*) fManger->GetObject("MuchCluster");
      fMuchPoints        = mcManager->InitBranch("MuchPoint");
      listMuchHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MuchPixelHitMatch"));
Administrator's avatar
Administrator committed
363
364
365
366
    }

    if (!fUseTOF) {
      fTofPoints         = 0;
367
      fTofHitDigiMatches = 0;
Administrator's avatar
Administrator committed
368
    }
369
370
    else {
      fTofPoints         = mcManager->InitBranch("TofPoint");
371
      fTofHitDigiMatches = static_cast<TClonesArray*>(fManger->GetObject("TofHitMatch"));
372
373
374
    }
  }
  else {
375
  }
376
377
  if (!fUseMVD) { listMvdHits = 0; }
  else {
Administrator's avatar
Administrator committed
378
    listMvdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHit"));
379
380
  }

Administrator's avatar
Administrator committed
381
382
  NMvdStations  = 0;
  NStsStations  = 0;
383
  NMuchStations = 0;
Administrator's avatar
Administrator committed
384
385
386
387
  NTrdStations  = 0;
  NTOFStation   = 0;
  NStation      = 0;

Sergei Zharko's avatar
Sergei Zharko committed
388
  // TODO: Replace algo initialization in the constructor (S.Zharko)
Administrator's avatar
Administrator committed
389
  algo = &algo_static;
390

Sergei Zharko's avatar
Sergei Zharko committed
391
  /**************************
392
393
   ** Field initialization **
   **************************/
Sergei Zharko's avatar
Sergei Zharko committed
394

395
396
397
398
  auto fieldGetterFcn = [](const double(&inPos)[3], double(&outB)[3]) {
    CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);
  };
  fpInitManager->SetFieldFunction(fieldGetterFcn);
399

Sergei Zharko's avatar
Sergei Zharko committed
400
  /***************************
401
402
   ** Target initialization **
   ***************************/
403

404
405
  auto& target = CbmKF::Instance()->vTargets[0];
  fpInitManager->SetTargetPosition(target.x, target.y, target.z);
Administrator's avatar
Administrator committed
406

Sergei Zharko's avatar
Sergei Zharko committed
407
  /*********************************
408
409
   ** Target field initialization **
   *********************************/
Administrator's avatar
Administrator committed
410

411
412
413
  fpInitManager->InitTargetField(2.5);


Sergei Zharko's avatar
Sergei Zharko committed
414
  /**************************************
415
416
417
418
419
420
   **                                  **
   ** STATIONS GEOMETRY INITIALIZATION **
   **                                  **
   **************************************/


Sergei Zharko's avatar
Sergei Zharko committed
421
  /***************************************************
422
423
   ** Active tracking detector subsystems selection **
   ***************************************************/
Administrator's avatar
Administrator committed
424

425
  fActiveTrackingDetectorIDs.insert(L1DetectorID::kSts);
Sergei Zharko's avatar
Sergei Zharko committed
426
  if (fUseMVD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMvd); }
427
  if (fUseMUCH) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kMuch); }
Sergei Zharko's avatar
Sergei Zharko committed
428
429
  if (fUseTRD) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTrd); }
  if (fUseTOF) { fActiveTrackingDetectorIDs.insert(L1DetectorID::kTof); }
430
431
  fpInitManager->SetActiveDetectorIDs(fActiveTrackingDetectorIDs);

Sergei Zharko's avatar
Sergei Zharko committed
432
  /*********************************************************************
433
434
435
436
   ** Counting numbers of stations for different detector subsystems  **
   *********************************************************************/

  CbmMuchGeoScheme* fGeoScheme = CbmMuchGeoScheme::Instance();
Sergei Zharko's avatar
Sergei Zharko committed
437

438
  /*** MuCh ***/
Administrator's avatar
Administrator committed
439
  if (fUseMUCH) {
440
441
442
443
    /// Save old global file and folder pointer to avoid messing with FairRoot
    TFile* oldFile     = gFile;
    TDirectory* oldDir = gDirectory;

Sergey Gorbunov's avatar
Sergey Gorbunov committed
444
    TFile* file         = new TFile(fMuchDigiFile, "READ");
Administrator's avatar
Administrator committed
445
446
    TObjArray* stations = (TObjArray*) file->Get("stations");
    fGeoScheme->Init(stations, 0);
447
448
449
450
451
    for (int iStation = 0; iStation < fGeoScheme->GetNStations(); iStation++) {
      const CbmMuchStation* station = fGeoScheme->GetStation(iStation);
      int nLayers                   = station->GetNLayers();
      NMuchStations += nLayers;
    }
452
453
454
455

    /// Restore old global file and folder pointer to avoid messing with FairRoot
    gFile      = oldFile;
    gDirectory = oldDir;
456
  }
Administrator's avatar
Administrator committed
457

458
  /*** TRD ***/
Administrator's avatar
Administrator committed
459
460
461
  if (fUseTRD) {
    Int_t layerCounter  = 0;
    TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
462
    for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
Administrator's avatar
Administrator committed
463
464
465
466
467
468
469
470
471
472
473
      TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
      if (TString(topNode->GetName()).Contains("trd")) {
        TObjArray* layers = topNode->GetNodes();
        for (Int_t iLayer = 0; iLayer < layers->GetEntriesFast(); iLayer++) {
          TGeoNode* layer = static_cast<TGeoNode*>(layers->At(iLayer));
          if (TString(layer->GetName()).Contains("layer")) { layerCounter++; }
        }
      }
    }
    NTrdStations = layerCounter;

474
475
    if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (fMissingHits)) { NTrdStations = NTrdStations - 1; }
  }
Administrator's avatar
Administrator committed
476

477
  /*** ToF ***/
478
479
  vector<float> TofStationZ;
  vector<int> TofStationN;
480
481
  if (fUseTOF) {
    NTOFStation = fTofDigiBdfPar->GetNbTrackingStations();
482
483
484
485
486
487
488
489
490
491
492
493
    TofStationZ.resize(NTOFStation, 0);
    TofStationN.resize(NTOFStation, 0);

    for (int iSmType = 0; iSmType < fTofDigiBdfPar->GetNbSmTypes(); iSmType++) {
      for (int iSm = 0; iSm < fTofDigiBdfPar->GetNbSm(iSmType); iSm++) {
        for (int iRpc = 0; iRpc < fTofDigiBdfPar->GetNbRpc(iSmType); iRpc++) {
          int iAddr                = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType);
          int station              = fTofDigiBdfPar->GetTrackingStation(iSmType, iSm, iRpc);
          CbmTofCell* fChannelInfo = fDigiPar->GetCell(iAddr);
          if (NULL == fChannelInfo) break;
          float z = fChannelInfo->GetZ();
          float x = fChannelInfo->GetX();
494
          //float y = fChannelInfo->GetY();
Valentina Akishina's avatar
Valentina Akishina committed
495
496
497

          if (station < 0) continue;

498
499
500
501
502
503
504
505
506
507
508
509
          if (fMissingHits) {
            if ((x > 20) && (z > 270) && (station == 1)) station = 2;
            if (z > 400) continue;
          }
          TofStationZ[station] += z;
          TofStationN[station] += 1;
        }
      }
    }
    for (Int_t i = 0; i < NTOFStation; i++)
      TofStationZ[i] = TofStationZ[i] / TofStationN[i];
  }
Administrator's avatar
Administrator committed
510

511
  /*** MVD and STS ***/
Administrator's avatar
Administrator committed
512
513
  CbmStsSetup* stsSetup = CbmStsSetup::Instance();
  if (!stsSetup->IsInit()) { stsSetup->Init(nullptr); }
514
515
516
  if (!stsSetup->IsModuleParsInit()) { stsSetup->SetModuleParameters(fStsParSetModule); }
  if (!stsSetup->IsSensorParsInit()) { stsSetup->SetSensorParameters(fStsParSetSensor); }
  if (!stsSetup->IsSensorCondInit()) { stsSetup->SetSensorConditions(fStsParSetSensorCond); }
517

Administrator's avatar
Administrator committed
518
  NMvdStations = (fUseMVD) ? CbmKF::Instance()->vMvdMaterial.size() : 0;
519
  NStsStations = CbmStsSetup::Instance()->GetNofStations();
520
  NStation     = NMvdStations + NStsStations + NMuchStations + NTrdStations + NTOFStation;
Sergei Zharko's avatar
Sergei Zharko committed
521

522
523
524
525
526
527
  // Provide crosscheck number of stations for the fpInitManager
  fpInitManager->SetStationsNumberCrosscheck(L1DetectorID::kMvd, NMvdStations);
  fpInitManager->SetStationsNumberCrosscheck(L1DetectorID::kSts, NStsStations);
  fpInitManager->SetStationsNumberCrosscheck(L1DetectorID::kMuch, NMuchStations);
  fpInitManager->SetStationsNumberCrosscheck(L1DetectorID::kTrd, NTrdStations);
  fpInitManager->SetStationsNumberCrosscheck(L1DetectorID::kTof, NTOFStation);
Sergei Zharko's avatar
Sergei Zharko committed
528

529
  {
Administrator's avatar
Administrator committed
530
    if (fSTAPDataMode % 2 == 1) {  // 1,3
531
      LOG(warn) << "CbmL1::Init: geo vector was removed, currently data cannot be written to a text-file";
532
      // TODO: Rewrite parameters i/o into L1InitManager (S.Zharko, 12.05.2022)
533
      //WriteSTAPGeoData(geo);
534
535
536
537
538
539
540
541
    };
    //if(fSTAPDataMode >= 2){ // 2,3
    //  int ind2, ind = geo.size();
    //  ReadSTAPGeoData(geo, ind2);
    //  if (ind2 != ind)  LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful.";
    //};
  }

Administrator's avatar
Administrator committed
542
  if (fSTAPDataMode >= 2) {  // 2,3
543
544
    LOG(fatal) << "CbmL1::Init: geo vector was removed, currently data cannot be read from a text-file. "
               << "Please, run CbmL1 task with STAPDataMode option < 2";
545
    // TODO: Rewrite parameters i/o into L1InitManager (S.Zharko, 12.05.2022)
546
547
548
549
550
    //int ind2, ind = geo.size();
    //ReadSTAPGeoData(geo, ind2);
    //if (ind2 != ind)
    //  LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful.";
  }
Sergei Zharko's avatar
Sergei Zharko committed
551
552

  /***************************************
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
   ** Stations geometry initialization  **
   ***************************************/

  /*** MVD stations info ***/
  for (int iSt = 0; iSt < NMvdStations; ++iSt) {  // NOTE: example using in-stack defined objects
    CbmMvdDetector* mvdDetector     = CbmMvdDetector::Instance();
    CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();

    CbmKFTube& t     = CbmKF::Instance()->vMvdMaterial[iSt];
    auto stationInfo = L1BaseStationInfo(L1DetectorID::kMvd, iSt);
    stationInfo.SetStationType(1);
    // MVD // TODO: to be exchanged with specific flags (timeInfo, fieldInfo etc.) (S.Zh.)
    stationInfo.SetTimeInfo(0);
    stationInfo.SetTimeResolution(1000.);
    stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
    stationInfo.SetZ(t.z);
    auto thickness = t.dz;
    auto radLength = t.RadLength;
    stationInfo.SetMaterial(thickness, radLength);
    stationInfo.SetXmax(t.R);
    stationInfo.SetYmax(t.R);
    stationInfo.SetRmin(t.r);
    stationInfo.SetRmax(t.R);
    fscal mvdFrontPhi   = 0;
    fscal mvdBackPhi    = TMath::Pi() / 2.;
    fscal mvdFrontSigma = mvdStationPar->GetXRes(iSt) / 10000;
    fscal mvdBackSigma  = mvdStationPar->GetYRes(iSt) / 10000;
    stationInfo.SetFrontBackStripsGeometry(mvdFrontPhi, mvdFrontSigma, mvdBackPhi, mvdBackSigma);
    fpInitManager->AddStation(stationInfo);
    LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZdouble();
583
584
  }

585
586
  /*** STS stations info ***/
  for (int iSt = 0; iSt < NStsStations; ++iSt) {  // NOTE: example using smart pointers
Sergei Zharko's avatar
Sergei Zharko committed
587
    auto cbmSts      = CbmStsSetup::Instance()->GetStation(iSt);
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
    auto stationInfo = L1BaseStationInfo(L1DetectorID::kSts, iSt);
    stationInfo.SetStationType(0);  // STS
    stationInfo.SetTimeInfo(1);
    stationInfo.SetTimeResolution(5.);
    stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
    // Setup station geometry and material
    stationInfo.SetZ(cbmSts->GetZ());
    double stsXmax = cbmSts->GetXmax();
    double stsYmax = cbmSts->GetYmax();
    stationInfo.SetXmax(stsXmax);
    stationInfo.SetYmax(stsYmax);
    stationInfo.SetRmin(0);
    stationInfo.SetRmax(stsXmax > stsYmax ? stsXmax : stsYmax);
    auto thickness = cbmSts->GetSensorD();
    auto radLength = cbmSts->GetRadLength();
    stationInfo.SetMaterial(thickness, radLength);
    // Setup strips geometry
    fscal stsFrontPhi   = cbmSts->GetSensorRotation() + cbmSts->GetSensorStereoAngle(0) * TMath::Pi() / 180.;
    fscal stsBackPhi    = cbmSts->GetSensorRotation() + cbmSts->GetSensorStereoAngle(1) * TMath::Pi() / 180.;
    fscal stsFrontSigma = cbmSts->GetSensorPitch(0) / sqrt(12);
    fscal stsBackSigma  = stsFrontSigma;
    stationInfo.SetFrontBackStripsGeometry(stsFrontPhi, stsFrontSigma, stsBackPhi, stsBackSigma);
    fpInitManager->AddStation(stationInfo);
    LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZdouble();
  }
Sergei Zharko's avatar
Sergei Zharko committed
613

614
615
616
617
618
619
620
621
622
623
624
625
626
627
  /*** MuCh stations info ***/
  for (int iSt = 0; iSt < NMuchStations; ++iSt) {
    int muchStationID       = iSt / 3;
    int muchLayerID         = iSt % 3;
    CbmMuchStation* station = (CbmMuchStation*) fGeoScheme->GetStation(muchStationID);
    CbmMuchLayer* layer     = station->GetLayer(muchLayerID);

    auto stationInfo = L1BaseStationInfo(L1DetectorID::kMuch, iSt);
    stationInfo.SetStationType(2);
    stationInfo.SetTimeInfo(1);
    stationInfo.SetTimeResolution(3.9);
    stationInfo.SetFieldStatus(0);
    stationInfo.SetZ(layer->GetZ());
    auto thickness = layer->GetDz();
Sergei Zharko's avatar
Sergei Zharko committed
628
    auto radLength = 0.;  // Why 0??? (S.Zharko)
629
630
631
632
633
634
635
636
637
638
639
640
641
    stationInfo.SetMaterial(thickness, radLength);
    stationInfo.SetXmax(100.);
    stationInfo.SetYmax(100.);
    stationInfo.SetRmin(10.);
    stationInfo.SetRmax(100.);
    fscal muchFrontPhi   = 0;
    fscal muchBackPhi    = TMath::Pi() / 2.;
    fscal muchFrontSigma = 0.35;
    fscal muchBackSigma  = 0.35;
    stationInfo.SetFrontBackStripsGeometry(muchFrontPhi, muchFrontSigma, muchBackPhi, muchBackSigma);
    fpInitManager->AddStation(stationInfo);
    LOG(info) << "- MuCh station " << iSt << " at z = " << stationInfo.GetZdouble();
  }
642

643
644
645
646
647
  /*** TRD stations info ***/
  for (int iSt = 0; iSt < NTrdStations; ++iSt) {
    int skip = iSt;  // temporary solution to remove TRD with id == 1 wrom mCBM
    if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (fMissingHits)) {
      if (iSt > 0) { skip++; }
Sergei Zharko's avatar
Sergei Zharko committed
648
    }
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
    int trdModuleID          = fTrdDigiPar->GetModuleId(skip);
    CbmTrdParModDigi* module = (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(trdModuleID);
    auto stationInfo         = L1BaseStationInfo(L1DetectorID::kTrd, skip);
    int stationType          = (iSt == 1 || iSt == 3) ? 6 : 3;
    stationInfo.SetStationType(stationType);
    stationInfo.SetTimeInfo(1);
    stationInfo.SetTimeResolution(10.);
    stationInfo.SetFieldStatus(0);
    stationInfo.SetZ(module->GetZ());
    auto thickness = 2. * module->GetSizeZ();
    auto radLength = 1.6;
    stationInfo.SetMaterial(thickness, radLength);
    stationInfo.SetXmax(module->GetSizeX());
    stationInfo.SetYmax(module->GetSizeY());
    stationInfo.SetRmin(0.);
    stationInfo.SetRmax(2. * module->GetSizeX());  // TODO: Why multiplied with 2.?
    fscal trdFrontPhi   = 0;
    fscal trdBackPhi    = TMath::Pi() / 2.;
    fscal trdFrontSigma = 0.15;
    fscal trdBackSigma  = 0.15;
    stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdFrontSigma, trdBackPhi, trdBackSigma);
    fpInitManager->AddStation(stationInfo);
    LOG(info) << "- TRD station " << iSt << " at z = " << stationInfo.GetZdouble();
  }
Sergei Zharko's avatar
Sergei Zharko committed
673

674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
  /*** TOF stations info ***/
  for (int iSt = 0; iSt < NTOFStation; ++iSt) {
    auto stationInfo = L1BaseStationInfo(L1DetectorID::kTof, iSt);
    stationInfo.SetStationType(4);
    stationInfo.SetTimeInfo(1);
    stationInfo.SetTimeResolution(0.075);
    stationInfo.SetFieldStatus(0);
    stationInfo.SetZ(TofStationZ[iSt]);
    auto thickness = 10.;
    auto radLength = 2.;
    stationInfo.SetMaterial(thickness, radLength);
    stationInfo.SetXmax(20.);
    stationInfo.SetYmax(20.);
    stationInfo.SetRmin(0.);
    stationInfo.SetRmax(150.);
    fscal tofFrontPhi   = 0;
    fscal tofBackPhi    = TMath::Pi() / 2.;
    fscal tofFrontSigma = 0.42;
    fscal tofBackSigma  = 0.23;
    stationInfo.SetFrontBackStripsGeometry(tofFrontPhi, tofFrontSigma, tofBackPhi, tofBackSigma);
    fpInitManager->AddStation(stationInfo);
    LOG(info) << "- TOF station " << iSt << " at z = " << stationInfo.GetZdouble();
  }
Sergei Zharko's avatar
Sergei Zharko committed
697

Sergei Zharko's avatar
Sergei Zharko committed
698
  /****************************************
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
   **                                    **
   ** TRACKING ITERATIONS INITIALIZATION **
   **                                    **
   ****************************************/

  // TODO: Need to provide a selection: default iterations input (these hardcoded ones), input from file or input
  //       from runing macro (S.Zharko)
  auto trackingIterFastPrim = L1CAIteration("FastPrimIter");
  trackingIterFastPrim.SetTrackChi2Cut(10.f);
  trackingIterFastPrim.SetTripletChi2Cut(23.4450f);  // = 7.815 * 3;  // prob = 0.05
  trackingIterFastPrim.SetDoubletChi2Cut(7.56327f);  // = 1.3449 * 2.f / 3.f;  // prob = 0.1
  trackingIterFastPrim.SetPickGather(3.0f);
  trackingIterFastPrim.SetPickNeighbour(5.0f);
  trackingIterFastPrim.SetMaxInvMom(1.0 / 0.5);
  trackingIterFastPrim.SetMaxSlopePV(1.1f);
  trackingIterFastPrim.SetMaxSlope(2.748f);
  trackingIterFastPrim.SetMaxDZ(0);
  trackingIterFastPrim.SetTargetPosSigmaXY(1, 1);
  trackingIterFastPrim.SetMinLevelTripletStart(0);
  trackingIterFastPrim.SetPrimary(true);

  auto trackingIterAllPrim = L1CAIteration("AllPrimIter");
  trackingIterAllPrim.SetTrackChi2Cut(10.f);
  trackingIterAllPrim.SetTripletChi2Cut(23.4450f);
  trackingIterAllPrim.SetDoubletChi2Cut(7.56327f);
  trackingIterAllPrim.SetPickGather(4.0f);
  trackingIterAllPrim.SetPickNeighbour(5.0f);
  trackingIterAllPrim.SetMaxInvMom(1.0 / 0.05);
  trackingIterAllPrim.SetMaxSlopePV(1.1f);
  trackingIterAllPrim.SetMaxSlope(2.748f);
  trackingIterAllPrim.SetMaxDZ(0.1);
  trackingIterAllPrim.SetTargetPosSigmaXY(1, 1);
  trackingIterAllPrim.SetMinLevelTripletStart(0);
  trackingIterAllPrim.SetPrimary(true);

  auto trackingIterFastPrim2 = L1CAIteration("FastPrim2Iter");
  trackingIterFastPrim2.SetTrackChi2Cut(10.f);
  trackingIterFastPrim2.SetTripletChi2Cut(21.1075f);
  trackingIterFastPrim2.SetDoubletChi2Cut(7.56327f);
  trackingIterFastPrim2.SetPickGather(3.0f);
  trackingIterFastPrim2.SetPickNeighbour(5.0f);
  trackingIterFastPrim2.SetMaxInvMom(1.0 / 0.5);
  trackingIterFastPrim2.SetMaxSlopePV(1.1f);
  trackingIterFastPrim2.SetMaxSlope(2.748f);
  trackingIterFastPrim2.SetMaxDZ(0);
  trackingIterFastPrim2.SetTargetPosSigmaXY(5, 5);
  trackingIterFastPrim2.SetMinLevelTripletStart(0);
  trackingIterFastPrim2.SetPrimary(true);

  auto trackingIterAllSec = L1CAIteration("AllSecIter");
  trackingIterAllSec.SetTrackChi2Cut(10.f);
  trackingIterAllSec.SetTripletChi2Cut(18.7560f);  // = 6.252 * 3;  // prob = 0.1
  trackingIterAllSec.SetDoubletChi2Cut(7.56327f);
  trackingIterAllSec.SetPickGather(4.0f);
  trackingIterAllSec.SetPickNeighbour(5.0f);
  trackingIterAllSec.SetMaxInvMom(1.0 / 0.1);
  trackingIterAllSec.SetMaxSlopePV(1.5f);
  trackingIterAllSec.SetMaxSlope(2.748f);
  trackingIterAllSec.SetMaxDZ(0.1);
  trackingIterAllSec.SetTargetPosSigmaXY(10, 10);
  trackingIterAllSec.SetMinLevelTripletStart(1);
  trackingIterAllSec.SetPrimary(false);

  auto trackingIterFastPrimJump = L1CAIteration("FastPrimJumpIter");
  trackingIterFastPrimJump.SetTrackChi2Cut(10.f);
  trackingIterFastPrimJump.SetTripletChi2Cut(21.1075f);  // prob = 0.01
  trackingIterFastPrimJump.SetDoubletChi2Cut(7.56327f);
  trackingIterFastPrimJump.SetPickGather(3.0f);
  trackingIterFastPrimJump.SetPickNeighbour(5.0f);
  trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.5);
  trackingIterFastPrimJump.SetMaxSlopePV(1.1f);
  trackingIterFastPrimJump.SetMaxSlope(2.748f);
  trackingIterFastPrimJump.SetMaxDZ(0);
  trackingIterFastPrimJump.SetTargetPosSigmaXY(5, 5);
  trackingIterFastPrimJump.SetMinLevelTripletStart(0);
  trackingIterFastPrimJump.SetPrimary(true);

  auto trackingIterAllPrimJump = L1CAIteration("AllPrimJumpIter");
  trackingIterAllPrimJump.SetTrackChi2Cut(10.f);
  trackingIterAllPrimJump.SetTripletChi2Cut(18.7560f);
  trackingIterAllPrimJump.SetDoubletChi2Cut(7.56327f);
  trackingIterAllPrimJump.SetPickGather(4.0f);
  trackingIterAllPrimJump.SetPickNeighbour(5.0f);
  trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.1);
  trackingIterAllPrimJump.SetMaxSlopePV(1.1f);
  trackingIterAllPrimJump.SetMaxSlope(2.748f);
  trackingIterAllPrimJump.SetMaxDZ(0.1);
  trackingIterAllPrimJump.SetTargetPosSigmaXY(5, 5);
  trackingIterAllPrimJump.SetMinLevelTripletStart(0);
  trackingIterAllPrimJump.SetPrimary(true);

  auto trackingIterAllSecJump = L1CAIteration("AllSecJumpIter");
  trackingIterAllSecJump.SetTrackChi2Cut(10.f);
  trackingIterAllSecJump.SetTripletChi2Cut(21.1075f);
  trackingIterAllSecJump.SetDoubletChi2Cut(7.56327f);
  trackingIterAllSecJump.SetPickGather(4.0f);
  trackingIterAllSecJump.SetPickNeighbour(5.0f);
  trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.1);
  trackingIterAllSecJump.SetMaxSlopePV(1.5f);
  trackingIterAllSecJump.SetMaxSlope(2.748f);
  trackingIterAllSecJump.SetMaxDZ(0.1);
  trackingIterAllSecJump.SetMinLevelTripletStart(1);
  trackingIterAllSecJump.SetTargetPosSigmaXY(10, 10);

  auto trackingIterAllPrimE = L1CAIteration("AllPrimEIter");
  trackingIterAllPrimE.SetTrackChi2Cut(10.f);
  trackingIterAllPrimE.SetTripletChi2Cut(23.4450f);
  trackingIterAllPrimE.SetDoubletChi2Cut(7.56327f);
  trackingIterAllPrimE.SetPickGather(4.0f);
  trackingIterAllPrimE.SetPickNeighbour(5.0f);
  trackingIterAllPrimE.SetMaxInvMom(1.0 / 0.05);
  trackingIterAllPrimE.SetMaxSlopePV(1.1f);
  trackingIterAllPrimE.SetMaxSlope(2.748f);
  trackingIterAllPrimE.SetMaxDZ(0.1);
  trackingIterAllPrimE.SetTargetPosSigmaXY(1, 1);
  trackingIterAllPrimE.SetMinLevelTripletStart(0);
  trackingIterAllPrimE.SetPrimary(true);

  auto trackingIterAllSecE = L1CAIteration("AllSecEIter");
  trackingIterAllSecE.SetTrackChi2Cut(10.f);
  trackingIterAllSecE.SetTripletChi2Cut(18.7560f);
  trackingIterAllSecE.SetDoubletChi2Cut(7.56327f);
  trackingIterAllSecE.SetPickGather(4.0f);
  trackingIterAllSecE.SetPickNeighbour(5.0f);
  trackingIterAllSecE.SetMaxInvMom(1.0 / 0.05);
  trackingIterAllSecE.SetMaxSlopePV(1.5f);
  trackingIterAllSecE.SetMaxSlope(2.748f);
  trackingIterAllSecE.SetMaxDZ(0.1);
  trackingIterAllSecE.SetMinLevelTripletStart(1);
  trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);

  // Select default track finder
  if (fTrackingMode == L1Algo::TrackingMode::kMcbm) {
    trackingIterAllPrim.SetMaxInvMom(1. / 0.1);
    trackingIterAllPrimE.SetMaxInvMom(1. / 0.1);
    trackingIterAllSecE.SetMaxInvMom(1. / 0.1);

    trackingIterFastPrim.SetMaxInvMom(1.0 / 0.3);
    trackingIterAllSec.SetMaxInvMom(1.0 / 0.3);
    trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.3);
    trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.3);
    trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.3);

    fpInitManager->SetCAIterationsNumberCrosscheck(4);
    // Initialize CA track finder iterations sequence
    fpInitManager->PushBackCAIteration(trackingIterFastPrim);
    fpInitManager->PushBackCAIteration(trackingIterAllPrim);
    fpInitManager->PushBackCAIteration(trackingIterAllPrimJump);
    fpInitManager->PushBackCAIteration(trackingIterAllSec);
  }
  else {
    fpInitManager->SetCAIterationsNumberCrosscheck(4);
    // Initialize CA track finder iterations sequence
    fpInitManager->PushBackCAIteration(trackingIterFastPrim);
    fpInitManager->PushBackCAIteration(trackingIterAllPrim);
    fpInitManager->PushBackCAIteration(trackingIterAllPrimJump);
    fpInitManager->PushBackCAIteration(trackingIterAllSec);
    //fpInitManager->PushBackCAIteration(trackingIterAllPrimE);
    //fpInitManager->PushBackCAIteration(trackingIterAllSecE);
    //fpInitManager->PushBackCAIteration(trackingIterFastPrimJump);
    //fpInitManager->PushBackCAIteration(trackingIterFastPrim2);
    //fpInitManager->PushBackCAIteration(trackingIterAllSecJump);
  }
Sergei Zharko's avatar
Sergei Zharko committed
862

Sergei Zharko's avatar
Sergei Zharko committed
863
  /**********************
864
865
   ** Set special cuts **
   **********************/
Sergei Zharko's avatar
Sergei Zharko committed
866

867
868
869
  fpInitManager->SetGhostSuppression(fGhostSuppression);
  fpInitManager->SetTrackingLevel(fTrackingLevel);
  fpInitManager->SetMomentumCutOff(fMomentumCutOff);
870

871
  /**********************/
872

873
  algo->Init(fUseHitErrors, fTrackingMode, fMissingHits);
874

Sergei Zharko's avatar
Sergei Zharko committed
875
  /*********************************
876
877
878
879
   ** Material map initialization **
   *********************************/

  // TODO: Move it to the L1InitManager (S.Zharko, 15.05.2022)
880
  algo->fRadThick.reset(algo->GetNstations());
881

Administrator's avatar
Administrator committed
882
  // Read STS  MVD TRD MuCh ToF Radiation Thickness table
883
  TString stationName = "Radiation Thickness [%], Station";
Administrator's avatar
Administrator committed
884
885
  if (fUseMVD) {
    if (fMvdMatBudgetFileName != "") {
886
887
888
      /// Save old global file and folder pointer to avoid messing with FairRoot
      TFile* oldFile     = gFile;
      TDirectory* oldDir = gDirectory;
Valentina Akishina's avatar
Valentina Akishina committed
889
      TFile* rlFile      = new TFile(fMvdMatBudgetFileName);
890
      cout << "MVD Material budget file is " << fMvdMatBudgetFileName << ".\n";
891
      for (int j = 0, iSta = 0; iSta < algo->GetNstationsBeforePipe(); iSta++, j++) {
892
893
894
        TString stationNameMvd = stationName;
        stationNameMvd += j;
        TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameMvd);
Administrator's avatar
Administrator committed
895
        if (!hStaRadLen) {
896
          cout << "L1: incorrect " << fMvdMatBudgetFileName << " file. No " << stationNameMvd << "\n";
897
          exit(1);
Administrator's avatar
Administrator committed
898
        }
899
900
        const int NBins  = hStaRadLen->GetNbinsX();            // should be same in Y
        const float RMax = hStaRadLen->GetXaxis()->GetXmax();  // should be same as min
Administrator's avatar
Administrator committed
901
        algo->fRadThick[iSta].SetBins(NBins, RMax);
902
        algo->fRadThick[iSta].table.resize(NBins);
903
        double sumRF = 0., nRF = 0.;
Administrator's avatar
Administrator committed
904
        for (int iB = 0; iB < NBins; iB++) {
905
          algo->fRadThick[iSta].table[iB].resize(NBins);
Administrator's avatar
Administrator committed
906
          for (int iB2 = 0; iB2 < NBins; iB2++) {
907
            algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
908
            // Correction for holes in material map
909
            if (algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
Administrator's avatar
Administrator committed
910
              if (iB2 > 0 && iB2 < NBins - 1)
911
912
                algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
                                                                  0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
913
            // Correction for the incorrect harcoded value of RadThick of MVD stations
914
            if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = 0.0015;
915
            //              algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
916
917
            sumRF += algo->fRadThick[iSta].table[iB][iB2];
            nRF++;
918
919
          }
        }
920
        cout << " iSta " << iSta << " average RadThick in the material map " << sumRF / nRF << endl;
921
922
923
      }
      rlFile->Close();
      rlFile->Delete();
924
925
926
      /// Restore old global file and folder pointer to avoid messing with FairRoot
      gFile      = oldFile;
      gDirectory = oldDir;
927
928
    }
    else {
Administrator's avatar
Administrator committed
929
930
      LOG(warn) << "No MVD material budget file is found. Homogenious budget "
                   "will be used";
931
      for (int iSta = 0; iSta < algo->GetNstationsBeforePipe(); iSta++) {
Administrator's avatar
Administrator committed
932
933
934
935
        algo->fRadThick[iSta].SetBins(1,
                                      100);  // mvd should be in +-100 cm square
        algo->fRadThick[iSta].table.resize(1);
        algo->fRadThick[iSta].table[0].resize(1);
936
        algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
Administrator's avatar
Administrator committed
937
      }
938
939
940
    }
  }
  if (fStsMatBudgetFileName != "") {
941
942
943
    /// Save old global file and folder pointer to avoid messing with FairRoot
    TFile* oldFile     = gFile;
    TDirectory* oldDir = gDirectory;
Valentina Akishina's avatar
Valentina Akishina committed
944
    TFile* rlFile      = new TFile(fStsMatBudgetFileName);
945
    cout << "STS Material budget file is " << fStsMatBudgetFileName << ".\n";
Sergei Zharko's avatar
Sergei Zharko committed
946
947
    for (int j = 1, iSta = algo->GetNstationsBeforePipe(); iSta < (algo->GetNstationsBeforePipe() + NStsStations);
         iSta++, j++) {
948
949
950
      TString stationNameSts = stationName;
      stationNameSts += j;
      TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
Administrator's avatar
Administrator committed
951
      if (!hStaRadLen) {
952
        cout << "L1: incorrect " << fStsMatBudgetFileName << " file. No " << stationNameSts << "\n";
953
954
        exit(1);
      }
955
956
      const int NBins  = hStaRadLen->GetNbinsX();            // should be same in Y
      const float RMax = hStaRadLen->GetXaxis()->GetXmax();  // should be same as min
Administrator's avatar
Administrator committed
957
      algo->fRadThick[iSta].SetBins(NBins, RMax);
958
      algo->fRadThick[iSta].table.resize(NBins);
959
      double sumRF = 0., nRF = 0.;
Administrator's avatar
Administrator committed
960
      for (int iB = 0; iB < NBins; iB++) {
961
        algo->fRadThick[iSta].table[iB].resize(NBins);
Administrator's avatar
Administrator committed
962
        for (int iB2 = 0; iB2 < NBins; iB2++) {
963
          algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
964
965
          if (algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
            algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
Administrator's avatar
Administrator committed
966
          // cout << " iSta " << iSta << " iB " << iB << "iB2 " << iB2 << " RadThick " << algo->fRadThick[iSta].table[iB][iB2] << endl;
967
968
          sumRF += algo->fRadThick[iSta].table[iB][iB2];
          nRF++;
969
970
        }
      }
971
      cout << " iSta " << iSta << " average RadThick in the material map " << sumRF / nRF << endl;
972
973
974
    }
    rlFile->Close();
    rlFile->Delete();
975
976
977
    /// Restore old global file and folder pointer to avoid messing with FairRoot
    gFile      = oldFile;
    gDirectory = oldDir;
978
979
  }
  else {
Administrator's avatar
Administrator committed
980
981
    LOG(warn) << "No STS material budget file is found. Homogenious budget "
                 "will be used";
982
    for (int iSta = algo->GetNstationsBeforePipe(); iSta < (algo->GetNstationsBeforePipe() + NStsStations); iSta++) {
983
984
985
      algo->fRadThick[iSta].SetBins(1, 100);
      algo->fRadThick[iSta].table.resize(1);
      algo->fRadThick[iSta].table[0].resize(1);
986
      algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
987
988
    }
  }
Administrator's avatar
Administrator committed
989
990

  if (fUseMUCH)
991
    if (fMuchMatBudgetFileName != "") {
992
993
994
      /// Save old global file and folder pointer to avoid messing with FairRoot
      TFile* oldFile     = gFile;
      TDirectory* oldDir = gDirectory;
Valentina Akishina's avatar
Valentina Akishina committed
995
      TFile* rlFile      = new TFile(fMuchMatBudgetFileName);
996
997
      cout << "Much Material budget file is " << fMuchMatBudgetFileName << ".\n";
      for (int j = 1, iSta = (NStsStations + NMvdStations); iSta < (NStsStations + NMvdStations + NMuchStations);
Administrator's avatar
Administrator committed
998
999
1000
           iSta++, j++) {
        TString stationNameSts = stationName;
        stationNameSts += j;