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