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 <TObjArray.h>
  41 #include <TObject.h>
  42 #include <TObjString.h>
  43
  44 #include <TPad.h>
  45 #include <TProfile.h>
  46 #include <TProfile2D.h>
  47 #include <TROOT.h>
  48 #include <TChain.h>
  49
  50 #include "AliLog.h"
  51 #include "AliTRDcluster.h"
  52 #include "AliESDHeader.h"
  53 #include "AliESDRun.h"
  54 #include "AliESDtrack.h"
  55 #include "AliExternalTrackParam.h"
  56 #include "AliTRDgeometry.h"
  57 #include "AliTRDpadPlane.h"
  58 #include "AliTRDSimParam.h"
  59 #include "AliTRDseedV1.h"
  60 #include "AliTRDtrackV1.h"
  61 #include "AliTRDtrackerV1.h"
  62 #include "AliTRDReconstructor.h"
  63 #include "AliTrackReference.h"
  64 #include "AliTrackPointArray.h"
  65 #include "AliTracker.h"
  66 #include "TTreeStream.h"
  67
  68 #include "info/AliTRDtrackInfo.h"
  69 #include "info/AliTRDeventInfo.h"
  70 #include "AliTRDcheckDET.h"
  71
  72 #include <cstdio>
  73 #include <iostream>
  74
  75 ClassImp(AliTRDcheckDET)
  76
  77 //_______________________________________________________
  78 AliTRDcheckDET::AliTRDcheckDET():
  79   AliTRDrecoTask()
  80   ,fEventInfo(NULL)
  81   ,fTriggerNames(NULL)
  82   ,fReconstructor(NULL)
  83   ,fGeo(NULL)
  84   ,fFlags(0)
  85 {
  86   //
  87   // Default constructor
  88   //
  89   SetNameTitle("checkDET", "Basic TRD data checker");
  90 }
  91
  92 //_______________________________________________________
  93 AliTRDcheckDET::AliTRDcheckDET(char* name):
  94   AliTRDrecoTask(name, "Basic TRD data checker")
  95   ,fEventInfo(NULL)
  96   ,fTriggerNames(NULL)
  97   ,fReconstructor(NULL)
  98   ,fGeo(NULL)
  99   ,fFlags(0)
 100 {
 101   //
 102   // Default constructor
 103   //
 104   DefineInput(2, AliTRDeventInfo::Class());
 105
 106   fReconstructor = new AliTRDReconstructor;
 107   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
 108   fGeo = new AliTRDgeometry;
 109   InitFunctorList();
 110 }
 111
 112
 113 //_______________________________________________________
 114 AliTRDcheckDET::~AliTRDcheckDET(){
 115   //
 116   // Destructor
 117   //
 118   if(fTriggerNames) delete fTriggerNames;
 119   delete fReconstructor;
 120   delete fGeo;
 121 }
 122
 123 //_______________________________________________________
 124 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
 125   //
 126   // Connect the Input data with the task
 127   //
 128   AliTRDrecoTask::ConnectInputData(opt);
 129   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
 130 }
 131
 132 //_______________________________________________________
 133 void AliTRDcheckDET::UserCreateOutputObjects(){
 134   //
 135   // Create Output Objects
 136   //
 137   OpenFile(1,"RECREATE");
 138   fContainer = Histos();
 139   if(!fTriggerNames) fTriggerNames = new TMap();
 140 }
 141
 142 //_______________________________________________________
 143 void AliTRDcheckDET::UserExec(Option_t *opt){
 144   //
 145   // Execution function
 146   // Filling TRD quality histos
 147   //
 148   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 149   AliTRDrecoTask::UserExec(opt);
 150   Int_t nTracks = 0;            // Count the number of tracks per event
 151   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 152   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 153   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
 154   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
 155   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 156     if(!fTracks->UncheckedAt(iti)) continue;
 157     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 158     if(!fTrackInfo->GetTrack()) continue;
 159     nTracks++;
 160   }
 161
 162   if(nTracks){
 163     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
 164     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
 165   }
 166   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 167     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 168     // also set the label for both histograms
 169     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
 170     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 171     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
 172     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 173   }
 174   PostData(1, fContainer);
 175 }
 176
 177
 178 //_______________________________________________________
 179 Bool_t AliTRDcheckDET::PostProcess(){
 180   //
 181   // Do Postprocessing (for the moment set the number of Reference histograms)
 182   //
 183
 184   TH1 * h = NULL;
 185
 186   // Calculate of the trigger clusters purity
 187   h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
 188   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
 189   h1->Divide(h);
 190   Float_t purities[20], val = 0;
 191   TString triggernames[20];
 192   Int_t nTriggerClasses = 0;
 193   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
 194     if((val = h1->GetBinContent(ibin))){
 195       purities[nTriggerClasses] = val;
 196       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
 197       nTriggerClasses++;
 198     }
 199   }
 200   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
 201   TAxis *ax = h->GetXaxis();
 202   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
 203     h->Fill(itrg, purities[itrg]);
 204     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
 205   }
 206   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
 207   h->GetYaxis()->SetRangeUser(0,1);
 208
 209   // track status
 210   h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
 211   Float_t ok = h->GetBinContent(1);
 212   Int_t nerr = h->GetNbinsX();
 213   for(Int_t ierr=nerr; ierr--;){
 214     h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
 215   }
 216   h->SetBinContent(1, 0.);
 217
 218   // tracklet status
 219
 220   TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
 221   for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
 222     h=dynamic_cast<TH1F*>(arr->At(ily));
 223     Float_t okB = h->GetBinContent(1);
 224     Int_t nerrB = h->GetNbinsX();
 225     for(Int_t ierr=nerrB; ierr--;){
 226       h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
 227     }
 228     h->SetBinContent(1, 0.);
 229   }
 230
 231   fNRefFigures = 17;
 232
 233   return kTRUE;
 234 }
 235
 236 //_______________________________________________________
 237 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
 238   //
 239   // Setting Reference Figures
 240   //
 241   enum FigureType_t{
 242     kFigNclustersTrack,
 243     kFigNclustersTracklet,
 244     kFigNtrackletsTrack,
 245     kFigNTrackletsP,
 246     kFigNtrackletsCross,
 247     kFigNtrackletsFindable,
 248     kFigNtracksEvent,
 249     kFigNtracksSector,
 250     kFigTrackStatus,
 251     kFigTrackletStatus,
 252     kFigChi2,
 253     kFigPH,
 254     kFigChargeCluster,
 255     kFigChargeTracklet,
 256     kFigNeventsTrigger,
 257     kFigNeventsTriggerTracks,
 258     kFigTriggerPurity
 259   };
 260   gPad->SetLogy(0);
 261   gPad->SetLogx(0);
 262   TH1 *h = NULL; TObjArray *arr=NULL;
 263   TLegend *leg = NULL;
 264   Bool_t kFIRST(1);
 265   switch(ifig){
 266   case kFigNclustersTrack:
 267     (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
 268     PutTrendValue("NClustersTrack", h->GetMean());
 269     PutTrendValue("NClustersTrackRMS", h->GetRMS());
 270     return kTRUE;
 271   case kFigNclustersTracklet:
 272     (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
 273     PutTrendValue("NClustersTracklet", h->GetMean());
 274     PutTrendValue("NClustersTrackletRMS", h->GetRMS());
 275     return kTRUE;
 276   case kFigNtrackletsTrack:
 277     h=MakePlotNTracklets();
 278     PutTrendValue("NTrackletsTrack", h->GetMean());
 279     PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
 280     return kTRUE;
 281   case kFigNTrackletsP:
 282     MakePlotnTrackletsVsP();
 283     return kTRUE;
 284   case kFigNtrackletsCross:
 285     h = (TH1F*)fContainer->FindObject("hNtlsCross");
 286     if(!MakeBarPlot(h, kRed)) break;
 287     PutTrendValue("NTrackletsCross", h->GetMean());
 288     PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
 289     return kTRUE;
 290   case kFigNtrackletsFindable:
 291     h = (TH1F*)fContainer->FindObject("hNtlsFindable");
 292     if(!MakeBarPlot(h, kGreen)) break;
 293     PutTrendValue("NTrackletsFindable", h->GetMean());
 294     PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
 295     return kTRUE;
 296   case kFigNtracksEvent:
 297     (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
 298     PutTrendValue("NTracksEvent", h->GetMean());
 299     PutTrendValue("NTracksEventRMS", h->GetRMS());
 300     return kTRUE;
 301   case kFigNtracksSector:
 302     h = (TH1F*)fContainer->FindObject("hNtrksSector");
 303     if(!MakeBarPlot(h, kGreen)) break;
 304     PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
 305     return kTRUE;
 306   case kFigTrackStatus:
 307     if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
 308     h->GetXaxis()->SetRangeUser(0.5, -1);
 309     h->GetYaxis()->CenterTitle();
 310     h->Draw("c");
 311     PutTrendValue("TrackStatus", h->Integral());
 312     gPad->SetLogy(0);
 313     return kTRUE;
 314   case kFigTrackletStatus:
 315     if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
 316     leg = new TLegend(.68, .7, .97, .97);
 317     leg->SetBorderSize(0);leg->SetFillStyle(0);
 318     leg->SetHeader("TRD layer");
 319     for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 320       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
 321       if(kFIRST){
 322         h->Draw("pl");
 323         h->GetXaxis()->SetRangeUser(0.5, -1);
 324         h->GetYaxis()->CenterTitle();
 325         kFIRST = kFALSE;
 326       } else h->Draw("samepl");
 327       leg->AddEntry(h, Form("ly = %d", ily), "l");
 328       PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
 329     }
 330     leg->Draw();
 331     gPad->SetLogy(0);
 332     return kTRUE;
 333   case kFigChi2:
 334     MakePlotChi2();
 335     return kTRUE;
 336   case kFigPH:
 337     gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
 338     MakePlotPulseHeight();
 339     gPad->SetLogy(0);
 340     return kTRUE;
 341   case kFigChargeCluster:
 342     (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
 343     gPad->SetLogy(1);
 344     PutTrendValue("ChargeCluster", h->GetMaximumBin());
 345     PutTrendValue("ChargeClusterRMS", h->GetRMS());
 346     return kTRUE;
 347   case kFigChargeTracklet:
 348     (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
 349     PutTrendValue("ChargeTracklet", h->GetMaximumBin());
 350     PutTrendValue("ChargeTrackletRMS", h->GetRMS());
 351     return kTRUE;
 352   case kFigNeventsTrigger:
 353     ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
 354     return kTRUE;
 355   case kFigNeventsTriggerTracks:
 356     ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
 357     return kTRUE;
 358   case kFigTriggerPurity:
 359     if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
 360     break;
 361   default:
 362     break;
 363   }
 364   AliInfo(Form("Reference plot [%d] missing result", ifig));
 365   return kFALSE;
 366 }
 367
 368 //_______________________________________________________
 369 TObjArray *AliTRDcheckDET::Histos(){
 370   //
 371   // Create QA histograms
 372   //
 373
 374   if(fContainer) return fContainer;
 375
 376   fContainer = new TObjArray(20);
 377   //fContainer->SetOwner(kTRUE);
 378
 379   // Register Histograms
 380   TH1 * h = NULL;
 381   TAxis *ax = NULL;
 382   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
 383     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
 384     h->GetXaxis()->SetTitle("N_{clusters}");
 385     h->GetYaxis()->SetTitle("Entries");
 386   } else h->Reset();
 387   fContainer->AddAt(h, kNclustersTrack);
 388
 389   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
 390     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
 391     h->GetXaxis()->SetTitle("N_{clusters}");
 392     h->GetYaxis()->SetTitle("Entries");
 393   } else h->Reset();
 394   fContainer->AddAt(h, kNclustersTracklet);
 395
 396   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
 397     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
 398     h->GetXaxis()->SetTitle("N^{tracklet}");
 399     h->GetYaxis()->SetTitle("freq. [%]");
 400   } else h->Reset();
 401   fContainer->AddAt(h, kNtrackletsTrack);
 402
 403   if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
 404     h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
 405     h->GetXaxis()->SetTitle("N^{tracklet}");
 406     h->GetYaxis()->SetTitle("freq. [%]");
 407   }
 408   fContainer->AddAt(h, kNtrackletsSTA);
 409
 410   // Binning for momentum dependent tracklet Plots
 411   const Int_t kNp(14);
 412   Float_t P=0.1;
 413   Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1];
 414   for(Int_t i=0;i<kNp+1; i++,P=TMath::Exp(i*.15)-1.) binsP[i]=P;
 415   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
 416   if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
 417     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
 418     h = new TH2F("hNtlsBAR",
 419     "N_{tracklets} / track;p [GeV/c];N^{tracklet};freq. [%]",
 420     kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt);
 421   }
 422   fContainer->AddAt(h, kNtrackletsBAR);
 423
 424   //
 425   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
 426     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
 427     h->GetXaxis()->SetTitle("n_{row cross}");
 428     h->GetYaxis()->SetTitle("freq. [%]");
 429   } else h->Reset();
 430   fContainer->AddAt(h, kNtrackletsCross);
 431
 432   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
 433     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
 434     h->GetXaxis()->SetTitle("r [a.u]");
 435     h->GetYaxis()->SetTitle("Entries");
 436   } else h->Reset();
 437   fContainer->AddAt(h, kNtrackletsFindable);
 438
 439   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
 440     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
 441     h->GetXaxis()->SetTitle("N_{tracks}");
 442     h->GetYaxis()->SetTitle("Entries");
 443   } else h->Reset();
 444   fContainer->AddAt(h, kNtracksEvent);
 445
 446   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
 447     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
 448     h->GetXaxis()->SetTitle("sector");
 449     h->GetYaxis()->SetTitle("freq. [%]");
 450   } else h->Reset();
 451   fContainer->AddAt(h, kNtracksSector);
 452
 453   if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
 454     const Int_t nerr = 7;
 455     h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
 456     const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
 457     ax = h->GetXaxis();
 458     for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 459     h->SetYTitle("Relative Error to Good [%]");
 460   }
 461   fContainer->AddAt(h, kTrackStatus);
 462
 463   TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
 464   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
 465   fContainer->AddAt(arr, kTrackletStatus);
 466   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 467     if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
 468       const Int_t nerr = 8;
 469       h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
 470       h->SetLineColor(ily+1);
 471       const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
 472       ax = h->GetXaxis();
 473       for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 474       h->SetYTitle("Relative Error to Good [%]");
 475     } else h->Reset();
 476     arr->AddAt(h, ily);
 477   }
 478
 479   // <PH> histos
 480   arr = new TObjArray(2);
 481   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
 482   fContainer->AddAt(arr, kPH);
 483   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
 484     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
 485     h->GetXaxis()->SetTitle("Time / 100ns");
 486     h->GetYaxis()->SetTitle("<PH> [a.u]");
 487   } else h->Reset();
 488   arr->AddAt(h, 0);
 489   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
 490     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
 491   else h->Reset();
 492   arr->AddAt(h, 1);
 493
 494   // Chi2 histos
 495   if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
 496     h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
 497     h->SetXTitle("ndf");
 498     h->SetYTitle("#chi^{2}/ndf");
 499     h->SetZTitle("entries");
 500   } else h->Reset();
 501   fContainer->AddAt(h, kChi2);
 502
 503   if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
 504     h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
 505     h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
 506     h->GetYaxis()->SetTitle("Entries");
 507   }else h->Reset();
 508   fContainer->AddAt(h, kChargeCluster);
 509
 510   if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
 511     h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
 512     h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
 513     h->GetYaxis()->SetTitle("Entries");
 514   }else h->Reset();
 515   fContainer->AddAt(h, kChargeTracklet);
 516
 517
 518   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 519     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 520   else h->Reset();
 521   printf("Histos Adding \n");
 522
 523   fContainer->AddAt(h, kNeventsTrigger);
 524
 525   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 526     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 527   else h->Reset();
 528   fContainer->AddAt(h, kNeventsTriggerTracks);
 529
 530   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
 531     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
 532     h->GetXaxis()->SetTitle("Trigger Cluster");
 533     h->GetYaxis()->SetTitle("freq.");
 534   } else h->Reset();
 535   fContainer->AddAt(h, kTriggerPurity);
 536
 537   return fContainer;
 538 }
 539
 540 /*
 541 * Plotting Functions
 542 */
 543
 544 //_______________________________________________________
 545 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
 546 {
 547 //
 548 // Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
 549 //   PROL - track prolongation failure
 550 //   PROP - track propagation failure
 551 //   AJST - crossing sectors failure
 552 //   SNP  - too large bending
 553 //   TINI - tracklet initialization failure
 554 //   UPDT - track position/covariance update failure
 555 //
 556 // Performance plot looks as below:
 557 //Begin_Html
 558 //<img src="TRD/trackStatus.gif">
 559 //End_Html
 560 //
 561   if(track) fkTrack = track;
 562   if(!fkTrack){
 563     AliWarning("No Track defined.");
 564     return NULL;
 565   }
 566   TH1 *h = NULL;
 567   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
 568     AliWarning("No Histogram defined.");
 569     return NULL;
 570   }
 571   h->Fill(fkTrack->GetStatusTRD());
 572   return h;
 573 }
 574
 575 //_______________________________________________________
 576 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
 577 {
 578 //
 579 // Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
 580 //   Geom   -
 581 //   Bound  - tracklet too close to chamber walls
 582 //   NoCl   - no clusters in the track roads
 583 //   NoAttach - fail to attach clusters
 584 //   NoClTr - fail to use clusters for fit
 585 //   NoFit  - tracklet fit failled
 586 //   Chi2   - chi2 tracklet-track over threshold
 587 //
 588 // Performance plot looks as below:
 589 //Begin_Html
 590 //<img src="TRD/trackletStatus.gif">
 591 //End_Html
 592 //
 593   if(track) fkTrack = track;
 594   if(!fkTrack){
 595     AliWarning("No Track defined.");
 596     return NULL;
 597   }
 598   TObjArray *arr =NULL;
 599   if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
 600     AliWarning("Histograms not defined.");
 601     return NULL;
 602   }
 603
 604   TH1 *h = NULL;
 605   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 606     if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
 607       AliWarning(Form("Missing histo for layer %d.", ily));
 608       continue;
 609     }
 610     h->Fill(fkTrack->GetStatusTRD(ily));
 611   }
 612   return h;
 613 }
 614
 615 //_______________________________________________________
 616 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
 617   //
 618   // Plot the mean number of clusters per tracklet
 619   //
 620   if(track) fkTrack = track;
 621   if(!fkTrack){
 622     AliWarning("No Track defined.");
 623     return NULL;
 624   }
 625   TH1 *h = NULL;
 626   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
 627     AliWarning("No Histogram defined.");
 628     return NULL;
 629   }
 630   AliTRDseedV1 *tracklet = NULL;
 631   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 632     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 633     h->Fill(tracklet->GetN2());
 634   }
 635   return h;
 636 }
 637
 638 //_______________________________________________________
 639 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
 640   //
 641   // Plot the number of clusters in one track
 642   //
 643   if(track) fkTrack = track;
 644   if(!fkTrack){
 645     AliWarning("No Track defined.");
 646     return NULL;
 647   }
 648   TH1 *h = NULL;
 649   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
 650     AliWarning("No Histogram defined.");
 651     return NULL;
 652   }
 653
 654   Int_t nclusters = 0;
 655   AliTRDseedV1 *tracklet = NULL;
 656   AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
 657   Double_t momentumRec = par->P();
 658   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 659     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 660     Int_t n(tracklet->GetN());
 661     nclusters += n;
 662     if(DebugLevel() > 2){
 663       Int_t crossing = Int_t(tracklet->IsRowCross());
 664       Int_t detector = tracklet->GetDetector();
 665       Float_t theta = TMath::ATan(tracklet->GetZref(1));
 666       Float_t phi = TMath::ATan(tracklet->GetYref(1));
 667       Float_t momentumMC = 0.;
 668       Int_t pdg = 0;
 669       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
 670       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
 671       if(fkMC){
 672         if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
 673         pdg = fkMC->GetPDG();
 674       }
 675       (*DebugStream()) << "NClustersTrack"
 676         << "Detector="  << detector
 677         << "crossing="  << crossing
 678         << "momentumMC="<< momentumMC
 679         << "momentumRec=" << momentumRec
 680         << "pdg="                               << pdg
 681         << "theta="                     << theta
 682         << "phi="                               << phi
 683         << "kinkIndex=" << kinkIndex
 684         << "TPCncls="           << nclsTPC
 685         << "TRDncls="   << n
 686         << "\n";
 687     }
 688   }
 689   h->Fill(nclusters);
 690   return h;
 691 }
 692
 693
 694 //_______________________________________________________
 695 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
 696   //
 697   // Plot the number of tracklets
 698   //
 699   if(track) fkTrack = track;
 700   if(!fkTrack){
 701     AliWarning("No Track defined.");
 702     return NULL;
 703   }
 704   TH1 *h = NULL, *hSta = NULL; TH2 *hBarrel = NULL;
 705   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
 706     AliWarning("No Histogram defined.");
 707     return NULL;
 708   }
 709   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
 710   h->Fill(nTracklets);
 711   if(!fkESD) return h;
 712   Int_t status = fkESD->GetStatus();
 713
 714 /*  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);*/
 715   Double_t p = 0.;
 716   Int_t method = -1;    // to distinguish between stand alone and full barrel tracks in the debugging
 717   if((status & AliESDtrack::kTRDin) != 0){
 718     method = 1;
 719     // Full Barrel Track: Save momentum dependence
 720     if(!(hBarrel = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsBAR)))){
 721       AliWarning("Method: Barrel.  Histogram not processed!");
 722       return NULL;
 723     }
 724     AliExternalTrackParam *par(fkTrack->GetTrackIn());
 725     if(!par){
 726       AliError("Input track params missing");
 727       return NULL;
 728     }
 729     hBarrel->Fill(par->P(), nTracklets);
 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       return NULL;
 736     }
 737     hSta->Fill(nTracklets);
 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(0); 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
