PWG1/TRD/AliTRDresolution.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercialf purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  TRD tracking resolution                                               //
  21 //
  22 // The class performs resolution and residual studies
  23 // of the TRD tracks for the following quantities :
  24 //   - spatial position (y, [z])
  25 //   - angular (phi) tracklet
  26 //   - momentum at the track level
  27 //
  28 // The class has to be used for regular detector performance checks using the official macros:
  29 //   - $ALICE_ROOT/TRD/qaRec/run.C
  30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
  31 //
  32 // For stand alone usage please refer to the following example:
  33 // {
  34 //   gSystem->Load("libANALYSIS.so");
  35 //   gSystem->Load("libTRDqaRec.so");
  36 //   AliTRDresolution *res = new AliTRDresolution();
  37 //   //res->SetMCdata();
  38 //   //res->SetVerbose();
  39 //   //res->SetVisual();
  40 //   res->Load();
  41 //   if(!res->PostProcess()) return;
  42 //   res->GetRefFigure(0);
  43 // }
  44 //
  45 //  Authors:                                                              //
  46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
  47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  48 //                                                                        //
  49 ////////////////////////////////////////////////////////////////////////////
  50
  51 #include <TSystem.h>
  52
  53 #include <TROOT.h>
  54 #include <TObjArray.h>
  55 #include <TH3.h>
  56 #include <TH2.h>
  57 #include <TH1.h>
  58 #include <TF1.h>
  59 #include <TCanvas.h>
  60 #include <TGaxis.h>
  61 #include <TBox.h>
  62 #include <TLegend.h>
  63 #include <TGraphErrors.h>
  64 #include <TGraphAsymmErrors.h>
  65 #include <TLinearFitter.h>
  66 #include <TMath.h>
  67 #include <TMatrixT.h>
  68 #include <TVectorT.h>
  69 #include <TTreeStream.h>
  70 #include <TGeoManager.h>
  71 #include <TDatabasePDG.h>
  72
  73 #include "AliPID.h"
  74 #include "AliLog.h"
  75 #include "AliESDtrack.h"
  76 #include "AliMathBase.h"
  77 #include "AliTrackPointArray.h"
  78
  79 #include "AliTRDresolution.h"
  80 #include "AliTRDgeometry.h"
  81 #include "AliTRDpadPlane.h"
  82 #include "AliTRDcluster.h"
  83 #include "AliTRDseedV1.h"
  84 #include "AliTRDtrackV1.h"
  85 #include "AliTRDReconstructor.h"
  86 #include "AliTRDrecoParam.h"
  87 #include "AliTRDpidUtil.h"
  88 #include "AliTRDinfoGen.h"
  89
  90 #include "info/AliTRDclusterInfo.h"
  91
  92 ClassImp(AliTRDresolution)
  93
  94 UChar_t const AliTRDresolution::fgNproj[kNviews] = {
  95   2, 2, 5, 5, 5,
  96   2, 5, 11, 11, 11
  97 };
  98 Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
  99      "Charge"
 100     ,"Cluster2Track"
 101     ,"Tracklet2Track"
 102     ,"Tracklet2TRDin"
 103     ,"Tracklet2TRDout"
 104     ,"Cluster2MC"
 105     ,"Tracklet2MC"
 106     ,"TRDin2MC"
 107     ,"TRDout2MC"
 108     ,"TRD2MC"
 109 };
 110 Char_t const * AliTRDresolution::fgParticle[11]={
 111   " p bar", " K -", " #pi -", " #mu -", " e -",
 112   " No PID",
 113   " e +", " #mu +", " #pi +", " K +", " p",
 114 };
 115
 116 // Configure segmentation for y resolution/residuals
 117 Int_t const AliTRDresolution::fgkNresYsegm[3] = {
 118   AliTRDgeometry::kNsector
 119  ,AliTRDgeometry::kNsector*AliTRDgeometry::kNstack
 120  ,AliTRDgeometry::kNdet
 121 };
 122 Char_t const *AliTRDresolution::fgkResYsegmName[3] = {
 123   "Sector", "Stack", "Detector"};
 124
 125
 126 //________________________________________________________
 127 AliTRDresolution::AliTRDresolution()
 128   :AliTRDrecoTask()
 129   ,fSegmentLevel(0)
 130   ,fIdxPlot(0)
 131   ,fIdxFrame(0)
 132   ,fPtThreshold(1.)
 133   ,fDyRange(1.5)
 134   ,fDBPDG(NULL)
 135   ,fGraphS(NULL)
 136   ,fGraphM(NULL)
 137   ,fCl(NULL)
 138   ,fMCcl(NULL)
 139 /*  ,fTrklt(NULL)
 140   ,fMCtrklt(NULL)*/
 141 {
 142   //
 143   // Default constructor
 144   //
 145   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
 146   SetSegmentationLevel();
 147   memset(fXcorr, 0, AliTRDgeometry::kNdet*30*sizeof(Float_t));
 148 }
 149
 150 //________________________________________________________
 151 AliTRDresolution::AliTRDresolution(char* name)
 152   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
 153   ,fSegmentLevel(0)
 154   ,fIdxPlot(0)
 155   ,fIdxFrame(0)
 156   ,fPtThreshold(1.)
 157   ,fDyRange(1.5)
 158   ,fDBPDG(NULL)
 159   ,fGraphS(NULL)
 160   ,fGraphM(NULL)
 161   ,fCl(NULL)
 162   ,fMCcl(NULL)
 163 /*  ,fTrklt(NULL)
 164   ,fMCtrklt(NULL)*/
 165 {
 166   //
 167   // Default constructor
 168   //
 169
 170   InitFunctorList();
 171   SetSegmentationLevel();
 172   memset(fXcorr, 0, AliTRDgeometry::kNdet*30*sizeof(Float_t));
 173
 174   DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
 175   DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
 176 /*  DefineOutput(kTrkltToTrk, TObjArray::Class()); // tracklet2track
 177   DefineOutput(kTrkltToMC, TObjArray::Class()); // tracklet2mc*/
 178 }
 179
 180 //________________________________________________________
 181 AliTRDresolution::~AliTRDresolution()
 182 {
 183   //
 184   // Destructor
 185   //
 186
 187   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 188   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 189   if(fCl){fCl->Delete(); delete fCl;}
 190   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 191 /*  if(fTrklt){fTrklt->Delete(); delete fTrklt;}
 192   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}*/
 193 }
 194
 195
 196 //________________________________________________________
 197 void AliTRDresolution::UserCreateOutputObjects()
 198 {
 199   // spatial resolution
 200
 201   AliTRDrecoTask::UserCreateOutputObjects();
 202   InitExchangeContainers();
 203 }
 204
 205 //________________________________________________________
 206 void AliTRDresolution::InitExchangeContainers()
 207 {
 208 // Init containers for subsequent tasks (AliTRDclusterResolution)
 209
 210   fCl = new TObjArray(200);
 211   fCl->SetOwner(kTRUE);
 212   fMCcl = new TObjArray();
 213   fMCcl->SetOwner(kTRUE);
 214 /*  fTrklt = new TObjArray();
 215   fTrklt->SetOwner(kTRUE);
 216   fMCtrklt = new TObjArray();
 217   fMCtrklt->SetOwner(kTRUE);*/
 218   PostData(kClToTrk, fCl);
 219   PostData(kClToMC, fMCcl);
 220 /*  PostData(kTrkltToTrk, fTrklt);
 221   PostData(kTrkltToMC, fMCtrklt);*/
 222 }
 223
 224 //________________________________________________________
 225 void AliTRDresolution::UserExec(Option_t *opt)
 226 {
 227   //
 228   // Execution part
 229   //
 230
 231   fCl->Delete();
 232   fMCcl->Delete();
 233   if(!HasLoadCorrection()) LoadCorrection("correction.lst");
 234   AliTRDrecoTask::UserExec(opt);
 235 }
 236
 237 //________________________________________________________
 238 Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt) const
 239 {
 240 // Helper function to calculate pulls in the yz plane
 241 // using proper tilt rotation
 242 // Uses functionality defined by AliTRDseedV1.
 243
 244   Double_t t2(tilt*tilt);
 245
 246   // rotate along pad
 247   Double_t cc[3];
 248   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
 249   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
 250   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
 251   // do sqrt
 252   Double_t sqr[3]={0., 0., 0.};
 253   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
 254   Double_t invsqr[3]={0., 0., 0.};
 255   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
 256   Double_t tmp(dyz[0]);
 257   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
 258   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
 259   return kTRUE;
 260 }
 261
 262 //________________________________________________________
 263 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
 264 {
 265   //
 266   // Plots the charge distribution
 267   //
 268
 269   if(track) fkTrack = track;
 270   if(!fkTrack){
 271     AliDebug(4, "No Track defined.");
 272     return NULL;
 273   }
 274   TObjArray *arr = NULL;
 275   if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCharge)))){
 276     AliWarning("No output container defined.");
 277     return NULL;
 278   }
 279   TH3S* h = NULL;
 280
 281   AliTRDseedV1 *fTracklet = NULL;
 282   AliTRDcluster *c = NULL;
 283   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 284     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 285     if(!fTracklet->IsOK()) continue;
 286     Float_t x0 = fTracklet->GetX0();
 287     Float_t dqdl, dl;
 288     for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
 289       if(!(c = fTracklet->GetClusters(itb))){
 290         if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
 291       }
 292       dqdl = fTracklet->GetdQdl(itb, &dl);
 293       if(dqdl<1.e-5) continue;
 294       dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
 295       (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dqdl);
 296     }
 297
 298 //     if(!HasMCdata()) continue;
 299 //     UChar_t s;
 300 //     Float_t pt0, y0, z0, dydx0, dzdx0;
 301 //     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
 302
 303   }
 304   return h;
 305 }
 306
 307
 308 //________________________________________________________
 309 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 310 {
 311   //
 312   // Plot the cluster distributions
 313   //
 314
 315   if(track) fkTrack = track;
 316   if(!fkTrack){
 317     AliDebug(4, "No Track defined.");
 318     return NULL;
 319   }
 320   TObjArray *arr = NULL;
 321   if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCluster)))){
 322     AliWarning("No output container defined.");
 323     return NULL;
 324   }
 325   ULong_t status = fkESD ? fkESD->GetStatus():0;
 326
 327   Int_t sgm[3];
 328   Double_t covR[7], cov[3], dy[2], dz[2];
 329   Float_t pt, x0, y0, z0, dydx, dzdx, tilt(0.);
 330   const AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 331   AliTRDseedV1 *fTracklet(NULL);  TObjArray *clInfoArr(NULL);
 332   // do LINEAR track refit if asked by the user
 333   // it is the user responsibility to check if B=0T
 334   Float_t param[10];  memset(param, 0, 10*sizeof(Float_t));
 335   Int_t np(0), nrc(0); AliTrackPoint clusters[300];
 336   if(HasTrackRefit()){
 337     Bool_t kPrimary(kFALSE);
 338     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 339       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 340       if(!fTracklet->IsOK()) continue;
 341       x0   = fTracklet->GetX0();
 342       tilt = fTracklet->GetTilt();
 343       AliTRDcluster *c = NULL;
 344       fTracklet->ResetClusterIter(kFALSE);
 345       while((c = fTracklet->PrevCluster())){
 346         Float_t xyz[3] = {c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()], c->GetY(), c->GetZ()};
 347         clusters[np].SetCharge(tilt);
 348         clusters[np].SetClusterType(0);
 349         clusters[np].SetVolumeID(ily);
 350         clusters[np].SetXYZ(xyz);
 351         np++;
 352       }
 353       if(fTracklet->IsRowCross()){
 354         Float_t xcross(0.), zcross(0.);
 355         if(fTracklet->GetEstimatedCrossPoint(xcross, zcross)){
 356           clusters[np].SetCharge(tilt);
 357           clusters[np].SetClusterType(1);
 358           clusters[np].SetVolumeID(ily);
 359           clusters[np].SetXYZ(xcross, 0., zcross);
 360           np++;
 361           nrc++;
 362         }
 363       }
 364       if(fTracklet->IsPrimary()) kPrimary = kTRUE;
 365     }
 366     if(kPrimary){
 367       clusters[np].SetCharge(tilt);
 368       clusters[np].SetClusterType(1);
 369       clusters[np].SetVolumeID(-1);
 370       clusters[np].SetXYZ(0., 0., 0.);
 371       np++;
 372     }
 373     if(!FitTrack(np, clusters, param)) {
 374       AliDebug(1, "Linear track Fit failed.");
 375       return NULL;
 376     }
 377     if(HasTrackSelection()){
 378       Bool_t kReject(kFALSE);
 379       if(fkTrack->GetNumberOfTracklets() != AliTRDgeometry::kNlayer)  kReject = kTRUE;
 380       if(!kReject && !UseTrack(np, clusters, param)) kReject = kTRUE;
 381       if(kReject){
 382         AliDebug(1, "Reject track for residuals analysis.");
 383         return NULL;
 384       }
 385     }
 386   }
 387   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 388     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 389     if(!fTracklet->IsOK()) continue;
 390     x0 = fTracklet->GetX0();
 391     pt = fTracklet->GetPt();
 392     sgm[2] = fTracklet->GetDetector();
 393     sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 394     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 395
 396     // retrive the track angle with the chamber
 397     if(HasTrackRefit()){
 398       Float_t par[3];
 399       if(!FitTracklet(ily, np, clusters, param, par)) continue;
 400       dydx = par[2];//param[3];
 401       dzdx = param[4];
 402       y0   = par[1] + dydx * (x0 - par[0]);//param[1] + dydx * (x0 - param[0]);
 403       z0   = param[2] + dzdx * (x0 - param[0]);
 404     } else {
 405       y0   = fTracklet->GetYref(0);
 406       z0   = fTracklet->GetZref(0);
 407       dydx = fTracklet->GetYref(1);
 408       dzdx = fTracklet->GetZref(1);
 409     }
 410     /*printf("RC[%c] Primary[%c]\n"
 411            "  Fit : y0[%f] z0[%f] dydx[%f] dzdx[%f]\n"
 412            "  Ref:  y0[%f] z0[%f] dydx[%f] dzdx[%f]\n", fTracklet->IsRowCross()?'y':'n', fTracklet->IsPrimary()?'y':'n', y0, z0, dydx, dzdx
 413            ,fTracklet->GetYref(0),fTracklet->GetZref(0),fTracklet->GetYref(1),fTracklet->GetZref(1));*/
 414     tilt = fTracklet->GetTilt();
 415     fTracklet->GetCovRef(covR);
 416     Double_t t2(tilt*tilt)
 417             ,corr(1./(1. + t2))
 418             ,cost(TMath::Sqrt(corr));
 419     AliTRDcluster *c = NULL;
 420     fTracklet->ResetClusterIter(kFALSE);
 421     while((c = fTracklet->PrevCluster())){
 422       Float_t xc = c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
 423       Float_t yc = c->GetY();
 424       Float_t zc = c->GetZ();
 425       Float_t dx = x0 - xc;
 426       Float_t yt = y0 - dx*dydx;
 427       Float_t zt = z0 - dx*dzdx;
 428       dy[0] = yc-yt; dz[0]= zc-zt;
 429
 430       // rotate along pad
 431       dy[1] = cost*(dy[0] - dz[0]*tilt);
 432       dz[1] = cost*(dz[0] + dy[0]*tilt);
 433       if(pt>fPtThreshold && c->IsInChamber()) ((TH3S*)arr->At(0))->Fill(dydx, dy[1], sgm[fSegmentLevel]);
 434
 435       // tilt rotation of covariance for clusters
 436       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
 437       cov[0] = (sy2+t2*sz2)*corr;
 438       cov[1] = tilt*(sz2 - sy2)*corr;
 439       cov[2] = (t2*sy2+sz2)*corr;
 440       // sum with track covariance
 441       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
 442       Double_t dyz[2]= {dy[1], dz[1]};
 443       Pulls(dyz, cov, tilt);
 444       ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 445
 446       // Get z-position with respect to anode wire
 447       Int_t istk = geo->GetStack(c->GetDetector());
 448       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
 449       Float_t row0 = pp->GetRow0();
 450       Float_t d  = row0 - zt + pp->GetAnodeWireOffset();
 451       d -= ((Int_t)(2 * d)) / 2.0;
 452       if (d > 0.25) d  = 0.5 - d;
 453
 454       AliTRDclusterInfo *clInfo(NULL);
 455       clInfo = new AliTRDclusterInfo;
 456       clInfo->SetCluster(c);
 457       Float_t covF[] = {cov[0], cov[1], cov[2]};
 458       clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
 459       clInfo->SetResolution(dy[1]);
 460       clInfo->SetAnisochronity(d);
 461       clInfo->SetDriftLength(dx);
 462       clInfo->SetTilt(tilt);
 463       if(fCl) fCl->Add(clInfo);
 464       else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 465
 466       if(DebugLevel()>=1){
 467         if(!clInfoArr){
 468           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 469           clInfoArr->SetOwner(kFALSE);
 470         }
 471         clInfoArr->Add(clInfo);
 472       }
 473     }
 474     if(DebugLevel()>=1 && clInfoArr){
 475       (*DebugStream()) << "cluster"
 476         <<"status="  << status
 477         <<"clInfo.=" << clInfoArr
 478         << "\n";
 479       clInfoArr->Clear();
 480     }
 481   }
 482   if(clInfoArr) delete clInfoArr;
 483
 484   return (TH3S*)arr->At(0);
 485 }
 486
 487
 488 //________________________________________________________
 489 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 490 {
 491 // Plot normalized residuals for tracklets to track.
 492 //
 493 // We start from the result that if X=N(|m|, |Cov|)
 494 // BEGIN_LATEX
 495 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
 496 // END_LATEX
 497 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
 498 // reference position.
 499   if(track) fkTrack = track;
 500   if(!fkTrack){
 501     AliDebug(4, "No Track defined.");
 502     return NULL;
 503   }
 504   TObjArray *arr = NULL;
 505   if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrack ))){
 506     AliWarning("No output container defined.");
 507     return NULL;
 508   }
 509
 510   Int_t sgm[3];
 511   Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
 512   Double_t pt, phi, tht, x, dx, dy[2], dz[2];
 513   AliTRDseedV1 *fTracklet(NULL);
 514   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
 515     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
 516     if(!fTracklet->IsOK()) continue;
 517     sgm[2] = fTracklet->GetDetector();
 518     sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 519     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 520     x   = fTracklet->GetX();
 521     dx  = fTracklet->GetX0() - x;
 522     pt  = fTracklet->GetPt();
 523     phi = fTracklet->GetYref(1);
 524     tht = fTracklet->GetZref(1);
 525     // compute dy and dz
 526     dy[0]= fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
 527     dz[0]= fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
 528     Double_t tilt(fTracklet->GetTilt())
 529             ,t2(tilt*tilt)
 530             ,corr(1./(1. + t2))
 531             ,cost(TMath::Sqrt(corr));
 532     Bool_t rc(fTracklet->IsRowCross());
 533
 534     // calculate residuals using tilt rotation
 535     dy[1]= cost*(dy[0] - dz[0]*tilt);
 536     dz[1]= cost*(dz[0] + dy[0]*tilt);
 537     ((TH3S*)arr->At(0))->Fill(phi, dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
 538     ((TH3S*)arr->At(2))->Fill(tht, dz[1], rc);
 539
 540     // compute covariance matrix
 541     fTracklet->GetCovAt(x, cov);
 542     fTracklet->GetCovRef(covR);
 543     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 544     Double_t dyz[2]= {dy[1], dz[1]};
 545     Pulls(dyz, cov, tilt);
 546     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 547     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
 548
 549     Double_t dphi((phi-fTracklet->GetYfit(1))/(1-phi*fTracklet->GetYfit(1)));
 550     Double_t dtht((tht-fTracklet->GetZfit(1))/(1-tht*fTracklet->GetZfit(1)));
 551     ((TH2I*)arr->At(4))->Fill(phi, TMath::ATan(dphi));
 552
 553     if(DebugLevel()>=1){
 554       UChar_t err(fTracklet->GetErrorMsg());
 555       (*DebugStream()) << "tracklet"
 556         <<"pt="  << pt
 557         <<"phi=" << phi
 558         <<"tht=" << tht
 559         <<"det=" << sgm[2]
 560         <<"dy0="  << dy[0]
 561         <<"dz0="  << dz[0]
 562         <<"dy="  << dy[1]
 563         <<"dz="  << dz[1]
 564         <<"dphi="<< dphi
 565         <<"dtht="<< dtht
 566         <<"dyp=" << dyz[0]
 567         <<"dzp=" << dyz[1]
 568         <<"rc="  << rc
 569         <<"err=" << err
 570         << "\n";
 571     }
 572   }
 573   return (TH2I*)arr->At(0);
 574 }
 575
 576
 577 //________________________________________________________
 578 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
 579 {
 580 // Store resolution/pulls of Kalman before updating with the TRD information
 581 // at the radial position of the first tracklet. The following points are used
 582 // for comparison
 583 //  - the (y,z,snp) of the first TRD tracklet
 584 //  - the (y, z, snp, tgl, pt) of the MC track reference
 585 //
 586 // Additionally the momentum resolution/pulls are calculated for usage in the
 587 // PID calculation.
 588
 589   if(track) fkTrack = track;
 590   if(!fkTrack){
 591     AliDebug(4, "No Track defined.");
 592     return NULL;
 593   }
 594   TObjArray *arr = NULL;
 595   if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackIn))){
 596     AliWarning("No output container defined.");
 597     return NULL;
 598   }
 599   AliExternalTrackParam *tin = NULL;
 600   if(!(tin = fkTrack->GetTrackIn())){
 601     AliWarning("Track did not entered TRD fiducial volume.");
 602     return NULL;
 603   }
 604   TH1 *h = NULL;
 605
 606   Double_t x = tin->GetX();
 607   AliTRDseedV1 *fTracklet = NULL;
 608   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 609     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 610     break;
 611   }
 612   if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
 613     AliWarning("Tracklet did not match Track.");
 614     return NULL;
 615   }
 616   Int_t sgm[3];
 617   sgm[2] = fTracklet->GetDetector();
 618   sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 619   sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 620   Double_t tilt(fTracklet->GetTilt())
 621           ,t2(tilt*tilt)
 622           ,corr(1./(1. + t2))
 623           ,cost(TMath::Sqrt(corr));
 624   Bool_t rc(fTracklet->IsRowCross());
 625
 626   const Int_t kNPAR(5);
 627   Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
 628   Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
 629   Double_t cov[3]; fTracklet->GetCovAt(x, cov);
 630
 631   // define sum covariances
 632   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
 633   Double_t *pc = &covR[0], *pp = &parR[0];
 634   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
 635     PAR(ir) = (*pp);
 636     for(Int_t ic = 0; ic<=ir; ic++,pc++){
 637       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
 638     }
 639   }
 640   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
 641   //COV.Print(); PAR.Print();
 642
 643   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
 644   Double_t dy[2]={parR[0] - fTracklet->GetY(), 0.}
 645           ,dz[2]={parR[1] - fTracklet->GetZ(), 0.}
 646           ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
 647   // calculate residuals using tilt rotation
 648   dy[1] = cost*(dy[0] - dz[0]*tilt);
 649   dz[1] = cost*(dz[0] + dy[0]*tilt);
 650
 651   if(1./PAR[4]>fPtThreshold) ((TH3S*)arr->At(0))->Fill(fTracklet->GetYref(1), dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
 652   ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
 653   ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
 654
 655   Double_t dyz[2] = {dy[1], dz[1]};
 656   Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
 657   Pulls(dyz, cc, tilt);
 658   ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 659   ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
 660
 661
 662
 663   // register reference histo for mini-task
 664   h = (TH2I*)arr->At(0);
 665
 666   if(DebugLevel()>=2){
 667     (*DebugStream()) << "trackIn"
 668       << "x="       << x
 669       << "P="       << &PAR
 670       << "C="       << &COV
 671       << "\n";
 672
 673     Double_t y = fTracklet->GetY();
 674     Double_t z = fTracklet->GetZ();
 675     (*DebugStream()) << "trackletIn"
 676       << "y="       << y
 677       << "z="       << z
 678       << "Vy="      << cov[0]
 679       << "Cyz="     << cov[1]
 680       << "Vz="      << cov[2]
 681       << "\n";
 682   }
 683
 684
 685   if(!HasMCdata()) return h;
 686   UChar_t s;
 687   Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
 688   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
 689   Int_t pdg = fkMC->GetPDG(),
 690         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 691         sign(0);
 692   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 693   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 694   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 695
 696   // translate to reference radial position
 697   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
 698   Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
 699   //Fill MC info
 700   TVectorD PARMC(kNPAR);
 701   PARMC[0]=y0; PARMC[1]=z0;
 702   PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
 703   PARMC[4]=1./pt0;
 704
 705 //   TMatrixDSymEigen eigen(COV);
 706 //   TVectorD evals = eigen.GetEigenValues();
 707 //   TMatrixDSym evalsm(kNPAR);
 708 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
 709 //   TMatrixD evecs = eigen.GetEigenVectors();
 710 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
 711
 712   // fill histos
 713   if(!(arr = (TObjArray*)fContainer->At(kMCtrackIn))) {
 714     AliWarning("No MC container defined.");
 715     return h;
 716   }
 717
 718   // y resolution/pulls
 719   if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
 720   ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)), (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
 721   // z resolution/pulls
 722   ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
 723   ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
 724   // phi resolution/snp pulls
 725   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
 726   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
 727   // theta resolution/tgl pulls
 728   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
 729   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
 730   // pt resolution\\1/pt pulls\\p resolution/pull
 731   ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
 732   ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
 733
 734   Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
 735   p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
 736   ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
 737 //   Float_t sp =
 738 //     p*p*PAR[4]*PAR[4]*COV(4,4)
 739 //    +2.*PAR[3]*COV(3,4)/PAR[4]
 740 //    +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
 741 //   if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
 742
 743   // fill debug for MC
 744   if(DebugLevel()>=3){
 745     (*DebugStream()) << "trackInMC"
 746       << "P="   << &PARMC
 747       << "\n";
 748   }
 749   return h;
 750 }
 751
 752 //________________________________________________________
 753 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
 754 {
 755 // Store resolution/pulls of Kalman after last update with the TRD information
 756 // at the radial position of the first tracklet. The following points are used
 757 // for comparison
 758 //  - the (y,z,snp) of the first TRD tracklet
 759 //  - the (y, z, snp, tgl, pt) of the MC track reference
 760 //
 761 // Additionally the momentum resolution/pulls are calculated for usage in the
 762 // PID calculation.
 763
 764   if(track) fkTrack = track;
 765   if(!fkTrack){
 766     AliDebug(4, "No Track defined.");
 767     return NULL;
 768   }
 769   TObjArray *arr = NULL;
 770   if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackOut))){
 771     AliWarning("No output container defined.");
 772     return NULL;
 773   }
 774   AliExternalTrackParam *tout = NULL;
 775   if(!(tout = fkTrack->GetTrackOut())){
 776     AliDebug(2, "Track did not exit TRD.");
 777     return NULL;
 778   }
 779   TH1 *h(NULL);
 780
 781   Double_t x = tout->GetX();
 782   AliTRDseedV1 *fTracklet(NULL);
 783   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 784     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 785     break;
 786   }
 787   if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
 788     AliWarning("Tracklet did not match Track position.");
 789     return NULL;
 790   }
 791   Int_t sgm[3];
 792   sgm[2] = fTracklet->GetDetector();
 793   sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 794   sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 795   Double_t tilt(fTracklet->GetTilt())
 796           ,t2(tilt*tilt)
 797           ,corr(1./(1. + t2))
 798           ,cost(TMath::Sqrt(corr));
 799   Bool_t rc(fTracklet->IsRowCross());
 800
 801   const Int_t kNPAR(5);
 802   Double_t parR[kNPAR]; memcpy(parR, tout->GetParameter(), kNPAR*sizeof(Double_t));
 803   Double_t covR[3*kNPAR]; memcpy(covR, tout->GetCovariance(), 3*kNPAR*sizeof(Double_t));
 804   Double_t cov[3]; fTracklet->GetCovAt(x, cov);
 805
 806   // define sum covariances
 807   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
 808   Double_t *pc = &covR[0], *pp = &parR[0];
 809   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
 810     PAR(ir) = (*pp);
 811     for(Int_t ic = 0; ic<=ir; ic++,pc++){
 812       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
 813     }
 814   }
 815   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
 816   //COV.Print(); PAR.Print();
 817
 818   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
 819   Double_t dy[3]={parR[0] - fTracklet->GetY(), 0., 0.}
 820           ,dz[3]={parR[1] - fTracklet->GetZ(), 0., 0.}
 821           ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
 822   // calculate residuals using tilt rotation
 823   dy[1] = cost*(dy[0] - dz[0]*tilt);
 824   dz[1] = cost*(dz[0] + dy[0]*tilt);
 825
 826   if(1./PAR[4]>fPtThreshold) ((TH3S*)arr->At(0))->Fill(fTracklet->GetYref(1), 1.e2*dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]); // scale to fit general residual range !!!
 827   ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
 828   ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
 829
 830   Double_t dyz[2] = {dy[1], dz[1]};
 831   Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
 832   Pulls(dyz, cc, tilt);
 833   ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 834   ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
 835
 836   // register reference histo for mini-task
 837   h = (TH2I*)arr->At(0);
 838
 839   if(DebugLevel()>=2){
 840     (*DebugStream()) << "trackOut"
 841       << "x="       << x
 842       << "P="       << &PAR
 843       << "C="       << &COV
 844       << "\n";
 845
 846     Double_t y = fTracklet->GetY();
 847     Double_t z = fTracklet->GetZ();
 848     (*DebugStream()) << "trackletOut"
 849       << "y="       << y
 850       << "z="       << z
 851       << "Vy="      << cov[0]
 852       << "Cyz="     << cov[1]
 853       << "Vz="      << cov[2]
 854       << "\n";
 855   }
 856
 857
 858   if(!HasMCdata()) return h;
 859   UChar_t s;
 860   Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
 861   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
 862   Int_t pdg = fkMC->GetPDG(),
 863         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 864         sign(0);
 865   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 866   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 867   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 868
 869   // translate to reference radial position
 870   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
 871   Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
 872   //Fill MC info
 873   TVectorD PARMC(kNPAR);
 874   PARMC[0]=y0; PARMC[1]=z0;
 875   PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
 876   PARMC[4]=1./pt0;
 877
 878 //   TMatrixDSymEigen eigen(COV);
 879 //   TVectorD evals = eigen.GetEigenValues();
 880 //   TMatrixDSym evalsm(kNPAR);
 881 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
 882 //   TMatrixD evecs = eigen.GetEigenVectors();
 883 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
 884
 885   // fill histos
 886   if(!(arr = (TObjArray*)fContainer->At(kMCtrackOut))){
 887     AliWarning("No MC container defined.");
 888     return h;
 889   }
 890   // y resolution/pulls
 891   if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
 892   ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)), (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
 893   // z resolution/pulls
 894   ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
 895   ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
 896   // phi resolution/snp pulls
 897   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
 898   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
 899   // theta resolution/tgl pulls
 900   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
 901   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
 902   // pt resolution\\1/pt pulls\\p resolution/pull
 903   ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
 904   ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
 905
 906   Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
 907   p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
 908   ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
 909 //   Float_t sp =
 910 //     p*p*PAR[4]*PAR[4]*COV(4,4)
 911 //    +2.*PAR[3]*COV(3,4)/PAR[4]
 912 //    +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
 913 //   if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
 914
 915   // fill debug for MC
 916   if(DebugLevel()>=3){
 917     (*DebugStream()) << "trackOutMC"
 918       << "P="   << &PARMC
 919       << "\n";
 920   }
 921   return h;
 922 }
 923
 924 //________________________________________________________
 925 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 926 {
 927   //
 928   // Plot MC distributions
 929   //
 930
 931   if(!HasMCdata()){
 932     AliDebug(2, "No MC defined. Results will not be available.");
 933     return NULL;
 934   }
 935   if(track) fkTrack = track;
 936   if(!fkTrack){
 937     AliDebug(4, "No Track defined.");
 938     return NULL;
 939   }
 940   if(!fContainer){
 941     AliWarning("No output container defined.");
 942     return NULL;
 943   }
 944   // retriev track characteristics
 945   Int_t pdg = fkMC->GetPDG(),
 946         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 947         sign(0),
 948         sgm[3],
 949         label(fkMC->GetLabel());
 950   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 951   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 952   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 953
 954   TObjArray *arr(NULL);TH1 *h(NULL);
 955   UChar_t s;
 956   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
 957   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 958   Double_t covR[7]/*, cov[3]*/;
 959
 960   if(DebugLevel()>=3){
 961     TVectorD dX(12), dY(12), dZ(12), vPt(12), dPt(12), cCOV(12*15);
 962     fkMC->PropagateKalman(&dX, &dY, &dZ, &vPt, &dPt, &cCOV);
 963     (*DebugStream()) << "MCkalman"
 964       << "pdg=" << pdg
 965       << "dx="  << &dX
 966       << "dy="  << &dY
 967       << "dz="  << &dZ
 968       << "pt="  << &vPt
 969       << "dpt=" << &dPt
 970       << "cov=" << &cCOV
 971       << "\n";
 972   }
 973   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 974   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
 975   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 976     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
 977        !fTracklet->IsOK())*/ continue;
 978
 979     sgm[2] = fTracklet->GetDetector();
 980     sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 981     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 982     Double_t tilt(fTracklet->GetTilt())
 983             ,t2(tilt*tilt)
 984             ,corr(1./(1. + t2))
 985             ,cost(TMath::Sqrt(corr));
 986     x0  = fTracklet->GetX0();
 987     //radial shift with respect to the MC reference (radial position of the pad plane)
 988     x= fTracklet->GetX();
 989     Bool_t rc(fTracklet->IsRowCross());
 990     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
 991     xAnode  = fTracklet->GetX0();
 992
 993     // MC track position at reference radial position
 994     dx  = x0 - x;
 995     if(DebugLevel()>=4){
 996       (*DebugStream()) << "MC"
 997         << "det="     << sgm[2]
 998         << "pdg="     << pdg
 999         << "sgn="     << sign
1000         << "pt="      << pt0
1001         << "x="       << x0
1002         << "y="       << y0
1003         << "z="       << z0
1004         << "dydx="    << dydx0
1005         << "dzdx="    << dzdx0
1006         << "\n";
1007     }
1008     Float_t ymc = y0 - dx*dydx0;
1009     Float_t zmc = z0 - dx*dzdx0;
1010     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
1011
1012     // Kalman position at reference radial position
1013     dx = xAnode - x;
1014     dydx = fTracklet->GetYref(1);
1015     dzdx = fTracklet->GetZref(1);
1016     dzdl = fTracklet->GetTgl();
1017     y  = fTracklet->GetYref(0) - dx*dydx;
1018     dy = y - ymc;
1019     z  = fTracklet->GetZref(0) - dx*dzdx;
1020     dz = z - zmc;
1021     pt = TMath::Abs(fTracklet->GetPt());
1022     fTracklet->GetCovRef(covR);
1023
1024     arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
1025     // y resolution/pulls
1026     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1027     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
1028     // z resolution/pulls
1029     ((TH3S*)arr->At(2))->Fill(dzdx0, dz, 0);
1030     ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
1031     // phi resolution/ snp pulls
1032     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
1033     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
1034     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
1035     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
1036     // theta resolution/ tgl pulls
1037     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
1038               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
1039     ((TH2I*)arr->At(6))->Fill(dzdl0,
1040     TMath::ATan(dtgl));
1041     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
1042     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
1043     Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
1044              p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
1045     ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
1046     ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
1047     ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
1048
1049     // Fill Debug stream for Kalman track
1050     if(DebugLevel()>=4){
1051       (*DebugStream()) << "MCtrack"
1052         << "pt="      << pt
1053         << "x="       << x
1054         << "y="       << y
1055         << "z="       << z
1056         << "dydx="    << dydx
1057         << "dzdx="    << dzdx
1058         << "s2y="     << covR[0]
1059         << "s2z="     << covR[2]
1060         << "\n";
1061     }
1062
1063     // recalculate tracklet based on the MC info
1064     AliTRDseedV1 tt(*fTracklet);
1065     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
1066     tt.SetZref(1, dzdx0);
1067     tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
1068     tt.Fit(1);
1069     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
1070     dydx = tt.GetYfit(1);
1071     dx = x0 - x;
1072     ymc = y0 - dx*dydx0;
1073     zmc = z0 - dx*dzdx0;
1074     dy = y-ymc;
1075     dz = z-zmc;
1076     Float_t dphi  = (dydx - dydx0);
1077     dphi /= (1.- dydx*dydx0);
1078
1079     // add tracklet residuals for y and dydx
1080     arr = (TObjArray*)fContainer->At(kMCtracklet);
1081
1082     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1083     if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
1084     ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
1085     if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
1086     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
1087
1088     // Fill Debug stream for tracklet
1089     if(DebugLevel()>=4){
1090       Float_t s2y = tt.GetS2Y();
1091       Float_t s2z = tt.GetS2Z();
1092       (*DebugStream()) << "MCtracklet"
1093         << "rc="    << rc
1094         << "x="     << x
1095         << "y="     << y
1096         << "z="     << z
1097         << "dydx="  << dydx
1098         << "s2y="   << s2y
1099         << "s2z="   << s2z
1100         << "\n";
1101     }
1102
1103     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
1104     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
1105     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
1106
1107     arr = (TObjArray*)fContainer->At(kMCcluster);
1108     AliTRDcluster *c = NULL;
1109     tt.ResetClusterIter(kFALSE);
1110     while((c = tt.PrevCluster())){
1111       Float_t  q = TMath::Abs(c->GetQ());
1112       x = c->GetX()+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; y = c->GetY();z = c->GetZ();
1113       dx = x0 - x;
1114       ymc= y0 - dx*dydx0;
1115       zmc= z0 - dx*dzdx0;
1116       dy = cost*(y - ymc - tilt*(z-zmc));
1117       dz = cost*(z - zmc + tilt*(y-ymc));
1118
1119       // Fill Histograms
1120       if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){
1121         ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1122         ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
1123       }
1124
1125       // Fill calibration container
1126       Float_t d = zr0 - zmc;
1127       d -= ((Int_t)(2 * d)) / 2.0;
1128       if (d > 0.25) d  = 0.5 - d;
1129       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
1130       clInfo->SetCluster(c);
1131       clInfo->SetMC(pdg, label);
1132       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
1133       clInfo->SetResolution(dy);
1134       clInfo->SetAnisochronity(d);
1135       clInfo->SetDriftLength(dx);
1136       clInfo->SetTilt(tilt);
1137       if(fMCcl) fMCcl->Add(clInfo);
1138       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
1139       if(DebugLevel()>=5){
1140         if(!clInfoArr){
1141           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
1142           clInfoArr->SetOwner(kFALSE);
1143         }
1144         clInfoArr->Add(clInfo);
1145       }
1146     }
1147     // Fill Debug Tree
1148     if(DebugLevel()>=5 && clInfoArr){
1149       (*DebugStream()) << "MCcluster"
1150         <<"clInfo.=" << clInfoArr
1151         << "\n";
1152       clInfoArr->Clear();
1153     }
1154   }
1155   if(clInfoArr) delete clInfoArr;
1156   return h;
1157 }
1158
1159
1160 //________________________________________________________
1161 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1162 {
1163   //
1164   // Get the reference figures
1165   //
1166
1167   Float_t xy[4] = {0., 0., 0., 0.};
1168   if(!gPad){
1169     AliWarning("Please provide a canvas to draw results.");
1170     return kFALSE;
1171   }
1172   Int_t selection[100], n(0), selStart(0); //
1173   Int_t ly0(0), dly(5);
1174   //Int_t ly0(1), dly(2); // used for SA
1175   TList *l(NULL); TVirtualPad *pad(NULL);
1176   TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
1177   switch(ifig){
1178   case 0: // charge resolution
1179     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1180     ((TVirtualPad*)l->At(0))->cd();
1181     ga=((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0));
1182     if(ga->GetN()) ga->Draw("apl");
1183     ((TVirtualPad*)l->At(1))->cd();
1184     g = ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0));
1185     if(g->GetN()) g->Draw("apl");
1186     break;
1187   case 1: // cluster2track residuals
1188     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1189     xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1190     pad = (TVirtualPad*)l->At(0); pad->cd();
1191     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1192     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1193     if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1194     pad=(TVirtualPad*)l->At(1); pad->cd();
1195     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1196     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1197     if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1198     return kTRUE;
1199   case 2: // cluster2track residuals
1200     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1201     xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1202     pad = (TVirtualPad*)l->At(0); pad->cd();
1203     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1204     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1205     if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1206     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-0.5; xy[3] = 2.5;
1207     pad=(TVirtualPad*)l->At(1); pad->cd();
1208     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1209     if(!GetGraphArray(xy, kCluster, 1, 1)) break;
1210     return kTRUE;
1211   case 3: // kTrack y
1212     gPad->Divide(3, 2, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1213     xy[0] = -.3; xy[1] = -20.; xy[2] = .3; xy[3] = 100.;
1214     ((TVirtualPad*)l->At(0))->cd();
1215     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1216     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1217
1218     ((TVirtualPad*)l->At(1))->cd();
1219     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1220     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1221
1222     ((TVirtualPad*)l->At(2))->cd();
1223     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1224     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1225
1226     ((TVirtualPad*)l->At(3))->cd();
1227     selStart=fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1228     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1229
1230     ((TVirtualPad*)l->At(4))->cd();
1231     selStart=fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1232     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1233
1234     ((TVirtualPad*)l->At(5))->cd();
1235     selStart=2*fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1236     if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1237     return kTRUE;
1238   case 4: // kTrack z
1239     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1240
1241     xy[0] = -1.; xy[1] = -150.; xy[2] = 1.; xy[3] = 1000.;
1242     ((TVirtualPad*)l->At(0))->cd();
1243     selection[0]=1;
1244     if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1245
1246     xy[0] = -1.; xy[1] = -1500.; xy[2] = 1.; xy[3] = 10000.;
1247     ((TVirtualPad*)l->At(1))->cd();
1248     selection[0]=0;
1249     if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1250
1251     return kTRUE;
1252   case 5: // kTrack  pulls
1253     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1254
1255     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1256     ((TVirtualPad*)l->At(0))->cd();
1257     if(!GetGraphArray(xy, kTrack, 1, 1)) break;
1258
1259     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1260     ((TVirtualPad*)l->At(1))->cd();
1261     if(!GetGraphArray(xy, kTrack, 3, 1)) break;
1262     return kTRUE;
1263   case 6: // kTrack  phi
1264     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1265     if(GetGraph(&xy[0], kTrack , 4)) return kTRUE;
1266     break;
1267   case 7: // kTrackIn y
1268     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1269     xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1270     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1271     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1272     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1273     if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1274     pad=((TVirtualPad*)l->At(1)); pad->cd();
1275     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1276     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1277     if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1278     return kTRUE;
1279   case 8: // kTrackIn y
1280     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1281     xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1282     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1283     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1284     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1285     if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1286     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1287     pad=((TVirtualPad*)l->At(1)); pad->cd();
1288     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1289     if(!GetGraphArray(xy, kTrackIn, 1, 1)) break;
1290     return kTRUE;
1291   case 9: // kTrackIn z
1292     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1293     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1294     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1295     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1296     selection[0]=1;
1297     if(!GetGraphArray(xy, kTrackIn, 2, 1, 1, selection)) break;
1298     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1299     pad = ((TVirtualPad*)l->At(1)); pad->cd();
1300     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1301     if(!GetGraphArray(xy, kTrackIn, 3, 1)) break;
1302     return kTRUE;
1303   case 10: // kTrackIn phi
1304     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1305     if(GetGraph(&xy[0], kTrackIn, 4)) return kTRUE;
1306     break;
1307   case 11: // kTrackOut y
1308     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1309     xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1310     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1311     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1312     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1313     if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1314     pad=((TVirtualPad*)l->At(1)); pad->cd();
1315     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1316     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1317     if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1318     return kTRUE;
1319   case 12: // kTrackOut y
1320     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1321     xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1322     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1323     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1324     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1325     if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1326     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1327     pad=((TVirtualPad*)l->At(1)); pad->cd();
1328     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1329     if(!GetGraphArray(xy, kTrackOut, 1, 1)) break;
1330     return kTRUE;
1331   case 13: // kTrackOut z
1332     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1333     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1334     pad = ((TVirtualPad*)l->At(0)); pad->cd();
1335     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1336     if(!GetGraphArray(xy, kTrackOut, 2, 1)) break;
1337     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1338     pad = ((TVirtualPad*)l->At(1)); pad->cd();
1339     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1340     if(!GetGraphArray(xy, kTrackOut, 3, 1)) break;
1341     return kTRUE;
1342   case 14: // kTrackOut phi
1343     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1344     if(GetGraph(&xy[0], kTrackOut, 4)) return kTRUE;
1345     break;
1346   case 15: // kMCcluster
1347     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1348     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
1349     ((TVirtualPad*)l->At(0))->cd();
1350     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1351     if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1352     ((TVirtualPad*)l->At(1))->cd();
1353     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1354     if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1355     return kTRUE;
1356   case 16: // kMCcluster
1357     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1358     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
1359     ((TVirtualPad*)l->At(0))->cd();
1360     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1361     if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1362     ((TVirtualPad*)l->At(1))->cd();
1363     xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1364     if(!GetGraphArray(xy, kMCcluster, 1, 1)) break;
1365     return kTRUE;
1366   case 17: //kMCtracklet [y]
1367     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1368     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1369     ((TVirtualPad*)l->At(0))->cd();
1370     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1371     if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1372     ((TVirtualPad*)l->At(1))->cd();
1373     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1374     if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1375     return kTRUE;
1376   case 18: //kMCtracklet [y]
1377     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1378     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1379     ((TVirtualPad*)l->At(0))->cd();
1380     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1381     if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1382     ((TVirtualPad*)l->At(1))->cd();
1383     xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1384     if(!GetGraphArray(xy, kMCtracklet, 1, 1)) break;
1385     return kTRUE;
1386   case 19: //kMCtracklet [z]
1387     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1388     xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
1389     ((TVirtualPad*)l->At(0))->cd();
1390     if(!GetGraphArray(xy, kMCtracklet, 2)) break;
1391     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1392     ((TVirtualPad*)l->At(1))->cd();
1393     if(!GetGraphArray(xy, kMCtracklet, 3)) break;
1394     return kTRUE;
1395   case 20: //kMCtracklet [phi]
1396     xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
1397     if(!GetGraph(&xy[0], kMCtracklet, 4)) break;
1398     return kTRUE;
1399   case 21: //kMCtrack [y] ly [0]
1400     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1401     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1402     ((TVirtualPad*)l->At(0))->cd();
1403     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*0.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1404     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1405     ((TVirtualPad*)l->At(1))->cd();
1406     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*0.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1407     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1408     return kTRUE;
1409   case 22: //kMCtrack [y] ly [1]
1410     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1411     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1412     ((TVirtualPad*)l->At(0))->cd();
1413     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*1.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1414     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1415     ((TVirtualPad*)l->At(1))->cd();
1416     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*1.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1417     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1418     return kTRUE;
1419   case 23: //kMCtrack [y] ly [2]
1420     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1421     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1422     ((TVirtualPad*)l->At(0))->cd();
1423     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*2.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1424     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1425     ((TVirtualPad*)l->At(1))->cd();
1426     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*2.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1427     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1428     return kTRUE;
1429   case 24: //kMCtrack [y] ly [3]
1430     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1431     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1432     ((TVirtualPad*)l->At(0))->cd();
1433     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*3.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1434     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1435     ((TVirtualPad*)l->At(1))->cd();
1436     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*3.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1437     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1438     return kTRUE;
1439   case 25: //kMCtrack [y] ly [4]
1440     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1441     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1442     ((TVirtualPad*)l->At(0))->cd();
1443     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1444     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1445     ((TVirtualPad*)l->At(1))->cd();
1446     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1447     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1448     return kTRUE;
1449   case 26: //kMCtrack [y] ly [5]
1450     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1451     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1452     ((TVirtualPad*)l->At(0))->cd();
1453     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1454     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1455     ((TVirtualPad*)l->At(1))->cd();
1456     selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1457     if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1458     return kTRUE;
1459   case 27: //kMCtrack [y pulls]
1460     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1461     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 5.5;
1462     ((TVirtualPad*)l->At(0))->cd();
1463     selStart=0; for(n=0; n<6; n++) selection[n]=selStart+n;
1464     if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1465     ((TVirtualPad*)l->At(1))->cd();
1466     selStart=6; for(n=0; n<6; n++) selection[n]=selStart+n;
1467     if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1468     return kTRUE;
1469   case 28: //kMCtrack [z]
1470     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1471     xy[0]=-1.; xy[1]=-1500.; xy[2]=1.; xy[3] =6000.;
1472     ((TVirtualPad*)l->At(0))->cd();
1473     if(!GetGraphArray(xy, kMCtrack, 2)) break;
1474     xy[0] = -1.; xy[1] = -1.5; xy[2] = 1.; xy[3] = 5.;
1475     ((TVirtualPad*)l->At(1))->cd();
1476     if(!GetGraphArray(xy, kMCtrack, 3)) break;
1477     return kTRUE;
1478   case 29: //kMCtrack [phi/snp]
1479     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1480     xy[0]=-.2; xy[1]=-0.5; xy[2]=.2; xy[3] =10.;
1481     ((TVirtualPad*)l->At(0))->cd();
1482     if(!GetGraphArray(xy, kMCtrack, 4)) break;
1483     xy[0] = -.2; xy[1] = -1.5; xy[2] = .2; xy[3] = 5.;
1484     ((TVirtualPad*)l->At(1))->cd();
1485     if(!GetGraphArray(xy, kMCtrack, 5)) break;
1486     return kTRUE;
1487   case 30: //kMCtrack [theta/tgl]
1488     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1489     xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
1490     ((TVirtualPad*)l->At(0))->cd();
1491     if(!GetGraphArray(xy, kMCtrack, 6)) break;
1492     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1493     ((TVirtualPad*)l->At(1))->cd();
1494     if(!GetGraphArray(xy, kMCtrack, 7)) break;
1495     return kTRUE;
1496   case 31: //kMCtrack [pt]
1497     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1498     pad = (TVirtualPad*)l->At(0); pad->cd();
1499     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1500     // pi selection
1501     n=0;
1502     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1503       selection[n++] = il*11 + 2; // pi-
1504       selection[n++] = il*11 + 8; // pi+
1505     }
1506     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1507     //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1508     if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#pi#pm")) break;
1509     pad->Modified(); pad->Update(); pad->SetLogx();
1510     pad = (TVirtualPad*)l->At(1); pad->cd();
1511     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1512     // mu selection
1513     n=0;
1514     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1515       selection[n++] = il*11 + 3; // mu-
1516       selection[n++] = il*11 + 7; // mu+
1517     }
1518     if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#mu#pm")) break;
1519     pad->Modified(); pad->Update(); pad->SetLogx();
1520     return kTRUE;
1521   case 32: //kMCtrack [pt]
1522     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1523     pad = (TVirtualPad*)l->At(0); pad->cd();
1524     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1525     // p selection
1526     n=0;
1527     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1528       selection[n++] = il*11 + 0; // p bar
1529       selection[n++] = il*11 + 10; // p
1530     }
1531     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1532     //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1533     if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "p&p bar")) break;
1534     pad->Modified(); pad->Update(); pad->SetLogx();
1535     pad = (TVirtualPad*)l->At(1); pad->cd();
1536     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1537     // e selection
1538     n=0;
1539     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1540       selection[n++] = il*11 + 4; // e-
1541       selection[n++] = il*11 + 6; // e+
1542     }
1543     xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1544     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1545     if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "e#pm")) break;
1546     pad->Modified(); pad->Update(); pad->SetLogx();
1547     return kTRUE;
1548   case 33: //kMCtrack [1/pt] pulls
1549     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1550     //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 4.5; // SA
1551     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1552     pad = (TVirtualPad*)l->At(0); pad->cd();
1553     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1554     // pi selection
1555     n=0;
1556     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1557       selection[n++] = il*11 + 2; // pi-
1558       selection[n++] = il*11 + 8; // pi+
1559     }
1560     if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#pi#pm")) break;
1561     pad = (TVirtualPad*)l->At(1); pad->cd();
1562     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1563     // mu selection
1564     n=0;
1565     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1566       selection[n++] = il*11 + 3; // mu-
1567       selection[n++] = il*11 + 7; // mu+
1568     }
1569     if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#mu#pm")) break;
1570     return kTRUE;
1571   case 34: //kMCtrack [1/pt] pulls
1572     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1573     pad = (TVirtualPad*)l->At(0); pad->cd();
1574     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1575     // p selection
1576     n=0;
1577     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1578       selection[n++] = il*11 + 0; // p bar
1579       selection[n++] = il*11 + 10; // p
1580     }
1581     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1582     //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 6.; // SA
1583     if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "p & p bar")) break;
1584     pad = (TVirtualPad*)l->At(1); pad->cd();
1585     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1586     // e selection
1587     n=0;
1588     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1589       selection[n++] = il*11 + 4; // e-
1590       selection[n++] = il*11 + 6; // e+
1591     }
1592     xy[0] = 0.; xy[1] = -2.; xy[2] = 2.; xy[3] = 4.5;
1593     if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "e#pm")) break;
1594     return kTRUE;
1595   case 35: //kMCtrack [p]
1596     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1597     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1598     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1599     pad = (TVirtualPad*)l->At(0); pad->cd();
1600     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1601     // pi selection
1602     n=0;
1603     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1604       selection[n++] = il*11 + 2; // pi-
1605       selection[n++] = il*11 + 8; // pi+
1606     }
1607     if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#pi#pm")) break;
1608     pad->Modified(); pad->Update(); pad->SetLogx();
1609     pad = (TVirtualPad*)l->At(1); pad->cd();
1610     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1611     // mu selection
1612     n=0;
1613     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1614       selection[n++] = il*11 + 3; // mu-
1615       selection[n++] = il*11 + 7; // mu+
1616     }
1617     if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#mu#pm")) break;
1618     pad->Modified(); pad->Update(); pad->SetLogx();
1619     return kTRUE;
1620   case 36: //kMCtrack [p]
1621     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1622     pad = (TVirtualPad*)l->At(0); pad->cd();
1623     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1624     // p selection
1625     n=0;
1626     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1627       selection[n++] = il*11 + 0; // p bar
1628       selection[n++] = il*11 + 10; // p
1629     }
1630     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1631     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.; // SA
1632     if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "p & p bar")) break;
1633     pad->Modified(); pad->Update(); pad->SetLogx();
1634     pad = (TVirtualPad*)l->At(1); pad->cd();
1635     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1636     // e selection
1637     n=0;
1638     for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1639       selection[n++] = il*11 + 4; // e-
1640       selection[n++] = il*11 + 6; // e+
1641     }
1642     xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1643     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1644     if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "e#pm")) break;
1645     pad->Modified(); pad->Update(); pad->SetLogx();
1646     return kTRUE;
1647   case 37: // kMCtrackIn [y]
1648     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1649     xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1650     ((TVirtualPad*)l->At(0))->cd();
1651     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1652     if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1653     ((TVirtualPad*)l->At(1))->cd();
1654     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1655     if(!GetGraphArray(&xy[0], kMCtrackIn, 0, 1, n, selection)) break;
1656     return kTRUE;
1657   case 38: // kMCtrackIn [y]
1658     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1659     xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1660     ((TVirtualPad*)l->At(0))->cd();
1661     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1662     if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1663     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1664     ((TVirtualPad*)l->At(1))->cd();
1665     if(!GetGraphArray(xy, kMCtrackIn, 1, 1)) break;
1666     return kTRUE;
1667   case 39: // kMCtrackIn [z]
1668     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1669     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
1670     ((TVirtualPad*)l->At(0))->cd();
1671     if(!GetGraphArray(xy, kMCtrackIn, 2, 1)) break;
1672     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1673     ((TVirtualPad*)l->At(1))->cd();
1674     if(!GetGraphArray(xy, kMCtrackIn, 3, 1)) break;
1675     return kTRUE;
1676   case 40: // kMCtrackIn [phi|snp]
1677     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1678     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1679     ((TVirtualPad*)l->At(0))->cd();
1680     if(!GetGraph(&xy[0], kMCtrackIn, 4)) break;
1681     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1682     ((TVirtualPad*)l->At(1))->cd();
1683     if(!GetGraph(&xy[0], kMCtrackIn, 5)) break;
1684     return kTRUE;
1685   case 41: // kMCtrackIn [theta|tgl]
1686     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1687     xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1688     ((TVirtualPad*)l->At(0))->cd();
1689     if(!GetGraph(&xy[0], kMCtrackIn, 6)) break;
1690     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1691     ((TVirtualPad*)l->At(1))->cd();
1692     if(!GetGraph(&xy[0], kMCtrackIn, 7)) break;
1693     return kTRUE;
1694   case 42: // kMCtrackIn [pt]
1695     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1696     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1697     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.; // SA
1698     pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
1699     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1700     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1701     if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1702     pad = (TVirtualPad*)l->At(1); pad->cd(); pad->SetLogx();
1703     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1704     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1705     if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1706     return kTRUE;
1707   case 43: //kMCtrackIn [1/pt] pulls
1708     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1709     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1710     pad = (TVirtualPad*)l->At(0); pad->cd();
1711     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1712     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1713     if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1714     pad = (TVirtualPad*)l->At(1); pad->cd();
1715     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1716     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1717     if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1718     return kTRUE;
1719   case 44: // kMCtrackIn [p]
1720     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1721     //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1722     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1723     pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1724     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1725     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1726     if(!GetGraphArray(xy, kMCtrackIn, 10, 1, n, selection)) break;
1727     pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1728     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1729     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1730     if(!GetGraphArray(xy, kMCtrackIn, 10, 1,  n, selection)) break;
1731     return kTRUE;
1732   case 45: // kMCtrackOut [y]
1733     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1734     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1735     ((TVirtualPad*)l->At(0))->cd();
1736     selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1737     if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1738     ((TVirtualPad*)l->At(1))->cd();
1739     selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1740     if(!GetGraphArray(&xy[0], kMCtrackOut, 0, 1, n, selection)) break;
1741     return kTRUE;
1742   case 46: // kMCtrackOut [y]
1743     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1744     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1745     ((TVirtualPad*)l->At(0))->cd();
1746     selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1747     if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1748     xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1749     ((TVirtualPad*)l->At(1))->cd();
1750     if(!GetGraphArray(xy, kMCtrackOut, 1, 1)) break;
1751     return kTRUE;
1752   case 47: // kMCtrackOut [z]
1753     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1754     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =1500.;
1755     ((TVirtualPad*)l->At(0))->cd();
1756     if(!GetGraphArray(xy, kMCtrackOut, 2, 1)) break;
1757     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1758     ((TVirtualPad*)l->At(1))->cd();
1759     if(!GetGraphArray(xy, kMCtrackOut, 3, 1)) break;
1760     return kTRUE;
1761   case 48: // kMCtrackOut [phi|snp]
1762     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1763     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1764     ((TVirtualPad*)l->At(0))->cd();
1765     if(!GetGraph(&xy[0], kMCtrackOut, 4)) break;
1766     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1767     ((TVirtualPad*)l->At(1))->cd();
1768     if(!GetGraph(&xy[0], kMCtrackOut, 5)) break;
1769     return kTRUE;
1770   case 49: // kMCtrackOut [theta|tgl]
1771     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1772     xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1773     ((TVirtualPad*)l->At(0))->cd();
1774     if(!GetGraph(&xy[0], kMCtrackOut, 6)) break;
1775     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 15.;
1776     ((TVirtualPad*)l->At(1))->cd();
1777     if(!GetGraph(&xy[0], kMCtrackOut, 7)) break;
1778     return kTRUE;
1779   case 50: // kMCtrackOut [pt]
1780     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1781     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1782     pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
1783     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1784     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1785     if(!GetGraphArray(xy, kMCtrackOut, 8, 1, n, selection)) break;
1786     pad = (TVirtualPad*)l->At(1); pad->cd();pad->SetLogx();
1787     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1788     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1789     if(!GetGraphArray(xy, kMCtrackOut, 8, 1, n, selection)) break;
1790     return kTRUE;
1791   case 51: //kMCtrackOut [1/pt] pulls
1792     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1793     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1794     pad = (TVirtualPad*)l->At(0); pad->cd();
1795     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1796     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1797     if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1798     pad = (TVirtualPad*)l->At(1); pad->cd();
1799     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1800     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1801     if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1802     return kTRUE;
1803   case 52: // kMCtrackOut [p]
1804     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1805     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1806     pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1807     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1808     n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1809     if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1810     pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1811     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1812     n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1813     if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1814     return kTRUE;
1815   }
1816   AliWarning(Form("Reference plot [%d] missing result", ifig));
1817   return kFALSE;
1818 }
1819
1820 //________________________________________________________
1821 void AliTRDresolution::MakeSummary()
1822 {
1823 // Build summary plots
1824
1825   if(!fGraphS || !fGraphM){
1826     AliError("Missing results");
1827     return;
1828   }
1829   Float_t xy[4] = {0., 0., 0., 0.};
1830   Float_t range[2];
1831   TH2 *h2 = new TH2I("h2SF", "", 20, -.2, .2, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5);
1832   h2->GetXaxis()->CenterTitle();
1833   h2->GetYaxis()->CenterTitle();
1834   h2->GetZaxis()->CenterTitle();h2->GetZaxis()->SetTitleOffset(1.4);
1835
1836   Int_t ih2(0), iSumPlot(0);
1837   TCanvas *cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Cluster & Tracklet Resolution", 1024, 768);
1838   cOut->Divide(3,2, 2.e-3, 2.e-3, kYellow-7);
1839   TVirtualPad *p(NULL);
1840
1841   p=cOut->cd(1);
1842   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1843   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1844   h2->SetTitle(Form("Cluster-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1845   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphS->At(kCluster))->At(0), h2);
1846   GetRange(h2, 1, range);
1847   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1848   h2->Draw("colz");
1849   h2->SetContour(7);
1850
1851   p=cOut->cd(2);
1852   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1853   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1854   h2->SetTitle(Form("Cluster-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1855   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphM->At(kCluster))->At(0), h2);
1856   GetRange(h2, 0, range);
1857   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1858   h2->Draw("colz");
1859   h2->SetContour(7);
1860
1861   p=cOut->cd(3);
1862   p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1863   xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1864   if(!GetGraphArray(xy, kCluster, 1, 1)) return;
1865
1866   p=cOut->cd(4);
1867   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1868   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1869   h2->SetTitle(Form("Tracklet-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1870   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphS->At(kTrack))->At(0), h2);
1871   GetRange(h2, 1, range);
1872   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1873   h2->Draw("colz");
1874   h2->SetContour(7);
1875
1876   p=cOut->cd(5);
1877   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1878   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1879   h2->SetTitle(Form("Tracklet-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1880   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphM->At(kTrack))->At(0), h2);
1881   GetRange(h2, 0, range);
1882   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1883   h2->Draw("colz");
1884   h2->SetContour(7);
1885
1886   p=cOut->cd(6);
1887   p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1888   xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1889   if(!GetGraphArray(xy, kTrack, 1, 1)) return;
1890
1891   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1892
1893   if(!HasMCdata() ||
1894     (!fGraphS->At(kMCcluster) || !fGraphM->At(kMCcluster) ||
1895      !fGraphS->At(kMCtracklet) || !fGraphM->At(kMCtracklet))){
1896     delete cOut;
1897     return;
1898   }
1899   cOut->Clear(); cOut->SetName(Form("TRDsummary%s_%d", GetName(), iSumPlot++));
1900   cOut->Divide(3, 2, 2.e-3, 2.e-3, kBlue-10);
1901
1902   p=cOut->cd(1);
1903   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1904   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1905   h2->SetTitle(Form("Cluster-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1906   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphS->At(kMCcluster))->At(0), h2);
1907   GetRange(h2, 1, range);
1908   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1909   h2->Draw("colz");
1910   h2->SetContour(7);
1911
1912   p=cOut->cd(2);
1913   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1914   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1915   h2->SetContour(7);
1916   h2->SetTitle(Form("Cluster-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1917   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphM->At(kMCcluster))->At(0), h2);
1918   GetRange(h2, 0, range);
1919   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1920   h2->Draw("colz");
1921   h2->SetContour(7);
1922
1923   p=cOut->cd(3);
1924   p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1925   xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1926   if(!GetGraphArray(xy, kMCcluster, 1, 1)) AliWarning("Missing MC cluster plot.");
1927
1928   p=cOut->cd(4);
1929   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1930   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1931   h2->SetContour(7);
1932   h2->SetTitle(Form("Tracklet-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1933   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphS->At(kMCtracklet))->At(0), h2);
1934   GetRange(h2, 1, range);
1935   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1936   h2->Draw("colz");
1937   h2->SetContour(7);
1938
1939   p=cOut->cd(5);
1940   p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1941   h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1942   h2->SetContour(7);
1943   h2->SetTitle(Form("Tracklet-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1944   MakeSummaryPlot((TObjArray*)  ((TObjArray*)fGraphM->At(kMCtracklet))->At(0), h2);
1945   GetRange(h2, 0, range);
1946   h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1947   h2->Draw("colz");
1948   h2->SetContour(7);
1949
1950   p=cOut->cd(6);
1951   p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1952   xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1953   GetGraphArray(xy, kMCtracklet, 1, 1);
1954
1955   cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1956   delete cOut;
1957 }
1958
1959 //________________________________________________________
1960 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1961 {
1962 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1963 // The "range" vector should be initialized with 2 elements
1964 // Option "mod" can be any of
1965 //   - 0 : gaussian like distribution
1966 //   - 1 : tailed distribution
1967
1968   Int_t nx(h2->GetNbinsX())
1969       , ny(h2->GetNbinsY())
1970       , n(nx*ny);
1971   Double_t *data=new Double_t[n];
1972   for(Int_t ix(1), in(0); ix<=nx; ix++){
1973     for(Int_t iy(1); iy<=ny; iy++)
1974       data[in++] = h2->GetBinContent(ix, iy);
1975   }
1976   Double_t mean, sigm;
1977   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1978
1979   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1980   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1981   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1982   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1983   h1.FillN(n,data,0);
1984   delete [] data;
1985
1986   switch(mod){
1987   case 0:// gaussian distribution
1988   {
1989     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1990     h1.Fit(&fg, "QN");
1991     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1992     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1993     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1994     break;
1995   }
1996   case 1:// tailed distribution
1997   {
1998     Int_t bmax(h1.GetMaximumBin());
1999     Int_t jBinMin(1), jBinMax(100);
2000     for(Int_t ibin(bmax); ibin--;){
2001       if(h1.GetBinContent(ibin)<1.){
2002         jBinMin=ibin; break;
2003       }
2004     }
2005     for(Int_t ibin(bmax); ibin++;){
2006       if(h1.GetBinContent(ibin)<1.){
2007         jBinMax=ibin; break;
2008       }
2009     }
2010     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
2011     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
2012     break;
2013   }
2014   }
2015
2016   return;
2017 }
2018
2019 //________________________________________________________
2020 void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
2021 {
2022 // Core functionality for MakeSummary function.
2023
2024   h2->Reset();
2025   Double_t x,y;
2026   TGraphErrors *g(NULL); TAxis *ax(h2->GetXaxis());
2027   for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
2028     g=(TGraphErrors*)a->At(iseg);
2029     for(Int_t in(0); in<g->GetN(); in++){
2030       g->GetPoint(in, x, y);
2031       h2->SetBinContent(ax->FindBin(x), iseg+1, y);
2032     }
2033   }
2034 }
2035
2036
2037 //________________________________________________________
2038 Bool_t AliTRDresolution::PostProcess()
2039 {
2040 // Fit, Project, Combine, Extract values from the containers filled during execution
2041
2042   /*fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));*/
2043   if (!fContainer) {
2044     AliError("ERROR: list not available");
2045     return kFALSE;
2046   }
2047
2048   // define general behavior parameters
2049   const Color_t fgColorS[11]={
2050   kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
2051   kGray,
2052   kRed, kViolet, kMagenta+1, kOrange-3, kOrange
2053   };
2054   const Color_t fgColorM[11]={
2055   kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
2056   kBlack,
2057   kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
2058   };
2059   const Marker_t fgMarker[11]={
2060   30, 30, 26, 25, 24,
2061   28,
2062   20, 21, 22, 29, 29
2063   };
2064
2065   TGraph *gm= NULL, *gs= NULL;
2066   if(!fGraphS && !fGraphM){
2067     TObjArray *aM(NULL), *aS(NULL);
2068     Int_t n = fContainer->GetEntriesFast();
2069     fGraphS = new TObjArray(n); fGraphS->SetOwner();
2070     fGraphM = new TObjArray(n); fGraphM->SetOwner();
2071     for(Int_t ig(0), nc(0); ig<n; ig++){
2072       fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
2073       fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
2074
2075       for(Int_t ic=0; ic<fgNproj[ig]; ic++, nc++){
2076         AliDebug(2, Form("building G[%d] P[%d] N[%d]", ig, ic, fNcomp[nc]));
2077         if(fNcomp[nc]>1){
2078           TObjArray *agS(NULL), *agM(NULL);
2079           aS->AddAt(agS = new TObjArray(fNcomp[nc]), ic);
2080           aM->AddAt(agM = new TObjArray(fNcomp[nc]), ic);
2081           for(Int_t is=fNcomp[nc]; is--;){
2082             agS->AddAt(gs = new TGraphErrors(), is);
2083             Int_t is0(is%11), il0(is/11);
2084             gs->SetMarkerStyle(fgMarker[is0]);
2085             gs->SetMarkerColor(fgColorS[is0]);
2086             gs->SetLineColor(fgColorS[is0]);
2087             gs->SetLineStyle(il0);gs->SetLineWidth(2);
2088             gs->SetName(Form("s_%d_%02d_%02d", ig, ic, is));
2089
2090             agM->AddAt(gm = new TGraphErrors(), is);
2091             gm->SetMarkerStyle(fgMarker[is0]);
2092             gm->SetMarkerColor(fgColorM[is0]);
2093             gm->SetLineColor(fgColorM[is0]);
2094             gm->SetLineStyle(il0);gm->SetLineWidth(2);
2095             gm->SetName(Form("m_%d_%02d_%02d", ig, ic, is));
2096             // this is important for labels in the legend
2097             if(ic==0) {
2098               gs->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2099               gm->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2100             } else if(ic==1) {
2101               gs->SetTitle(Form("%s Ly[%d]", is%2 ?"z":"y", is/2));
2102               gm->SetTitle(Form("%s Ly[%d]", is%2?"z":"y", is/2));
2103             } else if(ic==2||ic==3) {
2104               gs->SetTitle(Form("%s Ly[%d]", is%2 ?"RC":"no RC", is/2));
2105               gm->SetTitle(Form("%s Ly[%d]", is%2?"RC":"no RC", is/2));
2106             } else if(ic<=7) {
2107               gs->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2108               gm->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2109             } else {
2110               gs->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2111               gm->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2112             }
2113           }
2114         } else {
2115           aS->AddAt(gs = new TGraphErrors(), ic);
2116           gs->SetMarkerStyle(23);
2117           gs->SetMarkerColor(kRed);
2118           gs->SetLineColor(kRed);
2119           gs->SetNameTitle(Form("s_%d_%02d", ig, ic), "sigma");
2120
2121           aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
2122           gm->SetLineColor(kBlack);
2123           gm->SetMarkerStyle(7);
2124           gm->SetMarkerColor(kBlack);
2125           gm->SetNameTitle(Form("m_%d_%02d", ig, ic), "mean");
2126         }
2127       }
2128     }
2129   }
2130
2131 /*  printf("\n\n\n"); fGraphS->ls();
2132   printf("\n\n\n"); fGraphM->ls();*/
2133
2134
2135   // DEFINE MODELS
2136   // simple gauss
2137   TF1 fg("fGauss", "gaus", -.5, .5);
2138   // Landau for charge resolution
2139   TF1 fch("fClCh", "landau", 0., 1000.);
2140   // Landau for e+- pt resolution
2141   TF1 fpt("fPt", "landau", -0.1, 0.2);
2142
2143   //PROCESS EXPERIMENTAL DISTRIBUTIONS
2144   // Charge resolution
2145   //Process3DL(kCharge, 0, &fl);
2146   // Clusters residuals
2147   Process3D(kCluster, 0, &fg, 1.e4);
2148   Process3Dlinked(kCluster, 1, &fg);
2149   fNRefFigures = 3;
2150   // Tracklet residual/pulls
2151   Process3D(kTrack , 0, &fg, 1.e4);
2152   Process3Dlinked(kTrack , 1, &fg);
2153   Process3D(kTrack , 2, &fg, 1.e4);
2154   Process3D(kTrack , 3, &fg);
2155   Process2D(kTrack , 4, &fg, 1.e3);
2156   fNRefFigures = 7;
2157   // TRDin residual/pulls
2158   Process3D(kTrackIn, 0, &fg, 1.e4);
2159   Process3Dlinked(kTrackIn, 1, &fg);
2160   Process3D(kTrackIn, 2, &fg, 1.e4);
2161   Process3D(kTrackIn, 3, &fg);
2162   Process2D(kTrackIn, 4, &fg, 1.e3);
2163   fNRefFigures = 11;
2164   // TRDout residual/pulls
2165   Process3D(kTrackOut, 0, &fg, 1.e3); // scale to fit - see PlotTrackOut
2166   Process3Dlinked(kTrackOut, 1, &fg);
2167   Process3D(kTrackOut, 2, &fg, 1.e4);
2168   Process3D(kTrackOut, 3, &fg);
2169   Process2D(kTrackOut, 4, &fg, 1.e3);
2170   fNRefFigures = 15;
2171
2172   if(!HasMCdata()) return kTRUE;
2173
2174
2175   //PROCESS MC RESIDUAL DISTRIBUTIONS
2176
2177   // CLUSTER Y RESOLUTION/PULLS
2178   Process3D(kMCcluster, 0, &fg, 1.e4);
2179   Process3Dlinked(kMCcluster, 1, &fg, 1.);
2180   fNRefFigures = 17;
2181
2182   // TRACKLET RESOLUTION/PULLS
2183   Process3D(kMCtracklet, 0, &fg, 1.e4); // y
2184   Process3Dlinked(kMCtracklet, 1, &fg, 1.);   // y pulls
2185   Process3D(kMCtracklet, 2, &fg, 1.e4); // z
2186   Process3D(kMCtracklet, 3, &fg, 1.);   // z pulls
2187   Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
2188   fNRefFigures = 21;
2189
2190   // TRACK RESOLUTION/PULLS
2191   Process3Darray(kMCtrack, 0, &fg, 1.e4);   // y
2192   Process3DlinkedArray(kMCtrack, 1, &fg);   // y PULL
2193   Process3Darray(kMCtrack, 2, &fg, 1.e4);   // z
2194   Process3Darray(kMCtrack, 3, &fg);         // z PULL
2195   Process2Darray(kMCtrack, 4, &fg, 1.e3);   // phi
2196   Process2Darray(kMCtrack, 5, &fg);         // snp PULL
2197   Process2Darray(kMCtrack, 6, &fg, 1.e3);   // theta
2198   Process2Darray(kMCtrack, 7, &fg);         // tgl PULL
2199   Process3Darray(kMCtrack, 8, &fg, 1.e2);   // pt resolution
2200   Process3Darray(kMCtrack, 9, &fg);         // 1/pt pulls
2201   Process3Darray(kMCtrack, 10, &fg, 1.e2);  // p resolution
2202   fNRefFigures+=16;
2203
2204   // TRACK TRDin RESOLUTION/PULLS
2205   Process3D(kMCtrackIn, 0, &fg, 1.e4);// y resolution
2206   Process3Dlinked(kMCtrackIn, 1, &fg);      // y pulls
2207   Process3D(kMCtrackIn, 2, &fg, 1.e4);// z resolution
2208   Process3D(kMCtrackIn, 3, &fg);      // z pulls
2209   Process2D(kMCtrackIn, 4, &fg, 1.e3);// phi resolution
2210   Process2D(kMCtrackIn, 5, &fg);      // snp pulls
2211   Process2D(kMCtrackIn, 6, &fg, 1.e3);// theta resolution
2212   Process2D(kMCtrackIn, 7, &fg);      // tgl pulls
2213   Process3D(kMCtrackIn, 8, &fg, 1.e2);// pt resolution
2214   Process3D(kMCtrackIn, 9, &fg);      // 1/pt pulls
2215   Process3D(kMCtrackIn, 10, &fg, 1.e2);// p resolution
2216   fNRefFigures+=8;
2217
2218   // TRACK TRDout RESOLUTION/PULLS
2219   Process3D(kMCtrackOut, 0, &fg, 1.e4);// y resolution
2220   Process3Dlinked(kMCtrackOut, 1, &fg);      // y pulls
2221   Process3D(kMCtrackOut, 2, &fg, 1.e4);// z resolution
2222   Process3D(kMCtrackOut, 3, &fg);      // z pulls
2223   Process2D(kMCtrackOut, 4, &fg, 1.e3);// phi resolution
2224   Process2D(kMCtrackOut, 5, &fg);      // snp pulls
2225   Process2D(kMCtrackOut, 6, &fg, 1.e3);// theta resolution
2226   Process2D(kMCtrackOut, 7, &fg);      // tgl pulls
2227   Process3D(kMCtrackOut, 8, &fg, 1.e2);// pt resolution
2228   Process3D(kMCtrackOut, 9, &fg);      // 1/pt pulls
2229   Process3D(kMCtrackOut, 10, &fg, 1.e2);// p resolution
2230   fNRefFigures+=8;
2231
2232   return kTRUE;
2233 }
2234
2235
2236 //________________________________________________________
2237 void AliTRDresolution::Terminate(Option_t *opt)
2238 {
2239   AliTRDrecoTask::Terminate(opt);
2240   if(HasPostProcess()) PostProcess();
2241 }
2242
2243 //________________________________________________________
2244 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2245 {
2246 // Helper function to avoid duplication of code
2247 // Make first guesses on the fit parameters
2248
2249   // find the intial parameters of the fit !! (thanks George)
2250   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2251   Double_t sum = 0.;
2252   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2253   f->SetParLimits(0, 0., 3.*sum);
2254   f->SetParameter(0, .9*sum);
2255   Double_t rms = h->GetRMS();
2256   f->SetParLimits(1, -rms, rms);
2257   f->SetParameter(1, h->GetMean());
2258
2259   f->SetParLimits(2, 0., 2.*rms);
2260   f->SetParameter(2, rms);
2261   if(f->GetNpar() <= 4) return;
2262
2263   f->SetParLimits(3, 0., sum);
2264   f->SetParameter(3, .1*sum);
2265
2266   f->SetParLimits(4, -.3, .3);
2267   f->SetParameter(4, 0.);
2268
2269   f->SetParLimits(5, 0., 1.e2);
2270   f->SetParameter(5, 2.e-1);
2271 }
2272
2273 //________________________________________________________
2274 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand)
2275 {
2276 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2277 //  - y reziduals/pulls
2278
2279   TObjArray *arr = new TObjArray(2);
2280   arr->SetName(name); arr->SetOwner();
2281   TH1 *h(NULL); char hname[100], htitle[300];
2282
2283   // tracklet resolution/pull in y direction
2284   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
2285   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2286   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2287     Int_t nybins=fgkNresYsegm[fSegmentLevel];
2288     if(expand) nybins*=2;
2289     h = new TH3S(hname, htitle,
2290                  48, -.48, .48,            // phi
2291                  60, -fDyRange, fDyRange,  // dy
2292                  nybins, -0.5, nybins-0.5);// segment
2293   } else h->Reset();
2294   arr->AddAt(h, 0);
2295   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
2296   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2297   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2298     h = new TH3S(hname, htitle, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2299   } else h->Reset();
2300   arr->AddAt(h, 1);
2301
2302   return arr;
2303 }
2304
2305 //________________________________________________________
2306 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2307 {
2308 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2309 //  - y reziduals/pulls
2310 //  - z reziduals/pulls
2311 //  - phi reziduals
2312   TObjArray *arr = BuildMonitorContainerCluster(name, expand);
2313   arr->Expand(5);
2314   TH1 *h(NULL); char hname[100], htitle[300];
2315
2316   // tracklet resolution/pull in z direction
2317   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
2318   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
2319   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2320     h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
2321   } else h->Reset();
2322   arr->AddAt(h, 2);
2323   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
2324   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
2325   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2326     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2327     h->GetZaxis()->SetBinLabel(1, "no RC");
2328     h->GetZaxis()->SetBinLabel(2, "RC");
2329   } else h->Reset();
2330   arr->AddAt(h, 3);
2331
2332   // tracklet to track phi resolution
2333   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
2334   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
2335   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2336     h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
2337   } else h->Reset();
2338   arr->AddAt(h, 4);
2339
2340   return arr;
2341 }
2342
2343 //________________________________________________________
2344 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2345 {
2346 // Build performance histograms for AliExternalTrackParam.vs MC
2347 //  - y resolution/pulls
2348 //  - z resolution/pulls
2349 //  - phi resolution, snp pulls
2350 //  - theta resolution, tgl pulls
2351 //  - pt resolution, 1/pt pulls, p resolution
2352
2353   TObjArray *arr = BuildMonitorContainerTracklet(name);
2354   arr->Expand(11);
2355   TH1 *h(NULL); char hname[100], htitle[300];
2356   TAxis *ax(NULL);
2357
2358   // snp pulls
2359   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
2360   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
2361   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2362     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2363   } else h->Reset();
2364   arr->AddAt(h, 5);
2365
2366   // theta resolution
2367   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
2368   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2369   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2370     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2371   } else h->Reset();
2372   arr->AddAt(h, 6);
2373   // tgl pulls
2374   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
2375   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
2376   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2377     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2378   } else h->Reset();
2379   arr->AddAt(h, 7);
2380
2381   const Int_t kNpt(14);
2382   const Int_t kNdpt(150);
2383   const Int_t kNspc = 2*AliPID::kSPECIES+1;
2384   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
2385   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2386   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2387   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
2388   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
2389
2390   // Pt resolution
2391   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
2392   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2393   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2394     h = new TH3S(hname, htitle,
2395                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2396     ax = h->GetZaxis();
2397     for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2398   } else h->Reset();
2399   arr->AddAt(h, 8);
2400   // 1/Pt pulls
2401   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
2402   snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
2403   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2404     h = new TH3S(hname, htitle,
2405                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2406     ax = h->GetZaxis();
2407     for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2408   } else h->Reset();
2409   arr->AddAt(h, 9);
2410   // P resolution
2411   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
2412   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2413   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2414     h = new TH3S(hname, htitle,
2415                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2416     ax = h->GetZaxis();
2417     for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2418   } else h->Reset();
2419   arr->AddAt(h, 10);
2420
2421   return arr;
2422 }
2423
2424
2425 //________________________________________________________
2426 TObjArray* AliTRDresolution::Histos()
2427 {
2428   //
2429   // Define histograms
2430   //
2431
2432   if(fContainer) return fContainer;
2433
2434   fContainer  = new TObjArray(kNviews);
2435   //fContainer->SetOwner(kTRUE);
2436   TH1 *h(NULL);
2437   TObjArray *arr(NULL);
2438
2439   // binnings for plots containing momentum or pt
2440   const Int_t kNpt(14), kNphi(48), kNdy(60);
2441   Float_t lPhi=-.48, lDy=-.3, lPt=0.1;
2442   Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
2443   for(Int_t i=0; i<kNphi+1; i++,lPhi+=.02) binsPhi[i]=lPhi;
2444   for(Int_t i=0; i<kNdy+1; i++,lDy+=.01) binsDy[i]=lDy;
2445   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2446
2447   // cluster to track residuals/pulls
2448   fContainer->AddAt(arr = new TObjArray(2), kCharge);
2449   arr->SetName("Charge");
2450   if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
2451     h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
2452     h->GetXaxis()->SetTitle("dx/dx_{0}");
2453     h->GetYaxis()->SetTitle("x_{d} [cm]");
2454     h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
2455   } else h->Reset();
2456   arr->AddAt(h, 0);
2457
2458   // cluster to track residuals/pulls
2459   fContainer->AddAt(BuildMonitorContainerCluster("Cl"), kCluster);
2460   // tracklet to TRD track
2461   fContainer->AddAt(BuildMonitorContainerTracklet("Trk", kTRUE), kTrack);
2462   // tracklet to TRDin
2463   fContainer->AddAt(BuildMonitorContainerTracklet("TrkIN", kTRUE), kTrackIn);
2464   // tracklet to TRDout
2465   fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2466
2467
2468   // Resolution histos
2469   if(!HasMCdata()) return fContainer;
2470
2471   // cluster resolution
2472   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
2473
2474   // tracklet resolution
2475   fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
2476
2477   // track resolution
2478   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
2479   arr->SetName("MCtrk");
2480   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
2481
2482   // TRDin TRACK RESOLUTION
2483   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkIN"), kMCtrackIn);
2484
2485   // TRDout TRACK RESOLUTION
2486   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
2487
2488   return fContainer;
2489 }
2490
2491 //________________________________________________________
2492 Bool_t AliTRDresolution::Load(const Char_t *file, const Char_t *dir)
2493 {
2494 // Custom load function. Used to guess the segmentation level of the data.
2495
2496   if(!AliTRDrecoTask::Load(file, dir)) return kFALSE;
2497
2498   // look for cluster residual plot - always available
2499   TH3S* h3((TH3S*)((TObjArray*)fContainer->At(kClToTrk))->At(0));
2500   Int_t segmentation(h3->GetNbinsZ()/2);
2501   if(segmentation==fgkNresYsegm[0]){ // default segmentation. Nothing to do
2502     return kTRUE;
2503   } else if(segmentation==fgkNresYsegm[1]){ // stack segmentation.
2504     SetSegmentationLevel(1);
2505   } else if(segmentation==fgkNresYsegm[2]){ // detector segmentation.
2506     SetSegmentationLevel(2);
2507   } else {
2508     AliError(Form("Unknown segmentation [%d].", h3->GetNbinsZ()));
2509     return kFALSE;
2510   }
2511
2512   AliDebug(2, Form("Segmentation set to level \"%s\"", fgkResYsegmName[fSegmentLevel]));
2513   return kTRUE;
2514 }
2515
2516 //________________________________________________________
2517 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
2518 {
2519 // Robust function to process sigma/mean for 2D plot dy(x)
2520 // For each x bin a gauss fit is performed on the y projection and the range
2521 // with the minimum chi2/ndf is choosen
2522
2523   if(!h2) {
2524     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
2525     return kFALSE;
2526   }
2527   if(!Int_t(h2->GetEntries())){
2528     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
2529     return kFALSE;
2530   }
2531   if(!g || !g[0]|| !g[1]) {
2532     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
2533     return kFALSE;
2534   }
2535   // prepare
2536   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
2537   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
2538   TF1 f("f", "gaus", ymin, ymax);
2539   Int_t n(0);
2540   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2541   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2542   TH1D *h(NULL);
2543   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
2544   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
2545
2546
2547   // do actual loop
2548   Float_t chi2OverNdf(0.);
2549   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
2550     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
2551     ymin = ay->GetXmin(); ymax = ay->GetXmax();
2552
2553     h = h2->ProjectionY("py", ix, ix);
2554     if((n=(Int_t)h->GetEntries())<stat){
2555       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
2556       continue;
2557     }
2558     // looking for a first order mean value
2559     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
2560     h->Fit(&f, "QNW");
2561     chi2OverNdf = f.GetChisquare()/f.GetNDF();
2562     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
2563     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
2564     ey  = f.GetParError(1);
2565     sy = f.GetParameter(2); esy = f.GetParError(2);
2566 //     // looking for the best chi2/ndf value
2567 //     while(ymin<y0 && ymax>y1){
2568 //       f.SetParameter(1, y);
2569 //       f.SetParameter(2, sy);
2570 //       h->Fit(&f, "QNW", "", y0, y1);
2571 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
2572 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
2573 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
2574 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
2575 //         sy = f.GetParameter(2); esy = f.GetParError(2);
2576 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
2577 //       }
2578 //       y0-=dy; y1+=dy;
2579 //     }
2580     g[0]->SetPoint(np, x, y);
2581     g[0]->SetPointError(np, ex, ey);
2582     g[1]->SetPoint(np, x, sy);
2583     g[1]->SetPointError(np, ex, esy);
2584     np++;
2585   }
2586   return kTRUE;
2587 }
2588
2589
2590 //________________________________________________________
2591 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2592 {
2593   //
2594   // Do the processing
2595   //
2596
2597   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
2598   Int_t n = 0;
2599   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2600   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2601   if(Int_t(h2->GetEntries())){
2602     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2603   } else {
2604     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2605     fIdxPlot++;
2606     return kTRUE;
2607   }
2608
2609   const Int_t kINTEGRAL=1;
2610   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2611     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2612     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2613     TH1D *h = h2->ProjectionY(pn, abin, bbin);
2614     if((n=(Int_t)h->GetEntries())<150){
2615       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2616       continue;
2617     }
2618     h->Fit(f, "QN");
2619     Int_t ip = g[0]->GetN();
2620     AliDebug(4, Form("  x_%d[%f] range[%d %d] stat[%d] M[%f] Sgm[%f]", ip, x, abin, bbin, n, f->GetParameter(1), f->GetParameter(2)));
2621     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2622     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2623     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2624     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
2625 /*
2626     g[0]->SetPoint(ip, x, k*h->GetMean());
2627     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2628     g[1]->SetPoint(ip, x, k*h->GetRMS());
2629     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2630   }
2631   fIdxPlot++;
2632   return kTRUE;
2633 }
2634
2635 //________________________________________________________
2636 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
2637 {
2638   //
2639   // Do the processing
2640   //
2641
2642   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2643
2644   // retrive containers
2645   TH2I *h2(NULL);
2646   if(idx<0){
2647     if(!(h2= (TH2I*)(fContainer->At(plot)))) return kFALSE;
2648   } else{
2649     TObjArray *a0(NULL);
2650     if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2651     if(!(h2=(TH2I*)a0->At(idx))) return kFALSE;
2652   }
2653   if(Int_t(h2->GetEntries())){
2654     AliDebug(2, Form("p[%d] idx[%d] : h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
2655   } else {
2656     AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2657     return kFALSE;
2658   }
2659
2660   TGraphErrors *g[2];
2661   if(gidx<0) gidx=idx;
2662   if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
2663
2664   if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
2665
2666   return Process(h2, f, k, g);
2667 }
2668
2669 //________________________________________________________
2670 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2671 {
2672   //
2673   // Do the processing
2674   //
2675
2676   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2677
2678   // retrive containers
2679   TH3S *h3(NULL);
2680   if(idx<0){
2681     if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2682   } else{
2683     TObjArray *a0(NULL);
2684     if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2685     if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2686   }
2687   if(Int_t(h3->GetEntries())){
2688     AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2689   } else {
2690     AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2691     return kFALSE;
2692   }
2693
2694   TObjArray *gm, *gs;
2695   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2696   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2697   TGraphErrors *g[2];
2698
2699   TAxis *az = h3->GetZaxis();
2700   for(Int_t iz(0); iz<gm->GetEntriesFast(); iz++){
2701     if(!(g[0] = (TGraphErrors*)gm->At(iz))) return kFALSE;
2702     if(!(g[1] = (TGraphErrors*)gs->At(iz))) return kFALSE;
2703     az->SetRange(iz+1, iz+1);
2704     if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2705   }
2706
2707   return kTRUE;
2708 }
2709
2710
2711 //________________________________________________________
2712 Bool_t AliTRDresolution::Process3Dlinked(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2713 {
2714   //
2715   // Do the processing
2716   //
2717
2718   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2719
2720   // retrive containers
2721   TH3S *h3(NULL);
2722   if(idx<0){
2723     if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2724   } else{
2725     TObjArray *a0(NULL);
2726     if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2727     if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2728   }
2729   if(Int_t(h3->GetEntries())){
2730     AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2731   } else {
2732     AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2733     return kFALSE;
2734   }
2735
2736   TObjArray *gm, *gs;
2737   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2738   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2739   TGraphErrors *g[2];
2740
2741   if(!(g[0] = (TGraphErrors*)gm->At(0))) return kFALSE;
2742   if(!(g[1] = (TGraphErrors*)gs->At(0))) return kFALSE;
2743   if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2744
2745   if(!(g[0] = (TGraphErrors*)gm->At(1))) return kFALSE;
2746   if(!(g[1] = (TGraphErrors*)gs->At(1))) return kFALSE;
2747   if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2748
2749   return kTRUE;
2750 }
2751
2752
2753 //________________________________________________________
2754 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2755 {
2756   //
2757   // Do the processing
2758   //
2759
2760   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2761
2762   // retrive containers
2763   TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
2764   if(!h3) return kFALSE;
2765   AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2766
2767   TGraphAsymmErrors *gm;
2768   TGraphErrors *gs;
2769   if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
2770   if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
2771
2772   Float_t x(0.), r(0.), mpv(0.), xM(0.), xm(0.);
2773   TAxis *ay = h3->GetYaxis();
2774   for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
2775     ay->SetRange(iy, iy);
2776     x = ay->GetBinCenter(iy);
2777     TH2F *h2=(TH2F*)h3->Project3D("zx");
2778     TAxis *ax = h2->GetXaxis();
2779     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
2780       r = ax->GetBinCenter(ix);
2781       TH1D *h1 = h2->ProjectionY("py", ix, ix);
2782       if(h1->Integral()<50) continue;
2783       h1->Fit(f, "QN");
2784
2785       GetLandauMpvFwhm(f, mpv, xm, xM);
2786       Int_t ip = gm->GetN();
2787       gm->SetPoint(ip, x, k*mpv);
2788       gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2789       gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2790       gs->SetPointError(ip, 0., 0.);
2791     }
2792   }
2793
2794   return kTRUE;
2795 }
2796
2797 //________________________________________________________
2798 Bool_t AliTRDresolution::Process2Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2799 {
2800   //
2801   // Do the processing
2802   //
2803
2804   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2805
2806   // retrive containers
2807   TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2808   if(!arr) return kFALSE;
2809   AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2810
2811   TObjArray *gm, *gs;
2812   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2813   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2814
2815   TGraphErrors *g[2]; TH2I *h2(NULL); TObjArray *a0(NULL);
2816   for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2817     if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2818
2819     if(!(h2 = (TH2I*)a0->At(idx))) return kFALSE;
2820     if(Int_t(h2->GetEntries())){
2821       AliDebug(4, Form("   idx[%d] h[%s] %s", ia, h2->GetName(), h2->GetTitle()));
2822     } else {
2823       AliDebug(2, Form("   idx[%d] : Missing entries.", ia));
2824       continue;
2825     }
2826
2827     if(!(g[0] = (TGraphErrors*)gm->At(ia))) return kFALSE;
2828     if(!(g[1] = (TGraphErrors*)gs->At(ia))) return kFALSE;
2829     if(!Process(h2, f, k, g)) return kFALSE;
2830   }
2831
2832   return kTRUE;
2833 }
2834
2835 //________________________________________________________
2836 Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2837 {
2838   //
2839   // Do the processing
2840   //
2841
2842   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2843   //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2844
2845   // retrive containers
2846   TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2847   if(!arr) return kFALSE;
2848   AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2849
2850   TObjArray *gm, *gs;
2851   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2852   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2853
2854   TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2855   Int_t in(0);
2856   for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2857     if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2858
2859     if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2860      if(Int_t(h3->GetEntries())){
2861        AliDebug(4, Form("   idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2862     } else {
2863       AliDebug(2, Form("   idx[%d] : Missing entries.", ia));
2864       continue;
2865     }
2866     TAxis *az = h3->GetZaxis();
2867     for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
2868       if(in >= gm->GetEntriesFast()) break;
2869       if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2870       if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2871       az->SetRange(iz, iz);
2872       if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2873     }
2874   }
2875   AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2876
2877   return kTRUE;
2878 }
2879
2880 //________________________________________________________
2881 Bool_t AliTRDresolution::Process3DlinkedArray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2882 {
2883   //
2884   // Do the processing
2885   //
2886
2887   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2888   //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2889
2890   // retrive containers
2891   TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2892   if(!arr) return kFALSE;
2893   AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2894
2895   TObjArray *gm, *gs;
2896   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2897   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2898
2899   TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2900   Int_t in(0);
2901   for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2902     if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2903     if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2904     if(Int_t(h3->GetEntries())){
2905       AliDebug(4, Form("   idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2906     } else {
2907       AliDebug(2, Form("   idx[%d] : Missing entries.", ia));
2908       continue;
2909     }
2910     if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2911     if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2912     if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2913     in++;
2914
2915     if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2916     if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2917     if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2918     in++;
2919   }
2920   AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2921
2922   return kTRUE;
2923 }
2924
2925 //________________________________________________________
2926 Bool_t AliTRDresolution::GetGraph(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, const Char_t *explain)
2927 {
2928   //
2929   // Get the graphs
2930   //
2931
2932   if(!fGraphS || !fGraphM) return kFALSE;
2933   // axis titles look up
2934   Int_t nref = 0;
2935   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
2936   UChar_t jdx = idx<0?0:idx;
2937   for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2938   Char_t **at = fAxTitle[nref];
2939
2940   // build legends if requiered
2941   TLegend *leg(NULL);
2942   if(kLEG){
2943     leg=new TLegend(.65, .6, .95, .9);
2944     leg->SetBorderSize(0);
2945     leg->SetFillStyle(0);
2946   }
2947   // build frame
2948   TH1S *h1(NULL);
2949   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), Form("%s %s;%s;%s", at[0], explain?explain:"", at[1], at[2]), 2, bb[0], bb[2]);
2950   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2951   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2952   // axis range
2953   TAxis *ax = h1->GetXaxis();
2954   ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2955   ax = h1->GetYaxis();
2956   ax->SetRangeUser(bb[1], bb[3]);
2957   ax->CenterTitle(); ax->SetTitleOffset(1.4);
2958   h1->Draw();
2959   // bounding box
2960   TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2961   b->SetFillStyle(3002);b->SetLineColor(0);
2962   b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2963   b->Draw();
2964
2965   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2966   if(!gm) return kFALSE;
2967   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2968   if(!gs) return kFALSE;
2969
2970   Int_t n(0), nPlots(0);
2971   if((n=gm->GetN())) {
2972     nPlots++;
2973     gm->Draw("pl"); if(leg) leg->AddEntry(gm, gm->GetTitle(), "pl");
2974     PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2975     PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2976   }
2977
2978   if((n=gs->GetN())){
2979     nPlots++;
2980     gs->Draw("pl"); if(leg) leg->AddEntry(gs, gs->GetTitle(), "pl");
2981     gs->Sort(&TGraph::CompareY);
2982     PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2983     PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2984     gs->Sort(&TGraph::CompareX);
2985   }
2986   if(!nPlots) return kFALSE;
2987   if(leg) leg->Draw();
2988   return kTRUE;
2989 }
2990
2991 //________________________________________________________
2992 Bool_t AliTRDresolution::GetGraphArray(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, Int_t n, Int_t *sel, const Char_t *explain)
2993 {
2994   //
2995   // Get the graphs
2996   //
2997
2998   if(!fGraphS || !fGraphM) return kFALSE;
2999
3000   // axis titles look up
3001   Int_t nref(0);
3002   for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
3003   nref+=idx;
3004   Char_t **at = fAxTitle[nref];
3005
3006   // build legends if requiered
3007   TLegend *legM(NULL), *legS(NULL);
3008   if(kLEG){
3009     legM=new TLegend(.35, .6, .65, .9);
3010     legM->SetHeader("Mean");
3011     legM->SetBorderSize(0);
3012     legM->SetFillStyle(0);
3013     legS=new TLegend(.65, .6, .95, .9);
3014     legS->SetHeader("Sigma");
3015     legS->SetBorderSize(0);
3016     legS->SetFillStyle(0);
3017   }
3018   // build frame
3019   TH1S *h1(NULL);
3020   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), Form("%s %s;%s;%s", at[0], explain?explain:"", at[1], at[2]), 2, bb[0], bb[2]);
3021   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
3022   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
3023   // axis range
3024   TAxis *ax = h1->GetXaxis();
3025   ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
3026   ax = h1->GetYaxis();
3027   ax->SetRangeUser(bb[1], bb[3]);
3028   ax->CenterTitle(); ax->SetTitleOffset(1.4);
3029   h1->Draw();
3030
3031   TGraphErrors *gm(NULL), *gs(NULL);
3032   TObjArray *a0(NULL), *a1(NULL);
3033   a0 = (TObjArray*)((TObjArray*)fGraphM->At(ip))->At(idx);
3034   a1 = (TObjArray*)((TObjArray*)fGraphS->At(ip))->At(idx);
3035   if(!n) n=a0->GetEntriesFast();
3036   AliDebug(4, Form("Graph : Ref[%d] Title[%s] Limits{x[%f %f] y[%f %f]} Comp[%d] Selection[%c]", nref, at[0], bb[0], bb[2], bb[1], bb[3], n, sel ? 'y' : 'n'));
3037   Int_t nn(0), nPlots(0);
3038   for(Int_t is(0), is0(0); is<n; is++){
3039     is0 = sel ? sel[is] : is;
3040     if(!(gs =  (TGraphErrors*)a1->At(is0))) return kFALSE;
3041     if(!(gm =  (TGraphErrors*)a0->At(is0))) return kFALSE;
3042
3043     if((nn=gs->GetN())){
3044       nPlots++;
3045       gs->Draw("pc");
3046       if(legS){
3047         //printf("LegEntry %s [%s]%s\n", at[0], gs->GetName(), gs->GetTitle());
3048         legS->AddEntry(gs, gs->GetTitle(), "pl");
3049       }
3050       gs->Sort(&TGraph::CompareY);
3051       PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[0]);
3052       PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[nn-1]);
3053       gs->Sort(&TGraph::CompareX);
3054     }
3055     if(gm->GetN()){
3056       nPlots++;
3057       gm->Draw("pc");
3058       if(legM){
3059         //printf("LegEntry %s [%s]%s\n", at[0], gm->GetName(), gm->GetTitle());
3060         legM->AddEntry(gm, gm->GetTitle(), "pl");
3061       }
3062       PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
3063       PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
3064     }
3065   }
3066   if(!nPlots) return kFALSE;
3067   if(kLEG){
3068     legM->Draw();
3069     legS->Draw();
3070   }
3071   return kTRUE;
3072 }
3073
3074 //____________________________________________________________________
3075 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
3076 {
3077 //
3078 // Fit track with a staight line using the "np" clusters stored in the array "points".
3079 // The following particularities are stored in the clusters from points:
3080 //   1. pad tilt as cluster charge
3081 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
3082 // The parameters of the straight line fit are stored in the array "param" in the following order :
3083 //     param[0] - x0 reference radial position
3084 //     param[1] - y0 reference r-phi position @ x0
3085 //     param[2] - z0 reference z position @ x0
3086 //     param[3] - slope dy/dx
3087 //     param[4] - slope dz/dx
3088 //
3089 // Attention :
3090 // Function should be used to refit tracks for B=0T
3091 //
3092
3093   if(np<40){
3094     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
3095     return kFALSE;
3096   }
3097   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
3098
3099   Double_t x0(0.);
3100   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
3101   x0/=Float_t(np);
3102
3103   Double_t x, y, z, dx, tilt(0.);
3104   for(Int_t ip(0); ip<np; ip++){
3105     x = points[ip].GetX(); z = points[ip].GetZ();
3106     dx = x - x0;
3107     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3108   }
3109   if(zfitter.Eval() != 0) return kFALSE;
3110
3111   Double_t z0    = zfitter.GetParameter(0);
3112   Double_t dzdx  = zfitter.GetParameter(1);
3113   for(Int_t ip(0); ip<np; ip++){
3114     if(points[ip].GetClusterType()) continue;
3115     x    = points[ip].GetX();
3116     dx   = x - x0;
3117     y    = points[ip].GetY();
3118     z    = points[ip].GetZ();
3119     tilt = points[ip].GetCharge();
3120     y -= tilt*(-dzdx*dx + z - z0);
3121     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3122     yfitter.AddPoint(&dx, y, 1.);
3123   }
3124   if(yfitter.Eval() != 0) return kFALSE;
3125   Double_t y0   = yfitter.GetParameter(0);
3126   Double_t dydx = yfitter.GetParameter(1);
3127
3128   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3129   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
3130   return kTRUE;
3131 }
3132
3133 //____________________________________________________________________
3134 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3135 {
3136 //
3137 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3138 // See function FitTrack for the data stored in the "clusters" array
3139
3140 // The parameters of the straight line fit are stored in the array "param" in the following order :
3141 //     par[0] - x0 reference radial position
3142 //     par[1] - y0 reference r-phi position @ x0
3143 //     par[2] - slope dy/dx
3144 //
3145 // Attention :
3146 // Function should be used to refit tracks for B=0T
3147 //
3148
3149   TLinearFitter yfitter(2, "pol1");
3150
3151   // grep data for tracklet
3152   Double_t x0(0.), x[60], y[60], dy[60];
3153   Int_t nly(0);
3154   for(Int_t ip(0); ip<np; ip++){
3155     if(points[ip].GetClusterType()) continue;
3156     if(points[ip].GetVolumeID() != ly) continue;
3157     Float_t xt(points[ip].GetX())
3158            ,yt(param[1] + param[3] * (xt - param[0]));
3159     x[nly] = xt;
3160     y[nly] = points[ip].GetY();
3161     dy[nly]= y[nly]-yt;
3162     x0    += xt;
3163     nly++;
3164   }
3165   if(nly<10){
3166     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3167     return kFALSE;
3168   }
3169   // set radial reference for fit
3170   x0 /= Float_t(nly);
3171
3172   // find tracklet core
3173   Double_t mean(0.), sig(1.e3);
3174   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3175
3176   // simple cluster error parameterization
3177   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3178
3179   // fit tracklet core
3180   for(Int_t jly(0); jly<nly; jly++){
3181     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3182     Double_t dx(x[jly]-x0);
3183     yfitter.AddPoint(&dx, y[jly], 1.);
3184   }
3185   if(yfitter.Eval() != 0) return kFALSE;
3186   par[0] = x0;
3187   par[1] = yfitter.GetParameter(0);
3188   par[2] = yfitter.GetParameter(1);
3189   return kTRUE;
3190 }
3191
3192 //____________________________________________________________________
3193 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3194 {
3195 //
3196 // Global selection mechanism of tracksbased on cluster to fit residuals
3197 // The parameters are the same as used ni function FitTrack().
3198
3199   const Float_t kS(0.6), kM(0.2);
3200   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3201   Float_t dy, dz, s, m;
3202   for(Int_t ip(0); ip<np; ip++){
3203     if(points[ip].GetClusterType()) continue;
3204     Float_t x0(points[ip].GetX())
3205           ,y0(param[1] + param[3] * (x0 - param[0]))
3206           ,z0(param[2] + param[4] * (x0 - param[0]));
3207     dy=points[ip].GetY() - y0; h.Fill(dy);
3208     dz=points[ip].GetZ() - z0;
3209   }
3210   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3211   fg.SetParameter(1, 0.);
3212   fg.SetParameter(2, 2.e-2);
3213   h.Fit(&fg, "QN");
3214   m=fg.GetParameter(1); s=fg.GetParameter(2);
3215   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3216   return kTRUE;
3217 }
3218
3219 //________________________________________________________
3220 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3221 {
3222   //
3223   // Get the most probable value and the full width half mean
3224   // of a Landau distribution
3225   //
3226
3227   const Float_t dx = 1.;
3228   mpv = f->GetParameter(1);
3229   Float_t fx, max = f->Eval(mpv);
3230
3231   xm = mpv - dx;
3232   while((fx = f->Eval(xm))>.5*max){
3233     if(fx>max){
3234       max = fx;
3235       mpv = xm;
3236     }
3237     xm -= dx;
3238   }
3239
3240   xM += 2*(mpv - xm);
3241   while((fx = f->Eval(xM))>.5*max) xM += dx;
3242 }
3243
3244
3245 //________________________________________________________
3246 void AliTRDresolution::SetSegmentationLevel(Int_t l)
3247 {
3248 // Setting the segmentation level to "l"
3249   fSegmentLevel = l;
3250
3251   UShort_t const lNcomp[kNprojs] = {
3252     1,  1, //2,
3253     fgkNresYsegm[fSegmentLevel], 2, //2,
3254     2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3255     2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3256     2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3257   // MC
3258     fgkNresYsegm[fSegmentLevel], 2,          //2,
3259     fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3260     fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
3261     fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
3262     6*fgkNresYsegm[fSegmentLevel], 6*2, 6*2, 6*2, 6, 6, 6, 6, 6*11, 6*11, 6*11  //11
3263   };
3264   memcpy(fNcomp, lNcomp, kNprojs*sizeof(UShort_t));
3265
3266   Char_t const *lAxTitle[kNprojs][4] = {
3267     // Charge
3268     {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
3269   ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
3270     // Clusters to Kalman
3271   ,{"Cluster2Track residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3272   ,{"Cluster2Track  YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3273     // TRD tracklet to Kalman fit
3274   ,{"Tracklet2Track Y residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3275   ,{"Tracklet2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3276   ,{"Tracklet2Track Z residuals", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3277   ,{"Tracklet2Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3278   ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3279     // TRDin 2 first TRD tracklet
3280   ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3281   ,{"Tracklet2Track YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3282   ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3283   ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3284   ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3285     // TRDout 2 first TRD tracklet
3286   ,{"Tracklet2Track Y residuals @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3287   ,{"Tracklet2Track YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3288   ,{"Tracklet2Track Z residuals @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3289   ,{"Tracklet2Track Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3290   ,{"Tracklet2Track Phi residuals @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3291     // MC cluster
3292   ,{"MC Cluster Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3293   ,{"MC Cluster YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3294     // MC tracklet
3295   ,{"MC Tracklet Y resolution", "tg(#phi)", "y [#mum]",  "#sigma_{y}[#mum]"}
3296   ,{"MC Tracklet YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3297   ,{"MC Tracklet Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3298   ,{"MC Tracklet Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3299   ,{"MC Tracklet Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3300     // MC track TRDin
3301   ,{"MC Y resolution @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3302   ,{"MC YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3303   ,{"MC Z resolution @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3304   ,{"MC Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3305   ,{"MC #Phi resolution @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3306   ,{"MC SNP pulls @ TRDin", "tg(#phi)", "SNP", "#sigma_{snp}"}
3307   ,{"MC #Theta resolution @ TRDin", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3308   ,{"MC TGL pulls @ TRDin", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3309   ,{"MC P_{t} resolution @ TRDin", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
3310   ,{"MC 1/P_{t} pulls @ TRDin", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
3311   ,{"MC P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3312     // MC track TRDout
3313   ,{"MC Y resolution @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3314   ,{"MC YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3315   ,{"MC Z resolution @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3316   ,{"MC Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3317   ,{"MC #Phi resolution @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3318   ,{"MC SNP pulls @ TRDout", "tg(#phi)", "SNP", "#sigma_{snp}"}
3319   ,{"MC #Theta resolution @ TRDout", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3320   ,{"MC TGL pulls @ TRDout", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3321   ,{"MC P_{t} resolution @ TRDout", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
3322   ,{"MC 1/P_{t} pulls @ TRDout", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
3323   ,{"MC P resolution @ TRDout", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3324     // MC track in TRD
3325   ,{"MC Track Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3326   ,{"MC Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3327   ,{"MC Track Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3328   ,{"MC Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3329   ,{"MC Track #Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3330   ,{"MC Track SNP pulls", "tg(#phi)", "SNP", "#sigma_{snp}"}
3331   ,{"MC Track #Theta resolution", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3332   ,{"MC Track TGL pulls", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3333   ,{"MC P_{t} resolution", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
3334   ,{"MC 1/P_{t} pulls", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
3335   ,{"MC P resolution", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
3336   };
3337   memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));
3338 }
3339
3340 #include "TFile.h"
3341 //________________________________________________________
3342 Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3343 {
3344   if(!file){
3345     AliWarning("Use cluster position as in reconstruction.");
3346     SetLoadCorrection();
3347     return kTRUE;
3348   }
3349   TDirectory *cwd(gDirectory);
3350   TString fileList;
3351   FILE *filePtr = fopen(file, "rt");
3352   if(!filePtr){
3353     AliError(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3354     SetLoadCorrection();
3355     return kFALSE;
3356   }
3357   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3358   while(fileList.Gets(filePtr)){
3359     if(!TFile::Open(fileList.Data())) {
3360       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3361       continue;
3362     } else AliInfo(Form("\"%s\"", fileList.Data()));
3363
3364     TTree *tSys = (TTree*)gFile->Get("tSys");
3365     h2->SetDirectory(gDirectory); h2->Reset("ICE");
3366     tSys->Draw("det:t>>h2", "dx", "goff");
3367     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3368       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3369     }
3370     h2->SetDirectory(cwd);
3371     gFile->Close();
3372   }
3373   cwd->cd();
3374
3375   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
3376     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3377     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3378     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3379     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3380       printf("%03d|", idet);
3381       for(Int_t it(0); it<30; it++) printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3382       printf("\n");
3383     }
3384   }
3385   SetLoadCorrection();
3386   return kTRUE;
3387 }