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