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