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