Introducing hyperons in ESD (Yu.Belikov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Jul 2003 14:28:14 +0000 (14:28 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Jul 2003 14:28:14 +0000 (14:28 +0000)
21 files changed:
ITS/AliCascadeVertex.cxx
ITS/AliCascadeVertex.h
ITS/AliCascadeVertexer.cxx
ITS/AliCascadeVertexer.h
ITS/AliV0vertex.cxx
ITS/AliV0vertex.h
ITS/AliV0vertexer.cxx
ITS/AliV0vertexer.h
STEER/AliESD.cxx
STEER/AliESD.h
STEER/AliESDcascade.cxx [new file with mode: 0644]
STEER/AliESDcascade.h [new file with mode: 0644]
STEER/AliESDtest.C
STEER/AliESDtrack.cxx
STEER/AliESDtrack.h
STEER/AliESDv0.cxx [new file with mode: 0644]
STEER/AliESDv0.h [new file with mode: 0644]
STEER/AliESDv0Analysis.C [new file with mode: 0644]
STEER/AliKalmanTrack.h
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

index cdb0508..3671dfa 100644 (file)
@@ -18,9 +18,7 @@
 //
 //    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
 //-------------------------------------------------------------------------
-#include <Riostream.h>
 #include <TMath.h>
-#include <TPDGCode.h>
 
 #include "AliCascadeVertex.h"
 #include "AliITStrackV2.h"
 
 ClassImp(AliCascadeVertex)
 
-AliCascadeVertex::AliCascadeVertex() : TObject() {
-  //--------------------------------------------------------------------
-  // Default constructor  (Xi-)
-  //--------------------------------------------------------------------
-  fPdgCode=kXiMinus;
-  fEffMass=1.32131;
-  fChi2=1.e+33;
-  fPos[0]=fPos[1]=fPos[2]=0.;
-  fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.;
-}
-
-
 
 inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){
   // determinant 2x2
@@ -55,15 +41,14 @@ inline Double_t det (Double_t a00,Double_t a01,Double_t a02,
 }
 
 
-
 AliCascadeVertex::AliCascadeVertex(const AliV0vertex &v,const AliITStrackV2 &t) {
   //--------------------------------------------------------------------
   // Main constructor
   //--------------------------------------------------------------------
   fPdgCode=kXiMinus;
 
-  fV0lab[0]=v.GetNlabel(); fV0lab[1]=v.GetPlabel();
-  fBachLab=t.GetLabel(); 
+  fV0idx[0]=v.GetNindex(); fV0idx[1]=v.GetPindex();
+  fBachIdx=t.GetLabel(); 
 
   //Trivial estimation of the vertex parameters
   Double_t pt, phi, x, par[5];
@@ -122,187 +107,4 @@ AliCascadeVertex::AliCascadeVertex(const AliV0vertex &v,const AliITStrackV2 &t)
 
 }
 
-/*
-Double_t AliCascadeVertex::ChangeMassHypothesis(Double_t &v0q, Int_t code) {
-  //--------------------------------------------------------------------
-  // This function changes the mass hypothesis for this cascade
-  // and returns the "kinematical quality" of this hypothesis
-  // together with the "quality" of associated V0 (argument v0q) 
-  //--------------------------------------------------------------------
-  Double_t nmass=0.13957, pmass=0.93827, des0=0.9437-0.1723; 
-  Double_t bmass=0.13957, mass =1.3213,  des =1.1243-0.1970;
-
-  fPdgCode=code;
-
-  switch (code) {
-  case 213: 
-       bmass=0.93827; 
-       break;
-  case kXiMinus:
-       break;
-  case kXiPlusBar:
-       nmass=0.93827; pmass=0.13957; des0=-des0; 
-       des=-des;
-       break;
-  case kOmegaMinus: 
-       bmass=0.49368; mass=1.67245; des=1.1355-0.5369;
-       break;
-  case kOmegaPlusBar: 
-       nmass=0.93827; pmass=0.13957; des0=-des0; 
-       bmass=0.49368; mass=1.67245; des=0.5369-1.1355;
-       break;
-  default:
-       cerr<<"AliCascadeVertex::ChangeMassHypothesis: ";
-       cerr<<"Invalide PDG code !  Assuming XiMinus's...\n";
-       fPdgCode=kXiMinus;
-    break;
-  }
-
-  Double_t pxn=fV0mom[0][0], pyn=fV0mom[0][1], pzn=fV0mom[0][2];
-  Double_t pxp=fV0mom[1][0], pyp=fV0mom[1][1], pzp=fV0mom[1][2];
-  Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
-  Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
-  Double_t px0=pxn+pxp, py0=pyn+pyp, pz0=pzn+pzp;
-  Double_t p0=TMath::Sqrt(px0*px0 + py0*py0 + pz0*pz0);
-
-  Double_t gamma0=(en+ep)/1.11568, betagamma0=p0/1.11568;
-  Double_t pln=(pxn*px0 + pyn*py0 + pzn*pz0)/p0;
-  Double_t plp=(pxp*px0 + pyp*py0 + pzp*pz0)/p0;
-  Double_t plps=gamma0*plp - betagamma0*ep;
-
-  Double_t diff0=2*gamma0*plps + betagamma0*des0;
-
-
-  v0q=plp-pln-diff0;
-
-
-  Double_t pxb=fBachMom[0], pyb=fBachMom[1], pzb=fBachMom[2]; 
-
-  Double_t e0=TMath::Sqrt(1.11568*1.11568 + p0*p0);
-  Double_t eb=TMath::Sqrt(bmass*bmass + pxb*pxb + pyb*pyb + pzb*pzb);
-  Double_t pxl=px0+pxb, pyl=py0+pyb, pzl=pz0+pzb;
-  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
-  
-  fEffMass=TMath::Sqrt((e0+eb)*(e0+eb) - pl*pl);
-
-  Double_t gamma=(e0+eb)/mass, betagamma=pl/mass;
-  Double_t pl0=(px0*pxl + py0*pyl + pz0*pzl)/pl;
-  Double_t plb=(pxb*pxl + pyb*pyl + pzb*pzl)/pl;
-  Double_t pl0s=gamma*pl0 - betagamma*e0;
-
-  Double_t diff=2*gamma*pl0s + betagamma*des;
-
-  return (pl0-plb-diff);
-}
-*/
-
-Double_t AliCascadeVertex::ChangeMassHypothesis(Double_t &v0q, Int_t code) {
-  //--------------------------------------------------------------------
-  // This function changes the mass hypothesis for this cascade
-  // and returns the "kinematical quality" of this hypothesis
-  // together with the "quality" of associated V0 (argument v0q) 
-  //--------------------------------------------------------------------
-  Double_t nmass=0.13957, pmass=0.93827, ps0=0.101; 
-  Double_t bmass=0.13957, mass =1.3213,  ps =0.139;
-
-  fPdgCode=code;
-
-  switch (code) {
-  case 213: 
-       bmass=0.93827; 
-       break;
-  case kXiMinus:
-       break;
-  case kXiPlusBar:
-       nmass=0.93827; pmass=0.13957; 
-       break;
-  case kOmegaMinus: 
-       bmass=0.49368; mass=1.67245; ps=0.211;
-       break;
-  case kOmegaPlusBar: 
-       nmass=0.93827; pmass=0.13957; 
-       bmass=0.49368; mass=1.67245; ps=0.211;
-       break;
-  default:
-       cerr<<"AliCascadeVertex::ChangeMassHypothesis: ";
-       cerr<<"Invalide PDG code !  Assuming XiMinus's...\n";
-       fPdgCode=kXiMinus;
-    break;
-  }
-
-  Double_t pxn=fV0mom[0][0], pyn=fV0mom[0][1], pzn=fV0mom[0][2];
-  Double_t pxp=fV0mom[1][0], pyp=fV0mom[1][1], pzp=fV0mom[1][2];
-  Double_t px0=pxn+pxp, py0=pyn+pyp, pz0=pzn+pzp;
-  Double_t p0=TMath::Sqrt(px0*px0 + py0*py0 + pz0*pz0);
-
-  Double_t e0=TMath::Sqrt(1.11568*1.11568 + p0*p0);
-  Double_t beta0=p0/e0;
-  Double_t pln=(pxn*px0 + pyn*py0 + pzn*pz0)/p0;
-  Double_t plp=(pxp*px0 + pyp*py0 + pzp*pz0)/p0;
-  Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
-
-  Double_t a=(plp-pln)/(plp+pln);
-  a -= (pmass*pmass-nmass*nmass)/(1.11568*1.11568);
-  a = 0.25*beta0*beta0*1.11568*1.11568*a*a + pt2;
-
-
-  v0q=a - ps0*ps0;
-
-
-  Double_t pxb=fBachMom[0], pyb=fBachMom[1], pzb=fBachMom[2]; 
-
-  Double_t eb=TMath::Sqrt(bmass*bmass + pxb*pxb + pyb*pyb + pzb*pzb);
-  Double_t pxl=px0+pxb, pyl=py0+pyb, pzl=pz0+pzb;
-  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
-  
-  fEffMass=TMath::Sqrt((e0+eb)*(e0+eb) - pl*pl);
-
-  Double_t beta=pl/(e0+eb);
-  Double_t pl0=(px0*pxl + py0*pyl + pz0*pzl)/pl;
-  Double_t plb=(pxb*pxl + pyb*pyl + pzb*pzl)/pl;
-  pt2=p0*p0 - pl0*pl0;
-
-  a=(pl0-plb)/(pl0+plb);
-  a -= (1.11568*1.11568-bmass*bmass)/(mass*mass);
-  a = 0.25*beta*beta*mass*mass*a*a + pt2;
-
-  return (a - ps*ps);
-}
-
-void 
-AliCascadeVertex::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-  //--------------------------------------------------------------------
-  // This function returns the cascade momentum (global)
-  //--------------------------------------------------------------------
-  px=fV0mom[0][0]+fV0mom[1][0]+fBachMom[0]; 
-  py=fV0mom[0][1]+fV0mom[1][1]+fBachMom[1]; 
-  pz=fV0mom[0][2]+fV0mom[1][2]+fBachMom[2]; 
-}
-
-void AliCascadeVertex::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
-  //--------------------------------------------------------------------
-  // This function returns cascade position (global)
-  //--------------------------------------------------------------------
-  x=fPos[0]; 
-  y=fPos[1]; 
-  z=fPos[2]; 
-}
-
-Double_t AliCascadeVertex::GetD(Double_t x0, Double_t y0, Double_t z0) const {
-  //--------------------------------------------------------------------
-  // This function returns the cascade impact parameter
-  //--------------------------------------------------------------------
-
-  Double_t x=fPos[0],y=fPos[1],z=fPos[2];
-  Double_t px=fV0mom[0][0]+fV0mom[1][0]+fBachMom[0];
-  Double_t py=fV0mom[0][1]+fV0mom[1][1]+fBachMom[1];
-  Double_t pz=fV0mom[0][2]+fV0mom[1][2]+fBachMom[2];
-
-  Double_t dx=(y0-y)*pz - (z0-z)*py; 
-  Double_t dy=(x0-x)*pz - (z0-z)*px;
-  Double_t dz=(x0-x)*py - (y0-y)*px;
-  Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
-
-  return d;
-}
 
