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

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov, I. Rostovtseva (2016)
 *  
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  Read Event information to L1
 *
 *====================================================================
 */

22
#include "CbmEvent.h"
23
#include "CbmKF.h"
Administrator's avatar
Administrator committed
24
#include "CbmL1.h"
25
#include "CbmMCDataObject.h"
26
27
#include "CbmMatch.h"
#include "CbmMuchGeoScheme.h"
Administrator's avatar
Administrator committed
28
#include "CbmMuchPixelHit.h"
29
#include "CbmMuchPoint.h"
30
#include "CbmStsAddress.h"
31
32
33
34
#include "CbmStsCluster.h"
#include "CbmStsDigi.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
35
36
#include "CbmStsSetup.h"
#include "CbmTofDigiBdfPar.h"
Administrator's avatar
Administrator committed
37
#include "CbmTofHit.h"
38
#include "CbmTofPoint.h"
Administrator's avatar
Administrator committed
39
40
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
41

42
43
#include "FairMCEventHeader.h"

44
45
46
#include "L1Algo/L1Algo.h"
#include "L1AlgoInputData.h"

47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
//#include "CbmMvdHitMatch.h"


#include "TDatabasePDG.h"
#include "TRandom.h"

#include <iostream>

using std::cout;
using std::endl;
using std::find;

//#define MVDIDEALHITS
//#define STSIDEALHITS


63
64
static bool compareZ(const int& a, const int& b)
{
Administrator's avatar
Administrator committed
65
  //        return (CbmL1::Instance()->vMCPoints[a].z < CbmL1::Instance()->vMCPoints[b].z);
66
67
68
69
  const CbmL1* l1 = CbmL1::Instance();
  return l1->Get_Z_vMCPoint(a) < l1->Get_Z_vMCPoint(b);
}

70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
/// Local structure for a hit, containing both measured and MC information.
/// The structure is used to sort hits before writing them into L1 input arrays
///
struct TmpHit {
  int iStripF;             ///> index of front strip
  int iStripB;             ///> index of back strip
  int iStation;            ///> index of station
  int ExtIndex;            ///> index of hit in the external TClonesArray array (NOTE: negative for MVD)
  double u_front, u_back;  ///> positions of strips
  double x, y, z;          ///> position of hit (Cortesian coordinates)
  double xmc, ymc, p;      ///> mc position of hit, total momentum
  double tx, ty;           ///> mc slopes of the mc point
  double dx, dy, dxy;      ///> hit position errors in Cortesian coordinates
  double du, dv;           ///> hit position errors in strips coordinates
  int iMC;                 ///> index of MCPoint in the vMCPoints array
Sergei Zharko's avatar
Sergei Zharko committed
85
86
87
  double time = 0.;        ///> time of the hit
  double dt   = 1.e10;     ///> time error of the hit
  int Det;
88
  int id;
89
  int track;
90

Sergei Zharko's avatar
Sergei Zharko committed
91
  /// Provides comparison of two hits.
92
93
94
95
96
97
  /// If two hits belong to different stations,
  /// the smallest hit belongs to the station with the smallest index. Otherwise, the smallest hit
  /// has the smallest y coordinate
  /// \param  a  Left hit
  /// \param  b  Right hit
  /// \return boolean: true - the left hit is smaller then the right one
98
99
100
  static bool Compare(const TmpHit& a, const TmpHit& b)
  {
    return (a.iStation < b.iStation) || ((a.iStation == b.iStation) && (a.y < b.y));
101
  }
102

103
104
105
106
107
108
  /// Creates a hit from the CbmL1MCPoint object
  /// \param point  constant reference to the input MC-point
  /// \param det
  /// \param nTmpHits
  /// \param nStripF
  /// \param ip
Sergei Zharko's avatar
Sergei Zharko committed
109
  /// \param NStrips
110
111
  /// \param st       reference to the station info object
  // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
112
  void CreateHitFromPoint(const CbmL1MCPoint& point, int det, int nTmpHits, int nStripF, int ip, int& NStrips,
113
                          const L1Station& st)
114
115
116
117
118
119
  {
    ExtIndex = 0;
    Det      = det;
    id       = nTmpHits;  //tmpHits.size();
    iStation = point.iStation;

Valentina Akishina's avatar
Valentina Akishina committed
120
121
    dt   = st.dt[0];
    time = point.time + gRandom->Gaus(0, dt);
122
123
124
125
126
127
128
129
130

    iStripF = nStripF + ip;  //firstDetStrip + ip;
    iStripB = iStripF;
    if (NStrips <= iStripF) { NStrips = iStripF + 1; }

    float f_sigma = 0.1;  // sqrt(st.frontInfo.sigma2[0]);
    //1.;
    float b_sigma = 0.1;  // sqrt(st.backInfo.sigma2[0]);//1.;

Valentina Akishina's avatar
Valentina Akishina committed
131
132
133
134
135
    dx  = f_sigma;
    dy  = b_sigma;
    dxy = 0;
    du  = f_sigma;
    dv  = b_sigma;
136

Valentina Akishina's avatar
Valentina Akishina committed
137
138
    x = point.x;  // + gRandom->Gaus(0, dx);
    y = point.y;  // + gRandom->Gaus(0, dy);
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    z = point.z;

    xmc = point.x;
    ymc = point.y;

    track = point.ID;

    p = point.p;

    u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
    u_back  = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
    iMC     = ip;
  }

Sergei Zharko's avatar
Sergei Zharko committed
153
  /// Sets randomized position and time of the hit
154
155
156
157
158
  /// The positions are smeared within predefined errors dx, dy, dt; z coordinate
  /// of the hit is known precisely
  /// \param point  constant reference to the input MC-point
  /// \param st     reference to the station info object
  // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
159
  void SetHitFromPoint(const CbmL1MCPoint& point, const L1Station& st)
160
161
162
163
164
165
166
167
  {
    x       = 0.5 * (point.xIn + point.xOut) + gRandom->Gaus(0, dx);
    y       = 0.5 * (point.yIn + point.yOut) + gRandom->Gaus(0, dy);
    z       = 0.5 * (point.zIn + point.zOut);
    time    = point.time + gRandom->Gaus(0, dt);
    u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
    u_back  = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
  }
168
169
};

