HLT/exa/evalrowhough.C

   1 // $Id$
   2
   3 struct GoodTrack
   4 {
   5   Int_t label;
   6   Double_t eta;
   7   Int_t code;
   8   Double_t px,py,pz;
   9   Double_t x,y,z;
  10   Int_t nhits;
  11   Int_t sector;
  12 };
  13 typedef struct GoodTrack GoodTrack;
  14
  15 //Histograms
  16 TNtuple *fNtuppel;
  17 TH1F *fPtRes;
  18 TH1F *fNGoodTracksPt;
  19 TH1F *fNFoundTracksPt;
  20 TH1F *fNFakeTracksPt;
  21 TH1F *fTrackEffPt;
  22 TH1F *fFakeTrackEffPt;
  23 TH1F *fNGoodTracksEta;
  24 TH1F *fNFoundTracksEta;
  25 TH1F *fNFakeTracksEta;
  26 TH1F *fTrackEffEta;
  27 TH1F *fFakeTrackEffEta;
  28 TH1F *fNGoodTracksPhi;
  29 TH1F *fNFoundTracksPhi;
  30 TH1F *fNFakeTracksPhi;
  31 TH1F *fTrackEffPhi;
  32 TH1F *fFakeTrackEffPhi;
  33 TProfile *fFakesPt;
  34 TProfile *fFakesPhi;
  35 TProfile *fSlicesPt;
  36 TProfile *fSlicesPhi;
  37 TProfile2D *fFakesPtvsPhi;
  38 TNtuple *fNtupleRes;
  39 TH1F *fPtRes2;
  40 TH1F *fEtaRes;
  41 TH1F *fPsiRes;
  42
  43 TH1F *fNGoodTracksRad;
  44 TH1F *fNFoundTracksRad;
  45 TH1F *fNGoodTracksZ;
  46 TH1F *fNFoundTracksZ;
  47 TH1F *fRadEff;
  48 TH1F *fZEff;
  49
  50 void evalrowhough(Char_t *path="./fitter",Char_t *offlinepath="./",Int_t ievent=0, Float_t ptmin=0.4)
  51 {
  52   CreateHistos(25,0.1,4.1);
  53
  54   Char_t fname[1024];
  55   sprintf(fname,"%s/tracks_%d.raw",path,ievent);
  56   AliL3FileHandler *tfile = new AliL3FileHandler();
  57   if(!tfile->SetBinaryInput(fname)){
  58     LOG(AliL3Log::kError,"AliL3Evaluation::Setup","File Open")
  59       <<"Inputfile "<<fname<<" does not exist"<<ENDLOG;
  60     return;
  61   }
  62
  63   AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
  64   tfile->Binary2TrackArray(fTracks);
  65   //fTracks->QSort();
  66   tfile->CloseBinaryInput();
  67   delete tfile;
  68
  69   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
  70     {
  71       AliL3HoughTrack *track = fTracks->GetCheckedTrack(i);
  72       if(!track) continue;
  73
  74       track->SetEta(fTracks->GetCheckedTrack(i)->GetPseudoRapidity());
  75       UInt_t *ids = track->GetHitNumbers();
  76       Int_t slice = (ids[0]>>25)&0x7f;
  77       track->SetSlice(slice);
  78       cout<<"Track with label "<<ids<<" "<<track->GetNHits()<<" "<<track->GetMCid()<<" "<<ids[0]<<" "<<track->GetWeight()<<" "<<track->GetSlice()<<endl;
  79     }
  80
  81   Char_t filename[1024];
  82   sprintf(filename,"%s/good_tracks_tpc_%d",offlinepath,ievent);
  83   ifstream in(filename);
  84   if(!in)
  85     {
  86       cerr<<"AliL3Evaluate::GetGoodParticles : Problems opening file :"<<filename<<endl;
  87       return;
  88     }
  89
  90   Int_t MaxTracks=20000;
  91   Int_t fGoodGen=0;
  92   fGoodTracks = new GoodTrack[MaxTracks];
  93
  94   while (in>>fGoodTracks[fGoodGen].label>>fGoodTracks[fGoodGen].code>>
  95          fGoodTracks[fGoodGen].px>>fGoodTracks[fGoodGen].py>>fGoodTracks[fGoodGen].pz>>
  96          fGoodTracks[fGoodGen].x>>fGoodTracks[fGoodGen].y>>fGoodTracks[fGoodGen].z)
  97     {
  98       fGoodTracks[fGoodGen].nhits=-1;
  99       fGoodTracks[fGoodGen].sector=-1;
 100       fGoodGen++;
 101       if (fGoodGen==MaxTracks)
 102         {
 103           cerr<<"AliL3Evaluate::GetGoodParticles : Too many good tracks !\n";
 104           break;
 105         }
 106     }
 107
 108   Int_t fGoodGen2 = 0;
 109   Int_t fGoodFound = 0;
 110   Int_t fGoodFound2 = 0;
 111   Float_t fMinGoodPt = ptmin;
 112   Float_t fMaxGoodPt = 99999.;
 113
 114   if(!fGoodTracks)
 115     {
 116       cerr<<"AliL3Evaluate::FillEffHistos : No good tracks"<<endl;
 117       return;
 118     }
 119   cout<<"Comparing "<<fGoodGen<<" good tracks ..."<<endl;
 120   for(Int_t i=0; i<fGoodGen; i++)
 121     {
 122       //cout<<"Checking particle "<<i<<endl;
 123       Float_t ptg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py);
 124       Float_t pg = TMath::Sqrt(fGoodTracks[i].px*fGoodTracks[i].px + fGoodTracks[i].py*fGoodTracks[i].py + fGoodTracks[i].pz*fGoodTracks[i].pz);
 125       if(ptg < fMinGoodPt || ptg > fMaxGoodPt) continue;
 126       Float_t pzg=fGoodTracks[i].pz;
 127       Float_t dipangle=TMath::ATan2(pzg,ptg)*180./TMath::Pi();
 128       Float_t etag=0.5 * TMath::Log((pg+pzg)/(pg-pzg));
 129       Float_t phig=TMath::ATan2(fGoodTracks[i].py,fGoodTracks[i].px);
 130       Float_t rg = TMath::Sqrt(fGoodTracks[i].x*fGoodTracks[i].x + fGoodTracks[i].y*fGoodTracks[i].y);
 131       Float_t zg = TMath::Abs(fGoodTracks[i].z);
 132       if (phig<0) phig+=2*TMath::Pi();
 133
 134       fNGoodTracksPt->Fill(ptg);
 135       fNGoodTracksEta->Fill(etag);
 136       fNGoodTracksPhi->Fill(phig);
 137       if((zg != 0) || (rg != 0)) {
 138         fNGoodTracksRad->Fill(rg);
 139         fNGoodTracksZ->Fill(zg);
 140       }
 141       fGoodGen2++;
 142       Int_t found = 0;
 143       Int_t nslices = 1;
 144       Int_t org_slice;
 145       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 146         {
 147           AliL3HoughTrack *track = (AliL3HoughTrack *)fTracks->GetCheckedTrack(k);
 148           if(!track) continue;
 149           //      Int_t nHits = track->GetNumberOfPoints();
 150           //      if(nHits < fMinPointsOnTrack) break;
 151           Int_t nHits = 1;
 152           Int_t tracklabel;
 153           tracklabel = track->GetMCid();
 154
 155           if(TMath::Abs(tracklabel) != fGoodTracks[i].label) continue;
 156           if(found == 0) {
 157             found=1;
 158             Float_t pt=track->GetPt();
 159             org_slice = track->GetSlice();
 160             Float_t eta = track->GetEta();
 161             Float_t phi = track->GetPsi();
 162             //cout<<eta<<" "<<phi<<" "<<etag<<" "<<phig<<" "<<track->GetEtaIndex()<<endl;
 163             if(tracklabel == fGoodTracks[i].label)
 164               {
 165                 fGoodFound++;
 166
 167                 fNFoundTracksPt->Fill(ptg);
 168                 fNFoundTracksEta->Fill(etag);
 169                 fNFoundTracksPhi->Fill(phig);
 170                 if((zg != 0) || (rg != 0)) {
 171                   fNFoundTracksRad->Fill(rg);
 172                   fNFoundTracksZ->Fill(zg);
 173                 }
 174                 fNtuppel->Fill(ptg,pt,nHits);
 175                 fPtRes->Fill((pt-ptg)/ptg*100.);
 176                 fEtaRes->Fill(eta-etag);
 177                 fPsiRes->Fill(phi-phig);
 178               }
 179             else
 180               {
 181                 fNFakeTracksPt->Fill(ptg);
 182                 fNFakeTracksEta->Fill(etag);
 183                 fNFakeTracksPhi->Fill(phig);
 184               }
 185             //fPtRes->Fill((pt-ptg)/ptg*100.);
 186             //fNtuppel->Fill(ptg,pt,nHits);
 187           }
 188           else {
 189             found++;
 190             if(track->GetSlice() != org_slice)
 191               nslices = 2;
 192           }
 193         }
 194       if(!found)
 195         //cout<<"Track "<<fGoodTracks[i].label<<" was not found"<<endl;
 196         continue;
 197       else {
 198         fFakesPt->Fill(ptg,(Double_t)found-1.0);
 199         fFakesPhi->Fill(phig,(Double_t)found-1.0);
 200         fSlicesPt->Fill(ptg,(Double_t)nslices);
 201         fSlicesPhi->Fill(phig,(Double_t)nslices);
 202         //Float_t local_phi = phi-(phi/)
 203         //fFakesPtvsPhi->Fill(1.0/ptg,
 204       }
 205       fGoodFound2 += found;
 206     }
 207
 208   cout<<fGoodFound<<"("<<fGoodFound2<<") tracks found out of "<<fGoodGen2<<" generated tracks"<<endl;
 209   cout<<" The overall efficiency is "<<(Float_t)fGoodFound/(Float_t)fGoodGen2<<endl;
 210
 211   CalcEffHistos();
 212   plotptres(ptmin);
 213   Write2File();
 214 }
 215
 216 void CreateHistos(Int_t nbin,Float_t xlow,Float_t xup)
 217 {
 218   //Create the histograms
 219
 220   fNtuppel = new TNtuple("fNtuppel","Pt resolution","pt_gen:pt_found:nHits");
 221   fNtuppel->SetDirectory(0);
 222   fPtRes = new TH1F("fPtRes","Relative Pt resolution",30,-10.,10.);
 223   fNGoodTracksPt = new TH1F("fNGoodTracksPt","Good tracks vs pt",nbin,xlow,xup);
 224   fNFoundTracksPt = new TH1F("fNFoundTracksPt","Found tracks vs pt",nbin,xlow,xup);
 225   fNFakeTracksPt = new TH1F("fNFakeTracksPt","Fake tracks vs pt",nbin,xlow,xup);
 226   fTrackEffPt = new TH1F("fTrackEffPt","Tracking efficiency vs pt",nbin,xlow,xup);
 227   fFakeTrackEffPt = new TH1F("fFakeTrackEffPt","Efficiency for fake tracks vs pt",nbin,xlow,xup);
 228
 229   fNGoodTracksEta = new TH1F("fNGoodTracksEta","Good tracks vs eta",nbin,-1,1);
 230   fNFoundTracksEta = new TH1F("fNFoundTracksEta","Found tracks vs eta",nbin,-1,1);
 231   fNFakeTracksEta = new TH1F("fNFakeTracksEta","Fake tracks vs eta",nbin,-1,1);
 232   fTrackEffEta = new TH1F("fTrackEffEta","Tracking efficienct vs eta",nbin,-1,1);
 233   fFakeTrackEffEta = new TH1F("fFakeTrackEffEta","Efficiency for fake tracks vs eta",nbin,-1,1);
 234
 235
 236   fNGoodTracksPhi = new TH1F("fNGoodTracksPhi","Good tracks vs phi",nbin,0,6.28);
 237   fNFoundTracksPhi = new TH1F("fNFoundTracksPhi","Found tracks vs phi",nbin,0,6.28);
 238   fNFakeTracksPhi = new TH1F("fNFakeTracksPhi","Fake tracks vs phi",nbin,0,6.28);
 239   fTrackEffPhi = new TH1F("fTrackEffPhi","Tracking efficienct vs phi",nbin,0,6.28);
 240   fFakeTrackEffPhi = new TH1F("fFakeTrackEffPhi","Efficiency for fake tracks vs phi",nbin,0,6.28);
 241   fFakesPt = new TProfile("fFakesPt","Number of ghosts vs pt",nbin,xlow,xup,0,100);
 242   fFakesPhi = new TProfile("fFakesPhi","Number of ghosts vs phi",nbin,0,6.28,0,100);
 243   fSlicesPt = new TProfile("fSlicesPt","Number of slices vs pt",nbin,xlow,xup,0,100);
 244   fSlicesPhi = new TProfile("fSlicesPhi","Number of slices vs phi",nbin,0,6.28,0,100);
 245   fPtRes2 = new TH1F("fPtRes2","Relative Pt resolution (Pt>2 GeV)",20,-100.,100.);
 246   fEtaRes = new TH1F("fEtaRes","Eta resolution",30,-0.1,0.1);
 247   fPsiRes = new TH1F("fPsiRes","Psi resolution",30,-0.1,0.1);
 248
 249   fNGoodTracksRad = new TH1F("fNGoodTracksRad","Good tracks vs distance(x,y)",1,0,0.2);
 250   fNFoundTracksRad = new TH1F("fNFoundTracksRad","Found tracks vs distance(x,y)",1,0,0.2);
 251   fRadEff = new TH1F("fRadEff","Tracking efficiency vs distance(x,y)",1,0,0.2);
 252   fNGoodTracksZ = new TH1F("fNGoodTracksZ","Good tracks vs distance(z)",1,0,0.2);
 253   fNFoundTracksZ = new TH1F("fNFoundTracksZ","Found tracks vs distance(z)",1,0,0.2);
 254   fZEff = new TH1F("fZEff","Tracking efficiency vs distance(z)",1,0,0.2);
 255 }
 256
 257 void Write2File()
 258 {
 259   //Write histograms to file:
 260
 261   TFile *of = TFile::Open("hough_efficiencies.root","RECREATE");
 262   if(!of->IsOpen())
 263     {
 264       cout<<"Problems opening rootfile"<<endl;
 265       return;
 266     }
 267
 268   of->cd();
 269   fNtuppel->Write();
 270   fPtRes->Write();
 271   fNGoodTracksPt->Write();
 272   fNFoundTracksPt->Write();
 273   fNFakeTracksPt->Write();
 274   fTrackEffPt->Write();
 275   fFakeTrackEffPt->Write();
 276   fNGoodTracksEta->Write();
 277   fNFoundTracksEta->Write();
 278   fNFakeTracksEta->Write();
 279   fTrackEffEta->Write();
 280   fFakeTrackEffEta->Write();
 281
 282   fNGoodTracksPhi->Write();
 283   fNFoundTracksPhi->Write();
 284   fNFakeTracksPhi->Write();
 285   fTrackEffPhi->Write();
 286   fFakeTrackEffPhi->Write();
 287
 288   fFakesPt->Write();
 289   fFakesPhi->Write();
 290   fSlicesPt->Write();
 291   fSlicesPhi->Write();
 292
 293   fPtRes2->Write();
 294   fEtaRes->Write();
 295   fPsiRes->Write();
 296
 297   fNGoodTracksRad->Write();
 298   fNFoundTracksRad->Write();
 299   fRadEff->Write();
 300   fNGoodTracksZ->Write();
 301   fNFoundTracksZ->Write();
 302   fZEff->Write();
 303
 304   of->Close();
 305 }
 306
 307 void CalcEffHistos()
 308 {
 309   fNFoundTracksPt->Sumw2(); fNGoodTracksPt->Sumw2();
 310   fTrackEffPt->Divide(fNFoundTracksPt,fNGoodTracksPt,1,1.,"b");
 311   fFakeTrackEffPt->Divide(fNFakeTracksPt,fNGoodTracksPt,1,1.,"b");
 312   fTrackEffPt->SetMaximum(1.4);
 313   fTrackEffPt->SetXTitle("P_{T} [GeV]");
 314   fTrackEffPt->SetLineWidth(2);
 315   fFakeTrackEffPt->SetFillStyle(3013);
 316   fTrackEffPt->SetLineColor(4);
 317   fFakeTrackEffPt->SetFillColor(2);
 318
 319   fNFoundTracksEta->Sumw2(); fNGoodTracksEta->Sumw2();
 320   fTrackEffEta->Divide(fNFoundTracksEta,fNGoodTracksEta,1,1.,"b");
 321   fFakeTrackEffEta->Divide(fNFakeTracksEta,fNGoodTracksEta,1,1.,"b");
 322   fTrackEffEta->SetMaximum(1.4);
 323   fTrackEffEta->SetXTitle("#eta");
 324   fTrackEffEta->SetLineWidth(2);
 325   fFakeTrackEffEta->SetFillStyle(3013);
 326   fTrackEffEta->SetLineColor(4);
 327   fFakeTrackEffEta->SetFillColor(2);
 328
 329   fNFoundTracksPhi->Sumw2(); fNGoodTracksPhi->Sumw2();
 330   fTrackEffPhi->Divide(fNFoundTracksPhi,fNGoodTracksPhi,1,1.,"b");
 331   fFakeTrackEffPhi->Divide(fNFakeTracksPhi,fNGoodTracksPhi,1,1.,"b");
 332   fTrackEffPhi->SetMaximum(1.4);
 333   fTrackEffPhi->SetXTitle("#psi [rad]");
 334   fTrackEffPhi->SetLineWidth(2);
 335   fFakeTrackEffPhi->SetFillStyle(3013);
 336   fTrackEffPhi->SetLineColor(4);
 337   fFakeTrackEffPhi->SetFillColor(2);
 338
 339   fEtaRes->SetXTitle("#eta");
 340   fPsiRes->SetXTitle("#psi [rad]");
 341
 342   fRadEff->Divide(fNFoundTracksRad,fNGoodTracksRad,1,1.,"b");
 343   fRadEff->SetXTitle("Distance [cm]");
 344   fZEff->Divide(fNFoundTracksZ,fNGoodTracksZ,1,1.,"b");
 345   fZEff->SetXTitle("Distance [cm]");
 346 }
 347
 348 void plotptres(Float_t ptmin)
 349 {
 350   //Plot the relative pt resolution vs pt.
 351
 352   const Int_t n=8;
 353
 354   Double_t hltx[9] = {0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,2.0};
 355   Double_t hltxerr[n];
 356   Double_t hltavx[n];
 357   Double_t hlty[n];
 358   Double_t hltyerr[n];
 359   Char_t string[1024];
 360
 361   for(Int_t i=0; i<n; i++)
 362     {
 363       TH1F *hist = new TH1F("hist","",100,-50.*(float)(i+1),50.*(float)(i+1));
 364       sprintf(string,"pt_gen > %f && pt_gen <= %f",hltx[i],hltx[i+1]);
 365       fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>hist",string,"goff");
 366       TF1 *f1 = new TF1("f1","gaus");
 367       hist->Fit("f1","E");
 368       hlty[i] = f1->GetParameter(2);
 369       hltyerr[i] = f1->GetParError(2);
 370       //hlty[i] = hist->GetRMS();
 371       //hltyerr[i] = 0;
 372       hltavx[i] = (hltx[i]+hltx[i+1])/2.0;
 373       hltxerr[i] = (hltx[i+1]-hltx[i])/2.0;
 374       cout<<string<<" "<<i<<" "<<hlty[i]<<" "<<hltyerr[i]<<endl;
 375       delete f1;
 376       delete hist;
 377     }
 378
 379   sprintf(string,"pt_gen > %f",2.0);
 380   fNtuppel->Draw("(pt_found-pt_gen)/pt_gen*100>>fPtRes2",string,"goff");
 381   TF1 *f1 = new TF1("f1","gaus");
 382   fPtRes2->Fit("f1","E");
 383   fPtRes2->SetXTitle("#Delta P_{t} / P_{t} [%]");
 384   fPtRes2->SetYTitle("#Delta N");
 385   delete f1;
 386
 387   TGraphErrors *gr1 = new TGraphErrors(n,hltavx,hlty,hltxerr,hltyerr);
 388   TCanvas *c1 = new TCanvas("c1","",1);
 389   fPtRes = c1->DrawFrame(0.1,0,2,15);
 390   gr1->SetTitle("");
 391   gr1->Draw("P");
 392   gr1->SetMarkerStyle(20);
 393   gr1->SetLineWidth(2);
 394   gr1->SetMarkerSize(1.3);
 395   f1 = new TF1("f1","pol1",ptmin,2.0);
 396   gr1->Fit("f1","ER");
 397   fPtRes->SetXTitle("p_{t} [GeV]");
 398   fPtRes->SetYTitle("#Delta P_{t} / P_{t} [%]");
 399   delete f1;
 400 }