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 <TMath.h>
  67 #include <TMatrixT.h>
  68 #include <TVectorT.h>
  69 #include "TTreeStream.h"
  70 #include "TGeoManager.h"
  71
  72 #include "AliAnalysisManager.h"
  73 #include "AliTrackReference.h"
  74 #include "AliTrackPointArray.h"
  75 #include "AliCDBManager.h"
  76
  77 #include "AliTRDcalibDB.h"
  78 #include "AliTRDCommonParam.h"
  79 #include "AliTRDSimParam.h"
  80 #include "AliTRDgeometry.h"
  81 #include "AliTRDpadPlane.h"
  82 #include "AliTRDcluster.h"
  83 #include "AliTRDseedV1.h"
  84 #include "AliTRDtrackV1.h"
  85 #include "AliTRDtrackerV1.h"
  86 #include "AliTRDReconstructor.h"
  87 #include "AliTRDrecoParam.h"
  88
  89 #include "info/AliTRDclusterInfo.h"
  90 #include "info/AliTRDtrackInfo.h"
  91 #include "AliTRDresolution.h"
  92
  93 ClassImp(AliTRDresolution)
  94 UChar_t AliTRDresolution::fNElements[kNhistos] = {
  95   2, 5, 5,
  96   2, 5, 10, 2, 10
  97 };
  98 Char_t *AliTRDresolution::fAxTitle[32][3] = {
  99   // ESD
 100   {"tg(#phi)", "PULL: #mu_{y}^{cl}", "PULL: #sigma_{y}^{cl}"}
 101  ,{"tg(#phi)", "PULL: #mu_{y}^{trklt}", "PULL: #sigma_{y}^{trklt}"}
 102  ,{"tg(#phi)", "PULL: #mu_{#phi}^{trklt}", "PULL: #sigma_{#phi}^{trklt}"}
 103   // MC cluster
 104  ,{"tg(#phi)", "#mu_{y}^{cl} [#mum]", "#sigma_{y}^{cl} [#mum]"}
 105  ,{"tg(#phi)", "PULL: #mu_{y}^{cl}", "PULL: #sigma_{y}^{cl}"}
 106   // MC tracklet
 107  ,{"tg(#phi)", "#mu_{y}^{trklt} [#mum]", "#sigma_{y}^{trklt} [#mum]"}
 108  ,{"tg(#phi)", "PULL: #mu_{y}^{trklt}", "PULL: #sigma_{y}^{trklt}"}
 109  ,{"tg(#theta)", "#mu_{z}^{trklt} [#mum]", "#sigma_{z}^{trklt} [#mum]"}
 110  ,{"tg(#theta)", "PULL: #mu_{z}^{trklt}", "PULL: #sigma_{z}^{trklt}"}
 111  ,{"tg(#phi)", "#mu_{#phi}^{trklt} [mrad]", "#sigma_{#phi}^{trklt} [mrad]"}
 112   // MC track TPC
 113  ,{"tg(#phi)", "#mu_{y}^{TPC trk} [#mum]", "#sigma_{y}^{TPC trk} [#mum]"}
 114  ,{"tg(#phi)", "PULL: #mu_{y}^{TPC trk}", "PULL: #sigma_{y}^{TPC trk}"}
 115  ,{"tg(#theta)", "#mu_{z}^{TPC trk} [#mum]", "#sigma_{z}^{TPC trk} [#mum]"}
 116  ,{"tg(#theta)", "PULL: #mu_{z}^{TPC trk}", "PULL: #sigma_{z}^{TPC trk}"}
 117  ,{"tg(#phi)", "#mu_{#phi}^{TPC trk} [mrad]", "#sigma_{#phi}^{TPC trk} [mrad]"}
 118  ,{"tg(#phi)", "PULL: #mu_{snp}^{TPC trk}", "PULL: #sigma_{snp}^{TPC trk}"}
 119  ,{"tg(#theta)", "#mu_{#theta}^{TPC trk} [mrad]", "#sigma_{#theta}^{TPC trk} [mrad]"}
 120  ,{"tg(#theta)", "PULL: #mu_{tgl}^{TPC trk}", "PULL: #sigma_{tgl}^{TPC trk}"}
 121  ,{"p_{t}^{MC} [GeV/c]", "#mu^{TPC trk}(#Deltap_{t}/p_{t}^{MC}) [%]", "#sigma^{TPC trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
 122  ,{"1/p_{t}^{MC} [c/GeV]", "PULL: #mu_{1/p_{t}}^{TPC trk}", "PULL: #sigma_{1/p_{t}}^{TPC trk}"}
 123   // MC track HMPID
 124  ,{"tg(#theta)", "#mu_{z}^{trk} [#mum]", "#sigma_{z}^{trk} [#mum]"}
 125  ,{"tg(#theta)", "PULL: #mu_{z}^{trk}", "PULL: #sigma_{z}^{trk}"}
 126   // MC track in TRD
 127  ,{"tg(#phi)", "#mu_{y}^{trk} [#mum]", "#sigma_{y}^{trk} [#mum]"}
 128  ,{"tg(#phi)", "PULL: #mu_{y}^{trk}", "PULL: #sigma_{y}^{trk}"}
 129  ,{"tg(#theta)", "#mu_{z}^{Trk} [#mum]", "#sigma_{z}^{Trk} [#mum]"}
 130  ,{"tg(#theta)", "PULL: #mu_{z}^{Trk}", "PULL: #sigma_{z}^{Trk}"}
 131  ,{"tg(#phi)", "#mu_{#phi}^{Trk} [mrad]", "#sigma_{#phi}^{Trk} [mrad]"}
 132  ,{"tg(#phi)", "PULL: #mu_{snp}^{Trk}", "PULL: #sigma_{snp}^{Trk}"}
 133  ,{"tg(#theta)", "#mu_{#theta}^{Trk} [mrad]", "#sigma_{#theta}^{Trk} [mrad]"}
 134  ,{"tg(#theta)", "PULL: #mu_{tgl}^{Trk}", "PULL: #sigma_{tgl}^{Trk}"}
 135  ,{"p_{t}^{MC} [GeV/c]", "#mu^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]", "#sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
 136  ,{"1/p_{t}^{MC} [c/GeV]", "PULL: #mu_{1/p_{t}}^{Trk}", "PULL: #sigma_{1/p_{t}}^{Trk}"}
 137 };
 138
 139 //________________________________________________________
 140 AliTRDresolution::AliTRDresolution()
 141   :AliTRDrecoTask("resolution", "Spatial and momentum TRD resolution checker")
 142   ,fStatus(0)
 143   ,fReconstructor(0x0)
 144   ,fGeo(0x0)
 145   ,fGraphS(0x0)
 146   ,fGraphM(0x0)
 147   ,fCl(0x0)
 148   ,fTrklt(0x0)
 149   ,fMCcl(0x0)
 150   ,fMCtrklt(0x0)
 151 {
 152   fReconstructor = new AliTRDReconstructor();
 153   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 154   fGeo = new AliTRDgeometry();
 155
 156   InitFunctorList();
 157
 158   DefineOutput(1, TObjArray::Class()); // cluster2track
 159   DefineOutput(2, TObjArray::Class()); // tracklet2track
 160   DefineOutput(3, TObjArray::Class()); // cluster2mc
 161   DefineOutput(4, TObjArray::Class()); // tracklet2mc
 162 }
 163
 164 //________________________________________________________
 165 AliTRDresolution::~AliTRDresolution()
 166 {
 167   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 168   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 169   delete fGeo;
 170   delete fReconstructor;
 171   if(gGeoManager) delete gGeoManager;
 172   if(fCl){fCl->Delete(); delete fCl;}
 173   if(fTrklt){fTrklt->Delete(); delete fTrklt;}
 174   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 175   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
 176 }
 177
 178
 179 //________________________________________________________
 180 void AliTRDresolution::CreateOutputObjects()
 181 {
 182   // spatial resolution
 183   OpenFile(0, "RECREATE");
 184
 185   fContainer = Histos();
 186
 187   fCl = new TObjArray();
 188   fCl->SetOwner(kTRUE);
 189   fTrklt = new TObjArray();
 190   fTrklt->SetOwner(kTRUE);
 191   fMCcl = new TObjArray();
 192   fMCcl->SetOwner(kTRUE);
 193   fMCtrklt = new TObjArray();
 194   fMCtrklt->SetOwner(kTRUE);
 195 }
 196
 197 //________________________________________________________
 198 void AliTRDresolution::Exec(Option_t *opt)
 199 {
 200   fCl->Delete();
 201   fTrklt->Delete();
 202   fMCcl->Delete();
 203   fMCtrklt->Delete();
 204
 205   AliTRDrecoTask::Exec(opt);
 206
 207   PostData(1, fCl);
 208   PostData(2, fTrklt);
 209   PostData(3, fMCcl);
 210   PostData(4, fMCtrklt);
 211 }
 212
 213 //________________________________________________________
 214 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 215 {
 216   if(track) fTrack = track;
 217   if(!fTrack){
 218     AliWarning("No Track defined.");
 219     return 0x0;
 220   }
 221   TObjArray *arr = 0x0;
 222   if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
 223     AliWarning("No output container defined.");
 224     return 0x0;
 225   }
 226
 227   Double_t cov[3];
 228   Float_t x0, y0, z0, dy, dydx, dzdx;
 229   AliTRDseedV1 *fTracklet = 0x0;
 230   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 231     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 232     if(!fTracklet->IsOK()) continue;
 233     x0 = fTracklet->GetX0();
 234
 235     // retrive the track angle with the chamber
 236     y0   = fTracklet->GetYref(0);
 237     z0   = fTracklet->GetZref(0);
 238     dydx = fTracklet->GetYref(1);
 239     dzdx = fTracklet->GetZref(1);
 240     fTracklet->GetCovRef(cov);
 241     Float_t tilt = fTracklet->GetTilt();
 242     AliTRDcluster *c = 0x0;
 243     fTracklet->ResetClusterIter(kFALSE);
 244     while((c = fTracklet->PrevCluster())){
 245       Float_t xc = c->GetX();
 246       Float_t yc = c->GetY();
 247       Float_t zc = c->GetZ();
 248       Float_t dx = x0 - xc;
 249       Float_t yt = y0 - dx*dydx;
 250       Float_t zt = z0 - dx*dzdx;
 251       yc -= tilt*(zc-zt); // tilt correction
 252       dy = yt - yc;
 253
 254       Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
 255       Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2;
 256       ((TH2I*)arr->At(0))->Fill(dydx, dy);
 257       ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] + sx2 + sy2));
 258
 259       if(fDebugLevel>=1){
 260         // Get z-position with respect to anode wire
 261         //AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 262         Int_t istk = fGeo->GetStack(c->GetDetector());
 263         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 264         Float_t row0 = pp->GetRow0();
 265         Float_t d  =  row0 - zt + pp->GetAnodeWireOffset();
 266         d -= ((Int_t)(2 * d)) / 2.0;
 267         if (d > 0.25) d  = 0.5 - d;
 268
 269 /*        AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 270         fCl->Add(clInfo);
 271         clInfo->SetCluster(c);
 272         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
 273         clInfo->SetResolution(dy);
 274         clInfo->SetAnisochronity(d);
 275         clInfo->SetDriftLength(dx);
 276         (*fDebugStream) << "ClusterResiduals"
 277           <<"clInfo.=" << clInfo
 278           << "\n";*/
 279       }
 280     }
 281   }
 282   return (TH2I*)arr->At(0);
 283 }
 284
 285
 286 //________________________________________________________
 287 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 288 {
 289 // Plot normalized residuals for tracklets to track.
 290 //
 291 // We start from the result that if X=N(|m|, |Cov|)
 292 // BEGIN_LATEX
 293 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
 294 // END_LATEX
 295 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
 296 // reference position.
 297   if(track) fTrack = track;
 298   if(!fTrack){
 299     AliWarning("No Track defined.");
 300     return 0x0;
 301   }
 302   TObjArray *arr = 0x0;
 303   if(!(arr = (TObjArray*)fContainer->At(kTracklet))){
 304     AliWarning("No output container defined.");
 305     return 0x0;
 306   }
 307
 308   Double_t cov[3], covR[3]/*, sqr[3], inv[3]*/;
 309   Float_t x, dx, dy, dz;
 310   AliTRDseedV1 *fTracklet = 0x0;
 311   for(Int_t il=AliTRDgeometry::kNlayer; il--;){
 312     if(!(fTracklet = fTrack->GetTracklet(il))) continue;
 313     if(!fTracklet->IsOK()) continue;
 314     x    = fTracklet->GetX();
 315     dx   = fTracklet->GetX0() - x;
 316     // compute dy^2 and dz^2
 317     dy   = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
 318     dz   = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
 319     // compute covariance matrix
 320     fTracklet->GetCovAt(x, cov);
 321     fTracklet->GetCovRef(covR);
 322     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 323 /*  // Correct PULL calculation by considering off
 324     // diagonal elements in the covariance matrix
 325     // compute square root matrix
 326     if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
 327     if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
 328     Double_t y = sqr[0]*dy+sqr[1]*dz;
 329     Double_t z = sqr[1]*dy+sqr[2]*dz;
 330     ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
 331
 332     ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
 333     ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
 334     ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
 335     if(!fTracklet->IsRowCross()) continue;
 336     ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
 337     ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
 338   }
 339
 340
 341   return (TH2I*)arr->At(0);
 342 }
 343
 344
 345 //________________________________________________________
 346 TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
 347 {
 348   if(track) fTrack = track;
 349   if(!fTrack){
 350     AliWarning("No Track defined.");
 351     return 0x0;
 352   }
 353   AliExternalTrackParam *tin = 0x0;
 354   if(!(tin = fTrack->GetTrackLow())){
 355     AliWarning("Track did not entered TRD fiducial volume.");
 356     return 0x0;
 357   }
 358   TH1 *h = 0x0;
 359
 360   Double_t x = tin->GetX();
 361   AliTRDseedV1 *tracklet = 0x0;
 362   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 363     if(!(tracklet = fTrack->GetTracklet(ily))) continue;
 364     break;
 365   }
 366   if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
 367     AliWarning("Tracklet did not match TRD entrance.");
 368     return 0x0;
 369   }
 370   const Int_t kNPAR(5);
 371   Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
 372   Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
 373   Double_t cov[3]; tracklet->GetCovAt(x, cov);
 374
 375   // define sum covariances
 376   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
 377   Double_t *pc = &covR[0], *pp = &parR[0];
 378   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
 379     PAR(ir) = (*pp);
 380     for(Int_t ic = 0; ic<=ir; ic++,pc++){
 381       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
 382     }
 383   }
 384   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
 385   //COV.Print(); PAR.Print();
 386
 387   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
 388   Double_t dy = parR[0] - tracklet->GetY();
 389   TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
 390   ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
 391   ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
 392   if(tracklet->IsRowCross()){
 393     Double_t dz = parR[1] - tracklet->GetZ();
 394     ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
 395     ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
 396   }
 397   Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1));  ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
 398
 399
 400   // register reference histo for mini-task
 401   h = (TH2I*)arr->At(0);
 402
 403   if(fDebugLevel>=1){
 404     (*fDebugStream) << "trackIn"
 405       << "x="       << x
 406       << "P="       << &PAR
 407       << "C="       << &COV
 408       << "\n";
 409
 410     Double_t y = tracklet->GetY();
 411     Double_t z = tracklet->GetZ();
 412     (*fDebugStream) << "trackletIn"
 413       << "y="       << y
 414       << "z="       << z
 415       << "Vy="      << cov[0]
 416       << "Cyz="     << cov[1]
 417       << "Vz="      << cov[2]
 418       << "\n";
 419   }
 420
 421
 422   if(!HasMCdata()) return h;
 423   UChar_t s;
 424   Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
 425   if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
 426   // translate to reference radial position
 427   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
 428   //Fill MC info
 429   TVectorD PARMC(kNPAR);
 430   PARMC[0]=y0; PARMC[1]=z0;
 431   PARMC[2]=dydx0/TMath::Sqrt(1.+dydx0*dydx0); PARMC[3]=dzdx0;
 432   PARMC[4]=1./pt0;
 433
 434 //   TMatrixDSymEigen eigen(COV);
 435 //   TVectorD evals = eigen.GetEigenValues();
 436 //   TMatrixDSym evalsm(kNPAR);
 437 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
 438 //   TMatrixD evecs = eigen.GetEigenVectors();
 439 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
 440
 441   // fill histos
 442   arr = (TObjArray*)fContainer->At(kMCtrackTPC);
 443   // y resolution/pulls
 444   ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
 445   ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
 446   // z resolution/pulls
 447   ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
 448   ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
 449   // phi resolution/snp pulls
 450   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
 451   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
 452   // theta resolution/tgl pulls
 453   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
 454   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
 455   // pt resolution/1/pt pulls
 456   ((TH2I*)arr->At(8))->Fill(pt0, 1.-PARMC[4]/PAR[4]);
 457   ((TH2I*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)));
 458
 459   // fill debug for MC
 460   if(fDebugLevel>=1){
 461     (*fDebugStream) << "trackInMC"
 462       << "P="   << &PARMC
 463       << "\n";
 464   }
 465   return h;
 466 }
 467
 468 //________________________________________________________
 469 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 470 {
 471   if(!HasMCdata()){
 472     AliWarning("No MC defined. Results will not be available.");
 473     return 0x0;
 474   }
 475   if(track) fTrack = track;
 476   if(!fTrack){
 477     AliWarning("No Track defined.");
 478     return 0x0;
 479   }
 480   TObjArray *arr = 0x0;
 481   TH1 *h = 0x0;
 482   UChar_t s;
 483   Int_t pdg = fMC->GetPDG(), det=-1;
 484   Int_t label = fMC->GetLabel();
 485   Double_t xAnode, x, y, z, pt, dydx, dzdx;
 486   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 487   Double_t covR[3]/*, cov[3]*/;
 488
 489   if(fDebugLevel>=1){
 490     Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
 491     fMC->PropagateKalman(DX, DY, DZ, DPt, COV);
 492     (*fDebugStream) << "MCkalman"
 493       << "pdg="  << pdg
 494       << "dx0="  << DX[0]
 495       << "dx1="  << DX[1]
 496       << "dx2="  << DX[2]
 497       << "dy0="  << DY[0]
 498       << "dy1="  << DY[1]
 499       << "dy2="  << DY[2]
 500       << "dz0="  << DZ[0]
 501       << "dz1="  << DZ[1]
 502       << "dz2="  << DZ[2]
 503       << "dpt0=" << DPt[0]
 504       << "dpt1=" << DPt[1]
 505       << "dpt2=" << DPt[2]
 506       << "\n";
 507   }
 508
 509   AliTRDseedV1 *fTracklet = 0x0;
 510   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 511     if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
 512        !fTracklet->IsOK())*/ continue;
 513
 514     det = fTracklet->GetDetector();
 515     x0  = fTracklet->GetX0();
 516     //radial shift with respect to the MC reference (radial position of the pad plane)
 517     x= fTracklet->GetX();
 518     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
 519     xAnode  = fTracklet->GetX0();
 520
 521     // MC track position at reference radial position
 522     dx  = x0 - x;
 523     if(fDebugLevel>=1){
 524       (*fDebugStream) << "MC"
 525         << "det="     << det
 526         << "pdg="     << pdg
 527         << "pt="      << pt0
 528         << "x="       << x0
 529         << "y="       << y0
 530         << "z="       << z0
 531         << "dydx="    << dydx0
 532         << "dzdx="    << dzdx0
 533         << "\n";
 534     }
 535     Float_t yt = y0 - dx*dydx0;
 536     Float_t zt = z0 - dx*dzdx0;
 537     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
 538
 539     // Kalman position at reference radial position
 540     dx = xAnode - x;
 541     y  = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
 542     dy = yt - y;
 543     z  = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
 544     dz = zt - z;
 545     dydx = fTracklet->GetYref(1);
 546     dzdx = fTracklet->GetTgl();
 547     pt = TMath::Abs(fTracklet->GetPt());
 548     fTracklet->GetCovRef(covR);
 549
 550     arr = (TObjArray*)fContainer->At(kMCtrack);
 551     // y resolution/pulls
 552     ((TH2I*)arr->At(0))->Fill(dydx0, dy);
 553     ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
 554     // z resolution/pulls
 555     ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
 556     ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
 557     // phi resolution/ snp pulls
 558     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
 559     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
 560     //TODO ((TH2I*)arr->At(5))->Fill(dydx0, );
 561     // theta resolution/ tgl pulls
 562     Double_t dtgl = (dzdx - dzdx0)/(1.- dzdx*dzdx0);
 563     ((TH2I*)arr->At(6))->Fill(dzdx0,
 564     TMath::ATan(dtgl));
 565     //TODO ((TH2I*)arr->At(7))->Fill(dydx0, );
 566     // pt resolution/ 1/pt pulls
 567     if(pdg!=kElectron && pdg!=kPositron){
 568       ((TH2I*)arr->At(8))->Fill(pt0, 1.-pt/pt0);
 569       //TODO ((TH2I*)arr->At(9))->Fill(1./pt0, (pt0/pt-1.)/TMath::Sqrt(covR[4]));
 570     }
 571     // Fill Debug stream for Kalman track
 572     if(fDebugLevel>=1){
 573       (*fDebugStream) << "MCtrack"
 574         << "pt="      << pt
 575         << "x="       << x
 576         << "y="       << y
 577         << "z="       << z
 578         << "dydx="    << dydx
 579         << "dzdx="    << dzdx
 580         << "s2y="     << covR[0]
 581         << "s2z="     << covR[2]
 582         << "\n";
 583     }
 584
 585     // recalculate tracklet based on the MC info
 586     AliTRDseedV1 tt(*fTracklet);
 587     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
 588     tt.SetZref(1, dzdx0);
 589     tt.Fit(kTRUE, kTRUE);
 590     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
 591     dydx = tt.GetYfit(1);
 592     dx = x0 - x;
 593     yt = y0 - dx*dydx0;
 594     zt = z0 - dx*dzdx0;
 595     Bool_t rc = tt.IsRowCross();
 596
 597     // add tracklet residuals for y and dydx
 598     arr = (TObjArray*)fContainer->At(kMCtracklet);
 599     if(!rc){
 600       dy    = yt-y;
 601
 602       Float_t dphi  = (dydx - dydx0);
 603       dphi /= 1.- dydx*dydx0;
 604
 605       ((TH2I*)arr->At(0))->Fill(dydx0, dy);
 606       if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
 607       ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
 608     } else {
 609       // add tracklet residuals for z
 610       dz = zt-z;
 611       ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
 612       if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(tt.GetS2Z()));
 613     }
 614
 615     // Fill Debug stream for tracklet
 616     if(fDebugLevel>=1){
 617       Float_t s2y = tt.GetS2Y();
 618       Float_t s2z = tt.GetS2Z();
 619       (*fDebugStream) << "MCtracklet"
 620         << "rc="    << rc
 621         << "x="     << x
 622         << "y="     << y
 623         << "z="     << z
 624         << "dydx="  << dydx
 625         << "s2y="   << s2y
 626         << "s2z="   << s2z
 627         << "\n";
 628     }
 629
 630     Int_t istk = AliTRDgeometry::GetStack(det);
 631     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 632     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 633     Float_t tilt = fTracklet->GetTilt();
 634     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
 635
 636     arr = (TObjArray*)fContainer->At(kMCcluster);
 637     AliTRDcluster *c = 0x0;
 638     fTracklet->ResetClusterIter(kFALSE);
 639     while((c = fTracklet->PrevCluster())){
 640       Float_t  q = TMath::Abs(c->GetQ());
 641       x = c->GetX(); y = c->GetY();
 642 //       Int_t col = c->GetPadCol();
 643 //       Int_t row = c->GetPadRow();
 644 //       Double_t cw = pp->GetColSize(col);
 645 //       Double_t y0 = pp->GetColPos(col) + 0.5 * cw;
 646 //       Double_t s2 = AliTRDcalibDB::Instance()->GetPRFWidth(det, col, row); s2 *= s2; s2 -= - 1.5e-1;
 647 //       y = c->GetYloc(y0, s2, cw); y-=(xAnode-x)*exb;
 648
 649       z = c->GetZ();
 650       dx = x0 - x;
 651       yt = y0 - dx*dydx0;
 652       zt = z0 - dx*dzdx0;
 653       dy = yt - (y - tilt*(z-zt));
 654       Float_t sx2 = dydx0*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
 655       Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2;
 656
 657       // Fill Histograms
 658       if(q>20. && q<250.){
 659         ((TH2I*)arr->At(0))->Fill(dydx0, dy);
 660         ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(sx2+sy2));
 661       }
 662
 663       // Fill calibration container
 664       Float_t d = zr0 - zt;
 665       d -= ((Int_t)(2 * d)) / 2.0;
 666       if (d > 0.25) d  = 0.5 - d;
 667       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 668       fMCcl->Add(clInfo);
 669       clInfo->SetCluster(c);
 670       clInfo->SetMC(pdg, label);
 671       clInfo->SetGlobalPosition(yt, zt, dydx0, dzdx0);
 672       clInfo->SetResolution(dy);
 673       clInfo->SetAnisochronity(d);
 674       clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5);
 675       //dx-.5*AliTRDgeometry::CamHght());
 676       clInfo->SetTilt(tilt);
 677
 678       // Fill Debug Tree
 679       if(fDebugLevel>=2){
 680         //clInfo->Print();
 681         (*fDebugStream) << "MCcluster"
 682           <<"clInfo.=" << clInfo
 683           << "\n";
 684       }
 685     }
 686   }
 687   return h;
 688 }
 689
 690
 691 //________________________________________________________
 692 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
 693 {
 694   Float_t xy[4] = {0., 0., 0., 0.};
 695   if(!gPad){
 696     AliWarning("Please provide a canvas to draw results.");
 697     return kFALSE;
 698   }
 699   TList *l = 0x0;
 700   switch(ifig){
 701   case kCluster:
 702     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 703     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 2000.;
 704     ((TVirtualPad*)l->At(0))->cd();
 705     if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
 706     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 707     ((TVirtualPad*)l->At(1))->cd();
 708     if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
 709     return kTRUE;
 710   case kTracklet:
 711     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 712     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
 713     ((TVirtualPad*)l->At(0))->cd();
 714     if(!GetGraphPlot(&xy[0], kTracklet, 0)) break;
 715     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 716     ((TVirtualPad*)l->At(1))->cd();
 717     if(!GetGraphPlot(&xy[0], kTracklet, 1)) break;
 718     return kTRUE;
 719   case 2: // kTracklet z
 720     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 721     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
 722     ((TVirtualPad*)l->At(0))->cd();
 723     if(!GetGraphPlot(&xy[0], kTracklet, 2)) break;
 724     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 725     ((TVirtualPad*)l->At(1))->cd();
 726     if(!GetGraphPlot(&xy[0], kTracklet, 3)) break;
 727     return kTRUE;
 728   case 3: // kTracklet phi
 729     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 15.;
 730     if(GetGraphPlot(&xy[0], kTracklet, 4)) return kTRUE;
 731     break;
 732   case 4: // kTrackTPC y
 733     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 734     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
 735     ((TVirtualPad*)l->At(0))->cd();
 736     if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
 737     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 738     ((TVirtualPad*)l->At(1))->cd();
 739     if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
 740     return kTRUE;
 741   case 5: // kTrackTPC z
 742     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 743     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
 744     ((TVirtualPad*)l->At(0))->cd();
 745     if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
 746     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 747     ((TVirtualPad*)l->At(1))->cd();
 748     if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
 749     return kTRUE;
 750   case 6: // kTrackTPC phi
 751     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 15.;
 752     if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
 753     break;
 754   case 7: // kMCcluster
 755     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 756     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
 757     ((TVirtualPad*)l->At(0))->cd();
 758     if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
 759     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
 760     ((TVirtualPad*)l->At(1))->cd();
 761     if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
 762     return kTRUE;
 763   case 8: //kMCtracklet [y]
 764     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 765     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
 766     ((TVirtualPad*)l->At(0))->cd();
 767     if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
 768     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
 769     ((TVirtualPad*)l->At(1))->cd();
 770     if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
 771     return kTRUE;
 772   case 9: //kMCtracklet [z]
 773     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 774     xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
 775     ((TVirtualPad*)l->At(0))->cd();
 776     if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
 777     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 778     ((TVirtualPad*)l->At(1))->cd();
 779     if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
 780     return kTRUE;
 781   case 10: //kMCtracklet [phi]
 782     xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
 783     if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
 784     return kTRUE;
 785   case 11: //kMCtrack [y]
 786     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 787     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =200.;
 788     ((TVirtualPad*)l->At(0))->cd();
 789     if(!GetGraphPlot(&xy[0], kMCtrack, 0)) break;
 790     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
 791     ((TVirtualPad*)l->At(1))->cd();
 792     if(!GetGraphPlot(&xy[0], kMCtrack, 1)) break;
 793     return kTRUE;
 794   case 12: //kMCtrack [z]
 795     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 796     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =2000.;
 797     ((TVirtualPad*)l->At(0))->cd();
 798     if(!GetGraphPlot(&xy[0], kMCtrack, 2)) break;
 799     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 800     ((TVirtualPad*)l->At(1))->cd();
 801     if(!GetGraphPlot(&xy[0], kMCtrack, 3)) break;
 802     return kTRUE;
 803   case 13: //kMCtrack [phi/snp]
 804     //gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 805     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =200.;
 806     //((TVirtualPad*)l->At(0))->cd();
 807     if(!GetGraphPlot(&xy[0], kMCtrack, 4)) break;
 808 //     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
 809 //     ((TVirtualPad*)l->At(1))->cd();
 810 //     if(!GetGraphPlot(&xy[0], kMCtrack, 3)) break;
 811     return kTRUE;
 812   case 14: //kMCtrack [theta/tgl]
 813     //gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 814     xy[0]=-1.; xy[1]=-50.; xy[2]=1.; xy[3] =200.;
 815     //((TVirtualPad*)l->At(0))->cd();
 816     if(!GetGraphPlot(&xy[0], kMCtrack, 6)) break;
 817 //     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
 818 //     ((TVirtualPad*)l->At(1))->cd();
 819 //     if(!GetGraphPlot(&xy[0], kMCtrack, 3)) break;
 820     return kTRUE;
 821   case 15: //kMCtrack [pt]
 822     xy[0] = 0.; xy[1] = -0.5; xy[2] = 15.; xy[3] = 5.5;
 823     if(!GetGraphPlot(&xy[0], kMCtrack, 2)) break;
 824     return kTRUE;
 825   case 16: // kMCtrackTPC [y]
 826     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 827     xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
 828     ((TVirtualPad*)l->At(0))->cd();
 829     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
 830     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
 831     ((TVirtualPad*)l->At(1))->cd();
 832     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
 833     return kTRUE;
 834   case 17: // kMCtrackTPC [z]
 835     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 836     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
 837     ((TVirtualPad*)l->At(0))->cd();
 838     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
 839     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 840     ((TVirtualPad*)l->At(1))->cd();
 841     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
 842     return kTRUE;
 843   case 18: // kMCtrackTPC [phi|snp]
 844     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 845     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
 846     ((TVirtualPad*)l->At(0))->cd();
 847     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
 848     //xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
 849     ((TVirtualPad*)l->At(1))->cd();
 850     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
 851     return kTRUE;
 852   case 19: // kMCtrackTPC [theta|tgl]
 853     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 854     xy[0]=-1.; xy[1]=-10.5; xy[2]=1.; xy[3] =20.5;
 855     ((TVirtualPad*)l->At(0))->cd();
 856     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
 857     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
 858     ((TVirtualPad*)l->At(1))->cd();
 859     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
 860     return kTRUE;
 861   case 20: // kMCtrackTPC [pt]
 862     gPad->Divide(2, 1); l=gPad->GetListOfPrimitives();
 863     xy[0] = 0.; xy[1] = -.5; xy[2] = 15.; xy[3] = 2.5;
 864     ((TVirtualPad*)l->At(0))->cd();
 865     if(!GetGraphPlot(xy, AliTRDresolution::kMCtrackTPC, 8)) break;
 866     xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =2.5;
 867     ((TVirtualPad*)l->At(1))->cd();
 868     if(!GetGraphPlot(xy, AliTRDresolution::kMCtrackTPC, 9)) break;
 869     return kTRUE;
 870   case 21:  // kMCtrackHMPID [z]
 871     return kTRUE;
 872   }
 873   AliInfo(Form("Reference plot [%d] missing result", ifig));
 874   return kFALSE;
 875 }
 876
 877
 878 //________________________________________________________
 879 Bool_t AliTRDresolution::PostProcess()
 880 {
 881   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 882   if (!fContainer) {
 883     Printf("ERROR: list not available");
 884     return kFALSE;
 885   }
 886   TGraphErrors *gm= 0x0, *gs= 0x0;
 887   if(!fGraphS && !fGraphM){
 888     TObjArray *aM(0x0), *aS(0x0);
 889     Int_t n = fContainer->GetEntriesFast();
 890     fGraphS = new TObjArray(n); fGraphS->SetOwner();
 891     fGraphM = new TObjArray(n); fGraphM->SetOwner();
 892     for(Int_t ig=0; ig<n; ig++){
 893       if(fNElements[ig]>1){
 894         fGraphM->AddAt(aM = new TObjArray(fNElements[ig]), ig);
 895         fGraphS->AddAt(aS = new TObjArray(fNElements[ig]), ig);
 896       } else {
 897         aM = fGraphM;aS = fGraphS;
 898       }
 899       for(Int_t ic=0; ic<fNElements[ig]; ic++){
 900         aS->AddAt(gs = new TGraphErrors(), fNElements[ig]>1?ic:ig);
 901         gs->SetMarkerStyle(23);
 902         gs->SetMarkerColor(kRed);
 903         gs->SetLineColor(kRed);
 904         gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
 905
 906         aM->AddAt(gm = new TGraphErrors(), fNElements[ig]>1?ic:ig);
 907         gm->SetLineColor(kBlue);
 908         gm->SetMarkerStyle(7);
 909         gm->SetMarkerColor(kBlue);
 910         gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
 911       }
 912     }
 913   }
 914
 915
 916   // DEFINE MODELS
 917   // simple gauss
 918   TF1 f("f1", "gaus", -.5, .5);
 919   // gauss on a constant background
 920   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 921   // gauss on a gauss background
 922   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 923
 924
 925   //PROCESS EXPERIMENTAL DISTRIBUTIONS
 926   // Clusters residuals
 927   Process(kCluster, 0, &f, 1.e4);
 928   Process(kCluster, 1, &f);
 929   fNRefFigures = 1;
 930   // Tracklet residual/pulls
 931   Process(kTracklet, 0, &f, 1.e4);
 932   Process(kTracklet, 1, &f);
 933   Process(kTracklet, 2, &f, 1.e4);
 934   Process(kTracklet, 3, &f);
 935   Process(kTracklet, 4, &f, 1.e3);
 936   fNRefFigures = 4;
 937   // TPC track residual/pulls
 938   Process(kTrackTPC, 0, &f, 1.e4);
 939   Process(kTrackTPC, 1, &f);
 940   Process(kTrackTPC, 2, &f, 1.e4);
 941   Process(kTrackTPC, 3, &f);
 942   Process(kTrackTPC, 4, &f, 1.e3);
 943   fNRefFigures = 7;
 944
 945   if(!HasMCdata()) return kTRUE;
 946
 947
 948   //PROCESS MC RESIDUAL DISTRIBUTIONS
 949
 950   // CLUSTER Y RESOLUTION/PULLS
 951   Process(kMCcluster, 0, &f, 1.e4);
 952   Process(kMCcluster, 1, &f);
 953   fNRefFigures = 8;
 954
 955   // TRACKLET RESOLUTION/PULLS
 956   Process(kMCtracklet, 0, &f, 1.e4); // y
 957   Process(kMCtracklet, 1, &f);       // y pulls
 958   Process(kMCtracklet, 2, &f, 1.e4); // z
 959   Process(kMCtracklet, 3, &f);       // z pulls
 960   Process(kMCtracklet, 4, &f, 1.e3); // phi
 961   fNRefFigures = 11;
 962
 963   // TRACK RESOLUTION/PULLS
 964   Process(kMCtrack, 0, &f, 1.e4);   // y
 965   Process(kMCtrack, 1, &f);         // y PULL
 966   Process(kMCtrack, 2, &f, 1.e4);   // z
 967   Process(kMCtrack, 3, &f);         // z PULL
 968   Process(kMCtrack, 4, &f, 1.e3);   // phi
 969   //Process(kMCtrack, 5, &f);         // snp PULL
 970   Process(kMCtrack, 6, &f, 1.e3);   // theta
 971   //Process(kMCtrack, 7, &f);         // tgl PULL
 972   Process(kMCtrack, 8, &f, 1.e2);   // pt resolution
 973   //Process(kMCtrack, 9, &f);         // 1/pt pulls
 974   fNRefFigures = 16;
 975
 976
 977   // TRACK TPC RESOLUTION/PULLS
 978   Process(kMCtrackTPC, 0, &f, 1.e4);// y resolution
 979   Process(kMCtrackTPC, 1, &f);      // y pulls
 980   Process(kMCtrackTPC, 2, &f, 1.e4);// z resolution
 981   Process(kMCtrackTPC, 3, &f);      // z pulls
 982   Process(kMCtrackTPC, 4, &f, 1.e3);// phi resolution
 983   Process(kMCtrackTPC, 5, &f);      // snp pulls
 984   Process(kMCtrackTPC, 6, &f, 1.e3);// theta resolution
 985   Process(kMCtrackTPC, 7, &f);      // tgl pulls
 986   Process(kMCtrackTPC, 8, &f, 1.e2);// pt resolution
 987   Process(kMCtrackTPC, 9, &f);      // 1/pt pulls
 988   fNRefFigures = 21;
 989   return kTRUE;
 990
 991   // TRACK HMPID RESOLUTION/PULLS
 992   Process(kMCtrackHMPID, 0, &f, 1.e4); // z towards TOF
 993   Process(kMCtrackHMPID, 1, &f);       // z towards TOF
 994   fNRefFigures = 22;
 995
 996   return kTRUE;
 997 }
 998
 999
