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