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