1264
1265
1266
1267   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1268   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1269   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1270     g[itl] = new TGraphErrors();
1271     g[itl]->SetLineColor(color[itl]);
1272     g[itl]->SetMarkerColor(color[itl]);
1273     g[itl]->SetMarkerStyle(20 + itl);
1274   }
1275
1276   TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
1277   TAxis *ax(hBar->GetXaxis());
1278   Int_t np(ax->GetNbins());
1279   for(Int_t ipBin(1); ipBin<np; ipBin++){
1280     h = hBar->ProjectionY("npBin", ipBin, ipBin);
1281     h->Scale(100./h->Integral());
1282     Float_t p(ax->GetBinCenter(ipBin));
1283     Float_t dp(ax->GetBinWidth(ipBin));
1284     for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1285       g[itl]->SetPoint(ipBin-1, p, h->GetBinContent(itl+1));
1286       g[itl]->SetPointError(ipBin-1, dp/3.46, h->GetBinError(itl+1));
1287     }
1288   }
1289
1290   TLegend *leg = new TLegend(0.5, 0.75, 0.9, 0.9);
1291   leg->SetBorderSize(0);
1292   leg->SetHeader("Tracklet/Track");
1293   leg->SetFillStyle(0);
1294   h = hBar->ProjectionX("npxBin"); h->Reset();
1295   h->GetYaxis()->SetRangeUser(0., 90.);
1296   h->GetXaxis()->SetMoreLogLabels();
1297   h->SetYTitle("Prob. [%]");
1298   h->Draw("p");
1299   for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1300     g[itl]->Draw("pc");
1301     leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
1302   }
1303
1304   leg->Draw();
1305   gPad->SetLogx();
1306 }
1307
1308 //________________________________________________________
1309 TH1* AliTRDcheckDET::MakePlotPulseHeight(){
1310   //
1311   // Create Plot of the Pluse Height Spectrum
1312   //
1313   TH1 *h, *h1, *h2;
1314   TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1315   h = (TH1F*)arr->At(0);
1316   h->SetMarkerStyle(24);
1317   h->SetMarkerColor(kBlack);
1318   h->SetLineColor(kBlack);
1319   h->Draw("e1");
1320   // Trending for the pulse height: plateau value, slope and timebin of the maximum
1321   TLinearFitter fit(1,"pol1");
1322   Double_t time = 0.;
1323   for(Int_t itime = 10; itime <= 20; itime++){
1324     time = static_cast<Double_t>(itime);
1325     fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
1326   }
1327   fit.Eval();
1328   Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1329   Double_t slope = fit.GetParameter(1);
1330   PutTrendValue("PHplateau", plateau);
1331   PutTrendValue("PHslope", slope);
1332   PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1333   AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1334 //   copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1335   h1 = (TH1F *)arr->At(1);
1336   h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1337   for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
1338     h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1339   h2->SetMarkerStyle(22);
1340   h2->SetMarkerColor(kBlue);
1341   h2->SetLineColor(kBlue);
1342   h2->Draw("e1same");
1343   gPad->Update();
1344 //   create axis according to the histogram dimensions of the original second histogram
1345   TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1346                     gPad->GetUymax(),
1347                     gPad->GetUxmax(),
1348                     gPad->GetUymax(),
1349                     -0.08, 4.88, 510,"-L");
1350   axis->SetLineColor(kBlue);
1351   axis->SetLabelColor(kBlue);
1352   axis->SetTextColor(kBlue);
1353   axis->SetTitle("x_{0}-x_{c} [cm]");
1354   axis->Draw();
1355   return h1;
1356 }
1357
1358 //________________________________________________________
1359 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1360   //
1361   // Draw nice bar plots
1362   //
1363   if(!histo->GetEntries()) return kFALSE;
1364   histo->Scale(100./histo->Integral());
1365   histo->SetFillColor(color);
1366   histo->SetBarOffset(.2);
1367   histo->SetBarWidth(.6);
1368   histo->Draw("bar1");
1369   return kTRUE;
1370 }