[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),
 fKStarOut(0.0),
 fKStarSide(0.0),
 fKStarLong(0.0),
 fKStarCompNotCalc(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),
 fKStarOut(0.0),
 fKStarSide(0.0),
 fKStarLong(0.0),
 fKStarCompNotCalc(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),
 fKStarOut(0.0),
 fKStarSide(0.0),
 fKStarLong(0.0),
 fKStarCompNotCalc(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::GetKStarOut()
{
    CalculateKStarComp();
    return fKStarOut;
}
/************************************************************************/
Double_t AliAODPair::GetKStarSide()
{
    CalculateKStarComp();
    return fKStarSide;
}
/************************************************************************/
Double_t AliAODPair::GetKStarLong()
{
    CalculateKStarComp();
    return fKStarLong;
}
/************************************************************************/

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