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