Administrator's avatar
Administrator committed
170
void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, float& /*TsOverlap*/, int& FstHitinTs,
171
                      bool& areDataLeft, CbmEvent* event)
172
{
173
174
  if (fVerbose >= 10) cout << "ReadEvent: start." << endl;

175
176
  areDataLeft = false;  // no data left after reading the sub-timeslice

Sergey Gorbunov's avatar
Sergey Gorbunov committed
177
178
  fData_->Clear();

Administrator's avatar
Administrator committed
179
  // clear arrays for next event
Sergei Zharko's avatar
Sergei Zharko committed
180
181
182
183
184
185
186
187
188
  vMCPoints.clear();               /* <CbmL1MCPoint> */
  vMCPoints_in_Time_Slice.clear(); /* <int>          */
  vMCTracks.clear();               /* <CbmL1MCTrack> */
  vStsHits.clear();                /* <CbmL1Hit>     */
  vRTracks.clear();                /* <CbmL1Track>   */
  vHitMCRef.clear();               /* <int>: indexes of MC-points in vMCPoints (by index of algo->vStsHits) */
  vHitStore.clear();               /* <CbmL1HitStore> */
  dFEI2vMCPoints.clear();          /* dFEI vs MC-point index: dFEI = index * 10000 + fileID + eventID * 0.0001 */
  dFEI2vMCTracks.clear();          /* dFEI vs MC-track index: dFEI = index * 10000 + fileID + eventID * 0.0001 */
Sergey Gorbunov's avatar
Sergey Gorbunov committed
189

190
  if (fVerbose >= 10) cout << "ReadEvent: clear is done." << endl;
Administrator's avatar
Administrator committed
191

Sergey Gorbunov's avatar
Sergey Gorbunov committed
192
193
194
195
196
197
198
  L1Vector<TmpHit> tmpHits("CbmL1ReadEvent::tmpHits");

  {  // reserve enough space for hits
    int nHitsTotal = 0;
    if (listMvdHits) nHitsTotal += listMvdHits->GetEntriesFast();
    Int_t nEntSts = 0;
    if (listStsHits) {
199
      if (event) { nEntSts = event->GetNofData(ECbmDataType::kStsHit); }
Sergey Gorbunov's avatar
Sergey Gorbunov committed
200
201
202
203
204
205
206
207
208
209
210
      else {
        nEntSts = listStsHits->GetEntriesFast();
      }
      if (nEntSts < 0) nEntSts = 0;  // GetNofData() can return -1;
    }
    nHitsTotal += nEntSts;
    if (fMuchPixelHits) nHitsTotal += fMuchPixelHits->GetEntriesFast();
    if (listTrdHits) nHitsTotal += listTrdHits->GetEntriesFast();
    if (fTofHits) nHitsTotal += fTofHits->GetEntriesFast();
    tmpHits.reserve(nHitsTotal);
  }
211

Administrator's avatar
Administrator committed
212
  // -- produce Sts hits from space points --
213

Administrator's avatar
Administrator committed
214
  for (int i = 0; i < NStation; i++) {
215

216
    fData_->StsHitsStartIndex[i] = static_cast<L1HitIndex_t>(-1);
217
218
219
    fData_->StsHitsStopIndex[i]  = 0;
  }

220
221
222
223
224
225
226

  nMvdPoints      = 0;
  int nStsPoints  = 0;
  int nTrdPoints  = 0;
  int nMuchPoints = 0;
  int nTofPoints  = 0;

Administrator's avatar
Administrator committed
227
228
  // get MVD hits
  Int_t nMvdHits  = 0;
229
  Int_t nMuchHits = 0;
Administrator's avatar
Administrator committed
230
231
232
  Int_t nTrdHits  = 0;
  Int_t nTofHits  = 0;
  // get STS hits
233
  int nStsHits = 0;
Administrator's avatar
Administrator committed
234

235
236
237
238
  int firstTrdPoint  = 0;
  int firstStsPoint  = 0;
  int firstMuchPoint = 0;
  int firstTofPoint  = 0;
Administrator's avatar
Administrator committed
239

240
241
242
243
244
245
246
  /*
   * MC hits and tracks gathering: START
   *
   * If performance is studied, i.e. fPerformanse > 0, MC-information is used for
   * efficiencies calculation. Otherwise this 
   */

Administrator's avatar
Administrator committed
247
  if (fPerformance) {
248
    // Fill vMCTracks and dFEI2vMCTracks
249
    Fill_vMCTracks();
250
251

    // Fill vMCPoints, vMCPoints_in_Time_Slice and dFEI2vMCPoints
Sergey Gorbunov's avatar
Sergey Gorbunov committed
252
253
254
255
    vMCPoints.clear();
    vMCPoints.reserve(5 * vMCTracks.size() * NStation);
    vMCPoints_in_Time_Slice.clear();
    vMCPoints_in_Time_Slice.reserve(vMCPoints.capacity());
Administrator's avatar
Administrator committed
256

257
    for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) {
Administrator's avatar
Administrator committed
258
      Int_t iFile  = set_it->first;
259
260
      Int_t iEvent = set_it->second;

261
      // TODO: Probably to put this code in L1Algo interfaces (S.Zharko)
262
      if (fMvdPoints) {
Administrator's avatar
Administrator committed
263
        Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent);
264
        double maxDeviation     = 0;
Administrator's avatar
Administrator committed
265
        for (Int_t iMC = 0; iMC < nMvdPointsInEvent; iMC++) {
266
          CbmL1MCPoint MC;
Administrator's avatar
Administrator committed
267
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) {
268
            MC.iStation     = -1;
Sergei Zharko's avatar
Sergei Zharko committed
269
            const L1Station* sta = algo->GetStations().begin();
270
271
            double bestDist = 1.e20;
            for (Int_t iSt = 0; iSt < NMvdStations; iSt++) {
272
273
274
              // use z_in since z_out is sometimes very wrong
              // due to a problem in transport
              double d = (MC.zIn - sta[iSt].z[0]);
275
276
277
278
              if (fabs(d) < fabs(bestDist)) {
                bestDist    = d;
                MC.iStation = iSt;
              }
279
            }
280
            assert(MC.iStation >= 0);
281
            if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
Administrator's avatar
Administrator committed
282
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
283
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
284
            assert(trk_it != dFEI2vMCTracks.end());
285
            MC.ID = trk_it->second;
286
287
            vMCTracks[MC.ID].Points.push_back_no_warning(vMCPoints.size());
            dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC), vMCPoints.size()));
