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   InitFunctorList();
  98 }
  99
 100 //________________________________________________________
 101 AliTRDtrackingResolution::~AliTRDtrackingResolution()
 102 {
 103   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
 104   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
 105   delete fGeo;
 106   delete fReconstructor;
 107   if(gGeoManager) delete gGeoManager;
 108 }
 109
 110
 111 //________________________________________________________
 112 void AliTRDtrackingResolution::CreateOutputObjects()
 113 {
 114   // spatial resolution
 115   OpenFile(0, "RECREATE");
 116
 117   fContainer = Histos();
 118 }
 119
 120 // //________________________________________________________
 121 // void AliTRDtrackingResolution::Exec(Option_t *)
 122 // {
 123 //   // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
 124 //   // angular Resolution: res = Tracklet angle - TrackRef Angle
 125 //
 126 //   Int_t nTrackInfos = fTracks->GetEntriesFast();
 127 //   if(fDebugLevel>=2 && nTrackInfos){
 128 //     printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
 129 //   }
 130 //   const Int_t kNLayers = AliTRDgeometry::kNlayer;
 131 //   Int_t pdg, nly, ncrs;
 132 //   Double_t p, dy, theta/*, dphi, dymc, dzmc, dphimc*/;
 133 //   Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
 134 //   Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
 135 //
 136 //   AliTRDpadPlane *pp = 0x0;
 137 //   AliTrackPoint tr[kNLayers], tk[kNLayers];
 138 //   AliExternalTrackParam *fOp = 0x0;
 139 //   AliTRDtrackV1 *fTrack = 0x0;
 140 //   AliTRDtrackInfo *fInfo = 0x0;
 141 //   for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
 142 //     // check if ESD and MC-Information are available
 143 //     if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
 144 //     if(!(fTrack = fInfo->GetTrack())) continue;
 145 //     if(!(fOp = fInfo->GetOuterParam())) continue;
 146 //     pdg = fInfo->GetPDG();
 147 //     nly = 0; ncrs = 0; theta = 0.;
 148 //
 149 //     if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
 150 //
 151 //     p = fOp->P();
 152 //     Int_t npts = 0;
 153 //     memset(fP, 0, kNLayers*sizeof(Float_t));
 154 //     memset(fX, 0, kNLayers*sizeof(Float_t));
 155 //     memset(fY, 0, kNLayers*sizeof(Float_t));
 156 //     memset(fZ, 0, kNLayers*sizeof(Float_t));
 157 //     memset(fPhi, 0, kNLayers*sizeof(Float_t));
 158 //     memset(fTheta, 0, kNLayers*sizeof(Float_t));
 159 //     memset(fLayerMap, 0, kNLayers*sizeof(Bool_t));
 160 //     memset(fMCMap, 0, kNLayers*sizeof(Bool_t));
 161 //     AliTRDseedV1 *fTracklet = 0x0;
 162 //     for(Int_t iplane = 0; iplane < kNLayers; iplane++){
 163 //       if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
 164 //       if(!fTracklet->IsOK()) continue;
 165 //
 166 //       // Book point arrays
 167 //       fLayerMap[iplane] = kTRUE;
 168 //       tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
 169 //       tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
 170 //       npts++;
 171 //
 172 //       if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
 173 //
 174 //       // define reference values
 175 //       fP[iplane]   = p;
 176 //       fX[iplane]   = fTracklet->GetX0();
 177 //       fY[iplane]   = fTracklet->GetYref(0);
 178 //       fZ[iplane]   = fTracklet->GetZref(0);
 179 //       fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1));
 180 //       fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1));
 181 //
 182 //
 183 //       // RESOLUTION (compare to MC)
 184 //       if(HasMCdata()){
 185 //         if(fInfo->GetNTrackRefs() >= 2){
 186 //           Double_t pmc, ymc, zmc, phiMC, thetaMC;
 187 //           if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){
 188 //             fMCMap[iplane] = kTRUE;
 189 //             fP[iplane]     = pmc;
 190 //             fY[iplane]     = ymc;
 191 //             fZ[iplane]     = zmc;
 192 //             fPhi[iplane]   = phiMC;
 193 //             fTheta[iplane] = thetaMC;
 194 //           }
 195 //         }
 196 //       }
 197 //       Float_t phi   = fPhi[iplane]*TMath::RadToDeg();
 198 //       theta += fTheta[iplane]; nly++;
 199 //       if(fTracklet->GetNChange()) ncrs++;
 200 //
 201 //       // Do clusters residuals
 202 //       if(!fTracklet->Fit(kFALSE)) continue;
 203 //       AliTRDcluster *c = 0x0;
 204 //       for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 205 //         if(!(c = fTracklet->GetClusters(ic))) continue;
 206 //
 207 //         dy = fTracklet->GetYat(c->GetX()) - c->GetY();
 208 //         ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
 209 //
 210 //         if(fDebugLevel>=2){
 211 //           Float_t q = c->GetQ();
 212 //           // Get z-position with respect to anode wire
 213 //           AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 214 //           Int_t det = c->GetDetector();
 215 //           Float_t x = c->GetX();
 216 //           Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
 217 //           Int_t stack = fGeo->GetStack(det);
 218 //           pp = fGeo->GetPadPlane(iplane, stack);
 219 //           Float_t row0 = pp->GetRow0();
 220 //           Float_t d  =  row0 - z + simParam->GetAnodeWireOffset();
 221 //           d -= ((Int_t)(2 * d)) / 2.0;
 222 //           if (d > 0.25) d  = 0.5 - d;
 223 //
 224 //           (*fDebugStream) << "ResidualClusters"
 225 //             << "ly="   << iplane
 226 //             << "stk="  << stack
 227 //             << "pdg="  << pdg
 228 //             << "phi="  << fPhi[iplane]
 229 //             << "tht="  << fTheta[iplane]
 230 //             << "q="    << q
 231 //             << "x="    << x
 232 //             << "z="    << z
 233 //             << "d="    << d
 234 //             << "dy="   << dy
 235 //             << "\n";
 236 //         }
 237 //       }
 238 //       pp = 0x0;
 239 //     }
 240 //     if(nly) theta /= nly;
 241 //     if(fDebugLevel>=1){
 242 //       (*fDebugStream) << "TrackStatistics"
 243 //         << "nly="   << nly
 244 //         << "ncrs="  << ncrs
 245 //         << "tht="   << theta
 246 //         << "\n";
 247 //     }
 248 //
 249 //
 250 // //     // this protection we might drop TODO
 251 // //     if(fTrack->GetNumberOfTracklets() < 6) continue;
 252 // //
 253 // //     AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
 254 // //     Int_t iref = 0;
 255 // //     for(Int_t ip=0; ip<kNLayers; ip++){
 256 // //       if(!fLayerMap[ip]) continue;
 257 // //       fTracklet = fTrack->GetTracklet(ip);
 258 // //       // recalculate fit based on the new tilt correction
 259 // //       fTracklet->Fit();
 260 // //
 261 // //       dy = fTracklet->GetYfit(0) - tr[iref].GetY();
 262 // //       ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
 263 // //
 264 // //       dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
 265 // //       ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
 266 // //
 267 // //       if(HasMCdata()){
 268 // //         dymc = fY[ip] - tr[iref].GetY();
 269 // //         ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
 270 // //
 271 // //         dzmc = fZ[ip] - tr[iref].GetZ();
 272 // //         ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
 273 // //
 274 // //         dphimc = fPhi[ip] - fTracklet->GetYfit(1);
 275 // //         ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
 276 // //       }
 277 // //
 278 // //       iref++;
 279 // //
 280 // //       if(fDebugLevel>=1){
 281 // //         (*fDebugStream) << "RiemannTrack"
 282 // //           << "ly="    << ip
 283 // //           << "mc="    << fMCMap[ip]
 284 // //           << "p="     << fP[ip]
 285 // //           << "phi="   << fPhi[ip]
 286 // //           << "tht="   << fTheta[ip]
 287 // //           << "dy="    << dy
 288 // //           << "dphi="  << dphi
 289 // //           << "dymc="  << dymc
 290 // //           << "dzmc="  << dzmc
 291 // //           << "dphimc="<< dphimc
 292 // //           << "\n";
 293 // //       }
 294 // //     }
 295 //
 296 // //  if(!gGeoManager) TGeoManager::Import("geometry.root");
 297 // //     AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
 298 // //     for(Int_t ip=0; ip<nc; ip++){
 299 // //       dy = cl[ip].GetY() - tr[ip].GetY();
 300 // //      ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
 301 // //       dz = cl[ip].GetZ() - tr[ip].GetZ();
 302 // //       if(fDebugLevel>=1){
 303 // //         (*fDebugStream) << "KalmanTrack"
 304 // //           << "dy="                  << dy
 305 // //           << "dz="                  << dz
 306 // // /*          << "phi="                     << phi
 307 // //           << "theta="             << theta
 308 // //           << "dphi="              << dphi*/
 309 // //           << "\n";
 310 // //       }
 311 // //     }
 312 //
 313 //
 314 //   }
 315 //   PostData(0, fContainer);
 316 // }
 317
 318 //________________________________________________________
 319 TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
 320 {
 321   if(track) fTrack = track;
 322   if(!fTrack){
 323     AliWarning("No Track defined.");
 324     return 0x0;
 325   }
 326   TH1 *h = 0x0;
 327   if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
 328     AliWarning("No output histogram defined.");
 329     return 0x0;
 330   }
 331
 332   Int_t pdg = (HasMCdata() && fMC) ? fMC->GetPDG() : 0;
 333   UChar_t s;
 334   Float_t x0, y0, z0, dy, dydx, dzdx;
 335   AliTRDseedV1 *fTracklet = 0x0;
 336   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 337     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 338     if(!fTracklet->IsOK()) continue;
 339     x0 = fTracklet->GetX0();
 340
 341     // retrive the track angle with the chamber
 342     if(HasMCdata() && fMC){
 343       if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 344     }else{
 345       y0   = fTracklet->GetYref(0);
 346       z0   = fTracklet->GetZref(0);
 347       dydx = fTracklet->GetYref(1);
 348       dzdx = fTracklet->GetZref(1);
 349     }
 350
 351     AliTRDseedV1 trklt(*fTracklet);
 352     if(!trklt.Fit(kFALSE)) continue;
 353
 354     AliTRDcluster *c = 0x0;
 355     fTracklet->ResetClusterIter(kFALSE);
 356     while((c = fTracklet->PrevCluster())){
 357       dy = trklt.GetYat(c->GetX()) - c->GetY();
 358       h->Fill(dydx, dy);
 359
 360       if(fDebugLevel>=1){
 361         Float_t q = c->GetQ();
 362         // Get z-position with respect to anode wire
 363         AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
 364         Int_t det = c->GetDetector();
 365         Float_t x = c->GetX();
 366         Float_t z = z0-(x0-x)*dzdx;
 367         Int_t istk = fGeo->GetStack(det);
 368         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 369         Float_t row0 = pp->GetRow0();
 370         Float_t d  =  row0 - z + simParam->GetAnodeWireOffset();
 371         d -= ((Int_t)(2 * d)) / 2.0;
 372         if (d > 0.25) d  = 0.5 - d;
 373
 374         (*fDebugStream) << "ClusterResidual"
 375           << "ly="   << ily
 376           << "stk="  << istk
 377           << "pdg="  << pdg
 378           << "dydx=" << dydx
 379           << "dzdx=" << dzdx
 380           << "q="    << q
 381           << "x="    << x
 382           << "z="    << z
 383           << "d="    << d
 384           << "dy="   << dy
 385           << "\n";
 386       }
 387     }
 388   }
 389   return h;
 390 }
 391
 392 //________________________________________________________
 393 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
 394 {
 395   if(!HasMCdata()){
 396     AliWarning("No MC defined. Results will not be available.");
 397     return 0x0;
 398   }
 399   if(track) fTrack = track;
 400   if(!fTrack){
 401     AliWarning("No Track defined.");
 402     return 0x0;
 403   }
 404   TH1 *h = 0x0;
 405   if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
 406     AliWarning("No output histogram defined.");
 407     return 0x0;
 408   }
 409   if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
 410     AliWarning("No output histogram defined.");
 411     return 0x0;
 412   }
 413   if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
 414     AliWarning("No output histogram defined.");
 415     return 0x0;
 416   }
 417   if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
 418     AliWarning("No output histogram defined.");
 419     return 0x0;
 420   }
 421   //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
 422   UChar_t s;
 423   Int_t pdg = fMC->GetPDG(), det=-1;
 424   Float_t x0, y0, z0, dy, dydx, dzdx;
 425   AliTRDseedV1 *fTracklet = 0x0;
 426   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
 427     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
 428     if(!fTracklet->IsOK()) continue;
 429     //printf("process layer[%d]\n", ily);
 430     // retrive the track position and direction within the chamber
 431     det = fTracklet->GetDetector();
 432     x0  = fTracklet->GetX0();
 433     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue;
 434
 435     // recalculate tracklet based on the MC info
 436     AliTRDseedV1 tt(*fTracklet);
 437     tt.SetZref(0, z0);
 438     tt.SetZref(1, dzdx);
 439     if(!tt.Fit(kFALSE)) continue;
 440     //tt.Update();
 441     dy = tt.GetYfit(0) - y0;
 442     Float_t dphi   = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
 443     Float_t dz = 100.;
 444     Bool_t  cross = fTracklet->GetNChange();
 445     if(cross){
 446       Double_t *xyz = tt.GetCrossXYZ();
 447       dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
 448       ((TH2I*)fContainer->At(kTrackletZResolution))->Fill(dzdx, dz);
 449     }
 450
 451     // Fill Histograms
 452     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
 453     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
 454
 455     // Fill Debug stream
 456     if(fDebugLevel>=1){
 457       Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
 458       (*fDebugStream) << "TrkltResolution"
 459         << "det="                 << det
 460         << "pdg="     << pdg
 461         << "p="       << p
 462         << "ymc="     << y0
 463         << "zmc="     << z0
 464         << "dydx="    << dydx
 465         << "dzdx="    << dzdx
 466         << "cross="   << cross
 467         << "dy="                  << dy
 468         << "dz="                  << dz
 469         << "dphi="              << dphi
 470         << "\n";
 471     }
 472
 473     Int_t istk = AliTRDgeometry::GetStack(det);
 474     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
 475     Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
 476     Float_t tilt = fTracklet->GetTilt();
 477
 478     AliTRDcluster *c = 0x0;
 479     fTracklet->ResetClusterIter(kFALSE);
 480     while((c = fTracklet->PrevCluster())){
 481       Float_t  q = TMath::Abs(c->GetQ());
 482       Float_t xc = c->GetX();
 483       Float_t yc = c->GetY();
 484       Float_t zc = c->GetZ();
 485       Float_t dx = x0 - xc;
 486       Float_t yt = y0 - dx*dydx;
 487       Float_t zt = z0 - dx*dzdx;
 488       dy = yt - (yc - tilt*(zc-zt));
 489
 490       // Fill Histograms
 491       if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
 492
 493       // Fill Debug Tree
 494       if(fDebugLevel>=1){
 495         Float_t d = zr0 - zt;
 496         d -= ((Int_t)(2 * d)) / 2.0;
 497         if (d > 0.25) d  = 0.5 - d;
 498         Int_t label = fMC->GetLabel();
 499         (*fDebugStream) << "ClusterResolution"
 500           << "det="  << det
 501           << "pdg="  << pdg
 502           << "dydx=" << dydx
 503           << "dzdx=" << dzdx
 504           << "q="    << q
 505           << "d="    << d
 506           << "dy="   << dy
 507           << "xc="   << xc
 508           << "yc="   << yc
 509           << "zc="   << zc
 510           << "yt="   << yt
 511           << "zt="   << zt
 512           << "lbl="   << label
 513           << "\n";
 514       }
 515     }
 516   }
 517   return h;
 518 }
 519
 520
 521 //________________________________________________________
 522 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
 523 {
 524   TAxis *ax = 0x0;
 525   TGraphErrors *g = 0x0;
 526   switch(ifig){
 527   case kClusterResidual:
 528     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 529     g->Draw("apl");
 530     ax = g->GetHistogram()->GetYaxis();
 531     ax->SetRangeUser(-.5, 1.);
 532     ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
 533     ax = g->GetHistogram()->GetXaxis();
 534     ax->SetTitle("tg(#phi)");
 535     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 536     g->Draw("pl");
 537     return;
 538   case kClusterResolution:
 539     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 540     ax = g->GetHistogram()->GetYaxis();
 541     ax->SetRangeUser(-.5, 1.);
 542     ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
 543     ax = g->GetHistogram()->GetXaxis();
 544     ax->SetTitle("tg(#phi)");
 545     g->Draw("apl");
 546     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 547     g->Draw("pl");
 548     return;
 549   case kTrackletYResolution:
 550     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 551     ax = g->GetHistogram()->GetYaxis();
 552     ax->SetRangeUser(-.5, 1.);
 553     ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
 554     ax = g->GetHistogram()->GetXaxis();
 555     ax->SetTitle("tg(#phi)");
 556     g->Draw("apl");
 557     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 558     g->Draw("pl");
 559     return;
 560   case kTrackletZResolution:
 561     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 562     ax = g->GetHistogram()->GetYaxis();
 563     ax->SetRangeUser(-.5, 1.);
 564     ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
 565     ax = g->GetHistogram()->GetXaxis();
 566     ax->SetTitle("tg(#theta)");
 567     g->Draw("apl");
 568     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 569     g->Draw("pl");
 570     return;
 571   case kTrackletAngleResolution:
 572     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
 573     ax = g->GetHistogram()->GetYaxis();
 574     ax->SetRangeUser(-.05, .2);
 575     ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
 576     ax = g->GetHistogram()->GetXaxis();
 577     ax->SetTitle("tg(#phi)");
 578     g->Draw("apl");
 579     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
 580     g->Draw("pl");
 581     return;
 582   default:
 583     AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
 584     return;
 585   }
 586   AliInfo(Form("Reference plot [%d] missing result", ifig));
 587 }
 588
 589
 590 //________________________________________________________
 591 Bool_t AliTRDtrackingResolution::PostProcess()
 592 {
 593   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 594   if (!fContainer) {
 595     Printf("ERROR: list not available");
 596     return kFALSE;
 597   }
 598   fNRefFigures = fContainer->GetEntriesFast();
 599   if(!fGraphS){
 600     fGraphS = new TObjArray(fNRefFigures);
 601     fGraphS->SetOwner();
 602   }
 603   if(!fGraphM){
 604     fGraphM = new TObjArray(fNRefFigures);
 605     fGraphM->SetOwner();
 606   }
 607
 608   TH2I *h2 = 0x0;
 609   TH1D *h = 0x0;
 610   TGraphErrors *gm = 0x0, *gs = 0x0;
 611
 612   // define models
 613   TF1 f("f1", "gaus", -.5, .5);
 614
 615   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
 616
 617   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
 618
 619   TCanvas *c = 0x0;
 620   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
 621   char opt[5];
 622   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
 623
 624
 625   //PROCESS RESIDUAL DISTRIBUTIONS
 626
 627   // Clusters residuals
 628   h2 = (TH2I *)(fContainer->At(kClusterResidual));
 629   gm = new TGraphErrors(h2->GetNbinsX());
 630   gm->SetLineColor(kBlue);
 631   gm->SetMarkerStyle(7);
 632   gm->SetMarkerColor(kBlue);
 633   gm->SetNameTitle("clm", "");
 634   fGraphM->AddAt(gm, kClusterResidual);
 635   gs = new TGraphErrors(h2->GetNbinsX());
 636   gs->SetLineColor(kRed);
 637   gs->SetMarkerStyle(23);
 638   gs->SetMarkerColor(kRed);
 639   gs->SetNameTitle("cls", "");
 640   fGraphS->AddAt(gs, kClusterResidual);
 641   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 642     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 643     h = h2->ProjectionY("py", ibin, ibin);
 644     AdjustF1(h, &fc);
 645
 646     if(IsVisual()){c->cd(); c->SetLogy();}
 647     h->Fit(&fc, opt, "", -0.5, 0.5);
 648     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 649
 650     gm->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(1));
 651     gm->SetPointError(ibin - 1, 0., 10.*fc.GetParError(1));
 652     gs->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(2));
 653     gs->SetPointError(ibin - 1, 0., 10.*fc.GetParError(2));
 654   }
 655
 656
 657   //PROCESS RESOLUTION DISTRIBUTIONS
 658
 659   if(HasMCdata()){
 660     // cluster y resolution
 661     h2 = (TH2I*)fContainer->At(kClusterResolution);
 662     gm = new TGraphErrors(h2->GetNbinsX());
 663     gm->SetLineColor(kBlue);
 664     gm->SetMarkerStyle(7);
 665     gm->SetMarkerColor(kBlue);
 666     gm->SetNameTitle("clym", "");
 667     fGraphM->AddAt(gm, kClusterResolution);
 668     gs = new TGraphErrors(h2->GetNbinsX());
 669     gs->SetLineColor(kRed);
 670     gs->SetMarkerStyle(23);
 671     gs->SetMarkerColor(kRed);
 672     gs->SetNameTitle("clys", "");
 673     fGraphS->AddAt(gs, kClusterResolution);
 674     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 675       h = h2->ProjectionY("py", iphi, iphi);
 676       if(h->GetEntries()<100.) continue;
 677       AdjustF1(h, &f);
 678
 679       if(IsVisual()){c->cd(); c->SetLogy();}
 680       h->Fit(&f, opt, "", -0.5, 0.5);
 681       if(IsVerbose()){
 682         printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
 683       }
 684       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 685
 686       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 687       Int_t jphi = iphi -1;
 688       gm->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*f.GetParameter(1));
 689       gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
 690       gs->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*f.GetParameter(2));
 691       gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
 692     }
 693
 694     // tracklet y resolution
 695     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
 696     gm = new TGraphErrors(h2->GetNbinsX());
 697     gm->SetLineColor(kBlue);
 698     gm->SetMarkerStyle(7);
 699     gm->SetMarkerColor(kBlue);
 700     gm->SetNameTitle("trkltym", "");
 701     fGraphM->AddAt(gm, kTrackletYResolution);
 702     gs = new TGraphErrors(h2->GetNbinsX());
 703     gs->SetLineColor(kRed);
 704     gs->SetMarkerStyle(23);
 705     gs->SetMarkerColor(kRed);
 706     gs->SetNameTitle("trkltys", "");
 707     fGraphS->AddAt(gs, kTrackletYResolution);
 708     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 709       h = h2->ProjectionY("py", iphi, iphi);
 710       AdjustF1(h, &fb);
 711
 712       if(IsVisual()){c->cd(); c->SetLogy();}
 713       h->Fit(&fb, opt, "", -0.5, 0.5);
 714       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 715
 716       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 717       Int_t jphi = iphi -1;
 718       gm->SetPoint(jphi, phi, 10.*fb.GetParameter(1));
 719       gm->SetPointError(jphi, 0., 10.*fb.GetParError(1));
 720       gs->SetPoint(jphi, phi, 10.*fb.GetParameter(2));
 721       gs->SetPointError(jphi, 0., 10.*fb.GetParError(2));
 722     }
 723
 724     // tracklet z resolution
 725     h2 = (TH2I*)fContainer->At(kTrackletZResolution);
 726     gm = new TGraphErrors(h2->GetNbinsX());
 727     gm->SetLineColor(kBlue);
 728     gm->SetMarkerStyle(7);
 729     gm->SetMarkerColor(kBlue);
 730     gm->SetNameTitle("trkltzm", "");
 731     fGraphM->AddAt(gm, kTrackletZResolution);
 732     gs = new TGraphErrors(h2->GetNbinsX());
 733     gs->SetLineColor(kRed);
 734     gs->SetMarkerStyle(23);
 735     gs->SetMarkerColor(kRed);
 736     gs->SetNameTitle("trkltzs", "");
 737     fGraphS->AddAt(gs, kTrackletZResolution);
 738     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 739       h = h2->ProjectionY("py", iphi, iphi);
 740       AdjustF1(h, &fb);
 741
 742       if(IsVisual()){c->cd(); c->SetLogy();}
 743       h->Fit(&fb, opt, "", -0.5, 0.5);
 744       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 745
 746       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 747       Int_t jphi = iphi -1;
 748       gm->SetPoint(jphi, phi, 10.*fb.GetParameter(1));
 749       gm->SetPointError(jphi, 0., 10.*fb.GetParError(1));
 750       gs->SetPoint(jphi, phi, 10.*fb.GetParameter(2));
 751       gs->SetPointError(jphi, 0., 10.*fb.GetParError(2));
 752     }
 753
 754     // tracklet phi resolution
 755     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
 756     gm = new TGraphErrors(h2->GetNbinsX());
 757     gm->SetLineColor(kBlue);
 758     gm->SetMarkerStyle(7);
 759     gm->SetMarkerColor(kBlue);
 760     gm->SetNameTitle("trkltym", "");
 761     fGraphM->AddAt(gm, kTrackletAngleResolution);
 762     gs = new TGraphErrors(h2->GetNbinsX());
 763     gs->SetLineColor(kRed);
 764     gs->SetMarkerStyle(23);
 765     gs->SetMarkerColor(kRed);
 766     gs->SetNameTitle("trkltys", "");
 767     fGraphS->AddAt(gs, kTrackletAngleResolution);
 768     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 769       h = h2->ProjectionY("py", iphi, iphi);
 770
 771       if(IsVisual()){c->cd(); c->SetLogy();}
 772       h->Fit(&f, opt, "", -0.5, 0.5);
 773       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
 774
 775       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 776       Int_t jphi = iphi -1;
 777       gm->SetPoint(jphi, phi, f.GetParameter(1));
 778       gm->SetPointError(jphi, 0., f.GetParError(1));
 779       gs->SetPoint(jphi, phi, f.GetParameter(2));
 780       gs->SetPointError(jphi, 0., f.GetParError(2));
 781     }
 782   }
 783   if(c) delete c;
 784
 785   return kTRUE;
 786 }
 787
 788
 789 //________________________________________________________
 790 void AliTRDtrackingResolution::Terminate(Option_t *)
 791 {
 792   if(fDebugStream){
 793     delete fDebugStream;
 794     fDebugStream = 0x0;
 795     fDebugLevel = 0;
 796   }
 797   if(HasPostProcess()) PostProcess();
 798 }
 799
 800 //________________________________________________________
 801 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
 802 {
 803 // Helper function to avoid duplication of code
 804 // Make first guesses on the fit parameters
 805
 806   // find the intial parameters of the fit !! (thanks George)
 807   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
 808   Double_t sum = 0.;
 809   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
 810   f->SetParLimits(0, 0., 3.*sum);
 811   f->SetParameter(0, .9*sum);
 812
 813   f->SetParLimits(1, -.2, .2);
 814   f->SetParameter(1, -0.1);
 815
 816   f->SetParLimits(2, 0., 4.e-1);
 817   f->SetParameter(2, 2.e-2);
 818   if(f->GetNpar() <= 4) return;
 819
 820   f->SetParLimits(3, 0., sum);
 821   f->SetParameter(3, .1*sum);
 822
 823   f->SetParLimits(4, -.3, .3);
 824   f->SetParameter(4, 0.);
 825
 826   f->SetParLimits(5, 0., 1.e2);
 827   f->SetParameter(5, 2.e-1);
 828 }
 829
 830 //________________________________________________________
 831 TObjArray* AliTRDtrackingResolution::Histos()
 832 {
 833   if(fContainer) return fContainer;
 834
 835   fContainer  = new TObjArray(5);
 836
 837   TH1 *h = 0x0;
 838   // cluster to tracklet residuals [2]
 839   fContainer->AddAt(h = new TH2I("fYCl", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5), kClusterResidual);
 840   h->GetXaxis()->SetTitle("tg(#phi)");
 841   h->GetYaxis()->SetTitle("#Delta y [cm]");
 842   h->GetZaxis()->SetTitle("entries");
 843 //   // tracklet to Riemann fit residuals [2]
 844 //   fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
 845 //   fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
 846 //   fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
 847 //   fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
 848
 849   // Resolution histos
 850   if(HasMCdata()){
 851     // cluster y resolution [0]
 852     fContainer->AddAt(h = new TH2I("fYClMC", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5), kClusterResolution);
 853     h->GetXaxis()->SetTitle("tg(#phi)");
 854     h->GetYaxis()->SetTitle("#Delta y [cm]");
 855     h->GetZaxis()->SetTitle("entries");
 856     // tracklet y resolution [0]
 857     fContainer->AddAt(h = new TH2I("fYTrkltMC", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5), kTrackletYResolution);
 858     h->GetXaxis()->SetTitle("tg(#phi)");
 859     h->GetYaxis()->SetTitle("#Delta y [cm]");
 860     h->GetZaxis()->SetTitle("entries");
 861     // tracklet y resolution [0]
 862     fContainer->AddAt(h = new TH2I("fZTrkltMC", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5), kTrackletZResolution);
 863     h->GetXaxis()->SetTitle("tg(#theta)");
 864     h->GetYaxis()->SetTitle("#Delta z [cm]");
 865     h->GetZaxis()->SetTitle("entries");
 866     // tracklet angular resolution [1]
 867     fContainer->AddAt(h = new TH2I("fPhiTrkltMC", "Tracklet Resolution (Angular)", 31, -.48, .48, 100, -10., 10.), kTrackletAngleResolution);
 868     h->GetXaxis()->SetTitle("tg(#phi)");
 869     h->GetYaxis()->SetTitle("#Delta #phi [deg]");
 870     h->GetZaxis()->SetTitle("entries");
 871
 872 //     // Riemann track resolution [y, z, angular]
 873 //     fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
 874 //     fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
 875 //     fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
 876 //
 877 //     Kalman track resolution [y, z, angular]
 878 //     fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
 879 //     fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
 880 //     fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
 881   }
 882   return fContainer;
 883 }
 884
 885
 886 //________________________________________________________
 887 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 888 {
 889
 890   fReconstructor->SetRecoParam(r);
 891 }