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

//Simplified TParticle class
#include "AliHBTParticle.h"
#include <TParticle.h>
#include <TClass.h>

ClassImp(AliHBTParticle)

Int_t AliHBTParticle::fgDebug = 0;
//______________________________________________________________________________
AliHBTParticle::AliHBTParticle():  
 fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
 fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
{//empty particle
}
//______________________________________________________________________________

AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
               Double_t px, Double_t py, Double_t pz, Double_t etot,
               Double_t vx, Double_t vy, Double_t vz, Double_t time):  
  fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
  fCalcMass(0), 
  fPx(px), fPy(py),fPz(pz),fE(etot), 
  fVx(vx), fVy(vy),fVz(vz),fVt(time)
{
//mormal constructor
  SetPdgCode(pdg);
  if (GetPDG()) {
     fCalcMass    = GetPDG()->Mass();
  } else {
     Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
     if (a2 >= 0) fCalcMass =  TMath::Sqrt(a2);
     else         fCalcMass = -TMath::Sqrt(-a2);
  }
}
//______________________________________________________________________________

AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx, 
                               Double_t px, Double_t py, Double_t pz, Double_t etot,
                               Double_t vx, Double_t vy, Double_t vz, Double_t time):
  fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
  fCalcMass(0), 
  fPx(px), fPy(py),fPz(pz),fE(etot), 
  fVx(vx), fVy(vy),fVz(vz),fVt(time)
{
//mormal constructor
  SetPdgCode(pdg,prob);
  if (GetPDG()) {
     fCalcMass    = GetPDG()->Mass();
  } else {
     Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
     if (a2 >= 0) fCalcMass =  TMath::Sqrt(a2);
     else         fCalcMass = -TMath::Sqrt(-a2);
  }
}
//______________________________________________________________________________
AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
   fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
   fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
   fCalcMass(in.GetCalcMass()),
   fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()), 
   fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T())
{
 //Copy constructor
 for(Int_t i = 0; i<fNPids; i++)
  {
    fPids[i] =  in.fPids[i];
    fPidProb[i] = in.fPidProb[i];
  }
}

//______________________________________________________________________________
AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
   fPdgIdx(0), fIdxInEvent(idx),
   fNPids(0),fPids(0x0),fPidProb(0x0),
   fCalcMass(p.GetCalcMass()),
   fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
   fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
{
 //all copied in the initialization
 SetPdgCode(p.GetPdgCode());
}
//______________________________________________________________________________

void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
{
  SetPIDprobability(pdg,prob);
  fPdgIdx = GetPidSlot(pdg);
}

//______________________________________________________________________________
void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
{
//Sets another pdg code and corresponding probabilty
//Ids are set in decreasing order
//Check if total prbaility is not ivercoming unity is performed
//in case, warning is printed
  if (fgDebug) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);

  Float_t totprob = 0.0;//sums up probabilities
  Int_t idx = GetPidSlot(pdg);
  Int_t i;
  if (idx > -1) 
   {
     fPidProb[idx] = prob;
     for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
     if (totprob > 1.0)
       {
         Warning("SetPIDprobability","Total probability greater than UNITY: %f",totprob);
       }
     if (fgDebug) 
      {
        Info("SetPIDprobability","Current Total probability: %f",totprob);
      }
     return;
   }
    
  Int_t currentpid = GetPdgCode();
  fNPids++;
  Float_t* aPidProbNew = new Float_t[fNPids];
  Int_t* aPidsNew = new Int_t[fNPids];
  
  for (i = 0; i < fNPids-1;i++)//find a slot
   {
     if ( fPidProb[i] > prob)
      {
        if (fgDebug>4) Info("SetPID","Copying entry %d",i);
        aPidProbNew[i] = fPidProb[i];
        aPidsNew[i] = fPids[i];
        totprob+=fPidProb[i];
      }
     else break;
   }

  if (fgDebug>4) Info("SetPID","Setting new PID on entry %d",i);
  aPidProbNew[i] = prob;
  aPidsNew[i] = pdg;
  totprob+=prob;
  

  for (Int_t j = fNPids-1; j > i ;i--)//copy rest of old araays 
   {
     if (fgDebug>4) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
     aPidProbNew[j] = fPidProb[j-1];
     aPidsNew[j] = fPids[j-1];
     totprob+=fPidProb[j-1];
   }

  if (totprob > 1.0)
   {
     Warning("SetId","Total probability is greater than 1 !!!");
     Print();
   }
  delete [] fPidProb;
  delete [] fPids;
  
  fPidProb = aPidProbNew;
  fPids = aPidsNew;
  
  fPdgIdx = GetPidSlot(currentpid);
  if (fPdgIdx == -1) fPdgIdx = 0;
}
//______________________________________________________________________________