288
289
290
291
292
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
            nMvdPoints++;
          }
        }
Sergey Gorbunov's avatar
Sergey Gorbunov committed
293
        if (fVerbose > 2) { LOG(info) << "CbmL1ReadEvent: max deviation of Mvd points " << maxDeviation; }
294
        // ensure that the nominal station positions are not far from the sensors
295
        // assert(fabs(maxDeviation) < 1.);
296
297
      }

298
      firstStsPoint = vMCPoints.size();
299
      if (fStsPoints) {
300
301
        Int_t nMC           = fStsPoints->Size(iFile, iEvent);
        double maxDeviation = 0;
302
303
304
        for (Int_t iMC = 0; iMC < nMC; iMC++) {
          CbmL1MCPoint MC;
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) {
305
            MC.iStation     = -1;
Sergei Zharko's avatar
Sergei Zharko committed
306
            const L1Station* sta = algo->GetStations().begin() + NMvdStations;
307
            double bestDist = 1.e20;
308
            for (Int_t iSt = 0; iSt < NStsStations; iSt++) {
309
310
311
              // use z_in since z_out is sometimes very wrong
              // due to a problem in transport
              double d = (MC.zIn - sta[iSt].z[0]);
312
313
314
315
              if (fabs(d) < fabs(bestDist)) {
                bestDist    = d;
                MC.iStation = NMvdStations + iSt;
              }
316
317
            }
            assert(MC.iStation >= 0);
318
            if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
319
320
321
322
323
324
325
326
327
328
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            assert(trk_it != dFEI2vMCTracks.end());
            MC.ID = trk_it->second;
            vMCTracks[MC.ID].Points.push_back_no_warning(vMCPoints.size());
            dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints), vMCPoints.size()));
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
            nStsPoints++;
          }
Administrator's avatar
Administrator committed
329
        }
Sergey Gorbunov's avatar
Sergey Gorbunov committed
330
        if (fVerbose > 2) { LOG(info) << "CbmL1ReadEvent: max deviation of Sts points " << maxDeviation; }
331
        // ensure that the nominal station positions are not far from the sensors
332
        //  assert(fabs(maxDeviation) < 1.);
Administrator's avatar
Administrator committed
333
334
      }

335
      firstMuchPoint = vMCPoints.size();
Administrator's avatar
Administrator committed
336
      if (fMuchPoints) {
337
338
        Int_t nMC = fMuchPoints->Size(iFile, iEvent);
        for (Int_t iMC = 0; iMC < nMC; iMC++) {
339
          CbmL1MCPoint MC;
Administrator's avatar
Administrator committed
340
341
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) {
            MC.iStation    = -1;
342
            const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations;
343
344
345
346
            for (Int_t iSt = 0; iSt < NMuchStations; iSt++) {
              if (MC.z > sta[iSt].z[0] - 2.5) { MC.iStation = NMvdStations + NStsStations + iSt; }
            }
            assert(MC.iStation >= 0);
Administrator's avatar
Administrator committed
347
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
348
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
349
            assert(trk_it != dFEI2vMCTracks.end());
350
            MC.ID = trk_it->second;
351
352
353
            vMCTracks[MC.ID].Points.push_back_no_warning(vMCPoints.size());
            dFEI2vMCPoints.insert(
              DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints), vMCPoints.size()));
