CbmRoot
CbmTrdSetTracksPidModWkn.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- CbmTrdSetTracksPidModWkn source file -----
3 // ----- -----
4 // -------------------------------------------------------------------------
6 
7 #include "CbmMCTrack.h"
8 #include "CbmTrackMatch.h"
9 #include "CbmTrdHit.h"
10 #include "CbmTrdPoint.h"
11 #include "CbmTrdTrack.h"
12 
13 #include "FairRootManager.h"
14 
15 #include "TClonesArray.h"
16 #include "TMath.h"
17 #include "TStopwatch.h"
18 #include <iomanip>
19 #include <iostream>
20 
21 #ifdef HAVE_SSE
22 #include "P4_F32vec4.h"
23 #else
24 #include "PSEUDO_F32vec4.h"
25 #error NoSseFound
26 #endif // HAVE_SSE
27 
28 #include <cmath>
29 #include <vector>
30 using std::cout;
31 using std::endl;
32 
33 // ----- Default constructor -------------------------------------------
35  : CbmTrdSetTracksPidModWkn("TrdSetTracksPidModWkn", "") {}
36 // -------------------------------------------------------------------------
37 
38 // ----- Standard constructor ------------------------------------------
40  const char*)
41  : FairTask(name)
42  , fTrackArray(NULL)
43  , fTrdHitArray(NULL)
44  , fnSet(0)
45  , fdegWkn(0)
46  , fk1(0)
47  , fwkn0(0)
48  , fEmp(0)
49  , fXi(0)
50  , fWmin(0)
51  , fWmax(0)
52  , fDiff(0)
53  , fSISType("sis300") {}
54 // -------------------------------------------------------------------------
55 
56 
57 // ----- Destructor ----------------------------------------------------
59 // -------------------------------------------------------------------------
60 
61 
62 // ----- Public method Init (abstract in base class) --------------------
64 
65  // Get and check FairRootManager
66  FairRootManager* ioman = FairRootManager::Instance();
67  if (!ioman) {
68  cout << "-E- CbmTrdSetTracksPidModWkn::Init: "
69  << "RootManager not instantised!" << endl;
70  return kFATAL;
71  }
72 
73  // Get TrdTrack array
74  fTrackArray = (TClonesArray*) ioman->GetObject("TrdTrack");
75  if (!fTrackArray) {
76  cout << "-E- CbmTrdSetTracksPidModWkn::Init: No TrdTrack array!" << endl;
77  return kERROR;
78  }
79 
80  // Get TrdHit array
81  fTrdHitArray = (TClonesArray*) ioman->GetObject("TrdHit");
82  if (!fTrdHitArray) {
83  cout << "-E- CbmTrdSetTracksPidModWkn::Init: No TrdHit array!" << endl;
84  return kERROR;
85  }
86 
87  SetParameters();
88  // fmom = new TH1D("(p_{rec}-p_{mc})/p_{mc}","",400, -1, 1);
89  return kSUCCESS;
90 }
91 // -------------------------------------------------------------------------
92 
93 
94 // ----- Public method Exec --------------------------------------------
96 
97  TStopwatch timer;
98  timer.Start();
99 
100  if (!fTrackArray) return;
101 
102  float ti = 0, mom = 0;
103 
104  std::vector<float> eLossVector; // vector for energy losses
105 
106  fvec* mWkn = new fvec[fnSet];
107  fvec resWkn = 0;
108  fvec numTr = 0;
109 
110  // int nTracks_SIMD = fvecLen;
111  int NHits = 0;
112  int iV = 0;
113 
114  unsigned short nTracks = fTrackArray->GetEntriesFast();
115  unsigned short nTrue = 0;
116 
117  for (unsigned short itr = 0; itr < nTracks; itr++) {
118  CbmTrdTrack* pTr = (CbmTrdTrack*) fTrackArray->At(itr);
119  if ((pTr->GetNofHits()) > fnSet - 1)
120  nTrue++;
121  else
122  pTr->SetPidWkn(-2.);
123  }
124 
125  for (unsigned short itrack = 0; itrack < nTrue; itrack++) {
126  eLossVector.clear();
127 
128  CbmTrdTrack* pTrack = (CbmTrdTrack*) fTrackArray->At(itrack);
129  NHits = pTrack->GetNofHits();
130  numTr[iV] = itrack;
131  mom = fabs(1.f / pTrack->GetParamFirst()->GetQp());
132 
133  fEmp = 0.0002598 * mom * mom * mom - 0.008862 * mom * mom + 0.1176 * mom
134  + 0.9129;
135  fXi = 0.00008938 * mom * mom * mom - 0.003022 * mom * mom + 0.03999 * mom
136  + 0.5292;
137 
138  for (Int_t iHit = 0; iHit < NHits; iHit++) {
139  Int_t TRDindex = pTrack->GetHitIndex(iHit);
140  CbmTrdHit* trdHit = (CbmTrdHit*) fTrdHitArray->At(TRDindex);
141  eLossVector.push_back((trdHit->GetELoss()) * 1000000);
142  }
143 
144  for (unsigned int jVec = 0; jVec < eLossVector.size(); jVec++)
145  eLossVector[jVec] = (eLossVector[jVec] - fEmp) / fXi - 0.225;
146 
147  sort(eLossVector.begin(), eLossVector.end());
148 
149  for (unsigned int jVec = 0; jVec < eLossVector.size(); jVec++)
150  eLossVector[jVec] = TMath::LandauI(eLossVector[jVec]);
151 
152  for (Int_t iHit = 0; iHit < fnSet; iHit++)
153  mWkn[iHit][iV] = eLossVector[NHits - fnSet + iHit];
154 
155  iV++;
156 
157  if (!((iV == fvecLen) || ((itrack == (nTrue - 1)) && (iV != 0)))) continue;
158 
159  fvec sumWkn = 0;
160  for (Int_t iHit = 0; iHit < fnSet; iHit++) {
161  ti = iHit + 1;
162  fvec g1 = (ti - 1) / fnSet - mWkn[iHit];
163  fvec g2 = ti / fnSet - mWkn[iHit];
164  fvec s1 = 1, s2 = 1;
165  for (int iPow = 0; iPow < 5; iPow++) {
166  s1 *= g1;
167  s2 *= g2;
168  }
169  sumWkn += s1 - s2;
170  }
171 
172  // resWkn = -fwkn0*sumWkn;
173  // resWkn = 2*(-fwkn0*sumWkn-fWmin)/fDiff-1; //[-1;1]
174  resWkn = (-fwkn0 * sumWkn - fWmin) / fDiff; //[0;1]
175 
176 
177  for (Int_t iHit = 0; iHit < iV; iHit++) {
178  CbmTrdTrack* pTrack2 = (CbmTrdTrack*) fTrackArray->At(numTr[iHit]);
179  pTrack2->SetPidWkn(resWkn[iHit]);
180  }
181 
182  iV = 0;
183  }
184 
185  delete[] mWkn;
186 
187  timer.Stop();
188  static int nEv = 0;
189  nEv++;
190  static Double_t rtime = 0;
191  rtime += timer.RealTime();
192  static Double_t ctime = 0;
193  ctime += timer.CpuTime();
194  cout << endl << endl;
195  cout << "--CbmTrdSetTracksPidModWkn--" << endl;
196  cout << "Real time " << (rtime / Double_t(nEv)) << " s, CPU time "
197  << (ctime / Double_t(nEv)) << "s" << endl
198  << endl;
199 }
200 
201 
203  if (fSISType == "sis300") {
204  fnSet = 5; // number of the layers with TR
205  fdegWkn = 4; // statistics degree
206  }
207 
208  if (fSISType == "sis100") {
209  fnSet = 2; // number of the layers with TR
210  fdegWkn = 2; // statistics degree
211  }
212 
213  fk1 = fdegWkn + 1;
214  fwkn0 = pow(fnSet, 0.5 * fdegWkn) / fk1;
215  fWmin = 1 / (pow(2, fdegWkn) * pow(fnSet, fdegWkn / 2) * (fdegWkn + 1)),
216  fWmax = pow(fnSet, fdegWkn / 2) / (fdegWkn + 1), fDiff = fWmax - fWmin;
217 }
218 
219 
220 // ----- Public method Finish ------------------------------------------
222  //fmom->Write();
223 }
224 // -------------------------------------------------------------------------
225 
CbmTrdSetTracksPidModWkn::fDiff
float fDiff
Definition: CbmTrdSetTracksPidModWkn.h:84
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
CbmTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmTrack.h:53
CbmTrdHit::GetELoss
Double_t GetELoss() const
Definition: CbmTrdHit.h:79
CbmTrdHit
data class for a reconstructed Energy-4D measurement in the TRD
Definition: CbmTrdHit.h:35
CbmTrdSetTracksPidModWkn::CbmTrdSetTracksPidModWkn
CbmTrdSetTracksPidModWkn()
Definition: CbmTrdSetTracksPidModWkn.cxx:34
CbmTrdSetTracksPidModWkn::Finish
virtual void Finish()
Definition: CbmTrdSetTracksPidModWkn.cxx:221
CbmTrdSetTracksPidModWkn::fwkn0
float fwkn0
Definition: CbmTrdSetTracksPidModWkn.h:81
CbmTrackMatch.h
CbmTrdSetTracksPidModWkn::fTrackArray
TClonesArray * fTrackArray
Definition: CbmTrdSetTracksPidModWkn.h:75
CbmTrdSetTracksPidModWkn.h
CbmTrdSetTracksPidModWkn::SetParameters
void SetParameters()
Definition: CbmTrdSetTracksPidModWkn.cxx:202
CbmTrdSetTracksPidModWkn::fk1
float fk1
Definition: CbmTrdSetTracksPidModWkn.h:81
CbmTrdSetTracksPidModWkn::fTrdHitArray
TClonesArray * fTrdHitArray
Definition: CbmTrdSetTracksPidModWkn.h:76
CbmTrdSetTracksPidModWkn::fXi
float fXi
Definition: CbmTrdSetTracksPidModWkn.h:83
fvecLen
const int fvecLen
Definition: L1/vectors/P4_F32vec4.h:251
CbmTrdSetTracksPidModWkn
Definition: CbmTrdSetTracksPidModWkn.h:37
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmTrdHit.h
Class for hits in TRD detector.
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmTrdSetTracksPidModWkn::~CbmTrdSetTracksPidModWkn
virtual ~CbmTrdSetTracksPidModWkn()
Definition: CbmTrdSetTracksPidModWkn.cxx:58
CbmTrdSetTracksPidModWkn::Init
virtual InitStatus Init()
Definition: CbmTrdSetTracksPidModWkn.cxx:63
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmTrdSetTracksPidModWkn::fWmax
float fWmax
Definition: CbmTrdSetTracksPidModWkn.h:84
CbmTrdTrack
Definition: CbmTrdTrack.h:22
CbmTrdTrack::SetPidWkn
void SetPidWkn(Double_t pid)
Definition: CbmTrdTrack.h:41
CbmTrdSetTracksPidModWkn::fEmp
float fEmp
Definition: CbmTrdSetTracksPidModWkn.h:82
CbmTrdSetTracksPidModWkn::fnSet
int fnSet
Definition: CbmTrdSetTracksPidModWkn.h:79
CbmMCTrack.h
PSEUDO_F32vec4.h
CbmTrdSetTracksPidModWkn::fdegWkn
float fdegWkn
Definition: CbmTrdSetTracksPidModWkn.h:80
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
CbmTrdSetTracksPidModWkn::Exec
virtual void Exec(Option_t *opt)
Definition: CbmTrdSetTracksPidModWkn.cxx:95
CbmTrdPoint.h
CbmTrdSetTracksPidModWkn::fWmin
float fWmin
Definition: CbmTrdSetTracksPidModWkn.h:84
CbmTrdSetTracksPidModWkn::fSISType
std::string fSISType
Definition: CbmTrdSetTracksPidModWkn.h:91
CbmTrdTrack.h