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