index 32eaa9d..cc3e549 100644 (file)
@@ -10,8 +10,7 @@
 //    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
 //-------------------------------------------------------------------------
 
-#include <TObject.h>
-#include "AliPDG.h"
+#include "AliESDcascade.h"
 
 class AliITStrackV2;
 class AliV0vertex;
@@ -21,48 +20,11 @@ class AliV0vertex;
 #define kOmegaMinus    3334
 #define kOmegaPlusBar -3334
 
-class AliCascadeVertex : public TObject {
+class AliCascadeVertex : public AliESDcascade {
 public:
-  AliCascadeVertex();
+  AliCascadeVertex():AliESDcascade(){;}
   AliCascadeVertex(const AliV0vertex &vtx, const AliITStrackV2 &trk);
 
-  Double_t ChangeMassHypothesis(Double_t &v0q, Int_t code=kXiMinus); 
-
-  Int_t GetPdgCode() const {return fPdgCode;}
-  Double_t GetEffMass() const {return fEffMass;}
-  Double_t GetChi2() const {return fChi2;}
-  void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
-  void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const;
-  Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
-
-  void GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-     px=fV0mom[0][0]; py=fV0mom[0][1]; pz=fV0mom[0][2];
-  }
-  Int_t GetNlabel() const {return fV0lab[0];}
-  void GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-     px=fV0mom[1][0]; py=fV0mom[1][1]; pz=fV0mom[1][2];
-  }
-  Int_t GetPlabel() const {return fV0lab[1];}
-  void GetBPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-     px=fBachMom[0]; py=fBachMom[1]; pz=fBachMom[2];
-  }
-  Int_t GetBlabel() const {return fBachLab;}
-
-private: 
-  Int_t fPdgCode;           // reconstructed cascade type (PDG code)
-  Double_t fEffMass;        // reconstructed cascade effective mass
-  Double_t fChi2;           // chi2 value
-  Double_t fPos[3];         // cascade vertex position (global)
-  Double_t fPosCov[6];      // covariance matrix of the vertex position
-
-  Int_t fV0lab[2];          // labels of the V0 daughter tracks
-  Double_t fV0mom[2][3];    // V0 daughters' momenta (global)
-  Double_t fV0momCov[6];    // covariance matrix of the V0 momentum.
-
-  Int_t fBachLab;           // label of the bachelor track
-  Double_t fBachMom[3];      // bachelor momentum (global)
-  Double_t fBachMomCov[6];   // covariance matrix of the bachelor momentum.
-
   ClassDef(AliCascadeVertex,1)   // reconstructed cascade vertex
 };
 
index 1159193..377d6cd 100644 (file)
 //
 //    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
 //-------------------------------------------------------------------------
-#include <Riostream.h>
 #include <TObjArray.h>
-#include <TPDGCode.h>
 #include <TTree.h>
 
+#include "AliESD.h"
+#include "AliESDv0.h"
+#include "AliESDcascade.h"
 #include "AliCascadeVertex.h"
 #include "AliCascadeVertexer.h"
 #include "AliITStrackV2.h"
 
 ClassImp(AliCascadeVertexer)
 
