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