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