TRD/qaRec/AliTRDresolution.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: AliTRDresolution.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 //   AliTRDresolution *res = new AliTRDresolution();
  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 <TPDGCode.h>
  56 #include <TObjArray.h>
  57 #include <TH2.h>
  58 #include <TH1.h>
  59 #include <TF1.h>
  60 #include <TCanvas.h>
  61 #include <TBox.h>
  62 #include <TProfile.h>
  63 #include <TGraphErrors.h>
  64 #include <TMath.h>
  65 #include "TTreeStream.h"
  66 #include "TGeoManager.h"
  67
  68 #include "AliAnalysisManager.h"
  69 #include "AliTrackReference.h"
  70 #include "AliTrackPointArray.h"
  71 #include "AliCDBManager.h"
  72
  73 #include "AliTRDSimParam.h"
  74 #include "AliTRDgeometry.h"
  75 #include "AliTRDpadPlane.h"
  76 #include "AliTRDcluster.h"
  77 #include "AliTRDseedV1.h"
  78 #include "AliTRDtrackV1.h"
  79 #include "AliTRDtrackerV1.h"
  80 #include "AliTRDReconstructor.h"
  81 #include "AliTRDrecoParam.h"
  82
  83 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
  84 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  85 #include "AliTRDresolution.h"
  86
  87 ClassImp(AliTRDresolution)
  88
  89 //________________________________________________________
  90 AliTRDresolution::AliTRDresolution()
  91   :AliTRDrecoTask("Resolution", "Spatial and Momentum Resolution")
  92   ,fStatus(0)
  93   ,fReconstructor(0x0)
  94   ,fGeo(0x0)
  95   ,fGraphS(0x0)
  96   ,fGraphM(0x0)
  97   ,fCl(0x0)
  98   ,fTrklt(0x0)
  99   ,fMCcl(0x0)
 100   ,fMCtrklt(0x0)
 101 {
 102   fReconstructor = new AliTRDReconstructor();
 103   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 104   fGeo = new AliTRDgeometry();
 105
 106   InitFunctorList();
 107
 108   DefineOutput(1, TObjArray::Class()); // cluster2track
 109   DefineOutput(2, TObjArray::Class()); // tracklet2track
 110   DefineOutput(3, TObjArray::Class()); // cluster2mc
 111   DefineOutput(4, TObjArray::Class()); // tracklet2mc
 112 }
 113
 114 //________________________________________________________
 115 AliTRDresolution::~AliTRDresolution()
 116 {
 117   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 118   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 119   delete fGeo;
 120   delete fReconstructor;
 121   if(gGeoManager) delete gGeoManager;
 122   if(fCl){fCl->Delete(); delete fCl;}
 123   if(fTrklt){fTrklt->Delete(); delete fTrklt;}
 124   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 125   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
 126 }
 127
 128
 129 //________________________________________________________
 130 void AliTRDresolution::CreateOutputObjects()
 131 {
 132   // spatial resolution
 133   OpenFile(0, "RECREATE");
 134
 135   fContainer = Histos();
 136
 137   fCl = new TObjArray();
 138   fCl->SetOwner(kTRUE);
 139   fTrklt = new TObjArray();
 140   fTrklt->SetOwner(kTRUE);
 141   fMCcl = new TObjArray();
 142   fMCcl->SetOwner(kTRUE);
 143   fMCtrklt = new TObjArray();
 144   fMCtrklt->SetOwner(kTRUE);
 145 }
 146
 147 //________________________________________________________
 148 void AliTRDresolution::Exec(Option_t *opt)
 149 {
 150   fCl->Delete();
 151   fTrklt->Delete();
 152   fMCcl->Delete();
 153   fMCtrklt->Delete();
 154
 155   AliTRDrecoTask::Exec(opt);
 156
 157   PostData(1, fCl);
 158   PostData(2, fTrklt);
 159   PostData(3, fMCcl);
 160   PostData(4, fMCtrklt);
 161 }
 162
 163 //________________________________________________________
 164 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 165 {
 166   if(track) fTrack = track;
 167   if(!fTrack){
 168     AliWarning("No Track defined.");
 169     return 0x0;
 170   }
 171   TH1 *h = 0x0;
 172   if(!(h = ((TH2I*)fContainer->At(kCluster)))){
 173     AliWarning("No output histogram defined.");
 174     return 0x0;
 175   }
 176
 177   Double_t cov[3];
 178   Float_t x0, y0, z0, dy, dydx, dzdx;
 179   AliTRDseedV1 *fTracklet = 0x0;
 180   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 181     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 182     if(!fTracklet->IsOK()) continue;
 183     x0 = fTracklet->GetX0();
 184
 185     // retrive the track angle with the chamber
 186     y0   = fTracklet->GetYref(0);
 187     z0   = fTracklet->GetZref(0);
 188     dydx = fTracklet->GetYref(1);
 189     dzdx = fTracklet->GetZref(1);
 190     fTracklet->GetCovRef(cov);
 191     Float_t tilt = fTracklet->GetTilt();
 192     AliTRDcluster *c = 0x0;
 193     fTracklet->ResetClusterIter(kFALSE);
 194     while((c = fTracklet->PrevCluster())){
 195       Float_t xc = c->GetX();
 196       Float_t yc = c->GetY();
 197       Float_t zc = c->GetZ();
 198       Float_t dx = x0 - xc;
 199       Float_t yt = y0 - dx*dydx;
 200       Float_t zt = z0 - dx*dzdx;
 201       yc -= tilt*(zc-zt); // tilt correction
 202       dy = yt - yc;
 203
 204       h->Fill(dydx, dy/TMath::Sqrt(cov[0] + c->GetSigmaY2()));
 205
 206       if(fDebugLevel>=1){
 207         // Get z-position with respect to anode wire
 208         //AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 209         Int_t istk = fGeo->GetStack(c->GetDetector());
 210         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 211         Float_t row0 = pp->GetRow0();
 212         Float_t d  =  row0 - zt + pp->GetAnodeWireOffset();
 213         d -= ((Int_t)(2 * d)) / 2.0;
 214         if (d > 0.25) d  = 0.5 - d;
 215
 216 /*        AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 217         fCl->Add(clInfo);
 218         clInfo->SetCluster(c);
 219         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 220         clInfo->SetResolution(dy);
 221         clInfo->SetAnisochronity(d);
 222         clInfo->SetDriftLength(dx);
 223         (*fDebugStream) << "ClusterResiduals"
 224           <<"clInfo.=" << clInfo
 225           << "\n";*/
 226       }
 227     }
 228   }
 229   return h;
 230 }
 231
 232
 233 //________________________________________________________
 234 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 235 {
 236   if(track) fTrack = track;
 237   if(!fTrack){
 238     AliWarning("No Track defined.");
 239     return 0x0;
 240   }
 241   TH1 *h = 0x0;
 242   if(!(h = ((TH2I*)fContainer->At(kTrackletY)))){
 243     AliWarning("No output histogram defined.");
 244     return 0x0;
 245   }
 246
 247   Double_t cov[3], covR[3];
 248   Float_t x, y0, dx, dy, dydx;
 249   AliTRDseedV1 *fTracklet = 0x0;
 250   for(Int_t il=AliTRDgeometry::kNlayer; il--;){
 251     if(!(fTracklet = fTrack->GetTracklet(il))) continue;
 252     if(!fTracklet->IsOK()) continue;
 253     y0   = fTracklet->GetYref(0);
 254     dydx = fTracklet->GetYref(1);
 255     x    = fTracklet->GetX();
 256     dx   = fTracklet->GetX0() - x;
 257     dy   = y0-dx*dydx - fTracklet->GetY();
 258     fTracklet->GetCovAt(x, cov);
 259     fTracklet->GetCovRef(covR);
 260     h->Fill(dydx, dy/*/TMath::Sqrt(cov[0] + covR[0])*/);
 261   }
 262   return h;
 263 }
 264
 265 //________________________________________________________
 266 TH1* AliTRDresolution::PlotTrackletPhi(const AliTRDtrackV1 *track)
 267 {
 268   if(track) fTrack = track;
 269   if(!fTrack){
 270     AliWarning("No Track defined.");
 271     return 0x0;
 272   }
 273   TH1 *h = 0x0;
 274   if(!(h = ((TH2I*)fContainer->At(kTrackletPhi)))){
 275     AliWarning("No output histogram defined.");
 276     return 0x0;
 277   }
 278
 279   AliTRDseedV1 *tracklet = 0x0;
 280   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
 281     if(!(tracklet = fTrack->GetTracklet(il))) continue;
 282     if(!tracklet->IsOK()) continue;
 283     h->Fill(tracklet->GetYref(1), TMath::ATan(tracklet->GetYref(1))-TMath::ATan(tracklet->GetYfit(1)));
 284   }
 285   return h;
 286 }
 287
 288
 289 //________________________________________________________
 290 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 291 {
 292   if(!HasMCdata()){
 293     AliWarning("No MC defined. Results will not be available.");
 294     return 0x0;
 295   }
 296   if(track) fTrack = track;
 297   if(!fTrack){
 298     AliWarning("No Track defined.");
 299     return 0x0;
 300   }
 301   TH1 *h = 0x0;
 302   UChar_t s;
 303   Int_t pdg = fMC->GetPDG(), det=-1;
 304   Int_t label = fMC->GetLabel();
 305   Double_t x, y;
 306   Float_t p, pt, x0, y0, z0, dx, dy, dz, dydx, dzdx;
 307   Double_t covR[3], cov[3];
 308
 309   if(fDebugLevel>=1){
 310     Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
 311     fMC->PropagateKalman(DX, DY, DZ, DPt, COV);
 312     (*fDebugStream) << "MCkalman"
 313       << "pdg="  << pdg
 314       << "dx0="  << DX[0]
 315       << "dx1="  << DX[1]
 316       << "dx2="  << DX[2]
 317       << "dy0="  << DY[0]
 318       << "dy1="  << DY[1]
 319       << "dy2="  << DY[2]
 320       << "dz0="  << DZ[0]
 321       << "dz1="  << DZ[1]
 322       << "dz2="  << DZ[2]
 323       << "dpt0=" << DPt[0]
 324       << "dpt1=" << DPt[1]
 325       << "dpt2=" << DPt[2]
 326       << "\n";
 327   }
 328
 329   AliTRDseedV1 *fTracklet = 0x0;
 330   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 331     if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
 332        !fTracklet->IsOK())*/ continue;
 333
 334     det = fTracklet->GetDetector();
 335     x0  = fTracklet->GetX0();
 336     //radial shift with respect to the MC reference (radial position of the pad plane)
 337     x= fTracklet->GetX();
 338     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, pt, s)) continue;
 339     // MC track position at reference radial position
 340     dx  = x0 - x;
 341     if(fDebugLevel>=1){
 342       (*fDebugStream) << "MC"
 343         << "det="     << det
 344         << "pdg="     << pdg
 345         << "pt="      << pt
 346         << "dx="      << dx
 347         << "x0="      << x0
 348         << "y0="      << y0
 349         << "z0="      << z0
 350         << "dydx="    << dydx
 351         << "dzdx="    << dzdx
 352         << "\n";
 353     }
 354     Float_t yt = y0 - dx*dydx;
 355     Float_t zt = z0 - dx*dzdx;
 356     p = pt*(1.+dzdx*dzdx); // pt -> p
 357
 358     // add Kalman residuals for y, z and pt
 359     dx = fTracklet->GetX0() - x;
 360     Float_t yr = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
 361     dy = yt - yr;
 362     Float_t zr = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
 363     dz = zt - zr;
 364     Float_t tgl = fTracklet->GetTgl();
 365     Float_t dpt = pt - fTracklet->GetMomentum()/(1.+tgl*tgl);
 366     fTracklet->GetCovRef(covR);
 367
 368     ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx, dy);
 369     ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx, dy/TMath::Sqrt(covR[0]));
 370     if(ily==0){
 371       ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx, dz);
 372       ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
 373     } else if(ily==AliTRDgeometry::kNlayer-1) {
 374       ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx, dz);
 375       ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
 376     }
 377     if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt, dpt);
 378     // Fill Debug stream for Kalman track
 379     if(fDebugLevel>=1){
 380       Float_t dydxr = fTracklet->GetYref(1);
 381       (*fDebugStream) << "MCtrack"
 382         << "x="       << x
 383         << "dpt="     << dpt
 384         << "dy="      << dy
 385         << "dz="      << dz
 386         << "dydxr="   << dydxr
 387         << "dzdxr="   << tgl
 388         << "s2y="     << covR[0]
 389         << "s2z="     << covR[2]
 390         << "\n";
 391     }
 392
 393     // recalculate tracklet based on the MC info
 394     AliTRDseedV1 tt(*fTracklet);
 395     tt.SetZref(0, z0);
 396     tt.SetZref(1, dzdx);
 397     tt.Fit(kTRUE);
 398     x= tt.GetX(); // the true one
 399     dx = x0 - x;
 400     yt = y0 - dx*dydx;
 401     zt = z0 - dx*dzdx;
 402     Bool_t rc = tt.IsRowCross();
 403
 404     // add tracklet residuals for y and dydx
 405     Float_t dphi;
 406     tt.GetCovAt(x, cov);
 407     if(!rc){
 408       dy    = yt - tt.GetY();
 409       dphi  = (tt.GetYfit(1) - dydx);
 410       dphi /= 1.- tt.GetYfit(1)*dydx;
 411
 412       ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx, dy);
 413       if(cov[0]>0.) ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx, dy/TMath::Sqrt(cov[0]));
 414       ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx, dphi*TMath::RadToDeg());
 415     } else {
 416       // add tracklet residuals for z
 417       dz = tt.GetZ() - zt;
 418       ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx, dz);
 419       if(cov[2]>0.) ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx, dz/TMath::Sqrt(cov[2]));
 420     }
 421
 422     // Fill Debug stream for tracklet
 423     if(fDebugLevel>=1){
 424       (*fDebugStream) << "MCtracklet"
 425         << "rc="    << rc
 426         << "x="     << x
 427         << "dy="    << dy
 428         << "dz="    << dz
 429         << "dphi="  << dphi
 430         << "s2y="   << cov[0]
 431         << "s2z="   << cov[2]
 432         << "\n";
 433     }
 434
 435     Int_t istk = AliTRDgeometry::GetStack(det);
 436     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 437     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 438     Float_t tilt = fTracklet->GetTilt();
 439
 440     AliTRDcluster *c = 0x0;
 441     fTracklet->ResetClusterIter(kFALSE);
 442     while((c = fTracklet->PrevCluster())){
 443       Float_t  q = TMath::Abs(c->GetQ());
 444       AliTRDseedV1::GetClusterXY(c,x,y);
 445       //x = c->GetX(); y = c->GetY();
 446       Float_t xc = x;
 447       Float_t yc = y;
 448       Float_t zc = c->GetZ();
 449       dx = x0 - xc;
 450       yt = y0 - dx*dydx;
 451       zt = z0 - dx*dzdx;
 452       dy = yt - (yc - tilt*(zc-zt));
 453
 454       // Fill Histograms
 455       if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx, dy);
 456
 457       // Fill calibration container
 458       Float_t d = zr0 - zt;
 459       d -= ((Int_t)(2 * d)) / 2.0;
 460       if (d > 0.25) d  = 0.5 - d;
 461       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 462       fMCcl->Add(clInfo);
 463       clInfo->SetCluster(c);
 464       clInfo->SetMC(pdg, label);
 465       clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 466       clInfo->SetResolution(dy);
 467       clInfo->SetAnisochronity(d);
 468       clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
 469       clInfo->SetTilt(tilt);
 470
 471       // Fill Debug Tree
 472       if(fDebugLevel>=2){
 473         //clInfo->Print();
 474         (*fDebugStream) << "MCcluster"
 475           <<"clInfo.=" << clInfo
 476           << "\n";
 477       }
 478     }
 479   }
 480   return h;
 481 }
 482
 483
 484 //________________________________________________________
 485 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
 486 {
 487   Float_t y[2] = {0., 0.};
 488   TBox *b = 0x0;
 489   TAxis *ax = 0x0;
 490   TGraphErrors *g = 0x0;
 491   switch(ifig){
 492   case kCluster:
 493     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 494     g->Draw("apl");
 495     ax = g->GetHistogram()->GetYaxis();
 496     y[0] = -0.5; y[1] = 2.5;
 497     ax->SetRangeUser(y[0], y[1]);
 498     ax->SetTitle("Cluster-Track Pulls #sigma/#mu [mm]");
 499     ax = g->GetHistogram()->GetXaxis();
 500     ax->SetTitle("tg(#phi)");
 501     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 502     g->Draw("pl");
 503     b = new TBox(-.15, y[0], .15, y[1]);
 504     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 505     b->SetLineColor(0); b->Draw();
 506     return kTRUE;
 507   case kTrackletY:
 508     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 509     g->Draw("apl");
 510     ax = g->GetHistogram()->GetYaxis();
 511     ax->SetRangeUser(-.5, 3.);
 512     ax->SetTitle("Tracklet-Track Y-Pulls #sigma/#mu [mm]");
 513     ax = g->GetHistogram()->GetXaxis();
 514     ax->SetTitle("tg(#phi)");
 515     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 516     g->Draw("pl");
 517     b = new TBox(-.15, y[0], .15, y[1]);
 518     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 519     b->SetLineColor(0); b->Draw();
 520     return kTRUE;
 521   case kTrackletPhi:
 522     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 523     g->Draw("apl");
 524     ax = g->GetHistogram()->GetYaxis();
 525     ax->SetRangeUser(-.5, 2.);
 526     ax->SetTitle("Tracklet-Track Phi-Residuals #sigma/#mu [rad]");
 527     ax = g->GetHistogram()->GetXaxis();
 528     ax->SetTitle("tg(#phi)");
 529     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 530     g->Draw("pl");
 531     b = new TBox(-.15, y[0], .15, y[1]);
 532     b->SetFillStyle(3002);b->SetFillColor(kGreen);
 533     b->SetLineColor(0); b->Draw();
 534     return kTRUE;
 535   case kMCcluster:
 536     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 537     ax = g->GetHistogram()->GetYaxis();
 538     y[0] = -50.; y[1] = 600.;
 539     ax->SetRangeUser(y[0], y[1]);
 540     ax->SetTitle("Y_{cluster} #sigma/#mu [#mum]");
 541     ax = g->GetHistogram()->GetXaxis();
 542     ax->SetRangeUser(-.3, .3);
 543     ax->SetTitle("tg(#phi)");
 544     g->Draw("apl");
 545     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 546     g->Draw("pl");
 547     b = new TBox(-.15, y[0], .15, y[1]);
 548     b->SetFillStyle(3002);b->SetFillColor(kBlue);
 549     b->SetLineColor(0); b->Draw();
 550     return kTRUE;
 551   case kMCtrackletY:
 552     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 553     ax = g->GetHistogram()->GetYaxis();
 554     y[0] = -50.; y[1] = 250.;
 555     ax->SetRangeUser(y[0], y[1]);
 556     ax->SetTitle("Y_{tracklet} #sigma/#mu [#mum]");
 557     ax = g->GetHistogram()->GetXaxis();
 558     ax->SetRangeUser(-.2, .2);
 559     ax->SetTitle("tg(#phi)");
 560     g->Draw("apl");
 561     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 562     g->Draw("pl");
 563     b = new TBox(-.15, y[0], .15, y[1]);
 564     b->SetFillStyle(3002);b->SetFillColor(kBlue);
 565     b->SetLineColor(0); b->Draw();
 566     return kTRUE;
 567   case kMCtrackletZ:
 568     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 569     ax = g->GetHistogram()->GetYaxis();
 570     ax->SetRangeUser(-50., 700.);
 571     ax->SetTitle("Z_{tracklet} #sigma/#mu [#mum]");
 572     ax = g->GetHistogram()->GetXaxis();
 573     ax->SetTitle("tg(#theta)");
 574     g->Draw("apl");
 575     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 576     g->Draw("pl");
 577     return kTRUE;
 578   case kMCtrackletPhi:
 579     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 580     ax = g->GetHistogram()->GetYaxis();
 581     y[0] = -.05; y[1] = .2;
 582     ax->SetRangeUser(y[0], y[1]);
 583     ax->SetTitle("#Phi_{tracklet} #sigma/#mu [deg]");
 584     ax = g->GetHistogram()->GetXaxis();
 585     ax->SetTitle("tg(#phi)");
 586     g->Draw("apl");
 587     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 588     g->Draw("pl");
 589     return kTRUE;
 590   case kMCtrackY:
 591     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 592     ax = g->GetHistogram()->GetYaxis();
 593     y[0] = -50.; y[1] = 200.;
 594     ax->SetRangeUser(y[0], y[1]);
 595     ax->SetTitle("Y_{track} #sigma/#mu [#mum]");
 596     ax = g->GetHistogram()->GetXaxis();
 597     ax->SetRangeUser(-.2, .2);
 598     ax->SetTitle("tg(#phi)");
 599     g->Draw("apl");
 600     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 601     g->Draw("pl");
 602     b = new TBox(-.15, y[0], .15, y[1]);
 603     b->SetFillStyle(3002);b->SetFillColor(kBlue);
 604     b->SetLineColor(0); b->Draw();
 605     return kTRUE;
 606   case kMCtrackZIn:
 607   case kMCtrackZOut:
 608     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 609     ax = g->GetHistogram()->GetYaxis();
 610     ax->SetRangeUser(-500., 2000.);
 611     ax->SetTitle(Form("Z_{track}^{%s} #sigma/#mu [#mum]", ifig==kMCtrackZIn ? "in" : "out"));
 612     ax = g->GetHistogram()->GetXaxis();
 613     ax->SetTitle("tg(#theta)");
 614     g->Draw("apl");
 615     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 616     g->Draw("pl");
 617     return kTRUE;
 618   case kMCtrackPt:
 619     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 620     ax = g->GetHistogram()->GetYaxis();
 621     ax->SetRangeUser(-.5, 2.);
 622     ax->SetTitle("#epsilon_{P_{t}}^{track} / #mu [%]");
 623     ax = g->GetHistogram()->GetXaxis();
 624     ax->SetTitle("1/p_{t}");
 625     g->Draw("apl");
 626     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 627     g->Draw("pl");
 628     return kTRUE;
 629   case kMCtrackletYPull:
 630   case kMCtrackletZPull:
 631   case kMCtrackYPull:
 632   case kMCtrackZInPull:
 633   case kMCtrackZOutPull:
 634     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 635     ax = g->GetHistogram()->GetYaxis();
 636     ax->SetRangeUser(-.5, 2.);
 637     ax->SetTitle("MC Pulls");
 638     g->Draw("apl");
 639     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 640     g->Draw("pl");
 641     return kTRUE;
 642   }
 643   AliInfo(Form("Reference plot [%d] missing result", ifig));
 644   return kFALSE;
 645 }
 646
 647
 648 //________________________________________________________
 649 Bool_t AliTRDresolution::PostProcess()
 650 {
 651   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 652   if (!fContainer) {
 653     Printf("ERROR: list not available");
 654     return kFALSE;
 655   }
 656   fNRefFigures = fContainer->GetEntriesFast();
 657   TGraphErrors *gm= 0x0, *gs= 0x0;
 658   if(!fGraphS){
 659     fGraphS = new TObjArray(fNRefFigures);
 660     fGraphS->SetOwner();
 661     for(Int_t ig=0; ig<fNRefFigures; ig++){
 662       gs = new TGraphErrors();
 663       gs->SetLineColor(kRed);
 664       gs->SetMarkerStyle(23);
 665       gs->SetMarkerColor(kRed);
 666       gs->SetNameTitle(Form("s_%d", ig), "");
 667       fGraphS->AddAt(gs, ig);
 668     }
 669   }
 670   if(!fGraphM){
 671     fGraphM = new TObjArray(fNRefFigures);
 672     fGraphM->SetOwner();
 673     for(Int_t ig=0; ig<fNRefFigures; ig++){
 674       gm = new TGraphErrors();
 675       gm->SetLineColor(kBlue);
 676       gm->SetMarkerStyle(7);
 677       gm->SetMarkerColor(kBlue);
 678       gm->SetNameTitle(Form("m_%d", ig), "");
 679       fGraphM->AddAt(gm, ig);
 680     }
 681   }
 682
 683   // DEFINE MODELS
 684   // simple gauss
 685   TF1 f("f1", "gaus", -.5, .5);
 686   // gauss on a constant background
 687   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 688   // gauss on a gauss background
 689   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 690
 691
 692   //PROCESS EXPERIMENTAL DISTRIBUTIONS
 693
 694   // Clusters residuals
 695   Process(kCluster, &f);
 696
 697   // Tracklet y residuals
 698   Process(kTrackletY, &f);
 699
 700   // Tracklet phi residuals
 701   Process(kTrackletPhi, &f);
 702
 703   if(!HasMCdata()) return kTRUE;
 704
 705
 706   //PROCESS MC RESIDUAL DISTRIBUTIONS
 707
 708   // cluster y resolution
 709   Process(kMCcluster, &f, 1.e4);
 710
 711   // tracklet resolution
 712   Process(kMCtrackletY, &f, 1.e4); // y
 713   Process(kMCtrackletZ, &f, 1.e4); // z
 714   Process(kMCtrackletPhi, &f); // phi
 715
 716   // tracklet pulls
 717   Process(kMCtrackletYPull, &f); // y
 718   Process(kMCtrackletZPull, &f); // z
 719
 720   // track resolution
 721   Process(kMCtrackY, &f, 1.e4);    // y
 722   Process(kMCtrackZIn, &f, 1.e4);  // z towards TPC
 723   Process(kMCtrackZOut, &f, 1.e4); // z towards TOF
 724
 725   // track pulls
 726   Process(kMCtrackYPull, &f);    // y
 727   Process(kMCtrackZInPull, &f);  // z towards TPC
 728   Process(kMCtrackZOutPull, &f); // z towards TOF
 729
 730   // track Pt resolution
 731   TH2I *h2 = (TH2I*)fContainer->At(kMCtrackPt);
 732   TAxis *ax = h2->GetXaxis();
 733   gm = (TGraphErrors*)fGraphM->At(kMCtrackPt);
 734   gs = (TGraphErrors*)fGraphS->At(kMCtrackPt);
 735   TF1 fg("fg", "gaus", -1.5, 1.5);
 736   TF1 fl("fl", "landau", -4., 15.);
 737   TF1 fgl("fgl", "gaus(0)+landau(3)", -5., 20.);
 738   for(Int_t ip=1; ip<=ax->GetNbins(); ip++){
 739     TH1D *h = h2->ProjectionY("ppt", ip, ip);
 740     if(h->GetEntries()<70) continue;
 741
 742     h->Fit(&fg, "QN", "", -1.5, 1.5);
 743     fgl.SetParameter(0, fg.GetParameter(0));
 744     fgl.SetParameter(1, fg.GetParameter(1));
 745     fgl.SetParameter(2, fg.GetParameter(2));
 746     h->Fit(&fl, "QN", "", -4., 15.);
 747     fgl.SetParameter(3, fl.GetParameter(0));
 748     fgl.SetParameter(4, fl.GetParameter(1));
 749     fgl.SetParameter(5, fl.GetParameter(2));
 750
 751     h->Fit(&fgl, "NQ", "", -5., 20.);
 752
 753     Float_t invpt = ax->GetBinCenter(ip);
 754     Int_t ip = gm->GetN();
 755     gm->SetPoint(ip, invpt, fgl.GetParameter(1));
 756     gm->SetPointError(ip, 0., fgl.GetParError(1));
 757     gs->SetPoint(ip, invpt, fgl.GetParameter(2)*invpt);
 758     gs->SetPointError(ip, 0., fgl.GetParError(2));
 759     // fgl.GetParameter(4) // Landau MPV
 760     // fgl.GetParameter(5) // Landau Sigma
 761   }
 762
 763
 764   return kTRUE;
 765 }
 766
 767
 768 //________________________________________________________
 769 void AliTRDresolution::Terminate(Option_t *)
 770 {
 771   if(fDebugStream){
 772     delete fDebugStream;
 773     fDebugStream = 0x0;
 774     fDebugLevel = 0;
 775   }
 776   if(HasPostProcess()) PostProcess();
 777 }
 778
 779 //________________________________________________________
 780 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
 781 {
 782 // Helper function to avoid duplication of code
 783 // Make first guesses on the fit parameters
 784
 785   // find the intial parameters of the fit !! (thanks George)
 786   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 787   Double_t sum = 0.;
 788   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 789   f->SetParLimits(0, 0., 3.*sum);
 790   f->SetParameter(0, .9*sum);
 791
 792   f->SetParLimits(1, -.2, .2);
 793   f->SetParameter(1, -0.1);
 794
 795   f->SetParLimits(2, 0., 4.e-1);
 796   f->SetParameter(2, 2.e-2);
 797   if(f->GetNpar() <= 4) return;
 798
 799   f->SetParLimits(3, 0., sum);
 800   f->SetParameter(3, .1*sum);
 801
 802   f->SetParLimits(4, -.3, .3);
 803   f->SetParameter(4, 0.);
 804
 805   f->SetParLimits(5, 0., 1.e2);
 806   f->SetParameter(5, 2.e-1);
 807 }
 808
 809 //________________________________________________________
 810 TObjArray* AliTRDresolution::Histos()
 811 {
 812   if(fContainer) return fContainer;
 813
 814   fContainer  = new TObjArray(kNhistos);
 815   //fContainer->SetOwner(kTRUE);
 816
 817   TH1 *h = 0x0;
 818   // cluster to tracklet residuals [2]
 819   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
 820     h = new TH2I("hCl", "Clusters-Track Residuals", 21, -.33, .33, 100, -.5, .5);
 821     h->GetXaxis()->SetTitle("tg(#phi)");
 822     h->GetYaxis()->SetTitle("#Delta y [cm]");
 823     h->GetZaxis()->SetTitle("entries");
 824   } else h->Reset();
 825   fContainer->AddAt(h, kCluster);
 826
 827   // tracklet to track residuals [2]
 828   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
 829     h = new TH2I("hTrkltY", "Tracklets-Track Residuals (Y)", 21, -.33, .33, 100, -.5, .5);
 830     h->GetXaxis()->SetTitle("tg(#phi)");
 831     h->GetYaxis()->SetTitle("#Delta y [cm]");
 832     h->GetZaxis()->SetTitle("entries");
 833   } else h->Reset();
 834   fContainer->AddAt(h, kTrackletY);
 835
 836   // tracklet to track residuals angular [2]
 837   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
 838     h = new TH2I("hTrkltPhi", "Tracklets-Track Residuals (#Phi)", 21, -.33, .33, 100, -.5, .5);
 839     h->GetXaxis()->SetTitle("tg(#phi)");
 840     h->GetYaxis()->SetTitle("#Delta phi [#circ]");
 841     h->GetZaxis()->SetTitle("entries");
 842   } else h->Reset();
 843   fContainer->AddAt(h, kTrackletPhi);
 844
 845
 846   // Resolution histos
 847   if(!HasMCdata()) return fContainer;
 848
 849   // cluster y resolution [0]
 850   if(!(h = (TH2I*)gROOT->FindObject("hMCcl"))){
 851     h = new TH2I("hMCcl", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
 852     h->GetXaxis()->SetTitle("tg(#phi)");
 853     h->GetYaxis()->SetTitle("#Delta y [cm]");
 854     h->GetZaxis()->SetTitle("entries");
 855   } else h->Reset();
 856   fContainer->AddAt(h, kMCcluster);
 857
 858   // tracklet y resolution [0]
 859   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltY"))){
 860     h = new TH2I("hMCtrkltY", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.2, .2);
 861     h->GetXaxis()->SetTitle("tg(#phi)");
 862     h->GetYaxis()->SetTitle("#Delta y [cm]");
 863     h->GetZaxis()->SetTitle("entries");
 864   } else h->Reset();
 865   fContainer->AddAt(h, kMCtrackletY);
 866
 867   // tracklet y resolution [0]
 868   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltYPull"))){
 869     h = new TH2I("hMCtrkltYPull", "Tracklet Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
 870     h->GetXaxis()->SetTitle("tg(#phi)");
 871     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
 872     h->GetZaxis()->SetTitle("entries");
 873   } else h->Reset();
 874   fContainer->AddAt(h, kMCtrackletYPull);
 875
 876   // tracklet y resolution [0]
 877   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZ"))){
 878     h = new TH2I("hMCtrkltZ", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
 879     h->GetXaxis()->SetTitle("tg(#theta)");
 880     h->GetYaxis()->SetTitle("#Delta z [cm]");
 881     h->GetZaxis()->SetTitle("entries");
 882   } else h->Reset();
 883   fContainer->AddAt(h, kMCtrackletZ);
 884
 885   // tracklet y resolution [0]
 886   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZPull"))){
 887     h = new TH2I("hMCtrkltZPull", "Tracklet Pulls (Z)", 31, -.48, .48, 100, -3.5, 3.5);
 888     h->GetXaxis()->SetTitle("tg(#theta)");
 889     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
 890     h->GetZaxis()->SetTitle("entries");
 891   } else h->Reset();
 892   fContainer->AddAt(h, kMCtrackletZPull);
 893
 894   // tracklet angular resolution [1]
 895   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltPhi"))){
 896     h = new TH2I("hMCtrkltPhi", "Tracklet Resolution (#Phi)", 31, -.48, .48, 100, -10., 10.);
 897     h->GetXaxis()->SetTitle("tg(#phi)");
 898     h->GetYaxis()->SetTitle("#Delta #phi [deg]");
 899     h->GetZaxis()->SetTitle("entries");
 900   } else h->Reset();
 901   fContainer->AddAt(h, kMCtrackletPhi);
 902
 903   // Kalman track y resolution
 904   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkY"))){
 905     h = new TH2I("hMCtrkY", "Kalman Track Resolution (Y)", 31, -.48, .48, 100, -.2, .2);
 906     h->GetXaxis()->SetTitle("tg(#phi)");
 907     h->GetYaxis()->SetTitle("#Delta y [cm]");
 908     h->GetZaxis()->SetTitle("entries");
 909   } else h->Reset();
 910   fContainer->AddAt(h, kMCtrackY);
 911
 912   // Kalman track y resolution
 913   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYPull"))){
 914     h = new TH2I("hMCtrkYPull", "Kalman Track Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
 915     h->GetXaxis()->SetTitle("tg(#phi)");
 916     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
 917     h->GetZaxis()->SetTitle("entries");
 918   } else h->Reset();
 919   fContainer->AddAt(h, kMCtrackYPull);
 920
 921   // Kalman track Z resolution
 922   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZIn"))){
 923     h = new TH2I("hMCtrkZIn", "Kalman Track Resolution (Zin)", 20, -1., 1., 100, -1., 1.);
 924     h->GetXaxis()->SetTitle("tg(#theta)");
 925     h->GetYaxis()->SetTitle("#Delta z [cm]");
 926     h->GetZaxis()->SetTitle("entries");
 927   } else h->Reset();
 928   fContainer->AddAt(h, kMCtrackZIn);
 929
 930   // Kalman track Z resolution
 931   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOut"))){
 932     h = new TH2I("hMCtrkZOut", "Kalman Track Resolution (Zout)", 20, -1., 1., 100, -1., 1.);
 933     h->GetXaxis()->SetTitle("tg(#theta)");
 934     h->GetYaxis()->SetTitle("#Delta z [cm]");
 935     h->GetZaxis()->SetTitle("entries");
 936   } else h->Reset();
 937   fContainer->AddAt(h, kMCtrackZOut);
 938
 939   // Kalman track Z resolution
 940   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZInPull"))){
 941     h = new TH2I("hMCtrkZInPull", "Kalman Track Pulls (Zin)", 20, -1., 1., 100, -4.5, 4.5);
 942     h->GetXaxis()->SetTitle("tg(#theta)");
 943     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
 944     h->GetZaxis()->SetTitle("entries");
 945   } else h->Reset();
 946   fContainer->AddAt(h, kMCtrackZInPull);
 947
 948   // Kalman track Z resolution
 949   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOutPull"))){
 950     h = new TH2I("hMCtrkZOutPull", "Kalman Track Pulls (Zout)", 20, -1., 1., 100, -4.5, 4.5);
 951     h->GetXaxis()->SetTitle("tg(#theta)");
 952     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
 953     h->GetZaxis()->SetTitle("entries");
 954   } else h->Reset();
 955   fContainer->AddAt(h, kMCtrackZOutPull);
 956
 957   // Kalman track Pt resolution
 958   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPt"))){
 959     h = new TH2I("hMCtrkPt", "Kalman Track Resolution (Pt)", 100, 0., 2., 150, -5., 20.);
 960     h->GetXaxis()->SetTitle("1/p_{t}");
 961     h->GetYaxis()->SetTitle("#Delta p_{t} [GeV/c]");
 962     h->GetZaxis()->SetTitle("entries");
 963   } else h->Reset();
 964   fContainer->AddAt(h, kMCtrackPt);
 965
 966   return fContainer;
 967 }
 968
 969
 970 //________________________________________________________
 971 Bool_t AliTRDresolution::Process(ETRDresolutionPlot ip, TF1 *f, Float_t k)
 972 {
 973   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
 974   Bool_t kBUILD = kFALSE;
 975   if(!f){
 976     f = new TF1("f1", "gaus", -.5, .5);
 977     kBUILD = kTRUE;
 978   }
 979
 980   TH2I *h2 = 0x0;
 981   if(!(h2 = (TH2I *)(fContainer->At(ip)))) return kFALSE;
 982   TGraphErrors *gm = 0x0, *gs = 0x0;
 983   if(!(gm=(TGraphErrors*)fGraphM->At(ip))) return kFALSE;
 984   if(gm->GetN()) for(Int_t ip=gm->GetN(); ip--;) gm->RemovePoint(ip);
 985   if(!(gs=(TGraphErrors*)fGraphS->At(ip))) return kFALSE;
 986   if(gs->GetN()) for(Int_t ip=gs->GetN(); ip--;) gs->RemovePoint(ip);
 987   Char_t pn[10]; sprintf(pn, "p%02d", ip);
 988   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 989     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
 990     TH1D *h = h2->ProjectionY(pn, ibin, ibin);
 991     if(h->GetEntries()<100) continue;
 992     AdjustF1(h, f);
 993
 994     h->Fit(f, "QN");
 995
 996     Int_t ip = gm->GetN();
 997     gm->SetPoint(ip, x, k*f->GetParameter(1));
 998     gm->SetPointError(ip, 0., k*f->GetParError(1));
 999     gs->SetPoint(ip, x, k*f->GetParameter(2));
1000     gs->SetPointError(ip, 0., k*f->GetParError(2));
1001   }
1002
1003   if(kBUILD) delete f;
1004   return kTRUE;
1005 }
1006
1007 //________________________________________________________
1008 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
1009 {
1010
1011   fReconstructor->SetRecoParam(r);
1012 }