using AliCDBConnect task in test train of TRD
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcheckDET.cxx
index 8c003ac..f628f2c 100644 (file)
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
+#include <TArrayD.h>
 #include <TAxis.h>
 #include <TCanvas.h>
 #include <TFile.h>
 #include <TH1F.h>
 #include <TH1I.h>
 #include <TH2F.h>
+#include <TH3S.h>
+#include <TH3F.h>
 #include <TF1.h>
 #include <TGaxis.h>
 #include <TGraph.h>
@@ -37,6 +40,7 @@
 #include <TLinearFitter.h>
 #include <TMath.h>
 #include <TMap.h>
+#include <TProfile2D.h>
 #include <TObjArray.h>
 #include <TObject.h>
 #include <TObjString.h>
 
 #include "info/AliTRDtrackInfo.h"
 #include "info/AliTRDeventInfo.h"
+#include "AliTRDinfoGen.h"
 #include "AliTRDcheckDET.h"
+#include "AliTRDpwg1Helper.h"
 
 #include <cstdio>
 #include <iostream>
 
 ClassImp(AliTRDcheckDET)
 
+const Color_t AliTRDcheckDET::fkColorsCentrality[AliTRDeventInfo::kCentralityClasses] = {kTeal, kOrange, kGreen, kBlue ,kRed};
+
 //_______________________________________________________
 AliTRDcheckDET::AliTRDcheckDET():
   AliTRDrecoTask()
-  ,fEventInfo(NULL)
+  ,fCentralityClass(-1)
   ,fTriggerNames(NULL)
-  ,fReconstructor(NULL)
-  ,fGeo(NULL)
   ,fFlags(0)
 {
   //
   // Default constructor
   //
-  SetNameTitle("checkDET", "Basic TRD data checker");
+  SetNameTitle("TRDcheckDET", "Basic TRD data checker");
 }
 
 //_______________________________________________________
 AliTRDcheckDET::AliTRDcheckDET(char* name):
   AliTRDrecoTask(name, "Basic TRD data checker")
-  ,fEventInfo(NULL)
+  ,fCentralityClass(-1)
   ,fTriggerNames(NULL)
-  ,fReconstructor(NULL)
-  ,fGeo(NULL)
   ,fFlags(0)
 {
   //
   // Default constructor
   //
-  DefineInput(2, AliTRDeventInfo::Class());
-
-  fReconstructor = new AliTRDReconstructor;
-  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
-  fGeo = new AliTRDgeometry;
   InitFunctorList();
 }
 
@@ -116,8 +115,6 @@ AliTRDcheckDET::~AliTRDcheckDET(){
   // Destructor
   // 
   if(fTriggerNames) delete fTriggerNames;
-  delete fReconstructor;
-  delete fGeo;
 }
 
 //_______________________________________________________
@@ -125,9 +122,7 @@ void AliTRDcheckDET::UserCreateOutputObjects(){
   //
   // Create Output Objects
   //
-  if(!HasFunctorList()) InitFunctorList();
-  OpenFile(1,"RECREATE");
-  fContainer = Histos();
+  AliTRDrecoTask::UserCreateOutputObjects();
   if(!fTriggerNames) fTriggerNames = new TMap();
 }
 
@@ -137,226 +132,251 @@ void AliTRDcheckDET::UserExec(Option_t *opt){
   // Execution function
   // Filling TRD quality histos
   //
-
-  fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
-  if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return;        // For real data we select only physical events
+  AliDebug(2, Form("EventInfo[%p] Header[%p]", (void*)fEvent, (void*)(fEvent?fEvent->GetEventHeader():NULL)));
+  if(fEvent) fCentralityClass = fEvent->GetCentrality();
+  else fCentralityClass = -1;  // Assume pp
   AliTRDrecoTask::UserExec(opt);  
-  Int_t nTracks = 0;           // Count the number of tracks per event
-  Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
-  TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
-  AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
-  dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
+
+  TH1F *histo(NULL); AliTRDtrackInfo *fTrackInfo(NULL); Int_t nTracks(0);              // Count the number of tracks per event
   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
     if(!fTracks->UncheckedAt(iti)) continue;
-    AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
+    if(!(fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti)))) continue;
     if(!fTrackInfo->GetTrack()) continue;
     nTracks++;
   }
+  if(nTracks)
+    if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent)))) histo->Fill(nTracks);
+
+  if(!fEvent->GetEventHeader()) return; // For trigger statistics event header is essential
+  Int_t triggermask = fEvent->GetEventHeader()->GetTriggerMask();
+  TString triggername =  fEvent->GetRunInfo()->GetFiredTriggerClasses(triggermask);
+  AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
+  if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger)))) histo->Fill(triggermask);
 
   if(nTracks){
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
-    dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
+    if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks)))) histo->Fill(triggermask);
   }
   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);
+    if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))))
+      histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
+    if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))))
+      histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
   }
-  PostData(1, fContainer);
 }
 
 
 //_______________________________________________________
+//_______________________________________________________
 Bool_t AliTRDcheckDET::PostProcess(){
   //
   // Do Postprocessing (for the moment set the number of Reference histograms)
   //
   
-  TH1 * h = NULL;
-  
+  TH1F *h(NULL), *h1(NULL);
+
   // Calculate of the trigger clusters purity
-  h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
-  TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
-  h1->Divide(h);
-  Float_t purities[20], val = 0;
-  TString triggernames[20];
-  Int_t nTriggerClasses = 0;
-  for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
-    if((val = h1->GetBinContent(ibin))){
-      purities[nTriggerClasses] = val;
-      triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
-      nTriggerClasses++;
+  if((h  = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"))) &&
+     (h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks")))) {
+    h1->Divide(h);
+    Float_t purities[20], val = 0; memset(purities, 0, 20*sizeof(Float_t));
+    TString triggernames[20];
+    Int_t nTriggerClasses = 0;
+    for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
+      if((val = h1->GetBinContent(ibin))){
+        purities[nTriggerClasses] = val;
+        triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
+        nTriggerClasses++;
+      }
+    }
+
+    if((h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity)))){
+      TAxis *ax = h->GetXaxis();
+      for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
+        h->Fill(itrg, purities[itrg]);
+        ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
+      }
+      ax->SetRangeUser(-0.5, nTriggerClasses+.5);
+      h->GetYaxis()->SetRangeUser(0,1);
     }
   }
-  h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
-  TAxis *ax = h->GetXaxis();
-  for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
-    h->Fill(itrg, purities[itrg]);
-    ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
-  }
-  ax->SetRangeUser(-0.5, nTriggerClasses+.5);
-  h->GetYaxis()->SetRangeUser(0,1);
 
   // track status
-  h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
-  Float_t ok = h->GetBinContent(1);
-  Int_t nerr = h->GetNbinsX();
-  for(Int_t ierr=nerr; ierr--;){
-    h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
+  if((h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus)))){
+    Float_t ok = h->GetBinContent(1);
+    Int_t nerr = h->GetNbinsX();
+    for(Int_t ierr=nerr; ierr--;){
+      h->SetBinContent(ierr+1, ok>0.?1.e2*h->GetBinContent(ierr+1)/ok:0.);
+    }
+    h->SetBinContent(1, 0.);
   }
-  h->SetBinContent(1, 0.);
-
   // tracklet status
   
-  TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
-  for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
-    h=dynamic_cast<TH1F*>(arr->At(ily));
-    Float_t okB = h->GetBinContent(1);
-    Int_t nerrB = h->GetNbinsX();
-    for(Int_t ierr=nerrB; ierr--;){
-      h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
+  TObjArray *arr(NULL);
+  if(( arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus)))){
+    for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
+      if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
+      Float_t okB = h->GetBinContent(1);
+      Int_t nerrB = h->GetNbinsX();
+      for(Int_t ierr=nerrB; ierr--;){
+        h->SetBinContent(ierr+1, okB>0.?1.e2*h->GetBinContent(ierr+1)/okB:0.);
+      }
+      h->SetBinContent(1, 0.);
     }
-    h->SetBinContent(1, 0.);
   }
-
   fNRefFigures = 17;
 
   return kTRUE;
 }
 
