TRD/qaRec/AliTRDtrackingResolution.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercialf purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  TRD tracking resolution                                               //
  21 //
  22 // The class performs resolution and residual studies
  23 // of the TRD tracks for the following quantities :
  24 //   - spatial position (y, [z])
  25 //   - angular (phi) tracklet
  26 //   - momentum at the track level
  27 //
  28 // The class has to be used for regular detector performance checks using the official macros:
  29 //   - $ALICE_ROOT/TRD/qaRec/run.C
  30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
  31 //
  32 // For stand alone usage please refer to the following example:
  33 // {
  34 //   gSystem->Load("libANALYSIS.so");
  35 //   gSystem->Load("libTRDqaRec.so");
  36 //   AliTRDtrackingResolution *res = new AliTRDtrackingResolution();
  37 //   //res->SetMCdata();
  38 //   //res->SetVerbose();
  39 //   //res->SetVisual();
  40 //   res->Load("TRD.TaskResolution.root");
  41 //   if(!res->PostProcess()) return;
  42 //   res->GetRefFigure(0);
  43 // }
  44 //
  45 //  Authors:                                                              //
  46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
  47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  48 //                                                                        //
  49 ////////////////////////////////////////////////////////////////////////////
  50
  51 #include <cstring>
  52
  53 #include <TSystem.h>
  54 #include <TObjArray.h>
  55 #include <TH2.h>
  56 #include <TH1.h>
  57 #include <TF1.h>
  58 #include <TCanvas.h>
  59 #include <TBox.h>
  60 #include <TProfile.h>
  61 #include <TGraphErrors.h>
  62 #include <TMath.h>
  63 #include "TTreeStream.h"
  64 #include "TGeoManager.h"
  65
  66 #include "AliAnalysisManager.h"
  67 #include "AliTrackReference.h"
  68 #include "AliTrackPointArray.h"
  69 #include "AliCDBManager.h"
  70
  71 #include "AliTRDSimParam.h"
  72 #include "AliTRDgeometry.h"
  73 #include "AliTRDpadPlane.h"
  74 #include "AliTRDcluster.h"
  75 #include "AliTRDseedV1.h"
  76 #include "AliTRDtrackV1.h"
  77 #include "AliTRDtrackerV1.h"
  78 #include "AliTRDReconstructor.h"
  79 #include "AliTRDrecoParam.h"
  80
  81 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
  82 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  83 #include "AliTRDtrackingResolution.h"
  84
  85 ClassImp(AliTRDtrackingResolution)
  86
  87 //________________________________________________________
  88 AliTRDtrackingResolution::AliTRDtrackingResolution()
  89   :AliTRDrecoTask("Resolution", "Tracking Resolution")
  90   ,fStatus(0)
  91   ,fReconstructor(0x0)
  92   ,fGeo(0x0)
  93   ,fGraphS(0x0)
  94   ,fGraphM(0x0)
  95   ,fClResiduals(0x0)
  96   ,fClResolution(0x0)
  97   ,fTrkltResolution(0x0)
  98 {
  99   fReconstructor = new AliTRDReconstructor();
 100   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 101   fGeo = new AliTRDgeometry();
 102
 103   InitFunctorList();
 104
 105   DefineOutput(1+kClusterResidual, TObjArray::Class());
 106   DefineOutput(1+kClusterResolution, TObjArray::Class());
 107   DefineOutput(1+kTrackletYResolution, TObjArray::Class());
 108 }
 109
 110 //________________________________________________________
 111 AliTRDtrackingResolution::~AliTRDtrackingResolution()
 112 {
 113   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 114   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 115   delete fGeo;
 116   delete fReconstructor;
 117   if(gGeoManager) delete gGeoManager;
 118   if(fClResiduals){fClResiduals->Delete(); delete fClResiduals;}
 119   if(fClResolution){fClResolution->Delete(); delete fClResolution;}
 120   if(fTrkltResolution){fTrkltResolution->Delete(); delete fTrkltResolution;}
 121 }
 122
 123
 124 //________________________________________________________
 125 void AliTRDtrackingResolution::CreateOutputObjects()
 126 {
 127   // spatial resolution
 128   OpenFile(0, "RECREATE");
 129
 130   fContainer = Histos();
 131
 132   fClResiduals = new TObjArray();
 133   fClResiduals->SetOwner(kTRUE);
 134   fClResolution = new TObjArray();
 135   fClResolution->SetOwner(kTRUE);
 136   fTrkltResolution = new TObjArray();
 137   fTrkltResolution->SetOwner(kTRUE);
 138 }
 139
 140 //________________________________________________________
 141 void AliTRDtrackingResolution::Exec(Option_t *opt)
 142 {
 143   fClResiduals->Delete();
 144   fClResolution->Delete();
 145   fTrkltResolution->Delete();
 146
 147   AliTRDrecoTask::Exec(opt);
 148
 149   PostData(1+kClusterResidual, fClResiduals);
 150   PostData(1+kClusterResolution, fClResolution);
 151   PostData(1+kTrackletYResolution, fTrkltResolution);
 152 }
 153
 154 //________________________________________________________
 155 TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
 156 {
 157   if(track) fTrack = track;
 158   if(!fTrack){
 159     AliWarning("No Track defined.");
 160     return 0x0;
 161   }
 162   TH1 *h = 0x0;
 163   if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
 164     AliWarning("No output histogram defined.");
 165     return 0x0;
 166   }
 167
 168   Int_t pdg = (HasMCdata() && fMC) ? fMC->GetPDG() : 0;
 169   UChar_t s;
 170   Float_t x0, y0, z0, dy, dydx, dzdx;
 171   AliTRDseedV1 *fTracklet = 0x0;
 172   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 173     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 174     if(!fTracklet->IsOK()) continue;
 175     x0 = fTracklet->GetX0();
 176
 177     // retrive the track angle with the chamber
 178     if(HasMCdata() && fMC){
 179       if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 180     }else{
 181       y0   = fTracklet->GetYref(0);
 182       z0   = fTracklet->GetZref(0);
 183       dydx = fTracklet->GetYref(1);
 184       dzdx = fTracklet->GetZref(1);
 185     }
 186
 187     AliTRDseedV1 trklt(*fTracklet);
 188     if(!trklt.Fit(kFALSE)) continue;
 189
 190     AliTRDcluster *c = 0x0;
 191     fTracklet->ResetClusterIter(kFALSE);
 192     while((c = fTracklet->PrevCluster())){
 193       dy = trklt.GetYat(c->GetX()) - c->GetY();
 194       h->Fill(dydx, dy);
 195
 196       if(fDebugLevel>=1){
 197         Float_t q = c->GetQ();
 198         // Get z-position with respect to anode wire
 199         AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 200         Int_t det = c->GetDetector();
 201         Float_t x = c->GetX();
 202         Float_t z = z0-(x0-x)*dzdx;
 203         Int_t istk = fGeo->GetStack(det);
 204         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 205         Float_t row0 = pp->GetRow0();
 206         Float_t d  =  row0 - z + simParam->GetAnodeWireOffset();
 207         d -= ((Int_t)(2 * d)) / 2.0;
 208         if (d > 0.25) d  = 0.5 - d;
 209
 210         (*fDebugStream) << "ClusterResidual"
 211           << "ly="   << ily
 212           << "stk="  << istk
 213           << "pdg="  << pdg
 214           << "dydx=" << dydx
 215           << "dzdx=" << dzdx
 216           << "q="    << q
 217           << "x="    << x
 218           << "z="    << z
 219           << "d="    << d
 220           << "dy="   << dy
 221           << "\n";
 222
 223
 224       }
 225     }
 226   }
 227   return h;
 228 }
 229
 230 //________________________________________________________
 231 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
 232 {
 233   if(!HasMCdata()){
 234     AliWarning("No MC defined. Results will not be available.");
 235     return 0x0;
 236   }
 237   if(track) fTrack = track;
 238   if(!fTrack){
 239     AliWarning("No Track defined.");
 240     return 0x0;
 241   }
 242   TH1 *h = 0x0;
 243   if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
 244     AliWarning("No output histogram defined.");
 245     return 0x0;
 246   }
 247   if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
 248     AliWarning("No output histogram defined.");
 249     return 0x0;
 250   }
 251   if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
 252     AliWarning("No output histogram defined.");
 253     return 0x0;
 254   }
 255   if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
 256     AliWarning("No output histogram defined.");
 257     return 0x0;
 258   }
 259   //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
 260   UChar_t s;
 261   Int_t pdg = fMC->GetPDG(), det=-1;
 262   Float_t x0, y0, z0, dy, dydx, dzdx;
 263   AliTRDseedV1 *fTracklet = 0x0;
 264   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 265     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 266     if(!fTracklet->IsOK()) continue;
 267     //printf("process layer[%d]\n", ily);
 268     // retrive the track position and direction within the chamber
 269     det = fTracklet->GetDetector();
 270     x0  = fTracklet->GetX0();
 271     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 272
 273     // recalculate tracklet based on the MC info
 274     AliTRDseedV1 tt(*fTracklet);
 275     tt.SetZref(0, z0);
 276     tt.SetZref(1, dzdx);
 277     if(!tt.Fit(kFALSE)) continue;
 278     //tt.Update();
 279     dy = tt.GetYfit(0) - y0;
 280     Float_t dphi   = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
 281     Float_t dz = 100.;
 282     Bool_t  cross = fTracklet->GetNChange();
 283     if(cross){
 284       Double_t *xyz = tt.GetCrossXYZ();
 285       dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
 286       ((TH2I*)fContainer->At(kTrackletZResolution))->Fill(dzdx, dz);
 287     }
 288
 289     // Fill Histograms
 290     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
 291     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
 292
 293     // Fill Debug stream
 294     if(fDebugLevel>=1){
 295       Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
 296       (*fDebugStream) << "TrkltResolution"
 297         << "det="                 << det
 298         << "pdg="     << pdg
 299         << "p="       << p
 300         << "ymc="     << y0
 301         << "zmc="     << z0
 302         << "dydx="    << dydx
 303         << "dzdx="    << dzdx
 304         << "cross="   << cross
 305         << "dy="                  << dy
 306         << "dz="                  << dz
 307         << "dphi="              << dphi
 308         << "\n";
 309     }
 310
 311     Int_t istk = AliTRDgeometry::GetStack(det);
 312     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 313     Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
 314     Float_t tilt = fTracklet->GetTilt();
 315
 316     AliTRDcluster *c = 0x0;
 317     fTracklet->ResetClusterIter(kFALSE);
 318     while((c = fTracklet->PrevCluster())){
 319       Float_t  q = TMath::Abs(c->GetQ());
 320       Float_t xc = c->GetX();
 321       Float_t yc = c->GetY();
 322       Float_t zc = c->GetZ();
 323       Float_t dx = x0 - xc;
 324       Float_t yt = y0 - dx*dydx;
 325       Float_t zt = z0 - dx*dzdx;
 326       dy = yt - (yc - tilt*(zc-zt));
 327
 328       // Fill Histograms
 329       if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
 330
 331       // Fill Debug Tree
 332       if(fDebugLevel>=1){
 333         Float_t d = zr0 - zt;
 334         d -= ((Int_t)(2 * d)) / 2.0;
 335         if (d > 0.25) d  = 0.5 - d;
 336         Int_t label = fMC->GetLabel();
 337
 338         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 339         fClResolution->Add(clInfo);
 340         clInfo->SetCluster(c);
 341         clInfo->SetMC(pdg, label);
 342         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 343         clInfo->SetResolution(dy);
 344         clInfo->SetAnisochronity(d);
 345         clInfo->SetDriftLength(x0-xc);
 346         //clInfo->Print();
 347         (*fDebugStream) << "ClusterResolution"
 348           <<"clInfo.=" << clInfo
 349           << "\n";
 350       }
 351     }
 352   }
 353   return h;
 354 }
 355
 356
 357 //________________________________________________________
 358 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
 359 {
 360   TBox *b = 0x0;
 361   TAxis *ax = 0x0;
 362   TGraphErrors *g = 0x0;
 363   switch(ifig){
 364   case kClusterResidual:
 365     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 366     g->Draw("apl");
 367     ax = g->GetHistogram()->GetYaxis();
 368     ax->SetRangeUser(-.5, 1.);
 369     ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
 370     ax = g->GetHistogram()->GetXaxis();
 371     ax->SetTitle("tg(#phi)");
 372     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 373     g->Draw("pl");
 374     b = new TBox(-.15, -.5, .15, 1.);
 375     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 376     b->SetLineColor(0); b->Draw();
 377     return;
 378   case kClusterResolution:
 379     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 380     ax = g->GetHistogram()->GetYaxis();
 381     ax->SetRangeUser(-.5, 1.);
 382     ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
 383     ax = g->GetHistogram()->GetXaxis();
 384     ax->SetTitle("tg(#phi)");
 385     g->Draw("apl");
 386     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 387     g->Draw("pl");
 388     b = new TBox(-.15, -.5, .15, 1.);
 389     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 390     b->SetLineColor(0); b->Draw();
 391     return;
 392   case kTrackletYResolution:
 393     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 394     ax = g->GetHistogram()->GetYaxis();
 395     ax->SetRangeUser(-.5, 1.);
 396     ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
 397     ax = g->GetHistogram()->GetXaxis();
 398     ax->SetTitle("tg(#phi)");
 399     g->Draw("apl");
 400     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 401     g->Draw("pl");
 402     return;
 403   case kTrackletZResolution:
 404     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 405     ax = g->GetHistogram()->GetYaxis();
 406     ax->SetRangeUser(-.5, 1.);
 407     ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
 408     ax = g->GetHistogram()->GetXaxis();
 409     ax->SetTitle("tg(#theta)");
 410     g->Draw("apl");
 411     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 412     g->Draw("pl");
 413     return;
 414   case kTrackletAngleResolution:
 415     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 416     ax = g->GetHistogram()->GetYaxis();
 417     ax->SetRangeUser(-.05, .2);
 418     ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
 419     ax = g->GetHistogram()->GetXaxis();
 420     ax->SetTitle("tg(#phi)");
 421     g->Draw("apl");
 422     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 423     g->Draw("pl");
 424     return;
 425   default:
 426     AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
 427     return;
 428   }
 429   AliInfo(Form("Reference plot [%d] missing result", ifig));
 430 }
 431
 432
 433 //________________________________________________________
 434 Bool_t AliTRDtrackingResolution::PostProcess()
 435 {
 436   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 437   if (!fContainer) {
 438     Printf("ERROR: list not available");
 439     return kFALSE;
 440   }
 441   fNRefFigures = fContainer->GetEntriesFast();
 442   if(!fGraphS){
 443     fGraphS = new TObjArray(fNRefFigures);
 444     fGraphS->SetOwner();
 445   }
 446   if(!fGraphM){
 447     fGraphM = new TObjArray(fNRefFigures);
 448     fGraphM->SetOwner();
 449   }
 450
 451   TH2I *h2 = 0x0;
 452   TH1D *h = 0x0;
 453   TGraphErrors *gm = 0x0, *gs = 0x0;
 454
 455   // define models
 456   TF1 f("f1", "gaus", -.5, .5);
 457
 458   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 459
 460   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 461
 462   TCanvas *c = 0x0;
 463   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
 464   char opt[5];
 465   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
 466
 467
 468   //PROCESS RESIDUAL DISTRIBUTIONS
 469
 470   // Clusters residuals
 471   h2 = (TH2I *)(fContainer->At(kClusterResidual));
 472   gm = new TGraphErrors();
 473   gm->SetLineColor(kBlue);
 474   gm->SetMarkerStyle(7);
 475   gm->SetMarkerColor(kBlue);
 476   gm->SetNameTitle("clm", "");
 477   fGraphM->AddAt(gm, kClusterResidual);
 478   gs = new TGraphErrors();
 479   gs->SetLineColor(kRed);
 480   gs->SetMarkerStyle(23);
 481   gs->SetMarkerColor(kRed);
 482   gs->SetNameTitle("cls", "");
 483   fGraphS->AddAt(gs, kClusterResidual);
 484   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 485     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 486     h = h2->ProjectionY("py", ibin, ibin);
 487     if(h->GetEntries()<100) continue;
 488     AdjustF1(h, &f);
 489
 490     if(IsVisual()){c->cd(); c->SetLogy();}
 491     h->Fit(&f, opt, "", -0.5, 0.5);
 492     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 493
 494     Int_t ip = gm->GetN();
 495     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 496     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 497     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 498     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 499   }
 500
 501
 502   //PROCESS RESOLUTION DISTRIBUTIONS
 503
 504   if(HasMCdata()){
 505     // cluster y resolution
 506     h2 = (TH2I*)fContainer->At(kClusterResolution);
 507     gm = new TGraphErrors();
 508     gm->SetLineColor(kBlue);
 509     gm->SetMarkerStyle(7);
 510     gm->SetMarkerColor(kBlue);
 511     gm->SetNameTitle("clym", "");
 512     fGraphM->AddAt(gm, kClusterResolution);
 513     gs = new TGraphErrors();
 514     gs->SetLineColor(kRed);
 515     gs->SetMarkerStyle(23);
 516     gs->SetMarkerColor(kRed);
 517     gs->SetNameTitle("clys", "");
 518     fGraphS->AddAt(gs, kClusterResolution);
 519     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 520       h = h2->ProjectionY("py", iphi, iphi);
 521       if(h->GetEntries()<100) continue;
 522       AdjustF1(h, &f);
 523
 524       if(IsVisual()){c->cd(); c->SetLogy();}
 525       h->Fit(&f, opt, "", -0.5, 0.5);
 526       if(IsVerbose()){
 527         printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
 528       }
 529       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 530
 531       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 532       Int_t ip = gm->GetN();
 533       gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 534       gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 535       gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 536       gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 537     }
 538
 539     // tracklet y resolution
 540     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 541     gm = new TGraphErrors();
 542     gm->SetLineColor(kBlue);
 543     gm->SetMarkerStyle(7);
 544     gm->SetMarkerColor(kBlue);
 545     gm->SetNameTitle("trkltym", "");
 546     fGraphM->AddAt(gm, kTrackletYResolution);
 547     gs = new TGraphErrors();
 548     gs->SetLineColor(kRed);
 549     gs->SetMarkerStyle(23);
 550     gs->SetMarkerColor(kRed);
 551     gs->SetNameTitle("trkltys", "");
 552     fGraphS->AddAt(gs, kTrackletYResolution);
 553     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 554       h = h2->ProjectionY("py", iphi, iphi);
 555       if(h->GetEntries()<100) continue;
 556       AdjustF1(h, &fb);
 557
 558       if(IsVisual()){c->cd(); c->SetLogy();}
 559       h->Fit(&fb, opt, "", -0.5, 0.5);
 560       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 561
 562       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 563       Int_t ip = gm->GetN();
 564       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 565       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 566       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 567       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 568     }
 569
 570     // tracklet z resolution
 571     h2 = (TH2I*)fContainer->At(kTrackletZResolution);
 572     gm = new TGraphErrors();
 573     gm->SetLineColor(kBlue);
 574     gm->SetMarkerStyle(7);
 575     gm->SetMarkerColor(kBlue);
 576     gm->SetNameTitle("trkltzm", "");
 577     fGraphM->AddAt(gm, kTrackletZResolution);
 578     gs = new TGraphErrors();
 579     gs->SetLineColor(kRed);
 580     gs->SetMarkerStyle(23);
 581     gs->SetMarkerColor(kRed);
 582     gs->SetNameTitle("trkltzs", "");
 583     fGraphS->AddAt(gs, kTrackletZResolution);
 584     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 585       h = h2->ProjectionY("py", iphi, iphi);
 586       if(h->GetEntries()<100) continue;
 587       AdjustF1(h, &fb);
 588
 589       if(IsVisual()){c->cd(); c->SetLogy();}
 590       h->Fit(&fb, opt, "", -0.5, 0.5);
 591       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 592
 593       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 594       Int_t ip = gm->GetN();
 595       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 596       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 597       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 598       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 599     }
 600
 601     // tracklet phi resolution
 602     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 603     gm = new TGraphErrors();
 604     gm->SetLineColor(kBlue);
 605     gm->SetMarkerStyle(7);
 606     gm->SetMarkerColor(kBlue);
 607     gm->SetNameTitle("trkltym", "");
 608     fGraphM->AddAt(gm, kTrackletAngleResolution);
 609     gs = new TGraphErrors();
 610     gs->SetLineColor(kRed);
 611     gs->SetMarkerStyle(23);
 612     gs->SetMarkerColor(kRed);
 613     gs->SetNameTitle("trkltys", "");
 614     fGraphS->AddAt(gs, kTrackletAngleResolution);
 615     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 616       h = h2->ProjectionY("py", iphi, iphi);
 617       if(h->GetEntries()<100) continue;
 618
 619       if(IsVisual()){c->cd(); c->SetLogy();}
 620       h->Fit(&f, opt, "", -0.5, 0.5);
 621       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 622
 623       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 624       Int_t ip = gm->GetN();
 625       gm->SetPoint(ip, phi, f.GetParameter(1));
 626       gm->SetPointError(ip, 0., f.GetParError(1));
 627       gs->SetPoint(ip, phi, f.GetParameter(2));
 628       gs->SetPointError(ip, 0., f.GetParError(2));
 629     }
 630   }
 631   if(c) delete c;
 632
 633   return kTRUE;
 634 }
 635
 636
 637 //________________________________________________________
 638 void AliTRDtrackingResolution::Terminate(Option_t *)
 639 {
 640   if(fDebugStream){
 641     delete fDebugStream;
 642     fDebugStream = 0x0;
 643     fDebugLevel = 0;
 644   }
 645   if(HasPostProcess()) PostProcess();
 646 }
 647
 648 //________________________________________________________
 649 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
 650 {
 651 // Helper function to avoid duplication of code
 652 // Make first guesses on the fit parameters
 653
 654   // find the intial parameters of the fit !! (thanks George)
 655   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 656   Double_t sum = 0.;
 657   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 658   f->SetParLimits(0, 0., 3.*sum);
 659   f->SetParameter(0, .9*sum);
 660
 661   f->SetParLimits(1, -.2, .2);
 662   f->SetParameter(1, -0.1);
 663
 664   f->SetParLimits(2, 0., 4.e-1);
 665   f->SetParameter(2, 2.e-2);
 666   if(f->GetNpar() <= 4) return;
 667
 668   f->SetParLimits(3, 0., sum);
 669   f->SetParameter(3, .1*sum);
 670
 671   f->SetParLimits(4, -.3, .3);
 672   f->SetParameter(4, 0.);
 673
 674   f->SetParLimits(5, 0., 1.e2);
 675   f->SetParameter(5, 2.e-1);
 676 }
 677
 678 //________________________________________________________
 679 TObjArray* AliTRDtrackingResolution::Histos()
 680 {
 681   if(fContainer) return fContainer;
 682
 683   fContainer  = new TObjArray(5);
 684
 685   TH1 *h = 0x0;
 686   // cluster to tracklet residuals [2]
 687   fContainer->AddAt(h = new TH2I("fYCl", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5), kClusterResidual);
 688   h->GetXaxis()->SetTitle("tg(#phi)");
 689   h->GetYaxis()->SetTitle("#Delta y [cm]");
 690   h->GetZaxis()->SetTitle("entries");
 691 //   // tracklet to Riemann fit residuals [2]
 692 //   fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
 693 //   fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
 694 //   fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
 695 //   fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
 696
 697   // Resolution histos
 698   if(HasMCdata()){
 699     // cluster y resolution [0]
 700     fContainer->AddAt(h = new TH2I("fYClMC", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5), kClusterResolution);
 701     h->GetXaxis()->SetTitle("tg(#phi)");
 702     h->GetYaxis()->SetTitle("#Delta y [cm]");
 703     h->GetZaxis()->SetTitle("entries");
 704     // tracklet y resolution [0]
 705     fContainer->AddAt(h = new TH2I("fYTrkltMC", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5), kTrackletYResolution);
 706     h->GetXaxis()->SetTitle("tg(#phi)");
 707     h->GetYaxis()->SetTitle("#Delta y [cm]");
 708     h->GetZaxis()->SetTitle("entries");
 709     // tracklet y resolution [0]
 710     fContainer->AddAt(h = new TH2I("fZTrkltMC", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5), kTrackletZResolution);
 711     h->GetXaxis()->SetTitle("tg(#theta)");
 712     h->GetYaxis()->SetTitle("#Delta z [cm]");
 713     h->GetZaxis()->SetTitle("entries");
 714     // tracklet angular resolution [1]
 715     fContainer->AddAt(h = new TH2I("fPhiTrkltMC", "Tracklet Resolution (Angular)", 31, -.48, .48, 100, -10., 10.), kTrackletAngleResolution);
 716     h->GetXaxis()->SetTitle("tg(#phi)");
 717     h->GetYaxis()->SetTitle("#Delta #phi [deg]");
 718     h->GetZaxis()->SetTitle("entries");
 719
 720 //     // Riemann track resolution [y, z, angular]
 721 //     fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
 722 //     fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
 723 //     fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
 724 //
 725 //     Kalman track resolution [y, z, angular]
 726 //     fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
 727 //     fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
 728 //     fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
 729   }
 730   return fContainer;
 731 }
 732
 733
 734 //________________________________________________________
 735 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 736 {
 737
 738   fReconstructor->SetRecoParam(r);
 739 }