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;
 402
 403     if(k1!=0) tmp= k2/k1;
 404     else tmp=0;
 405
 406     return (tmp);
 407 }
 408 //_________________
 409 double AliFemtoPair::QSideCMS() const
 410 {
 411   // relative momentum side component in lab frame
 412     AliFemtoThreeVector tmp1 = fTrack1->FourMomentum().vect();
 413     AliFemtoThreeVector tmp2 = fTrack2->FourMomentum().vect();
 414
 415     double x1 = tmp1.x();  double y1 = tmp1.y();
 416     double x2 = tmp2.x();  double y2 = tmp2.y();
 417
 418     double xt = x1+x2;  double yt = y1+y2;
 419     double k1 = ::sqrt(xt*xt+yt*yt);
 420
 421     double tmp;
 422     if(k1!=0) tmp= 2.0*(x2*y1-x1*y2)/k1;
 423     else tmp=0;
 424
 425     return (tmp);
 426 }
 427
 428 //_________________________
 429 double AliFemtoPair::QLongCMS() const
 430 {
 431   // relative momentum component in lab frame
 432     AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 433     AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 434
 435     double dz = tmp1.z() - tmp2.z();
 436     double zz = tmp1.z() + tmp2.z();
 437
 438     double dt = tmp1.t() - tmp2.t();
 439     double tt = tmp1.t() + tmp2.t();
 440
 441     double beta = zz/tt;
 442     double gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
 443
 444     double temp = gamma*(dz - beta*dt);
 445     return (temp);
 446 }
 447
 448 //________________________________
 449 double AliFemtoPair::QOutPf() const
 450 {
 451   // relative momentum out component in pair frame
 452   AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 453   AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 454
 455   double dt = tmp1.t() - tmp2.t();
 456   double tt = tmp1.t() + tmp2.t();
 457
 458   double xt = tmp1.x() + tmp2.x();
 459   double yt = tmp1.y() + tmp2.y();
 460
 461   double k1 = ::sqrt(xt*xt + yt*yt);
 462   double bOut = k1/tt;
 463   double gOut = 1.0/TMath::Sqrt((1.-bOut)*(1.+bOut));
 464
 465   double temp = gOut*(QOutCMS() - bOut*dt);
 466   return (temp);
 467 }
 468
 469 //___________________________________
 470 double AliFemtoPair::QSidePf() const
 471 {
 472   // relative momentum side component in pair frame
 473
 474  return(this->QSideCMS());
 475 }
 476
 477 //___________________________________
 478
 479 double AliFemtoPair::QLongPf() const
 480 {
 481   // relative momentum long component in pair frame
 482
 483   return(this->QLongCMS());
 484 }
 485
 486 //___________________________________
 487 double AliFemtoPair::QOutBf(double /* beta */) const
 488 {
 489   // relative momentum out component
 490  return(this->QOutCMS());
 491 }
 492
 493 //___________________________________
 494
 495 double AliFemtoPair::QSideBf(double /* beta */) const
 496 {
 497   // relative momentum side component
 498  return(this->QSideCMS());
 499 }
 500
 501 //___________________________________
 502 double AliFemtoPair::QLongBf(double beta) const
 503 {
 504   // relative momentum long component
 505     AliFemtoLorentzVector tmp1 = fTrack1->FourMomentum();
 506     AliFemtoLorentzVector tmp2 = fTrack2->FourMomentum();
 507
 508     double dz = tmp1.z() - tmp2.z();
 509     double dt = tmp1.t() + tmp2.t();
 510
 511     double gamma = 1.0/::sqrt((1.-beta)*(1.+beta));
 512
 513     double temp = gamma*(dz - beta*dt);
 514     return (temp);
 515 }
 516
 517 double AliFemtoPair::Quality() const {
 518   // Calculate split quality of the pair
 519   unsigned long mapMask0 = 0xFFFFFF00;
 520   unsigned long mapMask1 = 0x1FFFFF;
 521   unsigned long padRow1To24Track1 = fTrack1->TopologyMap(0) & mapMask0;
 522   unsigned long padRow25To45Track1 = fTrack1->TopologyMap(1) & mapMask1;
 523   unsigned long padRow1To24Track2 = fTrack2->TopologyMap(0) & mapMask0;
 524   unsigned long padRow25To45Track2 = fTrack2->TopologyMap(1) & mapMask1;
 525   // AND logic
 526   unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
 527   unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
 528   // XOR logic
 529   unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
 530   unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
 531   unsigned long bitI;
 532   int ibits;
 533   int tQuality = 0;
 534   double normQual = 0.0;
 535   int tMaxQuality = fTrack1->NumberOfHits() + fTrack2->NumberOfHits();
 536   for (ibits=8;ibits<=31;ibits++) {
 537     bitI = 0;
 538     bitI |= 1UL<<(ibits);
 539     if ( onePad1To24 & bitI ) {
 540       tQuality++;
 541       continue;
 542     }
 543     else{
 544       if ( bothPads1To24 & bitI ) tQuality--;
 545     }
 546   }
 547   for (ibits=0;ibits<=20;ibits++) {
 548     bitI = 0;
 549     bitI |= 1UL<<(ibits);
 550     if ( onePad25To45 & bitI ) {
 551       tQuality++;
 552       continue;
 553     }
 554     else{
 555       if ( bothPads25To45 & bitI ) tQuality--;
 556     }
 557   }
 558   normQual = (double)tQuality/( (double) tMaxQuality );
 559   return ( normQual );
 560
 561 }
 562
 563 double AliFemtoPair::Quality2() const {
 564   // second implementation of split quality
 565   unsigned long mapMask0 = 0xFFFFFF00;
 566   unsigned long mapMask1 = 0x1FFFFF;
 567   unsigned long padRow1To24Track1 = fTrack1->TopologyMap(0) & mapMask0;
 568   unsigned long padRow25To45Track1 = fTrack1->TopologyMap(1) & mapMask1;
 569   unsigned long padRow1To24Track2 = fTrack2->TopologyMap(0) & mapMask0;
 570   unsigned long padRow25To45Track2 = fTrack2->TopologyMap(1) & mapMask1;
 571
 572   // AND logic
 573   //unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
 574   //unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
 575
 576   // XOR logic
 577   unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
 578   unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
 579   unsigned long bitI;
 580   int ibits;
 581   int tQuality = 0;
 582   double normQual = 0.0;
 583   int tMaxQuality = fTrack1->NumberOfHits() + fTrack2->NumberOfHits();
 584   for (ibits=8;ibits<=31;ibits++) {
 585     bitI = 0;
 586     bitI |= 1UL<<(ibits);
 587     if ( onePad1To24 & bitI ) {
 588       tQuality++;
 589       continue;
 590     }
 591     //else{
 592     //if ( bothPads1To24 & bitI ) tQuality--;
 593     //}
 594   }
 595   for (ibits=0;ibits<=20;ibits++) {
 596     bitI = 0;
 597     bitI |= 1UL<<(ibits);
 598     if ( onePad25To45 & bitI ) {
 599       tQuality++;
 600       continue;
 601     }
 602     //else{
 603     //if ( bothPads25To45 & bitI ) tQuality--;
 604     //}
 605   }
 606   normQual = (double)tQuality/( (double) tMaxQuality );
 607   return ( normQual );
 608
 609 }
 610
 611
 612 double AliFemtoPair::NominalTpcExitSeparation() const {
 613   // separation at exit from STAR TPC
 614   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->Track()->NominalTpcExitPoint();
 615   return (diff.Mag());
 616 }
 617
 618 double AliFemtoPair::NominalTpcEntranceSeparation() const {
 619   // separation at entrance to STAR TPC
 620   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->Track()->NominalTpcEntrancePoint();
 621   return (diff.Mag());
 622 }
 623
 624 // double AliFemtoPair::NominalTpcAverageSeparation() const {
 625 //   // average separation in STAR TPC
 626 //   AliFemtoThreeVector diff;
 627 //   double tAveSep = 0.0;
 628 //   int ipt = 0;
 629 //   if (fTrack1->fNominalPosSample && fTrack2->fNominalPosSample){
 630 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
 631 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
 632 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
 633 //       fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
 634 //       fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
 635 //       fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
 636 //       ipt<11
 637 //       ){
 638 //     //  for (int ipt=0; ipt<11; ipt++){
 639 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
 640 //     ipt++;
 641 //     tAveSep += diff.Mag();
 642 //   }
 643 //   tAveSep = tAveSep/(ipt+1.);
 644 //   return (tAveSep);}
 645 //   else return -1;
 646 // }
 647
 648 double AliFemtoPair::OpeningAngle() const {
 649   // opening angle
 650  return 57.296* fTrack1->FourMomentum().vect().Angle( fTrack2->FourMomentum().vect() );
 651 //   AliFemtoThreeVector p1 = fTrack1->FourMomentum().vect();
 652 //   AliFemtoThreeVector p2 = fTrack2->FourMomentum().vect();
 653 //   return 57.296*(p1.phi()-p2.phi());
 654 //   //double dAngInv = 57.296*acos((p1.dot(p2))/(p1.Mag()*p2.Mag()));
 655 //   //return (dAngInv);
 656 }
 657 //_________________
 658
 659
 660 double AliFemtoPair::KStarFlipped() const {
 661   // kstar with sign flipped
 662   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 663
 664   AliFmThreeVectorD qwe = tP1.vect();
 665   qwe *= -1.; // flip it
 666   tP1.SetVect(qwe);
 667
 668   AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 669   double tMass = abs(tSum);
 670   AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 671   double tGamma = tSum.e()/tMass;
 672   AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 673                               (tGammaBeta*tGammaBeta))*tGammaBeta;
 674   AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 675                       tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 676 //VP  tP1.vect() *= -1.; // unflip it
 677   return tK.vect().Mag();
 678 }
 679
 680 //double AliFemtoPair::CVK() const{
 681 //const AliFemtoLorentzVector& tP1 = fTrack1->FourMomentum();
 682 //AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 683 //double tMass = abs(tSum);
 684 //AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 685 //double tGamma = tSum.e()/tMass;
 686 //AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 687 //                    (tGammaBeta*tGammaBeta))*tGammaBeta;
 688 //AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 689 //            tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 690 //return (tK.vect())*tGammaBeta/tK.vect().Magnitude()/tGammaBeta.Magnitude();
 691 //}
 692
 693 double AliFemtoPair::CVKFlipped() const{
 694   // CVK with sign flipped
 695   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 696   AliFmThreeVectorD qwe = tP1.vect();
 697   qwe *= -1.; // flip it
 698   tP1.SetVect(qwe);
 699
 700   AliFemtoLorentzVector tSum = (tP1+fTrack2->FourMomentum());
 701   double tMass = abs(tSum);
 702   AliFmThreeVectorD tGammaBeta = (1./tMass)*tSum.vect();
 703   double tGamma = tSum.e()/tMass;
 704   AliFmThreeVectorD tLongMom  = ((tP1.vect()*tGammaBeta)/
 705                               (tGammaBeta*tGammaBeta))*tGammaBeta;
 706   AliFmLorentzVectorD tK(tGamma*tP1.e() - tP1.vect()*tGammaBeta,
 707                       tP1.vect() + (tGamma-1.)*tLongMom - tP1.e()*tGammaBeta);
 708 //VP  tP1.vect() *= -1.; // unflip it
 709   return (tK.vect())*tGammaBeta/tGamma;
 710 }
 711
 712 double AliFemtoPair::PInv() const{
 713   // invariant total momentum
 714   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 715   AliFemtoLorentzVector tP2 = fTrack2->FourMomentum();
 716   double tP = (tP1.px()+tP2.px())*(tP1.px()+tP2.px())+
 717               (tP1.py()+tP2.py())*(tP1.py()+tP2.py())+
 718               (tP1.pz()+tP2.pz())*(tP1.pz()+tP2.pz())-
 719               (tP1.e() -tP2.e() )*(tP1.e() -tP2.e() );
 720   return ::sqrt(fabs(tP));
 721 }
 722
 723 double AliFemtoPair::QInvFlippedXY() const{
 724   // qinv with X and Y flipped
 725   AliFemtoLorentzVector tP1 = fTrack1->FourMomentum();
 726   tP1.SetX(-1.*tP1.x());
 727   tP1.SetY(-1.*tP1.y());
 728   AliFemtoLorentzVector tDiff = (tP1-fTrack2->FourMomentum());
 729   return ( -1.* tDiff.m());
 730 }
 731
 732 void AliFemtoPair::CalcNonIdPar() const{ // fortran like function! faster?
 733   // Calculate generalized relative mometum
 734   // Use this instead of qXYZ() function when calculating
 735   // anything for non-identical particles
 736   fNonIdParNotCalculated=0;
 737   double px1 = fTrack1->FourMomentum().vect().x();
 738   double py1 = fTrack1->FourMomentum().vect().y();
 739   double pz1 = fTrack1->FourMomentum().vect().z();
 740   double pE1  = fTrack1->FourMomentum().e();
 741   double tParticle1Mass;
 742   if((pE1*pE1 - px1*px1 - py1*py1 - pz1*pz1)>0)
 743     tParticle1Mass = ::sqrt(pE1*pE1 - px1*px1 - py1*py1 - pz1*pz1);
 744   else
 745     tParticle1Mass = 0;
 746
 747   double px2 = fTrack2->FourMomentum().vect().x();
 748   double py2 = fTrack2->FourMomentum().vect().y();
 749   double pz2 = fTrack2->FourMomentum().vect().z();
 750   double pE2  = fTrack2->FourMomentum().e();
 751   double tParticle2Mass;
 752   if((pE2*pE2 - px2*px2 - py2*py2 - pz2*pz2)>0)
 753     tParticle2Mass = ::sqrt(pE2*pE2 - px2*px2 - py2*py2 - pz2*pz2);
 754   else
 755     tParticle2Mass = 0;
 756
 757   double tPx = px1+px2;
 758   double tPy = py1+py2;
 759   double tPz = pz1+pz2;
 760   double tPE = pE1+pE2;
 761
 762   double tPtrans = tPx*tPx + tPy*tPy;
 763   double tMtrans = tPE*tPE - tPz*tPz;
 764   double tPinv =   ::sqrt(tMtrans - tPtrans);
 765   tMtrans = ::sqrt(tMtrans);
 766   tPtrans = ::sqrt(tPtrans);
 767
 768   double tQinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
 769     (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
 770
 771   double tQ = (tParticle1Mass*tParticle1Mass - tParticle2Mass*tParticle2Mass)/tPinv;
 772   tQ = sqrt ( tQ*tQ - tQinvL);
 773
 774   fKStarCalc = tQ/2;
 775
 776   // ad 1) go to LCMS
 777   double beta = tPz/tPE;
 778   double gamma = tPE/tMtrans;
 779
 780   double pz1L = gamma * (pz1 - beta * pE1);
 781   double pE1L = gamma * (pE1 - beta * pz1);
 782
 783   // fill histogram for beam projection ( z - axis )
 784   fDKLong = pz1L;
 785
 786   // ad 2) rotation px -> tPt
 787   double px1R = (px1*tPx + py1*tPy)/tPtrans;
 788   double py1R = (-px1*tPy + py1*tPx)/tPtrans;
 789
 790   //fill histograms for side projection ( y - axis )
 791   fDKSide = py1R;
 792
 793   // ad 3) go from LCMS to CMS
 794   beta = tPtrans/tMtrans;
 795   gamma = tMtrans/tPinv;
 796
 797   double px1C = gamma * (px1R - beta * pE1L);
 798
 799   // fill histogram for out projection ( x - axis )
 800   fDKOut  = px1C;
 801
 802   fCVK = (fDKOut*tPtrans + fDKLong*tPz)/fKStarCalc/::sqrt(tPtrans*tPtrans+tPz*tPz);
 803 }
 804
 805
 806 /*void AliFemtoPair::calcNonIdParGlobal() const{ // fortran like function! faster?
 807   fNonIdParNotCalculatedGlobal=0;
 808   double px1 = fTrack1->Track()->PGlobal().x();
 809   double py1 = fTrack1->Track()->PGlobal().y();
 810   double pz1 = fTrack1->Track()->PGlobal().z();
 811   double tParticle1Mass =  fTrack1->FourMomentum().m2();
 812   double pE1  = ::sqrt(tParticle1Mass + px1*px1 + py1*py1 + pz1*pz1);
 813   tParticle1Mass = ::sqrt(tParticle1Mass);
 814
 815   double px2 = fTrack2->Track()->PGlobal().x();
 816   double py2 = fTrack2->Track()->PGlobal().y();
 817   double pz2 = fTrack2->Track()->PGlobal().z();
 818   double tParticle2Mass =  fTrack2->FourMomentum().m2();
 819   double pE2  = ::sqrt(tParticle2Mass + px2*px2 + py2*py2 + pz2*pz2);
 820   tParticle2Mass = ::sqrt(tParticle2Mass);
 821
 822   double Px = px1+px2;
 823   double Py = py1+py2;
 824   double Pz = pz1+pz2;
 825   double PE = pE1+pE2;
 826
 827   double Ptrans = Px*Px + Py*Py;
 828   double Mtrans = PE*PE - Pz*Pz;
 829   double Pinv =   ::sqrt(Mtrans - Ptrans);
 830   Mtrans = ::sqrt(Mtrans);
 831   Ptrans = ::sqrt(Ptrans);
 832
 833   double QinvL = (pE1-pE2)*(pE1-pE2) - (px1-px2)*(px1-px2) -
 834     (py1-py2)*(py1-py2) - (pz1-pz2)*(pz1-pz2);
 835
 836   double Q = (tParticle1Mass*tParticle1Mass - tParticle2Mass*tParticle2Mass)/Pinv;
 837   Q = sqrt ( Q*Q - QinvL);
 838
 839   kStarCalcGlobal = Q/2;
 840
 841   // ad 1) go to LCMS
 842   double beta = Pz/PE;
 843   double gamma = PE/Mtrans;
 844
 845   double pz1L = gamma * (pz1 - beta * pE1);
 846   double pE1L = gamma * (pE1 - beta * pz1);
 847
 848   // fill histogram for beam projection ( z - axis )
 849   fDKLongGlobal = pz1L;
 850
 851   // ad 2) rotation px -> Pt
 852   double px1R = (px1*Px + py1*Py)/Ptrans;
 853   double py1R = (-px1*Py + py1*Px)/Ptrans;
 854
 855   //fill histograms for side projection ( y - axis )
 856   fDKSideGlobal = py1R;
 857
 858   // ad 3) go from LCMS to CMS
 859   beta = Ptrans/Mtrans;
 860   gamma = Mtrans/Pinv;
 861
 862   double px1C = gamma * (px1R - beta * pE1L);
 863
 864   // fill histogram for out projection ( x - axis )
 865   fDKOutGlobal  = px1C;
 866
 867   fCVKGlobal = (fDKOutGlobal*Ptrans + fDKLongGlobal*Pz)/
 868     kStarCalcGlobal/::sqrt(Ptrans*Ptrans+Pz*Pz);
 869 }*/
 870
 871
 872
 873 // double AliFemtoPair::DcaInsideTpc() const{
 874 //   // dcs inside the STAR TPC
 875 //   double tMinDist=NominalTpcEntranceSeparation();
 876 //   double tExit = NominalTpcExitSeparation();
 877 //   tMinDist = (tExit>tMinDist) ? tMinDist : tExit;
 878 //   double tInsideDist;
 879 //   //tMinDist = 999.;
 880
 881 //   double rMin = 60.;
 882 //   double rMax = 190.;
 883 //   const AliFmPhysicalHelixD& tHelix1 = fTrack1->Helix();
 884 //   const AliFmPhysicalHelixD& tHelix2 = fTrack2->Helix();
 885 //   // --- One is a line and other one a helix
 886 //   //if (tHelix1.mSingularity != tHelix2.mSingularity) return -999.;
 887 //   // --- 2 lines : don't care right now
 888 //   //if (tHelix1.mSingularity)  return -999.;
 889 //   // --- 2 helix
 890 //   double dx = tHelix2.XCenter() - tHelix1.XCenter();
 891 //   double dy = tHelix2.YCenter() - tHelix1.YCenter();
 892 //   double dd = ::sqrt(dx*dx + dy*dy);
 893 //   double r1 = 1/tHelix1.Curvature();
 894 //   double r2 = 1/tHelix2.Curvature();
 895 //   double cosAlpha = (r1*r1 + dd*dd - r2*r2)/(2*r1*dd);
 896
 897 //   double x, y, r;
 898 //   double s;
 899 //   if (fabs(cosAlpha) < 1) {           // two solutions
 900 //     double sinAlpha = sin(acos(cosAlpha));
 901 //     x = tHelix1.XCenter() + r1*(cosAlpha*dx - sinAlpha*dy)/dd;
 902 //     y = tHelix1.YCenter() + r1*(sinAlpha*dx + cosAlpha*dy)/dd;
 903 //     r = ::sqrt(x*x+y*y);
 904 //     if( r > rMin &&  r < rMax &&
 905 //      fabs(atan2(y,x)-fTrack1->Track()->NominalTpcEntrancePoint().phi())< 0.5
 906 //      ){ // first solution inside
 907 //       s = tHelix1.PathLength(x, y);
 908 //       tInsideDist=tHelix2.Distance(tHelix1.At(s));
 909 //       if(tInsideDist<tMinDist) tMinDist = tInsideDist;
 910 //     }
 911 //     else{
 912 //       x = tHelix1.XCenter() + r1*(cosAlpha*dx + sinAlpha*dy)/dd;
 913 //       y = tHelix1.YCenter() + r1*(cosAlpha*dy - sinAlpha*dx)/dd;
 914 //       r = ::sqrt(x*x+y*y);
 915 //       if( r > rMin &&  r < rMax &&
 916 //        fabs(atan2(y,x)-fTrack1->Track()->NominalTpcEntrancePoint().phi())< 0.5
 917 //        ) {  // second solution inside
 918 //         s = tHelix1.PathLength(x, y);
 919 //         tInsideDist=tHelix2.Distance(tHelix1.At(s));
 920 //         if(tInsideDist<tMinDist) tMinDist = tInsideDist;
 921 //       }
 922 //     }
 923 //   }
 924 //   return tMinDist;
 925 // }
 926
 927 // void AliFemtoPair::CalcMergingPar() const{
 928 //   // Calculate merging factor for the pair in STAR TPC
 929 //   fMergingParNotCalculated=0;
 930
 931 //   double tDu, tDz;
 932 //   int tN = 0;
 933 //   fFracOfMergedRow = 0.;
 934 //   fWeightedAvSep =0.;
 935 //   double tDist;
 936 //   double tDistMax = 200.;
 937 //   for(int ti=0 ; ti<45 ; ti++){
 938 //     if(fTrack1->fSect[ti]==fTrack2->fSect[ti] && fTrack1->fSect[ti]!=-1){
 939 //       tDu = fabs(fTrack1->fU[ti]-fTrack2->fU[ti]);
 940 //       tDz = fabs(fTrack1->fZ[ti]-fTrack2->fZ[ti]);
 941 //       tN++;
 942 //       if(ti<13){
 943 //      fFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
 944 //      tDist = ::sqrt(tDu*tDu/fgMaxDuInner/fgMaxDuInner+
 945 //                   tDz*tDz/fgMaxDzInner/fgMaxDzInner);
 946 //      //fFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
 947 //       }
 948 //       else{
 949 //      fFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
 950 //      tDist = ::sqrt(tDu*tDu/fgMaxDuOuter/fgMaxDuOuter+
 951 //                   tDz*tDz/fgMaxDzOuter/fgMaxDzOuter);
 952 //      //fFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
 953 //       }
 954 //       if(tDist<tDistMax){
 955 //      fClosestRowAtDCA = ti+1;
 956 //      tDistMax = tDist;
 957 //       }
 958 //       fWeightedAvSep += tDist;
 959 //     }
 960 //   }
 961 //   if(tN>0){
 962 //     fWeightedAvSep /= tN;
 963 //     fFracOfMergedRow /= tN;
 964 //   }
 965 //   else{
 966 //     fClosestRowAtDCA = -1;
 967 //     fFracOfMergedRow = -1.;
 968 //     fWeightedAvSep = -1.;
 969 //   }
 970 // }
 971 // double AliFemtoPair::TpcExitSeparationTrackV0Pos() const {
 972 // //________________V0 daughters exit/entrance/average separation calc.
 973 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 974
 975 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->TpcV0PosExitPoint();
 976 //   return (diff.Mag());
 977 // }
 978
 979 // double AliFemtoPair::TpcEntranceSeparationTrackV0Pos() const {
 980 // //________________V0 daughters exit/entrance/average separation calc.
 981 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 982 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
 983 //   return (diff.Mag());
 984 // }
 985
 986 // double AliFemtoPair::TpcAverageSeparationTrackV0Pos() const {
 987 // //________________V0 daughters exit/entrance/average separation calc.
 988 // //_______1st part is a track 2nd is a V0 considering Pos daughter
 989 //   AliFemtoThreeVector diff;
 990 //   double tAveSep = 0.0;
 991 //   int ipt = 0;
 992 //   if (fTrack1->fNominalPosSample && fTrack2->fNominalPosSample){
 993 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
 994 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
 995 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
 996 //       fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
 997 //       fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
 998 //       fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
 999 //       (ipt<11)
1000 //       ){
1001 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
1002 //     ipt++;
1003 //     tAveSep += diff.Mag();
1004 //   }
1005 //   tAveSep = tAveSep/(ipt+1.);
1006 //   return (tAveSep);}
1007 //   else return -1;
1008 // }
1009 // double AliFemtoPair::TpcExitSeparationTrackV0Neg() const {
1010 // //_______1st part is a track 2nd is a V0 considering Neg daughter
1011 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcExitPoint() - fTrack2->TpcV0NegExitPoint();
1012 //   return (diff.Mag());
1013 // }
1014
1015 // double AliFemtoPair::TpcEntranceSeparationTrackV0Neg() const {
1016 // //_______1st part is a track 2nd is a V0 considering Neg daughter
1017 //   AliFemtoThreeVector diff = fTrack1->Track()->NominalTpcEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1018 //   return (diff.Mag());
1019 // }
1020
1021 // double AliFemtoPair::TpcAverageSeparationTrackV0Neg() const {
1022 // //_______1st part is a track 2nd is a V0 considering Neg daughter
1023 //   AliFemtoThreeVector diff;
1024 //   double tAveSep = 0.0;
1025 //   int ipt = 0;
1026 //   if (fTrack1->fNominalPosSample && fTrack2->fTpcV0NegPosSample){
1027 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
1028 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
1029 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
1030 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1031 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1032 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1033 //       (ipt<11)
1034 //       ){
1035 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1036 //     ipt++;
1037 //     tAveSep += diff.Mag();
1038 //   }
1039 //   tAveSep = tAveSep/(ipt+1.);
1040 //   return (tAveSep);}
1041 //   else return -1;
1042 // }
1043
1044 // double AliFemtoPair::TpcExitSeparationV0PosV0Pos() const {
1045 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1046 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosExitPoint() - fTrack2->TpcV0PosExitPoint();
1047 //   return (diff.Mag());
1048 // }
1049
1050 // double AliFemtoPair::TpcEntranceSeparationV0PosV0Pos() const {
1051 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1052 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
1053 //   return (diff.Mag());
1054 // }
1055 // double AliFemtoPair::TpcAverageSeparationV0PosV0Pos() const {
1056 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Pos daughter
1057 //   AliFemtoThreeVector diff;
1058 //   double tAveSep = 0.0;
1059 //   int ipt=0;
1060 //   if (fTrack1->fNominalPosSample && (fTrack2->fNominalPosSample)){
1061 //     while ((fabs(fTrack1->fNominalPosSample[ipt].x())<9999.) &&
1062 //      (fabs(fTrack1->fNominalPosSample[ipt].y())<9999.) &&
1063 //      (fabs(fTrack1->fNominalPosSample[ipt].z())<9999.) &&
1064 //      (fabs(fTrack2->fNominalPosSample[ipt].x())<9999.) &&
1065 //      (fabs(fTrack2->fNominalPosSample[ipt].y())<9999.) &&
1066 //      (fabs(fTrack2->fNominalPosSample[ipt].z())<9999.) &&
1067 //       (ipt<11)
1068 //      ){
1069 //       diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
1070 //       ipt++;
1071 //       tAveSep += diff.Mag();
1072 //     }
1073 //     tAveSep = tAveSep/(ipt+1);
1074 //     return (tAveSep);}
1075 //   else return -1;
1076 // }
1077
1078 // double AliFemtoPair::TpcExitSeparationV0PosV0Neg() const {
1079 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1080 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosExitPoint() - fTrack2->TpcV0NegExitPoint();
1081 //   return (diff.Mag());
1082 // }
1083
1084 // double AliFemtoPair::TpcEntranceSeparationV0PosV0Neg() const {
1085 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1086 //   AliFemtoThreeVector diff = fTrack1->TpcV0PosEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1087 //   return (diff.Mag());
1088 // }
1089 // double AliFemtoPair::TpcAverageSeparationV0PosV0Neg() const {
1090 // //_______1st part is a V0 considering Pos daughter 2nd is a V0 considering Neg daughter
1091 //   AliFemtoThreeVector diff;
1092 //   double tAveSep = 0.0;
1093 //   int ipt = 0;
1094 //   if (fTrack1->fNominalPosSample && fTrack2->fTpcV0NegPosSample){
1095 //   while (fabs(fTrack1->fNominalPosSample[ipt].x())<9999. &&
1096 //       fabs(fTrack1->fNominalPosSample[ipt].y())<9999. &&
1097 //       fabs(fTrack1->fNominalPosSample[ipt].z())<9999. &&
1098 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1099 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1100 //       fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1101 //       (ipt<11)
1102 //       ){
1103 //     diff = fTrack1->fNominalPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1104 //     ipt++;
1105 //     tAveSep += diff.Mag();
1106 //   }
1107 //   tAveSep = tAveSep/(ipt+1.);
1108 //   return (tAveSep);}
1109 //   else return -1;
1110 // }
1111 // double AliFemtoPair::TpcExitSeparationV0NegV0Pos() const {
1112 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1113 // // this is to check the upper case
1114 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegExitPoint() - fTrack2->TpcV0PosExitPoint();
1115 //   return (diff.Mag());
1116 // }
1117
1118 // double AliFemtoPair::TpcEntranceSeparationV0NegV0Pos() const {
1119 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1120 // // this is to check the upper case
1121 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegEntrancePoint() - fTrack2->TpcV0PosEntrancePoint();
1122 //   return (diff.Mag());
1123 // }
1124 // double AliFemtoPair::TpcAverageSeparationV0NegV0Pos() const {
1125 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Pos daughter
1126 // // this is to check the upper case
1127 //    AliFemtoThreeVector diff;
1128 //    double tAveSep = 0.0;
1129 //    int ipt = 0;
1130 //    if ( fTrack1->fTpcV0NegPosSample &&  fTrack2->fNominalPosSample){
1131 //      while (fabs(fTrack1->fTpcV0NegPosSample[ipt].x())<9999. &&
1132 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].y())<9999. &&
1133 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].z())<9999. &&
1134 //          fabs(fTrack2->fNominalPosSample[ipt].x())<9999. &&
1135 //          fabs(fTrack2->fNominalPosSample[ipt].y())<9999. &&
1136 //          fabs(fTrack2->fNominalPosSample[ipt].z())<9999. &&
1137 //          (ipt<11)
1138 //          ){
1139 //        diff = fTrack1->fTpcV0NegPosSample[ipt] - fTrack2->fNominalPosSample[ipt];
1140 //        ipt++;
1141 //        tAveSep += diff.Mag();
1142 //      }
1143 //      tAveSep = tAveSep/(ipt+1);
1144 //      return (tAveSep);}
1145 //      else return -1;
1146 // }
1147 // double AliFemtoPair::TpcExitSeparationV0NegV0Neg() const {
1148 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1149 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegExitPoint() - fTrack2->TpcV0NegExitPoint();
1150 //   return (diff.Mag());
1151 // }
1152
1153 // double AliFemtoPair::TpcEntranceSeparationV0NegV0Neg() const {
1154 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1155 //   AliFemtoThreeVector diff = fTrack1->TpcV0NegEntrancePoint() - fTrack2->TpcV0NegEntrancePoint();
1156 //   return (diff.Mag());
1157 // }
1158 // double AliFemtoPair::TpcAverageSeparationV0NegV0Neg() const {
1159 // //_______1st part is a V0 considering Neg daughter 2nd is a V0 considering Neg daughter
1160 //    AliFemtoThreeVector diff;
1161 //    double tAveSep = 0.0;
1162 //    int ipt=0;
1163 //    if (fTrack1->fTpcV0NegPosSample && fTrack2->fTpcV0NegPosSample){
1164 //      while (fabs(fTrack1->fTpcV0NegPosSample[ipt].x())<9999. &&
1165 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].y())<9999. &&
1166 //          fabs(fTrack1->fTpcV0NegPosSample[ipt].z())<9999. &&
1167 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].x())<9999. &&
1168 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].y())<9999. &&
1169 //          fabs(fTrack2->fTpcV0NegPosSample[ipt].z())<9999. &&
1170 //          (ipt<11)
1171 //          ){
1172 //        diff = fTrack1->fTpcV0NegPosSample[ipt] - fTrack2->fTpcV0NegPosSample[ipt];
1173 //        ipt++;
1174 //        tAveSep += diff.Mag();
1175 //      }
1176 //      tAveSep = tAveSep/(ipt+1);
1177 //      return (tAveSep);}
1178 //    else return -1;
1179 // }
1180
1181 // void AliFemtoPair::CalcMergingParFctn(short* tmpMergingParNotCalculatedFctn,
1182 //                                 float* tmpZ1,float* tmpU1,
1183 //                                 float* tmpZ2,float* tmpU2,
1184 //                                 int *tmpSect1,int *tmpSect2,
1185 //                                 double* tmpFracOfMergedRow,
1186 //                                 double* tmpClosestRowAtDCA
1187 //                                 ) const{
1188 // // calculate heper variables for merging
1189 //   tmpMergingParNotCalculatedFctn=0;
1190 //   double tDu, tDz;
1191 //   int tN = 0;
1192 //   *tmpFracOfMergedRow = 0.;
1193 //   *tmpClosestRowAtDCA = 0.;
1194 //   double tDist;
1195 //   double tDistMax = 100000000.;
1196 //   for(int ti=0 ; ti<45 ; ti++){
1197 //     if(tmpSect1[ti]==tmpSect2[ti] && tmpSect1[ti]!=-1){
1198 //      tDu = fabs(tmpU1[ti]-tmpU2[ti]);
1199 //      tDz = fabs(tmpZ1[ti]-tmpZ2[ti]);
1200 //      tN++;
1201 //       if(ti<13){
1202 //      *tmpFracOfMergedRow += (tDu<fgMaxDuInner && tDz<fgMaxDzInner);
1203 //      tDist = ::sqrt(tDu*tDu/fgMaxDuInner/fgMaxDuInner+
1204 //                   tDz*tDz/fgMaxDzInner/fgMaxDzInner);
1205 //       }
1206 //       else{
1207 //      *tmpFracOfMergedRow += (tDu<fgMaxDuOuter && tDz<fgMaxDzOuter);
1208 //      tDist = ::sqrt(tDu*tDu/fgMaxDuOuter/fgMaxDuOuter+
1209 //                   tDz*tDz/fgMaxDzOuter/fgMaxDzOuter);
1210 //      }
1211 //       if(tDist<tDistMax){
1212 //      fClosestRowAtDCA = ti+1;
1213 //      tDistMax = tDist;
1214 //       }
1215 //       //fWeightedAvSep += tDist; // now, wrong but not used
1216 //     }
1217 //   }
1218 //   if(tN>0){
1219 //     //fWeightedAvSep /= tN;
1220 //     *tmpFracOfMergedRow /= tN;
1221 //   }
1222 //   else{
1223 //     *tmpClosestRowAtDCA = -1;
1224 //     *tmpFracOfMergedRow = -1.;
1225 //     //fWeightedAvSep = -1.;
1226 //   }
1227 // }
1228