Adding static_cast to keep the compiler happy
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoParticle.cxx
CommitLineData
d0e92d9a 1///////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoParticle: main class halding all the necessary information //
4// about particle that is required during femtoscopic analysis //
5// This includes all the information about the quality of the track, //
6// its identification as well as track chracteristics with connection //
7// to the detector parts, e.g. entrance and exit points. //
8// //
9///////////////////////////////////////////////////////////////////////////
e5d53b83 10#include "AliFemtoKink.h"
d0e92d9a 11#include "AliFemtoParticle.h"
e5d53b83 12#include "AliFemtoXi.h"
67427ff7 13
d0e92d9a 14double AliFemtoParticle::fgPrimPimPar0= 9.05632e-01;
15double AliFemtoParticle::fgPrimPimPar1= -2.26737e-01;
16double AliFemtoParticle::fgPrimPimPar2= -1.03922e-01;
17double AliFemtoParticle::fgPrimPipPar0= 9.09616e-01;
18double AliFemtoParticle::fgPrimPipPar1= -9.00511e-02;
19double AliFemtoParticle::fgPrimPipPar2= -6.02940e-02;
20double AliFemtoParticle::fgPrimPmPar0= 0.;
21double AliFemtoParticle::fgPrimPmPar1= 0.;
22double AliFemtoParticle::fgPrimPmPar2= 0.;
23double AliFemtoParticle::fgPrimPpPar0= 0.;
24double AliFemtoParticle::fgPrimPpPar1= 0.;
25double AliFemtoParticle::fgPrimPpPar2= 0.;
67427ff7 26
27int TpcLocalTransform(AliFmThreeVectorD& xgl,
28 int& iSector,
29 int& iPadrow,
30 float& xlocal,
31 double& ttPhi);
32
33
34//_____________________
0215f606 35AliFemtoParticle::AliFemtoParticle() :
0b3bd1ac 36// fTpcV0NegPosSample(0),
37// fV0NegZ(0),
38// fV0NegU(0),
39// fV0NegSect(0),
0215f606 40 fTrack(0), fV0(0), fKink(0), fXi(0),
41 fFourMomentum(0),
42 fHelix(),
0b3bd1ac 43// fNominalTpcExitPoint(0),
44// fNominalTpcEntrancePoint(0),
0215f606 45 fHiddenInfo(0),
46 fPrimaryVertex(0),
47 fSecondaryVertex(0),
48 fHelixV0Pos(),
49 fTpcV0PosEntrancePoint(0),
50 fTpcV0PosExitPoint(0),
51 fHelixV0Neg(),
52 fTpcV0NegEntrancePoint(0),
53 fTpcV0NegExitPoint(0)
54{
d0e92d9a 55 // Default constructor
67427ff7 56 /* no-op for default */
57 // cout << "Created particle " << this << endl;
58}
59//_____________________
0215f606 60AliFemtoParticle::AliFemtoParticle(const AliFemtoParticle& aParticle):
0b3bd1ac 61// fTpcV0NegPosSample(0),
62// fV0NegZ(0),
63// fV0NegU(0),
64// fV0NegSect(0),
0215f606 65 fTrack(0), fV0(0), fKink(0), fXi(0),
66 fFourMomentum(0),
67 fHelix(),
0b3bd1ac 68// fNominalTpcExitPoint(0),
69// fNominalTpcEntrancePoint(0),
0215f606 70 fHiddenInfo(0),
71 fPrimaryVertex(0),
72 fSecondaryVertex(0),
73 fHelixV0Pos(),
74 fTpcV0PosEntrancePoint(0),
75 fTpcV0PosExitPoint(0),
76 fHelixV0Neg(),
77 fTpcV0NegEntrancePoint(0),
78 fTpcV0NegExitPoint(0)
79{
d0e92d9a 80 // Copy constructor
0215f606 81 if (aParticle.fTrack)
82 fTrack = new AliFemtoTrack(*aParticle.fTrack);
83 if (aParticle.fV0)
84 fV0 = new AliFemtoV0(*aParticle.fV0);
85 if (aParticle.fKink)
86 fKink = new AliFemtoKink(*aParticle.fKink);
87 if (aParticle.fXi)
88 fXi = new AliFemtoXi(*aParticle.fXi);
89 fFourMomentum = aParticle.fFourMomentum;
90 fHelix = aParticle.fHelix;
91
0b3bd1ac 92// for (int iter=0; iter<11; iter++)
93// fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
94
95// if (aParticle.fTpcV0NegPosSample) {
96// fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
97// for (int iter=0; iter<11; iter++)
98// fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
99// }
100
101// if (aParticle.fV0NegZ) {
102// fV0NegZ = (float *) malloc(sizeof(float) * 45);
103// for (int iter=0; iter<11; iter++)
104// fV0NegZ[iter] = aParticle.fV0NegZ[iter];
105// }
106// if (aParticle.fV0NegU) {
107// fV0NegU = (float *) malloc(sizeof(float) * 45);
108// for (int iter=0; iter<11; iter++)
109// fV0NegU[iter] = aParticle.fV0NegU[iter];
110// }
111// if (aParticle.fV0NegSect) {
112// fV0NegSect = (int *) malloc(sizeof(int) * 45);
113// for (int iter=0; iter<11; iter++)
114// fV0NegSect[iter] = aParticle.fV0NegSect[iter];
115// }
0215f606 116
117 fPrimaryVertex = aParticle.fPrimaryVertex;
118 fSecondaryVertex = aParticle.fSecondaryVertex;
119 CalculatePurity();
120 if(aParticle.fHiddenInfo){
d0e92d9a 121 fHiddenInfo= aParticle.HiddenInfo()->Clone();
0215f606 122 }
123
0b3bd1ac 124// fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
125// fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
0215f606 126
127 for (int iter=0; iter<6; iter++)
128 fPurity[iter] = aParticle.fPurity[iter];
129
130 fHelixV0Pos = aParticle.fHelixV0Pos;
131 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
132 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
133 fHelixV0Neg = aParticle.fHelixV0Neg;
134 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
135 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
136}
137//_____________________
67427ff7 138AliFemtoParticle::~AliFemtoParticle(){
139 // cout << "Issuing delete for AliFemtoParticle." << endl;
140
141 if (fTrack) delete fTrack;
142 if (fV0) {
0b3bd1ac 143// delete[] fTpcV0NegPosSample;
144// delete[] fV0NegZ;
145// delete[] fV0NegU;
146// delete[] fV0NegSect;
67427ff7 147 delete fV0;
148 }
149 if (fKink) delete fKink;
150 // cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
151 if (fHiddenInfo)
152 {
153 // cout << "Deleting HiddenInfo." << endl;
154 delete fHiddenInfo;
155 }
156 // cout << "Deleted particle " << this << endl;
157}
158//_____________________
0215f606 159AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) :
0b3bd1ac 160// fTpcV0NegPosSample(0),
161// fV0NegZ(0),
162// fV0NegU(0),
163// fV0NegSect(0),
0215f606 164 fTrack(0), fV0(0), fKink(0), fXi(0),
165 fFourMomentum(0),
166 fHelix(),
0b3bd1ac 167// fNominalTpcExitPoint(0),
168// fNominalTpcEntrancePoint(0),
0215f606 169 fHiddenInfo(0),
170 fPrimaryVertex(0),
171 fSecondaryVertex(0),
172 fHelixV0Pos(),
173 fTpcV0PosEntrancePoint(0),
174 fTpcV0PosExitPoint(0),
175 fHelixV0Neg(),
176 fTpcV0NegEntrancePoint(0),
177 fTpcV0NegExitPoint(0)
178{
d0e92d9a 179 // Constructor from normal track
67427ff7 180
181 // I know there is a better way to do this...
182 fTrack = new AliFemtoTrack(*hbtTrack);
183 AliFemtoThreeVector temp = hbtTrack->P();
184 fFourMomentum.setVect(temp);
185 double ener = ::sqrt(temp.mag2()+mass*mass);
186 fFourMomentum.setE(ener);
187// fMap[0] = hbtTrack->TopologyMap(0);
188 // fMap[1] = hbtTrack->TopologyMap(1);
189 // fNhits = hbtTrack->NHits();
190 fHelix = hbtTrack->Helix();
191 //CalculateNominalTpcExitAndEntrancePoints();
192
193
194 fPrimaryVertex.setX(0.);
195 fPrimaryVertex.setY(0.);
196 fPrimaryVertex.setZ(0.);
197 fSecondaryVertex.setX(0.);
198 fSecondaryVertex.setY(0.);
199 fSecondaryVertex.setZ(0.);
200 /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
201 CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
202 &fSecondaryVertex,
203 &fNominalTpcEntrancePoint,
204 &fNominalTpcExitPoint,
205 &mNominalPosSample[0],
206 &fZ[0],
207 &fU[0],
208 &fSect[0]);
209 */
210 CalculatePurity();
211 // ***
212 fHiddenInfo= 0;
213 if(hbtTrack->ValidHiddenInfo()){
d0e92d9a 214 fHiddenInfo= hbtTrack->GetHiddenInfo()->Clone();
67427ff7 215 }
216 // ***
217 // cout << "Created particle " << this << endl;
218
219}
220//_____________________
0215f606 221AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) :
0b3bd1ac 222// fTpcV0NegPosSample(0),
223// fV0NegZ(0),
224// fV0NegU(0),
225// fV0NegSect(0),
0215f606 226 fTrack(0), fV0(0), fKink(0), fXi(0),
227 fFourMomentum(0),
228 fHelix(),
0b3bd1ac 229// fNominalTpcExitPoint(0),
230// fNominalTpcEntrancePoint(0),
0215f606 231 fHiddenInfo(0),
232 fPrimaryVertex(0),
233 fSecondaryVertex(0),
234 fHelixV0Pos(),
235 fTpcV0PosEntrancePoint(0),
236 fTpcV0PosExitPoint(0),
237 fHelixV0Neg(),
238 fTpcV0NegEntrancePoint(0),
239 fTpcV0NegExitPoint(0)
240{
d0e92d9a 241 // Constructor from V0
67427ff7 242 fV0 = new AliFemtoV0(*hbtV0);
243 //fMap[0]= 0;
244 //fMap[1]= 0;
245 // I know there is a better way to do this...
d0e92d9a 246 AliFemtoThreeVector temp = hbtV0->MomV0();
67427ff7 247 fFourMomentum.setVect(temp);
248 double ener = ::sqrt(temp.mag2()+mass*mass);
249 fFourMomentum.setE(ener);
250 // Calculating TpcEntrancePoint for Positive V0 daugther
d0e92d9a 251 fPrimaryVertex = hbtV0->PrimaryVertex();
252 fSecondaryVertex = hbtV0->DecayVertexV0();
67427ff7 253 fHelixV0Pos = hbtV0->HelixPos();
254
0b3bd1ac 255// fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
256// fV0NegZ = new float[45];//for V0Neg
257// fV0NegU = new float[45];//for V0Neg
258// fV0NegSect = new int[45];//for V0Neg
259// CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
260// &fSecondaryVertex,
261// &fTpcV0PosEntrancePoint,
262// &fTpcV0PosExitPoint,
263// &fNominalPosSample[0],
264// &fZ[0],
265// &fU[0],&fSect[0]);
67427ff7 266 fHelixV0Neg = hbtV0->HelixNeg();
267
0b3bd1ac 268// CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
269// &fPrimaryVertex,
270// &fSecondaryVertex,
271// &fTpcV0NegEntrancePoint,
272// &fTpcV0NegExitPoint,
273// &fTpcV0NegPosSample[0],
274// &fV0NegZ[0],
275// &fV0NegU[0],&fV0NegSect[0]);
67427ff7 276
277 // ***
278 fHiddenInfo= 0;
279 if(hbtV0->ValidHiddenInfo()){
d0e92d9a 280 fHiddenInfo= hbtV0->GetHiddenInfo()->Clone();
67427ff7 281 }
282 // ***
283}
284//_____________________
0215f606 285AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) :
0b3bd1ac 286// fTpcV0NegPosSample(0),
287// fV0NegZ(0),
288// fV0NegU(0),
289// fV0NegSect(0),
0215f606 290 fTrack(0), fV0(0), fKink(0), fXi(0),
291 fFourMomentum(0),
292 fHelix(),
0b3bd1ac 293// fNominalTpcExitPoint(0),
294// fNominalTpcEntrancePoint(0),
0215f606 295 fHiddenInfo(0),
296 fPrimaryVertex(0),
297 fSecondaryVertex(0),
298 fHelixV0Pos(),
299 fTpcV0PosEntrancePoint(0),
300 fTpcV0PosExitPoint(0),
301 fHelixV0Neg(),
302 fTpcV0NegEntrancePoint(0),
303 fTpcV0NegExitPoint(0)
304{
d0e92d9a 305 // Constructor from Kink
67427ff7 306 fKink = new AliFemtoKink(*hbtKink);
307 // fMap[0]= 0;
308 //fMap[1]= 0;
309 // I know there is a better way to do this...
310 AliFemtoThreeVector temp = hbtKink->Parent().P();
311 fFourMomentum.setVect(temp);
312 double ener = ::sqrt(temp.mag2()+mass*mass);
313 fFourMomentum.setE(ener);
314}
315
316//_____________________
0215f606 317AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) :
0b3bd1ac 318// fTpcV0NegPosSample(0),
319// fV0NegZ(0),
320// fV0NegU(0),
321// fV0NegSect(0),
0215f606 322 fTrack(0), fV0(0), fKink(0), fXi(0),
323 fFourMomentum(0),
324 fHelix(),
0b3bd1ac 325// fNominalTpcExitPoint(0),
326// fNominalTpcEntrancePoint(0),
0215f606 327 fHiddenInfo(0),
328 fPrimaryVertex(0),
329 fSecondaryVertex(0),
330 fHelixV0Pos(),
331 fTpcV0PosEntrancePoint(0),
332 fTpcV0PosExitPoint(0),
333 fHelixV0Neg(),
334 fTpcV0NegEntrancePoint(0),
335 fTpcV0NegExitPoint(0)
336{
d0e92d9a 337 // Constructor from Xi
67427ff7 338 fXi = new AliFemtoXi(*hbtXi);
339 // fMap[0]= 0;
340 //fMap[1]= 0;
341 AliFemtoThreeVector temp;// = hbtXi->mMofXi;
342 fFourMomentum.setVect(temp);
343 double ener = ::sqrt(temp.mag2()+mass*mass);
344 fFourMomentum.setE(ener);
345 fHiddenInfo = 0;
346}
347//_____________________
0215f606 348AliFemtoParticle& AliFemtoParticle::operator=(const AliFemtoParticle& aParticle)
349{
d0e92d9a 350 // assignment operator
0215f606 351 if (this == &aParticle)
352 return *this;
353
354 if (aParticle.fTrack)
355 fTrack = new AliFemtoTrack(*aParticle.fTrack);
356 if (aParticle.fV0)
357 fV0 = new AliFemtoV0(*aParticle.fV0);
358 if (aParticle.fKink)
359 fKink = new AliFemtoKink(*aParticle.fKink);
360 if (aParticle.fXi)
361 fXi = new AliFemtoXi(*aParticle.fXi);
362 fFourMomentum = aParticle.fFourMomentum;
363 fHelix = aParticle.fHelix;
364
0b3bd1ac 365// for (int iter=0; iter<11; iter++)
366// fNominalPosSample[iter] = aParticle.fNominalPosSample[iter];
367
368// if (fTpcV0NegPosSample) delete fTpcV0NegPosSample;
369// if (aParticle.fTpcV0NegPosSample) {
370// fTpcV0NegPosSample = (AliFemtoThreeVector *) malloc(sizeof(AliFemtoThreeVector) * 11);
371// for (int iter=0; iter<11; iter++)
372// fTpcV0NegPosSample[iter] = aParticle.fTpcV0NegPosSample[iter];
373// }
374
375// if (fV0NegZ) delete fV0NegZ;
376// if (aParticle.fV0NegZ) {
377// fV0NegZ = (float *) malloc(sizeof(float) * 45);
378// for (int iter=0; iter<11; iter++)
379// fV0NegZ[iter] = aParticle.fV0NegZ[iter];
380// }
381// if (fV0NegU) delete fV0NegU;
382// if (aParticle.fV0NegU) {
383// fV0NegU = (float *) malloc(sizeof(float) * 45);
384// for (int iter=0; iter<11; iter++)
385// fV0NegU[iter] = aParticle.fV0NegU[iter];
386// }
387// if (fV0NegSect) delete fV0NegSect;
388// if (aParticle.fV0NegSect) {
389// fV0NegSect = (int *) malloc(sizeof(int) * 45);
390// for (int iter=0; iter<11; iter++)
391// fV0NegSect[iter] = aParticle.fV0NegSect[iter];
392// }
0215f606 393
394 fPrimaryVertex = aParticle.fPrimaryVertex;
395 fSecondaryVertex = aParticle.fSecondaryVertex;
396 CalculatePurity();
397 if(aParticle.fHiddenInfo){
d0e92d9a 398 fHiddenInfo= aParticle.fHiddenInfo->Clone();
0215f606 399 }
400
0b3bd1ac 401// fNominalTpcEntrancePoint = aParticle.fNominalTpcEntrancePoint;
402// fNominalTpcExitPoint = aParticle.fNominalTpcExitPoint;
403
0215f606 404 if (fHiddenInfo) delete fHiddenInfo;
405 if (aParticle.fHiddenInfo)
d0e92d9a 406 fHiddenInfo = aParticle.HiddenInfo()->Clone();
0215f606 407
408 for (int iter=0; iter<6; iter++)
409 fPurity[iter] = aParticle.fPurity[iter];
410
411 fHelixV0Pos = aParticle.fHelixV0Pos;
412 fTpcV0PosEntrancePoint = aParticle.fTpcV0PosEntrancePoint;
413 fTpcV0PosExitPoint = aParticle.fTpcV0PosExitPoint;
414 fHelixV0Neg = aParticle.fHelixV0Neg;
415 fTpcV0NegEntrancePoint = aParticle.fTpcV0NegEntrancePoint;
416 fTpcV0NegExitPoint = aParticle.fTpcV0NegExitPoint;
417
418 return *this;
419}
0b3bd1ac 420// //_____________________
421// const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
422// // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
423// // for now, we calculate Exit and Entrance points upon instantiation
424// return fNominalTpcExitPoint;
425// }
426// //_____________________
427// const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
428// // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
429// // for now, we calculate Exit and Entrance points upon instantiation
430// return fNominalTpcEntrancePoint;
431// }
67427ff7 432//_____________________
433void AliFemtoParticle::CalculatePurity(){
d0e92d9a 434 // Calculate additional parameterized purity
435
67427ff7 436 double tPt = fFourMomentum.perp();
437 // pi -
d0e92d9a 438 fPurity[0] = fgPrimPimPar0*(1.-exp((tPt-fgPrimPimPar1)/fgPrimPimPar2));
67427ff7 439 fPurity[0] *= fTrack->PidProbPion();
440 // pi+
d0e92d9a 441 fPurity[1] = fgPrimPipPar0*(1.-exp((tPt-fgPrimPipPar1)/fgPrimPipPar2));
67427ff7 442 fPurity[1] *= fTrack->PidProbPion();
443 // K-
444 fPurity[2] = fTrack->PidProbKaon();
445 // K+
446 fPurity[3] = fTrack->PidProbKaon();
447 // pbar
448 fPurity[4] = fTrack->PidProbProton();
449 // p
450 fPurity[5] = fTrack->PidProbProton();
451}
452
453double AliFemtoParticle::GetPionPurity()
454{
d0e92d9a 455 // Get full pion purity
67427ff7 456 if (fTrack->Charge()>0)
457 return fPurity[1];
458 else
459 return fPurity[0];
460}
461double AliFemtoParticle::GetKaonPurity()
462{
d0e92d9a 463 // Get full kaon purity
67427ff7 464 if (fTrack->Charge()>0)
465 return fPurity[3];
466 else
467 return fPurity[2];
468}
469double AliFemtoParticle::GetProtonPurity()
470{
d0e92d9a 471 // Get full proton purity
67427ff7 472 if (fTrack->Charge()>0)
473 return fPurity[5];
474 else
475 return fPurity[4];
476}
477
0b3bd1ac 478// void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
479// AliFemtoThreeVector* PrimVert,
480// AliFemtoThreeVector* SecVert,
481// AliFemtoThreeVector* tmpTpcEntrancePoint,
482// AliFemtoThreeVector* tmpTpcExitPoint,
483// AliFemtoThreeVector* tmpPosSample,
484// float* tmpZ,
485// float* tmpU,
486// int* tmpSect){
487// // this calculates the exit point of a secondary track,
488// // either through the endcap or through the Outer Field Cage
489// // We assume the track to start at tHelix.origin-PrimaryVertex
490// // it also calculates the entrance point of the secondary track,
491// // which is the point at which it crosses the
492// // inner field cage
493// // static AliFemtoThreeVector ZeroVec(0.,0.,0.);
494// AliFemtoThreeVector tZeroVec(0.,0.,0.);
495// // tZeroVec.setX(tHelix->origin().x()-PrimVert->x());
496// // tZeroVec.setY(tHelix->origin().y()-PrimVert->y());
497// // tZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
498// tZeroVec.setX(SecVert->x()-PrimVert->x());
499// tZeroVec.setY(SecVert->y()-PrimVert->y());
500// tZeroVec.setZ(SecVert->z()-PrimVert->z());
501// double dip, curv, phase;
502// int h;
503// curv = tHelix->Curvature();
504// dip = tHelix->DipAngle();
505// phase= tHelix->Phase();
506// h = tHelix->H();
67427ff7 507
0b3bd1ac 508// AliFmHelixD hel(curv,dip,phase,tZeroVec,h);
509
510// pairD candidates;
511// double sideLength; // this is how much length to go to leave through sides of TPC
512// double endLength; // this is how much length to go to leave through endcap of TPC
513// // figure out how far to go to leave through side...
514// candidates = hel.PathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
515// sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
516
517// static AliFemtoThreeVector tWestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
518// static AliFemtoThreeVector tEastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
519// static AliFemtoThreeVector tEndCapNormal(0.,0.,1.0);
520
521// endLength = hel.PathLength(tWestEnd,tEndCapNormal);
522// if (endLength < 0.0) endLength = hel.PathLength(tEastEnd,tEndCapNormal);
523
524// if (endLength < 0.0) cout <<
525// "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
526// << "Hey -- I cannot find an exit point out endcaps" << endl;
527// // OK, firstExitLength will be the shortest way out of the detector...
528// double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
529// // now then, let's return the POSITION at which particle leaves TPC...
530// *tmpTpcExitPoint = hel.At(firstExitLength);
531// // Finally, calculate the position at which the track crosses the inner field cage
532// candidates = hel.PathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
533
534// sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
535// // cout << "sideLength 2 ="<<sideLength << endl;
536// *tmpTpcEntrancePoint = hel.At(sideLength);
537// // This is the secure way !
538// if (IsNaN(tmpTpcEntrancePoint->x()) ||
539// IsNaN(tmpTpcEntrancePoint->y()) ||
540// IsNaN(tmpTpcEntrancePoint->z()) ){
541// cout << "tmpTpcEntrancePoint NAN"<< endl;
542// cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
543// tmpTpcEntrancePoint->setX(-9999.);
544// tmpTpcEntrancePoint->setY(-9999.);
545// tmpTpcEntrancePoint->setZ(-9999.);
546// }
67427ff7 547
0b3bd1ac 548// if (IsNaN(tmpTpcExitPoint->x()) ||
549// IsNaN(tmpTpcExitPoint->y()) ||
550// IsNaN(tmpTpcExitPoint->z()) ) {
551// // cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
552// // cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
553// // cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
554// // cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
555// tmpTpcExitPoint->setX(-9999.);
556// tmpTpcExitPoint->setY(-9999.);
557// tmpTpcExitPoint->setZ(-9999.);
558// }
559
560
561// // if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
562// // if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
563// // cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
564
565// // 03Oct00 - mal. OK, let's try something a little more
566// // along the lines of NA49 and E895 strategy.
567// // calculate the "nominal" position at N radii (say N=11)
568// // within the TPC, and for a pair cut
569// // use the average separation of these N
570// int irad = 0;
571// candidates = hel.PathLength(50.0);
572// sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
573// while (irad<11 && !IsNaN(sideLength)){
574// float radius = 50.0 + irad*15.0;
575// candidates = hel.PathLength(radius);
576// sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
577// tmpPosSample[irad] = hel.At(sideLength);
578// if(IsNaN(tmpPosSample[irad].x()) ||
579// IsNaN(tmpPosSample[irad].y()) ||
580// IsNaN(tmpPosSample[irad].z())
581// ){
582// cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
583// cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
584// tmpPosSample[irad] = AliFemtoThreeVector(-9999.,-9999.,-9999);
585// }
586// irad++;
587// if (irad<11){
588// float radius = 50.0 + irad*15.0;
589// candidates = hel.PathLength(radius);
590// sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
591// }
592// }
593// for (int i = irad; i<11; i++)
594// {
595// tmpPosSample[i] = AliFemtoThreeVector(-9999.,-9999.,-9999);
596// }
597
598// static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
599// 104,109.2,114.4,119.6,127.195,129.195,131.195,
600// 133.195,135.195,137.195,139.195,141.195,
601// 143.195,145.195,147.195,149.195,151.195,
602// 153.195,155.195,157.195,159.195,161.195,
603// 163.195,165.195,167.195,169.195,171.195,
604// 173.195,175.195,177.195,179.195,181.195,
605// 183.195,185.195,187.195,189.195};
606// int tRow,tSect,tOutOfBound;
607// double tLength,tPhi;
608// float tU;
609// AliFemtoThreeVector tPoint;
610// AliFmThreeVectorD tn(0,0,0);
611// AliFmThreeVectorD tr(0,0,0);
612// int ti =0;
613// // test to enter the loop
614// candidates = hel.PathLength(tRowRadius[ti]);
615// tLength = (candidates.first > 0) ? candidates.first : candidates.second;
616// if (IsNaN(tLength)){
617// cout <<"tLength Init tmp NAN" << endl;
618// cout <<"padrow number= "<<ti << "not reached" << endl;
619// cout << "*** DO NOT ENTER THE LOOP***" << endl;
620// tmpSect[ti]=-1;//sector
621// }
622// // end test
623// while(ti<45 && !IsNaN(tLength)){
624// candidates = hel.PathLength(tRowRadius[ti]);
625// tLength = (candidates.first > 0) ? candidates.first : candidates.second;
626// if (IsNaN(tLength)){
627// cout <<"tLength loop 1st NAN" << endl;
628// cout <<"padrow number= " << ti << " not reached" << endl;
629// cout << "*** THIS IS AN ERROR SHOULDN'T LOOP ***" << endl;
630// tmpSect[ti]=-1;//sector
631// }
632// tPoint = hel.At(tLength);
633// // Find which sector it is on
634// TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
635// if (IsNaN(tmpSect[ti])){
636// cout <<"***ERROR tmpSect"<< endl;
637// }
638// if (IsNaN(tRow)){
639// cout <<"***ERROR tRow"<< endl;
640// }
641// if (IsNaN(tU)){
642// cout <<"***ERROR tU"<< endl;
643// }
644// if (IsNaN(tPhi)){
645// cout <<"***ERROR tPhi"<< endl;
646// }
647// // calculate crossing plane
648// tn.setX(cos(tPhi));
649// tn.setY(sin(tPhi));
650// tr.setX(tRowRadius[ti]*cos(tPhi));
651// tr.setY(tRowRadius[ti]*sin(tPhi));
652// // find crossing point
653// tLength = hel.PathLength(tr,tn);
654// if (IsNaN(tLength)){
655// cout <<"tLength loop 2nd NAN" << endl;
656// cout <<"padrow number= " << ti << " not reached" << endl;
657// tmpSect[ti]=-2;//sector
658// }
659// tPoint = hel.At(tLength);
660// tmpZ[ti] = tPoint.z();
661// tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
662// if (IsNaN(tSect)){
663// cout <<"***ERROR tSect 2"<< endl;
664// }
665// if (IsNaN(tRow)){
666// cout <<"***ERROR tRow 2"<< endl;
667// }
668// if (IsNaN(tmpU[ti])){
669// cout <<"***ERROR tmpU[ti] 2"<< endl;
670// }
671// if (IsNaN(tPhi)){
672// cout <<"***ERROR tPhi 2 "<< endl;
673// }
674// if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
675// tmpSect[ti]=-2;
676// // cout << "missed once"<< endl;
677// }
678// else{
679// if(tmpSect[ti] != tSect){
680// // Try again on the other sector
681// tn.setX(cos(tPhi));
682// tn.setY(sin(tPhi));
683// tr.setX(tRowRadius[ti]*cos(tPhi));
684// tr.setY(tRowRadius[ti]*sin(tPhi));
685// // find crossing point
686// tLength = hel.PathLength(tr,tn);
687// tPoint = hel.At(tLength);
688// if (IsNaN(tLength)){
689// cout <<"tLength loop 3rd NAN" << endl;
690// cout <<"padrow number= "<< ti << " not reached" << endl;
691// tmpSect[ti]=-1;//sector
692// }
693// tmpZ[ti] = tPoint.z();
694// tmpSect[ti] = tSect;
695// tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
696// if (IsNaN(tSect)){
697// cout <<"***ERROR tSect 3"<< endl;
698// }
699// if (IsNaN(tRow)){
700// cout <<"***ERROR tRow 3"<< endl;
701// }
702// if (IsNaN(tmpU[ti])){
703// cout <<"***ERROR tmpU[ti] 3"<< endl;
704// }
705// if (IsNaN(tPhi)){
706// cout <<"***ERROR tPhi 3 "<< endl;
707// }
708// if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
709// tmpSect[ti]=-1;
710// }
711// }
712// }
713// if (IsNaN(tmpSect[ti])){
714// cout << "*******************ERROR***************************" << endl;
715// cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
716// cout << "*******************ERROR***************************" << endl;
717// }
718// if (IsNaN(tmpU[ti])){
719// cout << "*******************ERROR***************************" << endl;
720// cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
721// cout << "*******************ERROR***************************" << endl;
722// }
723// if (IsNaN(tmpZ[ti])){
724// cout << "*******************ERROR***************************" << endl;
725// cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
726// cout << "*******************ERROR***************************" << endl;
727// }
728// // If padrow ti not reached all other beyond are not reached
729// // in this case set sector to -1
730// if (tmpSect[ti]==-1){
731// for (int tj=ti; tj<45;tj++){
732// tmpSect[tj] = -1;
733// ti=45;
734// }
735// }
736// ti++;
737// if (ti<45){
738// candidates = hel.PathLength(tRowRadius[ti]);
739// tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
740// }
741// }
67427ff7 742//_____________________
743const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
744 return fTpcV0PosExitPoint;
745}
746//_____________________
747const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
748 return fTpcV0PosEntrancePoint;
749}
750//______________________
751const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
752 return fTpcV0NegExitPoint;
753}
754//_____________________
755const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
756 return fTpcV0NegEntrancePoint;
757}
758//______________________