Updated version of TOF PID for ESD (Yu.Belikov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Nov 2003 18:36:40 +0000 (18:36 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 13 Nov 2003 18:36:40 +0000 (18:36 +0000)
TOF/AliTOFComparison.C [new file with mode: 0644]
TOF/AliTOFpidESD.cxx
TOF/AliTOFpidESD.h

diff --git a/TOF/AliTOFComparison.C b/TOF/AliTOFComparison.C
new file mode 100644 (file)
index 0000000..dce1890
--- /dev/null
@@ -0,0 +1,346 @@
+/****************************************************************************
+ *      This macro estimates efficiency of matching with the TOF.           *
+ *      TOF "Good" tracks are those originating from the primary vertex,    *
+ *      being "good" in the ITS and having at least one digit in the TOF.   * 
+ *         (To get the list of "good" tracks one should first run           *
+ *          AliTPCComparison.C and AliITSComparisonV2.C macros)             *
+ *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
+ ****************************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+  #include <Riostream.h>
+  #include <fstream.h>
+
+  #include "AliRun.h"
+  #include "AliHeader.h"
+  #include "AliStack.h"
+  #include "AliRunLoader.h"
+  #include "AliLoader.h"
+
+  #include "AliESD.h"
+  #include "AliTOFdigit.h"
+
+  #include "TKey.h"
+  #include "TFile.h"
+  #include "TTree.h"
+  #include "TH1.h"
+  #include "TClonesArray.h"
+  #include "TStyle.h"
+  #include "TCanvas.h"
+  #include "TLine.h"
+  #include "TText.h"
+  #include "TParticle.h"
+#endif
+
+struct GoodTrackTOF {
+  Int_t lab;
+  Int_t code;
+  Float_t px,py,pz;
+  Float_t x,y,z;
+};
+
+extern AliRun *gAlice;
+
+Int_t AliTOFComparison() {
+   Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max);
+
+   cerr<<"Doing comparison...\n";
+
+   if (gAlice) { 
+     delete gAlice->GetRunLoader();
+     delete gAlice;//if everything was OK here it is already NULL
+     gAlice = 0x0;
+   }
+
+   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;   
+
+   {
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (rl == 0x0) {
+      cerr<<"Can not open session"<<endl;
+      return 1;
+   }
+   AliLoader* tofl = rl->GetLoader("TOFLoader");
+   if (tofl == 0x0) {
+      cerr<<"Can not get the TOF loader"<<endl;
+      return 2;
+   }
+   tofl->LoadDigits("read");
+
+   rl->GetEvent(0);
+
+   TTree *tofTree=tofl->TreeD();
+   if (!tofTree) {
+      cerr<<"Can't get the TOF cluster tree !\n";
+      return 3;
+   } 
+   TBranch *branch=tofTree->GetBranch("TOF");
+   if (!branch) { 
+      cerr<<"Can't get the branch with the TOF digits !\n";
+      return 4;
+   }
+   branch->SetAddress(&digits);
+
+   tofTree->GetEvent(0);
+
+   delete rl;
+   }
+
+   Int_t nd=digits->GetEntriesFast();
+   cerr<<"Number of digits: "<<nd<<endl;
+
+   const Int_t MAX=15000;
+   GoodTrackTOF gt[MAX];
+   Int_t ngood=0;
+   ifstream in("good_tracks_tof");
+   if (in) {
+      cerr<<"Reading good tracks...\n";
+      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++;
+         cerr<<ngood<<'\r';
+         if (ngood==MAX) {
+            cerr<<"Too many good tracks !\n";
+            break;
+         }
+      }
+      if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
+   } else {
+      cerr<<"Marking good tracks (this will take a while)...\n";
+      ngood=good_tracks_tof(gt,MAX);
+      ofstream out("good_tracks_tof");
+      if (out) {
+        for (Int_t ngd=0; ngd<ngood; ngd++)
+          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_tof) !\n";
+      out.close();
+   }
+
+   Double_t pmin=0.2;
+   Double_t pmax=4.0;
+
+   TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);    
+   TH1F *hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
+   TH1F *hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
+   TH1F *hgp=new TH1F("hgp","",30,pmin,pmax); //efficiency for good tracks
+   hgp->SetLineColor(4); hgp->SetLineWidth(2);
+   TH1F *hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
+   hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
+
+   TH1F *hgoo=new TH1F("hgoo","Good tracks",30,-1,1);    
+   TH1F *hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
+   TH1F *hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
+   TH1F *hgl=new TH1F("hgl","",30,-1,1); //efficiency for good tracks
+   hgl->SetLineColor(4); hgl->SetLineWidth(2);
+   TH1F *hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
+   hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
+
+   TFile *ef=TFile::Open("AliESDs.root");
+   if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+
+   TIter next(ef->GetListOfKeys());
+   TKey *key=0;
+   Int_t nev=0;
+   while ((key=(TKey*)next())!=0) {
+     cerr<<"Processing event number : "<<nev++<<endl;
+
+     AliESD *event=(AliESD*)key->ReadObj();
+     Int_t ntrk=event->GetNumberOfTracks();
+     cerr<<"Number of ESD tracks : "<<ntrk<<endl; 
+
+     Int_t matched=0;
+     Int_t mismatched=0;
+     for (Int_t i=0; i<ngood; i++) {
+         Int_t lab=gt[i].lab;
+         Double_t pxg=gt[i].px, pyg=gt[i].py, pzg=gt[i].pz;
+         Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+
+        if (ptg<0.1) continue;
+
+         Double_t tgl=pzg/ptg; //tan(lambda)
+
+         if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
+
+         Int_t j;
+        AliESDtrack *t=0;
+         for (j=0; j<ntrk; j++) {
+             AliESDtrack *tt=event->GetTrack(j);
+             if (lab!=TMath::Abs(tt->GetLabel())) continue;
+             t=tt;
+             //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
+             if (tt->GetTOFsignal() < 0) continue;
+             UInt_t idx=tt->GetTOFcluster();
+             if ((Int_t)idx>=nd) {
+              cerr<<"Wrong digit index ! "<<idx<<endl;
+               return 5;
+             }
+             AliTOFdigit *dig=(AliTOFdigit*)digits->UncheckedAt(idx);
+             Int_t *label=dig->GetTracks();
+             if (label[0]!=lab)
+             if (label[1]!=lab)
+              if (label[2]!=lab) {
+                 mismatched++; 
+                 if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); } 
+                 break;
+               }
+             if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
+             matched++;
+             break;
+         }
+         if (j==ntrk) {
+           cerr<<"Not matched: "<<lab<<"   ";
+            if (t) {
+               cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
+                  <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
+                  <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
+                  <<(t->GetStatus()&AliESDtrack::kTIME);
+           } else cerr<<"No ESD track !";
+            cerr<<endl;
+         }
+     }
+
+     cerr<<"Number of good tracks: "<<ngood<<endl;
+     cerr<<"Number of matched tracks: "<<matched<<endl;
+     cerr<<"Number of mismatched tracks: "<<mismatched<<endl;
+     if (ngood!=0) cerr<<"Efficiency: "<<Float_t(matched)/ngood<<endl;
+
+     hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
+     hgp->Divide(hfound,hgood,1,1.,"b");
+     hfp->Divide(hfake,hgood,1,1.,"b");
+     hgp->SetMaximum(1.4);
+     hgp->SetYTitle("Matching efficiency");
+     hgp->SetXTitle("Pt (GeV/c)");
+
+     hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
+     hgl->Divide(hfoun,hgoo,1,1.,"b");
+     hfl->Divide(hfak,hgoo,1,1.,"b");
+     hgl->SetMaximum(1.4);
+     hgl->SetYTitle("Matching efficiency");
+     hgl->SetXTitle("Tan(lambda)");
+
+     TCanvas *c1=new TCanvas("c1","",0,0,600,900);
+     c1->Divide(1,2);
+
+     c1->cd(1);
+
+     hgp->Draw();
+     hfp->Draw("histsame");
+     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");
+
+     c1->cd(2);
+
+     hgl->Draw();
+     hfl->Draw("histsame");
+     TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
+     line3->Draw("same");
+     TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
+     line4->Draw("same");
+
+     break;
+   }
+
+   return 0;
+}
+
+Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max) {
+   ifstream in("good_tracks_its");
+   if (!in) {
+     cerr<<"Can't get good_tracks_its !\n"; exit(11);
+   }   
+
+   AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
+   if (!rl) {
+       cerr<<"Can't start session !\n";
+       exit(1);
+   }
+
+   rl->GetEvent(0);
+
+
+   rl->LoadgAlice();
+   rl->LoadHeader();
+   Int_t np = rl->GetHeader()->GetNtrack();
+
+   Int_t *good=new Int_t[np];
+   Int_t k;
+   for (k=0; k<np; k++) good[k]=0;
+
+   AliLoader* tofl = rl->GetLoader("TOFLoader");
+   if (tofl == 0x0) {
+      cerr<<"Can not get the TOF loader"<<endl;
+      exit(2);
+   }
+   tofl->LoadDigits("read");
+
+   TTree *dTree=tofl->TreeD();
+   if (!dTree) {
+      cerr<<"Can't get the TOF cluster tree !\n";
+      exit(3);
+   } 
+
+   TBranch *branch=dTree->GetBranch("TOF");
+   if (!branch) { 
+     cerr<<"Can't get the branch with the TOF digits !\n";
+     exit(4);
+   }
+   TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
+   branch->SetAddress(&digits);
+   
+  dTree->GetEvent(0);
+  Int_t nd=digits->GetEntriesFast();
+  cerr<<"Number of digits: "<<nd<<endl;
+
+  for (Int_t i=0; i<nd; i++) {
+    AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
+    Int_t l0=d->GetTrack(0);
+       if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
+    Int_t l1=d->GetTrack(1);
+       if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
+    Int_t l2=d->GetTrack(2);
+       if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
+    if (l0>=0) good[l0]++; 
+    if (l1>=0) good[l1]++; 
+    if (l2>=0) good[l2]++;
+  }
+
+   
+  rl->LoadKinematics();
+  AliStack* stack = rl->Stack();
+  Int_t nt=0;
+  Double_t px,py,pz,x,y,z;
+  Int_t code,lab;
+  while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
+      if (good[lab] == 0) continue;
+      TParticle *p = (TParticle*)stack->Particle(lab);
+      if (p == 0x0) {
+         cerr<<"Can not get particle "<<lab<<endl;
+         exit(1);
+      }
+      if (TMath::Abs(p->Vx())>0.1) continue;
+      if (TMath::Abs(p->Vy())>0.1) continue;
+      if (TMath::Abs(p->Vz())>0.1) continue;
+
+      gt[nt].lab=lab;
+      gt[nt].code=p->GetPdgCode();
+//**** px py pz - in global coordinate system
+      gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
+      gt[nt].x=p->Vx(); gt[nt].y=p->Vy(); gt[nt].z=p->Vz();
+      nt++;
+      if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+  }
+
+  delete[] good;
+
+  rl->UnloadKinematics();
+  rl->UnloadHeader();
+  rl->UnloadgAlice();
+  delete rl;
+
+  return nt;
+}
index 2506b790e41c25fcd93f2f24e5a5d12f2b2fbbdc..1e1bf1d9a86f091d20928b8d4f0212ec70258710 100644 (file)
@@ -21,6 +21,7 @@
 #include "TFile.h"
 #include "TTree.h"
 #include "TClonesArray.h"
+#include "TError.h"
 
 #include "AliTOFpidESD.h"
 #include "AliESD.h"
 #include "AliTOFdigit.h"
 
 #include <TGeant3.h>
-#include <TVirtualMC.h>
 #include "../STRUCT/AliBODY.h"
 #include "../STRUCT/AliFRAMEv2.h"
 #include "AliTOFv2FHoles.h"
+#include "AliTOFConstants.h"
 
+#include <stdlib.h>
+
+class TVirtualMC;
+extern TVirtualMC *gMC;
 
 ClassImp(AliTOFpidESD)
 
@@ -40,15 +45,15 @@ static Int_t InitGeo() {
   //gSystem->Load("libgeant321");
   //new TGeant3("C++ Interface to Geant3");
 
-  AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
-  BODY->CreateGeometry();
+  AliBODY *body = new AliBODY("BODY", "Alice envelop");
+  body->CreateGeometry();
 
-  AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
-  FRAME->SetHoles(1);
-  FRAME->CreateGeometry();
+  AliFRAMEv2 *frame = new AliFRAMEv2("FRAME", "Space Frame");
+  frame->SetHoles(1);
+  frame->CreateGeometry();
 
-  AliTOF *TOF = new AliTOFv2FHoles("TOF", "TOF with Holes");
-  TOF->CreateGeometry();
+  AliTOF *tof = new AliTOFv2FHoles("TOF", "TOF with Holes");
+  tof->CreateGeometry();
 
   return 0;
 } 
@@ -58,9 +63,11 @@ AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) {
   //
   //  The main constructor
   //
-  if (InitGeo()) throw "AliTOFpidESD: can not initialize the geomtry !\n";
+  if (InitGeo()) throw "AliTOFpidESD: can not initialize the geometry !\n";
 
-  fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0;
+  fR=378.; 
+  fDy=AliTOFConstants::fgkXPad; fDz=AliTOFConstants::fgkZPad; 
+  fN=0; fEventN=0;
 
   fSigma=param[0];
   fRange=param[1];
@@ -68,17 +75,41 @@ AliTOFpidESD::AliTOFpidESD(Double_t *param) throw (const Char_t *) {
 }
 
 static Int_t DtoM(Int_t *dig, Float_t *g) {
-  const Int_t MAX=13;
-  Int_t lnam[MAX],lnum[MAX];
-
-  extern TVirtualMC *gMC;
+  const Int_t kMAX=13;
+  Int_t lnam[kMAX],lnum[kMAX];
 
   TGeant3 *geant3=(TGeant3*)gMC;
-  if (!geant3) {cerr<<"no geant3 found !\n"; return 1;}
+  if (!geant3) {::Info("DtoM","no geant3 found !"); return 1;}
 
   strncpy((Char_t*)(lnam+0),"ALIC",4); lnum[0]=1;
   strncpy((Char_t*)(lnam+1),"B077",4); lnum[1]=1;
 
+//4 padx  if z<=0 then ...
+if ((dig[1]==4)||(dig[1]==3)) dig[4]=AliTOFConstants::fgkNpadX-dig[4];
+else if (dig[1]==2) {
+   if (dig[2]>7) dig[4]=AliTOFConstants::fgkNpadX-dig[4];
+   else if (dig[2]==7) {
+      if (dig[3]==1) dig[4]=AliTOFConstants::fgkNpadX-dig[4];
+      else dig[4]+=1; 
+   } else dig[4]+=1; 
+} else dig[4]+=1;
+
+//3 padz
+if ((dig[1]==3)||(dig[1]==4)) dig[3]=AliTOFConstants::fgkNpadZ-dig[3];
+else dig[3]+=1;
+
+//2 strip
+if (dig[1]==0) dig[2]=AliTOFConstants::fgkNStripC-dig[2]; 
+else if (dig[1]==1) dig[2]=AliTOFConstants::fgkNStripB-dig[2];
+else dig[2]+=1; 
+
+//1 module
+dig[1]=AliTOFConstants::fgkNPlates-dig[1];
+
+//0 sector
+dig[0]+=(AliTOFConstants::fgkNSectors/2); dig[0]%=AliTOFConstants::fgkNSectors;
+dig[0]+=1;
+
   //sector
   switch (dig[0]) {
      case 1:
@@ -162,7 +193,7 @@ static Int_t DtoM(Int_t *dig, Float_t *g) {
          strncpy((Char_t*)(lnam+3),"BTO1",4); lnum[3]=1;
          break;
      default:
-         cerr<<"Wrong sector number : "<<dig[0]<<" !\n";
+         ::Info("DtoM","Wrong sector number : %d !",dig[0]);
          return 2;
   }
   //module
@@ -194,13 +225,13 @@ static Int_t DtoM(Int_t *dig, Float_t *g) {
         strncpy((Char_t*)(lnam+5),"FLTC",4); lnum[5]=0;
         break;
      default:
-        cerr<<"Wrong module number : "<<dig[1]<<" !\n";
+        ::Info("DtoM","Wrong module number : %d !",dig[1]);
         return 3;
   }
 
   //strip
-  if (dig[2]<1 || dig[2]>19) {
-     cerr<<"Wrong strip number : "<<dig[2]<<" !\n";
+  if (dig[2]<0/*1*/ || dig[2]>19) {
+     ::Info("DtoM","Wrong strip number : %d !",dig[2]);
      return 3;
   }
   strncpy((Char_t*)(lnam+6),"FSTR",4); lnum[6]=dig[2];
@@ -208,14 +239,14 @@ static Int_t DtoM(Int_t *dig, Float_t *g) {
 
   //divz
   if (dig[3]<1 || dig[3]>2) {
-     cerr<<"Wrong z-division number : "<<dig[3]<<" !\n";
+     ::Info("DtoM","Wrong z-division number : %d !",dig[3]);
      return 3;
   }
   strncpy((Char_t*)(lnam+8),"FSEZ",4); lnum[8]=dig[3];
 
   //divx
   if (dig[4]<1 || dig[4]>48) {
-     cerr<<"Wrong x-division number : "<<dig[4]<<" !\n";
+     ::Info("DtoM","Wrong x-division number : %d !",dig[4]);
      return 3;
   }
   strncpy((Char_t*)(lnam+9),"FSEX",4); lnum[9]=dig[4];
@@ -223,7 +254,9 @@ static Int_t DtoM(Int_t *dig, Float_t *g) {
 
   Gcvolu_t *gcvolu=geant3->Gcvolu(); gcvolu->nlevel=0;
   Int_t err=geant3->Glvolu(11,lnam,lnum);  //11-th level
-  if (err) {cout<<err<<" Wrong volume !\n"; return 3;}
+  if (err) {
+    ::Info("DtoM","Wrong volume !"); return 3;
+  }
   Float_t l[3]={0.,0.,0}; geant3->Gdtom(l,g,1);
   return 0;
 }
@@ -273,7 +306,7 @@ Int_t AliTOFpidESD::LoadClusters(const TFile *df) {
     h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; 
     h[3]=d->GetTdc(); h[4]=d->GetAdc();
 
-    AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks());
+    AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
     InsertCluster(cl);
   }  
 
@@ -314,7 +347,7 @@ Int_t AliTOFpidESD::LoadClusters(TTree *dTree) {
     h[1]=TMath::ATan2(g[1],g[0]); h[2]=g[2]; 
     h[3]=d->GetTdc(); h[4]=d->GetAdc();
 
-    AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks());
+    AliTOFcluster *cl=new AliTOFcluster(h,d->GetTracks(),i);
     InsertCluster(cl);
   }  
 
