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