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