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