+Int_t AliCascadeVertexer::V0sTracks2CascadeVertices(AliESD *event) {
+  //--------------------------------------------------------------------
+  // This function reconstructs cascade vertices
+  //      Adapted to the ESD by I.Belikov (Jouri.Belikov@cern.ch)
+  //--------------------------------------------------------------------
+
+   Int_t nV0=(Int_t)event->GetNumberOfV0s();
+   TObjArray vtcs(nV0);
+   Int_t i;
+   for (i=0; i<nV0; i++) {
+       const AliESDv0 *esdV0=event->GetV0(i);
+       vtcs.AddLast(new AliV0vertex(*esdV0));
+   }
+
+
+   Int_t ntr=(Int_t)event->GetNumberOfTracks();
+   TObjArray trks(ntr);
+   for (i=0; i<ntr; i++) {
+       AliESDtrack *esdtr=event->GetTrack(i);
+       AliITStrackV2 *iotrack=new AliITStrackV2(*esdtr);
+       iotrack->PropagateTo(3.,0.0023,65.19); iotrack->PropagateTo(2.5,0.,0.);
+       trks.AddLast(iotrack);
+   }   
+
+   Double_t massLambda=1.11568;
+   Int_t ncasc=0;
+
+   // Looking for the cascades...
+   for (i=0; i<nV0; i++) {
+      AliV0vertex *v=(AliV0vertex*)vtcs.UncheckedAt(i);
+      v->ChangeMassHypothesis(kLambda0); // the v0 must be Lambda 
+      if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue; 
+      if (v->GetD(0,0,0)<fDV0min) continue;
+      for (Int_t j=0; j<ntr; j++) {
+        AliITStrackV2 *b=(AliITStrackV2*)trks.UncheckedAt(j);
+
+         if (TMath::Abs(b->GetD())<fDBachMin) continue;
+         if (b->Get1Pt()<0.) continue;  // bachelor's charge 
+          
+        AliV0vertex v0(*v), *pv0=&v0;
+         AliITStrackV2 bt(*b), *pbt=&bt;
+
+         Double_t dca=PropagateToDCA(pv0,pbt);
+         if (dca > fDCAmax) continue;
+
+         AliCascadeVertex cascade(*pv0,*pbt);
+         if (cascade.GetChi2() > fChi2max) continue;
+
+        Double_t x,y,z; cascade.GetXYZ(x,y,z); 
+         Double_t r2=x*x + y*y; 
+         if (r2 > fRmax*fRmax) continue;   // condition on fiducial zone
+         if (r2 < fRmin*fRmin) continue;
+
+         {
+         Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
+         if (r2 > (x1*x1+y1*y1)) continue;
+         if (z*z > z1*z1) continue;
+         }
+
+        Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
+         Double_t p2=px*px+py*py+pz*pz;
+         Double_t cost=(x*px+y*py+z*pz)/TMath::Sqrt(p2*(r2+z*z));
+
+         if (cost<fCPAmax) continue; //condition on the cascade pointing angle 
+         //cascade.ChangeMassHypothesis(); //default is Xi
+
+         event->AddCascade(&cascade);
+
+         ncasc++;
+
+      }
+   }
+
+   // Looking for the anti-cascades...
+   for (i=0; i<nV0; i++) {
+      AliV0vertex *v=(AliV0vertex*)vtcs.UncheckedAt(i);
+      v->ChangeMassHypothesis(kLambda0Bar); //the v0 must be anti-Lambda 
+      if (TMath::Abs(v->GetEffMass()-massLambda)>fMassWin) continue; 
+      if (v->GetD(0,0,0)<fDV0min) continue;
+      for (Int_t j=0; j<ntr; j++) {
+        AliITStrackV2 *b=(AliITStrackV2*)trks.UncheckedAt(j);
+
+         if (TMath::Abs(b->GetD())<fDBachMin) continue;
+         if (b->Get1Pt()>0.) continue;  // bachelor's charge 
+          
+        AliV0vertex v0(*v), *pv0=&v0;
+         AliITStrackV2 bt(*b), *pbt=&bt;
+
+         Double_t dca=PropagateToDCA(pv0,pbt);
+         if (dca > fDCAmax) continue;
+
+         AliCascadeVertex cascade(*pv0,*pbt);
+         if (cascade.GetChi2() > fChi2max) continue;
+
+        Double_t x,y,z; cascade.GetXYZ(x,y,z); 
+         Double_t r2=x*x + y*y; 
+         if (r2 > fRmax*fRmax) continue;   // condition on fiducial zone
+         if (r2 < fRmin*fRmin) continue;
+
+         {
+         Double_t x1,y1,z1; pv0->GetXYZ(x1,y1,z1);
+         if (r2 > (x1*x1+y1*y1)) continue;
+         if (z*z > z1*z1) continue;
+         }
+
+        Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz);
+         Double_t p2=px*px+py*py+pz*pz;
+         Double_t cost=(x*px+y*py+z*pz)/TMath::Sqrt(p2*(r2+z*z));
+
+         if (cost<fCPAmax) continue; //condition on the cascade pointing angle 
+         //cascade.ChangeMassHypothesis(); //default is Xi
+
+         event->AddCascade(&cascade);
+
+         ncasc++;
+
+      }
+   }
+
+Info("V0sTracks2CascadeVertices","Number of reconstructed cascades: %d",ncasc);
+
+   trks.Delete();
+   vtcs.Delete();
+
+   return 0;
+}
+
 Int_t AliCascadeVertexer::
 V0sTracks2CascadeVertices(TTree *vTree,TTree *tTree, TTree *xTree) {
   //--------------------------------------------------------------------
@@ -160,7 +288,7 @@ V0sTracks2CascadeVertices(TTree *vTree,TTree *tTree, TTree *xTree) {
       }
    }
 
-   cerr<<"Number of reconstructed cascades: "<<ncasc<<endl;
+Info("V0sTracks2CascadeVertices","Number of reconstructed cascades: %d",ncasc);
 
    trks.Delete();
    vtxV0.Delete();
@@ -230,7 +358,7 @@ AliCascadeVertexer::PropagateToDCA(AliV0vertex *v, AliITStrackV2 *t) {
 
   x1=x1*cs1 + y1*sn1;
   if (!t->PropagateTo(x1,0.,0.)) {
-    cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n";
+    Error("PropagateToDCA","Propagation failed !");
     return 1.e+33;
   }  
 
index 3ce87ff..c14b5ff 100644 (file)
@@ -11,6 +11,7 @@
 
 #include "TObject.h"
 
+class AliESD;
 class TTree;
 class AliITStrackV2;
 class AliV0vertex;
@@ -22,6 +23,7 @@ public:
   AliCascadeVertexer(const Double_t cuts[8]);
   void SetCuts(const Double_t cuts[8]);
 
+  Int_t V0sTracks2CascadeVertices(AliESD *event);
   Int_t V0sTracks2CascadeVertices(TTree *v, TTree *t, TTree *x);
   Double_t PropagateToDCA(AliV0vertex *vtx, AliITStrackV2 *trk);
 
index f663e2f..9dd719a 100644 (file)
 //
 //     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
 //-------------------------------------------------------------------------
-#include <Riostream.h>
 #include <TMath.h>
-#include <TPDGCode.h>
 
 #include "AliV0vertex.h"
 #include "AliITStrackV2.h"
 
 ClassImp(AliV0vertex)
 
-AliV0vertex::AliV0vertex() : TObject() {
-  //--------------------------------------------------------------------
-  // Default constructor  (K0s)
-  //--------------------------------------------------------------------
-  fPdgCode=kK0Short;
-  fEffMass=0.497672;
-  fChi2=1.e+33;
-  fPos[0]=fPos[1]=fPos[2]=0.;
-  fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.;
-}
-
 AliV0vertex::AliV0vertex(const AliITStrackV2 &n, const AliITStrackV2 &p) {
   //--------------------------------------------------------------------
   // Main constructor
   //--------------------------------------------------------------------
   fPdgCode=kK0Short;
-  fNlab=n.GetLabel(); fPlab=p.GetLabel();
+  fNidx=n.GetLabel(); fPidx=p.GetLabel(); //indices in the array of ESD tracks
 
   //Trivial estimation of the vertex parameters
   Double_t pt, phi, x, par[5];
@@ -86,128 +73,4 @@ AliV0vertex::AliV0vertex(const AliITStrackV2 &n, const AliITStrackV2 &p) {
 
   fChi2=7.;   
 }
-/*
-Double_t AliV0vertex::ChangeMassHypothesis(Int_t code) {
-  //--------------------------------------------------------------------
-  // This function changes the mass hypothesis for this V0
-  // and returns the "kinematical quality" of this hypothesis 
-  //--------------------------------------------------------------------
-  Double_t nmass=0.13957, pmass=0.13957, mass=0.49767, des=0;
-
-  fPdgCode=code;
-
-  switch (code) {
-  case kLambda0:
-    nmass=0.13957; pmass=0.93827; mass=1.1157; des=0.9437-0.1723; break;
-  case kLambda0Bar:
-    pmass=0.13957; nmass=0.93827; mass=1.1157; des=0.1723-0.9437; break;
-  case kK0Short: 
-    break;
-  default:
-    cerr<<"AliV0vertex::ChangeMassHypothesis: ";
-    cerr<<"invalide PDG code ! Assuming K0s...\n";
-    fPdgCode=kK0Short;
-    break;
-  }
-
-  Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
-  Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
-
-  Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
-  Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
-  Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
-  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
-
-  fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
-
-  Double_t gamma=(en+ep)/mass, betagamma=pl/mass;
-  Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
-  Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
-  Double_t plps=gamma*plp - betagamma*ep;
-
-  Double_t diff=2*gamma*plps + betagamma*des;
-
-  return (plp-pln-diff);
-  
-}
-*/
-Double_t AliV0vertex::ChangeMassHypothesis(Int_t code) {
-  //--------------------------------------------------------------------
-  // This function changes the mass hypothesis for this V0
-  // and returns the "kinematical quality" of this hypothesis 
-  //--------------------------------------------------------------------
-  Double_t nmass=0.13957, pmass=0.13957, mass=0.49767, ps=0.206;
-
-  fPdgCode=code;
-
-  switch (code) {
-  case kLambda0:
-    nmass=0.13957; pmass=0.93827; mass=1.1157; ps=0.101; break;
-  case kLambda0Bar:
-    pmass=0.13957; nmass=0.93827; mass=1.1157; ps=0.101; break;
-  case kK0Short: 
-    break;
-  default:
-    cerr<<"AliV0vertex::ChangeMassHypothesis: ";
-    cerr<<"invalide PDG code ! Assuming K0s...\n";
-    fPdgCode=kK0Short;
-    break;
-  }
-
-  Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
-  Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
-
-  Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
-  Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
-  Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
-  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
-
-  fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
-
-  Double_t beta=pl/(en+ep);
-  Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
-  Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
-
-  Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
-
-  Double_t a=(plp-pln)/(plp+pln);
-  a -= (pmass*pmass-nmass*nmass)/(mass*mass);
-  a = 0.25*beta*beta*mass*mass*a*a + pt2;
-
-  return (a - ps*ps);
-  
-}
-
-void AliV0vertex::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-  //--------------------------------------------------------------------
-  // This function returns V0's momentum (global)
-  //--------------------------------------------------------------------
-  px=fNmom[0]+fPmom[0]; 
-  py=fNmom[1]+fPmom[1]; 
-  pz=fNmom[2]+fPmom[2]; 
-}
-
-void AliV0vertex::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
-  //--------------------------------------------------------------------
-  // This function returns V0's position (global)
-  //--------------------------------------------------------------------
-  x=fPos[0]; 
-  y=fPos[1]; 
-  z=fPos[2]; 
-}
 
