]> 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
This commit was generated by cvs2svn to compensate for changes in r18145,
[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$
17 * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki
18 * First version on CVS
19 *
20 * Revision 1.22 2003/09/02 17:58:32 perev
21 * gcc 3.2 updates + WarnOff
22 *
23 * Revision 1.21 2003/05/07 15:30:43 magestro
24 * Fixed bug related to finding merged hits (commit for Fabrice)
25 *
26 * Revision 1.20 2003/01/14 09:41:26 renault
27 * changes on average separation calculation, hit shared finder and memory optimisation
28 * for Z,U and Sectors variables.
29 *
30 * Revision 1.19 2002/12/12 17:01:49 kisiel
31 * Hidden Information handling and purity calculation
32 *
33 * Revision 1.18 2002/11/19 23:36:00 renault
34 * Enable calculation of exit/entrance separation for V0 daughters
35 *
36 * Revision 1.17 2001/12/14 23:11:30 fretiere
37 * 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
38 *
39 * Revision 1.16 2001/05/25 23:23:59 lisa
40 * Added in AliFemtoKink stuff
41 *
42 * Revision 1.15 2001/04/03 21:04:36 kisiel
43 *
44 *
45 * Changes needed to make the Theoretical code
46 * work. The main code is the ThCorrFctn directory.
47 * The most visible change is the addition of the
48 * HiddenInfo to AliFemtoPair.
49 *
50 * Revision 1.14 2000/10/05 23:09:05 lisa
51 * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
52 *
53 * Revision 1.13 2000/08/31 22:31:31 laue
54 * AliFemtoAnalysis: output changed (a little bit less)
55 * AliFemtoEvent: new version, members for reference mult added
56 * AliFemtoIOBinary: new IO for new AliFemtoEvent version
57 * AliFemtoTypes: TTree typedef to AliFemtoTTree added
58 * AliFemtoVertexAnalysis: overflow and underflow added
59 *
60 * Revision 1.12 2000/07/23 13:52:56 laue
61 * NominalExitPoint set to (-9999.,-9999.-9999.) if helix.at()
62 * returns nan (not a number).
63 *
64 * Revision 1.11 2000/07/19 17:18:48 laue
65 * Calculation of Entrance and Exit point added in AliFemtoParticle constructor
66 *
67 * Revision 1.10 2000/07/17 20:03:17 lisa
68 * Implemented tools for addressing and assessing trackmerging
69 *
70 * Revision 1.9 2000/07/16 21:38:23 laue
71 * AliFemtoCoulomb.cxx AliFemtoSectoredAnalysis.cxx : updated for standalone version
72 * AliFemtoV0.cc AliFemtoV0.h : some cast to prevent compiling warnings
73 * AliFemtoParticle.cc AliFemtoParticle.h : pointers fTrack,fV0 initialized to 0
74 * AliFemtoIOBinary.cc : some printouts in #ifdef STHBTDEBUG
75 * AliFemtoEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
76 * solution
77 *
78 * Revision 1.8 2000/05/03 17:44:43 laue
79 * AliFemtoEvent, AliFemtoTrack & AliFemtoV0 declared friend to AliFemtoIOBinary
80 * AliFemtoParticle updated for V0 pos,neg track Id
81 *
82 * Revision 1.7 2000/04/03 16:21:51 laue
83 * some include files changed
84 * Multi track cut added
85 *
86 * Revision 1.6 1999/12/11 15:58:29 lisa
87 * Add vertex decay position datum and accessor to AliFemtoParticle to allow pairwise cuts on seperation of V0s
88 *
89 * Revision 1.5 1999/09/17 22:38:02 lisa
90 * first full integration of V0s into AliFemto framework
91 *
92 * Revision 1.4 1999/09/01 19:04:53 lisa
93 * update Particle class AND add parity cf and Randys Coulomb correction
94 *
95 * Revision 1.3 1999/07/06 22:33:23 lisa
96 * Adjusted all to work in pro and new - dev itself is broken
97 *
98 * Revision 1.2 1999/06/29 17:50:27 fisyak
99 * formal changes to account new StEvent, does not complie yet
100 *
101 * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
102 * Installation of AliFemtoMaker
103 *
104 **************************************************************************/
105
106#include "Infrastructure/AliFemtoParticle.h"
107//#include "math_constants.h"
108#ifdef __CC5__
109 #include <math.h>
110#else
111 #include <cmath>
112#endif
113
114
115#include "TMath.h"
116using namespace TMath;
117
118double AliFemtoParticle::fPrimPimPar0= 9.05632e-01;
119double AliFemtoParticle::fPrimPimPar1= -2.26737e-01;
120double AliFemtoParticle::fPrimPimPar2= -1.03922e-01;
121double AliFemtoParticle::fPrimPipPar0= 9.09616e-01;
122double AliFemtoParticle::fPrimPipPar1= -9.00511e-02;
123double AliFemtoParticle::fPrimPipPar2= -6.02940e-02;
124double AliFemtoParticle::fPrimPmPar0= 0.;
125double AliFemtoParticle::fPrimPmPar1= 0.;
126double AliFemtoParticle::fPrimPmPar2= 0.;
127double AliFemtoParticle::fPrimPpPar0= 0.;
128double AliFemtoParticle::fPrimPpPar1= 0.;
129double AliFemtoParticle::fPrimPpPar2= 0.;
130
131int TpcLocalTransform(AliFmThreeVectorD& xgl,
132 int& iSector,
133 int& iPadrow,
134 float& xlocal,
135 double& ttPhi);
136
137
138//_____________________
139AliFemtoParticle::AliFemtoParticle() : fTrack(0), fV0(0), fKink(0), fHiddenInfo(0) {
140 /* no-op for default */
141 // cout << "Created particle " << this << endl;
142}
143//_____________________
144AliFemtoParticle::~AliFemtoParticle(){
145 // cout << "Issuing delete for AliFemtoParticle." << endl;
146
147 if (fTrack) delete fTrack;
148 if (fV0) {
149 delete[] fTpcV0NegPosSample;
150 delete[] fV0NegZ;
151 delete[] fV0NegU;
152 delete[] fV0NegSect;
153 delete fV0;
154 }
155 if (fKink) delete fKink;
156 // cout << "Trying to delete HiddenInfo: " << fHiddenInfo << endl;
157 if (fHiddenInfo)
158 {
159 // cout << "Deleting HiddenInfo." << endl;
160 delete fHiddenInfo;
161 }
162 // cout << "Deleted particle " << this << endl;
163}
164//_____________________
165AliFemtoParticle::AliFemtoParticle(const AliFemtoTrack* const hbtTrack,const double& mass) : fTrack(0), fV0(0), fKink(0), fHiddenInfo(0) {
166
167
168 // I know there is a better way to do this...
169 fTrack = new AliFemtoTrack(*hbtTrack);
170 AliFemtoThreeVector temp = hbtTrack->P();
171 fFourMomentum.setVect(temp);
172 double ener = ::sqrt(temp.mag2()+mass*mass);
173 fFourMomentum.setE(ener);
174// fMap[0] = hbtTrack->TopologyMap(0);
175 // fMap[1] = hbtTrack->TopologyMap(1);
176 // fNhits = hbtTrack->NHits();
177 fHelix = hbtTrack->Helix();
178 //CalculateNominalTpcExitAndEntrancePoints();
179
180
181 fPrimaryVertex.setX(0.);
182 fPrimaryVertex.setY(0.);
183 fPrimaryVertex.setZ(0.);
184 fSecondaryVertex.setX(0.);
185 fSecondaryVertex.setY(0.);
186 fSecondaryVertex.setZ(0.);
187 /* TO JA ODZNACZYLEM NIE WIEM DLACZEGO
188 CalculateTpcExitAndEntrancePoints(&fHelix,&fPrimaryVertex,
189 &fSecondaryVertex,
190 &fNominalTpcEntrancePoint,
191 &fNominalTpcExitPoint,
192 &mNominalPosSample[0],
193 &fZ[0],
194 &fU[0],
195 &fSect[0]);
196 */
197 CalculatePurity();
198 // ***
199 fHiddenInfo= 0;
200 if(hbtTrack->ValidHiddenInfo()){
201 fHiddenInfo= hbtTrack->getHiddenInfo()->getParticleHiddenInfo()->clone();
202 }
203 // ***
204 // cout << "Created particle " << this << endl;
205
206}
207//_____________________
208AliFemtoParticle::AliFemtoParticle(const AliFemtoV0* const hbtV0,const double& mass) : fTrack(0), fV0(0), fKink(0), fHiddenInfo(0) {
209 fV0 = new AliFemtoV0(*hbtV0);
210 //fMap[0]= 0;
211 //fMap[1]= 0;
212 // I know there is a better way to do this...
213 AliFemtoThreeVector temp = hbtV0->momV0();
214 fFourMomentum.setVect(temp);
215 double ener = ::sqrt(temp.mag2()+mass*mass);
216 fFourMomentum.setE(ener);
217 // Calculating TpcEntrancePoint for Positive V0 daugther
218 fPrimaryVertex = hbtV0->primaryVertex();
219 fSecondaryVertex = hbtV0->decayVertexV0();
220 fHelixV0Pos = hbtV0->HelixPos();
221
222 fTpcV0NegPosSample = new AliFemtoThreeVector[45];//for V0Neg
223 fV0NegZ = new float[45];//for V0Neg
224 fV0NegU = new float[45];//for V0Neg
225 fV0NegSect = new int[45];//for V0Neg
226 CalculateTpcExitAndEntrancePoints(&fHelixV0Pos,&fPrimaryVertex,
227 &fSecondaryVertex,
228 &fTpcV0PosEntrancePoint,
229 &fTpcV0PosExitPoint,
230 &fNominalPosSample[0],
231 &fZ[0],
232 &fU[0],&fSect[0]);
233 fHelixV0Neg = hbtV0->HelixNeg();
234
235 CalculateTpcExitAndEntrancePoints(&fHelixV0Neg,
236 &fPrimaryVertex,
237 &fSecondaryVertex,
238 &fTpcV0NegEntrancePoint,
239 &fTpcV0NegExitPoint,
240 &fTpcV0NegPosSample[0],
241 &fV0NegZ[0],
242 &fV0NegU[0],&fV0NegSect[0]);
243
244 // ***
245 fHiddenInfo= 0;
246 if(hbtV0->ValidHiddenInfo()){
247 fHiddenInfo= hbtV0->getHiddenInfo()->clone();
248 }
249 // ***
250}
251//_____________________
252AliFemtoParticle::AliFemtoParticle(const AliFemtoKink* const hbtKink,const double& mass) : fTrack(0), fV0(0), fHiddenInfo(0) {
253 fKink = new AliFemtoKink(*hbtKink);
254 // fMap[0]= 0;
255 //fMap[1]= 0;
256 // I know there is a better way to do this...
257 AliFemtoThreeVector temp = hbtKink->Parent().P();
258 fFourMomentum.setVect(temp);
259 double ener = ::sqrt(temp.mag2()+mass*mass);
260 fFourMomentum.setE(ener);
261}
262
263//_____________________
264AliFemtoParticle::AliFemtoParticle(const AliFemtoXi* const hbtXi, const double& mass) {
265 fXi = new AliFemtoXi(*hbtXi);
266 // fMap[0]= 0;
267 //fMap[1]= 0;
268 AliFemtoThreeVector temp;// = hbtXi->mMofXi;
269 fFourMomentum.setVect(temp);
270 double ener = ::sqrt(temp.mag2()+mass*mass);
271 fFourMomentum.setE(ener);
272 fHiddenInfo = 0;
273}
274//_____________________
275const AliFemtoThreeVector& AliFemtoParticle::NominalTpcExitPoint() const{
276 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
277 // for now, we calculate Exit and Entrance points upon instantiation
278 return fNominalTpcExitPoint;
279}
280//_____________________
281const AliFemtoThreeVector& AliFemtoParticle::NominalTpcEntrancePoint() const{
282 // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
283 // for now, we calculate Exit and Entrance points upon instantiation
284 return fNominalTpcEntrancePoint;
285}
286//_____________________
287void AliFemtoParticle::CalculatePurity(){
288 double tPt = fFourMomentum.perp();
289 // pi -
290 fPurity[0] = fPrimPimPar0*(1.-exp((tPt-fPrimPimPar1)/fPrimPimPar2));
291 fPurity[0] *= fTrack->PidProbPion();
292 // pi+
293 fPurity[1] = fPrimPipPar0*(1.-exp((tPt-fPrimPipPar1)/fPrimPipPar2));
294 fPurity[1] *= fTrack->PidProbPion();
295 // K-
296 fPurity[2] = fTrack->PidProbKaon();
297 // K+
298 fPurity[3] = fTrack->PidProbKaon();
299 // pbar
300 fPurity[4] = fTrack->PidProbProton();
301 // p
302 fPurity[5] = fTrack->PidProbProton();
303}
304
305double AliFemtoParticle::GetPionPurity()
306{
307 if (fTrack->Charge()>0)
308 return fPurity[1];
309 else
310 return fPurity[0];
311}
312double AliFemtoParticle::GetKaonPurity()
313{
314 if (fTrack->Charge()>0)
315 return fPurity[3];
316 else
317 return fPurity[2];
318}
319double AliFemtoParticle::GetProtonPurity()
320{
321 if (fTrack->Charge()>0)
322 return fPurity[5];
323 else
324 return fPurity[4];
325}
326
327void AliFemtoParticle::CalculateTpcExitAndEntrancePoints(AliFmPhysicalHelixD* tHelix,
328 AliFemtoThreeVector* PrimVert,
329 AliFemtoThreeVector* SecVert,
330 AliFemtoThreeVector* tmpTpcEntrancePoint,
331 AliFemtoThreeVector* tmpTpcExitPoint,
332 AliFemtoThreeVector* tmpPosSample,
333 float* tmpZ,
334 float* tmpU,
335 int* tmpSect){
336 // this calculates the exit point of a secondary track,
337 // either through the endcap or through the Outer Field Cage
338 // We assume the track to start at tHelix.origin-PrimaryVertex
339 // it also calculates the entrance point of the secondary track,
340 // which is the point at which it crosses the
341 // inner field cage
342 // static AliFemtoThreeVector ZeroVec(0.,0.,0.);
343 AliFemtoThreeVector ZeroVec(0.,0.,0.);
344// ZeroVec.setX(tHelix->origin().x()-PrimVert->x());
345// ZeroVec.setY(tHelix->origin().y()-PrimVert->y());
346// ZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
347 ZeroVec.setX(SecVert->x()-PrimVert->x());
348 ZeroVec.setY(SecVert->y()-PrimVert->y());
349 ZeroVec.setZ(SecVert->z()-PrimVert->z());
350 double dip, curv, phase;
351 int h;
352 curv = tHelix->curvature();
353 dip = tHelix->dipAngle();
354 phase= tHelix->phase();
355 h = tHelix->h();
356
357 AliFmHelixD hel(curv,dip,phase,ZeroVec,h);
358
359 pairD candidates;
360 double sideLength; // this is how much length to go to leave through sides of TPC
361 double endLength; // this is how much length to go to leave through endcap of TPC
362 // figure out how far to go to leave through side...
363 candidates = hel.pathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
364 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
365
366 static AliFemtoThreeVector WestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
367 static AliFemtoThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
368 static AliFemtoThreeVector EndCapNormal(0.,0.,1.0);
369
370 endLength = hel.pathLength(WestEnd,EndCapNormal);
371 if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
372
373 if (endLength < 0.0) cout <<
374 "AliFemtoParticle::CalculateTpcExitAndEntrancePoints(): "
375 << "Hey -- I cannot find an exit point out endcaps" << endl;
376 // OK, firstExitLength will be the shortest way out of the detector...
377 double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
378 // now then, let's return the POSITION at which particle leaves TPC...
379 *tmpTpcExitPoint = hel.at(firstExitLength);
380 // Finally, calculate the position at which the track crosses the inner field cage
381 candidates = hel.pathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
382
383 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
384 // cout << "sideLength 2 ="<<sideLength << endl;
385 *tmpTpcEntrancePoint = hel.at(sideLength);
386 // This is the secure way !
387 if (IsNaN(tmpTpcEntrancePoint->x()) ||
388 IsNaN(tmpTpcEntrancePoint->y()) ||
389 IsNaN(tmpTpcEntrancePoint->z()) ){
390 cout << "tmpTpcEntrancePoint NAN"<< endl;
391 cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
392 tmpTpcEntrancePoint->setX(-9999.);
393 tmpTpcEntrancePoint->setY(-9999.);
394 tmpTpcEntrancePoint->setZ(-9999.);
395 }
396
397 if (IsNaN(tmpTpcExitPoint->x()) ||
398 IsNaN(tmpTpcExitPoint->y()) ||
399 IsNaN(tmpTpcExitPoint->z()) ) {
400// cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
401// cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
402// cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
403// cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
404 tmpTpcExitPoint->setX(-9999.);
405 tmpTpcExitPoint->setY(-9999.);
406 tmpTpcExitPoint->setZ(-9999.);
407 }
408
409
410// if (IsNaN(tmpTpcExitPoint->x())) *tmpTpcExitPoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
411// if (IsNaN(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = AliFemtoThreeVector(-9999.,-9999.,-9999);
412 // cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
413
414 // 03Oct00 - mal. OK, let's try something a little more
415 // along the lines of NA49 and E895 strategy.
416 // calculate the "nominal" position at N radii (say N=11)
417 // within the TPC, and for a pair cut
418 // use the average separation of these N
419 int irad = 0;
420 candidates = hel.pathLength(50.0);
421 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
422 while (irad<11 && !IsNaN(sideLength)){
423 float radius = 50.0 + irad*15.0;
424 candidates = hel.pathLength(radius);
425 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
426 tmpPosSample[irad] = hel.at(sideLength);
427 if(IsNaN(tmpPosSample[irad].x()) ||
428 IsNaN(tmpPosSample[irad].y()) ||
429 IsNaN(tmpPosSample[irad].z())
430 ){
431 cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
432 cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
433 tmpPosSample[irad] = AliFemtoThreeVector(-9999.,-9999.,-9999);
434 }
435 irad++;
436 if (irad<11){
437 float radius = 50.0 + irad*15.0;
438 candidates = hel.pathLength(radius);
439 sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
440 }
441 }
442 for (int i = irad; i<11; i++)
443 {
444 tmpPosSample[i] = AliFemtoThreeVector(-9999.,-9999.,-9999);
445 }
446
447 static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
448 104,109.2,114.4,119.6,127.195,129.195,131.195,
449 133.195,135.195,137.195,139.195,141.195,
450 143.195,145.195,147.195,149.195,151.195,
451 153.195,155.195,157.195,159.195,161.195,
452 163.195,165.195,167.195,169.195,171.195,
453 173.195,175.195,177.195,179.195,181.195,
454 183.195,185.195,187.195,189.195};
455 int tRow,tSect,tOutOfBound;
456 double tLength,tPhi;
457 float tU;
458 AliFemtoThreeVector tPoint;
459 AliFmThreeVectorD tn(0,0,0);
460 AliFmThreeVectorD tr(0,0,0);
461 int ti =0;
462 // test to enter the loop
463 candidates = hel.pathLength(tRowRadius[ti]);
464 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
465 if (IsNaN(tLength)){
466 cout <<"tLength Init tmp NAN" << endl;
467 cout <<"padrow number= "<<ti << "not reached" << endl;
468 cout << "*** DO NOT ENTER THE LOOP***" << endl;
469 tmpSect[ti]=-1;//sector
470 }
471 // end test
472 while(ti<45 && !IsNaN(tLength)){
473 candidates = hel.pathLength(tRowRadius[ti]);
474 tLength = (candidates.first > 0) ? candidates.first : candidates.second;
475 if (IsNaN(tLength)){
476 cout <<"tLength loop 1st NAN" << endl;
477 cout <<"padrow number= " << ti << " not reached" << endl;
478 cout << "*** THIS IS AN ERROR SHOULDN'T LOOP ***" << endl;
479 tmpSect[ti]=-1;//sector
480 }
481 tPoint = hel.at(tLength);
482 // Find which sector it is on
483 TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
484 if (IsNaN(tmpSect[ti])){
485 cout <<"***ERROR tmpSect"<< endl;
486 }
487 if (IsNaN(tRow)){
488 cout <<"***ERROR tRow"<< endl;
489 }
490 if (IsNaN(tU)){
491 cout <<"***ERROR tU"<< endl;
492 }
493 if (IsNaN(tPhi)){
494 cout <<"***ERROR tPhi"<< endl;
495 }
496 // calculate crossing plane
497 tn.setX(cos(tPhi));
498 tn.setY(sin(tPhi));
499 tr.setX(tRowRadius[ti]*cos(tPhi));
500 tr.setY(tRowRadius[ti]*sin(tPhi));
501 // find crossing point
502 tLength = hel.pathLength(tr,tn);
503 if (IsNaN(tLength)){
504 cout <<"tLength loop 2nd NAN" << endl;
505 cout <<"padrow number= " << ti << " not reached" << endl;
506 tmpSect[ti]=-2;//sector
507 }
508 tPoint = hel.at(tLength);
509 tmpZ[ti] = tPoint.z();
510 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
511 if (IsNaN(tSect)){
512 cout <<"***ERROR tSect 2"<< endl;
513 }
514 if (IsNaN(tRow)){
515 cout <<"***ERROR tRow 2"<< endl;
516 }
517 if (IsNaN(tmpU[ti])){
518 cout <<"***ERROR tmpU[ti] 2"<< endl;
519 }
520 if (IsNaN(tPhi)){
521 cout <<"***ERROR tPhi 2 "<< endl;
522 }
523 if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
524 tmpSect[ti]=-2;
525 // cout << "missed once"<< endl;
526 }
527 else{
528 if(tmpSect[ti] != tSect){
529 // Try again on the other sector
530 tn.setX(cos(tPhi));
531 tn.setY(sin(tPhi));
532 tr.setX(tRowRadius[ti]*cos(tPhi));
533 tr.setY(tRowRadius[ti]*sin(tPhi));
534 // find crossing point
535 tLength = hel.pathLength(tr,tn);
536 tPoint = hel.at(tLength);
537 if (IsNaN(tLength)){
538 cout <<"tLength loop 3rd NAN" << endl;
539 cout <<"padrow number= "<< ti << " not reached" << endl;
540 tmpSect[ti]=-1;//sector
541 }
542 tmpZ[ti] = tPoint.z();
543 tmpSect[ti] = tSect;
544 tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
545 if (IsNaN(tSect)){
546 cout <<"***ERROR tSect 3"<< endl;
547 }
548 if (IsNaN(tRow)){
549 cout <<"***ERROR tRow 3"<< endl;
550 }
551 if (IsNaN(tmpU[ti])){
552 cout <<"***ERROR tmpU[ti] 3"<< endl;
553 }
554 if (IsNaN(tPhi)){
555 cout <<"***ERROR tPhi 3 "<< endl;
556 }
557 if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
558 tmpSect[ti]=-1;
559 }
560 }
561 }
562 if (IsNaN(tmpSect[ti])){
563 cout << "*******************ERROR***************************" << endl;
564 cout <<"AliFemtoParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
565 cout << "*******************ERROR***************************" << endl;
566 }
567 if (IsNaN(tmpU[ti])){
568 cout << "*******************ERROR***************************" << endl;
569 cout <<"AliFemtoParticle--Fctn tmpU=" << tmpU[ti] << endl;
570 cout << "*******************ERROR***************************" << endl;
571 }
572 if (IsNaN(tmpZ[ti])){
573 cout << "*******************ERROR***************************" << endl;
574 cout <<"AliFemtoParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
575 cout << "*******************ERROR***************************" << endl;
576 }
577 // If padrow ti not reached all other beyond are not reached
578 // in this case set sector to -1
579 if (tmpSect[ti]==-1){
580 for (int tj=ti; tj<45;tj++){
581 tmpSect[tj] = -1;
582 ti=45;
583 }
584 }
585 ti++;
586 if (ti<45){
587 candidates = hel.pathLength(tRowRadius[ti]);
588 tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
589 }
590}
591//_____________________
592const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosExitPoint() const{
593 return fTpcV0PosExitPoint;
594}
595//_____________________
596const AliFemtoThreeVector& AliFemtoParticle::TpcV0PosEntrancePoint() const{
597 return fTpcV0PosEntrancePoint;
598}
599//______________________
600const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegExitPoint() const{
601 return fTpcV0NegExitPoint;
602}
603//_____________________
604const AliFemtoThreeVector& AliFemtoParticle::TpcV0NegEntrancePoint() const{
605 return fTpcV0NegEntrancePoint;
606}
607//______________________