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) printf("Number of Histograms: %d\n", Histos()->GetEntries());
 194
 195   Int_t pdg;
 196   Double_t p, dy, dphi, dymc, dzmc, dphimc;
 197   Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
 198   Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
 199
 200   AliTRDpadPlane *pp = 0x0;
 201   AliTrackPoint tr[kNLayers], tk[kNLayers];
 202   AliExternalTrackParam *fOp = 0x0;
 203   AliTRDtrackV1 *fTrack = 0x0;
 204   AliTRDtrackInfo *fInfo = 0x0;
 205   if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
 206   for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
 207     // check if ESD and MC-Information are available
 208     if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
 209     if(!(fTrack = fInfo->GetTRDtrack())) continue;
 210     if(!(fOp = fInfo->GetOuterParam())) continue;
 211     pdg = fInfo->GetPDG();
 212
 213     if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
 214
 215     p = fOp->P();
 216     Int_t npts = 0;
 217     memset(fP, 0, kNLayers*sizeof(Float_t));
 218     memset(fX, 0, kNLayers*sizeof(Float_t));
 219     memset(fY, 0, kNLayers*sizeof(Float_t));
 220     memset(fZ, 0, kNLayers*sizeof(Float_t));
 221     memset(fPhi, 0, kNLayers*sizeof(Float_t));
 222     memset(fTheta, 0, kNLayers*sizeof(Float_t));
 223     memset(fLayerMap, 0, kNLayers*sizeof(Bool_t));
 224     memset(fMCMap, 0, kNLayers*sizeof(Bool_t));
 225     AliTRDseedV1 *fTracklet = 0x0;
 226     for(Int_t iplane = 0; iplane < kNLayers; iplane++){
 227       if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
 228       if(!fTracklet->IsOK()) continue;
 229
 230       // Book point arrays
 231       fLayerMap[iplane] = kTRUE;
 232       tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
 233       tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
 234       npts++;
 235
 236       if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
 237
 238       // define reference values
 239       fP[iplane]   = p;
 240       fX[iplane]   = fTracklet->GetX0();
 241       fY[iplane]   = fTracklet->GetYref(0);
 242       fZ[iplane]   = fTracklet->GetZref(0);
 243       fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1));
 244       fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1));
 245
 246
 247       // RESOLUTION (compare to MC)
 248       if(HasMCdata()){
 249         if(fInfo->GetNTrackRefs() >= 2){
 250           Double_t pmc, ymc, zmc, phiMC, thetaMC;
 251           if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){
 252             fMCMap[iplane] = kTRUE;
 253             fP[iplane]     = pmc;
 254             fY[iplane]     = ymc;
 255             fZ[iplane]     = zmc;
 256             fPhi[iplane]   = phiMC;
 257             fTheta[iplane] = thetaMC;
 258           }
 259         }
 260       }
 261       Float_t phi   = fPhi[iplane]*TMath::RadToDeg();
 262       //Float_t theta = fTheta[iplane]*TMath::RadToDeg();
 263
 264       // Do clusters residuals
 265       if(!fTracklet->Fit(kFALSE)) continue;
 266       AliTRDcluster *c = 0x0;
 267       for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 268         if(!(c = fTracklet->GetClusters(ic))) continue;
 269
 270         dy = fTracklet->GetYat(c->GetX()) - c->GetY();
 271         ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
 272
 273         if(fDebugLevel>=2){
 274           Float_t q = c->GetQ();
 275           // Get z-position with respect to anode wire
 276           AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 277           Int_t det = c->GetDetector();
 278           Float_t x = c->GetX();
 279           Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
 280           Int_t stack = fGeo->GetStack(det);
 281           pp = fGeo->GetPadPlane(iplane, stack);
 282           Float_t row0 = pp->GetRow0();
 283           Float_t d  =  row0 - z + simParam->GetAnodeWireOffset();
 284           d -= ((Int_t)(2 * d)) / 2.0;
 285           //if (d > 0.25) d  = 0.5 - d;
 286
 287           (*fDebugStream) << "ResidualClusters"
 288             << "ly="   << iplane
 289             << "stk="  << stack
 290             << "pdg="  << pdg
 291             << "phi="  << fPhi[iplane]
 292             << "tht="  << fTheta[iplane]
 293             << "q="    << q
 294             << "x="    << x
 295             << "z="    << z
 296             << "d="    << d
 297             << "dy="   << dy
 298             << "\n";
 299         }
 300       }
 301       pp = 0x0;
 302     }
 303
 304
 305     // this protection we might drop TODO
 306     if(fTrack->GetNumberOfTracklets() < 6) continue;
 307
 308     AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
 309     Int_t iref = 0;
 310     for(Int_t ip=0; ip<kNLayers; ip++){
 311       if(!fLayerMap[ip]) continue;
 312       fTracklet = fTrack->GetTracklet(ip);
 313       // recalculate fit based on the new tilt correction
 314       fTracklet->Fit();
 315
 316       dy = fTracklet->GetYfit(0) - tr[iref].GetY();
 317       ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
 318
 319       dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
 320       ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
 321
 322       if(HasMCdata()){
 323         dymc = fY[ip] - tr[iref].GetY();
 324         ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
 325
 326         dzmc = fZ[ip] - tr[iref].GetZ();
 327         ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
 328
 329         dphimc = fPhi[ip] - fTracklet->GetYfit(1);
 330         ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
 331       }
 332
 333       iref++;
 334
 335       if(fDebugLevel>=2){
 336         (*fDebugStream) << "RiemannTrack"
 337           << "ly="    << ip
 338           << "mc="    << fMCMap[ip]
 339           << "p="     << fP[ip]
 340           << "phi="   << fPhi[ip]
 341           << "tht="   << fTheta[ip]
 342           << "dy="    << dy
 343           << "dphi="  << dphi
 344           << "dymc="  << dymc
 345           << "dzmc="  << dzmc
 346           << "dphimc="<< dphimc
 347           << "\n";
 348       }
 349     }
 350
 351 //  if(!gGeoManager) TGeoManager::Import("geometry.root");
 352 //     AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
 353 //     for(Int_t ip=0; ip<nc; ip++){
 354 //       dy = cl[ip].GetY() - tr[ip].GetY();
 355 //      ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
 356 //       dz = cl[ip].GetZ() - tr[ip].GetZ();
 357 //       if(fDebugLevel>=1){
 358 //         (*fDebugStream) << "KalmanTrack"
 359 //           << "dy="             << dy
 360 //           << "dz="             << dz
 361 // /*          << "phi="                        << phi
 362 //           << "theta="                << theta
 363 //           << "dphi="         << dphi*/
 364 //           << "\n";
 365 //       }
 366 //     }
 367
 368
 369   }
 370   PostData(0, fContainer);
 371 }
 372
 373 //________________________________________________________
 374 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t */*opt*/)
 375 {
 376   switch(ifig){
 377   case 0:
 378     first = (Int_t)kGraphStart; last = first+3;
 379     break;
 380   case 1:
 381     first = (Int_t)kGraphStart+3; last = first+3;
 382     break;
 383   case 2:
 384     first = (Int_t)kGraphStart+6; last = first+3;
 385     break;
 386   default:
 387     first = (Int_t)kGraphStart; last = first;
 388     break;
 389   }
 390 }
 391
 392
 393 //________________________________________________________
 394 Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &p, Double_t &ymc, Double_t &zmc, Double_t &phi, Double_t &theta)
 395 {
 396
 397   AliTrackReference *fTrackRefs[2] = {0x0, 0x0},   *tempTrackRef = 0x0;
 398
 399   // check for 2 track ref where the radial position has a distance less than 3.7mm
 400   Int_t nFound = 0;
 401   for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
 402     if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
 403     if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
 404     if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
 405     fTrackRefs[nFound++] = tempTrackRef;
 406     if(nFound == 2) break;
 407   }
 408   if(nFound < 2){
 409     if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
 410     return kFALSE;
 411   }
 412   // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
 413
 414
 415   // RESOLUTION
 416   Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
 417   if(dx <= 0.){
 418     if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
 419     return kFALSE;
 420   }
 421   Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
 422   Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
 423   Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
 424   ymc =  fTrackRefs[1]->LocalY() - dydx*dx0;
 425   zmc =  fTrackRefs[1]->Z() - dzdx*dx0;
 426
 427   // recalculate tracklet based on the MC info
 428   tracklet->SetZref(0, zmc);
 429   tracklet->SetZref(1, -dzdx); // TODO
 430   tracklet->Fit();
 431   Double_t dy = tracklet->GetYfit(0) - ymc;
 432   Double_t dz = tracklet->GetZfit(0) - zmc;
 433
 434   //res_y *= 100; // in mm
 435   p = fTrackRefs[0]->P();
 436
 437   phi   = TMath::ATan(dydx);
 438   theta = TMath::ATan(dzdx);
 439   Double_t dphi   = TMath::ATan(tracklet->GetYfit(1)) - phi;
 440   if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
 441
 442   // Fill Histograms
 443   if(TMath::Abs(dx-3.7)<1.E-3){
 444     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(phi*TMath::RadToDeg(), dy);
 445     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
 446   }
 447   // Fill Debug Tree
 448   if(fDebugLevel>=2){
 449     Int_t iplane = tracklet->GetPlane();
 450     Int_t pdg = fInfo->GetPDG();
 451     (*fDebugStream) << "ResolutionTrklt"
 452       << "ly="            << iplane
 453       << "pdg="     << pdg
 454       << "p="       << p
 455       << "phi="                 << phi
 456       << "tht="           << theta
 457       << "ymc="     << ymc
 458       << "zmc="     << zmc
 459       << "dx="      << dx
 460       << "dy="            << dy
 461       << "dz="            << dz
 462       << "dphi="                << dphi
 463       << "\n";
 464   }
 465
 466   return kTRUE;
 467 }
 468
 469 //________________________________________________________
 470 Bool_t AliTRDtrackingResolution::PostProcess()
 471 {
 472   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 473   fNRefFigures = 0;
 474   if (!fContainer) {
 475     Printf("ERROR: list not available");
 476     return kFALSE;
 477   }
 478
 479   TH2I *h2 = 0x0;
 480   TH1D *h = 0x0;
 481   TGraphErrors *gm = 0x0, *gs = 0x0;
 482   TF1 f("f1", "gaus", -.5, .5);
 483   // define iterator over graphs
 484   Int_t jgraph = (Int_t)kGraphStart;
 485
 486   //PROCESS RESIDUAL DISTRIBUTIONS
 487
 488   // Clusters residuals
 489   // define model
 490   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 491   h2 = (TH2I *)(fContainer->At(kClusterYResidual));
 492   jgraph++; //skip the frame histo
 493   gm = (TGraphErrors*)fContainer->At(jgraph++);
 494   gs = (TGraphErrors*)fContainer->At(jgraph++);
 495   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 496     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 497     Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
 498     h = h2->ProjectionY("py", ibin, ibin);
 499     Fit(h, &fc);
 500     gm->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(1));
 501     gm->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(1));
 502     gs->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(2));
 503     gs->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(2));
 504   }
 505   fNRefFigures++;
 506
 507
 508   //PROCESS RESOLUTION DISTRIBUTIONS
 509   if(HasMCdata()){
 510     // tracklet y resolution
 511     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 512     jgraph++; //skip the frame histo
 513     gm = (TGraphErrors*)fContainer->At(jgraph++);
 514     gs = (TGraphErrors*)fContainer->At(jgraph++);
 515     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 516       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 517       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 518       h = h2->ProjectionY("py", iphi, iphi);
 519       Fit(h, &fc);
 520       Int_t jphi = iphi -1;
 521       gm->SetPoint(jphi, phi, 10.*f.GetParameter(1));
 522       gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
 523       gs->SetPoint(jphi, phi, 10.*f.GetParameter(2));
 524       gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
 525     }
 526     fNRefFigures++;
 527
 528     // tracklet phi resolution
 529     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 530     jgraph++; //skip the frame histo
 531     gm = (TGraphErrors*)fContainer->At(jgraph++);
 532     gs = (TGraphErrors*)fContainer->At(jgraph++);
 533     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 534       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 535       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 536       h = h2->ProjectionY("py", iphi, iphi);
 537       h->Fit(&f, "QN", "", -.5, .5);
 538       Int_t jphi = iphi -1;
 539       gm->SetPoint(jphi, phi, f.GetParameter(1));
 540       gm->SetPointError(jphi, 0., f.GetParError(1));
 541       gs->SetPoint(jphi, phi, f.GetParameter(2));
 542       gs->SetPointError(jphi, 0., f.GetParError(2));
 543     }
 544     fNRefFigures++;
 545   }
 546
 547   return kTRUE;
 548 }
 549
 550
 551 //________________________________________________________
 552 void AliTRDtrackingResolution::Terminate(Option_t *)
 553 {
 554   if(fDebugStream){
 555     delete fDebugStream;
 556     fDebugStream = 0x0;
 557     fDebugLevel = 0;
 558   }
 559   if(HasPostProcess()) PostProcess();
 560 }
 561
 562 //________________________________________________________
 563 void AliTRDtrackingResolution::Fit(TH1 *h, TF1 *f)
 564 {
 565 // Helper function to avoid duplication of code
 566 // Make first guesses on the fit parameters
 567
 568   // find the intial parameters of the fit !! (thanks George)
 569   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 570   Double_t sum = 0.;
 571   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 572   f->SetParLimits(0, 0., 3.*sum);
 573   f->SetParameter(0, .9*sum);
 574
 575   f->SetParLimits(1, -.1, .1);
 576   f->SetParameter(1, 0.);
 577
 578   f->SetParLimits(2, 0., 1.e-1);
 579   f->SetParameter(2, 2.e-2);
 580
 581   f->SetParLimits(3, 0., sum);
 582   f->SetParameter(3, .1*sum);
 583
 584   f->SetParLimits(4, -.3, .3);
 585   f->SetParameter(4, 0.);
 586
 587   f->SetParLimits(5, 0., 1.e2);
 588   f->SetParameter(5, 2.e-1);
 589
 590   h->Fit(f, "QN", "", -0.5, 0.5);
 591 }
 592
 593 //________________________________________________________
 594 TObjArray* AliTRDtrackingResolution::Histos()
 595 {
 596   if(!fContainer) fContainer  = new TObjArray(25);
 597   return fContainer;
 598 }
 599
 600
 601 //________________________________________________________
 602 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 603 {
 604
 605   fReconstructor->SetRecoParam(r);
 606 }