PWG1/TRD/AliTRDresolution.cxx

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