PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFemtoPair.cxx

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