]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/qaRec/AliTRDcheckDetector.cxx
Fixes to the Kalman filter (Markus)
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckDetector.cxx
index edeb60a8986395f704723379c49932a2f57ced77..9abe08a863a77d11bdb7553b9b4ad58f9cdee9d4 100644 (file)
@@ -2,19 +2,36 @@
 #include <TH1F.h>
 #include <TGraph.h>
 #include <TMath.h>
+#include <TMap.h>
 #include <TObjArray.h>
+#include <TObject.h>
+#include <TObjString.h>
 #include <TProfile.h>
+#include <TProfile2D.h>
 #include <TROOT.h>
 
+#include "AliLog.h"
 #include "AliTRDcluster.h"
+#include "AliESDHeader.h"
+#include "AliESDRun.h"
+#include "AliESDtrack.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
 #include "AliTRDseedV1.h"
 #include "AliTRDtrackV1.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDReconstructor.h"
+#include "AliTrackReference.h"
+#include "AliTrackPointArray.h"
+#include "AliTracker.h"
 #include "TTreeStream.h"
 
 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
+#include "AliTRDtrackInfo/AliTRDeventInfo.h"
 #include "AliTRDcheckDetector.h"
 
 #include <cstdio>
+#include <iostream>
 
 ////////////////////////////////////////////////////////////////////////////
 //                                                                        //
 
 //_______________________________________________________
 AliTRDcheckDetector::AliTRDcheckDetector():
-  AliTRDrecoTask("detChecker", "Basic Detector Checker")
-  ,fPHSdetector(0x0)
-  ,fPHSsector(0x0)
-  ,fQCLdetector(0x0)
-  ,fQCLsector(0x0)
-  ,fQTdetector(0x0)
-  ,fQTsector(0x0)
+  AliTRDrecoTask("DetChecker", "Basic Detector Checker")
+  ,fEventInfo(0x0)
+  ,fTriggerNames(0x0)
+  ,fReconstructor(0x0)
 {
   //
   // Default constructor
   //
+  DefineInput(1,AliTRDeventInfo::Class());
+  fReconstructor = new AliTRDReconstructor;
+  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
+  InitFunctorList();
 }
 
 //_______________________________________________________
@@ -48,12 +66,18 @@ AliTRDcheckDetector::~AliTRDcheckDetector(){
   //
   // Destructor
   // 
-  if(fPHSdetector) delete fPHSdetector;
-  if(fPHSsector) delete fPHSsector;
-  if(fQCLdetector) delete fQCLdetector;
-  if(fQCLsector) delete fQCLsector;
-  if(fQTdetector) delete fQTdetector;
-  if(fQTsector) delete fQTsector;
+  if(fEventInfo) delete fEventInfo;
+  if(fTriggerNames) delete fTriggerNames;
+  delete fReconstructor;
+}
+
+//_______________________________________________________
+void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
+  //
+  // Connect the Input data with the task
+  //
+  AliTRDrecoTask::ConnectInputData(opt);
+  fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
 }
 
 //_______________________________________________________
