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