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