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 <TMath.h>
  12   #include <TError.h>
  13   #include <Riostream.h>
  14   #include <TH1F.h>
  15   #include <TTree.h>
  16   #include <TParticle.h>
  17   #include <TCanvas.h>
  18   #include <TLine.h>
  19   #include <TText.h>
  20   #include <TBenchmark.h>
  21   #include <TStyle.h>
  22   #include <TFile.h>
  23   #include <TROOT.h>
  24
  25   #include "AliStack.h"
  26   #include "AliHeader.h"
  27   #include "AliTrackReference.h"
  28   #include "AliRunLoader.h"
  29   #include "AliRun.h"
  30   #include "AliESDEvent.h"
  31   #include "AliESDtrack.h"
  32
  33   #include "AliTOFcluster.h"
  34   #include "AliLoader.h"
  35
  36   #include "TClonesArray.h"
  37 #endif
  38
  39 Int_t GoodTracksTOF(const Char_t *dir=".");
  40
  41 extern AliRun *gAlice;
  42 extern TBenchmark *gBenchmark;
  43 extern TROOT *gROOT;
  44
  45 static Int_t allgood=0;
  46 static Int_t allmatched=0;
  47 static Int_t allmismatched=0;
  48
  49 Int_t AliTOFComparison(const Char_t *dir=".") {
  50    gBenchmark->Start("AliTOFComparison");
  51
  52    ::Info("AliTOFComparison.C","Doing comparison...");
  53
  54
  55    Double_t pmin=0.2;
  56    Double_t pmax=3.0;
  57
  58    TH1F *hgood=(TH1F*)gROOT->FindObject("hgood");
  59    if (!hgood) hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
  60
  61    TH1F *hfound=(TH1F*)gROOT->FindObject("hfound");
  62    if (!hfound) hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
  63
  64    TH1F *hfake=(TH1F*)gROOT->FindObject("hfake");
  65    if (!hfake) hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
  66
  67    TH1F *hgp=(TH1F*)gROOT->FindObject("hgp");
  68    if (!hgp) hgp=new TH1F("hgp","",30,pmin,pmax);
  69    hgp->SetLineColor(4); hgp->SetLineWidth(2);
  70
  71    TH1F *hfp=(TH1F*)gROOT->FindObject("hfp");
  72    if (!hfp) hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
  73    hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
  74
  75    TH1F *hgoo=(TH1F*)gROOT->FindObject("hgoo");
  76    if (!hgoo) hgoo=new TH1F("hgoo","Good tracks",30,-1,1);
  77
  78    TH1F *hfoun=(TH1F*)gROOT->FindObject("hfoun");
  79    if (!hfoun) hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
  80
  81    TH1F *hfak=(TH1F*)gROOT->FindObject("hfak");
  82    if (!hfak) hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
  83
  84    TH1F *hgl=(TH1F*)gROOT->FindObject("hgl");
  85    if (!hgl) hgl=new TH1F("hgl","",30,-1,1);
  86    hgl->SetLineColor(4); hgl->SetLineWidth(2);
  87
  88    TH1F *hfl=(TH1F*)gROOT->FindObject("hfl");
  89    if (!hfl) hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
  90    hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
  91
  92
  93
  94    Char_t fname[100];
  95    sprintf(fname,"%s/GoodTracksTOF.root",dir);
  96
  97    TFile *refFile=TFile::Open(fname,"old");
  98    if (!refFile || !refFile->IsOpen()) {
  99    ::Info("AliTOFComparison.C","Marking good tracks (will take a while)...");
 100      if (GoodTracksTOF(dir)) {
 101         ::Error("AliTOFComparison.C","Can't generate the reference file !");
 102         return 1;
 103      }
 104    }
 105    refFile=TFile::Open(fname,"old");
 106    if (!refFile || !refFile->IsOpen()) {
 107      ::Error("AliTOFComparison.C","Can't open the reference file !");
 108      return 1;
 109    }
 110
 111    TTree *tofTree=(TTree*)refFile->Get("tofTree");
 112    if (!tofTree) {
 113      ::Error("AliTOFComparison.C","Can't get the reference tree !");
 114      return 2;
 115    }
 116    TBranch *branch=tofTree->GetBranch("TOF");
 117    if (!branch) {
 118      ::Error("AliTOFComparison.C","Can't get the TOF branch !");
 119      return 3;
 120    }
 121    TClonesArray dummy("AliTrackReference",1000), *refs=&dummy;
 122    branch->SetAddress(&refs);
 123
 124
 125    if (gAlice) {
 126      delete AliRunLoader::Instance();
 127      delete gAlice;//if everything was OK here it is already NULL
 128      gAlice = 0x0;
 129    }
 130    sprintf(fname,"%s/galice.root",dir);
 131    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
 132    if (rl == 0x0) {
 133       cerr<<"Can not open session"<<endl;
 134       return 1;
 135    }
 136    AliLoader* tofl = rl->GetLoader("TOFLoader");
 137    if (tofl == 0x0) {
 138       cerr<<"Can not get the TOF loader"<<endl;
 139       return 2;
 140    }
 141    tofl->LoadRecPoints("read");
 142
 143
 144    sprintf(fname,"%s/AliESDs.root",dir);
 145    TFile *ef=TFile::Open(fname);
 146    if ((!ef)||(!ef->IsOpen())) {
 147       ::Error("AliTOFComparison.C","Can't open AliESDs.root !");
 148       delete rl;
 149       return 4;
 150    }
 151    AliESDEvent* event = new AliESDEvent();
 152    TTree* esdTree = (TTree*) ef->Get("esdTree");
 153    if (!esdTree) {
 154       ::Error("AliTOFComparison.C", "no ESD tree found");
 155       return 5;
 156    }
 157    event->ReadFromTree(esdTree);
 158
 159
 160
 161    //******* Loop over events *********
 162    Int_t e=0;
 163    while (esdTree->GetEvent(e)) {
 164      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
 165
 166      rl->GetEvent(e);
 167
 168      TTree *clsTree=tofl->TreeR();
 169      if (!clsTree) {
 170         cerr<<"Can't get the TOF cluster tree !\n";
 171         return 3;
 172      }
 173      TBranch *branch=clsTree->GetBranch("TOF");
 174      if (!branch) {
 175         cerr<<"Can't get the branch with the TOF digits !\n";
 176         return 4;
 177      }
 178      TClonesArray dummy("AliTOFcluster",10000), *clusters=&dummy;
 179      branch->SetAddress(&clusters);
 180
 181      clsTree->GetEvent(0);
 182      Int_t nd=clusters->GetEntriesFast();
 183      cerr<<"Number of the TOF clusters: "<<nd<<endl;
 184
 185
 186
 187      Int_t ntrk=event->GetNumberOfTracks();
 188      cerr<<"Number of ESD tracks : "<<ntrk<<endl;
 189
 190
 191      if (tofTree->GetEvent(e++)==0) {
 192         cerr<<"No reconstructable tracks !\n";
 193         continue;
 194      }
 195
 196      Int_t ngood=refs->GetEntriesFast();
 197
 198      Int_t matched=0;
 199      Int_t mismatched=0;
 200      for (Int_t i=0; i<ngood; i++) {
 201         AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(i);
 202         Int_t lab=ref->Label();
 203         Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
 204
 205         Double_t tgl=ref->Pz()/ptg; //tan(lambda)
 206
 207          if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
 208
 209          Int_t j;
 210          AliESDtrack *t=0;
 211          for (j=0; j<ntrk; j++) {
 212              AliESDtrack *tt=event->GetTrack(j);
 213              if (lab!=TMath::Abs(tt->GetLabel())) continue;
 214              t=tt;
 215              //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
 216              if (tt->GetTOFsignal() < 0) continue;
 217              UInt_t idx=tt->GetTOFcluster();
 218              if ((Int_t)idx>=nd) {
 219                cerr<<"Wrong cluster index ! "<<idx<<endl;
 220                return 5;
 221              }
 222              AliTOFcluster *cls=(AliTOFcluster*)clusters->At(idx);
 223              if (cls) {
 224                if (cls->GetLabel(0)!=lab)
 225                  if (cls->GetLabel(1)!=lab)
 226                    if (cls->GetLabel(2)!=lab) {
 227                      mismatched++;
 228                      if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); }
 229                      break;
 230                    }
 231                if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
 232                matched++;
 233                break;
 234              }
 235          }
 236          if (j==ntrk) {
 237             cerr<<"Not matched: "<<lab<<"   ";
 238             if (t) {
 239                cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
 240                    <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
 241                    <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
 242                    <<(t->GetStatus()&AliESDtrack::kTIME);
 243             } else cerr<<"No ESD track !";
 244             cerr<<endl;
 245          }
 246      }
 247      cout<<"Number of good tracks: "<<ngood<<endl;
 248      cout<<"Number of matched tracks: "<<matched<<endl;
 249      cout<<"Number of mismatched tracks: "<<mismatched<<endl;
 250
 251      allgood+=ngood; allmatched+=matched; allmismatched+=mismatched;
 252
 253      refs->Clear();
 254    } //***** End of the loop over events
 255
 256    delete event;
 257    delete esdTree;
 258    ef->Close();
 259
 260    delete tofTree;
 261    refFile->Close();
 262
 263    if (allgood!=0) cerr<<"\n\nEfficiency: "<<Float_t(allmatched)/allgood<<endl;
 264    cout<<"Total number of good tracks: "<<allgood<<endl;
 265    cout<<"Total number of matched tracks: "<<allmatched<<endl;
 266    cout<<"Total number of mismatched tracks: "<<allmismatched<<endl;
 267
 268    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
 269    if (!c1) {
 270       c1=new TCanvas("c1","",0,0,600,900);
 271       c1->Divide(1,2);
 272    }
 273    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 274    hgp->Divide(hfound,hgood,1,1.,"b");
 275    hfp->Divide(hfake,hgood,1,1.,"b");
 276    hgp->SetMaximum(1.4);
 277    hgp->SetYTitle("Matching efficiency");
 278    hgp->SetXTitle("Pt (GeV/c)");
 279
 280    hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
 281    hgl->Divide(hfoun,hgoo,1,1.,"b");
 282    hfl->Divide(hfak,hgoo,1,1.,"b");
 283    hgl->SetMaximum(1.4);
 284    hgl->SetYTitle("Matching efficiency");
 285    hgl->SetXTitle("Tan(lambda)");
 286
 287    c1->cd(1);
 288
 289    hgp->Draw();
 290    hfp->Draw("histsame");
 291    TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
 292    line1->Draw("same");
 293    TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
 294    line2->Draw("same");
 295
 296    c1->cd(2);
 297
 298    hgl->Draw();
 299    hfl->Draw("histsame");
 300    TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
 301    line3->Draw("same");
 302    TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
 303    line4->Draw("same");
 304
 305    c1->Update();
 306
 307    TFile fc("AliTOFComparison.root","RECREATE");
 308    c1->Write();
 309    fc.Close();
 310
 311    gBenchmark->Stop("AliTOFComparison");
 312    gBenchmark->Show("AliTOFComparison");
 313
 314    delete rl;
 315    return 0;
 316 }
 317
 318 Int_t GoodTracksTOF(const Char_t *dir) {
 319    if (gAlice) {
 320        delete AliRunLoader::Instance();
 321        delete gAlice;//if everything was OK here it is already NULL
 322        gAlice = 0x0;
 323    }
 324
 325    Char_t fname[100];
 326    sprintf(fname,"%s/galice.root",dir);
 327
 328    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
 329    if (!rl) {
 330       ::Error("GoodTracksTOF","Can't start session !");
 331       return 1;
 332    }
 333
 334    rl->LoadgAlice();
 335    rl->LoadHeader();
 336    rl->LoadKinematics();
 337
 338
 339    AliLoader* tofl = rl->GetLoader("TOFLoader");
 340    if (tofl == 0x0) {
 341       ::Error("GoodTracksTOF","Can not get the TOF loader !");
 342       delete rl;
 343       return 2;
 344    }
 345    tofl->LoadRecPoints("read");
 346
 347    Int_t nev=rl->GetNumberOfEvents();
 348    ::Info("GoodTracksTOF","Number of events : %d\n",nev);
 349
 350    sprintf(fname,"%s/GoodTracksITS.root",dir);
 351    TFile *itsFile=TFile::Open(fname);
 352    if ((!itsFile)||(!itsFile->IsOpen())) {
 353        ::Error("GoodTracksTOF","Can't open the GoodTracksITS.root !");
 354        delete rl;
 355        return 5;
 356    }
 357    TClonesArray dum("AliTrackReference",1000), *itsRefs=&dum;
 358    TTree *itsTree=(TTree*)itsFile->Get("itsTree");
 359    if (!itsTree) {
 360        ::Error("GoodTracksTOF","Can't get the ITS reference tree !");
 361        delete rl;
 362        return 6;
 363    }
 364    TBranch *itsBranch=itsTree->GetBranch("ITS");
 365    if (!itsBranch) {
 366       ::Error("GoodTracksTOF","Can't get the ITS reference branch !");
 367       delete rl;
 368       return 7;
 369    }
 370    itsBranch->SetAddress(&itsRefs);
 371
 372
 373    sprintf(fname,"%s/GoodTracksTOF.root",dir);
 374    TFile *tofFile=TFile::Open(fname,"recreate");
 375    TClonesArray dummy("AliTrackReference",1000), *tofRefs=&dummy;
 376    TTree tofTree("tofTree","Tree with info about the reconstructable TOF tracks");
 377    tofTree.Branch("TOF",&tofRefs);
 378
 379
 380    //********  Loop over generated events
 381    for (Int_t e=0; e<nev; e++) {
 382
 383      rl->GetEvent(e); tofFile->cd();
 384
 385      Int_t np = rl->GetHeader()->GetNtrack();
 386      cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
 387
 388      //******** Fill the "good" masks
 389      Int_t *good=new Int_t[np]; Int_t k; for (k=0; k<np; k++) good[k]=0;
 390
 391      TTree *cTree=tofl->TreeR();
 392      if (!cTree) {
 393         ::Error("GoodTracksTOF","Can't get the TOF cluster tree !");
 394         delete rl;
 395         return 8;
 396      }
 397      TBranch *branch=cTree->GetBranch("TOF");
 398      if (!branch) {
 399         ::Error("GoodTracksTOF","Can't get the branch with the TOF digits !");
 400         return 9;
 401      }
 402      TClonesArray dummy("AliTOFcluster",10000), *clusters=&dummy;
 403      branch->SetAddress(&clusters);
 404
 405      cTree->GetEvent(0);
 406      Int_t nd=clusters->GetEntriesFast();
 407
 408      for (Int_t i=0; i<nd; i++) {
 409        AliTOFcluster *c=(AliTOFcluster*)clusters->UncheckedAt(i);
 410        Int_t l0=c->GetLabel(0);
 411           if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
 412        Int_t l1=c->GetLabel(1);
 413           if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
 414        Int_t l2=c->GetLabel(2);
 415           if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
 416        if (l0>=0) good[l0]++;
 417        if (l1>=0) good[l1]++;
 418        if (l2>=0) good[l2]++;
 419      }
 420      clusters->Clear();
 421
 422
 423      //****** select tracks which are "good" enough
 424      AliStack* stack = rl->Stack();
 425
 426      itsTree->GetEvent(e);
 427      Int_t nk=itsRefs->GetEntriesFast();
 428
 429      Int_t nt=0;
 430      for (k=0; k<nk; k++) {
 431         AliTrackReference *itsRef=(AliTrackReference *)itsRefs->UncheckedAt(k);
 432         Int_t lab=itsRef->Label();
 433         if (good[lab] == 0) continue;
 434         TParticle *p = (TParticle*)stack->Particle(lab);
 435         if (p == 0x0) {
 436            cerr<<"Can not get particle "<<lab<<endl;
 437            continue;
 438         }
 439
 440         if (TMath::Abs(p->Vx())>0.1) continue;
 441         if (TMath::Abs(p->Vy())>0.1) continue;
 442         //if (TMath::Abs(p->Vz())>0.1) continue;
 443
 444         new((*tofRefs)[nt]) AliTrackReference(*itsRef);
 445         nt++;
 446      }
 447      itsRefs->Clear();
 448
 449      tofTree.Fill();
 450      tofRefs->Clear();
 451
 452      delete[] good;
 453
 454    } //*** end of the loop over generated events
 455
 456    tofTree.Write();
 457    tofFile->Close();
 458
 459    delete itsTree;
 460    itsFile->Close();
 461
 462    delete rl;
 463    return 0;
 464 }