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