+//__________________________________________________________
+///////////////////////////////////////////////////////////////////
+//
+// class AliHBTEvent
+//
+// This class is container for paticles coming from one event
+//
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////
+
+
#include "AliHBTEvent.h"
#include "AliHBTParticle.h"
-//_________________________________________________________________________
-///////////////////////////////////////////////////////////////////////////
-// //
-// class AliHBTEvent //
-// //
-// This class stores HBT perticles for one event //
-// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html //
-// //
-///////////////////////////////////////////////////////////////////////////
ClassImp(AliHBTEvent)
-const UInt_t AliHBTEvent::fgkInitEventSize = 10000;
-
+const UInt_t AliHBTEvent::fgkInitEventSize = 100;
/**************************************************************************/
-AliHBTEvent::AliHBTEvent():
- fSize(fgkInitEventSize),
- fParticles(new AliHBTParticle* [fSize]),
- fNParticles(0),
- fOwner(kTRUE)
+AliHBTEvent::AliHBTEvent()
{
- //default constructor
+ if(fgkInitEventSize<1)
+ {
+ Fatal("AliHBTEvent::AliHBTEvent()",
+ "fgkInitEventSize has a stiupid value (%d). Change it to positive number and recompile",
+ fgkInitEventSize);
+
+ }
+ fSize=fgkInitEventSize;
+ fParticles = new AliHBTParticle* [fSize];
+ fNParticles = 0;
+ fOwner = kTRUE;
}
/**************************************************************************/
AliHBTEvent::~AliHBTEvent()
{
- //destructor
this->Reset();//delete all particles
if(fParticles)
{
/**************************************************************************/
void AliHBTEvent::Reset()
{
- //deletes all particles from the event
+ //deletes all particles from the event
if(fParticles && fOwner)
{
for(Int_t i =0; i<fNParticles; i++)
}
fNParticles = 0;
}
-/**************************************************************************/
AliHBTParticle* AliHBTEvent::GetParticleSafely(Int_t n)
{
- //returns nth particle with range check
if( (n<0) || (fNParticles<=n) ) return 0x0;
else return fParticles[n];
void AliHBTEvent:: AddParticle(AliHBTParticle* hbtpart)
{
- //Adds new perticle to the event
+ //Adds new perticle to the event
if ( fNParticles+1 >= fSize) Expand(); //if there is no space in array, expand it
fParticles[fNParticles++] = hbtpart; //add a pointer
}
/**************************************************************************/
-void AliHBTEvent::AddParticle(TParticle* part)
+void AliHBTEvent::AddParticle(TParticle* part, Int_t idx)
{
- //Adds TParticle to event
- AddParticle( new AliHBTParticle(*part) );
+ AddParticle( new AliHBTParticle(*part,idx) );
}
/**************************************************************************/
void AliHBTEvent::
-AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+AddParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time)
{
- //adds particle to event
- AddParticle(new AliHBTParticle(pdg,px,py,pz,etot,vx,vy,vz,time) );
+ AddParticle(new AliHBTParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time) );
}
/**************************************************************************/
}
+
#ifndef ALIHBTEvent_H
#define ALIHBTEvent_H
-//_________________________________________________________________________
-///////////////////////////////////////////////////////////////////////////
-// //
-// class AliHBTEvent //
-// //
-// This class stores HBT perticles for one event //
-// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html //
-// //
-///////////////////////////////////////////////////////////////////////////
+//__________________________________________________________
+///////////////////////////////////////////////////////////////////
+//
+// class AliHBTEvent
+//
+// This class is container for paticles coming from one event
+//
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////
#include <TObject.h>
class AliHBTParticle;
class TParticle;
-
class AliHBTEvent: public TObject
{
public:
AliHBTEvent();
virtual ~AliHBTEvent();
-
+ const static UInt_t fgkInitEventSize; //initial number of the array
+ //if expanded, this size is used also
AliHBTParticle* GetParticle(Int_t n); //gets particle
AliHBTParticle* GetParticleSafely(Int_t n); //gets particle with index check
- void AddParticle(AliHBTParticle* hbtpart); //adds particle to the event
- void AddParticle(TParticle* part); //adds particle to the event
- void AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ void AddParticle(AliHBTParticle*); //adds particle to the event
+ void AddParticle(TParticle*, Int_t idx); //adds particle to the event
+ void AddParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time);
Int_t GetNumberOfParticles() const;
void Reset(); //deletes all entries
void SetOwner(Bool_t owns = kTRUE){ fOwner = owns; }
- Bool_t IsOwner() const {return fOwner;}
-
+ Bool_t IsOwner() {return fOwner;}
protected:
- Int_t fSize; //!current size of the array
AliHBTParticle ** fParticles; //!array of pointers to the particles
Int_t fNParticles; //!number of particles in Event
+ Int_t fSize; //!current size of the array
Bool_t fOwner; //flag if that event owns the
void Expand(); //expands the array if necessary
-
private:
- const static UInt_t fgkInitEventSize; //initial number of the array
- //if expanded, this size is used also
+
public:
ClassDef(AliHBTEvent,1)
};
//Simplified TParticle class
-
-
#include "AliHBTParticle.h"
-
#include <TParticle.h>
ClassImp(AliHBTParticle)
+Int_t AliHBTParticle::fgDebug = 0;
//______________________________________________________________________________
AliHBTParticle::AliHBTParticle():
- fPdgCode(0), fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
+ fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
{//empty particle
}
+//______________________________________________________________________________
-
+AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
+ Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time):
+ fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(0),
+ fPx(px), fPy(py),fPz(pz),fE(etot),
+ fVx(vx), fVy(vy),fVz(vz),fVt(time)
+{
+//mormal constructor
+ SetPdgCode(pdg);
+ if (GetPDG()) {
+ fCalcMass = GetPDG()->Mass();
+ } else {
+ Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
+ if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
+ else fCalcMass = -TMath::Sqrt(-a2);
+ }
+}
//______________________________________________________________________________
-AliHBTParticle::AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
- Double_t vx, Double_t vy, Double_t vz, Double_t time):
- fPdgCode(pdg),
+
+AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx,
+ Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time):
+ fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(0),
fPx(px), fPy(py),fPz(pz),fE(etot),
fVx(vx), fVy(vy),fVz(vz),fVt(time)
{
//mormal constructor
-
+ SetPdgCode(pdg,prob);
if (GetPDG()) {
fCalcMass = GetPDG()->Mass();
} else {
else fCalcMass = -TMath::Sqrt(-a2);
}
}
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
+ fPdgIdx(in.fPdgIdx), fIdxInEvent(in.fIdxInEvent),
+ fNPids(in.fNPids),fPids(new Int_t[fNPids]),fPidProb(new Float_t[fNPids]),
+ fCalcMass(in.GetCalcMass()),
+ fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()),
+ fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T())
+{
+ //Copy constructor
+ for(Int_t i = 0; i<fNPids; i++)
+ {
+ fPids[i] = in.fPids[i];
+ fPidProb[i] = in.fPidProb[i];
+ }
+}
//______________________________________________________________________________
-AliHBTParticle::AliHBTParticle(const TParticle &p):
- fPdgCode(p.GetPdgCode()),fCalcMass(p.GetCalcMass()),
+AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
+ fPdgIdx(0), fIdxInEvent(idx),
+ fNPids(0),fPids(0x0),fPidProb(0x0),
+ fCalcMass(p.GetCalcMass()),
fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()),
fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
{
//all copied in the initialization
-
+ SetPdgCode(p.GetPdgCode());
+}
+//______________________________________________________________________________
+
+void AliHBTParticle::SetPdgCode(Int_t pdg,Float_t prob)
+{
+ SetPIDprobability(pdg,prob);
+ fPdgIdx = GetPidSlot(pdg);
}
//______________________________________________________________________________
+void AliHBTParticle::SetPIDprobability(Int_t pdg, Float_t prob)
+{
+//Sets another pdg code and corresponding probabilty
+//Ids are set in decreasing order
+//Check if total prbaility is not ivercoming unity is performed
+//in case, warning is printed
+ if (fgDebug) Info("SetPIDprobability","Setting PID %d prob %f",pdg,prob);
+
+ Float_t totprob = 0.0;//sums up probabilities
+ Int_t idx = GetPidSlot(pdg);
+ Int_t i;
+ if (idx > -1)
+ {
+ fPidProb[idx] = prob;
+ for (i = 0; i < fNPids;i++) totprob+=fPidProb[i];
+ if (totprob > 1.0)
+ {
+ Warning("SetPIDprobability","Total probability greater than UNITY: %f",totprob);
+ }
+ if (fgDebug)
+ {
+ Info("SetPIDprobability","Current Total probability: %f",totprob);
+ }
+ return;
+ }
+
+ Int_t currentpid = GetPdgCode();
+ fNPids++;
+ Float_t* aPidProbNew = new Float_t[fNPids];
+ Int_t* aPidsNew = new Int_t[fNPids];
+
+ for (i = 0; i < fNPids-1;i++)//find a slot
+ {
+ if ( fPidProb[i] > prob)
+ {
+ if (fgDebug>4) Info("SetPID","Copying entry %d",i);
+ aPidProbNew[i] = fPidProb[i];
+ aPidsNew[i] = fPids[i];
+ totprob+=fPidProb[i];
+ }
+ else break;
+ }
+
+ if (fgDebug>4) Info("SetPID","Setting new PID on entry %d",i);
+ aPidProbNew[i] = prob;
+ aPidsNew[i] = pdg;
+ totprob+=prob;
+
+
+ for (Int_t j = fNPids-1; j > i ;i--)//copy rest of old araays
+ {
+ if (fgDebug>4) Info("SetPID","Copying from old entry %d to new entry %d",j-1,j);
+ aPidProbNew[j] = fPidProb[j-1];
+ aPidsNew[j] = fPids[j-1];
+ totprob+=fPidProb[j-1];
+ }
+
+ if (totprob > 1.0)
+ {
+ Warning("SetId","Total probability is greater than 1 !!!");
+ Print();
+ }
+ delete [] fPidProb;
+ delete [] fPids;
+
+ fPidProb = aPidProbNew;
+ fPids = aPidsNew;
+
+ fPdgIdx = GetPidSlot(currentpid);
+ if (fPdgIdx == -1) fPdgIdx = 0;
+}
+//______________________________________________________________________________
+
+Float_t AliHBTParticle::GetPIDprobability(Int_t pdg)
+{
+ Int_t idx = GetPidSlot(pdg);
+ if (idx < 0) return 0.0;//such pid was not specified for this particle
+ return fPidProb[idx];
+}
+//______________________________________________________________________________
+
const Char_t* AliHBTParticle::GetName() const
{
static char def[4] = "XXX";
- const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPdgCode);
+ const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
if (ap) return ap->GetName();
else return def;
}
//______________________________________________________________________________
+Int_t AliHBTParticle::GetPidSlot(Int_t pdg) const
+{
+ //returns position of the given PID in fPids (and fPidProb) array.
+ if (fPids == 0x0) return -1;
+ for (Int_t i = 0; i< fNPids; i++)
+ {
+ if (fPids[i] == pdg) return i;
+ }
+ return -1;
+}
+//______________________________________________________________________________
+void AliHBTParticle::Print() const
+{
+//prints information about particle
+ printf("____________________________________________________\n");
+ printf("Idx: %d PID: %d Name",fIdxInEvent,GetPdgCode());
+ TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(GetPdgCode());
+ if (pdgp)
+ {
+ printf("%s Mass %f\n",pdgp->GetName(),pdgp->Mass());
+ }
+ else
+ {
+ printf("Not known\n");
+ }
+
+ printf("Px: %+f Py: %+f Pz: %+f E: %+f Calculated Mass: %f Vx: %+f Vy: %+f Vz: %+f T: %+f",
+ Px(),Py(),Pz(),Energy(),GetCalcMass(),Vx(),Vy(),Vz(),T());
+ for (Int_t i = 0; i < fNPids; i++)
+ {
+ printf("# %d PID: %d Probability %f name",i,fPids[i],fPidProb[i]);
+ const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPids[i]);
+ if (ap)
+ {
+ printf("%s Mass %f\n",ap->GetName(),ap->Mass());
+ }
+ else
+ {
+ printf("Not known\n");
+ }
+ }
+}
public:
// ****** constructors and destructor
AliHBTParticle();
+ AliHBTParticle(const AliHBTParticle& in);
+
+ AliHBTParticle(Int_t pdg, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time);
- AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx, Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time);
- AliHBTParticle(const TParticle& p);
+ AliHBTParticle(const TParticle& p,Int_t idx);
virtual ~AliHBTParticle(){};
-
- Int_t GetPdgCode () const { return fPdgCode; }
+ void SetPIDprobability(Int_t pdg, Float_t prob = 1.0);
+ Float_t GetPIDprobability(Int_t pdg);
+
+ Int_t GetPdgCode () const { return (fPids)?fPids[fPdgIdx]:0;}
+ Int_t GetPid () const { return GetPdgCode();}
+ Float_t GetPidProb () const { return (fPidProb)?fPidProb[fPdgIdx]:0;}
+
+ Int_t GetUID () const { return fIdxInEvent;}
+
+ void SetPdgCode(Int_t pdg, Float_t prob = 1.0);
Double_t GetCalcMass () const { return fCalcMass; }
- Double_t GetMass () { return GetPDG()->Mass();}
+ Double_t GetMass () { return (GetPDG())?GetPDG()->Mass():-1.;}
- TParticlePDG* GetPDG (){return TDatabasePDG::Instance()->GetParticle(fPdgCode);}
+ TParticlePDG* GetPDG (){return TDatabasePDG::Instance()->GetParticle(GetPdgCode());}
Int_t Beauty () { return GetPDG()->Beauty(); }
Int_t Charm () { return GetPDG()->Charm(); }
void SetProductionVertex(const TLorentzVector& v)
{SetProductionVertex(v.X(),v.Y(),v.Z(),v.T());}
const Char_t* GetName() const;
+ void Print() const;
+
+ static void SetDebug(Int_t dbg=1){fgDebug=dbg;}
+ static Int_t GetDebug(){return fgDebug;}
+ static Int_t fgDebug; //debug printout level
protected:
-
- Int_t fPdgCode; // PDG code of the particle
+ Int_t GetPidSlot(Int_t pdg) const;//returns position of the given PID in fPids (and fPidProb) array.
+
+private:
+ Char_t fPdgIdx; // index of PDG code of the particle in fPids
+ Int_t fIdxInEvent; // index of a particle: the same particle can appear in the event
+ // many times with different pid's. Idx allows to check that they are the same particles
+ Char_t fNPids; // number of non-zero proboble Pids
+ Int_t *fPids; //[fNPids]
+ Float_t *fPidProb; //[fNPids]
Double_t fCalcMass; // Calculated mass
Double_t fPx; // x component of momentum
Double_t fVz; // z of production vertex
Double_t fVt; // t of production vertex
-public:
ClassDef(AliHBTParticle,1) // TParticle vertex particle information
};
if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
//if not take next partilce
- AliHBTParticle* part = new AliHBTParticle(*p);
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
Double_t y= iotrack->GetY();
Double_t z= iotrack->GetZ();
- AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), px, py , pz, tEtot, x, y, z, 0.);
+ AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i , px, py , pz, tEtot, x, y, z, 0.);
if(Pass(track)) { delete track;continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
//if not take next partilce
- AliHBTParticle* part = new AliHBTParticle(*p);
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
Double_t mass = p->GetMass();
Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
- AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
if(Pass(track))//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
{
if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
//if not take next partilce
- AliHBTParticle* part = new AliHBTParticle(*p);
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
particles->AddParticle(totalNevents,part);//put track and particle on the run
Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
- AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+ AliHBTParticle* part = new AliHBTParticle(gt.code, i, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
if(Pass(part)) continue;
Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
- AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ AliHBTParticle* track = new AliHBTParticle(gt.code, i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
if(Pass(track)) continue;
particles->AddParticle(currentEvent,part);
if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
//if not take next partilce
- AliHBTParticle* part = new AliHBTParticle(*p);
+ AliHBTParticle* part = new AliHBTParticle(*p,i);
if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
Double_t mass = p->GetMass();
Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
- AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
if(Pass(track))//check if meets all criteria of any of our cuts
//if it does not delete it and take next good track
{
}
/**************************************************************************/
-void AliHBTRun::AddParticle(Int_t event, TParticle* part)
+void AliHBTRun::AddParticle(Int_t event, TParticle* part, Int_t idx)
{
//if there is no event of this number, crate it and add to the collection
if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
- GetEvent(event)->AddParticle(part);
+ GetEvent(event)->AddParticle(part,idx);
}
/**************************************************************************/
-void AliHBTRun::AddParticle(Int_t event, Int_t pdg,
+void AliHBTRun::AddParticle(Int_t event, Int_t pdg, Int_t idx,
Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time)
{
//if there is no event of this number, crate it and add to the collection
if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
- GetEvent(event)->AddParticle(pdg,px,py,pz,etot,vx,vy,vz,time);
+ GetEvent(event)->AddParticle(pdg,idx,px,py,pz,etot,vx,vy,vz,time);
}
/**************************************************************************/
virtual ~AliHBTRun();
void AddParticle(Int_t event, AliHBTParticle* part); //inerface to AliHBTEvent::AddParticle(AliHBTParticle*)
- void AddParticle(Int_t event, TParticle* part);//inerface to AliHBTEvent::AddParticle(TParticle*)
+ void AddParticle(Int_t event, TParticle* part, Int_t idx);//inerface to AliHBTEvent::AddParticle(TParticle*)
//inerface to AliHBTEvent::AddParticle(Int_t.Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t)
- void AddParticle(Int_t event, Int_t pdg,
+ void AddParticle(Int_t event, Int_t pdg, Int_t idx,
Double_t px, Double_t py, Double_t pz, Double_t etot,
Double_t vx, Double_t vy, Double_t vz, Double_t time);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff