[u/mrichter/AliRoot.git] / ANALYSIS / AliAODPair.cxx

#include "AliAODPair.h"
//_________________________________________________________________________
///////////////////////////////////////////////////////////////////////////
//
// class AliAODPair
//
// class implements pair of particles and taking care of caluclation (almost)
// all of pair properties (Qinv, InvMass,...)
// 
// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
//
////////////////////////////////////////////////////////////////////////////

#include "AliVAODParticle.h"
#include "AliTrackPoints.h"

ClassImp(AliAODPair)

/************************************************************************/
AliAODPair::AliAODPair(Bool_t rev):
 fPart1(0x0),
 fPart2(0x0),
 fSwappedPair(0x0),
 fQSideLCMS(0.0),
 fQSideLCMSNotCalc(kTRUE),
 fQOutLCMS(0.0),
 fQOutLCMSNotCalc(kTRUE),
 fQLongLCMS(0.0),
 fQLongLCMSNotCalc(kTRUE),
 fQtLCMS(0.0),
 fQtLCMSNotCalc(kTRUE),
 fQt(0.0),
 fQtNotCalc(kTRUE),
 fQInv(0.0),
 fQInvNotCalc(kTRUE),
 fInvMass(0.0),
 fInvMassNotCalc(kTRUE),
 fKt(0.0),
 fKtNotCalc(kTRUE),
 fKStar(0.0),
 fKStarNotCalc(kTRUE),
 fPInv(0.0),
 fQSide(0.0),
 fOut(0.0),
 fQLong(0.0),
 fMt(0.0),
 fMtNotCalc(kTRUE),
 fInvMassSqr(0.0),
 fMassSqrNotCalc(kTRUE),
 fQInvL(0.0),
 fQInvLNotCalc(kTRUE),
 fAvarageDistance(0.0),
 fAvarageDistanceNotCalc(kTRUE),
 fPxSum(0.0),
 fPySum(0.0),
 fPzSum(0.0),
 fESum(0.0),
 fSumsNotCalc(kTRUE),
 fPxDiff(0.0),
 fPyDiff(0.0),
 fPzDiff(0.0),
 fEDiff(0.0),
 fDiffsNotCalc(kTRUE),
 fGammaLCMS(0.0),
 fGammaLCMSNotCalc(kTRUE),
 fChanged(kTRUE)
 {
//value of rev defines if it is Swapped
//if you pass kTRUE swpaped pair will NOT be created
//though you wont be able to get the swaped pair from this pair

  if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
  
 }
/************************************************************************/

AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
 fPart1(part1),
 fPart2(part2),
 fSwappedPair(0x0),
 fQSideLCMS(0.0),
 fQSideLCMSNotCalc(kTRUE),
 fQOutLCMS(0.0),
 fQOutLCMSNotCalc(kTRUE),
 fQLongLCMS(0.0),
 fQLongLCMSNotCalc(kTRUE),
 fQtLCMS(0.0),
 fQtLCMSNotCalc(kTRUE),
 fQInv(0.0),
 fQInvNotCalc(kTRUE),
 fInvMass(0.0),
 fInvMassNotCalc(kTRUE),
 fKt(0.0),
 fKtNotCalc(kTRUE),
 fKStar(0.0),
 fKStarNotCalc(kTRUE),
 fPInv(0.0),
 fQSide(0.0),
 fOut(0.0),
 fQLong(0.0),
 fMt(0.0),
 fMtNotCalc(kTRUE),
 fInvMassSqr(0.0),
 fMassSqrNotCalc(kTRUE),
 fQInvL(0.0),
 fQInvLNotCalc(kTRUE),
 fAvarageDistance(0.0),
 fAvarageDistanceNotCalc(kTRUE),
 fPxSum(0.0),
 fPySum(0.0),
 fPzSum(0.0),
 fESum(0.0),
 fSumsNotCalc(kTRUE),
 fPxDiff(0.0),
 fPyDiff(0.0),
 fPzDiff(0.0),
 fEDiff(0.0),
 fDiffsNotCalc(kTRUE),
 fGammaLCMS(0.0),
 fGammaLCMSNotCalc(kTRUE),
 fChanged(kTRUE)
 {
//value of rev defines if it is Swapped
//if you pass kTRUE swpaped pair will NOT be created
//though you wont be able to get the swaped pair from this pair

  if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
  
 }