Float_t AliHBTParticle::GetPIDprobability(Int_t pdg)
{
  Int_t idx = GetPidSlot(pdg);
  if (idx < 0) return 0.0;//such pid was not specified for this particle
  return fPidProb[idx];
}
//______________________________________________________________________________

const Char_t* AliHBTParticle::GetName() const 
{
   static char def[4] = "XXX";
   const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
   if (ap) return ap->GetName();
   else    return def;
}


//______________________________________________________________________________
Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
{
 //returns position of the given PID in fPids (and fPidProb) array.
 if (fPids == 0x0) return -1;
 for (Int_t i = 0; i< fNPids; i++)
  {
   if (fPids[i] == pdg) return i;
  }
 return -1;
}

//______________________________________________________________________________
void AliHBTParticle::Print() const
{
//prints information about particle
  printf("____________________________________________________\n");
  printf("Idx: %d  PID: %d  Name",fIdxInEvent,GetPdgCode());
  TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
  if (pdgp)
   {
     printf("%s Mass %f\n",pdgp->GetName(),pdgp->Mass());
   }
  else
   {
     printf("Not known\n");
   }
  
  printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
          Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());

  for (Int_t i = 0; i < fNPids; i++)
   {
     printf("# %d  PID: %d  Probability %f name ",i,fPids[i],fPidProb[i]);
     const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
     if (ap)
      {
        printf("%s Mass %f\n",ap->GetName(),ap->Mass());
      }
     else
      {
        printf("Not known\n");
      }
   }
}

//______________________________________________________________________________

//void AliHBTParticle::Streamer(TBuffer &b)
//{
//     // Stream all objects in the array to or from the I/O buffer.
//   UInt_t R__s, R__c;
//   Int_t i;
//   if (b.IsReading()) 
//    {
//      delete [] fPids;
//      delete [] fPidProb;
//      
//      Version_t v = b.ReadVersion(&R__s, &R__c);
//      if (v == 1)
//       {
//         AliHBTParticle::Class()->ReadBuffer(b, this);
//      }
//      else
//       {
//        TObject::Streamer(b);
//       b >> fPdgIdx;
//       b >> fIdxInEvent;
//       
//       b >> fNPids;
//       Int_t* fPids = new Int_t[fNPids];
//        Float_t* fPidProb = new Float_t[fNPids];
//        for (i = 0;i<fNPids;i++) 
//         {
//           b >> fPids[i];
//         }
//        for (i = 0;i<fNPids;i++) 
//        {
//          b >> fPidProb[i];
//         }
//        b >> fCalcMass;
//
//        b >> fPx;
//       b >> fPy;
//       b >> fPz;
//        b >> fE;
//
//       b >> fVx;
//        b >> fVy;
//        b >> fVz;
//       b >> fVt;
//       Info("Streamer","Read data");
//        Print();
//       }
//
//      b.CheckByteCount(R__s, R__c,AliHBTParticle::IsA());
//    } 
//  else 
//   {
//     R__c = b.WriteVersion(AliHBTParticle::IsA(), kTRUE);
//     TObject::Streamer(b);
//     Info("Streamer","Read data");
//     Print();
//
//     b << fPdgIdx;
//     b << fIdxInEvent;
//     b << fNPids;
//     for (i = 0;i<fNPids;i++) 
//      {
//        b << fPids[i];
//      }
//      {
//      {
//     for (i = 0;i<fNPids;i++) 
//     {
//        b << fPidProb[i];
//      }
//     b << fCalcMass;
//
//     b << fPx;
//     b << fPy;
//     b << fPz;
//     b << fE;
//
//     b << fVx;
//     b << fVy;
//     b << fVz;
//     b << fVt;
//
//    b.SetByteCount(R__c, kTRUE);
//   }
//}
Commit	Line	Data
1b446896	1	//Simplified TParticle class
1b446896	2	#include "AliHBTParticle.h"
1b446896	3	#include <TParticle.h>
ff3d1148	4	#include <TClass.h>
1b446896	5
	6	ClassImp(AliHBTParticle)
	7
8089a083	8	Int_t AliHBTParticle::fgDebug = 0;
1b446896	9	//______________________________________________________________________________
1b446896	10	AliHBTParticle::AliHBTParticle():
8089a083	11	fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
8089a083	12	fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
1b446896	13	{//empty particle
1b446896	14	}
8089a083	15	//______________________________________________________________________________
1b446896	16
8089a083	17	AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
	18	Double_t px, Double_t py, Double_t pz, Double_t etot,
	19	Double_t vx, Double_t vy, Double_t vz, Double_t time):
	20	fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
	21	fCalcMass(0),
	22	fPx(px), fPy(py),fPz(pz),fE(etot),
	23	fVx(vx), fVy(vy),fVz(vz),fVt(time)
	24	{
	25	//mormal constructor
	26	SetPdgCode(pdg);
	27	if (GetPDG()) {
	28	fCalcMass = GetPDG()->Mass();
	29	} else {
	30	Double_t a2 = fEfE -fPxfPx -fPyfPy -fPzfPz;
	31	if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
	32	else fCalcMass = -TMath::Sqrt(-a2);
	33	}
	34	}