@@ -61,128 +85,42 @@ void AliTRDcheckDetector::CreateOutputObjects(){
   //
   // Create Output Objects
   //
-  fContainer = new TObjArray();
-  TH1F *hNtrks = new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100);
-  TH1F *hNcls = new TH1F("hNcls", "Nr. of clusters per track", 30, 0, 180);
-  TH1F *hNtls = new TH1F("hNtls", "Nr. tracklets per track", 10, 0, 10);
-  TH1F *hNclTls = new TH1F("hNclTls","Mean Number of clusters per tracklet", 30, 0, 30);
-  TH1F *hChi2 = new TH1F("hChi2", "Chi2", 200, 0, 20);
-  TH1F *hChi2n = new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5);
-  TH1F *hSM = new TH1F("hSM", "Track Counts in Supermodule", 18, 0., 18.);
-  TH1F *hQcl = new TH1F("hQcl1", "Cluster charge", 200, -1200, 1200);
-  TProfile *hPH = new TProfile("hPH", "Average PH", 31, -0.5, 30.5);
-  TH1F * hQT = new TH1F("hQT", "Total Charge Deposit", 6000, 0, 6000);
-  TH1F * hEventsTrigger = new TH1F("hEventsTrigger", "Trigger Class", 20, 0, 20);
-  // Register Histograms
-  fContainer->Add(hNtrks);
-  fContainer->Add(hEventsTrigger);
-  fContainer->Add(hNcls);
-  fContainer->Add(hNtls);
-  fContainer->Add(hNclTls);
-  fContainer->Add(hChi2);
-  fContainer->Add(hChi2n);
-  fContainer->Add(hSM);
-  fContainer->Add(hPH);
-  fContainer->Add(hQcl);
-  fContainer->Add(hQT);
-
-  // Initialise the PHS, cluster charge and total charge deposit histograms for each detector
-  fPHSdetector = new TObjArray();
-  fPHSdetector->SetName("fPHSdetector");
-  fQCLdetector = new TObjArray();
-  fQCLdetector->SetName("fQCLdetector");
-  fQTdetector = new TObjArray();
-  fQTdetector->SetName("fQTdetector");
-  for(Int_t idet = 0; idet < kNDetectors; idet++){
-    fPHSdetector->Add(new TProfile(Form("hPHSdet%d", idet), Form("Average Pulse Height in Detector %d", idet), 31, -0.5, 30.5));
-    fQCLdetector->Add(new TH1F(Form("hQd%d",idet),Form("Qd%d",idet), 100,0,5000));
-    fQTdetector->Add(new TH1F(Form("hQTd%d", idet), Form("Qtotal%d", idet), 6000, 0, 6000));
-  }
-  fContainer->Add(fPHSdetector);
-  fContainer->Add(fQCLdetector);
-  fContainer->Add(fQTdetector);
-  
-  // Initialise the PHS, cluster charge and total charge deposit histograms for each sector
-  fPHSsector = new TObjArray();
-  fPHSsector->SetName("fPHSsector");
-  fQCLsector = new TObjArray();
-  fQCLsector->SetName("fQCLsector");
-  fQTsector = new TObjArray();
-  fQTsector->SetName("fQTsector");
-  for(Int_t isec = 0; isec < kNSectors; isec++){
-    fPHSsector->Add(new TProfile(Form("hPHSsec%d", isec), Form("Average Pulse Height in Sector %d", isec), 31, -0.5, 30.5));
-    fQCLsector->Add(new TH1F(Form("hQd%d",isec),Form("Qd%d",isec), 100,0,5000));
-    fQTsector->Add(new TH1F(Form("hQTd%d", isec), Form("Qtotal%d", isec), 6000, 0, 6000));
-  }
-  fContainer->Add(fPHSsector);
-  fContainer->Add(fQCLsector);
-  fContainer->Add(fQTsector);
+  OpenFile(0,"RECREATE");
+  fContainer = Histos();
+  if(!fTriggerNames) fTriggerNames = new TMap();
 }
 
 //_______________________________________________________
