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