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    int *good=new int[np];
  30    for (int ii=0; ii<np; ii++) good[ii]=0;
  31
  32    AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
  33    int ver=TPC->IsVersion();
  34    cerr<<"TPC version "<<ver<<" has been found !\n";
  35
  36    AliTPCD *digp= (AliTPCD*)file->Get(pname);
  37    if (digp!=0) TPC->SetDigParam(digp);
  38    else cerr<<"Warning: default TPC parameters will be used !\n";
  39
  40    int nrow_up=TPC->GetDigParam()->GetParam().GetNRowUp();
  41    int nrows=TPC->GetDigParam()->GetParam().GetNRowLow()+nrow_up;
  42    int zero=TPC->GetDigParam()->GetParam().GetZeroSup();
  43    int gap=int(0.125*nrows);
  44    int good_number=int(0.4*nrows);
  45
  46    switch (ver) {
  47    case 1:
  48       cerr<<"Making clusters...\n";
  49       TPC->Hits2Clusters();
  50       TClonesArray *clusters=TPC->Clusters();
  51       if (!clusters) {cerr<<"No clusters found !\n"; return;}
  52       int n=clusters->GetEntriesFast();
  53       cerr<<"Number of clusters "<<n<<"                                  \n";
  54
  55       cerr<<"Marking \"good\" tracks...                                  \n";
  56       for (int i=0; i<n; i++) {
  57           AliTPCcluster *c=(AliTPCcluster*)clusters->UncheckedAt(i);
  58           int lab=c->fTracks[0];
  59           if (lab<0) continue; //noise cluster
  60           lab=TMath::Abs(lab);
  61           int row=c->fPadRow;
  62           if (row==nrow_up-1    ) good[lab]|=0x1000;
  63           if (row==nrow_up-1-gap) good[lab]|=0x800;
  64           good[lab]++;
  65       }
  66
  67       break;
  68    case 2:
  69       cerr<<"Looking for clusters...\n";
  70       TPC->Digits2Clusters();
  71       TClonesArray *clusters=TPC->Clusters();
  72       if (!clusters) {cerr<<"No clusters found !\n"; return;}
  73       int n=clusters->GetEntriesFast();
  74       cerr<<"Number of clusters "<<n<<"                                  \n";
  75
  76       cerr<<"Marking \"good\" tracks...                                  \n";
  77       TTree *TD=gDirectory->Get(tname);
  78       TClonesArray *digits=TPC->Digits();
  79       TD->GetBranch("Digits")->SetAddress(&digits);
  80
  81       int *count = new int[np];
  82       int i;
  83       for (i=0; i<np; i++) count[i]=0;
  84       int sectors_by_rows=(int)TD->GetEntries();
  85       for (i=0; i<sectors_by_rows; i++) {
  86           if (!TD->GetEvent(i)) continue;
  87           int row;
  88           int ndigits=digits->GetEntriesFast();
  89           int j;
  90           for (j=0; j<ndigits; j++) {
  91               AliTPCdigit *dig = (AliTPCdigit*)digits->UncheckedAt(j);
  92               int idx0=dig->fTracks[0];
  93               int idx1=dig->fTracks[1];
  94               int idx2=dig->fTracks[2];
  95               row=dig->fPadRow;
  96               if (idx0>=0 && dig->fSignal>=zero) count[idx0]+=1;
  97               if (idx1>=0 && dig->fSignal>=zero) count[idx1]+=1;
  98               if (idx2>=0 && dig->fSignal>=zero) count[idx2]+=1;
  99           }
 100           for (j=0; j<np; j++) {
 101               if (count[j]>1) {
 102                  int ks;
 103                  if (row==nrow_up-1    ) good[j]|=0x1000;
 104                  if (row==nrow_up-1-gap) good[j]|=0x800;
 105                  good[j]++;
 106               }
 107               count[j]=0;
 108           }
 109       }
 110       delete[] count;
 111
 112       break;
 113    default:
 114       cerr<<"Invalid TPC version !\n";
 115       return;
 116    }
 117
 118    cerr<<"Looking for tracks...\n";
 119    TPC->Clusters2Tracks();
 120    int nt=0;
 121    TClonesArray *tracks=TPC->Tracks();
 122    if (tracks) nt=tracks->GetEntriesFast();
 123    cerr<<"Number of found tracks "<<nt<<endl;
 124
 125 /////////////////////////////////////////////////////////////////////////
 126    cerr<<"Doing comparison...\n";
 127    TH1F *hp=new TH1F("hp","PHI resolution",50,-100.,100.); hp->SetFillColor(4);
 128    TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-100,100); hl->SetFillColor(4);
 129    TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
 130    hpt->SetFillColor(2);
 131    TH1F *hd=new TH1F("hd","Impact parameter distribution ",30,0,25);
 132    hd->SetFillColor(6);
 133
 134    TH1F *hgood=new TH1F("hgood","Good tracks",20,0,10);
 135    TH1F *hfound=new TH1F("hfound","Found tracks",20,0,10);
 136    TH1F *hfake=new TH1F("hfake","Fake tracks",20,0,10);
 137    TH1F *hg=new TH1F("hg","",20,0,10); //efficiency for good tracks
 138    hg->SetLineColor(4); hg->SetLineWidth(2);
 139    TH1F *hf=new TH1F("hf","Efficiency for fake tracks",20,0,10);
 140    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
 141
 142    TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,500.);
 143    TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,500.);
 144
 145    for (int i=0; i<np; i++) {
 146       TParticle *p = (TParticle*)particles->UncheckedAt(i);
 147       if (p->GetFirstMother()>=0) continue;  //secondary particle
 148       if (good[i] < 0x1000+0x800+2+good_number) continue;
 149       Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py(),pzg=p->Pz();
 150       if (ptg<0.100) continue;
 151       if (fabs(pzg/ptg)>0.999) continue;
 152
 153       //cout<<i<<endl;
 154
 155       hgood->Fill(ptg);
 156       int found=0;
 157       for (int j=0; j<nt; j++) {
 158           AliTPCtrack *track=(AliTPCtrack*)tracks->UncheckedAt(j);
 159           int lab=track->GetLabel(nrows);
 160           if (fabs(lab)!=i) continue;
 161
 162           found=1;
 163
 164           Double_t xk=76.;
 165           track->PropagateTo(xk);
 166           xk-=0.11;
 167           track->PropagateTo(xk,42.7,2.27); //C
 168           xk-=26.;
 169           track->PropagateTo(xk,36.2,1.98e-3); //C02
 170           xk-=0.051;
 171           track->PropagateTo(xk,42.7,2.27); //C
 172           xk-=25.;
 173           track->PropagateTo(xk,36.7,1.29e-3);//Air
 174           xk-=0.4;                           //  +
 175           track->PropagateTo(xk,21.82,2.33);//ITS+beam_pipe+etc (approximately)
 176
 177           track->PropagateToVertex(); //comparison should be done at the vertex
 178
 179           if (lab==i) hfound->Fill(ptg);
 180           else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
 181           Double_t px,py,pz,pt=fabs(track->GetPt());track->GetPxPyPz(px,py,pz);
 182           Double_t phig=TMath::ATan(pyg/pxg);
 183           Double_t phi =TMath::ATan(py /px );
 184           hp->Fill((phi - phig)*1000.);
 185           Double_t lamg=TMath::ATan(pzg/ptg);
 186           Double_t lam =TMath::ATan(pz /pt );
 187           hl->Fill((lam - lamg)*1000.);
 188           hpt->Fill((pt - ptg)/ptg*100.);
 189           Double_t x=track->GetX(), y=track->GetY(), z=track->GetZ();
 190           hd->Fill(sqrt(x*x + y*y + z*z)*10.);
 191
 192           Double_t mom=track->GetP();
 193           Double_t dedx=track->GetdEdX(0.05,0.80)/
 194                         digp->GetParam().GetPadPitchLength();
 195
 196           hep->Fill(mom,dedx,1.);
 197           if (p->GetPdgCode()==211 || p->GetPdgCode()==-211)
 198             if (mom>0.4 && mom<0.5)
 199                    he->Fill(dedx,1.);
 200
 201           break;
 202       }
 203       if (!found) cerr<<"Track number "<<i<<" was not found !\n";
 204    }
 205
 206    delete[] good;
 207
 208    Stat_t ngood =hgood->GetEntries(); cerr<<"Good tracks "<<ngood<<endl;
 209    Stat_t nfound=hfound->GetEntries();
 210    if (ngood!=0)
 211       cerr<<"Integral efficiency is about "<<nfound/ngood*100.<<" %\n";
 212
 213    gStyle->SetOptStat(111110);
 214    gStyle->SetOptFit(1);
 215
 216    TCanvas *c1=new TCanvas("c1","",0,0,700,850);
 217
 218    TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
 219    p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
 220    hp->SetFillColor(4);  hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd();
 221
 222    TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
 223    p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
 224    hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd();
 225
 226    TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
 227    p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
 228    hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd();
 229
 230    TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
 231    p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
 232    hd->SetXTitle("(mm)"); hd->Draw(); c1->cd();
 233
 234    TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
 235    p5->SetFillColor(41); p5->SetFrameFillColor(10);
 236    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
 237    hg->Divide(hfound,hgood,1,1.,"b");
 238    hf->Divide(hfake,hgood,1,1.,"b");
 239    hg->SetMaximum(1.4);
 240    hg->SetYTitle("Tracking efficiency");
 241    hg->SetXTitle("Pt (GeV/c)");
 242    hg->Draw();
 243
 244    TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4);
 245    line1->Draw("same");
 246    TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4);
 247    line2->Draw("same");
 248
 249    hf->SetFillColor(1);
 250    hf->SetFillStyle(3013);
 251    hf->SetLineColor(2);
 252    hf->SetLineWidth(2);
 253    hf->Draw("histsame");
 254    TText *text = new TText(0.461176,0.248448,"Fake tracks");
 255    text->SetTextSize(0.05);
 256    text->Draw();
 257    text = new TText(0.453919,1.11408,"Good tracks");
 258    text->SetTextSize(0.05);
 259    text->Draw();
 260
 261
 262
 263    TCanvas *c2=new TCanvas("c2","",320,32,530,590);
 264
 265    TPad *p6=new TPad("p6","",0.,0.,1.,.5); p6->Draw();
 266    p6->cd(); p6->SetFillColor(42); p6->SetFrameFillColor(10);
 267    he->SetFillColor(2); he->SetFillStyle(3005);
 268    he->SetXTitle("Arbitrary Units");
 269    he->Fit("gaus"); c2->cd();
 270
 271    TPad *p7=new TPad("p7","",0.,0.5,1.,1.); p7->Draw();
 272    p7->cd(); p7->SetFillColor(42); p7->SetFrameFillColor(10);
 273    hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
 274    hep->Draw(); c1->cd();
 275
 276
 277 }
 278