CbmRichRingTrackAssignClosestD.cxx 7.29 KB
Newer Older
1
2
3
4
5
6
7
8
9
/**
* \file CbmRichRingTrackAssignClosestD.cxx
*
* \author Claudia Hoehne and Semen Lebedev
* \date 2007
**/

#include "CbmRichRingTrackAssignClosestD.h"

10
#include "CbmEvent.h"
Administrator's avatar
Administrator committed
11
#include "CbmGlobalTrack.h"
12
#include "CbmMCTrack.h"
13
#include "CbmRichRing.h"
Administrator's avatar
Administrator committed
14
#include "CbmTrdTrack.h"
15

16
17
#include "FairRootManager.h"
#include "FairTrackParam.h"
18
#include <Logger.h>
19
20
21
22

#include "TClonesArray.h"

#include <algorithm>
Administrator's avatar
Administrator committed
23
#include <iostream>
24
25
26
27
28
29
#include <vector>

using std::cout;
using std::endl;
using std::vector;

30
CbmRichRingTrackAssignClosestD::CbmRichRingTrackAssignClosestD() {}
Administrator's avatar
Administrator committed
31
32
33

CbmRichRingTrackAssignClosestD::~CbmRichRingTrackAssignClosestD() {}

34
35
36
37
void CbmRichRingTrackAssignClosestD::Init()
{
  FairRootManager* manager = FairRootManager::Instance();
  if (nullptr == manager) LOG(fatal) << "CbmRichRingTrackAssignClosestD::Init(): FairRootManager is nullptr.";
Administrator's avatar
Administrator committed
38

39
40
  fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
  if (fGlobalTracks == nullptr) LOG(fatal) << "CbmRichRingTrackAssignClosestD::Init(): No GlobalTrack.";
Administrator's avatar
Administrator committed
41

42
  fTrdTracks = (TClonesArray*) manager->GetObject("TrdTrack");
Administrator's avatar
Administrator committed
43
  //if (NULL == fTrdTracks) {Fatal("CbmRichRingTrackAssignClosestD::Init", "No TrdTrack array!");}
44
45
}

46
47
48
49
void CbmRichRingTrackAssignClosestD::DoAssign(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
{
  fEventNum++;
  if (fAlgorithmType == RingTrack) { DoAssignRingTrack(event, rings, richProj); }  // RingTrack algorithm
50

Administrator's avatar
Administrator committed
51
  else if (fAlgorithmType == TrackRing) {
52
    DoAssignTrackRing(event, rings, richProj);
Administrator's avatar
Administrator committed
53
  }  // TrackRing
54

Administrator's avatar
Administrator committed
55
  else if (fAlgorithmType == Combined) {
56
57
    DoAssignRingTrack(event, rings, richProj);
    DoAssignTrackRing(event, rings, richProj);
Administrator's avatar
Administrator committed
58
  }  // Combined
59
60
61

  Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
  Int_t nofRings  = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
62
63
64
  LOG(debug) << "CbmRichRingTrackAssignClosestD::DoAssign(): Event:" << fEventNum << " rings:" << nofRings
             << " ringsInTS:" << rings->GetEntriesFast() << " tracks:" << nofTracks
             << " tracksInTS:" << richProj->GetEntriesFast();
65
66
}

67
68
69
70
71
72
73
void CbmRichRingTrackAssignClosestD::DoAssignRingTrack(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
{

  const Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
  const Int_t nofRings  = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
  if (nofTracks <= 0 || nofRings <= 0) return;

Administrator's avatar
Administrator committed
74
75
76
77
78
79
80
81
82
83
  vector<Int_t> trackIndex;
  vector<Double_t> trackDist;
  trackIndex.resize(nofRings);
  trackDist.resize(nofRings);
  for (UInt_t i = 0; i < trackIndex.size(); i++) {
    trackIndex[i] = -1;
    trackDist[i]  = 999.;
  }

  for (Int_t iIter = 0; iIter < 4; iIter++) {
84
85
86
87
88
89
90
91
92
93
    for (Int_t iR0 = 0; iR0 < nofRings; iR0++) {
      Int_t iR = event ? event->GetIndex(ECbmDataType::kRichRing, iR0) : iR0;

      if (trackIndex[iR] != -1) continue;
      CbmRichRing* ring = static_cast<CbmRichRing*>(rings->At(iR));
      if (ring == nullptr) continue;
      if (ring->GetNofHits() < fMinNofHitsInRing) continue;

      Double_t xRing  = ring->GetCenterX();
      Double_t yRing  = ring->GetCenterY();
Administrator's avatar
Administrator committed
94
95
96
      Double_t rMin   = 999.;
      Int_t iTrackMin = -1;

97
98
99
100
      for (Int_t iT0 = 0; iT0 < nofTracks; iT0++) {
        Int_t iT = event ? event->GetIndex(ECbmDataType::kRichTrackProjection, iT0) : iT0;

        vector<Int_t>::iterator it = find(trackIndex.begin(), trackIndex.end(), iT);
Administrator's avatar
Administrator committed
101
102
        if (it != trackIndex.end()) continue;

103
104
105
        FairTrackParam* track = static_cast<FairTrackParam*>(richProj->At(iT));
        Double_t xTrack       = track->GetX();
        Double_t yTrack       = track->GetY();
Administrator's avatar
Administrator committed
106
107
108
        // no projection onto the photodetector plane
        if (xTrack == 0 && yTrack == 0) continue;

109
        if (fUseTrd && fTrdTracks != nullptr && !IsTrdElectron(iT)) continue;
Administrator's avatar
Administrator committed
110

111
        Double_t dist = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack) + (yRing - yTrack) * (yRing - yTrack));
Administrator's avatar
Administrator committed
112
113
114

        if (dist < rMin) {
          rMin      = dist;
115
          iTrackMin = iT;
Administrator's avatar
Administrator committed
116
117
        }
      }  // loop tracks
118
119
      trackIndex[iR] = iTrackMin;
      trackDist[iR]  = rMin;
Administrator's avatar
Administrator committed
120
121
122
123
124
125
126
127
128
    }  //loop rings

    for (UInt_t i1 = 0; i1 < trackIndex.size(); i1++) {
      for (UInt_t i2 = 0; i2 < trackIndex.size(); i2++) {
        if (i1 == i2) continue;
        if (trackIndex[i1] == trackIndex[i2] && trackIndex[i1] != -1) {
          if (trackDist[i1] >= trackDist[i2]) {
            trackDist[i1]  = 999.;
            trackIndex[i1] = -1;
129
130
          }
          else {
Administrator's avatar
Administrator committed
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
            trackDist[i2]  = 999.;
            trackIndex[i2] = -1;
          }
        }
      }
    }
  }  //iIter

  // fill global tracks
  for (UInt_t i = 0; i < trackIndex.size(); i++) {
    //		CbmRichRing* pRing = (CbmRichRing*)rings->At(i);
    // cout << "trackIndex[i]:" << trackIndex[i] << " trackDist[i]:" << trackDist[i] << " r:" << pRing->GetRadius() << " x:" << pRing->GetCenterX() << " y:" << pRing->GetCenterY()<< endl;
    if (trackIndex[i] == -1) continue;
    CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(trackIndex[i]);
    gTrack->SetRichRingIndex(i);
  }
147
148
}

