[u/mrichter/AliRoot.git] / HBTAN / AliHBTPair.cxx

#include "AliHBTPair.h"
#include "AliHBTParticle.h"

ClassImp(AliHBTPair)

/************************************************************************/
AliHBTPair::AliHBTPair(Bool_t rev):
 fPart1(0x0),
 fPart2(0x0),
 fSwapedPair(0x0),
 fQSideCMSLC(0.0),
 fQSideCMSLCNotCalc(kTRUE),
 fQOutCMSLC(0.0),
 fQOutCMSLCNotCalc(kTRUE),
 fQLongCMSLC(0.0),
 fQLongCMSLCNotCalc(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),
 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),
 fGammaCMSLC(0.0),
 fGammaCMSLCNotCalc(kTRUE),
 fChanged(kTRUE)
 {
//value of rev defines if it is Swaped
//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) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
  
 }
/************************************************************************/

AliHBTPair::AliHBTPair(AliHBTParticle* part1, AliHBTParticle* part2, Bool_t rev):
 fPart1(part1),
 fPart2(part2),
 fSwapedPair(0x0),
 fQSideCMSLC(0.0),
 fQSideCMSLCNotCalc(kTRUE),
 fQOutCMSLC(0.0),
 fQOutCMSLCNotCalc(kTRUE),
 fQLongCMSLC(0.0),
 fQLongCMSLCNotCalc(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),
 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),
 fGammaCMSLC(0.0),
 fGammaCMSLCNotCalc(kTRUE),
 fChanged(kTRUE)
 {
//value of rev defines if it is Swaped
//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) fSwapedPair = new AliHBTPair(part2,part1,kTRUE); //if false create swaped pair
  
 }
/************************************************************************/

Double_t AliHBTPair::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 AliHBTPair::GetQSideCMSLC()
{
  //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
 
  if (fQSideCMSLCNotCalc)
   {
    fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
    fQSideCMSLCNotCalc = kFALSE;
   }
  return fQSideCMSLC;
}
/************************************************************************/
Double_t AliHBTPair::GetQOutCMSLC()
{
 if(fQOutCMSLCNotCalc)
  {
   CalculateSums();
   CalculateDiffs();

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

      //beta=fPzSum/fESum;    // Longit. V == beta
      Double_t beta=fPzSum/fESum;
      Double_t gamma = GetGammaToCMSLC();
      
      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,GetGammaToCMSLC(),part1Pz,E1a);
//      LTR(part2Pz,E2,beta,GetGammaToCMSLC(),part2Pz,E2a);
      Double_t part1Pz=gamma*(fPart1->Pz()-beta*fPart1->Energy());
      Double_t part2Pz=gamma*(fPart2->Pz()-beta*fPart2->Energy());

      Double_t part1P2=part1Px*part1Px+part1Py*part1Py+part1Pz*part1Pz;
      Double_t part2P2=part2Px*part2Px+part2Py*part2Py+part2Pz*part2Pz;
      Double_t part1E=TMath::Sqrt(fPart1->GetMass()*fPart1->GetMass()+part1P2);
      Double_t part2E=TMath::Sqrt(fPart2->GetMass()*fPart2->GetMass()+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->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/(relmass*relmass);
      fQOutCMSLC=(part1Px-part2Px);//== id
      fQOutCMSLC=fQOutCMSLC-sumPx*hf; //sumPx == fPxSum ale po rotacji i transf
    }
   else
    {
      Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
      fQOutCMSLC = 0.5*k2/GetKt();
   // if (non-id)  fQOutCMSLC=fQOutCMSLC - sumPx*HF;
    }

    
   fQOutCMSLCNotCalc = kFALSE;
  }
 return fQOutCMSLC;
}
/************************************************************************/
Double_t AliHBTPair::GetQLongCMSLC()
{
 if (fQLongCMSLCNotCalc)
  {
    CalculateSums();
    CalculateDiffs();
    Double_t beta = fPzSum/fESum;
    fQLongCMSLC = GetGammaToCMSLC() * ( fPzDiff - beta*fEDiff );
    fQLongCMSLCNotCalc = kFALSE;
  }
 return fQLongCMSLC; 
}
/************************************************************************/
Double_t AliHBTPair::GetKt()
{
  if(fKtNotCalc)
   { 
     CalculateSums();
     fKt =  0.5*TMath::Hypot(fPxSum,fPySum);
     fKtNotCalc = kFALSE;
   }
  return fKt;
}
/************************************************************************/

Double_t AliHBTPair::GetKStar()
{
  if (fKStarNotCalc)
   { 
    CalculateSums();

    Double_t Ptrans = fPxSum*fPxSum + fPySum*fPySum;
    Double_t Mtrans = fESum*fESum - fPzSum*fPzSum;
    Double_t Pinv =   TMath::Sqrt(Mtrans - Ptrans);

    Double_t Q = (fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/Pinv;
    
    CalculateQInvL();
    
    Q = TMath::Sqrt( Q*Q - fQInvL);
    fKStar = Q/2.;
    fKStarNotCalc = kFALSE;
   }
  return fKStar;
}
/************************************************************************/

Double_t AliHBTPair::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 AliHBTPair::GetGammaToCMSLC()
{
  if(fGammaCMSLCNotCalc)
   {
     CalculateSums();
     Double_t beta = fPzSum/fESum;
     fGammaCMSLC = 1.0/TMath::Sqrt(1.0 - beta*beta);
     fGammaCMSLCNotCalc = kFALSE;
   }
  return fGammaCMSLC;
}
/************************************************************************/

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