BmnRoot
Loading...
Searching...
No Matches
BmnPidSRC.cxx
Go to the documentation of this file.
1#include "BmnPidSRC.h"
2
3static Float_t workTime = 0.0;
4
5using namespace std;
6using namespace TMath;
7
8Int_t nFoundAZ = 0;
9Int_t nFoundPZ = 0;
10
11
12 const Double_t AzUpBorder[8][14] = { // so called "pupochki"
13//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
14 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
15 {0, 1.11, 2.17, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
16 {0, 0, 0, 1.53, 2.15, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
17 {0, 0, 0, 0, 0, 1.55, 2.15, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
18 {0, 0, 0, 0, 0, 0, 0, 1.93, 0, 2.50, 0, 0, 0, 0}, // for Z=4
19 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.15, 2.30, 0, 0}, // for Z=5
20 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.00, 2.08, 0}, // for Z=6
21 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
22 };
23 const Double_t AzDownBorder[8][14] = {
24//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
25
26 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
27 {0, 0, 0.83, 1.83, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
28 {0, 0, 0, 1.01, 1.85, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
29 {0, 0, 0, 0, 0, 1.03, 1.85, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
30 {0, 0, 0, 0, 0, 0, 0, 1.45, 0, 2.06, 0, 0, 0, 0}, // for Z=4
31 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.91, 2.1, 0, 0}, // for Z=5
32 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.68, 1.92, 0}, // for Z=6
33 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
34 };
35
36 const Double_t ZAOutUpBorder[8][14] = {
37//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
38 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
39 {0, 1.3, 1.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
40 {0, 0, 0, 2.6, 2.6, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
41 {0, 0, 0, 0, 0, 3.3, 3.3, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
42 {0, 0, 0, 0, 0, 0, 0, 4.2, 0, 4.2, 0, 0, 0, 0}, // for Z=4
43 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.25, 5.25, 0, 0}, // for Z=5
44 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.3, 6.3, 0}, // for Z=6
45 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
46 };
47 const Double_t ZAOutDownBorder[8][14] = {
48//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
49 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
50 {0, 0.9, 0.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
51 {0, 0, 0, 2., 2., 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
52 {0, 0, 0, 0, 0, 2.7, 2.7, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
53 {0, 0, 0, 0, 0, 0, 0, 2.6, 0, 2.6, 0, 0, 0, 0}, // for Z=4
54 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.75, 4.75, 0, 0}, // for Z=5
55 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.7, 5.7, 0}, // for Z=6
56 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
57 };
58
59 const Double_t RigidityUpBorder[8][14] = {
60//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
61 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
62 {0, 0, 9.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
63 {0, 0, 0, 6.9, 8.9, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
64 {0, 0, 0, 0, 0, 6.9, 7, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
65 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=4
66 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=5
67 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8.7, 0}, // for Z=6
68 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
69 };
70 const Double_t RigidityDownBorder[8][14] = {
71//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
72
73 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
74 {0, 0, 7.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
75 {0, 0, 0, 5.5, 7.5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
76 {0, 0, 0, 0, 0, 5.5, 7.6, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
77 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=4
78 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=5
79 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7.5, 0}, // for Z=6
80 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
81 };
82 const Double_t ZOutUpBorder[8][14] = {
83//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
84 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
85 {0, 0, 1.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
86 {0, 0, 0, 2.6, 2.6, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
87 {0, 0, 0, 0, 0, 3.3, 3.3, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
88 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=4
89 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=5
90 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.3, 0}, // for Z=6
91 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
92 };
93 const Double_t ZOutDownBorder[8][14] = {
94//A 0 1 2 3 4 5 6 7 8 9 10 11 12 13
95 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=0
96 {0, 0, 1.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=1
97 {0, 0, 0, 2., 2., 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=2
98 {0, 0, 0, 0, 0, 2.7, 2.7, 0, 0, 0, 0, 0, 0, 0}, // for Z=3
99 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=4
100 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, // for Z=5
101 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.7, 0}, // for Z=6
102 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} // for Z=7
103 };
104
105BmnPidSRC::BmnPidSRC() {
106 fEventNo = 0;
107 fGlobalTracksArray = nullptr;
108 fDstEventHeader = nullptr;
109 fBmnDchTrack=nullptr;
110 fBmnTofHit=nullptr;
111 fKalman = nullptr;
112 fZout=0.0;
113 fGlobalTracksBranchName = "BmnGlobalTrack";
114 fDstEventHeaderBranchName = "DstEventHeader.";
115 fDchTrackBranchName = "BmnDchTrack";
116 fTofHitBranchName = "BmnTof700Hit";
117}
118
119InitStatus BmnPidSRC::Init() {
120
121 if (fVerbose > 1) cout << "=========================== Pid init started ====================" << endl;
122
123 //Get ROOT Manager
124 FairRootManager* ioman = FairRootManager::Instance();
125 if (NULL == ioman) Fatal("Init", "FairRootManager is not instantiated");
126
127 fGlobalTracksArray = (TClonesArray*) ioman->GetObject(fGlobalTracksBranchName); //in
128 if (!fGlobalTracksArray) {
129 cout << "BmnPidSRC::Init(): branch " << fGlobalTracksBranchName << " not found! Task will be deactivated" << endl;
130 SetActive(kFALSE);
131 return kERROR;
132 }
133
134 fBmnDchTrack = (TClonesArray*) ioman->GetObject(fDchTrackBranchName); //in
135 if (!fBmnDchTrack) {
136 cout << "BmnPidSRC::Init(): branch " << fDchTrackBranchName << " not found! Task will be deactivated" << endl;
137 SetActive(kFALSE);
138 return kERROR;
139 }
140
141 fBmnTofHit = (TClonesArray*) ioman->GetObject(fTofHitBranchName); //in
142 if (!fBmnTofHit) {
143 cout << "BmnPidSRC::Init(): branch " << fTofHitBranchName << " not found! Task will be deactivated" << endl;
144 SetActive(kFALSE);
145 return kERROR;
146 }
147
148 fDstEventHeader=(DstEventHeader*) ioman->GetObject(fDstEventHeaderBranchName);
149 if (!fDstEventHeader) {
150 cout << "BmnPidSRC::Init(): branch " << fDstEventHeaderBranchName << " not found! Task will be deactivated" << endl;
151 SetActive(kFALSE);
152 return kERROR;
153 }
154
155 //fField = FairRunAna::Instance()->GetField();
156 fKalman = new BmnKalmanFilter();
157
158 hPIDdch = new TH2F("PIDdch", "; total charge of the event Z; A/z", 500, 0, 10, 500, -2, 7);
159 hPIDgem = new TH2F("PIDgem", "; total charge of the event Z; P/q", 500, 0, 10, 500, 0, 15);
160
161 if (fVerbose > 1) cout << "=========================== Pid init finished ===================" << endl<<endl<<endl<<endl<<endl;
162
163 return kSUCCESS;
164}
165
166void BmnPidSRC::Exec(Option_t* opt) {
167 if (!IsActive())
168 return;
169 clock_t tStart = clock();
170
171 AzPID();
172
173// RigidityPID();
174
175 DrawPID();
176
177
178 if (fVerbose > 1) cout << "\n======================== Pid exec finished ======================" << endl;
179
180 clock_t tFinish = clock();
181 workTime += ((Float_t) (tFinish - tStart)) / CLOCKS_PER_SEC;
182}
183
184void BmnPidSRC::RigidityPID(){
185 if (fDstEventHeader->GetZ2out() > 0) {
186 fZout=sqrt(fDstEventHeader->GetZ2out());
187 // fZin=sqrt(fDstEventHeader->GetZ2in());
188
189 // Double_t fZoutTrue=fZout;//*0.797+0.235;
190 // cout << "fZoutTrue = " << fZoutTrue << endl;
191 for (int i=0; i<8; i++){
192 for (int j=0; j<14; j++){
193 if (RigidityDownBorder[i][j]==0) continue;
194 if (fZout>ZOutUpBorder[i][j] || fZout<ZOutDownBorder[i][j]) continue; // Charge cut for isotopes
195 for (int k=0; k<fGlobalTracksArray->GetEntriesFast(); k++){
196 BmnGlobalTrack* track = (BmnGlobalTrack*) fGlobalTracksArray->At(k);
197 //if (track->GetA() > -1) continue;
198 Double_t rigid = 1/track->GetParamLast()->GetQp();
199 //hPIDgem->Fill(fZout, rigid);
200 if (rigid>RigidityUpBorder[i][j] || rigid<RigidityDownBorder[i][j]) continue; // Rigidity cut
201 track->SetA(j);
202 track->SetZ(i);
203 //nFoundPZ++;
204 }
205 }
206
207 }
208 }
209
210}
211
212void BmnPidSRC::AzPID(){
213
214 if (fDstEventHeader->GetZ2out() > 0) {
215 //Double_t fZoutTrue=fZout;//*0.797+0.235;
216 // cout << "fZoutTrue = " << fZoutTrue << endl;
217 for (int i=0; i<8; i++){
218 for (int j=0; j<14; j++){
219
220 if (AzDownBorder[i][j]==0) continue; //unfounded isotopes cut
221
222 fZout=sqrt(fDstEventHeader->GetZ2out());
223 if (fZout>ZAOutUpBorder[i][j] || fZout<ZAOutDownBorder[i][j]) continue; // Charge cut for isotopes
224 //if(fGlobalTracksArray->GetEntriesFast()>1) continue; // no more than one track in entry
225 for (int k=0; k<fGlobalTracksArray->GetEntriesFast(); k++){
226
227 BmnGlobalTrack* track = (BmnGlobalTrack*) fGlobalTracksArray->At(k);
228
229 //Int_t index=track->GetTof2HitIndex();
230 //Int_t index1=track->GetDch1TrackIndex();
231 //if(index==-1 || index1==-1) continue; //DCH Tof cut
232 //fZin=sqrt(fDstEventHeader->GetZ2in());
233 //if (fZin>6.5 || fZin<5.5) continue; // Zin only Carbon
234
235 //for(int m=0; m<fBmnDchTrack->GetEntries(); m++){
236 //BmnDchTrack *dchTrack = (BmnDchTrack *)fBmnDchTrack->At(m);
237 //if(dchTrack->GetParamFirst()->GetZ()<550 || dchTrack->GetParamFirst()->GetZ()>650) continue; // global dch tracks
238 FairTrackParam parPrev(*(track->GetParamLast()));
239 fKalman->TGeoTrackPropagate(&parPrev, 750, 2212, NULL, NULL, kTRUE);
240 Double_t xdch=parPrev.GetX();
241 Double_t Txdch=parPrev.GetTx();
242 Double_t Az=(xdch - Txdch * 1470.88) * 0.026106 + 4.33385;
243 //hPIDdch->Fill(fZout, Az);
244 if (Az>AzUpBorder[i][j] || Az<AzDownBorder[i][j]) continue; // Az cut
245 track->SetA(j);
246 track->SetZ(i);
247 //nFoundAZ++;
248
249 }
250 }
251 }
252 }
253
254}
255
256void BmnPidSRC::DrawPID(){
257
258 if (fDstEventHeader->GetZ2out() > 0) {
259 fZout=sqrt(fDstEventHeader->GetZ2out());
260 // draw Gem-based PID
261 for (int k=0; k<fGlobalTracksArray->GetEntriesFast(); k++){
262 BmnGlobalTrack* track = (BmnGlobalTrack*) fGlobalTracksArray->At(k);
263 Double_t rigid = 1/track->GetParamLast()->GetQp();
264 hPIDgem->Fill(fZout, rigid);
265 }
266
267 // draw Dch-based PID
268 for(int m=0; m<fBmnDchTrack->GetEntries(); m++){
269 BmnDchTrack *dchTrack = (BmnDchTrack *)fBmnDchTrack->At(m);
270 if(dchTrack->GetParamFirst()->GetZ()<550 || dchTrack->GetParamFirst()->GetZ()>650) continue; // global dch tracks
271 FairTrackParam parPrev(*(dchTrack->GetParamFirst()));
272 Double_t xdch=parPrev.GetX();
273 Double_t Txdch=parPrev.GetTx();
274 Double_t Az=(xdch - Txdch * 1470.88) * 0.026106 + 3.508280;
275 hPIDdch->Fill(fZout, Az);
276 }
277 }
278}
279
280void BmnPidSRC::Finish() {
281 delete fKalman;
282 cout << "Work time of PID: " << workTime << endl;
283
284 TFile pidFile("pid.root", "RECREATE");
285
286 pidFile.Add(hPIDdch);
287 pidFile.Add(hPIDgem);
288
289 pidFile.Write();
290 pidFile.Close();
291
292 cout << "nFoundAZ = " << nFoundAZ << " nFoundPZ = " << nFoundPZ << endl;
293}
294
nFoundPZ
Int_t nFoundPZ
Definition BmnPidSRC.cxx:9
AzUpBorder
const Double_t AzUpBorder[8][14]
Definition BmnPidSRC.cxx:12
ZOutUpBorder
const Double_t ZOutUpBorder[8][14]
Definition BmnPidSRC.cxx:82
AzDownBorder
const Double_t AzDownBorder[8][14]
Definition BmnPidSRC.cxx:23
ZAOutUpBorder
const Double_t ZAOutUpBorder[8][14]
Definition BmnPidSRC.cxx:36
ZOutDownBorder
const Double_t ZOutDownBorder[8][14]
Definition BmnPidSRC.cxx:93
nFoundAZ
Int_t nFoundAZ
Definition BmnPidSRC.cxx:8
ZAOutDownBorder
const Double_t ZAOutDownBorder[8][14]
Definition BmnPidSRC.cxx:47
RigidityDownBorder
const Double_t RigidityDownBorder[8][14]
Definition BmnPidSRC.cxx:70
RigidityUpBorder
const Double_t RigidityUpBorder[8][14]
Definition BmnPidSRC.cxx:59
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
BmnGlobalTrack::SetZ
void SetZ(Int_t z)
Definition BmnGlobalTrack.h:95
BmnGlobalTrack::SetA
void SetA(Int_t a)
Definition BmnGlobalTrack.h:94
BmnKalmanFilter::TGeoTrackPropagate
BmnStatus TGeoTrackPropagate(FairTrackParam *par, Double_t zOut, Int_t pdg, vector< Double_t > *F, Double_t *length, Bool_t isField)
Definition BmnKalmanFilter.cxx:480
BmnPidSRC::DrawPID
void DrawPID()
Definition BmnPidSRC.cxx:256
BmnPidSRC::Exec
virtual void Exec(Option_t *opt)
Definition BmnPidSRC.cxx:166
BmnPidSRC::Finish
virtual void Finish()
Definition BmnPidSRC.cxx:280
BmnPidSRC::Init
virtual InitStatus Init()
Definition BmnPidSRC.cxx:119
BmnPidSRC::RigidityPID
void RigidityPID()
Definition BmnPidSRC.cxx:184
BmnPidSRC::AzPID
void AzPID()
Definition BmnPidSRC.cxx:212
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:72
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:76
DstEventHeader::GetZ2out
Double_t GetZ2out()
Definition DstEventHeader.h:69
std
STL namespace.