@@ -361,6 +394,16 @@ Int_t AliTOFpidESD::FindClusterIndex(Double_t z) const {
   return m;
 }
 
+static int cmp(const void *p1, const void *p2) {
+  AliESDtrack *t1=*((AliESDtrack**)p1);
+  AliESDtrack *t2=*((AliESDtrack**)p2);
+  Double_t c1[15]; t1->GetExternalCovariance(c1);
+  Double_t c2[15]; t2->GetExternalCovariance(c2);
+  if (c1[0]*c1[2] <c2[0]*c2[2]) return -1;
+  if (c1[0]*c1[2]==c2[0]*c2[2]) return 0;
+  return 1;
+}
+
 //_________________________________________________________________________
 Int_t AliTOFpidESD::MakePID(AliESD *event)
 {
@@ -368,23 +411,31 @@ Int_t AliTOFpidESD::MakePID(AliESD *event)
   //  This function calculates the "detector response" PID probabilities
   //                Just for a bare hint... 
 
-  static const Double_t masses[]={
+  static const Double_t kMasses[]={
     0.000511, 0.105658, 0.139570, 0.493677, 0.938272, 1.875613
   };
 
   Int_t ntrk=event->GetNumberOfTracks();
-  Int_t nmatch=0;
-  for (Int_t i=0; i<ntrk; i++) {
+  AliESDtrack **tracks=new AliESDtrack*[ntrk];
+
+  Int_t i;
+  for (i=0; i<ntrk; i++) {
     AliESDtrack *t=event->GetTrack(i);
+    tracks[i]=t;
+  }
+  qsort(tracks,ntrk,sizeof(AliESDtrack*),cmp);
+
+  Int_t nmatch=0;
+  for (i=0; i<ntrk; i++) {
+    AliESDtrack *t=tracks[i];
 
     if ((t->GetStatus()&AliESDtrack::kTRDout)==0) continue;
-    if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 
     Double_t x,par[5]; t->GetExternalParameters(x,par);
     Double_t cov[15]; t->GetExternalCovariance(cov);
 
-    Double_t dphi2=3*3*(cov[0] + fDy*fDy/12.)/fR/fR;
-    Double_t dz=3*TMath::Sqrt(cov[2] + fDz*fDz/12.);
+Double_t dphi=(5*TMath::Sqrt(cov[0]) + 0.5*fDy + 2.5*TMath::Abs(par[2]))/fR; 
+Double_t dz=5*TMath::Sqrt(cov[2]) + 0.5*fDz + 2.5*TMath::Abs(par[3]);
 
     Double_t phi=TMath::ATan2(par[0],x) + t->GetAlpha();
     if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
@@ -399,11 +450,11 @@ Int_t AliTOFpidESD::MakePID(AliESD *event)
       if (c->GetZ() > z+dz) break;
       if (c->IsUsed()) continue;
 
-      Double_t dphi=TMath::Abs(c->GetPhi()-phi);
-      if (dphi>TMath::Pi()) dphi-=TMath::Pi();
-      if (dphi*dphi > dphi2) continue;
+      Double_t dph=TMath::Abs(c->GetPhi()-phi);
+      if (dph>TMath::Pi()) dph-=TMath::Pi();
+      if (dph>dphi) continue;
 
-      Double_t d2=dphi*dphi*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
+      Double_t d2=dph*dph*fR*fR + (c->GetZ()-z)*(c->GetZ()-z);
       if (d2 > d2max) continue;
 
       d2max=d2;
@@ -418,34 +469,51 @@ Int_t AliTOFpidESD::MakePID(AliESD *event)
     nmatch++;
 
     AliTOFcluster *c=fClusters[index];
+    c->Use();
+
+    Double_t tof=50*c->GetTDC()+32; // in ps
+    t->SetTOFsignal(tof);
+    t->SetTOFcluster(c->GetIndex());
+
+    if ((t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 
     Double_t time[10]; t->GetIntegratedTimes(time);
-    Double_t tof=50*c->GetTDC(); // in ps
+
+    //track length correction
+    Double_t rc=TMath::Sqrt(c->GetR()*c->GetR() + c->GetZ()*c->GetZ());
+    Double_t rt=TMath::Sqrt(x*x + par[0]*par[0] + par[1]*par[1]);
+    Double_t dlt=rc-rt;
 
     Double_t p[10];
     Double_t mom=t->GetP();
     for (Int_t j=0; j<AliESDtrack::kSPECIES; j++) {
-      Double_t mass=masses[j];
-      Double_t dp_p=0.03;      //mean relative pt resolution;
-      Double_t sigma=dp_p*time[j]/(1.+ mom*mom/(mass*mass));
+      Double_t mass=kMasses[j];
+      Double_t dpp=0.01;      //mean relative pt resolution;
+      if (mom>0.5) dpp=0.01*mom;
+      Double_t sigma=dpp*time[j]/(1.+ mom*mom/(mass*mass));
       sigma=TMath::Sqrt(sigma*sigma + fSigma*fSigma);
+
+      time[j]+=dlt/3e-2*TMath::Sqrt(mom*mom+mass*mass)/mom;
+
       if (TMath::Abs(tof-time[j]) > fRange*sigma) {
        p[j]=TMath::Exp(-0.5*fRange*fRange);
         continue;
       }
       p[j]=TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sigma*sigma));
     }
-
-    //if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) continue;
-
-    t->SetTOFsignal(tof);
-    t->SetTOFcluster(index);
+    /*
+    if (c->GetLabel(0)!=TMath::Abs(t->GetLabel())) {
+       cerr<<"Wrong matching: "<<t->GetLabel()<<endl;
+       continue;
+    }
+    */
     t->SetTOFpid(p);
 
-    c->Use();
   }
 
   Info("MakePID","Number of matched ESD track: %d",nmatch);
 
+  delete[] tracks;
+
   return 0;
 }
index 98f3ee67293f52b8e6954aab2eea10c2da738f3f..c4c300bb2dc61137c17d6e16d373cc11e8222ff9 100644 (file)
@@ -15,9 +15,8 @@ class AliESD;
 class TFile;
 class TTree;
 
-const Int_t kMaxCluster=77777; //maximal number of the TOF clusters
-
 class AliTOFpidESD : public TObject {
+enum {kMaxCluster=77777}; //maximal number of the TOF clusters
 public:
   AliTOFpidESD(){fR=376.; fDy=2.5; fDz=3.5; fN=0; fEventN=0;}
   AliTOFpidESD(Double_t *param) throw (const Char_t *);
@@ -34,9 +33,10 @@ public:
 public:
   class AliTOFcluster {
   public:
-    AliTOFcluster(Double_t *h, Int_t *l) {
+    AliTOFcluster(Double_t *h, Int_t *l,Int_t idx) {
       fR=h[0]; fPhi=h[1]; fZ=h[2]; fTDC=h[3]; fADC=h[4];
       fLab[0]=l[0]; fLab[1]=l[1]; fLab[2]=l[2];
+      fIdx=idx;
     }
     void Use() {fADC=-fADC;}
 
@@ -47,17 +47,19 @@ public:
     Double_t GetADC() const {return TMath::Abs(fADC);}
     Int_t IsUsed() const {return (fADC<0) ? 1 : 0;}
     Int_t GetLabel(Int_t n) const {return fLab[n];}
+    Int_t GetIndex() const {return fIdx;}
   private:
-    Int_t fLab[3];
-    Double_t fR;
-    Double_t fPhi;
-    Double_t fZ;
-    Double_t fTDC;
-    Double_t fADC;
+    Int_t fLab[3]; //track labels
+    Double_t fR;   //r-coordinate
+    Double_t fPhi; //phi-coordinate
+    Double_t fZ;   //z-coordinate
+    Double_t fTDC; //TDC count
+    Double_t fADC; //ADC count
+    Int_t fIdx;    //index of this cluster
   };
 
 private:
-  Int_t InsertCluster(AliTOFcluster*);
+  Int_t InsertCluster(AliTOFcluster *c);
   Int_t FindClusterIndex(Double_t z) const;
 
   Int_t fEventN;          //event number