-//_______________________________________________________
-Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
+void AliTRDcheckDET::MakeSummary(){
   //
-  // Setting Reference Figures
+  // Create summary plots for TRD check DET
+  // This function creates 2 summary plots:
+  // - General Quantities
+  // - PHS
+  // The function will reuse GetRefFigure
   //
-  enum FigureType_t{
-    kFigNclustersTrack,
-    kFigNclustersTracklet,
-    kFigNtrackletsTrack,
-    kFigNTrackletsP,
-    kFigNtrackletsCross,
-    kFigNtrackletsFindable,
-    kFigNtracksEvent,
-    kFigNtracksSector,
-    kFigTrackStatus,
-    kFigTrackletStatus,
-    kFigChi2,
-    kFigPH,
-    kFigChargeCluster,
-    kFigChargeTracklet,
-    kFigNeventsTrigger,
-    kFigNeventsTriggerTracks,
-    kFigTriggerPurity
-  };
-  gPad->SetLogy(0);
-  gPad->SetLogx(0);
-  TH1 *h = NULL; TObjArray *arr=NULL;
-  TLegend *leg = NULL;
-  Bool_t kFIRST(1);
-  switch(ifig){
-  case kFigNclustersTrack:
-    (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
-    PutTrendValue("NClustersTrack", h->GetMean());
-    PutTrendValue("NClustersTrackRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNclustersTracklet:
-    (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
-    PutTrendValue("NClustersTracklet", h->GetMean());
-    PutTrendValue("NClustersTrackletRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNtrackletsTrack:
-    h=MakePlotNTracklets();
-    PutTrendValue("NTrackletsTrack", h->GetMean());
-    PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNTrackletsP:
-    MakePlotnTrackletsVsP();
-    return kTRUE;
-  case kFigNtrackletsCross:
-    h = (TH1F*)fContainer->FindObject("hNtlsCross");
-    if(!MakeBarPlot(h, kRed)) break;
-    PutTrendValue("NTrackletsCross", h->GetMean());
-    PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNtrackletsFindable:
-    h = (TH1F*)fContainer->FindObject("hNtlsFindable");
-    if(!MakeBarPlot(h, kGreen)) break;
-    PutTrendValue("NTrackletsFindable", h->GetMean());
-    PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNtracksEvent:
-    (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
-    PutTrendValue("NTracksEvent", h->GetMean());
-    PutTrendValue("NTracksEventRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNtracksSector:
-    h = (TH1F*)fContainer->FindObject("hNtrksSector");
-    if(!MakeBarPlot(h, kGreen)) break;
-    PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
-    return kTRUE;
-  case kFigTrackStatus:
-    if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
-    h->GetXaxis()->SetRangeUser(0.5, -1);
-    h->GetYaxis()->CenterTitle();
-    h->Draw("c");
-    PutTrendValue("TrackStatus", h->Integral());
-    gPad->SetLogy(0);
-    return kTRUE;
-  case kFigTrackletStatus:
-    if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
-    leg = new TLegend(.68, .7, .97, .97);
-    leg->SetBorderSize(0);leg->SetFillStyle(0);
-    leg->SetHeader("TRD layer");
-    for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
-      if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
-      if(kFIRST){
-        h->Draw("pl");
-        h->GetXaxis()->SetRangeUser(0.5, -1);
-        h->GetYaxis()->CenterTitle();
-        kFIRST = kFALSE;
-      } else h->Draw("samepl");
-      leg->AddEntry(h, Form("ly = %d", ily), "l");
-      PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
-    }
-    leg->Draw();
-    gPad->SetLogy(0);
-    return kTRUE;
-  case kFigChi2:
-    MakePlotChi2();
-    return kTRUE;
-  case kFigPH:
-    gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
-    MakePlotPulseHeight();
-    gPad->SetLogy(0);
-    return kTRUE;
-  case kFigChargeCluster:
-    (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
-    gPad->SetLogy(1);
-    PutTrendValue("ChargeCluster", h->GetMaximumBin());
-    PutTrendValue("ChargeClusterRMS", h->GetRMS());
-    return kTRUE;
-  case kFigChargeTracklet:
-    (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
-    PutTrendValue("ChargeTracklet", h->GetMaximumBin());
-    PutTrendValue("ChargeTrackletRMS", h->GetRMS());
-    return kTRUE;
-  case kFigNeventsTrigger:
-    ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
-    return kTRUE;
-  case kFigNeventsTriggerTracks:
-    ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
-    return kTRUE;
-  case kFigTriggerPurity: 
-    if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
-    break;
-  default:
-    break;
-  }
-  AliInfo(Form("Reference plot [%d] missing result", ifig));
-  return kFALSE;
+  
+  TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);
+  cOut->Divide(3,3);
+
+  // Create figures using GetRefFigure
+  cOut->cd(1); GetRefFigure(kFigNtracksEvent);  
+  cOut->cd(2); GetRefFigure(kFigNtracksSector);
+  cOut->cd(3); GetRefFigure(kFigNclustersTrack);
+  cOut->cd(4); GetRefFigure(kFigNclustersTracklet);
+  cOut->cd(5); GetRefFigure(kFigNtrackletsTrack);
+  cOut->cd(6); GetRefFigure(kFigNTrackletsP);
+  cOut->cd(7); GetRefFigure(kFigChargeCluster);
+  cOut->cd(8); GetRefFigure(kFigChargeTracklet);
+  cOut->SaveAs("TRD_TrackQuality.gif");
+  delete cOut;
+
+  // Second Plot: PHS
+  cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 512);
+  cOut->cd(); GetRefFigure(kFigPH);
+  cOut->SaveAs("TRD_PH.gif"); 
+  delete cOut;
+
+  // Third Plot: Mean Number of clusters as function of eta, phi and layer
+   cOut = new TCanvas(Form("summary%s3", GetName()), Form("Summary 3 for task %s", GetName()), 1024, 768);
+  cOut->cd(); MakePlotMeanClustersLayer();
+  cOut->SaveAs("TRD_MeanNclusters.gif"); 
+  delete cOut;
+
+}
+
+//_______________________________________________________
+Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
+       //
+       // Setting Reference Figures
+       //
+       gPad->SetLogy(0);
+       gPad->SetLogx(0);
+       TH1 *h = NULL; TObjArray *arr=NULL;
+       TLegend *leg = NULL;
+       Bool_t kFIRST(1);
+       switch(ifig){
+       case kFigNclustersTrack:
+               MakePlotNclustersTrack();
+               return kTRUE;
+       case kFigNclustersTracklet:
+               MakePlotNclustersTracklet();
+               return kTRUE;
+       case kFigNtrackletsTrack:
+               h=MakePlotNTracklets();
+               if(h){
+                       PutTrendValue("NTrackletsTrack", h->GetMean());
+                       PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
+               }
+               return kTRUE;
+       case kFigNTrackletsP:
+               MakePlotnTrackletsVsP();
+               return kTRUE;
+       case kFigNtrackletsCross:
+               h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsCross"), -1);
+               if(!MakeBarPlot(h, kRed)) break;
+               PutTrendValue("NTrackletsCross", h->GetMean());
+               PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
+               return kTRUE;
+       case kFigNtrackletsFindable:
+               h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsFindable"), -1);
+               if(!MakeBarPlot(h, kGreen)) break;
+               PutTrendValue("NTrackletsFindable", h->GetMean());
+               PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
+               return kTRUE;
+       case kFigNtracksEvent:
+               (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
+               PutTrendValue("NTracksEvent", h->GetMean());
+               PutTrendValue("NTracksEventRMS", h->GetRMS());
+               return kTRUE;
+       case kFigNtracksSector:
+               h = (TH1F*)fContainer->FindObject("hNtrksSector");
+               if(!MakeBarPlot(h, kGreen)) break;
+               PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
+               return kTRUE;
+       case kFigTrackStatus:
+               if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
+               h->GetXaxis()->SetRangeUser(0.5, -1);
+               h->GetYaxis()->CenterTitle();
+               h->Draw("c");
+               PutTrendValue("TrackStatus", h->Integral());
+               gPad->SetLogy(0);
+               return kTRUE;
+       case kFigTrackletStatus:
+               if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
+               leg = new TLegend(.68, .7, .97, .97);
+               leg->SetBorderSize(0);leg->SetFillStyle(0);
+               leg->SetHeader("TRD layer");
+               for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
+                       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
+                       if(kFIRST){
+                               h->Draw("pl");
+                               h->GetXaxis()->SetRangeUser(0.5, -1);
+                               h->GetYaxis()->CenterTitle();
+                               kFIRST = kFALSE;
+                       } else h->Draw("samepl");
+                       leg->AddEntry(h, Form("ly = %d", ily), "l");
+                       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
+               }
+               leg->Draw();
+               gPad->SetLogy(0);
+               return kTRUE;
+       case kFigChi2:
+               MakePlotChi2();
+               return kTRUE;
+       case kFigPH:
+               gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
+               MakePlotPulseHeight();
+               gPad->SetLogy(0);
+               return kTRUE;
+       case kFigChargeCluster:
+               h = ProjectCentrality((TH2*)fContainer->FindObject("hQcl"), -1);
+               h->Draw("c");
+               gPad->SetLogy(1);
+               PutTrendValue("ChargeCluster", h->GetMaximumBin());
+               PutTrendValue("ChargeClusterRMS", h->GetRMS());
+               return kTRUE;
+       case kFigChargeTracklet:
+               MakePlotTrackletCharge();
+               return kTRUE;
+       case kFigNeventsTrigger:
+               ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
+               return kTRUE;
+       case kFigNeventsTriggerTracks:
+               ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
+               return kTRUE;
+       case kFigTriggerPurity:
+               if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
+               break;
+       default:
+               break;
+       }
+       AliInfo(Form("Reference plot [%d] missing result", ifig));
+       return kFALSE;
 }
 
 //_______________________________________________________
@@ -368,67 +388,70 @@ TObjArray *AliTRDcheckDET::Histos(){
   if(fContainer) return fContainer;
   
   fContainer = new TObjArray(20);
-  //fContainer->SetOwner(kTRUE);
+  fContainer->SetOwner(kTRUE);
 
   // Register Histograms
   TH1 * h = NULL;
+  TH2 * h2 = NULL;      // Pointer for two dimensional histograms
+  TH3 * h3 = NULL;      // Pointer for tree dimensional histograms
   TAxis *ax = NULL;
-  if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
-    h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
-    h->GetXaxis()->SetTitle("N_{clusters}");
-    h->GetYaxis()->SetTitle("Entries");
-  } else h->Reset();
-  fContainer->AddAt(h, kNclustersTrack);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hNcls"))){
+    h2 = new TH2F("hNcls", "N_{clusters}/track;N_{clusters};Centrality;Entries", 181, -0.5, 180.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h2->Reset();
+  fContainer->AddAt(h2, kNclustersTrack);
+  TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
+  arr->SetOwner(kTRUE);  arr->SetName("clusters");
+  fContainer->AddAt(arr, kNclustersLayer);
+  for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
+    if(!(h = (TProfile2D *)gROOT->FindObject(Form("hNcl%d", ily)))){
+      h = new TProfile2D(Form("hNcl%d", ily), Form("Mean Number of clusters in Layer %d", ily), 100, -1.0, 1.0, 50, -1.1*TMath::Pi(), 1.1*TMath::Pi());
+      h->GetXaxis()->SetTitle("#eta");
+      h->GetYaxis()->SetTitle("#phi");
+    } else h->Reset();
+    arr->AddAt(h, ily);
+  }
 
-  if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
-    h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
-    h->GetXaxis()->SetTitle("N_{clusters}");
-    h->GetYaxis()->SetTitle("Entries");
-  } else h->Reset();
-  fContainer->AddAt(h, kNclustersTracklet);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hNclTls"))){
+    h2 = new TH2F("hNclTls","N_{clusters}/tracklet;N_{clusters};Entries", 51, -0.5, 50.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h2->Reset();
+  fContainer->AddAt(h2, kNclustersTracklet);
 
-  if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
-    h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
-    h->GetXaxis()->SetTitle("N^{tracklet}");
-    h->GetYaxis()->SetTitle("freq. [%]");
-  } else h->Reset();
-  fContainer->AddAt(h, kNtrackletsTrack);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hNtls"))){
+    h2 = new TH2F("hNtls", "N_{tracklets}/track;N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h2->Reset();
+  fContainer->AddAt(h2, kNtrackletsTrack);
 
-  if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
-    h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
-    h->GetXaxis()->SetTitle("N^{tracklet}");
-    h->GetYaxis()->SetTitle("freq. [%]");
+  if(!(h = (TH2F *)gROOT->FindObject("htlsSTA"))){
+    h = new TH2F("hNtlsSTA", "#splitline{N_{tracklets}/track}{Stand Alone};N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
   }
   fContainer->AddAt(h, kNtrackletsSTA);
 
   // Binning for momentum dependent tracklet Plots
   const Int_t kNp(30);
-  Float_t P=0.2;
-  Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1];
-  for(Int_t i=0;i<kNp+1; i++,P+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=P;
+  Float_t p=0.2;
+  Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1], binsCent[AliTRDeventInfo::kCentralityClasses+2];
+  for(Int_t i=0;i<kNp+1; i++,p+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=p;
   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
-  if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
+  for(Int_t icent = -1; icent < AliTRDeventInfo::kCentralityClasses + 1; icent++) binsCent[icent+1] = icent - 0.5;
+  if(!(h3 = (TH3F *)gROOT->FindObject("htlsBAR"))){
     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
-    h = new TH2F("hNtlsBAR", 
-    "N_{tracklets} / track;p [GeV/c];N^{tracklet};freq. [%]", 
-    kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt);
+    h3 = new TH3F("hNtlsBAR", 
+    "N_{tracklets}/track;p [GeV/c];N^{tracklet};freq. [%]",
+    kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt, AliTRDeventInfo::kCentralityClasses + 1, binsCent);
   }
-  fContainer->AddAt(h, kNtrackletsBAR);
+  fContainer->AddAt(h3, kNtrackletsBAR);
 
   // 
-  if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
-    h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
-    h->GetXaxis()->SetTitle("n_{row cross}");
-    h->GetYaxis()->SetTitle("freq. [%]");
-  } else h->Reset();
-  fContainer->AddAt(h, kNtrackletsCross);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsCross"))){
+    h2 = new TH2F("hNtlsCross", "N_{tracklets}^{cross}/track;n_{row cross};Centrality;freq.[%]", 7, -0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h2->Reset();
+  fContainer->AddAt(h2, kNtrackletsCross);
 
-  if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
-    h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
-    h->GetXaxis()->SetTitle("r [a.u]");
-    h->GetYaxis()->SetTitle("Entries");
-  } else h->Reset();
-  fContainer->AddAt(h, kNtrackletsFindable);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsFindable"))){
+    h2 = new TH2F("hNtlsFindable", "Found/Findable Tracklets;r[a.u];Centrality;Entries" , 101, -0.005, 1.005,  AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h2->Reset();
+  fContainer->AddAt(h2, kNtrackletsFindable);
 
   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
@@ -454,7 +477,7 @@ TObjArray *AliTRDcheckDET::Histos(){
   }
   fContainer->AddAt(h, kTrackStatus);
 
-  TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
+  arr = new TObjArray(AliTRDgeometry::kNlayer);
   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
   fContainer->AddAt(arr, kTrackletStatus);
   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
@@ -471,48 +494,38 @@ TObjArray *AliTRDcheckDET::Histos(){
   }
 
   // <PH> histos
-  arr = new TObjArray(2);
+  arr = new TObjArray(3);
   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
   fContainer->AddAt(arr, kPH);
-  if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
-    h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
-    h->GetXaxis()->SetTitle("Time / 100ns");
-    h->GetYaxis()->SetTitle("<PH> [a.u]");
-  } else h->Reset();
-  arr->AddAt(h, 0);
-  if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
-    h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
-  else h->Reset();
-  arr->AddAt(h, 1);
+  if(!(h3 = (TH3F *)gROOT->FindObject("hPHx"))){
+    h3 = new TH3F("hPHx", "<PH>(x);x[mm];Centrality;Charge[a.u.]", 31, -0.08, 4.88, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h3->Reset();
+  arr->AddAt(h3, 0);
+  if(!(h3 = (TH3F *)gROOT->FindObject("hPHt"))){
+    h3 = new TH3F("hPHt", "<PH>(t);time[100ns];Centrality;Charge[a.u.]", 31, -0.5, 30.5, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h3->Reset();
+  arr->AddAt(h3, 1);
 
   // Chi2 histos
-  if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
-    h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
-    h->SetXTitle("ndf");
-    h->SetYTitle("#chi^{2}/ndf");
-    h->SetZTitle("entries");
-  } else h->Reset();
-  fContainer->AddAt(h, kChi2);
+  if(!(h3 = (TH3S*)gROOT->FindObject("hChi2"))){
+    h3 = new TH3S("hChi2", "#chi^{2}/track;ndf;#chi^{2}/ndf;Centrality", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  } else h3->Reset();
+  fContainer->AddAt(h3, kChi2);
 
-  if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
-    h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
-    h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
-    h->GetYaxis()->SetTitle("Entries");
-  }else h->Reset();
-  fContainer->AddAt(h, kChargeCluster);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hQcl"))){
+    h2 = new TH2F("hQcl", "Q_{cluster};Charge[a.u.];Centrality;Entries", 200, 0, 1200, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  }else h2->Reset();
+  fContainer->AddAt(h2, kChargeCluster);
 
-  if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
-    h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
-    h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
-    h->GetYaxis()->SetTitle("Entries");
-  }else h->Reset();
-  fContainer->AddAt(h, kChargeTracklet);
+  if(!(h2 = (TH2F *)gROOT->FindObject("hQtrklt"))){
+    h2 = new TH2F("hQtrklt", "Q_{tracklet};Charge[a.u.];Centrality;Entries", 6000, 0, 6000, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
+  }else h2->Reset();
+  fContainer->AddAt(h2, kChargeTracklet);
 
 
   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
   else h->Reset();
-  printf("Histos Adding \n");
   
   fContainer->AddAt(h, kNeventsTrigger);
 
@@ -531,6 +544,18 @@ TObjArray *AliTRDcheckDET::Histos(){
   return fContainer;
 }
 
+//_______________________________________________________
+TH1 *AliTRDcheckDET::ProjectCentrality(TH2 *hIn, Int_t centralityBin){
+  //
+  // Project histogram to a given centrality Bin
+  //
+  if(!hIn) return NULL;
+  if(centralityBin >= AliTRDeventInfo::kCentralityClasses) centralityBin = -1;
+  Int_t binMin = centralityBin > -1 ? centralityBin + 2 : 0, 
+        binMax = centralityBin > -1 ? centralityBin + 2 : -1;
+  return hIn->ProjectionX(Form("%s_%d", hIn->GetName(), centralityBin), binMin, binMax);
+}
+
 /*
 * Plotting Functions
 */
@@ -554,7 +579,7 @@ TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
 //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
   TH1 *h = NULL;
@@ -586,7 +611,7 @@ TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
 //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
   TObjArray *arr =NULL;
@@ -613,18 +638,31 @@ TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
+  AliExternalTrackParam *par = fkTrack->GetTrackIn() ? fkTrack->GetTrackIn() : fkTrack->GetTrackOut();
+  TH2 *h = NULL;
+  TProfile2D *hlayer = NULL;
+  Double_t eta = 0., phi = 0.;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTracklet)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
   AliTRDseedV1 *tracklet = NULL;
+  TObjArray *histosLayer = dynamic_cast<TObjArray *>(fContainer->At(kNclustersLayer));
+  if(!histosLayer){
+    AliWarning("No Histograms for single layer defined");
+  }
   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
-    h->Fill(tracklet->GetN2());
+    h->Fill(tracklet->GetN2(), fCentralityClass);
+    if(histosLayer && par){
+      if((hlayer = dynamic_cast<TProfile2D *>(histosLayer->At(itl)))){
+        GetEtaPhiAt(par, tracklet->GetX0(), eta, phi);
+        hlayer->Fill(eta, phi, tracklet->GetN2());
+      }
+    }
   }
   return h;
 }
@@ -636,11 +674,11 @@ TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
+  TH2 *h = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTrack)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -648,6 +686,7 @@ TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
   Int_t nclusters = 0;
   AliTRDseedV1 *tracklet = NULL;
   AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
+  if(!par) return NULL;
   Double_t momentumRec = par->P();
   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
@@ -668,9 +707,10 @@ TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
       }
       (*DebugStream()) << "NClustersTrack"
         << "Detector="  << detector
+        << "Centrality="<< fCentralityClass
         << "crossing="  << crossing
         << "momentumMC="<< momentumMC
-        << "momentumRec=" << momentumRec
+        << "momentumRec="<< momentumRec
         << "pdg="                              << pdg
         << "theta="                    << theta
         << "phi="                              << phi
@@ -680,7 +720,7 @@ TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
         << "\n";
     }
   }
