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     if (ptrans > mtrans)
 331      {
 332        Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
 333        Error("GetKStar","Particle1:");
 334        fPart1->Print();
 335        Error("GetKStar","Particle2:");
 336        fPart2->Print();
 337        Error("GetKStar","");
 338
 339        fKStar = 10e5;
 340        fKStarNotCalc = kFALSE;
 341        return fKStar;
 342      }
 343     Double_t pinv =   TMath::Sqrt(mtrans - ptrans);
 344
 345     Double_t q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/pinv;
 346
 347     CalculateQInvL();
 348
 349     q = q*q - fQInvL;
 350     if ( q < 0)
 351      {
 352        Info("GetKStar","Sqrt of negative number q = %f",q);
 353        Error("GetKStar","Particle1:");
 354        fPart1->Print();
 355        Error("GetKStar","Particle2:");
 356        fPart2->Print();
 357        fKStar = 10e5;
 358        fKStarNotCalc = kFALSE;
 359        return fKStar;
 360      }
 361
 362     q = TMath::Sqrt(q);
 363     fKStar = q/2.;
 364     fKStarNotCalc = kFALSE;
 365    }
 366   return fKStar;
 367 }
 368 /************************************************************************/
 369
 370 Double_t AliHBTPair::GetQInv()
 371 {
 372 //returns Qinv
 373 //warning for non-id particles you want to use 2*KStar
 374   if(fQInvNotCalc)
 375    {
 376     CalculateQInvL();
 377     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
 378     fQInvNotCalc = kFALSE;
 379    }
 380   return fQInv;
 381 }
 382 /************************************************************************/
 383
 384 Double_t AliHBTPair::GetGammaToCMSLC()
 385 {
 386   //calculates gamma factor of the boost to CMSLC
 387   if(fGammaCMSLCNotCalc)
 388    {
 389      CalculateSums();
 390      Double_t beta = fPzSum/fESum;
 391      fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
 392      fGammaCMSLCNotCalc = kFALSE;
 393    }
 394   return fGammaCMSLC;
 395 }
 396 /************************************************************************/
 397
 398 Double_t AliHBTPair::GetMt()
 399 {
 400   //Calculates transverse mass of the pair
 401   if (fMtNotCalc)
 402    {
 403      CalculateSums();
 404      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
 405      fMtNotCalc = kFALSE;
 406    }
 407    return fMt;
 408 }
 409 /************************************************************************/
 410
 411 Double_t AliHBTPair::GetWeight()
 412 {
 413   //returns and buffers weight for this pair
 414   if (fWeightNotCalc)
 415    {
 416       fWeight = AliHBTWeights::Weight(this);
 417       fWeightNotCalc = kFALSE;
 418    }
 419   return fWeight;
 420 }
 421 /************************************************************************/
 422
 423 Double_t AliHBTPair::GetAvarageDistance()
 424 {
 425 //returns and buffers avarage distance between two tracks calculated
 426 // out of track points (see AliHBTTrackPoints class)
 427
 428   if (fAvarageDistanceNotCalc)
 429    {
 430      fAvarageDistance = AvDistance();
 431      fAvarageDistanceNotCalc = kFALSE;
 432    }
 433   return fAvarageDistance;
 434 }
 435 /************************************************************************/
 436
 437 Double_t AliHBTPair::AvDistance()
 438 {
 439   //returns avarage distance between two tracks in range
 440   //as defined in Track-Points of AliHBTParticle
 441   //returns negative value if error uccured f.g. tracks do not have track-points
 442   AliHBTTrackPoints* tpts1 = fPart1->GetTrackPoints();
 443   if ( tpts1 == 0x0)
 444    {//it could be simulated pair
 445 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
 446      return -1.0;
 447    }
 448
 449   AliHBTTrackPoints* tpts2 = fPart2->GetTrackPoints();
 450   if ( tpts2 == 0x0)
 451    {
 452 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
 453      return -1.0;
 454    }
 455
 456   return tpts1->AvarageDistance(*tpts2);
 457 }