PWG1/TRD/AliTRDresolution.cxx

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