-  h->Fill(nclusters);
+  h->Fill(nclusters, fCentralityClass);
   return h;
 }
 
@@ -692,16 +732,16 @@ TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL, *hSta = NULL; TH2 *hBarrel = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
+  TH2 *h = NULL, *hSta = NULL; TH3 *hBarrel = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsTrack)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
-  h->Fill(nTracklets);
+  h->Fill(nTracklets, fCentralityClass);
   if(!fkESD) return h;
   Int_t status = fkESD->GetStatus();
 
@@ -711,7 +751,7 @@ TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
   if((status & AliESDtrack::kTRDin) != 0){
     method = 1;
     // Full Barrel Track: Save momentum dependence
-    if(!(hBarrel = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsBAR)))){
+    if(!(hBarrel = dynamic_cast<TH3F *>(fContainer->At(kNtrackletsBAR)))){
       AliWarning("Method: Barrel.  Histogram not processed!");
       return NULL;
     }
@@ -720,45 +760,46 @@ TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
       AliError("Input track params missing");
       return NULL;
     }
-    hBarrel->Fill(par->P(), nTracklets);
+    p = par->P(); // p needed later in the debug streaming
+    hBarrel->Fill(p, nTracklets, fCentralityClass);
   } else {
     // Stand alone Track: momentum dependence not usefull
     method = 0;
-    if(!(hSta = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA)))) {
+    if(!(hSta = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsSTA)))) {
       AliWarning("Method: StandAlone.  Histogram not processed!");
       return NULL;
     }
