#endif
struct GoodTrackITS {
- Int_t event;
Int_t lab;
Int_t code;
Float_t px,py,pz;
Float_t x,y,z;
};
-Int_t AliITSComparisonV2() {
- Int_t good_tracks_its(GoodTrackITS *gt, Int_t max);
-
+Int_t AliITSComparisonV2(Int_t event=0) {
cerr<<"Doing comparison...\n";
- TFile *cf=TFile::Open("AliITSclustersV2.root");
- if (!cf->IsOpen()) {cerr<<"Can't open AliITSclustersV2.root !\n"; return 1;}
- AliITSgeom *geom=(AliITSgeom*)cf->Get("AliITSgeom");
- if (!geom) { cerr<<"Can't get the ITS geometry !\n"; return 2; }
- AliITStrackerV2 tracker(geom);
-
-// Load tracks
- TFile *tf=TFile::Open("AliITStracksV2.root");
- if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
- TObjArray tarray(2000);
- TTree *tracktree=(TTree*)tf->Get("TreeT_ITS_0");
- if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
- TBranch *tbranch=tracktree->GetBranch("tracks");
- Int_t nentr=(Int_t)tracktree->GetEntries(),i;
- for (i=0; i<nentr; i++) {
- AliITStrackV2 *iotrack=new AliITStrackV2;
- tbranch->SetAddress(&iotrack);
- tracktree->GetEvent(i);
-
- Int_t tpcLabel=iotrack->GetLabel();
- tracker.CookLabel(iotrack,0.);
- Int_t itsLabel=iotrack->GetLabel();
- if (itsLabel != tpcLabel) iotrack->SetLabel(-TMath::Abs(itsLabel));
- if (tpcLabel < 0) iotrack->SetLabel(-TMath::Abs(itsLabel));
- /*
- if (itsLabel==1234) {
- Int_t nc=iotrack->GetNumberOfClusters();
- for (Int_t k=0; k<nc; k++) {
- Int_t index=iotrack->GetClusterIndex(k);
- AliITSclusterV2 *c=tracker.GetCluster(index);
- cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
- }
- }
- */
- tarray.AddLast(iotrack);
+ const Int_t MAX=15000;
+ Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event);
+
+ Int_t nentr=0; TObjArray tarray(2000);
+ {/* Load tracks */
+ TFile *tf=TFile::Open("AliITStracksV2.root");
+ if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
+ char tname[100]; sprintf(tname,"TreeT_ITS_%d",event);
+ TTree *tracktree=(TTree*)tf->Get(tname);
+ if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ nentr=(Int_t)tracktree->GetEntries();
+ for (Int_t i=0; i<nentr; i++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tarray.AddLast(iotrack);
+ /*if (itsLabel != 1234) continue;
+ Int_t nc=iotrack->GetNumberOfClusters();
+ for (Int_t k=0; k<nc; k++) {
+ Int_t index=iotrack->GetClusterIndex(k);
+ AliITSclusterV2 *c=tracker.GetCluster(index);
+ cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
+ }*/
+ }
+ delete tracktree; //Thanks to Mariana Bondila
+ tf->Close();
}
- delete tracktree; //Thanks to Mariana Bondila
- tf->Close();
- delete geom; //Thanks to Mariana Bondila
- cf->Close();
-/////////////////////////////////////////////////////////////////////////
- const Int_t MAX=15000;
+ /* Generate a list of "good" tracks */
GoodTrackITS gt[MAX];
Int_t ngood=0;
ifstream in("good_tracks_its");
if (in) {
cerr<<"Reading good tracks...\n";
- while (in>>gt[ngood].event>>gt[ngood].lab>>gt[ngood].code>>
+ while (in>>gt[ngood].lab>>gt[ngood].code>>
gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
ngood++;
if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
} else {
cerr<<"Marking good tracks (this will take a while)...\n";
- ngood=good_tracks_its(gt,MAX);
+ ngood=good_tracks_its(gt,MAX,event);
ofstream out("good_tracks_its");
if (out) {
for (Int_t ngd=0; ngd<ngood; ngd++)
- out<<gt[ngd].event<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
+ out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
<<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
<<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
} else cerr<<"Can not open file (good_tracks_its) !\n";
return 0;
}
-Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) {
+Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event) {
if (gAlice) {delete gAlice; gAlice=0;}
TFile *file=TFile::Open("galice.root");
exit(5);
}
- Int_t np=gAlice->GetEvent(0);
+ Int_t np=gAlice->GetEvent(event);
Int_t *good=new Int_t[np];
Int_t k;
if (!cf->IsOpen()){
cerr<<"Can't open AliITSclustersV2.root !\n"; exit(6);
}
- TTree *cTree=(TTree*)cf->Get("TreeC_ITS_0");
+ char cname[100]; sprintf(cname,"TreeC_ITS_%d",event);
+ TTree *cTree=(TTree*)cf->Get(cname);
if (!cTree) {
cerr<<"Can't get cTree !\n"; exit(7);
}
}
Int_t nt=0;
Double_t px,py,pz,x,y,z;
- Int_t code,lab,event;
- while (in>>event>>lab>>code>>px>>py>>pz>>x>>y>>z) {
+ Int_t code,lab;
+ while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
if (good[lab] != 0x3F) continue;
TParticle *p = (TParticle*)gAlice->Particle(lab);
- gt[nt].event=event;
gt[nt].lab=lab;
gt[nt].code=p->GetPdgCode();
//**** px py pz - in global coordinate system
geom->Write();
TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
- //TTree *cTree=new TTree("cTree","ITS clusters");
TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters");
cTree->Branch("Clusters",&clusters);
cerr<<"Number of entries: "<<nentr<<endl;
+ Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ?
+
for (Int_t i=0; i<nentr; i++) {
points->Clear();
pTree->GetEvent(i);
nclusters+=ncl;
for (Int_t j=0; j<ncl; j++) {
AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
- Float_t lp[5];
+ //Float_t lp[5];
lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
lp[1]=p->GetZ()+zshift;
lp[2]=p->GetSigmaX2();
lp[3]=p->GetSigmaZ2();
lp[4]=p->GetQ();
- Int_t lab[6];
+ //Int_t lab[6];
lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
lab[3]=ndet;
new(cl[j]) AliITSclusterV2(lab,lp);
}
cTree->Fill(); clusters->Delete();
- points->Delete();
}
cTree->Write();
#ifndef __CINT__
#include <iostream.h>
+ #include "AliITSgeom.h"
#include "AliITStrackerV2.h"
#include "TFile.h"
TStopwatch timer;
AliITStrackerV2 tracker(geom);
+
+ //Double_t xyz[]={0.,0.,0.}; tracker.SetVertex(xyz); //primary vertex
+ //Int_t flag[]={1}; //some default flags
+ //flag[0]= 0; tracker.SetupFirstPass(flag); //no constraint
+ //flag[0]=-1; tracker.SetupSecondPass(flag); //skip second pass
+
Int_t rc=tracker.Clusters2Tracks(in,out);
timer.Stop(); timer.Print();
1.44e-4, 1.44e-4, 7.84e-6, 7.84e-6, 0.006889, 0.006889
};
- const Double_t kChi2PerCluster=7.;//10.;//7
- const Double_t kMaxChi2=15.;//20.; //15.
- const Double_t kMaxRoad=13.;
+ const Double_t kChi2PerCluster=7.;//5.;
+ const Double_t kMaxChi2=15.;//12;
+ const Double_t kMaxRoad=3.0;
const Double_t kSigmaYV=0.005e+0;
const Double_t kSigmaZV=0.010e+0;
ClassImp(AliITStrackV2)
-const Int_t kWARN=1;
+const Int_t kWARN=5;
//____________________________________________________________________________
AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) {
SetNumberOfClusters(0);
//SetConvConst(t.GetConvConst());
- fdEdx = 0.;
+ fdEdx = t.GetdEdx();
+ SetMass(t.GetMass());
+
fAlpha = t.GetAlpha();
if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
Int_t n=GetNumberOfClusters();
for (Int_t i=0; i<n; i++) fIndex[i]=t.fIndex[i];
}
-
+/*
//_____________________________________________________________________________
Int_t AliITStrackV2::Compare(const TObject *o) const {
//-----------------------------------------------------------------
// This function compares tracks according to the their curvature
//-----------------------------------------------------------------
- AliTPCtrack *t=(AliTPCtrack*)o;
+ AliITStrackV2 *t=(AliITStrackV2*)o;
Double_t co=TMath::Abs(t->Get1Pt());
Double_t c =TMath::Abs(Get1Pt());
if (c>co) return 1;
else if (c<co) return -1;
return 0;
}
+*/
+//_____________________________________________________________________________
+Int_t AliITStrackV2::Compare(const TObject *o) const {
+ //-----------------------------------------------------------------
+ // This function compares tracks according to the their curvature
+ //-----------------------------------------------------------------
+ AliITStrackV2 *t=(AliITStrackV2*)o;
+
+ Double_t p2=1./(Get1Pt()*Get1Pt());
+ Double_t b2=p2/(p2 + GetMass()*GetMass());
+ Double_t po2=1./(t->Get1Pt()*t->Get1Pt());
+ Double_t bo2=po2/(po2 + t->GetMass()*t->GetMass());
+ if (p2*b2>po2*bo2) return -1;
+ else if (p2*b2<po2*bo2) return 1;
+ return 0;
+}
//_____________________________________________________________________________
void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
}
//____________________________________________________________________________
-Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) {
+Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho) {
//------------------------------------------------------------------
//This function propagates a track to the minimal distance from the origin
//------------------------------------------------------------------
Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
- Propagate(fAlpha,xv,0.,0.,pm);
+ Propagate(fAlpha,xv,0.,0.);
return 0;
}
Double_t dx=xk-fX;
Double_t f1=fP2, f2=f1 + fP4*dx;
if (TMath::Abs(f2) >= 0.99999) {
- Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+ Int_t n=GetNumberOfClusters();
+ if (n>kWARN)
+ cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\n";
return 0;
}
r00+=fC00; r01+=fC10; r11+=fC11;
Double_t det=r00*r11 - r01*r01;
- if (TMath::Abs(det) < 1.e-10) {
+ if (TMath::Abs(det) < 1.e-30) {
Int_t n=GetNumberOfClusters();
- if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+ if (n>kWARN)
+ cerr<<n<<" AliKalmanTrack::GetPredictedChi2: Singular matrix !\n";
return 1e10;
}
Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
//_____________________________________________________________________________
Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
-Double_t x0, Double_t pm) const {
+Double_t x0) const {
//-----------------------------------------------------------------
// This function calculates a chi2 increment with a vertex contraint
//-----------------------------------------------------------------
//x0=0.;
if (x0!=0.) {
Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=pp2/(pp2 + pm*pm);
+ Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
v22 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
Double_t det=R.Determinant();
- if (TMath::Abs(det) < 1.e-25) {
+ if (TMath::Abs(det) < 1.e-30) {
Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Singular matrix !\n";
+ if (n>kWARN)
+ cerr<<n<<" AliITStrackV2::GetPredictedChi2: Singular matrix !\n";
return 1e10;
}
//____________________________________________________________________________
Int_t
-AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) {
+AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
//------------------------------------------------------------------
//This function propagates a track
//------------------------------------------------------------------
Double_t f1=fP2, f2=f1 + fP4*dx;
if (TMath::Abs(f2) >= 0.99999) {
Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+ if (n>kWARN)
+ cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\n";
return 0;
}
fX=x2;
Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + pm*pm);
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
//Multiple scattering******************
//x0=0.;
rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
if (x1 < x2) dE=-dE;
- fP4*=(1.- sqrt(p2+pm*pm)/p2*dE);
+ fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
}
if (!Invariant()) {cout<<"Propagate !\n"; return 0;}
Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
Double_t sf=fP2 + k20*dy + k21*dz;
- /*
- if (TMath::Abs(sf) >= 0.99999) {
- Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Filtering failed !\n";
- return 0;
- }
- */
+
fP0 += k00*dy + k01*dz;
fP1 += k10*dy + k11*dz;
fP2 = sf;
//------------------------------------------------------------------
// This function is for debugging purpose only
//------------------------------------------------------------------
+ Int_t n=GetNumberOfClusters();
+
//if (TMath::Abs(fP1)>11.5)
- //if (fP1*fP4<0) {cout<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
- if (TMath::Abs(fP2)>=1) {cout<<"fP2="<<fP2<<endl; return 0;}
-
- if (fC00<=0) {cout<<"fC00="<<fC00<<endl; return 0;}
- if (fC11<=0) {cout<<"fC11="<<fC11<<endl; return 0;}
- if (fC22<=0) {cout<<"fC22="<<fC22<<endl; return 0;}
- if (fC33<=0) {cout<<"fC33="<<fC33<<endl; return 0;}
- if (fC44<=0) {cout<<"fC44="<<fC44<<endl; return 0;}
+ //if (fP1*fP4<0) {
+ // if (n>kWARN) cout<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
+
+ if (TMath::Abs(fP2)>=1) {if (n>kWARN) cout<<"fP2="<<fP2<<endl; return 0;}
+
+ if (fC00<=0) {if (n>kWARN) cout<<"fC00="<<fC00<<endl; return 0;}
+ if (fC11<=0) {if (n>kWARN) cout<<"fC11="<<fC11<<endl; return 0;}
+ if (fC22<=0) {if (n>kWARN) cout<<"fC22="<<fC22<<endl; return 0;}
+ if (fC33<=0) {if (n>kWARN) cout<<"fC33="<<fC33<<endl; return 0;}
+ if (fC44<=0) {if (n>kWARN) cout<<"fC44="<<fC44<<endl; return 0;}
/*
TMatrixD m(5,5);
m(0,0)=fC00;
Double_t det=m.Determinant();
if (det <= 0) {
- cout<<" bad determinant "<<det<<endl;
- m.Print();
+ if (n>kWARN) { cout<<" bad determinant "<<det<<endl; m.Print(); }
return 0;
}
*/
}
//____________________________________________________________________________
-Int_t AliITStrackV2::Propagate(Double_t alp, Double_t xk,
-Double_t x0,Double_t rho,Double_t pm) {
+Int_t
+AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
//------------------------------------------------------------------
//This function propagates a track
//------------------------------------------------------------------
Double_t pp2=fP2*ca - cf*sa;
if (TMath::Abs(pp2) >= 0.99999) {
Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Rotation failed !\n";
+ if (n>kWARN)
+ cerr<<n<<" AliITStrackV2::Propagate: Rotation failed !\n";
return 0;
}
Double_t f1=fP2, f2=f1 + fP4*dx;
if (TMath::Abs(f2) >= 0.99999) {
Int_t n=GetNumberOfClusters();
- if (n>kWARN) cerr<<n<<" AliITStrackV2 warning: Propagation failed !\n";
+ if (n>kWARN)
+ cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
return 0;
}
fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=pp2/(pp2 + pm*pm);
+ Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
//Multiple scattering******************
//x0=0.;
rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
if (x1 < x2) dE=-dE;
- fP4*=(1.- sqrt(pp2+pm*pm)/pp2*dE);
+ fP4*=(1.- sqrt(pp2+GetMass()*GetMass())/pp2*dE);
}
if (!Invariant()) {
Double_t dy=fP0-yv, dz=fP1-zv;
Double_t r2=fX*fX+dy*dy;
Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + 0.14*0.14);
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
Double_t r2=fX*fX+fP0*fP0;
Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + 0.14*0.14);
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
AliITStrackV2():AliKalmanTrack(){}
AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *);
AliITStrackV2(const AliITStrackV2& t);
- Int_t
- PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
- void SetdEdx(Float_t dedx) {fdEdx=dedx;}
+ Int_t PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3);
+ Int_t Propagate(Double_t alpha, Double_t xr, Double_t x0, Double_t rho);
+ Int_t PropagateTo(Double_t xr,Double_t x0=21.82,Double_t rho=2.33);
+ Int_t Update(const Double_t *m, Double_t chi2, UInt_t i);
+ Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
+ Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
+ void SetdEdx(Double_t dedx) {fdEdx=dedx;}
+ void CookdEdx(Double_t low=0., Double_t up=1.) {}
void SetDetectorIndex(Int_t i) {SetLabel(i);}
+ void ResetCovariance();
+ void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
void *operator new(size_t s,void *p) { return p; }
void *operator new(size_t s) { return ::operator new(s); }
Int_t GetDetectorIndex() const {return GetLabel();}
Double_t GetX() const {return fX;}
Double_t GetAlpha()const {return fAlpha;}
- Float_t GetdEdx() const {return fdEdx;}
+ Double_t GetdEdx() const {return fdEdx;}
Double_t GetY() const {return fP0;}
Double_t GetZ() const {return fP1;}
Double_t GetSnp() const {return fP2;}
Double_t GetD() const;
Double_t GetSigmaY2() const {return fC00;}
Double_t GetSigmaZ2() const {return fC11;}
-
Int_t Compare(const TObject *o) const;
-
void GetExternalParameters(Double_t& xr, Double_t x[5]) const ;
void GetExternalCovariance(Double_t cov[15]) const ;
-
Int_t GetClusterIndex(Int_t i) const {return fIndex[i];}
Int_t GetGlobalXYZat(Double_t r,Double_t &x,Double_t &y,Double_t &z) const;
-
- Int_t Propagate(Double_t alpha,
- Double_t xr,Double_t x0,Double_t rho,Double_t pm=0.139);
-
Double_t GetPredictedChi2(const AliCluster *cluster) const;
- Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
-
- Double_t GetPredictedChi2(const AliCluster *cluster, Double_t *m,
- Double_t x0, Double_t pm=0.139) const;
- Int_t Update(const Double_t *m, Double_t chi2, UInt_t i);
- Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
- void ResetCovariance();
- void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
-
+ Double_t
+ GetPredictedChi2(const AliCluster *cluster, Double_t *m, Double_t x0) const;
Int_t Invariant() const;
-
- //protected:
-Int_t
-PropagateTo(Double_t xr,Double_t x0=21.82,Double_t rho=2.33,Double_t pm=0.139);
private:
Double_t fX; // X-coordinate of this track (reference plane)
//#define DEBUG
#ifdef DEBUG
-Int_t LAB=236;
+Int_t LAB=7;
#endif
-extern TRandom *gRandom;
-
AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
AliITStrackerV2::
-AliITStrackerV2(const AliITSgeom *geom, Int_t eventn) throw (const Char_t *) {
+AliITStrackerV2(const AliITSgeom *geom, Int_t eventn) throw (const Char_t *) :
+AliTracker() {
//--------------------------------------------------------------------
//This is the AliITStracker constructor
//It also reads clusters from a root file
//--------------------------------------------------------------------
-
- fEventN = eventn; //MI change add event number - used to generate identifier
- fYV=fZV=0.;
+ fEventN=eventn; //MI change add event number - used to generate identifier
AliITSgeom *g=(AliITSgeom*)geom;
cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
#endif
- AliITSdetector &det=fLayers[lay-1].GetDetector(c->GetDetectorIndex());
- Double_t r=det.GetR(); r=TMath::Sqrt(r*r + c->GetY()*c->GetY());
- if (r > TMath::Abs(c->GetZ())-0.5)
- fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
+ fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
}
clusters->Delete();
}
}
fI=kMaxLayer;
+
+ fPass=0;
+ fConstraint[0]=1; fConstraint[1]=0;
}
#ifdef DEBUG
return 2;
}
- in->cd();
- //
+ in->cd();
+
char tname[100];
- sprintf(tname,"TreeT_TPC_%d",fEventN);
- TTree *tpcTree=(TTree*)in->Get(tname);
-
- if (!tpcTree) {
- cerr<<"AliITStrackerV2::Clusters2Tracks() ";
- cerr<<"can't get a tree with TPC tracks !\n";
- return 3;
+ Int_t nentr=0; TObjArray itsTracks(10000);
+
+ {/* Read TPC tracks */
+ sprintf(tname,"TreeT_TPC_%d",fEventN);
+ TTree *tpcTree=(TTree*)in->Get(tname);
+ if (!tpcTree) {
+ cerr<<"AliITStrackerV2::Clusters2Tracks(): "
+ "can't get a tree with TPC tracks !\n";
+ return 3;
+ }
+ AliTPCtrack *itrack=new AliTPCtrack;
+ tpcTree->SetBranchAddress("tracks",&itrack);
+ nentr=(Int_t)tpcTree->GetEntries();
+ for (Int_t i=0; i<nentr; i++) {
+ tpcTree->GetEvent(i);
+ itsTracks.AddLast(new AliITStrackV2(*itrack));
+ }
+ delete tpcTree; //Thanks to Mariana Bondila
+ delete itrack;
}
-
- AliTPCtrack *itrack=new AliTPCtrack;
- tpcTree->SetBranchAddress("tracks",&itrack);
+ itsTracks.Sort();
out->cd();
- TTree itsTree("ITSf","Tree with ITS tracks");
+
+ sprintf(tname,"TreeT_ITS_%d",fEventN);
+ TTree itsTree(tname,"Tree with ITS tracks");
AliITStrackV2 *otrack=&fBestTrack;
itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
- Int_t ntrk=0;
-
- Int_t nentr=(Int_t)tpcTree->GetEntries();
- for (Int_t i=0; i<nentr; i++) {
-
- if (!tpcTree->GetEvent(i)) continue;
- Int_t tpcLabel=itrack->GetLabel(); //save the TPC track label
-
+ for (fPass=0; fPass<2; fPass++) {
+ Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
+ for (Int_t i=0; i<nentr; i++) {
+ if (i%10==0) cerr<<nentr-i<<" \r";
+ AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
+ if (t==0) continue; //this track has been already tracked
+ Int_t tpcLabel=t->GetLabel(); //save the TPC track label
#ifdef DEBUG
lbl=tpcLabel;
if (TMath::Abs(tpcLabel)!=LAB) continue;
cout<<tpcLabel<<" *****************\n";
#endif
+ try {
+ ResetTrackToFollow(*t);
+ } catch (const Char_t *msg) {
+ cerr<<msg<<endl;
+ continue;
+ }
+ ResetBestTrack();
- try {
- ResetTrackToFollow(AliITStrackV2(*itrack));
- } catch (const Char_t *msg) {
- cerr<<msg<<endl;
- continue;
- }
- ResetBestTrack();
-
- fYV=gRandom->Gaus(0.,kSigmaYV); fZV=gRandom->Gaus(0.,kSigmaZV);
-
- Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX();
- Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
- fTrackToFollow.Improve(x0,fYV,fZV);
-
- Double_t xk=80.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
-
- xk-=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk-=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk-=2.0;
- fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
- xk-=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk-=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
-
- xk=61.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //C02
-
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk -=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk -=0.5;
- fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex
- xk -=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
-
-
- for (FollowProlongation(); fI<kMaxLayer; fI++) {
- while (TakeNextProlongation()) FollowProlongation();
- }
+ Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX();
+ Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
+ if (constraint) fTrackToFollow.Improve(x0,GetY(),GetZ());
+
+ Double_t xk=80.;
+ fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
+
+ xk-=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk-=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk-=2.0;
+ fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
+ xk-=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk-=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+
+ xk=61.;
+ fTrackToFollow.PropagateTo(xk,0.,0.); //C02
+
+ xk -=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
+ xk -=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk -=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk -=0.5;
+ fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex
+ xk -=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk -=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk -=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
+
+
+ for (FollowProlongation(); fI<kMaxLayer; fI++) {
+ while (TakeNextProlongation()) FollowProlongation();
+ }
#ifdef DEBUG
cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
#endif
- if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
- cerr << ++ntrk << " \r";
- fBestTrack.SetLabel(tpcLabel);
- UseClusters(&fBestTrack);
- itsTree.Fill();
- }
+ if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
+ fBestTrack.SetLabel(tpcLabel);
+ fBestTrack.CookdEdx();
+ CookLabel(&fBestTrack,0.); //For comparison only
+ itsTree.Fill();
+ UseClusters(&fBestTrack);
+ delete itsTracks.RemoveAt(i);
+ }
+ }
}
- sprintf(tname,"TreeT_ITS_%d",fEventN);
- itsTree.Write(tname);
- savedir->cd();
- cerr<<"Number of TPC tracks: "<<nentr<<endl;
- cerr<<"Number of prolonged tracks: "<<ntrk<<endl;
- delete tpcTree; //Thanks to Mariana Bondila
+ itsTree.Write();
- delete itrack;
+ itsTracks.Delete();
+ savedir->cd();
+ cerr<<"Number of TPC tracks: "<<nentr<<endl;
+ cerr<<"Number of prolonged tracks: "<<itsTree.GetEntries()<<endl;
return 0;
}
Double_t rs=0.5*(fLayers[fI+1].GetR() + r);
Double_t ds=0.034; if (fI==3) ds=0.039;
Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
- fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,1*dx0r,dr);
+ fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,dx0r,dr);
}
//find intersection
Double_t x,y,z;
if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
- cerr<<"AliITStrackerV2::FollowProlongation: failed to estimate track !\n";
+ //cerr<<"AliITStrackerV2::FollowProlongation: "
+ //"failed to estimate track !\n";
break;
}
Double_t phi=TMath::ATan2(y,x);
Double_t d=layer.GetThickness(phi,z);
Int_t idet=layer.FindDetectorIndex(phi,z);
if (idet<0) {
- cerr<<"AliITStrackerV2::FollowProlongation: failed to find a detector !\n";
+ //cerr<<"AliITStrackerV2::FollowProlongation: "
+ //"failed to find a detector !\n";
break;
}
//propagate to the intersection
const AliITSdetector &det=layer.GetDetector(idet);
phi=det.GetPhi();
- if (!fTrackToFollow.Propagate(phi,det.GetR(),1*d/21.82*2.33,d*2.33)) {
- cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
+ if (!fTrackToFollow.Propagate(phi,det.GetR(),d/21.82*2.33,d*2.33)) {
+ //cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
break;
}
+ if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+1) break;
fTrackToFollow.SetDetectorIndex(idet);
//Select possible prolongations and store the current track estimation
track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[fI]);
if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
- if (dz > kMaxRoad/4) {
+ if (dz > kMaxRoad) {
//cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
break;
}
Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]);
if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
- if (dy > kMaxRoad/4) {
+ if (dy > kMaxRoad) {
//cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
break;
}
-
- Double_t zmin=track.GetZ() - dz;
+ Double_t zmin=track.GetZ() - dz;
Double_t zmax=track.GetZ() + dz;
Double_t ymin=track.GetY() + r*phi - dy;
Double_t ymax=track.GetY() + r*phi + dy;
AliITSlayer &layer=fLayers[fI];
AliITStrackV2 &t=fTracks[fI];
+ Int_t &constraint=fConstraint[fPass];
Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
if (t.GetDetectorIndex()!=idet) {
const AliITSdetector &det=layer.GetDetector(idet);
if (!t.Propagate(det.GetPhi(),det.GetR(),0.,0.)) {
- cerr<<"AliITStrackerV2::TakeNextProlongation: propagation failed !\n";
+ //cerr<<"AliITStrackerV2::TakeNextProlongation: "
+ //"propagation failed !\n";
continue;
}
t.SetDetectorIndex(idet);
//if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) {
if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
- cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
+ //cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
return 0;
}
- fTrackToFollow.Improve(d/21.82*2.33,fYV,fZV);
+ if (constraint) fTrackToFollow.Improve(d/21.82*2.33,GetY(),GetZ());
#ifdef DEBUG
cout<<"accepted lab="<<c->GetLabel(0)<<' '
//--------------------------------------------------------------------
fN=0;
fDetectors=0;
- for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
}
AliITStrackerV2::AliITSlayer::
fR=r; fPhiOffset=p; fZOffset=z;
fNladders=nl; fNdetectors=nd;
fDetectors=new AliITSdetector[fNladders*fNdetectors];
- for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
fN=0;
fI=0;
fI=i+1;
break;
}
+
return cluster;
}
//This function returns the thickness of this layer
//--------------------------------------------------------------------
//-pi<phi<+pi
- if (3 <fR&&fR<8 ) return 2.2*0.096;
+ if (3 <fR&&fR<8 ) return 2.5*0.096;
if (13<fR&&fR<26) return 1.1*0.088;
if (37<fR&&fR<41) return 1.1*0.085;
return 1.1*0.081;
class AliITStrackerV2 : public AliTracker {
public:
AliITStrackerV2():AliTracker(){}
- AliITStrackerV2(const AliITSgeom *geom, Int_t event=0) throw (const Char_t *);
-
+ AliITStrackerV2(const AliITSgeom *geom,Int_t event=0) throw (const Char_t *);
AliCluster *GetCluster(Int_t index) const;
Int_t Clusters2Tracks(const TFile *in, TFile *out);
Int_t PropagateBack(const TFile *in, TFile *out);
+ void SetupFirstPass(Int_t *flags, Double_t *cuts=0);
+ void SetupSecondPass(Int_t *flags, Double_t *cuts=0);
+
+ class AliITSdetector {
+ public:
+ AliITSdetector(){}
+ AliITSdetector(Double_t r,Double_t phi) {fR=r; fPhi=phi;}
+ void *operator new(size_t s,AliITSdetector *p) {return p;}
+ Double_t GetR() const {return fR;}
+ Double_t GetPhi() const {return fPhi;}
+ private:
+ Double_t fR; // polar coordinates
+ Double_t fPhi; // of this detector
+ };
+
+ class AliITSlayer {
+ public:
+ AliITSlayer();
+ AliITSlayer(Double_t r, Double_t p, Double_t z, Int_t nl, Int_t nd);
+ ~AliITSlayer();
+ Int_t InsertCluster(AliITSclusterV2 *c);
+ void SelectClusters(Double_t zmi,Double_t zma,Double_t ymi,Double_t yma);
+ const AliITSclusterV2 *GetNextCluster(Int_t &ci);
+ void *operator new(size_t s, AliITSlayer *p) {return p;}
+ Double_t GetR() const {return fR;}
+ AliITSclusterV2 *GetCluster(Int_t i) const {return fClusters[i];}
+ AliITSdetector &GetDetector(Int_t n) const { return fDetectors[n]; }
+ Int_t FindDetectorIndex(Double_t phi, Double_t z) const;
+ Double_t GetThickness(Double_t phi, Double_t z) const;
+ Int_t InRoad() const ;
+ Int_t GetNumberOfClusters() const {return fN;}
+ private:
+ Double_t fR; // mean radius of this layer
+ Double_t fPhiOffset; // offset of the first detector in Phi
+ Int_t fNladders; // number of ladders
+ Double_t fZOffset; // offset of the first detector in Z
+ Int_t fNdetectors; // detectors/ladder
+ AliITSdetector *fDetectors; // array of detectors
+ Int_t fN; // number of clusters
+ AliITSclusterV2 *fClusters[kMaxClusterPerLayer]; // pointers to clusters
+ Double_t fZmax; // edges
+ Double_t fYmin; // of the
+ Double_t fYmax; // "window"
+ Int_t fI; // index of the current cluster within the "window"
+ Int_t FindClusterIndex(Double_t z) const;
+ };
private:
- Int_t fEventN; //event number
+ void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
Double_t GetEffectiveThickness(Double_t phi, Double_t z) const;
-
void FollowProlongation();
Int_t TakeNextProlongation();
-
void ResetBestTrack() {
fBestTrack.~AliITStrackV2();
new(&fBestTrack) AliITStrackV2(fTrackToFollow);
}
-
void ResetTrackToFollow(const AliITStrackV2 &t) {
fTrackToFollow.~AliITStrackV2();
new(&fTrackToFollow) AliITStrackV2(t);
}
-
- //The two subclasses are public in order to access one from the another
- public:
-class AliITSdetector {
-private:
- Double_t fR; // polar coordinates
- Double_t fPhi; // of this detector
-
-public:
- AliITSdetector(){}
- AliITSdetector(Double_t r,Double_t phi) {fR=r; fPhi=phi;}
-
- void *operator new(size_t s,AliITSdetector *p) {return p;}
-
- Double_t GetR() const {return fR;}
- Double_t GetPhi() const {return fPhi;}
-};
-
-class AliITSlayer {
- Double_t fR; // mean radius of this layer
- Double_t fPhiOffset; // offset of the first detector in Phi
- Int_t fNladders; // number of ladders
- Double_t fZOffset; // offset of the first detector in Z
- Int_t fNdetectors; // detectors/ladder
- AliITSdetector *fDetectors; // array of detectors
-
- Int_t fN; // number of clusters
- AliITSclusterV2 *fClusters[kMaxClusterPerLayer]; // pointers to clusters
-
- Double_t fZmax; // edges
- Double_t fYmin; // of the
- Double_t fYmax; // "window"
- Int_t fI; // index of the current cluster within the "window"
-
- Int_t FindClusterIndex(Double_t z) const;
-
-public:
- AliITSlayer();
- AliITSlayer(Double_t r, Double_t p, Double_t z, Int_t nl, Int_t nd);
- ~AliITSlayer();
- Int_t InsertCluster(AliITSclusterV2 *c);
- void SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin,Double_t ymax);
- const AliITSclusterV2 *GetNextCluster(Int_t &ci);
-
- void *operator new(size_t s, AliITSlayer *p) {return p;}
-
- Double_t GetR() const {return fR;}
- AliITSclusterV2 *GetCluster(Int_t i) const {return fClusters[i];}
- AliITSdetector &GetDetector(Int_t n) const { return fDetectors[n]; }
- Int_t FindDetectorIndex(Double_t phi, Double_t z) const;
- Double_t GetThickness(Double_t phi, Double_t z) const;
-
- Int_t InRoad() const ;
- Int_t GetNumberOfClusters() const {return fN;}
-};
-
- private:
- Int_t fI; // index of the current layer
+ Int_t fEventN; //event number
+ Int_t fI; // index of the current layer
static AliITSlayer fLayers[kMaxLayer]; // ITS layers
- AliITStrackV2 fTracks[kMaxLayer]; // track estimations at the ITS layers
- AliITStrackV2 fBestTrack; // "best" track
- AliITStrackV2 fTrackToFollow; // followed track
-
- Double_t fYV; // Y-coordinate of the primary vertex
- Double_t fZV; // Z-coordinate of the primary vertex
+ AliITStrackV2 fTracks[kMaxLayer]; // track estimations at the ITS layers
+ AliITStrackV2 fBestTrack; // "best" track
+ AliITStrackV2 fTrackToFollow; // followed track
+ Int_t fPass; // current pass through the data
+ Int_t fConstraint[2]; // constraint flags
};
+inline void AliITStrackerV2::SetupFirstPass(Int_t *flags, Double_t *cuts) {
+ // This function sets up flags and cuts for the first tracking pass
+ //
+ // flags[0] - vertex constaint flag
+ // negative means "skip the pass"
+ // 0 means "no constraint"
+ // positive means "normal constraint"
+
+ fConstraint[0]=flags[0];
+ if (cuts==0) return;
+}
+
+inline void AliITStrackerV2::SetupSecondPass(Int_t *flags, Double_t *cuts) {
+ // This function sets up flags and cuts for the second tracking pass
+ //
+ // flags[0] - vertex constaint flag
+ // negative means "skip the pass"
+ // 0 means "no constraint"
+ // positive means "normal constraint"
+
+ fConstraint[1]=flags[0];
+ if (cuts==0) return;
+}
+
+inline void AliITStrackerV2::CookLabel(AliKalmanTrack *t,Float_t wrong) const {
+ //--------------------------------------------------------------------
+ //This function "cooks" a track label. If label<0, this track is fake.
+ //--------------------------------------------------------------------
+ Int_t tpcLabel=t->GetLabel();
+ if (tpcLabel<0) return;
+ AliTracker::CookLabel(t,wrong);
+ if (tpcLabel != t->GetLabel()) t->SetLabel(-tpcLabel);
+}
#endif
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+ #include <fstream.h>
+
+ #include "TH1.h"
+ #include "TFile.h"
+ #include "TTree.h"
+ #include "TObjArray.h"
+ #include "TStyle.h"
+ #include "TCanvas.h"
+ #include "TLine.h"
+ #include "TText.h"
+ #include "TParticle.h"
+ #include "TStopwatch.h"
+
+ #include "AliRun.h"
+ #include "AliPDG.h"
+ #include "AliV0vertex.h"
+#endif
+
+struct GoodVertex {
+ Int_t nlab,plab;
+ Int_t code;
+ Float_t px,py,pz;
+ Float_t x,y,z;
+};
+Int_t good_vertices(GoodVertex *gt, Int_t max);
+
+Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
+ cerr<<"Doing comparison...\n";
+
+ const Double_t V0window=0.05, V0width=0.015;
+ Double_t V0mass=0.5;
+ switch (code) {
+ case kK0Short: V0mass=0.497672; break;
+ case kLambda0: V0mass=1.115683; break;
+ case kLambda0Bar: V0mass=1.115683; break;
+ default: cerr<<"Invalid PDG code !\n"; return 1;
+ }
+
+ TStopwatch timer;
+
+ /*** Load reconstructed vertices ***/
+ TFile *vf=TFile::Open("AliV0vertices.root");
+ if (!vf->IsOpen()) {cerr<<"Can't open AliV0vertices.root !\n"; return 2;}
+ TObjArray varray(1000);
+ TTree *vTree=(TTree*)vf->Get("TreeV");
+ TBranch *branch=vTree->GetBranch("vertices");
+ Int_t nentr=(Int_t)vTree->GetEntries();
+ for (Int_t i=0; i<nentr; i++) {
+ AliV0vertex *iovertex=new AliV0vertex; branch->SetAddress(&iovertex);
+ vTree->GetEvent(i);
+ varray.AddLast(iovertex);
+ }
+
+ /*** Check if the file with the "good" vertices exists ***/
+ GoodVertex gv[1000];
+ Int_t ngood=0;
+ ifstream in("good_vertices");
+ if (in) {
+ cerr<<"Reading good vertices...\n";
+ while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
+ gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
+ gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
+ ngood++;
+ cerr<<ngood<<'\r';
+ if (ngood==1000) {
+ cerr<<"Too many good vertices !\n";
+ break;
+ }
+ }
+ if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
+ } else {
+ /*** generate a file with the "good" vertices ***/
+ cerr<<"Marking good vertices (this will take a while)...\n";
+ ngood=good_vertices(gv,1000);
+ ofstream out("good_vertices");
+ if (out) {
+ for (Int_t ngd=0; ngd<ngood; ngd++)
+ out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
+ gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
+ gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
+ } else cerr<<"Can not open file (good_vertices) !\n";
+ out.close();
+ }
+
+ vf->Close();
+
+
+ TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution
+ hp->SetXTitle("(mrad)"); hp->SetFillColor(2);
+ TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30);
+ hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013);
+ TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.);
+ hpt->SetXTitle("(%)"); hpt->SetFillColor(2);
+
+ TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res.
+ hx->SetXTitle("(mm)"); hx->SetFillColor(6);
+ TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.); //y res
+ hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013);
+ TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.); //z res.
+ hz->SetXTitle("(mm)"); hz->SetFillColor(6);
+
+
+ Double_t pmin=0.2, pmax=4.2; Int_t nchan=20;
+ TH1F *hgood=new TH1F("hgood","Good Vertices",nchan,pmin,pmax);
+ TH1F *hfound=new TH1F("hfound","Found Vertices",nchan,pmin,pmax);
+ TH1F *hfake=new TH1F("hfake","Fake Vertices",nchan,pmin,pmax);
+ TH1F *hg=new TH1F("hg","Efficiency for Good Vertices",nchan,pmin,pmax);
+ hg->SetLineColor(4); hg->SetLineWidth(2);
+ TH1F *hf=new TH1F("hf","Probability of Fake Vertices",nchan,pmin,pmax);
+ hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
+
+ Double_t mmin=V0mass-V0window, mmax=V0mass+V0window;
+ TH1F *v0s =new TH1F("v0s","V0s Effective Mass",40, mmin, mmax);
+ v0s->SetXTitle("(GeV/c)"); v0s->SetFillColor(6);
+
+
+ Double_t pxg=0.,pyg=0.,ptg=0.;
+ Int_t nlab=-1, plab=-1;
+ Int_t i;
+ for (i=0; i<nentr; i++) {
+ AliV0vertex *vertex=(AliV0vertex*)varray.UncheckedAt(i);
+ nlab=TMath::Abs(vertex->GetNlabel());
+ plab=TMath::Abs(vertex->GetPlabel());
+
+ vertex->ChangeMassHypothesis(code);
+
+ Double_t mass=vertex->GetEffMass();
+ v0s->Fill(mass);
+
+ Int_t j;
+ for (j=0; j<ngood; j++) {
+ if (gv[j].code != vertex->GetPdgCode()) continue;
+ if (gv[j].nlab == nlab)
+ if (gv[j].plab == plab) break;
+ }
+
+ if (TMath::Abs(mass-V0mass)>V0width) continue;
+
+ Double_t px,py,pz; vertex->GetPxPyPz(px,py,pz);
+ Double_t pt=TMath::Sqrt(px*px+py*py);
+
+ if (j==ngood) {
+ hfake->Fill(pt);
+ cerr<<"Fake vertex: ("<<nlab<<","<<plab<<")\n";
+ continue;
+ }
+
+ pxg=gv[j].px; pyg=gv[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px);
+ Double_t lamg=TMath::ATan2(gv[j].pz,ptg), lam=TMath::ATan2(pz,pt);
+ hp->Fill((phi - phig)*1000.);
+ hl->Fill((lam - lamg)*1000.);
+ hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.);
+
+ Double_t x,y,z; vertex->GetXYZ(x,y,z);
+ hx->Fill((x-gv[j].x)*10);
+ hy->Fill((y-gv[j].y)*10);
+ hz->Fill((z-gv[j].z)*10);
+
+ hfound->Fill(ptg);
+ gv[j].nlab=-1;
+
+ }
+ for (i=0; i<ngood; i++) {
+ if (gv[i].code != code) continue;
+ pxg=gv[i].px; pyg=gv[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ hgood->Fill(ptg);
+ nlab=gv[i].nlab; plab=gv[i].plab;
+ if (nlab < 0) continue;
+ cerr<<"Vertex ("<<nlab<<','<<plab<<") has not been found !\n";
+ }
+
+ varray.Delete();
+
+ Stat_t ng=hgood->GetEntries();
+ Stat_t nf=hfound->GetEntries();
+
+ cerr<<"Number of found vertices: "<<nentr<<" ("<<nf<<" in the peak)\n";
+ cerr<<"Number of good vertices: "<<ng<<endl;
+
+ if (ng!=0)
+ cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
+
+ gStyle->SetOptStat(111110);
+ gStyle->SetOptFit(1);
+
+ TCanvas *c1=new TCanvas("c1","",0,0,580,610);
+ c1->Divide(2,2);
+
+ c1->cd(1);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ //hp->Fit("gaus");
+ hp->Draw();
+ hl->Draw("same"); c1->cd();
+
+ c1->cd(2);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ //hpt->Fit("gaus"); c1->cd();
+ hpt->Draw(); c1->cd();
+
+ c1->cd(3);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ //hx->Fit("gaus");
+ hx->Draw();
+ hy->Draw("same"); c1->cd();
+
+ c1->cd(4);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ //hz->Fit("gaus");
+ hz->Draw();
+
+
+ TCanvas *c2=new TCanvas("c2","",600,0,580,610);
+ c2->Divide(1,2);
+
+ c2->cd(1);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
+ hg->Divide(hfound,hgood,1,1.,"b");
+ hf->Divide(hfake,hgood,1,1.,"b");
+ hg->SetMaximum(1.4);
+ hg->SetYTitle("V0 reconstruction efficiency");
+ hg->SetXTitle("Pt (GeV/c)");
+ hg->Draw();
+
+ TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
+ line1->Draw("same");
+ TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
+ line2->Draw("same");
+
+ hf->SetFillColor(1);
+ hf->SetFillStyle(3013);
+ hf->SetLineColor(2);
+ hf->SetLineWidth(2);
+ hf->Draw("histsame");
+ TText *text = new TText(0.461176,0.248448,"Fake vertices");
+ text->SetTextSize(0.05);
+ text->Draw();
+ text = new TText(0.453919,1.11408,"Good vertices");
+ text->SetTextSize(0.05);
+ text->Draw();
+
+
+ c2->cd(2);
+ gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
+ v0s->SetXTitle("(GeV/c)");
+ v0s->Fit("gaus","","",V0mass-V0width,V0mass+V0width);
+ Double_t max=v0s->GetMaximum();
+ TLine *line3 = new TLine(V0mass-V0width,0.,V0mass-V0width,max);
+ line3->Draw("same");
+ TLine *line4 = new TLine(V0mass+V0width,0.,V0mass+V0width,max);
+ line4->Draw("same");
+
+ timer.Stop(); timer.Print();
+
+ return 0;
+}
+
+Int_t good_vertices(GoodVertex *gv, Int_t max) {
+ Int_t nv=0;
+ /*** Get information about the cuts ***/
+ Double_t r2min=0.9*0.9;
+ Double_t r2max=2.9*2.9;
+
+ /*** Get labels of the "good" tracks ***/
+ Double_t dd; Int_t id, label[15000], ngt=0;
+ ifstream in("good_tracks_its");
+ if (!in) {
+ cerr<<"Can't open the file good_tracks_its \n";
+ return nv;
+ }
+ while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) {
+ ngt++;
+ if (ngt>=15000) {
+ cerr<<"Too many good ITS tracks !\n";
+ return nv;
+ }
+ }
+ if (!in.eof()) {
+ cerr<<"Read error (good_tracks_its) !\n";
+ return nv;
+ }
+
+ /*** Get an access to the kinematics ***/
+ if (gAlice) {delete gAlice; gAlice=0;}
+
+ TFile *file=TFile::Open("galice.root");
+ if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
+ if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
+ cerr<<"gAlice has not been found on galice.root !\n";
+ exit(5);
+ }
+
+ Int_t np=gAlice->GetEvent(0);
+ while (np--) {
+ cerr<<np<<'\r';
+ TParticle *p0=gAlice->Particle(np);
+
+ /*** only these V0s are "good" ***/
+ Int_t code=p0->GetPdgCode();
+ if (code!=kK0Short)
+ if (code!=kLambda0)
+ if (code!=kLambda0Bar) continue;
+
+ /*** daughter tracks should be "good" ***/
+ Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
+ if (nlab==plab) continue;
+ if (nlab<0) continue;
+ if (plab<0) continue;
+ Int_t i;
+ for (i=0; i<ngt; i++) if (label[i]==nlab) break;
+ if (i==ngt) continue;
+ for (i=0; i<ngt; i++) if (label[i]==plab) break;
+ if (i==ngt) continue;
+
+ /*** fiducial volume ***/
+ TParticle *p=gAlice->Particle(nlab);
+ Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y;
+ if (r2<r2min) continue;
+ if (r2>r2max) continue;
+
+ gv[nv].code=code;
+ gv[nv].plab=plab; gv[nv].nlab=nlab;
+ gv[nv].px=p0->Px(); gv[nv].py=p0->Py(); gv[nv].pz=p0->Pz();
+ gv[nv].x=x; gv[nv].y=y; gv[nv].z=z;
+ nv++;
+ }
+
+ delete gAlice; gAlice=0;
+
+ file->Close();
+
+ return nv;
+}
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+ #include "AliV0vertexer.h"
+ #include "TFile.h"
+ #include "TStopwatch.h"
+#endif
+
+Int_t AliV0FindVertices() {
+ cerr<<"Looking for V0 vertices...\n";
+
+ TFile *out=TFile::Open("AliV0vertices.root","new");
+ if (!out->IsOpen()) {cerr<<"Delete old AliV0vertices.root !\n"; return 1;}
+
+ TFile *in=TFile::Open("AliITStracksV2.root");
+ if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 2;}
+
+ Double_t cuts[]={33., // max. allowed chi2
+ 0.015,// min. allowed negative daughter's impact parameter
+ 0.015,// min. allowed positive daughter's impact parameter
+ 0.060,// max. allowed DCA between the daughter tracks
+ 0.997,// max. allowed cosine of V0's pointing angle
+ 0.9, // min. radius of the fiducial volume
+ 2.9 // max. radius of the fiducial volume
+ };
+ TStopwatch timer;
+ AliV0vertexer *vertexer=new AliV0vertexer(cuts);
+ Int_t rc=vertexer->Tracks2V0vertices(in,out);
+ delete vertexer;
+ timer.Stop(); timer.Print();
+
+ in->Close();
+ out->Close();
+
+ return rc;
+}
--- /dev/null
+/**************************************************************************
+ * 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 V0 vertex class
+//
+// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <iostream.h>
+#include <TMath.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();
+
+ //Trivial estimation of the vertex parameters
+ Double_t pt, phi, x, par[5];
+ Double_t alpha, cs, sn;
+
+ n.GetExternalParameters(x,par); alpha=n.GetAlpha();
+ pt=1./TMath::Abs(par[4]);
+ phi=TMath::ASin(par[2]) + alpha;
+ Double_t px1=pt*TMath::Cos(phi), py1=pt*TMath::Sin(phi), pz1=pt*par[3];
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x1=x*cs - par[0]*sn;
+ Double_t y1=x*sn + par[0]*cs;
+ Double_t z1=par[1];
+ Double_t sx1=sn*sn*n.GetSigmaY2(), sy1=cs*cs*n.GetSigmaY2();
+
+ p.GetExternalParameters(x,par); alpha=p.GetAlpha();
+ pt=1./TMath::Abs(par[4]);
+ phi=TMath::ASin(par[2]) + alpha;
+ Double_t px2=pt*TMath::Cos(phi), py2=pt*TMath::Sin(phi), pz2=pt*par[3];
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ Double_t x2=x*cs - par[0]*sn;
+ Double_t y2=x*sn + par[0]*cs;
+ Double_t z2=par[1];
+ Double_t sx2=sn*sn*p.GetSigmaY2(), sy2=cs*cs*p.GetSigmaY2();
+
+ Double_t sz1=n.GetSigmaZ2(), sz2=p.GetSigmaZ2();
+ Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
+ Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
+ Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
+ fPos[0]=wx1*x1 + wx2*x2; fPos[1]=wy1*y1 + wy2*y2; fPos[2]=wz1*z1 + wz2*z2;
+
+ //fPos[0]=0.5*(x1+x2); fPos[1]=0.5*(y1+y2); fPos[2]=0.5*(z1+z2);
+ fNmom[0]=px1; fNmom[1]=py1; fNmom[2]=pz1;
+ fPmom[0]=px2; fPmom[1]=py2; fPmom[2]=pz2;
+
+ Double_t e1=TMath::Sqrt(0.13957*0.13957 + px1*px1 + py1*py1 + pz1*pz1);
+ Double_t e2=TMath::Sqrt(0.13957*0.13957 + px2*px2 + py2*py2 + pz2*pz2);
+ fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)-
+ (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
+
+ fChi2=7.;
+}
+
+void AliV0vertex::ChangeMassHypothesis(Int_t code) {
+ //--------------------------------------------------------------------
+ // This function changes the mass hypothesis for this V0
+ //--------------------------------------------------------------------
+ Double_t nmass, pmass;
+
+ switch (code) {
+ case kLambda0:
+ nmass=0.13957; pmass=0.93827; break;
+ case kLambda0Bar:
+ pmass=0.13957; nmass=0.93827; break;
+ case kK0Short:
+ nmass=0.13957; pmass=0.13957; break;
+ default:
+ cerr<<"AliV0vertex::ChangeMassHypothesis: ";
+ cerr<<"invalide PDG code ! Assuming K0s...\n";
+ nmass=0.13957; pmass=0.13957; break;
+ }
+
+ Double_t px1=fNmom[0], py1=fNmom[1], pz1=fNmom[2];
+ Double_t px2=fPmom[0], py2=fPmom[1], pz2=fPmom[2];
+
+ Double_t e1=TMath::Sqrt(nmass*nmass + px1*px1 + py1*py1 + pz1*pz1);
+ Double_t e2=TMath::Sqrt(pmass*pmass + px2*px2 + py2*py2 + pz2*pz2);
+ fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)-
+ (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2));
+
+ fPdgCode=code;
+}
+
+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;
+}
--- /dev/null
+#ifndef ALIV0VERTEX_H
+#define ALIV0VERTEX_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+// V0 Vertex Class
+//
+// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+
+#include <TObject.h>
+#include "AliPDG.h"
+
+class AliITStrackV2;
+
+class AliV0vertex : public TObject {
+public:
+ AliV0vertex();
+ AliV0vertex(const AliITStrackV2 &neg, const AliITStrackV2 &pos);
+
+ void ChangeMassHypothesis(Int_t code=kLambda0);
+
+ 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;
+ 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
+};
+
+#endif
+
+
--- /dev/null
+/**************************************************************************
+ * 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 V0 vertexer class
+//
+// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//-------------------------------------------------------------------------
+#include <TFile.h>
+#include <TTree.h>
+#include <TObjArray.h>
+#include <iostream.h>
+
+#include "AliV0vertex.h"
+#include "AliV0vertexer.h"
+#include "AliITStrackV2.h"
+
+ClassImp(AliV0vertexer)
+
+Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
+ //--------------------------------------------------------------------
+ //This function reconstructs V0 vertices
+ //--------------------------------------------------------------------
+ TFile *in=(TFile*)inp;
+ TDirectory *savedir=gDirectory;
+
+ if (!in->IsOpen()) {
+ cerr<<"AliV0vertexer::Tracks2V0vertices(): ";
+ cerr<<"file with ITS tracks has not been open !\n";
+ return 1;
+ }
+
+ if (!out->IsOpen()) {
+ cerr<<"AliV0vertexer::Tracks2V0vertices(): ";
+ cerr<<"file for V0 vertices has not been open !\n";
+ return 2;
+ }
+
+ in->cd();
+
+ TTree *trkTree=(TTree*)in->Get("TreeT_ITS_0");
+ TBranch *branch=trkTree->GetBranch("tracks");
+ Int_t nentr=(Int_t)trkTree->GetEntries();
+
+ 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++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ branch->SetAddress(&iotrack);
+ trkTree->GetEvent(i);
+ if (iotrack->Get1Pt() > 0.) {nneg++; negtrks.AddLast(iotrack);}
+ else {npos++; postrks.AddLast(iotrack);}
+ }
+
+
+ out->cd();
+ TTree vtxTree("TreeV","Tree with V0 vertices");
+ AliV0vertex *ioVertex=0;
+ vtxTree.Branch("vertices","AliV0vertex",&ioVertex,32000,0);
+
+
+ for (i=0; i<nneg; i++) {
+ if (i%10==0) cerr<<nneg-i<<'\r';
+ AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
+ for (Int_t k=0; k<npos; k++) {
+ AliITStrackV2 *ptrk=(AliITStrackV2 *)postrks.UncheckedAt(k);
+
+ if (TMath::Abs(ntrk->GetD())<fDNmin) continue;
+ if (TMath::Abs(ptrk->GetD())<fDPmin) 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*px+y*py+z*pz)/TMath::Sqrt(p2*(r2+z*z));
+ if (cost<fCPAmax) continue;
+
+ //vertex.ChangeMassHypothesis(); //default is Lambda0
+
+ ioVertex=&vertex; vtxTree.Fill();
+
+ nvtx++;
+ }
+ }
+
+ cerr<<"Number of reconstructed V0 vertices: "<<nvtx<<endl;
+
+ vtxTree.Write();
+
+ negtrks.Delete();
+ postrks.Delete();
+
+ savedir->cd();
+
+ delete trkTree;
+
+ return 0;
+}
+
+
+static void External2Helix(const AliITStrackV2 *t, Double_t helix[6]) {
+ //--------------------------------------------------------------------
+ // External track parameters -> helix parameters
+ //--------------------------------------------------------------------
+ Double_t alpha,x,cs,sn;
+ t->GetExternalParameters(x,helix); alpha=t->GetAlpha();
+
+ cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
+ helix[5]=x*cs - helix[0]*sn; // x0
+ helix[0]=x*sn + helix[0]*cs; // y0
+//helix[1]= // z0
+ helix[2]=TMath::ASin(helix[2]) + alpha; // phi0
+//helix[3]= // tgl
+ helix[4]=helix[4]/t->GetConvConst(); // C
+}
+
+static void Evaluate(const Double_t *h, Double_t t,
+ Double_t r[3], //radius vector
+ Double_t g[3], //first defivatives
+ Double_t gg[3]) //second derivatives
+{
+ //--------------------------------------------------------------------
+ // Calculate position of a point on a track and some derivatives
+ //--------------------------------------------------------------------
+ Double_t phase=h[4]*t+h[2];
+ Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
+
+ r[0] = h[5] + (sn - h[6])/h[4];
+ r[1] = h[0] - (cs - h[7])/h[4];
+ r[2] = h[1] + h[3]*t;
+
+ g[0] = cs; g[1]=sn; g[2]=h[3];
+
+ gg[0]=-h[4]*sn; gg[1]=h[4]*cs; gg[2]=0.;
+}
+
+Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) {
+ //--------------------------------------------------------------------
+ // This function returns the DCA between two tracks
+ // The tracks will be moved to the point of DCA !
+ //--------------------------------------------------------------------
+ Double_t p1[8]; External2Helix(n,p1);
+ p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]);
+ Double_t p2[8]; External2Helix(p,p2);
+ p2[6]=TMath::Sin(p2[2]); p2[7]=TMath::Cos(p2[2]);
+
+
+ Double_t r1[3],g1[3],gg1[3]; Double_t t1=0.;
+ Evaluate(p1,t1,r1,g1,gg1);
+ Double_t r2[3],g2[3],gg2[3]; Double_t t2=0.;
+ Evaluate(p2,t2,r2,g2,gg2);
+
+ Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
+ Double_t dm=dx*dx + dy*dy + dz*dz;
+
+ Int_t max=27;
+ while (max--) {
+ Double_t gt1=-(dx*g1[0] + dy*g1[1] + dz*g1[2]);
+ Double_t gt2=+(dx*g2[0] + dy*g2[1] + dz*g2[2]);
+ Double_t h11=g1[0]*g1[0] - dx*gg1[0] +
+ g1[1]*g1[1] - dy*gg1[1] +
+ g1[2]*g1[2] - dz*gg1[2];
+ Double_t h22=g2[0]*g2[0] + dx*gg2[0] +
+ g2[1]*g2[1] + dy*gg2[1] +
+ g2[2]*g2[2] + dz*gg2[2];
+ Double_t h12=-(g1[0]*g2[0] + g1[1]*g2[1] + g1[2]*g2[2]);
+
+ Double_t det=h11*h22-h12*h12;
+
+ Double_t dt1,dt2;
+ if (TMath::Abs(det)<1.e-33) {
+ //(quasi)singular Hessian
+ dt1=-gt1; dt2=-gt2;
+ } else {
+ dt1=-(gt1*h22 - gt2*h12)/det;
+ dt2=-(h11*gt2 - h12*gt1)/det;
+ }
+
+ if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
+
+ //check delta(phase1) ?
+ //check delta(phase2) ?
+
+ if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
+ if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
+ if ((gt1*gt1+gt2*gt2) > 1.e-4)
+ cerr<<"AliV0vertexer::PropagateToDCA:"
+ " stopped at not a stationary point !\n";
+ Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
+ if (lmb < 0.)
+ cerr<<"AliV0vertexer::PropagateToDCA:"
+ " stopped at not a minimum !\n";
+ break;
+ }
+
+ Double_t dd=dm;
+ for (Int_t div=1 ; ; div*=2) {
+ Evaluate(p1,t1+dt1,r1,g1,gg1);
+ Evaluate(p2,t2+dt2,r2,g2,gg2);
+ dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
+ dd=dx*dx + dy*dy + dz*dz;
+ if (dd<dm) break;
+ dt1*=0.5; dt2*=0.5;
+ if (div>512) {
+ cerr<<"AliV0vertexer::PropagateToDCA: overshoot !\n"; break;
+ }
+ }
+ dm=dd;
+
+ t1+=dt1;
+ t2+=dt2;
+
+ }
+
+ if (max<=0) cerr<<"AliV0vertexer::PropagateToDCA: too many iterations !\n";
+
+ //propagate tracks to the points of DCA
+ Double_t cs=TMath::Cos(n->GetAlpha());
+ Double_t sn=TMath::Sin(n->GetAlpha());
+ Double_t x=r1[0]*cs + r1[1]*sn;
+ if (!n->PropagateTo(x,0.,0.)) {
+ //cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n";
+ return 1.e+33;
+ }
+
+ cs=TMath::Cos(p->GetAlpha());
+ sn=TMath::Sin(p->GetAlpha());
+ x=r2[0]*cs + r2[1]*sn;
+ if (!p->PropagateTo(x,0.,0.)) {
+ //cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n";
+ return 1.e+33;
+ }
+
+ return TMath::Sqrt(dm);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef ALIV0VERTEXER_H
+#define ALIV0VERTEXER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//------------------------------------------------------------------
+// V0 Vertexer Class
+//
+// Origin: Iouri Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch
+//------------------------------------------------------------------
+
+#include "TObject.h"
+
+class TFile;
+class AliITStrackV2;
+
+//_____________________________________________________________________________
+class AliV0vertexer : public TObject {
+public:
+ AliV0vertexer();
+ AliV0vertexer(const Double_t cuts[7]);
+ void SetCuts(const Double_t cuts[7]);
+
+ Int_t Tracks2V0vertices(const TFile *in, TFile *out);
+ Double_t PropagateToDCA(AliITStrackV2 *nt, AliITStrackV2 *pt);
+
+ void GetCuts(Double_t cuts[7]) const;
+
+private:
+ Double_t fChi2max; // maximal allowed chi2
+ Double_t fDNmin; // min. allowed negative daughter's impact parameter
+ Double_t fDPmin; // min. allowed positive daughter's impact parameter
+ Double_t fDCAmax; // maximal allowed DCA between the daughter tracks
+ Double_t fCPAmax; // maximal allowed cosine of V0's pointing angle
+ Double_t fRmin, fRmax; // max & min radii of the fiducial volume
+
+ ClassDef(AliV0vertexer,1) // V0 verterxer
+};
+
+inline AliV0vertexer::AliV0vertexer() {
+ fChi2max=33.;
+ fDNmin=0.015; fDPmin=0.015;
+ fDCAmax=0.01; fCPAmax=0.025;
+ fRmin=0.5; fRmax=2.5;
+}
+
+inline AliV0vertexer::AliV0vertexer(const Double_t cuts[7]) {
+ fChi2max=cuts[0];
+ fDNmin=cuts[1]; fDPmin=cuts[2];
+ fDCAmax=cuts[3]; fCPAmax=cuts[4];
+ fRmin=cuts[5]; fRmax=cuts[6];
+}
+
+inline void AliV0vertexer::SetCuts(const Double_t cuts[7]) {
+ fChi2max=cuts[0];
+ fDNmin=cuts[1]; fDPmin=cuts[2];
+ fDCAmax=cuts[3]; fCPAmax=cuts[4];
+ fRmin=cuts[5]; fRmax=cuts[6];
+}
+
+inline void AliV0vertexer::GetCuts(Double_t cuts[7]) const {
+ cuts[0]=fChi2max;
+ cuts[1]=fDNmin; cuts[2]=fDPmin;
+ cuts[3]=fDCAmax; cuts[4]=fCPAmax;
+ cuts[5]=fRmin; cuts[6]=fRmax;
+}
+
+#endif
+
+
#pragma link C++ class AliITSclusterV2+;
#pragma link C++ class AliITStrackV2+;
#pragma link C++ class AliITStrackerV2+;
-#pragma link C++ class AliITSVertex+;
+#pragma link C++ class AliV0vertex+;
+#pragma link C++ class AliV0vertexer+;
+#pragma link C++ class AliITSVertex+;
#endif
AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\
AliITSvtest.cxx \
AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \
- AliITSVertex.cxx
+ AliITSVertex.cxx \
+ AliV0vertex.cxx AliV0vertexer.cxx
# AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
# Fortran sources
AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\
AliITSvtest.cxx \
AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \
- AliITSVertex.cxx
+ AliITSVertex.cxx \
+ AliV0vertex.cxx AliV0vertexer.cxx
# AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
HDRS:= $(SRCS:.cxx=.h)
class AliKalmanTrack : public TObject {
public:
- AliKalmanTrack() { fLab=-3141593; fChi2=0; fN=0; }
- AliKalmanTrack(const AliKalmanTrack &t) {fLab=t.fLab;fChi2=t.fChi2;fN=t.fN;}
+ AliKalmanTrack() { fLab=-3141593; fChi2=0; fN=0; fMass=0.13957;}
+ AliKalmanTrack(const AliKalmanTrack &t) {
+ fLab=t.fLab; fChi2=t.fChi2; fN=t.fN; fMass=t.fMass;
+ }
virtual ~AliKalmanTrack(){};
void SetLabel(Int_t lab) {fLab=lab;}
Bool_t IsSortable() const {return kTRUE;}
Int_t GetLabel() const {return fLab;}
Double_t GetChi2() const {return fChi2;}
+ Double_t GetMass() const {return fMass;}
Int_t GetNumberOfClusters() const {return fN;}
virtual Int_t GetClusterIndex(Int_t i) const { //reserved for AliTracker
printf("AliKalmanTrack::GetClusterIndex(Int_t i) must be overloaded !\n");
virtual Int_t Compare(const TObject *o) const {return 0;}
- virtual void GetExternalParameters(Double_t &xr, Double_t x[5]) const {}
- virtual void GetExternalCovariance(Double_t cov[15]) const {}
+ virtual void GetExternalParameters(Double_t &xr, Double_t x[5]) const {;}
+ virtual void GetExternalCovariance(Double_t cov[15]) const {;}
- virtual Double_t GetPredictedChi2(const AliCluster *cluster) const {return 0;}
+ virtual Double_t GetPredictedChi2(const AliCluster *cluster) const {return 0.;}
virtual
- Int_t PropagateTo(Double_t xr,Double_t x0,Double_t rho,Double_t pm) {
- return 0;
- }
- virtual Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i) {
- return 0;
- }
+ Int_t PropagateTo(Double_t xr,Double_t x0,Double_t rho) {return 0;}
+ virtual Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i) {return 0;}
static void SetConvConst(Double_t cc) {fConvConst=cc;}
Double_t GetConvConst() const {return fConvConst;}
protected:
void SetChi2(Double_t chi2) {fChi2=chi2;}
+ void SetMass(Double_t mass) {fMass=mass;}
void SetNumberOfClusters(Int_t n) {fN=n;}
private:
Int_t fLab; // track label
Double_t fChi2; // total chi2 value for this track
+ Double_t fMass; // mass hypothesis
Int_t fN; // number of associated clusters
static Double_t fConvConst; //conversion constant cm -> GeV/c
ClassImp(AliTracker)
+Double_t AliTracker::fX;
+Double_t AliTracker::fY;
+Double_t AliTracker::fZ;
+
//__________________________________________________________________________
void AliTracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
//--------------------------------------------------------------------
class AliTracker {
public:
- AliTracker(){}
+ AliTracker() { fX=fY=fZ=0.; }
virtual ~AliTracker(){}
virtual Int_t Clusters2Tracks(const TFile *in, TFile *out)=0;
virtual Int_t PropagateBack(const TFile *in, TFile *out)=0;
+ static void SetVertex(Double_t *xyz) { fX=xyz[0]; fY=xyz[1]; fZ=xyz[2]; }
//protected:
virtual AliCluster *GetCluster(Int_t index) const=0;
virtual void UseClusters(const AliKalmanTrack *t, Int_t from=0) const;
virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
+ Double_t GetX() const {return fX;}
+ Double_t GetY() const {return fY;}
+ Double_t GetZ() const {return fZ;}
+
+private:
+ static Double_t fX; //X-coordinate of the primary vertex
+ static Double_t fY; //Y-coordinate of the primary vertex
+ static Double_t fZ; //Z-coordinate of the primary vertex
ClassDef(AliTracker,1) //abstract tracker
};
+/****************************************************************************
+ * Very important, delicate and rather obscure macro. *
+ * *
+ * Creates list of "trackable" tracks, *
+ * sorts tracks for matching with the ITS, *
+ * calculates efficiency, resolutions etc. *
+ * *
+ * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
+ * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de *
+ ****************************************************************************/
+
#ifndef __CINT__
#include "alles.h"
#include "AliTPCtracker.h"
#endif
struct GoodTrackTPC {
- Int_t fEventN; //event number
Int_t lab;
Int_t code;
Float_t px,py,pz;
Float_t x,y,z;
};
-Int_t AliTPCComparison(Int_t eventn=1) {
- Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn=1);
-
+Int_t AliTPCComparison(Int_t event=0) {
cerr<<"Doing comparison...\n";
- Int_t i;
- gBenchmark->Start("AliTPCComparison");
-
- TFile *cf=TFile::Open("AliTPCclusters.root");
- if (!cf->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 1;}
- AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
- if (!digp) { cerr<<"TPC parameters have not been found !\n"; return 2; }
+ const Int_t MAX=20000;
+ Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event);
+ gBenchmark->Start("AliTPCComparison");
- ///////////
- AliTPCtracker *tracker =0;
- TObjArray tarray(2000);
- AliTPCtrack *iotrack=0;
- Int_t nentr= 0;
- Int_t eventptr[1000];
+ Int_t nentr=0,i=0; TObjArray tarray(MAX);
+ { /*Load tracks*/
TFile *tf=TFile::Open("AliTPCtracks.root");
- TTree *tracktree=0;
- // Load clusters
- eventptr[0]=0;
- eventptr[1]=0;
- char tname[100]; //Y.B.
- for (Int_t event=0;event<eventn; event++){
- cf->cd();
- tracker = new AliTPCtracker(digp,event);
- tracker->LoadInnerSectors();
- tracker->LoadOuterSectors();
- //Y.B. char tname[100];
- if (eventn==-1) {
- sprintf(tname,"TreeT_TPC");
- }
- else {
- sprintf(tname,"TreeT_TPC_%d",event);
- }
+ if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
+
+ char tname[100]; sprintf(tname,"TreeT_TPC_%d",event);
+ TTree *tracktree=(TTree*)tf->Get(tname);
+ if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
- // Load tracks
- if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
- tracktree=(TTree*)tf->Get(tname);
- if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
- TBranch *tbranch=tracktree->GetBranch("tracks");
- Int_t nentr0=(Int_t)tracktree->GetEntries();
- nentr+=nentr0;
- for (i=0; i<nentr0; i++) {
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ nentr=(Int_t)tracktree->GetEntries();
+ AliTPCtrack *iotrack=0;
+ for (i=0; i<nentr; i++) {
iotrack=new AliTPCtrack;
tbranch->SetAddress(&iotrack);
tracktree->GetEvent(i);
- tracker->CookLabel(iotrack,0.1);
tarray.AddLast(iotrack);
- }
- eventptr[event+1] = nentr; //store end of the event
- delete tracker;
- delete tracktree; //Thanks to Mariana Bondila
- }
+ }
+ delete tracktree; //Thanks to Mariana Bondila
tf->Close();
- delete digp; //Thanks to Mariana Bondila
- cf->Close();
-
+ }
-/////////////////////////////////////////////////////////////////////////
- const Int_t MAX=45000;
+ /* Generate a list of "good" tracks */
GoodTrackTPC gt[MAX];
-
Int_t ngood=0;
ifstream in("good_tracks_tpc");
if (in) {
cerr<<"Reading good tracks...\n";
- while (in>>gt[ngood].fEventN>>gt[ngood].lab>>gt[ngood].code>>
+ while (in>>gt[ngood].lab>>gt[ngood].code>>
gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
ngood++;
if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
} else {
cerr<<"Marking good tracks (this will take a while)...\n";
- ngood=good_tracks_tpc(gt,MAX,eventn); //mi change
+ ngood=good_tracks_tpc(gt,MAX, event);
ofstream out("good_tracks_tpc");
if (out) {
for (Int_t ngd=0; ngd<ngood; ngd++)
- out<<gt[ngd].fEventN<<' '<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
+ out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '<<
gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
} else cerr<<"Can not open file (good_tracks_tpc) !\n";
out.close();
-
- cerr<<"Preparing tracks for matching with the ITS...\n";
- tarray.Sort();
- tf=TFile::Open("AliTPCtracks.root","recreate");
- tracktree=new TTree(tname,"Tree with \"forward\" TPC tracks");
- tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
- for (i=0; i<nentr; i++) {
- iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
- tracktree->Fill();
- }
- tracktree->Write();
- tf->Close();
}
- cerr<<"Number of good tracks : "<<ngood<<endl;
+
+
TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60);
hmpt->SetFillColor(6);
- AliTPCtrack *trk=(AliTPCtrack*)tarray.UncheckedAt(0);
- Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
- Double_t pmax=6.0+pmin;
-
- TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
- TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
- TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
- TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
+ TH1F *hgood=new TH1F("hgood","Good tracks",30,0.1,6.1);
+ TH1F *hfound=new TH1F("hfound","Found tracks",30,0.1,6.1);
+ TH1F *hfake=new TH1F("hfake","Fake tracks",30,0.1,6.1);
+ TH1F *hg=new TH1F("hg","",30,0.1,6.1); //efficiency for good tracks
hg->SetLineColor(4); hg->SetLineWidth(2);
- TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax);
+ TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.1,6.1);
hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
- TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,250.);
+ TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,400.);
hep->SetMarkerStyle(8);
hep->SetMarkerSize(0.4);
- Int_t mingood = ngood; //MI change
- Int_t * track_notfound = new Int_t[ngood];
- Int_t itrack_notfound =0;
- Int_t * track_fake = new Int_t[ngood];
- Int_t itrack_fake = 0;
- Int_t * track_multifound = new Int_t[ngood];
- Int_t * track_multifound_n = new Int_t[ngood];
- Int_t itrack_multifound =0;
+ //MI change
+ Int_t mingood=ngood;
+ Int_t track_notfound[MAX], itrack_notfound=0;
+ Int_t track_fake[MAX], itrack_fake=0;
+ Int_t track_multifound[MAX], track_multifound_n[MAX], itrack_multifound=0;
while (ngood--) {
Int_t lab=gt[ngood].lab,tlab=-1;
Float_t ptg=
TMath::Sqrt(gt[ngood].px*gt[ngood].px + gt[ngood].py*gt[ngood].py);
- if (ptg<pmin) continue;
+ if (ptg<1e-33) continue; // for those not crossing 0 pad row
hgood->Fill(ptg);
- Int_t ievent = gt[ngood].fEventN;
AliTPCtrack *track=0;
- // for (i=0; i<nentr; i++) {
- for (i=eventptr[ievent]; i<eventptr[ievent+1]; i++) {
-
+ for (i=0; i<nentr; i++) {
track=(AliTPCtrack*)tarray.UncheckedAt(i);
tlab=track->GetLabel();
if (lab==TMath::Abs(tlab)) break;
}
- // if (i==nentr) {
- if (i==eventptr[ievent+1]) {
-
- //cerr<<"Track "<<lab<<" was not found !\n";
+ if (i==nentr) {
track_notfound[itrack_notfound++]=lab;
continue;
}
Int_t micount=0;
Int_t mi;
AliTPCtrack * mitrack;
- // for (mi=0; mi<nentr; mi++) {
- for (mi=eventptr[ievent]; mi<eventptr[ievent+1]; mi++) {
-
+ for (mi=0; mi<nentr; mi++) {
mitrack=(AliTPCtrack*)tarray.UncheckedAt(mi);
if (lab==TMath::Abs(mitrack->GetLabel())) micount++;
}
if (micount>1) {
- //cout<<"Track no. "<<lab<<" found "<<micount<<" times\n";
track_multifound[itrack_multifound]=lab;
track_multifound_n[itrack_multifound]=micount;
itrack_multifound++;
}
-
- //
- Double_t xk=gt[ngood].x;//digp->GetPadRowRadii(0,0);
+ Double_t xk=gt[ngood].x;
track->PropagateTo(xk);
if (lab==tlab) hfound->Fill(ptg);
else {
track_fake[itrack_fake++]=lab;
hfake->Fill(ptg);
- //cerr<<lab<<" fake\n";
}
Double_t par[5]; track->GetExternalParameters(xk,par);
}
-
- Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<<ng<<endl;
- Stat_t nf=hfound->GetEntries();
- if (ng!=0)
- cerr<<"\n\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
-
- //MI change - addition
- cout<<"Total number of found tracks ="<<nentr<<"\n";
- cout<<"Total number of \"good\" tracks ="<<mingood<<"\n";
-
-
cout<<"\nList of Not found tracks :\n";
- // Int_t i;
for ( i = 0; i< itrack_notfound; i++){
cout<<track_notfound[i]<<"\t";
if ((i%5)==4) cout<<"\n";
}
cout<<"\nList of multiple found tracks :\n";
for ( i=0; i<itrack_multifound; i++) {
- cout<<"id. "<<track_multifound[i]<<" found - "<<track_multifound_n[i]<<"times\n";
+ cout<<"id. "<<track_multifound[i]
+ <<" found - "<<track_multifound_n[i]<<"times\n";
}
- cout<<"\n\n\n";
-
+ Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
+ if (ng!=0) cerr<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
+ cout<<"Total number of found tracks ="<<nentr<<endl;
+ cout<<"Total number of \"good\" tracks ="
+ <<mingood<<" (selected for comparison: "<<ng<<')'<<endl<<endl;
+
+
gStyle->SetOptStat(111110);
gStyle->SetOptFit(1);
hg->SetXTitle("Pt (GeV/c)");
hg->Draw();
- TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
+ TLine *line1 = new TLine(0.1,1.0,6.1,1.0); line1->SetLineStyle(4);
line1->Draw("same");
- TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
+ TLine *line2 = new TLine(0.1,0.9,6.1,0.9); line2->SetLineStyle(4);
line2->Draw("same");
hf->SetFillColor(1);
}
-Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) {
- //eventn - number of events in file
+Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event) {
+ Int_t nt=0;
TFile *file=TFile::Open("galice.root");
if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
exit(5);
}
- // Int_t np=gAlice->GetEvent(0); //MI change
+ Int_t np=gAlice->GetEvent(event);
AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
Int_t ver = TPC->IsVersion();
Int_t nrow_up=digp->GetNRowUp();
Int_t nrows=digp->GetNRowLow()+nrow_up;
Int_t zero=digp->GetZeroSup();
- Int_t gap=Int_t(0.125*nrows);
+ Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
Int_t good_number=Int_t(0.4*nrows);
-
- //MI change
- Int_t nt=0; //reset counter
- char treeName[100]; //declare event identifier
-
- for (Int_t event=0;event<eventn;event++){
-
- Int_t np=gAlice->GetEvent(event);
- Int_t *good=new Int_t[np];
- for (Int_t ii=0; ii<np; ii++) good[ii]=0;
-
-
- //MI change to be possible compile macro
- //definition out of the swith statemnet
- Int_t sectors_by_rows=0;
- TTree *TD=0;
- AliSimDigits da, *digits=&da;
- Int_t *count=0;
- switch (ver) {
- case 1:
- {
- TFile *cf=TFile::Open("AliTPCclusters.root");
- if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
- AliTPCClustersArray *ca=new AliTPCClustersArray;
- ca->Setup(digp);
- ca->SetClusterType("AliTPCcluster");
- //ca->ConnectTree("Segment Tree");
- ca->ConnectTree("TreeC_TPC_0");
- Int_t nrows=Int_t(ca->GetTree()->GetEntries());
- for (Int_t n=0; n<nrows; n++) {
- AliSegmentID *s=ca->LoadEntry(n);
- Int_t sec,row;
- digp->AdjustSectorRow(s->GetID(),sec,row);
- AliTPCClustersRow &clrow = *ca->GetRow(sec,row);
- Int_t ncl=clrow.GetArray()->GetEntriesFast();
- while (ncl--) {
- AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
- Int_t lab=c->GetLabel(0);
- if (lab<0) continue; //noise cluster
- lab=TMath::Abs(lab);
- if (sec>=digp->GetNInnerSector())
- if (row==nrow_up-1 ) good[lab]|=0x1000;
- if (sec>=digp->GetNInnerSector())
- if (row==nrow_up-1-gap) good[lab]|=0x800;
- good[lab]++;
- }
- ca->ClearRow(sec,row);
- }
- cf->Close();
- }
+ Int_t *good=new Int_t[np];
+ Int_t i;
+ for (i=0; i<np; i++) good[i]=0;
+
+
+ switch (ver) {
+ case 1:
+ {
+ char cname[100]; sprintf(cname,"TreeC_TPC_%d",event);
+ TFile *cf=TFile::Open("AliTPCclusters.root");
+ if (!cf->IsOpen()){cerr<<"Can't open AliTPCclusters.root !\n";exit(5);}
+ AliTPCClustersArray *pca=new AliTPCClustersArray, &ca=*pca;
+ ca.Setup(digp);
+ ca.SetClusterType("AliTPCcluster");
+ ca.ConnectTree(cname);
+ Int_t nrows=Int_t(ca.GetTree()->GetEntries());
+ for (Int_t n=0; n<nrows; n++) {
+ AliSegmentID *s=ca.LoadEntry(n);
+ Int_t sec,row;
+ digp->AdjustSectorRow(s->GetID(),sec,row);
+ AliTPCClustersRow &clrow = *ca.GetRow(sec,row);
+ Int_t ncl=clrow.GetArray()->GetEntriesFast();
+ while (ncl--) {
+ AliTPCcluster *c=(AliTPCcluster*)clrow[ncl];
+ Int_t lab=c->GetLabel(0);
+ if (lab<0) continue; //noise cluster
+ lab=TMath::Abs(lab);
+
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1) good[lab]|=0x4000;
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-gap) good[lab]|=0x1000;
+
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-shift) good[lab]|=0x2000;
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-gap-shift) good[lab]|=0x800;
+
+ good[lab]++;
+ }
+ ca.ClearRow(sec,row);
+ }
+ delete pca;
+ cf->Close();
+ }
break;
- case 2:
-
- sprintf(treeName,"TreeD_75x40_100x60_%d",event);
- TD=(TTree*)gDirectory->Get(treeName);
- TD->GetBranch("Segment")->SetAddress(&digits);
- count = new Int_t[np];
- Int_t i;
- for (i=0; i<np; i++) count[i]=0;
- sectors_by_rows=(Int_t)TD->GetEntries();
- for (i=0; i<sectors_by_rows; i++) {
- if (!TD->GetEvent(i)) continue;
- Int_t sec,row;
- digp->AdjustSectorRow(digits->GetID(),sec,row);
- cerr<<sec<<' '<<row<<" \r";
- digits->First();
- do { //Many thanks to J.Chudoba who noticed this
- Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
- Short_t dig = digits->GetDigit(it,ip);
- Int_t idx0=digits->GetTrackID(it,ip,0);
- Int_t idx1=digits->GetTrackID(it,ip,1);
- Int_t idx2=digits->GetTrackID(it,ip,2);
- if (idx0>=0 && dig>=zero) count[idx0]+=1;
- if (idx1>=0 && dig>=zero) count[idx1]+=1;
- if (idx2>=0 && dig>=zero) count[idx2]+=1;
+ case 2:
+ {
+ char dname[100]; sprintf(dname,"TreeD_75x40_100x60_%d",event);
+ TTree *TD=(TTree*)gDirectory->Get(dname);
+ AliSimDigits da, *digits=&da;
+ TD->GetBranch("Segment")->SetAddress(&digits);
+
+ Int_t *count = new Int_t[np];
+ Int_t i;
+ for (i=0; i<np; i++) count[i]=0;
+
+ Int_t sectors_by_rows=(Int_t)TD->GetEntries();
+ for (i=0; i<sectors_by_rows; i++) {
+ if (!TD->GetEvent(i)) continue;
+ Int_t sec,row;
+ digp->AdjustSectorRow(digits->GetID(),sec,row);
+ cerr<<sec<<' '<<row<<" \r";
+ digits->First();
+ do { //Many thanks to J.Chudoba who noticed this
+ Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
+ Short_t dig = digits->GetDigit(it,ip);
+ Int_t idx0=digits->GetTrackID(it,ip,0);
+ Int_t idx1=digits->GetTrackID(it,ip,1);
+ Int_t idx2=digits->GetTrackID(it,ip,2);
+ if (idx0>=0 && dig>=zero) count[idx0]+=1;
+ if (idx1>=0 && dig>=zero) count[idx1]+=1;
+ if (idx2>=0 && dig>=zero) count[idx2]+=1;
} while (digits->Next());
- for (Int_t j=0; j<np; j++) {
- if (count[j]>1) {
- if (sec>=digp->GetNInnerSector())
- if (row==nrow_up-1 ) good[j]|=0x1000;
- if (sec>=digp->GetNInnerSector())
- if (row==nrow_up-1-gap) good[j]|=0x800;
- good[j]++;
- }
- count[j]=0;
- }
- }
- delete[] count;
- delete TD; //Thanks to Mariana Bondila
- break;
- default:
- cerr<<"Invalid TPC version !\n";
- file->Close();
- exit(7);
- }
-
- TTree *TH=gAlice->TreeH();
- Int_t npart=(Int_t)TH->GetEntries();
-
- while (npart--) {
- AliTPChit *hit0=0;
-
- TPC->ResetHits();
- TH->GetEvent(npart);
- AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
- while (hit){
- if (hit->fQ==0.) break;
- hit = (AliTPChit*) TPC->NextHit();
- }
- if (hit) {
- hit0 = new AliTPChit(*hit); //Make copy of hit
- hit = hit0;
- }
- else continue;
- AliTPChit *hit1=(AliTPChit*)TPC->NextHit();
- if (hit1==0) continue;
- if (hit1->fQ != 0.) continue;
- Int_t i=hit->Track();
- TParticle *p = (TParticle*)gAlice->Particle(i);
-
- if (p->GetFirstMother()>=0) continue; //secondary particle
- if (good[i] < 0x1000+0x800+2+good_number) continue;
- if (p->Pt()<0.100) continue;
- if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
-
- gt[nt].lab=i;
- gt[nt].code=p->GetPdgCode();
- //**** px py pz - in global coordinate system, x y z - in local !
- gt[nt].px=hit->X(); gt[nt].py=hit->Y(); gt[nt].pz=hit->Z();
- Float_t cs,sn; digp->AdjustCosSin(hit1->fSector,cs,sn);
- gt[nt].fEventN=event; //MI change
- gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
- gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
- gt[nt].z = hit1->Z();
- nt++;
- if (hit0) delete hit0;
- cerr<<i<<" \r";
- if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+ for (Int_t j=0; j<np; j++) {
+ if (count[j]>1) {
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1 ) good[j]|=0x4000;
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-gap) good[j]|=0x1000;
+
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-shift) good[j]|=0x2000;
+ if (sec>=digp->GetNInnerSector())
+ if (row==nrow_up-1-gap-shift) good[j]|=0x800;
+ good[j]++;
+ }
+ count[j]=0;
+ }
+ }
+ delete[] count;
+ delete TD; //Thanks to Mariana Bondila
}
- delete[] good;
+ break;
+ default:
+ cerr<<"Invalid TPC version !\n";
+ file->Close();
+ exit(7);
}
+
+
+ /** select tracks which are "good" enough **/
+ for (i=0; i<np; i++) {
+ if ((good[i]&0x5000) != 0x5000)
+ if ((good[i]&0x2800) != 0x2800) continue;
+ if ((good[i]&0x7FF ) < good_number) continue;
+
+ TParticle *p = (TParticle*)gAlice->Particle(i);
+
+ if (p->Pt()<0.100) continue;
+ if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
+
+ Int_t j=p->GetFirstMother();
+ if (j>=0) {
+ TParticle *pp = (TParticle*)gAlice->Particle(j);
+ if (pp->GetFirstMother()>=0) continue;//only one decay is allowed
+ }
+
+ gt[nt].lab=i;
+ gt[nt].code=p->GetPdgCode();
+ gt[nt].px=0.; gt[nt].py=0.; gt[nt].pz=0.;
+ gt[nt].x=0.; gt[nt].z=0.; gt[nt].z=0.;
+ nt++;
+ if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+ cerr<<np-i<<" \r";
+ }
+
+ /** check if there is also information at the entrance of the TPC **/
+ TTree *TH=gAlice->TreeH(); np=(Int_t)TH->GetEntries();
+ for (i=0; i<np; i++) {
+ TPC->ResetHits();
+ TH->GetEvent(i);
+ AliTPChit *phit = (AliTPChit*)TPC->FirstHit(-1);
+ for ( ; phit; phit=(AliTPChit*)TPC->NextHit() ) {
+ if (phit->fQ !=0. ) continue;
+
+ Double_t px=phit->X(), py=phit->Y(), pz=phit->Z();
+
+ if ((phit=(AliTPChit*)TPC->NextHit())==0) break;
+ if (phit->fQ != 0.) continue;
+
+ Double_t x=phit->X(), y=phit->Y(), z=phit->Z();
+ if (TMath::Sqrt(x*x+y*y)>90.) continue;
+
+ Int_t j, lab=phit->Track();
+ for (j=0; j<nt; j++) {if (gt[j].lab==lab) break;}
+ if (j==nt) continue;
+
+ // (px,py,pz) - in global coordinate system, (x,y,z) - in local !
+ gt[j].px=px; gt[j].py=py; gt[j].pz=pz;
+ Float_t cs,sn; digp->AdjustCosSin(phit->fSector,cs,sn);
+ gt[j].x = x*cs + y*sn;
+ gt[j].y =-x*sn + y*cs;
+ gt[j].z = z;
+ }
+ cerr<<np-i<<" \r";
+ }
+
+ delete[] good;
+
delete gAlice; gAlice=0;
+
file->Close();
gBenchmark->Stop("AliTPCComparison");
gBenchmark->Show("AliTPCComparison");
#ifndef __CINT__
#include <iostream.h>
+ #include "AliTPCParam.h"
#include "AliTPCtracker.h"
#include "TFile.h"
TStopwatch timer;
+ Int_t rc=0;
for (Int_t i=0;i<eventn;i++){
printf("Processing event %d\n",i);
AliTPCtracker *tracker = new AliTPCtracker(par,i);
- Int_t rc=tracker->Clusters2Tracks(0,out);
+ //Double_t xyz[]={0.,0.,0.}; tracker->SetVertex(xyz); //primary vertex
+ rc=tracker->Clusters2Tracks(0,out);
delete tracker;
}
timer.Stop(); timer.Print();
in->Close();
out->Close();
- return 0;
+ return rc;
}
//-----------------------------------------------------------------
SetLabel(t.GetLabel());
SetChi2(0.);
+ SetMass(t.GetMass());
SetNumberOfClusters(0);
//SetConvConst(t.GetConvConst());
}
//_____________________________________________________________________________
-Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
-{
+Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
//-----------------------------------------------------------------
// This function propagates a track to a reference plane x=xk.
//-----------------------------------------------------------------
//Multiple scattering******************
Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fP0)*(y1-fP0)+(z1-fP1)*(z1-fP1));
Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + pm*pm);
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
//Double_t theta2=1.0259e-6*10*10/20/(beta2*p2)*d*rho;
Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
if (x1 < x2) dE=-dE;
cc=fP4;
- fP4*=(1.- sqrt(p2+pm*pm)/p2*dE);
+ fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
fP2+=fX*(fP4-cc);
return 1;
}
//_____________________________________________________________________________
-Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
+Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho)
{
//-----------------------------------------------------------------
// This function propagates tracks to the "vertex".
Double_t tgf=-fP2/(fP4*fP0 + sqrt(1-c*c));
Double_t snf=tgf/sqrt(1.+ tgf*tgf);
Double_t xv=(fP2+snf)/fP4;
- return PropagateTo(xv,x0,rho,pm);
+ return PropagateTo(xv,x0,rho);
}
//_____________________________________________________________________________
const Double_t cc[15], Double_t xr, Double_t alpha);
AliTPCtrack(const AliKalmanTrack& t, Double_t alpha);
AliTPCtrack(const AliTPCtrack& t);
- Int_t PropagateToVertex(
- Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
+ Int_t PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3);
Int_t Rotate(Double_t angle);
- void SetdEdx(Float_t dedx) {fdEdx=dedx;}
+ void SetdEdx(Double_t dedx) {fdEdx=dedx;}
Double_t GetX() const {return fX;}
Double_t GetAlpha() const {return fAlpha;}
//******************************************************
virtual Double_t GetPredictedChi2(const AliCluster *cluster) const;
- Int_t PropagateTo(Double_t xr,
- Double_t x0=28.94,Double_t rho=0.9e-3,Double_t pm=0.139);
+ Int_t PropagateTo(Double_t xr,Double_t x0=28.94,Double_t rho=0.9e-3);
Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i);
void ResetCovariance();
Double_t fAlpha; // Rotation angle the local (TPC sector)
// coordinate system and the global ALICE one.
- Double_t fdEdx; // dE/dx
+ Double_t fdEdx; // dE/dx
Double_t fP0; // Y-coordinate of a track
Double_t fP1; // Z-coordinate of a track
/*
$Log$
+Revision 1.14 2001/10/21 19:04:55 hristov
+Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
+
Revision 1.13 2001/07/23 12:01:30 hristov
Initialisation added
//_____________________________________________________________________________
AliTPCtracker::AliTPCtracker(const AliTPCParam *par, Int_t eventn):
-fkNIS(par->GetNInnerSector()/2),
-fkNOS(par->GetNOuterSector()/2)
+AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
{
//---------------------------------------------------------------------
// The main TPC tracker constructor
//------------------------------------------------------------------
delete[] fInnerSec;
delete[] fOuterSec;
- delete fSeeds;
+ fSeeds->Delete(); delete fSeeds;
}
//_____________________________________________________________________________
for (Int_t js=0; js < nl+nm+nu; js++) {
const AliTPCcluster *kcl;
Double_t x2, y2, z2;
- Double_t x3=0.,y3=0.;
+ Double_t x3=GetX(), y3=GetY(), z3=GetZ();
if (js<nl) {
const AliTPCRow& kr2=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
y2=xx2*sn+y2*cs;
}
- Double_t zz=z1 - z1/x1*(x1-x2);
+ Double_t zz=z1 - (z1-z3)/(x1-x3)*(x1-x2);
if (TMath::Abs(zz-z2)>5.) continue;
Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
if (TMath::Abs(x[3]) > 1.2) continue;
Double_t a=asin(x[2]);
Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
- if (TMath::Abs(zv)>10.) continue;
+ if (TMath::Abs(zv-z3)>10.) continue;
Double_t sy1=kr1[is]->GetSigmaY2(), sz1=kr1[is]->GetSigmaZ2();
Double_t sy2=kcl->GetSigmaY2(), sz2=kcl->GetSigmaZ2();
- Double_t sy3=100*0.025, sy=0.1, sz=0.1;
+ Double_t sy3=400*3./12., sy=0.1, sz=0.1;
Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
}
out->cd();
- TTree tracktree("TPCf","Tree with TPC tracks");
+
+ char tname[100];
+ sprintf(tname,"TreeT_TPC_%d",fEventN);
+ TTree tracktree(tname,"Tree with TPC tracks");
AliTPCtrack *iotrack=0;
tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
}
delete fSeeds->RemoveAt(i);
}
- UnloadOuterSectors();
+ //UnloadOuterSectors();
//tracking in inner sectors
LoadInnerSectors();
if (FollowProlongation(t)) {
if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
t.CookdEdx();
+ CookLabel(pt,0.1); //For comparison only
iotrack=pt;
tracktree.Fill();
UseClusters(&t,nc);
delete fSeeds->RemoveAt(i);
}
UnloadInnerSectors();
+ UnloadOuterSectors();
- char tname[100];
- if (fEventN==-1) {
- sprintf(tname,"TreeT_TPC");
- }
- else {
- sprintf(tname,"TreeT_TPC_%d",fEventN);
- }
-
-
- tracktree.Write(tname);
+ tracktree.Write();
cerr<<"Number of found tracks : "<<found<<endl;
for (i=nl; i<=nu; i++) dedx += fdEdxSample[i];
dedx /= (nu-nl+1);
SetdEdx(dedx);
+
+ //Very rough PID
+ Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
+
+ if (p<0.6) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
+ SetMass(0.93827); return;
+ }
+
+ if (p<1.2) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ SetMass(0.93827); return;
+ }
+
+ SetMass(0.13957); return;
+
}
Int_t fN; //number of loaded sectors
AliTPCSector *fSectors; //pointer to loaded sectors;
- Int_t fEventN; //event number
+ Int_t fEventN; //event number
AliTPCClustersArray fClustersArray; //array of TPC clusters
- TObjArray *fSeeds; //array of track seeds
+ TObjArray *fSeeds; //array of track seeds
};
#endif
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff