TPC/AliTPCTestTracking.C

   1 void AliTPCTestTracking() {
   2    const char *pname="Param1";
   3    const char *tname="TreeD0_Param1";
   4
   5 // Dynamically link some shared libs
   6    if (gClassTable->GetID("AliRun") < 0) {
   7       gROOT->LoadMacro("loadlibs.C");
   8       loadlibs();
   9    } else {
  10       delete gAlice;
  11       gAlice=0;
  12    }
  13
  14 // Connect the Root Galice file containing Geometry, Kine and Hits
  15    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  16    if (!file) file = new TFile("galice.root");
  17
  18 // Get AliRun object from file or create it if not on file
  19    if (!gAlice) {
  20       gAlice = (AliRun*)file->Get("gAlice");
  21       if (gAlice) printf("AliRun object found on file\n");
  22       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  23    }
  24
  25    gAlice->GetEvent(0);
  26
  27    TClonesArray *particles=gAlice->Particles();
  28    int np=particles->GetEntriesFast();
  29
  30    AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
  31    int ver=TPC->IsVersion();
  32    cerr<<"TPC version "<<ver<<" has been found !\n";
  33
  34    AliTPCD *digp= (AliTPCD*)file->Get(pname);
  35    if (digp!=0) TPC->SetDigParam(digp);
  36    else cerr<<"Warning: default TPC parameters will be used !\n";
  37
  38    int nrows=TPC->GetDigParam()->GetParam().GetNRowUp()-1;
  39    switch (ver) {
  40    case 1:
  41       cerr<<"Making clusters...\n";
  42       TPC->Hits2Clusters();
  43       TClonesArray *clusters=TPC->Clusters();
  44       if (!clusters) {cerr<<"No clusters found !\n"; return;}
  45       int n=clusters->GetEntriesFast();
  46       cerr<<"Number of clusters "<<n<<"                                  \n";
  47
  48       cerr<<"Marking \"good\" tracks...                                  \n";
  49       for (int i=0; i<n; i++) {
  50           AliTPCcluster *c=(AliTPCcluster*)clusters->UncheckedAt(i);
  51           int lab=c->fTracks[0];
  52           if (lab<0) continue; //noise cluster
  53           lab=TMath::Abs(lab);
  54           int sector=c->fSector-1, row=c->fPadRow-1;
  55           GParticle *p=(GParticle*)particles->UncheckedAt(lab);
  56           int ks;
  57           if (row==nrows)   {ks=p->GetKS()|0x1000; p->SetKS(ks);}
  58           if (row==nrows-8) {ks=p->GetKS()|0x800; p->SetKS(ks);}
  59           ks=p->GetKS()+1; p->SetKS(ks);
  60       }
  61
  62       break;
  63    case 2:
  64       cerr<<"Looking for clusters...\n";
  65       TPC->Digits2Clusters();
  66       TClonesArray *clusters=TPC->Clusters();
  67       if (!clusters) {cerr<<"No clusters found !\n"; return;}
  68       int n=clusters->GetEntriesFast();
  69       cerr<<"Number of clusters "<<n<<"                                  \n";
  70
  71       cerr<<"Marking \"good\" tracks...                                  \n";
  72       TTree *TD=gDirectory->Get(tname);
  73       TClonesArray *digits=TPC->Digits();
  74       TD->GetBranch("Digits")->SetAddress(&digits);
  75
  76       int *count = new int[np];
  77       int i;
  78       for (i=0; i<np; i++) count[i]=0;
  79       int sectors_by_rows=(int)TD->GetEntries();
  80       for (i=0; i<sectors_by_rows; i++) {
  81           if (!TD->GetEvent(i)) continue;
  82           int sec, row;
  83           int ndigits=digits->GetEntriesFast();
  84           int j;
  85           for (j=0; j<ndigits; j++) {
  86               AliTPCdigit *dig = (AliTPCdigit*)digits->UncheckedAt(j);
  87               int idx0=dig->fTracks[0];
  88               int idx1=dig->fTracks[1];
  89               int idx2=dig->fTracks[2];
  90               sec=dig->fSector-1; row=dig->fPadRow-1;
  91               if (idx0>=0 && dig->fSignal>0) count[idx0]+=1;
  92               if (idx1>=0 && dig->fSignal>0) count[idx1]+=1;
  93               if (idx2>=0 && dig->fSignal>0) count[idx2]+=1;
  94           }
  95           for (j=0; j<np; j++) {
  96               GParticle *p=(GParticle*)particles->UncheckedAt(j);
  97               if (count[j]>1) {
  98                  int ks;
  99                  if (row==nrows   ) {ks=p->GetKS()|0x1000; p->SetKS(ks);}
 100                  if (row==nrows-8 ) {ks=p->GetKS()|0x800;  p->SetKS(ks);}
 101                  ks=p->GetKS()+1; p->SetKS(ks);
 102               }
 103               count[j]=0;
 104           }
 105       }
 106       delete[] count;
 107
 108       break;
 109    default:
 110       cerr<<"Invalid TPC version !\n";
 111       return;
 112    }
 113
 114    cerr<<"Looking for tracks...\n";
 115    TPC->Clusters2Tracks();
 116    int nt=0;
 117    TClonesArray *tracks=TPC->Tracks();
 118    if (tracks) nt=tracks->GetEntriesFast();
 119    cerr<<"Number of found tracks "<<nt<<endl;
 120
 121 /////////////////////////////////////////////////////////////////////////
 122    cerr<<"Doing comparison...\n";
 123    TH1F *hp=new TH1F("hp","PHI resolution",50,-100.,100.); hp->SetFillColor(4);
 124    TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-100,100); hl->SetFillColor(4);
 125    TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
 126    hpt->SetFillColor(2);
 127    TH1F *hd=new TH1F("hd","Impact parameter distribution ",30,0,25);
 128    hd->SetFillColor(6);
 129
 130    TH1F *hgood=new TH1F("hgood","Good tracks",20,0,2);
 131    TH1F *hfound=new TH1F("hfound","Found tracks",20,0,2);
 132    TH1F *hfake=new TH1F("hfake","Fake tracks",20,0,2);
 133    TH1F *hg=new TH1F("hg","",20,0,2); //efficiency for good tracks
 134    hg->SetLineColor(4); hg->SetLineWidth(2);
 135    TH1F *hf=new TH1F("hf","Efficiency for fake tracks",20,0,2);
 136    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
 137
 138    for (int i=0; i<np; i++) {
 139       GParticle *p = (GParticle*)particles->UncheckedAt(i);
 140       if (p->GetParent()>=0) continue;  //secondary particle
 141       if (p->GetKS()<0x1000+0x800+2+30) continue;
 142       Double_t ptg=p->GetPT(),pxg=p->GetPx(),pyg=p->GetPy(),pzg=p->GetPz();
 143       if (ptg<0.100) continue;
 144       if (fabs(pzg/ptg)>0.999) continue;
 145
 146       //cout<<i<<endl;
 147
 148       hgood->Fill(ptg);
 149       int found=0;
 150       for (int j=0; j<nt; j++) {
 151           AliTPCtrack *track=(AliTPCtrack*)tracks->UncheckedAt(j);
 152           int lab=track->GetLab();
 153           if (fabs(lab)!=i) continue;
 154           //if (lab!=i) continue;
 155           found=1;
 156           Double_t xk=76.;
 157
 158           track->PropagateTo(xk);
 159           xk-=0.11;
 160           track->PropagateTo(xk,42.7,2.27); //C
 161           xk-=2.6;
 162           track->PropagateTo(xk,36.2,1.98e-3); //C02
 163           xk-=0.051;
 164           track->PropagateTo(xk,42.7,2.27); //C
 165
 166           xk-=0.4;
 167           track->PropagateTo(xk,21.82,2.33); //ITS+beam_pipe+etc (approximately)
 168
 169           track->PropagateToVertex(); //comparison should be done at the vertex
 170
 171           if (lab==i) hfound->Fill(ptg);
 172           else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
 173           Double_t px,py,pz,pt=fabs(track->GetPt());track->GetPxPyPz(px,py,pz);
 174           Double_t phig=TMath::ATan(pyg/pxg);
 175           Double_t phi =TMath::ATan(py /px );
 176           hp->Fill((phi - phig)*1000.);
 177           Double_t lamg=TMath::ATan(pzg/ptg);
 178           Double_t lam =TMath::ATan(pz /pt );
 179           hl->Fill((lam - lamg)*1000.);
 180           hpt->Fill((pt - ptg)/ptg*100.);
 181           Double_t x,y,z; track->GetXYZ(x,y,z);
 182           hd->Fill(sqrt(x*x + y*y + z*z)*10.);
 183           break;
 184       }
 185       if (!found) cerr<<"Track number "<<i<<" was not found !\n";
 186    }
 187    Stat_t ngood =hgood->GetEntries(); cerr<<"Good tracks "<<ngood<<endl;
 188    Stat_t nfound=hfound->GetEntries();
 189    if (ngood!=0)
 190       cerr<<"Integral efficiency is about "<<nfound/ngood*100.<<" %\n";
 191
 192    gStyle->SetOptStat(111110);
 193    gStyle->SetOptFit(1);
 194
 195    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
 196
 197    TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
 198    p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
 199    hp->SetFillColor(4);  hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd();
 200
 201    TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
 202    p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
 203    hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd();
 204
 205    TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
 206    p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
 207    hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd();
 208
 209    TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
 210    p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
 211    hd->SetXTitle("(mm)"); hd->Draw(); c1->cd();
 212
 213    TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
 214    /*p5->SetTopMargin(0.25);*/ p5->SetFillColor(41); p5->SetFrameFillColor(10);
 215    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 216    hg->Divide(hfound,hgood,1,1.,"b");
 217    hf->Divide(hfake,hgood,1,1.,"b");
 218    hg->SetMaximum(1.4);
 219    hg->SetYTitle("Tracking efficiency");
 220    hg->SetXTitle("Pt (GeV/c)");
 221    hg->Draw();
 222
 223    TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4);
 224    line1->Draw("same");
 225    TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4);
 226    line2->Draw("same");
 227
 228    hf->SetFillColor(1);
 229    hf->SetFillStyle(3013);
 230    hf->SetLineColor(2);
 231    hf->SetLineWidth(2);
 232    hf->Draw("histsame");
 233    TText *text = new TText(0.461176,0.248448,"Fake tracks");
 234    text->SetTextSize(0.05);
 235    text->Draw();
 236    text = new TText(0.453919,1.11408,"Good tracks");
 237    text->SetTextSize(0.05);
 238    text->Draw();
 239 }
 240