PWG1/TRD/AliTRDcheckDET.cxx

   1 ////////////////////////////////////////////////////////////////////////////
   2 //                                                                        //
   3 //                                                                        //
   4 //  Basic checks for tracking and detector performance                    //
   5 //
   6 //     For the moment (15.10.2009) the following checks are implemented    //
   7 //       - Number of clusters/track
   8 //       - Number of clusters/tracklet
   9 //       - Number of tracklets/track from different sources (Barrel, StandAlone)
  10 //       - Number of findable tracklets
  11 //       - Number of tracks per event and TRD sector
  12 //       - <PH>
  13 //       - Chi2 distribution for tracks
  14 //       - Charge distribution per cluster and tracklet
  15 //       - Number of events and tracks per trigger source
  16 //       - Trigger purity
  17 //       - Track and Tracklet propagation status
  18 //
  19 //  Authors:                                                              //
  20 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  21 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
  22 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  23 //                                                                        //
  24 ////////////////////////////////////////////////////////////////////////////
  25
  26 #include <TAxis.h>
  27 #include <TCanvas.h>
  28 #include <TFile.h>
  29 #include <TH1F.h>
  30 #include <TH1I.h>
  31 #include <TH2F.h>
  32 #include <TH3S.h>
  33 #include <TH3F.h>
  34 #include <TF1.h>
  35 #include <TGaxis.h>
  36 #include <TGraph.h>
  37 #include <TGraphErrors.h>
  38 #include <TLegend.h>
  39 #include <TLinearFitter.h>
  40 #include <TMath.h>
  41 #include <TMap.h>
  42 #include <TProfile2D.h>
  43 #include <TObjArray.h>
  44 #include <TObject.h>
  45 #include <TObjString.h>
  46
  47 #include <TPad.h>
  48 #include <TProfile.h>
  49 #include <TProfile2D.h>
  50 #include <TROOT.h>
  51 #include <TChain.h>
  52
  53 #include "AliLog.h"
  54 #include "AliTRDcluster.h"
  55 #include "AliESDHeader.h"
  56 #include "AliESDRun.h"
  57 #include "AliESDtrack.h"
  58 #include "AliExternalTrackParam.h"
  59 #include "AliTRDgeometry.h"
  60 #include "AliTRDpadPlane.h"
  61 #include "AliTRDSimParam.h"
  62 #include "AliTRDseedV1.h"
  63 #include "AliTRDtrackV1.h"
  64 #include "AliTRDtrackerV1.h"
  65 #include "AliTRDReconstructor.h"
  66 #include "AliTrackReference.h"
  67 #include "AliTrackPointArray.h"
  68 #include "AliTracker.h"
  69 #include "TTreeStream.h"
  70
  71 #include "info/AliTRDtrackInfo.h"
  72 #include "info/AliTRDeventInfo.h"
  73 #include "AliTRDinfoGen.h"
  74 #include "AliTRDcheckDET.h"
  75
  76 #include <cstdio>
  77 #include <iostream>
  78
  79 ClassImp(AliTRDcheckDET)
  80
  81 //_______________________________________________________
  82 AliTRDcheckDET::AliTRDcheckDET():
  83   AliTRDrecoTask()
  84   ,fCentralityClass(-1)
  85   ,fTriggerNames(NULL)
  86   ,fFlags(0)
  87 {
  88   //
  89   // Default constructor
  90   //
  91   SetNameTitle("TRDcheckDET", "Basic TRD data checker");
  92 }
  93
  94 //_______________________________________________________
  95 AliTRDcheckDET::AliTRDcheckDET(char* name):
  96   AliTRDrecoTask(name, "Basic TRD data checker")
  97   ,fCentralityClass(-1)
  98   ,fTriggerNames(NULL)
  99   ,fFlags(0)
 100 {
 101   //
 102   // Default constructor
 103   //
 104   InitFunctorList();
 105 }
 106
 107
 108 //_______________________________________________________
 109 AliTRDcheckDET::~AliTRDcheckDET(){
 110   //
 111   // Destructor
 112   //
 113   if(fTriggerNames) delete fTriggerNames;
 114 }
 115
 116 //_______________________________________________________
 117 void AliTRDcheckDET::UserCreateOutputObjects(){
 118   //
 119   // Create Output Objects
 120   //
 121   AliTRDrecoTask::UserCreateOutputObjects();
 122   if(!fTriggerNames) fTriggerNames = new TMap();
 123 }
 124
 125 //_______________________________________________________
 126 void AliTRDcheckDET::UserExec(Option_t *opt){
 127   //
 128   // Execution function
 129   // Filling TRD quality histos
 130   //
 131   AliDebug(2, Form("EventInfo[%p] Header[%p]", (void*)fEvent, (void*)(fEvent?fEvent->GetEventHeader():NULL)));
 132   if(fEvent) fCentralityClass = fEvent->GetCentrality();
 133   else fCentralityClass = -1;  // Assume pp
 134
 135   AliTRDrecoTask::UserExec(opt);
 136
 137   TH1F *histo(NULL); AliTRDtrackInfo *fTrackInfo(NULL); Int_t nTracks(0);               // Count the number of tracks per event
 138   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 139     if(!fTracks->UncheckedAt(iti)) continue;
 140     if(!(fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti)))) continue;
 141     if(!fTrackInfo->GetTrack()) continue;
 142     nTracks++;
 143   }
 144   if(nTracks)
 145     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent)))) histo->Fill(nTracks);
 146
 147   if(!fEvent->GetEventHeader()) return; // For trigger statistics event header is essential
 148   Int_t triggermask = fEvent->GetEventHeader()->GetTriggerMask();
 149   TString triggername =  fEvent->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 150   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
 151   if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger)))) histo->Fill(triggermask);
 152
 153   if(nTracks){
 154     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks)))) histo->Fill(triggermask);
 155   }
 156   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 157     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 158     // also set the label for both histograms
 159     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))))
 160       histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 161     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))))
 162       histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 163   }
 164 }
 165
 166
 167 //_______________________________________________________
 168 Bool_t AliTRDcheckDET::PostProcess(){
 169   //
 170   // Do Postprocessing (for the moment set the number of Reference histograms)
 171   //
 172
 173   TH1F *h(NULL), *h1(NULL);
 174
 175   // Calculate of the trigger clusters purity
 176   if((h  = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"))) &&
 177      (h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks")))) {
 178     h1->Divide(h);
 179     Float_t purities[20], val = 0; memset(purities, 0, 20*sizeof(Float_t));
 180     TString triggernames[20];
 181     Int_t nTriggerClasses = 0;
 182     for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
 183       if((val = h1->GetBinContent(ibin))){
 184         purities[nTriggerClasses] = val;
 185         triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
 186         nTriggerClasses++;
 187       }
 188     }
 189
 190     if((h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity)))){
 191       TAxis *ax = h->GetXaxis();
 192       for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
 193         h->Fill(itrg, purities[itrg]);
 194         ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
 195       }
 196       ax->SetRangeUser(-0.5, nTriggerClasses+.5);
 197       h->GetYaxis()->SetRangeUser(0,1);
 198     }
 199   }
 200
 201   // track status
 202   if((h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus)))){
 203     Float_t ok = h->GetBinContent(1);
 204     Int_t nerr = h->GetNbinsX();
 205     for(Int_t ierr=nerr; ierr--;){
 206       h->SetBinContent(ierr+1, ok>0.?1.e2*h->GetBinContent(ierr+1)/ok:0.);
 207     }
 208     h->SetBinContent(1, 0.);
 209   }
 210   // tracklet status
 211
 212   TObjArray *arr(NULL);
 213   if(( arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus)))){
 214     for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
 215       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
 216       Float_t okB = h->GetBinContent(1);
 217       Int_t nerrB = h->GetNbinsX();
 218       for(Int_t ierr=nerrB; ierr--;){
 219         h->SetBinContent(ierr+1, okB>0.?1.e2*h->GetBinContent(ierr+1)/okB:0.);
 220       }
 221       h->SetBinContent(1, 0.);
 222     }
 223   }
 224   fNRefFigures = 17;
 225
 226   return kTRUE;
 227 }
 228
 229 //_______________________________________________________
 230 void AliTRDcheckDET::MakeSummary(){
 231   //
 232   // Create summary plots for TRD check DET
 233   // This function creates 2 summary plots:
 234   // - General Quantities
 235   // - PHS
 236   // The function will reuse GetRefFigure
 237   //
 238
 239   TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);
 240   cOut->Divide(3,3);
 241
 242   // Create figures using GetRefFigure
 243   cOut->cd(1); GetRefFigure(kFigNtracksEvent);
 244   cOut->cd(2); GetRefFigure(kFigNtracksSector);
 245   cOut->cd(3); GetRefFigure(kFigNclustersTrack);
 246   cOut->cd(4); GetRefFigure(kFigNclustersTracklet);
 247   cOut->cd(5); GetRefFigure(kFigNtrackletsTrack);
 248   cOut->cd(6); GetRefFigure(kFigNTrackletsP);
 249   cOut->cd(7); GetRefFigure(kFigChargeCluster);
 250   cOut->cd(8); GetRefFigure(kFigChargeTracklet);
 251   cOut->SaveAs("TRD_TrackQuality.gif");
 252   delete cOut;
 253
 254   // Second Plot: PHS
 255   cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 512);
 256   cOut->cd(); GetRefFigure(kFigPH);
 257   cOut->SaveAs("TRD_PH.gif");
 258   delete cOut;
 259
 260   // Third Plot: Mean Number of clusters as function of eta, phi and layer
 261    cOut = new TCanvas(Form("summary%s3", GetName()), Form("Summary 3 for task %s", GetName()), 1024, 768);
 262   cOut->cd(); MakePlotMeanClustersLayer();
 263   cOut->SaveAs("TRD_MeanNclusters.gif");
 264   delete cOut;
 265
 266 }
 267
 268 //_______________________________________________________
 269 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
 270   //
 271   // Setting Reference Figures
 272   //
 273   gPad->SetLogy(0);
 274   gPad->SetLogx(0);
 275   TH1 *h = NULL; TObjArray *arr=NULL;
 276   TLegend *leg = NULL;
 277   Bool_t kFIRST(1);
 278   switch(ifig){
 279   case kFigNclustersTrack:
 280     (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
 281     PutTrendValue("NClustersTrack", h->GetMean());
 282     PutTrendValue("NClustersTrackRMS", h->GetRMS());
 283     return kTRUE;
 284   case kFigNclustersTracklet:
 285     (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
 286     PutTrendValue("NClustersTracklet", h->GetMean());
 287     PutTrendValue("NClustersTrackletRMS", h->GetRMS());
 288     return kTRUE;
 289   case kFigNtrackletsTrack:
 290     h=MakePlotNTracklets();
 291     if(h){
 292       PutTrendValue("NTrackletsTrack", h->GetMean());
 293       PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
 294     }
 295     return kTRUE;
 296   case kFigNTrackletsP:
 297     MakePlotnTrackletsVsP();
 298     return kTRUE;
 299   case kFigNtrackletsCross:
 300     h = (TH1F*)fContainer->FindObject("hNtlsCross");
 301     if(!MakeBarPlot(h, kRed)) break;
 302     PutTrendValue("NTrackletsCross", h->GetMean());
 303     PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
 304     return kTRUE;
 305   case kFigNtrackletsFindable:
 306     h = (TH1F*)fContainer->FindObject("hNtlsFindable");
 307     if(!MakeBarPlot(h, kGreen)) break;
 308     PutTrendValue("NTrackletsFindable", h->GetMean());
 309     PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
 310     return kTRUE;
 311   case kFigNtracksEvent:
 312     (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
 313     PutTrendValue("NTracksEvent", h->GetMean());
 314     PutTrendValue("NTracksEventRMS", h->GetRMS());
 315     return kTRUE;
 316   case kFigNtracksSector:
 317     h = (TH1F*)fContainer->FindObject("hNtrksSector");
 318     if(!MakeBarPlot(h, kGreen)) break;
 319     PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
 320     return kTRUE;
 321   case kFigTrackStatus:
 322     if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
 323     h->GetXaxis()->SetRangeUser(0.5, -1);
 324     h->GetYaxis()->CenterTitle();
 325     h->Draw("c");
 326     PutTrendValue("TrackStatus", h->Integral());
 327     gPad->SetLogy(0);
 328     return kTRUE;
 329   case kFigTrackletStatus:
 330     if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
 331     leg = new TLegend(.68, .7, .97, .97);
 332     leg->SetBorderSize(0);leg->SetFillStyle(0);
 333     leg->SetHeader("TRD layer");
 334     for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 335       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
 336       if(kFIRST){
 337         h->Draw("pl");
 338         h->GetXaxis()->SetRangeUser(0.5, -1);
 339         h->GetYaxis()->CenterTitle();
 340         kFIRST = kFALSE;
 341       } else h->Draw("samepl");
 342       leg->AddEntry(h, Form("ly = %d", ily), "l");
 343       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
 344     }
 345     leg->Draw();
 346     gPad->SetLogy(0);
 347     return kTRUE;
 348   case kFigChi2:
 349     MakePlotChi2();
 350     return kTRUE;
 351   case kFigPH:
 352     gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
 353     MakePlotPulseHeight();
 354     gPad->SetLogy(0);
 355     return kTRUE;
 356   case kFigChargeCluster:
 357     (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
 358     gPad->SetLogy(1);
 359     PutTrendValue("ChargeCluster", h->GetMaximumBin());
 360     PutTrendValue("ChargeClusterRMS", h->GetRMS());
 361     return kTRUE;
 362   case kFigChargeTracklet:
 363     (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
 364     PutTrendValue("ChargeTracklet", h->GetMaximumBin());
 365     PutTrendValue("ChargeTrackletRMS", h->GetRMS());
 366     return kTRUE;
 367   case kFigNeventsTrigger:
 368     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
 369     return kTRUE;
 370   case kFigNeventsTriggerTracks:
 371     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
 372     return kTRUE;
 373   case kFigTriggerPurity:
 374     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
 375     break;
 376   default:
 377     break;
 378   }
 379   AliInfo(Form("Reference plot [%d] missing result", ifig));
 380   return kFALSE;
 381 }
 382
 383 //_______________________________________________________
 384 TObjArray *AliTRDcheckDET::Histos(){
 385   //
 386   // Create QA histograms
 387   //
 388
 389   if(fContainer) return fContainer;
 390
 391   fContainer = new TObjArray(20);
 392   //fContainer->SetOwner(kTRUE);
 393
 394   // Register Histograms
 395   TH1 * h = NULL;
 396   TH2 * h2 = NULL;      // Pointer for two dimensional histograms
 397   TH3 * h3 = NULL;      // Pointer for tree dimensional histograms
 398   TAxis *ax = NULL;
 399   if(!(h2 = (TH2F *)gROOT->FindObject("hNcls"))){
 400     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);
 401   } else h2->Reset();
 402   fContainer->AddAt(h2, kNclustersTrack);
 403
 404   TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
 405   arr->SetOwner(kTRUE);  arr->SetName("clusters");
 406   fContainer->AddAt(arr, kNclustersLayer);
 407   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 408     if(!(h = (TProfile2D *)gROOT->FindObject(Form("hNcl%d", ily)))){
 409       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());
 410       h->GetXaxis()->SetTitle("#eta");
 411       h->GetYaxis()->SetTitle("#phi");
 412     } else h->Reset();
 413     arr->AddAt(h, ily);
 414   }
 415
 416   if(!(h2 = (TH2F *)gROOT->FindObject("hNclTls"))){
 417     h2 = new TH2F("hNclTls","N_{clusters}/tracklet;N_{clusters};Entries", 51, -0.5, 50.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 418   } else h2->Reset();
 419   fContainer->AddAt(h2, kNclustersTracklet);
 420
 421   if(!(h2 = (TH2F *)gROOT->FindObject("hNtls"))){
 422     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);
 423   } else h2->Reset();
 424   fContainer->AddAt(h2, kNtrackletsTrack);
 425
 426   if(!(h = (TH2F *)gROOT->FindObject("htlsSTA"))){
 427     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);
 428   }
 429   fContainer->AddAt(h, kNtrackletsSTA);
 430
 431   // Binning for momentum dependent tracklet Plots
 432   const Int_t kNp(30);
 433   Float_t p=0.2;
 434   Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1], binsCent[AliTRDeventInfo::kCentralityClasses+2];
 435   for(Int_t i=0;i<kNp+1; i++,p+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=p;
 436   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
 437   for(Int_t icent = -1; icent < AliTRDeventInfo::kCentralityClasses + 1; icent++) binsCent[icent+1] = icent - 0.5;
 438   if(!(h3 = (TH3F *)gROOT->FindObject("htlsBAR"))){
 439     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
 440     h3 = new TH3F("hNtlsBAR",
 441     "N_{tracklets}/track;p [GeV/c];N^{tracklet};freq. [%]",
 442     kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt, AliTRDeventInfo::kCentralityClasses + 1, binsCent);
 443   }
 444   fContainer->AddAt(h3, kNtrackletsBAR);
 445
 446   //
 447   if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsCross"))){
 448     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);
 449   } else h2->Reset();
 450   fContainer->AddAt(h2, kNtrackletsCross);
 451
 452   if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsFindable"))){
 453     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);
 454   } else h2->Reset();
 455   fContainer->AddAt(h2, kNtrackletsFindable);
 456
 457   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
 458     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
 459     h->GetXaxis()->SetTitle("N_{tracks}");
 460     h->GetYaxis()->SetTitle("Entries");
 461   } else h->Reset();
 462   fContainer->AddAt(h, kNtracksEvent);
 463
 464   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
 465     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
 466     h->GetXaxis()->SetTitle("sector");
 467     h->GetYaxis()->SetTitle("freq. [%]");
 468   } else h->Reset();
 469   fContainer->AddAt(h, kNtracksSector);
 470
 471   if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
 472     const Int_t nerr = 7;
 473     h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
 474     const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
 475     ax = h->GetXaxis();
 476     for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 477     h->SetYTitle("Relative Error to Good [%]");
 478   }
 479   fContainer->AddAt(h, kTrackStatus);
 480
 481   arr = new TObjArray(AliTRDgeometry::kNlayer);
 482   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
 483   fContainer->AddAt(arr, kTrackletStatus);
 484   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 485     if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
 486       const Int_t nerr = 8;
 487       h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
 488       h->SetLineColor(ily+1);
 489       const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
 490       ax = h->GetXaxis();
 491       for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 492       h->SetYTitle("Relative Error to Good [%]");
 493     } else h->Reset();
 494     arr->AddAt(h, ily);
 495   }
 496
 497   // <PH> histos
 498   arr = new TObjArray(3);
 499   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
 500   fContainer->AddAt(arr, kPH);
 501   if(!(h3 = (TH3F *)gROOT->FindObject("hPHx"))){
 502     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);
 503   } else h3->Reset();
 504   arr->AddAt(h3, 0);
 505   if(!(h3 = (TH3F *)gROOT->FindObject("hPHt"))){
 506     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);
 507   } else h3->Reset();
 508   arr->AddAt(h3, 1);
 509
 510   // Chi2 histos
 511   if(!(h3 = (TH3S*)gROOT->FindObject("hChi2"))){
 512     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);
 513   } else h3->Reset();
 514   fContainer->AddAt(h3, kChi2);
 515
 516   if(!(h2 = (TH2F *)gROOT->FindObject("hQcl"))){
 517     h2 = new TH2F("hQcl", "Q_{cluster};Charge[a.u.];Centrality;Entries", 200, 0, 1200, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 518   }else h2->Reset();
 519   fContainer->AddAt(h2, kChargeCluster);
 520
 521   if(!(h2 = (TH2F *)gROOT->FindObject("hQtrklt"))){
 522     h2 = new TH2F("hQtrklt", "Q_{tracklet};Charge[a.u.];Centrality;Entries", 6000, 0, 6000, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 523   }else h2->Reset();
 524   fContainer->AddAt(h2, kChargeTracklet);
 525
 526
 527   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 528     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 529   else h->Reset();
 530
 531   fContainer->AddAt(h, kNeventsTrigger);
 532
 533   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 534     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 535   else h->Reset();
 536   fContainer->AddAt(h, kNeventsTriggerTracks);
 537
 538   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
 539     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
 540     h->GetXaxis()->SetTitle("Trigger Cluster");
 541     h->GetYaxis()->SetTitle("freq.");
 542   } else h->Reset();
 543   fContainer->AddAt(h, kTriggerPurity);
 544
 545   return fContainer;
 546 }
 547
 548 /*
 549 * Plotting Functions
 550 */
 551
 552 //_______________________________________________________
 553 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
 554 {
 555 //
 556 // Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
 557 //   PROL - track prolongation failure
 558 //   PROP - track propagation failure
 559 //   AJST - crossing sectors failure
 560 //   SNP  - too large bending
 561 //   TINI - tracklet initialization failure
 562 //   UPDT - track position/covariance update failure
 563 //
 564 // Performance plot looks as below:
 565 //Begin_Html
 566 //<img src="TRD/trackStatus.gif">
 567 //End_Html
 568 //
 569   if(track) fkTrack = track;
 570   if(!fkTrack){
 571     AliDebug(4, "No Track defined.");
 572     return NULL;
 573   }
 574   TH1 *h = NULL;
 575   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
 576     AliWarning("No Histogram defined.");
 577     return NULL;
 578   }
 579   h->Fill(fkTrack->GetStatusTRD());
 580   return h;
 581 }
 582
 583 //_______________________________________________________
 584 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
 585 {
 586 //
 587 // Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
 588 //   Geom   -
 589 //   Bound  - tracklet too close to chamber walls
 590 //   NoCl   - no clusters in the track roads
 591 //   NoAttach - fail to attach clusters
 592 //   NoClTr - fail to use clusters for fit
 593 //   NoFit  - tracklet fit failled
 594 //   Chi2   - chi2 tracklet-track over threshold
 595 //
 596 // Performance plot looks as below:
 597 //Begin_Html
 598 //<img src="TRD/trackletStatus.gif">
 599 //End_Html
 600 //
 601   if(track) fkTrack = track;
 602   if(!fkTrack){
 603     AliDebug(4, "No Track defined.");
 604     return NULL;
 605   }
 606   TObjArray *arr =NULL;
 607   if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
 608     AliWarning("Histograms not defined.");
 609     return NULL;
 610   }
 611
 612   TH1 *h = NULL;
 613   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 614     if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
 615       AliWarning(Form("Missing histo for layer %d.", ily));
 616       continue;
 617     }
 618     h->Fill(fkTrack->GetStatusTRD(ily));
 619   }
 620   return h;
 621 }
 622
 623 //_______________________________________________________
 624 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
 625   //
 626   // Plot the mean number of clusters per tracklet
 627   //
 628   if(track) fkTrack = track;
 629   if(!fkTrack){
 630     AliDebug(4, "No Track defined.");
 631     return NULL;
 632   }
 633   AliExternalTrackParam *par = fkTrack->GetTrackIn() ? fkTrack->GetTrackIn() : fkTrack->GetTrackOut();
 634   TH2 *h = NULL;
 635   TProfile2D *hlayer = NULL;
 636   Double_t eta = 0., phi = 0.;
 637   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTracklet)))){
 638     AliWarning("No Histogram defined.");
 639     return NULL;
 640   }
 641   AliTRDseedV1 *tracklet = NULL;
 642   TObjArray *histosLayer = dynamic_cast<TObjArray *>(fContainer->At(kNclustersLayer));
 643   if(!histosLayer){
 644     AliWarning("No Histograms for single layer defined");
 645   }
 646   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 647     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 648     h->Fill(tracklet->GetN2(), fCentralityClass);
 649     if(histosLayer && par){
 650       if((hlayer = dynamic_cast<TProfile2D *>(histosLayer->At(itl)))){
 651         GetEtaPhiAt(par, tracklet->GetX0(), eta, phi);
 652         hlayer->Fill(eta, phi, tracklet->GetN2());
 653       }
 654     }
 655   }
 656   return h;
 657 }
 658
 659 //_______________________________________________________
 660 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
 661   //
 662   // Plot the number of clusters in one track
 663   //
 664   if(track) fkTrack = track;
 665   if(!fkTrack){
 666     AliDebug(4, "No Track defined.");
 667     return NULL;
 668   }
 669   TH2 *h = NULL;
 670   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTrack)))){
 671     AliWarning("No Histogram defined.");
 672     return NULL;
 673   }
 674
 675   Int_t nclusters = 0;
 676   AliTRDseedV1 *tracklet = NULL;
 677   AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
 678   if(!par) return NULL;
 679   Double_t momentumRec = par->P();
 680   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 681     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 682     Int_t n(tracklet->GetN());
 683     nclusters += n;
 684     if(DebugLevel() > 2){
 685       Int_t crossing = Int_t(tracklet->IsRowCross());
 686       Int_t detector = tracklet->GetDetector();
 687       Float_t theta = TMath::ATan(tracklet->GetZref(1));
 688       Float_t phi = TMath::ATan(tracklet->GetYref(1));
 689       Float_t momentumMC = 0.;
 690       Int_t pdg = 0;
 691       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
 692       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
 693       if(fkMC){
 694         if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
 695         pdg = fkMC->GetPDG();
 696       }
 697       (*DebugStream()) << "NClustersTrack"
 698         << "Detector="  << detector
 699         << "Centrality="<< fCentralityClass
 700         << "crossing="  << crossing
 701         << "momentumMC="<< momentumMC
 702         << "momentumRec="<< momentumRec
 703         << "pdg="                               << pdg
 704         << "theta="                     << theta
 705         << "phi="                               << phi
 706         << "kinkIndex=" << kinkIndex
 707         << "TPCncls="           << nclsTPC
 708         << "TRDncls="   << n
 709         << "\n";
 710     }
 711   }
 712   h->Fill(nclusters, fCentralityClass);
 713   return h;
 714 }
 715
 716
 717 //_______________________________________________________
 718 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
 719   //
 720   // Plot the number of tracklets
 721   //
 722   if(track) fkTrack = track;
 723   if(!fkTrack){
 724     AliDebug(4, "No Track defined.");
 725     return NULL;
 726   }
 727   TH2 *h = NULL, *hSta = NULL; TH3 *hBarrel = NULL;
 728   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsTrack)))){
 729     AliWarning("No Histogram defined.");
 730     return NULL;
 731   }
 732   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
 733   h->Fill(nTracklets, fCentralityClass);
 734   if(!fkESD) return h;
 735   Int_t status = fkESD->GetStatus();
 736
 737 /*  printf("in/out/refit/pid: TRD[%d|%d|%d|%d]\n", status &AliESDtrack::kTRDin ? 1 : 0, status &AliESDtrack::kTRDout ? 1 : 0, status &AliESDtrack::kTRDrefit ? 1 : 0, status &AliESDtrack::kTRDpid ? 1 : 0);*/
 738   Double_t p = 0.;
 739   Int_t method = -1;    // to distinguish between stand alone and full barrel tracks in the debugging
 740   if((status & AliESDtrack::kTRDin) != 0){
 741     method = 1;
 742     // Full Barrel Track: Save momentum dependence
 743     if(!(hBarrel = dynamic_cast<TH3F *>(fContainer->At(kNtrackletsBAR)))){
 744       AliWarning("Method: Barrel.  Histogram not processed!");
 745       return NULL;
 746     }
 747     AliExternalTrackParam *par(fkTrack->GetTrackIn());
 748     if(!par){
 749       AliError("Input track params missing");
 750       return NULL;
 751     }
 752     p = par->P(); // p needed later in the debug streaming
 753     hBarrel->Fill(p, nTracklets, fCentralityClass);
 754   } else {
 755     // Stand alone Track: momentum dependence not usefull
 756     method = 0;
 757     if(!(hSta = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsSTA)))) {
 758       AliWarning("Method: StandAlone.  Histogram not processed!");
 759       return NULL;
 760     }
 761     hSta->Fill(nTracklets, fCentralityClass);
 762   }
 763
 764   if(DebugLevel() > 2){
 765     AliTRDseedV1 *tracklet = NULL;
 766     AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 767     Int_t sector = -1, stack = -1, detector;
 768     for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 769       if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
 770       detector = tracklet->GetDetector();
 771       sector = geo->GetSector(detector);
 772       stack = geo->GetStack(detector);
 773       break;
 774     }
 775     (*DebugStream()) << "NTrackletsTrack"
 776       << "Sector="      << sector
 777       << "Stack="       << stack
 778       << "Centrality="  << fCentralityClass
 779       << "NTracklets="  << nTracklets
 780       << "Method="      << method
 781       << "p="           << p
 782       << "\n";
 783   }
 784   if(DebugLevel() > 3){
 785     AliTRDseedV1 *tracklet = NULL;
 786     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 787       if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
 788         (*DebugStream()) << "NTrackletsLayer"
 789         << "Layer=" << il
 790         << "p=" << p
 791         << "\n";
 792       }
 793     }
 794   }
 795   return h;
 796 }
 797
 798
 799 //_______________________________________________________
 800 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
 801   //
 802   // Plot the number of tracklets
 803   //
 804   if(track) fkTrack = track;
 805   if(!fkTrack){
 806     AliDebug(4, "No Track defined.");
 807     return NULL;
 808   }
 809   TH2 *h = NULL;
 810   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsCross)))){
 811     AliWarning("No Histogram defined.");
 812     return NULL;
 813   }
 814
 815   Int_t ncross = 0;
 816   AliTRDseedV1 *tracklet = NULL;
 817   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 818     if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
 819
 820     if(tracklet->IsRowCross()) ncross++;
 821   }
 822   h->Fill(ncross, fCentralityClass);
 823   return h;
 824 }
 825
 826 //_______________________________________________________
 827 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
 828   //
 829   // Plots the ratio of number of tracklets vs.
 830   // number of findable tracklets
 831   //
 832   // Findable tracklets are defined as track prolongation
 833   // to layer i does not hit the dead area +- epsilon
 834   //
 835   // In order to check whether tracklet hist active area in Layer i,
 836   // the track is refitted and the fitted position + an uncertainty
 837   // range is compared to the chamber border (also with a different
 838   // uncertainty)
 839   //
 840   // For the track fit two cases are distinguished:
 841   // If the track is a stand alone track (defined by the status bit
 842   // encoding, then the track is fitted with the tilted Rieman model
 843   // Otherwise the track is fitted with the Kalman fitter in two steps:
 844   // Since the track parameters are give at the outer point, we first
 845   // fit in direction inwards. Afterwards we fit again in direction outwards
 846   // to extrapolate the track to layers which are not reached by the track
 847   // For the Kalman model, the radial track points have to be shifted by
 848   // a distance epsilon in the direction that we want to fit
 849   //
 850   const Float_t epsilon = 0.01;   // dead area tolerance
 851   const Float_t epsilonR = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
 852   const Float_t deltaY = 0.7;    // Tolerance in the track position in y-direction
 853   const Float_t deltaZ = 7.0;    // Tolerance in the track position in z-direction (Padlength)
 854   Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 855
 856   if(track) fkTrack = track;
 857   if(!fkTrack){
 858     AliDebug(4, "No Track defined.");
 859     return NULL;
 860   }
 861   TH2 *h = NULL;
 862   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsFindable)))){
 863     AliWarning("No Histogram defined.");
 864     return NULL;
 865   }
 866   Int_t nFound = 0, nFindable = 0;
 867   Int_t stack = -1;
 868   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 869   Double_t y = 0., z = 0.;
 870   AliTRDseedV1 *tracklet = NULL;
 871   AliTRDpadPlane *pp;
 872   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 873     if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
 874       tracklet->SetReconstructor(AliTRDinfoGen::Reconstructor());
 875       nFound++;
 876     }
 877   }
 878   // 2 Different cases:
 879   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 880   // 2nd barrel track: here we propagate the track to the layers
 881   AliTrackPoint points[6];
 882   Float_t xyz[3];
 883   memset(xyz, 0, sizeof(Float_t) * 3);
 884   if(((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 885     // stand alone track
 886     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 887       xyz[0] = xAnode[il];
 888       points[il].SetXYZ(xyz);
 889     }
 890     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
 891   } else {
 892     // barrel track
 893     //
 894     // 2 Steps:
 895     // -> Kalman inwards
 896     // -> Kalman outwards
 897     AliTRDtrackV1 copyTrack(*fkTrack);  // Do Kalman on a (non-constant) copy of the track
 898     AliTrackPoint pointsInward[6], pointsOutward[6];
 899     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 900       // In order to avoid complications in the Kalman filter if the track points have the same radial
 901       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
 902       // in the direction we want to go
 903       // The track points have to be in reverse order for the Kalman Filter inwards
 904       xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
 905       pointsInward[il].SetXYZ(xyz);
 906       xyz[0] = xAnode[il] + epsilonR;
 907       pointsOutward[il].SetXYZ(xyz);
 908     }
 909     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 910       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
 911     // Kalman inwards
 912     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
 913     memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
 914     // Kalman outwards
 915     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
 916     memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 917   }
 918   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 919   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 920     y = points[il].GetY();
 921     z = points[il].GetZ();
 922     if((stack = geo->GetStack(z, il)) < 0) continue; // Not findable
 923     pp = geo->GetPadPlane(il, stack);
 924     ymin = pp->GetCol0() + epsilon;
 925     ymax = pp->GetColEnd() - epsilon;
 926     zmin = pp->GetRowEnd() + epsilon;
 927     zmax = pp->GetRow0() - epsilon;
 928     // ignore y-crossing (material)
 929     if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
 930       if(DebugLevel() > 3){
 931         Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
 932         Int_t hasTracklet = tracklet ? 1 : 0;
 933         (*DebugStream())   << "FindableTracklets"
 934           << "layer="     << il
 935           << "ytracklet=" << posTracklet[0]
 936           << "ytrack="    << y
 937           << "ztracklet=" << posTracklet[1]
 938           << "ztrack="    << z
 939           << "tracklet="  << hasTracklet
 940           << "\n";
 941       }
 942   }
 943
 944   h->Fill((nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1), fCentralityClass);
 945   AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 946   return h;
 947 }
 948
 949
 950 //_______________________________________________________
 951 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
 952   //
 953   // Plot the chi2 of the track
 954   //
 955   if(track) fkTrack = track;
 956   if(!fkTrack){
 957     AliDebug(4, "No Track defined.");
 958     return NULL;
 959   }
 960   TH3 *h = NULL;
 961   if(!(h = dynamic_cast<TH3S*>(fContainer->At(kChi2)))) {
 962     AliWarning("No Histogram defined.");
 963     return NULL;
 964   }
 965   Int_t n = fkTrack->GetNumberOfTracklets();
 966   if(!n) return NULL;
 967
 968   h->Fill(n, fkTrack->GetChi2()/n, fCentralityClass);
 969   return h;
 970 }
 971
 972
 973 //_______________________________________________________
 974 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
 975   //
 976   // Plot the average pulse height
 977   //
 978   if(track) fkTrack = track;
 979   if(!fkTrack){
 980     AliDebug(4, "No Track defined.");
 981     return NULL;
 982   }
 983   TH3F *h = NULL;
 984   if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
 985     AliWarning("No Histogram defined.");
 986     return NULL;
 987   }
 988   AliTRDseedV1 *tracklet = NULL;
 989   AliTRDcluster *c = NULL;
 990   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 991     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 992     Int_t crossing = Int_t(tracklet->IsRowCross());
 993     Int_t detector = tracklet->GetDetector();
 994     tracklet->ResetClusterIter();
 995     while((c = tracklet->NextCluster())){
 996       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
 997       Int_t localtime        = c->GetLocalTimeBin();
 998       Double_t absoluteCharge = TMath::Abs(c->GetQ());
 999       h->Fill(localtime, absoluteCharge, fCentralityClass);
1000       if(DebugLevel() > 3){
1001         Int_t inChamber = c->IsInChamber() ? 1 : 0;
1002         Double_t distance[2];
1003         GetDistanceToTracklet(distance, tracklet, c);
1004         Float_t theta = TMath::ATan(tracklet->GetZref(1));
1005         Float_t phi = TMath::ATan(tracklet->GetYref(1));
1006         AliExternalTrackParam *trdPar = fkTrack->GetTrackIn();
1007         Float_t momentumMC = 0, momentumRec = trdPar ? trdPar->P() : fkTrack->P(); // prefer Track Low
1008         Int_t pdg = 0;
1009         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1010         UShort_t tpcNCLS = fkESD ? fkESD->GetTPCncls() : 0;
1011         if(fkMC){
1012           if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
1013           pdg = fkMC->GetPDG();
1014         }
1015         (*DebugStream()) << "PHt"
1016           << "Detector="        << detector
1017           << "Centrality="<< fCentralityClass
1018           << "crossing="        << crossing
1019           << "inChamber=" << inChamber
1020           << "Timebin="         << localtime
1021           << "Charge="          << absoluteCharge
1022           << "momentumMC="      << momentumMC
1023           << "momentumRec="     << momentumRec
1024           << "pdg="                             << pdg
1025           << "theta="                   << theta
1026           << "phi="                             << phi
1027           << "kinkIndex="       << kinkIndex
1028           << "TPCncls="         << tpcNCLS
1029           << "dy="        << distance[0]
1030           << "dz="        << distance[1]
1031           << "c.="        << c
1032           << "\n";
1033       }
1034     }
1035   }
1036   return h;
1037 }
1038
1039 //_______________________________________________________
1040 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1041   //
1042   // Plots the average pulse height vs the distance from the anode wire
1043   // (plus const anode wire offset)
1044   //
1045   if(track) fkTrack = track;
1046   if(!fkTrack){
1047      AliDebug(4, "No Track defined.");
1048      return NULL;
1049   }
1050   TH3 *h = NULL;
1051   if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
1052     AliWarning("No Histogram defined.");
1053     return NULL;
1054   }
1055   AliTRDseedV1 *tracklet(NULL);
1056   AliTRDcluster *c(NULL);
1057   Double_t xd(0.), dqdl(0.);
1058   TVectorD vq(AliTRDseedV1::kNtb), vxd(AliTRDseedV1::kNtb), vdqdl(AliTRDseedV1::kNtb);
1059   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1060     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1061     Int_t det(tracklet->GetDetector());
1062     Bool_t rc(tracklet->IsRowCross());
1063     for(Int_t ic(0); ic<AliTRDseedV1::kNtb; ic++){
1064       Bool_t kFIRST(kFALSE);
1065       if(!(c = tracklet->GetClusters(ic))){
1066          if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1067       } else kFIRST=kTRUE;
1068       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1069       xd = tracklet->GetX0() - c->GetX(); vxd[ic] = xd;
1070       dqdl=tracklet->GetdQdl(ic); vdqdl[ic] = dqdl;
1071       vq[ic]=c->GetQ();
1072       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) vq[ic]+=c->GetQ();
1073       h->Fill(xd, dqdl, fCentralityClass);
1074     }
1075     if(DebugLevel() > 3){
1076       (*DebugStream()) << "PHx"
1077         << "det="  << det
1078         << "rc="   << rc
1079         << "xd="   << &vxd
1080         << "q="    << &vq
1081         << "dqdl=" << &vdqdl
1082         << "\n";
1083     }
1084   }
1085   return h;
1086 }
1087
1088 //_______________________________________________________
1089 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1090   //
1091   // Plot the cluster charge
1092   //
1093   if(track) fkTrack = track;
1094   if(!fkTrack){
1095     AliDebug(4, "No Track defined.");
1096     return NULL;
1097   }
1098   TH2 *h = NULL;
1099   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeCluster)))){
1100     AliWarning("No Histogram defined.");
1101     return NULL;
1102   }
1103   AliTRDseedV1 *tracklet = NULL;
1104   AliTRDcluster *c = NULL;
1105   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1106     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1107     for(Int_t ic(0); ic < AliTRDseedV1::kNtb; ic++){
1108       Bool_t kFIRST(kFALSE);
1109       if(!(c = tracklet->GetClusters(ic))) {
1110         if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1111       } else kFIRST = kTRUE;
1112       Float_t q(c->GetQ());
1113       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) q+=c->GetQ();
1114       h->Fill(q, fCentralityClass);
1115     }
1116   }
1117   return h;
1118 }
1119
1120 //_______________________________________________________
1121 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1122   //
1123   // Plot the charge deposit per chamber
1124   //
1125   if(track) fkTrack = track;
1126   if(!fkTrack){
1127     AliDebug(4, "No Track defined.");
1128     return NULL;
1129   }
1130   TH2 *h = NULL;
1131   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeTracklet)))){
1132     AliWarning("No Histogram defined.");
1133     return NULL;
1134   }
1135   AliTRDseedV1 *tracklet = NULL;
1136   AliTRDcluster *c = NULL;
1137   Double_t qTot = 0;
1138   Int_t nTracklets =fkTrack->GetNumberOfTracklets();
1139   for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1140     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1141     qTot = 0.;
1142     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1143       if(!(c = tracklet->GetClusters(ic))) continue;
1144       qTot += TMath::Abs(c->GetQ());
1145     }
1146     h->Fill(qTot, fCentralityClass);
1147     if(DebugLevel() > 3){
1148       Int_t crossing = (Int_t)tracklet->IsRowCross();
1149       Int_t detector = tracklet->GetDetector();
1150       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1151       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1152       Float_t momentum = 0.;
1153       Int_t pdg = 0;
1154       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1155       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1156       if(fkMC){
1157               if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1158         pdg = fkMC->GetPDG();
1159       }
1160       (*DebugStream()) << "ChargeTracklet"
1161         << "Detector="  << detector
1162         << "Centrality="<< fCentralityClass
1163         << "crossing="  << crossing
1164         << "momentum="  << momentum
1165         << "nTracklets="<< nTracklets
1166         << "pdg="                               << pdg
1167         << "theta="                     << theta
1168         << "phi="                               << phi
1169         << "kinkIndex=" << kinkIndex
1170         << "TPCncls="           << nclsTPC
1171         << "QT="        << qTot
1172         << "\n";
1173     }
1174   }
1175   return h;
1176 }
1177
1178 //_______________________________________________________
1179 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1180   //
1181   // Plot the number of tracks per Sector
1182   //
1183   if(track) fkTrack = track;
1184   if(!fkTrack){
1185     AliDebug(4, "No Track defined.");
1186     return NULL;
1187   }
1188   TH1 *h = NULL;
1189   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1190     AliWarning("No Histogram defined.");
1191     return NULL;
1192   }
1193
1194   // TODO we should compare with
1195   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1196
1197   AliTRDseedV1 *tracklet = NULL;
1198   Int_t sector = -1;
1199   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1200     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1201     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1202     break;
1203   }
1204   h->Fill(sector);
1205   return h;
1206 }
1207
1208
1209 //________________________________________________________
1210 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1211 {
1212   Float_t x = c->GetX();
1213   dist[0] = c->GetY() - tracklet->GetYat(x);
1214   dist[1] = c->GetZ() - tracklet->GetZat(x);
1215 }
1216
1217 //________________________________________________________
1218 void AliTRDcheckDET::GetEtaPhiAt(const AliExternalTrackParam *track, Double_t x, Double_t &eta, Double_t &phi){
1219   //
1220   // Get phi and eta at a given radial position
1221   //
1222   AliExternalTrackParam workpar(*track);
1223
1224   Double_t posLocal[3];
1225   Bool_t sucPos = workpar.GetXYZAt(x, fEvent->GetRunInfo()->GetMagneticField(), posLocal);
1226   Double_t sagPhi = sucPos ? TMath::ATan2(posLocal[1], posLocal[0]) : 0.;
1227   phi = sagPhi;
1228   eta = workpar.Eta();
1229 }
1230
1231
1232 //_______________________________________________________
1233 TH1* AliTRDcheckDET::MakePlotChi2() const
1234 {
1235 // Plot chi2/track normalized to number of degree of freedom
1236 // (tracklets) and compare with the theoretical distribution.
1237 //
1238 // Alex Bercuci <A.Bercuci@gsi.de>
1239
1240   return NULL;
1241
1242 /*  TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1243   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1244   f.SetParLimits(1,1, 1e100);
1245   TLegend *leg = new TLegend(.7,.7,.95,.95);
1246   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
1247   TH1D *h1 = NULL;
1248   Bool_t kFIRST = kTRUE;
1249   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1250     h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1251     if(h1->Integral()<50) continue;
1252     h1->Scale(1./h1->Integral());
1253     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1254     h1->SetLineColor(il);h1->SetLineStyle(2);
1255     f.SetParameter(1, .5*il);f.SetLineColor(il);
1256     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1257     leg->AddEntry(h1, Form("%d", il), "l");
1258     if(kFIRST){
1259       h1->GetXaxis()->SetRangeUser(0., 25.);
1260     }
1261     kFIRST = kFALSE;
1262   }
1263   leg->Draw();
1264   gPad->SetLogy();
1265   return h1;*/
1266 }
1267
1268
1269 //________________________________________________________
1270 TH1* AliTRDcheckDET::MakePlotNTracklets(){
1271   //
1272   // Make nice bar plot of the number of tracklets in each method
1273   //
1274   TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1275   TH1D *hBAR = tmp->ProjectionY();
1276   TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1277   TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1278   TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1279   leg->SetBorderSize(1);
1280   leg->SetFillColor(0);
1281
1282   Float_t scale = hCON->Integral();
1283   if(scale) hCON->Scale(100./scale);
1284   hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1285   hCON->SetBarWidth(0.2);
1286   hCON->SetBarOffset(0.6);
1287   hCON->SetStats(kFALSE);
1288   hCON->GetYaxis()->SetRangeUser(0.,40.);
1289   hCON->GetYaxis()->SetTitleOffset(1.2);
1290   hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1291   hCON->SetMaximum(55.);
1292
1293   if(scale) hBAR->Scale(100./scale);
1294   hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1295   hBAR->SetBarWidth(0.2);
1296   hBAR->SetBarOffset(0.2);
1297   hBAR->SetTitle("");
1298   hBAR->SetStats(kFALSE);
1299   hBAR->GetYaxis()->SetRangeUser(0.,40.);
1300   hBAR->GetYaxis()->SetTitleOffset(1.2);
1301   hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1302
1303   if(scale) hSTA->Scale(100./scale);
1304   hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1305   hSTA->SetBarWidth(0.2);
1306   hSTA->SetBarOffset(0.4);
1307   hSTA->SetTitle("");
1308   hSTA->SetStats(kFALSE);
1309   hSTA->GetYaxis()->SetRangeUser(0.,40.);
1310   hSTA->GetYaxis()->SetTitleOffset(1.2);
1311   hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1312   leg->Draw();
1313   gPad->Update();
1314   return hCON;
1315 }
1316
1317 //________________________________________________________
1318 void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1319   //
1320   // Plot abundance of tracks with number of tracklets as function of momentum
1321   //
1322
1323
1324
1325
1326   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1327   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1328   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1329     g[itl] = new TGraphErrors();
1330     g[itl]->SetLineColor(color[itl]);
1331     g[itl]->SetMarkerColor(color[itl]);
1332     g[itl]->SetMarkerStyle(20 + itl);
1333   }
1334
1335   TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
1336   TAxis *ax(hBar->GetXaxis());
1337   Int_t np(ax->GetNbins());
1338   for(Int_t ipBin(1); ipBin<np; ipBin++){
1339     h = hBar->ProjectionY("npBin", ipBin, ipBin);
1340     if(!Int_t(h->Integral())) continue;
1341     h->Scale(100./h->Integral());
1342     Float_t p(ax->GetBinCenter(ipBin));
1343     Float_t dp(ax->GetBinWidth(ipBin));
1344     Int_t ip(g[0]->GetN());
1345     for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1346       g[itl]->SetPoint(ip, p, h->GetBinContent(itl+1));
1347       g[itl]->SetPointError(ip, dp/3.46, h->GetBinError(itl+1));
1348     }
1349   }
1350
1351   TLegend *leg = new TLegend(0.76, 0.6, 1., 0.9);
1352   leg->SetBorderSize(0);
1353   leg->SetHeader("Tracklet/Track");
1354   leg->SetFillStyle(0);
1355   h = hBar->ProjectionX("npxBin"); h->Reset();
1356   h->SetTitle("");
1357   h->GetYaxis()->SetRangeUser(1., 99.);
1358   h->GetYaxis()->SetMoreLogLabels();
1359   h->GetYaxis()->CenterTitle();
1360   h->GetYaxis()->SetTitleOffset(1.2);
1361   h->SetYTitle("Prob. [%]");
1362   h->GetXaxis()->SetRangeUser(0.4, 12.);
1363   h->GetXaxis()->SetMoreLogLabels();
1364   h->GetXaxis()->CenterTitle();
1365   h->Draw("p");
1366   for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1367     g[itl]->Draw("pc");
1368     leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
1369   }
1370
1371   leg->Draw();
1372   gPad->SetLogx();gPad->SetLogy();
1373 }
1374
1375 //________________________________________________________
1376 Bool_t AliTRDcheckDET::MakePlotPulseHeight(){
1377   //
1378   // Create Plot of the Pluse Height Spectrum
1379   //
1380   TCanvas *output = gPad->GetCanvas();
1381   output->Divide(2);
1382   output->cd(1);
1383   TH1 *h, *h1, *h2;
1384   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1385   h = (TH1F*)arr->At(0);
1386   if((Int_t)h->GetEntries()){
1387     h->SetMarkerStyle(24);
1388     h->SetMarkerColor(kBlack);
1389     h->SetLineColor(kBlack);
1390     h->GetYaxis()->SetTitleOffset(1.5);
1391     h->Draw("e1");
1392     // Trending for the pulse height: plateau value, slope and timebin of the maximum
1393     TLinearFitter fit(1,"pol1");
1394     Double_t time = 0.;
1395     for(Int_t itime = 10; itime <= 20; itime++){
1396       time = Double_t(itime);
1397       Double_t err(h->GetBinError(itime + 1));
1398       if(err>1.e-10) fit.AddPoint(&time, h->GetBinContent(itime + 1), err);
1399     }
1400     if(!fit.Eval()){
1401       Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1402       Double_t slope = fit.GetParameter(1);
1403       PutTrendValue("PHplateau", plateau);
1404       PutTrendValue("PHslope", slope);
1405       PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1406       AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1407     }
1408   }
1409   //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1410   h1 = (TH1F *)arr->At(1);
1411   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1412   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
1413     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1414   h2->SetMarkerStyle(22);
1415   h2->SetMarkerColor(kBlue);
1416   h2->SetLineColor(kBlue);
1417   h2->Draw("e1same");
1418   gPad->Update();
1419
1420 //   create axis according to the histogram dimensions of the original second histogram
1421   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1422                     gPad->GetUymax(),
1423                     gPad->GetUxmax(),
1424                     gPad->GetUymax(),
1425                     -0.08, 4.88, 510,"-L");
1426   axis->SetLineColor(kBlue);
1427   axis->SetLabelColor(kBlue);
1428   axis->SetTextColor(kBlue);
1429   axis->SetTitle("x_{0}-x_{c} [cm]");
1430   axis->Draw();
1431
1432   output->cd(2);
1433   TH2 *ph2d = (TH2F *)arr->At(2);
1434   ph2d->GetYaxis()->SetTitleOffset(1.8);
1435   ph2d->SetStats(kFALSE);
1436   ph2d->Draw("colz");
1437   return kTRUE;
1438 }
1439
1440 //________________________________________________________
1441 void AliTRDcheckDET::MakePlotMeanClustersLayer(){
1442   //
1443   // Create Summary plot for the mean number of clusters per layer
1444   //
1445   TCanvas *output = gPad->GetCanvas();
1446   output->Divide(3,2);
1447   TObjArray *histos = (TObjArray *)fContainer->At(kNclustersLayer);
1448   if(!histos){
1449     AliWarning("Histos for each layer not found");
1450     return;
1451   }
1452   TProfile2D *hlayer = NULL;
1453   for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
1454     if(!(hlayer = dynamic_cast<TProfile2D *>(histos->At(ily)))) continue;
1455     output->cd(ily + 1);
1456     gPad->SetGrid(0,0);
1457     hlayer->Draw("colz");
1458   }
1459 }
1460
1461 //________________________________________________________
1462 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1463   //
1464   // Draw nice bar plots
1465   //
1466   if(!histo->GetEntries()) return kFALSE;
1467   histo->Scale(100./histo->Integral());
1468   histo->SetFillColor(color);
1469   histo->SetBarOffset(.2);
1470   histo->SetBarWidth(.6);
1471   histo->Draw("bar1");
1472   return kTRUE;
1473 }