1b446896	35	//______________________________________________________________________________
8089a083	36
	37	AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx,
	38	Double_t px, Double_t py, Double_t pz, Double_t etot,
	39	Double_t vx, Double_t vy, Double_t vz, Double_t time):
	40	fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
	41	fCalcMass(0),
35d14c47	42	fPx(px), fPy(py),fPz(pz),fE(etot),
35d14c47	43	fVx(vx), fVy(vy),fVz(vz),fVt(time)
1b446896	44	{
1b446896	45	//mormal constructor
8089a083	46	SetPdgCode(pdg,prob);
1b446896	47	if (GetPDG()) {
	48	fCalcMass = GetPDG()->Mass();
	49	} else {
	50	Double_t a2 = fEfE -fPxfPx -fPyfPy -fPzfPz;
	51	if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
	52	else fCalcMass = -TMath::Sqrt(-a2);
	53	}
	54	}
8089a083	55	//______________________________________________________________________________
	56	AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
	57	fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
	58	fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
	59	fCalcMass(in.GetCalcMass()),
	60	fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()),
	61	fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T())
	62	{
	63	//Copy constructor
	64	for(Int_t i = 0; i<fNPids; i++)
	65	{
	66	fPids[i] = in.fPids[i];
	67	fPidProb[i] = in.fPidProb[i];
	68	}
	69	}
1b446896	70
1b446896	71	//______________________________________________________________________________
8089a083	72	AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
	73	fPdgIdx(0), fIdxInEvent(idx),
	74	fNPids(0),fPids(0x0),fPidProb(0x0),
	75	fCalcMass(p.GetCalcMass()),
1b446896	76	fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()),
	77	fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
	78	{
	79	//all copied in the initialization
8089a083	80	SetPdgCode(p.GetPdgCode());
	81	}
	82	//______________________________________________________________________________
	83
	84	void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
	85	{
	86	SetPIDprobability(pdg,prob);
	87	fPdgIdx = GetPidSlot(pdg);
1b446896	88	}
	89
	90	//______________________________________________________________________________
8089a083	91	void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
	92	{
	93	//Sets another pdg code and corresponding probabilty
	94	//Ids are set in decreasing order
	95	//Check if total prbaility is not ivercoming unity is performed
	96	//in case, warning is printed
	97	if (fgDebug) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
	98
	99	Float_t totprob = 0.0;//sums up probabilities
	100	Int_t idx = GetPidSlot(pdg);
	101	Int_t i;
	102	if (idx > -1)
	103	{
	104	fPidProb[idx] = prob;
	105	for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
	106	if (totprob > 1.0)
	107	{
	108	Warning("SetPIDprobability","Total probability greater than UNITY: %f",totprob);
	109	}
	110	if (fgDebug)
	111	{
	112	Info("SetPIDprobability","Current Total probability: %f",totprob);
	113	}
	114	return;
	115	}
	116
	117	Int_t currentpid = GetPdgCode();
	118	fNPids++;
	119	Float_t* aPidProbNew = new Float_t[fNPids];
	120	Int_t* aPidsNew = new Int_t[fNPids];
	121
	122	for (i = 0; i < fNPids-1;i++)//find a slot
	123	{
	124	if ( fPidProb[i] > prob)
	125	{
	126	if (fgDebug>4) Info("SetPID","Copying entry %d",i);
	127	aPidProbNew[i] = fPidProb[i];
	128	aPidsNew[i] = fPids[i];
	129	totprob+=fPidProb[i];
	130	}
	131	else break;
	132	}
	133
	134	if (fgDebug>4) Info("SetPID","Setting new PID on entry %d",i);
	135	aPidProbNew[i] = prob;
	136	aPidsNew[i] = pdg;
	137	totprob+=prob;
	138
	139
	140	for (Int_t j = fNPids-1; j > i ;i--)//copy rest of old araays
	141	{
	142	if (fgDebug>4) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
	143	aPidProbNew[j] = fPidProb[j-1];
	144	aPidsNew[j] = fPids[j-1];
	145	totprob+=fPidProb[j-1];
	146	}
	147
	148	if (totprob > 1.0)
	149	{
	150	Warning("SetId","Total probability is greater than 1 !!!");
	151	Print();
	152	}
	153	delete [] fPidProb;
	154	delete [] fPids;
155
156	fPidProb = aPidProbNew;
157	fPids = aPidsNew;
158
159	fPdgIdx = GetPidSlot(currentpid);
160	if (fPdgIdx == -1) fPdgIdx = 0;
161	}
162	//______________________________________________________________________________
163
164	Float_t AliHBTParticle::GetPIDprobability(Int_t pdg)
165	{
166	Int_t idx = GetPidSlot(pdg);
167	if (idx < 0) return 0.0;//such pid was not specified for this particle
168	return fPidProb[idx];
169	}
170	//______________________________________________________________________________
171
1b446896	172	const Char_t* AliHBTParticle::GetName() const
	173	{
	174	static char def[4] = "XXX";
8089a083	175	const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
1b446896	176	if (ap) return ap->GetName();
	177	else return def;
	178	}
	179
	180
	181	//______________________________________________________________________________