-Double_t AliV0vertex::GetD(Double_t x0, Double_t y0, Double_t z0) const {
-  //--------------------------------------------------------------------
-  // This function returns V0's impact parameter
-  //--------------------------------------------------------------------
-  Double_t x=fPos[0],y=fPos[1],z=fPos[2];
-  Double_t px=fNmom[0]+fPmom[0];
-  Double_t py=fNmom[1]+fPmom[1];
-  Double_t pz=fNmom[2]+fPmom[2];
-
-  Double_t dx=(y0-y)*pz - (z0-z)*py; 
-  Double_t dy=(x0-x)*pz - (z0-z)*px;
-  Double_t dz=(x0-x)*py - (y0-y)*px;
-  Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
-  return d;
-}
index b18fca0..9391c6d 100644 (file)
 //    Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch 
 //-------------------------------------------------------------------------
 
-#include <TObject.h>
-#include <TPDGCode.h>
+#include <AliESDv0.h>
 
 class AliITStrackV2;
 
-class AliV0vertex : public TObject {
+class AliV0vertex : public AliESDv0 {
 public:
-  AliV0vertex();
+  AliV0vertex() : AliESDv0() {;}
+  AliV0vertex(const AliESDv0 &v) : AliESDv0(v) {;}
   AliV0vertex(const AliITStrackV2 &neg, const AliITStrackV2 &pos);
 
-  Double_t ChangeMassHypothesis(Int_t code=kK0Short); 
-
-  Int_t GetPdgCode() const {return fPdgCode;}
-  Double_t GetEffMass() const {return fEffMass;}
-  Double_t GetChi2() const {return fChi2;}
-  void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
-  void GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
-  void GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
-  void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const;
-  Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
-  Int_t GetNlabel() const {return fNlab;}
-  Int_t GetPlabel() const {return fPlab;}
-
-private: 
-  Int_t fPdgCode;           // reconstructed V0's type (PDG code)
-  Double_t fEffMass;        // reconstructed V0's effective mass
-  Double_t fChi2;           // V0's chi2 value
-  Double_t fPos[3];         // V0's position (global)
-  Double_t fPosCov[6];      // covariance matrix of the vertex position
-
-  Int_t fNlab;              // label of the negative daughter
-  Double_t fNmom[3];        // momentum of the negative daughter (global)
-  Double_t fNmomCov[6];     // covariance matrix of the negative daughter mom.
-
-  Int_t fPlab;              // label of the positive daughter
-  Double_t fPmom[3];        // momentum of the positive daughter (global)
-  Double_t fPmomCov[6];     // covariance matrix of the positive daughter mom.
-
-  ClassDef(AliV0vertex,1)   // reconstructed V0 vertex
+  ClassDef(AliV0vertex,1) // reconstructed V0 vertex
 };
 
-inline 
-void AliV0vertex::GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-px=fNmom[0]; py=fNmom[1]; pz=fNmom[2];
-}
-
-inline 
-void AliV0vertex::GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
-px=fPmom[0]; py=fPmom[1]; pz=fPmom[2];
-}
-
 #endif
 
 
index 49e2232..4eda554 100644 (file)
 //
 //     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
 //-------------------------------------------------------------------------
-#include <Riostream.h>
-#include <TPDGCode.h>
 #include <TObjArray.h>
 #include <TTree.h>
 
+#include "AliESD.h"
+#include "AliESDtrack.h"
 #include "AliITStrackV2.h"
 #include "AliV0vertex.h"
 #include "AliV0vertexer.h"
 
 ClassImp(AliV0vertexer)
 
