Corrected operator overloading (alpha)
[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  fTPCTrackPoints(0x0),fITSTrackPoints(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   fTPCTrackPoints(0x0),fITSTrackPoints(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   fTPCTrackPoints(0x0),fITSTrackPoints(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
74 AliAODParticle::AliAODParticle(const AliAODParticle& in):
75    AliVAODParticle(in),
76    fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
77    fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
78    fCalcMass(in.GetCalcMass()),
79    fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.E()), 
80    fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T()),
81    fTPCTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
82 {
83  //Copy constructor
84 // Info("AliAODParticle(const AliAODParticle& in)","");
85  for(Int_t i = 0; i<fNPids; i++)
86   {
87     fPids[i] =  in.fPids[i];
88     fPidProb[i] = in.fPidProb[i];
89   }
90  
91  if (in.fTPCTrackPoints)
92    fTPCTrackPoints = (AliTrackPoints*)in.fTPCTrackPoints->Clone();
93  if (in.fITSTrackPoints)
94    fITSTrackPoints = (AliTrackPoints*)in.fITSTrackPoints->Clone();
95  if (in.fClusterMap)
96    fClusterMap = (AliClusterMap*)in.fClusterMap->Clone();
97 }
98 //______________________________________________________________________________
99
100 AliAODParticle::AliAODParticle(const AliVAODParticle& in):
101    AliVAODParticle(in),
102    fPdgIdx(0), fIdxInEvent(in.GetUID()),
103    fNPids(0),fPids(0x0),fPidProb(0x0),
104    fCalcMass(-1.0),
105    fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.E()), 
106    fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T()),
107    fTPCTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
108 {
109  //Copy constructor
110 // Info("AliAODParticle(const AliVAODParticle& in)","");
111  for(Int_t i = 0; i<in.GetNumberOfPids(); i++)
112   {
113     SetPIDprobability(in.GetNthPid(i),in.GetNthPidProb(i));
114   }
115  SetPdgCode(in.GetPdgCode(),in.GetPidProb());
116  
117  AliTrackPoints* tpts = in.GetTPCTrackPoints();
118  if (tpts)  SetTPCTrackPoints((AliTrackPoints*)tpts->Clone());
119  
120  tpts = in.GetITSTrackPoints();  
121  if (tpts) SetITSTrackPoints((AliTrackPoints*)tpts->Clone());
122    
123  AliClusterMap* clmap = in.GetClusterMap();
124  if (clmap) SetClusterMap((AliClusterMap*)clmap->Clone());
125 }
126 //______________________________________________________________________________
127
128 AliAODParticle::AliAODParticle(const TParticle &p,Int_t idx):
129    fPdgIdx(0), fIdxInEvent(idx),
130    fNPids(0),fPids(0x0),fPidProb(0x0),
131    fCalcMass(p.GetCalcMass()),
132    fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
133    fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T()),
134    fTPCTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
135 {
136  //all copied in the initialization
137  SetPdgCode(p.GetPdgCode());
138 }
139 //______________________________________________________________________________
140
141 AliAODParticle::~AliAODParticle()
142 {
143 //dtor  
144   delete [] fPids;
145   delete [] fPidProb;
146   delete fTPCTrackPoints;
147   delete fITSTrackPoints;
148   delete fClusterMap;
149 }
150 //______________________________________________________________________________
151
152 void AliAODParticle::Clear(Option_t*)
153 {
154 //Must be implemented in order to store this object in Clones Array
155   delete [] fPids;
156   delete [] fPidProb;
157   delete fTPCTrackPoints;
158   delete fITSTrackPoints;
159   delete fClusterMap;
160   
161   fPids = 0x0;
162   fPidProb = 0x0;
163   fTPCTrackPoints = 0x0;
164   fITSTrackPoints = 0x0;
165   fClusterMap = 0x0;
166 }
167 //______________________________________________________________________________
168
169 AliVAODParticle& AliAODParticle::operator=(const AliVAODParticle& in)
170 {
171 //operator=
172 //  Info("operator=(const AliVAODParticle& in)","AliAODParticle");
173   
174   if (&in == this) return *this;
175
176   delete [] fPids;
177   delete [] fPidProb;
178   fPids = 0x0;
179   fPidProb = 0x0;
180   fNPids = 0;
181
182   Int_t npids = in.GetNumberOfPids();
183   for (Int_t i = 0; i < npids; i++)
184    {
185      SetPIDprobability(in.GetNthPid(i),in.GetNthPidProb(i));
186    }
187    
188   SetPdgCode(in.GetPdgCode(),in.GetPidProb());
189
190   SetUID(in.GetUID());
191
192   fCalcMass = in.Mass();
193   
194   fPx = in.Px();
195   fPy = in.Py();
196   fPz = in.Pz();
197   fE = in.E(); 
198   fVx = in.Vx();
199   fVy = in.Vy();
200   fVz = in.Vz();
201   fVt = in.T();
202   
203   delete fTPCTrackPoints;
204   AliTrackPoints* tpts = in.GetTPCTrackPoints();
205   fTPCTrackPoints = (tpts)?(AliTrackPoints*)tpts->Clone():0x0;
206
207   delete fITSTrackPoints;
208   tpts = in.GetITSTrackPoints();  
209   fITSTrackPoints = (tpts)?(AliTrackPoints*)tpts->Clone():0x0;
210   
211   delete fClusterMap;
212   AliClusterMap* incmap = in.GetClusterMap();
213   fClusterMap =  (incmap)?(AliClusterMap*)incmap->Clone():0x0;
214   
215   return *this;
216 }
217 //______________________________________________________________________________
218
219 AliAODParticle& AliAODParticle::operator=(const AliAODParticle& in)
220 {
221 //assigment operator
222 //  Info("operator=(const AliAODParticle& in)","AliAODParticle");
223   if (&in == this) return *this;
224   fNPids = in.fNPids;
225   delete [] fPids;
226   delete [] fPidProb;
227   fPids = new Int_t[fNPids];
228   fPidProb = new Float_t[fNPids];
229   for (Int_t i = 0; i < fNPids;i++)
230    {
231      fPids[i]    = in.fPids[i];
232      fPidProb[i] = in.fPidProb[i];
233    }
234    
235   fPdgIdx = in.fPdgIdx; 
236   fIdxInEvent = in.fIdxInEvent;  
237   fCalcMass = in.GetCalcMass();
238   fPx = in.Px();
239   fPy = in.Py();
240   fPz = in.Pz();
241   fE = in.E(); 
242   fVx = in.Vx();
243   fVy = in.Vy();
244   fVz = in.Vz();
245   fVt = in.T();
246   
247   delete fTPCTrackPoints;
248   fTPCTrackPoints = (in.fTPCTrackPoints)?(AliTrackPoints*)in.fTPCTrackPoints->Clone():0x0;
249
250   delete fITSTrackPoints;
251   fITSTrackPoints = (in.fITSTrackPoints)?(AliTrackPoints*)in.fITSTrackPoints->Clone():0x0;
252   
253   delete fClusterMap;
254   fClusterMap =  (in.fClusterMap)?(AliClusterMap*)in.fClusterMap->Clone():0x0;
255   
256   return *this;
257 }
258 //______________________________________________________________________________
259
260 void AliAODParticle::SetPdgCode(Int_t pdg,Float_t prob)
261 {
262 //Set PDG Code
263   SetPIDprobability(pdg,prob);
264   fPdgIdx = GetPidSlot(pdg);
265 }
266
267 //______________________________________________________________________________
268 void AliAODParticle::SetPIDprobability(Int_t pdg, Float_t prob)
269 {
270 //Sets another pdg code and corresponding probabilty
271 //Ids are set in decreasing order
272 //Check if total probability is not overcoming unity is performed
273 //in case, warning is printed
274   if (GetDebug() > 9) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
275
276   Float_t totprob = 0.0;//sums up probabilities
277   Int_t idx = GetPidSlot(pdg);
278   Int_t i;
279
280   if (idx > -1) 
281    {
282      fPidProb[idx] = prob;
283      for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
284      if (totprob > (1.0+0.000001))
285        {
286          Warning("SetPIDprobability","Total probability greater than unity (%f)",totprob);
287        }
288      if (GetDebug() > 9) 
289       {
290         Info("SetPIDprobability","Current Total probability: %f",totprob);
291       }
292      return;
293    }
294     
295   Int_t currentpid = GetPdgCode();
296   fNPids++;
297   Float_t* aPidProbNew = new Float_t[fNPids];
298   Int_t* aPidsNew = new Int_t[fNPids];
299   
300   for (i = 0; i < fNPids-1;i++)//find a slot
301    {
302      if ( fPidProb[i] > prob)
303       {
304         if (GetDebug()>9) Info("SetPID","Copying entry %d",i);
305         aPidProbNew[i] = fPidProb[i];
306         aPidsNew[i] = fPids[i];
307         totprob+=fPidProb[i];
308       }
309      else break;
310    }
311
312   if (GetDebug() > 9) Info("SetPID","Setting new PID on entry %d",i);
313   aPidProbNew[i] = prob;
314   aPidsNew[i] = pdg;
315   totprob+=prob;
316   
317
318   for (Int_t j = fNPids-1; j > i ;j--)//copy rest of old arays 
319    {
320      if (GetDebug() > 9) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
321      aPidProbNew[j] = fPidProb[j-1];
322      aPidsNew[j] = fPids[j-1];
323      totprob+=fPidProb[j-1];
324    }
325
326   delete [] fPidProb;
327   delete [] fPids;
328   
329   fPidProb = aPidProbNew;
330   fPids = aPidsNew;
331   
332   fPdgIdx = GetPidSlot(currentpid);
333   if (fPdgIdx == -1) fPdgIdx = 0;
334   
335   if (totprob > (1.0+0.000001))//space for numerical error
336    {
337      Warning("SetId","Total probability is greater than unity (%f)!!!",totprob);
338      Print();
339    }
340 }
341 //______________________________________________________________________________
342
343 Float_t AliAODParticle::GetPIDprobability(Int_t pdg) const
344 {
345 //Returns probability that this particle is the type of pdg
346   Int_t idx = GetPidSlot(pdg);
347   if (idx < 0) return 0.0;//such pid was not specified for this particle
348   return fPidProb[idx];
349 }
350 //______________________________________________________________________________
351
352 const Char_t* AliAODParticle::GetName() const 
353 {
354   //returns name of this particle 
355    static char def[4] = "XXX";
356    const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
357    if (ap) return ap->GetName();
358    else    return def;
359 }
360 //______________________________________________________________________________
361
362 Int_t AliAODParticle::GetPidSlot(Int_t pdg) const
363 {
364  //returns position of the given PID in fPids (and fPidProb) array.
365  if (fPids == 0x0) return -1;
366  for (Int_t i = 0; i< fNPids; i++)
367   {
368    if (fPids[i] == pdg) return i;
369   }
370  return -1;
371 }
372 //______________________________________________________________________________
373
374 Int_t AliAODParticle::GetNthPid(Int_t idx) const
375 {
376   //returns PID sitting on slot idx in fPids
377   if ( (idx < 0) || (idx >= fNPids) )
378    {
379      Error("GetNthPid","Out Of Bounds");
380      return 0;
381    }
382   return fPids[idx];
383 }
384 //______________________________________________________________________________
385
386 Float_t AliAODParticle::GetNthPidProb(Int_t idx) const
387 {
388   //returns PID sitting on slot idx in fPidProb
389   if ( (idx < 0) || (idx >= fNPids) )
390    {
391      Error("GetNthPid","Out Of Bounds");
392      return 0;
393    }
394   return fPidProb[idx];
395 }
396 //______________________________________________________________________________
397
398 void AliAODParticle::Print() const
399 {
400 //prints information about particle
401   printf("____________________________________________________\n");
402   printf("Idx: %d  PID: %d  Name: ",fIdxInEvent,GetPdgCode());
403   TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
404   if (pdgp)
405    {
406      printf("%s Mass: %f\n",pdgp->GetName(),pdgp->Mass());
407    }
408   else
409    {
410      printf("Not known\n");
411    }
412   
413   printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f\nVx: %+f Vy: %+f Vz: %+f T: %+f\n",
414           Px(),Py(),Pz(),E(),GetCalcMass(),Vx(),Vy(),Vz(),T());
415
416   for (Int_t i = 0; i < fNPids; i++)
417    {
418      printf("# %d  PID: %d  Probability %f name ",i,fPids[i],fPidProb[i]);
419      const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
420      if (ap)
421       {
422         printf("%s Mass %f\n",ap->GetName(),ap->Mass());
423       }
424      else
425       {
426         printf("Not known\n");
427       }
428    }
429 }
430
431 //______________________________________________________________________________
432
433 //void AliAODParticle::Streamer(TBuffer &b)
434 //{
435 //     // Stream all objects in the array to or from the I/O buffer.
436 //   UInt_t R__s, R__c;
437 //   Int_t i;
438 //   if (b.IsReading()) 
439 //    {
440 //      delete [] fPids;
441 //      delete [] fPidProb;
442 //      
443 //      Version_t v = b.ReadVersion(&R__s, &R__c);
444 //      if (v == 1)
445 //       {
446 //         AliAODParticle::Class()->ReadBuffer(b, this);
447 //      }
448 //      else
449 //       {
450 //        TObject::Streamer(b);
451 //       b >> fPdgIdx;
452 //       b >> fIdxInEvent;
453 //       
454 //       b >> fNPids;
455 //       Int_t* fPids = new Int_t[fNPids];
456 //        Float_t* fPidProb = new Float_t[fNPids];
457 //        for (i = 0;i<fNPids;i++) 
458 //         {
459 //           b >> fPids[i];
460 //         }
461 //        for (i = 0;i<fNPids;i++) 
462 //        {
463 //          b >> fPidProb[i];
464 //         }
465 //        b >> fCalcMass;
466 //
467 //        b >> fPx;
468 //       b >> fPy;
469 //       b >> fPz;
470 //        b >> fE;
471 //
472 //       b >> fVx;
473 //        b >> fVy;
474 //        b >> fVz;
475 //       b >> fVt;
476 //       Info("Streamer","Read data");
477 //        Print();
478 //       }
479 //
480 //      b.CheckByteCount(R__s, R__c,AliAODParticle::IsA());
481 //    } 
482 //  else 
483 //   {
484 //     R__c = b.WriteVersion(AliAODParticle::IsA(), kTRUE);
485 //     TObject::Streamer(b);
486 //     Info("Streamer","Read data");
487 //     Print();
488 //
489 //     b << fPdgIdx;
490 //     b << fIdxInEvent;
491 //     b << fNPids;
492 //     for (i = 0;i<fNPids;i++) 
493 //      {
494 //        b << fPids[i];
495 //      }
496 //      {
497 //      {
498 //     for (i = 0;i<fNPids;i++) 
499 //     {
500 //        b << fPidProb[i];
501 //      }
502 //     b << fCalcMass;
503 //
504 //     b << fPx;
505 //     b << fPy;
506 //     b << fPz;
507 //     b << fE;
508 //
509 //     b << fVx;
510 //     b << fVy;
511 //     b << fVz;
512 //     b << fVt;
513 //
514 //    b.SetByteCount(R__c, kTRUE);
515 //   }
516 //}