ANALYSIS/AliAODParticle.cxx

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