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