TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TAxis.h>
   2 #include <TCanvas.h>
   3 #include <TFile.h>
   4 #include <TH1F.h>
   5 #include <TGaxis.h>
   6 #include <TGraph.h>
   7 #include <TMath.h>
   8 #include <TMap.h>
   9 #include <TObjArray.h>
  10 #include <TObject.h>
  11 #include <TObjString.h>
  12 #include <TPad.h>
  13 #include <TProfile.h>
  14 #include <TProfile2D.h>
  15 #include <TROOT.h>
  16
  17 #include "AliLog.h"
  18 #include "AliTRDcluster.h"
  19 #include "AliESDHeader.h"
  20 #include "AliESDRun.h"
  21 #include "AliESDtrack.h"
  22 #include "AliTRDgeometry.h"
  23 #include "AliTRDpadPlane.h"
  24 #include "AliTRDSimParam.h"
  25 #include "AliTRDseedV1.h"
  26 #include "AliTRDtrackV1.h"
  27 #include "AliTRDtrackerV1.h"
  28 #include "AliTRDReconstructor.h"
  29 #include "AliTrackReference.h"
  30 #include "AliTrackPointArray.h"
  31 #include "AliTracker.h"
  32 #include "TTreeStream.h"
  33
  34 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  35 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  36 #include "AliTRDcheckDetector.h"
  37
  38 #include <cstdio>
  39 #include <iostream>
  40
  41 ////////////////////////////////////////////////////////////////////////////
  42 //                                                                        //
  43 //  Reconstruction QA                                                     //
  44 //                                                                        //
  45 //  Task doing basic checks for tracking and detector performance         //
  46 //                                                                        //
  47 //  Authors:                                                              //
  48 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  49 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  50 //                                                                        //
  51 ////////////////////////////////////////////////////////////////////////////
  52
  53 //_______________________________________________________
  54 AliTRDcheckDetector::AliTRDcheckDetector():
  55   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  56   ,fEventInfo(0x0)
  57   ,fTriggerNames(0x0)
  58   ,fReconstructor(0x0)
  59   ,fGeo(0x0)
  60 {
  61   //
  62   // Default constructor
  63   //
  64   DefineInput(1,AliTRDeventInfo::Class());
  65   fReconstructor = new AliTRDReconstructor;
  66   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  67   fGeo = new AliTRDgeometry;
  68   InitFunctorList();
  69 }
  70
  71 //_______________________________________________________
  72 AliTRDcheckDetector::~AliTRDcheckDetector(){
  73   //
  74   // Destructor
  75   //
  76   if(fTriggerNames) delete fTriggerNames;
  77   delete fReconstructor;
  78   delete fGeo;
  79 }
  80
  81 //_______________________________________________________
  82 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
  83   //
  84   // Connect the Input data with the task
  85   //
  86   AliTRDrecoTask::ConnectInputData(opt);
  87   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
  88 }
  89
  90 //_______________________________________________________
  91 void AliTRDcheckDetector::CreateOutputObjects(){
  92   //
  93   // Create Output Objects
  94   //
  95   OpenFile(0,"RECREATE");
  96   fContainer = Histos();
  97   if(!fTriggerNames) fTriggerNames = new TMap();
  98 }
  99
 100 //_______________________________________________________
 101 void AliTRDcheckDetector::Exec(Option_t *opt){
 102   //
 103   // Execution function
 104   // Filling TRD quality histos
 105   //
 106   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 107   AliTRDrecoTask::Exec(opt);
 108   Int_t nTracks = 0;            // Count the number of tracks per event
 109   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 110   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 111   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
 112   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
 113   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 114     if(!fTracks->UncheckedAt(iti)) continue;
 115     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 116     if(!fTrackInfo->GetTrack()) continue;
 117     nTracks++;
 118   }
 119   if(nTracks){
 120     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
 121     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
 122   }
 123   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 124     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 125     // also set the label for both histograms
 126     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 127     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 128     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 129     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 130   }
 131   PostData(0, fContainer);
 132 }
 133
 134 //_______________________________________________________
 135 void AliTRDcheckDetector::Terminate(Option_t *){
 136   //
 137   // Terminate function
 138   //
 139 }
 140
 141 //_______________________________________________________
 142 Bool_t AliTRDcheckDetector::PostProcess(){
 143   //
 144   // Do Postprocessing (for the moment set the number of Reference histograms)
 145   //
 146
 147   TH1 * histo = 0x0;
 148   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
 149   histo->GetXaxis()->SetTitle("Number of Tracks");
 150   histo->GetYaxis()->SetTitle("Events");
 151   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
 152   histo->GetXaxis()->SetTitle("Number of Clusters");
 153   histo->GetYaxis()->SetTitle("Entries");
 154   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
 155   histo->GetXaxis()->SetTitle("Number of Tracklets");
 156   histo->GetYaxis()->SetTitle("Entries");
 157   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTrackletsVsFindable));
 158   histo->GetXaxis()->SetTitle("Ratio Found/Findable Tracklets");
 159   histo->GetYaxis()->SetTitle("Entries");
 160   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
 161   histo->GetXaxis()->SetTitle("Number of Clusters");
 162   histo->GetYaxis()->SetTitle("Entries");
 163   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
 164   histo->GetXaxis()->SetTitle("#chi^2");
 165   histo->GetYaxis()->SetTitle("Entries");
 166   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 167   histo->GetXaxis()->SetTitle("Sector");
 168   histo->GetYaxis()->SetTitle("Number of Tracks");
 169   histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
 170   histo->GetXaxis()->SetTitle("Time / 100ns");
 171   histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
 172   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
 173   histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
 174   histo->GetYaxis()->SetTitle("Entries");
 175   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
 176   histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
 177   histo->GetYaxis()->SetTitle("Entries");
 178
 179   // Calculate the purity of the trigger clusters
 180   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 181   TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 182   histoTracks->Divide(histo);
 183   Float_t purities[20], val = 0;
 184   TString triggernames[20];
 185   Int_t nTriggerClasses = 0;
 186   for(Int_t ibin = 1; ibin <= histo->GetNbinsX(); ibin++){
 187     if((val = histoTracks->GetBinContent(ibin))){
 188       purities[nTriggerClasses] = val;
 189       triggernames[nTriggerClasses] = histoTracks->GetXaxis()->GetBinLabel(ibin);
 190       nTriggerClasses++;
 191     }
 192   }
 193   TH1F *hTriggerInf = new TH1F("fTriggerInf", "Trigger Information", TMath::Max(nTriggerClasses, 1), 0, TMath::Max(nTriggerClasses, 1));
 194   for(Int_t ibin = 1; ibin <= nTriggerClasses; ibin++){
 195     hTriggerInf->SetBinContent(ibin, purities[ibin-1]);
 196     hTriggerInf->GetXaxis()->SetBinLabel(ibin, triggernames[ibin-1].Data());
 197   }
 198   hTriggerInf->GetXaxis()->SetTitle("Trigger Cluster");
 199   hTriggerInf->GetYaxis()->SetTitle("Ratio");
 200   hTriggerInf->GetYaxis()->SetRangeUser(0,1);
 201 //      hTriggerInf->SetMarkerColor(kBlue);
 202 //      hTriggerInf->SetMarkerStyle(22);
 203   fContainer->Add(hTriggerInf);
 204   fNRefFigures = 11;
 205   return kTRUE;
 206 }
 207
 208 //_______________________________________________________
 209 Bool_t AliTRDcheckDetector::GetRefFigure(Int_t ifig){
 210   //
 211   // Setting Reference Figures
 212   //
 213   TH1 *h = 0x0, *h1 = 0x0, *h2 = 0x0;
 214   TGaxis *axis = 0x0;
 215   switch(ifig){
 216   case 0:
 217     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 218     break;
 219   case 1:
 220     ((TH1F*)fContainer->At(kNclustersHist))->Draw("pl");
 221     break;
 222   case 2:
 223     h = (TH1F*)fContainer->At(kNtrackletsHist);
 224     if(!h->GetEntries()) break;
 225     h->Scale(100./h->Integral());
 226     h->GetXaxis()->SetRangeUser(.5, 6.5);
 227     h->SetFillColor(kGreen);
 228     h->SetBarOffset(.2);
 229     h->SetBarWidth(.6);
 230     h->Draw("bar1");
 231     break;
 232   case 3:
 233     h = (TH1F*)fContainer->At(kNTrackletsVsFindable);
 234     if(!h->GetEntries()) break;
 235     h->Scale(100./h->Integral());
 236     h->GetXaxis()->SetRangeUser(0.005, 1.005);
 237     h->SetFillColor(kGreen);
 238     h->SetBarOffset(.2);
 239     h->SetBarWidth(.6);
 240     h->Draw("bar1");
 241     break;
 242   case 4:
 243     ((TH1F*)fContainer->At(kNclusterTrackletHist))->Draw("pc");
 244     break;
 245   case 5:
 246     ((TH1F*)fContainer->At(kChi2))->Draw("");
 247     break;
 248   case 6:
 249     h = (TH1F*)fContainer->At(kNTracksSectorHist);
 250     if(!h->GetEntries()) break;
 251     h->Scale(100./h->Integral());
 252     h->SetFillColor(kGreen);
 253     h->SetBarOffset(.2);
 254     h->SetBarWidth(.6);
 255     h->Draw("bar1");
 256     break;
 257   case 7:
 258     h = (TH1F*)fContainer->At(kPulseHeight);
 259     h->SetMarkerStyle(24);
 260     h->SetMarkerColor(kBlack);
 261     h->SetLineColor(kBlack);
 262     h->Draw("e1");
 263     // copy the second histogram in a new one with the same x-dimension as the phs with respect to time
 264     h1 = (TH1F *)fContainer->At(kPulseHeightDistance);
 265     h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
 266     for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
 267       h2->SetBinContent(ibin, h1->GetBinContent(ibin));
 268     h2->SetMarkerStyle(22);
 269     h2->SetMarkerColor(kBlue);
 270     h2->SetLineColor(kBlue);
 271     h2->Draw("e1same");
 272     gPad->Update();
 273     // create axis according to the histogram dimensions of the original second histogram
 274     axis = new TGaxis(gPad->GetUxmin(),
 275                       gPad->GetUymax(),
 276                       gPad->GetUxmax(),
 277                       gPad->GetUymax(),
 278                       -0.08, 4.88, 510,"-L");
 279     axis->SetLineColor(kBlue);
 280     axis->SetLabelColor(kBlue);
 281     axis->SetTextColor(kBlue);
 282     axis->SetTitle("x_{c}-x_{0} / cm");
 283     axis->Draw();
 284     break;
 285   case 8:
 286     ((TH1F*)fContainer->At(kClusterCharge))->Draw("c");
 287     break;
 288   case 9:
 289     ((TH1F*)fContainer->At(kChargeDeposit))->Draw("c");
 290     break;
 291   case 10:
 292     h=(TH1F*)fContainer->At(kPurity);
 293     h->SetBarOffset(.2);
 294     h->SetBarWidth(.6);
 295     h->SetFillColor(kGreen);
 296     h->Draw("bar1");
 297     break;
 298   default:
 299     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 300     break;
 301   }
 302   return kTRUE;
 303 }
 304
 305 //_______________________________________________________
 306 TObjArray *AliTRDcheckDetector::Histos(){
 307   //
 308   // Create QA histograms
 309   //
 310   if(fContainer) return fContainer;
 311
 312   fContainer = new TObjArray(25);
 313   //fContainer->SetOwner(kTRUE);
 314
 315   // Register Histograms
 316   TH1 * histptr = 0x0;
 317   if(!(histptr = (TH1F *)gROOT->FindObject("hNtrks")))
 318     histptr = new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100);
 319   else histptr->Reset();
 320   fContainer->AddAt(histptr, kNTracksEventHist);
 321   if(!(histptr = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 322     histptr = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 323   else histptr->Reset();
 324   fContainer->AddAt(histptr, kNEventsTriggerTracks);
 325   if(!(histptr = (TH1F *)gROOT->FindObject("hNcls")))
 326     histptr = new TH1F("hNcls", "Nr. of clusters per track", 181, -0.5, 180.5);
 327   else histptr->Reset();
 328   fContainer->AddAt(histptr, kNclustersHist);
 329   if(!(histptr = (TH1F *)gROOT->FindObject("hNtls")))
 330     histptr = new TH1F("hNtls", "Nr. tracklets per track", 7, -0.5, 6.5);
 331   else histptr->Reset();
 332   fContainer->AddAt(histptr, kNtrackletsHist);
 333   if(!(histptr = (TH1F *)gROOT->FindObject("hNtlsFindable")))
 334     histptr = new TH1F("hNtlsFindable", "Ratio of found/findable Tracklets" , 101, -0.005, 1.005);
 335   else histptr->Reset();
 336   fContainer->AddAt(histptr, kNTrackletsVsFindable);
 337   if(!(histptr = (TH1F *)gROOT->FindObject("hNclTls")))
 338     histptr = new TH1F("hNclTls","Mean Number of clusters per tracklet", 31, -0.5, 30.5);
 339   else histptr->Reset();
 340   fContainer->AddAt(histptr, kNclusterTrackletHist);
 341   if(!(histptr = (TH1F *)gROOT->FindObject("hChi2")))
 342     histptr = new TH1F("hChi2", "Chi2", 200, 0, 20);
 343   else histptr->Reset();
 344   fContainer->AddAt(histptr, kChi2);
 345   if(!(histptr = (TH1F *)gROOT->FindObject("hChi2n")))
 346     histptr = new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5);
 347   else histptr->Reset();
 348   fContainer->AddAt(histptr, kChi2Normalized);
 349   if(!(histptr = (TH1F *)gROOT->FindObject("hChi2n")))
 350     histptr = new TH1F("hSM", "Track Counts in Supermodule", 18, -0.5, 17.5);
 351   else histptr->Reset();
 352   fContainer->AddAt(histptr, kNTracksSectorHist);
 353   if(!(histptr = (TH1F *)gROOT->FindObject("hPHdetector")))
 354     histptr = new TProfile("hPHdetector", "Average PH", 31, -0.5, 30.5);
 355   else histptr->Reset();
 356   fContainer->AddAt(histptr, kPulseHeight);
 357   if(!(histptr = (TH1F *)gROOT->FindObject("hPHdistance")))
 358     histptr = new TProfile("hPHdistance", "Average PH", 31, -0.08, 4.88);
 359   else histptr->Reset();
 360   fContainer->AddAt(histptr, kPulseHeightDistance);
 361   if(!(histptr = (TH1F *)gROOT->FindObject("hQclDetector")))
 362     histptr = new TH1F("hQclDetector", "Cluster charge", 200, 0, 1200);
 363   else histptr->Reset();
 364   fContainer->AddAt(histptr, kClusterCharge);
 365   if(!(histptr = (TH1F *)gROOT->FindObject("hQTdetector")))
 366     histptr = new TH1F("hQTdetector", "Total Charge Deposit", 6000, 0, 6000);
 367   else histptr->Reset();
 368   fContainer->AddAt(histptr, kChargeDeposit);
 369   if(!(histptr = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 370     histptr = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 371   else histptr->Reset();
 372   fContainer->AddAt(histptr, kNEventsTrigger);
 373
 374   return fContainer;
 375 }
 376
 377 /*
 378 * Plotting Functions
 379 */
 380
 381 //_______________________________________________________
 382 TH1 *AliTRDcheckDetector::PlotMeanNClusters(const AliTRDtrackV1 *track){
 383   //
 384   // Plot the mean number of clusters per tracklet
 385   //
 386   if(track) fTrack = track;
 387   if(!fTrack){
 388     AliWarning("No Track defined.");
 389     return 0x0;
 390   }
 391   TH1 *h = 0x0;
 392   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclusterTrackletHist)))){
 393     AliWarning("No Histogram defined.");
 394     return 0x0;
 395   }
 396   AliTRDseedV1 *tracklet = 0x0;
 397   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 398     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 399     h->Fill(tracklet->GetN());
 400   }
 401   return h;
 402 }
 403
 404 //_______________________________________________________
 405 TH1 *AliTRDcheckDetector::PlotNClusters(const AliTRDtrackV1 *track){
 406   //
 407   // Plot the number of clusters in one track
 408   //
 409   if(track) fTrack = track;
 410   if(!fTrack){
 411     AliWarning("No Track defined.");
 412     return 0x0;
 413   }
 414   TH1 *h = 0x0;
 415   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersHist)))){
 416     AliWarning("No Histogram defined.");
 417     return 0x0;
 418   }
 419
 420   Int_t nclusters = 0;
 421   AliTRDseedV1 *tracklet = 0x0;
 422   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 423     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 424     nclusters += tracklet->GetN();
 425     if(fDebugLevel > 2){
 426       Int_t crossing = tracklet->GetNChange();
 427       AliTRDcluster *c = 0x0;
 428       for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 429         if(!(c = tracklet->GetClusters(itime))) continue;
 430         break;
 431       }
 432       Int_t detector = c->GetDetector();
 433       Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 434       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 435       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 436       Float_t momentum = 0.;
 437       Int_t pdg = 0;
 438       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 439       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 440       if(fMC){
 441         if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 442         pdg = fMC->GetPDG();
 443       }
 444       (*fDebugStream) << "NClusters"
 445         << "Detector="  << detector
 446         << "Sector="    << sector
 447         << "crossing="  << crossing
 448         << "momentum="  << momentum
 449         << "pdg="                               << pdg
 450         << "theta="                     << theta
 451         << "phi="                               << phi
 452         << "kinkIndex=" << kinkIndex
 453         << "TPCncls="           << TPCncls
 454         << "nclusters=" << nclusters
 455         << "\n";
 456     }
 457   }
 458   h->Fill(nclusters);
 459   return h;
 460 }
 461
 462 //_______________________________________________________
 463 TH1 *AliTRDcheckDetector::PlotChi2(const AliTRDtrackV1 *track){
 464   //
 465   // Plot the chi2 of the track
 466   //
 467   if(track) fTrack = track;
 468   if(!fTrack){
 469     AliWarning("No Track defined.");
 470     return 0x0;
 471   }
 472   TH1 *h = 0x0;
 473   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2)))){
 474     AliWarning("No Histogram defined.");
 475     return 0x0;
 476   }
 477   h->Fill(fTrack->GetChi2());
 478   return h;
 479 }
 480
 481 //_______________________________________________________
 482 TH1 *AliTRDcheckDetector::PlotNormalizedChi2(const AliTRDtrackV1 *track){
 483   //
 484   // Plot the chi2 of the track
 485   //
 486   if(track) fTrack = track;
 487   if(!fTrack){
 488     AliWarning("No Track defined.");
 489     return 0x0;
 490   }
 491   TH1 *h = 0x0;
 492   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2Normalized)))){
 493     AliWarning("No Histogram defined.");
 494     return 0x0;
 495   }
 496   Int_t nTracklets = 0;
 497   AliTRDseedV1 *tracklet = 0x0;
 498   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 499     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 500     nTracklets++;
 501   }
 502   h->Fill(fTrack->GetChi2()/nTracklets);
 503   return h;
 504 }
 505
 506
 507 //_______________________________________________________
 508 TH1 *AliTRDcheckDetector::PlotNTracklets(const AliTRDtrackV1 *track){
 509   //
 510   // Plot the number of tracklets
 511   //
 512   if(track) fTrack = track;
 513   if(!fTrack){
 514     AliWarning("No Track defined.");
 515     return 0x0;
 516   }
 517   TH1 *h = 0x0;
 518   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsHist)))){
 519     AliWarning("No Histogram defined.");
 520     return 0x0;
 521   }
 522   Int_t nTracklets = GetNTracklets(fTrack);
 523   h->Fill(nTracklets);
 524   if(fDebugLevel > 3){
 525     if(nTracklets == 1){
 526       // If we have one Tracklet, check in which layer this happens
 527       Int_t layer = -1;
 528       AliTRDseedV1 *tracklet = 0x0;
 529       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 530         if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
 531       }
 532       (*fDebugStream) << "PlotNTracklets"
 533         << "Layer=" << layer
 534         << "\n";
 535     }
 536   }
 537   return h;
 538 }
 539
 540 //_______________________________________________________
 541 TH1 *AliTRDcheckDetector::PlotTrackletsVsFindable(const AliTRDtrackV1 *track){
 542   //
 543   // Plots the ratio of number of tracklets vs.
 544   // number of findable tracklets
 545   //
 546   // Findable tracklets are defined as track prolongation
 547   // to layer i does not hit the dead area +- epsilon
 548   //
 549   // In order to check whether tracklet hist active area in Layer i,
 550   // the track is refitted and the fitted position + an uncertainty
 551   // range is compared to the chamber border (also with a different
 552   // uncertainty)
 553   //
 554   // For the track fit two cases are distinguished:
 555   // If the track is a stand alone track (defined by the status bit
 556   // encoding, then the track is fitted with the tilted Rieman model
 557   // Otherwise the track is fitted with the Kalman fitter in two steps:
 558   // Since the track parameters are give at the outer point, we first
 559   // fit in direction inwards. Afterwards we fit again in direction outwards
 560   // to extrapolate the track to layers which are not reached by the track
 561   // For the Kalman model, the radial track points have to be shifted by
 562   // a distance epsilon in the direction that we want to fit
 563   //
 564   const Float_t epsilon = 0.01;   // dead area tolerance
 565   const Float_t epsilon_R = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
 566   const Float_t delta_y = 0.7;    // Tolerance in the track position in y-direction
 567   const Float_t delta_z = 7.0;    // Tolerance in the track position in z-direction (Padlength)
 568   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 569
 570   if(track) fTrack = track;
 571   if(!fTrack){
 572     AliWarning("No Track defined.");
 573     return 0x0;
 574   }
 575   TH1 *h = 0x0;
 576   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTrackletsVsFindable)))){
 577     AliWarning("No Histogram defined.");
 578     return 0x0;
 579   }
 580   Int_t nFound = 0, nFindable = 0;
 581   Int_t stack = -1;
 582   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 583   Double_t y = 0., z = 0.;
 584   AliTRDseedV1 *tracklet = 0x0;
 585   AliTRDpadPlane *pp;
 586   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 587     if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
 588       tracklet->SetReconstructor(fReconstructor);
 589       nFound++;
 590     }
 591   }
 592   // 2 Different cases:
 593   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 594   // 2nd barrel track: here we propagate the track to the layers
 595   AliTrackPoint points[6];
 596   Float_t xyz[3];
 597   memset(xyz, 0, sizeof(Float_t) * 3);
 598   if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 599     // stand alone track
 600     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 601       xyz[0] = x_anode[il];
 602       points[il].SetXYZ(xyz);
 603     }
 604     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
 605   } else {
 606     // barrel track
 607     //
 608     // 2 Steps:
 609     // -> Kalman inwards
 610     // -> Kalman outwards
 611     AliTRDtrackV1 copy_track(*fTrack);  // Do Kalman on a (non-constant) copy of the track
 612     AliTrackPoint points_inward[6], points_outward[6];
 613     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 614       // In order to avoid complications in the Kalman filter if the track points have the same radial
 615       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
 616       // in the direction we want to go
 617       // The track points have to be in reverse order for the Kalman Filter inwards
 618       xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
 619       points_inward[il].SetXYZ(xyz);
 620       xyz[0] = x_anode[il] + epsilon_R;
 621       points_outward[il].SetXYZ(xyz);
 622     }
 623     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 624       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
 625     // Kalman inwards
 626     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kFALSE, 6, points_inward);
 627     memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
 628     // Kalman outwards
 629     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kTRUE, 6, points_inward);
 630     memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 631   }
 632   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 633     y = points[il].GetY();
 634     z = points[il].GetZ();
 635     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
 636     pp = fGeo->GetPadPlane(il, stack);
 637     ymin = pp->GetCol0() + epsilon;
 638     ymax = pp->GetColEnd() - epsilon;
 639     zmin = pp->GetRowEnd() + epsilon;
 640     zmax = pp->GetRow0() - epsilon;
 641     // ignore y-crossing (material)
 642     if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
 643       if(fDebugLevel > 3){
 644         Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetMeanz() : 0};
 645         Int_t hasTracklet = tracklet ? 1 : 0;
 646         (*fDebugStream)   << "GetFindableTracklets"
 647           << "layer="     << il
 648           << "ytracklet=" << pos_tracklet[0]
 649           << "ytrack="    << y
 650           << "ztracklet=" << pos_tracklet[1]
 651           << "ztrack="    << z
 652           << "tracklet="  << hasTracklet
 653           << "\n";
 654       }
 655   }
 656
 657   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
 658   if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 659   return h;
 660 }
 661
 662 //_______________________________________________________
 663 TH1 *AliTRDcheckDetector::PlotPulseHeight(const AliTRDtrackV1 *track){
 664   //
 665   // Plot the average pulse height
 666   //
 667   if(track) fTrack = track;
 668   if(!fTrack){
 669     AliWarning("No Track defined.");
 670     return 0x0;
 671   }
 672   TProfile *h = 0x0;
 673   if(!(h = dynamic_cast<TProfile *>(fContainer->At(kPulseHeight)))){
 674     AliWarning("No Histogram defined.");
 675     return 0x0;
 676   }
 677   AliTRDseedV1 *tracklet = 0x0;
 678   AliTRDcluster *c = 0x0;
 679   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 680     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 681     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 682       if(!(c = tracklet->GetClusters(itime))) continue;
 683       Int_t localtime        = c->GetLocalTimeBin();
 684       Double_t absolute_charge = TMath::Abs(c->GetQ());
 685       h->Fill(localtime, absolute_charge);
 686       if(fDebugLevel > 3){
 687         Int_t crossing = tracklet->GetNChange();
 688         Int_t detector = c->GetDetector();
 689         Double_t distance[2];
 690         GetDistanceToTracklet(distance, tracklet, c);
 691         Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 692         Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 693         Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 694         Float_t momentum = 0.;
 695         Int_t pdg = 0;
 696         Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 697         UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 698         if(fMC){
 699           if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 700           pdg = fMC->GetPDG();
 701         }
 702         (*fDebugStream) << "PulseHeight"
 703           << "Detector="        << detector
 704           << "Sector="          << sector
 705           << "crossing="        << crossing
 706           << "Timebin="         << localtime
 707           << "Charge="          << absolute_charge
 708           << "momentum="        << momentum
 709           << "pdg="                             << pdg
 710           << "theta="                   << theta
 711           << "phi="                             << phi
 712           << "kinkIndex="       << kinkIndex
 713           << "TPCncls="         << TPCncls
 714           << "dy="        << distance[0]
 715           << "dz="        << distance[1]
 716           << "c.="        << c
 717           << "\n";
 718       }
 719     }
 720   }
 721   return h;
 722 }
 723
 724 //_______________________________________________________
 725 TH1 *AliTRDcheckDetector::PlotPHSdistance(const AliTRDtrackV1 *track){
 726   //
 727   // Plots the average pulse height vs the distance from the anode wire
 728   // (plus const anode wire offset)
 729   //
 730   if(track) fTrack = track;
 731   if(!fTrack){
 732     AliWarning("No Track defined.");
 733     return 0x0;
 734   }
 735   TProfile *h = 0x0;
 736   if(!(h = dynamic_cast<TProfile *>(fContainer->At(kPulseHeightDistance)))){
 737     AliWarning("No Histogram defined.");
 738     return 0x0;
 739   }
 740   Float_t offset = .5*AliTRDgeometry::CamHght();
 741   AliTRDseedV1 *tracklet = 0x0;
 742   AliTRDcluster *c = 0x0;
 743   Double_t distance = 0;
 744   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 745     if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
 746     tracklet->ResetClusterIter();
 747     while((c = tracklet->NextCluster())){
 748       distance = tracklet->GetX0() - c->GetX() + offset;
 749       h->Fill(distance, TMath::Abs(c->GetQ()));
 750     }
 751   }
 752   return h;
 753 }
 754
 755 //_______________________________________________________
 756 TH1 *AliTRDcheckDetector::PlotClusterCharge(const AliTRDtrackV1 *track){
 757   //
 758   // Plot the cluster charge
 759   //
 760   if(track) fTrack = track;
 761   if(!fTrack){
 762     AliWarning("No Track defined.");
 763     return 0x0;
 764   }
 765   TH1 *h = 0x0;
 766   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kClusterCharge)))){
 767     AliWarning("No Histogram defined.");
 768     return 0x0;
 769   }
 770   AliTRDseedV1 *tracklet = 0x0;
 771   AliTRDcluster *c = 0x0;
 772   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 773     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 774     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 775       if(!(c = tracklet->GetClusters(itime))) continue;
 776       h->Fill(c->GetQ());
 777     }
 778   }
 779   return h;
 780 }
 781
 782 //_______________________________________________________
 783 TH1 *AliTRDcheckDetector::PlotChargeDeposit(const AliTRDtrackV1 *track){
 784   //
 785   // Plot the charge deposit per chamber
 786   //
 787   if(track) fTrack = track;
 788   if(!fTrack){
 789     AliWarning("No Track defined.");
 790     return 0x0;
 791   }
 792   TH1 *h = 0x0;
 793   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeDeposit)))){
 794     AliWarning("No Histogram defined.");
 795     return 0x0;
 796   }
 797   AliTRDseedV1 *tracklet = 0x0;
 798   AliTRDcluster *c = 0x0, *c1 = 0x0;    // c1 for the Debug Stream
 799   Double_t Qtot = 0;
 800   Int_t nTracklets = 0;
 801   if(fDebugLevel > 3)
 802     nTracklets = GetNTracklets(fTrack); // fill NTracklet to the Debug Stream
 803   for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
 804     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 805     Qtot = 0;
 806     c1 = 0x0;
 807     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 808       if(!(c = tracklet->GetClusters(itime))) continue;
 809       if(!c1) c1 = c;
 810       Qtot += TMath::Abs(c->GetQ());
 811     }
 812     h->Fill(Qtot);
 813     if(fDebugLevel > 3){
 814       Int_t crossing = tracklet->GetNChange();
 815       Int_t detector = c1->GetDetector();
 816       Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 817       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 818       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 819       Float_t momentum = 0.;
 820       Int_t pdg = 0;
 821       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 822       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 823       if(fMC){
 824               if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 825         pdg = fMC->GetPDG();
 826       }
 827       (*fDebugStream) << "ChargeDeposit"
 828         << "Detector="  << detector
 829         << "Sector="    << sector
 830         << "crossing="  << crossing
 831         << "momentum="  << momentum
 832         << "nTracklets="<< nTracklets
 833         << "pdg="                               << pdg
 834         << "theta="                     << theta
 835         << "phi="                               << phi
 836         << "kinkIndex=" << kinkIndex
 837         << "TPCncls="           << TPCncls
 838         << "QT="        << Qtot
 839         << "\n";
 840     }
 841   }
 842   return h;
 843 }
 844
 845 //_______________________________________________________
 846 TH1 *AliTRDcheckDetector::PlotTracksSector(const AliTRDtrackV1 *track){
 847   //
 848   // Plot the number of tracks per Sector
 849   //
 850   if(track) fTrack = track;
 851   if(!fTrack){
 852     AliWarning("No Track defined.");
 853     return 0x0;
 854   }
 855   TH1 *h = 0x0;
 856   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTracksSectorHist)))){
 857     AliWarning("No Histogram defined.");
 858     return 0x0;
 859   }
 860   AliTRDseedV1 *tracklet = 0x0;
 861   AliTRDcluster *c = 0x0;
 862   Int_t sector = -1;
 863   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 864     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 865     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 866       if(!(c = tracklet->GetClusters(itime))) continue;
 867       sector = static_cast<Int_t>(c->GetDetector()/AliTRDgeometry::kNdets);
 868     }
 869     break;
 870   }
 871   h->Fill(sector);
 872   return h;
 873 }
 874
 875 //_______________________________________________________
 876 Int_t AliTRDcheckDetector::GetNTracklets(const AliTRDtrackV1 *track){
 877   //
 878   // Count the number of tracklets per track
 879   //
 880   if(!track){
 881     AliError("No track");
 882     return 0;
 883   }
 884   Int_t nTracklets = 0;
 885   AliTRDseedV1 *tracklet = 0x0;
 886   for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 887     if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()) nTracklets++;
 888   }
 889   return nTracklets;
 890 }
 891
 892 //________________________________________________________
 893 void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
 894 {
 895
 896   fReconstructor->SetRecoParam(r);
 897 }
 898
 899 //________________________________________________________
 900 void AliTRDcheckDetector::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c){
 901   dist[0] = c->GetY() - (tracklet->GetYfit(0) + tracklet->GetYfit(1)*(c->GetX() - tracklet->GetX0()));
 902   dist[1] = c->GetZ() - tracklet->GetMeanz();
 903 }