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