Fix Coverity 10974-82
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoQATrackCut.cxx
CommitLineData
f9210de7 1/***************************************************************************
2 *
3 * $Id: AliFemtoQATrackCut.cxx 24360 2008-03-10 09:48:27Z akisiel $
4 *
5 *
6 ***************************************************************************
7 *
8 *
9 *
10 *
11 ***************************************************************************
12 *
13 * $Log$
14 * Revision 1.3 2007/05/22 09:01:42 akisiel
15 * Add the possibiloity to save cut settings in the ROOT file
16 *
17 * Revision 1.2 2007/05/21 10:38:25 akisiel
18 * More coding rule conformance
19 *
20 * Revision 1.1 2007/05/16 10:25:06 akisiel
21 * Making the directory structure of AliFemtoUser flat. All files go into one common directory
22 *
23 * Revision 1.4 2007/05/03 09:46:10 akisiel
24 * Fixing Effective C++ warnings
25 *
26 * Revision 1.3 2007/04/27 07:25:59 akisiel
27 * Make revisions needed for compilation from the main AliRoot tree
28 *
29 * Revision 1.1.1.1 2007/04/25 15:38:41 panos
30 * Importing the HBT code dir
31 *
32 * Revision 1.4 2007-04-03 16:00:08 mchojnacki
33 * Changes to iprove memory managing
34 *
35 * Revision 1.3 2007/03/13 15:30:03 mchojnacki
36 * adding reader for simulated data
37 *
38 * Revision 1.2 2007/03/08 14:58:03 mchojnacki
39 * adding some alice stuff
40 *
41 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
42 * First version on CVS
43 *
44 **************************************************************************/
45
46#include "AliFemtoQATrackCut.h"
47#include <cstdio>
48
49#ifdef __ROOT__
50ClassImp(AliFemtoQATrackCut)
51#endif
52
53
54// electron
55// 0.13 - 1.8
56// 0 7.594129e-02 8.256141e-03
57// 1 -5.535827e-01 8.170825e-02
58// 2 1.728591e+00 3.104210e-01
59// 3 -2.827893e+00 5.827802e-01
60// 4 2.503553e+00 5.736207e-01
61// 5 -1.125965e+00 2.821170e-01
62// 6 2.009036e-01 5.438876e-02
63
64// pion
65// 0.13 - 2.0
66// 0 1.063457e+00 8.872043e-03
67// 1 -4.222208e-01 2.534402e-02
68// 2 1.042004e-01 1.503945e-02
69
70// kaon
71// 0.18 - 2.0
72// 0 -7.289406e-02 1.686074e-03
73// 1 4.415666e-01 1.143939e-02
74// 2 -2.996790e-01 1.840964e-02
75// 3 6.704652e-02 7.783990e-03
76
77// proton
78// 0.26 - 2.0
79// 0 -3.730200e-02 2.347311e-03
80// 1 1.163684e-01 1.319316e-02
81// 2 8.354116e-02 1.997948e-02
82// 3 -4.608098e-02 8.336400e-03
83
84
85AliFemtoQATrackCut::AliFemtoQATrackCut() :
86 fCharge(0),
87 fLabel(0),
88 fStatus(0),
89 fminTPCclsF(0),
90 fminTPCncls(0),
91 fminITScls(0),
1e79fae4 92 fminTPCchiNdof(0),
f9210de7 93 fMaxTPCncls(1000),
94 fMaxITSchiNdof(1000.0),
95 fMaxTPCchiNdof(1000.0),
96 fMaxSigmaToVertex(1000.0),
97 fNTracksPassed(0),
98 fNTracksFailed(0),
99 fRemoveKinks(kFALSE),
100 fMostProbable(0),
1e79fae4 101 fTPCnclsExclusionSwitch(kFALSE),
102 fTPCchiNdofExclusionSwitch(kFALSE)
f9210de7 103{
104 // Default constructor
105 fNTracksPassed = fNTracksFailed = 0;
106 fCharge = 0; // takes both charges 0
107 fPt[0]=0.0; fPt[1] = 100.0;//100
108 fRapidity[0]=-2; fRapidity[1]=2;//-2 2
109 fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
110 fPidProbPion[0]=-1; fPidProbPion[1]=2;
111 fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
112 fPidProbProton[0]=-1;fPidProbProton[1]=2;
113 fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
114 fLabel=false;
115 fStatus=0;
116 fminTPCclsF=0;
117 fminITScls=0;
118 fTPCnclsExclusionSwitch = false;
119 fTPCnclsExclusion[0] = 0;
120 fTPCnclsExclusion[1] = 1000;
121 fTPCchiNdofExclusionSwitch = false;
122 fTPCchiNdofExclusion[0] = 0.0;
123 fTPCchiNdofExclusion[1] = 1000.0;
124}
125//------------------------------
126AliFemtoQATrackCut::~AliFemtoQATrackCut(){
127 /* noop */
128}
129//------------------------------
130bool AliFemtoQATrackCut::Pass(const AliFemtoTrack* track)
131{
132 // test the particle and return
133 // true if it meets all the criteria
134 // false if it doesn't meet at least one of the criteria
135 float tMost[5];
136
137 //cout<<"AliFemtoESD cut"<<endl;
138 //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
139 if (fStatus!=0)
140 {
141 //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
142 if ((track->Flags()&fStatus)!=fStatus)
143 {
144 // cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
145 return false;
146 }
147
148 }
149 if (fRemoveKinks) {
150 if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
151 return false;
152 }
153 if (fminTPCclsF>track->TPCnclsF())
154 {
155 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
156 return false;
157 }
158
159 // TPC number of clusters:
160 if (fTPCnclsExclusionSwitch) {
161 bool outTPCnclsExclusionZone[2];
162 outTPCnclsExclusionZone[0] = false;
163 outTPCnclsExclusionZone[1] = false;
164 if ( (fminTPCncls > track->TPCncls()) || (fTPCnclsExclusion[0] < track->TPCncls()) ) {
165 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
166 outTPCnclsExclusionZone[0] = true;
167 }
168 if ( (fMaxTPCncls < track->TPCncls()) || (fTPCnclsExclusion[1] > track->TPCncls()) ) {
169 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
170 outTPCnclsExclusionZone[1] = true;
171 }
172 if ( outTPCnclsExclusionZone[0] * outTPCnclsExclusionZone[1] ) { return false; }
173 }
174 else {
175 if (fminTPCncls > track->TPCncls()) {
176 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
177 return false;
178 }
179 if (fMaxTPCncls < track->TPCncls()) {
180 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
181 return false;
182 }
183 }
184
185 if (fminITScls>track->ITSncls())
186 {
187 //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
188 return false;
189 }
190
191 if (fMaxSigmaToVertex < track->SigmaToVertex()) {
192 return false;
193 }
194
195 if (track->ITSncls() > 0)
196 if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
197 return false;
198 }
199
200 // TPC chiNdof of tracks:
201 if (fTPCchiNdofExclusionSwitch && (track->TPCncls() > 0)) {
202 bool outTPCchiNdofExclusionZone[2];
203 outTPCchiNdofExclusionZone[0] = false;
204 outTPCchiNdofExclusionZone[1] = false;
205 if ( (fminTPCchiNdof > (track->TPCchi2()/track->TPCncls())) || (fTPCchiNdofExclusion[0] < (track->TPCchi2()/track->TPCncls())) ) {
206 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
207 outTPCchiNdofExclusionZone[0] = true;
208 }
209 if ( (fMaxTPCchiNdof < (track->TPCchi2()/track->TPCncls())) || (fTPCchiNdofExclusion[1] > (track->TPCchi2()/track->TPCncls())) ) {
210 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
211 outTPCchiNdofExclusionZone[1] = true;
212 }
213 if ( outTPCchiNdofExclusionZone[0] * outTPCchiNdofExclusionZone[1] ) { return false; }
214 }
215 else {
216 if (fminTPCchiNdof > (track->TPCchi2()/track->TPCncls())) {
217 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
218 return false;
219 }
220 if (fMaxTPCchiNdof < (track->TPCchi2()/track->TPCncls())) {
221 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
222 return false;
223 }
224 }
225
226 if (fLabel)
227 {
228 //cout<<"labels"<<endl;
229 if(track->Label()<0)
230 {
231 fNTracksFailed++;
232 // cout<<"No Go Through the cut"<<endl;
233 // cout<<fLabel<<" Label="<<track->Label()<<endl;
234 return false;
235 }
236 }
237 if (fCharge!=0)
238 {
239 //cout<<"AliFemtoESD cut ch "<<endl;
240 //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
241 if (track->Charge()!= fCharge)
242 {
243 fNTracksFailed++;
244 // cout<<"No Go Through the cut"<<endl;
245 // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
246 return false;
247 }
248 }
f6944668 249 float tEnergy = ::sqrt(track->P().Mag2()+fMass*fMass);
f9210de7 250 float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
251 float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
252 if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
253 {
254 fNTracksFailed++;
255 //cout<<"No Go Through the cut"<<endl;
256 //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
257 return false;
258 }
259 if ((tPt<fPt[0])||(tPt>fPt[1]))
260 {
261 fNTracksFailed++;
262 //cout<<"No Go Through the cut"<<endl;
263 //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
264 return false;
265 }
266// cout << "Track has pids: "
267// << track->PidProbElectron() << " "
268// << track->PidProbMuon() << " "
269// << track->PidProbPion() << " "
270// << track->PidProbKaon() << " "
271// << track->PidProbProton() << " "
272// << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
273
274
275 if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
276 {
277 fNTracksFailed++;
278 //cout<<"No Go Through the cut"<<endl;
279 //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
280 return false;
281 }
282 if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
283 {
284 fNTracksFailed++;
285 //cout<<"No Go Through the cut"<<endl;
286 //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
287 return false;
288 }
289 if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
290 {
291 fNTracksFailed++;
292 //cout<<"No Go Through the cut"<<endl;
293 //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
294 return false;
295 }
296 if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
297 {
298 fNTracksFailed++;
299 //cout<<"No Go Through the cut"<<endl;
300 //cout<<fPidProbProton[0]<<" < p ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
301 return false;
302 }
303 if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
304 {
305 fNTracksFailed++;
306 //cout<<"No Go Through the cut"<<endl;
307 //cout<<fPidProbMuon[0]<<" < mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
308 return false;
309 }
310
311 if (fMostProbable) {
f6944668 312 tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().Mag());
f9210de7 313 tMost[1] = 0.0;
f6944668 314 tMost[2] = track->PidProbPion()*PidFractionPion(track->P().Mag());
315 tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().Mag());
316 tMost[4] = track->PidProbProton()*PidFractionProton(track->P().Mag());
f9210de7 317 int imost=0;
318 float ipidmax = 0.0;
319 for (int ip=0; ip<5; ip++)
320 if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
321 if (imost != fMostProbable) return false;
322 }
323
324 // cout<<"Go Through the cut"<<endl;
325 // cout<<fLabel<<" Label="<<track->Label()<<endl;
326 // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
327 // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
328 //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
329 //cout<<fPidProbElectron[0]<<" < e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
330 //cout<<fPidProbPion[0]<<" < pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
331 //cout<<fPidProbKaon[0]<<" < k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
332 //cout<<fPidProbProton[0]<<" < p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
333 //cout<<fPidProbMuon[0]<<" < mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
334 fNTracksPassed++ ;
335 return true;
336
337
338}
339//------------------------------
340AliFemtoString AliFemtoQATrackCut::Report()
341{
342 // Prepare report from the execution
343 string tStemp;
344 char tCtemp[100];
345 sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
346 tStemp=tCtemp;
347 sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
348 tStemp+=tCtemp;
349 sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
350 tStemp+=tCtemp;
351 sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
352 tStemp+=tCtemp;
353 sprintf(tCtemp,"Number of tracks which passed:\t%ld Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
354 tStemp += tCtemp;
355 AliFemtoString returnThis = tStemp;
356 return returnThis;
357}
358TList *AliFemtoQATrackCut::ListSettings()
359{
360 // return a list of settings in a writable form
361 TList *tListSetttings = new TList();
362 char buf[200];
363 snprintf(buf, 200, "AliFemtoQATrackCut.mass=%f", this->Mass());
364 tListSetttings->AddLast(new TObjString(buf));
365
366 snprintf(buf, 200, "AliFemtoQATrackCut.charge=%i", fCharge);
367 tListSetttings->AddLast(new TObjString(buf));
368 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobpion.minimum=%f", fPidProbPion[0]);
369 tListSetttings->AddLast(new TObjString(buf));
370 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobpion.maximum=%f", fPidProbPion[1]);
371 tListSetttings->AddLast(new TObjString(buf));
372 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobkaon.minimum=%f", fPidProbKaon[0]);
373 tListSetttings->AddLast(new TObjString(buf));
374 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobkaon.maximum=%f", fPidProbKaon[1]);
375 tListSetttings->AddLast(new TObjString(buf));
376 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobproton.minimum=%f", fPidProbProton[0]);
377 tListSetttings->AddLast(new TObjString(buf));
378 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobproton.maximum=%f", fPidProbProton[1]);
379 tListSetttings->AddLast(new TObjString(buf));
380 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobelectron.minimum=%f", fPidProbElectron[0]);
381 tListSetttings->AddLast(new TObjString(buf));
382 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobelectron.maximum=%f", fPidProbElectron[1]);
383 tListSetttings->AddLast(new TObjString(buf));
384 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobMuon.minimum=%f", fPidProbMuon[0]);
385 tListSetttings->AddLast(new TObjString(buf));
386 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobMuon.maximum=%f", fPidProbMuon[1]);
387 tListSetttings->AddLast(new TObjString(buf));
388 snprintf(buf, 200, "AliFemtoQATrackCut.minimumtpcclusters=%i", fminTPCclsF);
389 tListSetttings->AddLast(new TObjString(buf));
390 snprintf(buf, 200, "AliFemtoQATrackCut.minimumitsclusters=%i", fminTPCclsF);
391 tListSetttings->AddLast(new TObjString(buf));
392 snprintf(buf, 200, "AliFemtoQATrackCut.pt.minimum=%f", fPt[0]);
393 tListSetttings->AddLast(new TObjString(buf));
394 snprintf(buf, 200, "AliFemtoQATrackCut.pt.maximum=%f", fPt[1]);
395 tListSetttings->AddLast(new TObjString(buf));
396 snprintf(buf, 200, "AliFemtoQATrackCut.rapidity.minimum=%f", fRapidity[0]);
397 tListSetttings->AddLast(new TObjString(buf));
398 snprintf(buf, 200, "AliFemtoQATrackCut.rapidity.maximum=%f", fRapidity[1]);
399 tListSetttings->AddLast(new TObjString(buf));
400 snprintf(buf, 200, "AliFemtoQATrackCut.removekinks=%i", fRemoveKinks);
401 tListSetttings->AddLast(new TObjString(buf));
402 snprintf(buf, 200, "AliFemtoQATrackCut.maxitschindof=%f", fMaxITSchiNdof);
403 tListSetttings->AddLast(new TObjString(buf));
404 snprintf(buf, 200, "AliFemtoQATrackCut.maxtpcchindof=%f", fMaxTPCchiNdof);
405 tListSetttings->AddLast(new TObjString(buf));
406 snprintf(buf, 200, "AliFemtoQATrackCut.maxsigmatovertex=%f", fMaxSigmaToVertex);
407 tListSetttings->AddLast(new TObjString(buf));
408 if (fMostProbable) {
409 if (fMostProbable == 2)
410 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Pion");
411 if (fMostProbable == 3)
412 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Kaon");
413 if (fMostProbable == 4)
414 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Proton");
415 tListSetttings->AddLast(new TObjString(buf));
416 }
417 return tListSetttings;
418}
419void AliFemtoQATrackCut::SetRemoveKinks(const bool& flag)
420{
421 fRemoveKinks = flag;
422}
423
424 // electron
425// 0.13 - 1.8
426// 0 7.594129e-02 8.256141e-03
427// 1 -5.535827e-01 8.170825e-02
428// 2 1.728591e+00 3.104210e-01
429// 3 -2.827893e+00 5.827802e-01
430// 4 2.503553e+00 5.736207e-01
431// 5 -1.125965e+00 2.821170e-01
432// 6 2.009036e-01 5.438876e-02
433float AliFemtoQATrackCut::PidFractionElectron(float mom) const
434{
435 // Provide a parameterized fraction of electrons dependent on momentum
436 if (mom<0.13) return 0.0;
437 if (mom>1.8) return 0.0;
438 return (7.594129e-02
439 -5.535827e-01*mom
440 +1.728591e+00*mom*mom
441 -2.827893e+00*mom*mom*mom
442 +2.503553e+00*mom*mom*mom*mom
443 -1.125965e+00*mom*mom*mom*mom*mom
444 +2.009036e-01*mom*mom*mom*mom*mom*mom);
445}
446
447// pion
448// 0.13 - 2.0
449// 0 1.063457e+00 8.872043e-03
450// 1 -4.222208e-01 2.534402e-02
451// 2 1.042004e-01 1.503945e-02
452float AliFemtoQATrackCut::PidFractionPion(float mom) const
453{
454 // Provide a parameterized fraction of pions dependent on momentum
455 if (mom<0.13) return 0.0;
456 if (mom>2.0) return 0.0;
457 return ( 1.063457e+00
458 -4.222208e-01*mom
459 +1.042004e-01*mom*mom);
460}
461
462// kaon
463// 0.18 - 2.0
464// 0 -7.289406e-02 1.686074e-03
465// 1 4.415666e-01 1.143939e-02
466// 2 -2.996790e-01 1.840964e-02
467// 3 6.704652e-02 7.783990e-03
468float AliFemtoQATrackCut::PidFractionKaon(float mom) const
469{
470 // Provide a parameterized fraction of kaons dependent on momentum
471 if (mom<0.18) return 0.0;
472 if (mom>2.0) return 0.0;
473 return (-7.289406e-02
474 +4.415666e-01*mom
475 -2.996790e-01*mom*mom
476 +6.704652e-02*mom*mom*mom);
477}
478
479// proton
480// 0.26 - 2.0
481// 0 -3.730200e-02 2.347311e-03
482// 1 1.163684e-01 1.319316e-02
483// 2 8.354116e-02 1.997948e-02
484// 3 -4.608098e-02 8.336400e-03
485float AliFemtoQATrackCut::PidFractionProton(float mom) const
486{
487 // Provide a parameterized fraction of protons dependent on momentum
488 if (mom<0.26) return 0.0;
489 if (mom>2.0) return 0.0;
490 return (-3.730200e-02
491 +1.163684e-01*mom
492 +8.354116e-02*mom*mom
493 -4.608098e-02*mom*mom*mom);
494}