PWG1/TRD/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();
  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 <TSystem.h>
  52
  53 #include <TROOT.h>
  54 #include <TObjArray.h>
  55 #include <TH3.h>
  56 #include <TH2.h>
  57 #include <TH1.h>
  58 #include <TF1.h>
  59 #include <TCanvas.h>
  60 #include <TGaxis.h>
  61 #include <TBox.h>
  62 #include <TLegend.h>
  63 #include <TGraphErrors.h>
  64 #include <TGraphAsymmErrors.h>
  65 #include <TMath.h>
  66 #include <TMatrixT.h>
  67 #include <TVectorT.h>
  68 #include <TTreeStream.h>
  69 #include <TGeoManager.h>
  70 #include <TDatabasePDG.h>
  71
  72 #include "AliPID.h"
  73 #include "AliESDtrack.h"
  74
  75 #include "AliTRDresolution.h"
  76 #include "AliTRDgeometry.h"
  77 #include "AliTRDpadPlane.h"
  78 #include "AliTRDcluster.h"
  79 #include "AliTRDseedV1.h"
  80 #include "AliTRDtrackV1.h"
  81 #include "AliTRDReconstructor.h"
  82 #include "AliTRDrecoParam.h"
  83 #include "AliTRDpidUtil.h"
  84
  85 #include "info/AliTRDclusterInfo.h"
  86
  87 ClassImp(AliTRDresolution)
  88
  89 UChar_t const AliTRDresolution::fgNhistos[kNviews] = {
  90   2, 2, 5, 5,
  91   2, 5, 12, 2, 14
  92 };
  93 UChar_t const AliTRDresolution::fgNproj[kNviews] = {
  94   2, 2, 5, 5,
  95   4, 7, 12, 2, 14
  96 };
  97 Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
  98      "Charge"
  99     ,"Cluster2Track"
 100     ,"Tracklet2Track"
 101     ,"Tracklet2TPC"
 102     ,"Cluster2MC"
 103     ,"Tracklet2MC"
 104     ,"TPC2MC"
 105     ,"TOF/HMPID2MC"
 106     ,"TRD2MC"
 107 };
 108 UChar_t const AliTRDresolution::fgNcomp[kNprojs] = {
 109   1, 1, //2,
 110   1, 1, //2,
 111   1, 1, 1, 1, 1, //5,
 112   1, 1, 1, 1, 1, 1, 1, //5,
 113 // MC
 114   1, 1, 1, 1,    //4,
 115   1, 1, 1, 1, 1, //5,
 116   1, 1, 1, 1, 1, 1, 1, 1, 11, 11, 11, 11, //12,
 117   1, 1, //2,
 118   1, 1, 1, 1, 1, 1, 1, 1, 66, 66, 66, 66, 66, 66 //14
 119 };
 120 Char_t const *AliTRDresolution::fgAxTitle[kNprojs][4] = {
 121   // Charge
 122   {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
 123  ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
 124   // Clusters to Kalman
 125  ,{"Cluster2Track residuals", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 126  ,{"Cluster2Track  pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 127   // TRD tracklet to Kalman fit
 128  ,{"Tracklet2Track Y residuals", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 129  ,{"Tracklet2Track Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 130  ,{"Tracklet2Track Z residuals", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
 131  ,{"Tracklet2Track Z pulls", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
 132  ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
 133   // TPC track 2 first TRD tracklet
 134  ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 135  ,{"Tracklet2Track Y pulls @ TRDin", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 136  ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
 137  ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
 138  ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
 139   // MC cluster
 140  ,{"MC Cluster Y resolution (p_{t}<1 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 141  ,{"MC Cluster Y resolution (1<p_{t}<2 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 142  ,{"MC Cluster Y resolution (p_{t}>3 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 143  ,{"MC Cluster Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 144   // MC tracklet
 145  ,{"MC Tracklet Y resolution (p_{t}<1 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 146  ,{"MC Tracklet Y resolution (1<p_{t}<2 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
 147  ,{"MC Tracklet Y resolution (p_{t}>3 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y}[#mum]"}
 148  ,{"MC Tracklet Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 149  ,{"MC Tracklet Cross Z resolution", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
 150  ,{"MC Tracklet Cross Z pulls", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
 151  ,{"MC Tracklet Phi resolution", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
 152   // MC track TPC
 153  ,{"Y resolution @ TRDin", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y}[#mum]"}
 154  ,{"Y pulls @ TRDin", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
 155  ,{"Z resolution @ TRDin", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
 156  ,{"Z pulls @ TRDin", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
 157  ,{"Phi resolution @ TRDin", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
 158  ,{"SNP pulls @ TRDin", "tg(#phi)", "#mu_{snp}", "#sigma_{snp}"}
 159  ,{"Theta resolution @ TRDin", "tg(#theta)", "#mu_{#theta} [mrad]", "#sigma_{#theta} [mrad]"}
 160  ,{"TGL pulls @ TRDin", "tg(#theta)", "#mu_{tgl}", "#sigma_{tgl}"}
 161  ,{"P_{t} resolution @ TRDin", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
 162  ,{"1/P_{t} pulls @ TRDin", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
 163  ,{"P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
 164  ,{"P pulls @ TRDin", "p^{MC} [GeV/c]", "1/p^{REC}-1/p^{MC}", "MC PULL: #sigma^{TPC}(#Deltap/#sigma_{p})"}
 165   // MC track TOF
 166  ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{TOF} [#mum]", "MC: #sigma_{z}^{TOF} [#mum]"}
 167  ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{TOF}", "MC PULL: #sigma_{z}^{TOF}"}
 168   // MC track in TRD
 169  ,{"TRD track MC Y resolution", "tg(#phi)", "#mu_{y}^{Trk} [#mum]", "#sigma_{y}^{Trk} [#mum]"}
 170  ,{"TRD track MC Y pulls", "tg(#phi)", "#mu_{y}^{Trk}", "#sigma_{y}^{Trk}"}
 171  ,{"TRD track MC Z resolution", "tg(#theta)", "#mu_{z}^{Trk} [#mum]", "#sigma_{z}^{Trk} [#mum]"}
 172  ,{"TRD track MC Z pulls", "tg(#theta)", "#mu_{z}^{Trk}", "#sigma_{z}^{Trk}"}
 173  ,{"TRD track MC Phi resolution", "tg(#phi)", "#mu_{#phi}^{Trk} [mrad]", "#sigma_{#phi}^{Trk} [mrad]"}
 174  ,{"TRD track MC SNP pulls", "tg(#phi)", "#mu_{snp}^{Trk}", "#sigma_{snp}^{Trk}"}
 175  ,{"TRD track MC Theta resolution", "tg(#theta)", "#mu_{#theta}^{Trk} [mrad]", "#sigma_{#theta}^{Trk} [mrad]"}
 176  ,{"TRD track MC TGL pulls", "tg(#theta)", "#mu_{tgl}^{Trk}", "#sigma_{tgl}^{Trk}"}
 177  ,{"P_{t} resolution TRD Layer", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
 178  ,{"1/P_{t} pulls TRD Layer", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
 179  ,{"P resolution TRD Layer", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
 180  ,{"[SA] P_{t} resolution TRD Layer", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
 181  ,{"[SA] 1/P_{t} pulls TRD Layer", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{Trk}"}
 182  ,{"[SA] P resolution TRD Layer", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{Trk}(#Deltap/p^{MC}) [%]"}
 183 };
 184
 185 //________________________________________________________
 186 AliTRDresolution::AliTRDresolution()
 187   :AliTRDrecoTask()
 188   ,fStatus(0)
 189   ,fIdxPlot(0)
 190   ,fIdxFrame(0)
 191   ,fReconstructor(NULL)
 192   ,fGeo(NULL)
 193   ,fDBPDG(NULL)
 194   ,fGraphS(NULL)
 195   ,fGraphM(NULL)
 196   ,fCl(NULL)
 197   ,fTrklt(NULL)
 198   ,fMCcl(NULL)
 199   ,fMCtrklt(NULL)
 200 {
 201   //
 202   // Default constructor
 203   //
 204   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
 205 }
 206
 207 AliTRDresolution::AliTRDresolution(char* name)
 208   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
 209   ,fStatus(0)
 210   ,fIdxPlot(0)
 211   ,fIdxFrame(0)
 212   ,fReconstructor(NULL)
 213   ,fGeo(NULL)
 214   ,fDBPDG(NULL)
 215   ,fGraphS(NULL)
 216   ,fGraphM(NULL)
 217   ,fCl(NULL)
 218   ,fTrklt(NULL)
 219   ,fMCcl(NULL)
 220   ,fMCtrklt(NULL)
 221 {
 222   //
 223   // Default constructor
 224   //
 225
 226   fReconstructor = new AliTRDReconstructor();
 227   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 228   fGeo = new AliTRDgeometry();
 229
 230   InitFunctorList();
 231
 232   DefineOutput(1, TObjArray::Class()); // cluster2track
 233   DefineOutput(2, TObjArray::Class()); // tracklet2track
 234   DefineOutput(3, TObjArray::Class()); // cluster2mc
 235   DefineOutput(4, TObjArray::Class()); // tracklet2mc
 236 }
 237
 238 //________________________________________________________
 239 AliTRDresolution::~AliTRDresolution()
 240 {
 241   //
 242   // Destructor
 243   //
 244
 245   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 246   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 247   delete fGeo;
 248   delete fReconstructor;
 249   if(gGeoManager) delete gGeoManager;
 250   if(fCl){fCl->Delete(); delete fCl;}
 251   if(fTrklt){fTrklt->Delete(); delete fTrklt;}
 252   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 253   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
 254 }
 255
 256
 257 //________________________________________________________
 258 void AliTRDresolution::UserCreateOutputObjects()
 259 {
 260   // spatial resolution
 261   OpenFile(1, "RECREATE");
 262
 263   fContainer = Histos();
 264
 265   fCl = new TObjArray();
 266   fCl->SetOwner(kTRUE);
 267   fTrklt = new TObjArray();
 268   fTrklt->SetOwner(kTRUE);
 269   fMCcl = new TObjArray();
 270   fMCcl->SetOwner(kTRUE);
 271   fMCtrklt = new TObjArray();
 272   fMCtrklt->SetOwner(kTRUE);
 273 }
 274
 275 //________________________________________________________
 276 void AliTRDresolution::UserExec(Option_t *opt)
 277 {
 278   //
 279   // Execution part
 280   //
 281
 282   fCl->Delete();
 283   fTrklt->Delete();
 284   fMCcl->Delete();
 285   fMCtrklt->Delete();
 286
 287   AliTRDrecoTask::UserExec(opt);
 288
 289   PostData(1, fCl);
 290   PostData(2, fTrklt);
 291   PostData(3, fMCcl);
 292   PostData(4, fMCtrklt);
 293 }
 294
 295 //________________________________________________________
 296 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
 297 {
 298   //
 299   // Plots the charge distribution
 300   //
 301
 302   if(track) fkTrack = track;
 303   if(!fkTrack){
 304     AliDebug(2, "No Track defined.");
 305     return NULL;
 306   }
 307   TObjArray *arr = NULL;
 308   if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
 309     AliWarning("No output container defined.");
 310     return NULL;
 311   }
 312   TH3S* h = NULL;
 313
 314   AliTRDseedV1 *fTracklet = NULL;
 315   AliTRDcluster *c = NULL;
 316   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 317     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 318     if(!fTracklet->IsOK()) continue;
 319     Float_t x0 = fTracklet->GetX0();
 320     Float_t dq, dl;
 321     for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
 322       if(!(c = fTracklet->GetClusters(itb))){
 323         if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
 324       }
 325       dq = fTracklet->GetdQdl(itb, &dl);
 326       dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
 327       (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
 328     }
 329
 330 //     if(!HasMCdata()) continue;
 331 //     UChar_t s;
 332 //     Float_t pt0, y0, z0, dydx0, dzdx0;
 333 //     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
 334
 335   }
 336   return h;
 337 }
 338
 339
 340 //________________________________________________________
 341 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 342 {
 343   //
 344   // Plot the cluster distributions
 345   //
 346
 347   if(track) fkTrack = track;
 348   if(!fkTrack){
 349     AliDebug(2, "No Track defined.");
 350     return NULL;
 351   }
 352   TObjArray *arr = NULL;
 353   if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
 354     AliWarning("No output container defined.");
 355     return NULL;
 356   }
 357   ULong_t status = fkESD ? fkESD->GetStatus():0;
 358
 359   Double_t covR[7], cov[3];
 360   Float_t x0, y0, z0, dy, dz, dydx, dzdx;
 361   AliTRDseedV1 *fTracklet(NULL);
 362   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 363     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 364     if(!fTracklet->IsOK()) continue;
 365     x0 = fTracklet->GetX0();
 366
 367     // retrive the track angle with the chamber
 368     y0   = fTracklet->GetYref(0);
 369     z0   = fTracklet->GetZref(0);
 370     dydx = fTracklet->GetYref(1);
 371     dzdx = fTracklet->GetZref(1);
 372     fTracklet->GetCovRef(covR);
 373     Float_t tilt(fTracklet->GetTilt());
 374     AliTRDcluster *c = NULL;
 375     fTracklet->ResetClusterIter(kFALSE);
 376     while((c = fTracklet->PrevCluster())){
 377       Float_t xc = c->GetX();
 378       Float_t yc = c->GetY();
 379       Float_t zc = c->GetZ();
 380       Float_t dx = x0 - xc;
 381       Float_t yt = y0 - dx*dydx;
 382       Float_t zt = z0 - dx*dzdx;
 383       // calculate residuals using tilt correction
 384 //       yc -= tilt*(zc-zt);
 385 //       dy = yt - yc;
 386
 387       // calculate residuals using correct covariance matrix
 388       cov[0] = c->GetSigmaY2();
 389       cov[1] = c->GetSigmaYZ();
 390       cov[2] = c->GetSigmaZ2();
 391       // do rotation
 392       Double_t sy2(cov[0]), sz2(cov[2]);
 393       Double_t t2 = tilt*tilt;
 394       Double_t correction = 1./(1. + t2);
 395       cov[0] = (sy2+t2*sz2)*correction;
 396       cov[1] = tilt*(sz2 - sy2)*correction;
 397       cov[2] = (t2*sy2+sz2)*correction;
 398       // do inversion
 399       Double_t r00=cov[0]+covR[0], r01=cov[1]+covR[1], r11=cov[2]+covR[2];
 400       Double_t det=r00*r11 - r01*r01;
 401       Double_t tmp=r00; r00=r11/det; r11=tmp/det;
 402       dy = (yc - yt)*TMath::Sqrt(r00);
 403       dz = (zc - zt)*TMath::Sqrt(r11);
 404
 405       ((TH2I*)arr->At(0))->Fill(dydx, dy/*, dz*/);
 406       ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
 407
 408       if(DebugLevel()>=2){
 409         // Get z-position with respect to anode wire
 410         Int_t istk = fGeo->GetStack(c->GetDetector());
 411         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 412         Float_t row0 = pp->GetRow0();
 413         Float_t d  =  row0 - zt + pp->GetAnodeWireOffset();
 414         d -= ((Int_t)(2 * d)) / 2.0;
 415         if (d > 0.25) d  = 0.5 - d;
 416
 417         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 418         fCl->Add(clInfo);
 419         clInfo->SetCluster(c);
 420         Float_t covF[] = {cov[0], cov[1], cov[2]};
 421         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
 422         clInfo->SetResolution(dy);
 423         clInfo->SetAnisochronity(d);
 424         clInfo->SetDriftLength(dx);
 425         clInfo->SetTilt(tilt);
 426         (*DebugStream()) << "ClusterREC"
 427           <<"status="  << status
 428           <<"clInfo.=" << clInfo
 429           << "\n";
 430       }
 431     }
 432   }
 433   return (TH2I*)arr->At(0);
 434 }
 435
 436
 437 //________________________________________________________
 438 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 439 {
 440 // Plot normalized residuals for tracklets to track.
 441 //
 442 // We start from the result that if X=N(|m|, |Cov|)
 443 // BEGIN_LATEX
 444 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
 445 // END_LATEX
 446 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
 447 // reference position.
 448   if(track) fkTrack = track;
 449   if(!fkTrack){
 450     AliDebug(2, "No Track defined.");
 451     return NULL;
 452   }
 453   TObjArray *arr = NULL;
 454   if(!(arr = (TObjArray*)fContainer->At(kTrackTRD ))){
 455     AliWarning("No output container defined.");
 456     return NULL;
 457   }
 458
 459   Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
 460   Float_t x, dx, dy, dz;
 461   AliTRDseedV1 *fTracklet = NULL;
 462   for(Int_t il=AliTRDgeometry::kNlayer; il--;){
 463     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
 464     if(!fTracklet->IsOK()) continue;
 465     x    = fTracklet->GetX();
 466     dx   = fTracklet->GetX0() - x;
 467     // compute dy^2 and dz^2
 468     dy   = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
 469     dz   = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
 470     // compute covariance matrix
 471     fTracklet->GetCovAt(x, cov);
 472     fTracklet->GetCovRef(covR);
 473     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 474 /*  // Correct PULL calculation by considering off
 475     // diagonal elements in the covariance matrix
 476     // compute square root matrix
 477     if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
 478     if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
 479     Double_t y = sqr[0]*dy+sqr[1]*dz;
 480     Double_t z = sqr[1]*dy+sqr[2]*dz;
 481     ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
 482
 483     ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
 484     ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
 485     ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
 486     if(!fTracklet->IsRowCross()) continue;
 487     ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
 488     ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
 489   }
 490
 491
 492   return (TH2I*)arr->At(0);
 493 }
 494
 495
 496 //________________________________________________________
 497 TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
 498 {
 499 // Store resolution/pulls of Kalman before updating with the TRD information
 500 // at the radial position of the first tracklet. The following points are used
 501 // for comparison
 502 //  - the (y,z,snp) of the first TRD tracklet
 503 //  - the (y, z, snp, tgl, pt) of the MC track reference
 504 //
 505 // Additionally the momentum resolution/pulls are calculated for usage in the
 506 // PID calculation.
 507
 508   if(track) fkTrack = track;
 509   if(!fkTrack){
 510     AliDebug(2, "No Track defined.");
 511     return NULL;
 512   }
 513   AliExternalTrackParam *tin = NULL;
 514   if(!(tin = fkTrack->GetTrackLow())){
 515     AliWarning("Track did not entered TRD fiducial volume.");
 516     return NULL;
 517   }
 518   TH1 *h = NULL;
 519
 520   Double_t x = tin->GetX();
 521   AliTRDseedV1 *tracklet = NULL;
 522   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 523     if(!(tracklet = fkTrack->GetTracklet(ily))) continue;
 524     break;
 525   }
 526   if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
 527     AliWarning("Tracklet did not match TRD entrance.");
 528     return NULL;
 529   }
 530   const Int_t kNPAR(5);
 531   Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
 532   Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
 533   Double_t cov[3]; tracklet->GetCovAt(x, cov);
 534
 535   // define sum covariances
 536   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
 537   Double_t *pc = &covR[0], *pp = &parR[0];
 538   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
 539     PAR(ir) = (*pp);
 540     for(Int_t ic = 0; ic<=ir; ic++,pc++){
 541       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
 542     }
 543   }
 544   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
 545   //COV.Print(); PAR.Print();
 546
 547   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
 548   Double_t dy = parR[0] - tracklet->GetY();
 549   TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
 550   ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
 551   ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
 552   if(tracklet->IsRowCross()){
 553     Double_t dz = parR[1] - tracklet->GetZ();
 554     ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
 555     ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
 556   }
 557   Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1));  ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
 558
 559
 560   // register reference histo for mini-task
 561   h = (TH2I*)arr->At(0);
 562
 563   if(DebugLevel()>=1){
 564     (*DebugStream()) << "trackIn"
 565       << "x="       << x
 566       << "P="       << &PAR
 567       << "C="       << &COV
 568       << "\n";
 569
 570     Double_t y = tracklet->GetY();
 571     Double_t z = tracklet->GetZ();
 572     (*DebugStream()) << "trackletIn"
 573       << "y="       << y
 574       << "z="       << z
 575       << "Vy="      << cov[0]
 576       << "Cyz="     << cov[1]
 577       << "Vz="      << cov[2]
 578       << "\n";
 579   }
 580
 581
 582   if(!HasMCdata()) return h;
 583   UChar_t s;
 584   Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
 585   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
 586   Int_t pdg = fkMC->GetPDG(),
 587         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 588         sign(0);
 589   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 590   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 591   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 592
 593   // translate to reference radial position
 594   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
 595   Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
 596   //Fill MC info
 597   TVectorD PARMC(kNPAR);
 598   PARMC[0]=y0; PARMC[1]=z0;
 599   PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
 600   PARMC[4]=1./pt0;
 601
 602 //   TMatrixDSymEigen eigen(COV);
 603 //   TVectorD evals = eigen.GetEigenValues();
 604 //   TMatrixDSym evalsm(kNPAR);
 605 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
 606 //   TMatrixD evecs = eigen.GetEigenVectors();
 607 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
 608
 609   // fill histos
 610   arr = (TObjArray*)fContainer->At(kMCtrackTPC);
 611   // y resolution/pulls
 612   ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
 613   ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
 614   // z resolution/pulls
 615   ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
 616   ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
 617   // phi resolution/snp pulls
 618   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
 619   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
 620   // theta resolution/tgl pulls
 621   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
 622   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
 623   // pt resolution\\1/pt pulls\\p resolution/pull
 624   ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
 625   ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
 626
 627   Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
 628   p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
 629   Float_t sp =
 630     p*p*PAR[4]*PAR[4]*COV(4,4)
 631    +2.*PAR[3]*COV(3,4)/PAR[4]
 632    +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
 633   ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
 634   if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
 635
 636   // fill debug for MC
 637   if(DebugLevel()>=1){
 638     (*DebugStream()) << "trackInMC"
 639       << "P="   << &PARMC
 640       << "\n";
 641   }
 642   return h;
 643 }
 644
 645 //________________________________________________________
 646 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 647 {
 648   //
 649   // Plot MC distributions
 650   //
 651
 652   if(!HasMCdata()){
 653     AliWarning("No MC defined. Results will not be available.");
 654     return NULL;
 655   }
 656   if(track) fkTrack = track;
 657   if(!fkTrack){
 658     AliDebug(2, "No Track defined.");
 659     return NULL;
 660   }
 661   // retriev track characteristics
 662   Int_t pdg = fkMC->GetPDG(),
 663         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 664         sign(0),
 665         det(-1),
 666         label(fkMC->GetLabel());
 667   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 668   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 669   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 670   Bool_t kBarrel = Bool_t(fkESD->GetStatus() & AliESDtrack::kTRDin);
 671
 672   TObjArray *arr(NULL);TH1 *h(NULL);
 673   UChar_t s;
 674   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
 675   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 676   Double_t covR[7]/*, cov[3]*/;
 677
 678   if(DebugLevel()>=1){
 679     TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
 680     fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
 681     (*DebugStream()) << "MCkalman"
 682       << "pdg=" << pdg
 683       << "dx="  << &dX
 684       << "dy="  << &dY
 685       << "dz="  << &dZ
 686       << "pt="  << &Pt
 687       << "dpt=" << &dPt
 688       << "cov=" << &cCOV
 689       << "\n";
 690   }
 691
 692   AliTRDReconstructor rec;
 693   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
 694   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 695     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
 696        !fTracklet->IsOK())*/ continue;
 697
 698     det = fTracklet->GetDetector();
 699     x0  = fTracklet->GetX0();
 700     //radial shift with respect to the MC reference (radial position of the pad plane)
 701     x= fTracklet->GetX();
 702     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
 703     xAnode  = fTracklet->GetX0();
 704
 705     // MC track position at reference radial position
 706     dx  = x0 - x;
 707     if(DebugLevel()>=1){
 708       (*DebugStream()) << "MC"
 709         << "det="     << det
 710         << "pdg="     << pdg
 711         << "sgn="     << sign
 712         << "barrel="  << kBarrel
 713         << "pt="      << pt0
 714         << "x="       << x0
 715         << "y="       << y0
 716         << "z="       << z0
 717         << "dydx="    << dydx0
 718         << "dzdx="    << dzdx0
 719         << "\n";
 720     }
 721     Float_t ymc = y0 - dx*dydx0;
 722     Float_t zmc = z0 - dx*dzdx0;
 723     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
 724
 725     // Kalman position at reference radial position
 726     dx = xAnode - x;
 727     dydx = fTracklet->GetYref(1);
 728     dzdx = fTracklet->GetZref(1);
 729     dzdl = fTracklet->GetTgl();
 730     y  = fTracklet->GetYref(0) - dx*dydx;
 731     dy = y - ymc;
 732     z  = fTracklet->GetZref(0) - dx*dzdx;
 733     dz = z - zmc;
 734     pt = TMath::Abs(fTracklet->GetPt());
 735     fTracklet->GetCovRef(covR);
 736
 737     arr = (TObjArray*)fContainer->At(kMCtrackTRD);
 738     // y resolution/pulls
 739     ((TH2I*)arr->At(0))->Fill(dydx0, dy);
 740     ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
 741     // z resolution/pulls
 742     ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
 743     ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
 744     // phi resolution/ snp pulls
 745     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
 746     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
 747     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
 748     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
 749     // theta resolution/ tgl pulls
 750     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
 751               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
 752     ((TH2I*)arr->At(6))->Fill(dzdl0,
 753     TMath::ATan(dtgl));
 754     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
 755     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
 756     Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
 757              p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
 758     if(kBarrel){
 759       ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., sign*sIdx);
 760       ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
 761       ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., sign*sIdx);
 762     } else {
 763       ((TH3S*)((TObjArray*)arr->At(11))->At(ily))->Fill(pt0, pt/pt0-1., sign*sIdx);
 764       ((TH3S*)((TObjArray*)arr->At(12))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
 765       ((TH3S*)((TObjArray*)arr->At(13))->At(ily))->Fill(p0, p/p0-1., sign*sIdx);
 766     }
 767
 768     // Fill Debug stream for Kalman track
 769     if(DebugLevel()>=1){
 770       (*DebugStream()) << "MCtrack"
 771         << "pt="      << pt
 772         << "x="       << x
 773         << "y="       << y
 774         << "z="       << z
 775         << "dydx="    << dydx
 776         << "dzdx="    << dzdx
 777         << "s2y="     << covR[0]
 778         << "s2z="     << covR[2]
 779         << "\n";
 780     }
 781
 782     // recalculate tracklet based on the MC info
 783     AliTRDseedV1 tt(*fTracklet);
 784     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
 785     tt.SetZref(1, dzdx0);
 786     tt.SetReconstructor(&rec);
 787     tt.Fit(kTRUE, kTRUE);
 788     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
 789     dydx = tt.GetYfit(1);
 790     dx = x0 - x;
 791     ymc = y0 - dx*dydx0;
 792     zmc = z0 - dx*dzdx0;
 793     Bool_t rc = tt.IsRowCross();
 794
 795     // add tracklet residuals for y and dydx
 796     arr = (TObjArray*)fContainer->At(kMCtracklet);
 797     if(!rc){
 798       dy    = y-ymc;
 799
 800       Float_t dphi  = (dydx - dydx0);
 801       dphi /= (1.- dydx*dydx0);
 802
 803       ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
 804       if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
 805       ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
 806     } else {
 807       // add tracklet residuals for z
 808       dz = z-zmc;
 809       ((TH2I*)arr->At(2))->Fill(dzdl0, dz);
 810       if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()));
 811     }
 812
 813     // Fill Debug stream for tracklet
 814     if(DebugLevel()>=1){
 815       Float_t s2y = tt.GetS2Y();
 816       Float_t s2z = tt.GetS2Z();
 817       (*DebugStream()) << "MCtracklet"
 818         << "rc="    << rc
 819         << "x="     << x
 820         << "y="     << y
 821         << "z="     << z
 822         << "dydx="  << dydx
 823         << "s2y="   << s2y
 824         << "s2z="   << s2z
 825         << "\n";
 826     }
 827
 828     Int_t istk = AliTRDgeometry::GetStack(det);
 829     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 830     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 831     Float_t tilt = fTracklet->GetTilt();
 832     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
 833
 834     arr = (TObjArray*)fContainer->At(kMCcluster);
 835     AliTRDcluster *c = NULL;
 836     fTracklet->ResetClusterIter(kFALSE);
 837     while((c = fTracklet->PrevCluster())){
 838       Float_t  q = TMath::Abs(c->GetQ());
 839       x = c->GetX(); y = c->GetY();z = c->GetZ();
 840       dx = x0 - x;
 841       ymc= y0 - dx*dydx0;
 842       zmc= z0 - dx*dzdx0;
 843       dy = (y - tilt*(z-zmc)) - ymc;
 844       // Fill Histograms
 845       if(q>20. && q<250.){
 846         ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
 847         ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2()));
 848       }
 849
 850       // Fill calibration container
 851       Float_t d = zr0 - zmc;
 852       d -= ((Int_t)(2 * d)) / 2.0;
 853       if (d > 0.25) d  = 0.5 - d;
 854       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 855       clInfo->SetCluster(c);
 856       clInfo->SetMC(pdg, label);
 857       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
 858       clInfo->SetResolution(dy);
 859       clInfo->SetAnisochronity(d);
 860       clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
 861       clInfo->SetTilt(tilt);
 862       fMCcl->Add(clInfo);
 863       if(DebugLevel()>=2){
 864         if(!clInfoArr) clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 865         clInfoArr->Add(clInfo);
 866       }
 867     }
 868     // Fill Debug Tree
 869     if(DebugLevel()>=2 && clInfoArr){
 870       (*DebugStream()) << "MCcluster"
 871         <<"clInfo.=" << clInfoArr
 872         << "\n";
 873       clInfoArr->Clear();
 874     }
 875   }
 876   if(clInfoArr) delete clInfoArr;
 877   return h;
 878 }
 879
 880
 881 //________________________________________________________
 882 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
 883 {
 884   //
 885   // Get the reference figures
 886   //
 887
 888   Float_t xy[4] = {0., 0., 0., 0.};
 889   if(!gPad){
 890     AliWarning("Please provide a canvas to draw results.");
 891     return kFALSE;
 892   }
 893   Int_t first(2),     // first particle species to be drawn
 894         nspecies(7); // last particle species to be drawn
 895   TList *l = NULL; TVirtualPad *pad=NULL;
 896   switch(ifig){
 897   case kCharge:
 898     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 899     ((TVirtualPad*)l->At(0))->cd();
 900     ((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0))->Draw("apl");
 901     ((TVirtualPad*)l->At(1))->cd();
 902     ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0))->Draw("apl");
 903     break;
 904   case kCluster:
 905     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 906     xy[0] = -.3; xy[1] = -200.; xy[2] = .3; xy[3] = 1000.;
 907     ((TVirtualPad*)l->At(0))->cd();
 908     if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
 909     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 910     ((TVirtualPad*)l->At(1))->cd();
 911     if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
 912     return kTRUE;
 913   case kTrackTRD :
 914     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 915     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
 916     ((TVirtualPad*)l->At(0))->cd();
 917     if(!GetGraphPlot(&xy[0], kTrackTRD , 0)) break;
 918     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 919     ((TVirtualPad*)l->At(1))->cd();
 920     if(!GetGraphPlot(&xy[0], kTrackTRD , 1)) break;
 921     return kTRUE;
 922   case 3: // kTrackTRD  z
 923     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 924     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
 925     ((TVirtualPad*)l->At(0))->cd();
 926     if(!GetGraphPlot(&xy[0], kTrackTRD , 2)) break;
 927     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 928     ((TVirtualPad*)l->At(1))->cd();
 929     if(!GetGraphPlot(&xy[0], kTrackTRD , 3)) break;
 930     return kTRUE;
 931   case 4: // kTrackTRD  phi
 932     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
 933     if(GetGraphPlot(&xy[0], kTrackTRD , 4)) return kTRUE;
 934     break;
 935   case 5: // kTrackTPC y
 936     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 937     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
 938     pad = ((TVirtualPad*)l->At(0)); pad->cd();
 939     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
 940     if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
 941     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 942     pad=((TVirtualPad*)l->At(1)); pad->cd();
 943     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
 944     if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
 945     return kTRUE;
 946   case 6: // kTrackTPC z
 947     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 948     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
 949     pad = ((TVirtualPad*)l->At(0)); pad->cd();
 950     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
 951     if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
 952     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 953     pad = ((TVirtualPad*)l->At(1)); pad->cd();
 954     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
 955     if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
 956     return kTRUE;
 957   case 7: // kTrackTPC phi
 958     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
 959     if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
 960     break;
 961   case 8: // kMCcluster pt <2 GeV/c
 962     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 963     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
 964     ((TVirtualPad*)l->At(0))->cd();
 965     if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
 966     ((TVirtualPad*)l->At(1))->cd();
 967     if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
 968     return kTRUE;
 969   case 9: // kMCcluster pt > 3 GeV/c
 970     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 971     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
 972     ((TVirtualPad*)l->At(0))->cd();
 973     if(!GetGraphPlot(&xy[0], kMCcluster, 2)) break;
 974     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 975     ((TVirtualPad*)l->At(1))->cd();
 976     if(!GetGraphPlot(&xy[0], kMCcluster, 3)) break;
 977     return kTRUE;
 978   case 10: //kMCtracklet [y] pt < 3 GeV/c
 979     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 980     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
 981     ((TVirtualPad*)l->At(0))->cd();
 982     if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
 983     ((TVirtualPad*)l->At(1))->cd();
 984     if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
 985     return kTRUE;
 986   case 11: //kMCtracklet [y] pt > 3 GeV/c
 987     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 988     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
 989     ((TVirtualPad*)l->At(0))->cd();
 990     if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
 991     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
 992     ((TVirtualPad*)l->At(1))->cd();
 993     if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
 994     return kTRUE;
 995   case 12: //kMCtracklet [z]
 996     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
 997     xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
 998     ((TVirtualPad*)l->At(0))->cd();
 999     if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
1000     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1001     ((TVirtualPad*)l->At(1))->cd();
1002     if(!GetGraphPlot(&xy[0], kMCtracklet, 5)) break;
1003     return kTRUE;
1004   case 13: //kMCtracklet [phi]
1005     xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
1006     if(!GetGraphPlot(&xy[0], kMCtracklet, 6)) break;
1007     return kTRUE;
1008   case 14: //kMCtrackTRD [y]
1009     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1010     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1011     ((TVirtualPad*)l->At(0))->cd();
1012     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 0)) break;
1013     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 3.5;
1014     ((TVirtualPad*)l->At(1))->cd();
1015     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 1)) break;
1016     return kTRUE;
1017   case 15: //kMCtrackTRD [z]
1018     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1019     xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.;
1020     ((TVirtualPad*)l->At(0))->cd();
1021     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break;
1022     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1023     ((TVirtualPad*)l->At(1))->cd();
1024     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break;
1025     return kTRUE;
1026   case 16: //kMCtrackTRD [phi/snp]
1027     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1028     xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.;
1029     ((TVirtualPad*)l->At(0))->cd();
1030     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break;
1031     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1032     ((TVirtualPad*)l->At(1))->cd();
1033     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break;
1034     return kTRUE;
1035   case 17: //kMCtrackTRD [theta/tgl]
1036     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1037     xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
1038     ((TVirtualPad*)l->At(0))->cd();
1039     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break;
1040     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1041     ((TVirtualPad*)l->At(1))->cd();
1042     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break;
1043     return kTRUE;
1044   case 18: //kMCtrackTRD [pt]
1045     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1046     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1047     pad = (TVirtualPad*)l->At(0); pad->cd();
1048     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1049     if(!GetGraphTrack(&xy[0], 8, first, nspecies)) break;
1050     pad->Modified(); pad->Update(); pad->SetLogx();
1051     pad = (TVirtualPad*)l->At(1); pad->cd();
1052     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1053     if(!GetGraphTrack(&xy[0], 8, 55+first, nspecies, kTRUE)) break;
1054     pad->Modified(); pad->Update(); pad->SetLogx();
1055     return kTRUE;
1056   case 19: //kMCtrackTRD [1/pt] pulls
1057     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1058     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1059     pad = (TVirtualPad*)l->At(0); pad->cd();
1060     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1061     if(!GetGraphTrack(&xy[0], 9, first, nspecies)) break;
1062     pad = (TVirtualPad*)l->At(1); pad->cd();
1063     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1064     if(!GetGraphTrack(&xy[0], 9, 55+first, nspecies, kTRUE)) break;
1065     return kTRUE;
1066   case 20: //kMCtrackTRD [p]
1067     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1068     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1069     pad = (TVirtualPad*)l->At(0); pad->cd();
1070     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1071     if(!GetGraphTrack(&xy[0], 10, first, nspecies)) break;
1072     pad->Modified(); pad->Update(); pad->SetLogx();
1073     pad = (TVirtualPad*)l->At(1); pad->cd();
1074     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1075     if(!GetGraphTrack(&xy[0], 10, 55+first, nspecies, kTRUE)) break;
1076     pad->Modified(); pad->Update(); pad->SetLogx();
1077     return kTRUE;
1078   case 21: //kMCtrackTRD - SA [pt]
1079     xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.;
1080     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1081     pad = (TVirtualPad*)l->At(0); pad->cd();
1082     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1083     if(!GetGraphTrack(&xy[0], 11, first, nspecies)) break;
1084     pad->SetLogx();
1085     pad = (TVirtualPad*)l->At(1); pad->cd();
1086     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1087     if(!GetGraphTrack(&xy[0], 11, 55+first, nspecies, kTRUE)) break;
1088     pad->SetLogx();
1089     return kTRUE;
1090   case 22: //kMCtrackTRD - SA [1/pt] pulls
1091     xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.;
1092     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1093     pad = (TVirtualPad*)l->At(0); pad->cd();
1094     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1095     if(!GetGraphTrack(&xy[0], 12, first, nspecies)) break;
1096     pad = (TVirtualPad*)l->At(1); pad->cd();
1097     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1098     if(!GetGraphTrack(&xy[0], 12, 55+first, nspecies, kTRUE)) break;
1099     return kTRUE;
1100   case 23: //kMCtrackTRD - SA [p]
1101     xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5;
1102     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1103     pad = (TVirtualPad*)l->At(0); pad->cd();
1104     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1105     if(!GetGraphTrack(&xy[0], 13, first, nspecies)) break;
1106     pad->SetLogx();
1107     pad = (TVirtualPad*)l->At(1); pad->cd();
1108     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1109     if(!GetGraphTrack(&xy[0], 13, 55+first, nspecies, kTRUE)) break;
1110     pad->SetLogx();
1111     return kTRUE;
1112   case 24: // kMCtrackTPC [y]
1113     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1114     xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
1115     ((TVirtualPad*)l->At(0))->cd();
1116     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
1117     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
1118     ((TVirtualPad*)l->At(1))->cd();
1119     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
1120     return kTRUE;
1121   case 25: // kMCtrackTPC [z]
1122     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1123     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
1124     ((TVirtualPad*)l->At(0))->cd();
1125     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
1126     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1127     ((TVirtualPad*)l->At(1))->cd();
1128     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
1129     return kTRUE;
1130   case 26: // kMCtrackTPC [phi|snp]
1131     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1132     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1133     ((TVirtualPad*)l->At(0))->cd();
1134     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
1135     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1136     ((TVirtualPad*)l->At(1))->cd();
1137     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
1138     return kTRUE;
1139   case 27: // kMCtrackTPC [theta|tgl]
1140     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1141     xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1142     ((TVirtualPad*)l->At(0))->cd();
1143     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
1144     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1145     ((TVirtualPad*)l->At(1))->cd();
1146     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
1147     return kTRUE;
1148   case 28: // kMCtrackTPC [pt]
1149     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1150     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1151     pad=(TVirtualPad*)l->At(0); pad->cd();
1152     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1153     if(!GetGraphTrackTPC(xy, 8, first, nspecies)) break;
1154     pad->SetLogx();
1155     xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =3.;
1156     pad = (TVirtualPad*)l->At(1); pad->cd();
1157     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1158     if(!GetGraphTrackTPC(xy, 9, first, nspecies, kTRUE)) break;
1159     return kTRUE;
1160   case 29: // kMCtrackTPC [p]
1161     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1162     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1163     pad = ((TVirtualPad*)l->At(0));pad->cd();
1164     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1165     if(!GetGraphTrackTPC(xy, 10, first, nspecies)) break;
1166     pad->SetLogx();
1167     xy[0]=0.; xy[1]=-0.5; xy[2]=8.; xy[3] =2.5;
1168     pad = ((TVirtualPad*)l->At(1)); pad->cd();
1169     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1170     if(!GetGraphTrackTPC(xy, 11, first, nspecies, kTRUE)) break;
1171     return kTRUE;
1172   case 30:  // kMCtrackTOF [z]
1173     return kTRUE;
1174   }
1175   AliWarning(Form("Reference plot [%d] missing result", ifig));
1176   return kFALSE;
1177 }
1178
1179 Char_t const *fgParticle[11]={
1180   " p bar", " K -", " #pi -", " #mu -", " e -",
1181   " No PID",
1182   " e +", " #mu +", " #pi +", " K +", " p",
1183 };
1184 const Color_t fgColorS[11]={
1185 kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
1186 kGray,
1187 kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1188 };
1189 const Color_t fgColorM[11]={
1190 kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1191 kBlack,
1192 kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
1193 };
1194 const Marker_t fgMarker[11]={
1195 30, 30, 26, 25, 24,
1196 28,
1197 20, 21, 22, 29, 29
1198 };
1199 //________________________________________________________
1200 Bool_t AliTRDresolution::PostProcess()
1201 {
1202   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1203   if (!fContainer) {
1204     AliError("ERROR: list not available");
1205     return kFALSE;
1206   }
1207   Int_t nc(0);
1208   TGraph *gm= NULL, *gs= NULL;
1209   if(!fGraphS && !fGraphM){
1210     TObjArray *aM(NULL), *aS(NULL);
1211     Int_t n = fContainer->GetEntriesFast();
1212     fGraphS = new TObjArray(n); fGraphS->SetOwner();
1213     fGraphM = new TObjArray(n); fGraphM->SetOwner();
1214     for(Int_t ig=0; ig<n; ig++){
1215       fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
1216       fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
1217
1218       for(Int_t ic=0; ic<fgNproj[ig]; ic++){
1219         //printf("g[%d] c[%d] n[%d]\n", ig, ic, fgNcomp[nc]);
1220         if(fgNcomp[nc]>1){
1221           TObjArray *agS(NULL), *agM(NULL);
1222           aS->AddAt(agS = new TObjArray(fgNcomp[nc]), ic);
1223           aM->AddAt(agM = new TObjArray(fgNcomp[nc]), ic);
1224           for(Int_t is=fgNcomp[nc]; is--;){
1225             agS->AddAt(gs = new TGraphErrors(), is);
1226             Int_t is0(is%11);
1227             gs->SetMarkerStyle(fgMarker[is0]);
1228             gs->SetMarkerColor(fgColorS[is0]);
1229             gs->SetLineColor(fgColorS[is0]);
1230             gs->SetNameTitle(Form("s_%d%02d%02d", ig, ic, is), fgParticle[is0]);
1231
1232             agM->AddAt(gm = new TGraphErrors(), is);
1233             gm->SetMarkerStyle(fgMarker[is0]);
1234             gm->SetMarkerColor(fgColorM[is0]);
1235             gm->SetLineColor(fgColorM[is0]);
1236             gm->SetNameTitle(Form("m_%d%02d%02d", ig, ic, is), fgParticle[is0]);
1237           }
1238         } else {
1239           aS->AddAt(gs = new TGraphErrors(), ic);
1240           gs->SetMarkerStyle(23);
1241           gs->SetMarkerColor(kRed);
1242           gs->SetLineColor(kRed);
1243           gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
1244
1245           aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
1246           gm->SetLineColor(kBlack);
1247           gm->SetMarkerStyle(7);
1248           gm->SetMarkerColor(kBlack);
1249           gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
1250         }
1251         nc+=fgNcomp[ic];
1252       }
1253     }
1254   }
1255
1256 /*  printf("\n\n\n"); fGraphS->ls();
1257   printf("\n\n\n"); fGraphM->ls();*/
1258
1259
1260   // DEFINE MODELS
1261   // simple gauss
1262   TF1 fg("fGauss", "gaus", -.5, .5);
1263   // Landau for charge resolution
1264   TF1 fch("fClCh", "landau", 0., 1000.);
1265   // Landau for e+- pt resolution
1266   TF1 fpt("fPt", "landau", -0.1, 0.2);
1267
1268   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1269   // Charge resolution
1270   //Process3DL(kCharge, 0, &fl);
1271   // Clusters residuals
1272   Process2D(kCluster, 0, &fg, 1.e4);
1273   Process2D(kCluster, 1, &fg);
1274   fNRefFigures = 1;
1275   // Tracklet residual/pulls
1276   Process2D(kTrackTRD , 0, &fg, 1.e4);
1277   Process2D(kTrackTRD , 1, &fg);
1278   Process2D(kTrackTRD , 2, &fg, 1.e4);
1279   Process2D(kTrackTRD , 3, &fg);
1280   Process2D(kTrackTRD , 4, &fg, 1.e3);
1281   fNRefFigures = 4;
1282   // TPC track residual/pulls
1283   Process2D(kTrackTPC, 0, &fg, 1.e4);
1284   Process2D(kTrackTPC, 1, &fg);
1285   Process2D(kTrackTPC, 2, &fg, 1.e4);
1286   Process2D(kTrackTPC, 3, &fg);
1287   Process2D(kTrackTPC, 4, &fg, 1.e3);
1288   fNRefFigures = 7;
1289
1290   if(!HasMCdata()) return kTRUE;
1291
1292
1293   //PROCESS MC RESIDUAL DISTRIBUTIONS
1294
1295   // CLUSTER Y RESOLUTION/PULLS
1296   Process3Drange(kMCcluster, 0, 0, &fg, 1.e4, 1, 1);
1297   Process3Drange(kMCcluster, 0, 1, &fg, 1.e4, 2, 2);
1298   Process3Drange(kMCcluster, 0, 2, &fg, 1.e4, 3, 12);
1299   Process2D(kMCcluster, 1, &fg, 1., 3);
1300   fNRefFigures = 10;
1301
1302   // TRACKLET RESOLUTION/PULLS
1303   Process3Drange(kMCtracklet, 0, 0, &fg, 1.e4, 1, 1); // y
1304   Process3Drange(kMCtracklet, 0, 1, &fg, 1.e4, 2, 2); // y
1305   Process3Drange(kMCtracklet, 0, 2, &fg, 1.e4, 3, 12); // y
1306   Process2D(kMCtracklet, 1, &fg, 1., 3);       // y pulls
1307   Process2D(kMCtracklet, 2, &fg, 1.e4, 4); // z
1308   Process2D(kMCtracklet, 3, &fg, 1., 5);       // z pulls
1309   Process2D(kMCtracklet, 4, &fg, 1.e3, 6); // phi
1310   fNRefFigures = 14;
1311
1312   // TRACK RESOLUTION/PULLS
1313   Process2D(kMCtrackTRD, 0, &fg, 1.e4);   // y
1314   Process2D(kMCtrackTRD, 1, &fg);         // y PULL
1315   Process2D(kMCtrackTRD, 2, &fg, 1.e4);   // z
1316   Process2D(kMCtrackTRD, 3, &fg);         // z PULL
1317   Process2D(kMCtrackTRD, 4, &fg, 1.e3);   // phi
1318   Process2D(kMCtrackTRD, 5, &fg);         // snp PULL
1319   Process2D(kMCtrackTRD, 6, &fg, 1.e3);   // theta
1320   Process2D(kMCtrackTRD, 7, &fg);         // tgl PULL
1321   Process4D(kMCtrackTRD, 8, &fg, 1.e2);   // pt resolution
1322   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 4);// pt resolution e1- @ L0
1323   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 6);// pt resolution e1+ @ L0
1324   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 55+4);// pt resolution e1- @ L5
1325   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 55+6);// pt resolution e1+ @ L5
1326   Process4D(kMCtrackTRD, 9, &fg);         // 1/pt pulls
1327   Process4D(kMCtrackTRD, 10, &fg, 1.e2);  // p resolution
1328   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 4);// p resolution e1- @ L0
1329   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 6);// p resolution e1+ @ L0
1330   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 55+4);// p resolution e1- @ L5
1331   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 55+6);// p resolution e1+ @ L5
1332   Process4D(kMCtrackTRD, 11, &fg, 1.e2);   // pt resolution stand alone
1333   Process4D(kMCtrackTRD, 12, &fg);         // 1/pt pulls  stand alone
1334   Process4D(kMCtrackTRD, 13, &fg, 1.e2);  // p resolution  stand alone
1335   fNRefFigures = 24;
1336
1337   // TRACK TPC RESOLUTION/PULLS
1338   Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
1339   Process2D(kMCtrackTPC, 1, &fg);      // y pulls
1340   Process2D(kMCtrackTPC, 2, &fg, 1.e4);// z resolution
1341   Process2D(kMCtrackTPC, 3, &fg);      // z pulls
1342   Process2D(kMCtrackTPC, 4, &fg, 1.e3);// phi resolution
1343   Process2D(kMCtrackTPC, 5, &fg);      // snp pulls
1344   Process2D(kMCtrackTPC, 6, &fg, 1.e3);// theta resolution
1345   Process2D(kMCtrackTPC, 7, &fg);      // tgl pulls
1346   Process3D(kMCtrackTPC, 8, &fg, 1.e2);// pt resolution
1347   Process3D(kMCtrackTPC, 9, &fg);      // 1/pt pulls
1348   Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
1349   Process3D(kMCtrackTPC, 11, &fg);      // p pulls
1350   fNRefFigures = 30;
1351
1352   // TRACK HMPID RESOLUTION/PULLS
1353   Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
1354   Process2D(kMCtrackTOF, 1, &fg);       // z towards TOF
1355   fNRefFigures = 31;
1356
1357   return kTRUE;
1358 }
1359
1360
1361 //________________________________________________________
1362 void AliTRDresolution::Terminate(Option_t *opt)
1363 {
1364   AliTRDrecoTask::Terminate(opt);
1365   if(HasPostProcess()) PostProcess();
1366 }
1367
1368 //________________________________________________________
1369 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1370 {
1371 // Helper function to avoid duplication of code
1372 // Make first guesses on the fit parameters
1373
1374   // find the intial parameters of the fit !! (thanks George)
1375   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1376   Double_t sum = 0.;
1377   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1378   f->SetParLimits(0, 0., 3.*sum);
1379   f->SetParameter(0, .9*sum);
1380   Double_t rms = h->GetRMS();
1381   f->SetParLimits(1, -rms, rms);
1382   f->SetParameter(1, h->GetMean());
1383
1384   f->SetParLimits(2, 0., 2.*rms);
1385   f->SetParameter(2, rms);
1386   if(f->GetNpar() <= 4) return;
1387
1388   f->SetParLimits(3, 0., sum);
1389   f->SetParameter(3, .1*sum);
1390
1391   f->SetParLimits(4, -.3, .3);
1392   f->SetParameter(4, 0.);
1393
1394   f->SetParLimits(5, 0., 1.e2);
1395   f->SetParameter(5, 2.e-1);
1396 }
1397
1398 //________________________________________________________
1399 TObjArray* AliTRDresolution::Histos()
1400 {
1401   //
1402   // Define histograms
1403   //
1404
1405   if(fContainer) return fContainer;
1406
1407   fContainer  = new TObjArray(kNviews);
1408   //fContainer->SetOwner(kTRUE);
1409   TH1 *h(NULL);
1410   TObjArray *arr(NULL);
1411
1412   // binnings for plots containing momentum or pt
1413   const Int_t kNpt(14), kNphi(48), kNdy(60);
1414   Float_t Phi=-.48, Dy=-.3, Pt=0.1;
1415   Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
1416   for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
1417   for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
1418   for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
1419
1420   // cluster to track residuals/pulls
1421   fContainer->AddAt(arr = new TObjArray(fgNhistos[kCharge]), kCharge);
1422   arr->SetName("Charge");
1423   if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
1424     h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
1425     h->GetXaxis()->SetTitle("dx/dx_{0}");
1426     h->GetYaxis()->SetTitle("x_{d} [cm]");
1427     h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
1428   } else h->Reset();
1429   arr->AddAt(h, 0);
1430
1431   // cluster to track residuals/pulls
1432   fContainer->AddAt(arr = new TObjArray(fgNhistos[kCluster]), kCluster);
1433   arr->SetName("Cl");
1434   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
1435     h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
1436     h->GetXaxis()->SetTitle("tg(#phi)");
1437     h->GetYaxis()->SetTitle("#Delta y [cm]");
1438     h->GetZaxis()->SetTitle("entries");
1439   } else h->Reset();
1440   arr->AddAt(h, 0);
1441   if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
1442     h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1443     h->GetXaxis()->SetTitle("tg(#phi)");
1444     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1445     h->GetZaxis()->SetTitle("entries");
1446   } else h->Reset();
1447   arr->AddAt(h, 1);
1448
1449   // tracklet to track residuals/pulls in y direction
1450   fContainer->AddAt(arr = new TObjArray(fgNhistos[kTrackTRD ]), kTrackTRD );
1451   arr->SetName("Trklt");
1452   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
1453     h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
1454     h->GetXaxis()->SetTitle("#tg(#phi)");
1455     h->GetYaxis()->SetTitle("#Delta y [cm]");
1456     h->GetZaxis()->SetTitle("entries");
1457   } else h->Reset();
1458   arr->AddAt(h, 0);
1459   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
1460     h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1461     h->GetXaxis()->SetTitle("#tg(#phi)");
1462     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1463     h->GetZaxis()->SetTitle("entries");
1464   } else h->Reset();
1465   arr->AddAt(h, 1);
1466   // tracklet to track residuals/pulls in z direction
1467   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
1468     h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1469     h->GetXaxis()->SetTitle("#tg(#theta)");
1470     h->GetYaxis()->SetTitle("#Delta z [cm]");
1471     h->GetZaxis()->SetTitle("entries");
1472   } else h->Reset();
1473   arr->AddAt(h, 2);
1474   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
1475     h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1476     h->GetXaxis()->SetTitle("#tg(#theta)");
1477     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1478     h->GetZaxis()->SetTitle("entries");
1479   } else h->Reset();
1480   arr->AddAt(h, 3);
1481   // tracklet to track phi residuals
1482   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
1483     h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1484     h->GetXaxis()->SetTitle("tg(#phi)");
1485     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1486     h->GetZaxis()->SetTitle("entries");
1487   } else h->Reset();
1488   arr->AddAt(h, 4);
1489
1490
1491   // tracklet to TPC track residuals/pulls in y direction
1492   fContainer->AddAt(arr = new TObjArray(fgNhistos[kTrackTPC]), kTrackTPC);
1493   arr->SetName("TrkTPC");
1494   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
1495     h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
1496     h->GetXaxis()->SetTitle("#tg(#phi)");
1497     h->GetYaxis()->SetTitle("#Delta y [cm]");
1498     h->GetZaxis()->SetTitle("entries");
1499   } else h->Reset();
1500   arr->AddAt(h, 0);
1501   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
1502     h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1503     h->GetXaxis()->SetTitle("#tg(#phi)");
1504     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1505     h->GetZaxis()->SetTitle("entries");
1506   } else h->Reset();
1507   arr->AddAt(h, 1);
1508   // tracklet to TPC track residuals/pulls in z direction
1509   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
1510     h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1511     h->GetXaxis()->SetTitle("#tg(#theta)");
1512     h->GetYaxis()->SetTitle("#Delta z [cm]");
1513     h->GetZaxis()->SetTitle("entries");
1514   } else h->Reset();
1515   arr->AddAt(h, 2);
1516   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
1517     h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1518     h->GetXaxis()->SetTitle("#tg(#theta)");
1519     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1520     h->GetZaxis()->SetTitle("entries");
1521   } else h->Reset();
1522   arr->AddAt(h, 3);
1523   // tracklet to TPC track phi residuals
1524   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
1525     h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1526     h->GetXaxis()->SetTitle("tg(#phi)");
1527     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1528     h->GetZaxis()->SetTitle("entries");
1529   } else h->Reset();
1530   arr->AddAt(h, 4);
1531
1532
1533   // Resolution histos
1534   if(!HasMCdata()) return fContainer;
1535
1536   // cluster y resolution [0]
1537   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCcluster]), kMCcluster);
1538   arr->SetName("McCl");
1539   if(!(h = (TH3S*)gROOT->FindObject("hMcCl"))){
1540     h = new TH3S("hMcCl",
1541     "Cluster Resolution;tg(#phi);#Delta y [cm];p_{t} [GeV/c]",
1542     kNphi, binsPhi, kNdy, binsDy, kNpt, binsPt);
1543   } else h->Reset();
1544   arr->AddAt(h, 0);
1545   if(!(h = (TH2I*)gROOT->FindObject("hMcClPull"))){
1546     h = new TH2I("hMcClPull", "Cluster Pulls", 48, -.48, .48, 100, -4.5, 4.5);
1547     h->GetXaxis()->SetTitle("tg(#phi)");
1548     h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
1549     h->GetZaxis()->SetTitle("p_{t} [GeV/c]");
1550   } else h->Reset();
1551   arr->AddAt(h, 1);
1552
1553
1554   // TRACKLET RESOLUTION
1555   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtracklet]), kMCtracklet);
1556   arr->SetName("McTrklt");
1557   // tracklet y resolution
1558   if(!(h = (TH3S*)gROOT->FindObject("hMcTrkltY"))){
1559     h = new TH3S("hMcTrkltY",
1560     "Tracklet Y Resolution;tg(#phi);#Delta y [cm];p_{t} [GeV/c]",
1561     kNphi, binsPhi, kNdy, binsDy, kNpt, binsPt);
1562   } else h->Reset();
1563   arr->AddAt(h, 0);
1564   // tracklet y pulls
1565   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltYPull"))){
1566     h = new TH2I("hMcTrkltYPull", "Tracklet Pulls (Y)", 48, -.48, .48, 100, -4.5, 4.5);
1567     h->GetXaxis()->SetTitle("tg(#phi)");
1568     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1569     h->GetZaxis()->SetTitle("entries");
1570   } else h->Reset();
1571   arr->AddAt(h, 1);
1572   // tracklet z resolution
1573   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZ"))){
1574     h = new TH2I("hMcTrkltZ", "Tracklet Resolution (Z)", 100, -1., 1., 100, -1., 1.);
1575     h->GetXaxis()->SetTitle("tg(#theta)");
1576     h->GetYaxis()->SetTitle("#Delta z [cm]");
1577     h->GetZaxis()->SetTitle("entries");
1578   } else h->Reset();
1579   arr->AddAt(h, 2);
1580   // tracklet z pulls
1581   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZPull"))){
1582     h = new TH2I("hMcTrkltZPull", "Tracklet Pulls (Z)", 100, -1., 1., 100, -3.5, 3.5);
1583     h->GetXaxis()->SetTitle("tg(#theta)");
1584     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1585     h->GetZaxis()->SetTitle("entries");
1586   } else h->Reset();
1587   arr->AddAt(h, 3);
1588   // tracklet phi resolution
1589   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltPhi"))){
1590     h = new TH2I("hMcTrkltPhi", "Tracklet Resolution (#Phi)", 48, -.48, .48, 100, -.15, .15);
1591     h->GetXaxis()->SetTitle("tg(#phi)");
1592     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1593     h->GetZaxis()->SetTitle("entries");
1594   } else h->Reset();
1595   arr->AddAt(h, 4);
1596
1597
1598   // KALMAN TRACK RESOLUTION
1599   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTRD]), kMCtrackTRD);
1600   arr->SetName("McTrkTRD");
1601   // Kalman track y resolution
1602   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkY"))){
1603     h = new TH2I("hMcTrkY", "Track Y Resolution", 48, -.48, .48, 100, -.2, .2);
1604     h->GetXaxis()->SetTitle("tg(#phi)");
1605     h->GetYaxis()->SetTitle("#Delta y [cm]");
1606     h->GetZaxis()->SetTitle("entries");
1607   } else h->Reset();
1608   arr->AddAt(h, 0);
1609   // Kalman track y pulls
1610   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkYPull"))){
1611     h = new TH2I("hMcTrkYPull", "Track Y Pulls", 48, -.48, .48, 100, -4., 4.);
1612     h->GetXaxis()->SetTitle("tg(#phi)");
1613     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1614     h->GetZaxis()->SetTitle("entries");
1615   } else h->Reset();
1616   arr->AddAt(h, 1);
1617   // Kalman track Z
1618   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZ"))){
1619     h = new TH2I("hMcTrkZ", "Track Z Resolution", 100, -1., 1., 100, -1., 1.);
1620     h->GetXaxis()->SetTitle("tg(#theta)");
1621     h->GetYaxis()->SetTitle("#Delta z [cm]");
1622     h->GetZaxis()->SetTitle("entries");
1623   } else h->Reset();
1624   arr->AddAt(h, 2);
1625   // Kalman track Z pulls
1626   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZPull"))){
1627     h = new TH2I("hMcTrkZPull", "Track Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1628     h->GetXaxis()->SetTitle("tg(#theta)");
1629     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1630     h->GetZaxis()->SetTitle("entries");
1631   } else h->Reset();
1632   arr->AddAt(h, 3);
1633   // Kalman track SNP
1634   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
1635     h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1636     h->GetXaxis()->SetTitle("tg(#phi)");
1637     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1638     h->GetZaxis()->SetTitle("entries");
1639   } else h->Reset();
1640   arr->AddAt(h, 4);
1641   // Kalman track SNP pulls
1642   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNPPull"))){
1643     h = new TH2I("hMcTrkSNPPull", "Track SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1644     h->GetXaxis()->SetTitle("tg(#phi)");
1645     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1646     h->GetZaxis()->SetTitle("entries");
1647   } else h->Reset();
1648   arr->AddAt(h, 5);
1649   // Kalman track TGL
1650   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
1651     h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1652     h->GetXaxis()->SetTitle("tg(#theta)");
1653     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1654     h->GetZaxis()->SetTitle("entries");
1655   } else h->Reset();
1656   arr->AddAt(h, 6);
1657   // Kalman track TGL pulls
1658   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGLPull"))){
1659     h = new TH2I("hMcTrkTGLPull", "Track TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1660     h->GetXaxis()->SetTitle("tg(#theta)");
1661     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1662     h->GetZaxis()->SetTitle("entries");
1663   } else h->Reset();
1664   arr->AddAt(h, 7);
1665
1666   const Int_t kNdpt(150);
1667   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1668   Float_t DPt=-.1, Spc=-5.5;
1669   Float_t binsSpc[kNspc+1], binsDPt[kNdpt+1];
1670   for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
1671   for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
1672   TObjArray *arr2 = NULL; TH3S* h3=NULL;
1673   // Kalman track Pt resolution
1674   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 8);
1675   arr2->SetName("Pt Resolution");
1676   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1677     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPt%d", il)))){
1678       h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1679     } else h3->Reset();
1680     arr2->AddAt(h3, il);
1681   }
1682   // Kalman track Pt pulls
1683   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 9);
1684   arr2->SetName("1/Pt Pulls");
1685   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1686     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPtPulls%d", il)))){
1687       h3 = new TH3S(Form("hMcTrkPtPulls%d", il),
1688       "Track 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES",
1689       kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1690     } else h3->Reset();
1691     arr2->AddAt(h3, il);
1692   }
1693   // Kalman track P resolution
1694   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10);
1695   arr2->SetName("P Resolution");
1696   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1697     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkP%d", il)))){
1698       h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1699     } else h3->Reset();
1700     arr2->AddAt(h3, il);
1701   }
1702   // TRD stand-alone track Pt resolution
1703   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 11);
1704   arr2->SetName("Pt Resolution [SA]");
1705   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1706     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkPt%d", il)))){
1707       h3 = new TH3S(Form("hMcSATrkPt%d", il),
1708       "Track Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES",
1709       kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1710     } else h3->Reset();
1711     arr2->AddAt(h3, il);
1712   }
1713   // TRD stand-alone track Pt pulls
1714   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 12);
1715   arr2->SetName("1/Pt Pulls [SA]");
1716   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1717     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkPtPulls%d", il)))){
1718       h3 = new TH3S(Form("hMcSATrkPtPulls%d", il),
1719       "Track 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES",
1720       kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1721     } else h3->Reset();
1722     arr2->AddAt(h3, il);
1723   }
1724   // TRD stand-alone track P resolution
1725   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 13);
1726   arr2->SetName("P Resolution [SA]");
1727   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1728     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkP%d", il)))){
1729       h3 = new TH3S(Form("hMcSATrkP%d", il),
1730       "Track P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES",
1731       kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1732     } else h3->Reset();
1733     arr2->AddAt(h3, il);
1734   }
1735
1736   // TPC TRACK RESOLUTION
1737   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTPC]), kMCtrackTPC);
1738   arr->SetName("McTrkTPC");
1739   // Kalman track Y
1740   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCY"))){
1741     h = new TH2I("hMcTrkTPCY", "Track[TPC] Y Resolution", 60, -.3, .3, 100, -.5, .5);
1742     h->GetXaxis()->SetTitle("tg(#phi)");
1743     h->GetYaxis()->SetTitle("#Delta y [cm]");
1744     h->GetZaxis()->SetTitle("entries");
1745   } else h->Reset();
1746   arr->AddAt(h, 0);
1747   // Kalman track Y pulls
1748   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCYPull"))){
1749     h = new TH2I("hMcTrkTPCYPull", "Track[TPC] Y Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1750     h->GetXaxis()->SetTitle("tg(#phi)");
1751     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1752     h->GetZaxis()->SetTitle("entries");
1753   } else h->Reset();
1754   arr->AddAt(h, 1);
1755   // Kalman track Z
1756   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZ"))){
1757     h = new TH2I("hMcTrkTPCZ", "Track[TPC] Z Resolution", 100, -1., 1., 100, -1., 1.);
1758     h->GetXaxis()->SetTitle("tg(#theta)");
1759     h->GetYaxis()->SetTitle("#Delta z [cm]");
1760     h->GetZaxis()->SetTitle("entries");
1761   } else h->Reset();
1762   arr->AddAt(h, 2);
1763   // Kalman track Z pulls
1764   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZPull"))){
1765     h = new TH2I("hMcTrkTPCZPull", "Track[TPC] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1766     h->GetXaxis()->SetTitle("tg(#theta)");
1767     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1768     h->GetZaxis()->SetTitle("entries");
1769   } else h->Reset();
1770   arr->AddAt(h, 3);
1771   // Kalman track SNP
1772   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){
1773     h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1774     h->GetXaxis()->SetTitle("tg(#phi)");
1775     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1776     h->GetZaxis()->SetTitle("entries");
1777   } else h->Reset();
1778   arr->AddAt(h, 4);
1779   // Kalman track SNP pulls
1780   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNPPull"))){
1781     h = new TH2I("hMcTrkTPCSNPPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1782     h->GetXaxis()->SetTitle("tg(#phi)");
1783     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1784     h->GetZaxis()->SetTitle("entries");
1785   } else h->Reset();
1786   arr->AddAt(h, 5);
1787   // Kalman track TGL
1788   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){
1789     h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1790     h->GetXaxis()->SetTitle("tg(#theta)");
1791     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1792     h->GetZaxis()->SetTitle("entries");
1793   } else h->Reset();
1794   arr->AddAt(h, 6);
1795   // Kalman track TGL pulls
1796   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGLPull"))){
1797     h = new TH2I("hMcTrkTPCTGLPull", "Track[TPC] TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1798     h->GetXaxis()->SetTitle("tg(#theta)");
1799     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1800     h->GetZaxis()->SetTitle("entries");
1801   } else h->Reset();
1802   arr->AddAt(h, 7);
1803   // Kalman track Pt resolution
1804   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
1805     h3 = new TH3S("hMcTrkTPCPt",
1806     "TRDin Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES",
1807     kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1808   } else h3->Reset();
1809   arr->AddAt(h3, 8);
1810   // Kalman track Pt pulls
1811   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPtPulls"))){
1812     h3 = new TH3S("hMcTrkTPCPtPulls",
1813     "Track[TPC] 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES",
1814     kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1815   } else h3->Reset();
1816   arr->AddAt(h3, 9);
1817   // Kalman track P resolution
1818   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){
1819     h3 = new TH3S("hMcTrkTPCP",
1820     "TRDin P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES",
1821     kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1822   } else h3->Reset();
1823   arr->AddAt(h3, 10);
1824   // Kalman track P pulls
1825   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPPulls"))){
1826     h3 = new TH3S("hMcTrkTPCPPulls",
1827     "TRDin P Pulls;p^{MC} [GeV/c];#Deltap/#sigma_{p};SPECIES",
1828     kNpt, 0., 12., 100, -5., 5., kNspc, -5.5, 5.5);
1829   } else h3->Reset();
1830   arr->AddAt(h3, 11);
1831
1832
1833
1834   // Kalman track Z resolution [TOF]
1835   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTOF]), kMCtrackTOF);
1836   arr->SetName("McTrkTOF");
1837   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZ"))){
1838     h = new TH2I("hMcTrkTOFZ", "Track[TOF] Z Resolution", 100, -1., 1., 100, -1., 1.);
1839     h->GetXaxis()->SetTitle("tg(#theta)");
1840     h->GetYaxis()->SetTitle("#Delta z [cm]");
1841     h->GetZaxis()->SetTitle("entries");
1842   } else h->Reset();
1843   arr->AddAt(h, 0);
1844   // Kalman track Z pulls
1845   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZPull"))){
1846     h = new TH2I("hMcTrkTOFZPull", "Track[TOF] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1847     h->GetXaxis()->SetTitle("tg(#theta)");
1848     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1849     h->GetZaxis()->SetTitle("entries");
1850   } else h->Reset();
1851   arr->AddAt(h, 1);
1852
1853   return fContainer;
1854 }
1855
1856 //________________________________________________________
1857 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1858 {
1859   //
1860   // Do the processing
1861   //
1862
1863   Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
1864   Int_t n = 0;
1865   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1866   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1867
1868   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1869     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1870     TH1D *h = h2->ProjectionY(pn, ibin, ibin);
1871     if(h->GetEntries()<100) continue;
1872     //AdjustF1(h, f);
1873
1874     h->Fit(f, "QN");
1875
1876     Int_t ip = g[0]->GetN();
1877     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1878     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1879     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1880     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1881
1882 /*
1883     g[0]->SetPoint(ip, x, k*h->GetMean());
1884     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1885     g[1]->SetPoint(ip, x, k*h->GetRMS());
1886     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1887   }
1888   fIdxPlot++;
1889   return kTRUE;
1890 }
1891
1892 //________________________________________________________
1893 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
1894 {
1895   //
1896   // Do the processing
1897   //
1898
1899   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1900
1901   // retrive containers
1902   TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1903   if(!h2) return kFALSE;
1904
1905   TGraphErrors *g[2];
1906   if(gidx<0) gidx=idx;
1907   if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1908
1909   if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1910
1911   return Process(h2, f, k, g);
1912 }
1913
1914 //________________________________________________________
1915 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1916 {
1917   //
1918   // Do the processing
1919   //
1920
1921   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1922
1923   // retrive containers
1924   TH3S *h3 = idx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1925   if(!h3) return kFALSE;
1926
1927   TObjArray *gm, *gs;
1928   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1929   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1930   TGraphErrors *g[2];
1931
1932   TAxis *az = h3->GetZaxis();
1933   for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1934     if(!(g[0] = (TGraphErrors*)gm->At(iz-1))) return kFALSE;
1935     if(!(g[1] = (TGraphErrors*)gs->At(iz-1))) return kFALSE;
1936     az->SetRange(iz, iz);
1937     if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1938   }
1939
1940   return kTRUE;
1941 }
1942
1943 //________________________________________________________
1944 Bool_t AliTRDresolution::Process3Drange(ETRDresolutionPlot plot, Int_t hidx, Int_t gidx, TF1 *f, Float_t k, Int_t zbin0, Int_t zbin1)
1945 {
1946   //
1947   // Do the processing
1948   //
1949
1950   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1951
1952   // retrive containers
1953   TH3S *h3 = hidx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(hidx);
1954   if(!h3) return kFALSE;
1955
1956   TGraphErrors *g[2];
1957   if(!(g[0] = (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1958   if(!(g[1] = (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1959
1960   TAxis *az = h3->GetZaxis();
1961   az->SetRange(zbin0, zbin1);
1962   if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1963   return kTRUE;
1964 }
1965
1966 //________________________________________________________
1967 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1968 {
1969   //
1970   // Do the processing
1971   //
1972
1973   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1974
1975   // retrive containers
1976   TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
1977   if(!h3) return kFALSE;
1978
1979   TGraphAsymmErrors *gm;
1980   TGraphErrors *gs;
1981   if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
1982   if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
1983
1984   Float_t x, r, mpv, xM, xm;
1985   TAxis *ay = h3->GetYaxis();
1986   for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
1987     ay->SetRange(iy, iy);
1988     x = ay->GetBinCenter(iy);
1989     TH2F *h2=(TH2F*)h3->Project3D("zx");
1990     TAxis *ax = h2->GetXaxis();
1991     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
1992       r = ax->GetBinCenter(ix);
1993       TH1D *h1 = h2->ProjectionY("py", ix, ix);
1994       if(h1->Integral()<50) continue;
1995       h1->Fit(f, "QN");
1996
1997       GetLandauMpvFwhm(f, mpv, xm, xM);
1998       Int_t ip = gm->GetN();
1999       gm->SetPoint(ip, x, k*mpv);
2000       gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2001       gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2002       gs->SetPointError(ip, 0., 0.);
2003     }
2004   }
2005
2006   return kTRUE;
2007 }
2008
2009 //________________________________________________________
2010 Bool_t AliTRDresolution::Process4D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t n)
2011 {
2012   //
2013   // Do the processing
2014   //
2015
2016   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2017   //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2018
2019   // retrive containers
2020   TObjArray *arr = (TObjArray*)((TObjArray*)(fContainer->At(plot)))->At(idx);
2021   if(!arr) return kFALSE;
2022
2023   TObjArray *gm, *gs;
2024   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2025   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2026
2027   TGraphErrors *g[2];
2028
2029   TH3S *h3(NULL);
2030   for(Int_t ix=0, in=0; ix<arr->GetEntriesFast(); ix++){
2031     if(!(h3 = (TH3S*)arr->At(ix))) return kFALSE;
2032     TAxis *az = h3->GetZaxis();
2033     //printf("  process ix[%d] bins[%d] in[%d]\n", ix, az->GetNbins(), in);
2034     for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
2035       if(n>=0 && n!=in) continue;
2036       if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2037       if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2038       //printf("    g0[%s] g1[%s]\n", g[0]->GetName(), g[1]->GetName());
2039       az->SetRange(iz, iz);
2040       if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2041     }
2042   }
2043
2044   return kTRUE;
2045 }
2046
2047 //________________________________________________________
2048 Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
2049 {
2050   //
2051   // Get the graphs
2052   //
2053
2054   if(!fGraphS || !fGraphM) return kFALSE;
2055
2056   //printf("plotting task[%d] gidx[%d]\n", ip, idx);
2057   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2058   if(!gm) return kFALSE;
2059   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2060   if(!gs) return kFALSE;
2061   //printf("gs[%s] gm[%s]\n", gs->GetName(), gm->GetName());
2062   gs->Draw("apl"); gm->Draw("pl");
2063   //return kTRUE;
2064   // titles look up
2065   Int_t nref = 0;
2066   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
2067   UChar_t jdx = idx<0?0:idx;
2068   for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2069   const Char_t **at = fgAxTitle[nref];
2070
2071   Int_t n(0);
2072   if((n=gm->GetN())) {
2073     PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2074     PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2075   }
2076
2077   if((n=gs->GetN())){
2078     gs->Sort(&TGraph::CompareY);
2079     PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2080     PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2081     gs->Sort(&TGraph::CompareX);
2082   }
2083
2084   // axis range
2085   TAxis *ax(NULL); TH1 *hf(NULL);
2086   hf = gs->GetHistogram();
2087   hf->SetTitle(at[0]);
2088   ax = hf->GetXaxis();
2089   ax->SetRangeUser(bb[0], bb[2]);
2090   ax->SetTitle(at[1]);ax->CenterTitle();
2091
2092   ax = hf->GetYaxis();
2093   ax->SetRangeUser(bb[1], bb[3]);
2094   ax->SetTitleOffset(1.1);
2095   ax->SetTitle(at[2]);ax->CenterTitle();
2096
2097   TGaxis *gax = NULL;
2098   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2099   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2100   //gax->SetVertical();
2101   gax->CenterTitle(); gax->SetTitleOffset(.7);
2102   gax->SetTitle(at[3]); gax->Draw();
2103
2104   // bounding box
2105   TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2106   b->SetFillStyle(3002);b->SetLineColor(0);
2107   b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2108   b->Draw();
2109   return kTRUE;
2110 }
2111
2112 //________________________________________________________
2113 Bool_t AliTRDresolution::GetGraphTrack(Float_t *bb, Int_t idx, Int_t first, Int_t n, Bool_t kLEG)
2114 {
2115   //
2116   // Get the graphs
2117   //
2118
2119   if(!fGraphS || !fGraphM) return kFALSE;
2120
2121   // axis titles look up
2122   Int_t nref(0);
2123   for(Int_t jp=0; jp<Int_t(kMCtrackTRD); jp++) nref+=fgNproj[jp];
2124   nref+=idx;
2125   const Char_t **at = fgAxTitle[nref];
2126   //printf("nref[%d] ax[%s] x[%f %f] y[%f %f]\n", nref, at[0], bb[0], bb[2], bb[1], bb[3]);
2127
2128   TLegend *legM(NULL), *legS(NULL);
2129   if(kLEG){
2130     legM=new TLegend(.35, .6, .65, .9);
2131     legM->SetHeader("Mean");
2132     legM->SetBorderSize(1);
2133     legM->SetFillColor(0);
2134     legS=new TLegend(.65, .6, .95, .9);
2135     legS->SetHeader("Sigma");
2136     legS->SetBorderSize(1);
2137     legS->SetFillColor(0);
2138   }
2139   Int_t layer(first/11);
2140   TH1S *h1(NULL);
2141   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), Form("%s %d;%s;%s", at[0], layer, at[1], at[2]), 2, bb[0], bb[2]);
2142   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2143   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2144   // axis range
2145   TAxis *ax = h1->GetXaxis();
2146   ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2147   ax = h1->GetYaxis();
2148   ax->SetRangeUser(bb[1], bb[3]);
2149   ax->CenterTitle(); ax->SetTitleOffset(1.4);
2150   h1->Draw();
2151 //   TGaxis *gax = NULL;
2152 //   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2153 //   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2154 //   //gax->SetVertical();
2155 //   gax->CenterTitle(); //gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
2156 //   gax->SetTitle(at[3]); gax->Draw();
2157
2158   TGraphErrors *gm = NULL, *gs = NULL;
2159   TObjArray *a0 = NULL, *a1 = NULL;
2160   a0 = (TObjArray*)((TObjArray*)fGraphM->At(kMCtrackTRD))->At(idx);
2161   a1 = (TObjArray*)((TObjArray*)fGraphS->At(kMCtrackTRD))->At(idx);
2162   Int_t nn(0), m(n/2);
2163   for(Int_t is=first, is0=0; is<first+n; is++, is0++){
2164     if(is0==m) continue; // no PID tracks
2165     if(is0==m-1||is0==m+1) continue; // electron tracks
2166     if(!(gs =  (TGraphErrors*)a1->At(is))) return kFALSE;
2167     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
2168
2169     if((nn=gs->GetN())){
2170       gs->Draw("pc"); if(legS) legS->AddEntry(gs, gs->GetTitle(), "pl");
2171       gs->Sort(&TGraph::CompareY);
2172       PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[0]);
2173       PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[nn-1]);
2174       gs->Sort(&TGraph::CompareX);
2175     }
2176     if(gm->GetN()){
2177       gm->Draw("pc");if(legM) legM->AddEntry(gm, gm->GetTitle(), "pl");
2178       PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gm->GetMean(2));
2179       PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gm->GetRMS(2));
2180     }
2181   }
2182   if(kLEG){legM->Draw();legS->Draw();}
2183   return kTRUE;
2184 }
2185
2186
2187 //________________________________________________________
2188 Bool_t AliTRDresolution::GetGraphTrackTPC(Float_t *bb, Int_t idx, Int_t ist, Int_t n, Bool_t kLEG)
2189 {
2190   //
2191   // Get the graphs
2192   //
2193
2194   if(!fGraphS || !fGraphM) return kFALSE;
2195
2196   // axis titles look up
2197   Int_t nref = 0;
2198   for(Int_t jp=0; jp<Int_t(kMCtrackTPC); jp++) nref+=fgNproj[jp];
2199   nref+=idx;
2200   const Char_t **at = fgAxTitle[nref];
2201
2202   TLegend *legM(NULL), *legS(NULL);
2203   if(kLEG){
2204     legM=new TLegend(.35, .6, .65, .9);
2205     legM->SetHeader("Mean");
2206     legM->SetBorderSize(1);
2207     legM->SetFillColor(0);
2208     legS=new TLegend(.65, .6, .95, .9);
2209     legS->SetHeader("Sigma");
2210     legS->SetBorderSize(1);
2211     legS->SetFillColor(0);
2212   }
2213   TH1S *h1(NULL);
2214   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), at[0], 2, bb[0], bb[2]);
2215   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2216   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2217   // axis range
2218   TAxis *ax = h1->GetXaxis();
2219   ax->SetTitle(at[1]);ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2220   ax = h1->GetYaxis();
2221   ax->SetRangeUser(bb[1], bb[3]);
2222   ax->SetTitleOffset(1.4);//ax->SetTitleSize(.06);
2223   ax->SetTitle(at[2]);ax->CenterTitle();
2224   h1->Draw();
2225 //   TGaxis *gax = NULL;
2226 //   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2227 //   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2228 //   //gax->SetVertical();
2229 //   gax->CenterTitle(); //gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
2230 //   gax->SetTitle(at[3]); gax->Draw();
2231
2232   Int_t nn(0), m(n/2);
2233   TGraphErrors *gm = NULL, *gs = NULL;
2234   TObjArray *a0 = NULL, *a1 = NULL;
2235   a0 = (TObjArray*)((TObjArray*)fGraphM->At(kMCtrackTPC))->At(idx);
2236   a1 = (TObjArray*)((TObjArray*)fGraphS->At(kMCtrackTPC))->At(idx);
2237   for(Int_t is=ist, is0=0; is<ist+n; is++, is0++){
2238     if(is0==m) continue; // no PID tracks
2239     //if(is0==m-1||is0==m+1) continue; // electron tracks
2240     if(!(gs =  (TGraphErrors*)a1->At(is))) return kFALSE;
2241     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
2242     if((nn=gs->GetN())){
2243       gs->Draw("pl");if(legS) legS->AddEntry(gs, gs->GetTitle(), "pl");
2244       gs->Sort(&TGraph::CompareY);
2245       PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[0]);
2246       PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[nn-1]);
2247       gs->Sort(&TGraph::CompareX);
2248     }
2249     if(gm->GetN()){
2250       gm->Draw("pl"); if(legM) legM->AddEntry(gm, gm->GetTitle(), "pl");
2251       PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gm->GetMean(2));
2252       PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gm->GetRMS(2));
2253     }
2254   }
2255   if(kLEG) {legM->Draw(); legS->Draw();}
2256   return kTRUE;
2257 }
2258
2259 //________________________________________________________
2260 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2261 {
2262   //
2263   // Get the most probable value and the full width half mean
2264   // of a Landau distribution
2265   //
2266
2267   const Float_t dx = 1.;
2268   mpv = f->GetParameter(1);
2269   Float_t fx, max = f->Eval(mpv);
2270
2271   xm = mpv - dx;
2272   while((fx = f->Eval(xm))>.5*max){
2273     if(fx>max){
2274       max = fx;
2275       mpv = xm;
2276     }
2277     xm -= dx;
2278   }
2279
2280   xM += 2*(mpv - xm);
2281   while((fx = f->Eval(xM))>.5*max) xM += dx;
2282 }
2283
2284
2285 //________________________________________________________
2286 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2287 {
2288
2289   fReconstructor->SetRecoParam(r);
2290 }