354
355
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
Administrator's avatar
Administrator committed
356
            nMuchPoints++;
357
          }
Administrator's avatar
Administrator committed
358
359
360
        }
      }

361
      firstTrdPoint = vMCPoints.size();
Administrator's avatar
Administrator committed
362
363
      if (fTrdPoints)
        for (Int_t iMC = 0; iMC < fTrdPoints->Size(iFile, iEvent); iMC++) {
364
          CbmL1MCPoint MC;
Administrator's avatar
Administrator committed
365
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) {
366
            MC.iStation    = -1;
367
            const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations;
368
369
370
            for (Int_t iSt = 0; iSt < NTrdStations; iSt++) {
              if (MC.z > sta[iSt].z[0] - 4.0) { MC.iStation = NMvdStations + NStsStations + NMuchStations + iSt; }
            }
371
            if (MC.iStation < 0) continue;
372
            assert(MC.iStation >= 0);
Administrator's avatar
Administrator committed
373
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
374
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
375
            assert(trk_it != dFEI2vMCTracks.end());
376
            MC.ID = trk_it->second;
377
378
379
            vMCTracks[MC.ID].Points.push_back_no_warning(vMCPoints.size());
            dFEI2vMCPoints.insert(
              DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints + nMuchPoints), vMCPoints.size()));
380
381
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
Administrator's avatar
Administrator committed
382
            nTrdPoints++;
383
          }
Administrator's avatar
Administrator committed
384
        }
385

Administrator's avatar
Administrator committed
386

387
      firstTofPoint = vMCPoints.size();
Administrator's avatar
Administrator committed
388
      if (fTofPoints) {
389

390
        vector<float> TofPointToTrackdZ[NTOFStation];
391

392
        TofPointToTrack.resize(NTOFStation);
Administrator's avatar
Administrator committed
393

394
        for (Int_t i = 0; i < NTOFStation; i++) {
395

396
397
398
          TofPointToTrack[i].resize(vMCTracks.size(), 0);
          TofPointToTrackdZ[i].resize(vMCTracks.size());
        }
399

400
        for (int iSt = 0; iSt < NTOFStation; iSt++)
Valentina Akishina's avatar
Valentina Akishina committed
401
          for (unsigned int i = 0; i < TofPointToTrackdZ[iSt].size(); i++) {
402
            TofPointToTrackdZ[iSt][i] = 100000;
Valentina Akishina's avatar
Valentina Akishina committed
403
404
            TofPointToTrack[iSt][i]   = -1;
          }
Administrator's avatar
Administrator committed
405

Valentina Akishina's avatar
Valentina Akishina committed
406

407
        for (Int_t iMC = 0; iMC < fTofPoints->Size(iFile, iEvent); iMC++) {
Administrator's avatar
Administrator committed
408

409
          CbmL1MCPoint MC;
Administrator's avatar
Administrator committed
410

411

412
413
414
415
416
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 4)) {
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            if (trk_it == dFEI2vMCTracks.end()) continue;
            Int_t IND_Track = trk_it->second;
Administrator's avatar
Administrator committed
417

418
            MC.iStation    = -1;
Sergei Zharko's avatar
Sergei Zharko committed
419
420
            const L1Station* sta =
              algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations + NTrdStations;
421
            for (Int_t iSt = 0; iSt < NTOFStation; iSt++)
Valentina Akishina's avatar
Valentina Akishina committed
422
              MC.iStation = (MC.z > sta[iSt].z[0] - 15)
423
424
                              ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
                              : MC.iStation;
425
426
427
            if (MC.iStation < 0) continue;
            assert(MC.iStation >= 0);
            int iTofSta = MC.iStation - (NMvdStations + NStsStations + NMuchStations + NTrdStations);
Valentina Akishina's avatar
Valentina Akishina committed
428
429
430
431
432
433

            if (iTofSta >= 0) {
              float dz = TofPointToTrackdZ[iTofSta][IND_Track];
              if (fabs(sta[MC.iStation].z[0] - MC.z) < dz) { TofPointToTrack[iTofSta][IND_Track] = iMC; }
              TofPointToTrackdZ[iTofSta][IND_Track] = fabs(sta[MC.iStation].z[0] - MC.z);
            }
434
435
          }
        }
Administrator's avatar
Administrator committed
436

