CbmRoot
CbmL1.cxx
Go to the documentation of this file.
1 /*
2  *====================================================================
3  *
4  * CBM Level 1 Reconstruction
5  *
6  * Authors: I.Kisel, S.Gorbunov
7  *
8  * e-mail : ikisel@kip.uni-heidelberg.de
9  *
10  *====================================================================
11  *
12  * L1 main program
13  *
14  *====================================================================
15  */
16 #include "CbmL1.h"
17 
18 #include "CbmKFVertex.h"
19 #include "CbmL1PFFitter.h"
20 
21 #include "CbmMCDataManager.h"
22 #include "L1Algo/L1Algo.h"
23 #include "L1Algo/L1Branch.h"
24 #include "L1Algo/L1Field.h"
25 #include "L1Algo/L1StsHit.h"
26 #include "L1AlgoInputData.h"
27 
28 
29 #include "CbmKF.h"
30 #include "CbmStsFindTracks.h"
31 #include "CbmStsParSetModule.h"
32 #include "CbmStsParSetSensor.h"
33 #include "CbmStsParSetSensorCond.h"
34 #include "CbmStsSetup.h"
35 #include "CbmStsStation.h"
36 #include "FairEventHeader.h"
37 #include "FairRunAna.h"
38 
39 #include "CbmTrdParModDigi.h" // for CbmTrdModule
40 #include "CbmTrdParSetDigi.h" // for CbmTrdParSetDigi
41 
42 #include "CbmTofCell.h"
43 #include "CbmTofDigiPar.h" // in tof/TofParam
44 
45 #include "TGeoArb8.h"
46 #include "TGeoBoolNode.h"
47 #include "TGeoCompositeShape.h"
48 #include "TGeoManager.h"
49 
50 #include "CbmMuchGeoScheme.h"
51 #include "CbmMuchModuleGem.h"
52 #include "CbmMuchPad.h"
53 #include "CbmMuchStation.h"
54 #include "TGeoManager.h"
55 #include "TGeoNode.h"
56 #include <list>
57 
58 #include "L1Event.h"
59 #include "TMatrixD.h"
60 #include "TROOT.h"
61 #include "TRandom3.h"
62 #include "TVector3.h"
63 #include "TVectorD.h"
64 #include <TFile.h>
65 
66 #include "CbmMvdDetector.h"
67 #include "CbmMvdStationPar.h"
68 
69 #include <fstream>
70 #include <iomanip>
71 #include <iostream>
72 #include <vector>
73 
74 using std::cout;
75 using std::endl;
76 using std::ios;
77 using std::vector;
78 
79 #include "KFTopoPerformance.h"
80 
81 ClassImp(CbmL1)
82 
83  static L1Algo algo_static _fvecalignment;
84 
85 //L1AlgoInputData* fData_static _fvecalignment;
86 
87 CbmL1* CbmL1::fInstance = 0;
88 
89 
90 CbmL1::CbmL1()
91  : FairTask("L1")
92  , algo(0)
93  , // for access to L1 Algorithm from L1::Instance
94  fDigiFile()
95  , fUseHitErrors(0)
96  , fmCBMmode(0)
97  , fGlobalMode(0)
98  , vRTracks()
99  , // reconstructed tracks
100  vFileEvent()
101  , vHitStore()
102  , vMCPoints()
103  , nMvdPoints(0)
104  , vMCPoints_in_Time_Slice()
105  , NStation(0)
106  , NMvdStations(0)
107  , NStsStations(0)
108  , NMuchStations(0)
109  , NTrdStations(0)
110  , NTOFStation(0)
111  , // number of detector stations (all\sts\mvd)
112  fPerformance(0)
113  , fSTAPDataMode(0)
114  , fSTAPDataDir("")
115  ,
116 
117  fTrackingLevel(2)
118  , // really doesn't used
119  fMomentumCutOff(0.1)
120  , // really doesn't used
121  fGhostSuppression(1)
122  , // really doesn't used
123  fUseMVD(0)
124  , // really doesn't used
125  fUseMUCH(0)
126  , fUseTRD(0)
127  , fUseTOF(0)
128  ,
129 
130  PrimVtx()
131  , fTimeSlice(nullptr)
132  , fEventList(nullptr)
133  , listStsDigi()
134  , fStsPoints(0)
135  , fMCTracks(0)
136  , fMvdPoints(0)
137  ,
138  //listMCTracks (0),
139  //listStsDigi(0),
140  listStsPts(0)
141  , listStsDigiMatch(0)
142  , listStsClusters(0)
143  , listStsHits(0)
144  , listStsHitMatch(0)
145  , listStsClusterMatch(0)
146  ,
147 
148  listMvdPts(0)
149  , listMvdHits(0)
150  , listMvdDigiMatches(0)
151  , listMvdHitMatches(0)
152  ,
153 
154  nMuchPoints(0)
155  , fMuchPoints(nullptr)
156  , listMuchHitMatches(nullptr)
157  , fDigiMatchesMuch(nullptr)
158  , fClustersMuch(nullptr)
159  , fMuchPixelHits(nullptr)
160  , fDigisMuch(nullptr)
161  , fTrdDigiPar(nullptr)
162  , fTrdModuleInfo(nullptr)
163  , fTrdPoints(nullptr)
164  , listTrdHits(nullptr)
165  , fTrdHitMatches(nullptr)
166  , fTofPoints(nullptr)
167  , fTofHitDigiMatches(nullptr)
168  , fTofHits(nullptr)
169  , fDigiPar(nullptr)
170  ,
171 
172  fPerfFile(nullptr)
173  , fHistoDir(nullptr)
174  , vStsHits()
175  , vMCTracks()
176  , vHitMCRef()
177  , dFEI2vMCPoints()
178  , dFEI2vMCTracks()
179  , fData(nullptr)
180  , histodir(nullptr)
181  , fFindParticlesMode()
182  , fStsMatBudgetFileName("")
183  , fMvdMatBudgetFileName("")
184  , fMuchMatBudgetFileName("")
185  , fTrdMatBudgetFileName("")
186  , fTofMatBudgetFileName("")
187  , fExtrapolateToTheEndOfSTS(false)
188  , fTimesliceMode(0)
189  , fTopoPerformance(nullptr)
190  , fEventEfficiency() {
191  if (!fInstance) fInstance = this;
192 }
193 
194 CbmL1::CbmL1(const char* name,
195  Int_t iVerbose,
196  Int_t _fPerformance,
197  int fSTAPDataMode_,
198  TString fSTAPDataDir_,
199  int findParticleMode_)
200  : FairTask(name, iVerbose)
201  , algo(0)
202  , // for access to L1 Algorithm from L1::Instance
203  fDigiFile()
204  , fUseHitErrors(0)
205  , fmCBMmode(0)
206  , fGlobalMode(0)
207  , vRTracks()
208  , // reconstructed tracks
209  vFileEvent()
210  , vHitStore()
211  , vMCPoints()
212  , nMvdPoints(0)
213  , vMCPoints_in_Time_Slice()
214  , NStation(0)
215  , NMvdStations(0)
216  , NStsStations(0)
217  , NMuchStations(0)
218  , NTrdStations(0)
219  , NTOFStation(0)
220  , // number of detector stations (all\sts\mvd)
221  fPerformance(_fPerformance)
222  , fSTAPDataMode(fSTAPDataMode_)
223  , fSTAPDataDir(fSTAPDataDir_)
224  ,
225 
226  fTrackingLevel(2)
227  , // really doesn't used
228  fMomentumCutOff(0.1)
229  , // really doesn't used
230  fGhostSuppression(1)
231  , // really doesn't used
232  fUseMVD(0)
233  , // really doesn't used
234  fUseMUCH(0)
235  , fUseTRD(0)
236  , fUseTOF(0)
237  ,
238 
239 
240  PrimVtx()
241  , fTimeSlice(nullptr)
242  , fEventList(nullptr)
243  , listStsDigi(0)
244  , fStsPoints(0)
245  , fMCTracks(0)
246  , fMvdPoints(NULL)
247  ,
248  //listMCTracks (0),
249 
250  listStsPts(0)
251  , listStsDigiMatch(0)
252  , listStsClusters(0)
253  , listStsHits(0)
254  , listStsHitMatch(0)
255  , listStsClusterMatch(0)
256  ,
257 
258  listMvdPts(0)
259  , listMvdHits(0)
260  , listMvdDigiMatches(0)
261  , listMvdHitMatches(0)
262  ,
263 
264  nMuchPoints(0)
265  , fMuchPoints(nullptr)
266  , listMuchHitMatches(nullptr)
267  , fDigiMatchesMuch(nullptr)
268  , fClustersMuch(nullptr)
269  , fMuchPixelHits(nullptr)
270  , fDigisMuch(nullptr)
271  , fTrdDigiPar(nullptr)
272  , fTrdModuleInfo(nullptr)
273  , fTrdPoints(nullptr)
274  , listTrdHits(nullptr)
275  , fTrdHitMatches(nullptr)
276  , fTofPoints(nullptr)
277  , fTofHitDigiMatches(nullptr)
278  , fTofHits(nullptr)
279  , fDigiPar(nullptr)
280  ,
281 
282  fPerfFile(nullptr)
283  , fHistoDir(nullptr)
284  , vStsHits()
285  , vMCTracks()
286  , vHitMCRef()
287  , dFEI2vMCPoints()
288  , dFEI2vMCTracks()
289  , fData(nullptr)
290  , histodir(nullptr)
291  , fFindParticlesMode(findParticleMode_)
292  , fStsMatBudgetFileName("")
293  , fMvdMatBudgetFileName("")
294  , fMuchMatBudgetFileName("")
295  , fTrdMatBudgetFileName("")
296  , fTofMatBudgetFileName("")
297  , fExtrapolateToTheEndOfSTS(false)
298  , fTimesliceMode(0)
299  , fTopoPerformance(nullptr)
300  , fEventEfficiency() {
301  if (!fInstance) fInstance = this;
302 }
303 
304 CbmL1::~CbmL1() {
305  if (fInstance == this) fInstance = 0;
306 }
307 
308 
309 void CbmL1::SetParContainers() {
310  FairRunAna* ana = FairRunAna::Instance();
311  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
312  fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
313  // fTrdDigiPar = (CbmTrdDigiPar*)(FairRunAna::Instance()->GetRuntimeDb()->getContainer("CbmTrdDigiPar"));
314  fTrdDigiPar = (CbmTrdParSetDigi*) (rtdb->getContainer("CbmTrdParSetDigi"));
315 
316  // Trigger loading of STS parameters
317  // If L1 runs standalone the parameters are not required by any STS task
318  // but they are needed by CbmStsSetup
319  fStsParSetModule =
320  dynamic_cast<CbmStsParSetModule*>(rtdb->getContainer("CbmStsParSetModule"));
321  fStsParSetSensor =
322  dynamic_cast<CbmStsParSetSensor*>(rtdb->getContainer("CbmStsParSetSensor"));
323  fStsParSetSensorCond = dynamic_cast<CbmStsParSetSensorCond*>(
324  rtdb->getContainer("CbmStsParSetSensorCond"));
325 }
326 
327 
328 InitStatus CbmL1::ReInit() {
329  SetParContainers();
330  return Init();
331 }
332 
333 InitStatus CbmL1::Init() {
334  fData = new L1AlgoInputData();
335 
336  if (fVerbose > 1) {
337  char y[20] = " [0;33;44m"; // yellow
338  char Y[20] = " [1;33;44m"; // yellow bold
339  char W[20] = " [1;37;44m"; // white bold
340  char o[20] = " [0m\n"; // color off
341  Y[0] = y[0] = W[0] = o[0] = 0x1B; // escape character
342 
343  cout << endl << endl;
344  cout << " " << W
345  << " "
346  << o;
347  cout << " " << W
348  << " ===////====================================================== "
349  << o;
350  cout << " " << W
351  << " = = "
352  << o;
353  cout << " " << W << " = " << Y << "L1 on-line finder"
354  << W << " = " << o;
355  cout << " " << W
356  << " = = "
357  << o;
358  cout << " " << W << " = " << W << "Cellular Automaton 3.1 Vector" << y
359  << " with " << W << "KF Quadro" << y << " technology" << W << " = "
360  << o;
361  cout << " " << W
362  << " = = "
363  << o;
364  cout << " " << W << " = " << y << "Designed for CBM collaboration" << W
365  << " = " << o;
366  cout << " " << W << " = " << y << "All rights reserved" << W
367  << " = " << o;
368  cout << " " << W
369  << " = = "
370  << o;
371  cout << " " << W
372  << " ========================================================////= "
373  << o;
374  cout << " " << W
375  << " "
376  << o;
377  cout << endl << endl;
378  }
379 
380  FairRootManager* fManger = FairRootManager::Instance();
381 
382  FairRunAna* Run = FairRunAna::Instance();
383  {
384  fUseMVD = 1;
385  CbmStsFindTracks* FindTask =
386  L1_DYNAMIC_CAST<CbmStsFindTracks*>(Run->GetTask("STSFindTracks"));
387  if (FindTask) fUseMVD = FindTask->MvdUsage();
388  }
389 
390  fUseMVD = 0;
391 
392  histodir = gROOT->mkdir("L1");
393 
394 
395  // turn on reconstruction in sub-detectors
396 
397  fUseMUCH = 0;
398  fUseTRD = 0;
399  fUseTOF = 0;
400 
401  if (fmCBMmode) {
402  fUseMUCH = 1;
403  fUseTRD = 1;
404  fUseTOF = 1;
405  }
406 
407 
408  if (fGlobalMode) {
409  fUseMUCH = 1;
410  fUseTRD = 1;
411  fUseTOF = 1;
412  }
413 
414 
415  fStsPoints = 0;
416  fMvdPoints = 0;
417  fMuchPoints = 0;
418  fTrdPoints = 0;
419  fTofPoints = 0;
420  fMCTracks = 0;
421 
422 
423  listStsClusters = 0;
424  listStsDigi.clear();
425  vFileEvent.clear();
426 
427 
428  if (fTimesliceMode) { // Time-slice mode selected
429  LOG(info) << GetName() << ": running in time-slice mode.";
430  fTimeSlice = NULL;
431  fTimeSlice = (CbmTimeSlice*) fManger->GetObject("TimeSlice.");
432  if (fTimeSlice == NULL)
433  LOG(fatal) << GetName() << ": No time slice branch in tree!";
434 
435  } //? time-slice mode
436 
437  else // event mode
438  LOG(info) << GetName() << ": running in event mode.";
439 
440 
441  listStsClusters =
442  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsCluster"));
443  listStsHitMatch =
444  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHitMatch"));
445  listStsClusterMatch =
446  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsClusterMatch"));
447 
448  listStsHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHit"));
449 
450  if (!fUseMUCH) {
451  fMuchPixelHits = 0;
452 
453  fDigisMuch = 0;
454  fDigiMatchesMuch = 0;
455  fClustersMuch = 0;
456 
457  fMuchPoints = 0;
458  listMuchHitMatches = 0;
459 
460  } else {
461 
462 
463  fMuchPixelHits = (TClonesArray*) fManger->GetObject("MuchPixelHit");
464  }
465 
466  if (!fUseTRD) {
467  fTrdPoints = 0;
468  fTrdHitMatches = 0;
469  fTrdPoints = 0;
470  listTrdHits = 0;
471 
472  } else {
473 
474  listTrdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("TrdHit"));
475  }
476 
477  if (!fUseTOF) {
478  fTofPoints = 0;
479  fTofHitDigiMatches = 0;
480  fTofHits = 0;
481 
482  } else {
483 
484 
485  fTofHits = (TClonesArray*) fManger->GetObject("TofHit");
486  }
487 
488  if (fPerformance) {
489  CbmMCDataManager* mcManager =
490  (CbmMCDataManager*) fManger->GetObject("MCDataManager");
491  if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!";
492 
493  fStsPoints = mcManager->InitBranch("StsPoint");
494  if (!fTimesliceMode) fMvdPoints = mcManager->InitBranch("MvdPoint");
495  fMCTracks = mcManager->InitBranch("MCTrack");
496  if (NULL == fStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
497  if (NULL == fMCTracks) LOG(fatal) << GetName() << ": No MCTrack data!";
498 
499  listStsPts = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsPoint"));
500 
501  if (fTimesliceMode) {
502  fEventList = (CbmMCEventList*) fManger->GetObject("MCEventList.");
503  if (NULL == fEventList)
504  LOG(fatal) << GetName() << ": No MCEventList data!";
505  }
506 
507  if (!fUseMVD) {
508  listMvdPts = 0;
509  listMvdHitMatches = 0;
510  } else {
511  listMvdPts =
512  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdPoint"));
513  listMvdDigiMatches =
514  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdDigiMatch"));
515  listMvdHitMatches =
516  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHitMatch"));
517 
518  if (!listMvdHitMatches && listMvdPts)
519  LOG(error)
520  << "No listMvdHitMatches provided, performance is not done correctly";
521  }
522 
523  if (!fUseTRD) {
524  fTrdPoints = 0;
525  fTrdHitMatches = 0;
526  fTrdPoints = 0;
527 
528  } else {
529  fTrdHitMatches = (TClonesArray*) fManger->GetObject("TrdHitMatch");
530  fTrdPoints = mcManager->InitBranch("TrdPoint");
531  }
532 
533  if (!fUseMUCH) {
534  fMuchPoints = 0;
535  listMuchHitMatches = 0;
536 
537  } else {
538 
539  fDigisMuch = (TClonesArray*) fManger->GetObject("MuchDigi");
540  fDigiMatchesMuch = (TClonesArray*) fManger->GetObject("MuchDigiMatch");
541  fClustersMuch = (TClonesArray*) fManger->GetObject("MuchCluster");
542  fMuchPoints = mcManager->InitBranch("MuchPoint");
543  listMuchHitMatches =
544  L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MuchPixelHitMatch"));
545  }
546 
547  if (!fUseTOF) {
548  fTofPoints = 0;
549  fTofHitDigiMatches = 0;
550  } else {
551 
552  fTofPoints = mcManager->InitBranch("TofPoint");
553  fTofHitDigiMatches =
554  static_cast<TClonesArray*>(fManger->GetObject("TofCalDigiMatch"));
555  }
556 
557  } else {
558  listMvdPts = 0;
559  listMvdHitMatches = 0;
560  fTrdPoints = 0;
561  fTrdHitMatches = 0;
562  fTrdPoints = 0;
563  fMuchPoints = 0;
564  listMuchHitMatches = 0;
565  fTofPoints = 0;
566  fTofHitDigiMatches = 0;
567  }
568  if (!fUseMVD) {
569  listMvdHits = 0;
570  } else {
571  listMvdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHit"));
572  }
573 
574  NMvdStations = 0;
575  NStsStations = 0;
576  NMuchStations = 0;
577  NTrdStations = 0;
578  NTOFStation = 0;
579  NStation = 0;
580 
581  if (fUseTOF) NTOFStation = 1;
582 
583  algo = &algo_static;
584 
585  vector<fscal> geo;
586  geo.clear();
587 
588  for (int i = 0; i < 3; i++) {
589  Double_t point[3] = {0, 0, 2.5 * i};
590  Double_t B[3] = {0, 0, 0};
591  if (CbmKF::Instance()->GetMagneticField())
592  CbmKF::Instance()->GetMagneticField()->GetFieldValue(point, B);
593  geo.push_back(2.5 * i);
594  geo.push_back(B[0]);
595  geo.push_back(B[1]);
596  geo.push_back(B[2]);
597  }
598 
599 
600  CbmMuchGeoScheme* fGeoScheme = CbmMuchGeoScheme::Instance();
601 
602  if (fUseMUCH) {
603  TFile* file = new TFile(fDigiFile);
604  TObjArray* stations = (TObjArray*) file->Get("stations");
605  fGeoScheme->Init(stations, 0);
606  for (int iStation = 0; iStation < fGeoScheme->GetNStations(); iStation++) {
607  const CbmMuchStation* station = fGeoScheme->GetStation(iStation);
608  int nLayers = station->GetNLayers();
609  NMuchStations += nLayers;
610  }
611  }
612 
613  // count TRD stations
614  if (fUseTRD) {
615  Int_t layerCounter = 0;
616  TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
617  for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast();
618  iTopNode++) {
619  TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
620  if (TString(topNode->GetName()).Contains("trd")) {
621  TObjArray* layers = topNode->GetNodes();
622  for (Int_t iLayer = 0; iLayer < layers->GetEntriesFast(); iLayer++) {
623  TGeoNode* layer = static_cast<TGeoNode*>(layers->At(iLayer));
624  if (TString(layer->GetName()).Contains("layer")) { layerCounter++; }
625  }
626  }
627  }
628  NTrdStations = layerCounter;
629  }
630 
631  // count ToF parameters
632 
633  float z_average = 0;
634 
635  float z_min = 100000;
636  float z_max = 0;
637 
638  float r_max = 0;
639  float r_min = 100000;
640 
641  if (fUseTOF) {
642 
643  Int_t nrOfCells = fDigiPar->GetNrOfModules();
644 
645  for (Int_t icell = 0; icell < nrOfCells; ++icell) {
646 
647  Int_t cellId = fDigiPar->GetCellId(icell);
648  CbmTofCell* fCellInfo = fDigiPar->GetCell(cellId);
649 
650  Double_t x = fCellInfo->GetX();
651  Double_t y = fCellInfo->GetY();
652  Double_t z = fCellInfo->GetZ();
653  // Double_t dx = fCellInfo->GetSizex();
654  // Double_t dy = fCellInfo->GetSizey();
655 
656  if (z < z_min) z_min = z;
657  if (z > z_max) z_max = z;
658 
659  if (x < r_min) r_min = x;
660  if (x > r_max) r_max = x;
661 
662  if (y < r_min) r_min = y;
663  if (y > r_max) r_max = y;
664 
665 
666  z_average += z;
667  }
668  z_average = z_average / nrOfCells;
669  }
670 
671  CbmStsSetup* stsSetup = CbmStsSetup::Instance();
672  if (!stsSetup->IsInit()) { stsSetup->Init(nullptr); }
673  if (!stsSetup->IsModuleParsInit()) {
674  stsSetup->SetModuleParameters(fStsParSetModule);
675  }
676  if (!stsSetup->IsSensorParsInit()) {
677  stsSetup->SetSensorParameters(fStsParSetSensor);
678  }
679  if (!stsSetup->IsSensorCondInit()) {
680  stsSetup->SetSensorConditions(fStsParSetSensorCond);
681  }
682 
683  NMvdStations = (fUseMVD) ? CbmKF::Instance()->vMvdMaterial.size() : 0;
684  NStsStations = CbmStsSetup::Instance()->GetNofStations();
685  NStation =
686  NMvdStations + NStsStations + NMuchStations + NTrdStations + NTOFStation;
687  geo.push_back(NStation);
688  geo.push_back(NMvdStations);
689  geo.push_back(NStsStations);
690 
691  // field
692  const int M = 5; // polinom order
693  const int N = 21; //(M+1)*(M+2)/2;
694 
695  // { // field at the z=0 plane
696  // const double Xmax = 10, Ymax = 10;
697  // const double z = 0.;
698  // double dx = 1.; // step for the field approximation
699  // double dy = 1.;
700  // if( dx > Xmax/N/2 ) dx = Xmax/N/4.;
701  // if( dy > Ymax/N/2 ) dy = Ymax/N/4.;
702  //
703  // TMatrixD A(N,N);
704  // TVectorD b0(N), b1(N), b2(N);
705  // for( int i=0; i<N; i++){
706  // for( int j=0; j<N; j++) A(i,j)==0.;
707  // b0(i)=b1(i)=b2(i) = 0.;
708  // }
709  // for( double x=-Xmax; x<=Xmax; x+=dx )
710  // for( double y=-Ymax; y<=Ymax; y+=dy )
711  // {
712  // double r = sqrt(fabs(x*x/Xmax/Xmax+y/Ymax*y/Ymax));
713  // if( r>1. ) continue;
714  // Double_t w = 1./(r*r+1);
715  // Double_t p[3] = { x, y, z};
716  // Double_t B[3] = {0.,0.,0.};
717  // CbmKF::Instance()->GetMagneticField()->GetFieldValue(p, B);
718  // TVectorD m(N);
719  // m(0)=1;
720  // for( int i=1; i<=M; i++){
721  // int k = (i-1)*(i)/2;
722  // int l = i*(i+1)/2;
723  // for( int j=0; j<i; j++ ) m(l+j) = x*m(k+j);
724  // m(l+i) = y*m(k+i-1);
725  // }
726  //
727  // TVectorD mt = m;
728  // for( int i=0; i<N; i++){
729  // for( int j=0; j<N;j++) A(i,j)+=w*m(i)*m(j);
730  // b0(i)+=w*B[0]*m(i);
731  // b1(i)+=w*B[1]*m(i);
732  // b2(i)+=w*B[2]*m(i);
733  // }
734  // }
735  // double det;
736  // A.Invert(&det);
737  // TVectorD c0 = A*b0, c1 = A*b1, c2 = A*b2;
738  //
739  // targetFieldSlice = new L1FieldSlice;
740  // for(int i=0; i<N; i++){
741  // targetFieldSlice->cx[i] = c0(i);
742  // targetFieldSlice->cy[i] = c1(i);
743  // targetFieldSlice->cz[i] = c2(i);
744  // }
745  //
746  // } // target field
747 
748  for (Int_t ist = 0; ist < NStation; ist++) {
749  double C[3][N];
750  double z = 0;
751  double Xmax = 0, Ymax = 0;
752  if (ist < NMvdStations) {
753 
754 
755  CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
756  CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
757 
758 
759  CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
760  geo.push_back(1);
761  geo.push_back(t.z);
762  geo.push_back(t.dz);
763  geo.push_back(t.r);
764  geo.push_back(t.R);
765  geo.push_back(t.RadLength);
766 
767  fscal f_phi = 0, f_sigma = mvdStationPar->GetXRes(ist) / 10000,
768  b_phi = 3.14159265358 / 2.,
769  b_sigma = mvdStationPar->GetYRes(ist) / 10000;
770  geo.push_back(f_phi);
771  geo.push_back(f_sigma);
772  geo.push_back(b_phi);
773  geo.push_back(b_sigma);
774  z = t.z;
775  Xmax = Ymax = t.R;
776  }
777 
778 
779  if (ist >= NMvdStations && ist < (NMvdStations + NStsStations)) {
780  CbmStsStation* station =
781  CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
782  geo.push_back(0);
783  geo.push_back(station->GetZ());
784  geo.push_back(station->GetSensorD());
785  geo.push_back(0);
786  geo.push_back(station->GetYmax() < station->GetXmax()
787  ? station->GetXmax()
788  : station->GetYmax());
789  geo.push_back(station->GetRadLength());
790 
791  double Pi = 3.14159265358;
792 
793  fscal f_phi, f_sigma, b_phi, b_sigma; // angle and sigma front/back side
794  f_phi = station->GetSensorRotation();
795  b_phi = f_phi;
796  f_phi += station->GetSensorStereoAngle(0) * Pi / 180.;
797  b_phi += station->GetSensorStereoAngle(1) * Pi / 180.;
798  f_sigma = station->GetSensorPitch(0) / TMath::Sqrt(12);
799  b_sigma = f_sigma;
800  //f_sigma *= cos(f_phi); // TODO: think about this
801  //b_sigma *= cos(b_phi);
802 
803  geo.push_back(f_phi);
804  geo.push_back(f_sigma);
805  geo.push_back(b_phi);
806  geo.push_back(b_sigma);
807  z = station->GetZ();
808 
809  Xmax = station->GetXmax();
810  Ymax = station->GetYmax();
811  }
812 
813  if ((ist < (NMvdStations + NStsStations + NMuchStations))
814  && (ist >= (NMvdStations + NStsStations))) {
815 
816  int iStation = (ist - NMvdStations - NStsStations) / 3;
817 
818 
819  CbmMuchStation* station =
820  (CbmMuchStation*) fGeoScheme->GetStation(iStation);
821 
822  CbmMuchLayer* layer =
823  station->GetLayer((ist - NMvdStations - NStsStations) % 3);
824 
825  // CbmMuchModuleGem* module = (CbmMuchModuleGem*) CbmMuchGeoScheme::Instance()->GetModule(0,0,0,0);
826 
827  // vector<CbmMuchPad*> pads = module->GetPads();
828 
829  z = layer->GetZ();
830 
831  geo.push_back(2);
832  geo.push_back(z);
833  geo.push_back(layer->GetDz());
834  geo.push_back(0);
835  geo.push_back(100); //station->GetRmax()
836  geo.push_back(0);
837 
838  fscal f_phi = 0, f_sigma = 0.1, b_phi = 3.14159265358 / 2., b_sigma = 0.1;
839  geo.push_back(f_phi);
840  geo.push_back(f_sigma);
841  geo.push_back(b_phi);
842  geo.push_back(b_sigma);
843 
844 
845  Xmax = 100; //station->GetRmax();
846  Ymax = 100; //station->GetRmax();
847  }
848 
849  // int num = 0;
850 
851  if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations))
852  && (ist >= (NMvdStations + NStsStations + NMuchStations))) {
853 
854  int num = ist - NMvdStations - NStsStations - NMuchStations;
855 
856  // if (num == 0) continue;//true_station = 0;
857  // if (!true_station) continue;
858 
859  // Int_t nrModules = fTrdDigiPar->GetNrOfModules();
860 
861  int ModuleId = fTrdDigiPar->GetModuleId(num);
862 
863  CbmTrdParModDigi* module =
864  (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(ModuleId);
865 
866  // if (!true_station[ist]) continue;
867 
868  if (num == 0 || num == 2 || num == 4) geo.push_back(3);
869  if (num == 1 || num == 3) geo.push_back(6);
870  geo.push_back(module->GetZ());
871 
872  geo.push_back(2 * module->GetSizeZ());
873  geo.push_back(0);
874  geo.push_back(2 * module->GetSizeX());
875  geo.push_back(10);
876 
877  fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2.,
878  b_sigma = 1 / 10000;
879  geo.push_back(f_phi);
880  geo.push_back(f_sigma);
881  geo.push_back(b_phi);
882  geo.push_back(b_sigma);
883  Xmax = Ymax = 20;
884  }
885 
886  if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations
887  + NTOFStation))
888  && (ist
889  >= (NMvdStations + NStsStations + NMuchStations + NTrdStations))) {
890 
891  geo.push_back(4);
892  geo.push_back(z_average);
893  geo.push_back(z_max - z_min);
894  geo.push_back(0);
895  geo.push_back(r_max);
896  geo.push_back(10);
897 
898  fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2.,
899  b_sigma = 1 / 10000;
900  geo.push_back(f_phi);
901  geo.push_back(f_sigma);
902  geo.push_back(b_phi);
903  geo.push_back(b_sigma);
904  Xmax = Ymax = 20;
905  }
906 
907  double dx = 1.; // step for the field approximation
908  double dy = 1.;
909 
910  if (dx > Xmax / N / 2) dx = Xmax / N / 4.;
911  if (dy > Ymax / N / 2) dy = Ymax / N / 4.;
912 
913  for (int i = 0; i < 3; i++)
914  for (int k = 0; k < N; k++)
915  C[i][k] = 0;
916  TMatrixD A(N, N);
917  TVectorD b0(N), b1(N), b2(N);
918  for (int i = 0; i < N; i++) {
919  for (int j = 0; j < N; j++)
920  A(i, j) = 0.;
921  b0(i) = b1(i) = b2(i) = 0.;
922  }
923 
924  if (CbmKF::Instance()->GetMagneticField()) {
925  for (double x = -Xmax; x <= Xmax; x += dx)
926  for (double y = -Ymax; y <= Ymax; y += dy) {
927  double r = sqrt(fabs(x * x / Xmax / Xmax + y / Ymax * y / Ymax));
928  if (r > 1.) continue;
929  Double_t w = 1. / (r * r + 1);
930  Double_t p[3] = {x, y, z};
931  Double_t B[3] = {0., 0., 0.};
932  CbmKF::Instance()->GetMagneticField()->GetFieldValue(p, B);
933  TVectorD m(N);
934  m(0) = 1;
935  for (int i = 1; i <= M; i++) {
936  int k = (i - 1) * (i) / 2;
937  int l = i * (i + 1) / 2;
938  for (int j = 0; j < i; j++)
939  m(l + j) = x * m(k + j);
940  m(l + i) = y * m(k + i - 1);
941  }
942 
943  TVectorD mt = m;
944  for (int i = 0; i < N; i++) {
945  for (int j = 0; j < N; j++)
946  A(i, j) += w * m(i) * m(j);
947  b0(i) += w * B[0] * m(i);
948  b1(i) += w * B[1] * m(i);
949  b2(i) += w * B[2] * m(i);
950  }
951  }
952  double det;
953  A.Invert(&det);
954  TVectorD c0 = A * b0, c1 = A * b1, c2 = A * b2;
955  for (int i = 0; i < N; i++) {
956  C[0][i] = c0(i);
957  C[1][i] = c1(i);
958  C[2][i] = c2(i);
959  }
960  }
961  geo.push_back(N);
962  for (int k = 0; k < 3; k++) {
963  for (int j = 0; j < N; j++) {
964  geo.push_back(C[k][j]);
965  }
966  }
967  }
968 
969  geo.push_back(fTrackingLevel);
970  geo.push_back(fMomentumCutOff);
971  geo.push_back(fGhostSuppression);
972 
973  {
974  if (fSTAPDataMode % 2 == 1) { // 1,3
975  WriteSTAPGeoData(geo);
976  };
977  //if(fSTAPDataMode >= 2){ // 2,3
978  // int ind2, ind = geo.size();
979  // ReadSTAPGeoData(geo, ind2);
980  // if (ind2 != ind) LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful.";
981  //};
982  }
983 
984  if (fSTAPDataMode >= 2) { // 2,3
985  int ind2, ind = geo.size();
986  ReadSTAPGeoData(geo, ind2);
987  if (ind2 != ind)
988  LOG(error) << "-E- CbmL1: Read geometry from file "
989  << fSTAPDataDir + "geo_algo.txt was NOT successful.";
990  }
991 
992  algo->Init(geo, fUseHitErrors, fmCBMmode);
993  geo.clear();
994 
995 
996  algo->fRadThick.resize(algo->NStations);
997 
998  // Read STS MVD TRD MuCh ToF Radiation Thickness table
999  TString stationName = "Radiation Thickness [%], Station";
1000  if (fUseMVD) {
1001  if (fMvdMatBudgetFileName != "") {
1002  TFile* oldfile = gFile;
1003  TFile* rlFile = new TFile(fMvdMatBudgetFileName);
1004  cout << "MVD Material budget file is " << fMvdMatBudgetFileName << ".\n";
1005  for (int j = 0, iSta = 0; iSta < algo->NMvdStations; iSta++, j++) {
1006  TString stationNameMvd = stationName;
1007  stationNameMvd += j;
1008  TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameMvd);
1009  if (!hStaRadLen) {
1010  cout << "L1: incorrect " << fMvdMatBudgetFileName << " file. No "
1011  << stationNameMvd << "\n";
1012  exit(1);
1013  }
1014  const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
1015  const float RMax =
1016  hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
1017  algo->fRadThick[iSta].SetBins(NBins, RMax);
1018  algo->fRadThick[iSta].table.resize(NBins);
1019 
1020  for (int iB = 0; iB < NBins; iB++) {
1021  algo->fRadThick[iSta].table[iB].resize(NBins);
1022  for (int iB2 = 0; iB2 < NBins; iB2++) {
1023  algo->fRadThick[iSta].table[iB][iB2] =
1024  0.01 * hStaRadLen->GetBinContent(iB, iB2);
1025  // Correction for holes in material map
1026  if (algo->fRadThick[iSta].table[iB][iB2]
1027  < algo->vStations[iSta].materialInfo.RadThick[0])
1028  if (iB2 > 0 && iB2 < NBins - 1)
1029  algo->fRadThick[iSta].table[iB][iB2] =
1030  TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
1031  0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
1032  // Correction for the incorrect harcoded value of RadThick of MVD stations
1033  if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015)
1034  algo->fRadThick[iSta].table[iB][iB2] = 0.0015;
1035  // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
1036  }
1037  }
1038  }
1039  rlFile->Close();
1040  rlFile->Delete();
1041  gFile = oldfile;
1042  } else {
1043  LOG(warn) << "No MVD material budget file is found. Homogenious budget "
1044  "will be used";
1045  for (int iSta = 0; iSta < algo->NMvdStations; iSta++) {
1046  cout << iSta << endl;
1047  algo->fRadThick[iSta].SetBins(1,
1048  100); // mvd should be in +-100 cm square
1049  algo->fRadThick[iSta].table.resize(1);
1050  algo->fRadThick[iSta].table[0].resize(1);
1051  algo->fRadThick[iSta].table[0][0] =
1052  algo->vStations[iSta].materialInfo.RadThick[0];
1053  }
1054  }
1055  }
1056  if (fStsMatBudgetFileName != "") {
1057  TFile* oldfile = gFile;
1058  TFile* rlFile = new TFile(fStsMatBudgetFileName);
1059  cout << "STS Material budget file is " << fStsMatBudgetFileName << ".\n";
1060  for (int j = 1, iSta = algo->NMvdStations;
1061  iSta < (algo->NMvdStations + NStsStations);
1062  iSta++, j++) {
1063  TString stationNameSts = stationName;
1064  stationNameSts += j;
1065  TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
1066  if (!hStaRadLen) {
1067  cout << "L1: incorrect " << fStsMatBudgetFileName << " file. No "
1068  << stationNameSts << "\n";
1069  exit(1);
1070  }
1071  const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
1072  const float RMax =
1073  hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
1074  algo->fRadThick[iSta].SetBins(NBins, RMax);
1075  algo->fRadThick[iSta].table.resize(NBins);
1076 
1077  for (int iB = 0; iB < NBins; iB++) {
1078  algo->fRadThick[iSta].table[iB].resize(NBins);
1079  for (int iB2 = 0; iB2 < NBins; iB2++) {
1080  algo->fRadThick[iSta].table[iB][iB2] =
1081  0.01 * hStaRadLen->GetBinContent(iB, iB2);
1082  if (algo->fRadThick[iSta].table[iB][iB2]
1083  < algo->vStations[iSta].materialInfo.RadThick[0])
1084  algo->fRadThick[iSta].table[iB][iB2] =
1085  algo->vStations[iSta].materialInfo.RadThick[0];
1086  // cout << " iSta " << iSta << " iB " << iB << "iB2 " << iB2 << " RadThick " << algo->fRadThick[iSta].table[iB][iB2] << endl;
1087  }
1088  }
1089  }
1090  rlFile->Close();
1091  rlFile->Delete();
1092  gFile = oldfile;
1093  } else {
1094  LOG(warn) << "No STS material budget file is found. Homogenious budget "
1095  "will be used";
1096  for (int iSta = algo->NMvdStations;
1097  iSta < (algo->NMvdStations + NStsStations);
1098  iSta++) {
1099  algo->fRadThick[iSta].SetBins(1, 100);
1100  algo->fRadThick[iSta].table.resize(1);
1101  algo->fRadThick[iSta].table[0].resize(1);
1102  algo->fRadThick[iSta].table[0][0] =
1103  algo->vStations[iSta].materialInfo.RadThick[0];
1104  }
1105  }
1106 
1107  if (fUseMUCH)
1108  if (fMuchMatBudgetFileName != "") {
1109  TFile* oldfile = gFile;
1110  TFile* rlFile = new TFile(fMuchMatBudgetFileName);
1111  cout << "Much Material budget file is " << fMuchMatBudgetFileName
1112  << ".\n";
1113  for (int j = 1, iSta = (NStsStations + NMvdStations);
1114  iSta < (NStsStations + NMvdStations + NMuchStations);
1115  iSta++, j++) {
1116  TString stationNameSts = stationName;
1117  stationNameSts += j;
1118  TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
1119  if (!hStaRadLen) {
1120  cout << "L1: incorrect " << fMuchMatBudgetFileName << " file. No "
1121  << stationNameSts << "\n";
1122  exit(1);
1123  }
1124 
1125 
1126  const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
1127  const float RMax =
1128  hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
1129  algo->fRadThick[iSta].SetBins(NBins, RMax);
1130  algo->fRadThick[iSta].table.resize(NBins);
1131 
1132  for (int iB = 0; iB < NBins; iB++) {
1133  algo->fRadThick[iSta].table[iB].resize(NBins);
1134  float hole = 0.15;
1135  for (int iB2 = 0; iB2 < NBins; iB2++) {
1136  algo->fRadThick[iSta].table[iB][iB2] =
1137  0.01 * hStaRadLen->GetBinContent(iB, iB2);
1138  // Correction for holes in material map
1139  //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
1140 
1141  if (iB2 > 0 && iB2 < NBins - 1)
1142  algo->fRadThick[iSta].table[iB][iB2] =
1143  TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
1144  0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
1145  // Correction for the incorrect harcoded value of RadThick of MVD stations
1146 
1147  if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015)
1148  hole = algo->fRadThick[iSta].table[iB][iB2];
1149  if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015)
1150  algo->fRadThick[iSta].table[iB][iB2] = hole;
1151  // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
1152  }
1153  }
1154  }
1155  rlFile->Close();
1156  rlFile->Delete();
1157  gFile = oldfile;
1158  } else {
1159  LOG(warn) << "No Much material budget file is found. Homogenious budget "
1160  "will be used";
1161  for (int iSta = (NStsStations + NMvdStations);
1162  iSta < (NStsStations + NMvdStations + NMuchStations);
1163  iSta++) {
1164  algo->fRadThick[iSta].SetBins(1, 100);
1165  algo->fRadThick[iSta].table.resize(1);
1166  algo->fRadThick[iSta].table[0].resize(1);
1167  algo->fRadThick[iSta].table[0][0] =
1168  algo->vStations[iSta].materialInfo.RadThick[0];
1169  }
1170  }
1171 
1172  if (fUseTRD)
1173  if (fTrdMatBudgetFileName != "") {
1174  TFile* oldfile = gFile;
1175  TFile* rlFile = new TFile(fTrdMatBudgetFileName);
1176  cout << "TRD Material budget file is " << fTrdMatBudgetFileName << ".\n";
1177  for (int j = 1, iSta = (NStsStations + NMvdStations + NMuchStations);
1178  iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations);
1179  iSta++, j++) {
1180  TString stationNameSts = stationName;
1181  stationNameSts += j;
1182  TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
1183  if (!hStaRadLen) {
1184  cout << "L1: incorrect " << fTrdMatBudgetFileName << " file. No "
1185  << stationNameSts << "\n";
1186  exit(1);
1187  }
1188 
1189 
1190  const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
1191  const float RMax =
1192  hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
1193  algo->fRadThick[iSta].SetBins(NBins, RMax);
1194  algo->fRadThick[iSta].table.resize(NBins);
1195 
1196  for (int iB = 0; iB < NBins; iB++) {
1197  algo->fRadThick[iSta].table[iB].resize(NBins);
1198  float hole = 0.15;
1199  for (int iB2 = 0; iB2 < NBins; iB2++) {
1200  algo->fRadThick[iSta].table[iB][iB2] =
1201  0.01 * hStaRadLen->GetBinContent(iB, iB2);
1202  // Correction for holes in material map
1203  //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
1204 
1205  if (iB2 > 0 && iB2 < NBins - 1)
1206  algo->fRadThick[iSta].table[iB][iB2] =
1207  TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
1208  0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
1209  // Correction for the incorrect harcoded value of RadThick of MVD stations
1210 
1211  if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015)
1212  hole = algo->fRadThick[iSta].table[iB][iB2];
1213  if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015)
1214  algo->fRadThick[iSta].table[iB][iB2] = hole;
1215  // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
1216  }
1217  }
1218  }
1219  rlFile->Close();
1220  rlFile->Delete();
1221  gFile = oldfile;
1222  } else {
1223  LOG(warn) << "No TRD material budget file is found. Homogenious budget "
1224  "will be used";
1225  for (int iSta = (NStsStations + NMvdStations + NMuchStations);
1226  iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations);
1227  iSta++) {
1228  algo->fRadThick[iSta].SetBins(1, 100);
1229  algo->fRadThick[iSta].table.resize(1);
1230  algo->fRadThick[iSta].table[0].resize(1);
1231  algo->fRadThick[iSta].table[0][0] =
1232  algo->vStations[iSta].materialInfo.RadThick[0];
1233  }
1234  }
1235 
1236  if (fUseTOF)
1237  if (fTofMatBudgetFileName != "") {
1238  TFile* oldfile = gFile;
1239  TFile* rlFile = new TFile(fTofMatBudgetFileName);
1240  cout << "TOF Material budget file is " << fTofMatBudgetFileName << ".\n";
1241  for (int j = 1,
1242  iSta =
1243  (NStsStations + NMvdStations + NMuchStations + NTrdStations);
1244  iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations
1245  + NTOFStation);
1246  iSta++, j++) {
1247  TString stationNameSts = stationName;
1248  stationNameSts += j;
1249  TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts);
1250  if (!hStaRadLen) {
1251  cout << "L1: incorrect " << fTofMatBudgetFileName << " file. No "
1252  << stationNameSts << "\n";
1253  exit(1);
1254  }
1255 
1256 
1257  const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
1258  const float RMax =
1259  hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
1260  algo->fRadThick[iSta].SetBins(NBins, RMax);
1261  algo->fRadThick[iSta].table.resize(NBins);
1262 
1263  for (int iB = 0; iB < NBins; iB++) {
1264  algo->fRadThick[iSta].table[iB].resize(NBins);
1265  float hole = 0.0015;
1266  for (int iB2 = 0; iB2 < NBins; iB2++) {
1267  algo->fRadThick[iSta].table[iB][iB2] =
1268  0.01 * hStaRadLen->GetBinContent(iB, iB2);
1269  // Correction for holes in material map
1270  //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
1271 
1272  if (iB2 > 0 && iB2 < NBins - 1)
1273  algo->fRadThick[iSta].table[iB][iB2] =
1274  TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
1275  0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
1276  // Correction for the incorrect harcoded value of RadThick of MVD stations
1277 
1278  if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015)
1279  hole = algo->fRadThick[iSta].table[iB][iB2];
1280  if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015)
1281  algo->fRadThick[iSta].table[iB][iB2] = hole;
1282  // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
1283  }
1284  }
1285  }
1286  rlFile->Close();
1287  rlFile->Delete();
1288  gFile = oldfile;
1289  } else {
1290  LOG(warn) << "No TOF material budget file is found. Homogenious budget "
1291  "will be used";
1292  for (int iSta =
1293  (NStsStations + NMvdStations + NMuchStations + NTrdStations);
1294  iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations
1295  + NTOFStation);
1296  iSta++) {
1297  algo->fRadThick[iSta].SetBins(1, 100);
1298  algo->fRadThick[iSta].table.resize(1);
1299  algo->fRadThick[iSta].table[0].resize(1);
1300  algo->fRadThick[iSta].table[0][0] =
1301  algo->vStations[iSta].materialInfo.RadThick[0];
1302  }
1303  }
1304  return kSUCCESS;
1305 }
1306 
1307 
1308 void CbmL1::Exec(Option_t* /*option*/) {}
1309 
1310 void CbmL1::Reconstruct(CbmEvent* event) {
1311  static int nevent = 0;
1312  vFileEvent.clear();
1313 
1314  if (fTimesliceMode && fPerformance) {
1315 
1316  int nofEvents = fEventList->GetNofEvents();
1317  for (int iE = 0; iE < nofEvents; iE++) {
1318 
1319  int fileId = fEventList->GetFileIdByIndex(iE);
1320  int eventId = fEventList->GetEventIdByIndex(iE);
1321  vFileEvent.insert(DFSET::value_type(fileId, eventId));
1322  }
1323 
1324  } else {
1325  Int_t iFile = FairRunAna::Instance()->GetEventHeader()->GetInputFileId();
1326  Int_t iEvent = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber();
1327  vFileEvent.insert(DFSET::value_type(iFile, iEvent));
1328  }
1329 
1330 
1331  if (fVerbose > 1)
1332  cout << endl << "CbmL1::Exec event " << ++nevent << " ..." << endl << endl;
1333 #ifdef _OPENMP
1334  omp_set_num_threads(1);
1335 #endif
1336  // repack data
1337 
1338  fData->Clear();
1339 
1340  if (fSTAPDataMode >= 2) { // 2,3
1341  fData->ReadHitsFromFile(fSTAPDataDir.Data(), 1, fVerbose);
1342 
1343  algo->SetData(fData->GetStsHits(),
1344  fData->GetStsStrips(),
1345  fData->GetStsStripsB(),
1346  fData->GetStsZPos(),
1347  fData->GetSFlag(),
1348  fData->GetSFlagB(),
1349  fData->GetStsHitsStartIndex(),
1350  fData->GetStsHitsStopIndex());
1351  } else
1352  ReadEvent(fData, event);
1353 
1354  if (0) { // correct hits on MC // dbg
1355  TRandom3 random;
1356  vector<int> sF, sB, zP;
1357  sF.clear();
1358  sB.clear();
1359  zP.clear();
1360  for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) {
1361  L1StsHit& h = const_cast<L1StsHit&>((*algo->vStsHits)[iH]);
1362 #ifdef USE_EVENT_NUMBER
1363  h.n = -1;
1364 #endif
1365  if (vStsHits[iH].mcPointIds.size() == 0) continue;
1366 
1367  const CbmL1MCPoint& mcp = vMCPoints[vStsHits[iH].mcPointIds[0]];
1368 
1369 #ifdef USE_EVENT_NUMBER
1370  h.n = mcp.event;
1371 #endif
1372  const int ista = (*algo->vSFlag)[h.f] / 4;
1373  const L1Station& sta = algo->vStations[ista];
1374  if (std::find(sF.begin(), sF.end(), h.f)
1375  != sF.end()) { // separate strips
1376 
1377  const_cast<vector<unsigned char>*>(algo->vSFlag)
1378  ->push_back((*algo->vSFlag)[h.f]);
1379 
1380  h.f = (*algo->vStsStrips).size();
1381 
1382  const_cast<std::vector<L1Strip>*>(algo->vStsStrips)
1383  ->push_back(L1Strip());
1384  }
1385  sF.push_back(h.f);
1386  if (std::find(sB.begin(), sB.end(), h.b) != sB.end()) {
1387  const_cast<vector<unsigned char>*>(algo->vSFlagB)
1388  ->push_back((*algo->vSFlagB)[h.b]);
1389  h.b = (*algo->vStsStripsB).size();
1390 
1391  const_cast<std::vector<L1Strip>*>(algo->vStsStripsB)
1392  ->push_back(L1Strip());
1393  }
1394  sB.push_back(h.b);
1395  if (std::find(zP.begin(), zP.end(), h.iz)
1396  != zP.end()) { // TODO why do we need it??gives prob=0
1397  h.iz = (*algo->vStsZPos).size();
1398  const_cast<std::vector<float>*>(algo->vStsZPos)->push_back(0);
1399  }
1400  zP.push_back(h.iz);
1401 
1402  const fscal idet = 1
1403  / (sta.xInfo.sin_phi * sta.yInfo.sin_phi
1404  - sta.xInfo.cos_phi * sta.yInfo.cos_phi)[0];
1405 
1406 #if 1 // GAUSS
1407  const_cast<std::vector<L1Strip>*>(algo->vStsStrips)->at(h.f) =
1408  idet * (+sta.yInfo.sin_phi[0] * mcp.x - sta.xInfo.cos_phi[0] * mcp.y)
1409  + random.Gaus(0, sqrt(sta.frontInfo.sigma2)[0]);
1410  const_cast<L1Strip&>((*algo->vStsStripsB)[h.b]) =
1411  idet * (-sta.yInfo.cos_phi[0] * mcp.x + sta.xInfo.sin_phi[0] * mcp.y)
1412  + random.Gaus(0, sqrt(sta.backInfo.sigma2)[0]);
1413 #else // UNIFORM
1414  (*algo->vStsStrips)[h.f] =
1415  idet * (+sta.yInfo.sin_phi[0] * mcp.x - sta.xInfo.cos_phi[0] * mcp.y)
1416  + random.Uniform(-sqrt(sta.frontInfo.sigma2)[0] * 3.5,
1417  sqrt(sta.frontInfo.sigma2)[0] * 3.5);
1418  (*algo->vStsStripsB)[h.b] =
1419  idet * (-sta.yInfo.cos_phi[0] * mcp.x + sta.xInfo.sin_phi[0] * mcp.y)
1420  + random.Uniform(-sqrt(sta.backInfo.sigma2)[0] * 3.5,
1421  sqrt(sta.backInfo.sigma2)[0] * 3.5);
1422 #endif
1423  const_cast<float&>((*algo->vStsZPos)[h.iz]) = mcp.z;
1424  }
1425  }
1426 
1427 
1428  if (fPerformance) {
1429  HitMatch();
1430  // calculate the max number of Hits\mcPoints on continuous(consecutive) stations
1431 
1432  for (vector<CbmL1MCTrack>::iterator it = vMCTracks.begin();
1433  it != vMCTracks.end();
1434  ++it)
1435  it->Init();
1436  }
1437 
1438  if (fSTAPDataMode % 2 == 1) { // 1,3
1439  WriteSTAPAlgoData();
1440  WriteSTAPPerfData();
1441  };
1442  if (fSTAPDataMode >= 2) { // 2,3
1443  //ReadSTAPAlgoData();
1444 
1445  ReadSTAPPerfData();
1446  };
1447 
1448  if ((fPerformance) && (fSTAPDataMode < 2)) { InputPerformance(); }
1449 
1450  // FieldApproxCheck();
1451  // FieldIntegralCheck();
1452 
1453  for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) {
1454 #ifdef USE_EVENT_NUMBER
1455  L1StsHit& h = const_cast<L1StsHit&>((*algo->vStsHits)[iH]);
1456  h.n = -1;
1457 #endif
1458  if (vStsHits[iH].mcPointIds.size() == 0) continue;
1459 #ifdef USE_EVENT_NUMBER
1460  const CbmL1MCPoint& mcp = vMCPoints[vStsHits[iH].mcPointIds[0]];
1461  h.n = mcp.event;
1462 #endif
1463  }
1464 
1465  for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin();
1466  i != vMCTracks.end();
1467  ++i) {
1468  CbmL1MCTrack& MC = *i;
1469 
1470  if (!MC.IsReconstructable()) continue;
1471  if (!(MC.ID >= 0)) continue;
1472 
1473  if (MC.StsHits.size() < 4) continue;
1474  vector<int> hitIndices(algo->NStations, -1);
1475 
1476  for (unsigned int iH = 0; iH < MC.StsHits.size(); iH++) {
1477  const int hitI = MC.StsHits[iH];
1478  CbmL1StsHit& hit = const_cast<CbmL1StsHit&>(vStsHits[hitI]);
1479 
1480  hit.event = MC.iEvent;
1481 
1482  // const int iStation = vMCPoints[hit.mcPointIds[0]].iStation;
1483  // hitIndices[iStation] = hitI;
1484  }
1485  }
1486 
1487 
1488  if (fVerbose > 1) cout << "L1 Track finder..." << endl;
1489  algo->CATrackFinder();
1490  // IdealTrackFinder();
1491 
1492 
1493  if (fVerbose > 1) cout << "L1 Track finder ok" << endl;
1494  // algo->L1KFTrackFitter( fExtrapolateToTheEndOfSTS );
1495 
1496 
1497  if (fmCBMmode || fGlobalMode) {
1498 
1499  L1FieldValue fB0 = algo->GetVtxFieldValue();
1500 
1501  if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001)
1502  && (fabs(fB0.z[0]) < 0.0000001))
1503  algo->KFTrackFitter_simple();
1504 
1505  else
1506  algo->L1KFTrackFitterMuch();
1507  } else {
1508 
1509  L1FieldValue fB0 = algo->GetVtxFieldValue();
1510 
1511  if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001)
1512  && (fabs(fB0.z[0]) < 0.0000001))
1513  algo->KFTrackFitter_simple();
1514 
1515  else
1516  algo->L1KFTrackFitter();
1517  }
1518 
1519 
1520  if (fVerbose > 1) cout << "L1 Track fitter ok" << endl;
1521 
1522  // save recontstructed tracks
1523  vRTracks.clear();
1524  int start_hit = 0;
1525 
1526  for (vector<L1Track>::iterator it = algo->vTracks.begin();
1527  it != (algo->vTracks.begin() + algo->NTracksIsecAll);
1528  it++) {
1529 
1530  CbmL1Track t;
1531  for (int i = 0; i < 7; i++)
1532  t.T[i] = it->TFirst[i];
1533  for (int i = 0; i < 21; i++)
1534  t.C[i] = it->CFirst[i];
1535  for (int i = 0; i < 7; i++)
1536  t.TLast[i] = it->TLast[i];
1537  for (int i = 0; i < 21; i++)
1538  t.CLast[i] = it->CLast[i];
1539  for (int k = 0; k < 7; k++)
1540  t.Tpv[k] = it->Tpv[k];
1541  for (int k = 0; k < 21; k++)
1542  t.Cpv[k] = it->Cpv[k];
1543  t.chi2 = it->chi2;
1544  t.NDF = it->NDF;
1545  //t.T[4] = it->Momentum;
1546  t.StsHits.clear();
1547  t.fTrackTime = it->fTrackTime;
1548 
1549  for (int i = 0; i < it->NHits; i++) {
1550  int start_hit1 = start_hit;
1551 
1552  if (algo->vRecoHits[start_hit1] > vStsHits.size() - 1)
1553  start_hit++;
1554  else
1555  t.StsHits.push_back(algo->vRecoHits[start_hit++]);
1556  }
1557  t.mass = 0.1395679; // pion mass
1558  t.is_electron = 0;
1559 
1560  t.SetId(vRTracks.size());
1561  CbmL1Track* prtra = &t;
1562 
1563  int indd = 0;
1564 
1565  for (vector<int>::iterator ih = (prtra->StsHits).begin();
1566  ih != (prtra->StsHits).end();
1567  ++ih) {
1568  if ((*ih) > int(vStsHits.size() - 1)) {
1569  indd = 1;
1570  break;
1571  }
1572  int nMCPoints = vStsHits[*ih].mcPointIds.size();
1573  for (int iP = 0; iP < nMCPoints; iP++) {
1574  int iMC = vStsHits[*ih].mcPointIds[iP];
1575  if (iMC > int(vMCPoints.size() - 1)) {
1576  // cout << " iMC " << iMC << " vMCPoints.size() " << vMCPoints.size() << endl;
1577  indd = 1;
1578  }
1579  }
1580  }
1581 
1582  if (indd == 1) continue;
1583 
1584  vRTracks.push_back(t);
1585  }
1586  // output performance
1587  if (fPerformance) {
1588  if (fVerbose > 1) cout << "Performance..." << endl;
1589  //HitMatch();
1590  TrackMatch();
1591  }
1592 
1593 
1594  if (fPerformance) {
1595  EfficienciesPerformance();
1596  HistoPerformance();
1597  TrackFitPerformance();
1598  // TimeHist();
1600  }
1601  if (fVerbose > 1) cout << "End of L1" << endl;
1602 
1603  static bool ask = 0;
1604  char symbol;
1605  if (ask) {
1606  std::cout << std::endl << "L1 run (any/r/q) > ";
1607  do {
1608  std::cin.get(symbol);
1609  if (symbol == 'r') ask = false;
1610  if ((symbol == 'e') || (symbol == 'q')) exit(0);
1611  } while (symbol != '\n');
1612  }
1613 }
1614 
1615 // ----- Finish CbmStsFitPerformanceTask task -----------------------------
1616 void CbmL1::Finish() {
1617  TDirectory* curr = gDirectory;
1618  TFile* currentFile = gFile;
1619 
1620  // Open output file and write histograms
1621  TFile* outfile = new TFile("L1_histo.root", "RECREATE");
1622  outfile->cd();
1623  writedir2current(histodir);
1624  outfile->Close();
1625  outfile->Delete();
1626 
1627  gFile = currentFile;
1628  gDirectory = curr;
1629 }
1630 
1631 
1632 void CbmL1::writedir2current(TObject* obj) {
1633  if (!obj->IsFolder())
1634  obj->Write();
1635  else {
1636  TDirectory* cur = gDirectory;
1637  TDirectory* sub = cur->mkdir(obj->GetName());
1638  sub->cd();
1639  TList* listSub = (L1_DYNAMIC_CAST<TDirectory*>(obj))->GetList();
1640  TIter it(listSub);
1641  while (TObject* obj1 = it())
1642  writedir2current(obj1);
1643  cur->cd();
1644  }
1645 }
1646 
1648 
1649 void CbmL1::IdealTrackFinder() {
1650  algo->vTracks.clear();
1651  algo->vRecoHits.clear();
1652 
1653  for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin();
1654  i != vMCTracks.end();
1655  ++i) {
1656  CbmL1MCTrack& MC = *i;
1657 
1658  if (!MC.IsReconstructable()) continue;
1659  if (!(MC.ID >= 0)) continue;
1660 
1661  if (MC.StsHits.size() < 4) continue;
1662 
1663  L1Track algoTr;
1664  algoTr.NHits = 0;
1665 
1666  vector<int> hitIndices(algo->NStations, -1);
1667 
1668  for (unsigned int iH = 0; iH < MC.StsHits.size(); iH++) {
1669  const int hitI = MC.StsHits[iH];
1670  const CbmL1StsHit& hit = vStsHits[hitI];
1671  const int iStation = vMCPoints[hit.mcPointIds[0]].iStation;
1672 
1673  if (iStation >= 0) hitIndices[iStation] = hitI;
1674  }
1675 
1676 
1677  for (int iH = 0; iH < algo->NStations; iH++) {
1678  const int hitI = hitIndices[iH];
1679  if (hitI < 0) continue;
1680 
1681  // algo->vRecoHits.push_back(hitI);
1682  algoTr.NHits++;
1683  }
1684 
1685 
1686  if (algoTr.NHits < 3) continue;
1687 
1688  for (int iH = 0; iH < algo->NStations; iH++) {
1689  const int hitI = hitIndices[iH];
1690  if (hitI < 0) continue;
1691  algo->vRecoHits.push_back(hitI);
1692  }
1693 
1694 
1695  algoTr.Momentum = MC.p;
1696  algoTr.TFirst[0] = MC.x;
1697  algoTr.TFirst[1] = MC.y;
1698  algoTr.TFirst[2] = MC.px / MC.pz;
1699  algoTr.TFirst[3] = MC.py / MC.pz;
1700  algoTr.TFirst[4] = MC.q / MC.p;
1701  algoTr.TFirst[5] = MC.z;
1702 
1703  algo->vTracks.push_back(algoTr);
1704  }
1705  algo->NTracksIsecAll = algo->vTracks.size();
1706 }; // void CbmL1::IdealTrackFinder()
1707 
1708 
1710 
1711 void CbmL1::WriteSTAPGeoData(const vector<float>& geo_) {
1712  // write geo in file
1713  TString fgeo_name = fSTAPDataDir + "geo_algo.txt";
1714  ofstream fgeo(fgeo_name);
1715  fgeo.setf(ios::scientific, ios::floatfield);
1716  fgeo.precision(20);
1717  int size = geo_.size();
1718  for (int i = 0; i < size; i++) {
1719  fgeo << geo_[i] << endl;
1720  };
1721  fgeo.close();
1722  cout << "-I- CbmL1: Geometry data has been written in " << fgeo_name << endl;
1723 } // void CbmL1::WriteSTAPGeoData(void* geo_, int size)
1724 
1725 
1726 void CbmL1::WriteSTAPAlgoData() // must be called after ReadEvent
1727 {
1728  // write algo data in file
1729  static int vNEvent = 1;
1730  fstream fadata;
1731 
1732  TString fadata_name = fSTAPDataDir + "data_algo.txt";
1733  // if ( vNEvent <= maxNEvent ) {
1734  if (1) {
1735 
1736  if (vNEvent == 1)
1737  fadata.open(fadata_name, fstream::out); // begin new file
1738  else
1739  fadata.open(fadata_name, fstream::out | fstream::app);
1740 
1741  fadata << "Event:"
1742  << " ";
1743  fadata << vNEvent << endl;
1744  // write vStsStrips
1745  int n = (*algo->vStsStrips).size(); // number of elements
1746  fadata << n << endl;
1747  for (int i = 0; i < n; i++) {
1748  fadata << (*algo->vStsStrips)[i] << endl;
1749  };
1750  if (fVerbose >= 4)
1751  cout << "vStsStrips[" << n << "]"
1752  << " have been written." << endl;
1753  // write vStsStripsB
1754  n = (*algo->vStsStripsB).size();
1755  fadata << n << endl;
1756  for (int i = 0; i < n; i++) {
1757  fadata << (*algo->vStsStripsB)[i] << endl;
1758  };
1759  if (fVerbose >= 4)
1760  cout << "vStsStripsB[" << n << "]"
1761  << " have been written." << endl;
1762  // write vStsZPos
1763  n = (*algo->vStsZPos).size();
1764  fadata << n << endl;
1765  for (int i = 0; i < n; i++) {
1766  fadata << (*algo->vStsZPos)[i] << endl;
1767  };
1768  if (fVerbose >= 4)
1769  cout << "vStsZPos[" << n << "]"
1770  << " have been written." << endl;
1771  // write vSFlag
1772  n = (*algo->vSFlag).size();
1773  fadata << n << endl;
1774  unsigned char element;
1775  for (int i = 0; i < n; i++) {
1776  element = (*algo->vSFlag)[i];
1777  fadata << static_cast<int>(element) << endl;
1778  };
1779  if (fVerbose >= 4)
1780  cout << "vSFlag[" << n << "]"
1781  << " have been written." << endl;
1782  // write vSFlagB
1783  n = (*algo->vSFlagB).size();
1784  fadata << n << endl;
1785  for (int i = 0; i < n; i++) {
1786  element = (*algo->vSFlagB)[i];
1787  fadata << static_cast<int>(element) << endl;
1788  };
1789  if (fVerbose >= 4)
1790  cout << "vSFlagB[" << n << "]"
1791  << " have been written." << endl;
1792  // write vStsHits
1793  n = (*algo->vStsHits).size();
1794  fadata << n << endl;
1795  for (int i = 0; i < n; i++) {
1796  fadata << static_cast<int>((*algo->vStsHits)[i].f) << " ";
1797  fadata << static_cast<int>((*algo->vStsHits)[i].b) << " ";
1798 #ifdef USE_EVENT_NUMBER
1799  fadata << static_cast<unsigned short int>((*algo->vStsHits)[i].n) << " ";
1800 #endif
1801  fadata << static_cast<int>((*algo->vStsHits)[i].iz) << " ";
1802  // fadata << (*algo->vStsHits)[i].time << endl;
1803 
1804  fadata << (*algo->vStsHits)[i].t_reco << endl;
1805  };
1806  if (fVerbose >= 4)
1807  cout << "vStsHits[" << n << "]"
1808  << " have been written." << endl;
1809  // write StsHitsStartIndex and StsHitsStopIndex
1810  n = 20;
1811  for (int i = 0; i < n; i++) {
1812  if (int(algo->MaxNStations) + 1 > i)
1813  fadata << algo->StsHitsStartIndex[i] << endl;
1814  else
1815  fadata << 0 << endl;
1816  };
1817  for (int i = 0; i < n; i++) {
1818  if (int(algo->MaxNStations) + 1 > i)
1819  fadata << algo->StsHitsStopIndex[i] << endl;
1820  else
1821  fadata << 0 << endl;
1822  };
1823 
1824 
1825  fadata.close();
1826  }
1827  cout << "-I- CbmL1: CATrackFinder data for event number " << vNEvent
1828  << " have been written in file " << fadata_name << endl;
1829  vNEvent++;
1830 } // void CbmL1::WriteSTAPAlgoData()
1831 
1832 
1833 void CbmL1::WriteSTAPPerfData() // must be called after ReadEvent
1834 {
1835  fstream fpdata;
1836  fpdata << setprecision(8);
1837 
1838  static int vNEvent = 1;
1839 
1840  TString fpdata_name = fSTAPDataDir + "data_perfo.txt";
1841  // write data for performance in file
1842  // if ( vNEvent <= maxNEvent ) {
1843  if (1) {
1844 
1845  if (vNEvent == 1)
1846  fpdata.open(fpdata_name, fstream::out); // begin new file
1847  else
1848  fpdata.open(fpdata_name, fstream::out | fstream::app);
1849 
1850  fpdata << "Event: ";
1851  fpdata << vNEvent << endl;
1852  // write vMCPoints
1853  Int_t n = vMCPoints.size(); // number of elements
1854  fpdata << n << endl;
1855  for (int i = 0; i < n; i++) {
1856  fpdata << vMCPoints[i].xIn << " ";
1857  fpdata << vMCPoints[i].yIn << " ";
1858  fpdata << vMCPoints[i].zIn << " ";
1859  fpdata << vMCPoints[i].pxIn << " ";
1860  fpdata << vMCPoints[i].pyIn << " ";
1861  fpdata << vMCPoints[i].pzIn << " " << endl;
1862  fpdata << vMCPoints[i].xOut << " ";
1863  fpdata << vMCPoints[i].yOut << " ";
1864  fpdata << vMCPoints[i].zOut << " ";
1865  fpdata << vMCPoints[i].pxOut << " ";
1866  fpdata << vMCPoints[i].pyOut << " ";
1867  fpdata << vMCPoints[i].pzOut << " " << endl;
1868 
1869  fpdata << vMCPoints[i].p << " ";
1870  fpdata << vMCPoints[i].q << " ";
1871  fpdata << vMCPoints[i].mass << " ";
1872  fpdata << vMCPoints[i].time << " ";
1873 
1874  fpdata << vMCPoints[i].pdg << " ";
1875  fpdata << vMCPoints[i].ID << " ";
1876  fpdata << vMCPoints[i].mother_ID << " ";
1877  fpdata << vMCPoints[i].iStation << endl;
1878 
1879  const int nhits = vMCPoints[i].hitIds.size();
1880  fpdata << nhits << endl << " ";
1881  for (int k = 0; k < nhits; k++) {
1882  fpdata << vMCPoints[i].hitIds[k] << " ";
1883  };
1884  fpdata << endl;
1885  };
1886  if (fVerbose >= 4)
1887  cout << "vMCPoints[" << n << "]"
1888  << " have been written." << endl;
1889 
1890  // write vMCTracks . without Points
1891  n = vMCTracks.size(); // number of elements
1892  fpdata << n << endl;
1893  for (int i = 0; i < n; i++) {
1894  fpdata << vMCTracks[i].x << " ";
1895  fpdata << vMCTracks[i].y << " ";
1896  fpdata << vMCTracks[i].z << " ";
1897  fpdata << vMCTracks[i].px << " ";
1898  fpdata << vMCTracks[i].py << " ";
1899  fpdata << vMCTracks[i].pz << " ";
1900  fpdata << vMCTracks[i].p << " ";
1901  fpdata << vMCTracks[i].q << " ";
1902  fpdata << vMCTracks[i].mass << " ";
1903  fpdata << vMCTracks[i].time << " ";
1904 
1905  fpdata << vMCTracks[i].pdg << " ";
1906  fpdata << vMCTracks[i].ID << " ";
1907  fpdata << vMCTracks[i].mother_ID << endl;
1908 
1909  int nhits = vMCTracks[i].StsHits.size();
1910  fpdata << " " << nhits << endl << " ";
1911  for (int k = 0; k < nhits; k++) {
1912  fpdata << vMCTracks[i].StsHits[k] << " ";
1913  };
1914  fpdata << endl;
1915 
1916  const int nPoints = vMCTracks[i].Points.size();
1917  fpdata << nPoints << endl << " ";
1918  for (int k = 0; k < nPoints; k++) {
1919  fpdata << vMCTracks[i].Points[k] << " ";
1920  };
1921  fpdata << endl;
1922 
1923  fpdata << vMCTracks[i].nMCContStations << " ";
1924  fpdata << vMCTracks[i].nHitContStations << " ";
1925  fpdata << vMCTracks[i].maxNStaMC << " ";
1926  fpdata << vMCTracks[i].maxNSensorMC << " ";
1927  fpdata << vMCTracks[i].maxNStaHits << " ";
1928  fpdata << vMCTracks[i].nStations << endl;
1929  };
1930  if (fVerbose >= 4)
1931  cout << "vMCTracks[" << n << "]"
1932  << " have been written." << endl;
1933 
1934  // write vHitMCRef
1935  n = vHitMCRef.size(); // number of elements
1936  fpdata << n << endl;
1937  for (int i = 0; i < n; i++) {
1938  fpdata << vHitMCRef[i] << endl;
1939  };
1940  if (fVerbose >= 4)
1941  cout << "vHitMCRef[" << n << "]"
1942  << " have been written." << endl;
1943 
1944  // write vHitStore
1945  n = vHitStore.size(); // number of elements
1946  fpdata << n << endl;
1947  for (int i = 0; i < n; i++) {
1948  fpdata << vHitStore[i].ExtIndex << " ";
1949  fpdata << vHitStore[i].iStation << " ";
1950 
1951  fpdata << vHitStore[i].x << " ";
1952  fpdata << vHitStore[i].y << endl;
1953  };
1954  if (fVerbose >= 4)
1955  cout << "vHitStore[" << n << "]"
1956  << " have been written." << endl;
1957 
1958  // write vStsHits
1959  n = vStsHits.size(); // number of elements
1960  fpdata << n << endl;
1961  for (int i = 0; i < n; i++) {
1962  fpdata << vStsHits[i].hitId << " ";
1963  fpdata << vStsHits[i].extIndex << endl;
1964 
1965  const int nPoints = vStsHits[i].mcPointIds.size();
1966  fpdata << nPoints << endl << " ";
1967  for (int k = 0; k < nPoints; k++) {
1968  fpdata << vStsHits[i].mcPointIds[k] << " ";
1969  };
1970  fpdata << endl;
1971  };
1972  if (fVerbose >= 4)
1973  cout << "vStsHits[" << n << "]"
1974  << " have been written." << endl;
1975 
1976  fpdata.close();
1977  }
1978  cout << "-I- CbmL1: Data for performance of event number " << vNEvent
1979  << " have been written in file " << fpdata_name << endl;
1980  vNEvent++;
1981 } // void CbmL1::WriteSTAPPerfData()
1982 
1983 istream& CbmL1::eatwhite(istream& is) // skip spaces
1984 {
1985  char c;
1986  while (is.get(c)) {
1987  if (isspace(c) == 0) {
1988  is.putback(c);
1989  break;
1990  }
1991  }
1992  return is;
1993 }
1994 
1995 //void CbmL1::ReadSTAPGeoData(vector<float> geo_, int &size)
1996 //void CbmL1::ReadSTAPGeoData(vector<fscal> geo_, int &size)
1997 void CbmL1::ReadSTAPGeoData(vector<fscal>& geo_, int& size) {
1998  TString fgeo_name = fSTAPDataDir + "geo_algo.txt";
1999  ifstream fgeo(fgeo_name);
2000 
2001  cout << "-I- CbmL1: Read geometry from file " << fgeo_name << endl;
2002  int i;
2003  for (i = 0; !fgeo.eof(); i++) {
2004  fscal tmp;
2005  fgeo >> tmp >> eatwhite;
2006  cout << " geo_[" << i << "]=" << geo_[i] << " tmp= " << tmp << endl;
2007  geo_[i] = tmp;
2008  };
2009  size = i;
2010  fgeo.close();
2011 } // void CbmL1::ReadSTAPGeoData(void* geo_, int &size)
2012 
2013 void CbmL1::ReadSTAPAlgoData() {
2014  static int nEvent = 1;
2015  static fstream fadata;
2016  TString fadata_name = fSTAPDataDir + "data_algo.txt";
2017  // if (nEvent <= maxNEvent){
2018  if (1) {
2019  if (nEvent == 1) fadata.open(fadata_name, fstream::in);
2020 
2021  if (algo->vStsHits)
2022  const_cast<std::vector<L1StsHit>*>(algo->vStsHits)->clear();
2023  if (algo->vStsStrips)
2024  const_cast<std::vector<L1Strip>*>(algo->vStsStrips)->clear();
2025  if (algo->vStsStripsB)
2026  const_cast<std::vector<L1Strip>*>(algo->vStsStripsB)->clear();
2027  if (algo->vStsZPos)
2028  const_cast<std::vector<float>*>(algo->vStsZPos)->clear();
2029  if (algo->vSFlag) const_cast<vector<unsigned char>*>(algo->vSFlag)->clear();
2030  if (algo->vSFlagB)
2031  const_cast<vector<unsigned char>*>(algo->vSFlagB)->clear();
2032 
2033  // check correct position in file
2034  char s[] = "Event: ";
2035  int nEv;
2036  fadata >> s;
2037  fadata >> nEv;
2038  if (nEv != nEvent)
2039  cout << "-E- CbmL1: Can't read event number " << nEvent << " from file "
2040  << fadata_name << endl;
2041 
2042  int n; // number of elements
2043  // read algo->vStsStrips
2044  fadata >> n;
2045  cout << " n " << n << endl;
2046  for (int i = 0; i < n; i++) {
2047  fscal element;
2048  fadata >> element;
2049  const_cast<std::vector<L1Strip>*>(algo->vStsStrips)->push_back(element);
2050  }
2051  if (fVerbose >= 4)
2052  cout << "vStsStrips[" << n << "]"
2053  << " have been read." << endl;
2054  // read algo->vStsStripsB
2055  fadata >> n;
2056  for (int i = 0; i < n; i++) {
2057  fscal element;
2058  fadata >> element;
2059  const_cast<std::vector<L1Strip>*>(algo->vStsStripsB)->push_back(element);
2060  }
2061  if (fVerbose >= 4)
2062  cout << "vStsStripsB[" << n << "]"
2063  << " have been read." << endl;
2064  // read algo->vStsZPos
2065  fadata >> n;
2066  for (int i = 0; i < n; i++) {
2067  fscal element;
2068  fadata >> element;
2069  const_cast<std::vector<float>*>(algo->vStsZPos)->push_back(element);
2070  }
2071  if (fVerbose >= 4)
2072  cout << "vStsZPos[" << n << "]"
2073  << " have been read." << endl;
2074  // read algo->vSFlag
2075  fadata >> n;
2076  for (int i = 0; i < n; i++) {
2077  int element;
2078  fadata >> element;
2079  const_cast<vector<unsigned char>*>(algo->vSFlag)
2080  ->push_back(static_cast<unsigned char>(element));
2081  }
2082  if (fVerbose >= 4)
2083  cout << "vSFlag[" << n << "]"
2084  << " have been read." << endl;
2085  // read algo->vSFlagB
2086  fadata >> n;
2087  for (int i = 0; i < n; i++) {
2088  int element;
2089  fadata >> element;
2090  const_cast<std::vector<L1Strip>*>(algo->vStsStripsB)
2091  ->push_back(static_cast<unsigned char>(element));
2092  }
2093  if (fVerbose >= 4)
2094  cout << "vSFlagB[" << n << "]"
2095  << " have been read." << endl;
2096  // read algo->vStsHits
2097  fadata >> n;
2098  int element_f; // for convert
2099  int element_b;
2100  int element_n;
2101  int element_iz;
2102  float time;
2103  for (int i = 0; i < n; i++) {
2104  L1StsHit element;
2105  fadata >> element_f >> element_b >> element_n >> element_iz >> time;
2106  element.f = static_cast<THitI>(element_f);
2107  element.b = static_cast<THitI>(element_b);
2108  element.iz = static_cast<TZPosI>(element_iz);
2109 
2110  element.t_reco = time;
2111  const_cast<std::vector<L1StsHit>*>(algo->vStsHits)->push_back(element);
2112  }
2113  if (fVerbose >= 4)
2114  cout << "vStsHits[" << n << "]"
2115  << " have been read." << endl;
2116  // read StsHitsStartIndex and StsHitsStopIndex
2117  n = 20;
2118  for (int i = 0; i < n; i++) {
2119  int tmp;
2120  fadata >> tmp;
2121  if (int(algo->MaxNStations) + 1 > i)
2122  (const_cast<unsigned int&>(algo->StsHitsStartIndex[i]) = tmp);
2123  }
2124  for (int i = 0; i < n; i++) {
2125  int tmp;
2126  fadata >> tmp;
2127  if (int(algo->MaxNStations) + 1 > i)
2128  (const_cast<unsigned int&>(algo->StsHitsStopIndex[i]) = tmp);
2129  }
2130 
2131  cout << "-I- CbmL1: CATrackFinder data for event " << nEvent
2132  << " has been read from file " << fadata_name << " successfully."
2133  << endl;
2134  }
2135  nEvent++;
2136 } // void CbmL1::ReadSTAPAlgoData()
2137 
2138 void CbmL1::ReadSTAPPerfData() {
2139  static int nEvent = 1;
2140  static fstream fpdata;
2141  TString fpdata_name = fSTAPDataDir + "data_perfo.txt";
2142  // if (nEvent <= maxNEvent){
2143  if (1) {
2144  if (nEvent == 1) { fpdata.open(fpdata_name, fstream::in); };
2145 
2146  vMCPoints.clear();
2147  vMCTracks.clear();
2148  vHitMCRef.clear();
2149  vHitStore.clear();
2150  vStsHits.clear();
2151  dFEI2vMCPoints.clear();
2152  dFEI2vMCTracks.clear();
2153  // check if it is right position in file
2154  char s[] = "EVENT: "; // buffer
2155  int nEv = 0; // event number
2156  fpdata >> s;
2157  fpdata >> nEv;
2158 
2159  if (nEv != nEvent)
2160  cout << "-E- CbmL1: Performance: can't read event number " << nEvent
2161  << " from file "
2162  << "data_perfo.txt" << endl;
2163  // vMCPoints
2164  int n; // number of elements
2165  fpdata >> n;
2166  for (int i = 0; i < n; i++) {
2167  CbmL1MCPoint element;
2168 
2169  fpdata >> element.xIn;
2170  fpdata >> element.yIn;
2171  fpdata >> element.zIn;
2172  fpdata >> element.pxIn;
2173  fpdata >> element.pyIn;
2174  fpdata >> element.pzIn;
2175 
2176  fpdata >> element.xOut;
2177  fpdata >> element.yOut;
2178  fpdata >> element.zOut;
2179  fpdata >> element.pxOut;
2180  fpdata >> element.pyOut;
2181  fpdata >> element.pzOut;
2182 
2183  fpdata >> element.p;
2184  fpdata >> element.q;
2185  fpdata >> element.mass;
2186  fpdata >> element.time;
2187 
2188  fpdata >> element.pdg;
2189  fpdata >> element.ID;
2190  fpdata >> element.mother_ID;
2191  fpdata >> element.iStation;
2192 
2193  int nhits;
2194  fpdata >> nhits;
2195  for (int k = 0; k < nhits; k++) {
2196  int helement;
2197  fpdata >> helement;
2198  element.hitIds.push_back(helement);
2199  };
2200 
2201  vMCPoints.push_back(element);
2202  };
2203  if (fVerbose >= 4)
2204  cout << "vMCPoints[" << n << "]"
2205  << " have been read." << endl;
2206 
2207  // vMCTracks . without Points
2208  fpdata >> n;
2209  for (int i = 0; i < n; i++) {
2210  CbmL1MCTrack element;
2211 
2212  fpdata >> element.x;
2213  fpdata >> element.y;
2214  fpdata >> element.z;
2215  fpdata >> element.px;
2216  fpdata >> element.py;
2217  fpdata >> element.pz;
2218  fpdata >> element.p;
2219  fpdata >> element.q;
2220  fpdata >> element.mass;
2221  fpdata >> element.time;
2222 
2223  fpdata >> element.pdg;
2224  fpdata >> element.ID;
2225  fpdata >> element.mother_ID;
2226 
2227  int nhits;
2228  fpdata >> nhits;
2229  for (int k = 0; k < nhits; k++) {
2230  int helement;
2231  fpdata >> helement;
2232  element.StsHits.push_back(helement);
2233  };
2234  fpdata >> nhits;
2235  for (int k = 0; k < nhits; k++) {
2236  int helement;
2237  fpdata >> helement;
2238  element.Points.push_back(helement);
2239  };
2240 
2241  fpdata >> element.nMCContStations;
2242  fpdata >> element.nHitContStations;
2243  fpdata >> element.maxNStaMC;
2244  fpdata >> element.maxNSensorMC;
2245  fpdata >> element.maxNStaHits;
2246  fpdata >> element.nStations;
2247 
2248  element.CalculateIsReconstructable();
2249  vMCTracks.push_back(element);
2250  };
2251  if (fVerbose >= 4)
2252  cout << "vMCTracks[" << n << "]"
2253  << " have been read." << endl;
2254 
2255  // vHitMCRef
2256  fpdata >> n;
2257  for (int i = 0; i < n; i++) {
2258  int element;
2259  fpdata >> element;
2260  vHitMCRef.push_back(element);
2261  };
2262  if (fVerbose >= 4)
2263  cout << "vHitMCRef[" << n << "]"
2264  << " have been read." << endl;
2265 
2266  // vHitStore
2267  fpdata >> n;
2268  for (int i = 0; i < n; i++) {
2269  CbmL1HitStore element;
2270  fpdata >> element.ExtIndex;
2271  fpdata >> element.iStation;
2272 
2273  fpdata >> element.x;
2274  fpdata >> element.y;
2275  vHitStore.push_back(element);
2276  };
2277  if (fVerbose >= 4)
2278  cout << "vHitStore[" << n << "]"
2279  << " have been read." << endl;
2280 
2281  // vStsHits
2282  fpdata >> n;
2283  for (int i = 0; i < n; i++) {
2284  CbmL1StsHit element;
2285  fpdata >> element.hitId;
2286  fpdata >> element.extIndex;
2287 
2288  int nPoints;
2289  fpdata >> nPoints;
2290  for (int k = 0; k < nPoints; k++) {
2291  int id;
2292  fpdata >> id;
2293  element.mcPointIds.push_back(id);
2294  };
2295  vStsHits.push_back(element);
2296  };
2297  if (fVerbose >= 4)
2298  cout << "vStsHits[" << n << "]"
2299  << " have been read." << endl;
2300 
2301 
2302  // if (nEvent == maxNEvent) { // file open on begin of all work class and close at end
2303  // fpdata.close();
2304  // cout << " -I- Performance: data read from file " << "data_perfo.txt" << " successfully"<< endl;
2305  // }
2306  cout << "-I- CbmL1: L1Performance data for event " << nEvent
2307  << " has been read from file " << fpdata_name << " successfully."
2308  << endl;
2309 
2310  } // if (nEvent <= maxNEvent)
2311  nEvent++;
2312 } // void CbmL1::ReadSTAPPerfData()
2313 
2314 void CbmL1::WriteSIMDKFData() {
2315  static bool first = 1;
2316 
2318  if (first) {
2319  FairField* dMF = CbmKF::Instance()->GetMagneticField();
2320 
2321  fstream FileGeo;
2322  FileGeo.open("geo.dat", ios::out);
2323 
2324  fstream FieldCheck;
2325  FieldCheck.open("field.dat", ios::out);
2326 
2327  Double_t bfg[3], rfg[3];
2328 
2329  rfg[0] = 0.;
2330  rfg[1] = 0.;
2331  rfg[2] = 0.;
2332  dMF->GetFieldValue(rfg, bfg);
2333  FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " "
2334  << endl;
2335 
2336  rfg[0] = 0.;
2337  rfg[1] = 0.;
2338  rfg[2] = 2.5;
2339  dMF->GetFieldValue(rfg, bfg);
2340  FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " "
2341  << endl;
2342 
2343  rfg[0] = 0.;
2344  rfg[1] = 0.;
2345  rfg[2] = 5.0;
2346  dMF->GetFieldValue(rfg, bfg);
2347  FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " "
2348  << endl
2349  << endl;
2350  FileGeo << NStation << endl;
2351 
2352  const int M = 5; // polinom order
2353  const int N = (M + 1) * (M + 2) / 2;
2354 
2355  for (Int_t ist = 0; ist < NStation; ist++) {
2356  fscal f_phi, f_sigma, b_phi, b_sigma;
2357 
2358  double C[3][N];
2359  double z = 0;
2360  double Xmax, Ymax;
2361  if (ist < NMvdStations) {
2362  CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
2363  f_phi = 0;
2364  f_sigma = 5.e-4;
2365  b_phi = 3.14159265358 / 2.;
2366  b_sigma = 5.e-4;
2367  z = t.z;
2368  Xmax = Ymax = t.R;
2369  } else {
2370  CbmStsStation* station =
2371  CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
2372  f_phi = station->GetSensorRotation();
2373  b_phi = f_phi;
2374  double Pi = 3.14159265358;
2375  f_phi += station->GetSensorStereoAngle(0) * Pi / 180.;
2376  b_phi += station->GetSensorStereoAngle(1) * Pi / 180.;
2377  f_sigma = station->GetSensorPitch(0) / TMath::Sqrt(12);
2378  b_sigma = f_sigma;
2379  //f_sigma *= cos(f_phi); // TODO: think about this
2380  //b_sigma *= cos(b_phi);
2381  z = station->GetZ();
2382 
2383  Xmax = station->GetXmax();
2384  Ymax = station->GetYmax();
2385  }
2386 
2387  double dx = 1.; // step for the field approximation
2388  double dy = 1.;
2389 
2390  if (dx > Xmax / N / 2) dx = Xmax / N / 4.;
2391  if (dy > Ymax / N / 2) dy = Ymax / N / 4.;
2392 
2393  for (int i = 0; i < 3; i++)
2394  for (int k = 0; k < N; k++)
2395  C[i][k] = 0;
2396  TMatrixD A(N, N);
2397  TVectorD b0(N), b1(N), b2(N);
2398  for (int i = 0; i < N; i++) {
2399  for (int j = 0; j < N; j++)
2400  A(i, j) = 0.;
2401  b0(i) = b1(i) = b2(i) = 0.;
2402  }
2403  for (double x = -Xmax; x <= Xmax; x += dx)
2404  for (double y = -Ymax; y <= Ymax; y += dy) {
2405  double r = sqrt(fabs(x * x / Xmax / Xmax + y / Ymax * y / Ymax));
2406  if (r > 1.) continue;
2407  Double_t w = 1. / (r * r + 1);
2408  Double_t p[3] = {x, y, z};
2409  Double_t B[3] = {0., 0., 0.};
2410  CbmKF::Instance()->GetMagneticField()->GetFieldValue(p, B);
2411  TVectorD m(N);
2412  m(0) = 1;
2413  for (int i = 1; i <= M; i++) {
2414  int k = (i - 1) * (i) / 2;
2415  int l = i * (i + 1) / 2;
2416  for (int j = 0; j < i; j++)
2417  m(l + j) = x * m(k + j);
2418  m(l + i) = y * m(k + i - 1);
2419  }
2420 
2421  TVectorD mt = m;
2422  for (int i = 0; i < N; i++) {
2423  for (int j = 0; j < N; j++)
2424  A(i, j) += w * m(i) * m(j);
2425  b0(i) += w * B[0] * m(i);
2426  b1(i) += w * B[1] * m(i);
2427  b2(i) += w * B[2] * m(i);
2428  }
2429  }
2430  double det;
2431  A.Invert(&det);
2432  TVectorD c0 = A * b0, c1 = A * b1, c2 = A * b2;
2433  for (int i = 0; i < N; i++) {
2434  C[0][i] = c0(i);
2435  C[1][i] = c1(i);
2436  C[2][i] = c2(i);
2437  }
2438 
2439  double c_f = cos(f_phi);
2440  double s_f = sin(f_phi);
2441  double c_b = cos(b_phi);
2442  double s_b = sin(b_phi);
2443 
2444  double det_m = c_f * s_b - s_f * c_b;
2445  det_m *= det_m;
2446  // double C00 = ( s_b*s_b*f_sigma*f_sigma + s_f*s_f*b_sigma*b_sigma )/det_m;
2447  // double C10 =-( s_b*c_b*f_sigma*f_sigma + s_f*c_f*b_sigma*b_sigma )/det_m;
2448  // double C11 = ( c_b*c_b*f_sigma*f_sigma + c_f*c_f*b_sigma*b_sigma )/det_m;
2449 
2450  // float delta_x = sqrt(C00);
2451  // float delta_y = sqrt(C11);
2452  FileGeo << " " << ist << " ";
2453  if (ist < NMvdStations) {
2454  CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[ist];
2455  FileGeo << t.z << " ";
2456  FileGeo << t.dz << " ";
2457  FileGeo << t.RadLength << " ";
2458  } else if (ist < (NStsStations + NMvdStations)) {
2459  CbmStsStation* station =
2460  CbmStsSetup::Instance()->GetStation(ist - NMvdStations);
2461  FileGeo << station->GetZ() << " ";
2462  FileGeo << station->GetSensorD() << " ";
2463  FileGeo << station->GetRadLength() << " ";
2464  }
2465  FileGeo << f_sigma << " ";
2466  FileGeo << b_sigma << " ";
2467  FileGeo << f_phi << " ";
2468  FileGeo << b_phi << endl;
2469  FileGeo << " " << N << endl;
2470  FileGeo << " ";
2471  for (int ik = 0; ik < N; ik++)
2472  FileGeo << C[0][ik] << " ";
2473  FileGeo << endl;
2474  FileGeo << " ";
2475  for (int ik = 0; ik < N; ik++)
2476  FileGeo << C[1][ik] << " ";
2477  FileGeo << endl;
2478  FileGeo << " ";
2479  for (int ik = 0; ik < N; ik++)
2480  FileGeo << C[2][ik] << " ";
2481  FileGeo << endl;
2482  }
2483  FileGeo.close();
2484  }
2485 
2487 
2488  static int TrNumber = 0;
2489  fstream Tracks, McTracksCentr, McTracksIn, McTracksOut;
2490  if (first) {
2491  Tracks.open("tracks.dat", fstream::out);
2492  McTracksCentr.open("mctrackscentr.dat", fstream::out);
2493  McTracksIn.open("mctracksin.dat", fstream::out);
2494  McTracksOut.open("mctracksout.dat", fstream::out);
2495  first = 0;
2496  } else {
2497  Tracks.open("tracks.dat", fstream::out | fstream::app);
2498  McTracksCentr.open("mctrackscentr.dat", fstream::out | fstream::app);
2499  McTracksIn.open("mctracksin.dat", fstream::out | fstream::app);
2500  McTracksOut.open("mctracksout.dat", fstream::out | fstream::app);
2501  }
2502 
2503  for (vector<CbmL1Track>::iterator RecTrack = vRTracks.begin();
2504  RecTrack != vRTracks.end();
2505  ++RecTrack) {
2506  if (RecTrack->IsGhost()) continue;
2507 
2508  CbmL1MCTrack* MCTrack = RecTrack->GetMCTrack();
2509  if (!(MCTrack->IsPrimary())) continue;
2510 
2511  int NHits = (RecTrack->StsHits).size();
2512  float x[20], y[20], z[20];
2513  int st[20];
2514  int jHit = 0;
2515  for (int iHit = 0; iHit < NHits; iHit++) {
2516  CbmL1HitStore& h = vHitStore[RecTrack->StsHits[iHit]];
2517  st[jHit] = h.iStation;
2518  if (h.ExtIndex < 0) {
2519  CbmMvdHit* MvdH = (CbmMvdHit*) listMvdHits->At(-h.ExtIndex - 1);
2520  x[jHit] = MvdH->GetX();
2521  y[jHit] = MvdH->GetY();
2522  z[jHit] = MvdH->GetZ();
2523  jHit++;
2524  } else {
2525  CbmStsHit* StsH = (CbmStsHit*) listStsHits->At(h.ExtIndex);
2526  x[jHit] = StsH->GetX();
2527  y[jHit] = StsH->GetY();
2528  z[jHit] = StsH->GetZ();
2529  jHit++;
2530  }
2531  }
2532 
2533  Tracks << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " "
2534  << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " "
2535  << MCTrack->pz << " " << MCTrack->q << " " << NHits << endl;
2536 
2537  float AngleXAxis = 0, AngleYAxis = 0;
2538  for (int i = 0; i < NHits; i++)
2539  Tracks << " " << st[i] << " " << x[i] << " " << y[i] << " " << z[i]
2540  << " " << AngleXAxis << " " << AngleYAxis << endl;
2541  Tracks << endl;
2542 
2543  int NMCPoints = (MCTrack->Points).size();
2544 
2545  McTracksIn << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " "
2546  << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " "
2547  << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl;
2548  McTracksOut << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " "
2549  << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " "
2550  << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl;
2551  McTracksCentr << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " "
2552  << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py
2553  << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints
2554  << endl;
2555 
2556  for (int iPoint = 0; iPoint < NMCPoints; iPoint++) {
2557  CbmL1MCPoint& MCPoint = vMCPoints[MCTrack->Points[iPoint]];
2558  McTracksIn << " " << MCPoint.iStation << " " << MCPoint.xIn << " "
2559  << MCPoint.yIn << " " << MCPoint.zIn << " " << MCPoint.pxIn
2560  << " " << MCPoint.pyIn << " " << MCPoint.pzIn << endl;
2561  McTracksOut << " " << MCPoint.iStation << " " << MCPoint.xOut << " "
2562  << MCPoint.yOut << " " << MCPoint.zOut << " " << MCPoint.pxOut
2563  << " " << MCPoint.pyOut << " " << MCPoint.pzOut << endl;
2564  McTracksCentr << " " << MCPoint.iStation << " " << MCPoint.x << " "
2565  << MCPoint.y << " " << MCPoint.z << " " << MCPoint.px << " "
2566  << MCPoint.py << " " << MCPoint.pz << endl;
2567  }
2568  McTracksIn << endl;
2569  McTracksOut << endl;
2570  McTracksCentr << endl;
2571 
2572  TrNumber++;
2573  }
2574 }
CbmMvdDetector::Instance
static CbmMvdDetector * Instance()
Definition: CbmMvdDetector.cxx:47
L1Algo::vStsZPos
const vector< fscal > * vStsZPos
Definition: L1Algo.h:342
CbmL1::fTofHits
TClonesArray * fTofHits
Definition: CbmL1.h:314
CbmTofCell::GetZ
Double_t GetZ() const
Definition: CbmTofCell.h:38
CbmL1::fStsParSetModule
CbmStsParSetModule * fStsParSetModule
Definition: CbmL1.h:372
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmL1::fUseHitErrors
bool fUseHitErrors
Definition: CbmL1.h:122
CbmStsStation::GetZ
Double_t GetZ() const
Definition: CbmStsStation.h:127
L1Algo::vStsHits
const vector< L1StsHit > * vStsHits
Definition: L1Algo.h:344
CbmL1MCPoint::q
double q
Definition: CbmL1MCPoint.h:59
CbmL1::listTrdHits
TClonesArray * listTrdHits
Definition: CbmL1.h:308
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition: CbmMCDataManager.cxx:137
CbmL1::fSTAPDataMode
int fSTAPDataMode
Definition: CbmL1.h:249
CbmL1MCPoint::mass
double mass
Definition: CbmL1MCPoint.h:59
CbmKFTube::z
Double_t z
Definition: CbmKFMaterial.h:92
fscal
float fscal
Definition: L1/vectors/P4_F32vec4.h:250
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:136
CbmMCEventList::GetFileIdByIndex
Int_t GetFileIdByIndex(UInt_t index)
File number by index @value File number for event at given index in list.
Definition: CbmMCEventList.cxx:106
CbmMvdDetector.h
CbmL1::SetParContainers
void SetParContainers()
Definition: CbmL1.cxx:309
CbmL1TrackPar::NDF
int NDF
Definition: CbmL1TrackPar.h:18
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
CbmMCDataManager.h
CbmL1::fMomentumCutOff
Double_t fMomentumCutOff
Definition: CbmL1.h:253
CbmL1Track::StsHits
vector< int > StsHits
Definition: CbmL1Track.h:67
CbmL1::InputPerformance
void InputPerformance()
Definition: CbmL1Performance.cxx:2267
CbmMuchGeoScheme
Definition: CbmMuchGeoScheme.h:43
L1Algo.h
CbmMuchStation::GetLayer
CbmMuchLayer * GetLayer(Int_t iLayer) const
Definition: CbmMuchStation.h:58
_fvecalignment
#define _fvecalignment
Definition: L1/vectors/P4_F32vec4.h:254
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmStsParSetSensor.h
CbmL1::fInstance
static CbmL1 * fInstance
Definition: CbmL1.h:354
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmL1::fEventList
CbmMCEventList * fEventList
Definition: CbmL1.h:263
CbmStsStation::GetYmax
Double_t GetYmax() const
Definition: CbmStsStation.h:119
CbmKF.h
CbmMvdDetector::GetParameterFile
CbmMvdStationPar * GetParameterFile()
Definition: CbmMvdDetector.h:111
CbmL1MCPoint::x
double x
Definition: CbmL1MCPoint.h:56
CbmL1::vStsHits
vector< CbmL1StsHit > vStsHits
Used data = Repacked input data.
Definition: CbmL1.h:335
CbmL1MCTrack::q
double q
Definition: CbmL1MCTrack.h:31
CbmL1::fDigiMatchesMuch
TClonesArray * fDigiMatchesMuch
Definition: CbmL1.h:296
CbmL1::algo
L1Algo * algo
Definition: CbmL1.h:119
L1AlgoInputData::GetStsHitsStopIndex
const THitI * GetStsHitsStopIndex() const
Definition: L1AlgoInputData.h:55
CbmL1MCPoint
Definition: CbmL1MCPoint.h:23
CbmL1MCPoint::time
double time
Definition: CbmL1MCPoint.h:59
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmKF::GetMagneticField
FairField * GetMagneticField()
Definition: CbmKF.h:88
CbmL1HitStore::iStation
int iStation
Definition: CbmL1.h:97
CbmStsSetup.h
L1Station
Definition: L1Station.h:9
L1Field.h
CbmL1::WriteSTAPGeoData
void WriteSTAPGeoData(const vector< float > &geo)
STandAlone Package service-functions.
Definition: CbmL1.cxx:1711
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
L1Algo::GetVtxFieldValue
const L1FieldValue & GetVtxFieldValue() const
Definition: L1Algo.h:453
L1AlgoInputData
Definition: L1AlgoInputData.h:15
CbmL1TrackPar::is_electron
bool is_electron
Definition: CbmL1TrackPar.h:20
CbmTofDigiPar::GetNrOfModules
Int_t GetNrOfModules()
Definition: CbmTofDigiPar.h:44
CbmL1::fTrdPoints
CbmMCDataArray * fTrdPoints
Definition: CbmL1.h:307
CbmL1::NMuchStations
int NMuchStations
Definition: CbmL1.h:244
L1FieldValue::z
fvec z
Definition: L1Field.h:17
CbmL1::listStsDigi
vector< CbmStsDigi > listStsDigi
MC event list (all)
Definition: CbmL1.h:265
CbmL1MCPoint::y
double y
Definition: CbmL1MCPoint.h:56
CbmL1::fTofMatBudgetFileName
TString fTofMatBudgetFileName
Definition: CbmL1.h:363
CbmStsSetup::SetSensorConditions
UInt_t SetSensorConditions(CbmStsParSetSensorCond *conds)
Set sensor conditions from parameter container.
Definition: CbmStsSetup.cxx:470
L1Algo::StsHitsStartIndex
const THitI * StsHitsStartIndex
Definition: L1Algo.h:358
CbmL1MCTrack::nStations
int nStations
Definition: CbmL1MCTrack.h:131
CbmL1::listStsHitMatch
TClonesArray * listStsHitMatch
Definition: CbmL1.h:277
CbmL1::fData
L1AlgoInputData * fData
Definition: CbmL1.h:350
CbmMvdStationPar::GetXRes
Double_t GetXRes(Int_t stationNumber) const
Definition: CbmMvdStationPar.cxx:130
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmL1HitStore::x
double x
Definition: CbmL1.h:98
L1AlgoInputData::GetStsHits
const vector< L1StsHit > & GetStsHits() const
Definition: L1AlgoInputData.h:48
L1Event.h
CbmL1::fTimeSlice
CbmTimeSlice * fTimeSlice
Input data.
Definition: CbmL1.h:262
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmL1Track::Tpv
double Tpv[7]
Definition: CbmL1Track.h:63
CbmL1::listMvdHits
TClonesArray * listMvdHits
Definition: CbmL1.h:288
CbmL1MCPoint::event
int event
Definition: CbmL1MCPoint.h:64
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition: CbmStsSetup.cxx:293
L1Algo::vStsStrips
const vector< L1Strip > * vStsStrips
Definition: L1Algo.h:340
CbmStsStation.h
CbmStsSetup::IsSensorParsInit
Bool_t IsSensorParsInit() const
Initialisation status for sensor parameters.
Definition: CbmStsSetup.h:142
CbmL1MCTrack::y
double y
Definition: CbmL1MCTrack.h:31
L1AlgoInputData::GetStsStrips
const vector< L1Strip > & GetStsStrips() const
Definition: L1AlgoInputData.h:49
CbmL1MCTrack::x
double x
Definition: CbmL1MCTrack.h:31
CbmL1MCTrack::maxNStaMC
int maxNStaMC
Definition: CbmL1MCTrack.h:127
L1Algo::KFTrackFitter_simple
void KFTrackFitter_simple()
Track fitting procedures.
Definition: L1TrackFitter.cxx:39
L1Algo::StsHitsStopIndex
const THitI * StsHitsStopIndex
Definition: L1Algo.h:359
CbmL1::listMvdHitMatches
TClonesArray * listMvdHitMatches
Definition: CbmL1.h:290
CbmStsSetup::SetSensorParameters
UInt_t SetSensorParameters(CbmStsParSetSensor *parSet)
Set sensor parameters from parameter container.
Definition: CbmStsSetup.cxx:486
CbmL1Track::CLast
double CLast[21]
Definition: CbmL1Track.h:66
CbmL1HitStore
Definition: CbmL1.h:94
CbmL1.h
CbmMuchGeoScheme::Init
void Init(TObjArray *stations, Int_t flag)
Definition: CbmMuchGeoScheme.cxx:121
L1StsHit
Definition: L1StsHit.h:12
CbmKF::vMvdMaterial
std::vector< CbmKFTube > vMvdMaterial
Definition: CbmKF.h:73
CbmL1MCTrack::px
double px
Definition: CbmL1MCTrack.h:31
CbmStsFindTracks::MvdUsage
Bool_t MvdUsage() const
Definition: CbmStsFindTracks.h:65
CbmTofDigiPar.h
CbmTrdParModDigi::GetSizeX
Double_t GetSizeX() const
Definition: CbmTrdParModDigi.h:109
CbmL1::listStsPts
TClonesArray * listStsPts
Definition: CbmL1.h:272
TZPosI
unsigned short int TZPosI
Definition: L1StsHit.h:9
CbmL1::fDigiPar
CbmTofDigiPar * fDigiPar
Definition: CbmL1.h:315
CbmKFTube::R
Double_t R
Definition: CbmKFMaterial.h:93
CbmL1::fDigiFile
TString fDigiFile
Definition: CbmL1.h:121
L1Algo::vStsStripsB
const vector< L1Strip > * vStsStripsB
Definition: L1Algo.h:341
CbmL1::dFEI2vMCTracks
DFEI2I dFEI2vMCTracks
Definition: CbmL1.h:345
CbmL1::fMvdMatBudgetFileName
TString fMvdMatBudgetFileName
Definition: CbmL1.h:360
CbmTofDigiPar::GetCell
CbmTofCell * GetCell(Int_t i)
Definition: CbmTofDigiPar.h:48
CbmL1::fTrackingLevel
Int_t fTrackingLevel
Definition: CbmL1.h:252
CbmL1::fMvdPoints
CbmMCDataArray * fMvdPoints
Definition: CbmL1.h:268
CbmStsSetup::GetNofStations
Int_t GetNofStations() const
Definition: CbmStsSetup.h:90
CbmL1::fTofHitDigiMatches
TClonesArray * fTofHitDigiMatches
Definition: CbmL1.h:313
CbmStsParSetModule.h
CbmL1MCTrack::pdg
int pdg
Definition: CbmL1MCTrack.h:32
CbmL1Track::TLast
double TLast[7]
Definition: CbmL1Track.h:66
CbmL1::NMvdStations
int NMvdStations
Definition: CbmL1.h:244
CbmL1::listMvdPts
TClonesArray * listMvdPts
Definition: CbmL1.h:287
L1StsHit::b
TStripI b
Definition: L1StsHit.h:15
CbmMvdHit
Definition: CbmMvdHit.h:29
CbmL1MCTrack::nMCContStations
int nMCContStations
Definition: CbmL1MCTrack.h:125
CbmStsParSetSensorCond.h
CbmStsStation::GetRadLength
Double_t GetRadLength() const
Definition: CbmStsStation.h:60
CbmL1TrackPar::T
double T[7]
Definition: CbmL1TrackPar.h:17
CbmL1MCTrack::nHitContStations
int nHitContStations
Definition: CbmL1MCTrack.h:126
CbmL1::fGhostSuppression
Bool_t fGhostSuppression
Definition: CbmL1.h:254
CbmL1MCPoint::pz
double pz
Definition: CbmL1MCPoint.h:56
L1Strip
Definition: L1Strip.h:12
CbmL1::NTrdStations
int NTrdStations
Definition: CbmL1.h:244
CbmStsStation::GetSensorStereoAngle
Double_t GetSensorStereoAngle(Int_t iSide) const
Definition: CbmStsStation.cxx:166
L1AlgoInputData::ReadHitsFromFile
bool ReadHitsFromFile(const char work_dir[100], const int maxNEvent, const int iVerbose)
Definition: L1AlgoInputData.cxx:67
CbmL1MCTrack::py
double py
Definition: CbmL1MCTrack.h:31
L1Algo::NTracksIsecAll
unsigned int NTracksIsecAll
Definition: L1Algo.h:364
CbmStsStation::GetSensorD
Double_t GetSensorD() const
Definition: CbmStsStation.h:72
CbmL1::NStsStations
int NStsStations
Definition: CbmL1.h:244
THitI
unsigned int THitI
Definition: L1StsHit.h:8
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition: CbmStsHit.h:31
L1FieldValue::y
fvec y
Definition: L1Field.h:17
CbmStsParSetSensorCond
Parameters container for CbmStsParSensorCond.
Definition: CbmStsParSetSensorCond.h:30
L1Algo::L1KFTrackFitterMuch
void L1KFTrackFitterMuch()
Definition: L1TrackFitter.cxx:775
L1Algo::vSFlagB
const vector< unsigned char > * vSFlagB
Definition: L1Algo.h:351
L1Algo::NStations
int NStations
Definition: L1Algo.h:333
CbmKF::Instance
static CbmKF * Instance()
Definition: CbmKF.h:39
CbmL1MCTrack::maxNSensorMC
int maxNSensorMC
Definition: CbmL1MCTrack.h:128
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition: CbmMCEventList.h:90
CbmMuchGeoScheme::GetNStations
Int_t GetNStations() const
Definition: CbmMuchGeoScheme.h:79
CbmL1::fStsParSetSensor
CbmStsParSetSensor * fStsParSetSensor
Definition: CbmL1.h:370
CbmL1StsHit::hitId
int hitId
Definition: CbmL1StsHit.h:29
CbmL1::TrackMatch
void TrackMatch()
Reconstruction Performance.
Definition: CbmL1Performance.cxx:60
CbmL1::WriteSIMDKFData
void WriteSIMDKFData()
SIMD KF Banchmark service-functions.
Definition: CbmL1.cxx:2314
CbmL1MCTrack::CalculateIsReconstructable
void CalculateIsReconstructable()
Definition: CbmL1MCTrack.cxx:233
CbmTrdParSet::GetModuleId
virtual Int_t GetModuleId(Int_t i) const
Definition: CbmTrdParSet.cxx:32
CbmL1::fStsMatBudgetFileName
TString fStsMatBudgetFileName
Definition: CbmL1.h:359
CbmL1::fMuchPoints
CbmMCDataArray * fMuchPoints
Definition: CbmL1.h:294
CbmStsParSetSensor
Parameters container for CbmStsParSensor.
Definition: CbmStsParSetSensor.h:30
CbmL1::dFEI2vMCPoints
DFEI2I dFEI2vMCPoints
Definition: CbmL1.h:344
CbmL1MCPoint::yIn
double yIn
Definition: CbmL1MCPoint.h:57
CbmL1MCTrack::p
double p
Definition: CbmL1MCTrack.h:31
CbmL1::fClustersMuch
TClonesArray * fClustersMuch
Definition: CbmL1.h:297
CbmL1::fTimesliceMode
int fTimesliceMode
Definition: CbmL1.h:365
CbmStsSetup::IsSensorCondInit
Bool_t IsSensorCondInit() const
Initialisation status for sensor conditions.
Definition: CbmStsSetup.h:136
L1Track::NHits
unsigned char NHits
Definition: L1Track.h:22
CbmKFMaterial::RadLength
Double_t RadLength
Definition: CbmKFMaterial.h:31
CbmTofCell::GetX
Double_t GetX() const
Definition: CbmTofCell.h:36
CbmStsFindTracks
Definition: CbmStsFindTracks.h:34
CbmL1::ReadSTAPPerfData
void ReadSTAPPerfData()
Definition: CbmL1.cxx:2138
CbmMvdStationPar::GetYRes
Double_t GetYRes(Int_t stationNumber) const
Definition: CbmMvdStationPar.cxx:142
CbmL1MCTrack::Points
vector< int > Points
Definition: CbmL1MCTrack.h:33
CbmL1MCTrack::mother_ID
int mother_ID
Definition: CbmL1MCTrack.h:32
CbmTofCell
Definition: CbmTofCell.h:8
CbmL1MCTrack::z
double z
Definition: CbmL1MCTrack.h:31
CbmMuchLayer::GetDz
Double_t GetDz()
Definition: CbmMuchLayer.cxx:77
CbmTrdParModDigi
Definition of chamber gain conversion for one TRD module.
Definition: CbmTrdParModDigi.h:14
CbmL1::listMvdDigiMatches
TClonesArray * listMvdDigiMatches
Definition: CbmL1.h:289
CbmMuchGeoScheme::Instance
static CbmMuchGeoScheme * Instance()
Definition: CbmMuchGeoScheme.cxx:113
CbmL1MCTrack::IsPrimary
bool IsPrimary() const
Definition: CbmL1MCTrack.h:104
CbmL1::HistoPerformance
void HistoPerformance()
Definition: CbmL1Performance.cxx:612
CbmL1MCTrack::ID
int ID
Definition: CbmL1MCTrack.h:32
CbmL1MCPoint::iStation
int iStation
Definition: CbmL1MCPoint.h:61
CbmL1MCPoint::xOut
double xOut
Definition: CbmL1MCPoint.h:58
CbmMvdStationPar.h
CbmTrdParModDigi.h
CbmL1MCPoint::p
double p
Definition: CbmL1MCPoint.h:59
CbmL1PFFitter.h
CbmL1::fUseTRD
Bool_t fUseTRD
Definition: CbmL1.h:255
CbmL1::ReadEvent
void ReadEvent(L1AlgoInputData *, CbmEvent *event=NULL)
Read information about hits, mcPoints and mcTracks into L1 classes.
Definition: CbmL1ReadEvent.cxx:92
L1Algo::vSFlag
const vector< unsigned char > * vSFlag
Definition: L1Algo.h:350
CbmL1::Exec
void Exec(Option_t *option)
Definition: CbmL1.cxx:1308
L1Algo::vTracks
L1Vector< L1Track > vTracks
Definition: L1Algo.h:355
CbmStsStation::GetSensorPitch
Double_t GetSensorPitch(Int_t iSide) const
Definition: CbmStsStation.cxx:151
CbmL1::fUseMVD
Bool_t fUseMVD
Definition: CbmL1.h:255
CbmTrdParSet::GetModulePar
virtual const CbmTrdParMod * GetModulePar(Int_t detId) const
Definition: CbmTrdParSet.cxx:45
CbmL1MCTrack::iEvent
int iEvent
Definition: CbmL1MCTrack.h:32
CbmL1StsHit::event
int event
Definition: CbmL1StsHit.h:39
CbmL1MCPoint::pyIn
double pyIn
Definition: CbmL1MCPoint.h:57
L1Algo
Definition: L1Algo.h:79
CbmL1MCPoint::hitIds
vector< int > hitIds
Definition: CbmL1MCPoint.h:73
CbmKFTube::dz
Double_t dz
Definition: CbmKFMaterial.h:92
CbmMCEventList
Container class for MC events with number, file and start time.
Definition: CbmMCEventList.h:38
CbmL1::ReadSTAPAlgoData
void ReadSTAPAlgoData()
Definition: CbmL1.cxx:2013
CbmL1::fStsPoints
CbmMCDataArray * fStsPoints
Definition: CbmL1.h:266
CbmStsSetup::IsModuleParsInit
Bool_t IsModuleParsInit() const
Initialisation status for module parameters.
Definition: CbmStsSetup.h:130
L1StsHit::t_reco
float t_reco
Definition: L1StsHit.h:17
CbmTrdParModDigi::GetZ
Double_t GetZ() const
Definition: CbmTrdParModDigi.h:114
CbmL1::fSTAPDataDir
TString fSTAPDataDir
Definition: CbmL1.h:250
CbmL1::IdealTrackFinder
void IdealTrackFinder()
--— Ideal Tracking --------------------------—
Definition: CbmL1.cxx:1649
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmStsSetup::Init
Bool_t Init(const char *geometryFile=nullptr)
Initialise the setup.
Definition: CbmStsSetup.cxx:201
CbmL1::fUseMUCH
Bool_t fUseMUCH
Definition: CbmL1.h:255
CbmTofCell.h
first
bool first
Definition: LKFMinuit.cxx:143
CbmL1::fDigisMuch
TClonesArray * fDigisMuch
Definition: CbmL1.h:300
CbmKFTube::r
Double_t r
Definition: CbmKFMaterial.h:93
CbmL1::Finish
void Finish()
Definition: CbmL1.cxx:1616
CbmStsSetup
Class representing the top level of the STS setup.
Definition: CbmStsSetup.h:39
CbmL1MCTrack::maxNStaHits
int maxNStaHits
Definition: CbmL1MCTrack.h:129
CbmL1StsHit::extIndex
int extIndex
Definition: CbmL1StsHit.h:33
CbmL1Track
Definition: CbmL1Track.h:33
ClassImp
ClassImp(CbmL1) static L1Algo algo_static _fvecalignment
L1Algo::Init
void Init(const vector< fscal > &geo, const bool UseHitErrors, const bool mCBMmode)
Definition: L1Algo.cxx:5
CbmL1::~CbmL1
~CbmL1()
Definition: CbmL1.cxx:304
CbmL1MCPoint::zIn
double zIn
Definition: CbmL1MCPoint.h:57
CbmTrdParSetDigi.h
CbmL1::fUseTOF
Bool_t fUseTOF
Definition: CbmL1.h:255
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
CbmL1::NStation
int NStation
Definition: CbmL1.h:244
CbmMvdDetector
Definition: CbmMvdDetector.h:39
CbmL1HitStore::ExtIndex
int ExtIndex
Definition: CbmL1.h:96
CbmTrdParModDigi::GetSizeZ
Double_t GetSizeZ() const
Definition: CbmTrdParModDigi.h:111
CbmStsSetup::GetStation
CbmStsStation * GetStation(Int_t stationId) const
Definition: CbmStsSetup.h:97
CbmMuchPad.h
CbmL1MCPoint::pxOut
double pxOut
Definition: CbmL1MCPoint.h:58
CbmL1::TrackFitPerformance
void TrackFitPerformance()
Definition: CbmL1Performance.cxx:1543
L1Track
Definition: L1Track.h:20
CbmL1MCPoint::z
double z
Definition: CbmL1MCPoint.h:56
CbmL1::ReInit
virtual InitStatus ReInit()
Definition: CbmL1.cxx:328
CbmL1::EfficienciesPerformance
void EfficienciesPerformance()
Definition: CbmL1Performance.cxx:293
L1AlgoInputData::GetSFlag
const L1Vector< unsigned char > & GetSFlag() const
Definition: L1AlgoInputData.h:52
CbmL1MCPoint::ID
int ID
Definition: CbmL1MCPoint.h:60
CbmL1::fTrdDigiPar
CbmTrdParSetDigi * fTrdDigiPar
Definition: CbmL1.h:304
L1Algo::vRecoHits
L1Vector< THitI > vRecoHits
Definition: L1Algo.h:356
CbmL1::NTOFStation
int NTOFStation
Definition: CbmL1.h:245
CbmTimeSlice
Bookkeeping of time-slice content.
Definition: CbmTimeSlice.h:29
L1Algo::NMvdStations
int NMvdStations
Definition: L1Algo.h:334
CbmL1StsHit
Definition: CbmL1StsHit.h:8
CbmMCEventList::GetEventIdByIndex
Int_t GetEventIdByIndex(UInt_t index)
Event number by index @value Event number for event at given index in list.
Definition: CbmMCEventList.cxx:77
L1AlgoInputData.h
CbmKFTube
Definition: CbmKFMaterial.h:77
CbmL1::writedir2current
static void writedir2current(TObject *obj)
Definition: CbmL1.cxx:1632
CbmL1MCPoint::mother_ID
int mother_ID
Definition: CbmL1MCPoint.h:60
CbmL1::Init
virtual InitStatus Init()
Definition: CbmL1.cxx:333
CbmL1::CbmL1
CbmL1()
Definition: CbmL1.cxx:90
L1Algo::L1KFTrackFitter
void L1KFTrackFitter()
Definition: L1TrackFitter.cxx:335
CbmL1MCTrack::time
double time
Definition: CbmL1MCTrack.h:31
CbmL1::vRTracks
vector< CbmL1Track > vRTracks
Definition: CbmL1.h:125
CbmL1MCPoint::pyOut
double pyOut
Definition: CbmL1MCPoint.h:58
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofDigiPar
Definition: CbmTofDigiPar.h:18
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
L1StsHit.h
CbmStsSetup::SetModuleParameters
UInt_t SetModuleParameters(CbmStsParSetModule *modulePars)
Set module parameters from parameter container.
Definition: CbmStsSetup.cxx:455
L1AlgoInputData::GetStsStripsB
const vector< L1Strip > & GetStsStripsB() const
Definition: L1AlgoInputData.h:50
CbmL1HitStore::y
double y
Definition: CbmL1.h:98
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofDigiPar::GetCellId
Int_t GetCellId(Int_t i)
Definition: CbmTofDigiPar.h:45
CbmL1MCTrack::pz
double pz
Definition: CbmL1MCTrack.h:31
CbmMvdStationPar
Definition: CbmMvdStationPar.h:22
L1Track::TFirst
fscal TFirst[7]
Definition: L1Track.h:25
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:137
CbmL1StsHit::mcPointIds
vector< int > mcPointIds
Definition: CbmL1StsHit.h:35
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition: CbmEvent.h:30
CbmL1::vHitMCRef
vector< int > vHitMCRef
Definition: CbmL1.h:339
CbmStsParSetModule
Parameters container for CbmStsParModule.
Definition: CbmStsParSetModule.h:30
CbmMuchGeoScheme::GetStation
CbmMuchStation * GetStation(Int_t iStation) const
Definition: CbmMuchGeoScheme.cxx:209
L1AlgoInputData::GetStsZPos
const vector< fscal > & GetStsZPos() const
Definition: L1AlgoInputData.h:51
CbmTrdParSetDigi
Definition: CbmTrdParSetDigi.h:15
CbmMuchStation::GetNLayers
Int_t GetNLayers() const
Definition: CbmMuchStation.h:50
L1AlgoInputData::Clear
void Clear()
Definition: L1AlgoInputData.h:79
CbmMuchLayer::GetZ
Double_t GetZ() const
Definition: CbmMuchLayer.h:56
CbmL1::fStsParSetSensorCond
CbmStsParSetSensorCond * fStsParSetSensorCond
Definition: CbmL1.h:371
CbmL1MCTrack::StsHits
vector< int > StsHits
Definition: CbmL1MCTrack.h:34
L1Branch.h
CbmKFTrackInterface::SetId
void SetId(int id)
Definition: CbmKFTrackInterface.h:68
CbmL1MCPoint::pzOut
double pzOut
Definition: CbmL1MCPoint.h:58
CbmL1::fMuchPixelHits
TClonesArray * fMuchPixelHits
Definition: CbmL1.h:299
CbmL1MCPoint::px
double px
Definition: CbmL1MCPoint.h:56
CbmL1MCPoint::pdg
int pdg
Definition: CbmL1MCPoint.h:60
CbmTofCell::GetY
Double_t GetY() const
Definition: CbmTofCell.h:37
CbmL1MCPoint::yOut
double yOut
Definition: CbmL1MCPoint.h:58
CbmL1MCTrack
Definition: CbmL1MCTrack.h:29
CbmL1TrackPar::C
double C[21]
Definition: CbmL1TrackPar.h:17
L1Algo::MaxNStations
@ MaxNStations
Definition: L1Algo.h:331
CbmL1::vMCPoints
vector< CbmL1MCPoint > vMCPoints
Definition: CbmL1.h:197
CbmStsStation::GetXmax
Double_t GetXmax() const
Definition: CbmStsStation.h:117
CbmL1::fmCBMmode
bool fmCBMmode
Definition: CbmL1.h:123
CbmMuchGeoScheme.h
CbmL1TrackPar::mass
double mass
Definition: CbmL1TrackPar.h:19
CbmL1Track::fTrackTime
double fTrackTime
Definition: CbmL1Track.h:71
CbmL1MCPoint::pxIn
double pxIn
Definition: CbmL1MCPoint.h:57
CbmL1::listStsClusters
TClonesArray * listStsClusters
Definition: CbmL1.h:275
CbmL1::fTrdHitMatches
TClonesArray * fTrdHitMatches
Definition: CbmL1.h:309
CbmMuchLayer
Definition: CbmMuchLayer.h:21
L1AlgoInputData::GetStsHitsStartIndex
const THitI * GetStsHitsStartIndex() const
Definition: L1AlgoInputData.h:54
CbmL1::vFileEvent
DFSET vFileEvent
Definition: CbmL1.h:127
CbmL1MCTrack::mass
double mass
Definition: CbmL1MCTrack.h:31
CbmL1
Definition: CbmL1.h:113
CbmL1::HitMatch
void HitMatch()
Input Performance.
Definition: CbmL1ReadEvent.cxx:1504
L1Algo::fRadThick
vector< L1Material > fRadThick
Definition: L1Algo.h:337
L1Track::Momentum
float Momentum
Definition: L1Track.h:24
CbmL1::listStsClusterMatch
TClonesArray * listStsClusterMatch
Definition: CbmL1.h:278
CbmL1MCPoint::xIn
double xIn
Definition: CbmL1MCPoint.h:57
CbmL1::ReadSTAPGeoData
void ReadSTAPGeoData(vector< float > &geo, int &size)
Definition: CbmL1.cxx:1997
CbmL1::eatwhite
static std::istream & eatwhite(std::istream &is)
Definition: CbmL1.cxx:1983
CbmMuchStation.h
CbmL1MCPoint::py
double py
Definition: CbmL1MCPoint.h:56
CbmL1::vHitStore
vector< CbmL1HitStore > vHitStore
Definition: CbmL1.h:180
L1FieldValue
Definition: L1Field.h:11
CbmL1::vMCTracks
vector< CbmL1MCTrack > vMCTracks
Definition: CbmL1.h:337
L1StsHit::f
TStripI f
Definition: L1StsHit.h:15
CbmL1::listStsHits
TClonesArray * listStsHits
Definition: CbmL1.h:276
L1StsHit::iz
TZPosI iz
Definition: L1StsHit.h:24
CbmStsFindTracks.h
CbmKFVertex.h
CbmL1::WriteSTAPPerfData
void WriteSTAPPerfData()
Definition: CbmL1.cxx:1833
CbmL1::Reconstruct
void Reconstruct(CbmEvent *event=NULL)
Definition: CbmL1.cxx:1310
CbmL1::fMCTracks
CbmMCDataArray * fMCTracks
Definition: CbmL1.h:267
CbmL1::WriteSTAPAlgoData
void WriteSTAPAlgoData()
Definition: CbmL1.cxx:1726
CbmL1::fPerformance
Int_t fPerformance
Definition: CbmL1.h:247
CbmL1::fMuchMatBudgetFileName
TString fMuchMatBudgetFileName
Definition: CbmL1.h:361
CbmL1MCPoint::zOut
double zOut
Definition: CbmL1MCPoint.h:58
CbmStsStation
Class representing a station of the StsSystem.
Definition: CbmStsStation.h:29
L1AlgoInputData::GetSFlagB
const L1Vector< unsigned char > & GetSFlagB() const
Definition: L1AlgoInputData.h:53
CbmStsSetup::IsInit
Bool_t IsInit() const
Initialisation status for sensor parameters.
Definition: CbmStsSetup.h:124
CbmStsStation::GetSensorRotation
Double_t GetSensorRotation() const
Definition: CbmStsStation.h:99
CbmL1::fGlobalMode
bool fGlobalMode
Definition: CbmL1.h:124
CbmMuchModuleGem.h
L1Algo::CATrackFinder
void CATrackFinder()
The main procedure - find tracks.
Definition: L1CATrackFinder.cxx:1786
L1Algo::SetData
void SetData(const vector< L1StsHit > &StsHits_, const vector< L1Strip > &StsStrips_, const vector< L1Strip > &StsStripsB_, const vector< fscal > &StsZPos_, const vector< unsigned char > &SFlag_, const vector< unsigned char > &SFlagB_, const THitI *StsHitsStartIndex_, const THitI *StsHitsStopIndex_)
Definition: L1Algo.cxx:179
CbmL1MCTrack::IsReconstructable
bool IsReconstructable() const
Definition: CbmL1MCTrack.h:105
CbmL1::fTofPoints
CbmMCDataArray * fTofPoints
Definition: CbmL1.h:312
CbmL1Track::Cpv
double Cpv[21]
Definition: CbmL1Track.h:63
L1FieldValue::x
fvec x
Definition: L1Field.h:15
CbmL1MCPoint::pzIn
double pzIn
Definition: CbmL1MCPoint.h:57
CbmL1::listMuchHitMatches
TClonesArray * listMuchHitMatches
Definition: CbmL1.h:295
CbmL1TrackPar::chi2
double chi2
Definition: CbmL1TrackPar.h:17
CbmL1::histodir
TDirectory * histodir
Definition: CbmL1.h:352
CbmL1::fTrdMatBudgetFileName
TString fTrdMatBudgetFileName
Definition: CbmL1.h:362