HBTAN/AliHBTPair.cxx

   1 #include "AliHBTPair.h"
   2 //_________________________________________________________________________
   3 ///////////////////////////////////////////////////////////////////////////
   4 //
   5 // class AliHBTPair
   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 "AliHBTParticle.h"
  15 #include "AliHBTWeights.h"
  16 #include "AliHBTTrackPoints.h"
  17
  18 ClassImp(AliHBTPair)
  19
  20 /************************************************************************/
  21 AliHBTPair::AliHBTPair(Bool_t rev):
  22  fPart1(0x0),
  23  fPart2(0x0),
  24  fSwapedPair(0x0),
  25  fQSideCMSLC(0.0),
  26  fQSideCMSLCNotCalc(kTRUE),
  27  fQOutCMSLC(0.0),
  28  fQOutCMSLCNotCalc(kTRUE),
  29  fQLongCMSLC(0.0),
  30  fQLongCMSLCNotCalc(kTRUE),
  31  fQInv(0.0),
  32  fQInvNotCalc(kTRUE),
  33  fInvMass(0.0),
  34  fInvMassNotCalc(kTRUE),
  35  fKt(0.0),
  36  fKtNotCalc(kTRUE),
  37  fKStar(0.0),
  38  fKStarNotCalc(kTRUE),
  39  fPInv(0.0),
  40  fQSide(0.0),
  41  fOut(0.0),
  42  fQLong(0.0),
  43  fMt(0.0),
  44  fMtNotCalc(kTRUE),
  45  fInvMassSqr(0.0),
  46  fMassSqrNotCalc(kTRUE),
  47  fQInvL(0.0),
  48  fQInvLNotCalc(kTRUE),
  49  fWeight(0.0),
  50  fWeightNotCalc(kTRUE),
  51  fAvarageDistance(0.0),
  52  fAvarageDistanceNotCalc(kTRUE),
  53  fPxSum(0.0),
  54  fPySum(0.0),
  55  fPzSum(0.0),
  56  fESum(0.0),
  57  fSumsNotCalc(kTRUE),
  58  fPxDiff(0.0),
  59  fPyDiff(0.0),
  60  fPzDiff(0.0),
  61  fEDiff(0.0),
  62  fDiffsNotCalc(kTRUE),
  63  fGammaCMSLC(0.0),
  64  fGammaCMSLCNotCalc(kTRUE),
  65  fChanged(kTRUE)
  66  {
  67 //value of rev defines if it is Swaped
  68 //if you pass kTRUE swpaped pair will NOT be created
  69 //though you wont be able to get the swaped pair from this pair
  70
  71   if(!rev) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
  72
  73  }
  74 /************************************************************************/
  75
  76 AliHBTPair::AliHBTPair(AliHBTParticle* part1, AliHBTParticle* part2, Bool_t rev):
  77  fPart1(part1),
  78  fPart2(part2),
  79  fSwapedPair(0x0),
  80  fQSideCMSLC(0.0),
  81  fQSideCMSLCNotCalc(kTRUE),
  82  fQOutCMSLC(0.0),
  83  fQOutCMSLCNotCalc(kTRUE),
  84  fQLongCMSLC(0.0),
  85  fQLongCMSLCNotCalc(kTRUE),
  86  fQInv(0.0),
  87  fQInvNotCalc(kTRUE),
  88  fInvMass(0.0),
  89  fInvMassNotCalc(kTRUE),
  90  fKt(0.0),
  91  fKtNotCalc(kTRUE),
  92  fKStar(0.0),
  93  fKStarNotCalc(kTRUE),
  94  fPInv(0.0),
  95  fQSide(0.0),
  96  fOut(0.0),
  97  fQLong(0.0),
  98  fMt(0.0),
  99  fMtNotCalc(kTRUE),
 100  fInvMassSqr(0.0),
 101  fMassSqrNotCalc(kTRUE),
 102  fQInvL(0.0),
 103  fQInvLNotCalc(kTRUE),
 104  fWeight(0.0),
 105  fWeightNotCalc(kTRUE),
 106  fAvarageDistance(0.0),
 107  fAvarageDistanceNotCalc(kTRUE),
 108  fPxSum(0.0),
 109  fPySum(0.0),
 110  fPzSum(0.0),
 111  fESum(0.0),
 112  fSumsNotCalc(kTRUE),
 113  fPxDiff(0.0),
 114  fPyDiff(0.0),
 115  fPzDiff(0.0),
 116  fEDiff(0.0),
 117  fDiffsNotCalc(kTRUE),
 118  fGammaCMSLC(0.0),
 119  fGammaCMSLCNotCalc(kTRUE),
 120  fChanged(kTRUE)
 121  {
 122 //value of rev defines if it is Swaped
 123 //if you pass kTRUE swpaped pair will NOT be created
 124 //though you wont be able to get the swaped pair from this pair
 125
 126   if(!rev) fSwapedPair = new AliHBTPair(part2,part1,kTRUE); //if false create swaped pair
 127
 128  }
 129 /************************************************************************/
 130 AliHBTPair::AliHBTPair(const AliHBTPair& in):
 131  TObject(in),
 132  fPart1(0x0),
 133  fPart2(0x0),
 134  fSwapedPair(0x0),
 135  fQSideCMSLC(0.0),
 136  fQSideCMSLCNotCalc(kTRUE),
 137  fQOutCMSLC(0.0),
 138  fQOutCMSLCNotCalc(kTRUE),
 139  fQLongCMSLC(0.0),
 140  fQLongCMSLCNotCalc(kTRUE),
 141  fQInv(0.0),
 142  fQInvNotCalc(kTRUE),
 143  fInvMass(0.0),
 144  fInvMassNotCalc(kTRUE),
 145  fKt(0.0),
 146  fKtNotCalc(kTRUE),
 147  fKStar(0.0),
 148  fKStarNotCalc(kTRUE),
 149  fPInv(0.0),
 150  fQSide(0.0),
 151  fOut(0.0),
 152  fQLong(0.0),
 153  fMt(0.0),
 154  fMtNotCalc(kTRUE),
 155  fInvMassSqr(0.0),
 156  fMassSqrNotCalc(kTRUE),
 157  fQInvL(0.0),
 158  fQInvLNotCalc(kTRUE),
 159  fWeight(0.0),
 160  fWeightNotCalc(kTRUE),
 161  fAvarageDistance(0.0),
 162  fAvarageDistanceNotCalc(kTRUE),
 163  fPxSum(0.0),
 164  fPySum(0.0),
 165  fPzSum(0.0),
 166  fESum(0.0),
 167  fSumsNotCalc(kTRUE),
 168  fPxDiff(0.0),
 169  fPyDiff(0.0),
 170  fPzDiff(0.0),
 171  fEDiff(0.0),
 172  fDiffsNotCalc(kTRUE),
 173  fGammaCMSLC(0.0),
 174  fGammaCMSLCNotCalc(kTRUE),
 175  fChanged(kTRUE)
 176 {
 177  //cpy constructor
 178  in.Copy(*this);
 179 }
 180 /************************************************************************/
 181
 182 AliHBTPair& AliHBTPair::operator=(const AliHBTPair& in)
 183 {
 184  //Assigment operator
 185  in.Copy(*this);
 186  return *this;
 187 }
 188 /************************************************************************/
 189
 190 Double_t AliHBTPair::GetInvMass()
 191 {
 192 //Returns qinv value for a pair
 193   if(fInvMassNotCalc)
 194    {
 195      CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
 196
 197      if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
 198      else fInvMass = TMath::Sqrt(fInvMassSqr);
 199
 200      fInvMassNotCalc = kFALSE;
 201    }
 202   return fInvMass;
 203 }
 204 /************************************************************************/
 205
 206 Double_t AliHBTPair::GetQSideCMSLC()
 207 {
 208 //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 209
 210   if (fQSideCMSLCNotCalc)
 211    {
 212     fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
 213     fQSideCMSLCNotCalc = kFALSE;
 214    }
 215   return fQSideCMSLC;
 216 }
 217 /************************************************************************/
 218
 219 Double_t AliHBTPair::GetQOutCMSLC()
 220 {
 221  //caculates Qout in Center Of Mass Longitudionally Co-Moving
 222  if(fQOutCMSLCNotCalc)
 223   {
 224    CalculateSums();
 225    CalculateDiffs();
 226
 227    if (fPart1->GetMass() != fPart2->GetMass())
 228     {
 229 /*
 230       //STAR algorithm
 231       Double_t beta  = fPzSum/fESum;
 232       Double_t gamma = GetGammaToCMSLC();
 233       Double_t el = gamma * (fPart1->Energy() - beta * fPart1->Pz());
 234       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
 235       beta  = 2.0*GetKt()/GetMt();
 236       gamma = GetMt()/GetQInv();
 237       fQOutCMSLC = gamma * (x - beta * el);
 238 */
 239
 240       //beta=fPzSum/fESum;    // Longit. V == beta
 241       Double_t beta=fPzSum/fESum;
 242       Double_t gamma = GetGammaToCMSLC();
 243
 244       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
 245       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
 246
 247 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
 248 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
 249       Double_t tmp;
 250       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
 251       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
 252       Double_t part1Px = tmp;
 253
 254       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
 255       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
 256       Double_t part2Px = tmp;
 257
 258
 259 //      LTR(part1Pz,E1,beta,GetGammaToCMSLC(),part1Pz,E1a);
 260 //      LTR(part2Pz,E2,beta,GetGammaToCMSLC(),part2Pz,E2a);
 261       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->Energy());
 262       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->Energy());
 263
 264       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
 265       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
 266       Double_t part1E=TMath::Sqrt(fPart1->GetMass()*fPart1->GetMass()+part1P2);
 267       Double_t part2E=TMath::Sqrt(fPart2->GetMass()*fPart2->GetMass()+part2P2);
 268       Double_t sumE=part1E+part2E;
 269       Double_t sumPx=part1Px+part2Px;
 270       Double_t sumPy=part1Py+part2Py;
 271       Double_t sumPZ=part1Pz+part2Pz;
 272       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
 273
 274       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
 275       Double_t hf = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/(relmass*relmass);
 276       fQOutCMSLC=(part1Px-part2Px);//== id
 277       fQOutCMSLC=fQOutCMSLC-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
 278     }
 279    else
 280     {
 281       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
 282       fQOutCMSLC = 0.5*k2/GetKt();
 283    // if (non-id)  fQOutCMSLC=fQOutCMSLC - sumPx*HF;
 284     }
 285
 286
 287    fQOutCMSLCNotCalc = kFALSE;
 288   }
 289  return fQOutCMSLC;
 290 }
 291 /************************************************************************/
 292
 293 Double_t AliHBTPair::GetQLongCMSLC()
 294 {
 295  //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
 296  if (fQLongCMSLCNotCalc)
 297   {
 298     CalculateSums();
 299     CalculateDiffs();
 300     Double_t beta = fPzSum/fESum;
 301     fQLongCMSLC = GetGammaToCMSLC() * ( fPzDiff - beta*fEDiff );
 302     fQLongCMSLCNotCalc = kFALSE;
 303   }
 304  return fQLongCMSLC;
 305 }
 306 /************************************************************************/
 307
 308 Double_t AliHBTPair::GetKt()
 309 {
 310  //calculates the evarage momentum of the pair
 311   if(fKtNotCalc)
 312    {
 313      CalculateSums();
 314      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
 315      fKtNotCalc = kFALSE;
 316    }
 317   return fKt;
 318 }
 319 /************************************************************************/
 320
 321 Double_t AliHBTPair::GetKStar()
 322 {
 323   //calculates invariant velocity difference
 324   if (fKStarNotCalc)
 325    {
 326     CalculateSums();
 327
 328     Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
 329     Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
 330     Double_t pinv =   TMath::Sqrt(mtrans - ptrans);
 331
 332     Double_t q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/pinv;
 333
 334     CalculateQInvL();
 335
 336     q = q*q - fQInvL;
 337     if ( q < 0)
 338      {
 339         Info("GetKStar","q = %f",q);
 340         fPart1->Print();
 341         fPart2->Print();
 342         q = TMath::Abs(q);
 343      }
 344
 345     q = TMath::Sqrt(q);
 346     fKStar = q/2.;
 347     fKStarNotCalc = kFALSE;
 348    }
 349   return fKStar;
 350 }
 351 /************************************************************************/
 352
 353 Double_t AliHBTPair::GetQInv()
 354 {
 355 //returns Qinv
 356 //warning for non-id particles you want to use 2*KStar
 357   if(fQInvNotCalc)
 358    {
 359     CalculateQInvL();
 360     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
 361     fQInvNotCalc = kFALSE;
 362    }
 363   return fQInv;
 364 }
 365 /************************************************************************/
 366
 367 Double_t AliHBTPair::GetGammaToCMSLC()
 368 {
 369   //calculates gamma factor of the boost to CMSLC
 370   if(fGammaCMSLCNotCalc)
 371    {
 372      CalculateSums();
 373      Double_t beta = fPzSum/fESum;
 374      fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
 375      fGammaCMSLCNotCalc = kFALSE;
 376    }
 377   return fGammaCMSLC;
 378 }
 379 /************************************************************************/
 380
 381 Double_t AliHBTPair::GetMt()
 382 {
 383   //Calculates transverse mass of the pair
 384   if (fMtNotCalc)
 385    {
 386      CalculateSums();
 387      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
 388      fMtNotCalc = kFALSE;
 389    }
 390    return fMt;
 391 }
 392 /************************************************************************/
 393
 394 Double_t AliHBTPair::GetWeight()
 395 {
 396   //returns and buffers weight for this pair
 397   if (fWeightNotCalc)
 398    {
 399       fWeight = AliHBTWeights::Weight(this);
 400       fWeightNotCalc = kFALSE;
 401    }
 402   return fWeight;
 403 }
 404 /************************************************************************/
 405
 406 Double_t AliHBTPair::GetAvarageDistance()
 407 {
 408 //returns and buffers avarage distance between two tracks calculated
 409 // out of track points (see AliHBTTrackPoints class)
 410
 411   if (fAvarageDistanceNotCalc)
 412    {
 413      fAvarageDistance = AvDistance();
 414      fAvarageDistanceNotCalc = kFALSE;
 415    }
 416   return fAvarageDistance;
 417 }
 418 /************************************************************************/
 419
 420 Double_t AliHBTPair::AvDistance()
 421 {
 422   //returns avarage distance between two tracks in range
 423   //as defined in Track-Points of AliHBTParticle
 424   //returns negative value if error uccured f.g. tracks do not have track-points
 425   AliHBTTrackPoints* tpts1 = fPart1->GetTrackPoints();
 426   if ( tpts1 == 0x0)
 427    {//it could be simulated pair
 428 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
 429      return -1.0;
 430    }
 431
 432   AliHBTTrackPoints* tpts2 = fPart2->GetTrackPoints();
 433   if ( tpts2 == 0x0)
 434    {
 435 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
 436      return -1.0;
 437    }
 438
 439   return tpts1->AvarageDistance(*tpts2);
 440 }