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