QA classes from Dave
[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),
92 fMaxTPCncls(1000),
93 fMaxITSchiNdof(1000.0),
94 fMaxTPCchiNdof(1000.0),
95 fMaxSigmaToVertex(1000.0),
96 fNTracksPassed(0),
97 fNTracksFailed(0),
98 fRemoveKinks(kFALSE),
99 fMostProbable(0),
100 fTPCnclsExclusionSwitch(kFALSE)
101{
102 // Default constructor
103 fNTracksPassed = fNTracksFailed = 0;
104 fCharge = 0; // takes both charges 0
105 fPt[0]=0.0; fPt[1] = 100.0;//100
106 fRapidity[0]=-2; fRapidity[1]=2;//-2 2
107 fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
108 fPidProbPion[0]=-1; fPidProbPion[1]=2;
109 fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
110 fPidProbProton[0]=-1;fPidProbProton[1]=2;
111 fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
112 fLabel=false;
113 fStatus=0;
114 fminTPCclsF=0;
115 fminITScls=0;
116 fTPCnclsExclusionSwitch = false;
117 fTPCnclsExclusion[0] = 0;
118 fTPCnclsExclusion[1] = 1000;
119 fTPCchiNdofExclusionSwitch = false;
120 fTPCchiNdofExclusion[0] = 0.0;
121 fTPCchiNdofExclusion[1] = 1000.0;
122}
123//------------------------------
124AliFemtoQATrackCut::~AliFemtoQATrackCut(){
125 /* noop */
126}
127//------------------------------
128bool AliFemtoQATrackCut::Pass(const AliFemtoTrack* track)
129{
130 // test the particle and return
131 // true if it meets all the criteria
132 // false if it doesn't meet at least one of the criteria
133 float tMost[5];
134
135 //cout<<"AliFemtoESD cut"<<endl;
136 //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
137 if (fStatus!=0)
138 {
139 //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
140 if ((track->Flags()&fStatus)!=fStatus)
141 {
142 // cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
143 return false;
144 }
145
146 }
147 if (fRemoveKinks) {
148 if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
149 return false;
150 }
151 if (fminTPCclsF>track->TPCnclsF())
152 {
153 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
154 return false;
155 }
156
157 // TPC number of clusters:
158 if (fTPCnclsExclusionSwitch) {
159 bool outTPCnclsExclusionZone[2];
160 outTPCnclsExclusionZone[0] = false;
161 outTPCnclsExclusionZone[1] = false;
162 if ( (fminTPCncls > track->TPCncls()) || (fTPCnclsExclusion[0] < track->TPCncls()) ) {
163 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
164 outTPCnclsExclusionZone[0] = true;
165 }
166 if ( (fMaxTPCncls < track->TPCncls()) || (fTPCnclsExclusion[1] > track->TPCncls()) ) {
167 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
168 outTPCnclsExclusionZone[1] = true;
169 }
170 if ( outTPCnclsExclusionZone[0] * outTPCnclsExclusionZone[1] ) { return false; }
171 }
172 else {
173 if (fminTPCncls > track->TPCncls()) {
174 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
175 return false;
176 }
177 if (fMaxTPCncls < track->TPCncls()) {
178 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
179 return false;
180 }
181 }
182
183 if (fminITScls>track->ITSncls())
184 {
185 //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
186 return false;
187 }
188
189 if (fMaxSigmaToVertex < track->SigmaToVertex()) {
190 return false;
191 }
192
193 if (track->ITSncls() > 0)
194 if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
195 return false;
196 }
197
198 // TPC chiNdof of tracks:
199 if (fTPCchiNdofExclusionSwitch && (track->TPCncls() > 0)) {
200 bool outTPCchiNdofExclusionZone[2];
201 outTPCchiNdofExclusionZone[0] = false;
202 outTPCchiNdofExclusionZone[1] = false;
203 if ( (fminTPCchiNdof > (track->TPCchi2()/track->TPCncls())) || (fTPCchiNdofExclusion[0] < (track->TPCchi2()/track->TPCncls())) ) {
204 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
205 outTPCchiNdofExclusionZone[0] = true;
206 }
207 if ( (fMaxTPCchiNdof < (track->TPCchi2()/track->TPCncls())) || (fTPCchiNdofExclusion[1] > (track->TPCchi2()/track->TPCncls())) ) {
208 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
209 outTPCchiNdofExclusionZone[1] = true;
210 }
211 if ( outTPCchiNdofExclusionZone[0] * outTPCchiNdofExclusionZone[1] ) { return false; }
212 }
213 else {
214 if (fminTPCchiNdof > (track->TPCchi2()/track->TPCncls())) {
215 //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
216 return false;
217 }
218 if (fMaxTPCchiNdof < (track->TPCchi2()/track->TPCncls())) {
219 //cout<<" No go because TPC Number of Cls"<<fMaxTPCclsF<< " "<<track->TPCnclsF()<<endl;
220 return false;
221 }
222 }
223
224 if (fLabel)
225 {
226 //cout<<"labels"<<endl;
227 if(track->Label()<0)
228 {
229 fNTracksFailed++;
230 // cout<<"No Go Through the cut"<<endl;
231 // cout<<fLabel<<" Label="<<track->Label()<<endl;
232 return false;
233 }
234 }
235 if (fCharge!=0)
236 {
237 //cout<<"AliFemtoESD cut ch "<<endl;
238 //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
239 if (track->Charge()!= fCharge)
240 {
241 fNTracksFailed++;
242 // cout<<"No Go Through the cut"<<endl;
243 // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
244 return false;
245 }
246 }
247 float tEnergy = ::sqrt(track->P().mag2()+fMass*fMass);
248 float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
249 float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
250 if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
251 {
252 fNTracksFailed++;
253 //cout<<"No Go Through the cut"<<endl;
254 //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
255 return false;
256 }
257 if ((tPt<fPt[0])||(tPt>fPt[1]))
258 {
259 fNTracksFailed++;
260 //cout<<"No Go Through the cut"<<endl;
261 //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
262 return false;
263 }
264// cout << "Track has pids: "
265// << track->PidProbElectron() << " "
266// << track->PidProbMuon() << " "
267// << track->PidProbPion() << " "
268// << track->PidProbKaon() << " "
269// << track->PidProbProton() << " "
270// << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
271
272
273 if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
274 {
275 fNTracksFailed++;
276 //cout<<"No Go Through the cut"<<endl;
277 //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
278 return false;
279 }
280 if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
281 {
282 fNTracksFailed++;
283 //cout<<"No Go Through the cut"<<endl;
284 //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
285 return false;
286 }
287 if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
288 {
289 fNTracksFailed++;
290 //cout<<"No Go Through the cut"<<endl;
291 //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
292 return false;
293 }
294 if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
295 {
296 fNTracksFailed++;
297 //cout<<"No Go Through the cut"<<endl;
298 //cout<<fPidProbProton[0]<<" < p ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
299 return false;
300 }
301 if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
302 {
303 fNTracksFailed++;
304 //cout<<"No Go Through the cut"<<endl;
305 //cout<<fPidProbMuon[0]<<" < mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
306 return false;
307 }
308
309 if (fMostProbable) {
310 tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().mag());
311 tMost[1] = 0.0;
312 tMost[2] = track->PidProbPion()*PidFractionPion(track->P().mag());
313 tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().mag());
314 tMost[4] = track->PidProbProton()*PidFractionProton(track->P().mag());
315 int imost=0;
316 float ipidmax = 0.0;
317 for (int ip=0; ip<5; ip++)
318 if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
319 if (imost != fMostProbable) return false;
320 }
321
322 // cout<<"Go Through the cut"<<endl;
323 // cout<<fLabel<<" Label="<<track->Label()<<endl;
324 // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
325 // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
326 //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
327 //cout<<fPidProbElectron[0]<<" < e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
328 //cout<<fPidProbPion[0]<<" < pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
329 //cout<<fPidProbKaon[0]<<" < k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
330 //cout<<fPidProbProton[0]<<" < p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
331 //cout<<fPidProbMuon[0]<<" < mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
332 fNTracksPassed++ ;
333 return true;
334
335
336}
337//------------------------------
338AliFemtoString AliFemtoQATrackCut::Report()
339{
340 // Prepare report from the execution
341 string tStemp;
342 char tCtemp[100];
343 sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
344 tStemp=tCtemp;
345 sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
346 tStemp+=tCtemp;
347 sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
348 tStemp+=tCtemp;
349 sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
350 tStemp+=tCtemp;
351 sprintf(tCtemp,"Number of tracks which passed:\t%ld Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
352 tStemp += tCtemp;
353 AliFemtoString returnThis = tStemp;
354 return returnThis;
355}
356TList *AliFemtoQATrackCut::ListSettings()
357{
358 // return a list of settings in a writable form
359 TList *tListSetttings = new TList();
360 char buf[200];
361 snprintf(buf, 200, "AliFemtoQATrackCut.mass=%f", this->Mass());
362 tListSetttings->AddLast(new TObjString(buf));
363
364 snprintf(buf, 200, "AliFemtoQATrackCut.charge=%i", fCharge);
365 tListSetttings->AddLast(new TObjString(buf));
366 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobpion.minimum=%f", fPidProbPion[0]);
367 tListSetttings->AddLast(new TObjString(buf));
368 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobpion.maximum=%f", fPidProbPion[1]);
369 tListSetttings->AddLast(new TObjString(buf));
370 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobkaon.minimum=%f", fPidProbKaon[0]);
371 tListSetttings->AddLast(new TObjString(buf));
372 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobkaon.maximum=%f", fPidProbKaon[1]);
373 tListSetttings->AddLast(new TObjString(buf));
374 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobproton.minimum=%f", fPidProbProton[0]);
375 tListSetttings->AddLast(new TObjString(buf));
376 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobproton.maximum=%f", fPidProbProton[1]);
377 tListSetttings->AddLast(new TObjString(buf));
378 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobelectron.minimum=%f", fPidProbElectron[0]);
379 tListSetttings->AddLast(new TObjString(buf));
380 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobelectron.maximum=%f", fPidProbElectron[1]);
381 tListSetttings->AddLast(new TObjString(buf));
382 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobMuon.minimum=%f", fPidProbMuon[0]);
383 tListSetttings->AddLast(new TObjString(buf));
384 snprintf(buf, 200, "AliFemtoQATrackCut.pidprobMuon.maximum=%f", fPidProbMuon[1]);
385 tListSetttings->AddLast(new TObjString(buf));
386 snprintf(buf, 200, "AliFemtoQATrackCut.minimumtpcclusters=%i", fminTPCclsF);
387 tListSetttings->AddLast(new TObjString(buf));
388 snprintf(buf, 200, "AliFemtoQATrackCut.minimumitsclusters=%i", fminTPCclsF);
389 tListSetttings->AddLast(new TObjString(buf));
390 snprintf(buf, 200, "AliFemtoQATrackCut.pt.minimum=%f", fPt[0]);
391 tListSetttings->AddLast(new TObjString(buf));
392 snprintf(buf, 200, "AliFemtoQATrackCut.pt.maximum=%f", fPt[1]);
393 tListSetttings->AddLast(new TObjString(buf));
394 snprintf(buf, 200, "AliFemtoQATrackCut.rapidity.minimum=%f", fRapidity[0]);
395 tListSetttings->AddLast(new TObjString(buf));
396 snprintf(buf, 200, "AliFemtoQATrackCut.rapidity.maximum=%f", fRapidity[1]);
397 tListSetttings->AddLast(new TObjString(buf));
398 snprintf(buf, 200, "AliFemtoQATrackCut.removekinks=%i", fRemoveKinks);
399 tListSetttings->AddLast(new TObjString(buf));
400 snprintf(buf, 200, "AliFemtoQATrackCut.maxitschindof=%f", fMaxITSchiNdof);
401 tListSetttings->AddLast(new TObjString(buf));
402 snprintf(buf, 200, "AliFemtoQATrackCut.maxtpcchindof=%f", fMaxTPCchiNdof);
403 tListSetttings->AddLast(new TObjString(buf));
404 snprintf(buf, 200, "AliFemtoQATrackCut.maxsigmatovertex=%f", fMaxSigmaToVertex);
405 tListSetttings->AddLast(new TObjString(buf));
406 if (fMostProbable) {
407 if (fMostProbable == 2)
408 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Pion");
409 if (fMostProbable == 3)
410 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Kaon");
411 if (fMostProbable == 4)
412 snprintf(buf, 200, "AliFemtoQATrackCut.mostprobable=%s", "Proton");
413 tListSetttings->AddLast(new TObjString(buf));
414 }
415 return tListSetttings;
416}
417void AliFemtoQATrackCut::SetRemoveKinks(const bool& flag)
418{
419 fRemoveKinks = flag;
420}
421
422 // electron
423// 0.13 - 1.8
424// 0 7.594129e-02 8.256141e-03
425// 1 -5.535827e-01 8.170825e-02
426// 2 1.728591e+00 3.104210e-01
427// 3 -2.827893e+00 5.827802e-01
428// 4 2.503553e+00 5.736207e-01
429// 5 -1.125965e+00 2.821170e-01
430// 6 2.009036e-01 5.438876e-02
431float AliFemtoQATrackCut::PidFractionElectron(float mom) const
432{
433 // Provide a parameterized fraction of electrons dependent on momentum
434 if (mom<0.13) return 0.0;
435 if (mom>1.8) return 0.0;
436 return (7.594129e-02
437 -5.535827e-01*mom
438 +1.728591e+00*mom*mom
439 -2.827893e+00*mom*mom*mom
440 +2.503553e+00*mom*mom*mom*mom
441 -1.125965e+00*mom*mom*mom*mom*mom
442 +2.009036e-01*mom*mom*mom*mom*mom*mom);
443}
444
445// pion
446// 0.13 - 2.0
447// 0 1.063457e+00 8.872043e-03
448// 1 -4.222208e-01 2.534402e-02
449// 2 1.042004e-01 1.503945e-02
450float AliFemtoQATrackCut::PidFractionPion(float mom) const
451{
452 // Provide a parameterized fraction of pions dependent on momentum
453 if (mom<0.13) return 0.0;
454 if (mom>2.0) return 0.0;
455 return ( 1.063457e+00
456 -4.222208e-01*mom
457 +1.042004e-01*mom*mom);
458}
459
460// kaon
461// 0.18 - 2.0
462// 0 -7.289406e-02 1.686074e-03
463// 1 4.415666e-01 1.143939e-02
464// 2 -2.996790e-01 1.840964e-02
465// 3 6.704652e-02 7.783990e-03
466float AliFemtoQATrackCut::PidFractionKaon(float mom) const
467{
468 // Provide a parameterized fraction of kaons dependent on momentum
469 if (mom<0.18) return 0.0;
470 if (mom>2.0) return 0.0;
471 return (-7.289406e-02
472 +4.415666e-01*mom
473 -2.996790e-01*mom*mom
474 +6.704652e-02*mom*mom*mom);
475}
476
477// proton
478// 0.26 - 2.0
479// 0 -3.730200e-02 2.347311e-03
480// 1 1.163684e-01 1.319316e-02
481// 2 8.354116e-02 1.997948e-02
482// 3 -4.608098e-02 8.336400e-03
483float AliFemtoQATrackCut::PidFractionProton(float mom) const
484{
485 // Provide a parameterized fraction of protons dependent on momentum
486 if (mom<0.26) return 0.0;
487 if (mom>2.0) return 0.0;
488 return (-3.730200e-02
489 +1.163684e-01*mom
490 +8.354116e-02*mom*mom
491 -4.608098e-02*mom*mom*mom);
492}