PWG1/TRD/AliTRDresolution.cxx

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