TOF/AliTOFComparison.C

   1 /****************************************************************************
   2  *      This macro estimates efficiency of matching with the TOF.           *
   3  *      TOF "Good" tracks are those originating from the primary vertex,    *
   4  *      being "good" in the ITS and having at least one digit in the TOF.   *
   5  *         (To get the list of "good" tracks one should first run           *
   6  *          AliTPCComparison.C and AliITSComparisonV2.C macros)             *
   7  *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
   8  ****************************************************************************/
   9
  10 #if !defined(__CINT__) || defined(__MAKECINT__)
  11   #include <Riostream.h>
  12   #include <fstream.h>
  13
  14   #include "AliRun.h"
  15   #include "AliHeader.h"
  16   #include "AliStack.h"
  17   #include "AliRunLoader.h"
  18   #include "AliLoader.h"
  19
  20   #include "AliESD.h"
  21   #include "AliTOFdigit.h"
  22
  23   #include "TKey.h"
  24   #include "TFile.h"
  25   #include "TTree.h"
  26   #include "TH1.h"
  27   #include "TClonesArray.h"
  28   #include "TStyle.h"
  29   #include "TCanvas.h"
  30   #include "TLine.h"
  31   #include "TText.h"
  32   #include "TParticle.h"
  33 #endif
  34
  35 struct GoodTrackTOF {
  36   Int_t lab;
  37   Int_t code;
  38   Float_t px,py,pz;
  39   Float_t x,y,z;
  40 };
  41
  42 extern AliRun *gAlice;
  43
  44 Int_t AliTOFComparison() {
  45    Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max);
  46
  47    cerr<<"Doing comparison...\n";
  48
  49    if (gAlice) {
  50      delete gAlice->GetRunLoader();
  51      delete gAlice;//if everything was OK here it is already NULL
  52      gAlice = 0x0;
  53    }
  54
  55    TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
  56
  57    {
  58    AliRunLoader *rl = AliRunLoader::Open("galice.root");
  59    if (rl == 0x0) {
  60       cerr<<"Can not open session"<<endl;
  61       return 1;
  62    }
  63    AliLoader* tofl = rl->GetLoader("TOFLoader");
  64    if (tofl == 0x0) {
  65       cerr<<"Can not get the TOF loader"<<endl;
  66       return 2;
  67    }
  68    tofl->LoadDigits("read");
  69
  70    rl->GetEvent(0);
  71
  72    TTree *tofTree=tofl->TreeD();
  73    if (!tofTree) {
  74       cerr<<"Can't get the TOF cluster tree !\n";
  75       return 3;
  76    }
  77    TBranch *branch=tofTree->GetBranch("TOF");
  78    if (!branch) {
  79       cerr<<"Can't get the branch with the TOF digits !\n";
  80       return 4;
  81    }
  82    branch->SetAddress(&digits);
  83
  84    tofTree->GetEvent(0);
  85
  86    delete rl;
  87    }
  88
  89    Int_t nd=digits->GetEntriesFast();
  90    cerr<<"Number of digits: "<<nd<<endl;
  91
  92    const Int_t MAX=15000;
  93    GoodTrackTOF gt[MAX];
  94    Int_t ngood=0;
  95    ifstream in("good_tracks_tof");
  96    if (in) {
  97       cerr<<"Reading good tracks...\n";
  98       while (in>>gt[ngood].lab>>gt[ngood].code>>
  99                  gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
 100                  gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
 101          ngood++;
 102          cerr<<ngood<<'\r';
 103          if (ngood==MAX) {
 104             cerr<<"Too many good tracks !\n";
 105             break;
 106          }
 107       }
 108       if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
 109    } else {
 110       cerr<<"Marking good tracks (this will take a while)...\n";
 111       ngood=good_tracks_tof(gt,MAX);
 112       ofstream out("good_tracks_tof");
 113       if (out) {
 114         for (Int_t ngd=0; ngd<ngood; ngd++)
 115           out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
 116              <<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
 117              <<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
 118       } else cerr<<"Can not open file (good_tracks_tof) !\n";
 119       out.close();
 120    }
 121
 122    Double_t pmin=0.2;
 123    Double_t pmax=4.0;
 124
 125    TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
 126    TH1F *hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
 127    TH1F *hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
 128    TH1F *hgp=new TH1F("hgp","",30,pmin,pmax); //efficiency for good tracks
 129    hgp->SetLineColor(4); hgp->SetLineWidth(2);
 130    TH1F *hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
 131    hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
 132
 133    TH1F *hgoo=new TH1F("hgoo","Good tracks",30,-1,1);
 134    TH1F *hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
 135    TH1F *hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
 136    TH1F *hgl=new TH1F("hgl","",30,-1,1); //efficiency for good tracks
 137    hgl->SetLineColor(4); hgl->SetLineWidth(2);
 138    TH1F *hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
 139    hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
 140
 141    TFile *ef=TFile::Open("AliESDs.root");
 142    if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
 143
 144    TIter next(ef->GetListOfKeys());
 145    TKey *key=0;
 146    Int_t nev=0;
 147    while ((key=(TKey*)next())!=0) {
 148      cerr<<"Processing event number : "<<nev++<<endl;
 149
 150      AliESD *event=(AliESD*)key->ReadObj();
 151      Int_t ntrk=event->GetNumberOfTracks();
 152      cerr<<"Number of ESD tracks : "<<ntrk<<endl;
 153
 154      Int_t matched=0;
 155      Int_t mismatched=0;
 156      for (Int_t i=0; i<ngood; i++) {
 157          Int_t lab=gt[i].lab;
 158          Double_t pxg=gt[i].px, pyg=gt[i].py, pzg=gt[i].pz;
 159          Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
 160
 161          if (ptg<0.1) continue;
 162
 163          Double_t tgl=pzg/ptg; //tan(lambda)
 164
 165          if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
 166
 167          Int_t j;
 168          AliESDtrack *t=0;
 169          for (j=0; j<ntrk; j++) {
 170              AliESDtrack *tt=event->GetTrack(j);
 171              if (lab!=TMath::Abs(tt->GetLabel())) continue;
 172              t=tt;
 173              //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
 174              if (tt->GetTOFsignal() < 0) continue;
 175              UInt_t idx=tt->GetTOFcluster();
 176              if ((Int_t)idx>=nd) {
 177                cerr<<"Wrong digit index ! "<<idx<<endl;
 178                return 5;
 179              }
 180              AliTOFdigit *dig=(AliTOFdigit*)digits->UncheckedAt(idx);
 181              Int_t *label=dig->GetTracks();
 182              if (label[0]!=lab)
 183              if (label[1]!=lab)
 184                if (label[2]!=lab) {
 185                  mismatched++;
 186                  if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); }
 187                  break;
 188                }
 189              if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
 190              matched++;
 191              break;
 192          }
 193          if (j==ntrk) {
 194             cerr<<"Not matched: "<<lab<<"   ";
 195             if (t) {
 196                cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
 197                    <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
 198                    <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
 199                    <<(t->GetStatus()&AliESDtrack::kTIME);
 200             } else cerr<<"No ESD track !";
 201             cerr<<endl;
 202          }
 203      }
 204
 205      cerr<<"Number of good tracks: "<<ngood<<endl;
 206      cerr<<"Number of matched tracks: "<<matched<<endl;
 207      cerr<<"Number of mismatched tracks: "<<mismatched<<endl;
 208      if (ngood!=0) cerr<<"Efficiency: "<<Float_t(matched)/ngood<<endl;
 209
 210      hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 211      hgp->Divide(hfound,hgood,1,1.,"b");
 212      hfp->Divide(hfake,hgood,1,1.,"b");
 213      hgp->SetMaximum(1.4);
 214      hgp->SetYTitle("Matching efficiency");
 215      hgp->SetXTitle("Pt (GeV/c)");
 216
 217      hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
 218      hgl->Divide(hfoun,hgoo,1,1.,"b");
 219      hfl->Divide(hfak,hgoo,1,1.,"b");
 220      hgl->SetMaximum(1.4);
 221      hgl->SetYTitle("Matching efficiency");
 222      hgl->SetXTitle("Tan(lambda)");
 223
 224      TCanvas *c1=new TCanvas("c1","",0,0,600,900);
 225      c1->Divide(1,2);
 226
 227      c1->cd(1);
 228
 229      hgp->Draw();
 230      hfp->Draw("histsame");
 231      TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
 232      line1->Draw("same");
 233      TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
 234      line2->Draw("same");
 235
 236      c1->cd(2);
 237
 238      hgl->Draw();
 239      hfl->Draw("histsame");
 240      TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
 241      line3->Draw("same");
 242      TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
 243      line4->Draw("same");
 244
 245      break;
 246    }
 247
 248    return 0;
 249 }
 250
 251 Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max) {
 252    ifstream in("good_tracks_its");
 253    if (!in) {
 254      cerr<<"Can't get good_tracks_its !\n"; exit(11);
 255    }
 256
 257    AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
 258    if (!rl) {
 259        cerr<<"Can't start session !\n";
 260        exit(1);
 261    }
 262
 263    rl->GetEvent(0);
 264
 265
 266    rl->LoadgAlice();
 267    rl->LoadHeader();
 268    Int_t np = rl->GetHeader()->GetNtrack();
 269
 270    Int_t *good=new Int_t[np];
 271    Int_t k;
 272    for (k=0; k<np; k++) good[k]=0;
 273
 274    AliLoader* tofl = rl->GetLoader("TOFLoader");
 275    if (tofl == 0x0) {
 276       cerr<<"Can not get the TOF loader"<<endl;
 277       exit(2);
 278    }
 279    tofl->LoadDigits("read");
 280
 281    TTree *dTree=tofl->TreeD();
 282    if (!dTree) {
 283       cerr<<"Can't get the TOF cluster tree !\n";
 284       exit(3);
 285    }
 286
 287    TBranch *branch=dTree->GetBranch("TOF");
 288    if (!branch) {
 289      cerr<<"Can't get the branch with the TOF digits !\n";
 290      exit(4);
 291    }
 292    TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
 293    branch->SetAddress(&digits);
 294
 295   dTree->GetEvent(0);
 296   Int_t nd=digits->GetEntriesFast();
 297   cerr<<"Number of digits: "<<nd<<endl;
 298
 299   for (Int_t i=0; i<nd; i++) {
 300     AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
 301     Int_t l0=d->GetTrack(0);
 302        if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
 303     Int_t l1=d->GetTrack(1);
 304        if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
 305     Int_t l2=d->GetTrack(2);
 306        if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
 307     if (l0>=0) good[l0]++;
 308     if (l1>=0) good[l1]++;
 309     if (l2>=0) good[l2]++;
 310   }
 311
 312
 313   rl->LoadKinematics();
 314   AliStack* stack = rl->Stack();
 315   Int_t nt=0;
 316   Double_t px,py,pz,x,y,z;
 317   Int_t code,lab;
 318   while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
 319       if (good[lab] == 0) continue;
 320       TParticle *p = (TParticle*)stack->Particle(lab);
 321       if (p == 0x0) {
 322          cerr<<"Can not get particle "<<lab<<endl;
 323          exit(1);
 324       }
 325       if (TMath::Abs(p->Vx())>0.1) continue;
 326       if (TMath::Abs(p->Vy())>0.1) continue;
 327       if (TMath::Abs(p->Vz())>0.1) continue;
 328
 329       gt[nt].lab=lab;
 330       gt[nt].code=p->GetPdgCode();
 331 //**** px py pz - in global coordinate system
 332       gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
 333       gt[nt].x=p->Vx(); gt[nt].y=p->Vy(); gt[nt].z=p->Vz();
 334       nt++;
 335       if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
 336   }
 337
 338   delete[] good;
 339
 340   rl->UnloadKinematics();
 341   rl->UnloadHeader();
 342   rl->UnloadgAlice();
 343   delete rl;
 344
 345   return nt;
 346 }