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