-void AliTRDcheckDetector::Exec(Option_t *){
+void AliTRDcheckDetector::Exec(Option_t *opt){
   //
   // Execution function
   // Filling TRD quality histos
   //
+  if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return;        // For real data we select only physical events
+  AliTRDrecoTask::Exec(opt);  
   Int_t nTracks = 0;           // Count the number of tracks per event
-  Int_t triggermask = 0;
-  AliTRDtrackInfo *fTrackInfo = 0x0;
-  AliTRDtrackV1 *fTRDtrack = 0x0;
-  AliTRDseedV1 *fTracklet = 0x0;
-  AliTRDcluster *fTRDcluster = 0x0;
+  Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
+  TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
+  if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
+  dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
-    fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
-    if(!fTrackInfo || !(fTRDtrack = fTrackInfo->GetTRDtrack())) continue;
-    if(!triggermask) triggermask = fTrackInfo->GetTriggerCluster();
-    Int_t nclusters = fTRDtrack->GetNumberOfClusters();
-    Int_t ntracklets = fTRDtrack->GetNumberOfTracklets();
-    Float_t chi2 = fTRDtrack->GetChi2();
-    // Fill Histograms
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist))->Fill(nclusters);
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist))->Fill(ntracklets);
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2))->Fill(chi2);
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2Normalized))->Fill(chi2/static_cast<Float_t>(ntracklets));
-    // now loop over single tracklets
-    if(ntracklets > 2){
-      Int_t sector = -1;
-      for(Int_t ilayer = 0; ilayer < kNLayers; ilayer++){
-        if(!(fTracklet = fTRDtrack->GetTracklet(ilayer)) || !fTracklet->IsOK()) continue;
-        Float_t nClustersTracklet = fTracklet->GetN();
-        dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist))->Fill(nClustersTracklet);
-        if(nClustersTracklet){
-          Float_t Qtot = 0;
-          Int_t detector = -1;
-          for(Int_t itb = 0; itb < kNTimebins; itb++){
-            if(!(fTRDcluster = fTracklet->GetClusters(itb))) continue;
-            Int_t localtime        = fTRDcluster->GetLocalTimeBin();
-            Double_t clusterCharge = fTRDcluster->GetQ();
-            detector               = fTRDcluster->GetDetector();
-            sector                 = static_cast<Int_t>(detector/kNDetectorsSector);
-            Double_t absolute_charge = TMath::Abs(clusterCharge);
-            Qtot += absolute_charge;
-            // Fill Histograms
-            // 1st. PHS
-            dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight))->Fill(localtime, absolute_charge);
-            dynamic_cast<TProfile *>(fPHSdetector->UncheckedAt(detector))->Fill(localtime, absolute_charge);
-            dynamic_cast<TProfile *>(fPHSsector->UncheckedAt(sector))->Fill(localtime, absolute_charge);
-            // 2nd. cluster charge
-            dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge))->Fill(absolute_charge);
-            dynamic_cast<TH1F *>(fQCLdetector->UncheckedAt(detector))->Fill(absolute_charge);
-            dynamic_cast<TH1F *>(fQCLsector->UncheckedAt(sector))->Fill(absolute_charge);
-          }
-          // Fill total cluster charge histograms
-          dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit))->Fill(Qtot);
-          dynamic_cast<TH1F *>(fQTdetector->UncheckedAt(detector))->Fill(Qtot);
-          dynamic_cast<TH1F *>(fQTsector->UncheckedAt(sector))->Fill(Qtot);
-        }
-      }
-      dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist))->Fill(sector);
-    }
+    if(!fTracks->UncheckedAt(iti)) continue;
+    AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
+    if(!fTrackInfo->GetTrack()) continue;
     nTracks++;
   }
-  dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
-       dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
+  if(nTracks){
+    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
+    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
+  }
+  if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
+    fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
+    // also set the label for both histograms
+    TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
+    histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
+    histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
+    histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
+  }
   PostData(0, fContainer);
 }
 