-    hSta->Fill(nTracklets);
+    hSta->Fill(nTracklets, fCentralityClass);
   }
 
   if(DebugLevel() > 2){
     AliTRDseedV1 *tracklet = NULL;
-    Int_t sector = -1;
+    AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
+    Int_t sector = -1, stack = -1, detector;
     for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
       if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
-      sector = fGeo->GetSector(tracklet->GetDetector());
+      detector = tracklet->GetDetector();
+      sector = geo->GetSector(detector);
+      stack = geo->GetStack(detector);
       break;
     }
     (*DebugStream()) << "NTrackletsTrack"
       << "Sector="      << sector
+      << "Stack="       << stack
+      << "Centrality="  << fCentralityClass
       << "NTracklets="  << nTracklets
       << "Method="      << method
       << "p="           << p
       << "\n";
   }
   if(DebugLevel() > 3){
-    if(nTracklets == 1){
-      // If we have one Tracklet, check in which layer this happens
-      Int_t layer = -1;
-      AliTRDseedV1 *tracklet = NULL;
-      for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
-        if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
-      }
-      if(layer >= 0){
+    AliTRDseedV1 *tracklet = NULL;
+    for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+      if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
         (*DebugStream()) << "NTrackletsLayer"
-          << "Layer=" << layer
-          << "p=" << p
-          << "\n";
+        << "Layer=" << il
+        << "p=" << p
+        << "\n";
       }
     }
   }
@@ -773,11 +814,11 @@ TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
+  TH2 *h = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsCross)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -789,7 +830,7 @@ TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
 
     if(tracklet->IsRowCross()) ncross++;
   }
-  h->Fill(ncross);
+  h->Fill(ncross, fCentralityClass);
   return h;
 }
 
@@ -825,11 +866,11 @@ TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
  
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
+  TH2 *h = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsFindable)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -841,7 +882,7 @@ TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
   AliTRDpadPlane *pp;  
   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
     if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
