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