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

#include "AliHBTPair.h"
#include "AliHBTParticle.h"
#include "AliHBTLLWeights.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),
 fLLWeight(0.0),
 ffLLWeightNotCalc(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),
 fLLWeight(0.0),
 ffLLWeightNotCalc(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;
}
/************************************************************************/

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