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