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