PWGPP/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 #include <TStyle.h>
  53 #include <TROOT.h>
  54 #include <TObjArray.h>
  55 #include <TH3.h>
  56 #include <TH2.h>
  57 #include <TH1.h>
  58 #include <THnSparse.h>
  59 #include <TF1.h>
  60 #include <TCanvas.h>
  61 #include <TGaxis.h>
  62 #include <TBox.h>
  63 #include <TLegend.h>
  64 #include <TGraphErrors.h>
  65 #include <TGraphAsymmErrors.h>
  66 #include <TLinearFitter.h>
  67 #include <TMath.h>
  68 #include <TMatrixT.h>
  69 #include <TVectorT.h>
  70 #include <TTreeStream.h>
  71 #include <TGeoManager.h>
  72 #include <TDatabasePDG.h>
  73
  74 #include "AliPID.h"
  75 #include "AliLog.h"
  76 #include "AliESDtrack.h"
  77 #include "AliMathBase.h"
  78 #include "AliTrackPointArray.h"
  79
  80 #include "AliTRDresolution.h"
  81 #include "AliTRDgeometry.h"
  82 #include "AliTRDtransform.h"
  83 #include "AliTRDpadPlane.h"
  84 #include "AliTRDcluster.h"
  85 #include "AliTRDseedV1.h"
  86 #include "AliTRDtrackV1.h"
  87 #include "AliTRDReconstructor.h"
  88 #include "AliTRDrecoParam.h"
  89 #include "AliTRDpidUtil.h"
  90 #include "AliTRDinfoGen.h"
  91 #include "AliAnalysisManager.h"
  92 #include "info/AliTRDclusterInfo.h"
  93 #include "info/AliTRDeventInfo.h"
  94
  95 ClassImp(AliTRDresolution)
  96 ClassImp(AliTRDresolution::AliTRDresolutionProjection)
  97
  98 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
  99   Int_t(kNbunchCross)/*bc*/,
 100   180/*phi*/,
 101   50/*eta*/,
 102   50/*dy*/,
 103   40/*dphi*/,
 104   50/*dz*/,
 105   3/*chg*species*/,
 106   kNpt/*pt*/
 107 };  //! no of bins/projection
 108 Double_t const AliTRDresolution::fgkMin[kNdim] = {
 109   -1.5,
 110   -TMath::Pi(),
 111   -1.,
 112   -1.5,
 113   -10.,
 114   -2.5,
 115   -1.5,
 116   -0.5
 117 };    //! low limits for projections
 118 Double_t const AliTRDresolution::fgkMax[kNdim] = {
 119   1.5,
 120   TMath::Pi(),
 121   1.,
 122   1.5,
 123   10.,
 124   2.5,
 125   1.5,
 126   kNpt-0.5
 127 };    //! high limits for projections
 128 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
 129   "bunch cross",
 130   "#phi [rad]",
 131   "#eta",
 132   "#Deltay [cm]",
 133   "#Delta#phi [deg]",
 134   "#Deltaz [cm]",
 135   "chg*spec*rc",
 136   "bin_p_{t}"
 137 };  //! title of projection
 138
 139 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
 140     "Det2Cluster"
 141     ,"Cluster2Track"
 142     ,"Tracklet2Track"
 143     ,"Tracklet2TRDin"
 144     ,"Cluster2MC"
 145     ,"Tracklet2MC"
 146     ,"TRDin2MC"
 147     ,"TRD2MC"
 148 //    ,"Tracklet2TRDout"
 149 //    ,"TRDout2MC"
 150 };
 151 Float_t AliTRDresolution::fgPtBin[25];
 152
 153 //________________________________________________________
 154 AliTRDresolution::AliTRDresolution()
 155   :AliTRDrecoTask()
 156   ,fSteer(0)
 157   ,fIdxFrame(0)
 158   ,fPtThreshold(.3)
 159   ,fDyRange(0.75)
 160   ,fBCbinTOF(0)
 161   ,fBCbinFill(0)
 162   ,fBsign(kFALSE)
 163   ,fProj(NULL)
 164   ,fDBPDG(NULL)
 165   ,fCl(NULL)
 166   ,fMCcl(NULL)
 167 {
 168   //
 169   // Default constructor
 170   //
 171   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
 172   MakePtSegmentation();
 173   SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
 174 }
 175
 176 //________________________________________________________
 177 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
 178   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
 179   ,fSteer(0)
 180   ,fIdxFrame(0)
 181   ,fPtThreshold(.3)
 182   ,fDyRange(0.75)
 183   ,fBCbinTOF(0)
 184   ,fBCbinFill(0)
 185   ,fBsign(kFALSE)
 186   ,fProj(NULL)
 187   ,fDBPDG(NULL)
 188   ,fCl(NULL)
 189   ,fMCcl(NULL)
 190 {
 191   //
 192   // Default constructor
 193   //
 194
 195   InitFunctorList();
 196   MakePtSegmentation();
 197   SetProcesses(kTRUE, kTRUE, kTRUE, kTRUE);
 198   if(xchange){
 199     SetUseExchangeContainers();
 200     DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
 201     DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
 202   }
 203 }
 204
 205 //________________________________________________________
 206 AliTRDresolution::~AliTRDresolution()
 207 {
 208   //
 209   // Destructor
 210   //
 211   if (AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
 212   if(fProj){fProj->Delete(); delete fProj;}
 213   if(fCl){fCl->Delete(); delete fCl;}
 214   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 215 }
 216
 217
 218 //________________________________________________________
 219 void AliTRDresolution::UserCreateOutputObjects()
 220 {
 221   // spatial resolution
 222
 223   AliTRDrecoTask::UserCreateOutputObjects();
 224   if(UseExchangeContainers()) InitExchangeContainers();
 225 }
 226
 227 //________________________________________________________
 228 void AliTRDresolution::InitExchangeContainers()
 229 {
 230 // Init containers for subsequent tasks (AliTRDclusterResolution)
 231
 232   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
 233   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
 234   PostData(kClToTrk, fCl);
 235   PostData(kClToMC, fMCcl);
 236 }
 237
 238 //________________________________________________________
 239 void AliTRDresolution::UserExec(Option_t *opt)
 240 {
 241   //
 242   // Execution part
 243   //
 244
 245   if(fCl) fCl->Delete();
 246   if(fMCcl) fMCcl->Delete();
 247   AliTRDrecoTask::UserExec(opt);
 248 }
 249
 250 //________________________________________________________
 251 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
 252 {
 253 // Helper function to calculate pulls in the yz plane
 254 // using proper tilt rotation
 255 // Uses functionality defined by AliTRDseedV1.
 256
 257   return kTRUE;
 258 /*
 259   Double_t t2(tilt*tilt);
 260   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
 261
 262   // rotate along pad
 263   Double_t cc[3];
 264   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
 265   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
 266   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
 267   // do sqrt
 268   Double_t sqr[3]={0., 0., 0.};
 269   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
 270   Double_t invsqr[3]={0., 0., 0.};
 271   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
 272   Double_t tmp(dyz[0]);
 273   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
 274   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
 275   return kTRUE;
 276 */
 277 }
 278
 279 //________________________________________________________
 280 TH1* AliTRDresolution::DetCluster(const TObjArray *cls)
 281 {
 282   //
 283   // Plot the cluster distributions
 284   //
 285
 286   if(cls) fkClusters = cls;
 287   if(!fkClusters){
 288     AliDebug(4, "No Clusters defined.");
 289     return NULL;
 290   }
 291   Int_t ncl(0);
 292   if(!(ncl = fkClusters->GetEntriesFast())){
 293     AliDebug(1, "No RecPoints defined.");
 294     return NULL;
 295   }
 296   Int_t det(-1);
 297   AliTRDcluster *cl(NULL);
 298   for(Int_t icl(0); icl<ncl; icl++){
 299     if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
 300     det = cl->GetDetector(); break;
 301   }
 302   if(det<0){
 303     AliDebug(1, "No useful clusters defined.");
 304     return NULL;
 305   }
 306   THnSparse *H(NULL);
 307   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kDetector)))){
 308     AliWarning("No output container defined.");
 309     return NULL;
 310   }
 311   Int_t ly(AliTRDgeometry::GetLayer(det)),
 312         stk(AliTRDgeometry::GetStack(det));
 313   Double_t val[kNdim],
 314            alpha((0.5+AliTRDgeometry::GetSector(det))*AliTRDgeometry::GetAlpha()),
 315            cs(TMath::Cos(alpha)),
 316            sn(TMath::Sin(alpha));
 317   for(Int_t icl(0); icl<ncl; icl++){
 318     if(!(cl = (AliTRDcluster*)(*fkClusters)[icl])) continue;
 319     val[kBC]  = ly;
 320     val[kPhi] = TMath::ATan2(cl->GetX()*sn + cl->GetY()*cs, cl->GetX()*cs - cl->GetY()*sn);
 321     val[kEta] = (5-stk)*16-cl->GetPadRow()-1-(stk<3?4:0);
 322     val[kYrez]= fEvent->GetMultiplicity();
 323     val[4]    = TMath::Abs(cl->GetQ());
 324     val[5]    = cl->IsFivePad()?1:cl->GetNPads();
 325     H->Fill(val);
 326   }
 327   return NULL;
 328 }
 329
 330 //________________________________________________________
 331 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 332 {
 333   //
 334   // Plot the cluster distributions
 335   //
 336
 337   if(track) fkTrack = track;
 338   if(!fkTrack){
 339     AliDebug(4, "No Track defined.");
 340     return NULL;
 341   }
 342   if(fkESD && TMath::Abs(fkESD->GetTOFbc()) > 1){
 343     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
 344     return NULL;
 345   }
 346   if(fkESD && fPt<fPtThreshold){
 347     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
 348     return NULL;
 349   }
 350   THnSparse *H(NULL);
 351   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
 352     AliWarning("No output container defined.");
 353     return NULL;
 354   }
 355
 356   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 357   Double_t val[kNdim+2],
 358            alpha(0.), cs(-2.), sn(0.); //Float_t exb, vd, t0, s2, dl, dt;
 359   TObjArray     *clInfoArr(NULL);
 360   AliTRDseedV1  *fTracklet(NULL);
 361   AliTRDcluster *c(NULL), *cc(NULL);
 362   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 363     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 364     if(!fTracklet->IsOK()) continue;
 365     //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
 366     val[kBC]  = ily;
 367     if(cs<-1.){
 368       alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
 369       cs    = TMath::Cos(alpha);
 370       sn    = TMath::Sin(alpha);
 371     }
 372     val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
 373     Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
 374     val[kEta] = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
 375     val[kPt]  = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
 376     Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
 377     Int_t row0(-1);
 378     Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
 379     fTracklet->ResetClusterIter(kTRUE);
 380     while((c = fTracklet->NextCluster())){
 381       Float_t xc(c->GetX()),
 382               q(TMath::Abs(c->GetQ()));
 383       if(row0<0) row0 = c->GetPadRow();
 384
 385       val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
 386       val[kPrez] = fTracklet->GetX0()-xc;
 387       val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
 388       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 389       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 390       if(ic) val[kZrez] /= (ic*q);
 391       val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
 392       val[kNdim]   = fEvent->GetMultiplicity();
 393       val[kNdim+1] = c->IsFivePad()?1:c->GetNPads();
 394       H->Fill(val);
 395 /*      // tilt rotation of covariance for clusters
 396       Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
 397       cov[0] = (sy2+t2*sz2)*corr;
 398       cov[1] = tilt*(sz2 - sy2)*corr;
 399       cov[2] = (t2*sy2+sz2)*corr;
 400       // sum with track covariance
 401       cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
 402       Double_t dyz[2]= {dy[1], dz[1]};
 403       Pulls(dyz, cov, tilt);*/
 404
 405       // Get z-position with respect to anode wire
 406       Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
 407               zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
 408       Int_t istk = geo->GetStack(c->GetDetector());
 409       AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
 410       Float_t rowZ = pp->GetRow0();
 411       Float_t d  = rowZ - zt + pp->GetAnodeWireOffset();
 412       d -= ((Int_t)(2 * d)) / 2.0;
 413       if (d > 0.25) d  = 0.5 - d;
 414
 415       AliTRDclusterInfo *clInfo(NULL);
 416       clInfo = new AliTRDclusterInfo;
 417       clInfo->SetCluster(c);
 418       //Float_t covF[] = {cov[0], cov[1], cov[2]};
 419       clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
 420       clInfo->SetResolution(val[kYrez]);
 421       clInfo->SetAnisochronity(d);
 422       clInfo->SetDriftLength(val[kZrez]);
 423       clInfo->SetTilt(fTracklet->GetTilt());
 424       if(fCl) fCl->Add(clInfo);
 425       //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 426
 427       if(DebugLevel()>=2){
 428         if(!clInfoArr){
 429           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 430           clInfoArr->SetOwner(kFALSE);
 431         }
 432         clInfoArr->Add(clInfo);
 433       }
 434     }
 435     if(DebugLevel()>=2 && clInfoArr){
 436       ULong_t status = fkESD->GetStatus();
 437       (*DebugStream()) << "cluster"
 438         <<"status="  << status
 439         <<"clInfo.=" << clInfoArr
 440         << "\n";
 441       clInfoArr->Clear();
 442     }
 443   }
 444   if(clInfoArr) delete clInfoArr;
 445
 446   if(!track) return NULL;
 447   // special care for EVE usage
 448   TH1 *h(NULL);
 449   if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 450   return H->Projection(kYrez);
 451 }
 452
 453
 454 //________________________________________________________
 455 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
 456 {
 457 // Plot normalized residuals for tracklets to track.
 458 //
 459 // We start from the result that if X=N(|m|, |Cov|)
 460 // BEGIN_LATEX
 461 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
 462 // END_LATEX
 463 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
 464 // reference position.
 465   if(track) fkTrack = track;
 466   if(!fkTrack){
 467     AliDebug(4, "No Track defined.");
 468     return NULL;
 469   }
 470   if(fkESD && TMath::Abs(fkESD->GetTOFbc())>1){
 471     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
 472     //return NULL;
 473   }
 474   THnSparse *H(NULL);
 475   if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
 476     AliWarning("No output container defined.");
 477     return NULL;
 478   }
 479
 480   const Int_t ndim(kNdim+8);
 481   Double_t val[ndim],
 482            alpha(0.), cs(-2.), sn(0.);
 483   Float_t sz[AliTRDseedV1::kNtb], pos[AliTRDseedV1::kNtb];
 484   Int_t ntbGap[AliTRDseedV1::kNtb];
 485   AliTRDseedV1 *fTracklet(NULL);
 486   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
 487     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
 488     if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()) continue;
 489     val [kBC] = il;
 490     if(cs<-1.){
 491       alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha();
 492       cs    = TMath::Cos(alpha);
 493       sn    = TMath::Sin(alpha);
 494     }
 495     val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
 496     Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());//fTracklet->GetTgl();
 497     val[kEta] = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
 498
 499     val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fkTrack->Charge();// fSpecies;
 500     val[kPt]  = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fTracklet->GetMomentum());
 501     Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
 502              dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
 503              dydx(fTracklet->GetYfit(1)),
 504              tilt(fTracklet->GetTilt());
 505     // correct for tilt rotation
 506     val[kYrez] = dyt - dzt*tilt;
 507     val[kZrez] = fTracklet->IsRowCross()?(dzt + dyt*tilt):(fTracklet->GetdQdl()*3.e-4-1.5);
 508     dydx+= tilt*fTracklet->GetZref(1);
 509     val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
 510     val[kNdim] = fEvent?fEvent->GetMultiplicity():0;
 511     val[kNdim+1] = 1.e2*fTracklet->GetTBoccupancy()/AliTRDseedV1::kNtb;
 512     Int_t n = fTracklet->GetChargeGaps(sz, pos, ntbGap);
 513     val[kNdim+2] = 0.; for(Int_t igap(0); igap<n; igap++) val[kNdim+2] += sz[igap];
 514     for(Int_t ifill(0); ifill<3; ifill++){ val[kNdim+3+ifill]=0.;val[kNdim+4+ifill]=0.;}
 515     for(Int_t igap(0), ifill(0); igap<n&&ifill<2; igap++){
 516       if(ntbGap[igap]<2) continue;
 517       val[kNdim+3+ifill] = sz[igap];
 518       val[kNdim+4+ifill] = pos[igap];
 519       ifill++;
 520     }
 521     H->Fill(val);
 522
 523 //     // compute covariance matrix
 524 //     fTracklet->GetCovAt(x, cov);
 525 //     fTracklet->GetCovRef(covR);
 526 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 527 //     Double_t dyz[2]= {dy[1], dz[1]};
 528 //     Pulls(dyz, cov, tilt);
 529 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 530 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
 531
 532     if(DebugLevel()>=3){
 533       Bool_t rc(fTracklet->IsRowCross());
 534       UChar_t err(fTracklet->GetErrorMsg());
 535       Double_t x(fTracklet->GetX()),
 536                pt(fTracklet->GetPt()),
 537                yt(fTracklet->GetYref(0)),
 538                zt(fTracklet->GetZref(0)),
 539                phi(fTracklet->GetYref(1)),
 540                tht(fTracklet->GetZref(1));
 541       Int_t ncl(fTracklet->GetN()), det(fTracklet->GetDetector());
 542       (*DebugStream()) << "tracklet"
 543         <<"pt="  << pt
 544         <<"x="   << x
 545         <<"yt="  << yt
 546         <<"zt="  << zt
 547         <<"phi=" << phi
 548         <<"tht=" << tht
 549         <<"det=" << det
 550         <<"n="   << ncl
 551         <<"dy0=" << dyt
 552         <<"dz0=" << dzt
 553         <<"dy="  << val[kYrez]
 554         <<"dz="  << val[kZrez]
 555         <<"dphi="<< val[kPrez]
 556         <<"dQ  ="<< val[kNdim]
 557         <<"rc="  << rc
 558         <<"err=" << err
 559         << "\n";
 560     }
 561   }
 562   if(!track) return NULL;
 563   // special care for EVE usage
 564   TH1 *h(NULL);
 565   if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 566   return H->Projection(kYrez);
 567 }
 568
 569
 570 //________________________________________________________
 571 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
 572 {
 573 // Store resolution/pulls of Kalman before updating with the TRD information
 574 // at the radial position of the first tracklet. The following points are used
 575 // for comparison
 576 //  - the (y,z,snp) of the first TRD tracklet
 577 //  - the (y, z, snp, tgl, pt) of the MC track reference
 578 //
 579 // Additionally the momentum resolution/pulls are calculated for usage in the
 580 // PID calculation.
 581   //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
 582
 583   if(track) fkTrack = track;
 584   if(!fkTrack){
 585     AliDebug(4, "No Track defined.");
 586     return NULL;
 587   }
 588   //fkTrack->Print();
 589   TH1 *h(NULL); // EVE projection
 590   // check container
 591   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
 592   if(!H){
 593     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
 594     return NULL;
 595   }
 596   // check input track status
 597   AliExternalTrackParam *tin(NULL);
 598   if(!(tin = fkTrack->GetTrackIn())){
 599     AliError("Track did not entered TRD fiducial volume.");
 600     return NULL;
 601   }
 602   // check first tracklet
 603   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
 604   if(!fTracklet){
 605     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
 606     if(!track) return NULL;
 607     // special care for EVE usage
 608     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 609     return H->Projection(kYrez);
 610   }
 611   if(!fTracklet->IsOK() || !fTracklet->IsChmbGood()){
 612     AliDebug(3, "Tracklet or Chamber not OK. Skip track.");
 613     if(!track) return NULL;
 614     // special care for EVE usage
 615     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 616     return H->Projection(kYrez);
 617   }
 618   // check radial position
 619   Double_t x = tin->GetX();
 620   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
 621     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
 622     if(!track) return NULL;
 623     // special care for EVE usage
 624     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 625     return H->Projection(kYrez);
 626   }
 627 //  printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
 628
 629   Int_t bc(fkESD?fkESD->GetTOFbc()/2:0);
 630   const Double_t *parR(tin->GetParameter());
 631   Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
 632             phit(fTracklet->GetYfit(1)),
 633             tilt(fTracklet->GetTilt()),
 634             norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
 635
 636   // correct for tilt rotation
 637   Double_t dy  = dyt - dzt*tilt,
 638            dz  = dzt + dyt*tilt,
 639            dx  = dy/(parR[2]*norm-parR[3]*norm*tilt);
 640   phit       += tilt*parR[3];
 641   Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
 642
 643   Double_t val[kNdim+3];
 644   val[kBC]          = bc==0?0:(bc<0?-1.:1.);
 645   Double_t alpha = (0.5+AliTRDgeometry::GetSector(fTracklet->GetDetector()))*AliTRDgeometry::GetAlpha(),
 646            cs    = TMath::Cos(alpha),
 647            sn    = TMath::Sin(alpha);
 648   val[kPhi] = TMath::ATan2(fTracklet->GetX()*sn + fTracklet->GetY()*cs, fTracklet->GetX()*cs - fTracklet->GetY()*sn);
 649   Float_t tgl = fTracklet->GetZ()/fTracklet->GetX()/TMath::Sqrt(1.+fTracklet->GetY()*fTracklet->GetY()/fTracklet->GetX()/fTracklet->GetX());
 650   val[kEta] = -TMath::Log(TMath::Tan(0.5 *  (0.5*TMath::Pi() - TMath::ATan(tgl))));
 651   val[kYrez]        = dy;
 652   val[kZrez]        = fTracklet->IsRowCross()?dz:(fTracklet->GetdQdl()*5.e-4 - 2.5);
 653   val[kPrez]        = dphi*TMath::RadToDeg();
 654   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies; //printf("val[%d] = %4.2f\n", kSpeciesChgRC, val[kSpeciesChgRC]);
 655   val[kPt]          = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fPt);
 656   val[kNdim]        = GetPtBin(fTracklet->GetMomentum());
 657   val[kNdim+1]      = dx;
 658   val[kNdim+2]      = fEvent?fEvent->GetBunchFill():0;
 659   H->Fill(val);
 660   if(DebugLevel()>=3){
 661     (*DebugStream()) << "trackIn"
 662       <<"tracklet.="  << fTracklet
 663       <<"trackIn.="   << tin
 664       << "\n";
 665   }
 666
 667   if(!track) return NULL;
 668   // special care for EVE usage
 669     if((h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 670     return H->Projection(kYrez);
 671 }
 672
 673 /*
 674 //________________________________________________________
 675 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
 676 {
 677 // Store resolution/pulls of Kalman after last update with the TRD information
 678 // at the radial position of the first tracklet. The following points are used
 679 // for comparison
 680 //  - the (y,z,snp) of the first TRD tracklet
 681 //  - the (y, z, snp, tgl, pt) of the MC track reference
 682 //
 683 // Additionally the momentum resolution/pulls are calculated for usage in the
 684 // PID calculation.
 685
 686   if(track) fkTrack = track;
 687   return NULL;
 688 }
 689 */
 690 //________________________________________________________
 691 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 692 {
 693   //
 694   // Plot MC distributions
 695   //
 696
 697   if(!HasMCdata()) return NULL;
 698   if(track) fkTrack = track;
 699   if(!fkTrack){
 700     AliDebug(4, "No Track defined.");
 701     return NULL;
 702   }
 703   Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
 704
 705   THnSparse *H(NULL);
 706   if(!fContainer){
 707     AliWarning("No output container defined.");
 708     return NULL;
 709   }
 710   // retriev track characteristics
 711   Int_t pdg = fkMC->GetPDG(),
 712         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 713         sign(0),
 714 //        sgm[3],
 715         label(fkMC->GetLabel());
 716 //        fSegmentLevel(0);
 717   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 718   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 719   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 720
 721   TH1 *h(NULL);
 722   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 723   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
 724   UChar_t s;
 725   Double_t x, y, z, pt, dydx, dzdx, dzdl;
 726   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 727   Double_t covR[7]/*, cov[3]*/;
 728
 729 /*  if(DebugLevel()>=3){
 730     // get first detector
 731     Int_t det = -1;
 732     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 733       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 734       det = fTracklet->GetDetector();
 735       break;
 736     }
 737     if(det>=0){
 738       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
 739       Double_t m(-1.);
 740       m = fkTrack->GetMass();
 741       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
 742         (*DebugStream()) << "MCkalman"
 743           << "pdg=" << pdg
 744           << "det=" << det
 745           << "x="   << &X
 746           << "y="   << &Y
 747           << "z="   << &Z
 748           << "dx="  << &dX
 749           << "dy="  << &dY
 750           << "dz="  << &dZ
 751           << "pt="  << &vPt
 752           << "dpt=" << &dPt
 753           << "bgt=" << &budget
 754           << "cov=" << &cCOV
 755           << "\n";
 756       }
 757     }
 758   }*/
 759   AliTRDcluster *c(NULL);
 760   Double_t val[kNdim+1];
 761   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 762     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
 763        !fTracklet->IsOK())*/ continue;
 764
 765     x= x0 = fTracklet->GetX();
 766     Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
 767     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue;
 768
 769     // MC track position at reference radial position
 770     dx  = x0 - x;
 771     Float_t ymc = y0 - dx*dydx0;
 772     Float_t zmc = z0 - dx*dzdx0;
 773     //phi -= TMath::Pi();
 774
 775     val[kBC]  = ily;
 776     val[kPhi] = phi;
 777     val[kEta] = eta;
 778     val[kSpeciesChgRC]= rc?0.:sign*sIdx;
 779     val[kPt]  = pt0<0.8?0:(pt0<1.5?1:2);//GetPtBin(pt0);
 780     Double_t tilt(fTracklet->GetTilt());
 781 //             ,t2(tilt*tilt)
 782 //             ,corr(1./(1. + t2))
 783 //             ,cost(TMath::Sqrt(corr));
 784
 785     AliExternalTrackParam *tin(fkTrack->GetTrackIn());
 786     if(ily==0 && tin){ // trackIn residuals
 787       // check radial position
 788       if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
 789       else{
 790         val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
 791         val[kYrez]        = tin->GetY()-ymc;
 792         val[kZrez]        = rc?(tin->GetZ()-zmc):(fTracklet->GetdQdl()*1.8e-4 - 0.9);
 793         val[kPrez]        = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
 794         if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
 795       }
 796     }
 797     //if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
 798
 799     // track residuals
 800     dydx = fTracklet->GetYref(1);
 801     dzdx = fTracklet->GetZref(1);
 802     dzdl = fTracklet->GetTgl();
 803     y  = fTracklet->GetYref(0);
 804     dy = y - ymc;
 805     z  = fTracklet->GetZref(0);
 806     dz = z - zmc;
 807     pt = TMath::Abs(fTracklet->GetPt());
 808     fTracklet->GetCovRef(covR);
 809
 810     val[kYrez] = dy;
 811     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
 812     val[kZrez] = dz;
 813     val[kNdim] = 1.e2*(pt/pt0-1.);
 814     if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
 815 /*      // theta resolution/ tgl pulls
 816       Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
 817                 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
 818       ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
 819       ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
 820       // pt resolution  \\ 1/pt pulls \\ p resolution for PID
 821       Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
 822               p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
 823       ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
 824       ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
 825       ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
 826
 827     // Fill Debug stream for MC track
 828     if(DebugLevel()>=4){
 829       Int_t det(fTracklet->GetDetector());
 830       (*DebugStream()) << "MC"
 831         << "det="     << det
 832         << "pdg="     << pdg
 833         << "sgn="     << sign
 834         << "pt="      << pt0
 835         << "x="       << x0
 836         << "y="       << y0
 837         << "z="       << z0
 838         << "dydx="    << dydx0
 839         << "dzdx="    << dzdx0
 840         << "\n";
 841
 842       // Fill Debug stream for Kalman track
 843       (*DebugStream()) << "MCtrack"
 844         << "pt="      << pt
 845         << "x="       << x
 846         << "y="       << y
 847         << "z="       << z
 848         << "dydx="    << dydx
 849         << "dzdx="    << dzdx
 850         << "s2y="     << covR[0]
 851         << "s2z="     << covR[2]
 852         << "\n";
 853     }
 854
 855     // tracklet residuals
 856     dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
 857     dzdx = fTracklet->GetZfit(1);
 858     y  = fTracklet->GetYfit(0);
 859     dy = y - ymc;
 860     z  = fTracklet->GetZfit(0);
 861     dz = z - zmc;
 862     val[kYrez] = dy - dz*tilt;
 863     val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
 864     val[kZrez] = rc?(dz + dy*tilt):(fTracklet->GetdQdl()*3.e-4 - 1.5);
 865 //      val[kNdim] = pt/pt0-1.;
 866     if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
 867
 868
 869     // Fill Debug stream for tracklet
 870     if(DebugLevel()>=4){
 871       Float_t s2y = fTracklet->GetS2Y();
 872       Float_t s2z = fTracklet->GetS2Z();
 873       (*DebugStream()) << "MCtracklet"
 874         << "rc="    << rc
 875         << "x="     << x
 876         << "y="     << y
 877         << "z="     << z
 878         << "dydx="  << dydx
 879         << "s2y="   << s2y
 880         << "s2z="   << s2z
 881         << "\n";
 882     }
 883
 884     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
 885     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 886     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
 887
 888     H = (THnSparse*)fContainer->At(kMCcluster);
 889     val[kPt]  = TMath::ATan(dydx0)*TMath::RadToDeg();
 890     //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
 891     Int_t row0(-1);
 892     Float_t padCorr(tilt*fTracklet->GetPadLength());
 893     fTracklet->ResetClusterIter(kTRUE);
 894     while((c = fTracklet->NextCluster())){
 895       if(row0<0) row0 = c->GetPadRow();
 896       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
 897       y = c->GetY()  + padCorr*(c->GetPadRow() - row0);
 898       z = c->GetZ();
 899       dx = x0 - x;
 900       ymc= y0 - dx*dydx0;
 901       zmc= z0 - dx*dzdx0;
 902       dy = y - ymc;
 903       dz = z - zmc;
 904       val[kYrez] = dy - dz*tilt;
 905       val[kPrez] = dx;
 906       val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t  q(TMath::Abs(c->GetQ()));
 907       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 908       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
 909       if(ic) val[kZrez] /= (ic*q);
 910       if(H) H->Fill(val);
 911
 912
 913       // Fill calibration container
 914       Float_t d = zr0 - zmc;
 915       d -= ((Int_t)(2 * d)) / 2.0;
 916       if (d > 0.25) d  = 0.5 - d;
 917       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 918       clInfo->SetCluster(c);
 919       clInfo->SetMC(pdg, label);
 920       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
 921       clInfo->SetResolution(dy);
 922       clInfo->SetAnisochronity(d);
 923       clInfo->SetDriftLength(dx);
 924       clInfo->SetTilt(tilt);
 925       if(fMCcl) fMCcl->Add(clInfo);
 926       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 927       if(DebugLevel()>=5){
 928         if(!clInfoArr){
 929           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 930           clInfoArr->SetOwner(kFALSE);
 931         }
 932         clInfoArr->Add(clInfo);
 933       }
 934     }
 935     // Fill Debug Tree
 936     if(DebugLevel()>=5 && clInfoArr){
 937       (*DebugStream()) << "MCcluster"
 938         <<"clInfo.=" << clInfoArr
 939         << "\n";
 940       clInfoArr->Clear();
 941     }
 942   }
 943   if(clInfoArr) delete clInfoArr;
 944   return h;
 945 }
 946
 947
 948 //__________________________________________________________________________
 949 Int_t AliTRDresolution::GetPtBin(Float_t pt)
 950 {
 951 // Find pt bin according to local pt segmentation
 952   Int_t ipt(-1);
 953   while(ipt<24){
 954     if(pt<fgPtBin[ipt+1]) break;
 955     ipt++;
 956   }
 957   return ipt;
 958 }
 959
 960 //________________________________________________________
 961 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
 962 {
 963 // return mean number of entries/bin of histogram "h"
 964 // if option "opt" is given the following values are accepted:
 965 //   "<" : consider only entries less than "cut"
 966 //   ">" : consider only entries greater than "cut"
 967
 968   //Int_t dim(h->GetDimension());
 969   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
 970   Double_t sum(0.); Int_t n(0);
 971   for(Int_t ix(1); ix<=nbx; ix++)
 972     for(Int_t iy(1); iy<=nby; iy++)
 973       for(Int_t iz(1); iz<=nbz; iz++){
 974         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
 975         else{
 976           if(strcmp(opt, "<")==0) {
 977             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 978           } else if(strcmp(opt, ">")==0){
 979             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 980           } else {sum += h->GetBinContent(ix, iy, iz); n++;}
 981         }
 982       }
 983   return n>0?sum/n:0.;
 984 }
 985
 986 //________________________________________________________
 987 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
 988 {
 989 // Get range on Z axis such to avoid outliers
 990
 991   Double_t cnt[20000], c, m, s;
 992   Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
 993   for(Int_t ix(1); ix<=nx; ix++){
 994     for(Int_t iy(1); iy<=ny; iy++){
 995       if((c = h2->GetBinContent(ix, iy))<10) continue;
 996       cnt[nc++] = c;
 997       if(nc==20000) break;
 998     }
 999     if(nc==20000) break;
1000   }
1001   AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
1002   min = m-s; max = m+2.*s;
1003 }
1004
1005 //________________________________________________________
1006 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1007 {
1008   //
1009   // Get the reference figures
1010   //
1011
1012   if(!gPad){
1013     AliWarning("Please provide a canvas to draw results.");
1014     return kFALSE;
1015   }
1016 /*  Int_t selection[100], n(0), selStart(0); //
1017   Int_t ly0(0), dly(5);
1018   TList *l(NULL); TVirtualPad *pad(NULL); */
1019   switch(ifig){
1020   case 0:
1021     break;
1022   }
1023   AliWarning(Form("Reference plot [%d] missing result", ifig));
1024   return kFALSE;
1025 }
1026
1027
1028 //________________________________________________________
1029 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
1030 {
1031 // Build pt segments
1032   for(Int_t j(0); j<=24; j++){
1033     pt0+=(TMath::Exp(j*j*dpt)-1.);
1034     fgPtBin[j]=pt0;
1035   }
1036 }
1037
1038 //________________________________________________________
1039 void AliTRDresolution::MakeSummary()
1040 {
1041 // Build summary plots
1042
1043   if(!fProj){
1044     AliError("Missing results");
1045     return;
1046   }
1047   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1048   TObjArray *arr(NULL); TH2 *h2(NULL);
1049
1050   // cluster resolution
1051   // define palette
1052   gStyle->SetPalette(1);
1053   const Int_t nClViews(9);
1054   const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1055   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1056   const Float_t vClSignMin[2] = {2.6e2, 4.4e2},
1057                 vClSignMax[2] = {4.4e2, 6.2e2},
1058                 vClMin[nClViews] = {3.2e2, vClSignMin[Int_t(fBsign)], vClSignMin[Int_t(!fBsign)], 0., 0., 0., 0., 0., 3.2e2},
1059                 vClMax[nClViews] = {5.e2, vClSignMax[Int_t(fBsign)], vClSignMax[Int_t(!fBsign)], 0., 0., 0., 0., 0., 5.e2};
1060   const Int_t nTrkltViews(20);
1061   const Char_t *vTrkltName[nTrkltViews] = {
1062     "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1063     "TrkltRCZ",
1064     "TrkltPh", "TrkltPhn", "TrkltPhp",
1065     "TrkltQ", "TrkltQn", "TrkltQp",
1066     "TrkltQS", "TrkltQSn", "TrkltQSp",
1067     "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1068 //    "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0",  // electrons low pt
1069 /*    "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1",  // muons low pt
1070     "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2"  // pions low pt*/
1071   };
1072   const UChar_t vTrkltOpt[nTrkltViews] = {
1073     0, 0, 0, 1, 1, 1,
1074     0,
1075     0, 0, 0,
1076     0, 0, 0,
1077     0, 0, 0,
1078     0, 0, 2, 2
1079   };
1080   const Int_t nTrkInViews(5);
1081   const Char_t *vTrkInName[nTrkInViews][6] = {
1082     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1083     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1084     {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1085     {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1086     {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1087     //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1088   };
1089   const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1090   const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1091   const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1092   const Char_t *ttt[6] = {"#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltay) [cm]", "#sigma(#Deltaz) [cm]", "#sigma(#Delta#phi) [deg]", "#sigma(#Delta#phi) [deg]"};
1093
1094   const Int_t nTrkViews(27);
1095   const Char_t *vTrkName[nTrkViews] = {
1096     "TrkY", "TrkYn", "TrkYp",
1097     "TrkPh", "TrkPhn", "TrkPhp",
1098     "TrkDPt", "TrkDPtn", "TrkDPtp",
1099     "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0",  // electrons low pt
1100     "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1",  // muons low pt
1101     "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2"  // pions low pt
1102   };
1103   const Char_t *typName[] = {"", "MC"};
1104   const Int_t nx(2048), ny(1536);
1105
1106   if((arr = (TObjArray*)fProj->At(kDetector))){
1107     cOut = new TCanvas(Form("%s_DetOccupancy", GetName()), "Detector performance", 2*nx, 2*ny);
1108     cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1109     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1110       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1111         p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1112         if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d%dEn", ily, icen)))) continue;
1113         SetNormZ(h2, 1, 11, 1, -1, 10.);
1114         SetRangeZ(h2, 0., 150, -25.);
1115         h2->GetZaxis()->SetTitle("Rel. Det. Occup. [%]");
1116         h2->GetZaxis()->CenterTitle();
1117         h2->Draw("colz");
1118         MakeDetectorPlot(ily, "pad");
1119       }
1120     }
1121     cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1122     cOut = new TCanvas(Form("%s_DetCharge", GetName()), "Detector performance", nx, ny);
1123     cOut->Divide(3,2, 1.e-5, 1.e-5);
1124     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1125       p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1126       if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d_2D", ily)))) continue;
1127       SetNormZ(h2, 1, 11, 1, -1, 10.);
1128       SetRangeZ(h2, 0., 45., -15.);
1129       h2->GetZaxis()->SetTitle("Rel. Mean(q) [%]");
1130       h2->GetZaxis()->CenterTitle();
1131       h2->Draw("colz");
1132       MakeDetectorPlot(ily, "pad");
1133     }
1134     cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1135   }
1136   for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1137     if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
1138       for(Int_t iview(0); iview<nClViews; iview++){
1139         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1140         cOut->Divide(3,2, 1.e-5, 1.e-5);
1141         Int_t nplot(0);
1142         for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1143           p=cOut->cd(ily+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1144           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], ily)))) continue;
1145           nplot++;
1146           if(vClOpt[iview]==0) h2->Draw("colz");
1147           else if(vClOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [#mum]", vClMin[iview], vClMax[iview], 1.e4);
1148           if(iview<5) MakeDetectorPlot(ily);
1149         }
1150         if(nplot==AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1151         else delete cOut;
1152       }
1153     }
1154     // tracklet systematic
1155     if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
1156       for(Int_t iview(0); iview<nTrkltViews; iview++){
1157         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1158         cOut->Divide(3,2, 1.e-5, 1.e-5);
1159         Int_t nplot(0);
1160         for(Int_t iplot(0); iplot<6; iplot++){
1161           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1162           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1163           nplot++;
1164           if(vTrkltOpt[iview]==0) h2->Draw("colz");
1165           else if (vTrkltOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .4);
1166           else if (vTrkltOpt[iview]==2) DrawSigma(h2, "#sigma(occupancy) [%]", 10.5, 15.);
1167           MakeDetectorPlot(iplot);
1168         }
1169         if(nplot==6){
1170           cOut->Modified();cOut->Update();
1171           cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1172         }
1173         delete cOut;
1174       }
1175     }
1176     // trackIn systematic
1177     if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
1178       for(Int_t iview(0); iview<nTrkInViews; iview++){
1179         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1180         cOut->Divide(3,2, 1.e-5, 1.e-5);
1181         Int_t nplot(0);
1182         for(Int_t iplot(0); iplot<6; iplot++){
1183           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1184           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){
1185             AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot]));
1186             continue;
1187           }
1188           nplot++;
1189           if(vTrkInOpt[iview]==0) h2->Draw("colz");
1190           else DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1191           MakeDetectorPlot(0);
1192         }
1193         if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1194         else delete cOut;
1195       }
1196       // species
1197       const Float_t zmin[] = {1., 20., 25., 10., 5.},
1198                     zmax[] = {10., 38., 43., 28., 14.};
1199       cOut = new TCanvas(Form("%s_%sTrkInSpc", GetName(), typName[ityp]), "Track IN PID", Int_t(1.5*nx), Int_t(1.5*ny));
1200       cOut->Divide(5,3, 1.e-5, 1.e-5);
1201       Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1202       for(Int_t ich(0), ipad(1); ich<3; ich++){
1203         TH2 *h2s(NULL);
1204         if(!(h2s = (TH2*)arr->FindObject(Form("H%sTrkInY%sEn", typName[ityp], chName[ich])))) {
1205           AliInfo(Form("Missing H%sTrkIn%sEn", typName[ityp], chName[ich]));
1206           continue;
1207         }
1208         for(Int_t ispec(0); ispec<AliPID::kSPECIES; ispec++){
1209           p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1210           if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInY%s%dEn", typName[ityp], chName[ich], ispec)))) {
1211             AliInfo(Form("Missing H%sTrkIn%s%dEn", typName[ityp], chName[ich], ispec));
1212             continue;
1213           }
1214           nplot++;
1215           h2->Divide(h2, h2s, 1.e2);
1216           h2->SetContour(9);
1217           h2->GetZaxis()->SetRangeUser(zmin[ispec], zmax[ispec]);
1218           h2->GetZaxis()->SetTitle("Rel. Abundancy [%]");h2->GetZaxis()->CenterTitle();
1219           TString tit(h2->GetTitle()); h2->SetTitle(Form("%s :: Relative Abundancy", ((*(tit.Tokenize("::")))[0])->GetName()));
1220           h2->Draw("colz");
1221           MakeDetectorPlot(0);
1222         }
1223       }
1224       if(nplot==15) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1225       else delete cOut;
1226     }
1227   }
1228   // track MC systematic
1229   if((arr = (TObjArray*)fProj->At(kMCtrack))) {
1230     for(Int_t iview(0); iview<nTrkViews; iview++){
1231       cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1232       cOut->Divide(3,2, 1.e-5, 1.e-5);
1233       Int_t nplot(0);
1234       for(Int_t iplot(0); iplot<6; iplot++){
1235         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1236         if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1237         h2->Draw("colz"); nplot++;
1238       }
1239       if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1240       else delete cOut;
1241     }
1242   }
1243
1244
1245   gStyle->SetPalette(1);
1246 }
1247
1248 //________________________________________________________
1249 void AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1250 {
1251   // Draw error bars scaled with "scale" instead of content values
1252   //use range [m,M] if limits are specified
1253
1254   if(!h2) return;
1255   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1256   TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1257                 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1258                 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1259                 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1260   h2e->SetContour(9);
1261   TAxis *az(h2e->GetZaxis());
1262   if(M>m) az->SetRangeUser(m, M);
1263   az->CenterTitle();
1264   az->SetTitleOffset(1.5);
1265   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1266     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1267       if(h2->GetBinContent(ix, iy)<-100.) continue;
1268       Float_t v(scale*h2->GetBinError(ix, iy));
1269       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1270       h2e->SetBinContent(ix, iy, v);
1271     }
1272   }
1273   h2e->Draw("colz");
1274 }
1275
1276 //________________________________________________________
1277 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1278 {
1279 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1280 // The "range" vector should be initialized with 2 elements
1281 // Option "mod" can be any of
1282 //   - 0 : gaussian like distribution
1283 //   - 1 : tailed distribution
1284
1285   Int_t nx(h2->GetNbinsX())
1286       , ny(h2->GetNbinsY())
1287       , n(nx*ny);
1288   Double_t *data=new Double_t[n];
1289   for(Int_t ix(1), in(0); ix<=nx; ix++){
1290     for(Int_t iy(1); iy<=ny; iy++)
1291       data[in++] = h2->GetBinContent(ix, iy);
1292   }
1293   Double_t mean, sigm;
1294   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1295
1296   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1297   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1298   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1299   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1300   h1.FillN(n,data,0);
1301   delete [] data;
1302
1303   switch(mod){
1304   case 0:// gaussian distribution
1305   {
1306     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1307     h1.Fit(&fg, "QN");
1308     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1309     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1310     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1311     break;
1312   }
1313   case 1:// tailed distribution
1314   {
1315     Int_t bmax(h1.GetMaximumBin());
1316     Int_t jBinMin(1), jBinMax(100);
1317     for(Int_t ibin(bmax); ibin--;){
1318       if(h1.GetBinContent(ibin)<1.){
1319         jBinMin=ibin; break;
1320       }
1321     }
1322     for(Int_t ibin(bmax); ibin++;){
1323       if(h1.GetBinContent(ibin)<1.){
1324         jBinMax=ibin; break;
1325       }
1326     }
1327     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1328     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1329     break;
1330   }
1331   }
1332
1333   return;
1334 }
1335
1336 //________________________________________________________
1337 Bool_t AliTRDresolution::MakeProjectionDetector()
1338 {
1339 // Analyse cluster
1340   const Int_t kNcontours(9);
1341   const Int_t kNstat(100);
1342   if(fProj && fProj->At(kDetector)) return kTRUE;
1343   if(!fContainer){
1344     AliError("Missing data container.");
1345     return kFALSE;
1346   }
1347   THnSparse *H(NULL);
1348   if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1349     AliInfo(Form("Missing/Wrong data @ hDet2Cluster."));
1350     return kTRUE;
1351   }
1352   Int_t ndim(H->GetNdimensions());
1353   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1354   TAxis *aa[kNdim], *an(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1355   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1356   if(ndim < 5) aa[4] = new TAxis(1, -0.5, 0.5);
1357   Int_t nPad(1);
1358   if(ndim > 5){
1359     an = H->GetAxis(5);
1360     nPad = kNpads+1;
1361   }
1362   // build list of projections
1363   const Int_t nsel(8*AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1364   // define rebinning strategy
1365   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1366   //const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1367   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1368   AliTRDresolutionProjection hp[kDetNproj];  TObjArray php(kDetNproj);
1369   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1370   for(Int_t ipad(0); ipad<nPad; ipad++){
1371     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1372       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1373         isel++; // new selection
1374         hp[ih].Build(Form("HDet%d%d%d", ily, icen, ipad),
1375                     Form("Detectors :: Ly[%d] Cen[%s] Pads[%s]", ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1376                     kEta, kPhi, 4, aa);
1377         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1378         hp[ih].SetShowRange(10., 55.);
1379         php.AddLast(&hp[ih++]); np[isel]++;
1380       }
1381     }
1382   }
1383   AliInfo(Form("Build %3d 3D projections.", ih));
1384
1385   Int_t ly(0), cen(0), npad(0);
1386   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1387     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1388     ly = coord[kBC]-1;    // layer selection
1389     cen = coord[kYrez]-1; // centrality selection
1390     npad = 0;             // no. of pads selection
1391     if(an) npad = TMath::Min(coord[5]-1, Int_t(kNpads));
1392     isel = npad*AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer+cen*AliTRDgeometry::kNlayer+ly;
1393     ((AliTRDresolutionProjection*)php.At(isel))->Increment(coord, v);
1394     //Int_t ioff=isel;for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1395   }
1396   TObjArray *arr(NULL);
1397   fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1398
1399   TH2 *h2(NULL);  Int_t jh(0);
1400   for(; ih--; ){
1401     if(!hp[ih].fH) continue;
1402     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1403     arr->AddAt(h2, jh++);
1404     if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].fH->GetName())))) continue;
1405     arr->AddAt(h2, jh++);
1406   }
1407   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1408   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1409     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1410       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, 0)))){
1411         for(Int_t ipad(1); ipad<nPad; ipad++){
1412           if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, ipad)))){
1413             (*pr0)+=(*pr1);
1414           }
1415         }
1416         pr0->fH->SetNameTitle(Form("HDet%d%d", ily, icen), Form("Detectors :: Ly[%d] Cen[%s]", ily, cenName[icen]));
1417         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1418         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
1419         if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))) (*pr1)+=(*pr0);
1420       }
1421     }
1422     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))){
1423       pr0->fH->SetNameTitle(Form("HDet%d", ily), Form("Detectors :: Ly[%d]", ily));
1424       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1425     }
1426   }
1427   AliInfo(Form("Done %3d 2D projections.", jh));
1428   return kTRUE;
1429 }
1430
1431 //________________________________________________________
1432 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1433 {
1434 // Analyse cluster
1435   const Int_t kNcontours(9);
1436   const Int_t kNstat(100);
1437   Int_t cidx=mc?kMCcluster:kCluster;
1438   if(fProj && fProj->At(cidx)) return kTRUE;
1439   if(!fContainer){
1440     AliError("Missing data container.");
1441     return kFALSE;
1442   }
1443   const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1444   THnSparse *H(NULL);
1445   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1446     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1447     return kFALSE;
1448   }
1449   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimCl));
1450   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1451   TAxis *aa[kNdim], *as(NULL), *apt(NULL), *acen(NULL), *anp(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1452   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1453   if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
1454   if(ndim > Int_t(kSpeciesChgRC)) as  = H->GetAxis(kSpeciesChgRC);
1455   if(ndim > Int_t(kNdim)) acen  = H->GetAxis(kNdim);
1456   if(ndim > Int_t(kNdim)+1) anp  = H->GetAxis(kNdim+1);
1457   // calculate size depending on debug level
1458   const Int_t nCh(apt?2:1);
1459   const Int_t nCen(acen?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1460   const Int_t nNpad(anp?(Int_t(kNpads)+1):1);
1461
1462   // build list of projections
1463   const Int_t nsel(AliTRDgeometry::kNlayer*kNcharge*(kNpads+1)*AliTRDeventInfo::kCentralityClasses);
1464   // define rebinning strategy
1465   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1466   AliTRDresolutionProjection hp[kClNproj];  TObjArray php(kClNproj);
1467   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1468   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1469   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1470   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1471     for(Int_t ich(0); ich<nCh; ich++){
1472       for(Int_t icen(0); icen<nCen; icen++){
1473         for(Int_t ipad(0); ipad<nNpad; ipad++){
1474           isel++; // new selection
1475           hp[ih].Build(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1476                       Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1477                       kEta, kPhi, kYrez, aa);
1478           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1479           php.AddLast(&hp[ih++]); np[isel]++;
1480           if(!debug) break;
1481           hp[ih].Build(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1482                       Form("Clusters[%c] :: Q Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1483                       kEta, kPhi, kSpeciesChgRC, aa);
1484           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1485           hp[ih].SetShowRange(24., 33.);
1486           php.AddLast(&hp[ih++]); np[isel]++;
1487           hp[ih].Build(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1488                       Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1489                       kPrez, kZrez, kYrez, aa);
1490           php.AddLast(&hp[ih++]); np[isel]++;
1491           hp[ih].Build(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1492                         Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s] Pads[%s]", chSgn[ich], ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1493                         kPrez, kPt, kYrez, aa);
1494           php.AddLast(&hp[ih++]); np[isel]++;
1495         }
1496       }
1497     }
1498   }
1499   AliInfo(Form("Build %3d 3D projections.", ih));
1500
1501   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1502   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1), ioff(0), cen(0), npad(0);
1503   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1504     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1505     ly = coord[kBC]-1;
1506     // RC selection
1507     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1508
1509     // charge selection
1510     ch = 0; // [-] track
1511     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1512     cen = 0; // centrality selection
1513     if(acen) cen = coord[kNdim]-1;
1514     npad = 0;             // no. of pads selection
1515     if(anp) npad = TMath::Min(coord[kNdim+1]-1, Int_t(kNpads));
1516
1517     if(debug){
1518       isel = ly*nCh*nCen*nNpad
1519             +ch*nCen*nNpad
1520             +cen*nNpad
1521             +npad;
1522       ioff=isel*4;
1523     } else {
1524       isel = ly; ioff = isel;
1525     }
1526     if(ioff>=ih){
1527       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1528       return kFALSE;
1529     }
1530     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
1531       AliError(Form("Missing projection %d", ioff));
1532       return kFALSE;
1533     }
1534     if(strcmp(pr0->fH->GetName(), Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ch], ly, cen, npad))!=0){
1535       AliError(Form("Projection mismatch :: request[H%sClY%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ly, cen, npad, pr0->fH->GetName()));
1536       return kFALSE;
1537     }
1538     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1539   }
1540   TObjArray *arr(NULL);
1541   fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1542
1543   TH2 *h2(NULL);  Int_t jh(0);
1544   for(; ih--; ){
1545     if(!hp[ih].fH) continue;
1546     if(strchr(hp[ih].fH->GetName(), 'Q')){
1547       if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1548       arr->AddAt(h2, jh++);
1549       if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].fH->GetName())))) continue;
1550       arr->AddAt(h2, jh++);
1551     } else {
1552       if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1553     }
1554   }
1555   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1556     for(Int_t ich(0); ich<nCh; ich++){
1557       for(Int_t icen(0); icen<nCen; icen++){
1558         /*!dy*/
1559         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1560           for(Int_t ipad(1); ipad<nNpad; ipad++){
1561             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1562             (*pr0)+=(*pr1);
1563           }
1564           pr0->fH->SetNameTitle(Form("H%sClY%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1565           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1566           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1567         }
1568         /*!Q*/
1569         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1570           for(Int_t ipad(1); ipad<nNpad; ipad++){
1571             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1572             (*pr0)+=(*pr1);
1573           }
1574           pr0->fH->SetNameTitle(Form("H%sClQ%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: Q Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1575           if((h2 = pr0->Projection2D(kNstat, kNcontours, 2, kFALSE))) arr->AddAt(h2, jh++);
1576           if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
1577           pr0->fH->SetName(Form("H%sClQS%c%d%d", mc?"MC":"", chName[ich], ily, icen));
1578           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1579           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1580         }
1581         /*!YXTC*/
1582         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1583           for(Int_t ipad(1); ipad<nNpad; ipad++){
1584             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1585             (*pr0)+=(*pr1);
1586           }
1587           pr0->fH->SetNameTitle(Form("H%sClYXTC%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,TC) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1588           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1589           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1590         }
1591         /*!YXPh*/
1592         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1593           for(Int_t ipad(1); ipad<nNpad; ipad++){
1594             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1595             (*pr0)+=(*pr1);
1596           }
1597           pr0->fH->SetNameTitle(Form("H%sClYXPh%c%d%d", mc?"MC":"", chName[ich], ily, icen), Form("Clusters[%c] :: #Deltay(x,#Phi) Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1598           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1599           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1600         }
1601       } // end centrality integration
1602       /*!dy*/
1603       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1604         pr0->fH->SetNameTitle(Form("H%sClY%c%d", mc?"MC":"", chName[ich], ily),
1605                               Form("Clusters[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
1606         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1607         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1608       }
1609       /*!Q*/
1610       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1611         pr0->fH->SetNameTitle(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], ily),
1612                               Form("Clusters[%c]:: Q Ly[%d]", chSgn[ich], ily));
1613         if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1614         pr0->fH->SetName(Form("H%sClQS%c%d", mc?"MC":"", chName[ich], ily));
1615         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1616         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1617       }
1618       /*!YXTC*/
1619       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1620         pr0->fH->SetNameTitle(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], ily),
1621                               Form("Clusters[%c]:: #Deltay(x,TC) Ly[%d]", chSgn[ich], ily));
1622         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1623         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1624       }
1625       /*!YXPh*/
1626       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1627         pr0->fH->SetNameTitle(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], ily),
1628                               Form("Clusters[%c]:: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], ily));
1629         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1630         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1631       }
1632     } // end charge integration
1633     /*!dy*/
1634     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1635       pr0->fH->SetNameTitle(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1636       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1637     }
1638     /*!YXPh*/
1639     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1640       pr0->fH->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1641       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1642     }
1643
1644   }
1645   AliInfo(Form("Done %3d 2D projections.", jh));
1646   return kTRUE;
1647 }
1648
1649 //________________________________________________________
1650 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1651 {
1652 // Analyse tracklet
1653   const Int_t kNcontours(9);
1654   const Int_t kNstat(30);
1655   const Int_t kNstatQ(30);
1656   Int_t cidx=mc?kMCtracklet:kTracklet;
1657   if(fProj && fProj->At(cidx)) return kTRUE;
1658   if(!fContainer){
1659     AliError("Missing data container.");
1660     return kFALSE;
1661   }
1662   const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1663   THnSparse *H(NULL);
1664   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1665     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1666     return kFALSE;
1667   }
1668   const Int_t mdim(kNdim+8);
1669   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimTrklt));
1670   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1671   TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1672   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1673   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1674   if(ndim > Int_t(kPt))           ap = H->GetAxis(kPt);           // init pt selection
1675   if(ndim > Int_t(kNdim))         ac = H->GetAxis(kNdim);         // init centrality selection
1676   // calculate size depending on debug level
1677   const Int_t nCen(debug?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1678   const Int_t nPt(debug?Int_t(kNpt):1);
1679   const Int_t nSpc(1);//ndim>kNdimTrklt?fgkNbins[kSpeciesChgRC]:1);
1680   const Int_t nCh(debug?Int_t(kNcharge):1);
1681
1682   // build list of projections
1683   const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*kNpt*(AliPID::kSPECIES*kNcharge + 1));
1684   // define rebinning strategy
1685   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1686   AliTRDresolutionProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1687   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1688   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1689   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
1690   const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
1691 //  const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1692   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1693   for(Int_t icen(0); icen<nCen; icen++){
1694     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1695       for(Int_t ipt(0); ipt<nPt; ipt++){
1696         for(Int_t isp(0); isp<nSpc; isp++){
1697           for(Int_t ich(0); ich<nCh; ich++){
1698             isel++; // new selection
1699             hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1700                         Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1701                         kEta, kPhi, kYrez, aa);
1702             //hp[ih].SetShowRange(-0.1,0.1);
1703             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1704             php.AddLast(&hp[ih++]); np[isel]++;
1705             hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1706                         Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1707                         kEta, kPhi, kPrez, aa);
1708             //hp[ih].SetShowRange(-0.5,0.5);
1709             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1710             php.AddLast(&hp[ih++]); np[isel]++;
1711             hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1712                         Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1713                         kEta, kPhi, kZrez, aa);
1714             hp[ih].SetShowRange(1.,2.3);
1715             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1716             php.AddLast(&hp[ih++]); np[isel]++;
1717             hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1718                         Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1719                         kEta, kPhi, kNdim+1, aa);
1720             hp[ih].SetShowRange(30., 70.);
1721             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1722             php.AddLast(&hp[ih++]); np[isel]++;
1723           }
1724         }
1725         if(ndim==kNdimTrklt) continue;
1726
1727         isel++; // new selection
1728         hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1729                     Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1730                     kEta, kPhi, kZrez, aa);
1731   //      hp[ih].SetShowRange(-0.1,0.1);
1732         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1733         php.AddLast(&hp[ih++]); np[isel]++;
1734 /*        hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1735                     Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1736                     kEta, kPhi, kYrez, aa);
1737         //hp[ih].SetShowRange(-0.1,0.1);
1738         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1739         php.AddLast(&hp[ih++]); np[isel]++;
1740         hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1741                     Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1742                     kEta, kPhi, kPrez, aa);
1743         //hp[ih].SetShowRange(-0.1,0.1);
1744         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1745         php.AddLast(&hp[ih++]); np[isel]++;
1746         hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1747                     Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1748                     kEta, kPhi, kNdim, aa);
1749         //hp[ih].SetShowRange(-0.1,0.1);
1750         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1751         php.AddLast(&hp[ih++]); np[isel]++;*/
1752       }
1753     }
1754   }
1755   AliInfo(Form("Build %3d 3D projections.", ih));
1756
1757   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1758   Int_t ly(0), ch(0), sp(2), rcBin(as?as->FindBin(0.):-1), pt(0), cen(0), ioff(0), jspc(nSpc*nCh+1), kspc(nSpc*nCh*4+1);
1759   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1760     v = H->GetBinContent(ib, coord);
1761     if(v<1.) continue;
1762     ly = coord[kBC]-1; // layer selection
1763     // charge selection
1764     ch = 0; sp=0;// [e-] track [dafault]
1765     if(rcBin>0){ // debug mode in which species/charge are also saved
1766       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1767       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1768       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1769     }
1770     // pt selection
1771     pt = 0; // low pt
1772     if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
1773     // centrality selection
1774     cen = 0; // default
1775     if(ac) cen = coord[kNdim]-1;
1776     // global selection
1777     if(ndim==kNdimTrklt){
1778       ioff = ly*4;
1779       isel = ly;
1780     } else {
1781       isel = cen*AliTRDgeometry::kNlayer*nPt*jspc+ly*nPt*jspc+pt*jspc; isel+=sp<0?(nSpc*nCh):ch;
1782       ioff = cen*AliTRDgeometry::kNlayer*nPt*kspc+ly*nPt*kspc+pt*kspc; ioff+=sp<0?((nSpc*nCh)*4):(ch*4);
1783     }
1784     if(ioff>=ih){
1785       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1786       return kFALSE;
1787     }
1788     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
1789       AliError(Form("Missing projection %d", ioff));
1790       return kFALSE;
1791     }
1792     if(sp>=0){
1793       if(strcmp(pr0->fH->GetName(), Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen))!=0){
1794         AliError(Form("Projection mismatch :: request[H%sTrkltY%c%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen, pr0->fH->GetName()));
1795         return kFALSE;
1796       }
1797     } else {
1798       if(strcmp(pr0->fH->GetName(), Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[pt], ly, cen))!=0){
1799         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c%d%d] found[%s]", mc?"MC":"", ptName[pt], ly, cen, pr0->fH->GetName()));
1800         return kFALSE;
1801       }
1802     }
1803     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1804   }
1805   TObjArray *arr(NULL);
1806   fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1807
1808   TH2 *h2(NULL); Int_t jh(0);
1809   for(; ih--; ){
1810     if(!hp[ih].fH) continue;
1811     Int_t mid(0), nstat(kNstat);
1812     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1813     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1814     arr->AddAt(h2, jh++);
1815   }
1816   // build combined performance plots
1817   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1818     for(Int_t ich(0); ich<nCh; ich++){
1819       for(Int_t icen(0); icen<nCen; icen++){
1820         for(Int_t ipt(0); ipt<nPt; ipt++){
1821           /*!dy*/
1822           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1823             for(Int_t isp(1); isp<nSpc; isp++){
1824               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1825               (*pr0)+=(*pr1);
1826             }
1827             pr0->fH->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1828                                       Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1829             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1830             if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1831           }
1832           /*!dphi*/
1833           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1834             for(Int_t isp(1); isp<nSpc; isp++){
1835               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1836               (*pr0)+=(*pr1);
1837             }
1838             pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1839                                       Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1840             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1841             if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1842           }
1843           /*!dQ/dl*/
1844           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1845             for(Int_t isp(1); isp<nSpc; isp++){
1846               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1847               (*pr0)+=(*pr1);
1848             }
1849             pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1850                                       Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1851             if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1852             pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1853                                       Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1854             pr0->SetShowRange(2.4, 5.1);
1855             if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1856             if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1857           }
1858           /*!TB occupancy*/
1859           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1860             for(Int_t isp(1); isp<nSpc; isp++){
1861               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1862               (*pr0)+=(*pr1);
1863             }
1864             pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1865                                       Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1866             if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1867             if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))) (*pr1)+=(*pr0);
1868           }
1869         } // end pt integration
1870         /*!dy*/
1871         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1872           pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1873                                 Form("Tracklets[%c]:: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1874           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1875           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1876         }
1877         /*!dphi*/
1878         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1879           pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1880                                 Form("Tracklets[%c]:: #Delta#phi Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1881           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1882           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1883         }
1884         /*!dQ/dl*/
1885         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1886           pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1887                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1888           pr0->SetShowRange(1.,2.3);
1889           if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1890           pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1891                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1892           pr0->SetShowRange(2.4,5.1);
1893           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1894           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1895         }
1896         /*!TB occupancy*/
1897         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1898           pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1899                                 Form("Tracklets[%c]:: OccupancyTB Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1900           if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1901           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1902         }
1903       } // end centrality integration
1904       /*!dy*/
1905       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1906         pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], ily));
1907         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1908         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1909       }
1910       /*!dphi*/
1911       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1912         pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], ily));
1913         pr0->SetShowRange(-.9,.9);
1914         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1915         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1916       }
1917       /*!dQ/dl*/
1918       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1919         pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1920         pr0->SetShowRange(1.,2.3);
1921         if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1922         pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1923         pr0->SetShowRange(2.4,5.1);
1924         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1925         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1926       }
1927       /*!TB occupancy*/
1928       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1929         pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], ily));
1930         if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1931         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))) (*pr1)+=(*pr0);
1932       }
1933     } // end charge integration
1934     /*!dy*/
1935     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1936       pr0->fH->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: #Deltay Ly[%d]", ily));
1937       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1938     }
1939     /*!dphi*/
1940     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1941       pr0->fH->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", ily), Form("Tracklets :: #Delta#phi Ly[%d]", ily));
1942       pr0->SetShowRange(-.45,.45);
1943       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1944     }
1945     /*!dQdl*/
1946     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1947       pr0->fH->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
1948       pr0->SetShowRange(1.,2.3);
1949       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1950       pr0->fH->SetName(Form("H%sTrkltQS%d", mc?"MC":"", ily));
1951       pr0->SetShowRange(2.4,5.1);
1952       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1953     }
1954     /*!TB occupancy*/
1955     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1956       pr0->fH->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", ily), Form("Tracklets :: OccupancyTB Ly[%d]", ily));
1957       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1958     }
1959
1960     /*! Row Cross processing*/
1961     for(Int_t icen(0); icen<nCen; icen++){
1962       /*!RC dz*/
1963       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, icen)))){
1964         for(Int_t ipt(0); ipt<kNpt; ipt++){
1965           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen)))) continue;
1966           (*pr0)+=(*pr1);
1967         }
1968         pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", ily, cenName[icen]));
1969         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1970         if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))) (*pr1)+=(*pr0);
1971       }
1972     } // end centrality integration for row cross
1973     /*!RC dz*/
1974     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))){
1975       pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", ily));
1976       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1977     }
1978   } // end layer loop
1979   AliInfo(Form("Done %3d 2D projections.", jh));
1980   return kTRUE;
1981 }
1982
1983 //________________________________________________________
1984 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1985 {
1986 // Analyse track in
1987   const Int_t kNcontours(9);
1988   const Int_t kNstat(30);
1989   Int_t cidx=mc?kMCtrackIn:kTrackIn;
1990   if(fProj && fProj->At(cidx)) return kTRUE;
1991   if(!fContainer){
1992     AliError("Missing data container.");
1993     return kFALSE;
1994   }
1995   const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC"};
1996   THnSparse *H(NULL);
1997   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1998     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1999     return kFALSE;
2000   }
2001
2002   const Int_t mdim(kNdim+3);
2003   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
2004   Int_t ndim(H->GetNdimensions());
2005   TAxis *aa[mdim], *as(NULL), *ap(NULL), *ax(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
2006   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2007   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
2008   if(ndim > Int_t(kPt))           ap = H->GetAxis(kPt);
2009   if(ndim > Int_t(kNdim)+1)       ax = H->GetAxis(kNdim+1);
2010   if(ndim > Int_t(kNdim)+2)      abf = H->GetAxis(kNdim+2);
2011   //AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
2012   const Int_t nPt(ndim>Int_t(kNdimTrkIn)?Int_t(kNpt):1);
2013
2014   // build list of projections
2015   const Int_t nsel(kNpt*(AliPID::kSPECIES*kNcharge + 1));
2016   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2017   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2018   const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
2019   // define rebinning strategy
2020   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
2021   AliTRDresolutionProjection hp[kMCTrkInNproj]; TObjArray php(kMCTrkInNproj+2);
2022   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2023   // define list of projections
2024   for(Int_t ipt(0); ipt<nPt; ipt++){
2025     for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2026       for(Int_t ich(0); ich<kNcharge; ich++){
2027         isel++; // new selection
2028         hp[ih].Build(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2029                      Form("TrackIn[%s%c]:: #Deltay{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2030                      kEta, kPhi, kYrez, aa);
2031         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2032         php.AddLast(&hp[ih++]); np[isel]++;
2033         hp[ih].Build(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2034                      Form("TrackIn[%s%c]:: #Delta#phi{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2035                      kEta, kPhi, kPrez, aa);
2036         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2037         php.AddLast(&hp[ih++]); np[isel]++;
2038         hp[ih].Build(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2039                      Form("TrackIn[%s%c]:: dQdl {%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2040                      kEta, kPhi, kZrez, aa);
2041         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2042         php.AddLast(&hp[ih++]); np[isel]++;
2043         if(!ax) continue;
2044         hp[ih].Build(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2045                      Form("TrackIn[%s%c]:: #Deltax{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2046                      kEta, kPhi, kNdim+1, aa);
2047         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2048         php.AddLast(&hp[ih++]); np[isel]++;
2049       }
2050     }
2051     isel++; // RC tracks
2052     hp[ih].Build(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt]),
2053                   Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
2054                   kEta, kPhi, kZrez, aa);
2055     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2056     php.AddLast(&hp[ih++]); np[isel]++;
2057     hp[ih].Build(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt]),
2058                   Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
2059                   kEta, kPhi, kYrez, aa);
2060     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2061     php.AddLast(&hp[ih++]); np[isel]++;
2062     hp[ih].Build(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt]),
2063                   Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
2064                   kEta, kPhi, kPrez, aa);
2065     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2066     php.AddLast(&hp[ih++]); np[isel]++;
2067     if(!ax) continue;
2068     hp[ih].Build(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt]),
2069                   Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
2070                   kEta, kPhi, kNdim+1, aa);
2071     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2072     php.AddLast(&hp[ih++]); np[isel]++;
2073   }
2074   AliInfo(Form("Build %3d 3D projections.", ih));
2075
2076   // fill projections
2077   Int_t ch(0), pt(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0);
2078   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
2079   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2080     v = H->GetBinContent(ib, coord);
2081     if(v<1.) continue;
2082     if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
2083     if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
2084     // charge selection
2085     ch = 0; sp=1;// [-] track
2086     if(rcBin>0){ // debug mode in which species are also saved
2087       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2088       if(sp>=AliPID::kSPECIES){
2089         AliDebug(2, Form("Wrong SpeciesIndex[%d]. Rescale", sp));
2090         sp = AliPID::kSPECIES-1;
2091       }
2092       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
2093       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2094     }
2095     // pt selection
2096     pt = 0; // low pt
2097     if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
2098     // global selection
2099     isel = pt*(AliPID::kSPECIES*kNcharge+1); isel+=sp<0?(AliPID::kSPECIES*kNcharge):(sp*kNcharge+ch);
2100     ioff = isel*(ax?4:3);
2101     if(ioff>=ih){
2102       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
2103       return kFALSE;
2104     }
2105     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
2106       AliError(Form("Missing projection %d", ioff));
2107       return kFALSE;
2108     }
2109     if(sp>=0){
2110       if(strcmp(pr0->fH->GetName(), Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ch], ptName[pt], sp))!=0){
2111         AliError(Form("Projection mismatch :: request[H%sTrkInY%c%c%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, pr0->fH->GetName()));
2112         return kFALSE;
2113       }
2114     } else {
2115       if(strcmp(pr0->fH->GetName(), Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[pt]))!=0){
2116         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c] found[%s]", mc?"MC":"", ptName[pt], pr0->fH->GetName()));
2117         return kFALSE;
2118       }
2119     }
2120     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
2121   }
2122   TObjArray *arr(NULL);
2123   fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
2124
2125   TH2 *h2(NULL); Int_t jh(0);
2126   for(; ih--; ){
2127     if(!hp[ih].fH) continue;
2128     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2129     arr->AddAt(h2, jh++);
2130   }
2131   // build combined performance plots
2132   // combine up the tree of projections
2133   AliTRDresolutionProjection xlow[2], specY[kNcharge*AliPID::kSPECIES], specPh[kNcharge*AliPID::kSPECIES], specQ[kNcharge*AliPID::kSPECIES];
2134   for(Int_t ich(0); ich<kNcharge; ich++){
2135     // PID dependency - summation over pt
2136     for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2137       /*!dy*/
2138       Int_t idx(ich*AliPID::kSPECIES+isp);
2139       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2140         specY[idx] = (*pr0);
2141         specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2142         specY[idx].fH->SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp),
2143                               Form("TrackIn[%s%c]:: #Deltay", AliPID::ParticleLatexName(isp), chSgn[ich]));
2144         for(Int_t ipt(1); ipt<nPt; ipt++){
2145           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2146           specY[idx]+=(*pr1);
2147         }
2148         php.AddLast(&specY[idx]);
2149         if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2150         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", specY[idx].fH->GetName())))) arr->AddAt(h2, jh++);
2151         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specY[idx];
2152       }
2153       /*!dphi*/
2154       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2155         specPh[idx] = (*pr0);
2156         specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2157         specPh[idx].fH->SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp),
2158                               Form("TrackIn[%s%c]:: #Delta#phi", AliPID::ParticleLatexName(isp), chSgn[ich]));
2159         specPh[idx].SetShowRange(-1.5, 1.5);
2160         for(Int_t ipt(1); ipt<nPt; ipt++){
2161           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2162           specPh[idx]+=(*pr1);
2163         }
2164         php.AddLast(&specPh[idx]);
2165         if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2166         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specPh[idx];
2167       }
2168       /*!dQdl*/
2169       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2170         specQ[idx] = (*pr0);
2171         specQ[idx].SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2172         specQ[idx].fH->SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp),
2173                               Form("TrackIn[%s%c]:: dQdl", AliPID::ParticleLatexName(isp), chSgn[ich]));
2174         specQ[idx].SetShowRange(-2.2, -1.75);
2175         specQ[idx].fH->GetZaxis()->SetTitle("dQdl [a.u.]");
2176         for(Int_t ipt(1); ipt<nPt; ipt++){
2177           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2178           specQ[idx]+=(*pr1);
2179         }
2180         php.AddLast(&specQ[idx]);
2181         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2182         specQ[idx].fH->SetName(Form("H%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp));
2183         specQ[idx].SetShowRange(-1.75, -1.25);
2184         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2185         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specQ[idx];
2186       }
2187     } // end PID loop for pt integration
2188
2189     // pt dependency - summation over PID
2190     for(Int_t ipt(0); ipt<nPt; ipt++){
2191       /*!dy*/
2192       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2193         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2194           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2195           (*pr0)+=(*pr1);
2196         }
2197         pr0->fH->SetNameTitle(Form("H%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2198                                   Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
2199         pr0->SetShowRange(-0.3, 0.3);
2200         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2201         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2202       }
2203       /*!dphi*/
2204       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2205         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2206           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2207           (*pr0)+=(*pr1);
2208         }
2209         pr0->fH->SetNameTitle(Form("H%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2210                                   Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
2211         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2212         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2213       }
2214       /*!dx*/
2215       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2216         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2217           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2218           (*pr0)+=(*pr1);
2219         }
2220         pr0->fH->SetNameTitle(Form("H%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2221                                   Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
2222         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2223         if(!ipt){
2224           xlow[ich] = (*pr0);
2225           xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5),
2226                                  Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2227           php.AddLast(&xlow[ich]);
2228         }
2229         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2230       }
2231     } // end pt loop for PID integration
2232
2233     /*!dy*/
2234     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2235       pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", chName[ich]),
2236                             Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2237       pr0->SetShowRange(-0.3, 0.3);
2238       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2239       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2240       if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2241     }
2242     /*!dy high pt*/
2243     if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[2], 0)))){
2244       if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[2], 0)))){
2245         (*pr0)+=(*pr1);
2246         pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", ptName[2]), Form("TrackIn :: #Deltay{%s}", ptCut[2]));
2247         pr0->SetShowRange(-0.3, 0.3);
2248         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2249       }
2250     }
2251     /*!dphi*/
2252     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2253       pr0->fH->SetNameTitle(Form("H%sTrkInPh%c", mc?"MC":"", chName[ich]),
2254                             Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2255       pr0->SetShowRange(-1., 1.);
2256       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2257       if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2258     }
2259     /*!dx*/
2260     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2261       pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", chName[ich]),
2262                             Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2263       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2264       if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2265     }
2266     /*!dx low pt*/
2267     if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 5)))){
2268       if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5)))){
2269         (*pr0)+=(*pr1);
2270         pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", ptName[0]), Form("TrackIn :: #Deltax{%s}", ptCut[0]));
2271         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2272       }
2273     }
2274   } // end charge loop
2275
2276   for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2277     /*!dy*/
2278     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))){
2279       pr0->fH->SetNameTitle(Form("H%sTrkInY%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Deltay", AliPID::ParticleLatexName(isp)));
2280       pr0->SetShowRange(-0.3, 0.3);
2281       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2282       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2283     }
2284     /*!dphi*/
2285     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))){
2286       pr0->fH->SetNameTitle(Form("H%sTrkInPh%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Delta#phi", AliPID::ParticleLatexName(isp)));
2287       pr0->SetShowRange(-1., 1.);
2288       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2289     }
2290     /*!dQdl*/
2291     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))){
2292       pr0->fH->SetNameTitle(Form("H%sTrkInQ%d", mc?"MC":"", isp), Form("TrackIn[%s] :: dQdl", AliPID::ParticleLatexName(isp)));
2293       pr0->SetShowRange(-2.2, -1.75);
2294       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2295       pr0->fH->SetName(Form("H%sTrkInQS%d", mc?"MC":"", isp));
2296       pr0->SetShowRange(-1.7, -1.25);
2297       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2298     }
2299   } // end PID processing
2300
2301   /*!dy*/
2302   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2303     pr0->fH->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
2304     pr0->SetShowRange(-0.3, 0.3);
2305     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2306     if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2307   }
2308   /*!dphi*/
2309   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2310     pr0->fH->SetNameTitle(Form("H%sTrkInPh", mc?"MC":""), "TrackIn :: #Delta#phi");
2311     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2312   }
2313   /*!dx*/
2314   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2315     pr0->fH->SetNameTitle(Form("H%sTrkInX", mc?"MC":""), "TrackIn :: #Deltax");
2316     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2317   }
2318
2319   // Row Cross processing
2320   /*!RC dz*/
2321   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[0])))){
2322     for(Int_t ipt(0); ipt<kNpt; ipt++){
2323       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt])))) continue;
2324       (*pr0)+=(*pr1);
2325     }
2326     pr0->fH->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
2327     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2328   }
2329   /*!RC dy*/
2330   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[0])))){
2331     for(Int_t ipt(0); ipt<kNpt; ipt++){
2332       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt])))) continue;
2333       (*pr0)+=(*pr1);
2334     }
2335     pr0->fH->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay");
2336     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2337   }
2338   /*!RC dphi*/
2339   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[0])))){
2340     for(Int_t ipt(0); ipt<kNpt; ipt++){
2341       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt])))) continue;
2342       (*pr0)+=(*pr1);
2343     }
2344     pr0->fH->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi");
2345     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2346   }
2347   /*!RC dx*/
2348   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[0])))){
2349     for(Int_t ipt(0); ipt<kNpt; ipt++){
2350       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt])))) continue;
2351       (*pr0)+=(*pr1);
2352     }
2353     pr0->fH->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax");
2354     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2355   }
2356   AliInfo(Form("Done %3d 2D projections.", jh));
2357   return kTRUE;
2358 }
2359
2360
2361 //________________________________________________________
2362 Bool_t AliTRDresolution::MakeProjectionTrack()
2363 {
2364 // Analyse tracklet
2365   const Int_t kNcontours(9);
2366   const Int_t kNstat(30);
2367   Int_t cidx(kMCtrack);
2368   if(fProj && fProj->At(cidx)) return kTRUE;
2369   if(!fContainer){
2370     AliError("Missing data container.");
2371     return kFALSE;
2372   }
2373   THnSparse *H(NULL);
2374   if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2375     AliError("Missing/Wrong data @ hTRD2MC.");
2376     return kFALSE;
2377   }
2378   Int_t ndim(H->GetNdimensions());
2379   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2380   TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2381   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2382   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2383   if(ndim > kPt) ap = H->GetAxis(kPt);
2384
2385   // build list of projections
2386   const Int_t nsel(AliTRDgeometry::kNlayer*kNpt*AliPID::kSPECIES*7);//, npsel(3);
2387   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2388   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2389   const Char_t *ptCut[kNpt] = {"p_{t}[GeV/c]<0.8", "0.8<=p_{t}[GeV/c]<1.5", "p_{t}[GeV/c]>=1.5"};
2390   // define rebinning strategy
2391   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
2392   AliTRDresolutionProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2393   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2394   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2395     for(Int_t ipt(0); ipt<kNpt; ipt++){
2396       for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2397         for(Int_t ich(0); ich<kNcharge; ich++){
2398           isel++; // new selection
2399           hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2400                        Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2401                        kEta, kPhi, kYrez, aa);
2402           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2403           php.AddLast(&hp[ih++]); np[isel]++;
2404           hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2405                        Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2406                        kEta, kPhi, kPrez, aa);
2407           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2408           php.AddLast(&hp[ih++]); np[isel]++;
2409           hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2410                        Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2411                        kEta, kPhi, kNdim, aa);
2412           hp[ih].SetShowRange(0.,10.);
2413           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2414           php.AddLast(&hp[ih++]); np[isel]++;
2415         }
2416       }
2417       isel++; // new selection
2418       hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], ily),
2419                     Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily),
2420                     kEta, kPhi, kZrez, aa);
2421       hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2422       php.AddLast(&hp[ih++]); np[isel]++;
2423     }
2424   }
2425
2426   Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2427   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2428     v = H->GetBinContent(ib, coord);
2429     if(v<1.) continue;
2430     ly = coord[kBC]-1; // layer selection
2431     // charge selection
2432     ch=0; sp=2;// [pi-] track [dafault]
2433     if(rcBin>0){ // debug mode in which species are also saved
2434       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2435       if(coord[kSpeciesChgRC] > rcBin) ch = 1;        // [+] track
2436       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2437     }
2438     // pt selection
2439     pt = 0; // low pt [default]
2440     if(ap) pt = coord[kPt]-1;
2441     // global selection
2442     Int_t ioff = ly*kNpt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2443     isel = ly*kNpt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2444     AliDebug(4, Form("SELECTION[%d] :: ch[%c] pt[%c] sp[%d] ly[%d]\n", np[isel], ch==2?'Z':chName[ch], ptName[pt], sp, ly));
2445     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
2446   }
2447   TObjArray *arr(NULL);
2448   fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2449
2450   TH2 *h2(NULL); Int_t jh(0);
2451   for(; ih--; ){
2452     if(!hp[ih].fH) continue;
2453     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2454     arr->AddAt(h2, jh++);
2455   }
2456
2457   // combine up the tree of projections
2458   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
2459   //Int_t iproj(0), jproj(0);
2460   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2461     for(Int_t ich(0); ich<kNcharge; ich++){
2462       for(Int_t ipt(0); ipt<kNpt; ipt++){
2463         /*!dy*/
2464         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2465           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2466             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2467             (*pr0)+=(*pr1);
2468           }
2469           AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2470           pr0->fH->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], ily),
2471                                     Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2472           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2473           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2474         }
2475         /*!dphi*/
2476         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2477           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2478             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2479             (*pr0)+=(*pr1);
2480           }
2481           AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2482           pr0->fH->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], ily),
2483                                     Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2484           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2485           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2486         }
2487
2488         /*!dpt/pt*/
2489         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2490           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2491             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2492             (*pr0)+=(*pr1);
2493           }
2494           AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2495           pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], ily),
2496                                     Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2497           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2498           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2499         }
2500       }
2501       /*!dy*/
2502       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2503         pr0->fH->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], ily),
2504                               Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
2505         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2506         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2507       }
2508       /*!dphi*/
2509       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2510         pr0->fH->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], ily),
2511                               Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], ily));
2512         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2513         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2514       }
2515       /*!dpt/pt*/
2516       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2517         pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], ily),
2518                               Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], ily));
2519         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2520         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2521       }
2522     }
2523     /*!dy*/
2524     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2525       pr0->fH->SetNameTitle(Form("HMCTrkY%d", ily), Form("Tracks :: #Deltay Ly[%d]", ily));
2526       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2527     }
2528     /*!dphi*/
2529     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2530       pr0->fH->SetNameTitle(Form("HMCTrkPh%d", ily), Form("Tracks :: #Delta#phi Ly[%d]", ily));
2531       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2532     }
2533     /*!dpt/pt*/
2534     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2535       pr0->fH->SetNameTitle(Form("HMCTrkDPt%d", ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", ily));
2536       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2537     }
2538   }
2539   AliInfo(Form("Done %3d 2D projections.", jh));
2540   return kTRUE;
2541 }
2542
2543 //________________________________________________________
2544 Bool_t AliTRDresolution::PostProcess()
2545 {
2546 // Fit, Project, Combine, Extract values from the containers filled during execution
2547
2548   if (!fContainer) {
2549     AliError("ERROR: list not available");
2550     return kFALSE;
2551   }
2552   if(!fProj){
2553     AliInfo("Building array of projections ...");
2554     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
2555   }
2556
2557   //PROCESS EXPERIMENTAL DISTRIBUTIONS
2558   // Clusters detector
2559   if(HasProcessDetector() && !MakeProjectionDetector()) return kFALSE;
2560   // Clusters residuals
2561   if(HasProcessCluster() && !MakeProjectionCluster()) return kFALSE;
2562   fNRefFigures = 3;
2563   // Tracklet residual/pulls
2564   if(HasProcessTrklt() && !MakeProjectionTracklet()) return kFALSE;
2565   fNRefFigures = 7;
2566   // TRDin residual/pulls
2567   if(HasProcessTrkIn() && !MakeProjectionTrackIn()) return kFALSE;
2568   fNRefFigures = 11;
2569
2570   if(!HasMCdata()) return kTRUE;
2571   //PROCESS MC RESIDUAL DISTRIBUTIONS
2572
2573   // CLUSTER Y RESOLUTION/PULLS
2574   if(!MakeProjectionCluster(kTRUE)) return kFALSE;
2575   fNRefFigures = 17;
2576
2577   // TRACKLET RESOLUTION/PULLS
2578   if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2579   fNRefFigures = 21;
2580
2581   // TRACK RESOLUTION/PULLS
2582   if(!MakeProjectionTrack()) return kFALSE;
2583   fNRefFigures+=16;
2584
2585   // TRACK TRDin RESOLUTION/PULLS
2586   if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2587   fNRefFigures+=8;
2588
2589   return kTRUE;
2590 }
2591
2592
2593 //________________________________________________________
2594 void AliTRDresolution::Terminate(Option_t *opt)
2595 {
2596   AliTRDrecoTask::Terminate(opt);
2597   if(HasPostProcess()) PostProcess();
2598 }
2599
2600 //________________________________________________________
2601 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2602 {
2603 // Helper function to avoid duplication of code
2604 // Make first guesses on the fit parameters
2605
2606   // find the intial parameters of the fit !! (thanks George)
2607   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2608   Double_t sum = 0.;
2609   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2610   f->SetParLimits(0, 0., 3.*sum);
2611   f->SetParameter(0, .9*sum);
2612   Double_t rms = h->GetRMS();
2613   f->SetParLimits(1, -rms, rms);
2614   f->SetParameter(1, h->GetMean());
2615
2616   f->SetParLimits(2, 0., 2.*rms);
2617   f->SetParameter(2, rms);
2618   if(f->GetNpar() <= 4) return;
2619
2620   f->SetParLimits(3, 0., sum);
2621   f->SetParameter(3, .1*sum);
2622
2623   f->SetParLimits(4, -.3, .3);
2624   f->SetParameter(4, 0.);
2625
2626   f->SetParLimits(5, 0., 1.e2);
2627   f->SetParameter(5, 2.e-1);
2628 }
2629
2630 //________________________________________________________
2631 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
2632 {
2633 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2634 //  - y reziduals/pulls
2635
2636   TObjArray *arr = new TObjArray(2);
2637   arr->SetName(name); arr->SetOwner();
2638   TH1 *h(NULL); char hname[100], htitle[300];
2639
2640   // tracklet resolution/pull in y direction
2641   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
2642   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
2643   Float_t rr = range<0.?fDyRange:range;
2644   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2645     Int_t nybins=50;
2646     if(expand) nybins*=2;
2647     h = new TH3S(hname, htitle,
2648                  48, -.48, .48,            // phi
2649                  60, -rr, rr,              // dy
2650                  nybins, -0.5, nybins-0.5);// segment
2651   } else h->Reset();
2652   arr->AddAt(h, 0);
2653   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
2654   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
2655   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2656     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2657   } else h->Reset();
2658   arr->AddAt(h, 1);
2659
2660   return arr;
2661 }
2662
2663 //________________________________________________________
2664 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2665 {
2666 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2667 //  - y reziduals/pulls
2668 //  - z reziduals/pulls
2669 //  - phi reziduals
2670   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
2671   arr->Expand(5);
2672   TH1 *h(NULL); char hname[100], htitle[300];
2673
2674   // tracklet resolution/pull in z direction
2675   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
2676   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
2677   if(!(h = (TH2S*)gROOT->FindObject(hname))){
2678     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
2679   } else h->Reset();
2680   arr->AddAt(h, 2);
2681   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
2682   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
2683   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2684     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2685     h->GetZaxis()->SetBinLabel(1, "no RC");
2686     h->GetZaxis()->SetBinLabel(2, "RC");
2687   } else h->Reset();
2688   arr->AddAt(h, 3);
2689
2690   // tracklet to track phi resolution
2691   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
2692   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
2693   Int_t nsgms=540;
2694   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2695     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
2696   } else h->Reset();
2697   arr->AddAt(h, 4);
2698
2699   return arr;
2700 }
2701
2702 //________________________________________________________
2703 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2704 {
2705 // Build performance histograms for AliExternalTrackParam.vs MC
2706 //  - y resolution/pulls
2707 //  - z resolution/pulls
2708 //  - phi resolution, snp pulls
2709 //  - theta resolution, tgl pulls
2710 //  - pt resolution, 1/pt pulls, p resolution
2711
2712   TObjArray *arr = BuildMonitorContainerTracklet(name);
2713   arr->Expand(11);
2714   TH1 *h(NULL); char hname[100], htitle[300];
2715   //TAxis *ax(NULL);
2716
2717   // snp pulls
2718   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
2719   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
2720   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2721     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2722   } else h->Reset();
2723   arr->AddAt(h, 5);
2724
2725   // theta resolution
2726   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
2727   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2728   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2729     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2730   } else h->Reset();
2731   arr->AddAt(h, 6);
2732   // tgl pulls
2733   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
2734   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
2735   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2736     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2737   } else h->Reset();
2738   arr->AddAt(h, 7);
2739
2740   const Int_t kNdpt(150);
2741   const Int_t kNspc = 2*AliPID::kSPECIES+1;
2742   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
2743   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2744   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2745   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
2746   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
2747
2748   // Pt resolution
2749   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
2750   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2751   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2752     h = new TH3S(hname, htitle,
2753                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2754     //ax = h->GetZaxis();
2755     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2756   } else h->Reset();
2757   arr->AddAt(h, 8);
2758   // 1/Pt pulls
2759   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
2760   snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
2761   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2762     h = new TH3S(hname, htitle,
2763                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2764     //ax = h->GetZaxis();
2765     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2766   } else h->Reset();
2767   arr->AddAt(h, 9);
2768   // P resolution
2769   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
2770   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2771   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2772     h = new TH3S(hname, htitle,
2773                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2774     //ax = h->GetZaxis();
2775     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2776   } else h->Reset();
2777   arr->AddAt(h, 10);
2778
2779   return arr;
2780 }
2781
2782
2783 //________________________________________________________
2784 TObjArray* AliTRDresolution::Histos()
2785 {
2786   //
2787   // Define histograms
2788   //
2789
2790   if(fContainer) return fContainer;
2791
2792   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
2793   THnSparse *H(NULL);
2794   const Int_t nhn(100); Char_t hn[nhn]; TString st;
2795
2796   //++++++++++++++++++++++
2797   // cluster to detector
2798   snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
2799   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2800     const Int_t mdim(6);
2801     const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]", "no. of pads"/*, "tb [10 Hz]"*/};
2802     const Int_t cldNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 50, kNpads};
2803     const Double_t cldMin[mdim]  = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0., 0.5},
2804                    cldMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200., kNpads+.5};
2805     st = "cluster proprieties;";
2806     // define minimum info to be saved in non debug mode
2807     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
2808     for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
2809     H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
2810   } else H->Reset();
2811   fContainer->AddAt(H, kDetector);
2812   //++++++++++++++++++++++
2813   // cluster to tracklet residuals/pulls
2814   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
2815   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2816     const Int_t mdim(10);
2817     const Char_t *clTitle[mdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]", "centrality", "no. of pads"};
2818     const Int_t clNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15, AliTRDeventInfo::kCentralityClasses, kNpads};
2819     const Double_t clMin[mdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45, -0.5, 0.5},
2820                    clMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45, AliTRDeventInfo::kCentralityClasses - 0.5, kNpads+.5};
2821     st = "cluster spatial&charge resolution;";
2822     // define minimum info to be saved in non debug mode
2823     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimCl);
2824     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2825     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2826   } else H->Reset();
2827   fContainer->AddAt(H, kCluster);
2828   //++++++++++++++++++++++
2829   // tracklet to TRD track
2830   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
2831   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2832     const Int_t mdim(kNdim+8);
2833     Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2834     Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2835     Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2836     Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2837     // set specific fields
2838     trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
2839     trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
2840     trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
2841     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2842     Int_t jdim(kNdim);
2843     trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
2844     jdim++; trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12; trMin[jdim] = 25.; trMax[jdim] = 85.;
2845 //    jdim++; trTitle[jdim]=StrDup("gap [cm]"); trNbins[jdim] = 80; trMin[jdim] = 0.1; trMax[jdim] = 2.1;
2846 //     jdim++; trTitle[jdim]=StrDup("size_{0} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2847 //     jdim++; trTitle[jdim]=StrDup("pos_{0} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2848 //     jdim++; trTitle[jdim]=StrDup("size_{1} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2849 //     jdim++; trTitle[jdim]=StrDup("pos_{1} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2850
2851     st = "tracklet spatial&charge resolution;";
2852     // define minimum info to be saved in non debug mode
2853     Int_t ndim=DebugLevel()>=1?(jdim+1):kNdimTrklt;
2854     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2855     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2856   } else H->Reset();
2857   fContainer->AddAt(H, kTracklet);
2858   //++++++++++++++++++++++
2859   // tracklet to TRDin
2860   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
2861   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2862     // set specific fields
2863     const Int_t mdim(kNdim+3);
2864     Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
2865     Int_t trinNbins[mdim];   memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
2866     Double_t trinMin[mdim];  memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
2867     Double_t trinMax[mdim];  memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
2868     trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1; trinMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trinMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2869     trinTitle[kNdim]=StrDup("p [GeV/c]"); trinNbins[kNdim] = 24; trinMin[kNdim] = -0.5; trinMax[kNdim] = 23.5;
2870     trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
2871     trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
2872     st = "r-#phi/z/angular residuals @ TRD entry;";
2873     // define minimum info to be saved in non debug mode
2874     Int_t ndim=DebugLevel()>=1?mdim:kNdim;//kNdimTrkIn;
2875     for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
2876     H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
2877   } else H->Reset();
2878   fContainer->AddAt(H, kTrackIn);
2879   // tracklet to TRDout
2880 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2881
2882
2883   // Resolution histos
2884   if(!HasMCdata()) return fContainer;
2885
2886   //++++++++++++++++++++++
2887   // cluster to TrackRef residuals/pulls
2888   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
2889   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2890     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
2891     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, Int_t(kNcharge)*AliPID::kSPECIES+1, 15};
2892     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, -AliPID::kSPECIES-0.5, -45},
2893                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, AliPID::kSPECIES+0.5, 45};
2894     st = "MC cluster spatial resolution;";
2895     // define minimum info to be saved in non debug mode
2896     Int_t ndim=DebugLevel()>=1?kNdim:4;
2897     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2898     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2899   } else H->Reset();
2900   fContainer->AddAt(H, kMCcluster);
2901   //++++++++++++++++++++++
2902   // tracklet to TrackRef
2903   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
2904   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2905     Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2906     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2907     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2908     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2909     // set specific fields
2910     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2911     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2912     trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
2913     trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
2914     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2915
2916     st = "MC tracklet spatial resolution;";
2917     // define minimum info to be saved in non debug mode
2918     Int_t ndim=DebugLevel()>=1?kNdim:4;
2919     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2920     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2921   } else H->Reset();
2922   fContainer->AddAt(H, kMCtracklet);
2923   //++++++++++++++++++++++
2924   // TRDin to TrackRef
2925   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
2926   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2927     st = "MC r-#phi/z/angular residuals @ TRD entry;";
2928     // set specific fields
2929     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2930     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2931     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2932     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2933     trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
2934     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2935     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2936     // define minimum info to be saved in non debug mode
2937     Int_t ndim=DebugLevel()>=1?kNdim:7;
2938     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
2939     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2940   } else H->Reset();
2941   fContainer->AddAt(H, kMCtrackIn);
2942   //++++++++++++++++++++++
2943   // track to TrackRef
2944   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
2945   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2946     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2947     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2948     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2949     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2950     // set specific fields
2951     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2952     trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
2953     trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
2954     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2955     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2956     trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
2957
2958     st = "MC track spatial&p_{t} resolution;";
2959     // define minimum info to be saved in non debug mode
2960     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
2961     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2962     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2963   } else H->Reset();
2964   fContainer->AddAt(H, kMCtrack);
2965
2966   return fContainer;
2967 }
2968
2969 //____________________________________________________________________
2970 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2971 {
2972 //
2973 // Fit track with a staight line using the "np" clusters stored in the array "points".
2974 // The following particularities are stored in the clusters from points:
2975 //   1. pad tilt as cluster charge
2976 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
2977 // The parameters of the straight line fit are stored in the array "param" in the following order :
2978 //     param[0] - x0 reference radial position
2979 //     param[1] - y0 reference r-phi position @ x0
2980 //     param[2] - z0 reference z position @ x0
2981 //     param[3] - slope dy/dx
2982 //     param[4] - slope dz/dx
2983 //
2984 // Attention :
2985 // Function should be used to refit tracks for B=0T
2986 //
2987
2988   if(np<40){
2989     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2990     return kFALSE;
2991   }
2992   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2993
2994   Double_t x0(0.);
2995   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2996   x0/=Float_t(np);
2997
2998   Double_t x, y, z, dx, tilt(0.);
2999   for(Int_t ip(0); ip<np; ip++){
3000     x = points[ip].GetX(); z = points[ip].GetZ();
3001     dx = x - x0;
3002     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3003   }
3004   if(zfitter.Eval() != 0) return kFALSE;
3005
3006   Double_t z0    = zfitter.GetParameter(0);
3007   Double_t dzdx  = zfitter.GetParameter(1);
3008   for(Int_t ip(0); ip<np; ip++){
3009     if(points[ip].GetClusterType()) continue;
3010     x    = points[ip].GetX();
3011     dx   = x - x0;
3012     y    = points[ip].GetY();
3013     z    = points[ip].GetZ();
3014     tilt = points[ip].GetCharge();
3015     y -= tilt*(-dzdx*dx + z - z0);
3016     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3017     yfitter.AddPoint(&dx, y, 1.);
3018   }
3019   if(yfitter.Eval() != 0) return kFALSE;
3020   Double_t y0   = yfitter.GetParameter(0);
3021   Double_t dydx = yfitter.GetParameter(1);
3022
3023   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3024   if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
3025   return kTRUE;
3026 }
3027
3028 //____________________________________________________________________
3029 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3030 {
3031 //
3032 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3033 // See function FitTrack for the data stored in the "clusters" array
3034
3035 // The parameters of the straight line fit are stored in the array "param" in the following order :
3036 //     par[0] - x0 reference radial position
3037 //     par[1] - y0 reference r-phi position @ x0
3038 //     par[2] - slope dy/dx
3039 //
3040 // Attention :
3041 // Function should be used to refit tracks for B=0T
3042 //
3043
3044   TLinearFitter yfitter(2, "pol1");
3045
3046   // grep data for tracklet
3047   Double_t x0(0.), x[60], y[60], dy[60];
3048   Int_t nly(0);
3049   for(Int_t ip(0); ip<np; ip++){
3050     if(points[ip].GetClusterType()) continue;
3051     if(points[ip].GetVolumeID() != ly) continue;
3052     Float_t xt(points[ip].GetX())
3053            ,yt(param[1] + param[3] * (xt - param[0]));
3054     x[nly] = xt;
3055     y[nly] = points[ip].GetY();
3056     dy[nly]= y[nly]-yt;
3057     x0    += xt;
3058     nly++;
3059   }
3060   if(nly<10){
3061     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3062     return kFALSE;
3063   }
3064   // set radial reference for fit
3065   x0 /= Float_t(nly);
3066
3067   // find tracklet core
3068   Double_t mean(0.), sig(1.e3);
3069   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3070
3071   // simple cluster error parameterization
3072   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3073
3074   // fit tracklet core
3075   for(Int_t jly(0); jly<nly; jly++){
3076     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3077     Double_t dx(x[jly]-x0);
3078     yfitter.AddPoint(&dx, y[jly], 1.);
3079   }
3080   if(yfitter.Eval() != 0) return kFALSE;
3081   par[0] = x0;
3082   par[1] = yfitter.GetParameter(0);
3083   par[2] = yfitter.GetParameter(1);
3084   return kTRUE;
3085 }
3086
3087 //____________________________________________________________________
3088 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3089 {
3090 //
3091 // Global selection mechanism of tracksbased on cluster to fit residuals
3092 // The parameters are the same as used ni function FitTrack().
3093
3094   const Float_t kS(0.6), kM(0.2);
3095   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3096   Float_t dy, dz, s, m;
3097   for(Int_t ip(0); ip<np; ip++){
3098     if(points[ip].GetClusterType()) continue;
3099     Float_t x0(points[ip].GetX())
3100           ,y0(param[1] + param[3] * (x0 - param[0]))
3101           ,z0(param[2] + param[4] * (x0 - param[0]));
3102     dy=points[ip].GetY() - y0; h.Fill(dy);
3103     dz=points[ip].GetZ() - z0;
3104   }
3105   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3106   fg.SetParameter(1, 0.);
3107   fg.SetParameter(2, 2.e-2);
3108   h.Fit(&fg, "QN");
3109   m=fg.GetParameter(1); s=fg.GetParameter(2);
3110   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3111   return kTRUE;
3112 }
3113
3114 //________________________________________________________
3115 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3116 {
3117   //
3118   // Get the most probable value and the full width half mean
3119   // of a Landau distribution
3120   //
3121
3122   const Float_t dx = 1.;
3123   mpv = f->GetParameter(1);
3124   Float_t fx, max = f->Eval(mpv);
3125
3126   xm = mpv - dx;
3127   while((fx = f->Eval(xm))>.5*max){
3128     if(fx>max){
3129       max = fx;
3130       mpv = xm;
3131     }
3132     xm -= dx;
3133   }
3134
3135   xM += 2*(mpv - xm);
3136   while((fx = f->Eval(xM))>.5*max) xM += dx;
3137 }
3138
3139
3140 // #include "TFile.h"
3141 // //________________________________________________________
3142 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3143 // {
3144 //   if(!file){
3145 //     AliWarning("Use cluster position as in reconstruction.");
3146 //     SetLoadCorrection();
3147 //     return kTRUE;
3148 //   }
3149 //   TDirectory *cwd(gDirectory);
3150 //   TString fileList;
3151 //   FILE *filePtr = fopen(file, "rt");
3152 //   if(!filePtr){
3153 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3154 //     SetLoadCorrection();
3155 //     return kFALSE;
3156 //   }
3157 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3158 //   while(fileList.Gets(filePtr)){
3159 //     if(!TFile::Open(fileList.Data())) {
3160 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3161 //       continue;
3162 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
3163 //
3164 //     TTree *tSys = (TTree*)gFile->Get("tSys");
3165 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
3166 //     tSys->Draw("det:t>>h2", "dx", "goff");
3167 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3168 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3169 //     }
3170 //     h2->SetDirectory(cwd);
3171 //     gFile->Close();
3172 //   }
3173 //   cwd->cd();
3174 //
3175 //   if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>=2){
3176 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3177 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3178 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3179 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3180 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3181 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3182 //       printf("%03d|", idet);
3183 //       fprintf(fout, "    {");
3184 //       for(Int_t it(0); it<30; it++){
3185 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3186 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3187 //       }
3188 //       printf("\n");
3189 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3190 //     }
3191 //     fprintf(fout, "  };\n");
3192 //   }
3193 //   SetLoadCorrection();
3194 //   return kTRUE;
3195 // }
3196
3197 //________________________________________________________
3198 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
3199   :TNamed()
3200   ,fH(NULL)
3201   ,fNrebin(0)
3202   ,fRebinX(NULL)
3203   ,fRebinY(NULL)
3204 {
3205   // constructor
3206   memset(fAx, 0, 3*sizeof(Int_t));
3207   memset(fRange, 0, 4*sizeof(Float_t));
3208 }
3209
3210 //________________________________________________________
3211 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
3212 {
3213   // destructor
3214   if(fH) delete fH;
3215 }
3216
3217 //________________________________________________________
3218 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
3219 {
3220 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
3221   if(!aa[ix] || !aa[iy] || !aa[iz]) return;
3222   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
3223   // check ax definiton to protect against older versions of the data
3224   if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){
3225     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));
3226     return;
3227   }
3228   if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){
3229     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));
3230     return;
3231   }
3232   if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){
3233     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));
3234     return;
3235   }
3236   SetNameTitle(n,t);
3237   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
3238     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
3239     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
3240     az->GetNbins(), az->GetXmin(), az->GetXmax());
3241   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
3242   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
3243   AliDebug(2, Form("H3(%s, %s) :: %s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]", n, t,
3244     ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
3245     ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
3246     az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));
3247 }
3248
3249 //________________________________________________________
3250 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator=(const AliTRDresolutionProjection& rhs)
3251 {
3252 // copy projections
3253   if(this == &rhs) return *this;
3254
3255   TNamed::operator=(rhs);
3256   if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}
3257   if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);
3258   memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));
3259   memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));
3260   if(fH) delete fH;
3261   if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));
3262   return *this;
3263 }
3264
3265 //________________________________________________________
3266 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator+=(const AliTRDresolutionProjection& other)
3267 {
3268 // increment projections
3269   if(!fH || !other.fH) return *this;
3270   AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));
3271   fH->Add(other.fH);
3272   return *this;
3273 }
3274
3275 //________________________________________________________
3276 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
3277 {
3278 // increment bin with value "v" pointed by general coord in "bin"
3279   if(!fH) return;
3280   AliDebug(4, Form("  %s[%2d]", fH->GetName(), Int_t(v)));
3281   fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));
3282 }
3283
3284 //________________________________________________________
3285 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)
3286 {
3287 // build the 2D projection and adjust binning
3288
3289   const Char_t *title[] = {"Mean", "#mu", "MPV"};
3290   if(!fH) return NULL;
3291   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
3292   TH2D *h2s(NULL), *hyx(NULL);
3293   if(!(h2s = (TH2D*)fH->Project3D("yx"))) return NULL;
3294   // save a copy of the original distribution
3295   if(!del) hyx = (TH2D*)h2s->Clone(Form("%sEn", fH->GetName()));
3296   Int_t irebin(0), dxBin(1), dyBin(1);
3297   while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
3298     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
3299     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
3300     irebin++;
3301   }
3302   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
3303   delete h2s;
3304
3305   // start projection
3306   TH1 *h(NULL); Int_t n(0);
3307   Float_t dz=(fRange[1]-fRange[1])/ncol;
3308   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
3309   Int_t nt(tokenTitle->GetEntriesFast());
3310   TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
3311             Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], nt>0?(*tokenTitle)[0]->GetName():"", nt>1?(*tokenTitle)[1]->GetName():""),
3312             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
3313   h2->SetContour(ncol);
3314   h2->GetZaxis()->CenterTitle();
3315   h2->GetZaxis()->SetTitleOffset(1.4);
3316   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
3317   AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));
3318   for(Int_t iy(0); iy<ny; iy++){
3319     for(Int_t ix(0); ix<nx; ix++){
3320       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
3321       Int_t ne((Int_t)h->Integral());
3322       if(ne<nstat/2){
3323         h2->SetBinContent(ix+1, iy+1, -999);
3324         h2->SetBinError(ix+1, iy+1, 1.);
3325         n++;
3326       }else{
3327         Float_t v(h->GetMean()), ve(h->GetRMS());
3328         if(mid==1){
3329           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
3330           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
3331           h->Fit(&fg, "WQ");
3332           v = fg.GetParameter(1); ve = fg.GetParameter(2);
3333         } else if (mid==2) {
3334           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
3335           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
3336           h->Fit(&fl, "WQ");
3337           v = fl.GetMaximumX(); ve = fl.GetParameter(2);
3338 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
3339           fgle.SetParameter(0, fl.GetParameter(0));
3340           fgle.SetParameter(1, fl.GetParameter(1));
3341           fgle.SetParameter(2, fl.GetParameter(2));
3342           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
3343           h->Fit(&fgle, "WQ");
3344           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
3345         }
3346         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
3347         else h2->SetBinContent(ix+1, iy+1, v);
3348         h2->SetBinError(ix+1, iy+1, ve);
3349       }
3350     }
3351   }
3352   if(h) delete h;
3353   if(n==nx*ny){delete h2; h2=NULL;} // clean empty projections
3354   return h2;
3355 }
3356
3357 //________________________________________________________
3358 void AliTRDresolution::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)
3359 {
3360 // Normalize histo content to the mean value in the range specified by bin ranges
3361 // [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.
3362 // Optionally a threshold "thr" can be specified to disregard entries with no meaning
3363
3364   Float_t s = 0., c=0.; Int_t is(0);
3365   for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){
3366     for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){
3367       if((c = h2->GetBinContent(ix, iy))<thr) continue;
3368       s += c; is++;
3369     }
3370   }
3371   s/= is?is:1;
3372   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){
3373     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){
3374       if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-100);
3375       else h2->SetBinContent(ix, iy, 100.*(c/s-1.));
3376     }
3377   }
3378 }
3379
3380 //________________________________________________________
3381 void AliTRDresolution::SetProcesses(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3382 {
3383 // steer processes
3384   if(det) SETBIT(fSteer, kDetector); else CLRBIT(fSteer, kDetector);
3385   if(cl)  SETBIT(fSteer, kCluster); else CLRBIT(fSteer, kCluster);
3386   if(trklt) SETBIT(fSteer, kTracklet); else CLRBIT(fSteer, kTracklet);
3387   if(trkin) SETBIT(fSteer, kTrackIn); else CLRBIT(fSteer, kTrackIn);
3388 }
3389
3390 //________________________________________________________
3391 void AliTRDresolution::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)
3392 {
3393 // Set range on Z axis such to avoid outliers
3394
3395   Float_t c(0.), dz(1.e-3*(max-min));
3396   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){
3397     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){
3398       if((c = h2->GetBinContent(ix, iy))<thr) continue;
3399       if(c<=min) h2->SetBinContent(ix, iy, min+dz);
3400     }
3401   }
3402   h2->GetZaxis()->SetRangeUser(min, max);
3403 }
3404
3405 //________________________________________________________
3406 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
3407 {
3408 // define rebinning strategy for this projection
3409   fNrebin = n;
3410   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
3411   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));
3412 }
3413
3414