/************************************************************************/
AliAODPair::AliAODPair(const AliAODPair& in):
 TObject(in),
 fPart1(0x0),
 fPart2(0x0),
 fSwappedPair(0x0),
 fQSideLCMS(0.0),
 fQSideLCMSNotCalc(kTRUE),
 fQOutLCMS(0.0),
 fQOutLCMSNotCalc(kTRUE),
 fQLongLCMS(0.0),
 fQLongLCMSNotCalc(kTRUE),
 fQtLCMS(0.0),
 fQtLCMSNotCalc(kTRUE),
 fQInv(0.0),
 fQInvNotCalc(kTRUE),
 fInvMass(0.0),
 fInvMassNotCalc(kTRUE),
 fKt(0.0),
 fKtNotCalc(kTRUE),
 fKStar(0.0),
 fKStarNotCalc(kTRUE),
 fPInv(0.0),
 fQSide(0.0),
 fOut(0.0),
 fQLong(0.0),
 fMt(0.0),
 fMtNotCalc(kTRUE),
 fInvMassSqr(0.0),
 fMassSqrNotCalc(kTRUE),
 fQInvL(0.0),
 fQInvLNotCalc(kTRUE),
 fAvarageDistance(0.0),
 fAvarageDistanceNotCalc(kTRUE),
 fPxSum(0.0),
 fPySum(0.0),
 fPzSum(0.0),
 fESum(0.0),
 fSumsNotCalc(kTRUE),
 fPxDiff(0.0),
 fPyDiff(0.0),
 fPzDiff(0.0),
 fEDiff(0.0),
 fDiffsNotCalc(kTRUE),
 fGammaLCMS(0.0),
 fGammaLCMSNotCalc(kTRUE),
 fChanged(kTRUE)
{
 //cpy constructor
 in.Copy(*this);
}
/************************************************************************/

AliAODPair& AliAODPair::operator=(const AliAODPair& in)
{
 //Assigment operator
 in.Copy(*this);
 return *this;
}
/************************************************************************/

Double_t AliAODPair::GetInvMass()
{
//Returns qinv value for a pair
  if(fInvMassNotCalc)
   {
     CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method 
     
     if(fInvMassSqr<0)  fInvMass = TMath::Sqrt(-fInvMassSqr);
     else fInvMass = TMath::Sqrt(fInvMassSqr); 
     
     fInvMassNotCalc = kFALSE;
   }
  return fInvMass;
}
/************************************************************************/

Double_t AliAODPair::GetQSideLCMS()
{
//return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 
  if (fQSideLCMSNotCalc)
   {
    fQSideLCMS = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
    fQSideLCMSNotCalc = kFALSE;
   }
  return fQSideLCMS;
}
/************************************************************************/

Double_t AliAODPair::GetQOutLCMS()
{
 //caculates Qout in Center Of Mass Longitudionally Co-Moving
 if(fQOutLCMSNotCalc)
  {
   CalculateSums();
   CalculateDiffs();

   if (fPart1->Mass() != fPart2->Mass())
    {
/*    
      //STAR algorithm
      Double_t beta  = fPzSum/fESum;
      Double_t gamma = GetGammaToLCMS();
      Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
      Double_t x  = ( fPart1->Px()*fPxSum + fPart1->Py()*fPySum) / ( 2.0*GetKt() );
      beta  = 2.0*GetKt()/GetMt();
      gamma = GetMt()/GetQInv();
      fQOutLCMS = gamma * (x - beta * el);
*/

      //beta=fPzSum/fESum;    // Longit. V == beta
      Double_t beta=fPzSum/fESum;
      Double_t gamma = GetGammaToLCMS();
      
      Double_t cosphi=fPxSum/(2.0*GetKt());  // cos(phi)
      Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
      
//      ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
//      ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
      Double_t tmp;
      tmp = fPart1->Px()*cosphi + fPart1->Py()*sinphi;
      Double_t part1Py = fPart1->Py()*cosphi - fPart1->Px()*sinphi;
      Double_t part1Px = tmp;

      tmp = fPart2->Px()*cosphi + fPart2->Py()*sinphi;
      Double_t part2Py = fPart2->Py()*cosphi - fPart2->Px()*sinphi;
      Double_t part2Px = tmp;
      
      
//      LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
//      LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
      Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->E());
      Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->E());

      Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
      Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
      Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
      Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
      Double_t sumE=part1E+part2E;
      Double_t sumPx=part1Px+part2Px;
      Double_t sumPy=part1Py+part2Py;
      Double_t sumPZ=part1Pz+part2Pz;
      Double_t sumP2=sumPx*sumPx+sumPy*sumPy+sumPZ*sumPZ;

      Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
      Double_t hf = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/(relmass*relmass);
      fQOutLCMS=(part1Px-part2Px);//== id
      fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
    }
   else
    {
      Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
      fQOutLCMS = 0.5*k2/GetKt();
   // if (non-id)  fQOutLCMS=fQOutLCMS - sumPx*HF;
    }

    
   fQOutLCMSNotCalc = kFALSE;
  }
 return fQOutLCMS;
}
/************************************************************************/

