PWGCF/FEMTOSCOPY/AliFemto/AliFemtoPair.cxx

   1 ///////////////////////////////////////////////////////////////////////////
   2 //                                                                       //
   3 // AliFemtoPair: the Pair object is passed to the PairCuts for           //
   4 // verification, and then to the AddRealPair and AddMixedPair methods of //
   5 // the Correlation Functions. It holds pair-specific variables like      //
   6 // relative momenta and has links to the particles and tracks that form  //
   7 // the pair.                                                             //
   8 //                                                                       //
   9 ///////////////////////////////////////////////////////////////////////////
  10 #include <TMath.h>
  11 #include "AliFemtoPair.h"
  12
  13 double AliFemtoPair::fgMaxDuInner = .8;
  14 double AliFemtoPair::fgMaxDzInner = 3.;
  15 double AliFemtoPair::fgMaxDuOuter = 1.4;
  16 double AliFemtoPair::fgMaxDzOuter = 3.2;
  17
  18
  19 AliFemtoPair::AliFemtoPair() :
  20   fTrack1(0), fTrack2(0),
  21   fPairAngleEP(0),
  22   fNonIdParNotCalculated(0),
  23   fDKSide(0),
  24   fDKOut(0),
  25   fDKLong(0),
  26   fCVK(0),
  27   fKStarCalc(0),
  28   fNonIdParNotCalculatedGlobal(0),
  29   fMergingParNotCalculated(0),
  30   fWeightedAvSep(0),
  31   fFracOfMergedRow(0),
  32   fClosestRowAtDCA(0),
  33   fMergingParNotCalculatedTrkV0Pos(0),
  34   fFracOfMergedRowTrkV0Pos(0),
  35   fClosestRowAtDCATrkV0Pos(0),
  36   fMergingParNotCalculatedTrkV0Neg(0),
  37   fFracOfMergedRowTrkV0Neg(0),
  38   fClosestRowAtDCATrkV0Neg(0),
  39   fMergingParNotCalculatedV0PosV0Neg(0),
  40   fFracOfMergedRowV0PosV0Neg(0),
  41   fClosestRowAtDCAV0PosV0Neg(0),
  42   fMergingParNotCalculatedV0NegV0Pos(0),
  43   fFracOfMergedRowV0NegV0Pos(0),
  44   fClosestRowAtDCAV0NegV0Pos(0),
  45   fMergingParNotCalculatedV0PosV0Pos(0),
  46   fFracOfMergedRowV0PosV0Pos(0),
  47   fClosestRowAtDCAV0PosV0Pos(0),
  48   fMergingParNotCalculatedV0NegV0Neg(0),
  49   fFracOfMergedRowV0NegV0Neg(0),
  50   fClosestRowAtDCAV0NegV0Neg(0)
  51 {
  52   // Default constructor
  53   fTrack1 = 0;
  54   fTrack2 = 0;
  55   SetDefaultHalfFieldMergingPar();
  56 }
  57
  58 AliFemtoPair::AliFemtoPair(AliFemtoParticle* a, AliFemtoParticle* b)
  59   : fTrack1(a), fTrack2(b),
  60   fPairAngleEP(0),
  61   fNonIdParNotCalculated(0),
  62   fDKSide(0),
  63   fDKOut(0),
  64   fDKLong(0),
  65   fCVK(0),
  66   fKStarCalc(0),
  67   fNonIdParNotCalculatedGlobal(0),
  68   fMergingParNotCalculated(0),
  69   fWeightedAvSep(0),
  70   fFracOfMergedRow(0),
  71   fClosestRowAtDCA(0),
  72   fMergingParNotCalculatedTrkV0Pos(0),
  73   fFracOfMergedRowTrkV0Pos(0),
  74   fClosestRowAtDCATrkV0Pos(0),
  75   fMergingParNotCalculatedTrkV0Neg(0),
  76   fFracOfMergedRowTrkV0Neg(0),
  77   fClosestRowAtDCATrkV0Neg(0),
  78   fMergingParNotCalculatedV0PosV0Neg(0),
  79   fFracOfMergedRowV0PosV0Neg(0),
  80   fClosestRowAtDCAV0PosV0Neg(0),
  81   fMergingParNotCalculatedV0NegV0Pos(0),
  82   fFracOfMergedRowV0NegV0Pos(0),
  83   fClosestRowAtDCAV0NegV0Pos(0),
  84   fMergingParNotCalculatedV0PosV0Pos(0),
  85   fFracOfMergedRowV0PosV0Pos(0),
  86   fClosestRowAtDCAV0PosV0Pos(0),
  87   fMergingParNotCalculatedV0NegV0Neg(0),
  88   fFracOfMergedRowV0NegV0Neg(0),
  89   fClosestRowAtDCAV0NegV0Neg(0)
  90 {
  91   // Construct a pair from two particles
  92   SetDefaultHalfFieldMergingPar();
  93 }
  94
  95 void AliFemtoPair::SetDefaultHalfFieldMergingPar(){
  96   fgMaxDuInner = 3;
  97   fgMaxDzInner = 4.;
  98   fgMaxDuOuter = 4.;
  99   fgMaxDzOuter = 6.;
 100 }
 101 void AliFemtoPair::SetDefaultFullFieldMergingPar(){
 102   // Set default TPC merging parameters for STAR TPC
 103   fgMaxDuInner = 0.8;
 104   fgMaxDzInner = 3.;
 105   fgMaxDuOuter = 1.4;
 106   fgMaxDzOuter = 3.2;
 107 }
 108 void AliFemtoPair::SetMergingPar(double aMaxDuInner, double aMaxDzInner,
 109                               double aMaxDuOuter, double aMaxDzOuter)
 110 {
 111   // Set TPC merging parameters for STAR TPC
 112   fgMaxDuInner = aMaxDuInner;
 113   fgMaxDzInner = aMaxDzInner;
 114   fgMaxDuOuter = aMaxDuOuter;
 115   fgMaxDzOuter = aMaxDzOuter;
 116 }
 117
 118 AliFemtoPair::~AliFemtoPair() {
 119   // Destructor
 120 /* no-op */
 121 }
 122
 123 AliFemtoPair::AliFemtoPair(const AliFemtoPair &aPair):
 124   fTrack1(0), fTrack2(0),
 125   fPairAngleEP(0),
 126   fNonIdParNotCalculated(0),
 127   fDKSide(0),
 128   fDKOut(0),
 129   fDKLong(0),
 130   fCVK(0),
 131   fKStarCalc(0),
 132   fNonIdParNotCalculatedGlobal(0),
 133   fMergingParNotCalculated(0),
 134   fWeightedAvSep(0),
 135   fFracOfMergedRow(0),
 136   fClosestRowAtDCA(0),
 137   fMergingParNotCalculatedTrkV0Pos(0),
 138   fFracOfMergedRowTrkV0Pos(0),
 139   fClosestRowAtDCATrkV0Pos(0),
 140   fMergingParNotCalculatedTrkV0Neg(0),
 141   fFracOfMergedRowTrkV0Neg(0),
 142   fClosestRowAtDCATrkV0Neg(0),
 143   fMergingParNotCalculatedV0PosV0Neg(0),
 144   fFracOfMergedRowV0PosV0Neg(0),
 145   fClosestRowAtDCAV0PosV0Neg(0),
 146   fMergingParNotCalculatedV0NegV0Pos(0),
 147   fFracOfMergedRowV0NegV0Pos(0),
 148   fClosestRowAtDCAV0NegV0Pos(0),
 149   fMergingParNotCalculatedV0PosV0Pos(0),
 150   fFracOfMergedRowV0PosV0Pos(0),
 151   fClosestRowAtDCAV0PosV0Pos(0),
 152   fMergingParNotCalculatedV0NegV0Neg(0),
 153   fFracOfMergedRowV0NegV0Neg(0),
 154   fClosestRowAtDCAV0NegV0Neg(0)
 155 {
 156   // Copy constructor
 157   fTrack1 = aPair.fTrack1;
 158   fTrack2 = aPair.fTrack2;
 159
 160   fNonIdParNotCalculated = aPair.fNonIdParNotCalculated;
 161   fDKSide = aPair.fDKSide;
 162   fDKOut = aPair.fDKOut;
 163   fDKLong = aPair.fDKLong;
 164   fCVK = aPair.fCVK;
 165   fKStarCalc = aPair.fKStarCalc;
 166
 167   fNonIdParNotCalculatedGlobal = aPair.fNonIdParNotCalculatedGlobal;
 168
 169   fMergingParNotCalculated = aPair.fMergingParNotCalculated;
 170   fWeightedAvSep = aPair.fWeightedAvSep;
 171   fFracOfMergedRow = aPair.fFracOfMergedRow;
 172   fClosestRowAtDCA = aPair.fClosestRowAtDCA;
 173
 174   fMergingParNotCalculatedTrkV0Pos = aPair.fMergingParNotCalculatedTrkV0Pos;
 175   fFracOfMergedRowTrkV0Pos = aPair.fFracOfMergedRowTrkV0Pos;
 176   fClosestRowAtDCATrkV0Pos = aPair.fClosestRowAtDCATrkV0Pos;
 177
 178   fMergingParNotCalculatedTrkV0Neg = aPair.fMergingParNotCalculatedTrkV0Neg;
 179   fFracOfMergedRowTrkV0Neg = aPair.fFracOfMergedRowTrkV0Neg;
 180   fClosestRowAtDCATrkV0Neg = aPair.fClosestRowAtDCATrkV0Neg;
 181
 182   fMergingParNotCalculatedV0PosV0Neg = aPair.fMergingParNotCalculatedV0PosV0Neg;
 183   fFracOfMergedRowV0PosV0Neg = aPair.fFracOfMergedRowV0PosV0Neg;
 184   fClosestRowAtDCAV0PosV0Neg = aPair.fClosestRowAtDCAV0PosV0Neg;
 185
 186   fMergingParNotCalculatedV0NegV0Pos = aPair.fMergingParNotCalculatedV0NegV0Pos;
 187   fFracOfMergedRowV0NegV0Pos = aPair.fFracOfMergedRowV0NegV0Pos;
 188   fClosestRowAtDCAV0NegV0Pos = aPair.fClosestRowAtDCAV0NegV0Pos;
 189
 190   fMergingParNotCalculatedV0PosV0Pos = aPair.fMergingParNotCalculatedV0PosV0Pos;
 191   fFracOfMergedRowV0PosV0Pos = aPair.fFracOfMergedRowV0PosV0Pos;
 192   fClosestRowAtDCAV0PosV0Pos = aPair.fClosestRowAtDCAV0PosV0Pos;
 193
 194   fMergingParNotCalculatedV0NegV0Neg = aPair.fMergingParNotCalculatedV0NegV0Neg;
 195   fFracOfMergedRowV0NegV0Neg = aPair.fFracOfMergedRowV0NegV0Neg;
 196   fClosestRowAtDCAV0NegV0Neg = aPair.fClosestRowAtDCAV0NegV0Neg;
 197 }
 198
 199 AliFemtoPair& AliFemtoPair::operator=(const AliFemtoPair &aPair)
 200 {
 201   // Assignment operator
 202   if (this == &aPair)
 203     return *this;
 204
 205   fTrack1 = aPair.fTrack1;
 206   fTrack2 = aPair.fTrack2;
 207
 208   fNonIdParNotCalculated = aPair.fNonIdParNotCalculated;
 209   fDKSide = aPair.fDKSide;
 210   fDKOut = aPair.fDKOut;
 211   fDKLong = aPair.fDKLong;
 212   fCVK = aPair.fCVK;
 213   fKStarCalc = aPair.fKStarCalc;
 214
 215   fNonIdParNotCalculatedGlobal = aPair.fNonIdParNotCalculatedGlobal;
 216
 217   fMergingParNotCalculated = aPair.fMergingParNotCalculated;
 218   fWeightedAvSep = aPair.fWeightedAvSep;
 219   fFracOfMergedRow = aPair.fFracOfMergedRow;
 220   fClosestRowAtDCA = aPair.fClosestRowAtDCA;
 221
 222   fMergingParNotCalculatedTrkV0Pos = aPair.fMergingParNotCalculatedTrkV0Pos;
 223   fFracOfMergedRowTrkV0Pos = aPair.fFracOfMergedRowTrkV0Pos;
 224   fClosestRowAtDCATrkV0Pos = aPair.fClosestRowAtDCATrkV0Pos;
 225
 226   fMergingParNotCalculatedTrkV0Neg = aPair.fMergingParNotCalculatedTrkV0Neg;
 227   fFracOfMergedRowTrkV0Neg = aPair.fFracOfMergedRowTrkV0Neg;
 228   fClosestRowAtDCATrkV0Neg = aPair.fClosestRowAtDCATrkV0Neg;
 229
 230   fMergingParNotCalculatedV0PosV0Neg = aPair.fMergingParNotCalculatedV0PosV0Neg;
 231   fFracOfMergedRowV0PosV0Neg = aPair.fFracOfMergedRowV0PosV0Neg;
 232   fClosestRowAtDCAV0PosV0Neg = aPair.fClosestRowAtDCAV0PosV0Neg;
 233
 234   fMergingParNotCalculatedV0NegV0Pos = aPair.fMergingParNotCalculatedV0NegV0Pos;
 235   fFracOfMergedRowV0NegV0Pos = aPair.fFracOfMergedRowV0NegV0Pos;
 236   fClosestRowAtDCAV0NegV0Pos = aPair.fClosestRowAtDCAV0NegV0Pos;
 237
 238   fMergingParNotCalculatedV0PosV0Pos = aPair.fMergingParNotCalculatedV0PosV0Pos;
 239   fFracOfMergedRowV0PosV0Pos = aPair.fFracOfMergedRowV0PosV0Pos;
 240   fClosestRowAtDCAV0PosV0Pos = aPair.fClosestRowAtDCAV0PosV0Pos;
 241
 242   fMergingParNotCalculatedV0NegV0Neg = aPair.fMergingParNotCalculatedV0NegV0Neg;
 243   fFracOfMergedRowV0NegV0Neg = aPair.fFracOfMergedRowV0NegV0Neg;
 244   fClosestRowAtDCAV0NegV0Neg = aPair.fClosestRowAtDCAV0NegV0Neg;
 245
 246   return *this;
 247 }
 248
 249 //________________________
 250 double AliFemtoPair::GetPairAngleEP() const
 251 {
 252         return fPairAngleEP;
 253 }
 254 //_________________
 255 double AliFemtoPair::MInv() const
 256 {
 257   // invariant mass
 258     double tInvariantMass = abs(fTrack1->FourMomentum() + fTrack2->FourMomentum());
 259     return (tInvariantMass);
 260 }
 261 //_________________
 262 double AliFemtoPair::KT() const
 263 {
 264   // transverse momentum
 265   double  tmp =
 266     (fTrack1->FourMomentum() + fTrack2->FourMomentum()).Perp();
 267   tmp *= .5;
 268
 269   return (tmp);
 270 }
 271 //_________________
 272 double AliFemtoPair::Rap() const
 273 {
 274   // longitudinal pair rapidity : Y = 0.5 ::log( E1 + E2 + pz1 + pz2 / E1 + E2 - pz1 - pz2 )
 275   double  tmp = 0.5 * log (
 276                            (fTrack1->FourMomentum().e() + fTrack2->FourMomentum().e() + fTrack1->FourMomentum().z() + fTrack2->FourMomentum().z()) /
 277                            (fTrack1->FourMomentum().e() + fTrack2->FourMomentum().e() - fTrack1->FourMomentum().z() - fTrack2->FourMomentum().z())
 278                            ) ;
 279   return (tmp);
 280 }
 281 //_________________
 282 double AliFemtoPair::EmissionAngle() const {
 283   // emission angle
 284   double pxTotal = this->FourMomentumSum().x();
 285   double pyTotal = this->FourMomentumSum().y();
 286   double angle = atan2(pyTotal,pxTotal)*180.0/3.1415926536;
 287   if (angle<0.0) angle+=360.0;
 288   return angle;
 289 }
 290 //_________________
 291 // get rid of ambiguously-named method fourMomentum() and replace it with
 292 // fourMomentumSum() and fourMomentumDiff() - mal 13feb2000
 293 AliFemtoLorentzVector AliFemtoPair::FourMomentumSum() const
 294 {
 295   // total momentum
 296   AliFemtoLorentzVector temp = fTrack1->FourMomentum()+fTrack2->FourMomentum();
 297   return temp;
 298 }
 299 AliFemtoLorentzVector AliFemtoPair::FourMomentumDiff() const
 300 {
 301   // momentum difference
 302   AliFemtoLorentzVector temp = fTrack1->FourMomentum()-fTrack2->FourMomentum();
 303   return temp;
 304 }
 305 //__________________________________
 306 void AliFemtoPair::QYKPCMS(double& qP, double& qT, double& q0) const
 307 {
 308   // Yano-Koonin-Podgoretskii Parametrisation in CMS
 309   ////
 310   // calculate momentum difference in source rest frame (= lab frame)
 311   ////
 312   AliFemtoLorentzVector l1 = fTrack1->FourMomentum() ;
 313   AliFemtoLorentzVector l2 = fTrack2->FourMomentum() ;
 314   AliFemtoLorentzVector  l ;
 315   // random ordering of the particles
 316   if ( rand()/(double)RAND_MAX > 0.50 )
 317     { l = l1-l2 ; }
 318   else
 319     { l = l2-l1 ; } ;
 320   // fill momentum differences into return variables
 321   qP = l.z() ;
 322   qT = l.vect().Perp() ;
 323   q0 = l.e() ;
 324 }
 325 //___________________________________
 326 void AliFemtoPair::QYKPLCMS(double& qP, double& qT, double& q0) const
 327 {
 328   // Yano-Koonin-Podgoretskii Parametrisation in LCMS
 329   ////
 330   //  calculate momentum difference in LCMS : frame where pz1 + pz2 = 0
 331   ////
 332   AliFemtoLorentzVector l1 = fTrack1->FourMomentum() ;
 333   AliFemtoLorentzVector l2 = fTrack2->FourMomentum() ;
 334   // determine beta to LCMS
 335   double beta = (l1.z()+l2.z()) / (l1.e()+l2.e()) ;
 336   double beta2 =  beta*beta ;
 337   // unfortunately STAR Class lib knows only boost(particle) not boost(beta) :(
 338   // -> create particle with velocity beta and mass 1.0
 339   // actually this is : dummyPz = ::sqrt( (dummyMass*dummyMass*beta2) / (1-beta2) ) ;
 340   double dummyPz = ::sqrt( (beta2) / (1-beta2) ) ;
 341   // boost in the correct direction
 342   if (beta>0.0) { dummyPz = -dummyPz; } ;
 343   // create dummy particle
 344   AliFemtoLorentzVector  l(0.0, 0.0, dummyPz) ;
 345   double dummyMass = 1.0 ;
 346   l.SetE(l.vect().MassHypothesis(dummyMass) );
 347   // boost particles along the beam into a frame with velocity beta
 348   AliFemtoLorentzVector l1boosted = l1.boost(l) ;
 349   AliFemtoLorentzVector l2boosted = l2.boost(l) ;
 350   // caculate the momentum difference with random ordering of the particle
 351   if ( rand()/(double)RAND_MAX >0.50)
 352     { l = l1boosted-l2boosted ; }
 353   else
 354     { l = l2boosted-l1boosted ;} ;
 355   // fill momentum differences into return variables
 356   qP = l.z() ;
 357   qT = l.vect().Perp() ;
 358   q0 = l.e() ;
 359 }
 360 //___________________________________
 361 // Yano-Koonin-Podgoretskii Parametrisation in pair rest frame
 362 void AliFemtoPair::QYKPPF(double& qP, double& qT, double& q0) const
 363 {
 364   ////
 365   //  calculate momentum difference in pair rest frame : frame where (pz1 + pz2, py1 + py2, px1 + px2) = (0,0,0)
 366   ////
 367   AliFemtoLorentzVector l1 = fTrack1->FourMomentum() ;
 368   AliFemtoLorentzVector l2 = fTrack2->FourMomentum() ;
 369   // the center of gravity of the pair travels with l
 370   AliFemtoLorentzVector  l = l1 + l2 ;
 371   l = -l ;
 372   l.SetE(-l.e()) ;
 373   // boost particles
 374   AliFemtoLorentzVector l1boosted = l1.boost(l) ;
 375   AliFemtoLorentzVector l2boosted = l2.boost(l) ;
 376   // caculate the momentum difference with random ordering of the particle
 377   if ( rand()/(double)RAND_MAX > 0.50)
 378     { l = l1boosted-l2boosted ; }
 379   else
 380     { l = l2boosted-l1boosted ;} ;
 381   // fill momentum differences into return variables
 382   qP = l.z();
 383   qT = l.vect().Perp();
 384   q0 = l.e();
 385 }
 386 //_________________
 387 double AliFemtoPair::QOutCMS() const
 388 {
 389   // relative momentum out component in lab frame
 390     AliFemtoThreeVector tmp1 = fTrack1->FourMomentum().vect();
 391     AliFemtoThreeVector tmp2 = fTrack2->FourMomentum().vect();
 392
 393     double dx = tmp1.x() - tmp2.x();
 394     double xt = tmp1.x() + tmp2.x();
 395
 396     double dy = tmp1.y() - tmp2.y();
 397     double yt = tmp1.y() + tmp2.y();
 398
 399     double k1 = (::sqrt(xt*xt+yt*yt));
 400     double k2 = (dx*xt+dy*yt);
 401     double tmp = k2/k1;
 402     return (tmp);
 403 }
 404 //_________________
 405 double AliFemtoPair::QSideCMS() const
 406 {
 407   // relative momentum side component in lab frame
 408     AliFemtoThreeVector tmp1 = fTrack1->FourMomentum().vect();
 409     AliFemtoThreeVector tmp2 = fTrack2->FourMomentum().vect();
 410
 411     double x1 = tmp1.x();  double y1 = tmp1.y();
 412     double x2 = tmp2.x();  double y2 = tmp2.y();
 413
 414     double xt = x1+x2;  double yt = y1+y2;
 415     double k1 = ::sqrt(xt*xt+yt*yt);
 416
 417     double tmp = 2.0*(x2*y1-x1*y2)/k1;
 418     return (tmp);
 419 }
 420
 421 //_________________________
 422 double AliFemtoPair::QLongCMS() const
 423 {
 424   // relative momentum component in lab frame
 425     AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 426     AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 427
 428     double dz = tmp1.z() - tmp2.z();
 429     double zz = tmp1.z() + tmp2.z();
 430
 431     double dt = tmp1.t() - tmp2.t();
 432     double tt = tmp1.t() + tmp2.t();
 433
 434     double beta = zz/tt;
 435     double gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
 436
 437     double temp = gamma*(dz - beta*dt);
 438     return (temp);
 439 }
 440
 441 //________________________________
 442 double AliFemtoPair::QOutPf() const
 443 {
 444   // relative momentum out component in pair frame
 445   AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 446   AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 447
 448   double dt = tmp1.t() - tmp2.t();
 449   double tt = tmp1.t() + tmp2.t();
 450
 451   double xt = tmp1.x() + tmp2.x();
 452   double yt = tmp1.y() + tmp2.y();
 453
 454   double k1 = ::sqrt(xt*xt + yt*yt);
 455   double bOut = k1/tt;
 456   double gOut = 1.0/TMath::Sqrt((1.-bOut)*(1.+bOut));
 457
 458   double temp = gOut*(QOutCMS() - bOut*dt);
 459   return (temp);
 460 }
 461
 462 //___________________________________
 463 double AliFemtoPair::QSidePf() const
 464 {
 465   // relative momentum side component in pair frame
 466
 467  return(this->QSideCMS());
 468 }
 469
 470 //___________________________________
 471
 472 double AliFemtoPair::QLongPf() const
 473 {
 474   // relative momentum long component in pair frame
 475
 476   return(this->QLongCMS());
 477 }
 478
 479 //___________________________________
 480 double AliFemtoPair::QOutBf(double /* beta */) const
 481 {
 482   // relative momentum out component
 483  return(this->QOutCMS());
 484 }
 485
 486 //___________________________________
 487
 488 double AliFemtoPair::QSideBf(double /* beta */) const
 489 {
 490   // relative momentum side component
 491  return(this->QSideCMS());
 492 }
 493
 494 //___________________________________
 495 double AliFemtoPair::QLongBf(double beta) const
 496 {
 497   // relative momentum long component
 498     AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 499     AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 500
 501     double dz = tmp1.z() - tmp2.z();
 502     double dt = tmp1.t() + tmp2.t();
 503
 504     double gamma = 1.0/::sqrt((1.-beta)*(1.+beta));
 505
 506     double temp = gamma*(dz - beta*dt);
 507     return (temp);
 508 }
 509
 510 double AliFemtoPair::Quality() const {
 511   // Calculate split quality of the pair
 512   unsigned long mapMask0 = 0xFFFFFF00;
 513   unsigned long mapMask1 = 0x1FFFFF;
 514   unsigned long padRow1To24Track1 = fTrack1->TopologyMap(0) & mapMask0;
 515   unsigned long padRow25To45Track1 = fTrack1->TopologyMap(1) & mapMask1;
 516   unsigned long padRow1To24Track2 = fTrack2->TopologyMap(0) & mapMask0;
 517   unsigned long padRow25To45Track2 = fTrack2->TopologyMap(1) & mapMask1;
 518   // AND logic
 519   unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
 520   unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
 521   // XOR logic
 522   unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
 523   unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
 524   unsigned long bitI;
 525   int ibits;
 526   int tQuality = 0;
 527   double normQual = 0.0;
 528   int tMaxQuality = fTrack1->NumberOfHits() + fTrack2->NumberOfHits();
 529   for (ibits=8;ibits<=31;ibits++) {
 530     bitI = 0;
 531     bitI |= 1UL<<(ibits);
 532     if ( onePad1To24 & bitI ) {
 533       tQuality++;
 534       continue;
 535     }
 536     else{
 537       if ( bothPads1To24 & bitI ) tQuality--;
 538     }
 539   }
 540   for (ibits=0;ibits<=20;ibits++) {
 541     bitI = 0;
 542     bitI |= 1UL<<(ibits);
 543     if ( onePad25To45 & bitI ) {
 544       tQuality++;
 545       continue;
 546     }
 547     else{
 548       if ( bothPads25To45 & bitI ) tQuality--;
 549     }
 550   }
 551   normQual = (double)tQuality/( (double) tMaxQuality );
 552   return ( normQual );
 553
 554 }
 555
 556 double AliFemtoPair::Quality2() const {
 557   // second implementation of split quality
 558   unsigned long mapMask0 = 0xFFFFFF00;
 559   unsigned long mapMask1 = 0x1FFFFF;
 560   unsigned long padRow1To24Track1 = fTrack1->TopologyMap(0) & mapMask0;
 561   unsigned long padRow25To45Track1 = fTrack1->TopologyMap(1) & mapMask1;
 562   unsigned long padRow1To24Track2 = fTrack2->TopologyMap(0) & mapMask0;
 563   unsigned long padRow25To45Track2 = fTrack2->TopologyMap(1) & mapMask1;
 564
 565   // AND logic
 566   //unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
 567   //unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
 568
 569   // XOR logic
 570   unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
 571   unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
 572   unsigned long bitI;
 573   int ibits;
 574   int tQuality = 0;
 575   double normQual = 0.0;
 576   int tMaxQuality = fTrack1->NumberOfHits() + fTrack2->NumberOfHits();
 577   for (ibits=8;ibits<=31;ibits++) {
 578     bitI = 0;
 579     bitI |= 1UL<<(ibits);
 580     if ( onePad1To24 & bitI ) {
 581       tQuality++;
 582       continue;
 583     }
 584     //else{
 585     //if ( bothPads1To24 & bitI ) tQuality--;
 586     //}
 587   }
 588   for (ibits=0;ibits<=20;ibits++) {
 589     bitI = 0;
 590     bitI |= 1UL<<(ibits);
 591     if ( onePad25To45 & bitI ) {
 592       tQuality++;
 593       continue;
 594     }
 595     //else{
 596     //if ( bothPads25To45 & bitI ) tQuality--;
 597     //}
 598   }
 599   normQual = (double)tQuality/( (double) tMaxQuality );
 600   return ( normQual );
 601
 602 }
 603
 604
 605 double AliFemtoPair::NominalTpcExitSeparation() const {
 606   // separation at exit from STAR TPC
 607   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->Track()->NominalTpcExitPoint();
 608   return (diff.Mag());
 609 }
 610
 611 double AliFemtoPair::NominalTpcEntranceSeparation() const {
 612   // separation at entrance to STAR TPC
 613   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->Track()->NominalTpcEntrancePoint();
 614   return (diff.Mag());
 615 }
 616
 617 // double AliFemtoPair::NominalTpcAverageSeparation() const {
 618 //   // average separation in STAR TPC
 619 //   AliFemtoThreeVector diff;
 620 //   double tAveSep = 0.0;
 621 //   int ipt = 0;
 622 //   if (fTrack1->fNominalPosSample && fTrack2->fNominalPosSample){
 623 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
 624 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
 625 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
 626 //       fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
 627 //       fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
 628 //       fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
 629 //       ipt<11
 630 //       ){
 631 //     //  for (int ipt=0; ipt<11; ipt++){
 632 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
 633 //     ipt++;
 634 //     tAveSep += diff.Mag();
 635 //   }
 636 //   tAveSep = tAveSep/(ipt+1.);
 637 //   return (tAveSep);}
 638 //   else return -1;
 639 // }
 640
 641 double AliFemtoPair::OpeningAngle() const {
 642   // opening angle
 643  return 57.296* fTrack1->FourMomentum().vect().Angle( fTrack2->FourMomentum().vect() );
 644 //   AliFemtoThreeVector p1 = fTrack1->FourMomentum().vect();
 645 //   AliFemtoThreeVector p2 = fTrack2->FourMomentum().vect();
 646 //   return 57.296*(p1.phi()-p2.phi());
 647 //   //double dAngInv = 57.296*acos((p1.dot(p2))/(p1.Mag()*p2.Mag()));
 648 //   //return (dAngInv);
 649 }
 650 //_________________
 651
 652
 653 double AliFemtoPair::KStarFlipped() const {
 654   // kstar with sign flipped
 655   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 656
 657   AliFmThreeVectorD qwe = tP1.vect();
 658   qwe *= -1.; // flip it
 659   tP1.SetVect(qwe);
 660
 661   AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 662   double tMass = abs(tSum);
 663   AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 664   double tGamma = tSum.e()/tMass;
 665   AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 666                               (tGammaBeta*tGammaBeta))*tGammaBeta;
 667   AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 668                       tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 669 //VP  tP1.vect() *= -1.; // unflip it
 670   return tK.vect().Mag();
 671 }
 672
 673 //double AliFemtoPair::CVK() const{
 674 //const AliFemtoLorentzVector& tP1 = fTrack1->FourMomentum();
 675 //AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 676 //double tMass = abs(tSum);
 677 //AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 678 //double tGamma = tSum.e()/tMass;
 679 //AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 680 //                    (tGammaBeta*tGammaBeta))*tGammaBeta;
 681 //AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 682 //            tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 683 //return (tK.vect())*tGammaBeta/tK.vect().Magnitude()/tGammaBeta.Magnitude();
 684 //}
 685
 686 double AliFemtoPair::CVKFlipped() const{
 687   // CVK with sign flipped
 688   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 689   AliFmThreeVectorD qwe = tP1.vect();
 690   qwe *= -1.; // flip it
 691   tP1.SetVect(qwe);
 692
 693   AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 694   double tMass = abs(tSum);
 695   AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 696   double tGamma = tSum.e()/tMass;
 697   AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 698                               (tGammaBeta*tGammaBeta))*tGammaBeta;
 699   AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 700                       tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 701 //VP  tP1.vect() *= -1.; // unflip it
 702   return (tK.vect())*tGammaBeta/tGamma;
 703 }
 704
 705 double AliFemtoPair::PInv() const{
 706   // invariant total momentum
 707   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 708   AliFemtoLorentzVector tP2 = fTrack2->FourMomentum();
 709   double tP = (tP1.px()+tP2.px())*(tP1.px()+tP2.px())+
 710               (tP1.py()+tP2.py())*(tP1.py()+tP2.py())+
 711               (tP1.pz()+tP2.pz())*(tP1.pz()+tP2.pz())-
 712               (tP1.e() -tP2.e() )*(tP1.e() -tP2.e() );
 713   return ::sqrt(fabs(tP));
 714 }
 715
 716 double AliFemtoPair::QInvFlippedXY() const{
 717   // qinv with X and Y flipped
 718   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 719   tP1.SetX(-1.*tP1.x());
 720   tP1.SetY(-1.*tP1.y());
 721   AliFemtoLorentzVector tDiff = (tP1-fTrack2->FourMomentum());
 722   return ( -1.* tDiff.m());
 723 }
 724
 725 void AliFemtoPair::CalcNonIdPar() const{ // fortran like function! faster?
 726   // Calculate generalized relative mometum
 727   // Use this instead of qXYZ() function when calculating
 728   // anything for non-identical particles
 729   fNonIdParNotCalculated=0;
 730   double px1 = fTrack1->FourMomentum().vect().x();
 731   double py1 = fTrack1->FourMomentum().vect().y();
 732   double pz1 = fTrack1->FourMomentum().vect().z();
 733   double pE1  = fTrack1->FourMomentum().e();
 734   double tParticle1Mass = ::sqrt(pE1*pE1 - px1*px1 - py1*py1 - pz1*pz1);
 735   double px2 = fTrack2->FourMomentum().vect().x();
 736   double py2 = fTrack2->FourMomentum().vect().y();
 737   double pz2 = fTrack2->FourMomentum().vect().z();
 738   double pE2  = fTrack2->FourMomentum().e();
 739   double tParticle2Mass = ::sqrt(pE2*pE2 - px2*px2 - py2*py2 - pz2*pz2);
 740
 741   double tPx = px1+px2;
 742   double tPy = py1+py2;
 743   double tPz = pz1+pz2;
 744   double tPE = pE1+pE2;
 745
 746   double tPtrans = tPx*tPx + tPy*tPy;
 747   double tMtrans = tPE*tPE - tPz*tPz;
 748   double tPinv =   ::sqrt(tMtrans - tPtrans);
 749   tMtrans = ::sqrt(tMtrans);
 750   tPtrans = ::sqrt(tPtrans);
 751
 752   double tQinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
 753     (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
 754
 755   double tQ = (tParticle1Mass*tParticle1Mass - tParticle2Mass*tParticle2Mass)/tPinv;
 756   tQ = sqrt ( tQ*tQ - tQinvL);
 757
 758   fKStarCalc = tQ/2;
 759
 760   // ad 1) go to LCMS
 761   double beta = tPz/tPE;
 762   double gamma = tPE/tMtrans;
 763
 764   double pz1L = gamma * (pz1 - beta * pE1);
 765   double pE1L = gamma * (pE1 - beta * pz1);
 766
 767   // fill histogram for beam projection ( z - axis )
 768   fDKLong = pz1L;
 769
 770   // ad 2) rotation px -> tPt
 771   double px1R = (px1*tPx + py1*tPy)/tPtrans;
 772   double py1R = (-px1*tPy + py1*tPx)/tPtrans;
 773
 774   //fill histograms for side projection ( y - axis )
 775   fDKSide = py1R;
 776
 777   // ad 3) go from LCMS to CMS
 778   beta = tPtrans/tMtrans;
 779   gamma = tMtrans/tPinv;
 780
 781   double px1C = gamma * (px1R - beta * pE1L);
 782
 783   // fill histogram for out projection ( x - axis )
 784   fDKOut  = px1C;
 785
 786   fCVK = (fDKOut*tPtrans + fDKLong*tPz)/fKStarCalc/::sqrt(tPtrans*tPtrans+tPz*tPz);
 787 }
 788
 789
 790 /*void AliFemtoPair::calcNonIdParGlobal() const{ // fortran like function! faster?
 791   fNonIdParNotCalculatedGlobal=0;
 792   double px1 = fTrack1->Track()->PGlobal().x();
 793   double py1 = fTrack1->Track()->PGlobal().y();
 794   double pz1 = fTrack1->Track()->PGlobal().z();
 795   double tParticle1Mass =  fTrack1->FourMomentum().m2();
 796   double pE1  = ::sqrt(tParticle1Mass + px1*px1 + py1*py1 + pz1*pz1);
 797   tParticle1Mass = ::sqrt(tParticle1Mass);
 798
 799   double px2 = fTrack2->Track()->PGlobal().x();
 800   double py2 = fTrack2->Track()->PGlobal().y();
 801   double pz2 = fTrack2->Track()->PGlobal().z();
 802   double tParticle2Mass =  fTrack2->FourMomentum().m2();
 803   double pE2  = ::sqrt(tParticle2Mass + px2*px2 + py2*py2 + pz2*pz2);
 804   tParticle2Mass = ::sqrt(tParticle2Mass);
 805
 806   double Px = px1+px2;
 807   double Py = py1+py2;
 808   double Pz = pz1+pz2;
 809   double PE = pE1+pE2;
 810
 811   double Ptrans = Px*Px + Py*Py;
 812   double Mtrans = PE*PE - Pz*Pz;
 813   double Pinv =   ::sqrt(Mtrans - Ptrans);
 814   Mtrans = ::sqrt(Mtrans);
 815   Ptrans = ::sqrt(Ptrans);
 816
 817   double QinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
 818     (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
 819
 820   double Q = (tParticle1Mass*tParticle1Mass - tParticle2Mass*tParticle2Mass)/Pinv;
 821   Q = sqrt ( Q*Q - QinvL);
 822
 823   kStarCalcGlobal = Q/2;
 824
 825   // ad 1) go to LCMS
 826   double beta = Pz/PE;
 827   double gamma = PE/Mtrans;
 828
 829   double pz1L = gamma * (pz1 - beta * pE1);
 830   double pE1L = gamma * (pE1 - beta * pz1);
 831
 832   // fill histogram for beam projection ( z - axis )
 833   fDKLongGlobal = pz1L;
 834
 835   // ad 2) rotation px -> Pt
 836   double px1R = (px1*Px + py1*Py)/Ptrans;
 837   double py1R = (-px1*Py + py1*Px)/Ptrans;
 838
 839   //fill histograms for side projection ( y - axis )
 840   fDKSideGlobal = py1R;
 841
 842   // ad 3) go from LCMS to CMS
 843   beta = Ptrans/Mtrans;
 844   gamma = Mtrans/Pinv;
 845
 846   double px1C = gamma * (px1R - beta * pE1L);
 847
 848   // fill histogram for out projection ( x - axis )
 849   fDKOutGlobal  = px1C;
 850
 851   fCVKGlobal = (fDKOutGlobal*Ptrans + fDKLongGlobal*Pz)/
 852     kStarCalcGlobal/::sqrt(Ptrans*Ptrans+Pz*Pz);
 853 }*/
 854
 855
 856
 857 // double AliFemtoPair::DcaInsideTpc() const{
 858 //   // dcs inside the STAR TPC
 859 //   double tMinDist=NominalTpcEntranceSeparation();
 860 //   double tExit = NominalTpcExitSeparation();
 861 //   tMinDist = (tExit>tMinDist) ? tMinDist : tExit;
 862 //   double tInsideDist;
 863 //   //tMinDist = 999.;
 864
 865 //   double rMin = 60.;
 866 //   double rMax = 190.;
 867 //   const AliFmPhysicalHelixD& tHelix1 = fTrack1->Helix();
 868 //   const AliFmPhysicalHelixD& tHelix2 = fTrack2->Helix();
 869 //   // --- One is a line and other one a helix
 870 //   //if (tHelix1.mSingularity != tHelix2.mSingularity) return -999.;
 871 //   // --- 2 lines : don't care right now
 872 //   //if (tHelix1.mSingularity)  return -999.;
 873 //   // --- 2 helix
 874 //   double dx = tHelix2.XCenter() - tHelix1.XCenter();
 875 //   double dy = tHelix2.YCenter() - tHelix1.YCenter();
 876 //   double dd = ::sqrt(dx*dx + dy*dy);
 877 //   double r1 = 1/tHelix1.Curvature();
 878 //   double r2 = 1/tHelix2.Curvature();
 879 //   double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
 880
 881 //   double x, y, r;
 882 //   double s;
 883 //   if (fabs(cosAlpha) < 1) {           // two solutions
 884 //     double sinAlpha = sin(acos(cosAlpha));
 885 //     x = tHelix1.XCenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
 886 //     y = tHelix1.YCenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
 887 //     r = ::sqrt(x*x+y*y);
 888 //     if( r > rMin &&  r < rMax &&
 889 //      fabs(atan2(y,x)-fTrack1->Track()->NominalTpcEntrancePoint().phi())< 0.5
 890 //      ){ // first solution inside
 891 //       s = tHelix1.PathLength(x, y);
 892 //       tInsideDist=tHelix2.Distance(tHelix1.At(s));
 893 //       if(tInsideDist<tMinDist) tMinDist = tInsideDist;
 894 //     }
 895 //     else{
 896 //       x = tHelix1.XCenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
 897 //       y = tHelix1.YCenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
 898 //       r = ::sqrt(x*x+y*y);
 899 //       if( r > rMin &&  r < rMax &&
 900 //        fabs(atan2(y,x)-fTrack1->Track()->NominalTpcEntrancePoint().phi())< 0.5
 901 //        ) {  // second solution inside
 902 //         s = tHelix1.PathLength(x, y);
 903 //         tInsideDist=tHelix2.Distance(tHelix1.At(s));
 904 //         if(tInsideDist<tMinDist) tMinDist = tInsideDist;
 905 //       }
 906 //     }
 907 //   }
 908 //   return tMinDist;
 909 // }
 910
 911 // void AliFemtoPair::CalcMergingPar() const{
 912 //   // Calculate merging factor for the pair in STAR TPC
 913 //   fMergingParNotCalculated=0;
 914
 915 //   double tDu, tDz;
 916 //   int tN = 0;
 917 //   fFracOfMergedRow = 0.;
 918 //   fWeightedAvSep =0.;
 919 //   double tDist;
 920 //   double tDistMax = 200.;
 921 //   for(int ti=0 ; ti<45 ; ti++){
 922 //     if(fTrack1->fSect[ti]==fTrack2->fSect[ti] && fTrack1->fSect[ti]!=-1){
 923 //       tDu = fabs(fTrack1->fU[ti]-fTrack2->fU[ti]);
 924 //       tDz = fabs(fTrack1->fZ[ti]-fTrack2->fZ[ti]);
 925 //       tN++;
 926 //       if(ti<13){
 927 //      fFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
 928 //      tDist = ::sqrt(tDu*tDu/fgMaxDuInner/fgMaxDuInner+
 929 //                   tDz*tDz/fgMaxDzInner/fgMaxDzInner);
 930 //      //fFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
 931 //       }
 932 //       else{
 933 //      fFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
 934 //      tDist = ::sqrt(tDu*tDu/fgMaxDuOuter/fgMaxDuOuter+
 935 //                   tDz*tDz/fgMaxDzOuter/fgMaxDzOuter);
 936 //      //fFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
 937 //       }
 938 //       if(tDist<tDistMax){
 939 //      fClosestRowAtDCA = ti+1;
 940 //      tDistMax = tDist;
 941 //       }
 942 //       fWeightedAvSep += tDist;
 943 //     }
 944 //   }
 945 //   if(tN>0){
 946 //     fWeightedAvSep /= tN;
 947 //     fFracOfMergedRow /= tN;
 948 //   }
 949 //   else{
 950 //     fClosestRowAtDCA = -1;
 951 //     fFracOfMergedRow = -1.;
 952 //     fWeightedAvSep = -1.;
 953 //   }
 954 // }
 955 // double AliFemtoPair::TpcExitSeparationTrackV0Pos() const {
 956 // //________________V0 daughters exit/entrance/average separation calc.
 957 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 958
 959 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->TpcV0PosExitPoint();
 960 //   return (diff.Mag());
 961 // }
 962
 963 // double AliFemtoPair::TpcEntranceSeparationTrackV0Pos() const {
 964 // //________________V0 daughters exit/entrance/average separation calc.
 965 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 966 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
 967 //   return (diff.Mag());
 968 // }
 969
 970 // double AliFemtoPair::TpcAverageSeparationTrackV0Pos() const {
 971 // //________________V0 daughters exit/entrance/average separation calc.
 972 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 973 //   AliFemtoThreeVector diff;
 974 //   double tAveSep = 0.0;
 975 //   int ipt = 0;
 976 //   if (fTrack1->fNominalPosSample && fTrack2->fNominalPosSample){
 977 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
 978 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
 979 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
 980 //       fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
 981 //       fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
 982 //       fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
 983 //       (ipt<11)
 984 //       ){
 985 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
 986 //     ipt++;
 987 //     tAveSep += diff.Mag();
 988 //   }
 989 //   tAveSep = tAveSep/(ipt+1.);
 990 //   return (tAveSep);}
 991 //   else return -1;
 992 // }
 993 // double AliFemtoPair::TpcExitSeparationTrackV0Neg() const {
 994 // //_______1st part is a track 2nd is a V0 considering Neg daughter
 995 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->TpcV0NegExitPoint();
 996 //   return (diff.Mag());
 997 // }
 998
 999 // double AliFemtoPair::TpcEntranceSeparationTrackV0Neg() const {
1000 // //_______1st part is a track 2nd is a V0 considering Neg daughter
1001 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1002 //   return (diff.Mag());
1003 // }
1004
1005 // double AliFemtoPair::TpcAverageSeparationTrackV0Neg() const {
1006 // //_______1st part is a track 2nd is a V0 considering Neg daughter
1007 //   AliFemtoThreeVector diff;
1008 //   double tAveSep = 0.0;
1009 //   int ipt = 0;
1010 //   if (fTrack1->fNominalPosSample && fTrack2->fTpcV0NegPosSample){
1011 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
1012 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
1013 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
1014 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1015 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1016 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1017 //       (ipt<11)
1018 //       ){
1019 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1020 //     ipt++;
1021 //     tAveSep += diff.Mag();
1022 //   }
1023 //   tAveSep = tAveSep/(ipt+1.);
1024 //   return (tAveSep);}
1025 //   else return -1;
1026 // }
1027
1028 // double AliFemtoPair::TpcExitSeparationV0PosV0Pos() const {
1029 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1030 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosExitPoint() - fTrack2->TpcV0PosExitPoint();
1031 //   return (diff.Mag());
1032 // }
1033
1034 // double AliFemtoPair::TpcEntranceSeparationV0PosV0Pos() const {
1035 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1036 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
1037 //   return (diff.Mag());
1038 // }
1039 // double AliFemtoPair::TpcAverageSeparationV0PosV0Pos() const {
1040 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1041 //   AliFemtoThreeVector diff;
1042 //   double tAveSep = 0.0;
1043 //   int ipt=0;
1044 //   if (fTrack1->fNominalPosSample && (fTrack2->fNominalPosSample)){
1045 //     while ((fabs(fTrack1->fNominalPosSample[ipt].x())<9999.) &&
1046 //      (fabs(fTrack1->fNominalPosSample[ipt].y())<9999.) &&
1047 //      (fabs(fTrack1->fNominalPosSample[ipt].z())<9999.) &&
1048 //      (fabs(fTrack2->fNominalPosSample[ipt].x())<9999.) &&
1049 //      (fabs(fTrack2->fNominalPosSample[ipt].y())<9999.) &&
1050 //      (fabs(fTrack2->fNominalPosSample[ipt].z())<9999.) &&
1051 //       (ipt<11)
1052 //      ){
1053 //       diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
1054 //       ipt++;
1055 //       tAveSep += diff.Mag();
1056 //     }
1057 //     tAveSep = tAveSep/(ipt+1);
1058 //     return (tAveSep);}
1059 //   else return -1;
1060 // }
1061
1062 // double AliFemtoPair::TpcExitSeparationV0PosV0Neg() const {
1063 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1064 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosExitPoint() - fTrack2->TpcV0NegExitPoint();
1065 //   return (diff.Mag());
1066 // }
1067
1068 // double AliFemtoPair::TpcEntranceSeparationV0PosV0Neg() const {
1069 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1070 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1071 //   return (diff.Mag());
1072 // }
1073 // double AliFemtoPair::TpcAverageSeparationV0PosV0Neg() const {
1074 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1075 //   AliFemtoThreeVector diff;
1076 //   double tAveSep = 0.0;
1077 //   int ipt = 0;
1078 //   if (fTrack1->fNominalPosSample && fTrack2->fTpcV0NegPosSample){
1079 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
1080 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
1081 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
1082 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1083 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1084 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1085 //       (ipt<11)
1086 //       ){
1087 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1088 //     ipt++;
1089 //     tAveSep += diff.Mag();
1090 //   }
1091 //   tAveSep = tAveSep/(ipt+1.);
1092 //   return (tAveSep);}
1093 //   else return -1;
1094 // }
1095 // double AliFemtoPair::TpcExitSeparationV0NegV0Pos() const {
1096 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1097 // // this is to check the upper case
1098 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegExitPoint() - fTrack2->TpcV0PosExitPoint();
1099 //   return (diff.Mag());
1100 // }
1101
1102 // double AliFemtoPair::TpcEntranceSeparationV0NegV0Pos() const {
1103 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1104 // // this is to check the upper case
1105 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
1106 //   return (diff.Mag());
1107 // }
1108 // double AliFemtoPair::TpcAverageSeparationV0NegV0Pos() const {
1109 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1110 // // this is to check the upper case
1111 //    AliFemtoThreeVector diff;
1112 //    double tAveSep = 0.0;
1113 //    int ipt = 0;
1114 //    if ( fTrack1->fTpcV0NegPosSample &&  fTrack2->fNominalPosSample){
1115 //      while (fabs(fTrack1->fTpcV0NegPosSample[ipt].x())<9999. &&
1116 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].y())<9999. &&
1117 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].z())<9999. &&
1118 //          fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
1119 //          fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
1120 //          fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
1121 //          (ipt<11)
1122 //          ){
1123 //        diff = fTrack1->fTpcV0NegPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
1124 //        ipt++;
1125 //        tAveSep += diff.Mag();
1126 //      }
1127 //      tAveSep = tAveSep/(ipt+1);
1128 //      return (tAveSep);}
1129 //      else return -1;
1130 // }
1131 // double AliFemtoPair::TpcExitSeparationV0NegV0Neg() const {
1132 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1133 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegExitPoint() - fTrack2->TpcV0NegExitPoint();
1134 //   return (diff.Mag());
1135 // }
1136
1137 // double AliFemtoPair::TpcEntranceSeparationV0NegV0Neg() const {
1138 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1139 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1140 //   return (diff.Mag());
1141 // }
1142 // double AliFemtoPair::TpcAverageSeparationV0NegV0Neg() const {
1143 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1144 //    AliFemtoThreeVector diff;
1145 //    double tAveSep = 0.0;
1146 //    int ipt=0;
1147 //    if (fTrack1->fTpcV0NegPosSample && fTrack2->fTpcV0NegPosSample){
1148 //      while (fabs(fTrack1->fTpcV0NegPosSample[ipt].x())<9999. &&
1149 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].y())<9999. &&
1150 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].z())<9999. &&
1151 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1152 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1153 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1154 //          (ipt<11)
1155 //          ){
1156 //        diff = fTrack1->fTpcV0NegPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1157 //        ipt++;
1158 //        tAveSep += diff.Mag();
1159 //      }
1160 //      tAveSep = tAveSep/(ipt+1);
1161 //      return (tAveSep);}
1162 //    else return -1;
1163 // }
1164
1165 // void AliFemtoPair::CalcMergingParFctn(short* tmpMergingParNotCalculatedFctn,
1166 //                                 float* tmpZ1,float* tmpU1,
1167 //                                 float* tmpZ2,float* tmpU2,
1168 //                                 int *tmpSect1,int *tmpSect2,
1169 //                                 double* tmpFracOfMergedRow,
1170 //                                 double* tmpClosestRowAtDCA
1171 //                                 ) const{
1172 // // calculate heper variables for merging
1173 //   tmpMergingParNotCalculatedFctn=0;
1174 //   double tDu, tDz;
1175 //   int tN = 0;
1176 //   *tmpFracOfMergedRow = 0.;
1177 //   *tmpClosestRowAtDCA = 0.;
1178 //   double tDist;
1179 //   double tDistMax = 100000000.;
1180 //   for(int ti=0 ; ti<45 ; ti++){
1181 //     if(tmpSect1[ti]==tmpSect2[ti] && tmpSect1[ti]!=-1){
1182 //      tDu = fabs(tmpU1[ti]-tmpU2[ti]);
1183 //      tDz = fabs(tmpZ1[ti]-tmpZ2[ti]);
1184 //      tN++;
1185 //       if(ti<13){
1186 //      *tmpFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
1187 //      tDist = ::sqrt(tDu*tDu/fgMaxDuInner/fgMaxDuInner+
1188 //                   tDz*tDz/fgMaxDzInner/fgMaxDzInner);
1189 //       }
1190 //       else{
1191 //      *tmpFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
1192 //      tDist = ::sqrt(tDu*tDu/fgMaxDuOuter/fgMaxDuOuter+
1193 //                   tDz*tDz/fgMaxDzOuter/fgMaxDzOuter);
1194 //      }
1195 //       if(tDist<tDistMax){
1196 //      fClosestRowAtDCA = ti+1;
1197 //      tDistMax = tDist;
1198 //       }
1199 //       //fWeightedAvSep += tDist; // now, wrong but not used
1200 //     }
1201 //   }
1202 //   if(tN>0){
1203 //     //fWeightedAvSep /= tN;
1204 //     *tmpFracOfMergedRow /= tN;
1205 //   }
1206 //   else{
1207 //     *tmpClosestRowAtDCA = -1;
1208 //     *tmpFracOfMergedRow = -1.;
1209 //     //fWeightedAvSep = -1.;
1210 //   }
1211 // }
1212