149
150
151
152
153
154
155
156
void CbmRichRingTrackAssignClosestD::DoAssignTrackRing(CbmEvent* event, TClonesArray* rings, TClonesArray* richProj)
{
  const Int_t nofTracks = event ? event->GetNofData(ECbmDataType::kRichTrackProjection) : richProj->GetEntriesFast();
  const Int_t nofRings  = event ? event->GetNofData(ECbmDataType::kRichRing) : rings->GetEntriesFast();
  if (nofTracks <= 0 || nofRings <= 0) return;
  for (Int_t iT0 = 0; iT0 < nofTracks; iT0++) {
    Int_t iT               = event ? event->GetIndex(ECbmDataType::kRichTrackProjection, iT0) : iT0;
    CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iT));
Administrator's avatar
Administrator committed
157
    // track already has rich segment
158
    if (gTrack == nullptr || gTrack->GetRichRingIndex() >= 0) continue;
Administrator's avatar
Administrator committed
159

160
161
162
    FairTrackParam* track = static_cast<FairTrackParam*>(richProj->At(iT));
    Double_t xTrack       = track->GetX();
    Double_t yTrack       = track->GetY();
Administrator's avatar
Administrator committed
163
    if (xTrack == 0 && yTrack == 0) continue;
164
    if (fUseTrd && fTrdTracks != nullptr && !IsTrdElectron(iT)) continue;
Administrator's avatar
Administrator committed
165
166
    Double_t rMin  = 999.;
    Int_t iRingMin = -1;
167
168
169
170
171
172
173
174
175
    for (Int_t iR0 = 0; iR0 < nofRings; iR0++) {
      Int_t iR          = event ? event->GetIndex(ECbmDataType::kRichRing, iR0) : iR0;
      CbmRichRing* ring = static_cast<CbmRichRing*>(rings->At(iR));
      if (ring == nullptr) continue;
      if (ring->GetNofHits() < fMinNofHitsInRing) continue;

      Double_t xRing = ring->GetCenterX();
      Double_t yRing = ring->GetCenterY();
      Double_t dist  = TMath::Sqrt((xRing - xTrack) * (xRing - xTrack) + (yRing - yTrack) * (yRing - yTrack));
Administrator's avatar
Administrator committed
176
177
      if (dist < rMin) {
        rMin     = dist;
178
        iRingMin = iR;
Administrator's avatar
Administrator committed
179
180
181
182
183
184
      }
    }  // loop rings
    if (iRingMin < 0) continue;
    //		CbmRichRing* pRing = (CbmRichRing*)rings->At(iRingMin);
    gTrack->SetRichRingIndex(iRingMin);
  }  //loop tracks
185
186
}

187
188
189
Bool_t CbmRichRingTrackAssignClosestD::IsTrdElectron(Int_t iTrack)
{
  CbmGlobalTrack* gTrack = static_cast<CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
Administrator's avatar
Administrator committed
190
191
  Int_t trdIndex         = gTrack->GetTrdTrackIndex();
  if (trdIndex == -1) return false;
192
  CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdIndex));
Administrator's avatar
Administrator committed
193
  if (NULL == trdTrack) return false;
194

Administrator's avatar
Administrator committed
195
  if (trdTrack->GetPidANN() > fTrdAnnCut) { return true; }
196

Administrator's avatar
Administrator committed
197
  return false;
198
}