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