-      tracklet->SetReconstructor(fReconstructor);
+      tracklet->SetReconstructor(AliTRDinfoGen::Reconstructor());
       nFound++;
     }
   }
@@ -885,11 +926,12 @@ TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
     memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
   }
+  AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
   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);
+    if((stack = geo->GetStack(z, il)) < 0) continue; // Not findable
+    pp = geo->GetPadPlane(il, stack);
     ymin = pp->GetCol0() + epsilon;
     ymax = pp->GetColEnd() - epsilon; 
     zmin = pp->GetRowEnd() + epsilon; 
@@ -910,7 +952,7 @@ TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
       }
   }
   
-  h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
+  h->Fill((nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1), fCentralityClass);
   AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
   return h;
 }
@@ -923,18 +965,18 @@ TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
+  TH3 *h = NULL;
+  if(!(h = dynamic_cast<TH3S*>(fContainer->At(kChi2)))) {
     AliWarning("No Histogram defined.");
     return NULL;
   }
   Int_t n = fkTrack->GetNumberOfTracklets();
   if(!n) return NULL;
 
-  h->Fill(n, fkTrack->GetChi2()/n);
+  h->Fill(n, fkTrack->GetChi2()/n, fCentralityClass);
   return h;
 }
 
@@ -946,11 +988,11 @@ TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TProfile *h = NULL;
-  if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
+  TH3F *h = NULL;
+  if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -965,31 +1007,36 @@ TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
       Int_t localtime        = c->GetLocalTimeBin();
       Double_t absoluteCharge = TMath::Abs(c->GetQ());
-      h->Fill(localtime, absoluteCharge);
+      h->Fill(localtime, absoluteCharge, fCentralityClass);
       if(DebugLevel() > 3){
+        Int_t inChamber = c->IsInChamber() ? 1 : 0;
         Double_t distance[2];
         GetDistanceToTracklet(distance, tracklet, c);
         Float_t theta = TMath::ATan(tracklet->GetZref(1));
         Float_t phi = TMath::ATan(tracklet->GetYref(1));
-        Float_t momentum = 0.;
+        AliExternalTrackParam *trdPar = fkTrack->GetTrackIn();
+        Float_t momentumMC = 0, momentumRec = trdPar ? trdPar->P() : fkTrack->P(); // prefer Track Low
         Int_t pdg = 0;
         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
-        UShort_t TPCncls = fkESD ? fkESD->GetTPCncls() : 0;
+        UShort_t tpcNCLS = fkESD ? fkESD->GetTPCncls() : 0;
         if(fkMC){
-          if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
+          if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
           pdg = fkMC->GetPDG();
         }
         (*DebugStream()) << "PHt"
           << "Detector="       << detector
+          << "Centrality="<< fCentralityClass
           << "crossing="       << crossing
+          << "inChamber=" << inChamber
           << "Timebin="                << localtime
           << "Charge="         << absoluteCharge
-          << "momentum="       << momentum
+          << "momentumMC="     << momentumMC
+          << "momentumRec="    << momentumRec
           << "pdg="                            << pdg
           << "theta="                  << theta
           << "phi="                            << phi
           << "kinkIndex="      << kinkIndex
-          << "TPCncls="                << TPCncls
+          << "TPCncls="                << tpcNCLS
           << "dy="        << distance[0]
           << "dz="        << distance[1]
           << "c.="        << c
@@ -1008,34 +1055,35 @@ TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
-    return NULL;
+     AliDebug(4, "No Track defined.");
+     return NULL;
   }
-  TProfile *h = NULL;
-  if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
+  TH3 *h = NULL;
+  if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
-  AliTRDseedV1 *tracklet = NULL;
-  AliTRDcluster *c = NULL;
+  AliTRDseedV1 *tracklet(NULL);
+  AliTRDcluster *c(NULL);
   Double_t xd(0.), dqdl(0.);
   TVectorD vq(AliTRDseedV1::kNtb), vxd(AliTRDseedV1::kNtb), vdqdl(AliTRDseedV1::kNtb);
   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
     Int_t det(tracklet->GetDetector());
     Bool_t rc(tracklet->IsRowCross());
-    Int_t jc(0);
-    for(Int_t ic(0); ic<AliTRDseedV1::kNclusters; ic++){
-      if(!(c = tracklet->GetClusters(ic))) continue;
+    for(Int_t ic(0); ic<AliTRDseedV1::kNtb; ic++){
+      Bool_t kFIRST(kFALSE);
+      if(!(c = tracklet->GetClusters(ic))){
+         if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
+      } else kFIRST=kTRUE;
       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
-      jc = (ic-AliTRDseedV1::kNtb);
-      if(jc<0) jc+=AliTRDseedV1::kNtb;
-      xd = tracklet->GetX0() - c->GetX(); vxd[jc] = xd;
-      dqdl=tracklet->GetdQdl(ic); vdqdl[jc] = dqdl;
-      h->Fill(xd, dqdl);
+      xd = tracklet->GetX0() - c->GetX(); vxd[ic] = xd;
+      dqdl=tracklet->GetdQdl(ic); vdqdl[ic] = dqdl;
+      vq[ic]=c->GetQ();
+      if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) vq[ic]+=c->GetQ();
+      h->Fill(xd, dqdl, fCentralityClass);
     }
     if(DebugLevel() > 3){
-      vq[jc]+=c->GetQ();
       (*DebugStream()) << "PHx"
         << "det="  << det
         << "rc="   << rc
@@ -1055,11 +1103,11 @@ TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
+  TH2 *h = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeCluster)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -1067,9 +1115,14 @@ TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
   AliTRDcluster *c = NULL;
   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
-    for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
-      if(!(c = tracklet->GetClusters(itime))) continue;
-      h->Fill(c->GetQ());
+    for(Int_t ic(0); ic < AliTRDseedV1::kNtb; ic++){
+      Bool_t kFIRST(kFALSE);
+      if(!(c = tracklet->GetClusters(ic))) {
+        if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
+      } else kFIRST = kTRUE;
+      Float_t q(c->GetQ());
+      if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) q+=c->GetQ();
+      h->Fill(q, fCentralityClass);
     }
   }
   return h;
@@ -1082,11 +1135,11 @@ TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
-  TH1 *h = NULL;
-  if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
+  TH2 *h = NULL;
+  if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeTracklet)))){
     AliWarning("No Histogram defined.");
     return NULL;
   }
@@ -1101,7 +1154,7 @@ TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
       if(!(c = tracklet->GetClusters(ic))) continue;
       qTot += TMath::Abs(c->GetQ());
     }
-    h->Fill(qTot);
+    h->Fill(qTot, fCentralityClass);
     if(DebugLevel() > 3){
       Int_t crossing = (Int_t)tracklet->IsRowCross();
       Int_t detector = tracklet->GetDetector();
@@ -1117,6 +1170,7 @@ TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
       }
       (*DebugStream()) << "ChargeTracklet"
         << "Detector="  << detector
+        << "Centrality="<< fCentralityClass
         << "crossing="  << crossing
         << "momentum=" << momentum
         << "nTracklets="<< nTracklets
@@ -1139,7 +1193,7 @@ TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
   //
   if(track) fkTrack = track;
   if(!fkTrack){
-    AliWarning("No Track defined.");
+    AliDebug(4, "No Track defined.");
     return NULL;
   }
   TH1 *h = NULL;
@@ -1164,13 +1218,6 @@ TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
 
 
 //________________________________________________________
-void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
-{
-
-  fReconstructor->SetRecoParam(r);
-}
-
-//________________________________________________________
 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
 {
   Float_t x = c->GetX();
@@ -1178,17 +1225,34 @@ void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const
   dist[1] = c->GetZ() - tracklet->GetZat(x);
 }
 
+//________________________________________________________
+void AliTRDcheckDET::GetEtaPhiAt(const AliExternalTrackParam *track, Double_t x, Double_t &eta, Double_t &phi){
+  //
+  // Get phi and eta at a given radial position
+  // 
+  AliExternalTrackParam workpar(*track);
+
+  Double_t posLocal[3];
+  Bool_t sucPos = workpar.GetXYZAt(x, fEvent->GetRunInfo()->GetMagneticField(), posLocal);
+  Double_t sagPhi = sucPos ? TMath::ATan2(posLocal[1], posLocal[0]) : 0.;
+  phi = sagPhi;
+  eta = workpar.Eta();
+}
+
 
 //_______________________________________________________
