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