From 7f6ddf589c191e48fb1ff2621a7dfa67bfaaa367 Mon Sep 17 00:00:00 2001 From: hristov Date: Thu, 8 Nov 2001 16:36:33 +0000 Subject: [PATCH] Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1) Possibility to pass the primary vertex position to the trackers (both for the TPC and the ITS) 2) Possibility to specify the number of tracking passes together with applying (or not applying) the vertex constraint (ITS only) 3) Possibility to make some use of partial PID provided by the TPC when doing tracking in the ITS (ITS only) 4) V0 reconstruction with a helix minimisation of the DCA. (new macros: AliV0FindVertices.C and AliV0Comparison.C) 4a) ( Consequence of the 4) ) All the efficiencies and resolutions are from now on calculated including *secondary*particles* too. (Don't be surprised by the drop in efficiency etc) --- ITS/AliITSComparisonV2.C | 92 +++---- ITS/AliITSFindClustersV2.C | 8 +- ITS/AliITSFindTracksV2.C | 7 + ITS/AliITSrecoV2.h | 6 +- ITS/AliITStrackV2.cxx | 107 ++++---- ITS/AliITStrackV2.h | 37 +-- ITS/AliITStrackerV2.cxx | 229 ++++++++--------- ITS/AliITStrackerV2.h | 158 ++++++------ ITS/AliV0Comparison.C | 336 +++++++++++++++++++++++++ ITS/AliV0FindVertices.C | 35 +++ ITS/AliV0vertex.cxx | 151 ++++++++++++ ITS/AliV0vertex.h | 54 +++++ ITS/AliV0vertexer.cxx | 283 +++++++++++++++++++++ ITS/AliV0vertexer.h | 70 ++++++ ITS/ITSLinkDef.h | 4 +- ITS/Makefile | 3 +- ITS/libITS.pkg | 3 +- STEER/AliKalmanTrack.h | 23 +- STEER/AliTracker.cxx | 4 + STEER/AliTracker.h | 11 +- TPC/AliTPCComparison.C | 486 +++++++++++++++++-------------------- TPC/AliTPCFindTracks.C | 7 +- TPC/AliTPCtrack.cxx | 12 +- TPC/AliTPCtrack.h | 10 +- TPC/AliTPCtracker.cxx | 53 ++-- TPC/AliTPCtracker.h | 4 +- 26 files changed, 1572 insertions(+), 621 deletions(-) create mode 100644 ITS/AliV0Comparison.C create mode 100644 ITS/AliV0FindVertices.C create mode 100644 ITS/AliV0vertex.cxx create mode 100644 ITS/AliV0vertex.h create mode 100644 ITS/AliV0vertexer.cxx create mode 100644 ITS/AliV0vertexer.h diff --git a/ITS/AliITSComparisonV2.C b/ITS/AliITSComparisonV2.C index ef036cfc491..2f82cc2489e 100644 --- a/ITS/AliITSComparisonV2.C +++ b/ITS/AliITSComparisonV2.C @@ -21,67 +21,51 @@ #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; iSetAddress(&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; kGetClusterIndex(k); - AliITSclusterV2 *c=tracker.GetCluster(index); - cout<GetLabel(0)<<' '<GetLabel(1)<<' '<GetLabel(2)<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; iSetAddress(&iotrack); + tracktree->GetEvent(i); + tarray.AddLast(iotrack); + /*if (itsLabel != 1234) continue; + Int_t nc=iotrack->GetNumberOfClusters(); + for (Int_t k=0; kGetClusterIndex(k); + AliITSclusterV2 *c=tracker.GetCluster(index); + cout<GetLabel(0)<<' '<GetLabel(1)<<' '<GetLabel(2)<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++; @@ -94,11 +78,11 @@ Int_t AliITSComparisonV2() { 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; ngdGetEvent(0); + Int_t np=gAlice->GetEvent(event); Int_t *good=new Int_t[np]; Int_t k; @@ -275,7 +259,8 @@ Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) { 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); } @@ -315,11 +300,10 @@ Int_t good_tracks_its(GoodTrackITS *gt, Int_t max) { } 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 diff --git a/ITS/AliITSFindClustersV2.C b/ITS/AliITSFindClustersV2.C index bc82fa92ded..6f9d1515ff0 100644 --- a/ITS/AliITSFindClustersV2.C +++ b/ITS/AliITSFindClustersV2.C @@ -82,7 +82,6 @@ Int_t AliITSFindClustersV2() { 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); @@ -105,6 +104,8 @@ Int_t AliITSFindClustersV2() { cerr<<"Number of entries: "<Clear(); pTree->GetEvent(i); @@ -117,13 +118,13 @@ Int_t AliITSFindClustersV2() { nclusters+=ncl; for (Int_t j=0; jUncheckedAt(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; @@ -143,7 +144,6 @@ Int_t AliITSFindClustersV2() { new(cl[j]) AliITSclusterV2(lab,lp); } cTree->Fill(); clusters->Delete(); - points->Delete(); } cTree->Write(); diff --git a/ITS/AliITSFindTracksV2.C b/ITS/AliITSFindTracksV2.C index aa3ef1dcbf6..00ac80a4ffa 100644 --- a/ITS/AliITSFindTracksV2.C +++ b/ITS/AliITSFindTracksV2.C @@ -1,5 +1,6 @@ #ifndef __CINT__ #include + #include "AliITSgeom.h" #include "AliITStrackerV2.h" #include "TFile.h" @@ -22,6 +23,12 @@ Int_t AliITSFindTracksV2() { 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(); diff --git a/ITS/AliITSrecoV2.h b/ITS/AliITSrecoV2.h index 93603a7389a..ef117753a9e 100644 --- a/ITS/AliITSrecoV2.h +++ b/ITS/AliITSrecoV2.h @@ -24,9 +24,9 @@ 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; diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index aa7ad013e23..3dfb99f838f 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -30,7 +30,7 @@ ClassImp(AliITStrackV2) -const Int_t kWARN=1; +const Int_t kWARN=5; //____________________________________________________________________________ AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) { @@ -42,7 +42,9 @@ 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(); @@ -88,19 +90,35 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) { Int_t n=GetNumberOfClusters(); for (Int_t i=0; iGet1Pt()); Double_t c =TMath::Abs(Get1Pt()); if (c>co) return 1; else if (cGet1Pt()*t->Get1Pt()); + Double_t bo2=po2/(po2 + t->GetMass()*t->GetMass()); + if (p2*b2>po2*bo2) return -1; + else if (p2*b2= 0.99999) { - Int_t n=GetNumberOfClusters(); - if (n>kWARN) cerr<kWARN) + cerr<4) cerr<kWARN) + cerr<kWARN) cerr<kWARN) + cerr<= 0.99999) { Int_t n=GetNumberOfClusters(); - if (n>kWARN) cerr<kWARN) + cerr<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<11.5) - //if (fP1*fP4<0) {cout<<"fP1*fP4="<=1) {cout<<"fP2="<kWARN) { cout<<" bad determinant "<= 0.99999) { Int_t n=GetNumberOfClusters(); - if (n>kWARN) cerr<kWARN) + cerr<= 0.99999) { Int_t n=GetNumberOfClusters(); - if (n>kWARN) cerr<kWARN) + cerr<GetLabel(0)!=LAB) continue; cout<GetY()<<' '<GetZ()<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(); } @@ -137,6 +131,9 @@ cout<GetY()<<' '<GetZ()<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; iGetEvent(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; iGetEvent(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; iGetLabel(); //save the TPC track label #ifdef DEBUG lbl=tpcLabel; if (TMath::Abs(tpcLabel)!=LAB) continue; cout<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-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: "<cd(); + cerr<<"Number of TPC tracks: "< 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; @@ -567,6 +575,7 @@ Int_t AliITStrackerV2::TakeNextProlongation() { 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]); @@ -585,7 +594,8 @@ Int_t AliITStrackerV2::TakeNextProlongation() { 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); @@ -617,10 +627,10 @@ cout<GetLabel(); + if (tpcLabel<0) return; + AliTracker::CookLabel(t,wrong); + if (tpcLabel != t->GetLabel()) t->SetLabel(-tpcLabel); +} #endif diff --git a/ITS/AliV0Comparison.C b/ITS/AliV0Comparison.C new file mode 100644 index 00000000000..53bc2f31216 --- /dev/null +++ b/ITS/AliV0Comparison.C @@ -0,0 +1,336 @@ +#ifndef __CINT__ + #include + #include + + #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; iSetAddress(&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<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; iGetNlabel()); + plab=TMath::Abs(vertex->GetPlabel()); + + vertex->ChangeMassHypothesis(code); + + Double_t mass=vertex->GetEffMass(); + v0s->Fill(mass); + + Int_t j; + for (j=0; jGetPdgCode()) 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: ("<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; iFill(ptg); + nlab=gv[i].nlab; plab=gv[i].plab; + if (nlab < 0) continue; + cerr<<"Vertex ("<GetEntries(); + Stat_t nf=hfound->GetEntries(); + + cerr<<"Number of found vertices: "<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<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; iParticle(nlab); + Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y; + if (r2r2max) 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; +} diff --git a/ITS/AliV0FindVertices.C b/ITS/AliV0FindVertices.C new file mode 100644 index 00000000000..fa10eb13773 --- /dev/null +++ b/ITS/AliV0FindVertices.C @@ -0,0 +1,35 @@ +#ifndef __CINT__ + #include + #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; +} diff --git a/ITS/AliV0vertex.cxx b/ITS/AliV0vertex.cxx new file mode 100644 index 00000000000..ca305e32c69 --- /dev/null +++ b/ITS/AliV0vertex.cxx @@ -0,0 +1,151 @@ +/************************************************************************** + * 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 +#include + +#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; +} diff --git a/ITS/AliV0vertex.h b/ITS/AliV0vertex.h new file mode 100644 index 00000000000..3de575a1331 --- /dev/null +++ b/ITS/AliV0vertex.h @@ -0,0 +1,54 @@ +#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 +#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 + + diff --git a/ITS/AliV0vertexer.cxx b/ITS/AliV0vertexer.cxx new file mode 100644 index 00000000000..077fb229047 --- /dev/null +++ b/ITS/AliV0vertexer.cxx @@ -0,0 +1,283 @@ +/************************************************************************** + * 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 +#include +#include +#include + +#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; iSetAddress(&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; iGetD())GetD()) 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 (costcd(); + + 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 (dd512) { + 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); +} + + + + + + + + + + + + + + + + diff --git a/ITS/AliV0vertexer.h b/ITS/AliV0vertexer.h new file mode 100644 index 00000000000..014418cfb2a --- /dev/null +++ b/ITS/AliV0vertexer.h @@ -0,0 +1,70 @@ +#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 + + diff --git a/ITS/ITSLinkDef.h b/ITS/ITSLinkDef.h index 30af2ac0cb5..151fa16bb50 100644 --- a/ITS/ITSLinkDef.h +++ b/ITS/ITSLinkDef.h @@ -125,6 +125,8 @@ #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 diff --git a/ITS/Makefile b/ITS/Makefile index ce7913ae180..55bc2f0a67e 100644 --- a/ITS/Makefile +++ b/ITS/Makefile @@ -42,7 +42,8 @@ SRCS = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \ 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 diff --git a/ITS/libITS.pkg b/ITS/libITS.pkg index 1b6506862ce..d9ffec9ef3d 100644 --- a/ITS/libITS.pkg +++ b/ITS/libITS.pkg @@ -31,7 +31,8 @@ SRCS = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \ 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) diff --git a/STEER/AliKalmanTrack.h b/STEER/AliKalmanTrack.h index 6fdcd8b927f..19dc66d1213 100644 --- a/STEER/AliKalmanTrack.h +++ b/STEER/AliKalmanTrack.h @@ -16,14 +16,17 @@ class AliCluster; 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"); @@ -32,28 +35,26 @@ public: 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 diff --git a/STEER/AliTracker.cxx b/STEER/AliTracker.cxx index c2d7f01e04f..8c3b66db752 100644 --- a/STEER/AliTracker.cxx +++ b/STEER/AliTracker.cxx @@ -26,6 +26,10 @@ ClassImp(AliTracker) +Double_t AliTracker::fX; +Double_t AliTracker::fY; +Double_t AliTracker::fZ; + //__________________________________________________________________________ void AliTracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const { //-------------------------------------------------------------------- diff --git a/STEER/AliTracker.h b/STEER/AliTracker.h index 0f13835f197..ce9769f86a9 100644 --- a/STEER/AliTracker.h +++ b/STEER/AliTracker.h @@ -16,15 +16,24 @@ class TFile; 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 }; diff --git a/TPC/AliTPCComparison.C b/TPC/AliTPCComparison.C index 5bef125d55c..056f67f20f7 100644 --- a/TPC/AliTPCComparison.C +++ b/TPC/AliTPCComparison.C @@ -1,87 +1,63 @@ +/**************************************************************************** + * 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;eventcd(); - 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; iGetBranch("tracks"); + nentr=(Int_t)tracktree->GetEntries(); + AliTPCtrack *iotrack=0; + for (i=0; iSetAddress(&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++; @@ -94,29 +70,18 @@ Int_t AliTPCComparison(Int_t eventn=1) { 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; ngdBranch("tracks","AliTPCtrack",&iotrack,32000,0); - for (i=0; iFill(); - } - tracktree->Write(); - tf->Close(); } - cerr<<"Number of good tracks : "<SetFillColor(4); TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4); @@ -125,54 +90,41 @@ Int_t AliTPCComparison(Int_t eventn=1) { 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.4SetMarkerStyle(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 (ptgFill(ptg); - Int_t ievent = gt[ngood].fEventN; AliTPCtrack *track=0; - // for (i=0; iGetLabel(); if (lab==TMath::Abs(tlab)) break; } - // if (i==nentr) { - if (i==eventptr[ievent+1]) { - - //cerr<<"Track "<GetLabel())) micount++; } if (micount>1) { - //cout<<"Track no. "<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<GetExternalParameters(xk,par); @@ -235,19 +181,7 @@ Int_t AliTPCComparison(Int_t eventn=1) { } - - Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<GetEntries(); - if (ng!=0) - cerr<<"\n\n\nIntegral efficiency is about "<GetEntries(), nf=hfound->GetEntries(); + if (ng!=0) cerr<<"\n\nIntegral efficiency is about "<SetOptStat(111110); gStyle->SetOptFit(1); @@ -295,9 +235,9 @@ Int_t AliTPCComparison(Int_t eventn=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); @@ -335,8 +275,8 @@ Int_t AliTPCComparison(Int_t eventn=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);} @@ -346,7 +286,7 @@ Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) { 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(); @@ -359,151 +299,171 @@ Int_t good_tracks_tpc(GoodTrackTPC *gt, Int_t max, Int_t eventn) { 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;eventGetEvent(event); - Int_t *good=new Int_t[np]; - for (Int_t ii=0; iiIsOpen()){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; nLoadEntry(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; iIsOpen()){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; nAdjustSectorRow(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; iGetEntries(); - for (i=0; iGetEvent(i)) continue; - Int_t sec,row; - digp->AdjustSectorRow(digits->GetID(),sec,row); - cerr<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; iGetEntries(); + for (i=0; iGetEvent(i)) continue; + Int_t sec,row; + digp->AdjustSectorRow(digits->GetID(),sec,row); + cerr<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; j1) { - 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<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; iParticle(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<TreeH(); np=(Int_t)TH->GetEntries(); + for (i=0; iResetHits(); + 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; jAdjustCosSin(phit->fSector,cs,sn); + gt[j].x = x*cs + y*sn; + gt[j].y =-x*sn + y*cs; + gt[j].z = z; + } + cerr<Close(); gBenchmark->Stop("AliTPCComparison"); gBenchmark->Show("AliTPCComparison"); diff --git a/TPC/AliTPCFindTracks.C b/TPC/AliTPCFindTracks.C index 7948dd348d1..a9ea84d0d95 100644 --- a/TPC/AliTPCFindTracks.C +++ b/TPC/AliTPCFindTracks.C @@ -1,5 +1,6 @@ #ifndef __CINT__ #include + #include "AliTPCParam.h" #include "AliTPCtracker.h" #include "TFile.h" @@ -20,10 +21,12 @@ Int_t AliTPCFindTracks(Int_t eventn=1) { TStopwatch timer; + Int_t rc=0; for (Int_t i=0;iClusters2Tracks(0,out); + //Double_t xyz[]={0.,0.,0.}; tracker->SetVertex(xyz); //primary vertex + rc=tracker->Clusters2Tracks(0,out); delete tracker; } timer.Stop(); timer.Print(); @@ -33,5 +36,5 @@ Int_t AliTPCFindTracks(Int_t eventn=1) { in->Close(); out->Close(); - return 0; + return rc; } diff --git a/TPC/AliTPCtrack.cxx b/TPC/AliTPCtrack.cxx index 36fd9f6c14a..b2fc9408cbb 100644 --- a/TPC/AliTPCtrack.cxx +++ b/TPC/AliTPCtrack.cxx @@ -60,6 +60,7 @@ AliTPCtrack::AliTPCtrack(const AliKalmanTrack& t,Double_t alpha) { //----------------------------------------------------------------- SetLabel(t.GetLabel()); SetChi2(0.); + SetMass(t.GetMass()); SetNumberOfClusters(0); //SetConvConst(t.GetConvConst()); @@ -170,8 +171,7 @@ Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const } //_____________________________________________________________________________ -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. //----------------------------------------------------------------- @@ -223,7 +223,7 @@ Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) //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; @@ -241,14 +241,14 @@ Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm) 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". @@ -257,7 +257,7 @@ Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm) 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); } //_____________________________________________________________________________ diff --git a/TPC/AliTPCtrack.h b/TPC/AliTPCtrack.h index 1626f60d490..8139b09e813 100644 --- a/TPC/AliTPCtrack.h +++ b/TPC/AliTPCtrack.h @@ -39,10 +39,9 @@ public: 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;} @@ -77,8 +76,7 @@ public: //****************************************************** 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(); @@ -87,7 +85,7 @@ private: 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 diff --git a/TPC/AliTPCtracker.cxx b/TPC/AliTPCtracker.cxx index 7538292ca1f..c7067762f9f 100644 --- a/TPC/AliTPCtracker.cxx +++ b/TPC/AliTPCtracker.cxx @@ -15,6 +15,9 @@ /* $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 @@ -66,8 +69,7 @@ Splitted from AliTPCtracking //_____________________________________________________________________________ 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 @@ -105,7 +107,7 @@ AliTPCtracker::~AliTPCtracker() { //------------------------------------------------------------------ delete[] fInnerSec; delete[] fOuterSec; - delete fSeeds; + fSeeds->Delete(); delete fSeeds; } //_____________________________________________________________________________ @@ -481,7 +483,7 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) { 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 (js5.) continue; Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); @@ -518,11 +520,11 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) { 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; @@ -618,7 +620,10 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) { } 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); @@ -645,7 +650,7 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) { } delete fSeeds->RemoveAt(i); } - UnloadOuterSectors(); + //UnloadOuterSectors(); //tracking in inner sectors LoadInnerSectors(); @@ -666,6 +671,7 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) { 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); @@ -676,17 +682,9 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) { 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 : "<