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   AliDebug(2, Form("%s[%d]", H->GetName(), H->GetNdimensions()));
1111   return kTRUE;
1112 }
1113
1114 //________________________________________________________
1115 Bool_t AliTRDresolution::MakeProjectionTrackIn()
1116 {
1117 // Analyse track in
1118
1119   Int_t cidx = kTrackIn;
1120   if(fProj && fProj->At(cidx)) return kTRUE;
1121   if(!fContainer){
1122     AliError("Missing data container.");
1123     return kFALSE;
1124   }
1125   THnSparse *H(NULL);
1126   if(!(H = (THnSparse*)fContainer->At(cidx))){
1127     AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1128     return kFALSE;
1129   }
1130
1131   Int_t  coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1132   TAxis //*abc(H->GetAxis(kBC)),
1133         *aphi(H->GetAxis(kPhi)),
1134         *aeta(H->GetAxis(kEta)),
1135         //*as(H->GetAxis(kSpeciesChgRC)),
1136         //*apt(H->GetAxis(kPt)),
1137         *ay(H->GetAxis(kYrez)),
1138         *az(H->GetAxis(kZrez)),
1139         *ap(H->GetAxis(kPrez));
1140   Int_t neta(aeta->GetNbins()), nphi(aphi->GetNbins());
1141   TH3I *h3[fgNproj[cidx]];
1142   h3[0] = new TH3I("h3TI0", Form("r-#phi residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ay->GetTitle()),
1143             neta, aeta->GetXmin(), aeta->GetXmax(),
1144             nphi, aphi->GetXmin(), aphi->GetXmax(),
1145             ay->GetNbins(), ay->GetXmin(), ay->GetXmax());
1146   h3[1] = (TH3I*)h3[0]->Clone("h3TI1"); h3[1]->SetTitle("r-#phi residuals for pos tracks");
1147   h3[2] = new TH3I("h3TI2", Form("z residuals for row cross;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), az->GetTitle()),
1148             neta, aeta->GetXmin(), aeta->GetXmax(),
1149             nphi, aphi->GetXmin(), aphi->GetXmax(),
1150             az->GetNbins(), az->GetXmin(), az->GetXmax());
1151   h3[3] = new TH3I("h3TI3", Form("angular residuals for neg tracks;%s;%s;%s", aeta->GetTitle(), aphi->GetTitle(), ap->GetTitle()),
1152             neta, aeta->GetXmin(), aeta->GetXmax(),
1153             nphi, aphi->GetXmin(), aphi->GetXmax(),
1154             ap->GetNbins(), ap->GetXmin(), ap->GetXmax());
1155   h3[4] = (TH3I*)h3[3]->Clone("h3TI4"); h3[4]->SetTitle("angular residuals for pos tracks");
1156   for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1157     v = H->GetBinContent(ib, coord);
1158     if(v<1.) continue;
1159     if(coord[kBC]>1) continue; // bunch cross cut
1160     // species selection
1161     if(coord[kSpeciesChgRC]<6){
1162       h3[0]->AddBinContent(
1163         h3[0]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1164       h3[3]->AddBinContent(
1165         h3[3]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1166     } else if(coord[kSpeciesChgRC]==6) {
1167       h3[2]->AddBinContent(
1168         h3[2]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1169     } else if(coord[kSpeciesChgRC]>6) {
1170       h3[1]->AddBinContent(
1171         h3[1]->GetBin(coord[kEta], coord[kPhi], coord[kYrez]), v);
1172       h3[4]->AddBinContent(
1173         h3[4]->GetBin(coord[kEta], coord[kPhi], coord[kPrez]), v);
1174     }
1175   }
1176   if(!fProj){
1177     AliInfo("Building array of projections ...");
1178     fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1179   }
1180   TObjArray *arr(NULL);
1181   fProj->AddAt(arr = new TObjArray(fgNproj[cidx]), cidx);
1182
1183   TH2F *h2(NULL);
1184   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++){
1185     TAxis *ax(h3[iproj]->GetZaxis());
1186     Float_t zm(ax->GetXmin()/3.), zM(ax->GetXmax()/3.), dz=(zM-zm)/50;
1187     arr->AddAt(h2 = new TH2F(Form("h2TI%d", iproj),
1188             Form("%s;%s;%s;%s", h3[iproj]->GetTitle(), aeta->GetTitle(), aphi->GetTitle(), ax->GetTitle()),
1189             neta, aeta->GetXmin(), aeta->GetXmax(),
1190             nphi, aphi->GetXmin(), aphi->GetXmax()), iproj);
1191     h2->SetContour(50);
1192     h2->GetZaxis()->CenterTitle();
1193     h2->GetZaxis()->SetRangeUser(zm, zM);
1194     zm+=dz; zM-=dz;
1195     for(Int_t iphi(1); iphi<=nphi; iphi++){
1196       for(Int_t ieta(1); ieta<=neta; ieta++){
1197         TH1 *h = h3[iproj]->ProjectionZ(Form("hy%d", iproj), ieta, ieta, iphi, iphi);
1198         h2->SetBinContent(ieta, iphi, GetMeanWithBoundary(h, zm, zM, dz));
1199       }
1200     }
1201   }
1202 //   h2[5] = (TH2F*)h2[0]->Clone("h25");
1203 //   h2[5]->SetTitle("Systematic shift between neg/pos tracks");
1204 //   h2[5]->SetZTitle("#Delta(#Delta^{-}y - #Delta^{+}y) [cm]"); h2[5]->Reset();
1205 //   h2[6] = (TH2F*)h2[1]->Clone("h26");
1206 //   h2[6]->SetTitle("Average shift of pos&neg tracks");
1207 //   h2[6]->SetZTitle("<#Delta^{-}y, #Delta^{+}y> [cm]"); h2[6]->Reset();
1208 //   for(Int_t iphi(1); iphi<=nphi; iphi++){
1209 //     for(Int_t ieta(1); ieta<=neta; ieta++){
1210 //       Float_t neg = h2[0]->GetBinContent(ieta, iphi),
1211 //               pos = h2[1]->GetBinContent(ieta, iphi);
1212 //       if(neg<-100 || pos<-100){
1213 //         h2[5]->SetBinContent(ieta, iphi, -999.);
1214 //         h2[6]->SetBinContent(ieta, iphi, -999.);
1215 //       } else {
1216 //         h2[5]->SetBinContent(ieta, iphi, neg-pos);
1217 //         h2[6]->SetBinContent(ieta, iphi, 0.5*(neg+pos));
1218 //       }
1219 //     }
1220 //   }
1221
1222   for(Int_t iproj(0); iproj<fgNproj[cidx]; iproj++) delete h3[iproj];
1223   return kTRUE;
1224 }
1225
1226
1227
1228 //________________________________________________________
1229 Bool_t AliTRDresolution::PostProcess()
1230 {
1231 // Fit, Project, Combine, Extract values from the containers filled during execution
1232
1233   if (!fContainer) {
1234     AliError("ERROR: list not available");
1235     return kFALSE;
1236   }
1237
1238   // DEFINE MODELS
1239   // simple gauss
1240   TF1 fg("fGauss", "gaus", -.5, .5);
1241   // Landau for charge resolution
1242   TF1 fch("fClCh", "landau", 0., 1000.);
1243   // Landau for e+- pt resolution
1244   TF1 fpt("fPt", "landau", -0.1, 0.2);
1245
1246   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1247   // Clusters residuals
1248   if(!MakeProjectionCluster()) return kFALSE;
1249   fNRefFigures = 3;
1250   // Tracklet residual/pulls
1251   if(!MakeProjectionTracklet()) return kFALSE;
1252   fNRefFigures = 7;
1253   // TRDin residual/pulls
1254   if(!MakeProjectionTrackIn()) return kFALSE;
1255   fNRefFigures = 11;
1256   // TRDout residual/pulls
1257 //  if(!MakeProjectionTrackOut()) return kFALSE;
1258   fNRefFigures = 15;
1259
1260   if(!HasMCdata()) return kTRUE;
1261
1262
1263   //PROCESS MC RESIDUAL DISTRIBUTIONS
1264
1265   // CLUSTER Y RESOLUTION/PULLS
1266 //  if(!MakeProjectionClusterMC()) return kFALSE;
1267   fNRefFigures = 17;
1268
1269   // TRACKLET RESOLUTION/PULLS
1270 //  if(!MakeProjectionTrackletMC()) return kFALSE;
1271   fNRefFigures = 21;
1272
1273   // TRACK RESOLUTION/PULLS
1274 /*  Process3Darray(kMCtrack, 0, &fg, 1.e4);   // y
1275   Process3DlinkedArray(kMCtrack, 1, &fg);   // y PULL
1276   Process3Darray(kMCtrack, 2, &fg, 1.e4);   // z
1277   Process3Darray(kMCtrack, 3, &fg);         // z PULL
1278   Process2Darray(kMCtrack, 4, &fg, 1.e3);   // phi
1279   Process2Darray(kMCtrack, 5, &fg);         // snp PULL
1280   Process2Darray(kMCtrack, 6, &fg, 1.e3);   // theta
1281   Process2Darray(kMCtrack, 7, &fg);         // tgl PULL
1282   Process3Darray(kMCtrack, 8, &fg, 1.e2);   // pt resolution
1283   Process3Darray(kMCtrack, 9, &fg);         // 1/pt pulls
1284   Process3Darray(kMCtrack, 10, &fg, 1.e2);  // p resolution*/
1285 //  if(!MakeProjectionTrackMC(kMCtrack)) return kFALSE;
1286   fNRefFigures+=16;
1287
1288   // TRACK TRDin RESOLUTION/PULLS
1289 //  if(!MakeProjectionTrackMC(kMCtrackIn)) return kFALSE;
1290   fNRefFigures+=8;
1291
1292   // TRACK TRDout RESOLUTION/PULLS
1293 //  if(!MakeProjectionTrackMC(kMCtrackOut)) return kFALSE;
1294   fNRefFigures+=8;
1295
1296   return kTRUE;
1297 }
1298
1299
1300 //________________________________________________________
1301 void AliTRDresolution::Terminate(Option_t *opt)
1302 {
1303   AliTRDrecoTask::Terminate(opt);
1304   if(HasPostProcess()) PostProcess();
1305 }
1306
1307 //________________________________________________________
1308 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1309 {
1310 // Helper function to avoid duplication of code
1311 // Make first guesses on the fit parameters
1312
1313   // find the intial parameters of the fit !! (thanks George)
1314   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1315   Double_t sum = 0.;
1316   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1317   f->SetParLimits(0, 0., 3.*sum);
1318   f->SetParameter(0, .9*sum);
1319   Double_t rms = h->GetRMS();
1320   f->SetParLimits(1, -rms, rms);
1321   f->SetParameter(1, h->GetMean());
1322
1323   f->SetParLimits(2, 0., 2.*rms);
1324   f->SetParameter(2, rms);
1325   if(f->GetNpar() <= 4) return;
1326
1327   f->SetParLimits(3, 0., sum);
1328   f->SetParameter(3, .1*sum);
1329
1330   f->SetParLimits(4, -.3, .3);
1331   f->SetParameter(4, 0.);
1332
1333   f->SetParLimits(5, 0., 1.e2);
1334   f->SetParameter(5, 2.e-1);
1335 }
1336
1337 //________________________________________________________
1338 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1339 {
1340 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1341 //  - y reziduals/pulls
1342
1343   TObjArray *arr = new TObjArray(2);
1344   arr->SetName(name); arr->SetOwner();
1345   TH1 *h(NULL); char hname[100], htitle[300];
1346
1347   // tracklet resolution/pull in y direction
1348   snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1349   snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1350   Float_t rr = range<0.?fDyRange:range;
1351   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1352     Int_t nybins=50;
1353     if(expand) nybins*=2;
1354     h = new TH3S(hname, htitle,
1355                  48, -.48, .48,            // phi
1356                  60, -rr, rr,              // dy
1357                  nybins, -0.5, nybins-0.5);// segment
1358   } else h->Reset();
1359   arr->AddAt(h, 0);
1360   snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1361   snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y  / #sigma_{y};#Delta z  / #sigma_{z}", GetNameId(), name, "Detector");
1362   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1363     h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1364   } else h->Reset();
1365   arr->AddAt(h, 1);
1366
1367   return arr;
1368 }
1369
1370 //________________________________________________________
1371 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1372 {
1373 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1374 //  - y reziduals/pulls
1375 //  - z reziduals/pulls
1376 //  - phi reziduals
1377   TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
1378   arr->Expand(5);
1379   TH1 *h(NULL); char hname[100], htitle[300];
1380
1381   // tracklet resolution/pull in z direction
1382   snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1383   snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1384   if(!(h = (TH2S*)gROOT->FindObject(hname))){
1385     h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1386   } else h->Reset();
1387   arr->AddAt(h, 2);
1388   snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1389   snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z  / #sigma_{z};row cross", GetNameId(), name);
1390   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1391     h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1392     h->GetZaxis()->SetBinLabel(1, "no RC");
1393     h->GetZaxis()->SetBinLabel(2, "RC");
1394   } else h->Reset();
1395   arr->AddAt(h, 3);
1396
1397   // tracklet to track phi resolution
1398   snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1399   snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1400   Int_t nsgms=540;
1401   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1402     h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1403   } else h->Reset();
1404   arr->AddAt(h, 4);
1405
1406   return arr;
1407 }
1408
1409 //________________________________________________________
1410 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1411 {
1412 // Build performance histograms for AliExternalTrackParam.vs MC
1413 //  - y resolution/pulls
1414 //  - z resolution/pulls
1415 //  - phi resolution, snp pulls
1416 //  - theta resolution, tgl pulls
1417 //  - pt resolution, 1/pt pulls, p resolution
1418
1419   TObjArray *arr = BuildMonitorContainerTracklet(name);
1420   arr->Expand(11);
1421   TH1 *h(NULL); char hname[100], htitle[300];
1422   //TAxis *ax(NULL);
1423
1424   // snp pulls
1425   snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1426   snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp  / #sigma_{snp};entries", GetNameId(), name);
1427   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1428     h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1429   } else h->Reset();
1430   arr->AddAt(h, 5);
1431
1432   // theta resolution
1433   snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1434   snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1435   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1436     h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1437   } else h->Reset();
1438   arr->AddAt(h, 6);
1439   // tgl pulls
1440   snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1441   snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl  / #sigma_{tgl};entries", GetNameId(), name);
1442   if(!(h = (TH2I*)gROOT->FindObject(hname))){
1443     h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1444   } else h->Reset();
1445   arr->AddAt(h, 7);
1446
1447   const Int_t kNdpt(150);
1448   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1449   Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1450   Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1451   for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1452   for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1453   for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1454
1455   // Pt resolution
1456   snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1457   snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1458   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1459     h = new TH3S(hname, htitle,
1460                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1461     //ax = h->GetZaxis();
1462     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1463   } else h->Reset();
1464   arr->AddAt(h, 8);
1465   // 1/Pt pulls
1466   snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1467   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);
1468   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1469     h = new TH3S(hname, htitle,
1470                  kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1471     //ax = h->GetZaxis();
1472     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1473   } else h->Reset();
1474   arr->AddAt(h, 9);
1475   // P resolution
1476   snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1477   snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1478   if(!(h = (TH3S*)gROOT->FindObject(hname))){
1479     h = new TH3S(hname, htitle,
1480                  kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1481     //ax = h->GetZaxis();
1482     //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1483   } else h->Reset();
1484   arr->AddAt(h, 10);
1485
1486   return arr;
1487 }
1488
1489
1490 //________________________________________________________
1491 TObjArray* AliTRDresolution::Histos()
1492 {
1493   //
1494   // Define histograms
1495   //
1496
1497   if(fContainer) return fContainer;
1498
1499   fContainer  = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1500   THnSparse *H(NULL);
1501   const Int_t nhn(100); Char_t hn[nhn]; TString st;
1502
1503   //++++++++++++++++++++++
1504   // cluster to track residuals/pulls
1505   snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1506   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1507     const Char_t *clTitle[5/*kNdim*/] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], "chg*Q/vd/angle", fgkTitle[kYrez]/*, "#Deltax [cm]", "Q</Q", "#Phi^{*} - ExB [deg]"*/};
1508     const Int_t clNbins[5/*kNdim*/]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], 51, fgkNbins[kYrez]/*, 33, 10, 30*/};
1509     const Double_t clMin[5/*kNdim*/]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], -102., fgkMin[kYrez]/3./*, 0., 0.1, -30*/},
1510                    clMax[5/*kNdim*/]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], 102., fgkMax[kYrez]/3./*, 3.3, 2.1, 30*/};
1511     st = "cluster spatial&charge resolution;";
1512     for(Int_t idim(0); idim<5/*kNdim*/; idim++){ st += clTitle[idim]; st+=";";}
1513     H = new THnSparseI(hn, st.Data(), kNdim, clNbins, clMin, clMax);
1514   } else H->Reset();
1515   fContainer->AddAt(H, kCluster);
1516   //++++++++++++++++++++++
1517   // tracklet to TRD track
1518   snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1519   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1520     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]*/};
1521     const Int_t trNbins[kNdim+1]   = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kSpeciesChgRC], fgkNbins[kYrez], fgkNbins[kZrez], fgkNbins[kPrez], 30/*, fgkNbins[kPt]*/};
1522     const Double_t trMin[kNdim+1]  = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kSpeciesChgRC], fgkMin[kYrez], fgkMin[kZrez], fgkMin[kPrez], 0./*, fgkMin[kPt]*/},
1523                    trMax[kNdim+1]  = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kSpeciesChgRC], fgkMax[kYrez], fgkMax[kZrez], fgkMax[kPrez], 3000./*, fgkMax[kPt]*/};
1524     st = "tracklet spatial&charge resolution;";
1525     for(Int_t idim(0); idim<kNdim+1; idim++){ st += trTitle[idim]; st+=";";}
1526     H = new THnSparseI(hn, st.Data(), kNdim+1, trNbins, trMin, trMax);
1527   } else H->Reset();
1528   fContainer->AddAt(H, kTracklet);
1529   //++++++++++++++++++++++
1530   // tracklet to TRDin
1531   snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1532   if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1533     st = "r-#phi/z/angular residuals @ TRD entry;";
1534     for(Int_t idim(0); idim<kNdim; idim++){ st += fgkTitle[idim]; st+=";";}
1535     H = new THnSparseI(hn, st.Data(), kNdim, fgkNbins, fgkMin, fgkMax);
1536   } else H->Reset();
1537   fContainer->AddAt(H, kTrackIn);
1538   // tracklet to TRDout
1539   fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1540
1541
1542   // Resolution histos
1543   if(!HasMCdata()) return fContainer;
1544
1545   // cluster resolution
1546   fContainer->AddAt(BuildMonitorContainerCluster("MCcl"),  kMCcluster);
1547
1548   // tracklet resolution
1549   fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1550
1551   // track resolution
1552   TObjArray *arr(NULL);
1553   fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1554   arr->SetName("MCtrk");
1555   for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1556
1557   // TRDin TRACK RESOLUTION
1558   fContainer->AddAt(H, kMCtrackIn);
1559
1560   // TRDout TRACK RESOLUTION
1561   fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1562
1563   return fContainer;
1564 }
1565
1566 //________________________________________________________
1567 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1568 {
1569 // Robust function to process sigma/mean for 2D plot dy(x)
1570 // For each x bin a gauss fit is performed on the y projection and the range
1571 // with the minimum chi2/ndf is choosen
1572
1573   if(!h2) {
1574     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
1575     return kFALSE;
1576   }
1577   if(!Int_t(h2->GetEntries())){
1578     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
1579     return kFALSE;
1580   }
1581   if(!g || !g[0]|| !g[1]) {
1582     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
1583     return kFALSE;
1584   }
1585   // prepare
1586   TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
1587   Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
1588   TF1 f("f", "gaus", ymin, ymax);
1589   Int_t n(0);
1590   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1591   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1592   TH1D *h(NULL);
1593   if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
1594   Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
1595
1596
1597   // do actual loop
1598   Float_t chi2OverNdf(0.);
1599   for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
1600     x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
1601     ymin = ay->GetXmin(); ymax = ay->GetXmax();
1602
1603     h = h2->ProjectionY("py", ix, ix);
1604     if((n=(Int_t)h->GetEntries())<stat){
1605       if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
1606       continue;
1607     }
1608     // looking for a first order mean value
1609     f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
1610     h->Fit(&f, "QNW");
1611     chi2OverNdf = f.GetChisquare()/f.GetNDF();
1612     printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
1613     y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
1614     ey  = f.GetParError(1);
1615     sy = f.GetParameter(2); esy = f.GetParError(2);
1616 //     // looking for the best chi2/ndf value
1617 //     while(ymin<y0 && ymax>y1){
1618 //       f.SetParameter(1, y);
1619 //       f.SetParameter(2, sy);
1620 //       h->Fit(&f, "QNW", "", y0, y1);
1621 //       printf("   range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
1622 //       if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
1623 //         chi2OverNdf = f.GetChisquare()/f.GetNDF();
1624 //         y  = f.GetParameter(1); ey  = f.GetParError(1);
1625 //         sy = f.GetParameter(2); esy = f.GetParError(2);
1626 //         printf("    set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
1627 //       }
1628 //       y0-=dy; y1+=dy;
1629 //     }
1630     g[0]->SetPoint(np, x, y);
1631     g[0]->SetPointError(np, ex, ey);
1632     g[1]->SetPoint(np, x, sy);
1633     g[1]->SetPointError(np, ex, esy);
1634     np++;
1635   }
1636   return kTRUE;
1637 }
1638
1639
1640 //________________________________________________________
1641 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1642 {
1643   //
1644   // Do the processing
1645   //
1646
1647   Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
1648   Int_t n = 0;
1649   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1650   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1651   if(Int_t(h2->GetEntries())){
1652     AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
1653   } else {
1654     AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
1655     fIdxPlot++;
1656     return kTRUE;
1657   }
1658
1659   const Int_t kINTEGRAL=1;
1660   for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
1661     Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
1662     Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
1663     TH1D *h = h2->ProjectionY(pn, abin, bbin);
1664     if((n=(Int_t)h->GetEntries())<150){
1665       AliDebug(4, Form("  x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
1666       continue;
1667     }
1668     h->Fit(f, "QN");
1669     Int_t ip = g[0]->GetN();
1670     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)));
1671     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1672     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1673     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1674     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1675 /*
1676     g[0]->SetPoint(ip, x, k*h->GetMean());
1677     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1678     g[1]->SetPoint(ip, x, k*h->GetRMS());
1679     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1680   }
1681   fIdxPlot++;
1682   return kTRUE;
1683 }
1684
1685
1686 //____________________________________________________________________
1687 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
1688 {
1689 //
1690 // Fit track with a staight line using the "np" clusters stored in the array "points".
1691 // The following particularities are stored in the clusters from points:
1692 //   1. pad tilt as cluster charge
1693 //   2. pad row cross or vertex constrain as fake cluster with cluster type 1
1694 // The parameters of the straight line fit are stored in the array "param" in the following order :
1695 //     param[0] - x0 reference radial position
1696 //     param[1] - y0 reference r-phi position @ x0
1697 //     param[2] - z0 reference z position @ x0
1698 //     param[3] - slope dy/dx
1699 //     param[4] - slope dz/dx
1700 //
1701 // Attention :
1702 // Function should be used to refit tracks for B=0T
1703 //
1704
1705   if(np<40){
1706     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
1707     return kFALSE;
1708   }
1709   TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
1710
1711   Double_t x0(0.);
1712   for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
1713   x0/=Float_t(np);
1714
1715   Double_t x, y, z, dx, tilt(0.);
1716   for(Int_t ip(0); ip<np; ip++){
1717     x = points[ip].GetX(); z = points[ip].GetZ();
1718     dx = x - x0;
1719     zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
1720   }
1721   if(zfitter.Eval() != 0) return kFALSE;
1722
1723   Double_t z0    = zfitter.GetParameter(0);
1724   Double_t dzdx  = zfitter.GetParameter(1);
1725   for(Int_t ip(0); ip<np; ip++){
1726     if(points[ip].GetClusterType()) continue;
1727     x    = points[ip].GetX();
1728     dx   = x - x0;
1729     y    = points[ip].GetY();
1730     z    = points[ip].GetZ();
1731     tilt = points[ip].GetCharge();
1732     y -= tilt*(-dzdx*dx + z - z0);
1733     Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
1734     yfitter.AddPoint(&dx, y, 1.);
1735   }
1736   if(yfitter.Eval() != 0) return kFALSE;
1737   Double_t y0   = yfitter.GetParameter(0);
1738   Double_t dydx = yfitter.GetParameter(1);
1739
1740   param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
1741   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);
1742   return kTRUE;
1743 }
1744
1745 //____________________________________________________________________
1746 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
1747 {
1748 //
1749 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
1750 // See function FitTrack for the data stored in the "clusters" array
1751
1752 // The parameters of the straight line fit are stored in the array "param" in the following order :
1753 //     par[0] - x0 reference radial position
1754 //     par[1] - y0 reference r-phi position @ x0
1755 //     par[2] - slope dy/dx
1756 //
1757 // Attention :
1758 // Function should be used to refit tracks for B=0T
1759 //
1760
1761   TLinearFitter yfitter(2, "pol1");
1762
1763   // grep data for tracklet
1764   Double_t x0(0.), x[60], y[60], dy[60];
1765   Int_t nly(0);
1766   for(Int_t ip(0); ip<np; ip++){
1767     if(points[ip].GetClusterType()) continue;
1768     if(points[ip].GetVolumeID() != ly) continue;
1769     Float_t xt(points[ip].GetX())
1770            ,yt(param[1] + param[3] * (xt - param[0]));
1771     x[nly] = xt;
1772     y[nly] = points[ip].GetY();
1773     dy[nly]= y[nly]-yt;
1774     x0    += xt;
1775     nly++;
1776   }
1777   if(nly<10){
1778     if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
1779     return kFALSE;
1780   }
1781   // set radial reference for fit
1782   x0 /= Float_t(nly);
1783
1784   // find tracklet core
1785   Double_t mean(0.), sig(1.e3);
1786   AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
1787
1788   // simple cluster error parameterization
1789   Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
1790
1791   // fit tracklet core
1792   for(Int_t jly(0); jly<nly; jly++){
1793     if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
1794     Double_t dx(x[jly]-x0);
1795     yfitter.AddPoint(&dx, y[jly], 1.);
1796   }
1797   if(yfitter.Eval() != 0) return kFALSE;
1798   par[0] = x0;
1799   par[1] = yfitter.GetParameter(0);
1800   par[2] = yfitter.GetParameter(1);
1801   return kTRUE;
1802 }
1803
1804 //____________________________________________________________________
1805 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
1806 {
1807 //
1808 // Global selection mechanism of tracksbased on cluster to fit residuals
1809 // The parameters are the same as used ni function FitTrack().
1810
1811   const Float_t kS(0.6), kM(0.2);
1812   TH1S h("h1", "", 100, -5.*kS, 5.*kS);
1813   Float_t dy, dz, s, m;
1814   for(Int_t ip(0); ip<np; ip++){
1815     if(points[ip].GetClusterType()) continue;
1816     Float_t x0(points[ip].GetX())
1817           ,y0(param[1] + param[3] * (x0 - param[0]))
1818           ,z0(param[2] + param[4] * (x0 - param[0]));
1819     dy=points[ip].GetY() - y0; h.Fill(dy);
1820     dz=points[ip].GetZ() - z0;
1821   }
1822   TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
1823   fg.SetParameter(1, 0.);
1824   fg.SetParameter(2, 2.e-2);
1825   h.Fit(&fg, "QN");
1826   m=fg.GetParameter(1); s=fg.GetParameter(2);
1827   if(s>kS || TMath::Abs(m)>kM) return kFALSE;
1828   return kTRUE;
1829 }
1830
1831 //________________________________________________________
1832 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
1833 {
1834   //
1835   // Get the most probable value and the full width half mean
1836   // of a Landau distribution
1837   //
1838
1839   const Float_t dx = 1.;
1840   mpv = f->GetParameter(1);
1841   Float_t fx, max = f->Eval(mpv);
1842
1843   xm = mpv - dx;
1844   while((fx = f->Eval(xm))>.5*max){
1845     if(fx>max){
1846       max = fx;
1847       mpv = xm;
1848     }
1849     xm -= dx;
1850   }
1851
1852   xM += 2*(mpv - xm);
1853   while((fx = f->Eval(xM))>.5*max) xM += dx;
1854 }
1855
1856
1857 // #include "TFile.h"
1858 // //________________________________________________________
1859 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
1860 // {
1861 //   if(!file){
1862 //     AliWarning("Use cluster position as in reconstruction.");
1863 //     SetLoadCorrection();
1864 //     return kTRUE;
1865 //   }
1866 //   TDirectory *cwd(gDirectory);
1867 //   TString fileList;
1868 //   FILE *filePtr = fopen(file, "rt");
1869 //   if(!filePtr){
1870 //     AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
1871 //     SetLoadCorrection();
1872 //     return kFALSE;
1873 //   }
1874 //   TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
1875 //   while(fileList.Gets(filePtr)){
1876 //     if(!TFile::Open(fileList.Data())) {
1877 //       AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
1878 //       continue;
1879 //     } else AliInfo(Form("\"%s\"", fileList.Data()));
1880 //
1881 //     TTree *tSys = (TTree*)gFile->Get("tSys");
1882 //     h2->SetDirectory(gDirectory); h2->Reset("ICE");
1883 //     tSys->Draw("det:t>>h2", "dx", "goff");
1884 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1885 //       for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
1886 //     }
1887 //     h2->SetDirectory(cwd);
1888 //     gFile->Close();
1889 //   }
1890 //   cwd->cd();
1891 //
1892 //   if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
1893 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1894 //     printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
1895 //     for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
1896 //     FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
1897 //     fprintf(fout, "  static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
1898 //     for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
1899 //       printf("%03d|", idet);
1900 //       fprintf(fout, "    {");
1901 //       for(Int_t it(0); it<30; it++){
1902 //         printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
1903 //         fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
1904 //       }
1905 //       printf("\n");
1906 //       fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
1907 //     }
1908 //     fprintf(fout, "  };\n");
1909 //   }
1910 //   SetLoadCorrection();
1911 //   return kTRUE;
1912 // }