437
438
        for (int iTofSta = 0; iTofSta < NTOFStation; iTofSta++)
          for (unsigned int iMC = 0; iMC < TofPointToTrack[iTofSta].size(); iMC++) {
Administrator's avatar
Administrator committed
439

Valentina Akishina's avatar
Valentina Akishina committed
440
            if (TofPointToTrack[iTofSta][iMC] == -1) continue;
Administrator's avatar
Administrator committed
441

442
            CbmL1MCPoint MC;
443

444
            if (!ReadMCPoint(&MC, TofPointToTrack[iTofSta][iMC], iFile, iEvent, 4)) {
Administrator's avatar
Administrator committed
445

446
              MC.iStation    = -1;
Sergei Zharko's avatar
Sergei Zharko committed
447
448
              const L1Station* sta =
                algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations + NTrdStations;
449
              for (Int_t iSt = 0; iSt < NTOFStation; iSt++)
Valentina Akishina's avatar
Valentina Akishina committed
450
                MC.iStation = (MC.z > sta[iSt].z[0] - 15)
451
452
453
                                ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
                                : MC.iStation;

Valentina Akishina's avatar
Valentina Akishina committed
454
455
              if (MC.iStation < 0) continue;
              TofPointToTrack[iTofSta][iMC] = vMCPoints.size();
456
              vMCTracks[iMC].Points.push_back_no_warning(vMCPoints.size());
Administrator's avatar
Administrator committed
457

458
459
460
461
462
463
464
465
466
467
              MC.ID = iMC;

              vMCPoints.push_back(MC);
              vMCPoints_in_Time_Slice.push_back(0);

              dFEI2vMCPoints.insert(DFEI2I::value_type(
                dFEI(iFile, iEvent, TofPointToTrack[iTofSta][iMC] + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints),
                vMCPoints.size() - 1));
              nTofPoints++;
            }
468
469
          }
      }
Administrator's avatar
Administrator committed
470
    }  //time_slice
471

Administrator's avatar
Administrator committed
472
    for (unsigned int iTr = 0; iTr < vMCTracks.size(); iTr++) {
473

474
      sort(vMCTracks[iTr].Points.begin(), vMCTracks[iTr].Points.end(), compareZ);
Administrator's avatar
Administrator committed
475
476

      if (vMCTracks[iTr].mother_ID >= 0) {
477
        Double_t dtrck          = dFEI(vMCTracks[iTr].iFile, vMCTracks[iTr].iEvent, vMCTracks[iTr].mother_ID);
Administrator's avatar
Administrator committed
478
        DFEI2I::iterator map_it = dFEI2vMCTracks.find(dtrck);
479
        if (map_it == dFEI2vMCTracks.end()) vMCTracks[iTr].mother_ID = -2;
Administrator's avatar
Administrator committed
480
481
482
483
        else
          vMCTracks[iTr].mother_ID = map_it->second;
      }
    }  //iTr
484
485
    if (fVerbose >= 10) cout << "Points in vMCTracks are sorted." << endl;

Administrator's avatar
Administrator committed
486
487
  }  //fPerformance

488
489
490
491
492
493
494
495
496
497
498
499
  /*
   * MC hits and tracks gathering: END
   */

  /*
   * Measured hits gathering: START
   *
   * In this section the measured hits from different detector subsystems are reformatted according to TmpHit structure
   * (NOTE: independent from particular detector design) and then pushed to the tmpHit vector. In the performance study mode
   * matching with MC points is done
   */

500
501
  int NStrips = 0;

502
  //
Administrator's avatar
Administrator committed
503
  // get MVD hits
504
  //
Administrator's avatar
Administrator committed
505
  if (listMvdHits) {
506

507
508
    int firstDetStrip = NStrips;

509
    for (int j = 0; j < listMvdHits->GetEntriesFast(); j++) {
510
511
      TmpHit th;
      {
Administrator's avatar
Administrator committed
512
513
        CbmMvdHit* mh = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(j));
        th.ExtIndex   = -(1 + j);
514
        th.id         = tmpHits.size();
Administrator's avatar
Administrator committed
515
        th.iStation   = mh->GetStationNr();
516
        th.iStripF    = firstDetStrip + j;
Sergey Gorbunov's avatar
Sergey Gorbunov committed
517
        th.iStripB    = th.iStripF;
518
        if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
Administrator's avatar
Administrator committed
519

520
521
522
        TVector3 pos, err;
        mh->Position(pos);
        mh->PositionError(err);
Administrator's avatar
Administrator committed
523

524
525
        th.dx = mh->GetDx();
        th.dy = mh->GetDy();
Administrator's avatar
Administrator committed
526
527
528

        th.du  = mh->GetDx();
        th.dv  = mh->GetDy();
529
530
531
532
        th.dxy = mh->GetDxy();

        th.x = pos.X();
        th.y = pos.Y();
533
        th.z = pos.Z();
534

535
        const L1Station& st = algo->GetStations()[th.iStation];
536
537
        th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
538
      }
539
540
      th.Det = 0;
      th.iMC = -1;
Administrator's avatar
Administrator committed
541
542
      if (fPerformance) {
        if (listMvdHitMatches) {
543
          CbmMatch* mvdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(listMvdHitMatches->At(j));
Administrator's avatar
Administrator committed
544
545
          if (mvdHitMatch->GetNofLinks() > 0)
            if (mvdHitMatch->GetLink(0).GetIndex() < nMvdPoints) {
546
              th.iMC = mvdHitMatch->GetLink(0).GetIndex();
547
548
549
#ifdef MVDIDEALHITS
//TODO
#endif
Administrator's avatar
Administrator committed
550
551
            }
        }
552
      }
Administrator's avatar
Administrator committed
553
      //if(  h.MC_Point >=0 ) // DEBUG !!!!
554
555
556
557
      {
        tmpHits.push_back(th);
        nMvdHits++;
      }
Administrator's avatar
Administrator committed
558
559
    }  // for j
  }    // if listMvdHits
