ANALYSIS/AliAODPair.cxx

   1 #include "AliAODPair.h"
   2 //_________________________________________________________________________
   3 ///////////////////////////////////////////////////////////////////////////
   4 //
   5 // class AliAODPair
   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 "AliVAODParticle.h"
  15 #include "AliTrackPoints.h"
  16
  17 ClassImp(AliAODPair)
  18
  19 /************************************************************************/
  20 AliAODPair::AliAODPair(Bool_t rev):
  21  fPart1(0x0),
  22  fPart2(0x0),
  23  fSwappedPair(0x0),
  24  fQSideLCMS(0.0),
  25  fQSideLCMSNotCalc(kTRUE),
  26  fQOutLCMS(0.0),
  27  fQOutLCMSNotCalc(kTRUE),
  28  fQLongLCMS(0.0),
  29  fQLongLCMSNotCalc(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  fAvarageDistance(0.0),
  49  fAvarageDistanceNotCalc(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  fGammaLCMS(0.0),
  61  fGammaLCMSNotCalc(kTRUE),
  62  fChanged(kTRUE)
  63  {
  64 //value of rev defines if it is Swapped
  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) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
  69
  70  }
  71 /************************************************************************/
  72
  73 AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
  74  fPart1(part1),
  75  fPart2(part2),
  76  fSwappedPair(0x0),
  77  fQSideLCMS(0.0),
  78  fQSideLCMSNotCalc(kTRUE),
  79  fQOutLCMS(0.0),
  80  fQOutLCMSNotCalc(kTRUE),
  81  fQLongLCMS(0.0),
  82  fQLongLCMSNotCalc(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  fAvarageDistance(0.0),
 102  fAvarageDistanceNotCalc(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  fGammaLCMS(0.0),
 114  fGammaLCMSNotCalc(kTRUE),
 115  fChanged(kTRUE)
 116  {
 117 //value of rev defines if it is Swapped
 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) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
 122
 123  }
 124 /************************************************************************/
 125 AliAODPair::AliAODPair(const AliAODPair& in):
 126  TObject(in),
 127  fPart1(0x0),
 128  fPart2(0x0),
 129  fSwappedPair(0x0),
 130  fQSideLCMS(0.0),
 131  fQSideLCMSNotCalc(kTRUE),
 132  fQOutLCMS(0.0),
 133  fQOutLCMSNotCalc(kTRUE),
 134  fQLongLCMS(0.0),
 135  fQLongLCMSNotCalc(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  fAvarageDistance(0.0),
 155  fAvarageDistanceNotCalc(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  fGammaLCMS(0.0),
 167  fGammaLCMSNotCalc(kTRUE),
 168  fChanged(kTRUE)
 169 {
 170  //cpy constructor
 171  in.Copy(*this);
 172 }
 173 /************************************************************************/
 174
 175 AliAODPair& AliAODPair::operator=(const AliAODPair& in)
 176 {
 177  //Assigment operator
 178  in.Copy(*this);
 179  return *this;
 180 }
 181 /************************************************************************/
 182
 183 Double_t AliAODPair::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
 199 Double_t AliAODPair::GetQSideLCMS()
 200 {
 201 //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 202
 203   if (fQSideLCMSNotCalc)
 204    {
 205     fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
 206     fQSideLCMSNotCalc = kFALSE;
 207    }
 208   return fQSideLCMS;
 209 }
 210 /************************************************************************/
 211
 212 Double_t AliAODPair::GetQOutLCMS()
 213 {
 214  //caculates Qout in Center Of Mass Longitudionally Co-Moving
 215  if(fQOutLCMSNotCalc)
 216   {
 217    CalculateSums();
 218    CalculateDiffs();
 219
 220    if (fPart1->Mass() != fPart2->Mass())
 221     {
 222 /*
 223       //STAR algorithm
 224       Double_t beta  = fPzSum/fESum;
 225       Double_t gamma = GetGammaToLCMS();
 226       Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
 227       Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
 228       beta  = 2.0*GetKt()/GetMt();
 229       gamma = GetMt()/GetQInv();
 230       fQOutLCMS = gamma * (x - beta * el);
 231 */
 232
 233       //beta=fPzSum/fESum;    // Longit. V == beta
 234       Double_t beta=fPzSum/fESum;
 235       Double_t gamma = GetGammaToLCMS();
 236
 237       Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
 238       Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
 239
 240 //      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
 241 //      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
 242       Double_t tmp;
 243       tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
 244       Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
 245       Double_t part1Px = tmp;
 246
 247       tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
 248       Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
 249       Double_t part2Px = tmp;
 250
 251
 252 //      LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
 253 //      LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
 254       Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
 255       Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());
 256
 257       Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
 258       Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
 259       Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
 260       Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
 261       Double_t sumE=part1E+part2E;
 262       Double_t sumPx=part1Px+part2Px;
 263       Double_t sumPy=part1Py+part2Py;
 264       Double_t sumPZ=part1Pz+part2Pz;
 265       Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;
 266
 267       Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
 268       Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
 269       fQOutLCMS=(part1Px-part2Px);//== id
 270       fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
 271     }
 272    else
 273     {
 274       Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
 275       fQOutLCMS = 0.5*k2/GetKt();
 276    // if (non-id)  fQOutLCMS=fQOutLCMS - sumPx*HF;
 277     }
 278
 279
 280    fQOutLCMSNotCalc = kFALSE;
 281   }
 282  return fQOutLCMS;
 283 }
 284 /************************************************************************/
 285
 286 Double_t AliAODPair::GetQLongLCMS()
 287 {
 288  //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
 289  if (fQLongLCMSNotCalc)
 290   {
 291     CalculateSums();
 292     CalculateDiffs();
 293     Double_t beta = fPzSum/fESum;
 294     fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
 295     fQLongLCMSNotCalc = kFALSE;
 296   }
 297  return fQLongLCMS;
 298 }
 299 /************************************************************************/
 300
 301 Double_t AliAODPair::GetKt()
 302 {
 303  //calculates the evarage momentum of the pair
 304   if(fKtNotCalc)
 305    {
 306      CalculateSums();
 307      fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
 308      fKtNotCalc = kFALSE;
 309    }
 310   return fKt;
 311 }
 312 /************************************************************************/
 313
 314 Double_t AliAODPair::GetKStar()
 315 {
 316   //calculates invariant velocity difference
 317   if (fKStarNotCalc)
 318    {
 319     CalculateSums();
 320
 321     Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
 322     Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
 323     if (ptrans > mtrans)
 324      {
 325        Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
 326        Error("GetKStar","Particle1:");
 327        fPart1->Print();
 328        Error("GetKStar","Particle2:");
 329        fPart2->Print();
 330        Error("GetKStar","");
 331
 332        fKStar = 10e5;
 333        fKStarNotCalc = kFALSE;
 334        return fKStar;
 335      }
 336     Double_t pinv =   TMath::Sqrt(mtrans - ptrans);
 337
 338     Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
 339
 340     CalculateQInvL();
 341
 342     q = q*q - fQInvL;
 343     if ( q < 0)
 344      {
 345        Info("GetKStar","Sqrt of negative number q = %f",q);
 346        Error("GetKStar","Particle1:");
 347        fPart1->Print();
 348        Error("GetKStar","Particle2:");
 349        fPart2->Print();
 350        fKStar = 10e5;
 351        fKStarNotCalc = kFALSE;
 352        return fKStar;
 353      }
 354
 355     q = TMath::Sqrt(q);
 356     fKStar = q/2.;
 357     fKStarNotCalc = kFALSE;
 358    }
 359   return fKStar;
 360 }
 361 /************************************************************************/
 362
 363 Double_t AliAODPair::GetQInv()
 364 {
 365 //returns Qinv
 366 //warning for non-id particles you want to use 2*KStar
 367   if(fQInvNotCalc)
 368    {
 369     CalculateQInvL();
 370     fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
 371     fQInvNotCalc = kFALSE;
 372    }
 373   return fQInv;
 374 }
 375 /************************************************************************/
 376
 377 Double_t AliAODPair::GetGammaToLCMS()
 378 {
 379   //calculates gamma factor of the boost to LCMS
 380   if(fGammaLCMSNotCalc)
 381    {
 382      CalculateSums();
 383      Double_t beta = fPzSum/fESum;
 384      fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
 385      fGammaLCMSNotCalc = kFALSE;
 386    }
 387   return fGammaLCMS;
 388 }
 389 /************************************************************************/
 390
 391 Double_t AliAODPair::GetMt()
 392 {
 393   //Calculates transverse mass of the pair
 394   if (fMtNotCalc)
 395    {
 396      CalculateSums();
 397      fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
 398      fMtNotCalc = kFALSE;
 399    }
 400    return fMt;
 401 }
 402 /************************************************************************/
 403
 404 Double_t AliAODPair::GetAvarageDistance()
 405 {
 406 //returns and buffers avarage distance between two tracks calculated
 407 // out of track points (see AliAODTrackPoints class)
 408
 409   if (fAvarageDistanceNotCalc)
 410    {
 411      fAvarageDistance = AvDistance();
 412      fAvarageDistanceNotCalc = kFALSE;
 413    }
 414   return fAvarageDistance;
 415 }
 416 /************************************************************************/
 417
 418 Double_t AliAODPair::AvDistance()
 419 {
 420   //returns avarage distance between two tracks in range
 421   //as defined in Track-Points of AliVAODParticle
 422   //returns negative value if error uccured f.g. tracks do not have track-points
 423   AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
 424   if ( tpts1 == 0x0)
 425    {//it could be simulated pair
 426 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
 427      return -1.0;
 428    }
 429
 430   AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
 431   if ( tpts2 == 0x0)
 432    {
 433 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
 434      return -1.0;
 435    }
 436
 437   return tpts1->AvarageDistance(*tpts2);
 438 }