CbmRoot
CbmStsWkn.cxx
Go to the documentation of this file.
1 
7 #include "CbmStsWkn.h"
8 #include "FairRootManager.h"
9 
10 #include "CbmGlobalTrack.h"
11 
12 #include "CbmStsCluster.h"
13 #include "CbmStsDigi.h"
14 #include "CbmStsHit.h"
15 #include "CbmStsTrack.h"
16 #include "TClonesArray.h"
17 #include "TMath.h"
18 #include "TString.h"
19 #include "TSystem.h"
20 
21 #include <cmath>
22 #include <iostream>
23 using std::cout;
24 using std::endl;
25 using std::fabs;
26 using std::map;
27 using std::vector;
28 
30  : fdegWkn(4), fEmp(2.4), fXi(0.5), fk1(fdegWkn + 1), fnSet(8), fwkn(-1) {
31  Init();
32 }
33 
35 
37  FairRootManager* ioman = FairRootManager::Instance();
38  if (ioman != nullptr) {
39  fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
40  fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
41  fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
42  fStsClusterArray = (TClonesArray*) ioman->GetObject("StsCluster");
43  fStsDigiArray = (TClonesArray*) ioman->GetObject("StsDigi");
44  }
45 }
46 
47 Double_t CbmStsWkn::GetStsWkn(Int_t StsTrackIndex) {
48  double dr = 1.;
49  std::vector<float> eLossVector;
50  std::vector<float> dEdxAllveto;
51 
52  CbmStsTrack* StsTrack = (CbmStsTrack*) fStsTracks->At(StsTrackIndex);
53  if (StsTrack == NULL) return fwkn;
54 
55  int nClustersWveto =
56  StsTrack->GetNofStsHits()
57  + StsTrack->GetNofStsHits(); //assume all clusters with veto
58  if (nClustersWveto < 8) return fwkn;
59 
60  //cout<<fnSet<<endl;
61  for (int iHit = 0; iHit < StsTrack->GetNofStsHits(); ++iHit) {
62  Int_t StsHitIndex = StsTrack->GetStsHitIndex(iHit);
63  CbmStsHit* stsHit = (CbmStsHit*) fStsHits->At(StsHitIndex);
64 
65  // ***********dx=dr is calculated from the track inclination
66  double x, y, z, xNext, yNext, zNext;
67  x = stsHit->GetX();
68  y = stsHit->GetY();
69  z = stsHit->GetZ();
70 
71  if (iHit != StsTrack->GetNofStsHits() - 1) {
72  CbmStsHit* stsHitNext =
73  (CbmStsHit*) fStsHits->At(StsTrack->GetStsHitIndex(iHit + 1));
74  xNext = stsHitNext->GetX();
75  yNext = stsHitNext->GetY();
76  zNext = stsHitNext->GetZ();
77  dr = sqrt((xNext - x) * (xNext - x) + (yNext - y) * (yNext - y)
78  + (zNext - z) * (zNext - z))
79  / (zNext - z); // if *300um, you get real reconstructed dr
80  } // else dr stay previous
81  // ***********End of dr calculation
82 
83  //************dE is defined as a total cluster charge
84  Int_t FrontClusterId = stsHit->GetFrontClusterId();
85  CbmStsCluster* frontCluster =
86  (CbmStsCluster*) fStsClusterArray->At(FrontClusterId);
87  Int_t BackClusterId = stsHit->GetBackClusterId();
88  CbmStsCluster* backCluster =
89  (CbmStsCluster*) fStsClusterArray->At(BackClusterId);
90 
91  if (!frontCluster || !backCluster) return fwkn;
92 
93  dEdxAllveto.push_back((frontCluster->GetCharge()) / dr);
94  dEdxAllveto.push_back((backCluster->GetCharge()) / dr);
95  }
96 
97  if (dEdxAllveto.size() != 0) {
98  unsigned int NSample = dEdxAllveto.size();
99 
100  for (unsigned int jVec = 0; jVec < NSample; jVec++)
101  dEdxAllveto[jVec] = dEdxAllveto[jVec] / 10000;
102 
103  for (unsigned int jVec = 0; jVec < NSample; jVec++)
104  dEdxAllveto[jVec] = (dEdxAllveto[jVec] - fEmp) / fXi - 0.225;
105 
106  sort(dEdxAllveto.begin(), dEdxAllveto.end());
107 
108  for (unsigned int jVec = 0; jVec < NSample; jVec++)
109  dEdxAllveto[jVec] = TMath::LandauI(dEdxAllveto[jVec]);
110 
111 
112  for (int iHit = 0; iHit < fnSet; iHit++)
113  eLossVector.push_back(dEdxAllveto[NSample - fnSet + iHit]);
114  }
115 
116  Double_t S = 0, ty = 0, ti = 0;
117 
118  for (Int_t i = 0; i < fnSet; i++) {
119  ty = eLossVector[i];
120  ti = i;
121  S += pow((ti - 1) / fnSet - ty, fk1) - pow(ti / fnSet - ty, fk1);
122  }
123  float fwkn0 = pow(fnSet, 0.5 * fdegWkn) / fk1;
124  Double_t result_wkn = -fwkn0 * S;
125 
126  return result_wkn;
127 }
128 
129 Double_t CbmStsWkn::GetStsWkn(const CbmStsTrack* StsTrack) {
130  double dr = 1.;
131  std::vector<float> eLossVector;
132  std::vector<float> dEdxAllveto;
133 
134  int nClustersWveto =
135  StsTrack->GetNofStsHits()
136  + StsTrack->GetNofStsHits(); //assume all clusters with veto
137  if (nClustersWveto < 8) return fwkn;
138 
139  //cout<<fnSet<<endl;
140  for (int iHit = 0; iHit < StsTrack->GetNofStsHits(); ++iHit) {
141  Int_t StsHitIndex = StsTrack->GetStsHitIndex(iHit);
142  CbmStsHit* stsHit = (CbmStsHit*) fStsHits->At(StsHitIndex);
143 
144  // ***********dx=dr is calculated from the track inclination
145  double x, y, z, xNext, yNext, zNext;
146  x = stsHit->GetX();
147  y = stsHit->GetY();
148  z = stsHit->GetZ();
149 
150  if (iHit != StsTrack->GetNofStsHits() - 1) {
151  CbmStsHit* stsHitNext =
152  (CbmStsHit*) fStsHits->At(StsTrack->GetStsHitIndex(iHit + 1));
153  xNext = stsHitNext->GetX();
154  yNext = stsHitNext->GetY();
155  zNext = stsHitNext->GetZ();
156  dr = sqrt((xNext - x) * (xNext - x) + (yNext - y) * (yNext - y)
157  + (zNext - z) * (zNext - z))
158  / (zNext - z); // if *300um, you get real reconstructed dr
159  } // else dr stay previous
160  // ***********End of dr calculation
161 
162  //************dE is defined as a total cluster charge
163  Int_t FrontClusterId = stsHit->GetFrontClusterId();
164  CbmStsCluster* frontCluster =
165  (CbmStsCluster*) fStsClusterArray->At(FrontClusterId);
166  Int_t BackClusterId = stsHit->GetBackClusterId();
167  CbmStsCluster* backCluster =
168  (CbmStsCluster*) fStsClusterArray->At(BackClusterId);
169 
170  if (!frontCluster || !backCluster) return fwkn;
171 
172  dEdxAllveto.push_back((frontCluster->GetCharge()) / dr);
173  dEdxAllveto.push_back((backCluster->GetCharge()) / dr);
174  }
175 
176  if (dEdxAllveto.size() != 0) {
177  unsigned int NSample = dEdxAllveto.size();
178 
179  for (unsigned int jVec = 0; jVec < NSample; jVec++)
180  dEdxAllveto[jVec] = dEdxAllveto[jVec] / 10000;
181 
182  for (unsigned int jVec = 0; jVec < NSample; jVec++)
183  dEdxAllveto[jVec] = (dEdxAllveto[jVec] - fEmp) / fXi - 0.225;
184 
185  sort(dEdxAllveto.begin(), dEdxAllveto.end());
186 
187  for (unsigned int jVec = 0; jVec < NSample; jVec++)
188  dEdxAllveto[jVec] = TMath::LandauI(dEdxAllveto[jVec]);
189 
190 
191  for (int iHit = 0; iHit < fnSet; iHit++)
192  eLossVector.push_back(dEdxAllveto[NSample - fnSet + iHit]);
193  }
194 
195  Double_t S = 0, ty = 0, ti = 0;
196 
197  for (Int_t i = 0; i < fnSet; i++) {
198  ty = eLossVector[i];
199  ti = i;
200  S += pow((ti - 1) / fnSet - ty, fk1) - pow(ti / fnSet - ty, fk1);
201  }
202  float fwkn0 = pow(fnSet, 0.5 * fdegWkn) / fk1;
203  Double_t result_wkn = -fwkn0 * S;
204 
205  return result_wkn;
206 }
207 
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmStsCluster
Data class for STS clusters.
Definition: CbmStsCluster.h:31
CbmStsWkn::fStsClusterArray
TClonesArray * fStsClusterArray
Definition: CbmStsWkn.h:56
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmStsWkn::~CbmStsWkn
virtual ~CbmStsWkn()
Definition: CbmStsWkn.cxx:34
CbmStsWkn::CbmStsWkn
CbmStsWkn()
Definition: CbmStsWkn.cxx:29
CbmStsWkn::fStsHits
TClonesArray * fStsHits
Definition: CbmStsWkn.h:55
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmStsWkn.h
CbmStsWkn::Init
void Init()
Definition: CbmStsWkn.cxx:36
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmStsWkn::fEmp
float fEmp
Definition: CbmStsWkn.h:47
CbmStsCluster::GetCharge
Double_t GetCharge() const
Get cluster charge @value Total cluster charge [e].
Definition: CbmStsCluster.h:55
CbmGlobalTrack.h
CbmStsWkn::GetStsWkn
Double_t GetStsWkn(Int_t StsTrackIndex)
Return Wkn value.
Definition: CbmStsWkn.cxx:47
CbmStsHit::GetFrontClusterId
Int_t GetFrontClusterId() const
Definition: CbmStsHit.h:93
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition: CbmStsHit.h:31
CbmStsDigi.h
CbmStsTrack.h
Data class for STS tracks.
CbmStsWkn
Definition: CbmStsWkn.h:16
ClassImp
ClassImp(CbmStsWkn)
CbmStsWkn::fXi
float fXi
Definition: CbmStsWkn.h:48
CbmStsWkn::fStsDigiArray
TClonesArray * fStsDigiArray
Definition: CbmStsWkn.h:57
CbmStsWkn::fStsTracks
TClonesArray * fStsTracks
Definition: CbmStsWkn.h:54
CbmStsWkn::fnSet
Int_t fnSet
Definition: CbmStsWkn.h:50
CbmStsWkn::fGlobalTracks
TClonesArray * fGlobalTracks
Definition: CbmStsWkn.h:53
CbmStsWkn::fwkn
Double_t fwkn
Definition: CbmStsWkn.h:51
CbmStsHit::GetBackClusterId
Int_t GetBackClusterId() const
Definition: CbmStsHit.h:69
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmStsWkn::fk1
float fk1
Definition: CbmStsWkn.h:49
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmStsWkn::fdegWkn
Int_t fdegWkn
Definition: CbmStsWkn.h:46
CbmStsTrack::GetStsHitIndex
Int_t GetStsHitIndex(Int_t iHit) const
Definition: CbmStsTrack.h:98
CbmStsCluster.h
Data class for STS clusters.
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsTrack::GetNofStsHits
Int_t GetNofStsHits() const
Definition: CbmStsTrack.h:90