+Int_t AliV0vertexer::Tracks2V0vertices(AliESD *event) {
+  //--------------------------------------------------------------------
+  //This function reconstructs V0 vertices
+  //--------------------------------------------------------------------
+
+   Int_t nentr=event->GetNumberOfTracks();
+
+   TObjArray negtrks(nentr/2);
+   TObjArray postrks(nentr/2);
+
+   Int_t nneg=0, npos=0, nvtx=0;
+
+   Int_t i;
+   for (i=0; i<nentr; i++) {
+     AliESDtrack *esd=event->GetTrack(i);
+     Int_t status=esd->GetStatus();
+
+  if ((status&AliESDtrack::kITSrefit)==0)
+     if ((status&AliESDtrack::kITSout)!=0 || (status&AliESDtrack::kITSin)==0)
+  continue;
+
+     AliITStrackV2 *iotrack=new AliITStrackV2(*esd);
+     iotrack->SetLabel(i);  // now it is the index in array of ESD tracks
+     iotrack->PropagateTo(3.,0.0023,65.19); 
+     iotrack->PropagateTo(2.5,0.,0.);
+
+     if (iotrack->Get1Pt() > 0.) {nneg++; negtrks.AddLast(iotrack);}
+     else {npos++; postrks.AddLast(iotrack);}
+   }   
+
+
+   for (i=0; i<nneg; i++) {
+      //if (i%10==0) cerr<<nneg-i<<'\r';
+      AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
+
+      if (TMath::Abs(ntrk->GetD(fX,fY))<fDPmin) continue;
+      if (TMath::Abs(ntrk->GetD(fX,fY))>fRmax) continue;
+
+      for (Int_t k=0; k<npos; k++) {
+         AliITStrackV2 *ptrk=(AliITStrackV2 *)postrks.UncheckedAt(k);
+
+         if (TMath::Abs(ptrk->GetD(fX,fY))<fDPmin) continue;
+         if (TMath::Abs(ptrk->GetD(fX,fY))>fRmax) continue;
+
+         if (TMath::Abs(ntrk->GetD(fX,fY))<fDNmin)
+         if (TMath::Abs(ptrk->GetD(fX,fY))<fDNmin) continue;
+
+
+         AliITStrackV2 nt(*ntrk), pt(*ptrk), *pnt=&nt, *ppt=&pt;
+
+         Double_t dca=PropagateToDCA(pnt,ppt);
+         if (dca > fDCAmax) continue;
+
+         AliV0vertex vertex(*pnt,*ppt);
+         if (vertex.GetChi2() > fChi2max) continue;
+        
+         /*  Think of something better here ! 
+         nt.PropagateToVertex(); if (TMath::Abs(nt.GetZ())<0.04) continue;
+         pt.PropagateToVertex(); if (TMath::Abs(pt.GetZ())<0.04) continue;
+        */
+
+         Double_t x,y,z; vertex.GetXYZ(x,y,z); 
+         Double_t r2=x*x + y*y; 
+         if (r2 > fRmax*fRmax) continue;
+         if (r2 < fRmin*fRmin) continue;
+
+         Double_t px,py,pz; vertex.GetPxPyPz(px,py,pz);
+         Double_t p2=px*px+py*py+pz*pz;
+         Double_t cost=((x-fX)*px + (y-fY)*py + (z-fZ)*pz)/
+               TMath::Sqrt(p2*((x-fX)*(x-fX) + (y-fY)*(y-fY) + (z-fZ)*(z-fZ)));
+
+        //if (cost < (5*fCPAmax-0.9-TMath::Sqrt(r2)*(fCPAmax-1))/4.1) continue;
+         if (cost < fCPAmax) continue;
+
+         //vertex.ChangeMassHypothesis(); //default is Lambda0 
+
+         event->AddV0(&vertex);
+
+         nvtx++;
+      }
+   }
+
+   Info("Tracks2V0vertices","Number of reconstructed V0 vertices: %d",nvtx);
+
+   negtrks.Delete();
+   postrks.Delete();
+
+   return 0;
+}
+
 Int_t AliV0vertexer::Tracks2V0vertices(TTree *tTree, TTree *vTree) {
   //--------------------------------------------------------------------
   //This function reconstructs V0 vertices
@@ -65,7 +155,7 @@ Int_t AliV0vertexer::Tracks2V0vertices(TTree *tTree, TTree *vTree) {
 
 
    for (i=0; i<nneg; i++) {
-      if (i%10==0) cerr<<nneg-i<<'\r';
+      //if (i%10==0) cerr<<nneg-i<<'\r';
       AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
 
       if (TMath::Abs(ntrk->GetD(fX,fY))<fDPmin) continue;
@@ -102,9 +192,9 @@ Int_t AliV0vertexer::Tracks2V0vertices(TTree *tTree, TTree *vTree) {
          Double_t px,py,pz; vertex.GetPxPyPz(px,py,pz);
          Double_t p2=px*px+py*py+pz*pz;
          Double_t cost=((x-fX)*px + (y-fY)*py + (z-fZ)*pz)/
-                TMath::Sqrt(p2*((x-fX)*(x-fX) + (y-fY)*(y-fY) + (z-fZ)*(z-fZ)));
+               TMath::Sqrt(p2*((x-fX)*(x-fX) + (y-fY)*(y-fY) + (z-fZ)*(z-fZ)));
 
-         //if (cost < (5*fCPAmax-0.9-TMath::Sqrt(r2)*(fCPAmax-1))/4.1) continue;
+        //if (cost < (5*fCPAmax-0.9-TMath::Sqrt(r2)*(fCPAmax-1))/4.1) continue;
          if (cost < fCPAmax) continue;
 
          //vertex.ChangeMassHypothesis(); //default is Lambda0 
@@ -115,7 +205,7 @@ Int_t AliV0vertexer::Tracks2V0vertices(TTree *tTree, TTree *vTree) {
       }
    }
 
-   cerr<<"Number of reconstructed V0 vertices: "<<nvtx<<endl;
+   Info("Tracks2V0vertices","Number of reconstructed V0 vertices: %d",nvtx);
 
    negtrks.Delete();
    postrks.Delete();
index 621bca7..0f5cd52 100644 (file)
@@ -13,6 +13,7 @@
 
 class TTree;
 class AliITStrackV2;
+class AliESD;
 
 //_____________________________________________________________________________
 class AliV0vertexer : public TObject {
@@ -22,6 +23,8 @@ public:
   void SetCuts(const Double_t cuts[7]);
   void SetVertex(Double_t *vtx) { fX=vtx[0]; fY=vtx[1]; fZ=vtx[2]; }
 
+  Int_t Tracks2V0vertices(AliESD *event);
+
   Int_t Tracks2V0vertices(TTree *in, TTree *out);
   Double_t PropagateToDCA(AliITStrackV2 *nt, AliITStrackV2 *pt);
 
index be2d805..ea915ff 100644 (file)
@@ -32,9 +32,9 @@ AliESD::AliESD():
   fRunNumber(0),
   fTrigger(0),
   fRecoVersion(0),
-  fTracks("AliESDtrack",15000)
-  //fV0s("AliV0vertex",200),
-  //fCascades("AliCascadeVertex",20)
+  fTracks("AliESDtrack",15000),
+  fV0s("AliESDv0",200),
+  fCascades("AliESDcascade",20)
 {
 }
 
index 066207a..eed56ba 100644 (file)
 #include "TObject.h"
 #include "TClonesArray.h"
 #include  "AliESDtrack.h"
+#include  "AliESDv0.h"
+#include  "AliESDcascade.h"
 
 class AliESD : public TObject {
 public:
   AliESD();
   virtual ~AliESD() {
     fTracks.Delete();
-    //fV0s.Delete();
-    //fCascades.Delete();
+    fV0s.Delete();
+    fCascades.Delete();
   }
 
   void SetEventNumber(Int_t n) {fEventNumber=n;}
@@ -35,13 +37,27 @@ public:
     new(fTracks[fTracks.GetEntriesFast()]) AliESDtrack(*t);
   }
 
+  AliESDv0 *GetV0(Int_t i) {
+    return (AliESDv0 *)fV0s.UncheckedAt(i);
+  }
+  void AddV0(const AliESDv0 *v) {
+    new(fV0s[fV0s.GetEntriesFast()]) AliESDv0(*v);
+  }
+
+  AliESDcascade *GetCascade(Int_t i) {
+    return (AliESDcascade *)fCascades.UncheckedAt(i);
+  }
+  void AddCascade(const AliESDcascade *c) {
+    new(fCascades[fCascades.GetEntriesFast()]) AliESDcascade(*c);
+  }
+
   Int_t  GetEventNumber() const {return fEventNumber;}
   Int_t  GetRunNumber() const {return fRunNumber;}
   Long_t GetTrigger() const {return fTrigger;}
   
   Int_t GetNumberOfTracks()   const {return fTracks.GetEntriesFast();}
-  //Int_t GetNumberOfV0s()      const {return fV0s.GetEntriesFast();}
-  //Int_t GetNumberOfCascades() const {return fCascades.GetEntriesFast();}
+  Int_t GetNumberOfV0s()      const {return fV0s.GetEntriesFast();}
+  Int_t GetNumberOfCascades() const {return fCascades.GetEntriesFast();}
   
 protected:
 
@@ -52,8 +68,8 @@ protected:
   Int_t        fRecoVersion;     // Version of reconstruction 
 
   TClonesArray  fTracks;         // ESD tracks
-  //TClonesArray  fV0s;            // V0 vertices
-  //TClonesArray  fCascades;       // Cascade vertices
+  TClonesArray  fV0s;            // V0 vertices
+  TClonesArray  fCascades;       // Cascade vertices
   
   ClassDef(AliESD,1)  //ESD class 
 };
diff --git a/STEER/AliESDcascade.cxx b/STEER/AliESDcascade.cxx
new file mode 100644 (file)
index 0000000..73a2801
--- /dev/null
@@ -0,0 +1,161 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+//               Implementation of the cascade vertex class
+//
+//    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
+//-------------------------------------------------------------------------
+#include <TMath.h>
+
+#include "AliESDcascade.h"
+#include "AliESDv0.h"
+
+ClassImp(AliESDcascade)
+
+AliESDcascade::AliESDcascade() : TObject() {
+  //--------------------------------------------------------------------
+  // Default constructor  (Xi-)
+  //--------------------------------------------------------------------
+  fPdgCode=kXiMinus;
+  fEffMass=1.32131;
+  fChi2=1.e+33;
+  fPos[0]=fPos[1]=fPos[2]=0.;
+  fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.;
+}
+
+inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){
+  // determinant 2x2
+  return a00*a11 - a01*a10;
+}
+
+inline Double_t det (Double_t a00,Double_t a01,Double_t a02,
+                     Double_t a10,Double_t a11,Double_t a12,
+                     Double_t a20,Double_t a21,Double_t a22) {
+  // determinant 3x3
+  return 
+  a00*det(a11,a12,a21,a22)-a01*det(a10,a12,a20,a22)+a02*det(a10,a11,a20,a21);
+}
+
+Double_t AliESDcascade::ChangeMassHypothesis(Double_t &v0q, Int_t code) {
+  //--------------------------------------------------------------------
+  // This function changes the mass hypothesis for this cascade
+  // and returns the "kinematical quality" of this hypothesis
+  // together with the "quality" of associated V0 (argument v0q) 
+  //--------------------------------------------------------------------
+  Double_t nmass=0.13957, pmass=0.93827, ps0=0.101; 
+  Double_t bmass=0.13957, mass =1.3213,  ps =0.139;
+
+  fPdgCode=code;
+
+  switch (code) {
+  case 213: 
+       bmass=0.93827; 
+       break;
+  case kXiMinus:
+       break;
+  case kXiPlusBar:
+       nmass=0.93827; pmass=0.13957; 
+       break;
+  case kOmegaMinus: 
+       bmass=0.49368; mass=1.67245; ps=0.211;
+       break;
+  case kOmegaPlusBar: 
+       nmass=0.93827; pmass=0.13957; 
+       bmass=0.49368; mass=1.67245; ps=0.211;
+       break;
+  default:
+       Info("AliCascadeVertex::ChangeMassHypothesis",
+            "Invalide PDG code !  Assuming XiMinus's...");
+       fPdgCode=kXiMinus;
+    break;
+  }
+
+  Double_t pxn=fV0mom[0][0], pyn=fV0mom[0][1], pzn=fV0mom[0][2];
+  Double_t pxp=fV0mom[1][0], pyp=fV0mom[1][1], pzp=fV0mom[1][2];
+  Double_t px0=pxn+pxp, py0=pyn+pyp, pz0=pzn+pzp;
+  Double_t p0=TMath::Sqrt(px0*px0 + py0*py0 + pz0*pz0);
+
+  Double_t e0=TMath::Sqrt(1.11568*1.11568 + p0*p0);
+  Double_t beta0=p0/e0;
+  Double_t pln=(pxn*px0 + pyn*py0 + pzn*pz0)/p0;
+  Double_t plp=(pxp*px0 + pyp*py0 + pzp*pz0)/p0;
+  Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
+
+  Double_t a=(plp-pln)/(plp+pln);
+  a -= (pmass*pmass-nmass*nmass)/(1.11568*1.11568);
+  a = 0.25*beta0*beta0*1.11568*1.11568*a*a + pt2;
+
+
+  v0q=a - ps0*ps0;
+
+
+  Double_t pxb=fBachMom[0], pyb=fBachMom[1], pzb=fBachMom[2]; 
+
+  Double_t eb=TMath::Sqrt(bmass*bmass + pxb*pxb + pyb*pyb + pzb*pzb);
+  Double_t pxl=px0+pxb, pyl=py0+pyb, pzl=pz0+pzb;
+  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
+  
+  fEffMass=TMath::Sqrt((e0+eb)*(e0+eb) - pl*pl);
+
+  Double_t beta=pl/(e0+eb);
+  Double_t pl0=(px0*pxl + py0*pyl + pz0*pzl)/pl;
+  Double_t plb=(pxb*pxl + pyb*pyl + pzb*pzl)/pl;
+  pt2=p0*p0 - pl0*pl0;
+
+  a=(pl0-plb)/(pl0+plb);
+  a -= (1.11568*1.11568-bmass*bmass)/(mass*mass);
+  a = 0.25*beta*beta*mass*mass*a*a + pt2;
+
+  return (a - ps*ps);
+}
+
+void 
+AliESDcascade::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+  //--------------------------------------------------------------------
+  // This function returns the cascade momentum (global)
+  //--------------------------------------------------------------------
+  px=fV0mom[0][0]+fV0mom[1][0]+fBachMom[0]; 
+  py=fV0mom[0][1]+fV0mom[1][1]+fBachMom[1]; 
+  pz=fV0mom[0][2]+fV0mom[1][2]+fBachMom[2]; 
+}
+
+void AliESDcascade::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
+  //--------------------------------------------------------------------
+  // This function returns cascade position (global)
+  //--------------------------------------------------------------------
+  x=fPos[0]; 
+  y=fPos[1]; 
+  z=fPos[2]; 
+}
+
+Double_t AliESDcascade::GetD(Double_t x0, Double_t y0, Double_t z0) const {
+  //--------------------------------------------------------------------
+  // This function returns the cascade impact parameter
+  //--------------------------------------------------------------------
+
+  Double_t x=fPos[0],y=fPos[1],z=fPos[2];
+  Double_t px=fV0mom[0][0]+fV0mom[1][0]+fBachMom[0];
+  Double_t py=fV0mom[0][1]+fV0mom[1][1]+fBachMom[1];
+  Double_t pz=fV0mom[0][2]+fV0mom[1][2]+fBachMom[2];
+
+  Double_t dx=(y0-y)*pz - (z0-z)*py; 
+  Double_t dy=(x0-x)*pz - (z0-z)*px;
+  Double_t dz=(x0-x)*py - (y0-y)*px;
+  Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
+
+  return d;
+}
+
diff --git a/STEER/AliESDcascade.h b/STEER/AliESDcascade.h
new file mode 100644 (file)
index 0000000..7abe482
--- /dev/null
@@ -0,0 +1,70 @@
+#ifndef ALIESDCASCADE_H
+#define ALIESDCASCADE_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                        ESD Cascade Vertex Class
+//
+//    Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include <TPDGCode.h>
+
+class AliESDtrack;
+class AliESDv0;
+
+#define kXiMinus       3312
+#define kXiPlusBar    -3312
+#define kOmegaMinus    3334
+#define kOmegaPlusBar -3334
+
+class AliESDcascade : public TObject {
+public:
+  AliESDcascade();
+
+  Double_t ChangeMassHypothesis(Double_t &v0q, Int_t code=kXiMinus); 
+
+  Int_t GetPdgCode() const {return fPdgCode;}
+  Double_t GetEffMass() const {return fEffMass;}
+  Double_t GetChi2() const {return fChi2;}
+  void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const;
+  Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
+
+  void GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+     px=fV0mom[0][0]; py=fV0mom[0][1]; pz=fV0mom[0][2];
+  }
+  Int_t GetNindex() const {return fV0idx[0];}
+  void GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+     px=fV0mom[1][0]; py=fV0mom[1][1]; pz=fV0mom[1][2];
+  }
+  Int_t GetPindex() const {return fV0idx[1];}
+  void GetBPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+     px=fBachMom[0]; py=fBachMom[1]; pz=fBachMom[2];
+  }
+  Int_t GetBindex() const {return fBachIdx;}
+
+protected: 
+  Int_t fPdgCode;           // reconstructed cascade type (PDG code)
+  Double_t fEffMass;        // reconstructed cascade effective mass
+  Double_t fChi2;           // chi2 value
+  Double_t fPos[3];         // cascade vertex position (global)
+  Double_t fPosCov[6];      // covariance matrix of the vertex position
+
+  Int_t fV0idx[2];          // indeices of the V0 daughter tracks
+  Double_t fV0mom[2][3];    // V0 daughters' momenta (global)
+  Double_t fV0momCov[6];    // covariance matrix of the V0 momentum.
+
+  Int_t fBachIdx;           // label of the bachelor track
+  Double_t fBachMom[3];     // bachelor momentum (global)
+  Double_t fBachMomCov[6];  // covariance matrix of the bachelor momentum.
+
+  ClassDef(AliESDcascade,1) // reconstructed cascade vertex
+};
+
+#endif
+
+
index ffaa41e..3ddf801 100644 (file)
@@ -33,6 +33,8 @@
   #include "AliITS.h"
   #include "AliITSgeom.h"
   #include "AliITStrackerV2.h"
