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