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