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
  92 #include "info/AliTRDclusterInfo.h"
  93
  94 ClassImp(AliTRDresolution)
  95 //ClassImp(AliTRDresolution::AliTRDresolutionProjection)
  96
  97 Int_t const   AliTRDresolution::fgkNbins[kNdim] = {
  98   Int_t(kNbunchCross)/*bc*/,
  99   180/*phi*/,
 100   50/*eta*/,
 101   50/*dy*/,
 102   40/*dphi*/,
 103   50/*dz*/,
 104   Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
 105   kNpt/*pt*/
 106 };  //! no of bins/projection
 107 Double_t const AliTRDresolution::fgkMin[kNdim] = {
 108   -0.5,
 109   -TMath::Pi(),
 110   -1.,
 111   -1.5,
 112   -10.,
 113   -2.5,
 114   -AliPID::kSPECIES-0.5,
 115   -0.5
 116 };    //! low limits for projections
 117 Double_t const AliTRDresolution::fgkMax[kNdim] = {
 118   Int_t(kNbunchCross)-0.5,
 119   TMath::Pi(),
 120   1.,
 121   1.5,
 122   10.,
 123   2.5,
 124   AliPID::kSPECIES+0.5,
 125   kNpt-0.5
 126 };    //! high limits for projections
 127 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
 128   "bunch cross",
 129   "#phi [rad]",
 130   "#eta",
 131   "#Deltay [cm]",
 132   "#Delta#phi [deg]",
 133   "#Deltaz [cm]",
 134   "chg*spec*rc",
 135   "bin_p_{t}"
 136 };  //! title of projection
 137
 138 Int_t const AliTRDresolution::fgkNproj[kNclasses] = {
 139   48, 72, 8, 5,
 140   2, 5, 11, 11, 11
 141 };
 142 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
 143     "Cluster2Track"
 144     ,"Tracklet2Track"
 145     ,"Tracklet2TRDin"
 146     ,"Tracklet2TRDout"
 147     ,"Cluster2MC"
 148     ,"Tracklet2MC"
 149     ,"TRDin2MC"
 150     ,"TRDout2MC"
 151     ,"TRD2MC"
 152 };
 153 Float_t AliTRDresolution::fgPtBin[kNpt+1];
 154
 155 //________________________________________________________
 156 AliTRDresolution::AliTRDresolution()
 157   :AliTRDrecoTask()
 158   ,fIdxPlot(0)
 159   ,fIdxFrame(0)
 160   ,fPtThreshold(1.)
 161   ,fDyRange(0.75)
 162   ,fProj(NULL)
 163   ,fDBPDG(NULL)
 164   ,fCl(NULL)
 165   ,fMCcl(NULL)
 166 {
 167   //
 168   // Default constructor
 169   //
 170   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
 171   MakePtSegmentation();
 172 }
 173
 174 //________________________________________________________
 175 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
 176   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
 177   ,fIdxPlot(0)
 178   ,fIdxFrame(0)
 179   ,fPtThreshold(1.)
 180   ,fDyRange(0.75)
 181   ,fProj(NULL)
 182   ,fDBPDG(NULL)
 183   ,fCl(NULL)
 184   ,fMCcl(NULL)
 185 {
 186   //
 187   // Default constructor
 188   //
 189
 190   InitFunctorList();
 191   MakePtSegmentation();
 192   if(xchange){
 193     SetUseExchangeContainers();
 194     DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
 195     DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
 196   }
 197 }
 198
 199 //________________________________________________________
 200 AliTRDresolution::~AliTRDresolution()
 201 {
 202   //
 203   // Destructor
 204   //
 205
 206   if(fProj){fProj->Delete(); delete fProj;}
 207   if(fCl){fCl->Delete(); delete fCl;}
 208   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
 209 }
 210
 211
 212 //________________________________________________________
 213 void AliTRDresolution::UserCreateOutputObjects()
 214 {
 215   // spatial resolution
 216
 217   AliTRDrecoTask::UserCreateOutputObjects();
 218   if(UseExchangeContainers()) InitExchangeContainers();
 219 }
 220
 221 //________________________________________________________
 222 void AliTRDresolution::InitExchangeContainers()
 223 {
 224 // Init containers for subsequent tasks (AliTRDclusterResolution)
 225
 226   fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
 227   fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
 228   PostData(kClToTrk, fCl);
 229   PostData(kClToMC, fMCcl);
 230 }
 231
 232 //________________________________________________________
 233 void AliTRDresolution::UserExec(Option_t *opt)
 234 {
 235   //
 236   // Execution part
 237   //
 238
 239   if(fCl) fCl->Delete();
 240   if(fMCcl) fMCcl->Delete();
 241   AliTRDrecoTask::UserExec(opt);
 242 }
 243
 244 //________________________________________________________
 245 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
 246 {
 247 // Helper function to calculate pulls in the yz plane
 248 // using proper tilt rotation
 249 // Uses functionality defined by AliTRDseedV1.
 250
 251   return kTRUE;
 252 /*
 253   Double_t t2(tilt*tilt);
 254   // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
 255
 256   // rotate along pad
 257   Double_t cc[3];
 258   cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
 259   cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
 260   cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
 261   // do sqrt
 262   Double_t sqr[3]={0., 0., 0.};
 263   if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
 264   Double_t invsqr[3]={0., 0., 0.};
 265   if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
 266   Double_t tmp(dyz[0]);
 267   dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
 268   dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
 269   return kTRUE;
 270 */
 271 }
 272
 273 //________________________________________________________
 274 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
 275 {
 276   //
 277   // Plot the cluster distributions
 278   //
 279
 280   if(track) fkTrack = track;
 281   if(!fkTrack){
 282     AliDebug(4, "No Track defined.");
 283     return NULL;
 284   }
 285   if(TMath::Abs(fkESD->GetTOFbc()) > 1){
 286     AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
 287     return NULL;
 288   }
 289   if(fPt<fPtThreshold){
 290     AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
 291     return NULL;
 292   }
 293   THnSparse *H(NULL);
 294   if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
 295     AliWarning("No output container defined.");
 296     return NULL;
 297   }
 298
 299   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 300   Double_t val[kNdim]; //Float_t exb, vd, t0, s2, dl, dt;
 301   TObjArray     *clInfoArr(NULL);
 302   AliTRDseedV1  *fTracklet(NULL);
 303   AliTRDcluster *c(NULL), *cc(NULL);
 304   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 305     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 306     if(!fTracklet->IsOK()) continue;
 307     //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
 308     val[kBC]  = ily;
 309     val[kPhi] = fPhi;
 310     val[kEta] = fEta;
 311     val[kPt]  = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
 312     Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
 313     Int_t row0(-1);
 314     Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
 315     fTracklet->ResetClusterIter(kTRUE);
 316     while((c = fTracklet->NextCluster())){
 317       Float_t xc(c->GetX()),
 318               q(c->GetQ());
 319       Int_t tb(c->GetLocalTimeBin());
 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);
 325       if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += cc->GetQ(); ic++;}
 326       if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += 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 //  return NULL;
 411   Double_t val[kNdim+1];
 412   AliTRDseedV1 *fTracklet(NULL);
 413   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
 414     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
 415     if(!fTracklet->IsOK()) continue;
 416     val [kBC] = il;
 417     val[kPhi] = fPhi;
 418     val[kEta] = fEta;
 419     val[kSpeciesChgRC]= fSpecies;
 420     val[kPt]  = GetPtBin(fTracklet->GetMomentum());
 421     Double_t dyt(fTracklet->GetYref(0) - fTracklet->GetYfit(0)),
 422              dzt(fTracklet->GetZref(0) - fTracklet->GetZfit(0)),
 423              dydx(fTracklet->GetYfit(1)),
 424              tilt(fTracklet->GetTilt());
 425     // correct for tilt rotation
 426     val[kYrez] = dyt - dzt*tilt;
 427     val[kZrez] = dzt + dyt*tilt;
 428     dydx+= tilt*fTracklet->GetZref(1);
 429     val[kPrez] = TMath::ATan((fTracklet->GetYref(1) - dydx)/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
 430     if(fTracklet->IsRowCross()){
 431       val[kSpeciesChgRC]= 0.;
 432 //      val[kPrez] = fkTrack->Charge(); // may be better defined
 433     }/* else {
 434       Float_t exb, vd, t0, s2, dl, dt;
 435       fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
 436       val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
 437     }*/
 438     val[kNdim] = fTracklet->GetdQdl();
 439     if(DebugLevel()>=1) H->Fill(val);
 440
 441 //     // compute covariance matrix
 442 //     fTracklet->GetCovAt(x, cov);
 443 //     fTracklet->GetCovRef(covR);
 444 //     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
 445 //     Double_t dyz[2]= {dy[1], dz[1]};
 446 //     Pulls(dyz, cov, tilt);
 447 //     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
 448 //     ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
 449
 450     if(DebugLevel()>=3){
 451       Bool_t rc(fTracklet->IsRowCross());
 452       UChar_t err(fTracklet->GetErrorMsg());
 453       Double_t x(fTracklet->GetX()),
 454                pt(fTracklet->GetPt()),
 455                yt(fTracklet->GetYref(0)),
 456                zt(fTracklet->GetZref(0)),
 457                phi(fTracklet->GetYref(1)),
 458                tht(fTracklet->GetZref(1));
 459       Int_t ncl(fTracklet->GetN()),
 460             det(fTracklet->GetDetector());
 461       (*DebugStream()) << "tracklet"
 462         <<"pt="  << pt
 463         <<"x="   << x
 464         <<"yt="  << yt
 465         <<"zt="  << zt
 466         <<"phi=" << phi
 467         <<"tht=" << tht
 468         <<"det=" << det
 469         <<"n="   << ncl
 470         <<"dy0=" << dyt
 471         <<"dz0=" << dzt
 472         <<"dy="  << val[kYrez]
 473         <<"dz="  << val[kZrez]
 474         <<"dphi="<< val[kPrez]
 475         <<"dQ  ="<< val[kNdim]
 476         <<"rc="  << rc
 477         <<"err=" << err
 478         << "\n";
 479     }
 480   }
 481   return NULL;//H->Projection(kEta, kPhi);
 482 }
 483
 484
 485 //________________________________________________________
 486 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
 487 {
 488 // Store resolution/pulls of Kalman before updating with the TRD information
 489 // at the radial position of the first tracklet. The following points are used
 490 // for comparison
 491 //  - the (y,z,snp) of the first TRD tracklet
 492 //  - the (y, z, snp, tgl, pt) of the MC track reference
 493 //
 494 // Additionally the momentum resolution/pulls are calculated for usage in the
 495 // PID calculation.
 496   //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
 497
 498   if(track) fkTrack = track;
 499   if(!fkTrack){
 500     AliDebug(4, "No Track defined.");
 501     return NULL;
 502   }
 503   //fkTrack->Print();
 504   // check container
 505   THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
 506   if(!H){
 507     AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
 508     return NULL;
 509   }
 510   // check input track status
 511   AliExternalTrackParam *tin(NULL);
 512   if(!(tin = fkTrack->GetTrackIn())){
 513     AliError("Track did not entered TRD fiducial volume.");
 514     return NULL;
 515   }
 516   // check first tracklet
 517   AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
 518   if(!fTracklet){
 519     AliDebug(3, "No Tracklet in ly[0]. Skip track.");
 520     return NULL;
 521   }
 522   // check radial position
 523   Double_t x = tin->GetX();
 524   if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
 525     AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
 526     return NULL;
 527   }
 528   //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
 529
 530   Int_t bc(TMath::Abs(fkESD->GetTOFbc())%2);
 531   const Double_t *parR(tin->GetParameter());
 532   Double_t dyt(parR[0] - fTracklet->GetYfit(0)), dzt(parR[1] - fTracklet->GetZfit(0)),
 533             phit(fTracklet->GetYfit(1)),
 534             tilt(fTracklet->GetTilt());
 535
 536   // correct for tilt rotation
 537   Double_t dy  = dyt - dzt*tilt,
 538            dz  = dzt + dyt*tilt;
 539   phit       += tilt*parR[3];
 540   Double_t dphi = TMath::ASin(parR[2])-TMath::ATan(phit);
 541
 542   Double_t val[kNdim];
 543   val[kBC]          = bc;
 544   val[kPhi]         = fPhi;
 545   val[kEta]         = fEta;
 546   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fSpecies;
 547   val[kPt]          = GetPtBin(fPt);
 548   val[kYrez]        = dy;
 549   val[kZrez]        = dz;
 550   val[kPrez]        = dphi*TMath::RadToDeg();
 551   H->Fill(val);
 552   if(DebugLevel()>=3){
 553     (*DebugStream()) << "trackIn"
 554       <<"tracklet.="  << fTracklet
 555       <<"trackIn.="   << tin
 556       << "\n";
 557   }
 558
 559   if(!HasMCdata()) return NULL; // H->Projection(kEta, kPhi);
 560   if(!(H = (THnSparseI*)fContainer->At(kMCtrackIn))) {
 561     AliError(Form("Missing container @ %d", Int_t(kMCtrackIn)));
 562     return NULL;
 563   }
 564
 565   // get MC info
 566   UChar_t s;
 567   Float_t pt0, eta, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
 568   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) return NULL;
 569   dyt = y0 - fTracklet->GetYfit(0);
 570   dzt = z0 - fTracklet->GetZfit(0);
 571   phit= fTracklet->GetYfit(1) + tilt*dzdx0;
 572   Float_t phi = TMath::ATan2(y0, x0);
 573   dy  = dyt - dzt*tilt;
 574   dz  = dzt + dyt*tilt;
 575   dphi= TMath::ASin(dydx0)-TMath::ATan(phit);
 576
 577   Int_t pdg = fkMC->GetPDG(),
 578         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 579         sign(0);
 580   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 581   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 582   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 583
 584
 585   val[kBC]          = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
 586   val[kPhi]         = phi;
 587   val[kEta]         = eta;
 588   val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:sign*(sIdx+1);
 589   val[kPt]          = GetPtBin(pt0);
 590   val[kYrez]        = dy;
 591   val[kZrez]        = dz;
 592   val[kPrez]        = dphi*TMath::RadToDeg();
 593   H->Fill(val);
 594
 595   return NULL; //H->Projection(kEta, kPhi);
 596 }
 597
 598 //________________________________________________________
 599 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
 600 {
 601 // Store resolution/pulls of Kalman after last update with the TRD information
 602 // at the radial position of the first tracklet. The following points are used
 603 // for comparison
 604 //  - the (y,z,snp) of the first TRD tracklet
 605 //  - the (y, z, snp, tgl, pt) of the MC track reference
 606 //
 607 // Additionally the momentum resolution/pulls are calculated for usage in the
 608 // PID calculation.
 609
 610   if(track) fkTrack = track;
 611   return NULL;
 612 }
 613
 614 //________________________________________________________
 615 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
 616 {
 617   //
 618   // Plot MC distributions
 619   //
 620
 621   if(!HasMCdata()){
 622     AliDebug(2, "No MC defined. Results will not be available.");
 623     return NULL;
 624   }
 625   if(track) fkTrack = track;
 626   if(!fkTrack){
 627     AliDebug(4, "No Track defined.");
 628     return NULL;
 629   }
 630   if(!fContainer){
 631     AliWarning("No output container defined.");
 632     return NULL;
 633   }
 634   // retriev track characteristics
 635   Int_t pdg = fkMC->GetPDG(),
 636         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
 637         sign(0),
 638         sgm[3],
 639         label(fkMC->GetLabel()),
 640         fSegmentLevel(0);
 641   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
 642   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
 643   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
 644
 645   TObjArray *arr(NULL);TH1 *h(NULL);
 646   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 647   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
 648   UChar_t s;
 649   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
 650   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
 651   Double_t covR[7]/*, cov[3]*/;
 652
 653   if(DebugLevel()>=3){
 654     // get first detector
 655     Int_t det = -1;
 656     for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 657       if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
 658       det = fTracklet->GetDetector();
 659       break;
 660     }
 661     if(det>=0){
 662       TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
 663       Double_t m(-1.);
 664       m = fkTrack->GetMass();
 665       if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
 666         (*DebugStream()) << "MCkalman"
 667           << "pdg=" << pdg
 668           << "det=" << det
 669           << "x="   << &X
 670           << "y="   << &Y
 671           << "z="   << &Z
 672           << "dx="  << &dX
 673           << "dy="  << &dY
 674           << "dz="  << &dZ
 675           << "pt="  << &vPt
 676           << "dpt=" << &dPt
 677           << "bgt=" << &budget
 678           << "cov=" << &cCOV
 679           << "\n";
 680       }
 681     }
 682   }
 683   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 684     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
 685        !fTracklet->IsOK())*/ continue;
 686
 687     sgm[2] = fTracklet->GetDetector();
 688     sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
 689     sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
 690     Double_t tilt(fTracklet->GetTilt())
 691             ,t2(tilt*tilt)
 692             ,corr(1./(1. + t2))
 693             ,cost(TMath::Sqrt(corr));
 694     x0  = fTracklet->GetX0();
 695     //radial shift with respect to the MC reference (radial position of the pad plane)
 696     x= fTracklet->GetX();
 697     Bool_t rc(fTracklet->IsRowCross()); Float_t eta;
 698     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, s)) continue;
 699     xAnode  = fTracklet->GetX0();
 700
 701     // MC track position at reference radial position
 702     dx  = x0 - x;
 703     if(DebugLevel()>=4){
 704       (*DebugStream()) << "MC"
 705         << "det="     << sgm[2]
 706         << "pdg="     << pdg
 707         << "sgn="     << sign
 708         << "pt="      << pt0
 709         << "x="       << x0
 710         << "y="       << y0
 711         << "z="       << z0
 712         << "dydx="    << dydx0
 713         << "dzdx="    << dzdx0
 714         << "\n";
 715     }
 716     Float_t ymc = y0 - dx*dydx0;
 717     Float_t zmc = z0 - dx*dzdx0;
 718     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
 719
 720     // Kalman position at reference radial position
 721     dx = xAnode - x;
 722     dydx = fTracklet->GetYref(1);
 723     dzdx = fTracklet->GetZref(1);
 724     dzdl = fTracklet->GetTgl();
 725     y  = fTracklet->GetYref(0) - dx*dydx;
 726     dy = y - ymc;
 727     z  = fTracklet->GetZref(0) - dx*dzdx;
 728     dz = z - zmc;
 729     pt = TMath::Abs(fTracklet->GetPt());
 730     fTracklet->GetCovRef(covR);
 731
 732     arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
 733     // y resolution/pulls
 734     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
 735     ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
 736     // z resolution/pulls
 737     ((TH2S*)arr->At(2))->Fill(dzdx0, dz);
 738     ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
 739     // phi resolution/ snp pulls
 740     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
 741     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
 742     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
 743     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
 744     // theta resolution/ tgl pulls
 745     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
 746               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
 747     ((TH2I*)arr->At(6))->Fill(dzdl0,
 748     TMath::ATan(dtgl));
 749     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
 750     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
 751     Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
 752              p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
 753     ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
 754     ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
 755     ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
 756
 757     // Fill Debug stream for Kalman track
 758     if(DebugLevel()>=4){
 759       (*DebugStream()) << "MCtrack"
 760         << "pt="      << pt
 761         << "x="       << x
 762         << "y="       << y
 763         << "z="       << z
 764         << "dydx="    << dydx
 765         << "dzdx="    << dzdx
 766         << "s2y="     << covR[0]
 767         << "s2z="     << covR[2]
 768         << "\n";
 769     }
 770
 771     // recalculate tracklet based on the MC info
 772     AliTRDseedV1 tt(*fTracklet);
 773     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
 774     tt.SetZref(1, dzdx0);
 775     tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
 776     tt.Fit(1);
 777     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
 778     dydx = tt.GetYfit(1);
 779     dx = x0 - x;
 780     ymc = y0 - dx*dydx0;
 781     zmc = z0 - dx*dzdx0;
 782     dy = y-ymc;
 783     dz = z-zmc;
 784     Float_t dphi  = (dydx - dydx0);
 785     dphi /= (1.- dydx*dydx0);
 786
 787     // add tracklet residuals for y and dydx
 788     arr = (TObjArray*)fContainer->At(kMCtracklet);
 789
 790     if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
 791     if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
 792     ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
 793     if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
 794     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
 795
 796     // Fill Debug stream for tracklet
 797     if(DebugLevel()>=4){
 798       Float_t s2y = tt.GetS2Y();
 799       Float_t s2z = tt.GetS2Z();
 800       (*DebugStream()) << "MCtracklet"
 801         << "rc="    << rc
 802         << "x="     << x
 803         << "y="     << y
 804         << "z="     << z
 805         << "dydx="  << dydx
 806         << "s2y="   << s2y
 807         << "s2z="   << s2z
 808         << "\n";
 809     }
 810
 811     AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
 812     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
 813     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
 814
 815     arr = (TObjArray*)fContainer->At(kMCcluster);
 816     AliTRDcluster *c = NULL;
 817     tt.ResetClusterIter(kFALSE);
 818     while((c = tt.PrevCluster())){
 819       Float_t  q = TMath::Abs(c->GetQ());
 820       x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()]; y = c->GetY();z = c->GetZ();
 821       dx = x0 - x;
 822       ymc= y0 - dx*dydx0;
 823       zmc= z0 - dx*dzdx0;
 824       dy = cost*(y - ymc - tilt*(z-zmc));
 825       dz = cost*(z - zmc + tilt*(y-ymc));
 826
 827       // Fill Histograms
 828       if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){
 829         ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
 830         ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
 831       }
 832
 833       // Fill calibration container
 834       Float_t d = zr0 - zmc;
 835       d -= ((Int_t)(2 * d)) / 2.0;
 836       if (d > 0.25) d  = 0.5 - d;
 837       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
 838       clInfo->SetCluster(c);
 839       clInfo->SetMC(pdg, label);
 840       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
 841       clInfo->SetResolution(dy);
 842       clInfo->SetAnisochronity(d);
 843       clInfo->SetDriftLength(dx);
 844       clInfo->SetTilt(tilt);
 845       if(fMCcl) fMCcl->Add(clInfo);
 846       else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
 847       if(DebugLevel()>=5){
 848         if(!clInfoArr){
 849           clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
 850           clInfoArr->SetOwner(kFALSE);
 851         }
 852         clInfoArr->Add(clInfo);
 853       }
 854     }
 855     // Fill Debug Tree
 856     if(DebugLevel()>=5 && clInfoArr){
 857       (*DebugStream()) << "MCcluster"
 858         <<"clInfo.=" << clInfoArr
 859         << "\n";
 860       clInfoArr->Clear();
 861     }
 862   }
 863   if(clInfoArr) delete clInfoArr;
 864   return h;
 865 }
 866
 867
 868 //__________________________________________________________________________
 869 Int_t AliTRDresolution::GetPtBin(Float_t pt)
 870 {
 871 // Find pt bin according to local pt segmentation
 872   Int_t ipt(-1);
 873   while(ipt<AliTRDresolution::kNpt){
 874     if(pt<fgPtBin[ipt+1]) break;
 875     ipt++;
 876   }
 877   return ipt;
 878 }
 879
 880 //________________________________________________________
 881 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
 882 {
 883 // return mean number of entries/bin of histogram "h"
 884 // if option "opt" is given the following values are accepted:
 885 //   "<" : consider only entries less than "cut"
 886 //   ">" : consider only entries greater than "cut"
 887
 888   //Int_t dim(h->GetDimension());
 889   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
 890   Double_t sum(0.); Int_t n(0);
 891   for(Int_t ix(1); ix<=nbx; ix++)
 892     for(Int_t iy(1); iy<=nby; iy++)
 893       for(Int_t iz(1); iz<=nbz; iz++){
 894         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
 895         else{
 896           if(strcmp(opt, "<")==0) {
 897             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 898           } else if(strcmp(opt, ">")==0){
 899             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
 900           } else {sum += h->GetBinContent(ix, iy, iz); n++;}
 901         }
 902       }
 903   return n>0?sum/n:0.;
 904 }
 905
 906 //________________________________________________________
 907 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
 908 {
 909   //
 910   // Get the reference figures
 911   //
 912
 913   if(!gPad){
 914     AliWarning("Please provide a canvas to draw results.");
 915     return kFALSE;
 916   }
 917 /*  Int_t selection[100], n(0), selStart(0); //
 918   Int_t ly0(0), dly(5);
 919   TList *l(NULL); TVirtualPad *pad(NULL); */
 920   switch(ifig){
 921   case 0:
 922     break;
 923   }
 924   AliWarning(Form("Reference plot [%d] missing result", ifig));
 925   return kFALSE;
 926 }
 927
 928
 929 //________________________________________________________
 930 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
 931 {
 932 // Build pt segments
 933   for(Int_t j(0); j<=kNpt; j++){
 934     pt0+=(TMath::Exp(j*j*dpt)-1.);
 935     fgPtBin[j]=pt0;
 936   }
 937 }
 938
 939 //________________________________________________________
 940 void AliTRDresolution::MakeSummary()
 941 {
 942 // Build summary plots
 943
 944   if(!fProj){
 945     AliError("Missing results");
 946     return;
 947   }
 948   TVirtualPad *p(NULL); TCanvas *cOut(NULL);
 949   TObjArray *arr(NULL); TH2 *h2(NULL);
 950
 951   // cluster resolution
 952   // define palette
 953   gStyle->SetPalette(1);
 954   const Int_t nClViews(8);
 955   const Char_t *vClName[nClViews] = {"HClY", "HClYn", "HClYp", "HClQn", "HClQp", "HClYXTCp", "HClYXTCn", "HClYXPh"};
 956   const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0};
 957   if((arr = (TObjArray*)fProj->At(kCluster))){
 958     for(Int_t iview(0); iview<nClViews; iview++){
 959       cOut = new TCanvas(Form("TRDsummary%s_Cl%02d", GetName(), iview), "Cluster Resolution", 1024, 768);
 960       cOut->Divide(3,2, 1.e-5, 1.e-5);
 961       Int_t nplot(0);
 962       for(Int_t iplot(0); iplot<6; iplot++){
 963         p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
 964         if(!(h2 = (TH2*)arr->FindObject(Form("%s%d_2D", vClName[iview], iplot)))) continue;
 965         nplot++;
 966         if(vClOpt[iview]==0) h2->Draw("colz");
 967         else if(vClOpt[iview]==1) DrawSigma(h2, 1.e4, 2.e2, 6.e2, "#sigma(#Deltay) [#mum]");
 968       }
 969       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 970       else delete cOut;
 971     }
 972   }
 973
 974   // tracklet systematic
 975   const Int_t nTrkltViews(10);
 976   const Char_t *vTrkltName[nTrkltViews] = {"HTrkltY", "HTrkltYn", "HTrkltYp", "HTrkltPhn", "HTrkltPhp", "HTrkltZ", "HTrkltQn", "HTrkltQp", "HTrkltPn", "HTrkltPp"};
 977   if((arr = (TObjArray*)fProj->At(kTracklet))){
 978     for(Int_t iview(0); iview<nTrkltViews; iview++){
 979       cOut = new TCanvas(Form("TRDsummary%s_Trklt%02d", GetName(), iview), "Tracklet Resolution", 1024, 768);
 980       cOut->Divide(3,2, 1.e-5, 1.e-5);
 981       Int_t nplot(0);
 982       for(Int_t iplot(0); iplot<6; iplot++){
 983         p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
 984         if(!(h2 = (TH2*)arr->FindObject(Form("%s%d_2D", vTrkltName[iview], iplot)))) continue;
 985         h2->Draw("colz"); nplot++;
 986       }
 987       if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
 988       else delete cOut;
 989     }
 990   }
 991   // trackIn systematic
 992   const Char_t *hname[] = {"HTrkInY", "HTrkInYn", "HTrkInYp", "HTrkInZ", "HTrkInPhn", "HTrkInPhp"};
 993   if((arr = (TObjArray*)fProj->At(kTrackIn))){
 994     cOut = new TCanvas(Form("TRDsummary%s_TrkIn", GetName()), "Track IN Resolution", 1024, 768);
 995     cOut->Divide(3,2, 1.e-5, 1.e-5);
 996     Int_t nplot(0);
 997     for(Int_t iplot(0); iplot<6; iplot++){
 998       p=cOut->cd(iplot+1);    p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
 999         if(!(h2 = (TH2*)arr->FindObject(Form("%s_2D", hname[iplot])))) continue;
1000         h2->Draw("colz"); nplot++;
1001     }
1002     if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1003     else delete cOut;
1004   }
1005   gStyle->SetPalette(1);
1006 }
1007
1008 //________________________________________________________
1009 void AliTRDresolution::DrawSigma(TH2 *h2, Float_t scale, Float_t m, Float_t M, const Char_t *title)
1010 {
1011   // Draw error bars scaled with "scale" instead of content values
1012   //use range [m,M] if limits are specified
1013
1014   if(!h2) return;
1015   TH2 *h2e = (TH2F*)h2->Clone(Form("%s_E", h2->GetName()));
1016   h2e->SetContour(10);
1017   if(M>m) h2e->GetZaxis()->SetRangeUser(m, M);
1018   if(title) h2e->GetZaxis()->SetTitle(title);
1019
1020   for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1021     for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1022       if(h2->GetBinContent(ix, iy)<-100.) continue;
1023       Float_t v(scale*h2->GetBinError(ix, iy));
1024       if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1025       h2e->SetBinContent(ix, iy, v);
1026     }
1027   }
1028   h2e->Draw("colz");
1029 }
1030
1031 //________________________________________________________
1032 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1033 {
1034 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1035 // The "range" vector should be initialized with 2 elements
1036 // Option "mod" can be any of
1037 //   - 0 : gaussian like distribution
1038 //   - 1 : tailed distribution
1039
1040   Int_t nx(h2->GetNbinsX())
1041       , ny(h2->GetNbinsY())
1042       , n(nx*ny);
1043   Double_t *data=new Double_t[n];
1044   for(Int_t ix(1), in(0); ix<=nx; ix++){
1045     for(Int_t iy(1); iy<=ny; iy++)
1046       data[in++] = h2->GetBinContent(ix, iy);
1047   }
1048   Double_t mean, sigm;
1049   AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1050
1051   range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1052   if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1053   AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1054   TH1S h1("h1SF0", "", 100, range[0], range[1]);
1055   h1.FillN(n,data,0);
1056   delete [] data;
1057
1058   switch(mod){
1059   case 0:// gaussian distribution
1060   {
1061     TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1062     h1.Fit(&fg, "QN");
1063     mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1064     range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1065     AliDebug(2, Form("     rangeG[%f %f]", range[0], range[1]));
1066     break;
1067   }
1068   case 1:// tailed distribution
1069   {
1070     Int_t bmax(h1.GetMaximumBin());
1071     Int_t jBinMin(1), jBinMax(100);
1072     for(Int_t ibin(bmax); ibin--;){
1073       if(h1.GetBinContent(ibin)<1.){
1074         jBinMin=ibin; break;
1075       }
1076     }
1077     for(Int_t ibin(bmax); ibin++;){
1078       if(h1.GetBinContent(ibin)<1.){
1079         jBinMax=ibin; break;
1080       }
1081     }
1082     range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1083     AliDebug(2, Form("     rangeT[%f %f]", range[0], range[1]));
1084     break;
1085   }
1086   }
1087
1088   return;
1089 }
1090
1091
1092 //________________________________________________________
1093 Bool_t AliTRDresolution::MakeProjectionCluster()
1094 {
1095 // Analyse cluster
1096   const Int_t kNcontours(9);
1097   const Int_t kNstat(300);
1098   Int_t cidx = kCluster;
1099   if(fProj && fProj->At(cidx)) return kTRUE;
1100   if(!fContainer){
1101     AliError("Missing data container.");
1102     return kFALSE;
1103   }
1104   THnSparse *H(NULL);
1105   if(!(H = (THnSparse*)fContainer->At(cidx))){
1106     AliError(Form("Missing/Wrong data @ %d.", cidx));
1107     return kFALSE;
1108   }
1109   Int_t ndim(H->GetNdimensions());
1110   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1111   TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1112   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1113   if(ndim > kPt) apt = H->GetAxis(kPt);
1114   if(ndim > kSpeciesChgRC) as  = H->GetAxis(kSpeciesChgRC);
1115   // build list of projections
1116   const Int_t nsel(12), npsel(5);
1117   // define rebinning strategy
1118   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1119   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1120   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1121   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1122     isel++; // new selection
1123     hp[ih].Build(Form("HClY%d", ily), Form("Clusters :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1124     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1125       php[isel][np[isel]++] = &hp[ih++];
1126     hp[ih].Build(Form("HClYn%d", ily), Form("Clusters[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1127     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1128       php[isel][np[isel]++] = &hp[ih++];
1129     hp[ih].Build(Form("HClQn%d", ily), Form("Clusters[-]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1130     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1131     hp[ih].SetShowRange(20., 40.);
1132       php[isel][np[isel]++] = &hp[ih++];
1133     hp[ih].Build(Form("HClYXTCn%d", ily), Form("Clusters[-]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1134 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1135       php[isel][np[isel]++] = &hp[ih++];
1136     hp[ih].Build(Form("HClYXPh%d", ily), Form("Clusters :: r-#phi(x,#Phi) residuals ly%d", ily), kPrez, kPt, kYrez, aa);
1137 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1138       php[isel][np[isel]++] = &hp[ih++];
1139     isel++; // new selection
1140       php[isel][np[isel]++] = &hp[ih-5]; // relink HClY
1141     hp[ih].Build(Form("HClYp%d", ily), Form("Clusters[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1142     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1143       php[isel][np[isel]++] = &hp[ih++];
1144     hp[ih].Build(Form("HClQp%d", ily), Form("Clusters[+]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1145     hp[ih].SetShowRange(20., 40.);
1146     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1147       php[isel][np[isel]++] = &hp[ih++];
1148     hp[ih].Build(Form("HClYXTCp%d", ily), Form("Clusters[+]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1149 //    hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1150       php[isel][np[isel]++] = &hp[ih++];
1151       php[isel][np[isel]++] = &hp[ih-4]; // relink HClYXPh
1152   }
1153
1154   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1);
1155   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1156     v = H->GetBinContent(ib, coord); if(v<1.) continue;
1157     ly = coord[kBC]-1;
1158     // RC selection
1159     if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1160
1161     // charge selection
1162     ch = 0; // [-] track
1163     if(chBin>0 && coord[kPt] > chBin) ch = 1;  // [+] track
1164
1165     isel = ly*2+ch;
1166     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1167   }
1168
1169   if(!fProj){
1170     AliInfo("Building array of projections ...");
1171     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1172   }
1173   TObjArray *arr(NULL);
1174   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1175
1176   TH2 *h2(NULL);
1177   for(; ih--; ){
1178     if(!hp[ih].fH) continue;
1179     Int_t mid(1), nstat(kNstat);
1180     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;}
1181     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1182     arr->AddAt(h2, ih);
1183   }
1184
1185   return kTRUE;
1186 }
1187
1188 //________________________________________________________
1189 Bool_t AliTRDresolution::MakeProjectionTracklet()
1190 {
1191 // Analyse tracklet
1192   const Int_t kNcontours(9);
1193   const Int_t kNstat(100);
1194   Int_t cidx = kTracklet;
1195   if(fProj && fProj->At(cidx)) return kTRUE;
1196   if(!fContainer){
1197     AliError("Missing data container.");
1198     return kFALSE;
1199   }
1200   THnSparse *H(NULL);
1201   if(!(H = (THnSparse*)fContainer->At(cidx))){
1202     AliError(Form("Missing/Wrong data @ %d.", cidx));
1203     return kFALSE;
1204   }
1205   Int_t ndim(H->GetNdimensions());
1206   Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1207   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1208   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1209   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1210   // build list of projections
1211   const Int_t nsel(18), npsel(6);
1212   // define rebinning strategy
1213   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1214   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1215   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1216   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1217     isel++; // new selection
1218     hp[ih].Build(Form("HTrkltY%d", ily), Form("Tracklets   :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1219     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1220       php[isel][np[isel]++] = &hp[ih++];
1221     hp[ih].Build(Form("HTrkltYn%d", ily), Form("Tracklets[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1222     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1223       php[isel][np[isel]++] = &hp[ih++];
1224     hp[ih].Build(Form("HTrkltPhn%d", ily), Form("Tracklets[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1225     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1226       php[isel][np[isel]++] = &hp[ih++];
1227     hp[ih].Build(Form("HTrkltPn%d", ily), Form("Tracklets[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1228     hp[ih].SetShowRange(6.,12.);
1229     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1230       php[isel][np[isel]++] = &hp[ih++];
1231     hp[ih].Build(Form("HTrkltYPn%d", ily), Form("Tracklets[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1232       php[isel][np[isel]++] = &hp[ih++];
1233     hp[ih].Build(Form("HTrkltQn%d", ily), Form("Tracklets[-]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1234     hp[ih].SetShowRange(700.,1100.);
1235     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1236       php[isel][np[isel]++] = &hp[ih++];
1237     isel++; // new selection
1238     php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1239     hp[ih].Build(Form("HTrkltYp%d", ily), Form("Tracklets[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1240     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1241       php[isel][np[isel]++] = &hp[ih++];
1242     hp[ih].Build(Form("HTrkltPhp%d", ily), Form("Tracklets[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1243     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1244       php[isel][np[isel]++] = &hp[ih++];
1245     hp[ih].Build(Form("HTrkltPp%d", ily), Form("Tracklets[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1246     hp[ih].SetShowRange(6.,12.);
1247     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1248       php[isel][np[isel]++] = &hp[ih++];
1249     hp[ih].Build(Form("HTrkltYPp%d", ily), Form("Tracklets[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1250       php[isel][np[isel]++] = &hp[ih++];
1251     hp[ih].Build(Form("HTrkltQp%d", ily), Form("Tracklets[+]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1252     hp[ih].SetShowRange(700.,1100.);
1253     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1254       php[isel][np[isel]++] = &hp[ih++];
1255     isel++; // new selection
1256     hp[ih].Build(Form("HTrkltZ%d", ily), Form("Tracklets[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1257     hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1258       php[isel][np[isel]++] = &hp[ih++];
1259   }
1260
1261   Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1262   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1263     v = H->GetBinContent(ib, coord);
1264     if(v<1.) continue;
1265     ly = coord[kBC]-1; // layer selection
1266     // charge selection
1267     ch = 0; // [-] track
1268     if(rcBin>0){ // debug mode in which species are also saved
1269       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1270       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1271     }
1272     isel = ly*3+ch;
1273     for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1274   }
1275
1276   if(!fProj){
1277     AliInfo("Building array of projections ...");
1278     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1279   }
1280   TObjArray *arr(NULL);
1281   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1282
1283   TH2 *h2(NULL);
1284   for(; ih--; ){
1285     if(!hp[ih].fH) continue;
1286     Int_t mid(0), nstat(kNstat);
1287     if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;}
1288     if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1289     arr->AddAt(h2, ih);
1290   }
1291   return kTRUE;
1292 }
1293
1294 //________________________________________________________
1295 Bool_t AliTRDresolution::MakeProjectionTrackIn()
1296 {
1297 // Analyse track in
1298
1299   const Int_t kNcontours(9);
1300   const Int_t kNstat(30);
1301
1302   Int_t cidx = kTrackIn;
1303   if(fProj && fProj->At(cidx)) return kTRUE;
1304   if(!fContainer){
1305     AliError("Missing data container.");
1306     return kFALSE;
1307   }
1308   THnSparse *H(NULL);
1309   if(!(H = (THnSparse*)fContainer->At(cidx))){
1310     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1311     return kFALSE;
1312   }
1313
1314   Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1315   Int_t ndim(H->GetNdimensions());
1316   TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1317   for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1318   if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1319   // build list of projections
1320   const Int_t nsel(3), npsel(4);
1321   // define rebinning strategy
1322   const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1323   const Int_t nPtPhi(2); Int_t rebinPtPhiX[nEtaPhi] = {1, 1}, rebinPtPhiY[nEtaPhi] = {2, 5};
1324   AliTRDresolutionProjection hp[fgkNproj[cidx]], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1325   Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1326   // define list of projections
1327   isel++;  // negative tracks
1328   hp[ih].Build("HTrkInY", "TrackIn :: r-#phi residuals", kEta, kPhi, kYrez, aa);
1329   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1330     php[isel][np[isel]++] = &hp[ih++];
1331   hp[ih].Build("HTrkInYn", "TrackIn[-]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1332   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1333     php[isel][np[isel]++] = &hp[ih++];
1334   hp[ih].Build("HTrkInPhn", "TrackIn[-]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1335   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1336     php[isel][np[isel]++] = &hp[ih++];
1337   hp[ih].Build("HTrkInYPn", "TrackIn[-]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1338   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1339     php[isel][np[isel]++] = &hp[ih++];
1340   isel++; // positive tracks
1341   php[isel][np[isel]++] = &hp[ih-4]; // relink first histo
1342   hp[ih].Build("HTrkInYp", "TrackIn[+]:: r-#phi residuals", kEta, kPhi, kYrez, aa);
1343   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1344     php[isel][np[isel]++] = &hp[ih++];
1345   hp[ih].Build("HTrkInPhp", "TrackIn[+]:: #Delta#phi residuals", kEta, kPhi, kPrez, aa);
1346   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1347     php[isel][np[isel]++] = &hp[ih++];
1348   hp[ih].Build("HTrkInYPp", "TrackIn[+]:: r-#phi/p_{t} residuals", kPt, kPhi, kYrez, aa);
1349   hp[ih].SetRebinStrategy(nPtPhi, rebinPtPhiX, rebinPtPhiY);
1350     php[isel][np[isel]++] = &hp[ih++];
1351   isel++; // RC tracks
1352   hp[ih].Build("HTrkInZ", "TrackIn[RC]:: z residuals", kEta, kPhi, kZrez, aa);
1353   hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1354     php[isel][np[isel]++] = &hp[ih++];
1355
1356   // fill projections
1357   Int_t ch(0), rcBin(as?as->FindBin(0.):-1);
1358   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1359     v = H->GetBinContent(ib, coord);
1360     if(v<1.) continue;
1361     if(coord[kBC]>1) continue; // bunch cross cut
1362     // charge selection
1363     ch = 0; // [-] track
1364     if(rcBin>0){ // debug mode in which species are also saved
1365       if(coord[kSpeciesChgRC] > rcBin) ch = 1;  // [+] track
1366       else if(coord[kSpeciesChgRC] == rcBin) ch = 2;  // [RC] track
1367     }
1368     for(Int_t jh(0); jh<np[ch]; jh++) php[ch][jh]->Increment(coord, v);
1369   }
1370   if(!fProj){
1371     AliInfo("Building array of projections ...");
1372     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1373   }
1374   TObjArray *arr(NULL);
1375   fProj->AddAt(arr = new TObjArray(fgkNproj[cidx]), cidx);
1376
1377   TH2 *h2(NULL);
1378   for(; ih--; ){
1379     if(!hp[ih].fH) continue;
1380     if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours/*, mid*/))) continue;
1381     arr->AddAt(h2, ih);
1382   }
1383   return kTRUE;
1384 }
1385
1386
1387
1388 //________________________________________________________
1389 Bool_t AliTRDresolution::PostProcess()
1390 {
1391 // Fit, Project, Combine, Extract values from the containers filled during execution
1392
1393   if (!fContainer) {
1394     AliError("ERROR: list not available");
1395     return kFALSE;
1396   }
1397
1398   // DEFINE MODELS
1399   // simple gauss
1400   TF1 fg("fGauss", "gaus", -.5, .5);
1401   // Landau for charge resolution
1402   TF1 fch("fClCh", "landau", 0., 1000.);
1403   // Landau for e+- pt resolution
1404   TF1 fpt("fPt", "landau", -0.1, 0.2);
1405
1406   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1407   // Clusters residuals
1408   if(!MakeProjectionCluster()) return kFALSE;
1409   fNRefFigures = 3;
1410   // Tracklet residual/pulls
1411   if(!MakeProjectionTracklet()) return kFALSE;
1412   fNRefFigures = 7;
1413   // TRDin residual/pulls
1414   if(!MakeProjectionTrackIn()) return kFALSE;
1415   fNRefFigures = 11;
1416   // TRDout residual/pulls
1417 //  if(!MakeProjectionTrackOut()) return kFALSE;
1418   fNRefFigures = 15;
1419
1420   if(!HasMCdata()) return kTRUE;
1421
1422
1423   //PROCESS MC RESIDUAL DISTRIBUTIONS
1424
1425   // CLUSTER Y RESOLUTION/PULLS
1426 //  if(!MakeProjectionClusterMC()) return kFALSE;
1427   fNRefFigures = 17;
1428
1429   // TRACKLET RESOLUTION/PULLS
1430 //  if(!MakeProjectionTrackletMC()) return kFALSE;
1431   fNRefFigures = 21;
1432
1433   // TRACK RESOLUTION/PULLS
1434 /*  Process3Darray(kMCtrack, 0, &fg, 1.e4);   // y
1435   Process3DlinkedArray(kMCtrack, 1, &fg);   // y PULL
1436   Process3Darray(kMCtrack, 2, &fg, 1.e4);   // z
1437   Process3Darray(kMCtrack, 3, &fg);         // z PULL
1438   Process2Darray(kMCtrack, 4, &fg, 1.e3);   // phi
1439   Process2Darray(kMCtrack, 5, &fg);         // snp PULL
1440   Process2Darray(kMCtrack, 6, &fg, 1.e3);   // theta
1441   Process2Darray(kMCtrack, 7, &fg);         // tgl PULL
1442   Process3Darray(kMCtrack, 8, &fg, 1.e2);   // pt resolution
1443   Process3Darray(kMCtrack, 9, &fg);         // 1/pt pulls
1444   Process3Darray(kMCtrack, 10, &fg, 1.e2);  // p resolution*/
1445 //  if(!MakeProjectionTrackMC(kMCtrack)) return kFALSE;
1446   fNRefFigures+=16;
1447
1448   // TRACK TRDin RESOLUTION/PULLS
1449 //  if(!MakeProjectionTrackMC(kMCtrackIn)) return kFALSE;
1450   fNRefFigures+=8;
1451
1452   // TRACK TRDout RESOLUTION/PULLS
1453 //  if(!MakeProjectionTrackMC(kMCtrackOut)) return kFALSE;
1454   fNRefFigures+=8;
1455
1456   return kTRUE;
1457 }
1458
1459
1460 //________________________________________________________
1461 void AliTRDresolution::Terminate(Option_t *opt)
1462 {
1463   AliTRDrecoTask::Terminate(opt);
1464   if(HasPostProcess()) PostProcess();
1465 }
1466
1467 //________________________________________________________
1468 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1469 {
1470 // Helper function to avoid duplication of code
1471 // Make first guesses on the fit parameters
1472
1473   // find the intial parameters of the fit !! (thanks George)
1474   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1475   Double_t sum = 0.;
1476   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1477   f->SetParLimits(0, 0., 3.*sum);
1478   f->SetParameter(0, .9*sum);
1479   Double_t rms = h->GetRMS();
1480   f->SetParLimits(1, -rms, rms);
1481   f->SetParameter(1, h->GetMean());
1482
1483   f->SetParLimits(2, 0., 2.*rms);
1484   f->SetParameter(2, rms);
1485   if(f->GetNpar() <= 4) return;
1486
1487   f->SetParLimits(3, 0., sum);
1488   f->SetParameter(3, .1*sum);
1489
1490   f->SetParLimits(4, -.3, .3);
1491   f->SetParameter(4, 0.);
1492
1493   f->SetParLimits(5, 0., 1.e2);
1494   f->SetParameter(5, 2.e-1);
1495 }
1496
1497 //________________________________________________________
1498 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1499 {
1500 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1501 //  - y reziduals/pulls
1502
1503   TObjArray *arr = new TObjArray(2);
1504   arr->SetName(name); arr->SetOwner();
1505   TH1 *h(NULL); char hname[100], htitle[300];
1506
1507   // tracklet resolution/pull in y direction
1508   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1509   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1510   Float_t rr = range<0.?fDyRange:range;
1511   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1512     Int_t nybins=50;
1513     if(expand) nybins*=2;
1514     h = new TH3S(hname, htitle,
1515                  48, -.48, .48,            // phi
1516                  60, -rr, rr,              // dy
1517                  nybins, -0.5, nybins-0.5);// segment
1518   } else h->Reset();
1519   arr->AddAt(h, 0);
1520   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1521   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1522   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1523     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1524   } else h->Reset();
1525   arr->AddAt(h, 1);
1526
1527   return arr;
1528 }
1529
1530 //________________________________________________________
1531 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1532 {
1533 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1534 //  - y reziduals/pulls
1535 //  - z reziduals/pulls
1536 //  - phi reziduals
1537   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
1538   arr->Expand(5);
1539   TH1 *h(NULL); char hname[100], htitle[300];
1540
1541   // tracklet resolution/pull in z direction
1542   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1543   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1544   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1545     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1546   } else h->Reset();
1547   arr->AddAt(h, 2);
1548   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1549   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1550   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1551     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1552     h->GetZaxis()->SetBinLabel(1, "no RC");
1553     h->GetZaxis()->SetBinLabel(2, "RC");
1554   } else h->Reset();
1555   arr->AddAt(h, 3);
1556
1557   // tracklet to track phi resolution
1558   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1559   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1560   Int_t nsgms=540;
1561   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1562     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1563   } else h->Reset();
1564   arr->AddAt(h, 4);
1565
1566   return arr;
1567 }
1568
1569 //________________________________________________________
1570 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1571 {
1572 // Build performance histograms for AliExternalTrackParam.vs MC
1573 //  - y resolution/pulls
1574 //  - z resolution/pulls
1575 //  - phi resolution, snp pulls
1576 //  - theta resolution, tgl pulls
1577 //  - pt resolution, 1/pt pulls, p resolution
1578
1579   TObjArray *arr = BuildMonitorContainerTracklet(name);
1580   arr->Expand(11);
1581   TH1 *h(NULL); char hname[100], htitle[300];
1582   //TAxis *ax(NULL);
1583
1584   // snp pulls
1585   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1586   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1587   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1588     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1589   } else h->Reset();
1590   arr->AddAt(h, 5);
1591
1592   // theta resolution
1593   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1594   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1595   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1596     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1597   } else h->Reset();
1598   arr->AddAt(h, 6);
1599   // tgl pulls
1600   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1601   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1602   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1603     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1604   } else h->Reset();
1605   arr->AddAt(h, 7);
1606
1607   const Int_t kNdpt(150);
1608   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1609   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1610   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1611   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1612   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1613   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1614
1615   // Pt resolution
1616   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1617   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1618   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1619     h = new TH3S(hname, htitle,
1620                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1621     //ax = h->GetZaxis();
1622     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1623   } else h->Reset();
1624   arr->AddAt(h, 8);
1625   // 1/Pt pulls
1626   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1627   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);
1628   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1629     h = new TH3S(hname, htitle,
1630                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1631     //ax = h->GetZaxis();
1632     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1633   } else h->Reset();
1634   arr->AddAt(h, 9);
1635   // P resolution
1636   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1637   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1638   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1639     h = new TH3S(hname, htitle,
1640                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1641     //ax = h->GetZaxis();
1642     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1643   } else h->Reset();
1644   arr->AddAt(h, 10);
1645
1646   return arr;
1647 }
1648
1649
1650 //________________________________________________________
1651 TObjArray* AliTRDresolution::Histos()
1652 {
1653   //
1654   // Define histograms
1655   //
1656
1657   if(fContainer) return fContainer;
1658
1659   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1660   THnSparse *H(NULL);
1661   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1662
1663   //++++++++++++++++++++++
1664   // cluster to tracklet residuals/pulls
1665   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1666   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1667     const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"};
1668     const Int_t clNbins[kNdim]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 33, 10, 30, 15};
1669     const Double_t clMin[kNdim]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, -2., -45},
1670                    clMax[kNdim]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 3.3, 2.1, 118., 45};
1671     st = "cluster spatial&charge resolution;";
1672     // define minimum info to be saved in non debug mode
1673     Int_t ndim=DebugLevel()>=1?kNdim:4;
1674     for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1675     H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1676   } else H->Reset();
1677   fContainer->AddAt(H, kCluster);
1678   //++++++++++++++++++++++
1679   // tracklet to TRD track
1680   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1681   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1682     Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1683     Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1684     Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1685     Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1686     // set specific fields
1687     trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1688     trTitle[kNdim]=StrDup("dq/dl [a.u.]"); trNbins[kNdim] = 30; trMin[kNdim] = 100.; trMax[kNdim] = 3100;
1689
1690     st = "tracklet spatial&charge resolution;";
1691     // define minimum info to be saved in non debug mode
1692     Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1693     for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1694     H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1695   } else H->Reset();
1696   fContainer->AddAt(H, kTracklet);
1697   //++++++++++++++++++++++
1698   // tracklet to TRDin
1699   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1700   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1701     st = "r-#phi/z/angular residuals @ TRD entry;";
1702     // define minimum info to be saved in non debug mode
1703     Int_t ndim=DebugLevel()>=1?kNdim:7;
1704     for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1705     H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, fgkMin, fgkMax);
1706   } else H->Reset();
1707   fContainer->AddAt(H, kTrackIn);
1708   // tracklet to TRDout
1709 //  fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1710
1711
1712   // Resolution histos
1713   if(!HasMCdata()) return fContainer;
1714
1715   // cluster resolution
1716   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1717
1718   // tracklet resolution
1719   fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1720
1721   // track resolution
1722   TObjArray *arr(NULL);
1723   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1724   arr->SetName("MCtrk");
1725   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1726
1727   // TRDin TRACK RESOLUTION
1728   fContainer->AddAt(H, kMCtrackIn);
1729
1730   // TRDout TRACK RESOLUTION
1731   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1732
1733   return fContainer;
1734 }
1735
1736 //________________________________________________________
1737 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1738 {
1739 // Robust function to process sigma/mean for 2D plot dy(x)
1740 // For each x bin a gauss fit is performed on the y projection and the range
1741 // with the minimum chi2/ndf is choosen
1742
1743   if(!h2) {
1744     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1745     return kFALSE;
1746   }
1747   if(!Int_t(h2->GetEntries())){
1748     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1749     return kFALSE;
1750   }
1751   if(!g || !g[0]|| !g[1]) {
1752     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1753     return kFALSE;
1754   }
1755   // prepare
1756   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1757   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1758   TF1 f("f", "gaus", ymin, ymax);
1759   Int_t n(0);
1760   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1761   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1762   TH1D *h(NULL);
1763   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1764   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1765
1766
1767   // do actual loop
1768   Float_t chi2OverNdf(0.);
1769   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1770     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1771     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1772
1773     h = h2->ProjectionY("py", ix, ix);
1774     if((n=(Int_t)h->GetEntries())<stat){
1775       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1776       continue;
1777     }
1778     // looking for a first order mean value
1779     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1780     h->Fit(&f, "QNW");
1781     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1782     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1783     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1784     ey  = f.GetParError(1);
1785     sy = f.GetParameter(2); esy = f.GetParError(2);
1786 //     // looking for the best chi2/ndf value
1787 //     while(ymin<y0 && ymax>y1){
1788 //       f.SetParameter(1, y);
1789 //       f.SetParameter(2, sy);
1790 //       h->Fit(&f, "QNW", "", y0, y1);
1791 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1792 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1793 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1794 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1795 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1796 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1797 //       }
1798 //       y0-=dy; y1+=dy;
1799 //     }
1800     g[0]->SetPoint(np, x, y);
1801     g[0]->SetPointError(np, ex, ey);
1802     g[1]->SetPoint(np, x, sy);
1803     g[1]->SetPointError(np, ex, esy);
1804     np++;
1805   }
1806   return kTRUE;
1807 }
1808
1809
1810 //________________________________________________________
1811 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1812 {
1813   //
1814   // Do the processing
1815   //
1816
1817   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1818   Int_t n = 0;
1819   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1820   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1821   if(Int_t(h2->GetEntries())){
1822     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1823   } else {
1824     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1825     fIdxPlot++;
1826     return kTRUE;
1827   }
1828
1829   const Int_t kINTEGRAL=1;
1830   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1831     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1832     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1833     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1834     if((n=(Int_t)h->GetEntries())<150){
1835       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1836       continue;
1837     }
1838     h->Fit(f, "QN");
1839     Int_t ip = g[0]->GetN();
1840     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)));
1841     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1842     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1843     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1844     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1845 /*
1846     g[0]->SetPoint(ip, x, k*h->GetMean());
1847     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1848     g[1]->SetPoint(ip, x, k*h->GetRMS());
1849     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1850   }
1851   fIdxPlot++;
1852   return kTRUE;
1853 }
1854
1855
1856 //____________________________________________________________________
1857 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
1858 {
1859 //
1860 // Fit track with a staight line using the "np" clusters stored in the array "points".
1861 // The following particularities are stored in the clusters from points:
1862 //   1. pad tilt as cluster charge
1863 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
1864 // The parameters of the straight line fit are stored in the array "param" in the following order :
1865 //     param[0] - x0 reference radial position
1866 //     param[1] - y0 reference r-phi position @ x0
1867 //     param[2] - z0 reference z position @ x0
1868 //     param[3] - slope dy/dx
1869 //     param[4] - slope dz/dx
1870 //
1871 // Attention :
1872 // Function should be used to refit tracks for B=0T
1873 //
1874
1875   if(np<40){
1876     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
1877     return kFALSE;
1878   }
1879   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
1880
1881   Double_t x0(0.);
1882   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
1883   x0/=Float_t(np);
1884
1885   Double_t x, y, z, dx, tilt(0.);
1886   for(Int_t ip(0); ip<np; ip++){
1887     x = points[ip].GetX(); z = points[ip].GetZ();
1888     dx = x - x0;
1889     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
1890   }
1891   if(zfitter.Eval() != 0) return kFALSE;
1892
1893   Double_t z0    = zfitter.GetParameter(0);
1894   Double_t dzdx  = zfitter.GetParameter(1);
1895   for(Int_t ip(0); ip<np; ip++){
1896     if(points[ip].GetClusterType()) continue;
1897     x    = points[ip].GetX();
1898     dx   = x - x0;
1899     y    = points[ip].GetY();
1900     z    = points[ip].GetZ();
1901     tilt = points[ip].GetCharge();
1902     y -= tilt*(-dzdx*dx + z - z0);
1903     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
1904     yfitter.AddPoint(&dx, y, 1.);
1905   }
1906   if(yfitter.Eval() != 0) return kFALSE;
1907   Double_t y0   = yfitter.GetParameter(0);
1908   Double_t dydx = yfitter.GetParameter(1);
1909
1910   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
1911   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);
1912   return kTRUE;
1913 }
1914
1915 //____________________________________________________________________
1916 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
1917 {
1918 //
1919 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
1920 // See function FitTrack for the data stored in the "clusters" array
1921
1922 // The parameters of the straight line fit are stored in the array "param" in the following order :
1923 //     par[0] - x0 reference radial position
1924 //     par[1] - y0 reference r-phi position @ x0
1925 //     par[2] - slope dy/dx
1926 //
1927 // Attention :
1928 // Function should be used to refit tracks for B=0T
1929 //
1930
1931   TLinearFitter yfitter(2, "pol1");
1932
1933   // grep data for tracklet
1934   Double_t x0(0.), x[60], y[60], dy[60];
1935   Int_t nly(0);
1936   for(Int_t ip(0); ip<np; ip++){
1937     if(points[ip].GetClusterType()) continue;
1938     if(points[ip].GetVolumeID() != ly) continue;
1939     Float_t xt(points[ip].GetX())
1940            ,yt(param[1] + param[3] * (xt - param[0]));
1941     x[nly] = xt;
1942     y[nly] = points[ip].GetY();
1943     dy[nly]= y[nly]-yt;
1944     x0    += xt;
1945     nly++;
1946   }
1947   if(nly<10){
1948     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
1949     return kFALSE;
1950   }
1951   // set radial reference for fit
1952   x0 /= Float_t(nly);
1953
1954   // find tracklet core
1955   Double_t mean(0.), sig(1.e3);
1956   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
1957
1958   // simple cluster error parameterization
1959   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
1960
1961   // fit tracklet core
1962   for(Int_t jly(0); jly<nly; jly++){
1963     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
1964     Double_t dx(x[jly]-x0);
1965     yfitter.AddPoint(&dx, y[jly], 1.);
1966   }
1967   if(yfitter.Eval() != 0) return kFALSE;
1968   par[0] = x0;
1969   par[1] = yfitter.GetParameter(0);
1970   par[2] = yfitter.GetParameter(1);
1971   return kTRUE;
1972 }
1973
1974 //____________________________________________________________________
1975 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
1976 {
1977 //
1978 // Global selection mechanism of tracksbased on cluster to fit residuals
1979 // The parameters are the same as used ni function FitTrack().
1980
1981   const Float_t kS(0.6), kM(0.2);
1982   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
1983   Float_t dy, dz, s, m;
1984   for(Int_t ip(0); ip<np; ip++){
1985     if(points[ip].GetClusterType()) continue;
1986     Float_t x0(points[ip].GetX())
1987           ,y0(param[1] + param[3] * (x0 - param[0]))
1988           ,z0(param[2] + param[4] * (x0 - param[0]));
1989     dy=points[ip].GetY() - y0; h.Fill(dy);
1990     dz=points[ip].GetZ() - z0;
1991   }
1992   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
1993   fg.SetParameter(1, 0.);
1994   fg.SetParameter(2, 2.e-2);
1995   h.Fit(&fg, "QN");
1996   m=fg.GetParameter(1); s=fg.GetParameter(2);
1997   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
1998   return kTRUE;
1999 }
2000
2001 //________________________________________________________
2002 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2003 {
2004   //
2005   // Get the most probable value and the full width half mean
2006   // of a Landau distribution
2007   //
2008
2009   const Float_t dx = 1.;
2010   mpv = f->GetParameter(1);
2011   Float_t fx, max = f->Eval(mpv);
2012
2013   xm = mpv - dx;
2014   while((fx = f->Eval(xm))>.5*max){
2015     if(fx>max){
2016       max = fx;
2017       mpv = xm;
2018     }
2019     xm -= dx;
2020   }
2021
2022   xM += 2*(mpv - xm);
2023   while((fx = f->Eval(xM))>.5*max) xM += dx;
2024 }
2025
2026
2027 // #include "TFile.h"
2028 // //________________________________________________________
2029 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
2030 // {
2031 //   if(!file){
2032 //     AliWarning("Use cluster position as in reconstruction.");
2033 //     SetLoadCorrection();
2034 //     return kTRUE;
2035 //   }
2036 //   TDirectory *cwd(gDirectory);
2037 //   TString fileList;
2038 //   FILE *filePtr = fopen(file, "rt");
2039 //   if(!filePtr){
2040 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
2041 //     SetLoadCorrection();
2042 //     return kFALSE;
2043 //   }
2044 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
2045 //   while(fileList.Gets(filePtr)){
2046 //     if(!TFile::Open(fileList.Data())) {
2047 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
2048 //       continue;
2049 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
2050 //
2051 //     TTree *tSys = (TTree*)gFile->Get("tSys");
2052 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
2053 //     tSys->Draw("det:t>>h2", "dx", "goff");
2054 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2055 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
2056 //     }
2057 //     h2->SetDirectory(cwd);
2058 //     gFile->Close();
2059 //   }
2060 //   cwd->cd();
2061 //
2062 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
2063 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2064 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
2065 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2066 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
2067 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
2068 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2069 //       printf("%03d|", idet);
2070 //       fprintf(fout, "    {");
2071 //       for(Int_t it(0); it<30; it++){
2072 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
2073 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
2074 //       }
2075 //       printf("\n");
2076 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
2077 //     }
2078 //     fprintf(fout, "  };\n");
2079 //   }
2080 //   SetLoadCorrection();
2081 //   return kTRUE;
2082 // }
2083
2084 //________________________________________________________
2085 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
2086   :fH(NULL)
2087   ,fNrebin(0)
2088   ,fRebinX(NULL)
2089   ,fRebinY(NULL)
2090 {
2091   // constructor
2092   memset(fAx, 0, 3*sizeof(Int_t));
2093   memset(fRange, 0, 4*sizeof(Float_t));
2094 }
2095
2096 //________________________________________________________
2097 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
2098 {
2099   // destructor
2100   if(fH) delete fH;
2101 }
2102
2103 //________________________________________________________
2104 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
2105 {
2106 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
2107   if(!aa[ix] || !aa[iy] || !aa[iz]) return;
2108   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
2109   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
2110     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
2111     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
2112     az->GetNbins(), az->GetXmin(), az->GetXmax());
2113   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
2114   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
2115 }
2116
2117 //________________________________________________________
2118 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
2119 {
2120 // increment bin with value "v" pointed by general coord in "bin"
2121   if(!fH) return;
2122   fH->AddBinContent(
2123         fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), v);
2124 }
2125
2126 //________________________________________________________
2127 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid)
2128 {
2129 // build the 2D projection and adjust binning
2130
2131   const Char_t *title[] = {"Mean", "#mu", "MPV"};
2132   if(!fH) return NULL;
2133   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
2134   TH2 *h2s = (TH2*)fH->Project3D("yx");
2135   Int_t irebin(0), dxBin(1), dyBin(1);
2136   while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
2137     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
2138     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
2139     irebin++;
2140   }
2141   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
2142   if(h2s) delete h2s;
2143
2144   // start projection
2145   TH1 *h(NULL);
2146   Float_t dz=(fRange[1]-fRange[1])/ncol;
2147   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
2148   TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
2149             Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""),
2150             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
2151   h2->SetContour(ncol);
2152   h2->GetZaxis()->CenterTitle();
2153   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
2154   //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny);
2155   for(Int_t iy(0); iy<ny; iy++){
2156     for(Int_t ix(0); ix<nx; ix++){
2157       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
2158       Int_t ne((Int_t)h->Integral());
2159       if(ne<nstat/2){
2160         h2->SetBinContent(ix+1, iy+1, -999);
2161         h2->SetBinError(ix+1, iy+1, 1.);
2162       }else{
2163         Float_t v(h->GetMean()), ve(h->GetRMS());
2164         if(mid==1){
2165           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
2166           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
2167           h->Fit(&fg, "WQ");
2168           v = fg.GetParameter(1); ve = fg.GetParameter(2);
2169         } else if (mid==2) {
2170           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
2171           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
2172           h->Fit(&fl, "WQ");
2173           v = fl.GetMaximumX(); ve = fl.GetParameter(2);
2174 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
2175           fgle.SetParameter(0, fl.GetParameter(0));
2176           fgle.SetParameter(1, fl.GetParameter(1));
2177           fgle.SetParameter(2, fl.GetParameter(2));
2178           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
2179           h->Fit(&fgle, "WQ");
2180           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
2181         }
2182         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
2183         else h2->SetBinContent(ix+1, iy+1, v);
2184         h2->SetBinError(ix+1, iy+1, ve);
2185       }
2186     }
2187   }
2188   if(h) delete h;
2189   return h2;
2190 }
2191
2192 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
2193 {
2194 // define rebinning strategy for this projection
2195   fNrebin = n;
2196   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
2197   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));
2198 }
2199
2200