TRD/qaRec/AliTRDtrackingResolution.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: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Reconstruction QA                                                     //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  24 //                                                                        //
  25 ////////////////////////////////////////////////////////////////////////////
  26
  27 #include <cstring>
  28
  29
  30 #include <TObjArray.h>
  31 #include <TH2.h>
  32 #include <TH1.h>
  33 #include <TF1.h>
  34 #include <TProfile.h>
  35 #include <TGraphErrors.h>
  36 #include <TMath.h>
  37 #include "TTreeStream.h"
  38 #include "TGeoManager.h"
  39
  40 #include "AliAnalysisManager.h"
  41 #include "AliTrackReference.h"
  42 #include "AliTrackPointArray.h"
  43 #include "AliCDBManager.h"
  44
  45 #include "AliTRDSimParam.h"
  46 #include "AliTRDgeometry.h"
  47 #include "AliTRDpadPlane.h"
  48 #include "AliTRDcluster.h"
  49 #include "AliTRDseedV1.h"
  50 #include "AliTRDtrackV1.h"
  51 #include "AliTRDtrackerV1.h"
  52 #include "AliTRDReconstructor.h"
  53 #include "AliTRDrecoParam.h"
  54
  55 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  56 #include "AliTRDtrackingResolution.h"
  57
  58 ClassImp(AliTRDtrackingResolution)
  59
  60 //________________________________________________________
  61 AliTRDtrackingResolution::AliTRDtrackingResolution()
  62   :AliTRDrecoTask("Resolution", "Tracking Resolution")
  63   ,fReconstructor(0x0)
  64   ,fGeo(0x0)
  65 {
  66   fReconstructor = new AliTRDReconstructor();
  67   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  68   fGeo = new AliTRDgeometry();
  69 }
  70
  71 //________________________________________________________
  72 AliTRDtrackingResolution::~AliTRDtrackingResolution()
  73 {
  74   delete fGeo;
  75   delete fReconstructor;
  76   if(gGeoManager) delete gGeoManager;
  77 }
  78
  79
  80 //________________________________________________________
  81 void AliTRDtrackingResolution::CreateOutputObjects()
  82 {
  83   // spatial resolution
  84   OpenFile(0, "RECREATE");
  85
  86   fContainer = Histos();
  87
  88   // cluster to tracklet residuals [2]
  89   fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -21., 21., 100, -.5, .5), kClusterYResidual);
  90   // tracklet to Riemann fit residuals [2]
  91   fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
  92   fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
  93   fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
  94   fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
  95
  96   // Resolution histos
  97   if(HasMCdata()){
  98     // tracklet resolution [0]
  99     fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 21, -21., 21., 100, -.5, .5), kTrackletYResolution);
 100     // tracklet angular resolution [1]
 101     fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackletAngleResolution);
 102
 103     // Riemann track resolution [y, z, angular]
 104     fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
 105     fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
 106     fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
 107
 108     // Kalman track resolution [y, z, angular]
 109     fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
 110     fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
 111     fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
 112   }
 113
 114   // CREATE GRAPHS for DISPLAY
 115
 116   // define iterator over graphs
 117   Int_t jgraph = (Int_t)kGraphStart;
 118   TH2I *h2 = (TH2I *)(fContainer->At(kClusterYResidual));
 119   // clusters tracklet residuals (mean-phi)
 120   TH1 *h = new TH1I("h", "", 100, -40., 40.);
 121   h->GetXaxis()->SetTitle("#Phi [deg]");
 122   h->GetYaxis()->SetTitle("Clusters Residuals : #sigma/#mu [mm]");
 123   h->GetYaxis()->SetRangeUser(-.05, 1.);
 124   fContainer->AddAt(h, jgraph++);
 125
 126   TGraphErrors *g = new TGraphErrors(h2->GetNbinsX());
 127   g->SetLineColor(kGreen);
 128   g->SetMarkerStyle(22);
 129   g->SetMarkerColor(kGreen);
 130   g->SetNameTitle("clm", "Residuals Clusters-Tracklet Mean");
 131   fContainer->AddAt(g, jgraph++);
 132
 133   // clusters tracklet residuals (sigma-phi)
 134   g = new TGraphErrors(h2->GetNbinsX());
 135   g->SetLineColor(kRed);
 136   g->SetMarkerStyle(23);
 137   g->SetMarkerColor(kRed);
 138   g->SetNameTitle("cls", "Residuals Clusters-Tracklet Sigma");
 139   fContainer->AddAt(g, jgraph++);
 140
 141   if(HasMCdata()){
 142     // tracklet y resolution
 143     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 144     h = new TH1I("h", "", 100, -40., 40.);
 145     h->GetXaxis()->SetTitle("#Phi [deg]");
 146     h->GetYaxis()->SetTitle("Tracklet Resolution : #sigma/#mu [mm]");
 147     h->GetYaxis()->SetRangeUser(-.05, 1.);
 148     fContainer->AddAt(h, jgraph++);
 149
 150     g = new TGraphErrors(h2->GetNbinsX());
 151     g->SetLineColor(kGreen);
 152     g->SetMarkerStyle(22);
 153     g->SetMarkerColor(kGreen);
 154     g->SetNameTitle("trkltym", "Resolution Tracklet Y Mean");
 155     fContainer->AddAt(g, jgraph++);
 156     g = new TGraphErrors(h2->GetNbinsX());
 157     g->SetLineColor(kRed);
 158     g->SetMarkerStyle(22);
 159     g->SetMarkerColor(kRed);
 160     g->SetNameTitle("trkltys", "Resolution Tracklet Y Sigma");
 161     fContainer->AddAt(g, jgraph++);
 162
 163     // tracklet phi resolution
 164     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 165     h = new TH1I("h", "", 100, -40., 40.);
 166     h->GetXaxis()->SetTitle("#Phi [deg]");
 167     h->GetYaxis()->SetTitle("Tracklet Angular Resolution : #sigma/#mu [deg]");
 168     h->GetYaxis()->SetRangeUser(-.05, .2);
 169     fContainer->AddAt(h, jgraph++);
 170
 171     g = new TGraphErrors(h2->GetNbinsX());
 172     g->SetLineColor(kGreen);
 173     g->SetMarkerStyle(22);
 174     g->SetMarkerColor(kGreen);
 175     g->SetNameTitle("trkltam", "Resolution Tracklet Y Mean");
 176     fContainer->AddAt(g, jgraph++);
 177     g = new TGraphErrors(h2->GetNbinsX());
 178     g->SetLineColor(kRed);
 179     g->SetMarkerStyle(22);
 180     g->SetMarkerColor(kRed);
 181     g->SetNameTitle("trkltas", "Angle Resolution Tracklet Sigma");
 182     fContainer->AddAt(g, jgraph++);
 183   }
 184 }
 185
 186 //________________________________________________________
 187 void AliTRDtrackingResolution::Exec(Option_t *)
 188 {
 189   // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
 190   // angular Resolution: res = Tracklet angle - TrackRef Angle
 191
 192   Int_t nTrackInfos = fTracks->GetEntriesFast();
 193   if(fDebugLevel>=2){
 194     printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
 195   }
 196
 197   Int_t pdg;
 198   Double_t p, dy, dphi, dymc, dzmc, dphimc;
 199   Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
 200   Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
 201
 202   AliTRDpadPlane *pp = 0x0;
 203   AliTrackPoint tr[kNLayers], tk[kNLayers];
 204   AliExternalTrackParam *fOp = 0x0;
 205   AliTRDtrackV1 *fTrack = 0x0;
 206   AliTRDtrackInfo *fInfo = 0x0;
 207   for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
 208     // check if ESD and MC-Information are available
 209     if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
 210     if(!(fTrack = fInfo->GetTRDtrack())) continue;
 211     if(!(fOp = fInfo->GetOuterParam())) continue;
 212     pdg = fInfo->GetPDG();
 213
 214     if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
 215
 216     p = fOp->P();
 217     Int_t npts = 0;
 218     memset(fP, 0, kNLayers*sizeof(Float_t));
 219     memset(fX, 0, kNLayers*sizeof(Float_t));
 220     memset(fY, 0, kNLayers*sizeof(Float_t));
 221     memset(fZ, 0, kNLayers*sizeof(Float_t));
 222     memset(fPhi, 0, kNLayers*sizeof(Float_t));
 223     memset(fTheta, 0, kNLayers*sizeof(Float_t));
 224     memset(fLayerMap, 0, kNLayers*sizeof(Bool_t));
 225     memset(fMCMap, 0, kNLayers*sizeof(Bool_t));
 226     AliTRDseedV1 *fTracklet = 0x0;
 227     for(Int_t iplane = 0; iplane < kNLayers; iplane++){
 228       if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
 229       if(!fTracklet->IsOK()) continue;
 230
 231       // Book point arrays
 232       fLayerMap[iplane] = kTRUE;
 233       tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
 234       tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
 235       npts++;
 236
 237       if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
 238
 239       // define reference values
 240       fP[iplane]   = p;
 241       fX[iplane]   = fTracklet->GetX0();
 242       fY[iplane]   = fTracklet->GetYref(0);
 243       fZ[iplane]   = fTracklet->GetZref(0);
 244       fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1));
 245       fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1));
 246
 247
 248       // RESOLUTION (compare to MC)
 249       if(HasMCdata()){
 250         if(fInfo->GetNTrackRefs() >= 2){
 251           Double_t pmc, ymc, zmc, phiMC, thetaMC;
 252           if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){
 253             fMCMap[iplane] = kTRUE;
 254             fP[iplane]     = pmc;
 255             fY[iplane]     = ymc;
 256             fZ[iplane]     = zmc;
 257             fPhi[iplane]   = phiMC;
 258             fTheta[iplane] = thetaMC;
 259           }
 260         }
 261       }
 262       Float_t phi   = fPhi[iplane]*TMath::RadToDeg();
 263       //Float_t theta = fTheta[iplane]*TMath::RadToDeg();
 264
 265       // Do clusters residuals
 266       if(!fTracklet->Fit(kFALSE)) continue;
 267       AliTRDcluster *c = 0x0;
 268       for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 269         if(!(c = fTracklet->GetClusters(ic))) continue;
 270
 271         dy = fTracklet->GetYat(c->GetX()) - c->GetY();
 272         ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
 273
 274         if(fDebugLevel>=2){
 275           Float_t q = c->GetQ();
 276           // Get z-position with respect to anode wire
 277           AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 278           Int_t det = c->GetDetector();
 279           Float_t x = c->GetX();
 280           Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
 281           Int_t stack = fGeo->GetStack(det);
 282           pp = fGeo->GetPadPlane(iplane, stack);
 283           Float_t row0 = pp->GetRow0();
 284           Float_t d  =  row0 - z + simParam->GetAnodeWireOffset();
 285           d -= ((Int_t)(2 * d)) / 2.0;
 286           //if (d > 0.25) d  = 0.5 - d;
 287
 288           (*fDebugStream) << "ResidualClusters"
 289             << "ly="   << iplane
 290             << "stk="  << stack
 291             << "pdg="  << pdg
 292             << "phi="  << fPhi[iplane]
 293             << "tht="  << fTheta[iplane]
 294             << "q="    << q
 295             << "x="    << x
 296             << "z="    << z
 297             << "d="    << d
 298             << "dy="   << dy
 299             << "\n";
 300         }
 301       }
 302       pp = 0x0;
 303     }
 304
 305
 306     // this protection we might drop TODO
 307     if(fTrack->GetNumberOfTracklets() < 6) continue;
 308
 309     AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
 310     Int_t iref = 0;
 311     for(Int_t ip=0; ip<kNLayers; ip++){
 312       if(!fLayerMap[ip]) continue;
 313       fTracklet = fTrack->GetTracklet(ip);
 314       // recalculate fit based on the new tilt correction
 315       fTracklet->Fit();
 316
 317       dy = fTracklet->GetYfit(0) - tr[iref].GetY();
 318       ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
 319
 320       dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
 321       ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
 322
 323       if(HasMCdata()){
 324         dymc = fY[ip] - tr[iref].GetY();
 325         ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
 326
 327         dzmc = fZ[ip] - tr[iref].GetZ();
 328         ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
 329
 330         dphimc = fPhi[ip] - fTracklet->GetYfit(1);
 331         ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
 332       }
 333
 334       iref++;
 335
 336       if(fDebugLevel>=2){
 337         (*fDebugStream) << "RiemannTrack"
 338           << "ly="    << ip
 339           << "mc="    << fMCMap[ip]
 340           << "p="     << fP[ip]
 341           << "phi="   << fPhi[ip]
 342           << "tht="   << fTheta[ip]
 343           << "dy="    << dy
 344           << "dphi="  << dphi
 345           << "dymc="  << dymc
 346           << "dzmc="  << dzmc
 347           << "dphimc="<< dphimc
 348           << "\n";
 349       }
 350     }
 351
 352 //  if(!gGeoManager) TGeoManager::Import("geometry.root");
 353 //     AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
 354 //     for(Int_t ip=0; ip<nc; ip++){
 355 //       dy = cl[ip].GetY() - tr[ip].GetY();
 356 //      ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
 357 //       dz = cl[ip].GetZ() - tr[ip].GetZ();
 358 //       if(fDebugLevel>=1){
 359 //         (*fDebugStream) << "KalmanTrack"
 360 //           << "dy="             << dy
 361 //           << "dz="             << dz
 362 // /*          << "phi="                        << phi
 363 //           << "theta="                << theta
 364 //           << "dphi="         << dphi*/
 365 //           << "\n";
 366 //       }
 367 //     }
 368
 369
 370   }
 371   PostData(0, fContainer);
 372 }
 373
 374 //________________________________________________________
 375 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt)
 376 {
 377   //sprintf(opt, "pl");
 378   switch(ifig){
 379   case 0:
 380     first = (Int_t)kGraphStart; last = first+3;
 381     break;
 382   case 1:
 383     first = (Int_t)kGraphStart+3; last = first+3;
 384     break;
 385   case 2:
 386     first = (Int_t)kGraphStart+6; last = first+3;
 387     break;
 388   default:
 389     first = (Int_t)kGraphStart; last = first;
 390     break;
 391   }
 392 }
 393
 394
 395 //________________________________________________________
 396 Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &p, Double_t &ymc, Double_t &zmc, Double_t &phi, Double_t &theta)
 397 {
 398
 399   AliTrackReference *fTrackRefs[2] = {0x0, 0x0};
 400   Float_t x0 = tracklet->GetX0();
 401
 402   // check for 2 track ref where the radial position has a distance less than 3.7mm
 403   Int_t nFound = 0;
 404   for(Int_t itr = 0; itr < AliTRDtrackInfo::kNTrackRefs; itr++){
 405     if(!(fTrackRefs[nFound] = fInfo->GetTrackRef(itr))) break;
 406
 407     if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, fTrackRefs[nFound]->LocalX());
 408     if(TMath::Abs(x0 - fTrackRefs[nFound]->LocalX()) > 3.7) continue;
 409     nFound++;
 410     if(nFound == 2) break;
 411   }
 412   if(nFound < 2){
 413     if(fDebugLevel>=3) printf("\t\tMissing track ref x0[%6.3f] ly[%d] nref[%d]\n", x0, tracklet->GetPlane(), fInfo->GetMCinfo()->GetNTrackRefs());
 414     return kFALSE;
 415   }
 416   // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
 417
 418
 419   // RESOLUTION
 420   Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
 421   if(dx <= 0.){
 422     if(fDebugLevel>=3) printf("\t\tTrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
 423     return kFALSE;
 424   }
 425   Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
 426   Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
 427   Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
 428   ymc =  fTrackRefs[1]->LocalY() - dydx*dx0;
 429   zmc =  fTrackRefs[1]->Z() - dzdx*dx0;
 430
 431   // recalculate tracklet based on the MC info
 432   tracklet->SetZref(0, zmc);
 433   tracklet->SetZref(1, -dzdx); // TODO
 434   tracklet->Fit();
 435   Double_t dy = tracklet->GetYfit(0) - ymc;
 436   Double_t dz = tracklet->GetZfit(0) - zmc;
 437
 438   //res_y *= 100; // in mm
 439   p = fTrackRefs[0]->P();
 440
 441   phi   = TMath::ATan(dydx);
 442   theta = TMath::ATan(dzdx);
 443   Double_t dphi   = TMath::ATan(tracklet->GetYfit(1)) - phi;
 444   if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
 445
 446   // Fill Histograms
 447   if(TMath::Abs(dx-3.7)<1.E-3){
 448     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(phi*TMath::RadToDeg(), dy);
 449     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
 450   }
 451   // Fill Debug Tree
 452   if(fDebugLevel>=2){
 453     Int_t iplane = tracklet->GetPlane();
 454     Int_t pdg = fInfo->GetPDG();
 455     (*fDebugStream) << "ResolutionTrklt"
 456       << "ly="            << iplane
 457       << "pdg="     << pdg
 458       << "p="       << p
 459       << "phi="                 << phi
 460       << "tht="           << theta
 461       << "ymc="     << ymc
 462       << "zmc="     << zmc
 463       << "dx="      << dx
 464       << "dy="            << dy
 465       << "dz="            << dz
 466       << "dphi="                << dphi
 467       << "\n";
 468
 469     Int_t det, stk;
 470     Float_t z0 = 0.;
 471     AliTRDpadPlane *pp = 0x0;
 472     AliTRDcluster *c = 0x0;
 473     tracklet->ResetClusterIter(kFALSE);
 474     while((c = tracklet->PrevCluster())){
 475       if(!pp){
 476         det = c->GetDetector();
 477         stk = fGeo->GetStack(det);
 478         pp = fGeo->GetPadPlane(iplane, stk);
 479         z0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
 480       }
 481       Float_t xc = c->GetX();
 482       dx = tracklet->GetX0() - xc;
 483       Float_t yt = ymc - dx*dydx;
 484       Float_t zt = zmc - dx*dzdx;
 485       Float_t yc = c->GetY() - TMath::Tan(tracklet->GetTilt()) * (c->GetZ() - zt);
 486       dy = yt - yc;
 487       Float_t d = z0 - zt;
 488       d -= ((Int_t)(2 * d)) / 2.0;
 489       (*fDebugStream) << "ResolutionClstr"
 490         << "det="                 << det
 491         << "ly="                  << iplane
 492         << "stk="                 << stk
 493         << "pdg="     << pdg
 494         << "p="       << p
 495         << "phi="                       << phi
 496         << "tht="                 << theta
 497         << "xc="       << xc
 498         << "yc="       << yc
 499         << "yt="       << yt
 500         << "zt="       << zt
 501         << "z0="       << z0
 502         << "d="       << d
 503         << "dy="                  << dy
 504         << "\n";
 505     }
 506   }
 507
 508   return kTRUE;
 509 }
 510
 511 //________________________________________________________
 512 Bool_t AliTRDtrackingResolution::PostProcess()
 513 {
 514   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 515   fNRefFigures = 0;
 516   if (!fContainer) {
 517     Printf("ERROR: list not available");
 518     return kFALSE;
 519   }
 520
 521   TH2I *h2 = 0x0;
 522   TH1D *h = 0x0;
 523   TGraphErrors *gm = 0x0, *gs = 0x0;
 524   TF1 f("f1", "gaus", -.5, .5);
 525   // define iterator over graphs
 526   Int_t jgraph = (Int_t)kGraphStart;
 527
 528   //PROCESS RESIDUAL DISTRIBUTIONS
 529
 530   // Clusters residuals
 531   // define model
 532   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 533   h2 = (TH2I *)(fContainer->At(kClusterYResidual));
 534   jgraph++; //skip the frame histo
 535   gm = (TGraphErrors*)fContainer->At(jgraph++);
 536   gs = (TGraphErrors*)fContainer->At(jgraph++);
 537   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 538     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 539     Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
 540     h = h2->ProjectionY("py", ibin, ibin);
 541     Fit(h, &fc);
 542     gm->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(1));
 543     gm->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(1));
 544     gs->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(2));
 545     gs->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(2));
 546   }
 547   fNRefFigures++;
 548
 549
 550   //PROCESS RESOLUTION DISTRIBUTIONS
 551   if(HasMCdata()){
 552     // tracklet y resolution
 553     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 554     jgraph++; //skip the frame histo
 555     gm = (TGraphErrors*)fContainer->At(jgraph++);
 556     gs = (TGraphErrors*)fContainer->At(jgraph++);
 557     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 558       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 559       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 560       h = h2->ProjectionY("py", iphi, iphi);
 561       Fit(h, &fc);
 562       Int_t jphi = iphi -1;
 563       gm->SetPoint(jphi, phi, 10.*f.GetParameter(1));
 564       gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
 565       gs->SetPoint(jphi, phi, 10.*f.GetParameter(2));
 566       gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
 567     }
 568     fNRefFigures++;
 569
 570     // tracklet phi resolution
 571     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 572     jgraph++; //skip the frame histo
 573     gm = (TGraphErrors*)fContainer->At(jgraph++);
 574     gs = (TGraphErrors*)fContainer->At(jgraph++);
 575     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 576       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 577       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 578       h = h2->ProjectionY("py", iphi, iphi);
 579       h->Fit(&f, "QN", "", -.5, .5);
 580       Int_t jphi = iphi -1;
 581       gm->SetPoint(jphi, phi, f.GetParameter(1));
 582       gm->SetPointError(jphi, 0., f.GetParError(1));
 583       gs->SetPoint(jphi, phi, f.GetParameter(2));
 584       gs->SetPointError(jphi, 0., f.GetParError(2));
 585     }
 586     fNRefFigures++;
 587   }
 588
 589   return kTRUE;
 590 }
 591
 592
 593 //________________________________________________________
 594 void AliTRDtrackingResolution::Terminate(Option_t *)
 595 {
 596   if(fDebugStream){
 597     delete fDebugStream;
 598     fDebugStream = 0x0;
 599     fDebugLevel = 0;
 600   }
 601   if(HasPostProcess()) PostProcess();
 602 }
 603
 604 //________________________________________________________
 605 void AliTRDtrackingResolution::Fit(TH1 *h, TF1 *f)
 606 {
 607 // Helper function to avoid duplication of code
 608 // Make first guesses on the fit parameters
 609
 610   // find the intial parameters of the fit !! (thanks George)
 611   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 612   Double_t sum = 0.;
 613   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 614   f->SetParLimits(0, 0., 3.*sum);
 615   f->SetParameter(0, .9*sum);
 616
 617   f->SetParLimits(1, -.1, .1);
 618   f->SetParameter(1, 0.);
 619
 620   f->SetParLimits(2, 0., 1.e-1);
 621   f->SetParameter(2, 2.e-2);
 622
 623   f->SetParLimits(3, 0., sum);
 624   f->SetParameter(3, .1*sum);
 625
 626   f->SetParLimits(4, -.3, .3);
 627   f->SetParameter(4, 0.);
 628
 629   f->SetParLimits(5, 0., 1.e2);
 630   f->SetParameter(5, 2.e-1);
 631
 632   h->Fit(f, "QN", "", -0.5, 0.5);
 633 }
 634
 635 //________________________________________________________
 636 TObjArray* AliTRDtrackingResolution::Histos()
 637 {
 638   if(!fContainer) fContainer  = new TObjArray(25);
 639   return fContainer;
 640 }
 641
 642
 643 //________________________________________________________
 644 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 645 {
 646
 647   fReconstructor->SetRecoParam(r);
 648 }