Double_t AliAODPair::GetQLongLCMS()
{
 //return Q Long in Central Of Mass System in Longitudialy Comoving Frame
 if (fQLongLCMSNotCalc)
  {
    CalculateSums();
    CalculateDiffs();
    Double_t beta = fPzSum/fESum;
    fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
    fQLongLCMSNotCalc = kFALSE;
  }
 return fQLongLCMS; 
}
/************************************************************************/

Double_t AliAODPair::GetQtLCMS()
{
 //returns Q transverse CMS longitudionally co-moving
 if (fQtLCMSNotCalc)
  {
    fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
    fQtLCMSNotCalc = kFALSE;
  }
 return fQtLCMS; 
}
/************************************************************************/

Double_t AliAODPair::GetQt()
{
 //returns Q transverse CMS longitudionally co-moving
 if (fQtNotCalc)
  {
    CalculateSums();
    CalculateDiffs();
    
    Double_t dotprod = fPxSum*fPxDiff + fPySum*fPyDiff + fPzSum*fPzDiff;
    Double_t klen =    fPxSum*fPxSum  + fPySum*fPySum  + fPzSum*fPzSum;
    klen = TMath::Sqrt(klen);
    Double_t qlen =    fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff;
    qlen = TMath::Sqrt(qlen);
    Double_t tmp = klen*qlen;
    if (tmp == 0.0)
     {
       fQt = 10e5;
       fQtNotCalc = kFALSE;
       return fQt;
     }
    Double_t cosopenangle = dotprod/tmp;
    Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
    
    fQt = sinopenangle*qlen;
    fQtNotCalc = kFALSE;
  }
 return fQt; 
}
/************************************************************************/

Double_t AliAODPair::GetKt()
{
 //calculates the evarage momentum of the pair
  if(fKtNotCalc)
   { 
     CalculateSums();
     fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
     fKtNotCalc = kFALSE;
   }
  return fKt;
}
/************************************************************************/

Double_t AliAODPair::GetKStar()
{
  //calculates invariant velocity difference
  if (fKStarNotCalc)
   { 
    CalculateSums();

    Double_t ptrans = fPxSum*fPxSum + fPySum*fPySum;
    Double_t mtrans = fESum*fESum - fPzSum*fPzSum;
    if (ptrans > mtrans)
     {
       Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
       Error("GetKStar","Particle1:");
       fPart1->Print();
       Error("GetKStar","Particle2:");
       fPart2->Print();
       Error("GetKStar","");
       
       fKStar = 10e5;
       fKStarNotCalc = kFALSE;
       return fKStar;
     }
    Double_t pinv =   TMath::Sqrt(mtrans - ptrans);

    Double_t q = (fPart1->Mass()*fPart1->Mass() - fPart2->Mass()*fPart2->Mass())/pinv;
    
    CalculateQInvL();
    
    q = q*q - fQInvL;
    if ( q < 0)
     {
       Info("GetKStar","Sqrt of negative number q = %f",q);
       Error("GetKStar","Particle1:");
       fPart1->Print();
       Error("GetKStar","Particle2:");
       fPart2->Print();
       fKStar = 10e5;
       fKStarNotCalc = kFALSE;
       return fKStar;
     }
     
    q = TMath::Sqrt(q);
    fKStar = q/2.;
    fKStarNotCalc = kFALSE;
   }
  return fKStar;
}
/************************************************************************/