560
561
  if (fVerbose >= 10) cout << "ReadEvent: mvd hits are gotten." << endl;

562
563
564
565
566
  if ((2 == fStsUseMcHit)) {  // SG!! create TRD hits from TRD points

    int firstDetStrip = NStrips;
    for (int ip = firstStsPoint; ip < firstStsPoint + nStsPoints; ip++) {
      const CbmL1MCPoint& p = vMCPoints[ip];
Valentina Akishina's avatar
Valentina Akishina committed
567
568
569
      //       int mcTrack           = p.ID;
      //       if (mcTrack < 0) continue;
      //       const CbmL1MCTrack& t = vMCTracks[mcTrack];
570
571
572
573
574
575
      //if (t.p < 1) continue;
      // if (t.Points.size() > 4) continue;
      // cout << "sts mc: station " << p.iStation - NMvdStations << " x " << p.x << " y " << p.y << " z " << p.z << " t "
      //            << p.time << " mc " << p.ID << " p " << p.p << endl;
      TmpHit th;
      int DetId = 1;
576
      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
577
578
579
580
581
      tmpHits.push_back(th);
      nStsHits++;
    }
  }

582
583
584
  //
  // Get STS hits
  //
585
  if (listStsHits && (2 != fStsUseMcHit)) {
586

Sergey Gorbunov's avatar
Sergey Gorbunov committed
587
    Int_t nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : listStsHits->GetEntriesFast());
Valentina Akishina's avatar
Valentina Akishina committed
588

589
    int firstDetStrip = NStrips;
Valentina Akishina's avatar
Valentina Akishina committed
590
591

    if (event) FstHitinTs = 0;
592

593
    for (Int_t j = FstHitinTs; j < nEntSts; j++) {
594
      Int_t hitIndex = 0;
Valentina Akishina's avatar
Valentina Akishina committed
595
      hitIndex       = (event ? event->GetIndex(ECbmDataType::kStsHit, j) : j);
Administrator's avatar
Administrator committed
596

597
      int hitIndexSort = 0;
598
      if (!fLegacyEventMode) hitIndexSort = StsIndex[hitIndex];
599
600
601
602
      else
        hitIndexSort = j;

      CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndexSort));
603
604
605

      TmpHit th;
      {
606
607
        CbmStsHit* mh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndexSort));
        th.ExtIndex   = hitIndexSort;
Administrator's avatar
Administrator committed
608
        th.Det        = 1;
609
610
        th.iStation =
          NMvdStations + CbmStsSetup::Instance()->GetStationNumber(mh->GetAddress());  //mh->GetStationNr() - 1;
611
612
613
614
615
        th.iStripF = firstDetStrip + mh->GetFrontClusterId();
        th.iStripB = firstDetStrip + mh->GetBackClusterId();

        if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
        if (NStrips <= th.iStripB) { NStrips = th.iStripB + 1; }
616
617

        //Get time
Administrator's avatar
Administrator committed
618
        th.time = mh->GetTime();
Sergey Gorbunov's avatar
Sergey Gorbunov committed
619
        th.dt   = mh->GetTimeError();
Administrator's avatar
Administrator committed
620

621
622
        if (!fLegacyEventMode) { th.id = nMvdHits + hitIndex; }
        else {
623
          th.id = tmpHits.size();
624
        }
625

626
        /// stop if reco TS ends and many hits left
627
        if (!event)
Valentina Akishina's avatar
Valentina Akishina committed
628
629
630
631
          if ((th.time > (TsStart + TsLength)) && ((nEntSts - hitIndex) > 300)) {
            areDataLeft = true;  // there are unprocessed data left in the time slice
            break;
          }
632

633
634
635
636
637
638
        TVector3 pos, err;
        mh->Position(pos);
        mh->PositionError(err);

        th.x = pos.X();
        th.y = pos.Y();
639
        th.z = pos.Z();
Administrator's avatar
Administrator committed
640
641
642

        th.dx  = mh->GetDx();
        th.dy  = mh->GetDy();
643
        th.dxy = mh->GetDxy();
Administrator's avatar
Administrator committed
644

645
646
647
        th.du = mh->GetDu();
        th.dv = mh->GetDv();

648
        const L1Station& st = algo->GetStations()[th.iStation];
649
650
        th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
651
652
653
654
655
      }
      th.iMC = -1;

      Int_t iMC = -1;

Administrator's avatar
Administrator committed
656
657
      if (fPerformance) {
        if (listStsClusterMatch) {
658
659
660
661
          const CbmMatch* frontClusterMatch =
            static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetFrontClusterId()));
          const CbmMatch* backClusterMatch =
            static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetBackClusterId()));
662
663
          CbmMatch stsHitMatch;

664
          for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) {
665
666
            const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);

667
            for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); iBackLink++) {
668
              const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
Administrator's avatar
Administrator committed
669
              if (frontLink == backLink) {
670
671
672
673
674
675
                stsHitMatch.AddLink(frontLink);
                stsHitMatch.AddLink(backLink);
              }
            }
          }