-TH1* AliTRDcheckDET::MakePlotChi2()
+TH1* AliTRDcheckDET::MakePlotChi2() const
 {
 // Plot chi2/track normalized to number of degree of freedom 
 // (tracklets) and compare with the theoretical distribution.
 // 
 // Alex Bercuci <A.Bercuci@gsi.de>
 
-  TH2S *h2 = (TH2S*)fContainer->At(kChi2);
+  return NULL;
+
+/*  TH2S *h2 = (TH2S*)fContainer->At(kChi2);
   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
+  f.SetParLimits(1,1, 1e100);
   TLegend *leg = new TLegend(.7,.7,.95,.95);
   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
   TH1D *h1 = NULL;
@@ -1209,56 +1273,254 @@ TH1* AliTRDcheckDET::MakePlotChi2()
   }
   leg->Draw();
   gPad->SetLogy();
-  return h1;
+  return h1;*/
 }
 
 
 //________________________________________________________
 TH1* AliTRDcheckDET::MakePlotNTracklets(){
-  //
-  // Make nice bar plot of the number of tracklets in each method
-  //
-  TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
-  TH1D *hBAR = tmp->ProjectionY();
-  TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
-  TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
-  TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
-  leg->SetBorderSize(1);
-  leg->SetFillColor(0);
-
-  Float_t scale = hCON->Integral();
-  hCON->Scale(100./scale);
-  hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
-  hCON->SetBarWidth(0.2);
-  hCON->SetBarOffset(0.6);
-  hCON->SetStats(kFALSE);
-  hCON->GetYaxis()->SetRangeUser(0.,40.);
-  hCON->GetYaxis()->SetTitleOffset(1.2);
-  hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
-  hCON->SetMaximum(55.);
-
-  hBAR->Scale(100./scale);
-  hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
-  hBAR->SetBarWidth(0.2);
-  hBAR->SetBarOffset(0.2);
-  hBAR->SetTitle("");
-  hBAR->SetStats(kFALSE);
-  hBAR->GetYaxis()->SetRangeUser(0.,40.);
-  hBAR->GetYaxis()->SetTitleOffset(1.2);
-  hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
-
-  hSTA->Scale(100./scale);
-  hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
-  hSTA->SetBarWidth(0.2);
-  hSTA->SetBarOffset(0.4);
-  hSTA->SetTitle("");
-  hSTA->SetStats(kFALSE);
-  hSTA->GetYaxis()->SetRangeUser(0.,40.);
-  hSTA->GetYaxis()->SetTitleOffset(1.2);
-  hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
-  leg->Draw();
-  gPad->Update();
-  return hCON;
+       //
+       // Make nice bar plot of the number of tracklets in each method
+       //
+       TH3 *hTracklets3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
+       if(!hTracklets3D){
+               AliError("Tracklet Histogram not found");
+               return NULL;
+       }
+       TH1 *hBAR = hTracklets3D->Project3D("y");
+       hBAR->SetName("hBAR");
+       hBAR->SetTitle("Number of Tracklets");
+       hBAR->Scale(1./hBAR->Integral());
+       hBAR->SetLineColor(kBlack);
+       hBAR->SetLineWidth(2);
+       hBAR->Draw();
+
+       // Draw also centrality-dependent plots
+       TH1 *hBarCent[AliTRDeventInfo::kCentralityClasses];
+       Int_t nHistsCentrality = 0;
+       for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+               hTracklets3D->GetZaxis()->SetRange(icent+2, icent+2);
+               hBarCent[icent] = hTracklets3D->Project3D("y");
+               if(!(hBarCent[icent] && hBarCent[icent]->GetEntries())){
+                       delete hBarCent[icent];
+                       hBarCent[icent] = NULL;
+                       continue;
+               }
+               hBarCent[icent]->SetName(Form("hBarCent_%d", icent));
+               hBarCent[icent]->SetTitle("Number of Tracklets");
+               hBarCent[icent]->Scale(1./hBarCent[icent]->Integral());
+               hBarCent[icent]->SetLineColor(fkColorsCentrality[icent]);
+               hBarCent[icent]->Draw("same");
+               nHistsCentrality++;
+       }
+       hTracklets3D->GetZaxis()->SetRange(0, hTracklets3D->GetNbinsZ());
+       AliInfo(Form("Number of Centrality Classes: %d", nHistsCentrality));
+       if(nHistsCentrality){
+               TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
+               leg->SetBorderSize(0);
+               leg->SetFillStyle(0);
+               for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++)
+                       if(hBarCent[icent]) leg->AddEntry(hBarCent[icent], Form("Centrality Class %d", icent), "l");
+               leg->Draw();
+               gPad->Update();
+       }
+       return hBAR;
+
+       /*
+         TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
+         TH1D *hBAR = tmp->ProjectionY();
+         TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
+         TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
+         TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
+         leg->SetBorderSize(1);
+         leg->SetFillColor(0);
+
+         Float_t scale = hCON->Integral();
+         if(scale) hCON->Scale(100./scale);
+         hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
+         hCON->SetBarWidth(0.2);
+         hCON->SetBarOffset(0.6);
+         hCON->SetStats(kFALSE);
+         hCON->GetYaxis()->SetRangeUser(0.,40.);
+         hCON->GetYaxis()->SetTitleOffset(1.2);
+         hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
+         hCON->SetMaximum(55.);
+
+         if(scale) hBAR->Scale(100./scale);
+         hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
+         hBAR->SetBarWidth(0.2);
+         hBAR->SetBarOffset(0.2);
+         hBAR->SetTitle("");
+         hBAR->SetStats(kFALSE);
+         hBAR->GetYaxis()->SetRangeUser(0.,40.);
+         hBAR->GetYaxis()->SetTitleOffset(1.2);
+         hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
+
+         if(scale) hSTA->Scale(100./scale);
+         hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
+         hSTA->SetBarWidth(0.2);
+         hSTA->SetBarOffset(0.4);
+         hSTA->SetTitle("");
+         hSTA->SetStats(kFALSE);
+         hSTA->GetYaxis()->SetRangeUser(0.,40.);
+         hSTA->GetYaxis()->SetTitleOffset(1.2);
+         hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
+         leg->Draw();
+         gPad->Update();
+         return hCON;
+        */
+}
+
+//________________________________________________________
+void AliTRDcheckDET::MakePlotNclustersTrack(){
+       //
+       // Plot number of clusters
+       // Put histos from all centrality classes into one pad
+       //
+       TH2 *hClusters = dynamic_cast<TH2 *>(fContainer->FindObject("hNcls"));
+       if(!hClusters){
+               AliError("Cluster histogram not found in the output");
+       }
+       TH1 *hAllCentrality = hClusters->ProjectionX("hNcls_all");
+       hAllCentrality->SetTitle("Number of clusters/track");
+       hAllCentrality->Scale(1./hAllCentrality->Integral());
+       hAllCentrality->SetLineColor(kBlack);
+       hAllCentrality->SetLineWidth(2);
+       hAllCentrality->GetYaxis()->SetRangeUser(0, 0.02);
+       hAllCentrality->Draw();
+       PutTrendValue("NClustersTrack", hAllCentrality->GetMean());
+       PutTrendValue("NClustersTrackRMS", hAllCentrality->GetRMS());
+
+       // Now look at single centrality classes
+       TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
+       Int_t nHistsCentrality = 0;
+       for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+               hProjCentral[icent] = hClusters->ProjectionX(Form("hNcls_%d", icent), icent+1, icent+1);
+               if(!hProjCentral[icent]->GetEntries()){
+                       delete hProjCentral[icent];
+                       hProjCentral[icent] = NULL;
+                       continue;
+               }
+               hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
+               hProjCentral[icent]->SetTitle("Number of clusters/track");
+               hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
+               hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.03);
+               hProjCentral[icent]->Draw("same");
+               nHistsCentrality++;
+       }
+       AliInfo(Form("%d centrality classes found", nHistsCentrality));
+       if(nHistsCentrality){
+               // Draw nice legend
+               TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
+               leg->SetBorderSize(0);
+               leg->SetFillStyle(0);
+               for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+                       if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
+               }
+               leg->Draw();
+               gPad->Update();
+       }
+}
+
+//________________________________________________________
+void AliTRDcheckDET::MakePlotNclustersTracklet(){
+       //
+       // Plot number of clusters for different centrality classes
+       //
+       TH2 *hClusters = dynamic_cast<TH2*>(fContainer->FindObject("hNclTls"));
+       if(!hClusters){
+               AliError("Histogram for the number of clusters per tracklet not available");
+               return;
+       }
+       TH1 *hAllCentrality = hClusters->ProjectionX("hNclsTls_all");
+       hAllCentrality->SetTitle("Number of clusters/track");
+       hAllCentrality->Scale(1./hAllCentrality->Integral());
+       hAllCentrality->SetLineColor(kBlack);
+       hAllCentrality->SetLineWidth(2);
+       hAllCentrality->GetYaxis()->SetRangeUser(0, 0.3);
+       hAllCentrality->Draw("pc");
+       PutTrendValue("NClustersTracklet", hAllCentrality->GetMean());
+       PutTrendValue("NClustersTrackletRMS", hAllCentrality->GetRMS());
+
+       // Now look at single centrality classes
+       TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
+       Int_t nHistsCentrality = 0;
+       for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+               hProjCentral[icent] = hClusters->ProjectionX(Form("hNclsTls_%d", icent), icent+1, icent+1);
+               if(!hProjCentral[icent]->GetEntries()){
+                       delete hProjCentral[icent];
+                       hProjCentral[icent] = NULL;
+                       continue;
+               }
+               hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
+               hProjCentral[icent]->SetTitle("Number of clusters/track");
+               hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
+               hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.3);
+               hProjCentral[icent]->Draw("pcsame");
+               nHistsCentrality++;
+       }
+       AliInfo(Form("%d centrality classes found", nHistsCentrality));
+       if(nHistsCentrality){
+               // Draw nice legend
+               TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
+               leg->SetBorderSize(0);
+               leg->SetFillStyle(0);
+               for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+                       if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
+               }
+               leg->Draw();
+               gPad->Update();
+       }
+}
+
+//________________________________________________________
+void AliTRDcheckDET::MakePlotTrackletCharge(){
+       //
+       // Draw tracklet charge for different centrality classes
+       //
+       TH2 *hQt = dynamic_cast<TH2*>(fContainer->FindObject("hQtrklt"));
+       if(!hQt){
+               AliError("Histogram for tracklet charges not found");
+               return;
+       }
+       // First Project all charhes
+       TH1 *hQtAll = hQt->ProjectionX("hQtAll");
+       hQtAll->SetTitle("Tracklet Charge");
+       Double_t scalefactor = 0.7 * hQtAll->Integral() / hQtAll->GetMaximum();
+       hQtAll->Scale(scalefactor/hQtAll->Integral());
+       hQtAll->GetYaxis()->SetRangeUser(0., 1.);
+       hQtAll->SetLineColor(kBlack);
+       hQtAll->SetLineWidth(2);
+       hQtAll->Draw("c");
+       PutTrendValue("ChargeTracklet", hQtAll->GetMaximumBin());
+       PutTrendValue("ChargeTrackletRMS", hQtAll->GetRMS());
+
+       TH1 *hQtCent[AliTRDeventInfo::kCentralityClasses];
+       Int_t nHistsCentrality = 0;
+       for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+               hQtCent[icent] = hQt->ProjectionX(Form("hQt_%d", icent), icent+1, icent+1);
+               if(!hQtCent[icent]->GetEntries()){
+                       delete hQtCent[icent];
+                       continue;
+               }
+               hQtCent[icent]->SetTitle("Tracklet Charge");
+               hQtCent[icent]->Scale(scalefactor/hQtCent[icent]->Integral());
+               hQtCent[icent]->GetYaxis()->SetRangeUser(0., 1.);
+               hQtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
+               hQtCent[icent]->Draw("csame");
+               nHistsCentrality++;
+       }
+       AliInfo(Form("%d centrality classes found", nHistsCentrality));
+       if(nHistsCentrality){
+               TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
+               leg->SetBorderSize(0);
+               leg->SetFillStyle(0);
+               for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+                       if(hQtCent[icent]) leg->AddEntry(hQtCent[icent], Form("Centrality Class %d", icent), "l");
+               }
+               leg->Draw();
+               gPad->Update();
+       }
 }
 
 //________________________________________________________
