HBTAN/AliHBTParticle.cxx

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