+  #include "AliV0vertexer.h"
+  #include "AliCascadeVertexer.h"
   #include "AliITSpidESD.h"
   #include "AliITSLoader.h"
 
@@ -113,6 +115,27 @@ Int_t AliESDtest(Int_t nev=1) {
    Double_t parITS[]={34.,0.15,10.};
    AliITSpidESD itsPID(parITS);
 
+   //An instance of the V0 finder
+   Double_t cuts[]={33,  // max. allowed chi2
+                    0.16,// min. allowed negative daughter's impact parameter 
+                    0.05,// min. allowed positive daughter's impact parameter 
+                    0.080,// max. allowed DCA between the daughter tracks
+                    0.998,// max. allowed cosine of V0's pointing angle
+                    0.9,  // min. radius of the fiducial volume
+                    2.9   // max. radius of the fiducial volume
+                   };
+   AliV0vertexer vtxer(cuts);
+
+   Double_t cts[]={33.,    // max. allowed chi2
+                    0.05,   // min. allowed V0 impact parameter 
+                    0.008,  // window around the Lambda mass 
+                    0.035,  // min. allowed bachelor's impact parameter 
+                    0.10,   // max. allowed DCA between a V0 and a track
+                    0.9985, //max. allowed cosine of the cascade pointing angle
+                    0.9,    // min. radius of the fiducial volume
+                    2.9     // max. radius of the fiducial volume
+                   };
+   AliCascadeVertexer cvtxer=AliCascadeVertexer(cts);
 
 /**** The TPC corner ********************/
 
@@ -214,9 +237,12 @@ Int_t AliESDtest(Int_t nev=1) {
      itsTracker.LoadClusters(itsTree);
      rc+=itsTracker.Clusters2Tracks(event);
 
+     rc+=vtxer.Tracks2V0vertices(event);            // V0 finding
+     rc+=cvtxer.V0sTracks2CascadeVertices(event);   // cascade finding
+
      rc+=itsTracker.PropagateBack(event); 
      itsTracker.UnloadClusters();
-     itsPID.MakePID(event);
+     //itsPID.MakePID(event);
      
      rc+=tpcTracker.PropagateBack(event);
      tpcTracker.UnloadClusters();
index 3e6f9bc..cbd6aa4 100644 (file)
@@ -103,6 +103,24 @@ Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags) {
 }
 
 //_______________________________________________________________________
+void AliESDtrack::GetExternalParametersAt(Double_t x, Double_t p[5]) const {
+  //---------------------------------------------------------------------
+  // This function returns external representation of the track parameters
+  // at the plane x
+  //---------------------------------------------------------------------
+  Double_t dx=x-fRx;
+  Double_t c=fRp[4]/AliKalmanTrack::GetConvConst();
+  Double_t f1=fRp[2], f2=f1 + c*dx;
+  Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);    
+
+  p[0]=fRp[0]+dx*(f1+f2)/(r1+r2);
+  p[1]=fRp[1]+dx*(f1+f2)/(f1*r2 + f2*r1)*fRp[3];
+  p[2]=fRp[2]+dx*c;
+  p[3]=fRp[3];
+  p[4]=fRp[4];
+}
+
+//_______________________________________________________________________
 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
   //---------------------------------------------------------------------
   // This function returns external representation of the track parameters
