ITS/AliITSComparisonV2.C

   1 /****************************************************************************
   2  *           Very important, delicate and rather obscure macro.             *
   3  *                                                                          *
   4  *               Creates list of "trackable" tracks,                        *
   5  *             calculates efficiency, resolutions etc.                      *
   6  *  The ESD tracks must be in an appropriate state (kITSin or kITSrefit)    *
   7  *                                                                          *
   8  *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
   9  ****************************************************************************/
  10
  11 #if !defined(__CINT__) || defined(__MAKECINT__)
  12   #include <TMath.h>
  13   #include <TError.h>
  14   #include <Riostream.h>
  15   #include <TH1F.h>
  16   #include <TH2F.h>
  17   #include <TTree.h>
  18   #include <TParticle.h>
  19   #include <TCanvas.h>
  20   #include <TLine.h>
  21   #include <TText.h>
  22   #include <TBenchmark.h>
  23   #include <TStyle.h>
  24   #include <TKey.h>
  25   #include <TROOT.h>
  26
  27   #include "AliStack.h"
  28   #include "AliHeader.h"
  29   #include "AliTrackReference.h"
  30   #include "AliRunLoader.h"
  31   #include "AliRun.h"
  32   #include "AliESD.h"
  33
  34   #include "AliITS.h"
  35   #include "AliITSgeom.h"
  36   #include "AliITStrackV2.h"
  37   #include "AliITSclusterV2.h"
  38   #include "AliITSLoader.h"
  39 #endif
  40
  41 Int_t GoodTracksITS(const Char_t *dir=".");
  42
  43 extern AliRun *gAlice;
  44 extern TBenchmark *gBenchmark;
  45 extern TROOT *gROOT;
  46
  47 static Int_t allgood=0;
  48 static Int_t allselected=0;
  49 static Int_t allfound=0;
  50
  51 Int_t AliITSComparisonV2
  52 (Float_t ptcutl=0.2, Float_t ptcuth=10., const Char_t *dir=".", Float_t ratio=0.0) {
  53    gBenchmark->Start("AliITSComparisonV2");
  54
  55    ::Info("AliITSComparisonV2.C","Doing comparison...");
  56
  57
  58    TH1F *hp=(TH1F*)gROOT->FindObject("hp");
  59    if (!hp) hp=new TH1F("hp","PHI resolution",50,-20.,20.);
  60    hp->SetFillColor(4);
  61
  62    TH1F *hl=(TH1F*)gROOT->FindObject("hl");
  63    if (!hl) hl=new TH1F("hl","LAMBDA resolution",50,-20,20);
  64    hl->SetFillColor(4);
  65
  66    TH1F *hpt=(TH1F*)gROOT->FindObject("hpt");
  67    if (!hpt) hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
  68    hpt->SetFillColor(2);
  69
  70    TH1F *hmpt=(TH1F*)gROOT->FindObject("hmpt");
  71    if (!hmpt)
  72       hmpt=new TH1F("hmpt","Transverse impact parameter",30,-300,300);
  73    hmpt->SetFillColor(6);
  74
  75    TH1F *hz=(TH1F*)gROOT->FindObject("hz");
  76    if (!hz) hz=new TH1F("hz","Longitudinal impact parameter",30,-300,300);
  77
  78
  79
  80    TH1F *hgood=(TH1F*)gROOT->FindObject("hgood");
  81    if (!hgood) hgood=new TH1F("hgood","Good tracks",30,0.2,6.1);
  82
  83    TH1F *hfound=(TH1F*)gROOT->FindObject("hfound");
  84    if (!hfound) hfound=new TH1F("hfound","Found tracks",30,0.2,6.1);
  85
  86    TH1F *hfake=(TH1F*)gROOT->FindObject("hfake");
  87    if (!hfake) hfake=new TH1F("hfake","Fake tracks",30,0.2,6.1);
  88
  89    TH1F *hg=(TH1F*)gROOT->FindObject("hg");
  90    if (!hg) hg=new TH1F("hg","Efficiency for good tracks",30,0.2,6.1);
  91    hg->SetLineColor(4); hg->SetLineWidth(2);
  92
  93    TH1F *hf=(TH1F*)gROOT->FindObject("hf");
  94    if (!hf) hf=new TH1F("hf","Efficiency for fake tracks",30,0.2,6.1);
  95    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
  96
  97    TH1F *he=(TH1F*)gROOT->FindObject("he");
  98    if (!he)
  99       he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
 100
 101    TH2F *hep=(TH2F*)gROOT->FindObject("hep");
 102    if (!hep) hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,400.);
 103    hep->SetMarkerStyle(8);
 104    hep->SetMarkerSize(0.4);
 105
 106
 107    Char_t fname[100];
 108    sprintf(fname,"%s/GoodTracksITS.root",dir);
 109
 110    TFile *refFile=TFile::Open(fname,"old");
 111    if (!refFile || !refFile->IsOpen()) {
 112    ::Info("AliITSComparisonV2.C","Marking good tracks (will take a while)...");
 113      if (GoodTracksITS(dir)) {
 114         ::Error("AliITSComparisonV2.C","Can't generate the reference file !");
 115         return 1;
 116      }
 117    }
 118    refFile=TFile::Open(fname,"old");
 119    if (!refFile || !refFile->IsOpen()) {
 120      ::Error("AliITSComparisonV2.C","Can't open the reference file !");
 121      return 1;
 122    }
 123
 124    TTree *itsTree=(TTree*)refFile->Get("itsTree");
 125    if (!itsTree) {
 126      ::Error("AliITSComparisonV2.C","Can't get the reference tree !");
 127      return 2;
 128    }
 129    TBranch *branch=itsTree->GetBranch("ITS");
 130    if (!branch) {
 131      ::Error("AliITSComparisonV2.C","Can't get the ITS branch !");
 132      return 3;
 133    }
 134    TClonesArray dummy("AliTrackReference",1000), *refs=&dummy;
 135    branch->SetAddress(&refs);
 136
 137
 138    sprintf(fname,"%s/AliESDs.root",dir);
 139    TFile *ef=TFile::Open(fname);
 140    if ((!ef)||(!ef->IsOpen())) {
 141       sprintf(fname,"%s/AliESDits.root",dir);
 142       ef=TFile::Open(fname);
 143       if ((!ef)||(!ef->IsOpen())) {
 144          ::Error("AliITSComparisonV2.C","Can't open AliESDits.root !");
 145          return 4;
 146       }
 147    }
 148    TKey *key=0;
 149    TIter next(ef->GetListOfKeys());
 150
 151
 152
 153    //******* Loop over events *********
 154    Int_t e=0;
 155    while ((key=(TKey*)next())!=0) {
 156      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
 157
 158      AliESD *event=(AliESD*)key->ReadObj();
 159
 160      Float_t field=event->GetMagneticField();
 161      AliKalmanTrack::SetConvConst(1000/0.299792458/field);
 162
 163      Int_t nentr=event->GetNumberOfTracks();
 164      allfound+=nentr;
 165
 166      if (itsTree->GetEvent(e++)==0) {
 167         cerr<<"No reconstructable tracks !\n";
 168         continue;
 169      }
 170
 171      Int_t ngood=refs->GetEntriesFast();
 172      allgood+=ngood;
 173
 174      const Int_t MAX=15000;
 175      Int_t notf[MAX], nnotf=0;
 176      Int_t fake[MAX], nfake=0;
 177      Int_t mult[MAX], numb[MAX], nmult=0;
 178      Int_t k;
 179      for (k=0; k<ngood; k++) {
 180         AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(k);
 181         Int_t lab=ref->Label(), tlab=-1;
 182         Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
 183
 184         if (ptg<ptcutl) continue;
 185         if (ptg>ptcuth) continue;
 186
 187         allselected++;
 188
 189         hgood->Fill(ptg);
 190
 191         AliESDtrack *esd=0;
 192         Int_t cnt=0;
 193         Int_t pipe=0;
 194         for (Int_t i=0; i<nentr; i++) {
 195            AliESDtrack *t=event->GetTrack(i);
 196            UInt_t status=t->GetStatus();
 197            UInt_t flags=AliESDtrack::kTPCin|AliESDtrack::kITSin;
 198
 199            if ((status&AliESDtrack::kITSrefit)==0) continue;
 200
 201            if ((status&flags)==status) pipe=1;
 202
 203            Int_t lbl=t->GetLabel();
 204            if (lab==TMath::Abs(lbl)) {
 205               if (cnt==0) {esd=t; tlab=lbl;}
 206               cnt++;
 207            }
 208         }
 209         if (cnt==0) {
 210            notf[nnotf++]=lab;
 211            continue;
 212         } else if (cnt>1){
 213            mult[nmult]=lab;
 214            numb[nmult]=cnt; nmult++;
 215         }
 216
 217         if (esd->GetITSFakeRatio()<ratio) tlab = TMath::Abs(tlab);
 218
 219         if (lab==tlab) hfound->Fill(ptg);
 220         else {
 221           fake[nfake++]=lab;
 222           hfake->Fill(ptg);
 223         }
 224
 225         AliITStrackV2 track(*esd);
 226         if (pipe!=0) {
 227            track.PropagateTo(3.,0.0028,65.19);
 228            track.PropagateToVertex();  // This is not "exactly" the vertex
 229         }
 230
 231         Double_t xv,par[5]; track.GetExternalParameters(xv,par);
 232         Float_t phi=TMath::ASin(par[2]) + track.GetAlpha();
 233         if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 234         if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 235         Float_t lam=TMath::ATan(par[3]);
 236         Float_t pt_1=TMath::Abs(par[4]);
 237
 238         Float_t phig=TMath::ATan2(ref->Py(),ref->Px());
 239         hp->Fill((phi - phig)*1000.);
 240
 241         Float_t lamg=TMath::ATan2(ref->Pz(),ptg);
 242         hl->Fill((lam - lamg)*1000.);
 243
 244         Double_t d=10000*track.GetD(ref->X(),ref->Y());
 245         hmpt->Fill(d);
 246
 247         Double_t z=10000*(track.GetZ()-ref->Z());
 248         hz->Fill(z);
 249
 250         hpt->Fill((pt_1 - 1/ptg)/(1/ptg)*100.);
 251
 252         Float_t mom=1./(pt_1*TMath::Cos(lam));
 253         Float_t dedx=track.GetdEdx();
 254         hep->Fill(mom,dedx,1.);
 255
 256         Int_t pdg=(Int_t)ref->GetLength();  //this is particle's PDG !
 257
 258         if (TMath::Abs(pdg)==211) //pions
 259            if (mom>0.4 && mom<0.5) he->Fill(dedx,1.);
 260
 261      }
 262
 263      cout<<"\nList of Not found tracks :\n";
 264      for (k=0; k<nnotf; k++){
 265        cout<<notf[k]<<"\t";
 266        if ((k%9)==8) cout<<"\n";
 267      }
 268      cout<<"\n\nList of fake  tracks :\n";
 269      for (k=0; k<nfake; k++){
 270        cout<<fake[k]<<"\t";
 271        if ((k%9)==8) cout<<"\n";
 272      }
 273      cout<<"\n\nList of multiple found tracks :\n";
 274      for (k=0; k<nmult; k++) {
 275          cout<<"id.   "<<mult[k]
 276              <<"     found - "<<numb[k]<<"times\n";
 277      }
 278      cout<<endl;
 279
 280      cout<<"Number of found tracks : "<<nentr<<endl;
 281      cout<<"Number of \"good\" tracks : "<<ngood<<endl;
 282
 283      refs->Clear();
 284      delete event;
 285    } //***** End of the loop over events
 286
 287    ef->Close();
 288
 289    delete itsTree;
 290    refFile->Close();
 291
 292    Stat_t ng=hgood->GetEntries(), nf=hfound->GetEntries();
 293    if (ng!=0) cout<<"\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
 294    cout<<"Total number selected of \"good\" tracks ="<<allselected<<endl<<endl;
 295    cout<<"Total number of found tracks ="<<allfound<<endl;
 296    cout<<"Total number of \"good\" tracks ="<<allgood<<endl;
 297    cout<<endl;
 298
 299    gStyle->SetOptStat(111110);
 300    gStyle->SetOptFit(1);
 301
 302    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
 303
 304    Int_t minc=33;
 305
 306    TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
 307    p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
 308    hp->SetFillColor(4);  hp->SetXTitle("(mrad)");
 309    if (hp->GetEntries()<minc) hp->Draw(); else hp->Fit("gaus"); c1->cd();
 310
 311    TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
 312    p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
 313    hl->SetXTitle("(mrad)");
 314    if (hl->GetEntries()<minc) hl->Draw(); else hl->Fit("gaus"); c1->cd();
 315
 316    TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
 317    p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
 318    hpt->SetXTitle("(%)");
 319    if (hpt->GetEntries()<minc) hpt->Draw(); else hpt->Fit("gaus"); c1->cd();
 320
 321    TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
 322    p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
 323    hmpt->SetXTitle("(micron)");
 324    if (hmpt->GetEntries()<minc) hmpt->Draw(); else hmpt->Fit("gaus");
 325    hz->Draw("same"); c1->cd();
 326
 327
 328    TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
 329    p5->SetFillColor(41); p5->SetFrameFillColor(10);
 330    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 331    hg->Divide(hfound,hgood,1,1.,"b");
 332    hf->Divide(hfake,hgood,1,1.,"b");
 333    hg->SetMaximum(1.4);
 334    hg->SetYTitle("Tracking efficiency");
 335    hg->SetXTitle("Pt (GeV/c)");
 336    hg->Draw();
 337
 338    TLine *line1 = new TLine(0.2,1.0,6.1,1.0); line1->SetLineStyle(4);
 339    line1->Draw("same");
 340    TLine *line2 = new TLine(0.2,0.9,6.1,0.9); line2->SetLineStyle(4);
 341    line2->Draw("same");
 342
 343    hf->SetFillColor(1);
 344    hf->SetFillStyle(3013);
 345    hf->SetLineColor(2);
 346    hf->SetLineWidth(2);
 347    hf->Draw("histsame");
 348    TText *text = new TText(0.461176,0.248448,"Fake tracks");
 349    text->SetTextSize(0.05);
 350    text->Draw();
 351    text = new TText(0.453919,1.11408,"Good tracks");
 352    text->SetTextSize(0.05);
 353    text->Draw();
 354
 355    TCanvas *c2=new TCanvas("c2","",320,32,530,590);
 356    TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw();
 357    p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
 358    he->SetFillColor(2); he->SetFillStyle(3005);
 359    he->SetXTitle("Arbitrary Units");
 360    if (he->GetEntries()<minc) he->Draw(); else he->Fit("gaus"); c2->cd();
 361
 362    TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw();
 363    p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
 364    hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
 365    hep->Draw(); c1->cd();
 366
 367    TFile fc("AliITSComparisonV2.root","RECREATE");
 368    c1->Write();
 369    c2->Write();
 370    fc.Close();
 371
 372    gBenchmark->Stop("AliITSComparisonV2");
 373    gBenchmark->Show("AliITSComparisonV2");
 374
 375    return 0;
 376 }
 377
 378
 379
 380 Int_t GoodTracksITS(const Char_t *dir) {
 381    if (gAlice) {
 382        delete gAlice->GetRunLoader();
 383        delete gAlice;//if everything was OK here it is already NULL
 384        gAlice = 0x0;
 385    }
 386
 387    Char_t fname[100];
 388    sprintf(fname,"%s/galice.root",dir);
 389
 390    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
 391    if (!rl) {
 392       ::Error("GoodTracksITS","Can't start session !");
 393       return 1;
 394    }
 395
 396    rl->LoadgAlice();
 397    rl->LoadHeader();
 398    rl->LoadKinematics();
 399
 400    AliITS *ITS=(AliITS*)rl->GetAliRun()->GetDetector("ITS");
 401    if (!ITS) {
 402       ::Error("GoodTracksITS","Can't get the ITS !");
 403       delete rl;
 404       return 2;
 405    }
 406    AliITSgeom *geom=ITS->GetITSgeom();
 407    if (!geom) {
 408       ::Error("GoodTracksITS","Can't get the ITS geometry !");
 409       delete rl;
 410       return 3;
 411    }
 412
 413    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 414    if (itsl == 0x0) {
 415        ::Error("GoodTracksITS","Can not find the ITSLoader");
 416        delete rl;
 417        return 4;
 418    }
 419    itsl->LoadRecPoints();
 420
 421
 422    Int_t nev=rl->GetNumberOfEvents();
 423    ::Info("GoodTracksITS","Number of events : %d\n",nev);
 424
 425    sprintf(fname,"%s/GoodTracksTPC.root",dir);
 426    TFile *tpcFile=TFile::Open(fname);
 427    if ((!tpcFile)||(!tpcFile->IsOpen())) {
 428        ::Error("GoodTracksITS","Can't open the GoodTracksTPC.root !");
 429        delete rl;
 430        return 5;
 431    }
 432    TClonesArray dum("AliTrackReference",1000), *tpcRefs=&dum;
 433    TTree *tpcTree=(TTree*)tpcFile->Get("tpcTree");
 434    if (!tpcTree) {
 435        ::Error("GoodTracksITS","Can't get the TPC reference tree !");
 436        delete rl;
 437        return 6;
 438    }
 439    TBranch *tpcBranch=tpcTree->GetBranch("TPC");
 440    if (!tpcBranch) {
 441       ::Error("GoodTracksITS","Can't get the TPC reference branch !");
 442       delete rl;
 443       return 7;
 444    }
 445    tpcBranch->SetAddress(&tpcRefs);
 446
 447    sprintf(fname,"%s/GoodTracksITS.root",dir);
 448    TFile *itsFile=TFile::Open(fname,"recreate");
 449    TClonesArray dummy("AliTrackReference",1000), *itsRefs=&dummy;
 450    TTree itsTree("itsTree","Tree with info about the reconstructable ITS tracks");
 451    itsTree.Branch("ITS",&itsRefs);
 452
 453    //********  Loop over generated events
 454    for (Int_t e=0; e<nev; e++) {
 455      Int_t k;
 456
 457      rl->GetEvent(e);  itsFile->cd();
 458
 459      Int_t np = rl->GetHeader()->GetNtrack();
 460      cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
 461
 462      //******** Fill the "good" masks
 463      Int_t *good=new Int_t[np]; for (k=0; k<np; k++) good[k]=0;
 464
 465      TTree *cTree=itsl->TreeR();
 466      if (!cTree) {
 467         ::Error("GoodTracksITS","Can't get the cluster tree !");
 468         delete rl;
 469         return 8;
 470      }
 471      TBranch *branch=cTree->GetBranch("Clusters");
 472      if (!branch) {
 473         ::Error("GoodTracksITS","Can't get the clusters branch !");
 474         delete rl;
 475         return 9;
 476      }
 477      TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
 478      branch->SetAddress(&clusters);
 479
 480      Int_t entr=(Int_t)cTree->GetEntries();
 481      for (k=0; k<entr; k++) {
 482          cTree->GetEvent(k);
 483          Int_t ncl=clusters->GetEntriesFast(); if (ncl==0) continue;
 484          Int_t lay,lad,det;  geom->GetModuleId(k,lay,lad,det);
 485          if (lay<1 || lay>6) {
 486             ::Error("GoodTracksITS","Wrong layer !");
 487             delete rl;
 488             return 10;
 489          }
 490          while (ncl--) {
 491             AliITSclusterV2 *pnt=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
 492             Int_t l0=pnt->GetLabel(0);
 493                if (l0>=np) {
 494 //               cerr<<"Wrong label: "<<l0<<endl;
 495                  continue;
 496                }
 497             Int_t l1=pnt->GetLabel(1);
 498                if (l1>=np) {
 499 //               cerr<<"Wrong label: "<<l1<<endl;
 500                  continue;
 501                }
 502             Int_t l2=pnt->GetLabel(2);
 503                if (l2>=np) {
 504 //               cerr<<"Wrong label: "<<l2<<endl;
 505                  continue;
 506                }
 507             Int_t mask=1<<(lay-1);
 508             if (l0>=0) good[l0]|=mask;
 509             if (l1>=0) good[l1]|=mask;
 510             if (l2>=0) good[l2]|=mask;
 511          }
 512          clusters->Clear();
 513      }
 514
 515
 516      //****** select tracks which are "good" enough
 517      AliStack* stack = rl->Stack();
 518
 519      tpcTree->GetEvent(e);
 520      Int_t nk=tpcRefs->GetEntriesFast();
 521      Int_t nt=0;
 522      for (k=0; k<nk; k++) {
 523         AliTrackReference *tpcRef=(AliTrackReference *)tpcRefs->UncheckedAt(k);
 524         Int_t lab=tpcRef->Label();
 525         if (good[lab] != 0x3F) continue;
 526         TParticle *p = (TParticle*)stack->Particle(lab);
 527         if (p == 0x0) {
 528            cerr<<"Can not get particle "<<lab<<endl;
 529            continue;
 530         }
 531
 532         AliTrackReference *ref=new((*itsRefs)[nt]) AliTrackReference(*tpcRef);
 533         ref->SetMomentum(p->Px(),p->Py(),p->Pz());
 534         ref->SetPosition(p->Vx(),p->Vy(),p->Vz());
 535         nt++;
 536      }
 537      tpcRefs->Clear();
 538
 539      itsTree.Fill();
 540      itsRefs->Clear();
 541
 542      delete[] good;
 543
 544    } //*** end of the loop over generated events
 545
 546    itsTree.Write();
 547    itsFile->Close();
 548
 549    delete tpcTree;
 550    tpcFile->Close();
 551
 552    delete rl;
 553    return 0;
 554 }
 555
 556