ITS/AliV0Comparison.C

   1 /****************************************************************************
   2  *           Very important, delicate and rather obscure macro.             *
   3  *                                                                          *
   4  *                  Creates list of "findable" V0s,                         *
   5  *             calculates efficiency, resolutions etc.                      *
   6  *                                                                          *
   7  *   Origin: I.Belikov, IReS, Strasbourg, 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 "AliRunLoader.h"
  17   #include "AliITSLoader.h"
  18
  19   #include "TH1.h"
  20   #include "TFile.h"
  21   #include "TTree.h"
  22   #include "TObjArray.h"
  23   #include "TStyle.h"
  24   #include "TCanvas.h"
  25   #include "TLine.h"
  26   #include "TText.h"
  27   #include "TParticle.h"
  28   #include "TStopwatch.h"
  29
  30   #include "AliRun.h"
  31   #include "AliPDG.h"
  32   #include "AliV0vertex.h"
  33 #endif
  34
  35 struct GoodVertex {
  36   Int_t nlab,plab;
  37   Int_t code;
  38   Float_t px,py,pz;
  39   Float_t x,y,z;
  40 };
  41 Int_t good_vertices(GoodVertex *gt, Int_t max);
  42
  43 extern AliRun *gAlice;
  44
  45 Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
  46    cerr<<"Doing comparison...\n";
  47
  48    TStopwatch timer;
  49
  50    const Double_t V0window=0.05;
  51    Double_t ptncut=0.13, ptpcut=0.13, kinecut=0.03;
  52    Double_t V0mass=0.497672, V0width=0.020;
  53    switch (code) {
  54    case kK0Short:
  55         break;
  56    case kLambda0:
  57         V0mass=1.115683; V0width=0.015; ptpcut=0.50; kinecut=0.002;
  58         break;
  59    case kLambda0Bar:
  60         V0mass=1.115683; V0width=0.015; ptncut=0.50; kinecut=0.002;
  61         break;
  62    default: cerr<<"Invalid PDG code !\n"; return 1;
  63    }
  64
  65    TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution
  66    hp->SetXTitle("(mrad)"); hp->SetFillColor(2);
  67    TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30);
  68    hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013);
  69    TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.);
  70    hpt->SetXTitle("(%)"); hpt->SetFillColor(2);
  71
  72    TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res.
  73    hx->SetXTitle("(mm)"); hx->SetFillColor(6);
  74    TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.);   //y res
  75    hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013);
  76    TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.);   //z res.
  77    hz->SetXTitle("(mm)"); hz->SetFillColor(6);
  78
  79
  80    Double_t pmin=0.2, pmax=4.2; Int_t nchan=20;
  81    TH1F *hgood=new TH1F("hgood","Good Vertices",nchan,pmin,pmax);
  82    TH1F *hfound=new TH1F("hfound","Found Vertices",nchan,pmin,pmax);
  83    TH1F *hfake=new TH1F("hfake","Fake Vertices",nchan,pmin,pmax);
  84    TH1F *hg=new TH1F("hg","Efficiency for Good Vertices",nchan,pmin,pmax);
  85    hg->SetLineColor(4); hg->SetLineWidth(2);
  86    TH1F *hf=new TH1F("hf","Probability of Fake Vertices",nchan,pmin,pmax);
  87    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
  88
  89    Double_t mmin=V0mass-V0window, mmax=V0mass+V0window;
  90    TH1F *v0s =new TH1F("v0s","V0s Effective Mass",40, mmin, mmax);
  91    v0s->SetXTitle("(GeV)");
  92    v0s->SetLineColor(4); v0s->SetLineWidth(4);
  93    TH1F *v0sf =new TH1F("v0sf","Fake V0s Effective Mass",40, mmin, mmax);
  94    v0sf->SetXTitle("(GeV)"); v0sf->SetFillColor(6);
  95
  96
  97    if (gAlice) {
  98       delete gAlice->GetRunLoader();
  99       delete gAlice;
 100       gAlice=0;
 101    }
 102    AliRunLoader *rl = AliRunLoader::Open("galice.root");
 103    if (!rl) {
 104        cerr<<"AliV0Comparison.C :Can't start sesion !\n";
 105        return 1;
 106    }
 107    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 108    if (itsl == 0x0) {
 109        cerr<<"AliV0Comparison.C : Can not find the ITSLoader\n";
 110        delete rl;
 111        return 2;
 112    }
 113
 114    /*** Load reconstructed vertices ***/
 115    TObjArray varray(1000);
 116    itsl->LoadV0s();
 117    TTree *vTree=itsl->TreeV0();
 118    TBranch *branch=vTree->GetBranch("vertices");
 119    Int_t nentr=(Int_t)vTree->GetEntries();
 120    for (Int_t i=0; i<nentr; i++) {
 121        AliV0vertex *iovertex=new AliV0vertex; branch->SetAddress(&iovertex);
 122        vTree->GetEvent(i);
 123        varray.AddLast(iovertex);
 124    }
 125
 126    /*** Check if the file with the "good" vertices exists ***/
 127    GoodVertex gv[1000];
 128    Int_t ngood=0;
 129    ifstream in("good_vertices");
 130    if (in) {
 131       cerr<<"Reading good vertices...\n";
 132       while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
 133                  gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
 134                  gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
 135          ngood++;
 136          cerr<<ngood<<'\r';
 137          if (ngood==1000) {
 138             cerr<<"Too many good vertices !\n";
 139             break;
 140          }
 141       }
 142       if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
 143    } else {
 144      /*** generate a file with the "good" vertices ***/
 145       cerr<<"Marking good vertices (this will take a while)...\n";
 146       ngood=good_vertices(gv,1000);
 147       ofstream out("good_vertices");
 148       if (out) {
 149          for (Int_t ngd=0; ngd<ngood; ngd++)
 150             out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
 151                  gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
 152                  gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
 153       } else cerr<<"Can not open file (good_vertices) !\n";
 154       out.close();
 155    }
 156
 157
 158    Double_t pxg=0.,pyg=0.,ptg=0.;
 159    Int_t nlab=-1, plab=-1;
 160    Int_t i;
 161    for (i=0; i<nentr; i++) {
 162        AliV0vertex *vertex=(AliV0vertex*)varray.UncheckedAt(i);
 163        nlab=TMath::Abs(vertex->GetNindex());
 164        plab=TMath::Abs(vertex->GetPindex());
 165
 166        /** Kinematical cuts **/
 167        Double_t pxn,pyn,pzn; vertex->GetNPxPyPz(pxn,pyn,pzn);
 168        Double_t ptn=TMath::Sqrt(pxn*pxn + pyn*pyn);
 169        if (ptn < ptncut) continue;
 170        Double_t pxp,pyp,pzp; vertex->GetPPxPyPz(pxp,pyp,pzp);
 171        Double_t ptp=TMath::Sqrt(pxp*pxp + pyp*pyp);
 172        if (ptp < ptpcut) continue;
 173        Double_t kine=vertex->ChangeMassHypothesis(code);
 174        //if (TMath::Abs(kine)>kinecut) continue;
 175
 176        Double_t mass=vertex->GetEffMass();
 177        v0s->Fill(mass);
 178        v0sf->Fill(mass);
 179
 180        Int_t j;
 181        for (j=0; j<ngood; j++) {
 182           if (gv[j].code != vertex->GetPdgCode()) continue;
 183           if (gv[j].nlab == nlab)
 184           if (gv[j].plab == plab) break;
 185        }
 186
 187        if (TMath::Abs(mass-V0mass)>V0width) continue;
 188
 189        Double_t px,py,pz; vertex->GetPxPyPz(px,py,pz);
 190        Double_t pt=TMath::Sqrt(px*px+py*py);
 191
 192        if (j==ngood) {
 193           hfake->Fill(pt);
 194           cerr<<"Fake vertex: ("<<nlab<<","<<plab<<")\n";
 195           continue;
 196        }
 197        v0sf->Fill(mass,-1);
 198
 199        pxg=gv[j].px; pyg=gv[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
 200        Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px);
 201        Double_t lamg=TMath::ATan2(gv[j].pz,ptg), lam=TMath::ATan2(pz,pt);
 202        hp->Fill((phi - phig)*1000.);
 203        hl->Fill((lam - lamg)*1000.);
 204        hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.);
 205
 206        Double_t x,y,z; vertex->GetXYZ(x,y,z);
 207        hx->Fill((x-gv[j].x)*10);
 208        hy->Fill((y-gv[j].y)*10);
 209        hz->Fill((z-gv[j].z)*10);
 210
 211        hfound->Fill(ptg);
 212        gv[j].nlab=-1;
 213
 214    }
 215    for (i=0; i<ngood; i++) {
 216       if (gv[i].code != code) continue;
 217       pxg=gv[i].px; pyg=gv[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
 218       hgood->Fill(ptg);
 219       nlab=gv[i].nlab; plab=gv[i].plab;
 220       if (nlab < 0) continue;
 221       cerr<<"Vertex ("<<nlab<<','<<plab<<") has not been found !\n";
 222    }
 223
 224    varray.Delete();
 225
 226    Stat_t ng=hgood->GetEntries();
 227    Stat_t nf=hfound->GetEntries();
 228
 229    cerr<<"Number of found vertices: "<<nentr<<" ("<<nf<<" in the peak)\n";
 230    cerr<<"Number of good vertices: "<<ng<<endl;
 231
 232    if (ng!=0)
 233       cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
 234
 235    gStyle->SetOptStat(111110);
 236    gStyle->SetOptFit(1);
 237
 238    TCanvas *c1=new TCanvas("c1","",0,0,580,610);
 239    c1->Divide(2,2);
 240
 241    c1->cd(1);
 242    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 243    //hp->Fit("gaus");
 244    hp->Draw();
 245    hl->Draw("same"); c1->cd();
 246
 247    c1->cd(2);
 248    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 249    //hpt->Fit("gaus"); c1->cd();
 250    hpt->Draw(); c1->cd();
 251
 252    c1->cd(3);
 253    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 254    //hx->Fit("gaus");
 255    hx->Draw();
 256    hy->Draw("same"); c1->cd();
 257
 258    c1->cd(4);
 259    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 260    //hz->Fit("gaus");
 261    hz->Draw();
 262
 263
 264    TCanvas *c2=new TCanvas("c2","",600,0,580,610);
 265    c2->Divide(1,2);
 266
 267    c2->cd(1);
 268    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 269    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 270    hg->Divide(hfound,hgood,1,1.,"b");
 271    hf->Divide(hfake,hgood,1,1.,"b");
 272    hg->SetMaximum(1.4);
 273    hg->SetYTitle("V0 reconstruction efficiency");
 274    hg->SetXTitle("Pt (GeV/c)");
 275    hg->Draw();
 276
 277    TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
 278    line1->Draw("same");
 279    TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
 280    line2->Draw("same");
 281
 282    hf->SetFillColor(1);
 283    hf->SetFillStyle(3013);
 284    hf->SetLineColor(2);
 285    hf->SetLineWidth(2);
 286    hf->Draw("histsame");
 287    TText *text = new TText(0.461176,0.248448,"Fake vertices");
 288    text->SetTextSize(0.05);
 289    text->Draw();
 290    text = new TText(0.453919,1.11408,"Good vertices");
 291    text->SetTextSize(0.05);
 292    text->Draw();
 293
 294
 295    c2->cd(2);
 296    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
 297    //v0s->Fit("gaus","","",V0mass-V0width,V0mass+V0width);
 298    v0s->Draw();
 299    v0sf->Draw("same");
 300    Double_t max=v0s->GetMaximum();
 301    TLine *line3 = new TLine(V0mass-V0width,0.,V0mass-V0width,max);
 302    line3->Draw("same");
 303    TLine *line4 = new TLine(V0mass+V0width,0.,V0mass+V0width,max);
 304    line4->Draw("same");
 305
 306    timer.Stop(); timer.Print();
 307
 308    delete rl;
 309
 310    return 0;
 311 }
 312
 313 Int_t good_vertices(GoodVertex *gv, Int_t max) {
 314    Int_t nv=0;
 315    /*** Get information about the cuts ***/
 316    Double_t r2min=0.9*0.9;
 317    Double_t r2max=2.9*2.9;
 318
 319    /*** Get labels of the "good" tracks ***/
 320    Double_t dd; Int_t id, label[15000], ngt=0;
 321    ifstream in("good_tracks_its");
 322    if (!in) {
 323      cerr<<"Can't open the file good_tracks_its \n";
 324      return nv;
 325    }
 326    while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) {
 327      ngt++;
 328      if (ngt>=15000) {
 329        cerr<<"Too many good ITS tracks !\n";
 330        return nv;
 331      }
 332    }
 333    if (!in.eof()) {
 334       cerr<<"Read error (good_tracks_its) !\n";
 335       return nv;
 336    }
 337
 338    /*** Get an access to the kinematics ***/
 339    AliRunLoader *rl =
 340         AliRunLoader::GetRunLoader(AliConfig::fgkDefaultEventFolderName);
 341    if (rl == 0x0) {
 342      ::Fatal("AliV0Comparison.C::good_vertices","Can not find Run Loader !");
 343    }
 344
 345    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
 346    if (itsl == 0x0) {
 347        cerr<<"AliITSComparisonV2.C : Can not find TPCLoader\n";
 348        delete rl;
 349        return 1;
 350    }
 351    rl->LoadgAlice();
 352    rl->LoadHeader();
 353    rl->LoadKinematics();
 354    Int_t np = rl->GetHeader()->GetNtrack();
 355
 356    while (np--) {
 357       cerr<<np<<'\r';
 358       TParticle *p0=gAlice->GetMCApp()->Particle(np);
 359
 360       /*** only these V0s are "good" ***/
 361       Int_t code=p0->GetPdgCode();
 362       if (code!=kK0Short)
 363       if (code!=kLambda0)
 364       if (code!=kLambda0Bar) continue;
 365
 366       /*** daughter tracks should be "good" ***/
 367       Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
 368       if (nlab==plab) continue;
 369       if (nlab<0) continue;
 370       if (plab<0) continue;
 371       Int_t i;
 372       for (i=0; i<ngt; i++) if (label[i]==nlab) break;
 373       if (i==ngt) continue;
 374       for (i=0; i<ngt; i++) if (label[i]==plab) break;
 375       if (i==ngt) continue;
 376
 377       /*** fiducial volume ***/
 378       TParticle *p=gAlice->GetMCApp()->Particle(nlab);
 379       Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y;
 380       if (r2<r2min) continue;
 381       if (r2>r2max) continue;
 382
 383       if (gAlice->GetMCApp()->Particle(plab)->GetPDG()->Charge() < 0.) {
 384          i=plab; plab=nlab; nlab=i;
 385       }
 386
 387       gv[nv].code=code;
 388       gv[nv].plab=plab; gv[nv].nlab=nlab;
 389       gv[nv].px=p0->Px(); gv[nv].py=p0->Py(); gv[nv].pz=p0->Pz();
 390       gv[nv].x=x; gv[nv].y=y; gv[nv].z=z;
 391       nv++;
 392    }
 393
 394    delete rl;
 395
 396    return nv;
 397 }