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?fEvent->GetMultiplicity():0.;
 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   if(!track) return NULL;
 945   // special care for EVE usage
 946   if(H && (h = (TH1*)gDirectory->Get(Form("%s_proj_%d", H->GetName(), kYrez)))) delete h;
 947   return H?H->Projection(kYrez):NULL;
 948 }
 949
 950
 951 //__________________________________________________________________________
 952 Int_t AliTRDresolution::GetPtBin(Float_t pt)
 953 {
 954 // Find pt bin according to local pt segmentation
 955   Int_t ipt(-1);
 956   while(ipt<24){
 957     if(pt<fgPtBin[ipt+1]) break;
 958     ipt++;
 959   }
 960   return ipt;
 961 }
 962
 963 //________________________________________________________
 964 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
 965 {
 966 // return mean number of entries/bin of histogram "h"
 967 // if option "opt" is given the following values are accepted:
 968 //   "<" : consider only entries less than "cut"
 969 //   ">" : consider only entries greater than "cut"
 970
 971   //Int_t dim(h->GetDimension());
 972   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
 973   Double_t sum(0.); Int_t n(0);
 974   for(Int_t ix(1); ix<=nbx; ix++)
 975     for(Int_t iy(1); iy<=nby; iy++)
 976       for(Int_t iz(1); iz<=nbz; iz++){
 977         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
 978         else{
 979           if(strcmp(opt, "<")==0) {
 980             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 981           } else if(strcmp(opt, ">")==0){
 982             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 983           } else {sum += h->GetBinContent(ix, iy, iz); n++;}
 984         }
 985       }
 986   return n>0?sum/n:0.;
 987 }
 988
 989 //________________________________________________________
 990 void AliTRDresolution::GetRangeZ(TH2 *h2, Float_t &min, Float_t &max)
 991 {
 992 // Get range on Z axis such to avoid outliers
 993
 994   Double_t cnt[20000], c, m, s;
 995   Int_t nx(h2->GetXaxis()->GetNbins()), ny(h2->GetYaxis()->GetNbins()), nc(0);
 996   for(Int_t ix(1); ix<=nx; ix++){
 997     for(Int_t iy(1); iy<=ny; iy++){
 998       if((c = h2->GetBinContent(ix, iy))<10) continue;
 999       cnt[nc++] = c;
1000       if(nc==20000) break;
1001     }
1002     if(nc==20000) break;
1003   }
1004   AliMathBase::EvaluateUni(nc, cnt, m, s, 0);
1005   min = m-s; max = m+2.*s;
1006 }
1007
1008 //________________________________________________________
1009 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1010 {
1011   //
1012   // Get the reference figures
1013   //
1014
1015   if(!gPad){
1016     AliWarning("Please provide a canvas to draw results.");
1017     return kFALSE;
1018   }
1019 /*  Int_t selection[100], n(0), selStart(0); //
1020   Int_t ly0(0), dly(5);
1021   TList *l(NULL); TVirtualPad *pad(NULL); */
1022   switch(ifig){
1023   case 0:
1024     break;
1025   }
1026   AliWarning(Form("Reference plot [%d] missing result", ifig));
1027   return kFALSE;
1028 }
1029
1030
1031 //________________________________________________________
1032 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
1033 {
1034 // Build pt segments
1035   for(Int_t j(0); j<=24; j++){
1036     pt0+=(TMath::Exp(j*j*dpt)-1.);
1037     fgPtBin[j]=pt0;
1038   }
1039 }
1040
1041 //________________________________________________________
1042 void AliTRDresolution::MakeSummary()
1043 {
1044 // Build summary plots
1045
1046   if(!fProj){
1047     AliError("Missing results");
1048     return;
1049   }
1050   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
1051   TObjArray *arr(NULL); TH2 *h2(NULL);
1052
1053   // cluster resolution
1054   // define palette
1055   gStyle->SetPalette(1);
1056   const Int_t nClViews(9);
1057   const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh"};
1058   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1};
1059   const Float_t vClSignMin[2] = {2.6e2, 4.4e2},
1060                 vClSignMax[2] = {4.4e2, 6.2e2},
1061                 vClMin[nClViews] = {3.2e2, vClSignMin[Int_t(fBsign)], vClSignMin[Int_t(!fBsign)], 0., 0., 0., 0., 0., 3.2e2},
1062                 vClMax[nClViews] = {5.e2, vClSignMax[Int_t(fBsign)], vClSignMax[Int_t(!fBsign)], 0., 0., 0., 0., 0., 5.e2};
1063   const Int_t nTrkltViews(20);
1064   const Char_t *vTrkltName[nTrkltViews] = {
1065     "TrkltY", "TrkltYn", "TrkltYp", "TrkltY", "TrkltYn", "TrkltYp",
1066     "TrkltRCZ",
1067     "TrkltPh", "TrkltPhn", "TrkltPhp",
1068     "TrkltQ", "TrkltQn", "TrkltQp",
1069     "TrkltQS", "TrkltQSn", "TrkltQSp",
1070     "TrkltTBn", "TrkltTBp", "TrkltTBn", "TrkltTBp"
1071 //    "TrkltYnl0", "TrkltYpl0", "TrkltPhnl0", "TrkltPhpl0", "TrkltQnl0", "TrkltQpl0",  // electrons low pt
1072 /*    "TrkltYnl1", "TrkltYpl1", "TrkltPhnl1", "TrkltPhpl1", "TrkltQnl1", "TrkltQpl1",  // muons low pt
1073     "TrkltYnl2", "TrkltYpl2", "TrkltPhnl2", "TrkltPhpl2", "TrkltQnl2", "TrkltQpl2"  // pions low pt*/
1074   };
1075   const UChar_t vTrkltOpt[nTrkltViews] = {
1076     0, 0, 0, 1, 1, 1,
1077     0,
1078     0, 0, 0,
1079     0, 0, 0,
1080     0, 0, 0,
1081     0, 0, 2, 2
1082   };
1083   const Int_t nTrkInViews(5);
1084   const Char_t *vTrkInName[nTrkInViews][6] = {
1085     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1086     {"TrkInY", "TrkInYn", "TrkInYp", "TrkInRCZ", "TrkInPhn", "TrkInPhp"},
1087     {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
1088     {"TrkInXnl", "TrkInXpl", "TrkInXl", "TrkInYnh", "TrkInYph", "TrkInYh"},
1089     {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"},
1090     //{"TrkInRCX", "TrkInRCY", "TrkInRCPh", "TrkInRCZl", "TrkInRCZi", "TrkInRCZh"}
1091   };
1092   const UChar_t vTrkInOpt[nTrkInViews] = {0, 1, 0, 0, 0};
1093   const Float_t min[6] = {0.15, 0.15, 0.15, 0.15, 0.5, 0.5};
1094   const Float_t max[6] = {0.6, 0.6, 0.6, 0.6, 2.3, 2.3};
1095   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]"};
1096
1097   const Int_t nTrkViews(27);
1098   const Char_t *vTrkName[nTrkViews] = {
1099     "TrkY", "TrkYn", "TrkYp",
1100     "TrkPh", "TrkPhn", "TrkPhp",
1101     "TrkDPt", "TrkDPtn", "TrkDPtp",
1102     "TrkYnl0", "TrkYpl0", "TrkPhnl0", "TrkPhpl0", "TrkDPtnl0", "TrkDPtpl0",  // electrons low pt
1103     "TrkYnl1", "TrkYpl1", "TrkPhnl1", "TrkPhpl1", "TrkDPtnl1", "TrkDPtpl1",  // muons low pt
1104     "TrkYnl2", "TrkYpl2", "TrkPhnl2", "TrkPhpl2", "TrkDPtnl2", "TrkDPtpl2"  // pions low pt
1105   };
1106   const Char_t *typName[] = {"", "MC"};
1107   const Int_t nx(2048), ny(1536);
1108
1109   if((arr = (TObjArray*)fProj->At(kDetector))){
1110     cOut = new TCanvas(Form("%s_DetOccupancy", GetName()), "Detector performance", 2*nx, 2*ny);
1111     cOut->Divide(AliTRDgeometry::kNlayer,AliTRDeventInfo::kCentralityClasses, 1.e-5, 1.e-5);
1112     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1113       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1114         p=cOut->cd(icen*AliTRDgeometry::kNlayer+ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1115         if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d%dEn", ily, icen)))) continue;
1116         SetNormZ(h2, 1, 11, 1, -1, 10.);
1117         SetRangeZ(h2, 0., 150, -25.);
1118         h2->GetZaxis()->SetTitle("Rel. Det. Occup. [%]");
1119         h2->GetZaxis()->CenterTitle();
1120         h2->Draw("colz");
1121         MakeDetectorPlot(ily, "pad");
1122       }
1123     }
1124     cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1125     cOut = new TCanvas(Form("%s_DetCharge", GetName()), "Detector performance", nx, ny);
1126     cOut->Divide(3,2, 1.e-5, 1.e-5);
1127     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1128       p=cOut->cd(ily+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1129       if(!(h2 = (TH2*)arr->FindObject(Form("HDet%d_2D", ily)))) continue;
1130       SetNormZ(h2, 1, 11, 1, -1, 10.);
1131       SetRangeZ(h2, 0., 45., -15.);
1132       h2->GetZaxis()->SetTitle("Rel. Mean(q) [%]");
1133       h2->GetZaxis()->CenterTitle();
1134       h2->Draw("colz");
1135       MakeDetectorPlot(ily, "pad");
1136     }
1137     cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1138   }
1139   for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
1140     if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
1141       for(Int_t iview(0); iview<nClViews; iview++){
1142         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vClName[iview], vClOpt[iview]), "Cluster Resolution", nx, ny);
1143         cOut->Divide(3,2, 1.e-5, 1.e-5);
1144         Int_t nplot(0);
1145         for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1146           p=cOut->cd(ily+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
1147           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], ily)))) continue;
1148           nplot++;
1149           if(vClOpt[iview]==0) h2->Draw("colz");
1150           else if(vClOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [#mum]", vClMin[iview], vClMax[iview], 1.e4);
1151           if(iview<5) MakeDetectorPlot(ily);
1152         }
1153         if(nplot==AliTRDgeometry::kNlayer) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1154         else delete cOut;
1155       }
1156     }
1157     // tracklet systematic
1158     if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
1159       for(Int_t iview(0); iview<nTrkltViews; iview++){
1160         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkltName[iview], vTrkltOpt[iview]), "Tracklet Resolution", nx, ny);
1161         cOut->Divide(3,2, 1.e-5, 1.e-5);
1162         Int_t nplot(0);
1163         for(Int_t iplot(0); iplot<6; iplot++){
1164           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1165           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
1166           nplot++;
1167           if(vTrkltOpt[iview]==0) h2->Draw("colz");
1168           else if (vTrkltOpt[iview]==1) DrawSigma(h2, "#sigma(#Deltay) [cm]", .15, .4);
1169           else if (vTrkltOpt[iview]==2) DrawSigma(h2, "#sigma(occupancy) [%]", 10.5, 15.);
1170           MakeDetectorPlot(iplot);
1171         }
1172         if(nplot==6){
1173           cOut->Modified();cOut->Update();
1174           cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1175         }
1176         delete cOut;
1177       }
1178     }
1179     // trackIn systematic
1180     if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
1181       for(Int_t iview(0); iview<nTrkInViews; iview++){
1182         cOut = new TCanvas(Form("%s_%s%s_%d", GetName(), typName[ityp], vTrkInName[iview][0], vTrkInOpt[iview]), "Track IN Resolution", nx, ny);
1183         cOut->Divide(3,2, 1.e-5, 1.e-5);
1184         Int_t nplot(0);
1185         for(Int_t iplot(0); iplot<6; iplot++){
1186           p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1187           if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))){
1188             AliInfo(Form("Missing H%s%s_2D", typName[ityp], vTrkInName[iview][iplot]));
1189             continue;
1190           }
1191           nplot++;
1192           if(vTrkInOpt[iview]==0) h2->Draw("colz");
1193           else DrawSigma(h2, ttt[iplot], min[iplot], max[iplot]);
1194           MakeDetectorPlot(0);
1195         }
1196         if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1197         else delete cOut;
1198       }
1199       // species
1200       const Float_t zmin[] = {1., 20., 25., 10., 5.},
1201                     zmax[] = {10., 38., 43., 28., 14.};
1202       cOut = new TCanvas(Form("%s_%sTrkInSpc", GetName(), typName[ityp]), "Track IN PID", Int_t(1.5*nx), Int_t(1.5*ny));
1203       cOut->Divide(5,3, 1.e-5, 1.e-5);
1204       Int_t nplot(0); const Char_t *chName[] = {"p", "n", ""};
1205       for(Int_t ich(0), ipad(1); ich<3; ich++){
1206         TH2 *h2s(NULL);
1207         if(!(h2s = (TH2*)arr->FindObject(Form("H%sTrkInY%sEn", typName[ityp], chName[ich])))) {
1208           AliInfo(Form("Missing H%sTrkIn%sEn", typName[ityp], chName[ich]));
1209           continue;
1210         }
1211         for(Int_t ispec(0); ispec<AliPID::kSPECIES; ispec++){
1212           p=cOut->cd(ipad++); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1213           if(!(h2 = (TH2*)arr->FindObject(Form("H%sTrkInY%s%dEn", typName[ityp], chName[ich], ispec)))) {
1214             AliInfo(Form("Missing H%sTrkIn%s%dEn", typName[ityp], chName[ich], ispec));
1215             continue;
1216           }
1217           nplot++;
1218           h2->Divide(h2, h2s, 1.e2);
1219           h2->SetContour(9);
1220           h2->GetZaxis()->SetRangeUser(zmin[ispec], zmax[ispec]);
1221           h2->GetZaxis()->SetTitle("Rel. Abundancy [%]");h2->GetZaxis()->CenterTitle();
1222           TString tit(h2->GetTitle()); h2->SetTitle(Form("%s :: Relative Abundancy", ((*(tit.Tokenize("::")))[0])->GetName()));
1223           h2->Draw("colz");
1224           MakeDetectorPlot(0);
1225         }
1226       }
1227       if(nplot==15) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1228       else delete cOut;
1229     }
1230   }
1231   // track MC systematic
1232   if((arr = (TObjArray*)fProj->At(kMCtrack))) {
1233     for(Int_t iview(0); iview<nTrkViews; iview++){
1234       cOut = new TCanvas(Form("%s_MC%s", GetName(), vTrkName[iview]), "Track Resolution", nx, ny);
1235       cOut->Divide(3,2, 1.e-5, 1.e-5);
1236       Int_t nplot(0);
1237       for(Int_t iplot(0); iplot<6; iplot++){
1238         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
1239         if(!(h2 = (TH2*)arr->FindObject(Form("HMC%s%d_2D", vTrkName[iview], iplot)))) continue;
1240         h2->Draw("colz"); nplot++;
1241       }
1242       if(nplot==6) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1243       else delete cOut;
1244     }
1245   }
1246
1247
1248   gStyle->SetPalette(1);
1249 }
1250
1251 //________________________________________________________
1252 void AliTRDresolution::DrawSigma(TH2 *h2, const Char_t *title, Float_t m, Float_t M, Float_t scale)
1253 {
1254   // Draw error bars scaled with "scale" instead of content values
1255   //use range [m,M] if limits are specified
1256
1257   if(!h2) return;
1258   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1259   TH2F *h2e = new TH2F(Form("%s_E", h2->GetName()),
1260                 Form("%s;%s;%s;%s", h2->GetTitle(), ax->GetTitle(), ay->GetTitle(), title),
1261                 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
1262                 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1263   h2e->SetContour(9);
1264   TAxis *az(h2e->GetZaxis());
1265   if(M>m) az->SetRangeUser(m, M);
1266   az->CenterTitle();
1267   az->SetTitleOffset(1.5);
1268   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1269     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1270       if(h2->GetBinContent(ix, iy)<-100.) continue;
1271       Float_t v(scale*h2->GetBinError(ix, iy));
1272       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1273       h2e->SetBinContent(ix, iy, v);
1274     }
1275   }
1276   h2e->Draw("colz");
1277 }
1278
1279 //________________________________________________________
1280 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1281 {
1282 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1283 // The "range" vector should be initialized with 2 elements
1284 // Option "mod" can be any of
1285 //   - 0 : gaussian like distribution
1286 //   - 1 : tailed distribution
1287
1288   Int_t nx(h2->GetNbinsX())
1289       , ny(h2->GetNbinsY())
1290       , n(nx*ny);
1291   Double_t *data=new Double_t[n];
1292   for(Int_t ix(1), in(0); ix<=nx; ix++){
1293     for(Int_t iy(1); iy<=ny; iy++)
1294       data[in++] = h2->GetBinContent(ix, iy);
1295   }
1296   Double_t mean, sigm;
1297   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1298
1299   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1300   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1301   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1302   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1303   h1.FillN(n,data,0);
1304   delete [] data;
1305
1306   switch(mod){
1307   case 0:// gaussian distribution
1308   {
1309     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1310     h1.Fit(&fg, "QN");
1311     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1312     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1313     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1314     break;
1315   }
1316   case 1:// tailed distribution
1317   {
1318     Int_t bmax(h1.GetMaximumBin());
1319     Int_t jBinMin(1), jBinMax(100);
1320     for(Int_t ibin(bmax); ibin--;){
1321       if(h1.GetBinContent(ibin)<1.){
1322         jBinMin=ibin; break;
1323       }
1324     }
1325     for(Int_t ibin(bmax); ibin++;){
1326       if(h1.GetBinContent(ibin)<1.){
1327         jBinMax=ibin; break;
1328       }
1329     }
1330     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1331     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1332     break;
1333   }
1334   }
1335
1336   return;
1337 }
1338
1339 //________________________________________________________
1340 Bool_t AliTRDresolution::MakeProjectionDetector()
1341 {
1342 // Analyse cluster
1343   const Int_t kNcontours(9);
1344   const Int_t kNstat(100);
1345   if(fProj && fProj->At(kDetector)) return kTRUE;
1346   if(!fContainer){
1347     AliError("Missing data container.");
1348     return kFALSE;
1349   }
1350   THnSparse *H(NULL);
1351   if(!(H = (THnSparse*)fContainer->FindObject("hDet2Cluster"))){
1352     AliInfo(Form("Missing/Wrong data @ hDet2Cluster."));
1353     return kTRUE;
1354   }
1355   Int_t ndim(H->GetNdimensions());
1356   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1357   TAxis *aa[kNdim], *an(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1358   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1359   if(ndim < 5) aa[4] = new TAxis(1, -0.5, 0.5);
1360   Int_t nPad(1);
1361   if(ndim > 5){
1362     an = H->GetAxis(5);
1363     nPad = kNpads+1;
1364   }
1365   // build list of projections
1366   const Int_t nsel(8*AliTRDgeometry::kNlayer*AliTRDeventInfo::kCentralityClasses);
1367   // define rebinning strategy
1368   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 2, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1369   //const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1370   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1371   AliTRDresolutionProjection hp[kDetNproj];  TObjArray php(kDetNproj);
1372   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1373   for(Int_t ipad(0); ipad<nPad; ipad++){
1374     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1375       for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1376         isel++; // new selection
1377         hp[ih].Build(Form("HDet%d%d%d", ily, icen, ipad),
1378                     Form("Detectors :: Ly[%d] Cen[%s] Pads[%s]", ily, cenName[icen], ipad?(ipad<kNpads?Form("%d", ipad+1):Form(">%d", kNpads)):"deconv"),
1379                     kEta, kPhi, 4, aa);
1380         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1381         hp[ih].SetShowRange(10., 55.);
1382         php.AddLast(&hp[ih++]); np[isel]++;
1383       }
1384     }
1385   }
1386   AliInfo(Form("Build %3d 3D projections.", ih));
1387
1388   Int_t ly(0), cen(0), npad(0);
1389   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1390     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1391     ly = coord[kBC]-1;    // layer selection
1392     cen = coord[kYrez]-1; // centrality selection
1393     npad = 0;             // no. of pads selection
1394     if(an) npad = TMath::Min(coord[5]-1, Int_t(kNpads));
1395     isel = npad*AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer+cen*AliTRDgeometry::kNlayer+ly;
1396     ((AliTRDresolutionProjection*)php.At(isel))->Increment(coord, v);
1397     //Int_t ioff=isel;for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1398   }
1399   TObjArray *arr(NULL);
1400   fProj->AddAt(arr = new TObjArray(kDetNproj), kDetector);
1401
1402   TH2 *h2(NULL);  Int_t jh(0);
1403   for(; ih--; ){
1404     if(!hp[ih].fH) continue;
1405     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1406     arr->AddAt(h2, jh++);
1407     if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].fH->GetName())))) continue;
1408     arr->AddAt(h2, jh++);
1409   }
1410   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1411   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1412     for(Int_t icen(0); icen<AliTRDeventInfo::kCentralityClasses; icen++){
1413       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, 0)))){
1414         for(Int_t ipad(1); ipad<nPad; ipad++){
1415           if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, icen, ipad)))){
1416             (*pr0)+=(*pr1);
1417           }
1418         }
1419         pr0->fH->SetNameTitle(Form("HDet%d%d", ily, icen), Form("Detectors :: Ly[%d] Cen[%s]", ily, cenName[icen]));
1420         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0, kFALSE))) arr->AddAt(h2, jh++);
1421         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
1422         if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))) (*pr1)+=(*pr0);
1423       }
1424     }
1425     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HDet%d%d%d", ily, 0, 0)))){
1426       pr0->fH->SetNameTitle(Form("HDet%d", ily), Form("Detectors :: Ly[%d]", ily));
1427       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1428     }
1429   }
1430   AliInfo(Form("Done %3d 2D projections.", jh));
1431   return kTRUE;
1432 }
1433
1434 //________________________________________________________
1435 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1436 {
1437 // Analyse cluster
1438   const Int_t kNcontours(9);
1439   const Int_t kNstat(100);
1440   Int_t cidx=mc?kMCcluster:kCluster;
1441   if(fProj && fProj->At(cidx)) return kTRUE;
1442   if(!fContainer){
1443     AliError("Missing data container.");
1444     return kFALSE;
1445   }
1446   const Char_t *projName[] = {"hCluster2Track", "hCluster2MC"};
1447   THnSparse *H(NULL);
1448   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1449     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1450     return kFALSE;
1451   }
1452   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimCl));
1453   Int_t coord[10]; memset(coord, 0, sizeof(Int_t) * 10); Double_t v = 0.;
1454   TAxis *aa[kNdim], *as(NULL), *apt(NULL), *acen(NULL), *anp(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1455   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1456   if(ndim > Int_t(kPt)) apt = H->GetAxis(kPt);
1457   if(ndim > Int_t(kSpeciesChgRC)) as  = H->GetAxis(kSpeciesChgRC);
1458   if(ndim > Int_t(kNdim)) acen  = H->GetAxis(kNdim);
1459   if(ndim > Int_t(kNdim)+1) anp  = H->GetAxis(kNdim+1);
1460   // calculate size depending on debug level
1461   const Int_t nCh(apt?2:1);
1462   const Int_t nCen(acen?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1463   const Int_t nNpad(anp?(Int_t(kNpads)+1):1);
1464
1465   // build list of projections
1466   const Int_t nsel(AliTRDgeometry::kNlayer*kNcharge*(kNpads+1)*AliTRDeventInfo::kCentralityClasses);
1467   // define rebinning strategy
1468   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1469   AliTRDresolutionProjection hp[kClNproj];  TObjArray php(kClNproj);
1470   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1471   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1472   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1473   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1474     for(Int_t ich(0); ich<nCh; ich++){
1475       for(Int_t icen(0); icen<nCen; icen++){
1476         for(Int_t ipad(0); ipad<nNpad; ipad++){
1477           isel++; // new selection
1478           hp[ih].Build(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1479                       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"),
1480                       kEta, kPhi, kYrez, aa);
1481           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1482           php.AddLast(&hp[ih++]); np[isel]++;
1483           if(!debug) break;
1484           hp[ih].Build(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1485                       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"),
1486                       kEta, kPhi, kSpeciesChgRC, aa);
1487           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1488           hp[ih].SetShowRange(24., 33.);
1489           php.AddLast(&hp[ih++]); np[isel]++;
1490           hp[ih].Build(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1491                       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"),
1492                       kPrez, kZrez, kYrez, aa);
1493           php.AddLast(&hp[ih++]); np[isel]++;
1494           hp[ih].Build(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, ipad),
1495                         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"),
1496                         kPrez, kPt, kYrez, aa);
1497           php.AddLast(&hp[ih++]); np[isel]++;
1498         }
1499       }
1500     }
1501   }
1502   AliInfo(Form("Build %3d 3D projections.", ih));
1503
1504   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1505   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1), ioff(0), cen(0), npad(0);
1506   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1507     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1508     ly = coord[kBC]-1;
1509     // RC selection
1510     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1511
1512     // charge selection
1513     ch = 0; // [-] track
1514     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1515     cen = 0; // centrality selection
1516     if(acen) cen = coord[kNdim]-1;
1517     npad = 0;             // no. of pads selection
1518     if(anp) npad = TMath::Min(coord[kNdim+1]-1, Int_t(kNpads));
1519
1520     if(debug){
1521       isel = ly*nCh*nCen*nNpad
1522             +ch*nCen*nNpad
1523             +cen*nNpad
1524             +npad;
1525       ioff=isel*4;
1526     } else {
1527       isel = ly; ioff = isel;
1528     }
1529     if(ioff>=ih){
1530       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1531       return kFALSE;
1532     }
1533     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
1534       AliError(Form("Missing projection %d", ioff));
1535       return kFALSE;
1536     }
1537     if(strcmp(pr0->fH->GetName(), Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ch], ly, cen, npad))!=0){
1538       AliError(Form("Projection mismatch :: request[H%sClY%c%d%d%d] found[%s]", mc?"MC":"", chName[ch], ly, cen, npad, pr0->fH->GetName()));
1539       return kFALSE;
1540     }
1541     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1542   }
1543   TObjArray *arr(NULL);
1544   fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1545
1546   TH2 *h2(NULL);  Int_t jh(0);
1547   for(; ih--; ){
1548     if(!hp[ih].fH) continue;
1549     if(strchr(hp[ih].fH->GetName(), 'Q')){
1550       if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours, 0, kFALSE))) continue;
1551       arr->AddAt(h2, jh++);
1552       if(!(h2 = (TH2*)gDirectory->Get(Form("%sEn", hp[ih].fH->GetName())))) continue;
1553       arr->AddAt(h2, jh++);
1554     } else {
1555       if((h2 = hp[ih].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1556     }
1557   }
1558   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1559     for(Int_t ich(0); ich<nCh; ich++){
1560       for(Int_t icen(0); icen<nCen; icen++){
1561         /*!dy*/
1562         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1563           for(Int_t ipad(1); ipad<nNpad; ipad++){
1564             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1565             (*pr0)+=(*pr1);
1566           }
1567           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]));
1568           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1569           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1570         }
1571         /*!Q*/
1572         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1573           for(Int_t ipad(1); ipad<nNpad; ipad++){
1574             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1575             (*pr0)+=(*pr1);
1576           }
1577           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]));
1578           if((h2 = pr0->Projection2D(kNstat, kNcontours, 2, kFALSE))) arr->AddAt(h2, jh++);
1579           if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
1580           pr0->fH->SetName(Form("H%sClQS%c%d%d", mc?"MC":"", chName[ich], ily, icen));
1581           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1582           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1583         }
1584         /*!YXTC*/
1585         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1586           for(Int_t ipad(1); ipad<nNpad; ipad++){
1587             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1588             (*pr0)+=(*pr1);
1589           }
1590           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]));
1591           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1592           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1593         }
1594         /*!YXPh*/
1595         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, icen, 0)))){
1596           for(Int_t ipad(1); ipad<nNpad; ipad++){
1597             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[1], ily, icen, ipad)))) continue;
1598             (*pr0)+=(*pr1);
1599           }
1600           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]));
1601           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1602           if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))) (*pr1)+=(*pr0);
1603         }
1604       } // end centrality integration
1605       /*!dy*/
1606       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1607         pr0->fH->SetNameTitle(Form("H%sClY%c%d", mc?"MC":"", chName[ich], ily),
1608                               Form("Clusters[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
1609         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1610         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1611       }
1612       /*!Q*/
1613       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1614         pr0->fH->SetNameTitle(Form("H%sClQ%c%d", mc?"MC":"", chName[ich], ily),
1615                               Form("Clusters[%c]:: Q Ly[%d]", chSgn[ich], ily));
1616         if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1617         pr0->fH->SetName(Form("H%sClQS%c%d", mc?"MC":"", chName[ich], ily));
1618         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1619         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClQ%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1620       }
1621       /*!YXTC*/
1622       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1623         pr0->fH->SetNameTitle(Form("H%sClYXTC%c%d", mc?"MC":"", chName[ich], ily),
1624                               Form("Clusters[%c]:: #Deltay(x,TC) Ly[%d]", chSgn[ich], ily));
1625         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1626         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXTC%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1627       }
1628       /*!YXPh*/
1629       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[ich], ily, 0, 0)))){
1630         pr0->fH->SetNameTitle(Form("H%sClYXPh%c%d", mc?"MC":"", chName[ich], ily),
1631                               Form("Clusters[%c]:: #Deltay(x,#Phi) Ly[%d]", chSgn[ich], ily));
1632         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1633         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))) (*pr1)+=(*pr0);
1634       }
1635     } // end charge integration
1636     /*!dy*/
1637     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClY%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1638       pr0->fH->SetNameTitle(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1639       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1640     }
1641     /*!YXPh*/
1642     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sClYXPh%c%d%d%d", mc?"MC":"", chName[0], ily, 0, 0)))){
1643       pr0->fH->SetNameTitle(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: #Deltay Ly[%d]", ily));
1644       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1645     }
1646
1647   }
1648   AliInfo(Form("Done %3d 2D projections.", jh));
1649   return kTRUE;
1650 }
1651
1652 //________________________________________________________
1653 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1654 {
1655 // Analyse tracklet
1656   const Int_t kNcontours(9);
1657   const Int_t kNstat(30);
1658   const Int_t kNstatQ(30);
1659   Int_t cidx=mc?kMCtracklet:kTracklet;
1660   if(fProj && fProj->At(cidx)) return kTRUE;
1661   if(!fContainer){
1662     AliError("Missing data container.");
1663     return kFALSE;
1664   }
1665   const Char_t *projName[] = {"hTracklet2Track", "hTracklet2MC"};
1666   THnSparse *H(NULL);
1667   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
1668     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
1669     return kFALSE;
1670   }
1671   const Int_t mdim(kNdim+8);
1672   Int_t ndim(H->GetNdimensions()); Bool_t debug(ndim>Int_t(kNdimTrklt));
1673   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
1674   TAxis *aa[mdim], *as(NULL), *ap(NULL), *ac(NULL); memset(aa, 0, sizeof(TAxis*) * mdim);
1675   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1676   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC); // init species/charge selection
1677   if(ndim > Int_t(kPt))           ap = H->GetAxis(kPt);           // init pt selection
1678   if(ndim > Int_t(kNdim))         ac = H->GetAxis(kNdim);         // init centrality selection
1679   // calculate size depending on debug level
1680   const Int_t nCen(debug?Int_t(AliTRDeventInfo::kCentralityClasses):1);
1681   const Int_t nPt(debug?Int_t(kNpt):1);
1682   const Int_t nSpc(1);//ndim>kNdimTrklt?fgkNbins[kSpeciesChgRC]:1);
1683   const Int_t nCh(debug?Int_t(kNcharge):1);
1684
1685   // build list of projections
1686   const Int_t nsel(AliTRDeventInfo::kCentralityClasses*AliTRDgeometry::kNlayer*kNpt*(AliPID::kSPECIES*kNcharge + 1));
1687   // define rebinning strategy
1688   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1689   AliTRDresolutionProjection hp[kTrkltNproj]; TObjArray php(kTrkltNproj);
1690   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1691   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
1692   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
1693   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"};
1694 //  const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"0-10%", "10-20%", "20-50%", "50-80%", "80-100%"};
1695   const Char_t *cenName[AliTRDeventInfo::kCentralityClasses] = {"2800-inf", "2100-2799", "1400-2099", "700-1399", "0-699"};
1696   for(Int_t icen(0); icen<nCen; icen++){
1697     for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1698       for(Int_t ipt(0); ipt<nPt; ipt++){
1699         for(Int_t isp(0); isp<nSpc; isp++){
1700           for(Int_t ich(0); ich<nCh; ich++){
1701             isel++; // new selection
1702             hp[ih].Build(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1703                         Form("Tracklets[%s%c]:: #Deltay{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1704                         kEta, kPhi, kYrez, aa);
1705             //hp[ih].SetShowRange(-0.1,0.1);
1706             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1707             php.AddLast(&hp[ih++]); np[isel]++;
1708             hp[ih].Build(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1709                         Form("Tracklets[%s%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1710                         kEta, kPhi, kPrez, aa);
1711             //hp[ih].SetShowRange(-0.5,0.5);
1712             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1713             php.AddLast(&hp[ih++]); np[isel]++;
1714             hp[ih].Build(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1715                         Form("Tracklets[%s%c]:: dQdl{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1716                         kEta, kPhi, kZrez, aa);
1717             hp[ih].SetShowRange(1.,2.3);
1718             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1719             php.AddLast(&hp[ih++]); np[isel]++;
1720             hp[ih].Build(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen),
1721                         Form("Tracklets[%s%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily, cenName[icen]),
1722                         kEta, kPhi, kNdim+1, aa);
1723             hp[ih].SetShowRange(30., 70.);
1724             hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1725             php.AddLast(&hp[ih++]); np[isel]++;
1726           }
1727         }
1728         if(ndim==kNdimTrklt) continue;
1729
1730         isel++; // new selection
1731         hp[ih].Build(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1732                     Form("Tracklets[RC]:: #Deltaz{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1733                     kEta, kPhi, kZrez, aa);
1734   //      hp[ih].SetShowRange(-0.1,0.1);
1735         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1736         php.AddLast(&hp[ih++]); np[isel]++;
1737 /*        hp[ih].Build(Form("H%sTrkltRCY%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1738                     Form("Tracklets[RC]:: #Deltay{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1739                     kEta, kPhi, kYrez, aa);
1740         //hp[ih].SetShowRange(-0.1,0.1);
1741         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1742         php.AddLast(&hp[ih++]); np[isel]++;
1743         hp[ih].Build(Form("H%sTrkltRCPh%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1744                     Form("Tracklets[RC]:: #Delta#phi{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1745                     kEta, kPhi, kPrez, aa);
1746         //hp[ih].SetShowRange(-0.1,0.1);
1747         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1748         php.AddLast(&hp[ih++]); np[isel]++;
1749         hp[ih].Build(Form("H%sTrkltRCQ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen),
1750                     Form("Tracklets[RC]:: dQdl{%s} Ly[%d] Cen[%s]", ptCut[ipt], ily, cenName[icen]),
1751                     kEta, kPhi, kNdim, aa);
1752         //hp[ih].SetShowRange(-0.1,0.1);
1753         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1754         php.AddLast(&hp[ih++]); np[isel]++;*/
1755       }
1756     }
1757   }
1758   AliInfo(Form("Build %3d 3D projections.", ih));
1759
1760   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
1761   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);
1762   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1763     v = H->GetBinContent(ib, coord);
1764     if(v<1.) continue;
1765     ly = coord[kBC]-1; // layer selection
1766     // charge selection
1767     ch = 0; sp=0;// [e-] track [dafault]
1768     if(rcBin>0){ // debug mode in which species/charge are also saved
1769       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1770       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1771       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1772     }
1773     // pt selection
1774     pt = 0; // low pt
1775     if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
1776     // centrality selection
1777     cen = 0; // default
1778     if(ac) cen = coord[kNdim]-1;
1779     // global selection
1780     if(ndim==kNdimTrklt){
1781       ioff = ly*4;
1782       isel = ly;
1783     } else {
1784       isel = cen*AliTRDgeometry::kNlayer*nPt*jspc+ly*nPt*jspc+pt*jspc; isel+=sp<0?(nSpc*nCh):ch;
1785       ioff = cen*AliTRDgeometry::kNlayer*nPt*kspc+ly*nPt*kspc+pt*kspc; ioff+=sp<0?((nSpc*nCh)*4):(ch*4);
1786     }
1787     if(ioff>=ih){
1788       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
1789       return kFALSE;
1790     }
1791     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
1792       AliError(Form("Missing projection %d", ioff));
1793       return kFALSE;
1794     }
1795     if(sp>=0){
1796       if(strcmp(pr0->fH->GetName(), Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ch], ptName[pt], sp, ly, cen))!=0){
1797         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()));
1798         return kFALSE;
1799       }
1800     } else {
1801       if(strcmp(pr0->fH->GetName(), Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[pt], ly, cen))!=0){
1802         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c%d%d] found[%s]", mc?"MC":"", ptName[pt], ly, cen, pr0->fH->GetName()));
1803         return kFALSE;
1804       }
1805     }
1806     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
1807   }
1808   TObjArray *arr(NULL);
1809   fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1810
1811   TH2 *h2(NULL); Int_t jh(0);
1812   for(; ih--; ){
1813     if(!hp[ih].fH) continue;
1814     Int_t mid(0), nstat(kNstat);
1815     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=kNstatQ;}
1816     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1817     arr->AddAt(h2, jh++);
1818   }
1819   // build combined performance plots
1820   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1821     for(Int_t ich(0); ich<nCh; ich++){
1822       for(Int_t icen(0); icen<nCen; icen++){
1823         for(Int_t ipt(0); ipt<nPt; ipt++){
1824           /*!dy*/
1825           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1826             for(Int_t isp(1); isp<nSpc; isp++){
1827               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1828               (*pr0)+=(*pr1);
1829             }
1830             pr0->fH->SetNameTitle(Form("H%sTrkltY%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1831                                       Form("Tracklets[%c]:: #Deltay{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1832             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1833             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);
1834           }
1835           /*!dphi*/
1836           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1837             for(Int_t isp(1); isp<nSpc; isp++){
1838               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1839               (*pr0)+=(*pr1);
1840             }
1841             pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1842                                       Form("Tracklets[%c]:: #Delta#phi{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1843             if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1844             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);
1845           }
1846           /*!dQ/dl*/
1847           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1848             for(Int_t isp(1); isp<nSpc; isp++){
1849               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1850               (*pr0)+=(*pr1);
1851             }
1852             pr0->fH->SetNameTitle(Form("H%sTrkltQ%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             if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1855             pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1856                                       Form("Tracklets[%c]:: dQdl{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1857             pr0->SetShowRange(2.4, 5.1);
1858             if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1859             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);
1860           }
1861           /*!TB occupancy*/
1862           if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], 0, ily, icen)))){
1863             for(Int_t isp(1); isp<nSpc; isp++){
1864               if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[ipt], isp, ily, icen)))) continue;
1865               (*pr0)+=(*pr1);
1866             }
1867             pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%c%d%d", mc?"MC":"", chName[ich], ptName[ipt], ily, icen),
1868                                       Form("Tracklets[%c]:: OccupancyTB{%s} Ly[%d] Cen[%s]", chSgn[ich], ptCut[ipt], ily, cenName[icen]));
1869             if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1870             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);
1871           }
1872         } // end pt integration
1873         /*!dy*/
1874         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1875           pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1876                                 Form("Tracklets[%c]:: #Deltay Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1877           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1878           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);
1879         }
1880         /*!dphi*/
1881         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1882           pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1883                                 Form("Tracklets[%c]:: #Delta#phi Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1884           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1885           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);
1886         }
1887         /*!dQ/dl*/
1888         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1889           pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1890                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1891           pr0->SetShowRange(1.,2.3);
1892           if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1893           pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1894                                 Form("Tracklets[%c]:: dQdl Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1895           pr0->SetShowRange(2.4,5.1);
1896           if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1897           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);
1898         }
1899         /*!TB occupancy*/
1900         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, icen)))){
1901           pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d%d", mc?"MC":"", chName[ich], ily, icen),
1902                                 Form("Tracklets[%c]:: OccupancyTB Ly[%d] Cen[%s]", chSgn[ich], ily, cenName[icen]));
1903           if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1904           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);
1905         }
1906       } // end centrality integration
1907       /*!dy*/
1908       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1909         pr0->fH->SetNameTitle(Form("H%sTrkltY%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Deltay Ly[%d]", chSgn[ich], ily));
1910         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1911         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);
1912       }
1913       /*!dphi*/
1914       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1915         pr0->fH->SetNameTitle(Form("H%sTrkltPh%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: #Delta#phi Ly[%d]", chSgn[ich], ily));
1916         pr0->SetShowRange(-.9,.9);
1917         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1918         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);
1919       }
1920       /*!dQ/dl*/
1921       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1922         pr0->fH->SetNameTitle(Form("H%sTrkltQ%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1923         pr0->SetShowRange(1.,2.3);
1924         if((h2 = pr0->Projection2D(kNstatQ, kNcontours, 2))) arr->AddAt(h2, jh++);
1925         pr0->fH->SetNameTitle(Form("H%sTrkltQS%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: dQdl Ly[%d]", chSgn[ich], ily));
1926         pr0->SetShowRange(2.4,5.1);
1927         if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1928         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);
1929       }
1930       /*!TB occupancy*/
1931       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[ich], ptName[0], 0, ily, 0)))){
1932         pr0->fH->SetNameTitle(Form("H%sTrkltTB%c%d", mc?"MC":"", chName[ich], ily), Form("Tracklets[%c] :: OccupancyTB Ly[%d]", chSgn[ich], ily));
1933         if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1934         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);
1935       }
1936     } // end charge integration
1937     /*!dy*/
1938     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltY%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1939       pr0->fH->SetNameTitle(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: #Deltay Ly[%d]", ily));
1940       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1941     }
1942     /*!dphi*/
1943     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltPh%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1944       pr0->fH->SetNameTitle(Form("H%sTrkltPh%d", mc?"MC":"", ily), Form("Tracklets :: #Delta#phi Ly[%d]", ily));
1945       pr0->SetShowRange(-.45,.45);
1946       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1947     }
1948     /*!dQdl*/
1949     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltQ%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1950       pr0->fH->SetNameTitle(Form("H%sTrkltQ%d", mc?"MC":"", ily), Form("Tracklets :: dQdl Ly[%d]", ily));
1951       pr0->SetShowRange(1.,2.3);
1952       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
1953       pr0->fH->SetName(Form("H%sTrkltQS%d", mc?"MC":"", ily));
1954       pr0->SetShowRange(2.4,5.1);
1955       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
1956     }
1957     /*!TB occupancy*/
1958     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltTB%c%c%d%d%d", mc?"MC":"", chName[0], ptName[0], 0, ily, 0)))){
1959       pr0->fH->SetNameTitle(Form("H%sTrkltTB%d", mc?"MC":"", ily), Form("Tracklets :: OccupancyTB Ly[%d]", ily));
1960       if((h2 = pr0->Projection2D(kNstat, kNcontours))) arr->AddAt(h2, jh++);
1961     }
1962
1963     /*! Row Cross processing*/
1964     for(Int_t icen(0); icen<nCen; icen++){
1965       /*!RC dz*/
1966       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, icen)))){
1967         for(Int_t ipt(0); ipt<kNpt; ipt++){
1968           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[ipt], ily, icen)))) continue;
1969           (*pr0)+=(*pr1);
1970         }
1971         pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d%d", mc?"MC":"", ily, icen), Form("Tracklets[RC]:: #Deltaz Ly[%d] Cen[%s]", ily, cenName[icen]));
1972         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1973         if(icen && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))) (*pr1)+=(*pr0);
1974       }
1975     } // end centrality integration for row cross
1976     /*!RC dz*/
1977     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkltRCZ%c%d%d", mc?"MC":"", ptName[0], ily, 0)))){
1978       pr0->fH->SetNameTitle(Form("H%sTrkltRCZ%d", mc?"MC":"", ily), Form("Tracklets[RC] :: #Deltaz Ly[%d]", ily));
1979       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
1980     }
1981   } // end layer loop
1982   AliInfo(Form("Done %3d 2D projections.", jh));
1983   return kTRUE;
1984 }
1985
1986 //________________________________________________________
1987 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1988 {
1989 // Analyse track in
1990   const Int_t kNcontours(9);
1991   const Int_t kNstat(30);
1992   Int_t cidx=mc?kMCtrackIn:kTrackIn;
1993   if(fProj && fProj->At(cidx)) return kTRUE;
1994   if(!fContainer){
1995     AliError("Missing data container.");
1996     return kFALSE;
1997   }
1998   const Char_t *projName[] = {"hTracklet2TRDin", "hTRDin2MC"};
1999   THnSparse *H(NULL);
2000   if(!(H = (THnSparse*)fContainer->FindObject(projName[Int_t(mc)]))){
2001     AliError(Form("Missing/Wrong data @ %s.", projName[Int_t(mc)]));
2002     return kFALSE;
2003   }
2004
2005   const Int_t mdim(kNdim+3);
2006   Int_t coord[mdim]; memset(coord, 0, sizeof(Int_t) * mdim); Double_t v = 0.;
2007   Int_t ndim(H->GetNdimensions());
2008   TAxis *aa[mdim], *as(NULL), *ap(NULL), *ax(NULL), *abf(NULL); memset(aa, 0, sizeof(TAxis*) * (mdim));
2009   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2010   if(ndim > Int_t(kSpeciesChgRC)) as = H->GetAxis(kSpeciesChgRC);
2011   if(ndim > Int_t(kPt))           ap = H->GetAxis(kPt);
2012   if(ndim > Int_t(kNdim)+1)       ax = H->GetAxis(kNdim+1);
2013   if(ndim > Int_t(kNdim)+2)      abf = H->GetAxis(kNdim+2);
2014   //AliInfo(Form("Using : Species[%c] Pt[%c] BunchFill[%c]", as?'y':'n', ap?'y':'n', abf?'y':'n'));
2015   const Int_t nPt(ndim>Int_t(kNdimTrkIn)?Int_t(kNpt):1);
2016
2017   // build list of projections
2018   const Int_t nsel(kNpt*(AliPID::kSPECIES*kNcharge + 1));
2019   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2020   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2021   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"};
2022   // define rebinning strategy
2023   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
2024   AliTRDresolutionProjection hp[kMCTrkInNproj]; TObjArray php(kMCTrkInNproj+2);
2025   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2026   // define list of projections
2027   for(Int_t ipt(0); ipt<nPt; ipt++){
2028     for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2029       for(Int_t ich(0); ich<kNcharge; ich++){
2030         isel++; // new selection
2031         hp[ih].Build(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2032                      Form("TrackIn[%s%c]:: #Deltay{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2033                      kEta, kPhi, kYrez, aa);
2034         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2035         php.AddLast(&hp[ih++]); np[isel]++;
2036         hp[ih].Build(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2037                      Form("TrackIn[%s%c]:: #Delta#phi{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2038                      kEta, kPhi, kPrez, aa);
2039         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2040         php.AddLast(&hp[ih++]); np[isel]++;
2041         hp[ih].Build(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2042                      Form("TrackIn[%s%c]:: dQdl {%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2043                      kEta, kPhi, kZrez, aa);
2044         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2045         php.AddLast(&hp[ih++]); np[isel]++;
2046         if(!ax) continue;
2047         hp[ih].Build(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp),
2048                      Form("TrackIn[%s%c]:: #Deltax{%s}", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt]),
2049                      kEta, kPhi, kNdim+1, aa);
2050         hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2051         php.AddLast(&hp[ih++]); np[isel]++;
2052       }
2053     }
2054     isel++; // RC tracks
2055     hp[ih].Build(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt]),
2056                   Form("TrackIn[RC]:: #Deltaz{%s}", ptCut[ipt]),
2057                   kEta, kPhi, kZrez, aa);
2058     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2059     php.AddLast(&hp[ih++]); np[isel]++;
2060     hp[ih].Build(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt]),
2061                   Form("TrackIn[RC]:: #Deltay{%s}", ptCut[ipt]),
2062                   kEta, kPhi, kYrez, aa);
2063     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2064     php.AddLast(&hp[ih++]); np[isel]++;
2065     hp[ih].Build(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt]),
2066                   Form("TrackIn[RC]:: #Delta#phi{%s}", ptCut[ipt]),
2067                   kEta, kPhi, kPrez, aa);
2068     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2069     php.AddLast(&hp[ih++]); np[isel]++;
2070     if(!ax) continue;
2071     hp[ih].Build(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt]),
2072                   Form("TrackIn[RC]:: #Deltax{%s}", ptCut[ipt]),
2073                   kEta, kPhi, kNdim+1, aa);
2074     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2075     php.AddLast(&hp[ih++]); np[isel]++;
2076   }
2077   AliInfo(Form("Build %3d 3D projections.", ih));
2078
2079   // fill projections
2080   Int_t ch(0), pt(0), sp(1), rcBin(as?as->FindBin(0.):-1), ioff(0);
2081   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
2082   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2083     v = H->GetBinContent(ib, coord);
2084     if(v<1.) continue;
2085     if(fBCbinTOF>0 && coord[kBC]!=fBCbinTOF) continue; // TOF bunch cross cut
2086     if(fBCbinFill>0 && abf && coord[kNdim+2]!=fBCbinTOF) continue; // Fill bunch cut
2087     // charge selection
2088     ch = 0; sp=1;// [-] track
2089     if(rcBin>0){ // debug mode in which species are also saved
2090       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2091       if(sp>=AliPID::kSPECIES){
2092         AliDebug(2, Form("Wrong SpeciesIndex[%d]. Rescale", sp));
2093         sp = AliPID::kSPECIES-1;
2094       }
2095       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
2096       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2097     }
2098     // pt selection
2099     pt = 0; // low pt
2100     if(ap) pt = TMath::Min(coord[kPt]-1, Int_t(kNpt)-1);
2101     // global selection
2102     isel = pt*(AliPID::kSPECIES*kNcharge+1); isel+=sp<0?(AliPID::kSPECIES*kNcharge):(sp*kNcharge+ch);
2103     ioff = isel*(ax?4:3);
2104     if(ioff>=ih){
2105       AliError(Form("Wrong selection %d [%3d]", ioff, ih));
2106       return kFALSE;
2107     }
2108     if(!(pr0=(AliTRDresolutionProjection*)php.At(ioff))) {
2109       AliError(Form("Missing projection %d", ioff));
2110       return kFALSE;
2111     }
2112     if(sp>=0 && ch<2){
2113       if(strcmp(pr0->fH->GetName(), Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ch], ptName[pt], sp))!=0){
2114         AliError(Form("Projection mismatch :: request[H%sTrkInY%c%c%d] found[%s]", mc?"MC":"", chName[ch], ptName[pt], sp, pr0->fH->GetName()));
2115         return kFALSE;
2116       }
2117     } else {
2118       if(strcmp(pr0->fH->GetName(), Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[pt]))!=0){
2119         AliError(Form("Projection mismatch :: request[H%sTrkltRCZ%c] found[%s]", mc?"MC":"", ptName[pt], pr0->fH->GetName()));
2120         return kFALSE;
2121       }
2122     }
2123     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
2124   }
2125   TObjArray *arr(NULL);
2126   fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
2127
2128   TH2 *h2(NULL); Int_t jh(0);
2129   for(; ih--; ){
2130     if(!hp[ih].fH) continue;
2131     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2132     arr->AddAt(h2, jh++);
2133   }
2134   // build combined performance plots
2135   // combine up the tree of projections
2136   AliTRDresolutionProjection xlow[2], specY[kNcharge*AliPID::kSPECIES], specPh[kNcharge*AliPID::kSPECIES], specQ[kNcharge*AliPID::kSPECIES];
2137   for(Int_t ich(0); ich<kNcharge; ich++){
2138     // PID dependency - summation over pt
2139     for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2140       /*!dy*/
2141       Int_t idx(ich*AliPID::kSPECIES+isp);
2142       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2143         specY[idx] = (*pr0);
2144         specY[idx].SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2145         specY[idx].fH->SetNameTitle(Form("H%sTrkInY%c%d", mc?"MC":"", chName[ich], isp),
2146                               Form("TrackIn[%s%c]:: #Deltay", AliPID::ParticleLatexName(isp), chSgn[ich]));
2147         for(Int_t ipt(1); ipt<nPt; ipt++){
2148           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2149           specY[idx]+=(*pr1);
2150         }
2151         php.AddLast(&specY[idx]);
2152         if((h2 = specY[idx].Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2153         if((h2 = (TH2*)gDirectory->Get(Form("%sEn", specY[idx].fH->GetName())))) arr->AddAt(h2, jh++);
2154         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specY[idx];
2155       }
2156       /*!dphi*/
2157       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2158         specPh[idx] = (*pr0);
2159         specPh[idx].SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2160         specPh[idx].fH->SetNameTitle(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[ich], isp),
2161                               Form("TrackIn[%s%c]:: #Delta#phi", AliPID::ParticleLatexName(isp), chSgn[ich]));
2162         specPh[idx].SetShowRange(-1.5, 1.5);
2163         for(Int_t ipt(1); ipt<nPt; ipt++){
2164           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2165           specPh[idx]+=(*pr1);
2166         }
2167         php.AddLast(&specPh[idx]);
2168         if((h2 = specPh[idx].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2169         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specPh[idx];
2170       }
2171       /*!dQdl*/
2172       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[0], isp)))){
2173         specQ[idx] = (*pr0);
2174         specQ[idx].SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp), "Sum over pt");
2175         specQ[idx].fH->SetNameTitle(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[ich], isp),
2176                               Form("TrackIn[%s%c]:: dQdl", AliPID::ParticleLatexName(isp), chSgn[ich]));
2177         specQ[idx].SetShowRange(-2.2, -1.75);
2178         specQ[idx].fH->GetZaxis()->SetTitle("dQdl [a.u.]");
2179         for(Int_t ipt(1); ipt<nPt; ipt++){
2180           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2181           specQ[idx]+=(*pr1);
2182         }
2183         php.AddLast(&specQ[idx]);
2184         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2185         specQ[idx].fH->SetName(Form("H%sTrkInQS%c%d", mc?"MC":"", chName[ich], isp));
2186         specQ[idx].SetShowRange(-1.75, -1.25);
2187         if((h2 = specQ[idx].Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2188         if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))) (*pr1)+=specQ[idx];
2189       }
2190     } // end PID loop for pt integration
2191
2192     // pt dependency - summation over PID
2193     for(Int_t ipt(0); ipt<nPt; ipt++){
2194       /*!dy*/
2195       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2196         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2197           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2198           (*pr0)+=(*pr1);
2199         }
2200         pr0->fH->SetNameTitle(Form("H%sTrkInY%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2201                                   Form("TrackIn[%c]:: #Deltay{%s}", chSgn[ich], ptCut[ipt]));
2202         pr0->SetShowRange(-0.3, 0.3);
2203         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2204         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2205       }
2206       /*!dphi*/
2207       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2208         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2209           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2210           (*pr0)+=(*pr1);
2211         }
2212         pr0->fH->SetNameTitle(Form("H%sTrkInPh%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2213                                   Form("TrackIn[%c]:: #Delta#phi{%s}", chSgn[ich], ptCut[ipt]));
2214         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2215         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2216       }
2217       /*!dx*/
2218       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], 0)))){
2219         for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2220           if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[ipt], isp)))) continue;
2221           (*pr0)+=(*pr1);
2222         }
2223         pr0->fH->SetNameTitle(Form("H%sTrkInX%c%c", mc?"MC":"", chName[ich], ptName[ipt]),
2224                                   Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[ipt]));
2225         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2226         if(!ipt){
2227           xlow[ich] = (*pr0);
2228           xlow[ich].SetNameTitle(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5),
2229                                  Form("TrackIn[%c]:: #Deltax{%s}", chSgn[ich], ptCut[0]));
2230           php.AddLast(&xlow[ich]);
2231         }
2232         if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))) (*pr1)+=(*pr0);
2233       }
2234     } // end pt loop for PID integration
2235
2236     /*!dy*/
2237     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2238       pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", chName[ich]),
2239                             Form("TrackIn[%c]:: #Deltay", chSgn[ich]));
2240       pr0->SetShowRange(-0.3, 0.3);
2241       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2242       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2243       if(ich && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2244     }
2245     /*!dy high pt*/
2246     if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[2], 0)))){
2247       if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[ich], ptName[2], 0)))){
2248         (*pr0)+=(*pr1);
2249         pr0->fH->SetNameTitle(Form("H%sTrkInY%c", mc?"MC":"", ptName[2]), Form("TrackIn :: #Deltay{%s}", ptCut[2]));
2250         pr0->SetShowRange(-0.3, 0.3);
2251         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2252       }
2253     }
2254     /*!dphi*/
2255     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2256       pr0->fH->SetNameTitle(Form("H%sTrkInPh%c", mc?"MC":"", chName[ich]),
2257                             Form("TrackIn[%c]:: #Delta#phi", chSgn[ich]));
2258       pr0->SetShowRange(-1., 1.);
2259       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2260       if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2261     }
2262     /*!dx*/
2263     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 0)))){
2264       pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", chName[ich]),
2265                             Form("TrackIn[%c]:: #Deltax", chSgn[ich]));
2266       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2267       if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))) (*pr1)+=(*pr0);
2268     }
2269     /*!dx low pt*/
2270     if(ich && (pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 5)))){
2271       if((pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[ich], ptName[0], 5)))){
2272         (*pr0)+=(*pr1);
2273         pr0->fH->SetNameTitle(Form("H%sTrkInX%c", mc?"MC":"", ptName[0]), Form("TrackIn :: #Deltax{%s}", ptCut[0]));
2274         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2275       }
2276     }
2277   } // end charge loop
2278
2279   for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2280     /*!dy*/
2281     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%d", mc?"MC":"", chName[0], isp)))){
2282       pr0->fH->SetNameTitle(Form("H%sTrkInY%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Deltay", AliPID::ParticleLatexName(isp)));
2283       pr0->SetShowRange(-0.3, 0.3);
2284       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2285       if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2286     }
2287     /*!dphi*/
2288     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%d", mc?"MC":"", chName[0], isp)))){
2289       pr0->fH->SetNameTitle(Form("H%sTrkInPh%d", mc?"MC":"", isp), Form("TrackIn[%s] :: #Delta#phi", AliPID::ParticleLatexName(isp)));
2290       pr0->SetShowRange(-1., 1.);
2291       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2292     }
2293     /*!dQdl*/
2294     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInQ%c%d", mc?"MC":"", chName[0], isp)))){
2295       pr0->fH->SetNameTitle(Form("H%sTrkInQ%d", mc?"MC":"", isp), Form("TrackIn[%s] :: dQdl", AliPID::ParticleLatexName(isp)));
2296       pr0->SetShowRange(-2.2, -1.75);
2297       if((h2 = pr0->Projection2D(kNstat, kNcontours, 2))) arr->AddAt(h2, jh++);
2298       pr0->fH->SetName(Form("H%sTrkInQS%d", mc?"MC":"", isp));
2299       pr0->SetShowRange(-1.7, -1.25);
2300       if((h2 = pr0->Projection2D(kNstat, kNcontours, 0))) arr->AddAt(h2, jh++);
2301     }
2302   } // end PID processing
2303
2304   /*!dy*/
2305   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInY%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2306     pr0->fH->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
2307     pr0->SetShowRange(-0.3, 0.3);
2308     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1, kFALSE))) arr->AddAt(h2, jh++);
2309     if((h2 = (TH2*)gDirectory->Get(Form("%sEn", pr0->fH->GetName())))) arr->AddAt(h2, jh++);
2310   }
2311   /*!dphi*/
2312   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInPh%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2313     pr0->fH->SetNameTitle(Form("H%sTrkInPh", mc?"MC":""), "TrackIn :: #Delta#phi");
2314     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2315   }
2316   /*!dx*/
2317   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInX%c%c%d", mc?"MC":"", chName[0], ptName[0], 0)))){
2318     pr0->fH->SetNameTitle(Form("H%sTrkInX", mc?"MC":""), "TrackIn :: #Deltax");
2319     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2320   }
2321
2322   // Row Cross processing
2323   /*!RC dz*/
2324   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[0])))){
2325     for(Int_t ipt(0); ipt<kNpt; ipt++){
2326       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCZ%c", mc?"MC":"", ptName[ipt])))) continue;
2327       (*pr0)+=(*pr1);
2328     }
2329     pr0->fH->SetNameTitle(Form("H%sTrkInRCZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
2330     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2331   }
2332   /*!RC dy*/
2333   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[0])))){
2334     for(Int_t ipt(0); ipt<kNpt; ipt++){
2335       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCY%c", mc?"MC":"", ptName[ipt])))) continue;
2336       (*pr0)+=(*pr1);
2337     }
2338     pr0->fH->SetNameTitle(Form("H%sTrkInRCY", mc?"MC":""), "TrackIn[RC]:: #Deltay");
2339     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2340   }
2341   /*!RC dphi*/
2342   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[0])))){
2343     for(Int_t ipt(0); ipt<kNpt; ipt++){
2344       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCPh%c", mc?"MC":"", ptName[ipt])))) continue;
2345       (*pr0)+=(*pr1);
2346     }
2347     pr0->fH->SetNameTitle(Form("H%sTrkInRCPh", mc?"MC":""), "TrackIn[RC]:: #Delta#phi");
2348     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2349   }
2350   /*!RC dx*/
2351   if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[0])))){
2352     for(Int_t ipt(0); ipt<kNpt; ipt++){
2353       if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("H%sTrkInRCX%c", mc?"MC":"", ptName[ipt])))) continue;
2354       (*pr0)+=(*pr1);
2355     }
2356     pr0->fH->SetNameTitle(Form("H%sTrkInRCX", mc?"MC":""), "TrackIn[RC]:: #Deltax");
2357     if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2358   }
2359   AliInfo(Form("Done %3d 2D projections.", jh));
2360   return kTRUE;
2361 }
2362
2363
2364 //________________________________________________________
2365 Bool_t AliTRDresolution::MakeProjectionTrack()
2366 {
2367 // Analyse tracklet
2368   const Int_t kNcontours(9);
2369   const Int_t kNstat(30);
2370   Int_t cidx(kMCtrack);
2371   if(fProj && fProj->At(cidx)) return kTRUE;
2372   if(!fContainer){
2373     AliError("Missing data container.");
2374     return kFALSE;
2375   }
2376   THnSparse *H(NULL);
2377   if(!(H = (THnSparse*)fContainer->FindObject("hTRD2MC"))){
2378     AliError("Missing/Wrong data @ hTRD2MC.");
2379     return kFALSE;
2380   }
2381   Int_t ndim(H->GetNdimensions());
2382   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
2383   TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
2384   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
2385   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
2386   if(ndim > kPt) ap = H->GetAxis(kPt);
2387
2388   // build list of projections
2389   const Int_t nsel(AliTRDgeometry::kNlayer*kNpt*AliPID::kSPECIES*7);//, npsel(3);
2390   const Char_t chName[kNcharge] = {'n', 'p'};const Char_t chSgn[kNcharge] = {'-', '+'};
2391   const Char_t ptName[kNpt] = {'l', 'i', 'h'};
2392   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"};
2393   // define rebinning strategy
2394   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
2395   AliTRDresolutionProjection hp[kTrkNproj]; TObjArray php(kTrkNproj);
2396   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
2397   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2398     for(Int_t ipt(0); ipt<kNpt; ipt++){
2399       for(Int_t isp(0); isp<AliPID::kSPECIES; isp++){
2400         for(Int_t ich(0); ich<kNcharge; ich++){
2401           isel++; // new selection
2402           hp[ih].Build(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2403                        Form("Tracks[%s%c]:: #Deltay{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2404                        kEta, kPhi, kYrez, aa);
2405           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2406           php.AddLast(&hp[ih++]); np[isel]++;
2407           hp[ih].Build(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2408                        Form("Tracks[%s%c]:: #Delta#phi{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2409                        kEta, kPhi, kPrez, aa);
2410           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2411           php.AddLast(&hp[ih++]); np[isel]++;
2412           hp[ih].Build(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily),
2413                        Form("Tracks[%s%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", AliPID::ParticleLatexName(isp), chSgn[ich], ptCut[ipt], ily),
2414                        kEta, kPhi, kNdim, aa);
2415           hp[ih].SetShowRange(0.,10.);
2416           hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2417           php.AddLast(&hp[ih++]); np[isel]++;
2418         }
2419       }
2420       isel++; // new selection
2421       hp[ih].Build(Form("HMCTrkZ%c%d", ptName[ipt], ily),
2422                     Form("Tracks[RC]:: #Deltaz{%s} Ly[%d]", ptCut[ipt], ily),
2423                     kEta, kPhi, kZrez, aa);
2424       hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
2425       php.AddLast(&hp[ih++]); np[isel]++;
2426     }
2427   }
2428
2429   Int_t ly(0), ch(0), pt(0), sp(2), rcBin(as?as->FindBin(0.):-1);
2430   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
2431     v = H->GetBinContent(ib, coord);
2432     if(v<1.) continue;
2433     ly = coord[kBC]-1; // layer selection
2434     // charge selection
2435     ch=0; sp=2;// [pi-] track [dafault]
2436     if(rcBin>0){ // debug mode in which species are also saved
2437       sp = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
2438       if(coord[kSpeciesChgRC] > rcBin) ch = 1;        // [+] track
2439       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
2440     }
2441     // pt selection
2442     pt = 0; // low pt [default]
2443     if(ap) pt = coord[kPt]-1;
2444     // global selection
2445     Int_t ioff = ly*kNpt*31+pt*31; ioff+=3*(sp<0?10:(sp*kNcharge+ch));
2446     isel = ly*kNpt*11+pt*11; isel+=sp<0?10:(sp*kNcharge+ch);
2447     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));
2448     for(Int_t jh(0); jh<np[isel]; jh++) ((AliTRDresolutionProjection*)php.At(ioff+jh))->Increment(coord, v);
2449   }
2450   TObjArray *arr(NULL);
2451   fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
2452
2453   TH2 *h2(NULL); Int_t jh(0);
2454   for(; ih--; ){
2455     if(!hp[ih].fH) continue;
2456     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
2457     arr->AddAt(h2, jh++);
2458   }
2459
2460   // combine up the tree of projections
2461   AliTRDresolutionProjection *pr0(NULL), *pr1(NULL);
2462   //Int_t iproj(0), jproj(0);
2463   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
2464     for(Int_t ich(0); ich<kNcharge; ich++){
2465       for(Int_t ipt(0); ipt<kNpt; ipt++){
2466         /*!dy*/
2467         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2468           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2469             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2470             (*pr0)+=(*pr1);
2471           }
2472           AliDebug(2, Form("Rename %s to HMCTrkY%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2473           pr0->fH->SetNameTitle(Form("HMCTrkY%c%c%d", chName[ich], ptName[ipt], ily),
2474                                     Form("Tracks[%c]:: #Deltay{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2475           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2476           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2477         }
2478         /*!dphi*/
2479         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2480           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2481             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2482             (*pr0)+=(*pr1);
2483           }
2484           AliDebug(2, Form("Rename %s to HMCTrkPh%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2485           pr0->fH->SetNameTitle(Form("HMCTrkPh%c%c%d", chName[ich], ptName[ipt], ily),
2486                                     Form("Tracks[%c]:: #Delta#phi{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2487           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2488           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2489         }
2490
2491         /*!dpt/pt*/
2492         if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], 0, ily)))){
2493           for(Int_t isp(1); isp<AliPID::kSPECIES; isp++){
2494             if(!(pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[ipt], isp, ily)))) continue;
2495             (*pr0)+=(*pr1);
2496           }
2497           AliDebug(2, Form("Rename %s to HMCTrkDPt%c%c%d", pr0->fH->GetName(), chName[ich], ptName[ipt], ily));
2498           pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%c%d", chName[ich], ptName[ipt], ily),
2499                                     Form("Tracks[%c]:: #Deltap_{t}/p_{t}{%s} Ly[%d]", chSgn[ich], ptCut[ipt], ily));
2500           if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2501           if(ipt && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2502         }
2503       }
2504       /*!dy*/
2505       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2506         pr0->fH->SetNameTitle(Form("HMCTrkY%c%d", chName[ich], ily),
2507                               Form("Tracks[%c]:: #Deltay Ly[%d]", chSgn[ich], ily));
2508         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2509         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2510       }
2511       /*!dphi*/
2512       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2513         pr0->fH->SetNameTitle(Form("HMCTrkPh%c%d", chName[ich], ily),
2514                               Form("Tracks[%c]:: #Delta#phi Ly[%d]", chSgn[ich], ily));
2515         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2516         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2517       }
2518       /*!dpt/pt*/
2519       if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[ich], ptName[0], 0, ily)))){
2520         pr0->fH->SetNameTitle(Form("HMCTrkDPt%c%d", chName[ich], ily),
2521                               Form("Tracks[%c]:: #Deltap_{t}/p_{t} Ly[%d]", chSgn[ich], ily));
2522         if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2523         if(ich==1 && (pr1 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))) (*pr1)+=(*pr0);
2524       }
2525     }
2526     /*!dy*/
2527     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkY%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2528       pr0->fH->SetNameTitle(Form("HMCTrkY%d", ily), Form("Tracks :: #Deltay Ly[%d]", ily));
2529       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2530     }
2531     /*!dphi*/
2532     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkPh%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2533       pr0->fH->SetNameTitle(Form("HMCTrkPh%d", ily), Form("Tracks :: #Delta#phi Ly[%d]", ily));
2534       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2535     }
2536     /*!dpt/pt*/
2537     if((pr0 = (AliTRDresolutionProjection*)php.FindObject(Form("HMCTrkDPt%c%c%d%d", chName[0], ptName[0], 0, ily)))){
2538       pr0->fH->SetNameTitle(Form("HMCTrkDPt%d", ily), Form("Tracks :: #Deltap_{t}/p_{t} Ly[%d]", ily));
2539       if((h2 = pr0->Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, jh++);
2540     }
2541   }
2542   AliInfo(Form("Done %3d 2D projections.", jh));
2543   return kTRUE;
2544 }
2545
2546 //________________________________________________________
2547 Bool_t AliTRDresolution::PostProcess()
2548 {
2549 // Fit, Project, Combine, Extract values from the containers filled during execution
2550
2551   if (!fContainer) {
2552     AliError("ERROR: list not available");
2553     return kFALSE;
2554   }
2555   if(!fProj){
2556     AliInfo("Building array of projections ...");
2557     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
2558   }
2559
2560   //PROCESS EXPERIMENTAL DISTRIBUTIONS
2561   // Clusters detector
2562   if(HasProcessDetector() && !MakeProjectionDetector()) return kFALSE;
2563   // Clusters residuals
2564   if(HasProcessCluster() && !MakeProjectionCluster()) return kFALSE;
2565   fNRefFigures = 3;
2566   // Tracklet residual/pulls
2567   if(HasProcessTrklt() && !MakeProjectionTracklet()) return kFALSE;
2568   fNRefFigures = 7;
2569   // TRDin residual/pulls
2570   if(HasProcessTrkIn() && !MakeProjectionTrackIn()) return kFALSE;
2571   fNRefFigures = 11;
2572
2573   if(!HasMCdata()) return kTRUE;
2574   //PROCESS MC RESIDUAL DISTRIBUTIONS
2575
2576   // CLUSTER Y RESOLUTION/PULLS
2577   if(!MakeProjectionCluster(kTRUE)) return kFALSE;
2578   fNRefFigures = 17;
2579
2580   // TRACKLET RESOLUTION/PULLS
2581   if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
2582   fNRefFigures = 21;
2583
2584   // TRACK RESOLUTION/PULLS
2585   if(!MakeProjectionTrack()) return kFALSE;
2586   fNRefFigures+=16;
2587
2588   // TRACK TRDin RESOLUTION/PULLS
2589   if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
2590   fNRefFigures+=8;
2591
2592   return kTRUE;
2593 }
2594
2595
2596 //________________________________________________________
2597 void AliTRDresolution::Terminate(Option_t *opt)
2598 {
2599   AliTRDrecoTask::Terminate(opt);
2600   if(HasPostProcess()) PostProcess();
2601 }
2602
2603 //________________________________________________________
2604 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2605 {
2606 // Helper function to avoid duplication of code
2607 // Make first guesses on the fit parameters
2608
2609   // find the intial parameters of the fit !! (thanks George)
2610   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2611   Double_t sum = 0.;
2612   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2613   f->SetParLimits(0, 0., 3.*sum);
2614   f->SetParameter(0, .9*sum);
2615   Double_t rms = h->GetRMS();
2616   f->SetParLimits(1, -rms, rms);
2617   f->SetParameter(1, h->GetMean());
2618
2619   f->SetParLimits(2, 0., 2.*rms);
2620   f->SetParameter(2, rms);
2621   if(f->GetNpar() <= 4) return;
2622
2623   f->SetParLimits(3, 0., sum);
2624   f->SetParameter(3, .1*sum);
2625
2626   f->SetParLimits(4, -.3, .3);
2627   f->SetParameter(4, 0.);
2628
2629   f->SetParLimits(5, 0., 1.e2);
2630   f->SetParameter(5, 2.e-1);
2631 }
2632
2633 //________________________________________________________
2634 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
2635 {
2636 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2637 //  - y reziduals/pulls
2638
2639   TObjArray *arr = new TObjArray(2);
2640   arr->SetName(name); arr->SetOwner();
2641   TH1 *h(NULL); char hname[100], htitle[300];
2642
2643   // tracklet resolution/pull in y direction
2644   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
2645   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
2646   Float_t rr = range<0.?fDyRange:range;
2647   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2648     Int_t nybins=50;
2649     if(expand) nybins*=2;
2650     h = new TH3S(hname, htitle,
2651                  48, -.48, .48,            // phi
2652                  60, -rr, rr,              // dy
2653                  nybins, -0.5, nybins-0.5);// segment
2654   } else h->Reset();
2655   arr->AddAt(h, 0);
2656   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
2657   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
2658   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2659     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2660   } else h->Reset();
2661   arr->AddAt(h, 1);
2662
2663   return arr;
2664 }
2665
2666 //________________________________________________________
2667 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2668 {
2669 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2670 //  - y reziduals/pulls
2671 //  - z reziduals/pulls
2672 //  - phi reziduals
2673   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
2674   arr->Expand(5);
2675   TH1 *h(NULL); char hname[100], htitle[300];
2676
2677   // tracklet resolution/pull in z direction
2678   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
2679   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
2680   if(!(h = (TH2S*)gROOT->FindObject(hname))){
2681     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
2682   } else h->Reset();
2683   arr->AddAt(h, 2);
2684   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
2685   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
2686   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2687     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2688     h->GetZaxis()->SetBinLabel(1, "no RC");
2689     h->GetZaxis()->SetBinLabel(2, "RC");
2690   } else h->Reset();
2691   arr->AddAt(h, 3);
2692
2693   // tracklet to track phi resolution
2694   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
2695   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
2696   Int_t nsgms=540;
2697   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2698     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
2699   } else h->Reset();
2700   arr->AddAt(h, 4);
2701
2702   return arr;
2703 }
2704
2705 //________________________________________________________
2706 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2707 {
2708 // Build performance histograms for AliExternalTrackParam.vs MC
2709 //  - y resolution/pulls
2710 //  - z resolution/pulls
2711 //  - phi resolution, snp pulls
2712 //  - theta resolution, tgl pulls
2713 //  - pt resolution, 1/pt pulls, p resolution
2714
2715   TObjArray *arr = BuildMonitorContainerTracklet(name);
2716   arr->Expand(11);
2717   TH1 *h(NULL); char hname[100], htitle[300];
2718   //TAxis *ax(NULL);
2719
2720   // snp pulls
2721   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
2722   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
2723   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2724     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2725   } else h->Reset();
2726   arr->AddAt(h, 5);
2727
2728   // theta resolution
2729   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
2730   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2731   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2732     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2733   } else h->Reset();
2734   arr->AddAt(h, 6);
2735   // tgl pulls
2736   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
2737   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
2738   if(!(h = (TH2I*)gROOT->FindObject(hname))){
2739     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2740   } else h->Reset();
2741   arr->AddAt(h, 7);
2742
2743   const Int_t kNdpt(150);
2744   const Int_t kNspc = 2*AliPID::kSPECIES+1;
2745   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
2746   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2747   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2748   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
2749   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
2750
2751   // Pt resolution
2752   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
2753   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2754   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2755     h = new TH3S(hname, htitle,
2756                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2757     //ax = h->GetZaxis();
2758     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2759   } else h->Reset();
2760   arr->AddAt(h, 8);
2761   // 1/Pt pulls
2762   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
2763   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);
2764   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2765     h = new TH3S(hname, htitle,
2766                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2767     //ax = h->GetZaxis();
2768     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2769   } else h->Reset();
2770   arr->AddAt(h, 9);
2771   // P resolution
2772   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
2773   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2774   if(!(h = (TH3S*)gROOT->FindObject(hname))){
2775     h = new TH3S(hname, htitle,
2776                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2777     //ax = h->GetZaxis();
2778     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2779   } else h->Reset();
2780   arr->AddAt(h, 10);
2781
2782   return arr;
2783 }
2784
2785
2786 //________________________________________________________
2787 TObjArray* AliTRDresolution::Histos()
2788 {
2789   //
2790   // Define histograms
2791   //
2792
2793   if(fContainer) return fContainer;
2794
2795   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
2796   THnSparse *H(NULL);
2797   const Int_t nhn(100); Char_t hn[nhn]; TString st;
2798
2799   //++++++++++++++++++++++
2800   // cluster to detector
2801   snprintf(hn, nhn, "h%s", fgPerformanceName[kDetector]);
2802   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2803     const Int_t mdim(6);
2804     const Char_t *cldTitle[mdim] = {"layer", fgkTitle[kPhi], "pad row", "centrality", "q [a.u.]", "no. of pads"/*, "tb [10 Hz]"*/};
2805     const Int_t cldNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], 76, AliTRDeventInfo::kCentralityClasses, 50, kNpads};
2806     const Double_t cldMin[mdim]  = {-0.5, fgkMin[kPhi], -0.5, -0.5, 0., 0.5},
2807                    cldMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], 75.5, AliTRDeventInfo::kCentralityClasses - 0.5, 1200., kNpads+.5};
2808     st = "cluster proprieties;";
2809     // define minimum info to be saved in non debug mode
2810     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimDet);
2811     for(Int_t idim(0); idim<ndim; idim++){ st += cldTitle[idim]; st+=";";}
2812     H = new THnSparseI(hn, st.Data(), ndim, cldNbins, cldMin, cldMax);
2813   } else H->Reset();
2814   fContainer->AddAt(H, kDetector);
2815   //++++++++++++++++++++++
2816   // cluster to tracklet residuals/pulls
2817   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
2818   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2819     const Int_t mdim(10);
2820     const Char_t *clTitle[mdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]", "centrality", "no. of pads"};
2821     const Int_t clNbins[mdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 30, 30, 15, AliTRDeventInfo::kCentralityClasses, kNpads};
2822     const Double_t clMin[mdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45, -0.5, 0.5},
2823                    clMax[mdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45, AliTRDeventInfo::kCentralityClasses - 0.5, kNpads+.5};
2824     st = "cluster spatial&charge resolution;";
2825     // define minimum info to be saved in non debug mode
2826     Int_t ndim=DebugLevel()>=1?mdim:Int_t(kNdimCl);
2827     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2828     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2829   } else H->Reset();
2830   fContainer->AddAt(H, kCluster);
2831   //++++++++++++++++++++++
2832   // tracklet to TRD track
2833   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
2834   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2835     const Int_t mdim(kNdim+8);
2836     Char_t *trTitle[mdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2837     Int_t trNbins[mdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2838     Double_t trMin[mdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2839     Double_t trMax[mdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2840     // set specific fields
2841     trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
2842     trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
2843     trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
2844     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2845     Int_t jdim(kNdim);
2846     trTitle[jdim]=StrDup("centrality"); trNbins[jdim] = AliTRDeventInfo::kCentralityClasses; trMin[jdim] = -.5; trMax[jdim] = AliTRDeventInfo::kCentralityClasses - 0.5;
2847     jdim++; trTitle[jdim]=StrDup("occupancy [%]"); trNbins[jdim] = 12; trMin[jdim] = 25.; trMax[jdim] = 85.;
2848 //    jdim++; trTitle[jdim]=StrDup("gap [cm]"); trNbins[jdim] = 80; trMin[jdim] = 0.1; trMax[jdim] = 2.1;
2849 //     jdim++; trTitle[jdim]=StrDup("size_{0} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2850 //     jdim++; trTitle[jdim]=StrDup("pos_{0} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2851 //     jdim++; trTitle[jdim]=StrDup("size_{1} [cm]"); trNbins[jdim] = 16; trMin[jdim] = 0.15; trMax[jdim] = 1.75;
2852 //     jdim++; trTitle[jdim]=StrDup("pos_{1} [cm]"); trNbins[jdim] = 10; trMin[jdim] = 0.; trMax[jdim] = 3.5;
2853
2854     st = "tracklet spatial&charge resolution;";
2855     // define minimum info to be saved in non debug mode
2856     Int_t ndim=DebugLevel()>=1?(jdim+1):kNdimTrklt;
2857     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2858     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2859   } else H->Reset();
2860   fContainer->AddAt(H, kTracklet);
2861   //++++++++++++++++++++++
2862   // tracklet to TRDin
2863   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
2864   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2865     // set specific fields
2866     const Int_t mdim(kNdim+3);
2867     Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
2868     Int_t trinNbins[mdim];   memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
2869     Double_t trinMin[mdim];  memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
2870     Double_t trinMax[mdim];  memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
2871     trinNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1; trinMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trinMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2872     trinTitle[kNdim]=StrDup("p [GeV/c]"); trinNbins[kNdim] = 24; trinMin[kNdim] = -0.5; trinMax[kNdim] = 23.5;
2873     trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
2874     trinTitle[kNdim+2]=StrDup("Fill Bunch"); trinNbins[kNdim+2]=3500; trinMin[kNdim+2]=-0.5; trinMax[kNdim+2]=3499.5;
2875     st = "r-#phi/z/angular residuals @ TRD entry;";
2876     // define minimum info to be saved in non debug mode
2877     Int_t ndim=DebugLevel()>=1?mdim:kNdim;//kNdimTrkIn;
2878     for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
2879     H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
2880   } else H->Reset();
2881   fContainer->AddAt(H, kTrackIn);
2882   // tracklet to TRDout
2883 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2884
2885
2886   // Resolution histos
2887   if(!HasMCdata()) return fContainer;
2888
2889   //++++++++++++++++++++++
2890   // cluster to TrackRef residuals/pulls
2891   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
2892   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2893     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
2894     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, Int_t(kNcharge)*AliPID::kSPECIES+1, 15};
2895     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, -AliPID::kSPECIES-0.5, -45},
2896                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, AliPID::kSPECIES+0.5, 45};
2897     st = "MC cluster spatial resolution;";
2898     // define minimum info to be saved in non debug mode
2899     Int_t ndim=DebugLevel()>=1?kNdim:4;
2900     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
2901     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
2902   } else H->Reset();
2903   fContainer->AddAt(H, kMCcluster);
2904   //++++++++++++++++++++++
2905   // tracklet to TrackRef
2906   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
2907   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2908     Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2909     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2910     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2911     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2912     // set specific fields
2913     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2914     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2915     trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
2916     trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
2917     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2918
2919     st = "MC tracklet spatial resolution;";
2920     // define minimum info to be saved in non debug mode
2921     Int_t ndim=DebugLevel()>=1?kNdim:4;
2922     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2923     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2924   } else H->Reset();
2925   fContainer->AddAt(H, kMCtracklet);
2926   //++++++++++++++++++++++
2927   // TRDin to TrackRef
2928   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
2929   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2930     st = "MC r-#phi/z/angular residuals @ TRD entry;";
2931     // set specific fields
2932     Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2933     Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2934     Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2935     trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
2936     trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
2937     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2938     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2939     // define minimum info to be saved in non debug mode
2940     Int_t ndim=DebugLevel()>=1?kNdim:7;
2941     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
2942     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2943   } else H->Reset();
2944   fContainer->AddAt(H, kMCtrackIn);
2945   //++++++++++++++++++++++
2946   // track to TrackRef
2947   snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
2948   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
2949     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
2950     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
2951     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
2952     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
2953     // set specific fields
2954     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
2955     trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
2956     trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
2957     trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
2958     trNbins[kSpeciesChgRC] = Int_t(kNcharge)*AliPID::kSPECIES+1;trMin[kSpeciesChgRC] = -AliPID::kSPECIES-0.5; trMax[kSpeciesChgRC] = AliPID::kSPECIES+0.5;
2959     trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
2960
2961     st = "MC track spatial&p_{t} resolution;";
2962     // define minimum info to be saved in non debug mode
2963     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
2964     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
2965     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
2966   } else H->Reset();
2967   fContainer->AddAt(H, kMCtrack);
2968
2969   return fContainer;
2970 }
2971
2972 //____________________________________________________________________
2973 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2974 {
2975 //
2976 // Fit track with a staight line using the "np" clusters stored in the array "points".
2977 // The following particularities are stored in the clusters from points:
2978 //   1. pad tilt as cluster charge
2979 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
2980 // The parameters of the straight line fit are stored in the array "param" in the following order :
2981 //     param[0] - x0 reference radial position
2982 //     param[1] - y0 reference r-phi position @ x0
2983 //     param[2] - z0 reference z position @ x0
2984 //     param[3] - slope dy/dx
2985 //     param[4] - slope dz/dx
2986 //
2987 // Attention :
2988 // Function should be used to refit tracks for B=0T
2989 //
2990
2991   if(np<40){
2992     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2993     return kFALSE;
2994   }
2995   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2996
2997   Double_t x0(0.);
2998   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2999   x0/=Float_t(np);
3000
3001   Double_t x, y, z, dx, tilt(0.);
3002   for(Int_t ip(0); ip<np; ip++){
3003     x = points[ip].GetX(); z = points[ip].GetZ();
3004     dx = x - x0;
3005     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3006   }
3007   if(zfitter.Eval() != 0) return kFALSE;
3008
3009   Double_t z0    = zfitter.GetParameter(0);
3010   Double_t dzdx  = zfitter.GetParameter(1);
3011   for(Int_t ip(0); ip<np; ip++){
3012     if(points[ip].GetClusterType()) continue;
3013     x    = points[ip].GetX();
3014     dx   = x - x0;
3015     y    = points[ip].GetY();
3016     z    = points[ip].GetZ();
3017     tilt = points[ip].GetCharge();
3018     y -= tilt*(-dzdx*dx + z - z0);
3019     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3020     yfitter.AddPoint(&dx, y, 1.);
3021   }
3022   if(yfitter.Eval() != 0) return kFALSE;
3023   Double_t y0   = yfitter.GetParameter(0);
3024   Double_t dydx = yfitter.GetParameter(1);
3025
3026   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3027   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);
3028   return kTRUE;
3029 }
3030
3031 //____________________________________________________________________
3032 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3033 {
3034 //
3035 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3036 // See function FitTrack for the data stored in the "clusters" array
3037
3038 // The parameters of the straight line fit are stored in the array "param" in the following order :
3039 //     par[0] - x0 reference radial position
3040 //     par[1] - y0 reference r-phi position @ x0
3041 //     par[2] - slope dy/dx
3042 //
3043 // Attention :
3044 // Function should be used to refit tracks for B=0T
3045 //
3046
3047   TLinearFitter yfitter(2, "pol1");
3048
3049   // grep data for tracklet
3050   Double_t x0(0.), x[60], y[60], dy[60];
3051   Int_t nly(0);
3052   for(Int_t ip(0); ip<np; ip++){
3053     if(points[ip].GetClusterType()) continue;
3054     if(points[ip].GetVolumeID() != ly) continue;
3055     Float_t xt(points[ip].GetX())
3056            ,yt(param[1] + param[3] * (xt - param[0]));
3057     x[nly] = xt;
3058     y[nly] = points[ip].GetY();
3059     dy[nly]= y[nly]-yt;
3060     x0    += xt;
3061     nly++;
3062   }
3063   if(nly<10){
3064     if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
3065     return kFALSE;
3066   }
3067   // set radial reference for fit
3068   x0 /= Float_t(nly);
3069
3070   // find tracklet core
3071   Double_t mean(0.), sig(1.e3);
3072   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3073
3074   // simple cluster error parameterization
3075   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3076
3077   // fit tracklet core
3078   for(Int_t jly(0); jly<nly; jly++){
3079     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3080     Double_t dx(x[jly]-x0);
3081     yfitter.AddPoint(&dx, y[jly], 1.);
3082   }
3083   if(yfitter.Eval() != 0) return kFALSE;
3084   par[0] = x0;
3085   par[1] = yfitter.GetParameter(0);
3086   par[2] = yfitter.GetParameter(1);
3087   return kTRUE;
3088 }
3089
3090 //____________________________________________________________________
3091 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
3092 {
3093 //
3094 // Global selection mechanism of tracksbased on cluster to fit residuals
3095 // The parameters are the same as used ni function FitTrack().
3096
3097   const Float_t kS(0.6), kM(0.2);
3098   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3099   Float_t dy, dz, s, m;
3100   for(Int_t ip(0); ip<np; ip++){
3101     if(points[ip].GetClusterType()) continue;
3102     Float_t x0(points[ip].GetX())
3103           ,y0(param[1] + param[3] * (x0 - param[0]))
3104           ,z0(param[2] + param[4] * (x0 - param[0]));
3105     dy=points[ip].GetY() - y0; h.Fill(dy);
3106     dz=points[ip].GetZ() - z0;
3107   }
3108   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3109   fg.SetParameter(1, 0.);
3110   fg.SetParameter(2, 2.e-2);
3111   h.Fit(&fg, "QN");
3112   m=fg.GetParameter(1); s=fg.GetParameter(2);
3113   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3114   return kTRUE;
3115 }
3116
3117 //________________________________________________________
3118 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3119 {
3120   //
3121   // Get the most probable value and the full width half mean
3122   // of a Landau distribution
3123   //
3124
3125   const Float_t dx = 1.;
3126   mpv = f->GetParameter(1);
3127   Float_t fx, max = f->Eval(mpv);
3128
3129   xm = mpv - dx;
3130   while((fx = f->Eval(xm))>.5*max){
3131     if(fx>max){
3132       max = fx;
3133       mpv = xm;
3134     }
3135     xm -= dx;
3136   }
3137
3138   xM += 2*(mpv - xm);
3139   while((fx = f->Eval(xM))>.5*max) xM += dx;
3140 }
3141
3142
3143 // #include "TFile.h"
3144 // //________________________________________________________
3145 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
3146 // {
3147 //   if(!file){
3148 //     AliWarning("Use cluster position as in reconstruction.");
3149 //     SetLoadCorrection();
3150 //     return kTRUE;
3151 //   }
3152 //   TDirectory *cwd(gDirectory);
3153 //   TString fileList;
3154 //   FILE *filePtr = fopen(file, "rt");
3155 //   if(!filePtr){
3156 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
3157 //     SetLoadCorrection();
3158 //     return kFALSE;
3159 //   }
3160 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
3161 //   while(fileList.Gets(filePtr)){
3162 //     if(!TFile::Open(fileList.Data())) {
3163 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
3164 //       continue;
3165 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
3166 //
3167 //     TTree *tSys = (TTree*)gFile->Get("tSys");
3168 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
3169 //     tSys->Draw("det:t>>h2", "dx", "goff");
3170 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3171 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
3172 //     }
3173 //     h2->SetDirectory(cwd);
3174 //     gFile->Close();
3175 //   }
3176 //   cwd->cd();
3177 //
3178 //   if(AliLog::GetDebugLevel("PWGPP", "AliTRDresolution")>=2){
3179 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3180 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
3181 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
3182 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
3183 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
3184 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
3185 //       printf("%03d|", idet);
3186 //       fprintf(fout, "    {");
3187 //       for(Int_t it(0); it<30; it++){
3188 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
3189 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
3190 //       }
3191 //       printf("\n");
3192 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
3193 //     }
3194 //     fprintf(fout, "  };\n");
3195 //   }
3196 //   SetLoadCorrection();
3197 //   return kTRUE;
3198 // }
3199
3200 //________________________________________________________
3201 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
3202   :TNamed()
3203   ,fH(NULL)
3204   ,fNrebin(0)
3205   ,fRebinX(NULL)
3206   ,fRebinY(NULL)
3207 {
3208   // constructor
3209   memset(fAx, 0, 3*sizeof(Int_t));
3210   memset(fRange, 0, 4*sizeof(Float_t));
3211 }
3212
3213 //________________________________________________________
3214 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
3215 {
3216   // destructor
3217   if(fH) delete fH;
3218 }
3219
3220 //________________________________________________________
3221 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
3222 {
3223 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
3224   if(!aa[ix] || !aa[iy] || !aa[iz]) return;
3225   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
3226   // check ax definiton to protect against older versions of the data
3227   if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){
3228     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));
3229     return;
3230   }
3231   if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){
3232     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));
3233     return;
3234   }
3235   if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){
3236     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));
3237     return;
3238   }
3239   SetNameTitle(n,t);
3240   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
3241     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
3242     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
3243     az->GetNbins(), az->GetXmin(), az->GetXmax());
3244   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
3245   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
3246   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,
3247     ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
3248     ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
3249     az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));
3250 }
3251
3252 //________________________________________________________
3253 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator=(const AliTRDresolutionProjection& rhs)
3254 {
3255 // copy projections
3256   if(this == &rhs) return *this;
3257
3258   TNamed::operator=(rhs);
3259   if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}
3260   if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);
3261   memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));
3262   memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));
3263   if(fH) delete fH;
3264   if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));
3265   return *this;
3266 }
3267
3268 //________________________________________________________
3269 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator+=(const AliTRDresolutionProjection& other)
3270 {
3271 // increment projections
3272   if(!fH || !other.fH) return *this;
3273   AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));
3274   fH->Add(other.fH);
3275   return *this;
3276 }
3277
3278 //________________________________________________________
3279 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
3280 {
3281 // increment bin with value "v" pointed by general coord in "bin"
3282   if(!fH) return;
3283   AliDebug(4, Form("  %s[%2d]", fH->GetName(), Int_t(v)));
3284   fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));
3285 }
3286
3287 //________________________________________________________
3288 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)
3289 {
3290 // build the 2D projection and adjust binning
3291
3292   const Char_t *title[] = {"Mean", "#mu", "MPV"};
3293   if(!fH) return NULL;
3294   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
3295   TH2D *h2s(NULL), *hyx(NULL);
3296   if(!(h2s = (TH2D*)fH->Project3D("yx"))) return NULL;
3297   // save a copy of the original distribution
3298   if(!del){
3299     hyx = (TH2D*)h2s->Clone();
3300     hyx->SetName(Form("%sEn", fH->GetName()));
3301   }
3302   Int_t irebin(0), dxBin(1), dyBin(1);
3303   while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
3304     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
3305     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
3306     irebin++;
3307   }
3308   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
3309   delete h2s;
3310
3311   // start projection
3312   TH1 *h(NULL); Int_t n(0);
3313   Float_t dz=(fRange[1]-fRange[1])/ncol;
3314   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
3315   Int_t nt(tokenTitle->GetEntriesFast());
3316   TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
3317             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():""),
3318             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
3319   h2->SetContour(ncol);
3320   h2->GetZaxis()->CenterTitle();
3321   h2->GetZaxis()->SetTitleOffset(1.4);
3322   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
3323   AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));
3324   for(Int_t iy(0); iy<ny; iy++){
3325     for(Int_t ix(0); ix<nx; ix++){
3326       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
3327       Int_t ne((Int_t)h->Integral());
3328       if(ne<nstat/2){
3329         h2->SetBinContent(ix+1, iy+1, -999);
3330         h2->SetBinError(ix+1, iy+1, 1.);
3331         n++;
3332       }else{
3333         Float_t v(h->GetMean()), ve(h->GetRMS());
3334         if(mid==1){
3335           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
3336           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
3337           h->Fit(&fg, "WQ");
3338           v = fg.GetParameter(1); ve = fg.GetParameter(2);
3339         } else if (mid==2) {
3340           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
3341           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
3342           h->Fit(&fl, "WQ");
3343           v = fl.GetMaximumX(); ve = fl.GetParameter(2);
3344 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
3345           fgle.SetParameter(0, fl.GetParameter(0));
3346           fgle.SetParameter(1, fl.GetParameter(1));
3347           fgle.SetParameter(2, fl.GetParameter(2));
3348           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
3349           h->Fit(&fgle, "WQ");
3350           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
3351         }
3352         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
3353         else h2->SetBinContent(ix+1, iy+1, v);
3354         h2->SetBinError(ix+1, iy+1, ve);
3355       }
3356     }
3357   }
3358   if(h) delete h;
3359   if(n==nx*ny){delete h2; h2=NULL;} // clean empty projections
3360   return h2;
3361 }
3362
3363 //________________________________________________________
3364 void AliTRDresolution::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)
3365 {
3366 // Normalize histo content to the mean value in the range specified by bin ranges
3367 // [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.
3368 // Optionally a threshold "thr" can be specified to disregard entries with no meaning
3369
3370   Float_t s = 0., c=0.; Int_t is(0);
3371   for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){
3372     for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){
3373       if((c = h2->GetBinContent(ix, iy))<thr) continue;
3374       s += c; is++;
3375     }
3376   }
3377   s/= is?is:1;
3378   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){
3379     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){
3380       if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-100);
3381       else h2->SetBinContent(ix, iy, 100.*(c/s-1.));
3382     }
3383   }
3384 }
3385
3386 //________________________________________________________
3387 void AliTRDresolution::SetProcesses(Bool_t det, Bool_t cl, Bool_t trklt, Bool_t trkin)
3388 {
3389 // steer processes
3390   if(det) SETBIT(fSteer, kDetector); else CLRBIT(fSteer, kDetector);
3391   if(cl)  SETBIT(fSteer, kCluster); else CLRBIT(fSteer, kCluster);
3392   if(trklt) SETBIT(fSteer, kTracklet); else CLRBIT(fSteer, kTracklet);
3393   if(trkin) SETBIT(fSteer, kTrackIn); else CLRBIT(fSteer, kTrackIn);
3394 }
3395
3396 //________________________________________________________
3397 void AliTRDresolution::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)
3398 {
3399 // Set range on Z axis such to avoid outliers
3400
3401   Float_t c(0.), dz(1.e-3*(max-min));
3402   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){
3403     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){
3404       if((c = h2->GetBinContent(ix, iy))<thr) continue;
3405       if(c<=min) h2->SetBinContent(ix, iy, min+dz);
3406     }
3407   }
3408   h2->GetZaxis()->SetRangeUser(min, max);
3409 }
3410
3411 //________________________________________________________
3412 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
3413 {
3414 // define rebinning strategy for this projection
3415   fNrebin = n;
3416   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
3417   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));
3418 }
3419
3420