PWG1/TRD/AliTRDresolution.cxx

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