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