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 AliHBTPair& AliHBTPair::operator=(const AliHBTPair& in)
 176 {
 177  //Assigment operator
 178  in.Copy(*this);
 179  return *this;
 180 }
 181 /************************************************************************/
 182
 183 Double_t AliHBTPair::GetInvMass()
 184 {
 185 //Returns qinv value for a pair
 186   if(fInvMassNotCalc)
 187    {
 188      CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
 189
 190      if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
 191      else fInvMass = TMath::Sqrt(fInvMassSqr);
 192
 193      fInvMassNotCalc = kFALSE;
 194    }
 195   return fInvMass;
 196 }
 197 /************************************************************************/
 198 Double_t AliHBTPair::GetQSideCMSLC()
 199 {
 200   //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 201
 202   if (fQSideCMSLCNotCalc)
 203    {
 204     fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
 205     fQSideCMSLCNotCalc = kFALSE;
 206    }
 207   return fQSideCMSLC;
 208 }
 209 /************************************************************************/
 210 Double_t AliHBTPair::GetQOutCMSLC()
 211 {
 212  //caculates Qout in Center Of Mass Longitudionally Co-Moving
 213  if(fQOutCMSLCNotCalc)
 214   {
 215    CalculateSums();
 216    CalculateDiffs();
 217
 218    if (fPart1->GetMass() != fPart2->GetMass())
 219     {
 220 /*
 221       //STAR algorithm
 222       Double_t beta  = fPzSum/fESum;
 223       Double_t gamma = GetGammaToCMSLC();
 224       Double_t el = gamma * (fPart1->Energy() - beta * fPart1->Pz());
 225       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
 226       beta  = 2.0*GetKt()/GetMt();
 227       gamma = GetMt()/GetQInv();
 228       fQOutCMSLC = gamma * (x - beta * el);
 229 */
 230
 231       //beta=fPzSum/fESum;    // Longit. V == beta
 232       Double_t beta=fPzSum/fESum;
 233       Double_t gamma = GetGammaToCMSLC();
 234
 235       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
 236       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
 237
 238 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
 239 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
 240       Double_t tmp;
 241       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
 242       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
 243       Double_t part1Px = tmp;
 244
 245       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
 246       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
 247       Double_t part2Px = tmp;
 248
 249
 250 //      LTR(part1Pz,E1,beta,GetGammaToCMSLC(),part1Pz,E1a);
 251 //      LTR(part2Pz,E2,beta,GetGammaToCMSLC(),part2Pz,E2a);
 252       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->Energy());
 253       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->Energy());
 254
 255       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
 256       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
 257       Double_t part1E=TMath::Sqrt(fPart1->GetMass()*fPart1->GetMass()+part1P2);
 258       Double_t part2E=TMath::Sqrt(fPart2->GetMass()*fPart2->GetMass()+part2P2);
 259       Double_t sumE=part1E+part2E;
 260       Double_t sumPx=part1Px+part2Px;
 261       Double_t sumPy=part1Py+part2Py;
 262       Double_t sumPZ=part1Pz+part2Pz;
 263       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
 264
 265       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
 266       Double_t hf = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/(relmass*relmass);
 267       fQOutCMSLC=(part1Px-part2Px);//== id
 268       fQOutCMSLC=fQOutCMSLC-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
 269     }
 270    else
 271     {
 272       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
 273       fQOutCMSLC = 0.5*k2/GetKt();
 274    // if (non-id)  fQOutCMSLC=fQOutCMSLC - sumPx*HF;
 275     }
 276
 277
 278    fQOutCMSLCNotCalc = kFALSE;
 279   }
 280  return fQOutCMSLC;
 281 }
 282 /************************************************************************/
 283 Double_t AliHBTPair::GetQLongCMSLC()
 284 {
 285  //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
 286  if (fQLongCMSLCNotCalc)
 287   {
 288     CalculateSums();
 289     CalculateDiffs();
 290     Double_t beta = fPzSum/fESum;
 291     fQLongCMSLC = GetGammaToCMSLC() * ( fPzDiff - beta*fEDiff );
 292     fQLongCMSLCNotCalc = kFALSE;
 293   }
 294  return fQLongCMSLC;
 295 }
 296 /************************************************************************/
 297 Double_t AliHBTPair::GetKt()
 298 {
 299  //calculates the evarage momentum of the pair
 300   if(fKtNotCalc)
 301    {
 302      CalculateSums();
 303      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
 304      fKtNotCalc = kFALSE;
 305    }
 306   return fKt;
 307 }
 308 /************************************************************************/
 309
 310 Double_t AliHBTPair::GetKStar()
 311 {
 312   //calculates invariant velocity difference
 313   if (fKStarNotCalc)
 314    {
 315     CalculateSums();
 316
 317     Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
 318     Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
 319     Double_t pinv =   TMath::Sqrt(mtrans - ptrans);
 320
 321     Double_t q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/pinv;
 322
 323     CalculateQInvL();
 324
 325     q = q*q - fQInvL;
 326     if ( q < 0)
 327      {
 328         Info("GetKStar","q = %f",q);
 329         fPart1->Print();
 330         fPart2->Print();
 331         q = TMath::Abs(q);
 332      }
 333
 334     q = TMath::Sqrt(q);
 335     fKStar = q/2.;
 336     fKStarNotCalc = kFALSE;
 337    }
 338   return fKStar;
 339 }
 340 /************************************************************************/
 341
 342 Double_t AliHBTPair::GetQInv()
 343 {
 344 //returns Qinv
 345 //warning for non-id particles you want to use 2*KStar
 346   if(fQInvNotCalc)
 347    {
 348     CalculateQInvL();
 349     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
 350     fQInvNotCalc = kFALSE;
 351    }
 352   return fQInv;
 353 }
 354 /************************************************************************/
 355
 356 Double_t AliHBTPair::GetGammaToCMSLC()
 357 {
 358   //calculates gamma factor of the boost to CMSLC
 359   if(fGammaCMSLCNotCalc)
 360    {
 361      CalculateSums();
 362      Double_t beta = fPzSum/fESum;
 363      fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
 364      fGammaCMSLCNotCalc = kFALSE;
 365    }
 366   return fGammaCMSLC;
 367 }
 368 /************************************************************************/
 369
 370 Double_t AliHBTPair::GetMt()
 371 {
 372   //Calculates transverse mass of the pair
 373   if (fMtNotCalc)
 374    {
 375      CalculateSums();
 376      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
 377      fMtNotCalc = kFALSE;
 378    }
 379    return fMt;
 380 }
 381 /************************************************************************/
 382
 383 Double_t AliHBTPair::GetWeight()
 384 {
 385   //returns and buffers weight for this pair
 386   if (fWeightNotCalc)
 387    {
 388       fWeight = AliHBTWeights::Weight(this);
 389       fWeightNotCalc = kFALSE;
 390    }
 391   return fWeight;
 392 }