Administrator's avatar
Administrator committed
676
          if (stsHitMatch.GetNofLinks() > 0) {
677
            Float_t bestWeight = 0.f;
Administrator's avatar
Administrator committed
678
679
680
            for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
              if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
                bestWeight   = stsHitMatch.GetLink(iLink).GetWeight();
681
682
                Int_t iFile  = stsHitMatch.GetLink(iLink).GetFile();
                Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry();
683
                //                 if(fLegacyEventMode) //TODO Fix the event number in links
Administrator's avatar
Administrator committed
684
                //                   iEvent+=1;
685
                Int_t iIndex            = stsHitMatch.GetLink(iLink).GetIndex() + nMvdPoints;
Administrator's avatar
Administrator committed
686
                Double_t dtrck          = dFEI(iFile, iEvent, iIndex);
687
                DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
Administrator's avatar
Administrator committed
688
                if (trk_it == dFEI2vMCPoints.end()) continue;
689
690
691
692
                iMC = trk_it->second;
              }
            }
          }
693
694
695
696
697
#ifdef STSIDEALHITS
          // TODO
          // CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex));
          // h.z = 0.5 * (point->GetZOut() + point->GetZIn());
#endif
698
        }
699
        else {
Administrator's avatar
Administrator committed
700
          iMC = sh->GetRefId();  // TODO1: don't need this with FairLinks
701
702
        }
      }  //fPerformance
703

Administrator's avatar
Administrator committed
704
705
706
      if (iMC > -1) {
        th.iMC = iMC;
        // th.track = iMC;
707
708
709
710
711
      }

      tmpHits.push_back(th);
      nStsHits++;

Administrator's avatar
Administrator committed
712
713
    }  // for j
  }    // if listStsHits
714

715
716
717
  //
  // Get MuCh hits
  //
718
719
720
721
722
723
  if ((2 == fMuchUseMcHit) && fUseMUCH) {  // SG!! create TRD hits from TRD points
    int firstDetStrip = NStrips;

    for (int ip = firstMuchPoint; ip < firstMuchPoint + nMuchPoints; ip++) {
      const CbmL1MCPoint& p = vMCPoints[ip];

Valentina Akishina's avatar
Valentina Akishina committed
724
725
726
      //       int mcTrack = p.ID;
      //       if (mcTrack < 0) continue;
      //       const CbmL1MCTrack& t = vMCTracks[mcTrack];
727
728
729
730
731
      //if (t.p < 1) continue;
      // if (t.Points.size() > 4) continue;

      TmpHit th;
      int DetId = 2;
732
      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
733
734
735
736
737
738
739

      tmpHits.push_back(th);
      nMuchHits++;
    }
  }

  if (fUseMUCH && fMuchPixelHits && (2 != fMuchUseMcHit)) {
740

741
    Int_t nEnt = fMuchPixelHits->GetEntriesFast();
Administrator's avatar
Administrator committed
742

743
    int firstDetStrip = NStrips;
Administrator's avatar
Administrator committed
744
745

    for (int j = 0; j < nEnt; j++) {
746
747
      TmpHit th;
      {
Administrator's avatar
Administrator committed
748

749
        CbmMuchPixelHit* mh = static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(j));
Administrator's avatar
Administrator committed
750

751
        th.ExtIndex = j;
Administrator's avatar
Administrator committed
752
        th.Det      = 2;
753
        th.id       = tmpHits.size();
754

Administrator's avatar
Administrator committed
755

756
757
        Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(mh->GetAddress());
        Int_t layerNumber   = CbmMuchGeoScheme::GetLayerIndex(mh->GetAddress());
Administrator's avatar
Administrator committed
758
759
760
761
762
763

        int DetId = stationNumber * 3 + layerNumber;


        th.iStation = DetId + NMvdStations + NStsStations;
        //Get time
764
        th.time = mh->GetTime() - 14.5;
Sergey Gorbunov's avatar
Sergey Gorbunov committed
765
        th.dt   = mh->GetTimeError();
766
767


Administrator's avatar
Administrator committed
768
        //   th.iSector  = 0;
769
770
771
        th.iStripF = firstDetStrip + j;
        th.iStripB = th.iStripF;
        if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
772
773

        th.x = mh->GetX();
774
        th.y = mh->GetY();
775
        th.z = mh->GetZ();
Administrator's avatar
Administrator committed
776
777
778

        th.dx  = mh->GetDx();
        th.dy  = mh->GetDy();
779
        th.dxy = 0;
Administrator's avatar
Administrator committed
780

781
782
783
        th.du = mh->GetDx();
        th.dv = mh->GetDy();

Administrator's avatar
Administrator committed
784

785
        const L1Station& st = algo->GetStations()[th.iStation];
786
787
        th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
788
      }
Administrator's avatar
Administrator committed
789
      th.iMC  = -1;
790
791
792
      int iMC = -1;


Administrator's avatar
Administrator committed
793
794
      if (fPerformance) {
        if (listMuchHitMatches) {
795
          CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(listMuchHitMatches->At(j));
Administrator's avatar
Administrator committed
796
797
798
799
800
801
802
803
804
805
806
807
808
809


          for (Int_t iLink = 0; iLink < matchHitMatch->GetNofLinks(); iLink++) {
            if (matchHitMatch->GetLink(iLink).GetIndex() < nMuchPoints) {
              iMC          = matchHitMatch->GetLink(iLink).GetIndex();
              Int_t iIndex = iMC + nMvdPoints + nStsPoints;

              Int_t iFile  = matchHitMatch->GetLink(0).GetFile();
              Int_t iEvent = matchHitMatch->GetLink(0).GetEntry();

              Double_t dtrck          = dFEI(iFile, iEvent, iIndex);
              DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
              if (trk_it == dFEI2vMCPoints.end()) continue;
              th.iMC = trk_it->second;
810
              if ((1 == fMuchUseMcHit) && (th.iMC > -1))
811
                th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
Administrator's avatar
Administrator committed
812
            }
813
          }
Administrator's avatar
Administrator committed
814
        }
815
816
      }

Administrator's avatar
Administrator committed
817
818
819
820
821
      tmpHits.push_back(th);
      nMuchHits++;

    }  // for j
  }    // if listMvdHits