8089a083	182	Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
	183	{
	184	//returns position of the given PID in fPids (and fPidProb) array.
	185	if (fPids == 0x0) return -1;
	186	for (Int_t i = 0; i< fNPids; i++)
	187	{
	188	if (fPids[i] == pdg) return i;
	189	}
	190	return -1;
	191	}
1b446896	192
8089a083	193	//______________________________________________________________________________
	194	void AliHBTParticle::Print() const
	195	{
	196	//prints information about particle
	197	printf("____________________________________________________\n");
	198	printf("Idx: %d PID: %d Name",fIdxInEvent,GetPdgCode());
	199	TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
	200	if (pdgp)
	201	{
	202	printf("%s Mass %f\n",pdgp->GetName(),pdgp->Mass());
	203	}
	204	else
	205	{
	206	printf("Not known\n");
	207	}
	208
ff3d1148	209	printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
8089a083	210	Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
1b446896	211
8089a083	212	for (Int_t i = 0; i < fNPids; i++)
8089a083	213	{
ff3d1148	214	printf("# %d PID: %d Probability %f name ",i,fPids[i],fPidProb[i]);
8089a083	215	const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
	216	if (ap)
	217	{
	218	printf("%s Mass %f\n",ap->GetName(),ap->Mass());
	219	}
	220	else
	221	{
	222	printf("Not known\n");
	223	}
	224	}
	225	}
ff3d1148	226
	227	//______________________________________________________________________________
	228
	229	//void AliHBTParticle::Streamer(TBuffer &b)
	230	//{
	231	// // Stream all objects in the array to or from the I/O buffer.
	232	// UInt_t R__s, R__c;
	233	// Int_t i;
	234	// if (b.IsReading())
	235	// {
	236	// delete [] fPids;
	237	// delete [] fPidProb;
	238	//
	239	// Version_t v = b.ReadVersion(&R__s, &R__c);
	240	// if (v == 1)
	241	// {
	242	// AliHBTParticle::Class()->ReadBuffer(b, this);
	243	// }
	244	// else
	245	// {
	246	// TObject::Streamer(b);
	247	// b >> fPdgIdx;
	248	// b >> fIdxInEvent;
	249	//
	250	// b >> fNPids;
	251	// Int_t* fPids = new Int_t[fNPids];
	252	// Float_t* fPidProb = new Float_t[fNPids];
	253	// for (i = 0;i<fNPids;i++)
	254	// {
	255	// b >> fPids[i];
	256	// }
	257	// for (i = 0;i<fNPids;i++)
	258	// {
	259	// b >> fPidProb[i];
	260	// }
	261	// b >> fCalcMass;
	262	//
	263	// b >> fPx;
	264	// b >> fPy;
	265	// b >> fPz;
	266	// b >> fE;
	267	//
	268	// b >> fVx;
	269	// b >> fVy;
	270	// b >> fVz;
	271	// b >> fVt;
	272	// Info("Streamer","Read data");
	273	// Print();
	274	// }
	275	//
	276	// b.CheckByteCount(R__s, R__c,AliHBTParticle::IsA());
	277	// }
	278	// else
	279	// {
	280	// R__c = b.WriteVersion(AliHBTParticle::IsA(), kTRUE);
	281	// TObject::Streamer(b);
	282	// Info("Streamer","Read data");
	283	// Print();
	284	//
	285	// b << fPdgIdx;
	286	// b << fIdxInEvent;
	287	// b << fNPids;
	288	// for (i = 0;i<fNPids;i++)
	289	// {
290	// b << fPids[i];
291	// }
292	// {
293	// {
294	// for (i = 0;i<fNPids;i++)
295	// {
296	// b << fPidProb[i];
297	// }
298	// b << fCalcMass;
299	//
300	// b << fPx;
301	// b << fPy;
302	// b << fPz;
303	// b << fE;
304	//
305	// b << fVx;
306	// b << fVy;
307	// b << fVz;
308	// b << fVt;
309	//
310	// b.SetByteCount(R__c, kTRUE);
311	// }
312	//}