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 <TList.h>
  32 #include <TH2.h>
  33 #include <TH1.h>
  34 #include <TF1.h>
  35 #include <TProfile.h>
  36 #include <TGraphErrors.h>
  37 #include <TMath.h>
  38 #include "TTreeStream.h"
  39 #include "TGeoManager.h"
  40
  41 #include "AliAnalysisManager.h"
  42 #include "AliTrackReference.h"
  43 #include "AliTrackPointArray.h"
  44 #include "AliCDBManager.h"
  45
  46 #include "AliTRDcluster.h"
  47 #include "AliTRDseedV1.h"
  48 #include "AliTRDtrackV1.h"
  49 #include "AliTRDtrackerV1.h"
  50 #include "AliTRDReconstructor.h"
  51 #include "AliTRDrecoParam.h"
  52
  53 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  54 #include "AliTRDtrackingResolution.h"
  55
  56 ClassImp(AliTRDtrackingResolution)
  57
  58 //________________________________________________________
  59 AliTRDtrackingResolution::AliTRDtrackingResolution(const char * name):
  60   AliAnalysisTask(name, "")
  61   ,fTracks(0x0)
  62   ,fHistos(0x0)
  63   ,fReconstructor(0x0)
  64   ,fHasMCdata(kTRUE)
  65   ,fDebugLevel(0)
  66   ,fDebugStream(0x0)
  67 {
  68   TGeoManager::Import("geometry.root");
  69   fReconstructor = new AliTRDReconstructor();
  70   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  71
  72   DefineInput(0, TObjArray::Class());
  73   DefineOutput(0, TList::Class());
  74 }
  75
  76 //________________________________________________________
  77 AliTRDtrackingResolution::~AliTRDtrackingResolution()
  78 {
  79   fHistos->Delete(); delete fHistos;
  80   delete fReconstructor;
  81 }
  82
  83 //________________________________________________________
  84 void AliTRDtrackingResolution::ConnectInputData(Option_t *){
  85   fTracks = dynamic_cast<TObjArray *>(GetInputData(0));
  86 }
  87
  88 //________________________________________________________
  89 void AliTRDtrackingResolution::CreateOutputObjects()
  90 {
  91   // spatial resolution
  92   OpenFile(0, "RECREATE");
  93
  94   fHistos = new TList();
  95
  96   // Resolution histos
  97   if(HasMCdata()){
  98     // tracklet resolution [0]
  99     fHistos->AddAt(new TH2I("fY", "", 21, -21., 21., 100, -.5, .5), 0);
 100     // tracklet angular resolution [1]
 101     fHistos->AddAt(new TH2I("fPhi", "", 21, -21., 21., 100, -10., 10.), 1);
 102   }
 103   // Residual histos
 104   // cluster to tracklet residuals [2]
 105   Int_t position = HasMCdata() ? 2 : 0;
 106   fHistos->AddAt(new TH2I("fYClRes", "", 21, -21., 21., 100, -.5, .5), position);
 107 }
 108
 109 //________________________________________________________
 110 void AliTRDtrackingResolution::Exec(Option_t *)
 111 {
 112   // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
 113   // angular Resolution: res = Tracklet angle - TrackRef Angle
 114
 115   Int_t nTrackInfos = fTracks->GetEntriesFast();
 116   if(fDebugLevel>=2) printf("Number of Histograms: %d\n", fHistos->GetEntries());
 117
 118   Double_t dy, dz;
 119   AliTrackPoint tr[kNLayers], tk[kNLayers];
 120   AliTRDtrackV1 *fTrack = 0x0;
 121   AliTRDtrackInfo *fInfo = 0x0;
 122   if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
 123   for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
 124     // check if ESD and MC-Information are available
 125     if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
 126     if(!(fTrack = fInfo->GetTRDtrack())) continue;
 127     if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
 128
 129
 130     Int_t npts = 0;
 131     AliTRDseedV1 *fTracklet = 0x0;
 132     for(Int_t iplane = 0; iplane < kNLayers; iplane++){
 133       if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
 134       if(!fTracklet->IsOK()) continue;
 135
 136       // Book point arrays
 137       tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
 138       tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
 139       npts++;
 140
 141       if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
 142
 143       Float_t phi = TMath::ATan(fTracklet->GetYref(1));
 144
 145       // RESOLUTION (compare to MC)
 146       if(HasMCdata()){
 147         if(fInfo->GetNTrackRefs() >= 2){
 148         Double_t phiMC;
 149         if(Resolution(fTracklet, fInfo, phiMC)) phi = phiMC;
 150         }
 151       }
 152
 153       // Do clusters residuals
 154       if(!fTracklet->Fit(kFALSE)) continue;
 155       Int_t histpos = HasMCdata() ? 2 : 0;
 156       AliTRDcluster *c = 0x0;
 157       for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 158         if(!(c = fTracklet->GetClusters(ic))) continue;
 159
 160         dy = fTracklet->GetYat(c->GetX()) - c->GetY();
 161         ((TH2I*)fHistos->At(histpos))->Fill(phi*TMath::RadToDeg(), dy);
 162         if(fDebugLevel>=1){
 163           Float_t q = c->GetQ();
 164           (*fDebugStream) << "ClsTrkltResidual"
 165             << "plane="         << iplane
 166             << "phi="     << phi
 167             << "q="       << q
 168             << "dy="      << dy
 169             << "\n";
 170         }
 171       }
 172     }
 173
 174
 175     // this protection we might drop TODO
 176     if(fTrack->GetNumberOfTracklets() < 6) continue;
 177
 178     AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
 179     for(Int_t ip=0; ip<npts; ip++){
 180       dy = tk[ip].GetY() - tr[ip].GetY();
 181       dz = tk[ip].GetZ() - tr[ip].GetZ();
 182       if(fDebugLevel>=1){
 183         (*fDebugStream) << "ResidualsRT"
 184           << "dy="                << dy
 185           << "dz="                << dz
 186 /*          << "phi="                   << phi
 187           << "theta="           << theta
 188           << "dphi="            << dphi*/
 189           << "\n";
 190       }
 191     }
 192
 193 //     AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
 194 //     for(Int_t ip=0; ip<nc; ip++){
 195 //       dy = cl[ip].GetY() - tr[ip].GetY();
 196 //       dz = cl[ip].GetZ() - tr[ip].GetZ();
 197 //       if(fDebugLevel>=1){
 198 //         (*fDebugStream) << "ResidualsKF"
 199 //           << "dy="             << dy
 200 //           << "dz="             << dz
 201 // /*          << "phi="                        << phi
 202 //           << "theta="                << theta
 203 //           << "dphi="         << dphi*/
 204 //           << "\n";
 205 //       }
 206 //     }
 207
 208
 209   }
 210   PostData(0, fHistos);
 211 }
 212
 213
 214 //________________________________________________________
 215 Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &phi)
 216 {
 217
 218   AliTrackReference *fTrackRefs[2] = {0x0, 0x0},   *tempTrackRef = 0x0;
 219
 220   // check for 2 track ref where the radial position has a distance less than 3.7mm
 221   Int_t nFound = 0;
 222   for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
 223     if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
 224     if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
 225     if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
 226     fTrackRefs[nFound++] = tempTrackRef;
 227     if(nFound == 2) break;
 228   }
 229   if(nFound < 2){
 230     if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
 231     return kFALSE;
 232   }
 233   // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
 234
 235
 236   // RESOLUTION
 237   Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
 238   if(dx <= 0.){
 239     if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
 240     return kFALSE;
 241   }
 242   Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
 243   Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
 244   Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
 245   Double_t ymc =  fTrackRefs[1]->LocalY() - dydx*dx0;
 246   Double_t zmc =  fTrackRefs[1]->Z() - dzdx*dx0;
 247
 248   // recalculate tracklet based on the MC info
 249   tracklet->SetZref(0, zmc);
 250   tracklet->SetZref(1, dzdx);
 251   tracklet->Fit();
 252   Double_t dy = tracklet->GetYfit(0) - ymc;
 253   Double_t dz = tracklet->GetZfit(0) - zmc;
 254
 255   //res_y *= 100; // in mm
 256   Double_t momentum = fTrackRefs[0]->P();
 257
 258   phi   = TMath::ATan(dydx);
 259   Double_t theta = TMath::ATan(dzdx);
 260   Double_t dphi   = TMath::ATan(tracklet->GetYfit(1)) - phi;
 261   if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
 262
 263   // Fill Histograms
 264   if(TMath::Abs(dx-3.7)<1.E-3){
 265     ((TH2I*)fHistos->At(0))->Fill(phi*TMath::RadToDeg(), dy);
 266     ((TH2I*)fHistos->At(1))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
 267   }
 268   // Fill Debug Tree
 269   if(fDebugLevel>=1){
 270     Int_t iplane = tracklet->GetPlane();
 271     (*fDebugStream) << "TrkltResolution"
 272       << "plane="               << iplane
 273       << "p="       << momentum
 274       << "dx="      << dx
 275       << "dy="            << dy
 276       << "dz="            << dz
 277       << "phi="                 << phi
 278       << "theta="               << theta
 279       << "dphi="                << dphi
 280       << "ymc="     << ymc
 281       << "zmc="     << zmc
 282       << "tracklet.="<<tracklet
 283       << "\n";
 284   }
 285
 286   return kTRUE;
 287 }
 288
 289
 290 //________________________________________________________
 291 void AliTRDtrackingResolution::Terminate(Option_t *)
 292 {
 293   if(fDebugStream) delete fDebugStream;
 294
 295   TH2I *h2 = 0x0;
 296   TH1D *h = 0x0;
 297   TF1 f("f1", "gaus", -.5, .5);
 298   if(HasMCdata()){
 299     //process distributions
 300     fHistos = dynamic_cast<TList*>(GetOutputData(0));
 301     if (!fHistos) {
 302       Printf("ERROR: list not available");
 303       return;
 304     }
 305
 306   // y resolution
 307     h2 = (TH2I*)fHistos->At(0);
 308     TGraphErrors *gm = new TGraphErrors(h2->GetNbinsX());
 309     gm->SetNameTitle("meany", "Mean dy");
 310     TGraphErrors *gs = new TGraphErrors(h2->GetNbinsX());
 311     gs->SetNameTitle("sigmy", "Sigma y");
 312     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 313       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 314       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 315       h = h2->ProjectionY("py", iphi, iphi);
 316       h->Fit(&f, "QN", "", -.5, .5);
 317       Int_t jphi = iphi -1;
 318       gm->SetPoint(jphi, phi, f.GetParameter(1));
 319       gm->SetPointError(jphi, 0., f.GetParError(1));
 320       gs->SetPoint(jphi, phi, f.GetParameter(2));
 321       gs->SetPointError(jphi, 0., f.GetParError(2));
 322     }
 323     fHistos->Add(gm);
 324     fHistos->Add(gs);
 325
 326     // phi resolution
 327     h2 = (TH2I*)fHistos->At(1);
 328     gm = new TGraphErrors(h2->GetNbinsX());
 329     gm->SetNameTitle("meanphi", "Mean Phi");
 330     gs = new TGraphErrors(h2->GetNbinsX());
 331     gs->SetNameTitle("sigmphi", "Sigma Phi");
 332     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
 333       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
 334       f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
 335       h = h2->ProjectionY("py", iphi, iphi);
 336       h->Fit(&f, "QN", "", -.5, .5);
 337       Int_t jphi = iphi -1;
 338       gm->SetPoint(jphi, phi, f.GetParameter(1));
 339       gm->SetPointError(jphi, 0., f.GetParError(1));
 340       gs->SetPoint(jphi, phi, f.GetParameter(2));
 341       gs->SetPointError(jphi, 0., f.GetParError(2));
 342     }
 343     fHistos->Add(gm);
 344     fHistos->Add(gs);
 345   }
 346
 347   // Fit clusters residuals
 348   Int_t position_residuals = fHasMCdata ? 2 : 0;
 349   h2 = (TH2I *)(fHistos->At(position_residuals));
 350   TGraphErrors *residuals_mean = new TGraphErrors(h2->GetNbinsX());
 351   residuals_mean->SetLineColor(kGreen);
 352   residuals_mean->SetMarkerStyle(22);
 353   residuals_mean->SetMarkerColor(kGreen);
 354   TGraphErrors *residuals_sigma = new TGraphErrors(h2->GetNbinsX());
 355   residuals_mean->SetNameTitle("residuals_mean", "Residuals Mean Phi");
 356   residuals_sigma->SetNameTitle("residuals_sigma", "Residuals Sigma Phi");
 357   residuals_sigma->SetLineColor(kRed);
 358   residuals_sigma->SetMarkerStyle(23);
 359   residuals_sigma->SetMarkerColor(kRed);
 360   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
 361     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
 362     Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
 363     h = h2->ProjectionY("py", ibin, ibin);
 364     h->Fit(&f, "QN", "", -0.5, 0.5);
 365     residuals_mean->SetPoint(ibin - 1, phi, f.GetParameter(1));
 366     residuals_mean->SetPointError(ibin - 1, dphi, f.GetParError(1));
 367     residuals_sigma->SetPoint(ibin - 1, phi, f.GetParameter(2));
 368     residuals_sigma->SetPointError(ibin - 1, dphi, f.GetParError(2));
 369   }
 370   fHistos->Add(residuals_mean);
 371   fHistos->Add(residuals_sigma);
 372 }
 373
 374 //________________________________________________________
 375 void AliTRDtrackingResolution::SetDebugLevel(Int_t level){
 376   fDebugLevel = level;
 377   if(!fDebugLevel) return;
 378   if(fDebugStream) return;
 379   fDebugStream = new TTreeSRedirector("TRD.Resolution.root");
 380 }
 381
 382 //________________________________________________________
 383 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
 384 {
 385   fReconstructor->SetRecoParam(r);
 386 }