Catching up to NewIO -> Particle stores all passible PID and their probabilities
[u/mrichter/AliRoot.git] / HBTAN / AliHBTParticle.cxx
1 //Simplified TParticle class
2 #include "AliHBTParticle.h"
3 #include <TParticle.h>
4
5 ClassImp(AliHBTParticle)
6
7 Int_t AliHBTParticle::fgDebug = 0;
8 //______________________________________________________________________________
9 AliHBTParticle::AliHBTParticle():  
10  fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
11  fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
12 {//empty particle
13 }
14 //______________________________________________________________________________
15
16 AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
17                Double_t px, Double_t py, Double_t pz, Double_t etot,
18                Double_t vx, Double_t vy, Double_t vz, Double_t time):  
19   fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
20   fCalcMass(0), 
21   fPx(px), fPy(py),fPz(pz),fE(etot), 
22   fVx(vx), fVy(vy),fVz(vz),fVt(time)
23 {
24 //mormal constructor
25   SetPdgCode(pdg);
26   if (GetPDG()) {
27      fCalcMass    = GetPDG()->Mass();
28   } else {
29      Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
30      if (a2 >= 0) fCalcMass =  TMath::Sqrt(a2);
31      else         fCalcMass = -TMath::Sqrt(-a2);
32   }
33 }
34 //______________________________________________________________________________
35
36 AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx, 
37                                Double_t px, Double_t py, Double_t pz, Double_t etot,
38                                Double_t vx, Double_t vy, Double_t vz, Double_t time):
39   fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
40   fCalcMass(0), 
41   fPx(px), fPy(py),fPz(pz),fE(etot), 
42   fVx(vx), fVy(vy),fVz(vz),fVt(time)
43 {
44 //mormal constructor
45   SetPdgCode(pdg,prob);
46   if (GetPDG()) {
47      fCalcMass    = GetPDG()->Mass();
48   } else {
49      Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
50      if (a2 >= 0) fCalcMass =  TMath::Sqrt(a2);
51      else         fCalcMass = -TMath::Sqrt(-a2);
52   }
53 }
54 //______________________________________________________________________________
55 AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
56    fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
57    fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
58    fCalcMass(in.GetCalcMass()),
59    fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()), 
60    fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T())
61 {
62  //Copy constructor
63  for(Int_t i = 0; i<fNPids; i++)
64   {
65     fPids[i] =  in.fPids[i];
66     fPidProb[i] = in.fPidProb[i];
67   }
68 }
69
70 //______________________________________________________________________________
71 AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
72    fPdgIdx(0), fIdxInEvent(idx),
73    fNPids(0),fPids(0x0),fPidProb(0x0),
74    fCalcMass(p.GetCalcMass()),
75    fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
76    fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
77 {
78  //all copied in the initialization
79  SetPdgCode(p.GetPdgCode());
80 }
81 //______________________________________________________________________________
82
83 void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
84 {
85   SetPIDprobability(pdg,prob);
86   fPdgIdx = GetPidSlot(pdg);
87 }
88
89 //______________________________________________________________________________
90 void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
91 {
92 //Sets another pdg code and corresponding probabilty
93 //Ids are set in decreasing order
94 //Check if total prbaility is not ivercoming unity is performed
95 //in case, warning is printed
96   if (fgDebug) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
97
98   Float_t totprob = 0.0;//sums up probabilities
99   Int_t idx = GetPidSlot(pdg);
100   Int_t i;
101   if (idx > -1) 
102    {
103      fPidProb[idx] = prob;
104      for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
105      if (totprob > 1.0)
106        {
107          Warning("SetPIDprobability","Total probability greater than UNITY: %f",totprob);
108        }
109      if (fgDebug) 
110       {
111         Info("SetPIDprobability","Current Total probability: %f",totprob);
112       }
113      return;
114    }
115     
116   Int_t currentpid = GetPdgCode();
117   fNPids++;
118   Float_t* aPidProbNew = new Float_t[fNPids];
119   Int_t* aPidsNew = new Int_t[fNPids];
120   
121   for (i = 0; i < fNPids-1;i++)//find a slot
122    {
123      if ( fPidProb[i] > prob)
124       {
125         if (fgDebug>4) Info("SetPID","Copying entry %d",i);
126         aPidProbNew[i] = fPidProb[i];
127         aPidsNew[i] = fPids[i];
128         totprob+=fPidProb[i];
129       }
130      else break;
131    }
132
133   if (fgDebug>4) Info("SetPID","Setting new PID on entry %d",i);
134   aPidProbNew[i] = prob;
135   aPidsNew[i] = pdg;
136   totprob+=prob;
137   
138
139   for (Int_t j = fNPids-1; j > i ;i--)//copy rest of old araays 
140    {
141      if (fgDebug>4) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
142      aPidProbNew[j] = fPidProb[j-1];
143      aPidsNew[j] = fPids[j-1];
144      totprob+=fPidProb[j-1];
145    }
146
147   if (totprob > 1.0)
148    {
149      Warning("SetId","Total probability is greater than 1 !!!");
150      Print();
151    }
152   delete [] fPidProb;
153   delete [] fPids;
154   
155   fPidProb = aPidProbNew;
156   fPids = aPidsNew;
157   
158   fPdgIdx = GetPidSlot(currentpid);
159   if (fPdgIdx == -1) fPdgIdx = 0;
160 }
161 //______________________________________________________________________________
162
163 Float_t AliHBTParticle::GetPIDprobability(Int_t pdg)
164 {
165   Int_t idx = GetPidSlot(pdg);
166   if (idx < 0) return 0.0;//such pid was not specified for this particle
167   return fPidProb[idx];
168 }
169 //______________________________________________________________________________
170
171 const Char_t* AliHBTParticle::GetName() const 
172 {
173    static char def[4] = "XXX";
174    const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
175    if (ap) return ap->GetName();
176    else    return def;
177 }
178
179
180 //______________________________________________________________________________
181 Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
182 {
183  //returns position of the given PID in fPids (and fPidProb) array.
184  if (fPids == 0x0) return -1;
185  for (Int_t i = 0; i< fNPids; i++)
186   {
187    if (fPids[i] == pdg) return i;
188   }
189  return -1;
190 }
191
192 //______________________________________________________________________________
193 void AliHBTParticle::Print() const
194 {
195 //prints information about particle
196   printf("____________________________________________________\n");
197   printf("Idx: %d  PID: %d  Name",fIdxInEvent,GetPdgCode());
198   TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
199   if (pdgp)
200    {
201      printf("%s Mass %f\n",pdgp->GetName(),pdgp->Mass());
202    }
203   else
204    {
205      printf("Not known\n");
206    }
207   
208   printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f Vx: %+f Vy: %+f Vz: %+f T: %+f",
209           Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
210
211   for (Int_t i = 0; i < fNPids; i++)
212    {
213      printf("# %d  PID: %d  Probability %f name",i,fPids[i],fPidProb[i]);
214      const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
215      if (ap)
216       {
217         printf("%s Mass %f\n",ap->GetName(),ap->Mass());
218       }
219      else
220       {
221         printf("Not known\n");
222       }
223    }
224 }