StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtPair.cc
1 /***************************************************************************
2  *
3  * StHbtPair.cc,v 1.23
4  *
5  * Author: Brian Laziuk, Yale University
6  * slightly modified by Mike Lisa
7  ***************************************************************************
8  *
9  * Description: part of STAR HBT Framework: StHbtMaker package
10  * the Pair object is passed to the PairCuts for verification, and
11  * then to the AddRealPair and AddMixedPair methods of the
12  * Correlation Functions
13  *
14  ***************************************************************************
15  * Revision 1.23 2002/09/25 19:23:25 rcwells
16  * Added const to emissionAngle()
17  *
18  * Revision 1.22 2002/04/22 22:48:11 laue
19  * corrected calculation of opening angle
20  **
21  * $Log: StHbtPair.cc,v $
22  * Revision 1.29 2009/09/23 04:43:41 fine
23  * Fix StLorentxVector ctor
24  *
25  * Revision 1.28 2009/09/23 00:51:21 jeromel
26  * Fix for StThreevector
27  *
28  * Revision 1.27 2003/09/02 17:58:32 perev
29  * gcc 3.2 updates + WarnOff
30  *
31  * Revision 1.26 2003/01/31 19:57:15 magestro
32  * Cleared up simple compiler warnings on i386_linux24
33  *
34  * Revision 1.25 2003/01/14 09:44:08 renault
35  * corrections on average separation calculation for tracks which doesn't cross
36  * all 45 padrows.
37  *
38  * Revision 1.24 2002/11/19 23:33:10 renault
39  * Enable average separation calculation for all combinaisons of
40  * V0 daughters and tracks
41  *
42  * Revision 1.21 2002/02/28 14:18:36 rcwells
43  * Added emissionAngle function to StHbtPair
44  *
45  * Revision 1.20 2001/12/14 23:11:30 fretiere
46  * 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
47  *
48  * Revision 1.19 2001/04/25 18:05:09 perev
49  * HPcorrs
50  *
51  * Revision 1.18 2001/04/03 21:04:36 kisiel
52  *
53  *
54  * Changes needed to make the Theoretical code
55  * work. The main code is the ThCorrFctn directory.
56  * The most visible change is the addition of the
57  * HiddenInfo to StHbtPair.
58  *
59  * Revision 1.17 2001/03/28 22:35:20 flierl
60  * changes and bugfixes in qYKP*
61  * add pairrapidity
62  *
63  * Revision 1.16 2001/02/15 19:23:00 rcwells
64  * Fixed sign in qSideCMS
65  *
66  * Revision 1.15 2001/01/22 22:56:41 laue
67  * Yano-Koonin-Podgoretskii Parametrisation added
68  *
69  * Revision 1.14 2000/12/11 21:44:30 rcwells
70  * Corrected qSideCMS function
71  *
72  * Revision 1.13 2000/10/26 16:09:16 lisa
73  * Added OpeningAngle PairCut class and method to StHbtPair
74  *
75  * Revision 1.12 2000/10/05 23:09:05 lisa
76  * Added kT-dependent radii to mixed-event simulator AND implemented AverageSeparation Cut and CorrFctn
77  *
78  * Revision 1.11 2000/07/17 20:03:16 lisa
79  * Implemented tools for addressing and assessing trackmerging
80  *
81  * Revision 1.10 2000/04/04 16:27:03 rcwells
82  * Removed an errant cout in StHbtPair.cc
83  *
84  * Revision 1.9 2000/04/04 16:13:09 lisa
85  * StHbtPair:quality() now returns normalized value (and so is double) and add a CorrFctn which looks at quality()
86  *
87  * Revision 1.8 2000/04/03 22:09:12 rcwells
88  * Add member function ... quality().
89  *
90  * Revision 1.7 2000/02/13 21:13:33 lisa
91  * changed ambiguous StHbtPair::fourMomentum() to fourMomentumSum() and fourMomentumDiff() and fixed related bug in QvecCorrFctn
92  *
93  * Revision 1.6 1999/07/29 16:16:34 lisa
94  * Selemons upgrade of StHbtPair class
95  *
96  * Revision 1.5 1999/07/22 18:49:10 lisa
97  * Implement idea of Fabrice to not create and delete StHbtPair all the time
98  *
99  * Revision 1.4 1999/07/12 18:57:05 lisa
100  * fixed small bug in fourMomentum method of StHbtPair
101  *
102  * Revision 1.3 1999/07/06 22:33:22 lisa
103  * Adjusted all to work in pro and new - dev itself is broken
104  *
105  * Revision 1.2 1999/06/29 17:50:27 fisyak
106  * formal changes to account new StEvent, does not complie yet
107  *
108  * Revision 1.1.1.1 1999/06/29 16:02:57 lisa
109  * Installation of StHbtMaker
110  *
111  **************************************************************************/
112 
113 #include "StHbtMaker/Infrastructure/StHbtPair.hh"
114 
115 double StHbtPair::mMaxDuInner = .8;
116 double StHbtPair::mMaxDzInner = 3.;
117 double StHbtPair::mMaxDuOuter = 1.4;
118 double StHbtPair::mMaxDzOuter = 3.2;
119 
120 
121 StHbtPair::StHbtPair(){
122  mTrack1 = 0;
123  mTrack2 = 0;
124  setDefaultHalfFieldMergingPar();
125 }
126 
127 StHbtPair::StHbtPair(StHbtParticle* a, StHbtParticle* b)
128  : mTrack1(a), mTrack2(b)
129 {
130  setDefaultHalfFieldMergingPar();
131 }
132 
133 void StHbtPair::setDefaultHalfFieldMergingPar(){
134  mMaxDuInner = 3;
135  mMaxDzInner = 4.;
136  mMaxDuOuter = 4.;
137  mMaxDzOuter = 6.;
138 }
139 void StHbtPair::setDefaultFullFieldMergingPar(){
140  mMaxDuInner = 0.8;
141  mMaxDzInner = 3.;
142  mMaxDuOuter = 1.4;
143  mMaxDzOuter = 3.2;
144 }
145 void StHbtPair::setMergingPar(double aMaxDuInner, double aMaxDzInner,
146  double aMaxDuOuter, double aMaxDzOuter){
147  mMaxDuInner = aMaxDuInner;
148  mMaxDzInner = aMaxDzInner;
149  mMaxDuOuter = aMaxDuOuter;
150  mMaxDzOuter = aMaxDzOuter;
151 };
152 
153 StHbtPair::~StHbtPair() {/* no-op */}
154 
155 //StHbtPair::StHbtPair(const StHbtPair &a) {/* missing */}
156 
157 //StHbtPair& StHbtPair::operator=(const StHbtPair &a)
158 
159 //_________________
160 double StHbtPair::mInv() const
161 {
162  double InvariantMass = abs(mTrack1->FourMomentum() + mTrack2->FourMomentum());
163  return (InvariantMass);
164 }
165 //_________________
166 double StHbtPair::kT() const
167 {
168 
169  double tmp =
170  (mTrack1->FourMomentum() + mTrack2->FourMomentum()).perp();
171  tmp *= .5;
172 
173  return (tmp);
174 }
175 //_________________
176 double StHbtPair::rap() const
177 {
178  // longitudinal pair rapidity : Y = 0.5 ::log( E1 + E2 + pz1 + pz2 / E1 + E2 - pz1 - pz2 )
179  double tmp = 0.5 * log (
180  (mTrack1->FourMomentum().e() + mTrack2->FourMomentum().e() + mTrack1->FourMomentum().z() + mTrack2->FourMomentum().z()) /
181  (mTrack1->FourMomentum().e() + mTrack2->FourMomentum().e() - mTrack1->FourMomentum().z() - mTrack2->FourMomentum().z())
182  ) ;
183  return (tmp);
184 }
185 //_________________
186 double StHbtPair::emissionAngle()const {
187  double pxTotal = this->fourMomentumSum().x();
188  double pyTotal = this->fourMomentumSum().y();
189  double angle = atan2(pyTotal,pxTotal)*180.0/3.1415926536;
190  if (angle<0.0) angle+=360.0;
191  return angle;
192 }
193 //_________________
194 // get rid of ambiguously-named method fourMomentum() and replace it with
195 // fourMomentumSum() and fourMomentumDiff() - mal 13feb2000
196 StHbtLorentzVector StHbtPair::fourMomentumSum() const
197 {
198  StHbtLorentzVector temp = mTrack1->FourMomentum()+mTrack2->FourMomentum();
199  return temp;
200 }
201 StHbtLorentzVector StHbtPair::fourMomentumDiff() const
202 {
203  StHbtLorentzVector temp = mTrack1->FourMomentum()-mTrack2->FourMomentum();
204  return temp;
205 }
206 //__________________________________
207 // Yano-Koonin-Podgoretskii Parametrisation in CMS
208 void StHbtPair::qYKPCMS(double& qP, double& qT, double& q0) const
209 {
211  // calculate momentum difference in source rest frame (= lab frame)
213  StHbtLorentzVector l1 = mTrack1->FourMomentum() ;
214  StHbtLorentzVector l2 = mTrack2->FourMomentum() ;
216  // random ordering of the particles
217  if ( rand()/(double)RAND_MAX > 0.50 )
218  { l = l1-l2 ; }
219  else
220  { l = l2-l1 ; } ;
221  // fill momentum differences into return variables
222  qP = l.z() ;
223  qT = l.vect().perp() ;
224  q0 = l.e() ;
225 }
226 //___________________________________
227 // Yano-Koonin-Podgoretskii Parametrisation in LCMS
228 void StHbtPair::qYKPLCMS(double& qP, double& qT, double& q0) const
229 {
231  // calculate momentum difference in LCMS : frame where pz1 + pz2 = 0
233  StHbtLorentzVector l1 = mTrack1->FourMomentum() ;
234  StHbtLorentzVector l2 = mTrack2->FourMomentum() ;
235  // determine beta to LCMS
236  double beta = (l1.z()+l2.z()) / (l1.e()+l2.e()) ;
237  double beta2 = beta*beta ;
238  // unfortunately STAR Class lib knows only boost(particle) not boost(beta) :(
239  // -> create particle with velocity beta and mass 1.0
240  // actually this is : dummyPz = ::sqrt( (dummyMass*dummyMass*beta2) / (1-beta2) ) ;
241  double dummyPz = ::sqrt( (beta2) / (1-beta2) ) ;
242  // boost in the correct direction
243  if (beta>0.0) { dummyPz = -dummyPz; } ;
244 
245  // create dummy particle
246  StHbtLorentzVector l(0.0, 0.0, dummyPz,0) ;
247  double dummyMass = 1.0 ;
248 
249  l.setZ(dummyPz);
250  l.setE(l.vect().massHypothesis(dummyMass) );
251 
252  // boost particles along the beam into a frame with velocity beta
253  StHbtLorentzVector l1boosted = l1.boost(l) ;
254  StHbtLorentzVector l2boosted = l2.boost(l) ;
255  // caculate the momentum difference with random ordering of the particle
256  if ( rand()/(double)RAND_MAX >0.50)
257  { l = l1boosted-l2boosted ; }
258  else
259  { l = l2boosted-l1boosted ;} ;
260  // fill momentum differences into return variables
261  qP = l.z() ;
262  qT = l.vect().perp() ;
263  q0 = l.e() ;
264 }
265 //___________________________________
266 // Yano-Koonin-Podgoretskii Parametrisation in pair rest frame
267 void StHbtPair::qYKPPF(double& qP, double& qT, double& q0) const
268 {
270  // calculate momentum difference in pair rest frame : frame where (pz1 + pz2, py1 + py2, px1 + px2) = (0,0,0)
272  StHbtLorentzVector l1 = mTrack1->FourMomentum() ;
273  StHbtLorentzVector l2 = mTrack2->FourMomentum() ;
274  // the center of gravity of the pair travels with l
275  StHbtLorentzVector l = l1 + l2 ;
276  l = -l ;
277  l.setE(-l.e()) ;
278  // boost particles
279  StHbtLorentzVector l1boosted = l1.boost(l) ;
280  StHbtLorentzVector l2boosted = l2.boost(l) ;
281  // caculate the momentum difference with random ordering of the particle
282  if ( rand()/(double)RAND_MAX > 0.50)
283  { l = l1boosted-l2boosted ; }
284  else
285  { l = l2boosted-l1boosted ;} ;
286  // fill momentum differences into return variables
287  qP = l.z();
288  qT = l.vect().perp();
289  q0 = l.e();
290 }
291 //_________________
292 double StHbtPair::qOutCMS() const
293 {
294  StHbtThreeVector tmp1 = mTrack1->FourMomentum().vect();
295  StHbtThreeVector tmp2 = mTrack2->FourMomentum().vect();
296 
297  double dx = tmp1.x() - tmp2.x();
298  double xt = tmp1.x() + tmp2.x();
299 
300  double dy = tmp1.y() - tmp2.y();
301  double yt = tmp1.y() + tmp2.y();
302 
303  double k1 = (::sqrt(xt*xt+yt*yt));
304  double k2 = (dx*xt+dy*yt);
305  double tmp = k2/k1;
306  return (tmp);
307 }
308 //_________________
309 double StHbtPair::qSideCMS() const
310 {
311  StHbtThreeVector tmp1 = mTrack1->FourMomentum().vect();
312  StHbtThreeVector tmp2 = mTrack2->FourMomentum().vect();
313 
314  double x1 = tmp1.x(); double y1 = tmp1.y();
315  double x2 = tmp2.x(); double y2 = tmp2.y();
316 
317  double xt = x1+x2; double yt = y1+y2;
318  double k1 = ::sqrt(xt*xt+yt*yt);
319 
320  double tmp = 2.0*(x2*y1-x1*y2)/k1;
321  return (tmp);
322 }
323 
324 //_________________________
325 double StHbtPair::qLongCMS() const
326 {
327  StHbtLorentzVector tmp1 = mTrack1->FourMomentum();
328  StHbtLorentzVector tmp2 = mTrack2->FourMomentum();
329 
330  double dz = tmp1.z() - tmp2.z();
331  double zz = tmp1.z() + tmp2.z();
332 
333  double dt = tmp1.t() - tmp2.t();
334  double tt = tmp1.t() + tmp2.t();
335 
336  double beta = zz/tt;
337  double gamma = 1.0/::sqrt(1.0 - beta*beta);
338 
339  double temp = gamma*(dz - beta*dt);
340  return (temp);
341 }
342 
343 //________________________________
344 double StHbtPair::qOutPf() const
345 {
346  StHbtLorentzVector tmp1 = mTrack1->FourMomentum();
347  StHbtLorentzVector tmp2 = mTrack2->FourMomentum();
348 
349  double dt = tmp1.t() - tmp2.t();
350  double tt = tmp1.t() + tmp2.t();
351 
352  double xt = tmp1.x() + tmp2.x();
353  double yt = tmp1.y() + tmp2.y();
354 
355  double k1 = ::sqrt(xt*xt + yt*yt);
356  double bOut = k1/tt;
357  double gOut = 1.0/::sqrt(1.0 - bOut*bOut);
358 
359  double temp = gOut*(this->qOutCMS() - bOut*dt);
360  return (temp);
361 }
362 
363 //___________________________________
364 double StHbtPair::qSidePf() const
365 {
366  return(this->qSideCMS());
367 }
368 
369 //___________________________________
370 
371 double StHbtPair::qLongPf() const
372 {
373  return(this->qLongCMS());
374 }
375 
376 //___________________________________
377 double StHbtPair::qOutBf(double beta) const
378 {
379  return(this->qOutCMS());
380 }
381 
382 //___________________________________
383 
384 double StHbtPair::qSideBf(double beta) const
385 {
386  return(this->qSideCMS());
387 }
388 
389 //___________________________________
390 double StHbtPair::qLongBf(double beta) const
391 {
392  StHbtLorentzVector tmp1 = mTrack1->FourMomentum();
393  StHbtLorentzVector tmp2 = mTrack2->FourMomentum();
394 
395  double dz = tmp1.z() - tmp2.z();
396  double dt = tmp1.t() + tmp2.t();
397 
398  double gamma = 1.0/::sqrt(1.0 - beta*beta);
399 
400  double temp = gamma*(dz - beta*dt);
401  return (temp);
402 }
403 
404 double StHbtPair::quality() const {
405  unsigned long mapMask0 = 0xFFFFFF00;
406  unsigned long mapMask1 = 0x1FFFFF;
407  unsigned long padRow1To24Track1 = mTrack1->TopologyMap(0) & mapMask0;
408  unsigned long padRow25To45Track1 = mTrack1->TopologyMap(1) & mapMask1;
409  unsigned long padRow1To24Track2 = mTrack2->TopologyMap(0) & mapMask0;
410  unsigned long padRow25To45Track2 = mTrack2->TopologyMap(1) & mapMask1;
411  // AND logic
412  unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
413  unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
414  // XOR logic
415  unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
416  unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
417  unsigned long bitI;
418  int ibits;
419  int Quality = 0;
420  double normQual = 0.0;
421  int MaxQuality = mTrack1->NumberOfHits() + mTrack2->NumberOfHits();
422  for (ibits=8;ibits<=31;ibits++) {
423  bitI = 0;
424  bitI |= 1UL<<(ibits);
425  if ( onePad1To24 & bitI ) {
426  Quality++;
427  continue;
428  }
429  else{
430  if ( bothPads1To24 & bitI ) Quality--;
431  }
432  }
433  for (ibits=0;ibits<=20;ibits++) {
434  bitI = 0;
435  bitI |= 1UL<<(ibits);
436  if ( onePad25To45 & bitI ) {
437  Quality++;
438  continue;
439  }
440  else{
441  if ( bothPads25To45 & bitI ) Quality--;
442  }
443  }
444  normQual = (double)Quality/( (double) MaxQuality );
445  return ( normQual );
446 
447 }
448 
449 double StHbtPair::quality2() const {
450  unsigned long mapMask0 = 0xFFFFFF00;
451  unsigned long mapMask1 = 0x1FFFFF;
452  unsigned long padRow1To24Track1 = mTrack1->TopologyMap(0) & mapMask0;
453  unsigned long padRow25To45Track1 = mTrack1->TopologyMap(1) & mapMask1;
454  unsigned long padRow1To24Track2 = mTrack2->TopologyMap(0) & mapMask0;
455  unsigned long padRow25To45Track2 = mTrack2->TopologyMap(1) & mapMask1;
456 
457  // AND logic
458  //unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
459  //unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
460 
461  // XOR logic
462  unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
463  unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
464  unsigned long bitI;
465  int ibits;
466  int Quality = 0;
467  double normQual = 0.0;
468  int MaxQuality = mTrack1->NumberOfHits() + mTrack2->NumberOfHits();
469  for (ibits=8;ibits<=31;ibits++) {
470  bitI = 0;
471  bitI |= 1UL<<(ibits);
472  if ( onePad1To24 & bitI ) {
473  Quality++;
474  continue;
475  }
476  //else{
477  //if ( bothPads1To24 & bitI ) Quality--;
478  //}
479  }
480  for (ibits=0;ibits<=20;ibits++) {
481  bitI = 0;
482  bitI |= 1UL<<(ibits);
483  if ( onePad25To45 & bitI ) {
484  Quality++;
485  continue;
486  }
487  //else{
488  //if ( bothPads25To45 & bitI ) Quality--;
489  //}
490  }
491  normQual = (double)Quality/( (double) MaxQuality );
492  return ( normQual );
493 
494 }
495 
496 
497 double StHbtPair::NominalTpcExitSeparation() const {
498  StHbtThreeVector diff = mTrack1->NominalTpcExitPoint() - mTrack2->NominalTpcExitPoint();
499  return (diff.mag());
500 }
501 
502 double StHbtPair::NominalTpcEntranceSeparation() const {
503  StHbtThreeVector diff = mTrack1->NominalTpcEntrancePoint() - mTrack2->NominalTpcEntrancePoint();
504  return (diff.mag());
505 }
506 
507 double StHbtPair::NominalTpcAverageSeparation() const {
508  StHbtThreeVector diff;
509  double AveSep = 0.0;
510  int ipt = 0;
511  if (mTrack1->mNominalPosSample && mTrack2->mNominalPosSample){
512  while (fabs(mTrack1->mNominalPosSample[ipt].x())<9999. &&
513  fabs(mTrack1->mNominalPosSample[ipt].y())<9999. &&
514  fabs(mTrack1->mNominalPosSample[ipt].z())<9999. &&
515  fabs(mTrack2->mNominalPosSample[ipt].x())<9999. &&
516  fabs(mTrack2->mNominalPosSample[ipt].y())<9999. &&
517  fabs(mTrack2->mNominalPosSample[ipt].z())<9999. &&
518  ipt<11
519  ){
520  // for (int ipt=0; ipt<11; ipt++){
521  diff = mTrack1->mNominalPosSample[ipt] - mTrack2->mNominalPosSample[ipt];
522  ipt++;
523  AveSep += diff.mag();
524  }
525  AveSep = AveSep/(ipt+1.);
526  return (AveSep);}
527  else return -1;
528 }
529 
530 double StHbtPair::OpeningAngle() const {
531  return 57.296* mTrack1->FourMomentum().vect().angle( mTrack2->FourMomentum().vect() );
532 // StHbtThreeVector p1 = mTrack1->FourMomentum().vect();
533 // StHbtThreeVector p2 = mTrack2->FourMomentum().vect();
534 // return 57.296*(p1.phi()-p2.phi());
535 // //double dAngInv = 57.296*acos((p1.dot(p2))/(p1.mag()*p2.mag()));
536 // //return (dAngInv);
537 }
538 //_________________
539 
540 
541 double StHbtPair::KStarFlipped() const {
542  StHbtLorentzVector tP1 = mTrack1->FourMomentum();
543 
544  StThreeVectorD qwe = tP1.vect();
545  qwe *= -1.; // flip it
546  tP1.setVect(qwe);
547 
548  StHbtLorentzVector tSum = (tP1+mTrack2->FourMomentum());
549  double tMass = abs(tSum);
550  StThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
551  double tGamma = tSum.e()/tMass;
552  StThreeVectorD tLongMom = ((tP1.vect()*tGammaBeta)/
553  (tGammaBeta*tGammaBeta))*tGammaBeta;
554  StLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
555  tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
556 //VP tP1.vect() *= -1.; // unflip it
557  return tK.vect().mag();
558 }
559 
560 //double StHbtPair::CVK() const{
561 //const StHbtLorentzVector& tP1 = mTrack1->FourMomentum();
562 //StHbtLorentzVector tSum = (tP1+mTrack2->FourMomentum());
563 //double tMass = abs(tSum);
564 //StThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
565 //double tGamma = tSum.e()/tMass;
566 //StThreeVectorD tLongMom = ((tP1.vect()*tGammaBeta)/
567 // (tGammaBeta*tGammaBeta))*tGammaBeta;
568 //StLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
569 // tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
570 //return (tK.vect())*tGammaBeta/tK.vect().magnitude()/tGammaBeta.magnitude();
571 //}
572 
573 double StHbtPair::CVKFlipped() const{
574  StHbtLorentzVector tP1 = mTrack1->FourMomentum();
575  StThreeVectorD qwe = tP1.vect();
576  qwe *= -1.; // flip it
577  tP1.setVect(qwe);
578 
579  StHbtLorentzVector tSum = (tP1+mTrack2->FourMomentum());
580  double tMass = abs(tSum);
581  StThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
582  double tGamma = tSum.e()/tMass;
583  StThreeVectorD tLongMom = ((tP1.vect()*tGammaBeta)/
584  (tGammaBeta*tGammaBeta))*tGammaBeta;
585  StLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
586  tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
587 //VP tP1.vect() *= -1.; // unflip it
588  return (tK.vect())*tGammaBeta/tGamma;
589 }
590 
591 double StHbtPair::pInv() const{
592  StHbtLorentzVector tP1 = mTrack1->FourMomentum();
593  StHbtLorentzVector tP2 = mTrack2->FourMomentum();
594  double tP = (tP1.px()+tP2.px())*(tP1.px()+tP2.px())+
595  (tP1.py()+tP2.py())*(tP1.py()+tP2.py())+
596  (tP1.pz()+tP2.pz())*(tP1.pz()+tP2.pz())-
597  (tP1.e() -tP2.e() )*(tP1.e() -tP2.e() );
598  return ::sqrt(fabs(tP));
599 }
600 
601 double StHbtPair::qInvFlippedXY() const{
602  StHbtLorentzVector tP1 = mTrack1->FourMomentum();
603  tP1.setX(-1.*tP1.x());
604  tP1.setY(-1.*tP1.y());
605  StHbtLorentzVector tDiff = (tP1-mTrack2->FourMomentum());
606  return ( -1.* tDiff.m());
607 }
608 
609 void StHbtPair::calcNonIdPar() const{ // fortran like function! faster?
610  mNonIdParNotCalculated=0;
611  double px1 = mTrack1->FourMomentum().vect().x();
612  double py1 = mTrack1->FourMomentum().vect().y();
613  double pz1 = mTrack1->FourMomentum().vect().z();
614  double pE1 = mTrack1->FourMomentum().e();
615  double Particle1Mass = ::sqrt(pE1*pE1 - px1*px1 - py1*py1 - pz1*pz1);
616  double px2 = mTrack2->FourMomentum().vect().x();
617  double py2 = mTrack2->FourMomentum().vect().y();
618  double pz2 = mTrack2->FourMomentum().vect().z();
619  double pE2 = mTrack2->FourMomentum().e();
620  double Particle2Mass = ::sqrt(pE2*pE2 - px2*px2 - py2*py2 - pz2*pz2);
621 
622  double Px = px1+px2;
623  double Py = py1+py2;
624  double Pz = pz1+pz2;
625  double PE = pE1+pE2;
626 
627  double Ptrans = Px*Px + Py*Py;
628  double Mtrans = PE*PE - Pz*Pz;
629  double Pinv = ::sqrt(Mtrans - Ptrans);
630  Mtrans = ::sqrt(Mtrans);
631  Ptrans = ::sqrt(Ptrans);
632 
633  double QinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
634  (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
635 
636  double Q = (Particle1Mass*Particle1Mass - Particle2Mass*Particle2Mass)/Pinv;
637  Q = sqrt ( Q*Q - QinvL);
638 
639  kStarCalc = Q/2;
640 
641  // ad 1) go to LCMS
642  double beta = Pz/PE;
643  double gamma = PE/Mtrans;
644 
645  double pz1L = gamma * (pz1 - beta * pE1);
646  double pE1L = gamma * (pE1 - beta * pz1);
647 
648  // fill histogram for beam projection ( z - axis )
649  mDKLong = pz1L;
650 
651  // ad 2) rotation px -> Pt
652  double px1R = (px1*Px + py1*Py)/Ptrans;
653  double py1R = (-px1*Py + py1*Px)/Ptrans;
654 
655  //fill histograms for side projection ( y - axis )
656  mDKSide = py1R;
657 
658  // ad 3) go from LCMS to CMS
659  beta = Ptrans/Mtrans;
660  gamma = Mtrans/Pinv;
661 
662  double px1C = gamma * (px1R - beta * pE1L);
663 
664  // fill histogram for out projection ( x - axis )
665  mDKOut = px1C;
666 
667  mCVK = (mDKOut*Ptrans + mDKLong*Pz)/kStarCalc/::sqrt(Ptrans*Ptrans+Pz*Pz);
668 }
669 
670 
671 void StHbtPair::calcNonIdParGlobal() const{ // fortran like function! faster?
672  mNonIdParNotCalculatedGlobal=0;
673  double px1 = mTrack1->Track()->PGlobal().x();
674  double py1 = mTrack1->Track()->PGlobal().y();
675  double pz1 = mTrack1->Track()->PGlobal().z();
676  double Particle1Mass = mTrack1->FourMomentum().m2();
677  double pE1 = ::sqrt(Particle1Mass + px1*px1 + py1*py1 + pz1*pz1);
678  Particle1Mass = ::sqrt(Particle1Mass);
679 
680  double px2 = mTrack2->Track()->PGlobal().x();
681  double py2 = mTrack2->Track()->PGlobal().y();
682  double pz2 = mTrack2->Track()->PGlobal().z();
683  double Particle2Mass = mTrack2->FourMomentum().m2();
684  double pE2 = ::sqrt(Particle2Mass + px2*px2 + py2*py2 + pz2*pz2);
685  Particle2Mass = ::sqrt(Particle2Mass);
686 
687  double Px = px1+px2;
688  double Py = py1+py2;
689  double Pz = pz1+pz2;
690  double PE = pE1+pE2;
691 
692  double Ptrans = Px*Px + Py*Py;
693  double Mtrans = PE*PE - Pz*Pz;
694  double Pinv = ::sqrt(Mtrans - Ptrans);
695  Mtrans = ::sqrt(Mtrans);
696  Ptrans = ::sqrt(Ptrans);
697 
698  double QinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
699  (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
700 
701  double Q = (Particle1Mass*Particle1Mass - Particle2Mass*Particle2Mass)/Pinv;
702  Q = sqrt ( Q*Q - QinvL);
703 
704  kStarCalcGlobal = Q/2;
705 
706  // ad 1) go to LCMS
707  double beta = Pz/PE;
708  double gamma = PE/Mtrans;
709 
710  double pz1L = gamma * (pz1 - beta * pE1);
711  double pE1L = gamma * (pE1 - beta * pz1);
712 
713  // fill histogram for beam projection ( z - axis )
714  mDKLongGlobal = pz1L;
715 
716  // ad 2) rotation px -> Pt
717  double px1R = (px1*Px + py1*Py)/Ptrans;
718  double py1R = (-px1*Py + py1*Px)/Ptrans;
719 
720  //fill histograms for side projection ( y - axis )
721  mDKSideGlobal = py1R;
722 
723  // ad 3) go from LCMS to CMS
724  beta = Ptrans/Mtrans;
725  gamma = Mtrans/Pinv;
726 
727  double px1C = gamma * (px1R - beta * pE1L);
728 
729  // fill histogram for out projection ( x - axis )
730  mDKOutGlobal = px1C;
731 
732  mCVKGlobal = (mDKOutGlobal*Ptrans + mDKLongGlobal*Pz)/
733  kStarCalcGlobal/::sqrt(Ptrans*Ptrans+Pz*Pz);
734 }
735 
736 
737 
738 double StHbtPair::dcaInsideTpc() const{
739 
740  double tMinDist=NominalTpcEntranceSeparation();
741  double tExit = NominalTpcExitSeparation();
742  tMinDist = (tExit>tMinDist) ? tMinDist : tExit;
743  double tInsideDist;
744  //tMinDist = 999.;
745 
746  double rMin = 60.;
747  double rMax = 190.;
748  const StPhysicalHelixD& tHelix1 = mTrack1->Helix();
749  const StPhysicalHelixD& tHelix2 = mTrack2->Helix();
750  // --- One is a line and other one a helix
751  //if (tHelix1.mSingularity != tHelix2.mSingularity) return -999.;
752  // --- 2 lines : don't care right now
753  //if (tHelix1.mSingularity) return -999.;
754  // --- 2 helix
755  double dx = tHelix2.xcenter() - tHelix1.xcenter();
756  double dy = tHelix2.ycenter() - tHelix1.ycenter();
757  double dd = ::sqrt(dx*dx + dy*dy);
758  double r1 = 1/tHelix1.curvature();
759  double r2 = 1/tHelix2.curvature();
760  double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
761 
762  double x, y, r;
763  double s;
764  if (fabs(cosAlpha) < 1) { // two solutions
765  double sinAlpha = sin(acos(cosAlpha));
766  x = tHelix1.xcenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
767  y = tHelix1.ycenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
768  r = ::sqrt(x*x+y*y);
769  if( r > rMin && r < rMax &&
770  fabs(atan2(y,x)-mTrack1->NominalTpcEntrancePoint().phi())< 0.5
771  ){ // first solution inside
772  s = tHelix1.pathLength(x, y);
773  tInsideDist=tHelix2.distance(tHelix1.at(s));
774  if(tInsideDist<tMinDist) tMinDist = tInsideDist;
775  }
776  else{
777  x = tHelix1.xcenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
778  y = tHelix1.ycenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
779  r = ::sqrt(x*x+y*y);
780  if( r > rMin && r < rMax &&
781  fabs(atan2(y,x)-mTrack1->NominalTpcEntrancePoint().phi())< 0.5
782  ) { // second solution inside
783  s = tHelix1.pathLength(x, y);
784  tInsideDist=tHelix2.distance(tHelix1.at(s));
785  if(tInsideDist<tMinDist) tMinDist = tInsideDist;
786  }
787  }
788  }
789  return tMinDist;
790 }
791 
792 void StHbtPair::calcMergingPar() const{
793  mMergingParNotCalculated=0;
794 
795  double tDu, tDz;
796  int tN = 0;
797  mFracOfMergedRow = 0.;
798  mWeightedAvSep =0.;
799  double tDist;
800  double tDistMax = 200.;
801  for(int ti=0 ; ti<45 ; ti++){
802  if(mTrack1->mSect[ti]==mTrack2->mSect[ti] && mTrack1->mSect[ti]!=-1){
803  tDu = fabs(mTrack1->mU[ti]-mTrack2->mU[ti]);
804  tDz = fabs(mTrack1->mZ[ti]-mTrack2->mZ[ti]);
805  tN++;
806  if(ti<13){
807  mFracOfMergedRow += (tDu<mMaxDuInner && tDz<mMaxDzInner);
808  tDist = ::sqrt(tDu*tDu/mMaxDuInner/mMaxDuInner+
809  tDz*tDz/mMaxDzInner/mMaxDzInner);
810  //mFracOfMergedRow += (tDu<mMaxDuInner && tDz<mMaxDzInner);
811  }
812  else{
813  mFracOfMergedRow += (tDu<mMaxDuOuter && tDz<mMaxDzOuter);
814  tDist = ::sqrt(tDu*tDu/mMaxDuOuter/mMaxDuOuter+
815  tDz*tDz/mMaxDzOuter/mMaxDzOuter);
816  //mFracOfMergedRow += (tDu<mMaxDuOuter && tDz<mMaxDzOuter);
817  }
818  if(tDist<tDistMax){
819  mClosestRowAtDCA = ti+1;
820  tDistMax = tDist;
821  }
822  mWeightedAvSep += tDist;
823  }
824  }
825  if(tN>0){
826  mWeightedAvSep /= tN;
827  mFracOfMergedRow /= tN;
828  }
829  else{
830  mClosestRowAtDCA = -1;
831  mFracOfMergedRow = -1.;
832  mWeightedAvSep = -1.;
833  }
834 }
835 //________________V0 daughters exit/entrance/average separation calc.
836 //_______1st part is a track 2nd is a V0 considering Pos daughter
837 double StHbtPair::TpcExitSeparationTrackV0Pos() const {
838  StHbtThreeVector diff = mTrack1->NominalTpcExitPoint() - mTrack2->TpcV0PosExitPoint();
839  return (diff.mag());
840 }
841 
842 double StHbtPair::TpcEntranceSeparationTrackV0Pos() const {
843  StHbtThreeVector diff = mTrack1->NominalTpcEntrancePoint() - mTrack2->TpcV0PosEntrancePoint();
844  return (diff.mag());
845 }
846 
847 double StHbtPair::TpcAverageSeparationTrackV0Pos() const {
848  StHbtThreeVector diff;
849  double AveSep = 0.0;
850  int ipt = 0;
851  if (mTrack1->mNominalPosSample && mTrack2->mNominalPosSample){
852  while (fabs(mTrack1->mNominalPosSample[ipt].x())<9999. &&
853  fabs(mTrack1->mNominalPosSample[ipt].y())<9999. &&
854  fabs(mTrack1->mNominalPosSample[ipt].z())<9999. &&
855  fabs(mTrack2->mNominalPosSample[ipt].x())<9999. &&
856  fabs(mTrack2->mNominalPosSample[ipt].y())<9999. &&
857  fabs(mTrack2->mNominalPosSample[ipt].z())<9999. &&
858  (ipt<11)
859  ){
860  diff = mTrack1->mNominalPosSample[ipt] - mTrack2->mNominalPosSample[ipt];
861  ipt++;
862  AveSep += diff.mag();
863  }
864  AveSep = AveSep/(ipt+1.);
865  return (AveSep);}
866  else return -1;
867 }
868 //_______1st part is a track 2nd is a V0 considering Neg daughter
869 double StHbtPair::TpcExitSeparationTrackV0Neg() const {
870  StHbtThreeVector diff = mTrack1->NominalTpcExitPoint() - mTrack2->TpcV0NegExitPoint();
871  return (diff.mag());
872 }
873 
874 double StHbtPair::TpcEntranceSeparationTrackV0Neg() const {
875  StHbtThreeVector diff = mTrack1->NominalTpcEntrancePoint() - mTrack2->TpcV0NegEntrancePoint();
876  return (diff.mag());
877 }
878 
879 double StHbtPair::TpcAverageSeparationTrackV0Neg() const {
880  StHbtThreeVector diff;
881  double AveSep = 0.0;
882  int ipt = 0;
883  if (mTrack1->mNominalPosSample && mTrack2->mTpcV0NegPosSample){
884  while (fabs(mTrack1->mNominalPosSample[ipt].x())<9999. &&
885  fabs(mTrack1->mNominalPosSample[ipt].y())<9999. &&
886  fabs(mTrack1->mNominalPosSample[ipt].z())<9999. &&
887  fabs(mTrack2->mTpcV0NegPosSample[ipt].x())<9999. &&
888  fabs(mTrack2->mTpcV0NegPosSample[ipt].y())<9999. &&
889  fabs(mTrack2->mTpcV0NegPosSample[ipt].z())<9999. &&
890  (ipt<11)
891  ){
892  diff = mTrack1->mNominalPosSample[ipt] - mTrack2->mTpcV0NegPosSample[ipt];
893  ipt++;
894  AveSep += diff.mag();
895  }
896  AveSep = AveSep/(ipt+1.);
897  return (AveSep);}
898  else return -1;
899 }
900 
901 //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
902 double StHbtPair::TpcExitSeparationV0PosV0Pos() const {
903  StHbtThreeVector diff = mTrack1->TpcV0PosExitPoint() - mTrack2->TpcV0PosExitPoint();
904  return (diff.mag());
905 }
906 
907 double StHbtPair::TpcEntranceSeparationV0PosV0Pos() const {
908  StHbtThreeVector diff = mTrack1->TpcV0PosEntrancePoint() - mTrack2->TpcV0PosEntrancePoint();
909  return (diff.mag());
910 }
911 double StHbtPair::TpcAverageSeparationV0PosV0Pos() const {
912  StHbtThreeVector diff;
913  double AveSep = 0.0;
914  int ipt=0;
915  if (mTrack1->mNominalPosSample && (mTrack2->mNominalPosSample)){
916  while ((fabs(mTrack1->mNominalPosSample[ipt].x())<9999.) &&
917  (fabs(mTrack1->mNominalPosSample[ipt].y())<9999.) &&
918  (fabs(mTrack1->mNominalPosSample[ipt].z())<9999.) &&
919  (fabs(mTrack2->mNominalPosSample[ipt].x())<9999.) &&
920  (fabs(mTrack2->mNominalPosSample[ipt].y())<9999.) &&
921  (fabs(mTrack2->mNominalPosSample[ipt].z())<9999.) &&
922  (ipt<11)
923  ){
924  diff = mTrack1->mNominalPosSample[ipt] - mTrack2->mNominalPosSample[ipt];
925  ipt++;
926  AveSep += diff.mag();
927  }
928  AveSep = AveSep/(ipt+1);
929  return (AveSep);}
930  else return -1;
931 }
932 
933 //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
934 double StHbtPair::TpcExitSeparationV0PosV0Neg() const {
935  StHbtThreeVector diff = mTrack1->TpcV0PosExitPoint() - mTrack2->TpcV0NegExitPoint();
936  return (diff.mag());
937 }
938 
939 double StHbtPair::TpcEntranceSeparationV0PosV0Neg() const {
940  StHbtThreeVector diff = mTrack1->TpcV0PosEntrancePoint() - mTrack2->TpcV0NegEntrancePoint();
941  return (diff.mag());
942 }
943 double StHbtPair::TpcAverageSeparationV0PosV0Neg() const {
944  StHbtThreeVector diff;
945  double AveSep = 0.0;
946  int ipt = 0;
947  if (mTrack1->mNominalPosSample && mTrack2->mTpcV0NegPosSample){
948  while (fabs(mTrack1->mNominalPosSample[ipt].x())<9999. &&
949  fabs(mTrack1->mNominalPosSample[ipt].y())<9999. &&
950  fabs(mTrack1->mNominalPosSample[ipt].z())<9999. &&
951  fabs(mTrack2->mTpcV0NegPosSample[ipt].x())<9999. &&
952  fabs(mTrack2->mTpcV0NegPosSample[ipt].y())<9999. &&
953  fabs(mTrack2->mTpcV0NegPosSample[ipt].z())<9999. &&
954  (ipt<11)
955  ){
956  diff = mTrack1->mNominalPosSample[ipt] - mTrack2->mTpcV0NegPosSample[ipt];
957  ipt++;
958  AveSep += diff.mag();
959  }
960  AveSep = AveSep/(ipt+1.);
961  return (AveSep);}
962  else return -1;
963 }
964 //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
965 // this is to check the upper case
966 double StHbtPair::TpcExitSeparationV0NegV0Pos() const {
967  StHbtThreeVector diff = mTrack1->TpcV0NegExitPoint() - mTrack2->TpcV0PosExitPoint();
968  return (diff.mag());
969 }
970 
971 double StHbtPair::TpcEntranceSeparationV0NegV0Pos() const {
972  StHbtThreeVector diff = mTrack1->TpcV0NegEntrancePoint() - mTrack2->TpcV0PosEntrancePoint();
973  return (diff.mag());
974 }
975 double StHbtPair::TpcAverageSeparationV0NegV0Pos() const {
976  StHbtThreeVector diff;
977  double AveSep = 0.0;
978  int ipt = 0;
979  if ( mTrack1->mTpcV0NegPosSample && mTrack2->mNominalPosSample){
980  while (fabs(mTrack1->mTpcV0NegPosSample[ipt].x())<9999. &&
981  fabs(mTrack1->mTpcV0NegPosSample[ipt].y())<9999. &&
982  fabs(mTrack1->mTpcV0NegPosSample[ipt].z())<9999. &&
983  fabs(mTrack2->mNominalPosSample[ipt].x())<9999. &&
984  fabs(mTrack2->mNominalPosSample[ipt].y())<9999. &&
985  fabs(mTrack2->mNominalPosSample[ipt].z())<9999. &&
986  (ipt<11)
987  ){
988  diff = mTrack1->mTpcV0NegPosSample[ipt] - mTrack2->mNominalPosSample[ipt];
989  ipt++;
990  AveSep += diff.mag();
991  }
992  AveSep = AveSep/(ipt+1);
993  return (AveSep);}
994  else return -1;
995 }
996 //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
997 double StHbtPair::TpcExitSeparationV0NegV0Neg() const {
998  StHbtThreeVector diff = mTrack1->TpcV0NegExitPoint() - mTrack2->TpcV0NegExitPoint();
999  return (diff.mag());
1000 }
1001 
1002 double StHbtPair::TpcEntranceSeparationV0NegV0Neg() const {
1003  StHbtThreeVector diff = mTrack1->TpcV0NegEntrancePoint() - mTrack2->TpcV0NegEntrancePoint();
1004  return (diff.mag());
1005 }
1006 double StHbtPair::TpcAverageSeparationV0NegV0Neg() const {
1007  StHbtThreeVector diff;
1008  double AveSep = 0.0;
1009  int ipt=0;
1010  if (mTrack1->mTpcV0NegPosSample && mTrack2->mTpcV0NegPosSample){
1011  while (fabs(mTrack1->mTpcV0NegPosSample[ipt].x())<9999. &&
1012  fabs(mTrack1->mTpcV0NegPosSample[ipt].y())<9999. &&
1013  fabs(mTrack1->mTpcV0NegPosSample[ipt].z())<9999. &&
1014  fabs(mTrack2->mTpcV0NegPosSample[ipt].x())<9999. &&
1015  fabs(mTrack2->mTpcV0NegPosSample[ipt].y())<9999. &&
1016  fabs(mTrack2->mTpcV0NegPosSample[ipt].z())<9999. &&
1017  (ipt<11)
1018  ){
1019  diff = mTrack1->mTpcV0NegPosSample[ipt] - mTrack2->mTpcV0NegPosSample[ipt];
1020  ipt++;
1021  AveSep += diff.mag();
1022  }
1023  AveSep = AveSep/(ipt+1);
1024  return (AveSep);}
1025  else return -1;
1026 }
1027 
1028 //________________end V0 daughters exit/entrance/average separation calc.
1029 void StHbtPair::CalcMergingParFctn(short* tmpMergingParNotCalculatedFctn,
1030  float* tmpZ1,float* tmpU1,
1031  float* tmpZ2,float* tmpU2,
1032  int *tmpSect1,int *tmpSect2,
1033  double* tmpFracOfMergedRow,
1034  double* tmpClosestRowAtDCA
1035  ) const{
1036  tmpMergingParNotCalculatedFctn=0;
1037  double tDu, tDz;
1038  int tN = 0;
1039  *tmpFracOfMergedRow = 0.;
1040  *tmpClosestRowAtDCA = 0.;
1041  double tDist;
1042  double tDistMax = 100000000.;
1043  for(int ti=0 ; ti<45 ; ti++){
1044  if(tmpSect1[ti]==tmpSect2[ti] && tmpSect1[ti]!=-1){
1045  tDu = fabs(tmpU1[ti]-tmpU2[ti]);
1046  tDz = fabs(tmpZ1[ti]-tmpZ2[ti]);
1047  tN++;
1048  if(ti<13){
1049  *tmpFracOfMergedRow += (tDu<mMaxDuInner && tDz<mMaxDzInner);
1050  tDist = ::sqrt(tDu*tDu/mMaxDuInner/mMaxDuInner+
1051  tDz*tDz/mMaxDzInner/mMaxDzInner);
1052  }
1053  else{
1054  *tmpFracOfMergedRow += (tDu<mMaxDuOuter && tDz<mMaxDzOuter);
1055  tDist = ::sqrt(tDu*tDu/mMaxDuOuter/mMaxDuOuter+
1056  tDz*tDz/mMaxDzOuter/mMaxDzOuter);
1057  }
1058  if(tDist<tDistMax){
1059  mClosestRowAtDCA = ti+1;
1060  tDistMax = tDist;
1061  }
1062  //mWeightedAvSep += tDist; // now, wrong but not used
1063  }
1064  }
1065  if(tN>0){
1066  //mWeightedAvSep /= tN;
1067  *tmpFracOfMergedRow /= tN;
1068  }
1069  else{
1070  *tmpClosestRowAtDCA = -1;
1071  *tmpFracOfMergedRow = -1.;
1072  //mWeightedAvSep = -1.;
1073  }
1074 }
1075 
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
double ycenter() const
x-center of circle in xy-plane
Definition: StHelix.cc:215
double xcenter() const
aziumth in xy-plane measured from ring center
Definition: StHelix.cc:207