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