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   ,fTrkltResiduals(0x0)
  97   ,fTrkltPhiResiduals(0x0)
  98   ,fClResolution(0x0)
  99   ,fTrkltResolution(0x0)
 100 {
 101   fReconstructor = new AliTRDReconstructor();
 102   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 103   fGeo = new AliTRDgeometry();
 104
 105   InitFunctorList();
 106
 107   DefineOutput(1+kClusterResidual, TObjArray::Class());
 108   DefineOutput(1+kTrackletYResidual, TObjArray::Class());
 109   DefineOutput(1+kTrackletPhiResidual, TObjArray::Class());
 110   DefineOutput(1+kClusterResolution, TObjArray::Class());
 111   DefineOutput(1+kTrackletYResolution, TObjArray::Class());
 112 }
 113
 114 //________________________________________________________
 115 AliTRDtrackingResolution::~AliTRDtrackingResolution()
 116 {
 117   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 118   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 119   delete fGeo;
 120   delete fReconstructor;
 121   if(gGeoManager) delete gGeoManager;
 122   if(fClResiduals){fClResiduals->Delete(); delete fClResiduals;}
 123   if(fTrkltResiduals){fTrkltResiduals->Delete(); delete fTrkltResiduals;}
 124   if(fTrkltPhiResiduals){fTrkltPhiResiduals->Delete(); delete fTrkltPhiResiduals;}
 125   if(fClResolution){fClResolution->Delete(); delete fClResolution;}
 126   if(fTrkltResolution){fTrkltResolution->Delete(); delete fTrkltResolution;}
 127 }
 128
 129
 130 //________________________________________________________
 131 void AliTRDtrackingResolution::CreateOutputObjects()
 132 {
 133   // spatial resolution
 134   OpenFile(0, "RECREATE");
 135
 136   fContainer = Histos();
 137
 138   fClResiduals = new TObjArray();
 139   fClResiduals->SetOwner(kTRUE);
 140   fTrkltResiduals = new TObjArray();
 141   fTrkltResiduals->SetOwner(kTRUE);
 142   fTrkltPhiResiduals = new TObjArray();
 143   fTrkltPhiResiduals->SetOwner(kTRUE);
 144   fClResolution = new TObjArray();
 145   fClResolution->SetOwner(kTRUE);
 146   fTrkltResolution = new TObjArray();
 147   fTrkltResolution->SetOwner(kTRUE);
 148 }
 149
 150 //________________________________________________________
 151 void AliTRDtrackingResolution::Exec(Option_t *opt)
 152 {
 153   fClResiduals->Delete();
 154   fTrkltResiduals->Delete();
 155   fTrkltPhiResiduals->Delete();
 156   fClResolution->Delete();
 157   fTrkltResolution->Delete();
 158
 159   AliTRDrecoTask::Exec(opt);
 160
 161   PostData(1+kClusterResidual, fClResiduals);
 162   PostData(1+kTrackletYResidual, fTrkltResiduals);
 163   PostData(1+kTrackletPhiResidual, fTrkltPhiResiduals);
 164   PostData(1+kClusterResolution, fClResolution);
 165   PostData(1+kTrackletYResolution, fTrkltResolution);
 166 }
 167
 168 //________________________________________________________
 169 TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
 170 {
 171   if(track) fTrack = track;
 172   if(!fTrack){
 173     AliWarning("No Track defined.");
 174     return 0x0;
 175   }
 176   TH1 *h = 0x0;
 177   if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
 178     AliWarning("No output histogram defined.");
 179     return 0x0;
 180   }
 181
 182   Float_t x0, y0, z0, dy, dydx, dzdx;
 183   AliTRDseedV1 *fTracklet = 0x0;
 184   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 185     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 186     if(!fTracklet->IsOK()) continue;
 187     x0 = fTracklet->GetX0();
 188
 189     // retrive the track angle with the chamber
 190     y0   = fTracklet->GetYref(0);
 191     z0   = fTracklet->GetZref(0);
 192     dydx = fTracklet->GetYref(1);
 193     dzdx = fTracklet->GetZref(1);
 194     Float_t tilt = fTracklet->GetTilt();
 195     AliTRDcluster *c = 0x0;
 196     fTracklet->ResetClusterIter(kFALSE);
 197     while((c = fTracklet->PrevCluster())){
 198       Float_t xc = c->GetX();
 199       Float_t yc = c->GetY();
 200       Float_t zc = c->GetZ();
 201       Float_t dx = x0 - xc;
 202       Float_t yt = y0 - dx*dydx;
 203       Float_t zt = z0 - dx*dzdx;
 204       dy = yt - (yc - tilt*(zc-zt));
 205
 206       //dy = trklt.GetYat(c->GetX()) - c->GetY();
 207       h->Fill(dydx, dy);
 208
 209       if(fDebugLevel>=1){
 210         // Get z-position with respect to anode wire
 211         AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 212         Int_t istk = fGeo->GetStack(c->GetDetector());
 213         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 214         Float_t row0 = pp->GetRow0();
 215         Float_t d  =  row0 - zt + simParam->GetAnodeWireOffset();
 216         d -= ((Int_t)(2 * d)) / 2.0;
 217         if (d > 0.25) d  = 0.5 - d;
 218
 219         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 220         fClResiduals->Add(clInfo);
 221         clInfo->SetCluster(c);
 222         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 223         clInfo->SetResolution(dy);
 224         clInfo->SetAnisochronity(d);
 225         clInfo->SetDriftLength(dx);
 226         (*fDebugStream) << "ClusterResiduals"
 227           <<"clInfo.=" << clInfo
 228           << "\n";
 229       }
 230     }
 231   }
 232   return h;
 233 }
 234
 235
 236 //________________________________________________________
 237 TH1* AliTRDtrackingResolution::PlotTrackletResiduals(const AliTRDtrackV1 *track)
 238 {
 239   if(track) fTrack = track;
 240   if(!fTrack){
 241     AliWarning("No Track defined.");
 242     return 0x0;
 243   }
 244   TH1 *h = 0x0;
 245   if(!(h = ((TH2I*)fContainer->At(kTrackletYResidual)))){
 246     AliWarning("No output histogram defined.");
 247     return 0x0;
 248   }
 249
 250   Float_t  dydx, dzdx, yref, zref, yfit, zfit, x0, xmean;
 251   AliTRDseedV1 *fTracklet = 0x0;
 252   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 253     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 254     if(!fTracklet->IsOK()) continue;
 255
 256     AliTRDseedV1 tt(*fTracklet);
 257     // tt.SetX0(fTracklet->GetX0()+2.4);
 258     if(!tt.Fit(kTRUE)) continue;
 259
 260     x0     = tt.GetX0();
 261     xmean  = x0;//tt.GetXref();
 262     // printf("xmean[%f]\n", xmean);
 263
 264     dydx   = tt.GetYref(1);
 265     yfit   = -tt.GetYfit(0)*xmean+tt.GetYfit(0);
 266     yref   = -tt.GetYref(0)*xmean+tt.GetYref(0);
 267
 268     dzdx   = tt.GetZref(1);
 269     zfit   = -tt.GetZfit(0)*xmean+tt.GetZfit(0);
 270     zref   = -tt.GetZref(0)*xmean+tt.GetZref(0);
 271
 272     h->Fill(dydx, yref-yfit);
 273
 274
 275 /*    AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 276     fClResiduals->Add(clInfo);
 277     clInfo->SetGlobalPosition(yref, zref, phiref, thetaref);
 278     clInfo->SetResolution(yref-yfit);
 279
 280     (*fDebugStream) << "TrackletResiduals"
 281                     <<"clInfo.=" << clInfo
 282                     << "\n";*/
 283     if(fDebugLevel>=1){
 284       (*fDebugStream) << "TrkltResiduals"
 285         << "yref="     << yref
 286         << "yfit="     << yfit
 287         << "dydx="     << dydx
 288         << "zref="     << zref
 289         << "zfit="     << zfit
 290         << "dzdx="     << dzdx
 291         << "xmean="    << xmean
 292         << "x0="       << x0
 293         << "\n";
 294     }
 295   }
 296   return h;
 297 }
 298
 299 //________________________________________________________
 300 TH1* AliTRDtrackingResolution::PlotTrackletPhiResiduals(const AliTRDtrackV1 *track)
 301 {
 302   if(track) fTrack = track;
 303   if(!fTrack){
 304     AliWarning("No Track defined.");
 305     return 0x0;
 306   }
 307   TH1 *h = 0x0;
 308   if(!(h = ((TH2I*)fContainer->At(kTrackletPhiResidual)))){
 309     AliWarning("No output histogram defined.");
 310     return 0x0;
 311   }
 312
 313   Float_t  dydx_ref, dydx_fit;
 314   AliTRDseedV1 *fTracklet = 0x0;
 315   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 316     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 317     if(!fTracklet->IsOK()) continue;
 318
 319     dydx_ref   = fTracklet->GetYref(1);
 320     dydx_fit   = fTracklet->GetYfit(1);
 321
 322     h->Fill(dydx_ref, dydx_ref-dydx_fit);
 323
 324     if(fDebugLevel>=1){
 325       (*fDebugStream) << "TrkltPhiResiduals"
 326         << "dydx_ref="     << dydx_ref
 327         << "dydx_fit="     << dydx_fit
 328         << "\n";
 329     }
 330   }
 331   return h;
 332 }
 333
 334
 335 //________________________________________________________
 336 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
 337 {
 338   if(!HasMCdata()){
 339     AliWarning("No MC defined. Results will not be available.");
 340     return 0x0;
 341   }
 342   if(track) fTrack = track;
 343   if(!fTrack){
 344     AliWarning("No Track defined.");
 345     return 0x0;
 346   }
 347   TH1 *h = 0x0;
 348   if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
 349     AliWarning("No output histogram defined.");
 350     return 0x0;
 351   }
 352   if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
 353     AliWarning("No output histogram defined.");
 354     return 0x0;
 355   }
 356   if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
 357     AliWarning("No output histogram defined.");
 358     return 0x0;
 359   }
 360   if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
 361     AliWarning("No output histogram defined.");
 362     return 0x0;
 363   }
 364   //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
 365   UChar_t s;
 366   Int_t pdg = fMC->GetPDG(), det=-1;
 367   Int_t label = fMC->GetLabel();
 368   Float_t x0, y0, z0, dy, dydx, dzdx;
 369   AliTRDseedV1 *fTracklet = 0x0;
 370   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 371     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 372     if(!fTracklet->IsOK()) continue;
 373     //printf("process layer[%d]\n", ily);
 374     // retrive the track position and direction within the chamber
 375     det = fTracklet->GetDetector();
 376     x0  = fTracklet->GetX0();
 377     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 378
 379     // recalculate tracklet based on the MC info
 380     AliTRDseedV1 tt(*fTracklet);
 381     tt.SetZref(0, z0);
 382     tt.SetZref(1, dzdx);
 383     if(!tt.Fit(kFALSE)) continue;
 384     //tt.Update();
 385     dy = tt.GetYfit(0) - y0;
 386     Float_t dphi   = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
 387     Float_t dz = 100.;
 388     Bool_t  cross = fTracklet->GetNChange();
 389     if(cross){
 390       Double_t *xyz = tt.GetCrossXYZ();
 391       dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
 392       ((TH2I*)fContainer->At(kTrackletZResolution))->Fill(dzdx, dz);
 393     }
 394
 395     // Fill Histograms
 396     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
 397     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
 398
 399     // Fill Debug stream
 400     if(fDebugLevel>=1){
 401       Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
 402       (*fDebugStream) << "TrkltResolution"
 403         << "det="                 << det
 404         << "pdg="     << pdg
 405         << "p="       << p
 406         << "ymc="     << y0
 407         << "zmc="     << z0
 408         << "dydx="    << dydx
 409         << "dzdx="    << dzdx
 410         << "cross="   << cross
 411         << "dy="                  << dy
 412         << "dz="                  << dz
 413         << "dphi="              << dphi
 414         << "\n";
 415     }
 416
 417     Int_t istk = AliTRDgeometry::GetStack(det);
 418     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 419     Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
 420     Float_t tilt = fTracklet->GetTilt();
 421
 422     AliTRDcluster *c = 0x0;
 423     fTracklet->ResetClusterIter(kFALSE);
 424     while((c = fTracklet->PrevCluster())){
 425       Float_t  q = TMath::Abs(c->GetQ());
 426       Float_t xc = c->GetX();
 427       Float_t yc = c->GetY();
 428       Float_t zc = c->GetZ();
 429       Float_t dx = x0 - xc;
 430       Float_t yt = y0 - dx*dydx;
 431       Float_t zt = z0 - dx*dzdx;
 432       dy = yt - (yc - tilt*(zc-zt));
 433
 434       // Fill Histograms
 435       if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
 436
 437       // Fill Debug Tree
 438       if(fDebugLevel>=1){
 439         Float_t d = zr0 - zt;
 440         d -= ((Int_t)(2 * d)) / 2.0;
 441         if (d > 0.25) d  = 0.5 - d;
 442
 443         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 444         fClResolution->Add(clInfo);
 445         clInfo->SetCluster(c);
 446         clInfo->SetMC(pdg, label);
 447         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 448         clInfo->SetResolution(dy);
 449         clInfo->SetAnisochronity(d);
 450         clInfo->SetDriftLength(x0-xc);
 451         //clInfo->Print();
 452         (*fDebugStream) << "ClusterResolution"
 453           <<"clInfo.=" << clInfo
 454           << "\n";
 455       }
 456     }
 457   }
 458   return h;
 459 }
 460
 461
 462 //________________________________________________________
 463 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
 464 {
 465   TBox *b = 0x0;
 466   TAxis *ax = 0x0;
 467   TGraphErrors *g = 0x0;
 468   switch(ifig){
 469   case kClusterResidual:
 470     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 471     g->Draw("apl");
 472     ax = g->GetHistogram()->GetYaxis();
 473     ax->SetRangeUser(-.5, 1.);
 474     ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
 475     ax = g->GetHistogram()->GetXaxis();
 476     ax->SetTitle("tg(#phi)");
 477     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 478     g->Draw("pl");
 479     b = new TBox(-.15, -.5, .15, 1.);
 480     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 481     b->SetLineColor(0); b->Draw();
 482     return;
 483   case kTrackletYResidual:
 484     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 485     g->Draw("apl");
 486     ax = g->GetHistogram()->GetYaxis();
 487     ax->SetRangeUser(-.5, 1.);
 488     ax->SetTitle("Tracklet Y Residuals #sigma/#mu [mm]");
 489     ax = g->GetHistogram()->GetXaxis();
 490     ax->SetTitle("tg(#phi)");
 491     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 492     g->Draw("pl");
 493     b = new TBox(-.15, -.5, .15, 1.);
 494     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 495     b->SetLineColor(0); b->Draw();
 496     return;
 497   case kTrackletPhiResidual:
 498     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 499     g->Draw("apl");
 500     ax = g->GetHistogram()->GetYaxis();
 501     ax->SetRangeUser(-.5, 1.);
 502     ax->SetTitle("Tracklet Phi Residuals #sigma/#mu [mm]");
 503     ax = g->GetHistogram()->GetXaxis();
 504     ax->SetTitle("tg(#phi)");
 505     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 506     g->Draw("pl");
 507     b = new TBox(-.15, -.5, .15, 1.);
 508     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 509     b->SetLineColor(0); b->Draw();
 510     return;
 511   case kClusterResolution:
 512     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 513     ax = g->GetHistogram()->GetYaxis();
 514     ax->SetRangeUser(-.5, 1.);
 515     ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
 516     ax = g->GetHistogram()->GetXaxis();
 517     ax->SetTitle("tg(#phi)");
 518     g->Draw("apl");
 519     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 520     g->Draw("pl");
 521     b = new TBox(-.15, -.5, .15, 1.);
 522     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 523     b->SetLineColor(0); b->Draw();
 524     return;
 525   case kTrackletYResolution:
 526     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 527     ax = g->GetHistogram()->GetYaxis();
 528     ax->SetRangeUser(-.5, 1.);
 529     ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
 530     ax = g->GetHistogram()->GetXaxis();
 531     ax->SetTitle("tg(#phi)");
 532     g->Draw("apl");
 533     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 534     g->Draw("pl");
 535     return;
 536   case kTrackletZResolution:
 537     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 538     ax = g->GetHistogram()->GetYaxis();
 539     ax->SetRangeUser(-.5, 1.);
 540     ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
 541     ax = g->GetHistogram()->GetXaxis();
 542     ax->SetTitle("tg(#theta)");
 543     g->Draw("apl");
 544     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 545     g->Draw("pl");
 546     return;
 547   case kTrackletAngleResolution:
 548     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 549     ax = g->GetHistogram()->GetYaxis();
 550     ax->SetRangeUser(-.05, .2);
 551     ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
 552     ax = g->GetHistogram()->GetXaxis();
 553     ax->SetTitle("tg(#phi)");
 554     g->Draw("apl");
 555     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 556     g->Draw("pl");
 557     return;
 558   default:
 559     AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
 560     return;
 561   }
 562   AliInfo(Form("Reference plot [%d] missing result", ifig));
 563 }
 564
 565
 566 //________________________________________________________
 567 Bool_t AliTRDtrackingResolution::PostProcess()
 568 {
 569   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 570   if (!fContainer) {
 571     Printf("ERROR: list not available");
 572     return kFALSE;
 573   }
 574   fNRefFigures = fContainer->GetEntriesFast();
 575   if(!fGraphS){
 576     fGraphS = new TObjArray(fNRefFigures);
 577     fGraphS->SetOwner();
 578   }
 579   if(!fGraphM){
 580     fGraphM = new TObjArray(fNRefFigures);
 581     fGraphM->SetOwner();
 582   }
 583
 584   TH2I *h2 = 0x0;
 585   TH1D *h = 0x0;
 586   TGraphErrors *gm = 0x0, *gs = 0x0;
 587
 588   // define models
 589   TF1 f("f1", "gaus", -.5, .5);
 590
 591   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 592
 593   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 594
 595   TCanvas *c = 0x0;
 596   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
 597   char opt[5];
 598   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
 599
 600
 601   //PROCESS RESIDUAL DISTRIBUTIONS
 602
 603   // Clusters residuals
 604   h2 = (TH2I *)(fContainer->At(kClusterResidual));
 605   gm = new TGraphErrors();
 606   gm->SetLineColor(kBlue);
 607   gm->SetMarkerStyle(7);
 608   gm->SetMarkerColor(kBlue);
 609   gm->SetNameTitle("clm", "");
 610   fGraphM->AddAt(gm, kClusterResidual);
 611   gs = new TGraphErrors();
 612   gs->SetLineColor(kRed);
 613   gs->SetMarkerStyle(23);
 614   gs->SetMarkerColor(kRed);
 615   gs->SetNameTitle("cls", "");
 616   fGraphS->AddAt(gs, kClusterResidual);
 617   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 618     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 619     h = h2->ProjectionY("py", ibin, ibin);
 620     if(h->GetEntries()<100) continue;
 621     AdjustF1(h, &f);
 622
 623     if(IsVisual()){c->cd(); c->SetLogy();}
 624     h->Fit(&f, opt, "", -0.5, 0.5);
 625     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 626
 627     Int_t ip = gm->GetN();
 628     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 629     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 630     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 631     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 632   }
 633
 634   // Tracklet y residuals
 635   h2 = (TH2I *)(fContainer->At(kTrackletYResidual));
 636   gm = new TGraphErrors();
 637   gm->SetLineColor(kBlue);
 638   gm->SetMarkerStyle(7);
 639   gm->SetMarkerColor(kBlue);
 640   gm->SetNameTitle("tktm", "");
 641   fGraphM->AddAt(gm, kTrackletYResidual);
 642   gs = new TGraphErrors();
 643   gs->SetLineColor(kRed);
 644   gs->SetMarkerStyle(23);
 645   gs->SetMarkerColor(kRed);
 646   gs->SetNameTitle("tkts", "");
 647   fGraphS->AddAt(gs, kTrackletYResidual);
 648   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 649     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 650     h = h2->ProjectionY("py", ibin, ibin);
 651     if(h->GetEntries()<100) continue;
 652     AdjustF1(h, &f);
 653
 654     if(IsVisual()){c->cd(); c->SetLogy();}
 655     h->Fit(&f, opt, "", -0.5, 0.5);
 656     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 657
 658     Int_t ip = gm->GetN();
 659     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 660     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 661     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 662     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 663   }
 664
 665   // Tracklet phi residuals
 666   h2 = (TH2I *)(fContainer->At(kTrackletPhiResidual));
 667   gm = new TGraphErrors();
 668   gm->SetLineColor(kBlue);
 669   gm->SetMarkerStyle(7);
 670   gm->SetMarkerColor(kBlue);
 671   gm->SetNameTitle("tktphim", "");
 672   fGraphM->AddAt(gm, kTrackletPhiResidual);
 673   gs = new TGraphErrors();
 674   gs->SetLineColor(kRed);
 675   gs->SetMarkerStyle(23);
 676   gs->SetMarkerColor(kRed);
 677   gs->SetNameTitle("tktphis", "");
 678   fGraphS->AddAt(gs, kTrackletPhiResidual);
 679   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 680     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 681     h = h2->ProjectionY("py", ibin, ibin);
 682     if(h->GetEntries()<100) continue;
 683     AdjustF1(h, &f);
 684
 685     if(IsVisual()){c->cd(); c->SetLogy();}
 686     h->Fit(&f, opt, "", -0.5, 0.5);
 687     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 688
 689     Int_t ip = gm->GetN();
 690     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 691     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 692     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 693     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 694   }
 695
 696
 697   //PROCESS RESOLUTION DISTRIBUTIONS
 698
 699   if(HasMCdata()){
 700     // cluster y resolution
 701     h2 = (TH2I*)fContainer->At(kClusterResolution);
 702     gm = new TGraphErrors();
 703     gm->SetLineColor(kBlue);
 704     gm->SetMarkerStyle(7);
 705     gm->SetMarkerColor(kBlue);
 706     gm->SetNameTitle("clym", "");
 707     fGraphM->AddAt(gm, kClusterResolution);
 708     gs = new TGraphErrors();
 709     gs->SetLineColor(kRed);
 710     gs->SetMarkerStyle(23);
 711     gs->SetMarkerColor(kRed);
 712     gs->SetNameTitle("clys", "");
 713     fGraphS->AddAt(gs, kClusterResolution);
 714     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 715       h = h2->ProjectionY("py", iphi, iphi);
 716       if(h->GetEntries()<100) continue;
 717       AdjustF1(h, &f);
 718
 719       if(IsVisual()){c->cd(); c->SetLogy();}
 720       h->Fit(&f, opt, "", -0.5, 0.5);
 721       if(IsVerbose()){
 722         printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
 723       }
 724       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 725
 726       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 727       Int_t ip = gm->GetN();
 728       gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 729       gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 730       gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 731       gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 732     }
 733
 734     // tracklet y resolution
 735     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 736     gm = new TGraphErrors();
 737     gm->SetLineColor(kBlue);
 738     gm->SetMarkerStyle(7);
 739     gm->SetMarkerColor(kBlue);
 740     gm->SetNameTitle("trkltym", "");
 741     fGraphM->AddAt(gm, kTrackletYResolution);
 742     gs = new TGraphErrors();
 743     gs->SetLineColor(kRed);
 744     gs->SetMarkerStyle(23);
 745     gs->SetMarkerColor(kRed);
 746     gs->SetNameTitle("trkltys", "");
 747     fGraphS->AddAt(gs, kTrackletYResolution);
 748     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 749       h = h2->ProjectionY("py", iphi, iphi);
 750       if(h->GetEntries()<100) continue;
 751       AdjustF1(h, &fb);
 752
 753       if(IsVisual()){c->cd(); c->SetLogy();}
 754       h->Fit(&fb, opt, "", -0.5, 0.5);
 755       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 756
 757       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 758       Int_t ip = gm->GetN();
 759       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 760       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 761       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 762       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 763     }
 764
 765     // tracklet z resolution
 766     h2 = (TH2I*)fContainer->At(kTrackletZResolution);
 767     gm = new TGraphErrors();
 768     gm->SetLineColor(kBlue);
 769     gm->SetMarkerStyle(7);
 770     gm->SetMarkerColor(kBlue);
 771     gm->SetNameTitle("trkltzm", "");
 772     fGraphM->AddAt(gm, kTrackletZResolution);
 773     gs = new TGraphErrors();
 774     gs->SetLineColor(kRed);
 775     gs->SetMarkerStyle(23);
 776     gs->SetMarkerColor(kRed);
 777     gs->SetNameTitle("trkltzs", "");
 778     fGraphS->AddAt(gs, kTrackletZResolution);
 779     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 780       h = h2->ProjectionY("py", iphi, iphi);
 781       if(h->GetEntries()<100) continue;
 782       AdjustF1(h, &fb);
 783
 784       if(IsVisual()){c->cd(); c->SetLogy();}
 785       h->Fit(&fb, opt, "", -0.5, 0.5);
 786       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 787
 788       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 789       Int_t ip = gm->GetN();
 790       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 791       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 792       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 793       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 794     }
 795
 796     // tracklet phi resolution
 797     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 798     gm = new TGraphErrors();
 799     gm->SetLineColor(kBlue);
 800     gm->SetMarkerStyle(7);
 801     gm->SetMarkerColor(kBlue);
 802     gm->SetNameTitle("trkltym", "");
 803     fGraphM->AddAt(gm, kTrackletAngleResolution);
 804     gs = new TGraphErrors();
 805     gs->SetLineColor(kRed);
 806     gs->SetMarkerStyle(23);
 807     gs->SetMarkerColor(kRed);
 808     gs->SetNameTitle("trkltys", "");
 809     fGraphS->AddAt(gs, kTrackletAngleResolution);
 810     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 811       h = h2->ProjectionY("py", iphi, iphi);
 812       if(h->GetEntries()<100) continue;
 813
 814       if(IsVisual()){c->cd(); c->SetLogy();}
 815       h->Fit(&f, opt, "", -0.5, 0.5);
 816       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 817
 818       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 819       Int_t ip = gm->GetN();
 820       gm->SetPoint(ip, phi, f.GetParameter(1));
 821       gm->SetPointError(ip, 0., f.GetParError(1));
 822       gs->SetPoint(ip, phi, f.GetParameter(2));
 823       gs->SetPointError(ip, 0., f.GetParError(2));
 824     }
 825   }
 826   if(c) delete c;
 827
 828   return kTRUE;
 829 }
 830
 831
 832 //________________________________________________________
 833 void AliTRDtrackingResolution::Terminate(Option_t *)
 834 {
 835   if(fDebugStream){
 836     delete fDebugStream;
 837     fDebugStream = 0x0;
 838     fDebugLevel = 0;
 839   }
 840   if(HasPostProcess()) PostProcess();
 841 }
 842
 843 //________________________________________________________
 844 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
 845 {
 846 // Helper function to avoid duplication of code
 847 // Make first guesses on the fit parameters
 848
 849   // find the intial parameters of the fit !! (thanks George)
 850   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 851   Double_t sum = 0.;
 852   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 853   f->SetParLimits(0, 0., 3.*sum);
 854   f->SetParameter(0, .9*sum);
 855
 856   f->SetParLimits(1, -.2, .2);
 857   f->SetParameter(1, -0.1);
 858
 859   f->SetParLimits(2, 0., 4.e-1);
 860   f->SetParameter(2, 2.e-2);
 861   if(f->GetNpar() <= 4) return;
 862
 863   f->SetParLimits(3, 0., sum);
 864   f->SetParameter(3, .1*sum);
 865
 866   f->SetParLimits(4, -.3, .3);
 867   f->SetParameter(4, 0.);
 868
 869   f->SetParLimits(5, 0., 1.e2);
 870   f->SetParameter(5, 2.e-1);
 871 }
 872
 873 //________________________________________________________
 874 TObjArray* AliTRDtrackingResolution::Histos()
 875 {
 876   if(fContainer) return fContainer;
 877
 878   fContainer  = new TObjArray(5);
 879
 880   TH1 *h = 0x0;
 881   // cluster to tracklet residuals [2]
 882   fContainer->AddAt(h = new TH2I("fYCl", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5), kClusterResidual);
 883   h->GetXaxis()->SetTitle("tg(#phi)");
 884   h->GetYaxis()->SetTitle("#Delta y [cm]");
 885   h->GetZaxis()->SetTitle("entries");
 886   // tracklet to track residuals [2]
 887   fContainer->AddAt(h = new TH2I("hTrkltYRez", "Tracklets", 21, -.33, .33, 100, -.5, .5), kTrackletYResidual);
 888   h->GetXaxis()->SetTitle("tg(#phi)");
 889   h->GetYaxis()->SetTitle("#Delta y [cm]");
 890   h->GetZaxis()->SetTitle("entries");
 891   // tracklet to track residuals angular [2]
 892   fContainer->AddAt(h = new TH2I("hTrkltPhiRez", "Tracklets", 21, -.33, .33, 100, -.5, .5), kTrackletPhiResidual);
 893   h->GetXaxis()->SetTitle("tg(#phi)");
 894   h->GetYaxis()->SetTitle("#Delta phi [#circ]");
 895   h->GetZaxis()->SetTitle("entries");
 896
 897
 898   // Resolution histos
 899   if(HasMCdata()){
 900     // cluster y resolution [0]
 901     fContainer->AddAt(h = new TH2I("fYClMC", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5), kClusterResolution);
 902     h->GetXaxis()->SetTitle("tg(#phi)");
 903     h->GetYaxis()->SetTitle("#Delta y [cm]");
 904     h->GetZaxis()->SetTitle("entries");
 905     // tracklet y resolution [0]
 906     fContainer->AddAt(h = new TH2I("fYTrkltMC", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5), kTrackletYResolution);
 907     h->GetXaxis()->SetTitle("tg(#phi)");
 908     h->GetYaxis()->SetTitle("#Delta y [cm]");
 909     h->GetZaxis()->SetTitle("entries");
 910     // tracklet y resolution [0]
 911     fContainer->AddAt(h = new TH2I("fZTrkltMC", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5), kTrackletZResolution);
 912     h->GetXaxis()->SetTitle("tg(#theta)");
 913     h->GetYaxis()->SetTitle("#Delta z [cm]");
 914     h->GetZaxis()->SetTitle("entries");
 915     // tracklet angular resolution [1]
 916     fContainer->AddAt(h = new TH2I("fPhiTrkltMC", "Tracklet Resolution (Angular)", 31, -.48, .48, 100, -10., 10.), kTrackletAngleResolution);
 917     h->GetXaxis()->SetTitle("tg(#phi)");
 918     h->GetYaxis()->SetTitle("#Delta #phi [deg]");
 919     h->GetZaxis()->SetTitle("entries");
 920
 921 //     // Riemann track resolution [y, z, angular]
 922 //     fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
 923 //     fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
 924 //     fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
 925 //
 926 //     Kalman track resolution [y, z, angular]
 927 //     fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
 928 //     fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
 929 //     fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
 930   }
 931   return fContainer;
 932 }
 933
 934
 935 //________________________________________________________
 936 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 937 {
 938
 939   fReconstructor->SetRecoParam(r);
 940 }