Double_t AliAODPair::GetQInv()
{
//returns Qinv 
//warning for non-id particles you want to use 2*KStar
  if(fQInvNotCalc)
   {
    CalculateQInvL();
    fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
    fQInvNotCalc = kFALSE;
   }
  return fQInv;
}
/************************************************************************/

Double_t AliAODPair::GetGammaToLCMS()
{
  //calculates gamma factor of the boost to LCMS
  if(fGammaLCMSNotCalc)
   {
     CalculateSums();
     Double_t beta = fPzSum/fESum;
     fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
     fGammaLCMSNotCalc = kFALSE;
   }
  return fGammaLCMS;
}
/************************************************************************/

Double_t AliAODPair::GetMt()
{
  //Calculates transverse mass of the pair
  if (fMtNotCalc)
   {
     CalculateSums();
     fMt = TMath::Sqrt(fESum*fESum - fPzSum*fPzSum);
     fMtNotCalc = kFALSE;
   }
   return fMt;
}
/************************************************************************/

Double_t AliAODPair::GetAvarageDistance()
{
//returns and buffers avarage distance between two tracks calculated 
// out of track points (see AliAODTrackPoints class)

  if (fAvarageDistanceNotCalc)
   {
     fAvarageDistance = AvDistance();
     fAvarageDistanceNotCalc = kFALSE;
   }
  return fAvarageDistance;
}
/************************************************************************/

Double_t AliAODPair::AvDistance()
{
  //returns avarage distance between two tracks in range 
  //as defined in Track-Points of AliVAODParticle
  //returns negative value if error uccured f.g. tracks do not have track-points
  AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
  if ( tpts1 == 0x0)
   {//it could be simulated pair
//     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
     return -1.0;
   }

  AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
  if ( tpts2 == 0x0)
   {
//     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
     return -1.0;
   }

  return tpts1->AvarageDistance(*tpts2);
}
Commit	Line	Data
78d7c6d3	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),
f422a7da	30	fQtLCMS(0.0),
f422a7da	31	fQtLCMSNotCalc(kTRUE),
fb629e2d	32	fQt(0.0),
fb629e2d	33	fQtNotCalc(kTRUE),
78d7c6d3	34	fQInv(0.0),
	35	fQInvNotCalc(kTRUE),
	36	fInvMass(0.0),
	37	fInvMassNotCalc(kTRUE),
	38	fKt(0.0),
	39	fKtNotCalc(kTRUE),
	40	fKStar(0.0),
	41	fKStarNotCalc(kTRUE),
	42	fPInv(0.0),
	43	fQSide(0.0),
	44	fOut(0.0),
	45	fQLong(0.0),
	46	fMt(0.0),
	47	fMtNotCalc(kTRUE),
	48	fInvMassSqr(0.0),
	49	fMassSqrNotCalc(kTRUE),
	50	fQInvL(0.0),
	51	fQInvLNotCalc(kTRUE),
	52	fAvarageDistance(0.0),
	53	fAvarageDistanceNotCalc(kTRUE),
	54	fPxSum(0.0),
	55	fPySum(0.0),
	56	fPzSum(0.0),
	57	fESum(0.0),
	58	fSumsNotCalc(kTRUE),
	59	fPxDiff(0.0),
	60	fPyDiff(0.0),
	61	fPzDiff(0.0),
	62	fEDiff(0.0),
	63	fDiffsNotCalc(kTRUE),
	64	fGammaLCMS(0.0),
	65	fGammaLCMSNotCalc(kTRUE),
	66	fChanged(kTRUE)
	67	{
	68	//value of rev defines if it is Swapped
	69	//if you pass kTRUE swpaped pair will NOT be created
	70	//though you wont be able to get the swaped pair from this pair
	71
	72	if(!rev) fSwappedPair = new AliAODPair(kTRUE); //if false create swaped pair
	73
	74	}
	75	/************************************************************************/
	76
	77	AliAODPair::AliAODPair(AliVAODParticle* part1, AliVAODParticle* part2, Bool_t rev):
	78	fPart1(part1),
	79	fPart2(part2),
	80	fSwappedPair(0x0),
	81	fQSideLCMS(0.0),
	82	fQSideLCMSNotCalc(kTRUE),
	83	fQOutLCMS(0.0),
	84	fQOutLCMSNotCalc(kTRUE),
	85	fQLongLCMS(0.0),
	86	fQLongLCMSNotCalc(kTRUE),