index 49bf6bb..ca5ddf6 100644 (file)
@@ -28,6 +28,7 @@ public:
   ULong_t GetStatus() const {return fFlags;}
   Int_t GetLabel() const {return fLabel;}
   Double_t GetAlpha() const {return fRalpha;}
+  void GetExternalParametersAt(Double_t x, Double_t p[5]) const;
   void GetExternalParameters(Double_t &x, Double_t p[5]) const;
   void GetExternalCovariance(Double_t cov[15]) const;
   Double_t GetIntegratedLength() const {return fTrackLength;}
diff --git a/STEER/AliESDv0.cxx b/STEER/AliESDv0.cxx
new file mode 100644 (file)
index 0000000..67e0ea3
--- /dev/null
@@ -0,0 +1,119 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-------------------------------------------------------------------------
+//               Implementation of the ESD V0 vertex class
+//
+//     Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <Riostream.h>
+#include <TMath.h>
+#include <TPDGCode.h>
+
+#include "AliESDv0.h"
+
+ClassImp(AliESDv0)
+
+AliESDv0::AliESDv0() : TObject() {
+  //--------------------------------------------------------------------
+  // Default constructor  (K0s)
+  //--------------------------------------------------------------------
+  fPdgCode=kK0Short;
+  fEffMass=0.497672;
+  fChi2=1.e+33;
+  fPos[0]=fPos[1]=fPos[2]=0.;
+  fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.;
+}
+
+Double_t AliESDv0::ChangeMassHypothesis(Int_t code) {
+  //--------------------------------------------------------------------
+  // This function changes the mass hypothesis for this V0
+  // and returns the "kinematical quality" of this hypothesis 
+  //--------------------------------------------------------------------
+  Double_t nmass=0.13957, pmass=0.13957, mass=0.49767, ps=0.206;
+
+  fPdgCode=code;
+
+  switch (code) {
+  case kLambda0:
+    nmass=0.13957; pmass=0.93827; mass=1.1157; ps=0.101; break;
+  case kLambda0Bar:
+    pmass=0.13957; nmass=0.93827; mass=1.1157; ps=0.101; break;
+  case kK0Short: 
+    break;
+  default:
+    cerr<<"AliV0vertex::ChangeMassHypothesis: ";
+    cerr<<"invalide PDG code ! Assuming K0s...\n";
+    fPdgCode=kK0Short;
+    break;
+  }
+
+  Double_t pxn=fNmom[0], pyn=fNmom[1], pzn=fNmom[2]; 
+  Double_t pxp=fPmom[0], pyp=fPmom[1], pzp=fPmom[2];
+
+  Double_t en=TMath::Sqrt(nmass*nmass + pxn*pxn + pyn*pyn + pzn*pzn);
+  Double_t ep=TMath::Sqrt(pmass*pmass + pxp*pxp + pyp*pyp + pzp*pzp);
+  Double_t pxl=pxn+pxp, pyl=pyn+pyp, pzl=pzn+pzp;
+  Double_t pl=TMath::Sqrt(pxl*pxl + pyl*pyl + pzl*pzl);
+
+  fEffMass=TMath::Sqrt((en+ep)*(en+ep)-pl*pl);
+
+  Double_t beta=pl/(en+ep);
+  Double_t pln=(pxn*pxl + pyn*pyl + pzn*pzl)/pl;
+  Double_t plp=(pxp*pxl + pyp*pyl + pzp*pzl)/pl;
+
+  Double_t pt2=pxp*pxp + pyp*pyp + pzp*pzp - plp*plp;
+
+  Double_t a=(plp-pln)/(plp+pln);
+  a -= (pmass*pmass-nmass*nmass)/(mass*mass);
+  a = 0.25*beta*beta*mass*mass*a*a + pt2;
+
+  return (a - ps*ps);
+  
+}
+
+void AliESDv0::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's momentum (global)
+  //--------------------------------------------------------------------
+  px=fNmom[0]+fPmom[0]; 
+  py=fNmom[1]+fPmom[1]; 
+  pz=fNmom[2]+fPmom[2]; 
+}
+
+void AliESDv0::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's position (global)
+  //--------------------------------------------------------------------
+  x=fPos[0]; 
+  y=fPos[1]; 
+  z=fPos[2]; 
+}
+
+Double_t AliESDv0::GetD(Double_t x0, Double_t y0, Double_t z0) const {
+  //--------------------------------------------------------------------
+  // This function returns V0's impact parameter
+  //--------------------------------------------------------------------
+  Double_t x=fPos[0],y=fPos[1],z=fPos[2];
+  Double_t px=fNmom[0]+fPmom[0];
+  Double_t py=fNmom[1]+fPmom[1];
+  Double_t pz=fNmom[2]+fPmom[2];
+
+  Double_t dx=(y0-y)*pz - (z0-z)*py; 
+  Double_t dy=(x0-x)*pz - (z0-z)*px;
+  Double_t dz=(x0-x)*py - (y0-y)*px;
+  Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz));
+  return d;
+}
diff --git a/STEER/AliESDv0.h b/STEER/AliESDv0.h
new file mode 100644 (file)
index 0000000..837dccf
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef ALIESDV0_H
+#define ALIESDV0_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//-------------------------------------------------------------------------
+//                          ESD V0 Vertex Class
+//
+//    Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include <TPDGCode.h>
+
+class AliESDtrack;
+
+class AliESDv0 : public TObject {
+public:
+  AliESDv0();
+
+  Double_t ChangeMassHypothesis(Int_t code=kK0Short); 
+
+  Int_t GetPdgCode() const {return fPdgCode;}
+  Double_t GetEffMass() const {return fEffMass;}
+  Double_t GetChi2() const {return fChi2;}
+  void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const;
+  void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const;
+  Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
+  Int_t GetNindex() const {return fNidx;}
+  Int_t GetPindex() const {return fPidx;}
+
+protected: 
+  Int_t fPdgCode;           // reconstructed V0's type (PDG code)
+  Double_t fEffMass;        // reconstructed V0's effective mass
+  Double_t fChi2;           // V0's chi2 value
+  Double_t fPos[3];         // V0's position (global)
+  Double_t fPosCov[6];      // covariance matrix of the vertex position
+
+  Int_t fNidx;              // index of the negative daughter
+  Double_t fNmom[3];        // momentum of the negative daughter (global)
+  Double_t fNmomCov[6];     // covariance matrix of the negative daughter mom.
+
+  Int_t fPidx;              // index of the positive daughter
+  Double_t fPmom[3];        // momentum of the positive daughter (global)
+  Double_t fPmomCov[6];     // covariance matrix of the positive daughter mom.
+
+  ClassDef(AliESDv0,1)      // ESD V0 vertex
+};
+
+inline 
+void AliESDv0::GetNPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+px=fNmom[0]; py=fNmom[1]; pz=fNmom[2];
+}
+
+inline 
+void AliESDv0::GetPPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const {
+px=fPmom[0]; py=fPmom[1]; pz=fPmom[2];
+}
+
+#endif
+
+
diff --git a/STEER/AliESDv0Analysis.C b/STEER/AliESDv0Analysis.C
new file mode 100644 (file)
index 0000000..101dbd8
--- /dev/null
@@ -0,0 +1,88 @@
+//********************************************************************
+//     Example (very naive for the moment) of the data analysis 
+//                    using the ESD classes.
+//       It demonstrates the idea of the "combined PID" 
+//            applied to the Lambda0 reconstruction. 
+//      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+//********************************************************************
+
+#if !defined( __CINT__) || defined(__MAKECINT__)
+  #include <Riostream.h>
+  #include "TKey.h"
+  #include "TFile.h"
+  #include "TH1F.h"
+  #include "TH2F.h"
+  #include "TCanvas.h"
+  #include "TStopwatch.h"
+  #include "TParticle.h"
+
+  #include "AliRun.h"
+
+  #include "AliESD.h"
+
+#endif
+
+extern AliRun *gAlice;
+
+Int_t AliESDv0Analysis(Int_t nev=1) { 
+   TH1F *hm=new TH1F("hm","Effective Mass",40,1.065,1.165);
+   hm->SetXTitle("Mass (GeV/c**2)");
+
+   TFile *ef=TFile::Open("AliESDs.root");
+   if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+
+   TStopwatch timer;
+   Int_t rc=0,n=0;
+   TKey *key=0;
+   TIter next(ef->GetListOfKeys());
+
+   //****** Tentative particle type "concentrations"
+   Double_t c[5]={0.0, 0.0, 0.1, 0.1, 0.1};
+
+   //******* The loop over events
+   while ((key=(TKey*)next())!=0) {
+
+     cerr<<"Processing event number : "<<n++<<endl;
+
+     AliESD *event=(AliESD*)key->ReadObj();
+
+     Int_t nv0=event->GetNumberOfV0s();
+     cerr<<"Number of ESD v0s : "<<nv0<<endl; 
+
+     while (nv0--) {
+       AliESDv0 *v0=event->GetV0(nv0);
+       Int_t pi=v0->GetPindex();
+       AliESDtrack *t=event->GetTrack(pi);
+       Int_t isProton=1;
+       if ((t->GetStatus()&AliESDtrack::kESDpid)!=0) {
+        Double_t r[10]; t->GetESDpid(r);
+         Double_t rcc=0.;
+         Int_t i;
+         for (i=0; i<AliESDtrack::kSPECIES; i++) rcc+=(c[i]*r[i]);
+         if (rcc==0.) continue;
+         //Here we apply Bayes' formula
+         Double_t w[10];
+         for (i=0; i<AliESDtrack::kSPECIES; i++) w[i]=c[i]*r[i]/rcc;
+
+         if (w[4]<w[3]) isProton=0;
+         if (w[4]<w[2]) isProton=0;
+         if (w[4]<w[1]) isProton=0;
+        if (w[4]<w[0]) isProton=0;
+       }
+       if (!isProton) continue;
+       v0->ChangeMassHypothesis(3122);
+       Double_t mass=v0->GetEffMass();
+       hm->Fill(mass);
+     } 
+     delete event;
+   }
+
+   timer.Stop(); timer.Print();
+
+   hm->Draw();
+
+   ef->Close();
+
+
+   return rc;
+}
index f58cf12..027ffec 100644 (file)
@@ -87,7 +87,7 @@ public:
 
   static void SetConvConst(Double_t cc) {fgConvConst=cc;}
   static void SetConvConst();
