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