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   AliTRDseedV1 *tracklet = 0x0;
 252   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 253     if(!(tracklet = fTrack->GetTracklet(il))) continue;
 254     if(!tracklet->IsOK()) continue;
 255     h->Fill(tracklet->GetYref(1), tracklet->GetYref(0)-tracklet->GetYfit(0));
 256   }
 257   return h;
 258
 259   // refit the track
 260 //   AliRieman fRim(fTrack->GetNumberOfClusters());
 261 //   Float_t x[AliTRDgeometry::kNlayer] = {-1., -1., -1., -1., -1., -1.}, y[AliTRDgeometry::kNlayer], dydx[AliTRDgeometry::kNlayer];
 262
 263 //     AliTRDcluster *c = 0x0;
 264 //     tracklet->ResetClusterIter(kFALSE);
 265 //     while((c = tracklet->PrevCluster())){
 266 //       Float_t xc = c->GetX();
 267 //       Float_t yc = c->GetY();
 268 //       Float_t zc = c->GetZ();
 269 //       Float_t zt = tracklet->GetZref(0) - (tracklet->GetX0()-xc)*tracklet->GetZref(1);
 270 //       yc -= tracklet->GetTilt()*(zc-zt);
 271 //       fRim.AddPoint(xc, yc, zc, 1, 10);
 272 //     }
 273 //     tracklet->Fit(kTRUE);
 274 //
 275 //     x[il] = tracklet->GetX0();
 276 //     y[il] = tracklet->GetYfit(0)-tracklet->GetYfit(1)*(tracklet->GetX0()-x[il]);
 277 //     dydx[il] = tracklet->GetYref(1);
 278
 279
 280 //   fRim.Update();
 281
 282 //   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 283 //     if(x[il] < 0.) continue;
 284 //     Float_t dy = y[il]-fRim.GetYat(x[il])/*/sigma_track*/;
 285 //     h->Fill(dydx[il], dy);
 286 //
 287 //     if(fDebugLevel>=1){
 288 //       Double_t sigmay = fRim.GetErrY(x[il]);
 289 //       Float_t yt = fRim.GetYat(x[il]);
 290 //       (*fDebugStream) << "TrkltResiduals"
 291 //         << "layer="  << il
 292 //         << "x="      << x[il]
 293 //         << "y="      << y[il]
 294 //         << "yt="     << yt
 295 //         << "dydx="   << dydx[il]
 296 //         << "dy="     << dy
 297 //         << "sigmay=" << sigmay
 298 //         << "\n";
 299 //     }
 300 //   }
 301 }
 302
 303 //________________________________________________________
 304 TH1* AliTRDtrackingResolution::PlotTrackletPhiResiduals(const AliTRDtrackV1 *track)
 305 {
 306   if(track) fTrack = track;
 307   if(!fTrack){
 308     AliWarning("No Track defined.");
 309     return 0x0;
 310   }
 311   TH1 *h = 0x0;
 312   if(!(h = ((TH2I*)fContainer->At(kTrackletPhiResidual)))){
 313     AliWarning("No output histogram defined.");
 314     return 0x0;
 315   }
 316
 317   AliTRDseedV1 *tracklet = 0x0;
 318   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 319     if(!(tracklet = fTrack->GetTracklet(il))) continue;
 320     if(!tracklet->IsOK()) continue;
 321     h->Fill(tracklet->GetYref(1), TMath::ATan(tracklet->GetYref(1))-TMath::ATan(tracklet->GetYfit(1)));
 322   }
 323   return h;
 324
 325 //   // refit the track
 326 //   AliRieman fRim(fTrack->GetNumberOfClusters());
 327 //   Float_t x[AliTRDgeometry::kNlayer] = {-1., -1., -1., -1., -1., -1.}, y[AliTRDgeometry::kNlayer], dydx[AliTRDgeometry::kNlayer];
 328 //   Float_t dydx_ref=0, dydx_fit=0, phiref=0, phifit=0, phidiff=0;
 329 //   AliTRDseedV1 *tracklet = 0x0;
 330 //   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 331 //     if(!(tracklet = fTrack->GetTracklet(il))) continue;
 332 //     if(!tracklet->IsOK()) continue;
 333 //     AliTRDcluster *c = 0x0;
 334 //     tracklet->ResetClusterIter(kFALSE);
 335 //     while((c = tracklet->PrevCluster())){
 336 //       Float_t xc = c->GetX();
 337 //       Float_t yc = c->GetY();
 338 //       Float_t zc = c->GetZ();
 339 //       Float_t zt = tracklet->GetZref(0) - (tracklet->GetX0()-xc)*tracklet->GetZref(1);
 340 //       yc -= tracklet->GetTilt()*(zc-zt);
 341 //       fRim.AddPoint(xc, yc, zc, 1, 10);
 342 //     }
 343 //     tracklet->Fit(kTRUE);
 344 //
 345 //     x[il] = tracklet->GetX0();
 346 //     y[il] = tracklet->GetYfit(0)-tracklet->GetYfit(1)*(tracklet->GetX0()-x[il]);
 347 //     dydx[il] = tracklet->GetYref(1);
 348 //   }
 349 //   fRim.Update();
 350 //
 351 //   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 352 //     if(x[il] < 0.) continue;
 353 //     if(!(tracklet = fTrack->GetTracklet(il))) continue;
 354 //     if(!tracklet->IsOK()) continue;
 355 //
 356 //     dydx_ref = fRim.GetDYat(x[il]);
 357 //     dydx_fit = tracklet->GetYfit(1);
 358 //
 359 //     phiref = TMath::ATan(dydx_ref);//*TMath::RadToDeg();
 360 //     phifit = TMath::ATan(dydx_fit);//*TMath::RadToDeg();
 361 //
 362 //     phidiff = phiref-phifit; /*/sigma_phi*/;
 363 //
 364 //     h->Fill(dydx_ref, phidiff);
 365 //
 366 //
 367 //     if(fDebugLevel>=1){
 368 //       (*fDebugStream) << "TrkltPhiResiduals"
 369 //         << "layer="  << il
 370 //         << "dydx_fit="   << dydx_fit
 371 //         << "dydx_ref="   << dydx_ref
 372 //         << "phiref="     << phiref
 373 //         << "phifit="     << phifit
 374 //         << "phidiff="     << phidiff
 375 //         << "\n";
 376 //     }
 377 //   }
 378 }
 379
 380
 381 //________________________________________________________
 382 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
 383 {
 384   if(!HasMCdata()){
 385     AliWarning("No MC defined. Results will not be available.");
 386     return 0x0;
 387   }
 388   if(track) fTrack = track;
 389   if(!fTrack){
 390     AliWarning("No Track defined.");
 391     return 0x0;
 392   }
 393   TH1 *h = 0x0;
 394   if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
 395     AliWarning("No output histogram defined.");
 396     return 0x0;
 397   }
 398   if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
 399     AliWarning("No output histogram defined.");
 400     return 0x0;
 401   }
 402   if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
 403     AliWarning("No output histogram defined.");
 404     return 0x0;
 405   }
 406   if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
 407     AliWarning("No output histogram defined.");
 408     return 0x0;
 409   }
 410   //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
 411   UChar_t s;
 412   Int_t pdg = fMC->GetPDG(), det=-1;
 413   Int_t label = fMC->GetLabel();
 414   Float_t x0, y0, z0, dx, dy, dydx, dzdx;
 415   AliTRDseedV1 *fTracklet = 0x0;
 416   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 417     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 418     if(!fTracklet->IsOK()) continue;
 419     //printf("process layer[%d]\n", ily);
 420     // retrive the track position and direction within the chamber
 421     det = fTracklet->GetDetector();
 422     x0  = fTracklet->GetX0();
 423     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 424
 425     // recalculate tracklet based on the MC info
 426     AliTRDseedV1 tt(*fTracklet);
 427     tt.SetZref(0, z0);
 428     tt.SetZref(1, dzdx);
 429     if(!tt.Fit(kTRUE)) continue;
 430     //tt.Update();
 431
 432     dx = 0.;//x0 - tt.GetXref();
 433     Float_t yt = y0 - dx*dydx;
 434     Float_t yf = tt.GetYfit(0) - dx*tt.GetYfit(1);
 435     dy = yf-yt;
 436     Float_t dphi   = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
 437     Float_t dz = 100.;
 438     Bool_t  cross = fTracklet->GetNChange();
 439     if(cross){
 440       Double_t *xyz = tt.GetCrossXYZ();
 441       dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
 442       ((TH2I*)fContainer->At(kTrackletZResolution))->Fill(dzdx, dz);
 443     }
 444
 445     // Fill Histograms
 446     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
 447     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
 448
 449     // Fill Debug stream
 450     if(fDebugLevel>=1){
 451       Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
 452       (*fDebugStream) << "TrkltResolution"
 453         << "det="                 << det
 454         << "pdg="     << pdg
 455         << "p="       << p
 456         << "ymc="     << y0
 457         << "zmc="     << z0
 458         << "dydx="    << dydx
 459         << "dzdx="    << dzdx
 460         << "cross="   << cross
 461         << "dy="                  << dy
 462         << "dz="                  << dz
 463         << "dphi="              << dphi
 464         << "\n";
 465     }
 466
 467     Int_t istk = AliTRDgeometry::GetStack(det);
 468     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 469     Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
 470     Float_t tilt = fTracklet->GetTilt();
 471
 472     AliTRDcluster *c = 0x0;
 473     fTracklet->ResetClusterIter(kFALSE);
 474     while((c = fTracklet->PrevCluster())){
 475       Float_t  q = TMath::Abs(c->GetQ());
 476       Float_t xc = c->GetX();
 477       Float_t yc = c->GetY();
 478       Float_t zc = c->GetZ();
 479       dx = x0 - xc;
 480       Float_t yt = y0 - dx*dydx;
 481       Float_t zt = z0 - dx*dzdx;
 482       dy = yt - (yc - tilt*(zc-zt));
 483
 484       // Fill Histograms
 485       if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
 486
 487       // Fill Debug Tree
 488       if(fDebugLevel>=1){
 489         Float_t d = zr0 - zt;
 490         d -= ((Int_t)(2 * d)) / 2.0;
 491         if (d > 0.25) d  = 0.5 - d;
 492
 493         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 494         fClResolution->Add(clInfo);
 495         clInfo->SetCluster(c);
 496         clInfo->SetMC(pdg, label);
 497         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 498         clInfo->SetResolution(dy);
 499         clInfo->SetAnisochronity(d);
 500         clInfo->SetDriftLength(x0-xc);
 501         //clInfo->Print();
 502         (*fDebugStream) << "ClusterResolution"
 503           <<"clInfo.=" << clInfo
 504           << "\n";
 505       }
 506     }
 507   }
 508   return h;
 509 }
 510
 511
 512 //________________________________________________________
 513 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
 514 {
 515   TBox *b = 0x0;
 516   TAxis *ax = 0x0;
 517   TGraphErrors *g = 0x0;
 518   switch(ifig){
 519   case kClusterResidual:
 520     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 521     g->Draw("apl");
 522     ax = g->GetHistogram()->GetYaxis();
 523     ax->SetRangeUser(-.5, 2.5);
 524     ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
 525     ax = g->GetHistogram()->GetXaxis();
 526     ax->SetTitle("tg(#phi)");
 527     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 528     g->Draw("pl");
 529     b = new TBox(-.15, -.5, .15, 1.);
 530     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 531     b->SetLineColor(0); b->Draw();
 532     return;
 533   case kTrackletYResidual:
 534     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 535     g->Draw("apl");
 536     ax = g->GetHistogram()->GetYaxis();
 537     ax->SetRangeUser(-.5, 3.);
 538     ax->SetTitle("Tracklet Y Residuals #sigma/#mu [mm]");
 539     ax = g->GetHistogram()->GetXaxis();
 540     ax->SetTitle("tg(#phi)");
 541     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 542     g->Draw("pl");
 543     b = new TBox(-.15, -.5, .15, 1.);
 544     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 545     b->SetLineColor(0); b->Draw();
 546     return;
 547   case kTrackletPhiResidual:
 548     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 549     g->Draw("apl");
 550     ax = g->GetHistogram()->GetYaxis();
 551     ax->SetRangeUser(-.5, 2.);
 552     ax->SetTitle("Tracklet Phi Residuals #sigma/#mu [rad]");
 553     ax = g->GetHistogram()->GetXaxis();
 554     ax->SetTitle("tg(#phi)");
 555     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 556     g->Draw("pl");
 557     b = new TBox(-.15, -.5, .15, 1.);
 558     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 559     b->SetLineColor(0); b->Draw();
 560     return;
 561   case kClusterResolution:
 562     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 563     ax = g->GetHistogram()->GetYaxis();
 564     ax->SetRangeUser(-.5, 1.);
 565     ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
 566     ax = g->GetHistogram()->GetXaxis();
 567     ax->SetTitle("tg(#phi)");
 568     g->Draw("apl");
 569     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 570     g->Draw("pl");
 571     b = new TBox(-.15, -.5, .15, 1.);
 572     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 573     b->SetLineColor(0); b->Draw();
 574     return;
 575   case kTrackletYResolution:
 576     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 577     ax = g->GetHistogram()->GetYaxis();
 578     ax->SetRangeUser(-.5, 1.);
 579     ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
 580     ax = g->GetHistogram()->GetXaxis();
 581     ax->SetTitle("tg(#phi)");
 582     g->Draw("apl");
 583     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 584     g->Draw("pl");
 585     return;
 586   case kTrackletZResolution:
 587     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 588     ax = g->GetHistogram()->GetYaxis();
 589     ax->SetRangeUser(-.5, 1.);
 590     ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
 591     ax = g->GetHistogram()->GetXaxis();
 592     ax->SetTitle("tg(#theta)");
 593     g->Draw("apl");
 594     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 595     g->Draw("pl");
 596     return;
 597   case kTrackletAngleResolution:
 598     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 599     ax = g->GetHistogram()->GetYaxis();
 600     ax->SetRangeUser(-.05, .2);
 601     ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
 602     ax = g->GetHistogram()->GetXaxis();
 603     ax->SetTitle("tg(#phi)");
 604     g->Draw("apl");
 605     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 606     g->Draw("pl");
 607     return;
 608   default:
 609     AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
 610     return;
 611   }
 612   AliInfo(Form("Reference plot [%d] missing result", ifig));
 613 }
 614
 615
 616 //________________________________________________________
 617 Bool_t AliTRDtrackingResolution::PostProcess()
 618 {
 619   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 620   if (!fContainer) {
 621     Printf("ERROR: list not available");
 622     return kFALSE;
 623   }
 624   fNRefFigures = fContainer->GetEntriesFast();
 625   if(!fGraphS){
 626     fGraphS = new TObjArray(fNRefFigures);
 627     fGraphS->SetOwner();
 628   }
 629   if(!fGraphM){
 630     fGraphM = new TObjArray(fNRefFigures);
 631     fGraphM->SetOwner();
 632   }
 633
 634   TH2I *h2 = 0x0;
 635   TH1D *h = 0x0;
 636   TGraphErrors *gm = 0x0, *gs = 0x0;
 637
 638   // define models
 639   TF1 f("f1", "gaus", -.5, .5);
 640
 641   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 642
 643   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 644
 645   TCanvas *c = 0x0;
 646   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
 647   char opt[5];
 648   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
 649
 650
 651   //PROCESS RESIDUAL DISTRIBUTIONS
 652
 653   // Clusters residuals
 654   h2 = (TH2I *)(fContainer->At(kClusterResidual));
 655   gm = new TGraphErrors();
 656   gm->SetLineColor(kBlue);
 657   gm->SetMarkerStyle(7);
 658   gm->SetMarkerColor(kBlue);
 659   gm->SetNameTitle("clm", "");
 660   fGraphM->AddAt(gm, kClusterResidual);
 661   gs = new TGraphErrors();
 662   gs->SetLineColor(kRed);
 663   gs->SetMarkerStyle(23);
 664   gs->SetMarkerColor(kRed);
 665   gs->SetNameTitle("cls", "");
 666   fGraphS->AddAt(gs, kClusterResidual);
 667   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 668     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 669     h = h2->ProjectionY("py", ibin, ibin);
 670     if(h->GetEntries()<100) continue;
 671     AdjustF1(h, &f);
 672
 673     if(IsVisual()){c->cd(); c->SetLogy();}
 674     h->Fit(&f, opt, "", -0.5, 0.5);
 675     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 676
 677     Int_t ip = gm->GetN();
 678     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 679     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 680     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 681     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 682   }
 683
 684   // Tracklet y residuals
 685   h2 = (TH2I *)(fContainer->At(kTrackletYResidual));
 686   gm = new TGraphErrors();
 687   gm->SetLineColor(kBlue);
 688   gm->SetMarkerStyle(7);
 689   gm->SetMarkerColor(kBlue);
 690   gm->SetNameTitle("tktm", "");
 691   fGraphM->AddAt(gm, kTrackletYResidual);
 692   gs = new TGraphErrors();
 693   gs->SetLineColor(kRed);
 694   gs->SetMarkerStyle(23);
 695   gs->SetMarkerColor(kRed);
 696   gs->SetNameTitle("tkts", "");
 697   fGraphS->AddAt(gs, kTrackletYResidual);
 698   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 699     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 700     h = h2->ProjectionY("py", ibin, ibin);
 701     if(h->GetEntries()<100) continue;
 702     AdjustF1(h, &f);
 703
 704     if(IsVisual()){c->cd(); c->SetLogy();}
 705     h->Fit(&f, opt, "", -0.5, 0.5);
 706     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 707
 708     Int_t ip = gm->GetN();
 709     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 710     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 711     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 712     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 713   }
 714
 715   // Tracklet phi residuals
 716   h2 = (TH2I *)(fContainer->At(kTrackletPhiResidual));
 717   gm = new TGraphErrors();
 718   gm->SetLineColor(kBlue);
 719   gm->SetMarkerStyle(7);
 720   gm->SetMarkerColor(kBlue);
 721   gm->SetNameTitle("tktphim", "");
 722   fGraphM->AddAt(gm, kTrackletPhiResidual);
 723   gs = new TGraphErrors();
 724   gs->SetLineColor(kRed);
 725   gs->SetMarkerStyle(23);
 726   gs->SetMarkerColor(kRed);
 727   gs->SetNameTitle("tktphis", "");
 728   fGraphS->AddAt(gs, kTrackletPhiResidual);
 729   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 730     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 731     h = h2->ProjectionY("py", ibin, ibin);
 732     if(h->GetEntries()<100) continue;
 733     AdjustF1(h, &f);
 734
 735     if(IsVisual()){c->cd(); c->SetLogy();}
 736     h->Fit(&f, opt, "", -0.5, 0.5);
 737     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 738
 739     Int_t ip = gm->GetN();
 740     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 741     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 742     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 743     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 744   }
 745
 746
 747   //PROCESS RESOLUTION DISTRIBUTIONS
 748
 749   if(HasMCdata()){
 750     // cluster y resolution
 751     h2 = (TH2I*)fContainer->At(kClusterResolution);
 752     gm = new TGraphErrors();
 753     gm->SetLineColor(kBlue);
 754     gm->SetMarkerStyle(7);
 755     gm->SetMarkerColor(kBlue);
 756     gm->SetNameTitle("clym", "");
 757     fGraphM->AddAt(gm, kClusterResolution);
 758     gs = new TGraphErrors();
 759     gs->SetLineColor(kRed);
 760     gs->SetMarkerStyle(23);
 761     gs->SetMarkerColor(kRed);
 762     gs->SetNameTitle("clys", "");
 763     fGraphS->AddAt(gs, kClusterResolution);
 764     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 765       h = h2->ProjectionY("py", iphi, iphi);
 766       if(h->GetEntries()<100) continue;
 767       AdjustF1(h, &f);
 768
 769       if(IsVisual()){c->cd(); c->SetLogy();}
 770       h->Fit(&f, opt, "", -0.5, 0.5);
 771       if(IsVerbose()){
 772         printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
 773       }
 774       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 775
 776       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 777       Int_t ip = gm->GetN();
 778       gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
 779       gm->SetPointError(ip, 0., 10.*f.GetParError(1));
 780       gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
 781       gs->SetPointError(ip, 0., 10.*f.GetParError(2));
 782     }
 783
 784     // tracklet y resolution
 785     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 786     gm = new TGraphErrors();
 787     gm->SetLineColor(kBlue);
 788     gm->SetMarkerStyle(7);
 789     gm->SetMarkerColor(kBlue);
 790     gm->SetNameTitle("trkltym", "");
 791     fGraphM->AddAt(gm, kTrackletYResolution);
 792     gs = new TGraphErrors();
 793     gs->SetLineColor(kRed);
 794     gs->SetMarkerStyle(23);
 795     gs->SetMarkerColor(kRed);
 796     gs->SetNameTitle("trkltys", "");
 797     fGraphS->AddAt(gs, kTrackletYResolution);
 798     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 799       h = h2->ProjectionY("py", iphi, iphi);
 800       if(h->GetEntries()<100) continue;
 801       AdjustF1(h, &fb);
 802
 803       if(IsVisual()){c->cd(); c->SetLogy();}
 804       h->Fit(&fb, opt, "", -0.5, 0.5);
 805       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 806
 807       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 808       Int_t ip = gm->GetN();
 809       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 810       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 811       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 812       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 813     }
 814
 815     // tracklet z resolution
 816     h2 = (TH2I*)fContainer->At(kTrackletZResolution);
 817     gm = new TGraphErrors();
 818     gm->SetLineColor(kBlue);
 819     gm->SetMarkerStyle(7);
 820     gm->SetMarkerColor(kBlue);
 821     gm->SetNameTitle("trkltzm", "");
 822     fGraphM->AddAt(gm, kTrackletZResolution);
 823     gs = new TGraphErrors();
 824     gs->SetLineColor(kRed);
 825     gs->SetMarkerStyle(23);
 826     gs->SetMarkerColor(kRed);
 827     gs->SetNameTitle("trkltzs", "");
 828     fGraphS->AddAt(gs, kTrackletZResolution);
 829     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 830       h = h2->ProjectionY("py", iphi, iphi);
 831       if(h->GetEntries()<100) continue;
 832       AdjustF1(h, &fb);
 833
 834       if(IsVisual()){c->cd(); c->SetLogy();}
 835       h->Fit(&fb, opt, "", -0.5, 0.5);
 836       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 837
 838       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 839       Int_t ip = gm->GetN();
 840       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
 841       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
 842       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
 843       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
 844     }
 845
 846     // tracklet phi resolution
 847     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 848     gm = new TGraphErrors();
 849     gm->SetLineColor(kBlue);
 850     gm->SetMarkerStyle(7);
 851     gm->SetMarkerColor(kBlue);
 852     gm->SetNameTitle("trkltym", "");
 853     fGraphM->AddAt(gm, kTrackletAngleResolution);
 854     gs = new TGraphErrors();
 855     gs->SetLineColor(kRed);
 856     gs->SetMarkerStyle(23);
 857     gs->SetMarkerColor(kRed);
 858     gs->SetNameTitle("trkltys", "");
 859     fGraphS->AddAt(gs, kTrackletAngleResolution);
 860     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 861       h = h2->ProjectionY("py", iphi, iphi);
 862       if(h->GetEntries()<100) continue;
 863
 864       if(IsVisual()){c->cd(); c->SetLogy();}
 865       h->Fit(&f, opt, "", -0.5, 0.5);
 866       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 867
 868       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 869       Int_t ip = gm->GetN();
 870       gm->SetPoint(ip, phi, f.GetParameter(1));
 871       gm->SetPointError(ip, 0., f.GetParError(1));
 872       gs->SetPoint(ip, phi, f.GetParameter(2));
 873       gs->SetPointError(ip, 0., f.GetParError(2));
 874     }
 875   }
 876   if(c) delete c;
 877
 878   return kTRUE;
 879 }
 880
 881
 882 //________________________________________________________
 883 void AliTRDtrackingResolution::Terminate(Option_t *)
 884 {
 885   if(fDebugStream){
 886     delete fDebugStream;
 887     fDebugStream = 0x0;
 888     fDebugLevel = 0;
 889   }
 890   if(HasPostProcess()) PostProcess();
 891 }
 892
 893 //________________________________________________________
 894 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
 895 {
 896 // Helper function to avoid duplication of code
 897 // Make first guesses on the fit parameters
 898
 899   // find the intial parameters of the fit !! (thanks George)
 900   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 901   Double_t sum = 0.;
 902   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 903   f->SetParLimits(0, 0., 3.*sum);
 904   f->SetParameter(0, .9*sum);
 905
 906   f->SetParLimits(1, -.2, .2);
 907   f->SetParameter(1, -0.1);
 908
 909   f->SetParLimits(2, 0., 4.e-1);
 910   f->SetParameter(2, 2.e-2);
 911   if(f->GetNpar() <= 4) return;
 912
 913   f->SetParLimits(3, 0., sum);
 914   f->SetParameter(3, .1*sum);
 915
 916   f->SetParLimits(4, -.3, .3);
 917   f->SetParameter(4, 0.);
 918
 919   f->SetParLimits(5, 0., 1.e2);
 920   f->SetParameter(5, 2.e-1);
 921 }
 922
 923 //________________________________________________________
 924 TObjArray* AliTRDtrackingResolution::Histos()
 925 {
 926   if(fContainer) return fContainer;
 927
 928   fContainer  = new TObjArray(7);
 929
 930   TH1 *h = 0x0;
 931   // cluster to tracklet residuals [2]
 932   if(!(h = (TH2I*)gROOT->FindObject("fYCl"))){
 933     h = new TH2I("fYCl", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5);
 934     h->GetXaxis()->SetTitle("tg(#phi)");
 935     h->GetYaxis()->SetTitle("#Delta y [cm]");
 936     h->GetZaxis()->SetTitle("entries");
 937   } else h->Reset();
 938   fContainer->AddAt(h, kClusterResidual);
 939
 940   // tracklet to track residuals [2]
 941   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYRez"))){
 942     h = new TH2I("hTrkltYRez", "Tracklets", 21, -.33, .33, 100, -.5, .5);
 943     h->GetXaxis()->SetTitle("tg(#phi)");
 944     h->GetYaxis()->SetTitle("#Delta y [cm]");
 945     h->GetZaxis()->SetTitle("entries");
 946   } else h->Reset();
 947   fContainer->AddAt(h, kTrackletYResidual);
 948
 949   // tracklet to track residuals angular [2]
 950   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhiRez"))){
 951     h = new TH2I("hTrkltPhiRez", "Tracklets", 21, -.33, .33, 100, -.5, .5);
 952     h->GetXaxis()->SetTitle("tg(#phi)");
 953     h->GetYaxis()->SetTitle("#Delta phi [#circ]");
 954     h->GetZaxis()->SetTitle("entries");
 955   } else h->Reset();
 956   fContainer->AddAt(h, kTrackletPhiResidual);
 957
 958
 959   // Resolution histos
 960   if(HasMCdata()){
 961     // cluster y resolution [0]
 962     if(!(h = (TH2I*)gROOT->FindObject("fYClMC"))){
 963       h = new TH2I("fYClMC", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
 964       h->GetXaxis()->SetTitle("tg(#phi)");
 965       h->GetYaxis()->SetTitle("#Delta y [cm]");
 966       h->GetZaxis()->SetTitle("entries");
 967     } else h->Reset();
 968     fContainer->AddAt(h, kClusterResolution);
 969
 970     // tracklet y resolution [0]
 971     if(!(h = (TH2I*)gROOT->FindObject("fYTrkltMC"))){
 972       h = new TH2I("fYTrkltMC", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
 973       h->GetXaxis()->SetTitle("tg(#phi)");
 974       h->GetYaxis()->SetTitle("#Delta y [cm]");
 975       h->GetZaxis()->SetTitle("entries");
 976     } else h->Reset();
 977     fContainer->AddAt(h, kTrackletYResolution);
 978
 979     // tracklet y resolution [0]
 980     if(!(h = (TH2I*)gROOT->FindObject("fZTrkltMC"))){
 981       h = new TH2I("fZTrkltMC", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
 982       h->GetXaxis()->SetTitle("tg(#theta)");
 983       h->GetYaxis()->SetTitle("#Delta z [cm]");
 984       h->GetZaxis()->SetTitle("entries");
 985     } else h->Reset();
 986     fContainer->AddAt(h, kTrackletZResolution);
 987
 988     // tracklet angular resolution [1]
 989     if(!(h = (TH2I*)gROOT->FindObject("fPhiTrkltMC"))){
 990       h = new TH2I("fPhiTrkltMC", "Tracklet Resolution (Angular)", 31, -.48, .48, 100, -10., 10.);
 991       h->GetXaxis()->SetTitle("tg(#phi)");
 992       h->GetYaxis()->SetTitle("#Delta #phi [deg]");
 993       h->GetZaxis()->SetTitle("entries");
 994     } else h->Reset();
 995     fContainer->AddAt(h, kTrackletAngleResolution);
 996   }
 997   return fContainer;
 998 }
 999
1000
1001 //________________________________________________________
1002 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
1003 {
1004
1005   fReconstructor->SetRecoParam(r);
1006 }