@@ -195,84 +133,614 @@ void AliTRDcheckDetector::Terminate(Option_t *){
 
 //_______________________________________________________
 Bool_t AliTRDcheckDetector::PostProcess(){
-       //
-       // Do Postprocessing (for the moment set the number of Reference histograms)
-       //
-       
-       TH1 * histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
-       histo->GetXaxis()->SetTitle("Number of Tracks");
-       histo->GetYaxis()->SetTitle("Events");
-       // Calculate percentage of events containing tracks
-       Int_t entries = (Int_t)histo->GetEntries();
-  Float_t percTracks = entries > 0 ? (entries - histo->GetBinContent(1))/Float_t(entries) : 0.;
-       printf("Percentage of Events containing TRD Tracks: %6.3f %%\n", 1.E2*percTracks);
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
-       histo->GetXaxis()->SetTitle("Number of Clusters");
-       histo->GetYaxis()->SetTitle("Entries");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
-       histo->GetXaxis()->SetTitle("Number of Tracklets");
-       histo->GetYaxis()->SetTitle("Entries");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
-       histo->GetXaxis()->SetTitle("Number of Clusters");
-       histo->GetYaxis()->SetTitle("Entries");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
-       histo->GetXaxis()->SetTitle("#chi^2");
-       histo->GetYaxis()->SetTitle("Entries");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
-       histo->GetXaxis()->SetTitle("Sector");
-       histo->GetYaxis()->SetTitle("Number of Tracks");
-       histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
-       histo->GetXaxis()->SetTitle("Time / 100ns");
-       histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
-       histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
-       histo->GetYaxis()->SetTitle("Entries");
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
-       histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
-       histo->GetYaxis()->SetTitle("Entries");
-       // Calculate the Percentage of events containing tracks as function of the trigger cluster
-       histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
-       Double_t nEvents = histo->Integral();
-       TGraph *percentages = new TGraph(histo->GetNbinsX());
-       for(Int_t ibin = 0; ibin < histo->GetNbinsX(); ibin++)
-               percentages->SetPoint(ibin, ibin, histo->GetBinContent(histo->FindBin(ibin))/nEvents);
-       percentages->GetXaxis()->SetTitle("Trigger Cluster");
-       percentages->GetYaxis()->SetTitle("#%Events");
-       percentages->SetMarkerColor(kBlue);
-       percentages->SetMarkerStyle(22);
-       fContainer->Add(percentages);
-       fNRefFigures = 10;
-       return kTRUE;
+  //
+  // Do Postprocessing (for the moment set the number of Reference histograms)
+  //
+  
+  TH1 * histo = 0x0;
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
+  histo->GetXaxis()->SetTitle("Number of Tracks");
+  histo->GetYaxis()->SetTitle("Events");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
+  histo->GetXaxis()->SetTitle("Number of Clusters");
+  histo->GetYaxis()->SetTitle("Entries");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
+  histo->GetXaxis()->SetTitle("Number of Tracklets");
+  histo->GetYaxis()->SetTitle("Entries");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
+  histo->GetXaxis()->SetTitle("Number of Clusters");
+  histo->GetYaxis()->SetTitle("Entries");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
+  histo->GetXaxis()->SetTitle("#chi^2");
+  histo->GetYaxis()->SetTitle("Entries");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
+  histo->GetXaxis()->SetTitle("Sector");
+  histo->GetYaxis()->SetTitle("Number of Tracks");
+  histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
+  histo->GetXaxis()->SetTitle("Time / 100ns");
+  histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
+  histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
+  histo->GetYaxis()->SetTitle("Entries");
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
+  histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
+  histo->GetYaxis()->SetTitle("Entries");
+  
+  // Calculate the purity of the trigger clusters
+  histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
+  TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
+  histoTracks->Divide(histo);
+  Float_t purities[20], val = 0;
+  TString triggernames[20];
+  Int_t nTriggerClasses = 0;
+  for(Int_t ibin = 1; ibin <= histo->GetNbinsX(); ibin++){
+    if((val = histoTracks->GetBinContent(ibin))){
+      purities[nTriggerClasses] = val;
+      triggernames[nTriggerClasses] = histoTracks->GetXaxis()->GetBinLabel(ibin);
+      nTriggerClasses++;
+    }
+  }
+  TH1F *hTriggerInf = new TH1F("fTriggerInf", "Trigger Information", TMath::Max(nTriggerClasses, 1), 0, TMath::Max(nTriggerClasses, 1));
+  for(Int_t ibin = 1; ibin <= nTriggerClasses; ibin++){
+    hTriggerInf->SetBinContent(ibin, purities[ibin-1]);
+    hTriggerInf->GetXaxis()->SetBinLabel(ibin, triggernames[ibin-1].Data());
+  }
+  hTriggerInf->GetXaxis()->SetTitle("Trigger Cluster");
+  hTriggerInf->GetYaxis()->SetTitle("Ratio");
+  hTriggerInf->GetYaxis()->SetRangeUser(0,1);
+//     hTriggerInf->SetMarkerColor(kBlue);
+//     hTriggerInf->SetMarkerStyle(22);
+  fContainer->Add(hTriggerInf);
+  fNRefFigures = 10;
+  return kTRUE;
+}
+
+//_______________________________________________________
+void AliTRDcheckDetector::GetRefFigure(Int_t ifig){
+  //
+  // Setting Reference Figures
+  //
+  TH1 *h = 0x0;
+  switch(ifig){
+  case 0:      
+    ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
+    break;
+  case 1:
+    ((TH1F*)fContainer->At(kNclustersHist))->Draw("pl");
+    break;
+  case 2:
+    h = (TH1F*)fContainer->At(kNtrackletsHist);
+    if(!h->GetEntries()) break;
+    h->Scale(100./h->Integral());
+    h->GetXaxis()->SetRangeUser(.5, 6.5);
+    h->SetFillColor(kGreen);
+    h->SetBarOffset(.2);
+    h->SetBarWidth(.6);
+    h->Draw("bar1");
+    break;
+  case 3:
+    ((TH1F*)fContainer->At(kNclusterTrackletHist))->Draw("pc");
+    break;
+  case 4:
+    ((TH1F*)fContainer->At(kChi2))->Draw("");
+    break;
+  case 5:
+    h = (TH1F*)fContainer->At(kNTracksSectorHist);
+    if(!h->GetEntries()) break;
+    h->Scale(100./h->Integral());
+    h->SetFillColor(kGreen);
+    h->SetBarOffset(.2);
+    h->SetBarWidth(.6);
+    h->Draw("bar1");
+    break;
+  case 6:
+    h = (TH1F*)fContainer->At(kPulseHeight);
+    h->SetMarkerStyle(24);
+    h->Draw("e1");
+    break;
+  case 7:
+    ((TH1F*)fContainer->At(kClusterCharge))->Draw("c");
+    break;
+  case 8:
+    ((TH1F*)fContainer->At(kChargeDeposit))->Draw("c");
+    break;
+  case 9: 
+    h=(TH1F*)fContainer->At(kPurity);
+    h->SetBarOffset(.2);
+    h->SetBarWidth(.6);
+    h->SetFillColor(kGreen);
+    h->Draw("bar1");
+    break;
+  default:
+    ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
+    break;
+  }
+}
+
+//_______________________________________________________
+TObjArray *AliTRDcheckDetector::Histos(){
+  //
+  // Create QA histograms
+  //
+  if(fContainer) return fContainer;
+  
+  fContainer = new TObjArray(25);
+  // Register Histograms
+  fContainer->AddAt(new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100), kNTracksEventHist);
+  fContainer->AddAt(new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100), kNEventsTriggerTracks);
+  fContainer->AddAt(new TH1F("hNcls", "Nr. of clusters per track", 181, -0.5, 180.5), kNclustersHist);
+  fContainer->AddAt(new TH1F("hNtls", "Nr. tracklets per track", 7, -0.5, 6.5), kNtrackletsHist);
+  fContainer->AddAt(new TH1F("hNtlsFindable", "Ratio of found/findable Tracklets" , 11, -0.05, 1.05), kNTrackletsVsFindable);
+  fContainer->AddAt(new TH1F("hNclTls","Mean Number of clusters per tracklet", 31, -0.5, 30.5), kNclusterTrackletHist);
+  fContainer->AddAt(new TH1F("hChi2", "Chi2", 200, 0, 20), kChi2);
+  fContainer->AddAt(new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5), kChi2Normalized);
+  fContainer->AddAt(new TH1F("hSM", "Track Counts in Supermodule", 18, -0.5, 17.5), kNTracksSectorHist);
+  // Detector signal on Detector-by-Detector basis
+  fContainer->AddAt(new TProfile("hPHdetector", "Average PH", 31, -0.5, 30.5), kPulseHeight);
+  fContainer->AddAt(new TH1F("hQclDetector", "Cluster charge", 200, 0, 1200), kClusterCharge);
+  fContainer->AddAt(new TH1F("hQTdetector", "Total Charge Deposit", 6000, 0, 6000), kChargeDeposit);
+  fContainer->AddAt(new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100), kNEventsTrigger);
+
+  return fContainer;
+}
+
+/*
+* Plotting Functions
+*/
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotMeanNClusters(const AliTRDtrackV1 *track){
+  //
+  // Plot the mean number of clusters per tracklet
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclusterTrackletHist)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  AliTRDseedV1 *tracklet = 0x0;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
+    h->Fill(tracklet->GetN());
+  }
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotNClusters(const AliTRDtrackV1 *track){
+  //
+  // Plot the number of clusters in one track
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersHist)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  
+  Int_t nclusters = 0;
+  AliTRDseedV1 *tracklet = 0x0;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
+    nclusters += tracklet->GetN();
+    if(fDebugLevel > 2){
+      Int_t crossing = tracklet->GetNChange();
+      AliTRDcluster *c = 0x0;
+      for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
+        if(!(c = tracklet->GetClusters(itime))) continue;
+        break;
+      }
+      Int_t detector = c->GetDetector();
+      Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
+      Float_t theta = TMath::ATan(tracklet->GetZfit(1));
+      Float_t phi = TMath::ATan(tracklet->GetYfit(1));
+      Float_t momentum = 0.;
+      Int_t pdg = 0;
+      Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
+      UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
+      if(fMC){
+        if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
+        pdg = fMC->GetPDG();
+      }
+      (*fDebugStream) << "NClusters"
+        << "Detector="  << detector
+        << "Sector="    << sector
+        << "crossing="  << crossing
+        << "momentum=" << momentum
+        << "pdg="                              << pdg
+        << "theta="                    << theta
+        << "phi="                              << phi
+        << "kinkIndex="        << kinkIndex
+        << "TPCncls="          << TPCncls
+        << "nclusters=" << nclusters
+        << "\n";
+    }
+  }
+  h->Fill(nclusters);
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotChi2(const AliTRDtrackV1 *track){
+  //
+  // Plot the chi2 of the track
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  h->Fill(fTrack->GetChi2());
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotNormalizedChi2(const AliTRDtrackV1 *track){
+  //
+  // Plot the chi2 of the track
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2Normalized)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  Int_t nTracklets = 0;
+  AliTRDseedV1 *tracklet = 0x0;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
+    nTracklets++;
+  }
+  h->Fill(fTrack->GetChi2()/nTracklets);
+  return h;
+}
+
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotNTracklets(const AliTRDtrackV1 *track){
+  //
+  // Plot the number of tracklets
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsHist)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  Int_t nTracklets = GetNTracklets(track);
+  h->Fill(nTracklets);
+  if(fDebugLevel > 3){
+    if(nTracklets == 1){
+      // If we have one Tracklet, check in which layer this happens
+      Int_t layer = -1;
+      AliTRDseedV1 *tracklet = 0x0;
+      for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+        if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
+      }
+      (*fDebugStream) << "PlotNTracklets"
+        << "Layer=" << layer
+        << "\n";
+    }
+  }
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotTrackletsVsFindable(const AliTRDtrackV1 *track){
+  //
+  // Plots the ratio of number of tracklets vs.
+  // number of findable tracklets
+  //
+  // Findable tracklets are defined as track prolongation
+  // to layer i does not hit the dead area +- epsilon
+  //
+  const Float_t epsilon = 0.01;   // dead area tolerance
+  const Float_t epsilon_R = 1;
+  Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTrackletsVsFindable)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  Int_t nFound = 0, nFindable = 0;
+  Int_t stack = -1;
+  Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
+  Double_t y = 0., z = 0.;
+  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDpadPlane *pp;  
+  AliTRDgeometry *fGeo = new AliTRDgeometry;
+  for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+    if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
+      tracklet->SetReconstructor(fReconstructor);
+      nFound++;
+      // printout tracklet
+    }
+  }
+  // 2 Different cases:
+  // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
+  // 2nd barrel track: here we propagate the track to the layers
+  AliTrackPoint points[6];
+  Float_t xyz[3];
+  memset(xyz, 0, sizeof(Float_t) * 3);
+  if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
+    // stand alone track
+    for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+      xyz[0] = x_anode[il];
+      points[il].SetXYZ(xyz);
+    }
+    AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
+  } else {
+    // barrel track
+    for(Int_t il = AliTRDgeometry::kNlayer; il--;){
+      // The track points have to be in reverse order for the Kalman Filter (since we move back)
+      xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
+      points[il].SetXYZ(xyz);
+    }
+    for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
+      printf("%d. X = %f\n", ipt, points[ipt].GetX());
+    AliTRDtrackerV1::FitKalman(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kFALSE, 6, points);
+    AliTrackPoint tempPointCont[AliTRDgeometry::kNlayer];
+    memcpy(tempPointCont, points, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
+    for(Int_t il = AliTRDgeometry::kNlayer; il--; ){
+      tempPointCont[il].GetXYZ(xyz);
+      points[AliTRDgeometry::kNlayer - il -1].SetXYZ(xyz);
+    }
+  }
+  for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+    y = points[il].GetY();
+    z = points[il].GetZ();
+    if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
+    pp = fGeo->GetPadPlane(il, stack);
+    ymin = pp->GetCol0() + epsilon;
+    ymax = pp->GetColEnd() - epsilon; 
+    zmin = pp->GetRowEnd() + epsilon; 
+    zmax = pp->GetRow0() - epsilon;
+    // ignore y-crossing (material)
+    if((z > zmin && z < zmax) && (y > ymin && y < ymax)) nFindable++;
+/*
+      if(fDebugLevel > 3){
+        Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetMeanz() : 0};
+        Int_t hasTracklet = tracklet ? 1 : 0;
+        (*fDebugStream)   << "GetFindableTracklets"
+          << "layer="     << il
+          << "ytracklet=" << pos_tracklet[0]
+          << "ytrack="    << y
+          << "ztracklet=" << pos_tracklet[1]
+          << "ztrack="    << z
+          << "tracklet="  << hasTracklet
+          << "\n";
+      }
+*/
+  }
+  
+  delete fGeo;
+  h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
+  if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotPulseHeight(const AliTRDtrackV1 *track){
+  //
+  // Plot the average pulse height
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TProfile *h = 0x0;
+  if(!(h = dynamic_cast<TProfile *>(fContainer->At(kPulseHeight)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDcluster *c = 0x0;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
+    for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
+      if(!(c = tracklet->GetClusters(itime))) continue;
+      Int_t localtime        = c->GetLocalTimeBin();
+      Double_t absolute_charge = TMath::Abs(c->GetQ());
+      h->Fill(localtime, absolute_charge);
+      if(fDebugLevel > 3){
+        Int_t crossing = tracklet->GetNChange();
+        Int_t detector = c->GetDetector();
+        Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
+        Float_t theta = TMath::ATan(tracklet->GetZfit(1));
+        Float_t phi = TMath::ATan(tracklet->GetYfit(1));
+        Float_t momentum = 0.;
+        Int_t pdg = 0;
+        Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
+        UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
+        if(fMC){
+          if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
+          pdg = fMC->GetPDG();
+        }
+        (*fDebugStream) << "PulseHeight"
+          << "Detector="       << detector
+          << "Sector="         << sector
+          << "crossing="       << crossing
+          << "Timebin="                << localtime
+          << "Charge="         << absolute_charge
+          << "momentum="       << momentum
+          << "pdg="                            << pdg
+          << "theta="                  << theta
+          << "phi="                            << phi
+          << "kinkIndex="      << kinkIndex
+          << "TPCncls="                << TPCncls
+          << "\n";
+      }
+    }
+  }
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotClusterCharge(const AliTRDtrackV1 *track){
+  //
+  // Plot the cluster charge
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kClusterCharge)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDcluster *c = 0x0;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
+    for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
+      if(!(c = tracklet->GetClusters(itime))) continue;
+      h->Fill(c->GetQ());
+    }
+  }
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotChargeDeposit(const AliTRDtrackV1 *track){
+  //
+  // Plot the charge deposit per chamber
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeDeposit)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDcluster *c = 0x0, *c1 = 0x0;   // c1 for the Debug Stream
+  Double_t Qtot = 0;
+  Int_t nTracklets = 0;
+  if(fDebugLevel > 3)
+    nTracklets = GetNTracklets(fTrack); // fill NTracklet to the Debug Stream
+  for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
+    Qtot = 0;
+    c1 = 0x0;
+    for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
+      if(!(c = tracklet->GetClusters(itime))) continue;
+      if(!c1) c1 = c;
+      Qtot += TMath::Abs(c->GetQ());
+    }
+    h->Fill(Qtot);
+    if(fDebugLevel > 3){
+      Int_t crossing = tracklet->GetNChange();
+      Int_t detector = c1->GetDetector();
+      Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
+      Float_t theta = TMath::ATan(tracklet->GetZfit(1));
+      Float_t phi = TMath::ATan(tracklet->GetYfit(1));
+      Float_t momentum = 0.;
+      Int_t pdg = 0;
+      Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
+      UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
+      if(fMC){
+             if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
+        pdg = fMC->GetPDG();
+      }
+      (*fDebugStream) << "ChargeDeposit"
+        << "Detector="  << detector
+        << "Sector="    << sector
+        << "crossing="  << crossing
+        << "momentum=" << momentum
+        << "nTracklets="<< nTracklets
+        << "pdg="                              << pdg
+        << "theta="                    << theta
+        << "phi="                              << phi
+        << "kinkIndex="        << kinkIndex
+        << "TPCncls="          << TPCncls
+        << "QT="        << Qtot
+        << "\n";
+    }
+  }
+  return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotTracksSector(const AliTRDtrackV1 *track){
+  //
+  // Plot the number of tracks per Sector
+  //
+  if(track) fTrack = track;
+  if(!fTrack){
+    AliWarning("No Track defined.");
+    return 0x0;
+  }
+  TH1 *h = 0x0;
+  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTracksSectorHist)))){
+    AliWarning("No Histogram defined.");
+    return 0x0;
+  }
+  AliTRDseedV1 *tracklet = 0x0;
+  AliTRDcluster *c = 0x0;
+  Int_t sector = -1;
+  for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
+    if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
+    for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
+      if(!(c = tracklet->GetClusters(itime))) continue;
+      sector = static_cast<Int_t>(c->GetDetector()/AliTRDgeometry::kNdets);
+    }
+    break;
+  }
+  h->Fill(sector);
+  return h;
 }
 
 //_______________________________________________________
