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