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