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