ANALYSIS/AliAODPair.cxx

   1 #include "AliAODPair.h"
   2 //_________________________________________________________________________
   3 ///////////////////////////////////////////////////////////////////////////
   4 //
   5 // class AliAODPair
   6 //
   7 // class implements pair of particles and taking care of caluclation (almost)
   8 // all of pair properties (Qinv, InvMass,...)
   9 //
  10 // more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
  11 //
  12 ////////////////////////////////////////////////////////////////////////////
  13
  14 #include "AliVAODParticle.h"
  15 #include "AliTrackPoints.h"
  16
  17 ClassImp(AliAODPair)
  18
  19 /************************************************************************/
  20 AliAODPair::AliAODPair(Bool_t rev):
  21  fPart1(0x0),
  22  fPart2(0x0),
  23  fSwappedPair(0x0),
  24  fQSideLCMS(0.0),
  25  fQSideLCMSNotCalc(kTRUE),
  26  fQOutLCMS(0.0),
  27  fQOutLCMSNotCalc(kTRUE),
  28  fQLongLCMS(0.0),
  29  fQLongLCMSNotCalc(kTRUE),
  30  fQtLCMS(0.0),
  31  fQtLCMSNotCalc(kTRUE),
  32  fQt(0.0),
  33  fQtNotCalc(kTRUE),
  34  fQInv(0.0),
  35  fQInvNotCalc(kTRUE),
  36  fInvMass(0.0),
  37  fInvMassNotCalc(kTRUE),
  38  fKt(0.0),
  39  fKtNotCalc(kTRUE),
  40  fKStar(0.0),
  41  fKStarNotCalc(kTRUE),
  42  fKStarOut(0.0),
  43  fKStarSide(0.0),
  44  fKStarLong(0.0),
  45  fKStarCompNotCalc(kTRUE),
  46  fPInv(0.0),
  47  fQSide(0.0),
  48  fOut(0.0),
  49  fQLong(0.0),
  50  fMt(0.0),
  51  fMtNotCalc(kTRUE),
  52  fInvMassSqr(0.0),
  53  fMassSqrNotCalc(kTRUE),
  54  fQInvL(0.0),
  55  fQInvLNotCalc(kTRUE),
  56  fAvarageDistance(0.0),
  57  fAvarageDistanceNotCalc(kTRUE),
  58  fPxSum(0.0),
  59  fPySum(0.0),
  60  fPzSum(0.0),
  61  fESum(0.0),
  62  fSumsNotCalc(kTRUE),
  63  fPxDiff(0.0),
  64  fPyDiff(0.0),
  65  fPzDiff(0.0),
  66  fEDiff(0.0),
  67  fDiffsNotCalc(kTRUE),
  68  fGammaLCMS(0.0),
  69  fGammaLCMSNotCalc(kTRUE),
  70  fChanged(kTRUE)
  71  {
  72 //value of rev defines if it is Swapped
  73 //if you pass kTRUE swpaped pair will NOT be created
  74 //though you wont be able to get the swaped pair from this pair
  75
  76   if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
  77
  78  }
  79 /************************************************************************/
  80
  81 AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
  82  fPart1(part1),
  83  fPart2(part2),
  84  fSwappedPair(0x0),
  85  fQSideLCMS(0.0),
  86  fQSideLCMSNotCalc(kTRUE),
  87  fQOutLCMS(0.0),
  88  fQOutLCMSNotCalc(kTRUE),
  89  fQLongLCMS(0.0),
  90  fQLongLCMSNotCalc(kTRUE),
  91  fQtLCMS(0.0),
  92  fQtLCMSNotCalc(kTRUE),
  93  fQInv(0.0),
  94  fQInvNotCalc(kTRUE),
  95  fInvMass(0.0),
  96  fInvMassNotCalc(kTRUE),
  97  fKt(0.0),
  98  fKtNotCalc(kTRUE),
  99  fKStar(0.0),
 100  fKStarNotCalc(kTRUE),
 101  fKStarOut(0.0),
 102  fKStarSide(0.0),
 103  fKStarLong(0.0),
 104  fKStarCompNotCalc(kTRUE),
 105  fPInv(0.0),
 106  fQSide(0.0),
 107  fOut(0.0),
 108  fQLong(0.0),
 109  fMt(0.0),
 110  fMtNotCalc(kTRUE),
 111  fInvMassSqr(0.0),
 112  fMassSqrNotCalc(kTRUE),
 113  fQInvL(0.0),
 114  fQInvLNotCalc(kTRUE),
 115  fAvarageDistance(0.0),
 116  fAvarageDistanceNotCalc(kTRUE),
 117  fPxSum(0.0),
 118  fPySum(0.0),
 119  fPzSum(0.0),
 120  fESum(0.0),
 121  fSumsNotCalc(kTRUE),
 122  fPxDiff(0.0),
 123  fPyDiff(0.0),
 124  fPzDiff(0.0),
 125  fEDiff(0.0),
 126  fDiffsNotCalc(kTRUE),
 127  fGammaLCMS(0.0),
 128  fGammaLCMSNotCalc(kTRUE),
 129  fChanged(kTRUE)
 130  {
 131 //value of rev defines if it is Swapped
 132 //if you pass kTRUE swpaped pair will NOT be created
 133 //though you wont be able to get the swaped pair from this pair
 134
 135   if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
 136
 137  }
 138 /************************************************************************/
 139 AliAODPair::AliAODPair(const AliAODPair& in):
 140  TObject(in),
 141  fPart1(0x0),
 142  fPart2(0x0),
 143  fSwappedPair(0x0),
 144  fQSideLCMS(0.0),
 145  fQSideLCMSNotCalc(kTRUE),
 146  fQOutLCMS(0.0),
 147  fQOutLCMSNotCalc(kTRUE),
 148  fQLongLCMS(0.0),
 149  fQLongLCMSNotCalc(kTRUE),
 150  fQtLCMS(0.0),
 151  fQtLCMSNotCalc(kTRUE),
 152  fQInv(0.0),
 153  fQInvNotCalc(kTRUE),
 154  fInvMass(0.0),
 155  fInvMassNotCalc(kTRUE),
 156  fKt(0.0),
 157  fKtNotCalc(kTRUE),
 158  fKStar(0.0),
 159  fKStarNotCalc(kTRUE),
 160  fKStarOut(0.0),
 161  fKStarSide(0.0),
 162  fKStarLong(0.0),
 163  fKStarCompNotCalc(kTRUE),
 164  fPInv(0.0),
 165  fQSide(0.0),
 166  fOut(0.0),
 167  fQLong(0.0),
 168  fMt(0.0),
 169  fMtNotCalc(kTRUE),
 170  fInvMassSqr(0.0),
 171  fMassSqrNotCalc(kTRUE),
 172  fQInvL(0.0),
 173  fQInvLNotCalc(kTRUE),
 174  fAvarageDistance(0.0),
 175  fAvarageDistanceNotCalc(kTRUE),
 176  fPxSum(0.0),
 177  fPySum(0.0),
 178  fPzSum(0.0),
 179  fESum(0.0),
 180  fSumsNotCalc(kTRUE),
 181  fPxDiff(0.0),
 182  fPyDiff(0.0),
 183  fPzDiff(0.0),
 184  fEDiff(0.0),
 185  fDiffsNotCalc(kTRUE),
 186  fGammaLCMS(0.0),
 187  fGammaLCMSNotCalc(kTRUE),
 188  fChanged(kTRUE)
 189 {
 190  //cpy constructor
 191  in.Copy(*this);
 192 }
 193 /************************************************************************/
 194
 195 AliAODPair& AliAODPair::operator=(const AliAODPair& in)
 196 {
 197  //Assigment operator
 198  in.Copy(*this);
 199  return *this;
 200 }
 201 /************************************************************************/
 202
 203 Double_t AliAODPair::GetInvMass()
 204 {
 205 //Returns qinv value for a pair
 206   if(fInvMassNotCalc)
 207    {
 208      CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
 209
 210      if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
 211      else fInvMass = TMath::Sqrt(fInvMassSqr);
 212
 213      fInvMassNotCalc = kFALSE;
 214    }
 215   return fInvMass;
 216 }
 217 /************************************************************************/
 218
 219 Double_t AliAODPair::GetQSideLCMS()
 220 {
 221 //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 222
 223   if (fQSideLCMSNotCalc)
 224    {
 225     fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
 226     fQSideLCMSNotCalc = kFALSE;
 227    }
 228   return fQSideLCMS;
 229 }
 230 /************************************************************************/
 231
 232 Double_t AliAODPair::GetQOutLCMS()
 233 {
 234  //caculates Qout in Center Of Mass Longitudionally Co-Moving
 235  if(fQOutLCMSNotCalc)
 236   {
 237    CalculateSums();
 238    CalculateDiffs();
 239
 240    if (fPart1->Mass() != fPart2->Mass())
 241     {
 242 /*
 243       //STAR algorithm
 244       Double_t beta  = fPzSum/fESum;
 245       Double_t gamma = GetGammaToLCMS();
 246       Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
 247       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
 248       beta  = 2.0*GetKt()/GetMt();
 249       gamma = GetMt()/GetQInv();
 250       fQOutLCMS = gamma * (x - beta * el);
 251 */
 252
 253       //beta=fPzSum/fESum;    // Longit. V == beta
 254       Double_t beta=fPzSum/fESum;
 255       Double_t gamma = GetGammaToLCMS();
 256
 257       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
 258       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
 259
 260 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
 261 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
 262       Double_t tmp;
 263       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
 264       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
 265       Double_t part1Px = tmp;
 266
 267       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
 268       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
 269       Double_t part2Px = tmp;
 270
 271
 272 //      LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
 273 //      LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
 274       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
 275       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
 276
 277       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
 278       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
 279       Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
 280       Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
 281       Double_t sumE=part1E+part2E;
 282       Double_t sumPx=part1Px+part2Px;
 283       Double_t sumPy=part1Py+part2Py;
 284       Double_t sumPZ=part1Pz+part2Pz;
 285       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
 286
 287       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
 288       Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
 289       fQOutLCMS=(part1Px-part2Px);//== id
 290       fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
 291     }
 292    else
 293     {
 294       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
 295       fQOutLCMS = 0.5*k2/GetKt();
 296    // if (non-id)  fQOutLCMS=fQOutLCMS - sumPx*HF;
 297     }
 298
 299
 300    fQOutLCMSNotCalc = kFALSE;
 301   }
 302  return fQOutLCMS;
 303 }
 304 /************************************************************************/
 305
 306 Double_t AliAODPair::GetQLongLCMS()
 307 {
 308  //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
 309  if (fQLongLCMSNotCalc)
 310   {
 311     CalculateSums();
 312     CalculateDiffs();
 313     Double_t beta = fPzSum/fESum;
 314     fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
 315     fQLongLCMSNotCalc = kFALSE;
 316   }
 317  return fQLongLCMS;
 318 }
 319 /************************************************************************/
 320
 321 Double_t AliAODPair::GetQtLCMS()
 322 {
 323  //returns Q transverse CMS longitudionally co-moving
 324  if (fQtLCMSNotCalc)
 325   {
 326     fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
 327     fQtLCMSNotCalc = kFALSE;
 328   }
 329  return fQtLCMS;
 330 }
 331 /************************************************************************/
 332
 333 Double_t AliAODPair::GetQt()
 334 {
 335  //returns Q transverse CMS longitudionally co-moving
 336  if (fQtNotCalc)
 337   {
 338     CalculateSums();
 339     CalculateDiffs();
 340
 341     Double_t dotprod = fPxSum*fPxDiff + fPySum*fPyDiff + fPzSum*fPzDiff;
 342     Double_t klen =    fPxSum*fPxSum  + fPySum*fPySum  + fPzSum*fPzSum;
 343     klen = TMath::Sqrt(klen);
 344     Double_t qlen =    fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff;
 345     qlen = TMath::Sqrt(qlen);
 346     Double_t tmp = klen*qlen;
 347     if (tmp == 0.0)
 348      {
 349        fQt = 10e5;
 350        fQtNotCalc = kFALSE;
 351        return fQt;
 352      }
 353     Double_t cosopenangle = dotprod/tmp;
 354     Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
 355
 356     fQt = sinopenangle*qlen;
 357     fQtNotCalc = kFALSE;
 358   }
 359  return fQt;
 360 }
 361 /************************************************************************/
 362
 363 Double_t AliAODPair::GetKt()
 364 {
 365  //calculates the evarage momentum of the pair
 366   if(fKtNotCalc)
 367    {
 368      CalculateSums();
 369      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
 370      fKtNotCalc = kFALSE;
 371    }
 372   return fKt;
 373 }
 374 /************************************************************************/
 375
 376 Double_t AliAODPair::GetKStar()
 377 {
 378   //calculates invariant velocity difference
 379   if (fKStarNotCalc)
 380    {
 381     CalculateSums();
 382
 383     Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
 384     Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
 385     if (ptrans > mtrans)
 386      {
 387        Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
 388        Error("GetKStar","Particle1:");
 389        fPart1->Print();
 390        Error("GetKStar","Particle2:");
 391        fPart2->Print();
 392        Error("GetKStar","");
 393
 394        fKStar = 10e5;
 395        fKStarNotCalc = kFALSE;
 396        return fKStar;
 397      }
 398     Double_t pinv =   TMath::Sqrt(mtrans - ptrans);
 399
 400     Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
 401
 402     CalculateQInvL();
 403
 404     q = q*q - fQInvL;
 405     if ( q < 0)
 406      {
 407        Info("GetKStar","Sqrt of negative number q = %f",q);
 408        Error("GetKStar","Particle1:");
 409        fPart1->Print();
 410        Error("GetKStar","Particle2:");
 411        fPart2->Print();
 412        fKStar = 10e5;
 413        fKStarNotCalc = kFALSE;
 414        return fKStar;
 415      }
 416
 417     q = TMath::Sqrt(q);
 418     fKStar = q/2.;
 419     fKStarNotCalc = kFALSE;
 420    }
 421   return fKStar;
 422 }
 423 /************************************************************************/
 424 Double_t AliAODPair::GetKStarOut()
 425 {
 426     CalculateKStarComp();
 427     return fKStarOut;
 428 }
 429 /************************************************************************/
 430 Double_t AliAODPair::GetKStarSide()
 431 {
 432     CalculateKStarComp();
 433     return fKStarSide;
 434 }
 435 /************************************************************************/
 436 Double_t AliAODPair::GetKStarLong()
 437 {
 438     CalculateKStarComp();
 439     return fKStarLong;
 440 }
 441 /************************************************************************/
 442
 443 Double_t AliAODPair::GetQInv()
 444 {
 445 //returns Qinv
 446 //warning for non-id particles you want to use 2*KStar
 447   if(fQInvNotCalc)
 448    {
 449     CalculateQInvL();
 450     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
 451     fQInvNotCalc = kFALSE;
 452    }
 453   return fQInv;
 454 }
 455 /************************************************************************/
 456
 457 Double_t AliAODPair::GetGammaToLCMS()
 458 {
 459   //calculates gamma factor of the boost to LCMS
 460   if(fGammaLCMSNotCalc)
 461    {
 462      CalculateSums();
 463      Double_t beta = fPzSum/fESum;
 464      fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
 465      fGammaLCMSNotCalc = kFALSE;
 466    }
 467   return fGammaLCMS;
 468 }
 469 /************************************************************************/
 470
 471 Double_t AliAODPair::GetMt()
 472 {
 473   //Calculates transverse mass of the pair
 474   if (fMtNotCalc)
 475    {
 476      CalculateSums();
 477      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
 478      fMtNotCalc = kFALSE;
 479    }
 480    return fMt;
 481 }
 482 /************************************************************************/
 483
 484 Double_t AliAODPair::GetAvarageDistance()
 485 {
 486 //returns and buffers avarage distance between two tracks calculated
 487 // out of track points (see AliAODTrackPoints class)
 488
 489   if (fAvarageDistanceNotCalc)
 490    {
 491      fAvarageDistance = AvDistance();
 492      fAvarageDistanceNotCalc = kFALSE;
 493    }
 494   return fAvarageDistance;
 495 }
 496 /************************************************************************/
 497
 498 Double_t AliAODPair::AvDistance()
 499 {
 500   //returns avarage distance between two tracks in range
 501   //as defined in Track-Points of AliVAODParticle
 502   //returns negative value if error uccured f.g. tracks do not have track-points
 503   AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
 504   if ( tpts1 == 0x0)
 505    {//it could be simulated pair
 506 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
 507      return -1.0;
 508    }
 509
 510   AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
 511   if ( tpts2 == 0x0)
 512    {
 513 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
 514      return -1.0;
 515    }
 516
 517   return tpts1->AvarageDistance(*tpts2);
 518 }