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