f422a7da	87	fQtLCMS(0.0),
f422a7da	88	fQtLCMSNotCalc(kTRUE),
78d7c6d3	89	fQInv(0.0),
	90	fQInvNotCalc(kTRUE),
	91	fInvMass(0.0),
	92	fInvMassNotCalc(kTRUE),
	93	fKt(0.0),
	94	fKtNotCalc(kTRUE),
	95	fKStar(0.0),
	96	fKStarNotCalc(kTRUE),
	97	fPInv(0.0),
	98	fQSide(0.0),
	99	fOut(0.0),
	100	fQLong(0.0),
	101	fMt(0.0),
	102	fMtNotCalc(kTRUE),
	103	fInvMassSqr(0.0),
	104	fMassSqrNotCalc(kTRUE),
	105	fQInvL(0.0),
	106	fQInvLNotCalc(kTRUE),
	107	fAvarageDistance(0.0),
	108	fAvarageDistanceNotCalc(kTRUE),
	109	fPxSum(0.0),
	110	fPySum(0.0),
	111	fPzSum(0.0),
	112	fESum(0.0),
	113	fSumsNotCalc(kTRUE),
	114	fPxDiff(0.0),
	115	fPyDiff(0.0),
	116	fPzDiff(0.0),
	117	fEDiff(0.0),
	118	fDiffsNotCalc(kTRUE),
	119	fGammaLCMS(0.0),
	120	fGammaLCMSNotCalc(kTRUE),
	121	fChanged(kTRUE)
	122	{
	123	//value of rev defines if it is Swapped
	124	//if you pass kTRUE swpaped pair will NOT be created
	125	//though you wont be able to get the swaped pair from this pair
	126
	127	if(!rev) fSwappedPair = new AliAODPair(part2,part1,kTRUE); //if false create swaped pair
	128
	129	}
	130	/************************************************************************/
	131	AliAODPair::AliAODPair(const AliAODPair& in):
	132	TObject(in),
	133	fPart1(0x0),
	134	fPart2(0x0),
	135	fSwappedPair(0x0),
	136	fQSideLCMS(0.0),
	137	fQSideLCMSNotCalc(kTRUE),
	138	fQOutLCMS(0.0),
	139	fQOutLCMSNotCalc(kTRUE),
	140	fQLongLCMS(0.0),
	141	fQLongLCMSNotCalc(kTRUE),