-void AliTRDcheckDetector::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt){
-       //
-       // Setting Reference Figures
-       //
-       switch(ifig){
-               case 0: first = last = kNTracksEventHist;
-                                               break;
-               case 1: first = last = kNclustersHist;
-                                               break;
-               case 2: first = last = kNtrackletsHist;
-                                               break;
-               case 3: first = last = kNclusterTrackletHist;
-                                               break;
-               case 4: first = last = kChi2;
-                                               break;
-               case 5: first = last = kNTracksSectorHist;
-                                               break;
-               case 6: first = last = kPulseHeight;
-                                               break;
-               case 7: first = last = kClusterCharge;
-                                               break;
-               case 8: first = last = kChargeDeposit;
-                                               break;
-               case 9: first = last = kPostProcessing;
-                                               break;
-               default: first = last = kNTracksEventHist;
-                                               break;
-       };
+Int_t AliTRDcheckDetector::GetNTracklets(const AliTRDtrackV1 *track){
+  //
+  // Count the number of tracklets per track
+  //
+  if(!track) return 0;
+  Int_t nTracklets = 0;
+  AliTRDseedV1 *tracklet = 0x0;
+  for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+    if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()) nTracklets++;
+  }
+  return nTracklets;
 }
 
+//________________________________________________________
+void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
+{
+
+  fReconstructor->SetRecoParam(r);
+}