StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtParticle.cc
1 /***************************************************************************
2  *
3  * $Id: StHbtParticle.cc,v 1.23 2015/11/02 20:12:58 perev Exp $
4  *
5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6  ***************************************************************************
7  *
8  * Description: part of STAR HBT Framework: StHbtMaker 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: StHbtParticle.cc,v $
17  * Revision 1.23 2015/11/02 20:12:58 perev
18  * isnan for new compiler
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 StHbtTrack, StHbtParticle, StHbtHiddenInfo, StHbtPair 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 StHbtKink 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 StHbtPair.
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  * StHbtAnalysis: output changed (a little bit less)
55  * StHbtEvent: new version, members for reference mult added
56  * StHbtIOBinary: new IO for new StHbtEvent version
57  * StHbtTypes: TTree typedef to StHbtTTree added
58  * StHbtVertexAnalysis: 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 StHbtParticle 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  * StHbtCoulomb.cxx StHbtSectoredAnalysis.cxx : updated for standalone version
72  * StHbtV0.cc StHbtV0.hh : some cast to prevent compiling warnings
73  * StHbtParticle.cc StHbtParticle.hh : pointers mTrack,mV0 initialized to 0
74  * StHbtIOBinary.cc : some printouts in #ifdef STHBTDEBUG
75  * StHbtEvent.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  * StHbtEvent, StHbtTrack & StHbtV0 declared friend to StHbtIOBinary
80  * StHbtParticle 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 StHbtParticle 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 StHbt 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 StHbtMaker
103  *
104  **************************************************************************/
105 
106 #include "StHbtMaker/Infrastructure/StHbtParticle.hh"
107 #include "math_constants.h"
108 #ifdef __CC5__
109  #include <math.h>
110 #else
111  #include <cmath>
112 #endif
113 
114 double StHbtParticle::mPrimPimPar0= 9.05632e-01;
115 double StHbtParticle::mPrimPimPar1= -2.26737e-01;
116 double StHbtParticle::mPrimPimPar2= -1.03922e-01;
117 double StHbtParticle::mPrimPipPar0= 9.09616e-01;
118 double StHbtParticle::mPrimPipPar1= -9.00511e-02;
119 double StHbtParticle::mPrimPipPar2= -6.02940e-02;
120 double StHbtParticle::mPrimPmPar0= 0.;
121 double StHbtParticle::mPrimPmPar1= 0.;
122 double StHbtParticle::mPrimPmPar2= 0.;
123 double StHbtParticle::mPrimPpPar0= 0.;
124 double StHbtParticle::mPrimPpPar1= 0.;
125 double StHbtParticle::mPrimPpPar2= 0.;
126 
127 int TpcLocalTransform(StThreeVectorD& xgl,
128  int& iSector,
129  int& iPadrow,
130  float& xlocal,
131  double& ttPhi);
132 
133 
134 //_____________________
135 StHbtParticle::StHbtParticle() : mTrack(0), mV0(0), mKink(0), mHiddenInfo(0) {
136  /* no-op for default */
137 }
138 //_____________________
139 StHbtParticle::~StHbtParticle(){
140  // cout << "Issuing delete for StHbtParticle." << endl;
141 
142  if (mTrack) delete mTrack;
143  if (mV0) {
144  delete[] mTpcV0NegPosSample;
145  delete[] mV0NegZ;
146  delete[] mV0NegU;
147  delete[] mV0NegSect;
148  delete mV0;
149  }
150  if (mKink) delete mKink;
151  // cout << "Trying to delete HiddenInfo: " << mHiddenInfo << endl;
152  if (mHiddenInfo)
153  {
154  // cout << "Deleting HiddenInfo." << endl;
155  delete mHiddenInfo;
156  }
157 }
158 //_____________________
159 StHbtParticle::StHbtParticle(const StHbtTrack* const hbtTrack,const double& mass) : mTrack(0), mV0(0), mKink(0), mHiddenInfo(0) {
160 
161 
162  // I know there is a better way to do this...
163  mTrack = new StHbtTrack(*hbtTrack);
164  StHbtThreeVector temp = hbtTrack->P();
165  mFourMomentum.setVect(temp);
166  double ener = ::sqrt(temp.mag2()+mass*mass);
167  mFourMomentum.setE(ener);
168  mMap[0] = hbtTrack->TopologyMap(0);
169  mMap[1] = hbtTrack->TopologyMap(1);
170  mNhits = hbtTrack->NHits();
171  mHelix = hbtTrack->Helix();
172  //CalculateNominalTpcExitAndEntrancePoints();
173 
174 
175  mPrimaryVertex.setX(0.);
176  mPrimaryVertex.setY(0.);
177  mPrimaryVertex.setZ(0.);
178  mSecondaryVertex.setX(0.);
179  mSecondaryVertex.setY(0.);
180  mSecondaryVertex.setZ(0.);
181 
182  CalculateTpcExitAndEntrancePoints(&mHelix,&mPrimaryVertex,
183  &mSecondaryVertex,
184  &mNominalTpcEntrancePoint,
185  &mNominalTpcExitPoint,
186  &mNominalPosSample[0],
187  &mZ[0],
188  &mU[0],
189  &mSect[0]);
190 
191  CalculatePurity();
192  // ***
193  mHiddenInfo= 0;
194  if(hbtTrack->ValidHiddenInfo()){
195  mHiddenInfo= hbtTrack->getHiddenInfo()->getParticleHiddenInfo();
196  }
197  // ***
198 
199 }
200 //_____________________
201 StHbtParticle::StHbtParticle(const StHbtV0* const hbtV0,const double& mass) : mTrack(0), mV0(0), mKink(0), mHiddenInfo(0) {
202  mV0 = new StHbtV0(*hbtV0);
203  mMap[0]= 0;
204  mMap[1]= 0;
205  // I know there is a better way to do this...
206  StHbtThreeVector temp = hbtV0->momV0();
207  mFourMomentum.setVect(temp);
208  double ener = ::sqrt(temp.mag2()+mass*mass);
209  mFourMomentum.setE(ener);
210  // Calculating TpcEntrancePoint for Positive V0 daugther
211  mPrimaryVertex = hbtV0->primaryVertex();
212  mSecondaryVertex = hbtV0->decayVertexV0();
213  mHelixV0Pos = hbtV0->HelixPos();
214 
215  mTpcV0NegPosSample = new StHbtThreeVector[45];//for V0Neg
216  mV0NegZ = new float[45];//for V0Neg
217  mV0NegU = new float[45];//for V0Neg
218  mV0NegSect = new int[45];//for V0Neg
219 
220  CalculateTpcExitAndEntrancePoints(&mHelixV0Pos,&mPrimaryVertex,
221  &mSecondaryVertex,
222  &mTpcV0PosEntrancePoint,
223  &mTpcV0PosExitPoint,
224  &mNominalPosSample[0],
225  &mZ[0],
226  &mU[0],&mSect[0]);
227 
228  mHelixV0Neg = hbtV0->HelixNeg();
229 
230  CalculateTpcExitAndEntrancePoints(&mHelixV0Neg,
231  &mPrimaryVertex,
232  &mSecondaryVertex,
233  &mTpcV0NegEntrancePoint,
234  &mTpcV0NegExitPoint,
235  &mTpcV0NegPosSample[0],
236  &mV0NegZ[0],
237  &mV0NegU[0],&mV0NegSect[0]);
238 
239  // ***
240  mHiddenInfo= 0;
241  if(hbtV0->ValidHiddenInfo()){
242  mHiddenInfo= hbtV0->getHiddenInfo()->clone();
243  }
244  // ***
245 }
246 //_____________________
247 StHbtParticle::StHbtParticle(const StHbtKink* const hbtKink,const double& mass) : mTrack(0), mV0(0), mHiddenInfo(0) {
248  mKink = new StHbtKink(*hbtKink);
249  mMap[0]= 0;
250  mMap[1]= 0;
251  // I know there is a better way to do this...
252  StHbtThreeVector temp = hbtKink->Parent().P();
253  mFourMomentum.setVect(temp);
254  double ener = ::sqrt(temp.mag2()+mass*mass);
255  mFourMomentum.setE(ener);
256 }
257 
258 //_____________________
259 StHbtParticle::StHbtParticle(const StHbtXi* const hbtXi, const double& mass) {
260  mXi = new StHbtXi(*hbtXi);
261  mMap[0]= 0;
262  mMap[1]= 0;
263  StHbtThreeVector temp;// = hbtXi->mMomXi;
264  mFourMomentum.setVect(temp);
265  double ener = ::sqrt(temp.mag2()+mass*mass);
266  mFourMomentum.setE(ener);
267  mHiddenInfo = 0;
268 }
269 //_____________________
270 const StHbtThreeVector& StHbtParticle::NominalTpcExitPoint() const{
271  // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
272  // for now, we calculate Exit and Entrance points upon instantiation
273  return mNominalTpcExitPoint;
274 }
275 //_____________________
276 const StHbtThreeVector& StHbtParticle::NominalTpcEntrancePoint() const{
277  // in future, may want to calculate this "on demand" only, sot this routine may get more sophisticated
278  // for now, we calculate Exit and Entrance points upon instantiation
279  return mNominalTpcEntrancePoint;
280 }
281 //_____________________
282 //***************WARNING*************************
283 // Gael includes this in a function with parameters to use it for
284 // tracks and V0 daughters GR 5 dec 02
285 //***************WARNING*************************
286 // void StHbtParticle::CalculateNominalTpcExitAndEntrancePoints(){
287 // // this calculates the "nominal" exit point of a track, either through the endcap or through the Outer Field Cage
288 // // "nominal" means the track is assumed to start at (0,0,0)
289 // // it also calculates the "nominal" entrance point of the track, which is the point at which it crosses the
290 // // inner field cage
291 // static StHbtThreeVector ZeroVec(0.,0.,0.);
292 // double dip, curv, phase;
293 // int h;
294 // curv = mHelix.curvature();
295 // dip = mHelix.dipAngle();
296 // phase= mHelix.phase();
297 // h = mHelix.h();
298 // StHelixD hel(curv,dip,phase,ZeroVec,h);
299 
300 // pairD candidates;
301 // double sideLength; // this is how much length to go to leave through sides of TPC
302 // double endLength; // this is how much length to go to leave through endcap of TPC
303 // // figure out how far to go to leave through side...
304 // candidates = hel.pathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
305 // sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
306 
307 // static StHbtThreeVector WestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
308 // static StHbtThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
309 // static StHbtThreeVector EndCapNormal(0.,0.,1.0);
310 
311 // endLength = hel.pathLength(WestEnd,EndCapNormal);
312 // if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
313 
314 // if (endLength < 0.0) cout << "StHbtParticle::CalculateNominalTpcExitAndEntrancePoints(): "
315 // << "Hey-- I cannot find an exit point out endcaps" << endl;
316 
317 // // OK, firstExitLength will be the shortest way out of the detector...
318 // double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
319 
320 // // now then, let's return the POSITION at which particle leaves TPC...
321 // mNominalTpcExitPoint = hel.at(firstExitLength);
322 
323 
324 // // Finally, calculate the position at which the track crosses the inner field cage
325 // candidates = hel.pathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
326 
327 // sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
328 // // if (sideLength < 0.0)
329 // // {
330 // // cout
331 // // << "no crossing with IFC"
332 // // << " curve=" << curv
333 // // << " candidates=" << candidates.first << " " << candidates.second
334 // // << "origin=" << mHelix.origin() << " "<< dip << " " << phase << " " << h << endl;
335 // // }
336 // // else
337 // // {
338 // // cout
339 // // << "does cross IFC"
340 // // << " curve=" << curv
341 // // << " candidates=" << candidates.first << " " << candidates.second
342 // // << "origin=" << mHelix.origin() << " "<< dip << " " << phase << " " << h << endl;
343 // // }
344 
345 
346 // // if (sideLength < 0.0)
347 // // {
348 // // if (phase > C_PI)
349 // // {
350 // // cout << "righto" << endl;
351 // // }
352 // // else
353 // // {
354 // // cout << "WRONGO!! 1 " << phase << endl;
355 // // }
356 // // }
357 // // else
358 // // {
359 // // if (phase > C_PI)
360 // // {
361 // // cout << "WRONGO!! 2 " << phase << endl;
362 // // }
363 // // else
364 // // {
365 // // cout << "righto " << endl;
366 // // }
367 // // }
368 
369 // mNominalTpcEntrancePoint = hel.at(sideLength);
370 
371 
372 // // This is the secure way !
373 // // if (::isnan(mNominalTpcEntrancePoint.x()) ||
374 // // ::isnan(mNominalTpcEntrancePoint.x()) ||
375 // // ::isnan(mNominalTpcEntrancePoint.x()) ) mNominalTpcEntrancePoint = StHbtThreeVector(-9999.,-9999.,-9999);
376 // // if (::isnan(mNominalTpcExitPoint.x()) ||
377 // // ::isnan(mNominalTpcExitPoint.x()) ||
378 // // ::isnan(mNominalTpcExitPoint.x()) ) mNominalTpcExitPoint = StHbtThreeVector(-9999.,-9999.,-9999);
379 
380 // // This is faster
381 // if (::isnan(mNominalTpcExitPoint.x())) mNominalTpcExitPoint = StHbtThreeVector(-9999.,-9999.,-9999);
382 
383 
384 // // 03Oct00 - mal. OK, let's try something a little more along the lines of NA49 and E895 strategy.
385 // // calculate the "nominal" position at N radii (say N=11) within the TPC, and for a pair cut
386 // // use the average separation of these N
387 // for (int irad=0; irad<11; irad++){
388 // float radius = 50.0 + irad*15.0;
389 // candidates = hel.pathLength(radius);
390 // sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
391 // mNominalPosSample[irad] = hel.at(sideLength);
392 // }
393 
394 
395 
396 // static double tSectToPhi[24]={2.,1.,0.,11.,10.,9.,8. ,7. ,6.,5.,4.,3.,
397 // 4.,5.,6., 7., 8.,9.,10.,11.,0.,1.,2.,3.};
398 // static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
399 // 104,109.2,114.4,119.6,127.195,129.195,131.195,
400 // 133.195,135.195,137.195,139.195,141.195,
401 // 143.195,145.195,147.195,149.195,151.195,
402 // 153.195,155.195,157.195,159.195,161.195,
403 // 163.195,165.195,167.195,169.195,171.195,
404 // 173.195,175.195,177.195,179.195,181.195,
405 // 183.195,185.195,187.195,189.195};
406 // int tRow,tSect,tOutOfBound;
407 // double tU,tLength,tPhi;
408 // StHbtThreeVector tPoint;
409 // StThreeVectorD tn(0,0,0);
410 // StThreeVectorD tr(0,0,0);
411 // for(int ti=0;ti<45;ti++){
412 // // Find which sector it is on
413 // candidates = hel.pathLength(tRowRadius[ti]);
414 // tLength = (candidates.first > 0) ? candidates.first : candidates.second;
415 // tPoint = hel.at(tLength);
416 // TpcLocalTransform(tPoint,mSect[ti],tRow,tU,tPhi);
417 // // calculate crossing plane
418 // //tPhi = tSectToPhi[mSect[ti]-1]*TMath::Pi()/6.;
419 // tn.setX(cos(tPhi));
420 // tn.setY(sin(tPhi));
421 // tr.setX(tRowRadius[ti]*cos(tPhi));
422 // tr.setY(tRowRadius[ti]*sin(tPhi));
423 // // find crossing point
424 // tLength = hel.pathLength(tr,tn);
425 // tPoint = hel.at(tLength);
426 // mZ[ti] = tPoint.z();
427 // tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,mU[ti],tPhi);
428 // if(tOutOfBound || (mSect[ti] == tSect && tRow!=(ti+1))){
429 // //cout << "Out of bound " << tOutOfBound2 << " " << tOutOfBound << " "
430 // // << tSect << " " << mSect[ti] << " "
431 // // << ti+1 << " " << tRow << " " << tRowRadius[ti] << " "
432 // // << tU << " " << mU[ti] << endl;
433 // mSect[ti]=-1;
434 // }
435 // else{
436 // if(mSect[ti] != tSect){
437 // // Try again on the other sector
438 // tn.setX(cos(tPhi));
439 // tn.setY(sin(tPhi));
440 // tr.setX(tRowRadius[ti]*cos(tPhi));
441 // tr.setY(tRowRadius[ti]*sin(tPhi));
442 // // find crossing point
443 // tLength = hel.pathLength(tr,tn);
444 // tPoint = hel.at(tLength);
445 // mZ[ti] = tPoint.z();
446 // mSect[ti] = tSect;
447 // tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,mU[ti],tPhi);
448 // if(tOutOfBound || tSect!= mSect[ti] || tRow!=(ti+1)){
449 // mSect[ti]=-1;
450 // //cout << "Twice bad : OutOfBound = " << tOutOfBound
451 // // << " SectOk = " << (tSect!= mSect[ti])
452 // // << " RowOk = " << (tRow!=(ti+1)) << endl;
453 // }
454 // }
455 // }
456 // }
457 // }
458 //_____________________
459 void StHbtParticle::CalculatePurity(){
460  double tPt = mFourMomentum.perp();
461  // pi -
462  mPurity[0] = mPrimPimPar0*(1.-exp((tPt-mPrimPimPar1)/mPrimPimPar2));
463  mPurity[0] *= mTrack->PidProbPion();
464  // pi+
465  mPurity[1] = mPrimPipPar0*(1.-exp((tPt-mPrimPipPar1)/mPrimPipPar2));
466  mPurity[1] *= mTrack->PidProbPion();
467  // K-
468  mPurity[2] = mTrack->PidProbKaon();
469  // K+
470  mPurity[3] = mTrack->PidProbKaon();
471  // pbar
472  mPurity[4] = mTrack->PidProbProton();
473  // p
474  mPurity[5] = mTrack->PidProbProton();
475 }
476 
477 double StHbtParticle::GetPionPurity()
478 {
479  if (mTrack->Charge()>0)
480  return mPurity[1];
481  else
482  return mPurity[0];
483 }
484 double StHbtParticle::GetKaonPurity()
485 {
486  if (mTrack->Charge()>0)
487  return mPurity[3];
488  else
489  return mPurity[2];
490 }
491 double StHbtParticle::GetProtonPurity()
492 {
493  if (mTrack->Charge()>0)
494  return mPurity[5];
495  else
496  return mPurity[4];
497 }
498 
499 void StHbtParticle::CalculateTpcExitAndEntrancePoints(StPhysicalHelixD* tHelix,
500  StHbtThreeVector* PrimVert,
501  StHbtThreeVector* SecVert,
502  StHbtThreeVector* tmpTpcEntrancePoint,
503  StHbtThreeVector* tmpTpcExitPoint,
504  StHbtThreeVector* tmpPosSample,
505  float* tmpZ,
506  float* tmpU,
507  int* tmpSect){
508  // this calculates the exit point of a secondary track,
509  // either through the endcap or through the Outer Field Cage
510  // We assume the track to start at tHelix.origin-PrimaryVertex
511  // it also calculates the entrance point of the secondary track,
512  // which is the point at which it crosses the
513  // inner field cage
514  // static StHbtThreeVector ZeroVec(0.,0.,0.);
515  StHbtThreeVector ZeroVec(0.,0.,0.);
516 // ZeroVec.setX(tHelix->origin().x()-PrimVert->x());
517 // ZeroVec.setY(tHelix->origin().y()-PrimVert->y());
518 // ZeroVec.setZ(tHelix->origin().z()-PrimVert->z());
519  ZeroVec.setX(SecVert->x()-PrimVert->x());
520  ZeroVec.setY(SecVert->y()-PrimVert->y());
521  ZeroVec.setZ(SecVert->z()-PrimVert->z());
522  double dip, curv, phase;
523  int h;
524  curv = tHelix->curvature();
525  dip = tHelix->dipAngle();
526  phase= tHelix->phase();
527  h = tHelix->h();
528 
529  StHelixD hel(curv,dip,phase,ZeroVec,h);
530 
531  pairD candidates;
532  double sideLength; // this is how much length to go to leave through sides of TPC
533  double endLength; // this is how much length to go to leave through endcap of TPC
534  // figure out how far to go to leave through side...
535  candidates = hel.pathLength(200.0); // bugfix MAL jul00 - 200cm NOT 2cm
536  sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
537 
538  static StHbtThreeVector WestEnd(0.,0.,200.); // bugfix MAL jul00 - 200cm NOT 2cm
539  static StHbtThreeVector EastEnd(0.,0.,-200.); // bugfix MAL jul00 - 200cm NOT 2cm
540  static StHbtThreeVector EndCapNormal(0.,0.,1.0);
541 
542  endLength = hel.pathLength(WestEnd,EndCapNormal);
543  if (endLength < 0.0) endLength = hel.pathLength(EastEnd,EndCapNormal);
544 
545  if (endLength < 0.0) cout <<
546  "StHbtParticle::CalculateTpcExitAndEntrancePoints(): "
547  << "Hey -- I cannot find an exit point out endcaps" << endl;
548  // OK, firstExitLength will be the shortest way out of the detector...
549  double firstExitLength = (endLength < sideLength) ? endLength : sideLength;
550  // now then, let's return the POSITION at which particle leaves TPC...
551  *tmpTpcExitPoint = hel.at(firstExitLength);
552  // Finally, calculate the position at which the track crosses the inner field cage
553  candidates = hel.pathLength(50.0); // bugfix MAL jul00 - 200cm NOT 2cm
554 
555  sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
556  // cout << "sideLength 2 ="<<sideLength << endl;
557  *tmpTpcEntrancePoint = hel.at(sideLength);
558  // This is the secure way !
559  if (::isnan(tmpTpcEntrancePoint->x()) ||
560  ::isnan(tmpTpcEntrancePoint->y()) ||
561  ::isnan(tmpTpcEntrancePoint->z()) ){
562  cout << "tmpTpcEntrancePoint NAN"<< endl;
563  cout << "tmpNominalTpcEntrancePoint = " <<tmpTpcEntrancePoint<< endl;
564  tmpTpcEntrancePoint->setX(-9999.);
565  tmpTpcEntrancePoint->setY(-9999.);
566  tmpTpcEntrancePoint->setZ(-9999.);
567  }
568 
569  if (::isnan(tmpTpcExitPoint->x()) ||
570  ::isnan(tmpTpcExitPoint->y()) ||
571  ::isnan(tmpTpcExitPoint->z()) ) {
572 // cout << "tmpTpcExitPoint NAN set at (-9999,-9999,-9999)"<< endl;
573 // cout << "tmpTpcExitPoint X= " <<tmpTpcExitPoint->x()<< endl;
574 // cout << "tmpTpcExitPoint Y= " <<tmpTpcExitPoint->y()<< endl;
575 // cout << "tmpTpcExitPoint Z= " <<tmpTpcExitPoint->z()<< endl;
576  tmpTpcExitPoint->setX(-9999.);
577  tmpTpcExitPoint->setY(-9999.);
578  tmpTpcExitPoint->setZ(-9999.);
579  }
580 
581 
582 // if (::isnan(tmpTpcExitPoint->x())) *tmpTpcExitPoint = StHbtThreeVector(-9999.,-9999.,-9999);
583 // if (::isnan(tmpTpcEntrancetPoint->x())) *tmpTpcEntrancePoint = StHbtThreeVector(-9999.,-9999.,-9999);
584  // cout << "tmpTpcEntrancePoint"<<*tmpTpcEntrancePoint << endl;
585 
586  // 03Oct00 - mal. OK, let's try something a little more
587  // along the lines of NA49 and E895 strategy.
588  // calculate the "nominal" position at N radii (say N=11)
589  // within the TPC, and for a pair cut
590  // use the average separation of these N
591  int irad = 0;
592  candidates = hel.pathLength(50.0);
593  sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
594  while (irad<11 && !::isnan(sideLength)){
595  float radius = 50.0 + irad*15.0;
596  candidates = hel.pathLength(radius);
597  sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
598  tmpPosSample[irad] = hel.at(sideLength);
599  if(::isnan(tmpPosSample[irad].x()) ||
600  ::isnan(tmpPosSample[irad].y()) ||
601  ::isnan(tmpPosSample[irad].z())
602  ){
603  cout << "tmpPosSample for radius=" << radius << " NAN"<< endl;
604  cout << "tmpPosSample=(" <<tmpPosSample[irad]<<")"<< endl;
605  tmpPosSample[irad] = StHbtThreeVector(-9999.,-9999.,-9999);
606  }
607  irad++;
608  if (irad<11){
609  float radius = 50.0 + irad*15.0;
610  candidates = hel.pathLength(radius);
611  sideLength = (candidates.first > 0) ? candidates.first : candidates.second;
612  }
613  }
614  for (int i = irad; i<11; i++)
615  {
616  tmpPosSample[i] = StHbtThreeVector(-9999.,-9999.,-9999);
617  }
618 
619  static float tRowRadius[45] = {60,64.8,69.6,74.4,79.2,84,88.8,93.6,98.8,
620  104,109.2,114.4,119.6,127.195,129.195,131.195,
621  133.195,135.195,137.195,139.195,141.195,
622  143.195,145.195,147.195,149.195,151.195,
623  153.195,155.195,157.195,159.195,161.195,
624  163.195,165.195,167.195,169.195,171.195,
625  173.195,175.195,177.195,179.195,181.195,
626  183.195,185.195,187.195,189.195};
627  int tRow,tSect,tOutOfBound;
628  double tLength,tPhi;
629  float tU;
630  StHbtThreeVector tPoint;
631  StThreeVectorD tn(0,0,0);
632  StThreeVectorD tr(0,0,0);
633  int ti =0;
634  // test to enter the loop
635  candidates = hel.pathLength(tRowRadius[ti]);
636  tLength = (candidates.first > 0) ? candidates.first : candidates.second;
637  if (::isnan(tLength)){
638  cout <<"tLength Init tmp NAN" << endl;
639  cout <<"padrow number= "<<ti << "not reached" << endl;
640  cout << "*** DO NOT ENTER THE LOOP***" << endl;
641  tmpSect[ti]=-1;//sector
642  }
643  // end test
644  while(ti<45 && !::isnan(tLength)){
645  candidates = hel.pathLength(tRowRadius[ti]);
646  tLength = (candidates.first > 0) ? candidates.first : candidates.second;
647  if (::isnan(tLength)){
648  cout <<"tLength loop 1st NAN" << endl;
649  cout <<"padrow number= " << ti << " not reached" << endl;
650  cout << "*** THIS IS AN ERROR SHOULDN'T LOOP ***" << endl;
651  tmpSect[ti]=-1;//sector
652  }
653  tPoint = hel.at(tLength);
654  // Find which sector it is on
655  TpcLocalTransform(tPoint,tmpSect[ti],tRow,tU,tPhi);
656  if (::isnan(tmpSect[ti])){
657  cout <<"***ERROR tmpSect"<< endl;
658  }
659  if (::isnan(tRow)){
660  cout <<"***ERROR tRow"<< endl;
661  }
662  if (::isnan(tU)){
663  cout <<"***ERROR tU"<< endl;
664  }
665  if (::isnan(tPhi)){
666  cout <<"***ERROR tPhi"<< endl;
667  }
668  // calculate crossing plane
669  tn.setX(cos(tPhi));
670  tn.setY(sin(tPhi));
671  tr.setX(tRowRadius[ti]*cos(tPhi));
672  tr.setY(tRowRadius[ti]*sin(tPhi));
673  // find crossing point
674  tLength = hel.pathLength(tr,tn);
675  if (::isnan(tLength)){
676  cout <<"tLength loop 2nd NAN" << endl;
677  cout <<"padrow number= " << ti << " not reached" << endl;
678  tmpSect[ti]=-2;//sector
679  }
680  tPoint = hel.at(tLength);
681  tmpZ[ti] = tPoint.z();
682  tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
683  if (::isnan(tSect)){
684  cout <<"***ERROR tSect 2"<< endl;
685  }
686  if (::isnan(tRow)){
687  cout <<"***ERROR tRow 2"<< endl;
688  }
689  if (::isnan(tmpU[ti])){
690  cout <<"***ERROR tmpU[ti] 2"<< endl;
691  }
692  if (::isnan(tPhi)){
693  cout <<"***ERROR tPhi 2 "<< endl;
694  }
695  if(tOutOfBound || (tmpSect[ti] == tSect && tRow!=(ti+1))){
696  tmpSect[ti]=-2;
697  // cout << "missed once"<< endl;
698  }
699  else{
700  if(tmpSect[ti] != tSect){
701  // Try again on the other sector
702  tn.setX(cos(tPhi));
703  tn.setY(sin(tPhi));
704  tr.setX(tRowRadius[ti]*cos(tPhi));
705  tr.setY(tRowRadius[ti]*sin(tPhi));
706  // find crossing point
707  tLength = hel.pathLength(tr,tn);
708  tPoint = hel.at(tLength);
709  if (::isnan(tLength)){
710  cout <<"tLength loop 3rd NAN" << endl;
711  cout <<"padrow number= "<< ti << " not reached" << endl;
712  tmpSect[ti]=-1;//sector
713  }
714  tmpZ[ti] = tPoint.z();
715  tmpSect[ti] = tSect;
716  tOutOfBound = TpcLocalTransform(tPoint,tSect,tRow,tmpU[ti],tPhi);
717  if (::isnan(tSect)){
718  cout <<"***ERROR tSect 3"<< endl;
719  }
720  if (::isnan(tRow)){
721  cout <<"***ERROR tRow 3"<< endl;
722  }
723  if (::isnan(tmpU[ti])){
724  cout <<"***ERROR tmpU[ti] 3"<< endl;
725  }
726  if (::isnan(tPhi)){
727  cout <<"***ERROR tPhi 3 "<< endl;
728  }
729  if(tOutOfBound || tSect!= tmpSect[ti] || tRow!=(ti+1)){
730  tmpSect[ti]=-1;
731  }
732  }
733  }
734  if (::isnan(tmpSect[ti])){
735  cout << "*******************ERROR***************************" << endl;
736  cout <<"StHbtParticle--Fctn tmpSect=" << tmpSect[ti] << endl;
737  cout << "*******************ERROR***************************" << endl;
738  }
739  if (::isnan(tmpU[ti])){
740  cout << "*******************ERROR***************************" << endl;
741  cout <<"StHbtParticle--Fctn tmpU=" << tmpU[ti] << endl;
742  cout << "*******************ERROR***************************" << endl;
743  }
744  if (::isnan(tmpZ[ti])){
745  cout << "*******************ERROR***************************" << endl;
746  cout <<"StHbtParticle--Fctn tmpZ=" << tmpZ[ti] << endl;
747  cout << "*******************ERROR***************************" << endl;
748  }
749  // If padrow ti not reached all other beyond are not reached
750  // in this case set sector to -1
751  if (tmpSect[ti]==-1){
752  for (int tj=ti; tj<45;tj++){
753  tmpSect[tj] = -1;
754  ti=45;
755  }
756  }
757  ti++;
758  if (ti<45){
759  candidates = hel.pathLength(tRowRadius[ti]);
760  tLength = (candidates.first > 0) ? candidates.first : candidates.second;}
761  }
762 }
763 //_____________________
764 const StHbtThreeVector& StHbtParticle::TpcV0PosExitPoint() const{
765  return mTpcV0PosExitPoint;
766 }
767 //_____________________
768 const StHbtThreeVector& StHbtParticle::TpcV0PosEntrancePoint() const{
769  return mTpcV0PosEntrancePoint;
770 }
771 //______________________
772 const StHbtThreeVector& StHbtParticle::TpcV0NegExitPoint() const{
773  return mTpcV0NegExitPoint;
774 }
775 //_____________________
776 const StHbtThreeVector& StHbtParticle::TpcV0NegEntrancePoint() const{
777  return mTpcV0NegEntrancePoint;
778 }
779 //______________________
StHelix::h
int h() const
y-center of circle in xy-plane
Definition: StHelix.hh:174
StHelix::phase
double phase() const
1/R in xy-plane
Definition: StHelix.hh:180
Generated by   doxygen 1.8.5