f422a7da	142	fQtLCMS(0.0),
f422a7da	143	fQtLCMSNotCalc(kTRUE),
78d7c6d3	144	fQInv(0.0),
	145	fQInvNotCalc(kTRUE),
	146	fInvMass(0.0),
	147	fInvMassNotCalc(kTRUE),
	148	fKt(0.0),
	149	fKtNotCalc(kTRUE),
	150	fKStar(0.0),
	151	fKStarNotCalc(kTRUE),
	152	fPInv(0.0),
	153	fQSide(0.0),
	154	fOut(0.0),
	155	fQLong(0.0),
	156	fMt(0.0),
	157	fMtNotCalc(kTRUE),
	158	fInvMassSqr(0.0),
	159	fMassSqrNotCalc(kTRUE),
	160	fQInvL(0.0),
	161	fQInvLNotCalc(kTRUE),
	162	fAvarageDistance(0.0),
	163	fAvarageDistanceNotCalc(kTRUE),
	164	fPxSum(0.0),
	165	fPySum(0.0),
	166	fPzSum(0.0),
	167	fESum(0.0),
	168	fSumsNotCalc(kTRUE),
	169	fPxDiff(0.0),
	170	fPyDiff(0.0),
	171	fPzDiff(0.0),
	172	fEDiff(0.0),
	173	fDiffsNotCalc(kTRUE),
	174	fGammaLCMS(0.0),
	175	fGammaLCMSNotCalc(kTRUE),
	176	fChanged(kTRUE)
	177	{
	178	//cpy constructor
	179	in.Copy(*this);
	180	}
	181	/************************************************************************/
	182
	183	AliAODPair& AliAODPair::operator=(const AliAODPair& in)
	184	{
	185	//Assigment operator
	186	in.Copy(*this);
	187	return *this;
	188	}
	189	/************************************************************************/
	190
	191	Double_t AliAODPair::GetInvMass()
	192	{
	193	//Returns qinv value for a pair
	194	if(fInvMassNotCalc)
	195	{
	196	CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
	197
	198	if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
	199	else fInvMass = TMath::Sqrt(fInvMassSqr);
	200
	201	fInvMassNotCalc = kFALSE;
	202	}
	203	return fInvMass;
	204	}
	205	/************************************************************************/
	206
	207	Double_t AliAODPair::GetQSideLCMS()
208	{
209	//return Q Side in Central Of Mass System in Longitudialy Comoving Frame
210
211	if (fQSideLCMSNotCalc)
212	{
213	fQSideLCMS = (fPart1->Px()fPart2->Py()-fPart2->Px()fPart1->Py())/GetKt();
214	fQSideLCMSNotCalc = kFALSE;
215	}
216	return fQSideLCMS;
217	}
218	/************************************************************************/
219
220	Double_t AliAODPair::GetQOutLCMS()
221	{
222	//caculates Qout in Center Of Mass Longitudionally Co-Moving
223	if(fQOutLCMSNotCalc)
224	{
225	CalculateSums();
226	CalculateDiffs();
227
228	if (fPart1->Mass() != fPart2->Mass())
229	{
230	/*
231	//STAR algorithm
232	Double_t beta = fPzSum/fESum;
233	Double_t gamma = GetGammaToLCMS();
234	Double_t el = gamma * (fPart1->E() - beta * fPart1->Pz());
235	Double_t x = ( fPart1->Px()fPxSum + fPart1->Py()fPySum) / ( 2.0*GetKt() );
236	beta = 2.0*GetKt()/GetMt();
237	gamma = GetMt()/GetQInv();
238	fQOutLCMS = gamma * (x - beta * el);
239	*/
240
241	//beta=fPzSum/fESum; // Longit. V == beta
242	Double_t beta=fPzSum/fESum;
243	Double_t gamma = GetGammaToLCMS();
244
245	Double_t cosphi=fPxSum/(2.0*GetKt()); // cos(phi)
246	Double_t sinphi=fPySum/(2.0*GetKt()); // sin(phi)
247
248	// ROTATE(part1Px,part1Py,SPHI,CPHI,part1Px,part1Py);//ROT8
249	// ROTATE(part2Px,part2Py,SPHI,CPHI,part2Px,part2Py);//ROT8
250	Double_t tmp;
251	tmp = fPart1->Px()cosphi + fPart1->Py()sinphi;
252	Double_t part1Py = fPart1->Py()cosphi - fPart1->Px()sinphi;
253	Double_t part1Px = tmp;
254
255	tmp = fPart2->Px()cosphi + fPart2->Py()sinphi;
256	Double_t part2Py = fPart2->Py()cosphi - fPart2->Px()sinphi;
257	Double_t part2Px = tmp;
258
259
260	// LTR(part1Pz,E1,beta,GetGammaToLCMS(),part1Pz,E1a);
261	// LTR(part2Pz,E2,beta,GetGammaToLCMS(),part2Pz,E2a);
262	Double_t part1Pz=gamma(fPart1->Pz()-betafPart1->E());
263	Double_t part2Pz=gamma(fPart2->Pz()-betafPart2->E());
264
265	Double_t part1P2=part1Pxpart1Px+part1Pypart1Py+part1Pz*part1Pz;
266	Double_t part2P2=part2Pxpart2Px+part2Pypart2Py+part2Pz*part2Pz;
267	Double_t part1E=TMath::Sqrt(fPart1->Mass()*fPart1->Mass()+part1P2);
268	Double_t part2E=TMath::Sqrt(fPart2->Mass()*fPart2->Mass()+part2P2);
269	Double_t sumE=part1E+part2E;
270	Double_t sumPx=part1Px+part2Px;
271	Double_t sumPy=part1Py+part2Py;
272	Double_t sumPZ=part1Pz+part2Pz;
273	Double_t sumP2=sumPxsumPx+sumPysumPy+sumPZ*sumPZ;
274
275	Double_t relmass=TMath::Sqrt(sumE*sumE-sumP2);
276	Double_t hf = (fPart1->Mass()fPart1->Mass() - fPart2->Mass()fPart2->Mass())/(relmass*relmass);
277	fQOutLCMS=(part1Px-part2Px);//== id
278	fQOutLCMS=fQOutLCMS-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
279	}
280	else
281	{
282	Double_t k2 = fPxSumfPxDiff+fPySumfPyDiff;
283	fQOutLCMS = 0.5*k2/GetKt();
284	// if (non-id) fQOutLCMS=fQOutLCMS - sumPx*HF;
285	}
286
287
288	fQOutLCMSNotCalc = kFALSE;
289	}
290	return fQOutLCMS;
291	}
292	/************************************************************************/
293
294	Double_t AliAODPair::GetQLongLCMS()
295	{
296	//return Q Long in Central Of Mass System in Longitudialy Comoving Frame
297	if (fQLongLCMSNotCalc)
298	{
299	CalculateSums();
300	CalculateDiffs();
301	Double_t beta = fPzSum/fESum;
302	fQLongLCMS = GetGammaToLCMS() * ( fPzDiff - beta*fEDiff );
303	fQLongLCMSNotCalc = kFALSE;
304	}
305	return fQLongLCMS;
306	}
307	/************************************************************************/
308
f422a7da	309	Double_t AliAODPair::GetQtLCMS()
	310	{
	311	//returns Q transverse CMS longitudionally co-moving
	312	if (fQtLCMSNotCalc)
	313	{
	314	fQtLCMS = TMath::Hypot(GetQOutLCMS(),GetQSideLCMS());
	315	fQtLCMSNotCalc = kFALSE;
	316	}
	317	return fQtLCMS;
	318	}
