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