JETAN/AliJetParticle.cxx

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