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