fb629e2d	319	/************************************************************************/
	320
	321	Double_t AliAODPair::GetQt()
	322	{
	323	//returns Q transverse CMS longitudionally co-moving
	324	if (fQtNotCalc)
	325	{
05eb9603	326	CalculateSums();
	327	CalculateDiffs();
	328
fb629e2d	329	Double_t dotprod = fPxSumfPxDiff + fPySumfPyDiff + fPzSum*fPzDiff;
	330	Double_t klen = fPxSumfPxSum + fPySumfPySum + fPzSum*fPzSum;
	331	klen = TMath::Sqrt(klen);
	332	Double_t qlen = fPxDifffPxDiff + fPyDifffPyDiff + fPzDiff*fPzDiff;
	333	qlen = TMath::Sqrt(qlen);
05eb9603	334	Double_t tmp = klen*qlen;
	335	if (tmp == 0.0)
	336	{
	337	fQt = 10e5;
	338	fQtNotCalc = kFALSE;
	339	return fQt;
	340	}
	341	Double_t cosopenangle = dotprod/tmp;
fb629e2d	342	Double_t sinopenangle = TMath::Sqrt(1.0 - cosopenangle*cosopenangle);
	343
	344	fQt = sinopenangle*qlen;
	345	fQtNotCalc = kFALSE;
	346	}
	347	return fQt;
	348	}
	349	/************************************************************************/
