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