]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFemtoParticle.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Pad size less then cell size + ideal geom in v2
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Infrastructure / AliFemtoParticle.cxx
CommitLineData
67427ff7 1/***************************************************************************
2 *
3 * $Id$
4 *
5 * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6 ***************************************************************************
7 *
8 * Description: part of STAR HBT Framework: AliFemtoMaker package
9 * Particle objects are part of the PicoEvent, which is what is
10 * stored in the EventMixingBuffers
11 * A Track object gets converted to a Particle object if it
12 * passes the ParticleCut of an Analysis
13 *
14 ***************************************************************************
15 *
16 * $Log$
0215f606 17 * Revision 1.1.1.1 2007/04/25 15:38:41 panos
18 * Importing the HBT code dir
19 *
67427ff7 20 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
21 * First version on CVS
22 *
23 * Revision 1.22 2003/09/02 17:58:32 perev
24 * gcc 3.2 updates + WarnOff
25 *
26 * Revision 1.21 2003/05/07 15:30:43 magestro
27 * Fixed bug related to finding merged hits (commit for Fabrice)
28 *
29 * Revision 1.20 2003/01/14 09:41:26 renault
30 * changes on average separation calculation, hit shared finder and memory optimisation
31 * for Z,U and Sectors variables.
32 *
33 * Revision 1.19 2002/12/12 17:01:49 kisiel
34 * Hidden Information handling and purity calculation
35 *
36 * Revision 1.18 2002/11/19 23:36:00 renault
37 * Enable calculation of exit/entrance separation for V0 daughters
38 *
39 * Revision 1.17 2001/12/14 23:11:30 fretiere
40 * Add class HitMergingCut. Add class fabricesPairCut = HitMerginCut + pair purity cuts. Add TpcLocalTransform function which convert to local tpc coord (not pretty). Modify AliFemtoTrack, AliFemtoParticle, AliFemtoHiddenInfo, AliFemtoPair to handle the hit information and cope with my code
41 *
42 * Revision 1.16 2001/05/25 23:23:59 lisa
43 * Added in AliFemtoKink stuff
44 *
45 * Revision 1.15 2001/04/03 21:04:36 kisiel
46 *
47 *
48 * Changes needed to make the Theoretical code
49 * work. The main code is the ThCorrFctn directory.
50 * The most visible change is the addition of the
51 * HiddenInfo to AliFemtoPair.
52 *
53 * Revision 1.14 2000/10/05 23:09:05 lisa
54 * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
55 *
56 * Revision 1.13 2000/08/31 22:31:31 laue
57 * AliFemtoAnalysis: output changed (a little bit less)
58 * AliFemtoEvent: new version, members for reference mult added
59 * AliFemtoIOBinary: new IO for new AliFemtoEvent version
60 * AliFemtoTypes: TTree typedef to AliFemtoTTree added
61 * AliFemtoVertexAnalysis: overflow and underflow added
62 *
63 * Revision 1.12 2000/07/23 13:52:56 laue
64 * NominalExitPoint set to (-9999.,-9999.-9999.) if helix.at()
65 * returns nan (not a number).
66 *
67 * Revision 1.11 2000/07/19 17:18:48 laue
68 * Calculation of Entrance and Exit point added in AliFemtoParticle constructor
69 *
70 * Revision 1.10 2000/07/17 20:03:17 lisa
71 * Implemented tools for addressing and assessing trackmerging
72 *
73 * Revision 1.9 2000/07/16 21:38:23 laue
74 * AliFemtoCoulomb.cxx AliFemtoSectoredAnalysis.cxx : updated for standalone version
75 * AliFemtoV0.cc AliFemtoV0.h : some cast to prevent compiling warnings
76 * AliFemtoParticle.cc AliFemtoParticle.h : pointers fTrack,fV0 initialized to 0
77 * AliFemtoIOBinary.cc : some printouts in #ifdef STHBTDEBUG
78 * AliFemtoEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
79 * solution
80 *
81 * Revision 1.8 2000/05/03 17:44:43 laue
82 * AliFemtoEvent, AliFemtoTrack & AliFemtoV0 declared friend to AliFemtoIOBinary
83 * AliFemtoParticle updated for V0 pos,neg track Id
84 *
85 * Revision 1.7 2000/04/03 16:21:51 laue
86 * some include files changed
87 * Multi track cut added
88 *
89 * Revision 1.6 1999/12/11 15:58:29 lisa
90 * Add vertex decay position datum and accessor to AliFemtoParticle to allow pairwise cuts on seperation of V0s
91 *
92 * Revision 1.5 1999/09/17 22:38:02 lisa
93 * first full integration of V0s into AliFemto framework
94 *
95 * Revision 1.4 1999/09/01 19:04:53 lisa
96 * update Particle class AND add parity cf and Randys Coulomb correction
97 *
98 * Revision 1.3 1999/07/06 22:33:23 lisa
99 * Adjusted all to work in pro and new - dev itself is broken
100 *
101 * Revision 1.2 1999/06/29 17:50:27 fisyak
102 * formal changes to account new StEvent, does not complie yet
103 *
104 * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
105 * Installation of AliFemtoMaker
106 *
107 **************************************************************************/
108
109#include "Infrastructure/AliFemtoParticle.h"
110//#include "math_constants.h"
111#ifdef __CC5__
112 #include <math.h>
113#else
114 #include <cmath>
115#endif
116
117
118#include "TMath.h"
119using namespace TMath;
120
121double AliFemtoParticle::fPrimPimPar0= 9.05632e-01;
122double AliFemtoParticle::fPrimPimPar1= -2.26737e-01;
123double AliFemtoParticle::fPrimPimPar2= -1.03922e-01;
124double AliFemtoParticle::fPrimPipPar0= 9.09616e-01;
125double AliFemtoParticle::fPrimPipPar1= -9.00511e-02;
126double AliFemtoParticle::fPrimPipPar2= -6.02940e-02;
127double AliFemtoParticle::fPrimPmPar0= 0.;
128double AliFemtoParticle::fPrimPmPar1= 0.;
129double AliFemtoParticle::fPrimPmPar2= 0.;
130double AliFemtoParticle::fPrimPpPar0= 0.;
131double AliFemtoParticle::fPrimPpPar1= 0.;
132double AliFemtoParticle::fPrimPpPar2= 0.;
133
134int TpcLocalTransform(AliFmThreeVectorD& xgl,
135 int& iSector,
136 int& iPadrow,
137 float& xlocal,
138 double& ttPhi);
139
140
141//_____________________
0215f606 142AliFemtoParticle::AliFemtoParticle() :
143 fTpcV0NegPosSample(0),
144 fV0NegZ(0),
145 fV0NegU(0),
146 fV0NegSect(0),
147 fTrack(0), fV0(0), fKink(0), fXi(0),
148 fFourMomentum(0),
149 fHelix(),
150 fNominalTpcExitPoint(0),
151 fNominalTpcEntrancePoint(0),
152 fHiddenInfo(0),
153 fPrimaryVertex(0),
154 fSecondaryVertex(0),
155 fHelixV0Pos(),
156 fTpcV0PosEntrancePoint(0),
157 fTpcV0PosExitPoint(0),
158 fHelixV0Neg(),
159 fTpcV0NegEntrancePoint(0),
160 fTpcV0NegExitPoint(0)
161{
67427ff7 162 /* no-op for default */
163 // cout << "Created particle " << this << endl;
164}
165//_____________________
0215f606 166AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
167 fTpcV0NegPosSample(0),
168 fV0NegZ(0),
169 fV0NegU(0),
170 fV0NegSect(0),
171 fTrack(0), fV0(0), fKink(0), fXi(0),
172 fFourMomentum(0),
173 fHelix(),
174 fNominalTpcExitPoint(0),
175 fNominalTpcEntrancePoint(0),
176 fHiddenInfo(0),
177 fPrimaryVertex(0),
178 fSecondaryVertex(0),
179 fHelixV0Pos(),
180 fTpcV0PosEntrancePoint(0),
181 fTpcV0PosExitPoint(0),
182 fHelixV0Neg(),
183 fTpcV0NegEntrancePoint(0),
184 fTpcV0NegExitPoint(0)
185{
186 if (aParticle.fTrack)
187 fTrack = new AliFemtoTrack(*aParticle.fTrack);
188 if (aParticle.fV0)
189 fV0 = new AliFemtoV0(*aParticle.fV0);
190 if (aParticle.fKink)
191 fKink = new AliFemtoKink(*aParticle.fKink);
192 if (aParticle.fXi)
193 fXi = new AliFemtoXi(*aParticle.fXi);
194 fFourMomentum = aParticle.fFourMomentum;
195 fHelix = aParticle.fHelix;
196
197 for (int iter=0; iter<11; iter++)
198 fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
199
200 if (aParticle.fTpcV0NegPosSample) {
201 fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
202 for (int iter=0; iter<11; iter++)
203 fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
204 }
205
206 if (aParticle.fV0NegZ) {
207 fV0NegZ = (float *) malloc(sizeof(float) * 45);
208 for (int iter=0; iter<11; iter++)
209 fV0NegZ[iter] = aParticle.fV0NegZ[iter];
210 }
211 if (aParticle.fV0NegU) {
212 fV0NegU = (float *) malloc(sizeof(float) * 45);
213 for (int iter=0; iter<11; iter++)
214 fV0NegU[iter] = aParticle.fV0NegU[iter];
215 }
216 if (aParticle.fV0NegSect) {
217 fV0NegSect = (int *) malloc(sizeof(int) * 45);
218 for (int iter=0; iter<11; iter++)
219 fV0NegSect[iter] = aParticle.fV0NegSect[iter];
220 }
221
222 fPrimaryVertex = aParticle.fPrimaryVertex;
223 fSecondaryVertex = aParticle.fSecondaryVertex;
224 CalculatePurity();
225 if(aParticle.fHiddenInfo){
226 fHiddenInfo= aParticle.HiddenInfo()->clone();
227 }
228
229 fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
230 fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
231
232 for (int iter=0; iter<6; iter++)
233 fPurity[iter] = aParticle.fPurity[iter];
234
235 fHelixV0Pos = aParticle.fHelixV0Pos;
236 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
237 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
238 fHelixV0Neg = aParticle.fHelixV0Neg;
239 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
240 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
241}
242//_____________________
67427ff7 243AliFemtoParticle::~AliFemtoParticle(){
244 // cout << "Issuing delete for AliFemtoParticle." << endl;
245
246 if (fTrack) delete fTrack;
247 if (fV0) {
248 delete[] fTpcV0NegPosSample;
249 delete[] fV0NegZ;
250 delete[] fV0NegU;
251 delete[] fV0NegSect;
252 delete fV0;
253 }
254 if (fKink) delete fKink;
255 // cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
256 if (fHiddenInfo)
257 {
258 // cout << "Deleting HiddenInfo." << endl;
259 delete fHiddenInfo;
260 }
261 // cout << "Deleted particle " << this << endl;
262}
263//_____________________
0215f606 264AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) :
265 fTpcV0NegPosSample(0),
266 fV0NegZ(0),
267 fV0NegU(0),
268 fV0NegSect(0),
269 fTrack(0), fV0(0), fKink(0), fXi(0),
270 fFourMomentum(0),
271 fHelix(),
272 fNominalTpcExitPoint(0),
273 fNominalTpcEntrancePoint(0),
274 fHiddenInfo(0),
275 fPrimaryVertex(0),
276 fSecondaryVertex(0),
277 fHelixV0Pos(),
278 fTpcV0PosEntrancePoint(0),
279 fTpcV0PosExitPoint(0),
280 fHelixV0Neg(),
281 fTpcV0NegEntrancePoint(0),
282 fTpcV0NegExitPoint(0)
283{
67427ff7 284
285
286 // I know there is a better way to do this...
287 fTrack = new AliFemtoTrack(*hbtTrack);
288 AliFemtoThreeVector temp = hbtTrack->P();
289 fFourMomentum.setVect(temp);
290 double ener = ::sqrt(temp.mag2()+mass*mass);
291 fFourMomentum.setE(ener);
292// fMap[0] = hbtTrack->TopologyMap(0);
293 // fMap[1] = hbtTrack->TopologyMap(1);
294 // fNhits = hbtTrack->NHits();
295 fHelix = hbtTrack->Helix();
296 //CalculateNominalTpcExitAndEntrancePoints();
297
298
299 fPrimaryVertex.setX(0.);
300 fPrimaryVertex.setY(0.);
301 fPrimaryVertex.setZ(0.);
302 fSecondaryVertex.setX(0.);
303 fSecondaryVertex.setY(0.);
304 fSecondaryVertex.setZ(0.);
305 /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
306 CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
307 &fSecondaryVertex,
308 &fNominalTpcEntrancePoint,
309 &fNominalTpcExitPoint,
310 &mNominalPosSample[0],
311 &fZ[0],
312 &fU[0],
313 &fSect[0]);
314 */
315 CalculatePurity();
316 // ***
317 fHiddenInfo= 0;
318 if(hbtTrack->ValidHiddenInfo()){
319 fHiddenInfo= hbtTrack->getHiddenInfo()->getParticleHiddenInfo()->clone();
320 }
321 // ***
322 // cout << "Created particle " << this << endl;
323
324}
325//_____________________
0215f606 326AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) :
327 fTpcV0NegPosSample(0),
328 fV0NegZ(0),
329 fV0NegU(0),
330 fV0NegSect(0),
331 fTrack(0), fV0(0), fKink(0), fXi(0),
332 fFourMomentum(0),
333 fHelix(),
334 fNominalTpcExitPoint(0),
335 fNominalTpcEntrancePoint(0),
336 fHiddenInfo(0),
337 fPrimaryVertex(0),
338 fSecondaryVertex(0),
339 fHelixV0Pos(),
340 fTpcV0PosEntrancePoint(0),
341 fTpcV0PosExitPoint(0),
342 fHelixV0Neg(),
343 fTpcV0NegEntrancePoint(0),
344 fTpcV0NegExitPoint(0)
345{
67427ff7 346 fV0 = new AliFemtoV0(*hbtV0);
347 //fMap[0]= 0;
348 //fMap[1]= 0;
349 // I know there is a better way to do this...
350 AliFemtoThreeVector temp = hbtV0->momV0();
351 fFourMomentum.setVect(temp);
352 double ener = ::sqrt(temp.mag2()+mass*mass);
353 fFourMomentum.setE(ener);
354 // Calculating TpcEntrancePoint for Positive V0 daugther
355 fPrimaryVertex = hbtV0->primaryVertex();
356 fSecondaryVertex = hbtV0->decayVertexV0();
357 fHelixV0Pos = hbtV0->HelixPos();
358
359 fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
360 fV0NegZ = new float[45];//for V0Neg
361 fV0NegU = new float[45];//for V0Neg
362 fV0NegSect = new int[45];//for V0Neg
363 CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
364 &fSecondaryVertex,
365 &fTpcV0PosEntrancePoint,
366 &fTpcV0PosExitPoint,
367 &fNominalPosSample[0],
368 &fZ[0],
369 &fU[0],&fSect[0]);
370 fHelixV0Neg = hbtV0->HelixNeg();
371
372 CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
373 &fPrimaryVertex,
374 &fSecondaryVertex,
375 &fTpcV0NegEntrancePoint,
376 &fTpcV0NegExitPoint,
377 &fTpcV0NegPosSample[0],
378 &fV0NegZ[0],
379 &fV0NegU[0],&fV0NegSect[0]);
380
381 // ***
382 fHiddenInfo= 0;
383 if(hbtV0->ValidHiddenInfo()){
384 fHiddenInfo= hbtV0->getHiddenInfo()->clone();
385 }
386 // ***
387}
388//_____________________
0215f606 389AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) :
390 fTpcV0NegPosSample(0),
391 fV0NegZ(0),
392 fV0NegU(0),
393 fV0NegSect(0),
394 fTrack(0), fV0(0), fKink(0), fXi(0),
395 fFourMomentum(0),
396 fHelix(),
397 fNominalTpcExitPoint(0),
398 fNominalTpcEntrancePoint(0),
399 fHiddenInfo(0),
400 fPrimaryVertex(0),
401 fSecondaryVertex(0),
402 fHelixV0Pos(),
403 fTpcV0PosEntrancePoint(0),
404 fTpcV0PosExitPoint(0),
405 fHelixV0Neg(),
406 fTpcV0NegEntrancePoint(0),
407 fTpcV0NegExitPoint(0)
408{
67427ff7 409 fKink = new AliFemtoKink(*hbtKink);
410 // fMap[0]= 0;
411 //fMap[1]= 0;
412 // I know there is a better way to do this...
413 AliFemtoThreeVector temp = hbtKink->Parent().P();
414 fFourMomentum.setVect(temp);
415 double ener = ::sqrt(temp.mag2()+mass*mass);
416 fFourMomentum.setE(ener);
417}
418
419//_____________________
0215f606 420AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
421 fTpcV0NegPosSample(0),
422 fV0NegZ(0),
423 fV0NegU(0),
424 fV0NegSect(0),
425 fTrack(0), fV0(0), fKink(0), fXi(0),
426 fFourMomentum(0),
427 fHelix(),
428 fNominalTpcExitPoint(0),
429 fNominalTpcEntrancePoint(0),
430 fHiddenInfo(0),
431 fPrimaryVertex(0),
432 fSecondaryVertex(0),
433 fHelixV0Pos(),
434 fTpcV0PosEntrancePoint(0),
435 fTpcV0PosExitPoint(0),
436 fHelixV0Neg(),
437 fTpcV0NegEntrancePoint(0),
438 fTpcV0NegExitPoint(0)
439{
67427ff7 440 fXi = new AliFemtoXi(*hbtXi);
441 // fMap[0]= 0;
442 //fMap[1]= 0;
443 AliFemtoThreeVector temp;// = hbtXi->mMofXi;
444 fFourMomentum.setVect(temp);
445 double ener = ::sqrt(temp.mag2()+mass*mass);
446 fFourMomentum.setE(ener);
447 fHiddenInfo = 0;
448}
449//_____________________
0215f606 450AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
451{
452 if (this == &aParticle)
453 return *this;
454
455 if (aParticle.fTrack)
456 fTrack = new AliFemtoTrack(*aParticle.fTrack);
457 if (aParticle.fV0)
458 fV0 = new AliFemtoV0(*aParticle.fV0);
459 if (aParticle.fKink)
460 fKink = new AliFemtoKink(*aParticle.fKink);
461 if (aParticle.fXi)
462 fXi = new AliFemtoXi(*aParticle.fXi);
463 fFourMomentum = aParticle.fFourMomentum;
464 fHelix = aParticle.fHelix;
465
466 for (int iter=0; iter<11; iter++)
467 fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
468
469 if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
470 if (aParticle.fTpcV0NegPosSample) {
471 fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
472 for (int iter=0; iter<11; iter++)
473 fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
474 }
475
476 if (fV0NegZ) delete fV0NegZ;
477 if (aParticle.fV0NegZ) {
478 fV0NegZ = (float *) malloc(sizeof(float) * 45);
479 for (int iter=0; iter<11; iter++)
480 fV0NegZ[iter] = aParticle.fV0NegZ[iter];
481 }
482 if (fV0NegU) delete fV0NegU;
483 if (aParticle.fV0NegU) {
484 fV0NegU = (float *) malloc(sizeof(float) * 45);
485 for (int iter=0; iter<11; iter++)
486 fV0NegU[iter] = aParticle.fV0NegU[iter];
487 }
488 if (fV0NegSect) delete fV0NegSect;
489 if (aParticle.fV0NegSect) {
490 fV0NegSect = (int *) malloc(sizeof(int) * 45);
491 for (int iter=0; iter<11; iter++)
492 fV0NegSect[iter] = aParticle.fV0NegSect[iter];
493 }
494
495 fPrimaryVertex = aParticle.fPrimaryVertex;
496 fSecondaryVertex = aParticle.fSecondaryVertex;
497 CalculatePurity();
498 if(aParticle.fHiddenInfo){
499 fHiddenInfo= aParticle.fHiddenInfo->clone();
500 }
501
502 fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
503 fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
504
505 if (fHiddenInfo) delete fHiddenInfo;
506 if (aParticle.fHiddenInfo)
507 fHiddenInfo = aParticle.HiddenInfo()->clone();
508
509 for (int iter=0; iter<6; iter++)
510 fPurity[iter] = aParticle.fPurity[iter];
511
512 fHelixV0Pos = aParticle.fHelixV0Pos;
513 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
514 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
515 fHelixV0Neg = aParticle.fHelixV0Neg;
516 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
517 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
518
519 return *this;
520}
521//_____________________
67427ff7 522const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
523 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
524 // for now, we calculate Exit and Entrance points upon instantiation
525 return fNominalTpcExitPoint;
526}
527//_____________________
528const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
529 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
530 // for now, we calculate Exit and Entrance points upon instantiation
531 return fNominalTpcEntrancePoint;
532}
533//_____________________
534void AliFemtoParticle::CalculatePurity(){
535 double tPt = fFourMomentum.perp();
536 // pi -
537 fPurity[0] = fPrimPimPar0*(1.-exp((tPt-fPrimPimPar1)/fPrimPimPar2));
538 fPurity[0] *= fTrack->PidProbPion();
539 // pi+
540 fPurity[1] = fPrimPipPar0*(1.-exp((tPt-fPrimPipPar1)/fPrimPipPar2));
541 fPurity[1] *= fTrack->PidProbPion();
542 // K-
543 fPurity[2] = fTrack->PidProbKaon();
544 // K+
545 fPurity[3] = fTrack->PidProbKaon();
546 // pbar
547 fPurity[4] = fTrack->PidProbProton();
548 // p
549 fPurity[5] = fTrack->PidProbProton();
550}
551
552double AliFemtoParticle::GetPionPurity()
553{
554 if (fTrack->Charge()>0)
555 return fPurity[1];
556 else
557 return fPurity[0];
558}
559double AliFemtoParticle::GetKaonPurity()
560{
561 if (fTrack->Charge()>0)
562 return fPurity[3];
563 else
564 return fPurity[2];
565}
566double AliFemtoParticle::GetProtonPurity()
567{
568 if (fTrack->Charge()>0)
569 return fPurity[5];
570 else
571 return fPurity[4];
572}
573
574void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
575 AliFemtoThreeVector* PrimVert,
576 AliFemtoThreeVector* SecVert,
577 AliFemtoThreeVector* tmpTpcEntrancePoint,
578 AliFemtoThreeVector* tmpTpcExitPoint,
579 AliFemtoThreeVector* tmpPosSample,
580 float* tmpZ,
581 float* tmpU,
582 int* tmpSect){
583 // this calculates the exit point of a secondary track,
584 // either through the endcap or through the Outer Field Cage
585 // We assume the track to start at tHelix.origin-PrimaryVertex
586 // it also calculates the entrance point of the secondary track,
587 // which is the point at which it crosses the
588 // inner field cage
589 // static AliFemtoThreeVector ZeroVec(0.,0.,0.);
590 AliFemtoThreeVector ZeroVec(0.,0.,0.);
591// ZeroVec.setX(tHelix->origin().x()-PrimVert->x());
592// ZeroVec.setY(tHelix->origin().y()-PrimVert->y());
593// ZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
594 ZeroVec.setX(SecVert->x()-PrimVert->x());
595 ZeroVec.setY(SecVert->y()-PrimVert->y());
596 ZeroVec.setZ(SecVert->z()-PrimVert->z());
597 double dip, curv, phase;
598 int h;
599 curv = tHelix->curvature();
600 dip = tHelix->dipAngle();
601 phase= tHelix->phase();
602 h = tHelix->h();
603
604 AliFmHelixD hel(curv,dip,phase,ZeroVec,h);
605
606 pairD candidates;
607 double sideLength; // this is how much length to go to leave through sides of TPC
608 double endLength; // this is how much length to go to leave through endcap of TPC
609 // figure out how far to go to leave through side...
610 candidates = hel.pathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
611 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
612
613 static AliFemtoThreeVector WestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
614 static AliFemtoThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
615 static AliFemtoThreeVector EndCapNormal(0.,0.,1.0);
616
617 endLength = hel.pathLength(WestEnd,EndCapNormal);
618 if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
619
620 if (endLength < 0.0) cout <<
621 "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
622 << "Hey -- I cannot find an exit point out endcaps" << endl;
623 // OK, firstExitLength will be the shortest way out of the detector...
624 double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
625 // now then, let's return the POSITION at which particle leaves TPC...
626 *tmpTpcExitPoint = hel.at(firstExitLength);
627 // Finally, calculate the position at which the track crosses the inner field cage
628 candidates = hel.pathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
629
630 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
631 // cout << "sideLength 2 ="<<sideLength << endl;
632 *tmpTpcEntrancePoint = hel.at(sideLength);
633 // This is the secure way !
634 if (IsNaN(tmpTpcEntrancePoint->x()) ||
635 IsNaN(tmpTpcEntrancePoint->y()) ||
636 IsNaN(tmpTpcEntrancePoint->z()) ){
637 cout << "tmpTpcEntrancePoint NAN"<< endl;
638 cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
639 tmpTpcEntrancePoint->setX(-9999.);
640 tmpTpcEntrancePoint->setY(-9999.);
641 tmpTpcEntrancePoint->setZ(-9999.);
642 }
643
644 if (IsNaN(tmpTpcExitPoint->x()) ||
645 IsNaN(tmpTpcExitPoint->y()) ||
646 IsNaN(tmpTpcExitPoint->z()) ) {
647// cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
648// cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
649// cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
650// cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
651 tmpTpcExitPoint->setX(-9999.);
652 tmpTpcExitPoint->setY(-9999.);
653 tmpTpcExitPoint->setZ(-9999.);
654 }
655
656
657// if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
658// if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
659 // cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
660
661 // 03Oct00 - mal. OK, let's try something a little more
662 // along the lines of NA49 and E895 strategy.
663 // calculate the "nominal" position at N radii (say N=11)
664 // within the TPC, and for a pair cut
665 // use the average separation of these N
666 int irad = 0;
667 candidates = hel.pathLength(50.0);
668 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
669 while (irad<11 && !IsNaN(sideLength)){
670 float radius = 50.0 + irad*15.0;
671 candidates = hel.pathLength(radius);
672 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
673 tmpPosSample[irad] = hel.at(sideLength);
674 if(IsNaN(tmpPosSample[irad].x()) ||
675 IsNaN(tmpPosSample[irad].y()) ||
676 IsNaN(tmpPosSample[irad].z())
677 ){
678 cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
679 cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
680 tmpPosSample[irad] = AliFemtoThreeVector(-9999.,-9999.,-9999);
681 }
682 irad++;
683 if (irad<11){
684 float radius = 50.0 + irad*15.0;
685 candidates = hel.pathLength(radius);
686 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
687 }
688 }
689 for (int i = irad; i<11; i++)
690 {
691 tmpPosSample[i] = AliFemtoThreeVector(-9999.,-9999.,-9999);
692 }
693
694 static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
695 104,109.2,114.4,119.6,127.195,129.195,131.195,
696 133.195,135.195,137.195,139.195,141.195,
697 143.195,145.195,147.195,149.195,151.195,
698 153.195,155.195,157.195,159.195,161.195,
699 163.195,165.195,167.195,169.195,171.195,
700 173.195,175.195,177.195,179.195,181.195,
701 183.195,185.195,187.195,189.195};
702 int tRow,tSect,tOutOfBound;
703 double tLength,tPhi;
704 float tU;
705 AliFemtoThreeVector tPoint;
706 AliFmThreeVectorD tn(0,0,0);
707 AliFmThreeVectorD tr(0,0,0);
708 int ti =0;
709 // test to enter the loop
710 candidates = hel.pathLength(tRowRadius[ti]);
711 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
712 if (IsNaN(tLength)){
713 cout <<"tLength Init tmp NAN" << endl;
714 cout <<"padrow number= "<<ti << "not reached" << endl;
715 cout << "*** DO NOT ENTER THE LOOP***" << endl;
716 tmpSect[ti]=-1;//sector
717 }
718 // end test
719 while(ti<45 && !IsNaN(tLength)){
720 candidates = hel.pathLength(tRowRadius[ti]);
721 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
722 if (IsNaN(tLength)){
723 cout <<"tLength loop 1st NAN" << endl;
724 cout <<"padrow number= " << ti << " not reached" << endl;
725 cout << "*** THIS IS AN ERROR SHOULDN'T LOOP ***" << endl;
726 tmpSect[ti]=-1;//sector
727 }
728 tPoint = hel.at(tLength);
729 // Find which sector it is on
730 TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
731 if (IsNaN(tmpSect[ti])){
732 cout <<"***ERROR tmpSect"<< endl;
733 }
734 if (IsNaN(tRow)){
735 cout <<"***ERROR tRow"<< endl;
736 }
737 if (IsNaN(tU)){
738 cout <<"***ERROR tU"<< endl;
739 }
740 if (IsNaN(tPhi)){
741 cout <<"***ERROR tPhi"<< endl;
742 }
743 // calculate crossing plane
744 tn.setX(cos(tPhi));
745 tn.setY(sin(tPhi));
746 tr.setX(tRowRadius[ti]*cos(tPhi));
747 tr.setY(tRowRadius[ti]*sin(tPhi));
748 // find crossing point
749 tLength = hel.pathLength(tr,tn);
750 if (IsNaN(tLength)){
751 cout <<"tLength loop 2nd NAN" << endl;
752 cout <<"padrow number= " << ti << " not reached" << endl;
753 tmpSect[ti]=-2;//sector
754 }
755 tPoint = hel.at(tLength);
756 tmpZ[ti] = tPoint.z();
757 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
758 if (IsNaN(tSect)){
759 cout <<"***ERROR tSect 2"<< endl;
760 }
761 if (IsNaN(tRow)){
762 cout <<"***ERROR tRow 2"<< endl;
763 }
764 if (IsNaN(tmpU[ti])){
765 cout <<"***ERROR tmpU[ti] 2"<< endl;
766 }
767 if (IsNaN(tPhi)){
768 cout <<"***ERROR tPhi 2 "<< endl;
769 }
770 if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
771 tmpSect[ti]=-2;
772 // cout << "missed once"<< endl;
773 }
774 else{
775 if(tmpSect[ti] != tSect){
776 // Try again on the other sector
777 tn.setX(cos(tPhi));
778 tn.setY(sin(tPhi));
779 tr.setX(tRowRadius[ti]*cos(tPhi));
780 tr.setY(tRowRadius[ti]*sin(tPhi));
781 // find crossing point
782 tLength = hel.pathLength(tr,tn);
783 tPoint = hel.at(tLength);
784 if (IsNaN(tLength)){
785 cout <<"tLength loop 3rd NAN" << endl;
786 cout <<"padrow number= "<< ti << " not reached" << endl;
787 tmpSect[ti]=-1;//sector
788 }
789 tmpZ[ti] = tPoint.z();
790 tmpSect[ti] = tSect;
791 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
792 if (IsNaN(tSect)){
793 cout <<"***ERROR tSect 3"<< endl;
794 }
795 if (IsNaN(tRow)){
796 cout <<"***ERROR tRow 3"<< endl;
797 }
798 if (IsNaN(tmpU[ti])){
799 cout <<"***ERROR tmpU[ti] 3"<< endl;
800 }
801 if (IsNaN(tPhi)){
802 cout <<"***ERROR tPhi 3 "<< endl;
803 }
804 if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
805 tmpSect[ti]=-1;
806 }
807 }
808 }
809 if (IsNaN(tmpSect[ti])){
810 cout << "*******************ERROR***************************" << endl;
811 cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
812 cout << "*******************ERROR***************************" << endl;
813 }
814 if (IsNaN(tmpU[ti])){
815 cout << "*******************ERROR***************************" << endl;
816 cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
817 cout << "*******************ERROR***************************" << endl;
818 }
819 if (IsNaN(tmpZ[ti])){
820 cout << "*******************ERROR***************************" << endl;
821 cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
822 cout << "*******************ERROR***************************" << endl;
823 }
824 // If padrow ti not reached all other beyond are not reached
825 // in this case set sector to -1
826 if (tmpSect[ti]==-1){
827 for (int tj=ti; tj<45;tj++){
828 tmpSect[tj] = -1;
829 ti=45;
830 }
831 }
832 ti++;
833 if (ti<45){
834 candidates = hel.pathLength(tRowRadius[ti]);
835 tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
836 }
837}
838//_____________________
839const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
840 return fTpcV0PosExitPoint;
841}
842//_____________________
843const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
844 return fTpcV0PosEntrancePoint;
845}
846//______________________
847const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
848 return fTpcV0NegExitPoint;
849}
850//_____________________
851const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
852 return fTpcV0NegEntrancePoint;
853}
854//______________________