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