1000 //________________________________________________________
1001 void AliTRDresolution::Terminate(Option_t *)
1002 {
1003   if(fDebugStream){
1004     delete fDebugStream;
1005     fDebugStream = 0x0;
1006     fDebugLevel = 0;
1007   }
1008   if(HasPostProcess()) PostProcess();
1009 }
1010
1011 //________________________________________________________
1012 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1013 {
1014 // Helper function to avoid duplication of code
1015 // Make first guesses on the fit parameters
1016
1017   // find the intial parameters of the fit !! (thanks George)
1018   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1019   Double_t sum = 0.;
1020   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1021   f->SetParLimits(0, 0., 3.*sum);
1022   f->SetParameter(0, .9*sum);
1023
1024   f->SetParLimits(1, -.2, .2);
1025   f->SetParameter(1, -0.1);
1026
1027   f->SetParLimits(2, 0., 4.e-1);
1028   f->SetParameter(2, 2.e-2);
1029   if(f->GetNpar() <= 4) return;
1030
1031   f->SetParLimits(3, 0., sum);
1032   f->SetParameter(3, .1*sum);
1033
1034   f->SetParLimits(4, -.3, .3);
1035   f->SetParameter(4, 0.);
1036
1037   f->SetParLimits(5, 0., 1.e2);
1038   f->SetParameter(5, 2.e-1);
1039 }
1040
1041 //________________________________________________________
1042 TObjArray* AliTRDresolution::Histos()
1043 {
1044   if(fContainer) return fContainer;
1045
1046   fContainer  = new TObjArray(kNhistos);
1047   //fContainer->SetOwner(kTRUE);
1048   TH1 *h = 0x0;
1049   TObjArray *arr = 0x0;
1050
1051   // cluster to track residuals/pulls
1052   fContainer->AddAt(arr = new TObjArray(fNElements[kCluster]), kCluster);
1053   arr->SetName("Cl");
1054   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
1055     h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
1056     h->GetXaxis()->SetTitle("tg(#phi)");
1057     h->GetYaxis()->SetTitle("#Delta y [cm]");
1058     h->GetZaxis()->SetTitle("entries");
1059   } else h->Reset();
1060   arr->AddAt(h, 0);
1061   if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
1062     h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1063     h->GetXaxis()->SetTitle("tg(#phi)");
1064     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1065     h->GetZaxis()->SetTitle("entries");
1066   } else h->Reset();
1067   arr->AddAt(h, 1);
1068
1069   // tracklet to track residuals/pulls in y direction
1070   fContainer->AddAt(arr = new TObjArray(fNElements[kTracklet]), kTracklet);
1071   arr->SetName("Trklt");
1072   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
1073     h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
1074     h->GetXaxis()->SetTitle("#tg(#phi)");
1075     h->GetYaxis()->SetTitle("#Delta y [cm]");
1076     h->GetZaxis()->SetTitle("entries");
1077   } else h->Reset();
1078   arr->AddAt(h, 0);
1079   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
1080     h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1081     h->GetXaxis()->SetTitle("#tg(#phi)");
1082     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1083     h->GetZaxis()->SetTitle("entries");
1084   } else h->Reset();
1085   arr->AddAt(h, 1);
1086   // tracklet to track residuals/pulls in z direction
1087   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
1088     h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 21, -1., 1., 100, -1.5, 1.5);
1089     h->GetXaxis()->SetTitle("#tg(#theta)");
1090     h->GetYaxis()->SetTitle("#Delta z [cm]");
1091     h->GetZaxis()->SetTitle("entries");
1092   } else h->Reset();
1093   arr->AddAt(h, 2);
1094   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
1095     h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 21, -1., 1., 100, -5.5, 5.5);
1096     h->GetXaxis()->SetTitle("#tg(#theta)");
1097     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1098     h->GetZaxis()->SetTitle("entries");
1099   } else h->Reset();
1100   arr->AddAt(h, 3);
1101   // tracklet to track phi residuals
1102   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
1103     h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1104     h->GetXaxis()->SetTitle("tg(#phi)");
1105     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1106     h->GetZaxis()->SetTitle("entries");
1107   } else h->Reset();
1108   arr->AddAt(h, 4);
1109
1110
1111   // tracklet to TPC track residuals/pulls in y direction
1112   fContainer->AddAt(arr = new TObjArray(fNElements[kTrackTPC]), kTrackTPC);
1113   arr->SetName("TrkTPC");
1114   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
1115     h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
1116     h->GetXaxis()->SetTitle("#tg(#phi)");
1117     h->GetYaxis()->SetTitle("#Delta y [cm]");
1118     h->GetZaxis()->SetTitle("entries");
1119   } else h->Reset();
1120   arr->AddAt(h, 0);
1121   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
1122     h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1123     h->GetXaxis()->SetTitle("#tg(#phi)");
1124     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1125     h->GetZaxis()->SetTitle("entries");
1126   } else h->Reset();
1127   arr->AddAt(h, 1);
1128   // tracklet to TPC track residuals/pulls in z direction
1129   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
1130     h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 21, -1., 1., 100, -1.5, 1.5);
1131     h->GetXaxis()->SetTitle("#tg(#theta)");
1132     h->GetYaxis()->SetTitle("#Delta z [cm]");
1133     h->GetZaxis()->SetTitle("entries");
1134   } else h->Reset();
1135   arr->AddAt(h, 2);
1136   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
1137     h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 21, -1., 1., 100, -5.5, 5.5);
1138     h->GetXaxis()->SetTitle("#tg(#theta)");
1139     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1140     h->GetZaxis()->SetTitle("entries");
1141   } else h->Reset();
1142   arr->AddAt(h, 3);
1143   // tracklet to TPC track phi residuals
1144   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
1145     h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1146     h->GetXaxis()->SetTitle("tg(#phi)");
1147     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1148     h->GetZaxis()->SetTitle("entries");
1149   } else h->Reset();
1150   arr->AddAt(h, 4);
1151
1152
1153   // Resolution histos
1154   if(!HasMCdata()) return fContainer;
1155
1156   // cluster y resolution [0]
1157   fContainer->AddAt(arr = new TObjArray(fNElements[kMCcluster]), kMCcluster);
1158   arr->SetName("McCl");
1159   if(!(h = (TH2I*)gROOT->FindObject("hMCcl"))){
1160     h = new TH2I("hMCcl", "Cluster Resolution", 31, -.48, .48, 100, -.3, .3);
1161     h->GetXaxis()->SetTitle("tg(#phi)");
1162     h->GetYaxis()->SetTitle("#Delta y [cm]");
1163     h->GetZaxis()->SetTitle("entries");
1164   } else h->Reset();
1165   arr->AddAt(h, 0);
1166   if(!(h = (TH2I*)gROOT->FindObject("hMCclPull"))){
1167     h = new TH2I("hMCclPull", "Cluster Pulls", 31, -.48, .48, 100, -4.5, 4.5);
1168     h->GetXaxis()->SetTitle("tg(#phi)");
1169     h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
1170     h->GetZaxis()->SetTitle("entries");
1171   } else h->Reset();
1172   arr->AddAt(h, 1);
1173
1174
1175   // TRACKLET RESOLUTION
1176   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtracklet]), kMCtracklet);
1177   arr->SetName("McTrklt");
1178   // tracklet y resolution
1179   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltY"))){
1180     h = new TH2I("hMCtrkltY", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.2, .2);
1181     h->GetXaxis()->SetTitle("tg(#phi)");
1182     h->GetYaxis()->SetTitle("#Delta y [cm]");
1183     h->GetZaxis()->SetTitle("entries");
1184   } else h->Reset();
1185   arr->AddAt(h, 0);
1186   // tracklet y pulls
1187   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltYPull"))){
1188     h = new TH2I("hMCtrkltYPull", "Tracklet Pulls (Y)", 31, -.48, .48, 100, -4.5, 4.5);
1189     h->GetXaxis()->SetTitle("tg(#phi)");
1190     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1191     h->GetZaxis()->SetTitle("entries");
1192   } else h->Reset();
1193   arr->AddAt(h, 1);
1194   // tracklet z resolution
1195   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZ"))){
1196     h = new TH2I("hMCtrkltZ", "Tracklet Resolution (Z)", 50, -1., 1., 100, -1., 1.);
1197     h->GetXaxis()->SetTitle("tg(#theta)");
1198     h->GetYaxis()->SetTitle("#Delta z [cm]");
1199     h->GetZaxis()->SetTitle("entries");
1200   } else h->Reset();
1201   arr->AddAt(h, 2);
1202   // tracklet z pulls
1203   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZPull"))){
1204     h = new TH2I("hMCtrkltZPull", "Tracklet Pulls (Z)", 31, -1., 1., 100, -3.5, 3.5);
1205     h->GetXaxis()->SetTitle("tg(#theta)");
1206     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1207     h->GetZaxis()->SetTitle("entries");
1208   } else h->Reset();
1209   arr->AddAt(h, 3);
1210   // tracklet phi resolution
1211   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltPhi"))){
1212     h = new TH2I("hMCtrkltPhi", "Tracklet Resolution (#Phi)", 31, -.48, .48, 100, -.15, .15);
1213     h->GetXaxis()->SetTitle("tg(#phi)");
1214     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1215     h->GetZaxis()->SetTitle("entries");
1216   } else h->Reset();
1217   arr->AddAt(h, 4);
1218
1219
1220   // KALMAN TRACK RESOLUTION
1221   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrack]), kMCtrack);
1222   arr->SetName("McTrack");
1223   // Kalman track y resolution
1224   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkY"))){
1225     h = new TH2I("hMCtrkY", "Track Y Resolution", 31, -.48, .48, 100, -.2, .2);
1226     h->GetXaxis()->SetTitle("tg(#phi)");
1227     h->GetYaxis()->SetTitle("#Delta y [cm]");
1228     h->GetZaxis()->SetTitle("entries");
1229   } else h->Reset();
1230   arr->AddAt(h, 0);
1231   // Kalman track y pulls
1232   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYPull"))){
1233     h = new TH2I("hMCtrkYPull", "Track Y Pulls", 31, -.48, .48, 100, -4., 4.);
1234     h->GetXaxis()->SetTitle("tg(#phi)");
1235     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1236     h->GetZaxis()->SetTitle("entries");
1237   } else h->Reset();
1238   arr->AddAt(h, 1);
1239   // Kalman track Z
1240   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZ"))){
1241     h = new TH2I("hMCtrkZ", "Track Z Resolution", 20, -1., 1., 100, -1., 1.);
1242     h->GetXaxis()->SetTitle("tg(#theta)");
1243     h->GetYaxis()->SetTitle("#Delta z [cm]");
1244     h->GetZaxis()->SetTitle("entries");
1245   } else h->Reset();
1246   arr->AddAt(h, 2);
1247   // Kalman track Z pulls
1248   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZPull"))){
1249     h = new TH2I("hMCtrkZPull", "Track Z Pulls", 20, -1., 1., 100, -4.5, 4.5);
1250     h->GetXaxis()->SetTitle("tg(#theta)");
1251     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1252     h->GetZaxis()->SetTitle("entries");
1253   } else h->Reset();
1254   arr->AddAt(h, 3);
1255   // Kalman track SNP
1256   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkSNP"))){
1257     h = new TH2I("hMCtrkSNP", "Track Phi Resolution", 20, -.3, .3, 100, -.02, .02);
1258     h->GetXaxis()->SetTitle("tg(#phi)");
1259     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1260     h->GetZaxis()->SetTitle("entries");
1261   } else h->Reset();
1262   arr->AddAt(h, 4);
1263   // Kalman track SNP pulls
1264   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkSNPPull"))){
1265     h = new TH2I("hMCtrkSNPPull", "Track SNP Pulls", 20, -.3, .3, 100, -4.5, 4.5);
1266     h->GetXaxis()->SetTitle("tg(#phi)");
1267     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1268     h->GetZaxis()->SetTitle("entries");
1269   } else h->Reset();
1270   arr->AddAt(h, 5);
1271   // Kalman track TGL
1272   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkTGL"))){
1273     h = new TH2I("hMCtrkTGL", "Track Theta Resolution", 20, -1., 1., 100, -.1, .1);
1274     h->GetXaxis()->SetTitle("tg(#theta)");
1275     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1276     h->GetZaxis()->SetTitle("entries");
1277   } else h->Reset();
1278   arr->AddAt(h, 6);
1279   // Kalman track TGL pulls
1280   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkTGLPull"))){
1281     h = new TH2I("hMCtrkTGLPull", "Track TGL  Pulls", 20, -1., 1., 100, -4.5, 4.5);
1282     h->GetXaxis()->SetTitle("tg(#theta)");
1283     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1284     h->GetZaxis()->SetTitle("entries");
1285   } else h->Reset();
1286   arr->AddAt(h, 7);
1287   // Kalman track Pt resolution
1288   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPt"))){
1289     h = new TH2I("hMCtrkPt", "Track Pt Resolution", 40, 0., 20., 150, -.15, .15);
1290     h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1291     h->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1292     h->GetZaxis()->SetTitle("entries");
1293   } else h->Reset();
1294   arr->AddAt(h, 8);
1295   // Kalman track Pt pulls
1296   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPtPulls"))){
1297     h = new TH2I("hMCtrkPtPulls", "Track 1/Pt Pulls", 40, 0., 2., 100, -4., 4.);
1298     h->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1299     h->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1300     h->GetZaxis()->SetTitle("entries");
1301   } else h->Reset();
1302   arr->AddAt(h, 9);
1303
1304
1305   // TPC TRACK RESOLUTION
1306   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTPC]), kMCtrackTPC);
1307   arr->SetName("McTrackTPC");
1308   // Kalman track Y
1309   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYIn"))){
1310     h = new TH2I("hMCtrkYIn", "Track[TPC] Y Resolution", 20, -.3, .3, 100, -.5, .5);
1311     h->GetXaxis()->SetTitle("tg(#phi)");
1312     h->GetYaxis()->SetTitle("#Delta y [cm]");
1313     h->GetZaxis()->SetTitle("entries");
1314   } else h->Reset();
1315   arr->AddAt(h, 0);
1316   // Kalman track Y pulls
1317   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYInPull"))){
1318     h = new TH2I("hMCtrkYInPull", "Track[TPC] Y Pulls", 20, -.3, .3, 100, -4.5, 4.5);
1319     h->GetXaxis()->SetTitle("tg(#phi)");
1320     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1321     h->GetZaxis()->SetTitle("entries");
1322   } else h->Reset();
1323   arr->AddAt(h, 1);
1324   // Kalman track Z
1325   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZIn"))){
1326     h = new TH2I("hMCtrkZIn", "Track[TPC] Z Resolution", 20, -1., 1., 100, -1., 1.);
1327     h->GetXaxis()->SetTitle("tg(#theta)");
1328     h->GetYaxis()->SetTitle("#Delta z [cm]");
1329     h->GetZaxis()->SetTitle("entries");
1330   } else h->Reset();
1331   arr->AddAt(h, 2);
1332   // Kalman track Z pulls
1333   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZInPull"))){
1334     h = new TH2I("hMCtrkZInPull", "Track[TPC] Z Pulls", 20, -1., 1., 100, -4.5, 4.5);
1335     h->GetXaxis()->SetTitle("tg(#theta)");
1336     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1337     h->GetZaxis()->SetTitle("entries");
1338   } else h->Reset();
1339   arr->AddAt(h, 3);
1340   // Kalman track SNP
1341   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkSNPIn"))){
1342     h = new TH2I("hMCtrkSNPIn", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -.02, .02);
1343     h->GetXaxis()->SetTitle("tg(#phi)");
1344     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1345     h->GetZaxis()->SetTitle("entries");
1346   } else h->Reset();
1347   arr->AddAt(h, 4);
1348   // Kalman track SNP pulls
1349   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkSNPInPull"))){
1350     h = new TH2I("hMCtrkSNPInPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1351     h->GetXaxis()->SetTitle("tg(#phi)");
1352     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1353     h->GetZaxis()->SetTitle("entries");
1354   } else h->Reset();
1355   arr->AddAt(h, 5);
1356   // Kalman track TGL
1357   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkTGLIn"))){
1358     h = new TH2I("hMCtrkTGLIn", "Track[TPC] Theta Resolution", 60, -1., 1., 100, -.1, .1);
1359     h->GetXaxis()->SetTitle("tg(#theta)");
1360     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1361     h->GetZaxis()->SetTitle("entries");
1362   } else h->Reset();
1363   arr->AddAt(h, 6);
1364   // Kalman track TGL pulls
1365   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkTGLInPull"))){
1366     h = new TH2I("hMCtrkTGLInPull", "Track[TPC] TGL  Pulls", 60, -1., 1., 100, -4.5, 4.5);
1367     h->GetXaxis()->SetTitle("tg(#theta)");
1368     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1369     h->GetZaxis()->SetTitle("entries");
1370   } else h->Reset();
1371   arr->AddAt(h, 7);
1372   // Kalman track Pt resolution
1373   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPtIn"))){
1374     h = new TH2I("hMCtrkPtIn", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.15, .15);
1375     h->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1376     h->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1377     h->GetZaxis()->SetTitle("entries");
1378   } else h->Reset();
1379   arr->AddAt(h, 8);
1380   // Kalman track Pt pulls
1381   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPtInPulls"))){
1382     h = new TH2I("hMCtrkPtInPulls", "Track[TPC] 1/Pt Pulls", 80, 0., 2., 100, -4., 4.);
1383     h->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1384     h->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1385     h->GetZaxis()->SetTitle("entries");
1386   } else h->Reset();
1387   arr->AddAt(h, 9);
1388
1389
1390
1391   // Kalman track Z resolution [OUT]
1392   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackHMPID]), kMCtrackHMPID);
1393   arr->SetName("McTrackHMPID");
1394   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOut"))){
1395     h = new TH2I("hMCtrkZOut", "Track[TOF] Z Resolution", 20, -1., 1., 100, -1., 1.);
1396     h->GetXaxis()->SetTitle("tg(#theta)");
1397     h->GetYaxis()->SetTitle("#Delta z [cm]");
1398     h->GetZaxis()->SetTitle("entries");
1399   } else h->Reset();
1400   arr->AddAt(h, 0);
1401   // Kalman track Z pulls
1402   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOutPull"))){
1403     h = new TH2I("hMCtrkZOutPull", "Track[TOF] Z Pulls", 20, -1., 1., 100, -4.5, 4.5);
1404     h->GetXaxis()->SetTitle("tg(#theta)");
1405     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1406     h->GetZaxis()->SetTitle("entries");
1407   } else h->Reset();
1408   arr->AddAt(h, 1);
1409
1410   return fContainer;
1411 }
1412
1413
1414 //________________________________________________________
1415 Bool_t AliTRDresolution::Process(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1416 {
1417   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1418   Bool_t kBUILD = kFALSE;
1419   if(!f){
1420     f = new TF1("f1", "gaus", -.5, .5);
1421     kBUILD = kTRUE;
1422   }
1423
1424   // retrive containers
1425   TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1426   if(!h2) return kFALSE;
1427   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(idx);
1428   if(!gm) return kFALSE;
1429   if(gm->GetN()) for(Int_t ip=gm->GetN(); ip--;) gm->RemovePoint(ip);
1430
1431   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(idx);
1432   if(!gs) return kFALSE;
1433   if(gs->GetN()) for(Int_t ip=gs->GetN(); ip--;) gs->RemovePoint(ip);
1434
1435   Char_t pn[10]; sprintf(pn, "p%02d%02d", plot, idx);
1436   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1437     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1438     TH1D *h = h2->ProjectionY(pn, ibin, ibin);
1439     if(h->GetEntries()<100) continue;
1440     AdjustF1(h, f);
1441
1442     h->Fit(f, "QN");
1443
1444     Int_t ip = gm->GetN();
1445     gm->SetPoint(ip, x, k*f->GetParameter(1));
1446     gm->SetPointError(ip, 0., k*f->GetParError(1));
1447     gs->SetPoint(ip, x, k*f->GetParameter(2));
1448     gs->SetPointError(ip, 0., k*f->GetParError(2));
1449   }
1450
1451   if(kBUILD) delete f;
1452   return kTRUE;
1453 }
1454
1455 //________________________________________________________
1456 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, TF1 *f, Float_t k)
1457 {
1458   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1459   Bool_t kBUILD = kFALSE;
1460   if(!f){
1461     f = new TF1("f1", "gaus", -.5, .5);
1462     kBUILD = kTRUE;
1463   }
1464
1465   TH3I *h3 = 0x0;
1466   if(!(h3 = (TH3I*)(fContainer->At(plot)))) return kFALSE;
1467   TGraphErrors *gm = 0x0, *gs = 0x0;
1468   if(!(gm=(TGraphErrors*)fGraphM->At(plot))) return kFALSE;
1469   if(gm->GetN()) for(Int_t ip=gm->GetN(); ip--;) gm->RemovePoint(ip);
1470   if(!(gs=(TGraphErrors*)fGraphS->At(plot))) return kFALSE;
1471   if(gs->GetN()) for(Int_t ip=gs->GetN(); ip--;) gs->RemovePoint(ip);
1472   Char_t pn[10]; sprintf(pn, "p%02d", plot);
1473   for(Int_t ibin = 1; ibin <= h3->GetNbinsX(); ibin++){
1474     Double_t x = h3->GetXaxis()->GetBinCenter(ibin);
1475     TH1D *h = h3->ProjectionZ(pn, ibin, ibin);
1476     if(h->GetEntries()<100) continue;
1477     AdjustF1(h, f);
1478
1479     h->Fit(f, "QN");
1480
1481     Int_t ip = gm->GetN();
1482     gm->SetPoint(ip, x, k*f->GetParameter(1));
1483     gm->SetPointError(ip, 0., k*f->GetParError(1));
1484     gs->SetPoint(ip, x, k*f->GetParameter(2));
1485     gs->SetPointError(ip, 0., k*f->GetParError(2));
1486   }
1487
1488   if(kBUILD) delete f;
1489   return kTRUE;
1490 }
1491
1492 //________________________________________________________
1493 Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
1494 {
1495   if(!fGraphS || !fGraphM) return kFALSE;
1496   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
1497   if(!gm) return kFALSE;
1498   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
1499   if(!gs) return kFALSE;
1500   gs->Draw("apl"); gm->Draw("pl");
1501
1502   //printf("bb[%f %f %f %f]\n", bb[0], bb[1], bb[2], bb[3]);
1503
1504   // axis range
1505   TAxis *ax = 0x0;
1506   ax = gs->GetHistogram()->GetYaxis();
1507   ax->SetRangeUser(bb[1], bb[3]);
1508   ax = gs->GetHistogram()->GetXaxis();
1509   ax->SetRangeUser(bb[0], bb[2]);
1510
1511   // axis titles
1512   Int_t nref = 0;
1513   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fNElements[jp];
1514   UChar_t jdx = idx<0?0:idx;
1515   for(Int_t jc=0; jc<TMath::Min(jdx,fNElements[ip]-1); jc++) nref++;
1516   Char_t **at = fAxTitle[nref];
1517   ax->SetTitle(*at);ax->CenterTitle();
1518   TGaxis *gax = 0x0;
1519   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+L");
1520   gax->SetLineColor(kBlue);gax->SetLineWidth(2);gax->SetTextColor(kBlue);
1521   //gax->SetVertical();
1522   gax->CenterTitle();
1523   gax->SetTitle(*(++at)); gax->Draw();
1524
1525   gax = new TGaxis(bb[0], bb[1], bb[0], bb[3], bb[1], bb[3], 510, "");
1526   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1527   /*gax->SetVertical();*/gax->CenterTitle();
1528   gax->SetTitle(*(++at)); gax->Draw();
1529   //return kTRUE;
1530
1531   // bounding box
1532   TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
1533   b->SetFillStyle(3002);b->SetLineColor(0);
1534   b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
1535   b->Draw();
1536   return kTRUE;
1537 }
1538
1539 // gSystem->Load("libANALYSIS.so")
1540 // gSystem->Load("libTRDqaRec.so")
1541 // AliTRDresolution res
1542 void AliTRDresolution::DumpAxTitle(Int_t ip, Int_t idx)
1543 {
1544   Int_t nref = 0;
1545   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fNElements[jp];
1546   UChar_t jdx = idx<0?0:idx;
1547   printf("ref entry* %d jdx[%d]\n", nref, jdx);
1548   for(Int_t jc=0; jc<TMath::Min(jdx,fNElements[ip]-1); jc++) nref++;
1549   printf("ref entry %d\n", nref);
1550   Char_t **at = fAxTitle[nref];
1551
1552   printf("%s || ", (*at));
1553   at++; printf("%s || ", (*at));
1554   at++; printf("%s\n", (*at));
1555 }
1556
1557
1558 //________________________________________________________
1559 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
1560 {
1561
1562   fReconstructor->SetRecoParam(r);
1563 }