@@ -1267,9 +1529,6 @@ void AliTRDcheckDET::MakePlotnTrackletsVsP(){
   // Plot abundance of tracks with number of tracklets as function of momentum
   //
 
-
-
-
   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
@@ -1279,7 +1538,12 @@ void AliTRDcheckDET::MakePlotnTrackletsVsP(){
     g[itl]->SetMarkerStyle(20 + itl);
   }
 
-  TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
+  TH3 *hBar3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
+  if(!hBar3D){
+         AliError("Histogram for the number of tracklets vs p not available");
+         return;
+  }
+  TH2 *hBar = (TH2 *)hBar3D->Project3D("yx");
   TAxis *ax(hBar->GetXaxis());
   Int_t np(ax->GetNbins());
   for(Int_t ipBin(1); ipBin<np; ipBin++){
@@ -1320,53 +1584,183 @@ void AliTRDcheckDET::MakePlotnTrackletsVsP(){
 }
 
 //________________________________________________________
-TH1* AliTRDcheckDET::MakePlotPulseHeight(){
+Bool_t AliTRDcheckDET::MakePlotPulseHeight(){
+       //
+       // Create Plot of the Pluse Height Spectrum
+       //
+       TCanvas *output = gPad->GetCanvas();
+       output->Divide(2);
+       output->cd(1);
+
+       TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
+       //TH3F *hPhx = dynamic_cast<TH3F *>(arr->At(0)),
+       TH3F *hPht = dynamic_cast<TH3F *>(arr->At(1));
+       if(!hPht) return kFALSE;
+       // Project centrality of the 2 histograms
+       //TH2 *hProjCentX = dynamic_cast<TH2 *>(hPhx->Project3D("yx")),
+       TH2 *hProjCentT = dynamic_cast<TH2 *>(hPht->Project3D("yx"));
+  if(!hProjCentT) return kFALSE;
+       //hProjCentX->SetName("hProjCentX");
+       hProjCentT->SetName("hProjCentT");
+       // Draw 2D histogram versus time on pad 2
+       output->cd(2);
+       hProjCentT->SetTitle("2D Pulse Height vs. Time");
+       hProjCentT->GetYaxis()->SetTitleOffset(1.8);
+  hProjCentT->GetYaxis()->SetRangeUser(0., 600.);
+       hProjCentT->SetStats(kFALSE);
+       hProjCentT->Draw("colz");
+
+  // Fit each slice of the pulse height spectrum with a Landau function
+  TGraphErrors *landaufit = new TGraphErrors;
+  landaufit->SetTitle();
+  landaufit->SetMarkerColor(kWhite);
+  landaufit->SetLineColor(kWhite);
+  landaufit->SetMarkerStyle(25);
+  TH1 *projection;
+  Int_t ntime = 0;
+  for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetLast(); it++){
+    projection = hProjCentT->ProjectionY("proj", it, it);
+    if(!projection->GetEntries()){
+      delete projection;
+      continue;  
+    }
+    TF1 fitfun("fitfun", "landau", 0, 1000);
+    fitfun.SetParLimits(0, 1e-1, 1e9); fitfun.SetParameter(0, 1000);
+    fitfun.SetParLimits(1, 10, 200); fitfun.SetParameter(1, 50);
+    fitfun.SetParLimits(2, 1e-1, 200); fitfun.SetParameter(2, 10);
+    projection->Fit(&fitfun, "QN", "", 0, 1000); 
+    landaufit->SetPoint(ntime, hProjCentT->GetXaxis()->GetBinCenter(it), fitfun.GetParameter(1));
+    landaufit->SetPointError(ntime, hProjCentT->GetXaxis()->GetBinWidth(it)/2., fitfun.GetParError(1));
+    ntime++;
+    delete projection;
+  }
+  landaufit->Draw("lpesame");
+
+       // Fill 1D PHS as function of time resp. radius (same binning of the 2 histograms)
+       //TH1 *hPhsX = new TH1F("hPhsX", "Average PH vs X", hProjCentT->GetXaxis()->GetNbins(), hProjCentT->GetXaxis()->GetXbins()->GetArray()),
+  Int_t nbinsT = hProjCentT->GetXaxis()->GetNbins();
+       TH1     *hPhsT = new TH1F("hPhsT", "Average Ph vs Time; Time (100 ns); Average Pulse Height (a.u.)", nbinsT, -0.5, nbinsT - 0.5),
+               *htmp = NULL;
+       /*for(Int_t irad = 1; irad <= hProjCentX->GetXaxis()->GetNbins(); irad++){
+               htmp = hProjCentX->ProjectionY("tmp", irad, irad);
+               hPhsX->SetBinContent(irad, htmp->GetMean());
+               hPhsX->SetBinError(irad, htmp->GetMeanError());
+               delete htmp;
+       }*/
+       for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetNbins(); it++){
+               htmp = hProjCentT->ProjectionY("tmp", it, it);
+               hPhsT->SetBinContent(it, htmp->GetMean());
+               hPhsT->SetBinError(it, htmp->GetMeanError());
+               delete htmp;
+       }
+       output->cd(1);
+       // Draw 1D histograms
+       if(hPhsT->GetEntries()){
+               hPhsT->SetMarkerStyle(24);
+               hPhsT->SetMarkerColor(kBlack);
+               hPhsT->SetLineColor(kBlack);
+               hPhsT->GetYaxis()->SetTitleOffset(1.5);
+               hPhsT->Draw("e1");
+               // Now fit the PHS with respect to time to get plateau and slope
+               // Trending for the pulse height: plateau value, slope and timebin of the maximum
+               TLinearFitter fit(1,"pol1");
+               Double_t time = 0.;
+               for(Int_t itime = 10; itime <= 20; itime++){
+                       time = Double_t(itime);
+                       Double_t err(hPhsT->GetBinError(itime + 1));
+                       if(err>1.e-10) fit.AddPoint(&time, hPhsT->GetBinContent(itime + 1), err);
+               }
+               if(!fit.Eval()){
+                       Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
+                       Double_t slope = fit.GetParameter(1);
+                       PutTrendValue("PHplateau", plateau);
+                       PutTrendValue("PHslope", slope);
+                       PutTrendValue("PHamplificationPeak", static_cast<Double_t>(hPhsT->GetMaximumBin()-1));
+                       AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(hPhsT->GetMaximumBin()-1)));
+               }
+       }
+       /*if(hPhsX->GetEntries()){
+               hPhsX->SetMarkerStyle(22);
+               hPhsX->SetMarkerColor(kBlue);
+               hPhsX->SetLineColor(kBlue);
+               hPhsX->GetYaxis()->SetTitleOffset(1.5);
+               hPhsX->Draw("e1same");
+               // create axis according to the histogram dimensions of the original second histogram
+               TGaxis *axis = new TGaxis(gPad->GetUxmin(),
+                               gPad->GetUymax(),
+                               gPad->GetUxmax(),
+                               gPad->GetUymax(),
+                               -0.08, 4.88, 510,"-L");
+               axis->SetLineColor(kBlue);
+               axis->SetLabelColor(kBlue);
+               axis->SetTextColor(kBlue);
+               axis->SetTitle("x_{0}-x_{c} [cm]");
+               axis->Draw();
+               gPad->Update();
+       }*/
+
+       // Centrality-dependent Pulse-Height Spectrum
+       TH1 *hPhtCent[AliTRDeventInfo::kCentralityClasses];
+       TH2 *hPtmp;
+       Int_t nHistsCentrality = 0;
+       for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+               hPht->GetZaxis()->SetRange(icent+2,icent+2);
+               hPtmp = dynamic_cast<TH2*>(hPht->Project3D("yx"));
+               if(!(hPtmp && hPtmp->GetEntries())){
+                       hPhtCent[icent] = NULL;
+                       continue;
+               }
+               hPhtCent[icent] = new TH1F(Form("hPhtCent_%d", icent), "Average Ph vs Time", hPtmp->GetNbinsX(), hPtmp->GetXaxis()->GetXbins()->GetArray());
+               for(Int_t it = 1; it <= hPtmp->GetNbinsX(); it++){
+                       htmp = hPtmp->ProjectionY("htmp", it, it);
+                       hPhtCent[icent]->SetBinContent(it, htmp->GetMean());
+                       hPhtCent[icent]->SetBinError(it, htmp->GetMeanError());
+                       delete htmp;
+               }
+               delete hPtmp;
+               hPhtCent[icent]->SetMarkerStyle(24);
+               hPhtCent[icent]->SetMarkerColor(fkColorsCentrality[icent]);
+               hPhtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
+               hPhtCent[icent]->Draw("e1same");
+               nHistsCentrality++;
+       }
+       hPht->GetZaxis()->SetRange(0, hPht->GetNbinsZ());
+       if(nHistsCentrality){
+               TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
+               leg->SetBorderSize(0);
+               leg->SetFillStyle(0);
+               for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
+                       if(hPhtCent[icent]) leg->AddEntry(hPhtCent[icent], Form("Centrality Class %d", icent), "p");
+               }
+               leg->Draw();
+               gPad->Update();
+       }
+
+       // delete 2D Projection of the PHS vs x since it is only used to calculate the 1D projection
+       //delete hProjCentX;
+
+       return kTRUE;
+}
+
+//________________________________________________________
+void AliTRDcheckDET::MakePlotMeanClustersLayer(){
   //
-  // Create Plot of the Pluse Height Spectrum
+  // Create Summary plot for the mean number of clusters per layer
   //
-  TH1 *h, *h1, *h2;
-  TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
-  h = (TH1F*)arr->At(0);
-  h->SetMarkerStyle(24);
-  h->SetMarkerColor(kBlack);
-  h->SetLineColor(kBlack);
-  h->Draw("e1");
-  // Trending for the pulse height: plateau value, slope and timebin of the maximum
-  TLinearFitter fit(1,"pol1");
-  Double_t time = 0.;
-  for(Int_t itime = 10; itime <= 20; itime++){
-    time = static_cast<Double_t>(itime);
-    fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
+  TCanvas *output = gPad->GetCanvas();
+  output->Divide(3,2);
+  TObjArray *histos = (TObjArray *)fContainer->At(kNclustersLayer);
+  if(!histos){
+    AliWarning("Histos for each layer not found");
+    return;
+  }
+  TProfile2D *hlayer = NULL;
+  for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
+    if(!(hlayer = dynamic_cast<TProfile2D *>(histos->At(ily)))) continue;
+    output->cd(ily + 1);
+    gPad->SetGrid(0,0);
+    hlayer->Draw("colz");
   }
-  fit.Eval();
-  Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
-  Double_t slope = fit.GetParameter(1);
-  PutTrendValue("PHplateau", plateau);
-  PutTrendValue("PHslope", slope);
-  PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
-  AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
-//   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
-  h1 = (TH1F *)arr->At(1);
-  h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
-  for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++) 
-    h2->SetBinContent(ibin, h1->GetBinContent(ibin));
-  h2->SetMarkerStyle(22);
-  h2->SetMarkerColor(kBlue);
-  h2->SetLineColor(kBlue);
-  h2->Draw("e1same");
-  gPad->Update();
-//   create axis according to the histogram dimensions of the original second histogram
-  TGaxis *axis = new TGaxis(gPad->GetUxmin(),
-                    gPad->GetUymax(),
-                    gPad->GetUxmax(),
-                    gPad->GetUymax(),
-                    -0.08, 4.88, 510,"-L");
-  axis->SetLineColor(kBlue);
-  axis->SetLabelColor(kBlue);
-  axis->SetTextColor(kBlue);
-  axis->SetTitle("x_{0}-x_{c} [cm]");
-  axis->Draw();
-  return h1;
 }
 
 //________________________________________________________