PWG1/TRD/AliTRDresolution.cxx

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