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