Track Points Added. Class needed for Anti-Merging cut
[u/mrichter/AliRoot.git] / HBTAN / AliHBTParticle.cxx
CommitLineData
1b446896 1#include "AliHBTParticle.h"
f5de2f09 2//___________________________________________________________
3/////////////////////////////////////////////////////////////
4//
5// class AliHBTParticle
6//
7// Ali HBT Particle: simplified class TParticle
8// Simplified in order to minimize the size of object
9// - we want to keep a lot of such a objects in memory
10// Additionaly adjusted for HBT Analysies purposes
11//
12/////////////////////////////////////////////////////////////
1b446896 13#include <TParticle.h>
ff3d1148 14#include <TClass.h>
f5de2f09 15#include "AliHBTTrackPoints.h"
1b446896 16
17ClassImp(AliHBTParticle)
18
8089a083 19Int_t AliHBTParticle::fgDebug = 0;
1b446896 20//______________________________________________________________________________
21AliHBTParticle::AliHBTParticle():
8089a083 22 fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
f5de2f09 23 fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0),
24 fTrackPoints(0x0)
1b446896 25{//empty particle
26}
8089a083 27//______________________________________________________________________________
1b446896 28
8089a083 29AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
30 Double_t px, Double_t py, Double_t pz, Double_t etot,
31 Double_t vx, Double_t vy, Double_t vz, Double_t time):
e665964c 32 fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
8089a083 33 fCalcMass(0),
34 fPx(px), fPy(py),fPz(pz),fE(etot),
f5de2f09 35 fVx(vx), fVy(vy),fVz(vz),fVt(time),
36 fTrackPoints(0x0)
8089a083 37{
38//mormal constructor
39 SetPdgCode(pdg);
40 if (GetPDG()) {
41 fCalcMass = GetPDG()->Mass();
42 } else {
43 Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
44 if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
45 else fCalcMass = -TMath::Sqrt(-a2);
46 }
47}
1b446896 48//______________________________________________________________________________
8089a083 49
50AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx,
51 Double_t px, Double_t py, Double_t pz, Double_t etot,
52 Double_t vx, Double_t vy, Double_t vz, Double_t time):
e665964c 53 fPdgIdx(0), fIdxInEvent(idx),fNPids(0),fPids(0x0),fPidProb(0x0),
8089a083 54 fCalcMass(0),
35d14c47 55 fPx(px), fPy(py),fPz(pz),fE(etot),
f5de2f09 56 fVx(vx), fVy(vy),fVz(vz),fVt(time),
57 fTrackPoints(0x0)
1b446896 58{
59//mormal constructor
8089a083 60 SetPdgCode(pdg,prob);
1b446896 61 if (GetPDG()) {
62 fCalcMass = GetPDG()->Mass();
63 } else {
64 Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
65 if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
66 else fCalcMass = -TMath::Sqrt(-a2);
67 }
68}
8089a083 69//______________________________________________________________________________
70AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
bed069a4 71 TObject(in),
8089a083 72 fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
73 fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
74 fCalcMass(in.GetCalcMass()),
75 fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()),
76 fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T())
77{
78 //Copy constructor
79 for(Int_t i = 0; i<fNPids; i++)
80 {
81 fPids[i] = in.fPids[i];
82 fPidProb[i] = in.fPidProb[i];
83 }
f5de2f09 84
85 fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)fTrackPoints->Clone():0x0;
86
8089a083 87}
1b446896 88
89//______________________________________________________________________________
8089a083 90AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
91 fPdgIdx(0), fIdxInEvent(idx),
92 fNPids(0),fPids(0x0),fPidProb(0x0),
93 fCalcMass(p.GetCalcMass()),
1b446896 94 fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()),
f5de2f09 95 fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T()),
96 fTrackPoints(0x0)
1b446896 97{
98 //all copied in the initialization
8089a083 99 SetPdgCode(p.GetPdgCode());
100}
101//______________________________________________________________________________
102
bed069a4 103AliHBTParticle::~AliHBTParticle()
104{
105//dtor
106 delete [] fPids;
107 delete [] fPidProb;
f5de2f09 108 delete fTrackPoints;
bed069a4 109}
110//______________________________________________________________________________
111
112AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
113{
114//assigment operator
115
116 fNPids = in.fNPids;
117 delete [] fPids;
118 delete [] fPidProb;
119 Int_t* fPids = new Int_t[fNPids];
120 Float_t* fPidProb = new Float_t[fNPids];
121 for (Int_t i = 0; i < fNPids;i++)
122 {
123 fPids[i] = in.fPids[i];
124 fPidProb[i] = in.fPidProb[i];
125 }
126
127 fPdgIdx = in.fPdgIdx;
128 fIdxInEvent = in.fIdxInEvent;
129 fCalcMass = in.GetCalcMass();
130 fPx = in.Px();
131 fPy = in.Py();
132 fPz = in.Pz();
133 fE = in.Energy();
134 fVx = in.Vx();
135 fVy = in.Vy();
136 fVz = in.Vz();
137 fVt = in.T();
138
f5de2f09 139 delete fTrackPoints;
140 fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)fTrackPoints->Clone():0x0;
141
bed069a4 142 return *this;
143}
144//______________________________________________________________________________
145
8089a083 146void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
147{
148 SetPIDprobability(pdg,prob);
149 fPdgIdx = GetPidSlot(pdg);
1b446896 150}
151
152//______________________________________________________________________________
8089a083 153void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
154{
155//Sets another pdg code and corresponding probabilty
156//Ids are set in decreasing order
157//Check if total prbaility is not ivercoming unity is performed
158//in case, warning is printed
e665964c 159 if (fgDebug > 9) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
8089a083 160
161 Float_t totprob = 0.0;//sums up probabilities
162 Int_t idx = GetPidSlot(pdg);
163 Int_t i;
164 if (idx > -1)
165 {
166 fPidProb[idx] = prob;
167 for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
e665964c 168 if (totprob > (1.0+0.000001))
8089a083 169 {
e665964c 170 Warning("SetPIDprobability","Total probability greater than unity (%f)",totprob);
8089a083 171 }
e665964c 172 if (fgDebug > 9)
8089a083 173 {
174 Info("SetPIDprobability","Current Total probability: %f",totprob);
175 }
176 return;
177 }
178
179 Int_t currentpid = GetPdgCode();
180 fNPids++;
181 Float_t* aPidProbNew = new Float_t[fNPids];
182 Int_t* aPidsNew = new Int_t[fNPids];
183
184 for (i = 0; i < fNPids-1;i++)//find a slot
185 {
186 if ( fPidProb[i] > prob)
187 {
e665964c 188 if (fgDebug>9) Info("SetPID","Copying entry %d",i);
8089a083 189 aPidProbNew[i] = fPidProb[i];
190 aPidsNew[i] = fPids[i];
191 totprob+=fPidProb[i];
192 }
193 else break;
194 }
195
e665964c 196 if (fgDebug > 9) Info("SetPID","Setting new PID on entry %d",i);
8089a083 197 aPidProbNew[i] = prob;
198 aPidsNew[i] = pdg;
199 totprob+=prob;
200
201
e665964c 202 for (Int_t j = fNPids-1; j > i ;j--)//copy rest of old araays
8089a083 203 {
e665964c 204 if (fgDebug > 9) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
8089a083 205 aPidProbNew[j] = fPidProb[j-1];
206 aPidsNew[j] = fPids[j-1];
207 totprob+=fPidProb[j-1];
208 }
209
8089a083 210 delete [] fPidProb;
211 delete [] fPids;
212
213 fPidProb = aPidProbNew;
214 fPids = aPidsNew;
215
216 fPdgIdx = GetPidSlot(currentpid);
217 if (fPdgIdx == -1) fPdgIdx = 0;
e665964c 218
219 if (totprob > (1.0+0.000001))//place for numerical error
220 {
221 Warning("SetId","Total probability is greater than unity (%f)!!!",totprob);
222 Print();
223 }
8089a083 224}
225//______________________________________________________________________________
226
227Float_t AliHBTParticle::GetPIDprobability(Int_t pdg)
228{
229 Int_t idx = GetPidSlot(pdg);
230 if (idx < 0) return 0.0;//such pid was not specified for this particle
231 return fPidProb[idx];
232}
233//______________________________________________________________________________
234
1b446896 235const Char_t* AliHBTParticle::GetName() const
236{
237 static char def[4] = "XXX";
8089a083 238 const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
1b446896 239 if (ap) return ap->GetName();
240 else return def;
241}
242
243
244//______________________________________________________________________________
8089a083 245Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
246{
247 //returns position of the given PID in fPids (and fPidProb) array.
248 if (fPids == 0x0) return -1;
249 for (Int_t i = 0; i< fNPids; i++)
250 {
251 if (fPids[i] == pdg) return i;
252 }
253 return -1;
254}
e665964c 255//______________________________________________________________________________
256
257Int_t AliHBTParticle::GetNthPid(Int_t idx) const
258{
259 //returns PID sitting on slot idx in fPids
260 if ( (idx < 0) || (idx >= fNPids) )
261 {
262 Error("GetNthPid","Out Of Bounds");
263 return 0;
264 }
265 return fPids[idx];
266}
267//______________________________________________________________________________
1b446896 268
e665964c 269Float_t AliHBTParticle::GetNthPidProb(Int_t idx) const
270{
271 //returns PID sitting on slot idx in fPidProb
272 if ( (idx < 0) || (idx >= fNPids) )
273 {
274 Error("GetNthPid","Out Of Bounds");
275 return 0;
276 }
277 return fPidProb[idx];
278}
8089a083 279//______________________________________________________________________________
e665964c 280
8089a083 281void AliHBTParticle::Print() const
282{
283//prints information about particle
284 printf("____________________________________________________\n");
e665964c 285 printf("Idx: %d PID: %d Name: ",fIdxInEvent,GetPdgCode());
8089a083 286 TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
287 if (pdgp)
288 {
e665964c 289 printf("%s Mass: %f\n",pdgp->GetName(),pdgp->Mass());
8089a083 290 }
291 else
292 {
293 printf("Not known\n");
294 }
295
ff3d1148 296 printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
8089a083 297 Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
1b446896 298
8089a083 299 for (Int_t i = 0; i < fNPids; i++)
300 {
ff3d1148 301 printf("# %d PID: %d Probability %f name ",i,fPids[i],fPidProb[i]);
8089a083 302 const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
303 if (ap)
304 {
305 printf("%s Mass %f\n",ap->GetName(),ap->Mass());
306 }
307 else
308 {
309 printf("Not known\n");
310 }
311 }
312}
ff3d1148 313
314//______________________________________________________________________________
315
316//void AliHBTParticle::Streamer(TBuffer &b)
317//{
318// // Stream all objects in the array to or from the I/O buffer.
319// UInt_t R__s, R__c;
320// Int_t i;
321// if (b.IsReading())
322// {
323// delete [] fPids;
324// delete [] fPidProb;
325//
326// Version_t v = b.ReadVersion(&R__s, &R__c);
327// if (v == 1)
328// {
329// AliHBTParticle::Class()->ReadBuffer(b, this);
330// }
331// else
332// {
333// TObject::Streamer(b);
334// b >> fPdgIdx;
335// b >> fIdxInEvent;
336//
337// b >> fNPids;
338// Int_t* fPids = new Int_t[fNPids];
339// Float_t* fPidProb = new Float_t[fNPids];
340// for (i = 0;i<fNPids;i++)
341// {
342// b >> fPids[i];
343// }
344// for (i = 0;i<fNPids;i++)
345// {
346// b >> fPidProb[i];
347// }
348// b >> fCalcMass;
349//
350// b >> fPx;
351// b >> fPy;
352// b >> fPz;
353// b >> fE;
354//
355// b >> fVx;
356// b >> fVy;
357// b >> fVz;
358// b >> fVt;
359// Info("Streamer","Read data");
360// Print();
361// }
362//
363// b.CheckByteCount(R__s, R__c,AliHBTParticle::IsA());
364// }
365// else
366// {
367// R__c = b.WriteVersion(AliHBTParticle::IsA(), kTRUE);
368// TObject::Streamer(b);
369// Info("Streamer","Read data");
370// Print();
371//
372// b << fPdgIdx;
373// b << fIdxInEvent;
374// b << fNPids;
375// for (i = 0;i<fNPids;i++)
376// {
377// b << fPids[i];
378// }
379// {
380// {
381// for (i = 0;i<fNPids;i++)
382// {
383// b << fPidProb[i];
384// }
385// b << fCalcMass;
386//
387// b << fPx;
388// b << fPy;
389// b << fPz;
390// b << fE;
391//
392// b << fVx;
393// b << fVy;
394// b << fVz;
395// b << fVt;
396//
397// b.SetByteCount(R__c, kTRUE);
398// }
399//}