-  Double_t GetConvConst() const {return fgConvConst;}
+  static Double_t GetConvConst() {return fgConvConst;}
 
   static void SetMagneticField(Double_t f) {// f - Magnetic field in T
     fgConvConst=100/0.299792458/f;
index c121df1..543354e 100644 (file)
@@ -65,6 +65,8 @@
 #pragma link C++ class  AliTrackReference+;
 #pragma link C++ class  AliESD+;
 #pragma link C++ class  AliESDtrack+;
+#pragma link C++ class  AliESDv0+;
+#pragma link C++ class  AliESDcascade+;
 #pragma link C++ class  AliESDvertex+;
 #pragma link C++ class  AliESDpid+;
 #pragma link C++ class  AliTrackMap+;
index a296811..262cd28 100644 (file)
@@ -16,7 +16,7 @@ AliMergeCombi.cxx AliMagFMaps.cxx AliFieldMap.cxx \
 AliGausCorr.cxx AliTrackReference.cxx AliESD.cxx \
 AliTrackMap.cxx AliTrackMapper.cxx AliCollisionGeometry.cxx \
 AliMemoryWatcher.cxx AliBarrelTrack.cxx \
-AliESDtrack.cxx AliESDvertex.cxx AliESDpid.cxx \
+AliESDtrack.cxx AliESDv0.cxx AliESDcascade.cxx AliESDvertex.cxx AliESDpid.cxx \
 AliRawReader.cxx
 
 HDRS:= $(SRCS:.cxx=.h)