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