822

823
  if ((2 == fTrdUseMcHit) && fTrdHitMatches && fUseTRD) {  // SG!! create TRD hits from TRD points
824
825
826
827
828
    int firstDetStrip = NStrips;

    for (int ip = firstTrdPoint; ip < firstTrdPoint + nTrdPoints; ip++) {
      const CbmL1MCPoint& p = vMCPoints[ip];

Valentina Akishina's avatar
Valentina Akishina committed
829
830
831
      //       int mcTrack = p.ID;
      //       if (mcTrack < 0) continue;
      //       const CbmL1MCTrack& t = vMCTracks[mcTrack];
832
833

      TmpHit th;
834
      int DetId = 3;
835
      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
836
837
838
839
      tmpHits.push_back(th);
      nTrdHits++;
    }
  }
Sergei Zharko's avatar
Sergei Zharko committed
840

841
842
843
  //
  // Get TRD hits
  //
844
  if (fUseTRD && listTrdHits && (2 != fTrdUseMcHit)) {
845
    int firstDetStrip = NStrips;
846
    vector<bool> mcUsed(nTrdPoints, 0);
847

848
    for (int j = 0; j < listTrdHits->GetEntriesFast(); j++) {
849
      TmpHit th;
Administrator's avatar
Administrator committed
850
851
852
853

      CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(j));
      th.ExtIndex   = j;
      th.Det        = 3;
Valentina Akishina's avatar
Valentina Akishina committed
854

855
      th.id = tmpHits.size();
856

857
858
      int sta = mh->GetPlaneId();

859
      if ((fTrackingMode == L1Algo::TrackingMode::kMcbm) && (sta > 1) && (fMissingHits)) { sta = sta - 1; }
860
861
862

      th.iStation = NMvdStations + sta + NStsStations + NMuchStations;

863
864
865
866
      //  if (mh->GetPlaneId()==0) continue;
      //  if (mh->GetPlaneId()==1) continue;
      //  if (mh->GetPlaneId()==2) continue;
      //  if (mh->GetPlaneId()==3) continue;
Administrator's avatar
Administrator committed
867
868

      th.time = mh->GetTime();
Sergey Gorbunov's avatar
Sergey Gorbunov committed
869
      th.dt   = mh->GetTimeError();
Administrator's avatar
Administrator committed
870
871

      //   th.iSector  = 0;
872
873
874
      th.iStripF = firstDetStrip + j;
      th.iStripB = th.iStripF;  //TrdHitsOnStationIndex[num+1][k];
      if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
Administrator's avatar
Administrator committed
875
876
877
878
879
880
881

      TVector3 pos, err;
      mh->Position(pos);
      mh->PositionError(err);

      th.x = pos.X();
      th.y = pos.Y();
882
      th.z = pos.Z();
Administrator's avatar
Administrator committed
883
884
885
886
887
888
889
890
891


      th.dx  = fabs(mh->GetDx());
      th.dy  = fabs(mh->GetDy());
      th.dxy = 0;

      th.du = fabs(mh->GetDx());
      th.dv = fabs(mh->GetDy());

892
      const L1Station& st = algo->GetStations()[th.iStation];
893
894
      th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
      th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
Administrator's avatar
Administrator committed
895
896

      th.iMC  = -1;
897
898
      int iMC = -1;

899
      if (fPerformance && fTrdHitMatches) {
900
901
902

        CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j));

Valentina Akishina's avatar
Valentina Akishina committed
903
        if (trdHitMatch->GetNofLinks() > 0)
904
905
906
          if (trdHitMatch->GetLink(0).GetIndex() < nTrdPoints) {
            iMC    = trdHitMatch->GetLink(0).GetIndex();
            th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
Valentina Akishina's avatar
Valentina Akishina committed
907
            //                                th.track = vMCPoints[th.iMC].ID;
908

909
910
            if ((1 == fTrdUseMcHit) && (th.iMC > -1))  //SG!!! replace hits by MC points

911
              th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
912
913

            if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {  //SG!!! replace hits by MC points
Valentina Akishina's avatar
Valentina Akishina committed
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936