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