f422a7da	350
78d7c6d3	351	Double_t AliAODPair::GetKt()
	352	{
	353	//calculates the evarage momentum of the pair
	354	if(fKtNotCalc)
	355	{
	356	CalculateSums();
	357	fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
	358	fKtNotCalc = kFALSE;
	359	}
	360	return fKt;
	361	}
	362	/************************************************************************/
	363
	364	Double_t AliAODPair::GetKStar()
	365	{
	366	//calculates invariant velocity difference
	367	if (fKStarNotCalc)
	368	{
	369	CalculateSums();
	370
	371	Double_t ptrans = fPxSumfPxSum + fPySumfPySum;
	372	Double_t mtrans = fESumfESum - fPzSumfPzSum;
	373	if (ptrans > mtrans)
	374	{
	375	Error("GetKStar","Tranverse momentum bigger than transverse mass. Not normal for on-shell particles");
	376	Error("GetKStar","Particle1:");
	377	fPart1->Print();
	378	Error("GetKStar","Particle2:");
	379	fPart2->Print();
	380	Error("GetKStar","");
	381
	382	fKStar = 10e5;
	383	fKStarNotCalc = kFALSE;
	384	return fKStar;
	385	}
	386	Double_t pinv = TMath::Sqrt(mtrans - ptrans);
	387
	388	Double_t q = (fPart1->Mass()fPart1->Mass() - fPart2->Mass()fPart2->Mass())/pinv;
	389
	390	CalculateQInvL();
	391
	392	q = q*q - fQInvL;
	393	if ( q < 0)
	394	{
	395	Info("GetKStar","Sqrt of negative number q = %f",q);
	396	Error("GetKStar","Particle1:");
	397	fPart1->Print();
	398	Error("GetKStar","Particle2:");
	399	fPart2->Print();
	400	fKStar = 10e5;
	401	fKStarNotCalc = kFALSE;
	402	return fKStar;
	403	}
	404
	405	q = TMath::Sqrt(q);
	406	fKStar = q/2.;
	407	fKStarNotCalc = kFALSE;
	408	}
	409	return fKStar;
	410	}
	411	/************************************************************************/
	412
	413	Double_t AliAODPair::GetQInv()
	414	{
415	//returns Qinv
416	//warning for non-id particles you want to use 2*KStar
417	if(fQInvNotCalc)
418	{
419	CalculateQInvL();
420	fQInv = TMath::Sqrt(TMath::Abs(fQInvL));
421	fQInvNotCalc = kFALSE;
422	}
423	return fQInv;
424	}
425	/************************************************************************/
426
427	Double_t AliAODPair::GetGammaToLCMS()
428	{
429	//calculates gamma factor of the boost to LCMS
430	if(fGammaLCMSNotCalc)
431	{
432	CalculateSums();
433	Double_t beta = fPzSum/fESum;
434	fGammaLCMS = 1.0/TMath::Sqrt(1.0 - beta*beta);
435	fGammaLCMSNotCalc = kFALSE;
436	}
437	return fGammaLCMS;
438	}
439	/************************************************************************/
440
441	Double_t AliAODPair::GetMt()
442	{
443	//Calculates transverse mass of the pair
444	if (fMtNotCalc)
445	{
446	CalculateSums();
447	fMt = TMath::Sqrt(fESumfESum - fPzSumfPzSum);
448	fMtNotCalc = kFALSE;
449	}
450	return fMt;
451	}
452	/************************************************************************/
453
454	Double_t AliAODPair::GetAvarageDistance()
455	{
456	//returns and buffers avarage distance between two tracks calculated
457	// out of track points (see AliAODTrackPoints class)
458
459	if (fAvarageDistanceNotCalc)
460	{
461	fAvarageDistance = AvDistance();
462	fAvarageDistanceNotCalc = kFALSE;
463	}
464	return fAvarageDistance;
465	}
466	/************************************************************************/
467
468	Double_t AliAODPair::AvDistance()
469	{
470	//returns avarage distance between two tracks in range
471	//as defined in Track-Points of AliVAODParticle
472	//returns negative value if error uccured f.g. tracks do not have track-points
473	AliTrackPoints* tpts1 = fPart1->GetTPCTrackPoints();
474	if ( tpts1 == 0x0)
475	{//it could be simulated pair
476	// Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
477	return -1.0;
478	}
479
480	AliTrackPoints* tpts2 = fPart2->GetTPCTrackPoints();
481	if ( tpts2 == 0x0)
482	{
483	// Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
484	return -1.0;
485	}
486
487	return tpts1->AvarageDistance(*tpts2);
488	}