PWGPP/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 <TArrayD.h>
  27 #include <TAxis.h>
  28 #include <TCanvas.h>
  29 #include <TFile.h>
  30 #include <TH1F.h>
  31 #include <TH1I.h>
  32 #include <TH2F.h>
  33 #include <TH3S.h>
  34 #include <TH3F.h>
  35 #include <TF1.h>
  36 #include <TGaxis.h>
  37 #include <TGraph.h>
  38 #include <TGraphErrors.h>
  39 #include <TLegend.h>
  40 #include <TLinearFitter.h>
  41 #include <TMath.h>
  42 #include <TMap.h>
  43 #include <TProfile2D.h>
  44 #include <TObjArray.h>
  45 #include <TObject.h>
  46 #include <TObjString.h>
  47
  48 #include <TPad.h>
  49 #include <TProfile.h>
  50 #include <TProfile2D.h>
  51 #include <TROOT.h>
  52 #include <TChain.h>
  53
  54 #include "AliLog.h"
  55 #include "AliTRDcluster.h"
  56 #include "AliESDHeader.h"
  57 #include "AliESDRun.h"
  58 #include "AliESDtrack.h"
  59 #include "AliExternalTrackParam.h"
  60 #include "AliTRDgeometry.h"
  61 #include "AliTRDpadPlane.h"
  62 #include "AliTRDSimParam.h"
  63 #include "AliTRDseedV1.h"
  64 #include "AliTRDtrackV1.h"
  65 #include "AliTRDtrackerV1.h"
  66 #include "AliTRDReconstructor.h"
  67 #include "AliTrackReference.h"
  68 #include "AliTrackPointArray.h"
  69 #include "AliTracker.h"
  70 #include "TTreeStream.h"
  71
  72 #include "info/AliTRDtrackInfo.h"
  73 #include "info/AliTRDeventInfo.h"
  74 #include "AliTRDinfoGen.h"
  75 #include "AliTRDcheckDET.h"
  76 #include "AliTRDpwgppHelper.h"
  77
  78 #include <cstdio>
  79 #include <iostream>
  80
  81 ClassImp(AliTRDcheckDET)
  82
  83 const Color_t AliTRDcheckDET::fkColorsCentrality[AliTRDeventInfo::kCentralityClasses] = {kTeal, kOrange, kGreen, kBlue ,kRed};
  84
  85 //_______________________________________________________
  86 AliTRDcheckDET::AliTRDcheckDET():
  87   AliTRDrecoTask()
  88   ,fCentralityClass(-1)
  89   ,fTriggerNames(NULL)
  90   ,fFlags(0)
  91 {
  92   //
  93   // Default constructor
  94   //
  95   SetNameTitle("TRDcheckDET", "Basic TRD data checker");
  96 }
  97
  98 //_______________________________________________________
  99 AliTRDcheckDET::AliTRDcheckDET(char* name):
 100   AliTRDrecoTask(name, "Basic TRD data checker")
 101   ,fCentralityClass(-1)
 102   ,fTriggerNames(NULL)
 103   ,fFlags(0)
 104 {
 105   //
 106   // Default constructor
 107   //
 108   InitFunctorList();
 109 }
 110
 111
 112 //_______________________________________________________
 113 AliTRDcheckDET::~AliTRDcheckDET(){
 114   //
 115   // Destructor
 116   //
 117   if(fTriggerNames) delete fTriggerNames;
 118 }
 119
 120 //_______________________________________________________
 121 void AliTRDcheckDET::UserCreateOutputObjects(){
 122   //
 123   // Create Output Objects
 124   //
 125   AliTRDrecoTask::UserCreateOutputObjects();
 126   if(!fTriggerNames) fTriggerNames = new TMap();
 127 }
 128
 129 //_______________________________________________________
 130 void AliTRDcheckDET::UserExec(Option_t *opt){
 131   //
 132   // Execution function
 133   // Filling TRD quality histos
 134   //
 135   AliDebug(2, Form("EventInfo[%p] Header[%p]", (void*)fEvent, (void*)(fEvent?fEvent->GetEventHeader():NULL)));
 136   if(fEvent) fCentralityClass = fEvent->GetCentrality();
 137   else fCentralityClass = -1;  // Assume pp
 138
 139   AliTRDrecoTask::UserExec(opt);
 140
 141   TH1F *histo(NULL); AliTRDtrackInfo *fTrackInfo(NULL); Int_t nTracks(0);               // Count the number of tracks per event
 142   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 143     if(!fTracks->UncheckedAt(iti)) continue;
 144     if(!(fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti)))) continue;
 145     if(!fTrackInfo->GetTrack()) continue;
 146     nTracks++;
 147   }
 148   if(nTracks)
 149     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent)))) histo->Fill(nTracks);
 150
 151   if(!fEvent->GetEventHeader()) return; // For trigger statistics event header is essential
 152   Int_t triggermask = fEvent->GetEventHeader()->GetTriggerMask();
 153   TString triggername =  fEvent->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 154   AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
 155   if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger)))) histo->Fill(triggermask);
 156
 157   if(nTracks){
 158     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks)))) histo->Fill(triggermask);
 159   }
 160   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 161     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 162     // also set the label for both histograms
 163     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))))
 164       histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 165     if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))))
 166       histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 167   }
 168 }
 169
 170
 171 //_______________________________________________________
 172 //_______________________________________________________
 173 Bool_t AliTRDcheckDET::PostProcess(){
 174   //
 175   // Do Postprocessing (for the moment set the number of Reference histograms)
 176   //
 177
 178   TH1F *h(NULL), *h1(NULL);
 179
 180   // Calculate of the trigger clusters purity
 181   if((h  = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"))) &&
 182      (h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks")))) {
 183     h1->Divide(h);
 184     Float_t purities[20], val = 0; memset(purities, 0, 20*sizeof(Float_t));
 185     TString triggernames[20];
 186     Int_t nTriggerClasses = 0;
 187     for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
 188       if((val = h1->GetBinContent(ibin))){
 189         purities[nTriggerClasses] = val;
 190         triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
 191         nTriggerClasses++;
 192       }
 193     }
 194
 195     if((h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity)))){
 196       TAxis *ax = h->GetXaxis();
 197       for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
 198         h->Fill(itrg, purities[itrg]);
 199         ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
 200       }
 201       ax->SetRangeUser(-0.5, nTriggerClasses+.5);
 202       h->GetYaxis()->SetRangeUser(0,1);
 203     }
 204   }
 205
 206   // track status
 207   if((h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus)))){
 208     Float_t ok = h->GetBinContent(1);
 209     Int_t nerr = h->GetNbinsX();
 210     for(Int_t ierr=nerr; ierr--;){
 211       h->SetBinContent(ierr+1, ok>0.?1.e2*h->GetBinContent(ierr+1)/ok:0.);
 212     }
 213     h->SetBinContent(1, 0.);
 214   }
 215   // tracklet status
 216
 217   TObjArray *arr(NULL);
 218   if(( arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus)))){
 219     for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
 220       if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
 221       Float_t okB = h->GetBinContent(1);
 222       Int_t nerrB = h->GetNbinsX();
 223       for(Int_t ierr=nerrB; ierr--;){
 224         h->SetBinContent(ierr+1, okB>0.?1.e2*h->GetBinContent(ierr+1)/okB:0.);
 225       }
 226       h->SetBinContent(1, 0.);
 227     }
 228   }
 229   fNRefFigures = 17;
 230
 231   return kTRUE;
 232 }
 233
 234 void AliTRDcheckDET::MakeSummary(){
 235   //
 236   // Create summary plots for TRD check DET
 237   // This function creates 2 summary plots:
 238   // - General Quantities
 239   // - PHS
 240   // The function will reuse GetRefFigure
 241   //
 242
 243   TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);
 244   cOut->Divide(3,3);
 245
 246   // Create figures using GetRefFigure
 247   cOut->cd(1); GetRefFigure(kFigNtracksEvent);
 248   cOut->cd(2); GetRefFigure(kFigNtracksSector);
 249   cOut->cd(3); GetRefFigure(kFigNclustersTrack);
 250   cOut->cd(4); GetRefFigure(kFigNclustersTracklet);
 251   cOut->cd(5); GetRefFigure(kFigNtrackletsTrack);
 252   cOut->cd(6); GetRefFigure(kFigNTrackletsP);
 253   cOut->cd(7); GetRefFigure(kFigChargeCluster);
 254   cOut->cd(8); GetRefFigure(kFigChargeTracklet);
 255   cOut->SaveAs("TRD_TrackQuality.gif");
 256   delete cOut;
 257
 258   // Second Plot: PHS
 259   cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 512);
 260   cOut->cd(); GetRefFigure(kFigPH);
 261   cOut->SaveAs("TRD_PH.gif");
 262   delete cOut;
 263
 264   // Third Plot: Mean Number of clusters as function of eta, phi and layer
 265    cOut = new TCanvas(Form("summary%s3", GetName()), Form("Summary 3 for task %s", GetName()), 1024, 768);
 266   cOut->cd(); MakePlotMeanClustersLayer();
 267   cOut->SaveAs("TRD_MeanNclusters.gif");
 268   delete cOut;
 269
 270 }
 271
 272 //_______________________________________________________
 273 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
 274         //
 275         // Setting Reference Figures
 276         //
 277         gPad->SetLogy(0);
 278         gPad->SetLogx(0);
 279         TH1 *h = NULL; TObjArray *arr=NULL;
 280         TLegend *leg = NULL;
 281         Bool_t kFIRST(1);
 282         switch(ifig){
 283         case kFigNclustersTrack:
 284                 MakePlotNclustersTrack();
 285                 return kTRUE;
 286         case kFigNclustersTracklet:
 287                 MakePlotNclustersTracklet();
 288                 return kTRUE;
 289         case kFigNtrackletsTrack:
 290                 h=MakePlotNTracklets();
 291                 if(h){
 292                         PutTrendValue("NTrackletsTrack", h->GetMean());
 293                         PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
 294                 }
 295                 return kTRUE;
 296         case kFigNTrackletsP:
 297                 MakePlotnTrackletsVsP();
 298                 return kTRUE;
 299         case kFigNtrackletsCross:
 300                 h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsCross"), -1);
 301                 if(!MakeBarPlot(h, kRed)) break;
 302                 PutTrendValue("NTrackletsCross", h->GetMean());
 303                 PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
 304                 return kTRUE;
 305         case kFigNtrackletsFindable:
 306                 h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsFindable"), -1);
 307                 if(!MakeBarPlot(h, kGreen)) break;
 308                 PutTrendValue("NTrackletsFindable", h->GetMean());
 309                 PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
 310                 return kTRUE;
 311         case kFigNtracksEvent:
 312                 (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
 313                 PutTrendValue("NTracksEvent", h->GetMean());
 314                 PutTrendValue("NTracksEventRMS", h->GetRMS());
 315                 return kTRUE;
 316         case kFigNtracksSector:
 317                 h = (TH1F*)fContainer->FindObject("hNtrksSector");
 318                 if(!MakeBarPlot(h, kGreen)) break;
 319                 PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
 320                 return kTRUE;
 321         case kFigTrackStatus:
 322                 if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
 323                 h->GetXaxis()->SetRangeUser(0.5, -1);
 324                 h->GetYaxis()->CenterTitle();
 325                 h->Draw("c");
 326                 PutTrendValue("TrackStatus", h->Integral());
 327                 gPad->SetLogy(0);
 328                 return kTRUE;
 329         case kFigTrackletStatus:
 330                 if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
 331                 leg = new TLegend(.68, .7, .97, .97);
 332                 leg->SetBorderSize(0);leg->SetFillStyle(0);
 333                 leg->SetHeader("TRD layer");
 334                 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 335                         if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
 336                         if(kFIRST){
 337                                 h->Draw("pl");
 338                                 h->GetXaxis()->SetRangeUser(0.5, -1);
 339                                 h->GetYaxis()->CenterTitle();
 340                                 kFIRST = kFALSE;
 341                         } else h->Draw("samepl");
 342                         leg->AddEntry(h, Form("ly = %d", ily), "l");
 343                         PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
 344                 }
 345                 leg->Draw();
 346                 gPad->SetLogy(0);
 347                 return kTRUE;
 348         case kFigChi2:
 349                 MakePlotChi2();
 350                 return kTRUE;
 351         case kFigPH:
 352                 gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
 353                 MakePlotPulseHeight();
 354                 gPad->SetLogy(0);
 355                 return kTRUE;
 356         case kFigChargeCluster:
 357                 h = ProjectCentrality((TH2*)fContainer->FindObject("hQcl"), -1);
 358                 h->Draw("c");
 359                 gPad->SetLogy(1);
 360                 PutTrendValue("ChargeCluster", h->GetMaximumBin());
 361                 PutTrendValue("ChargeClusterRMS", h->GetRMS());
 362                 return kTRUE;
 363         case kFigChargeTracklet:
 364                 MakePlotTrackletCharge();
 365                 return kTRUE;
 366         case kFigNeventsTrigger:
 367                 ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
 368                 return kTRUE;
 369         case kFigNeventsTriggerTracks:
 370                 ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
 371                 return kTRUE;
 372         case kFigTriggerPurity:
 373                 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
 374                 break;
 375         default:
 376                 break;
 377         }
 378         AliInfo(Form("Reference plot [%d] missing result", ifig));
 379         return kFALSE;
 380 }
 381
 382 //_______________________________________________________
 383 TObjArray *AliTRDcheckDET::Histos(){
 384   //
 385   // Create QA histograms
 386   //
 387
 388   if(fContainer) return fContainer;
 389
 390   fContainer = new TObjArray(20);
 391   fContainer->SetOwner(kTRUE);
 392
 393   // Register Histograms
 394   TH1 * h = NULL;
 395   TH2 * h2 = NULL;      // Pointer for two dimensional histograms
 396   TH3 * h3 = NULL;      // Pointer for tree dimensional histograms
 397   TAxis *ax = NULL;
 398   if(!(h2 = (TH2F *)gROOT->FindObject("hNcls"))){
 399     h2 = new TH2F("hNcls", "N_{clusters}/track;N_{clusters};Centrality;Entries", 181, -0.5, 180.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 400   } else h2->Reset();
 401   fContainer->AddAt(h2, kNclustersTrack);
 402
 403   TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
 404   arr->SetOwner(kTRUE);  arr->SetName("clusters");
 405   fContainer->AddAt(arr, kNclustersLayer);
 406   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 407     if(!(h = (TProfile2D *)gROOT->FindObject(Form("hNcl%d", ily)))){
 408       h = new TProfile2D(Form("hNcl%d", ily), Form("Mean Number of clusters in Layer %d", ily), 100, -1.0, 1.0, 50, -1.1*TMath::Pi(), 1.1*TMath::Pi());
 409       h->GetXaxis()->SetTitle("#eta");
 410       h->GetYaxis()->SetTitle("#phi");
 411     } else h->Reset();
 412     arr->AddAt(h, ily);
 413   }
 414
 415   if(!(h2 = (TH2F *)gROOT->FindObject("hNclTls"))){
 416     h2 = new TH2F("hNclTls","N_{clusters}/tracklet;N_{clusters};Entries", 51, -0.5, 50.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 417   } else h2->Reset();
 418   fContainer->AddAt(h2, kNclustersTracklet);
 419
 420   if(!(h2 = (TH2F *)gROOT->FindObject("hNtls"))){
 421     h2 = new TH2F("hNtls", "N_{tracklets}/track;N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 422   } else h2->Reset();
 423   fContainer->AddAt(h2, kNtrackletsTrack);
 424
 425   if(!(h = (TH2F *)gROOT->FindObject("htlsSTA"))){
 426     h = new TH2F("hNtlsSTA", "#splitline{N_{tracklets}/track}{Stand Alone};N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 427   }
 428   fContainer->AddAt(h, kNtrackletsSTA);
 429
 430   // Binning for momentum dependent tracklet Plots
 431   const Int_t kNp(30);
 432   Float_t p=0.2;
 433   Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1], binsCent[AliTRDeventInfo::kCentralityClasses+2];
 434   for(Int_t i=0;i<kNp+1; i++,p+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=p;
 435   for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
 436   for(Int_t icent = -1; icent < AliTRDeventInfo::kCentralityClasses + 1; icent++) binsCent[icent+1] = icent - 0.5;
 437   if(!(h3 = (TH3F *)gROOT->FindObject("htlsBAR"))){
 438     // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
 439     h3 = new TH3F("hNtlsBAR",
 440     "N_{tracklets}/track;p [GeV/c];N^{tracklet};freq. [%]",
 441     kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt, AliTRDeventInfo::kCentralityClasses + 1, binsCent);
 442   }
 443   fContainer->AddAt(h3, kNtrackletsBAR);
 444
 445   //
 446   if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsCross"))){
 447     h2 = new TH2F("hNtlsCross", "N_{tracklets}^{cross}/track;n_{row cross};Centrality;freq.[%]", 7, -0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 448   } else h2->Reset();
 449   fContainer->AddAt(h2, kNtrackletsCross);
 450
 451   if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsFindable"))){
 452     h2 = new TH2F("hNtlsFindable", "Found/Findable Tracklets;r[a.u];Centrality;Entries" , 101, -0.005, 1.005,  AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 453   } else h2->Reset();
 454   fContainer->AddAt(h2, kNtrackletsFindable);
 455
 456   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
 457     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
 458     h->GetXaxis()->SetTitle("N_{tracks}");
 459     h->GetYaxis()->SetTitle("Entries");
 460   } else h->Reset();
 461   fContainer->AddAt(h, kNtracksEvent);
 462
 463   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
 464     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
 465     h->GetXaxis()->SetTitle("sector");
 466     h->GetYaxis()->SetTitle("freq. [%]");
 467   } else h->Reset();
 468   fContainer->AddAt(h, kNtracksSector);
 469
 470   if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
 471     const Int_t nerr = 7;
 472     h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
 473     const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
 474     ax = h->GetXaxis();
 475     for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 476     h->SetYTitle("Relative Error to Good [%]");
 477   }
 478   fContainer->AddAt(h, kTrackStatus);
 479
 480   arr = new TObjArray(AliTRDgeometry::kNlayer);
 481   arr->SetOwner(kTRUE);  arr->SetName("TrackletStatus");
 482   fContainer->AddAt(arr, kTrackletStatus);
 483   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 484     if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
 485       const Int_t nerr = 8;
 486       h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
 487       h->SetLineColor(ily+1);
 488       const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
 489       ax = h->GetXaxis();
 490       for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
 491       h->SetYTitle("Relative Error to Good [%]");
 492     } else h->Reset();
 493     arr->AddAt(h, ily);
 494   }
 495
 496   // <PH> histos
 497   arr = new TObjArray(3);
 498   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
 499   fContainer->AddAt(arr, kPH);
 500   if(!(h3 = (TH3F *)gROOT->FindObject("hPHx"))){
 501     h3 = new TH3F("hPHx", "<PH>(x);x[mm];Charge[a.u.];Centrality", 31, -0.08, 4.88, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 502   } else h3->Reset();
 503   arr->AddAt(h3, 0);
 504   if(!(h3 = (TH3F *)gROOT->FindObject("hPHt"))){
 505     h3 = new TH3F("hPHt", "<PH>(t);time[100ns];Charge[a.u.];Centrality", 31, -0.5, 30.5, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 506   } else h3->Reset();
 507   arr->AddAt(h3, 1);
 508
 509   // Chi2 histos
 510   if(!(h3 = (TH3S*)gROOT->FindObject("hChi2"))){
 511     h3 = new TH3S("hChi2", "#chi^{2}/track;ndf;#chi^{2}/ndf;Centrality", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 512   } else h3->Reset();
 513   fContainer->AddAt(h3, kChi2);
 514
 515   if(!(h2 = (TH2F *)gROOT->FindObject("hQcl"))){
 516     h2 = new TH2F("hQcl", "Q_{cluster};Charge[a.u.];Centrality;Entries", 200, 0, 1200, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 517   }else h2->Reset();
 518   fContainer->AddAt(h2, kChargeCluster);
 519
 520   if(!(h2 = (TH2F *)gROOT->FindObject("hQtrklt"))){
 521     h2 = new TH2F("hQtrklt", "Q_{tracklet};Charge[a.u.];Centrality;Entries", 6000, 0, 6000, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
 522   }else h2->Reset();
 523   fContainer->AddAt(h2, kChargeTracklet);
 524
 525
 526   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 527     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 528   else h->Reset();
 529
 530   fContainer->AddAt(h, kNeventsTrigger);
 531
 532   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 533     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 534   else h->Reset();
 535   fContainer->AddAt(h, kNeventsTriggerTracks);
 536
 537   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
 538     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
 539     h->GetXaxis()->SetTitle("Trigger Cluster");
 540     h->GetYaxis()->SetTitle("freq.");
 541   } else h->Reset();
 542   fContainer->AddAt(h, kTriggerPurity);
 543
 544   return fContainer;
 545 }
 546
 547 //_______________________________________________________
 548 TH1 *AliTRDcheckDET::ProjectCentrality(TH2 *hIn, Int_t centralityBin){
 549   //
 550   // Project histogram to a given centrality Bin
 551   //
 552   if(!hIn) return NULL;
 553   if(centralityBin >= AliTRDeventInfo::kCentralityClasses) centralityBin = -1;
 554   Int_t binMin = centralityBin > -1 ? centralityBin + 2 : 0,
 555         binMax = centralityBin > -1 ? centralityBin + 2 : -1;
 556   return hIn->ProjectionX(Form("%s_%d", hIn->GetName(), centralityBin), binMin, binMax);
 557 }
 558
 559 /*
 560 * Plotting Functions
 561 */
 562
 563 //_______________________________________________________
 564 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
 565 {
 566 //
 567 // Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
 568 //   PROL - track prolongation failure
 569 //   PROP - track propagation failure
 570 //   AJST - crossing sectors failure
 571 //   SNP  - too large bending
 572 //   TINI - tracklet initialization failure
 573 //   UPDT - track position/covariance update failure
 574 //
 575 // Performance plot looks as below:
 576 //Begin_Html
 577 //<img src="TRD/trackStatus.gif">
 578 //End_Html
 579 //
 580   if(track) fkTrack = track;
 581   if(!fkTrack){
 582     AliDebug(4, "No Track defined.");
 583     return NULL;
 584   }
 585   TH1 *h = NULL;
 586   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
 587     AliWarning("No Histogram defined.");
 588     return NULL;
 589   }
 590   h->Fill(fkTrack->GetStatusTRD());
 591   return h;
 592 }
 593
 594 //_______________________________________________________
 595 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
 596 {
 597 //
 598 // Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
 599 //   Geom   -
 600 //   Bound  - tracklet too close to chamber walls
 601 //   NoCl   - no clusters in the track roads
 602 //   NoAttach - fail to attach clusters
 603 //   NoClTr - fail to use clusters for fit
 604 //   NoFit  - tracklet fit failled
 605 //   Chi2   - chi2 tracklet-track over threshold
 606 //
 607 // Performance plot looks as below:
 608 //Begin_Html
 609 //<img src="TRD/trackletStatus.gif">
 610 //End_Html
 611 //
 612   if(track) fkTrack = track;
 613   if(!fkTrack){
 614     AliDebug(4, "No Track defined.");
 615     return NULL;
 616   }
 617   TObjArray *arr =NULL;
 618   if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
 619     AliWarning("Histograms not defined.");
 620     return NULL;
 621   }
 622
 623   TH1 *h = NULL;
 624   for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
 625     if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
 626       AliWarning(Form("Missing histo for layer %d.", ily));
 627       continue;
 628     }
 629     h->Fill(fkTrack->GetStatusTRD(ily));
 630   }
 631   return h;
 632 }
 633
 634 //_______________________________________________________
 635 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
 636   //
 637   // Plot the mean number of clusters per tracklet
 638   //
 639   if(track) fkTrack = track;
 640   if(!fkTrack){
 641     AliDebug(4, "No Track defined.");
 642     return NULL;
 643   }
 644   AliExternalTrackParam *par = fkTrack->GetTrackIn() ? fkTrack->GetTrackIn() : fkTrack->GetTrackOut();
 645   TH2 *h = NULL;
 646   TProfile2D *hlayer = NULL;
 647   Double_t eta = 0., phi = 0.;
 648   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTracklet)))){
 649     AliWarning("No Histogram defined.");
 650     return NULL;
 651   }
 652   AliTRDseedV1 *tracklet = NULL;
 653   TObjArray *histosLayer = dynamic_cast<TObjArray *>(fContainer->At(kNclustersLayer));
 654   if(!histosLayer){
 655     AliWarning("No Histograms for single layer defined");
 656   }
 657   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 658     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 659     h->Fill(tracklet->GetN2(), fCentralityClass);
 660     if(histosLayer && par){
 661       if((hlayer = dynamic_cast<TProfile2D *>(histosLayer->At(itl)))){
 662         GetEtaPhiAt(par, tracklet->GetX0(), eta, phi);
 663         hlayer->Fill(eta, phi, tracklet->GetN2());
 664       }
 665     }
 666   }
 667   return h;
 668 }
 669
 670 //_______________________________________________________
 671 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
 672   //
 673   // Plot the number of clusters in one track
 674   //
 675   if(track) fkTrack = track;
 676   if(!fkTrack){
 677     AliDebug(4, "No Track defined.");
 678     return NULL;
 679   }
 680   TH2 *h = NULL;
 681   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTrack)))){
 682     AliWarning("No Histogram defined.");
 683     return NULL;
 684   }
 685
 686   Int_t nclusters = 0;
 687   AliTRDseedV1 *tracklet = NULL;
 688   AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
 689   if(!par) return NULL;
 690   Double_t momentumRec = par->P();
 691   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 692     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 693     Int_t n(tracklet->GetN());
 694     nclusters += n;
 695     if(DebugLevel() > 2){
 696       Int_t crossing = Int_t(tracklet->IsRowCross());
 697       Int_t detector = tracklet->GetDetector();
 698       Float_t theta = TMath::ATan(tracklet->GetZref(1));
 699       Float_t phi = TMath::ATan(tracklet->GetYref(1));
 700       Float_t momentumMC = 0.;
 701       Int_t pdg = 0;
 702       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
 703       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
 704       if(fkMC){
 705         if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
 706         pdg = fkMC->GetPDG();
 707       }
 708       (*DebugStream()) << "NClustersTrack"
 709         << "Detector="  << detector
 710         << "Centrality="<< fCentralityClass
 711         << "crossing="  << crossing
 712         << "momentumMC="<< momentumMC
 713         << "momentumRec="<< momentumRec
 714         << "pdg="                               << pdg
 715         << "theta="                     << theta
 716         << "phi="                               << phi
 717         << "kinkIndex=" << kinkIndex
 718         << "TPCncls="           << nclsTPC
 719         << "TRDncls="   << n
 720         << "\n";
 721     }
 722   }
 723   h->Fill(nclusters, fCentralityClass);
 724   return h;
 725 }
 726
 727
 728 //_______________________________________________________
 729 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
 730   //
 731   // Plot the number of tracklets
 732   //
 733   if(track) fkTrack = track;
 734   if(!fkTrack){
 735     AliDebug(4, "No Track defined.");
 736     return NULL;
 737   }
 738   TH2 *h = NULL, *hSta = NULL; TH3 *hBarrel = NULL;
 739   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsTrack)))){
 740     AliWarning("No Histogram defined.");
 741     return NULL;
 742   }
 743   Int_t nTracklets = fkTrack->GetNumberOfTracklets();
 744   h->Fill(nTracklets, fCentralityClass);
 745   if(!fkESD) return h;
 746   Int_t status = fkESD->GetStatus();
 747
 748 /*  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);*/
 749   Double_t p = 0.;
 750   Int_t method = -1;    // to distinguish between stand alone and full barrel tracks in the debugging
 751   if((status & AliESDtrack::kTRDin) != 0){
 752     method = 1;
 753     // Full Barrel Track: Save momentum dependence
 754     if(!(hBarrel = dynamic_cast<TH3F *>(fContainer->At(kNtrackletsBAR)))){
 755       AliWarning("Method: Barrel.  Histogram not processed!");
 756       return NULL;
 757     }
 758     AliExternalTrackParam *par(fkTrack->GetTrackIn());
 759     if(!par){
 760       AliError("Input track params missing");
 761       return NULL;
 762     }
 763     p = par->P(); // p needed later in the debug streaming
 764     hBarrel->Fill(p, nTracklets, fCentralityClass);
 765   } else {
 766     // Stand alone Track: momentum dependence not usefull
 767     method = 0;
 768     if(!(hSta = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsSTA)))) {
 769       AliWarning("Method: StandAlone.  Histogram not processed!");
 770       return NULL;
 771     }
 772     hSta->Fill(nTracklets, fCentralityClass);
 773   }
 774
 775   if(DebugLevel() > 2){
 776     AliTRDseedV1 *tracklet = NULL;
 777     AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 778     Int_t sector = -1, stack = -1, detector;
 779     for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 780       if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
 781       detector = tracklet->GetDetector();
 782       sector = geo->GetSector(detector);
 783       stack = geo->GetStack(detector);
 784       break;
 785     }
 786     (*DebugStream()) << "NTrackletsTrack"
 787       << "Sector="      << sector
 788       << "Stack="       << stack
 789       << "Centrality="  << fCentralityClass
 790       << "NTracklets="  << nTracklets
 791       << "Method="      << method
 792       << "p="           << p
 793       << "\n";
 794   }
 795   if(DebugLevel() > 3){
 796     AliTRDseedV1 *tracklet = NULL;
 797     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 798       if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
 799         (*DebugStream()) << "NTrackletsLayer"
 800         << "Layer=" << il
 801         << "p=" << p
 802         << "\n";
 803       }
 804     }
 805   }
 806   return h;
 807 }
 808
 809
 810 //_______________________________________________________
 811 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
 812   //
 813   // Plot the number of tracklets
 814   //
 815   if(track) fkTrack = track;
 816   if(!fkTrack){
 817     AliDebug(4, "No Track defined.");
 818     return NULL;
 819   }
 820   TH2 *h = NULL;
 821   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsCross)))){
 822     AliWarning("No Histogram defined.");
 823     return NULL;
 824   }
 825
 826   Int_t ncross = 0;
 827   AliTRDseedV1 *tracklet = NULL;
 828   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 829     if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
 830
 831     if(tracklet->IsRowCross()) ncross++;
 832   }
 833   h->Fill(ncross, fCentralityClass);
 834   return h;
 835 }
 836
 837 //_______________________________________________________
 838 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
 839   //
 840   // Plots the ratio of number of tracklets vs.
 841   // number of findable tracklets
 842   //
 843   // Findable tracklets are defined as track prolongation
 844   // to layer i does not hit the dead area +- epsilon
 845   //
 846   // In order to check whether tracklet hist active area in Layer i,
 847   // the track is refitted and the fitted position + an uncertainty
 848   // range is compared to the chamber border (also with a different
 849   // uncertainty)
 850   //
 851   // For the track fit two cases are distinguished:
 852   // If the track is a stand alone track (defined by the status bit
 853   // encoding, then the track is fitted with the tilted Rieman model
 854   // Otherwise the track is fitted with the Kalman fitter in two steps:
 855   // Since the track parameters are give at the outer point, we first
 856   // fit in direction inwards. Afterwards we fit again in direction outwards
 857   // to extrapolate the track to layers which are not reached by the track
 858   // For the Kalman model, the radial track points have to be shifted by
 859   // a distance epsilon in the direction that we want to fit
 860   //
 861   const Float_t epsilon = 0.01;   // dead area tolerance
 862   const Float_t epsilonR = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
 863   const Float_t deltaY = 0.7;    // Tolerance in the track position in y-direction
 864   const Float_t deltaZ = 7.0;    // Tolerance in the track position in z-direction (Padlength)
 865   Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 866
 867   if(track) fkTrack = track;
 868   if(!fkTrack){
 869     AliDebug(4, "No Track defined.");
 870     return NULL;
 871   }
 872   TH2 *h = NULL;
 873   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsFindable)))){
 874     AliWarning("No Histogram defined.");
 875     return NULL;
 876   }
 877   Int_t nFound = 0, nFindable = 0;
 878   Int_t stack = -1;
 879   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 880   Double_t y = 0., z = 0.;
 881   AliTRDseedV1 *tracklet = NULL;
 882   AliTRDpadPlane *pp;
 883   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 884     if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
 885       tracklet->SetReconstructor(AliTRDinfoGen::Reconstructor());
 886       nFound++;
 887     }
 888   }
 889   // 2 Different cases:
 890   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 891   // 2nd barrel track: here we propagate the track to the layers
 892   AliTrackPoint points[6];
 893   Float_t xyz[3];
 894   memset(xyz, 0, sizeof(Float_t) * 3);
 895   if(fkESD && ((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 896     // stand alone track
 897     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 898       xyz[0] = xAnode[il];
 899       points[il].SetXYZ(xyz);
 900     }
 901     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
 902   } else {
 903     // barrel track
 904     //
 905     // 2 Steps:
 906     // -> Kalman inwards
 907     // -> Kalman outwards
 908     AliTRDtrackV1 copyTrack(*fkTrack);  // Do Kalman on a (non-constant) copy of the track
 909     AliTrackPoint pointsInward[6], pointsOutward[6];
 910     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 911       // In order to avoid complications in the Kalman filter if the track points have the same radial
 912       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
 913       // in the direction we want to go
 914       // The track points have to be in reverse order for the Kalman Filter inwards
 915       xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
 916       pointsInward[il].SetXYZ(xyz);
 917       xyz[0] = xAnode[il] + epsilonR;
 918       pointsOutward[il].SetXYZ(xyz);
 919     }
 920     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 921       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
 922     // Kalman inwards
 923     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
 924     memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
 925     // Kalman outwards
 926     AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
 927     memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 928   }
 929   AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
 930   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 931     y = points[il].GetY();
 932     z = points[il].GetZ();
 933     if((stack = geo->GetStack(z, il)) < 0) continue; // Not findable
 934     pp = geo->GetPadPlane(il, stack);
 935     ymin = pp->GetCol0() + epsilon;
 936     ymax = pp->GetColEnd() - epsilon;
 937     zmin = pp->GetRowEnd() + epsilon;
 938     zmax = pp->GetRow0() - epsilon;
 939     // ignore y-crossing (material)
 940     if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
 941       if(DebugLevel() > 3){
 942         Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
 943         Int_t hasTracklet = tracklet ? 1 : 0;
 944         (*DebugStream())   << "FindableTracklets"
 945           << "layer="     << il
 946           << "ytracklet=" << posTracklet[0]
 947           << "ytrack="    << y
 948           << "ztracklet=" << posTracklet[1]
 949           << "ztrack="    << z
 950           << "tracklet="  << hasTracklet
 951           << "\n";
 952       }
 953   }
 954
 955   h->Fill((nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1), fCentralityClass);
 956   AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 957   return h;
 958 }
 959
 960
 961 //_______________________________________________________
 962 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
 963   //
 964   // Plot the chi2 of the track
 965   //
 966   if(track) fkTrack = track;
 967   if(!fkTrack){
 968     AliDebug(4, "No Track defined.");
 969     return NULL;
 970   }
 971   TH3 *h = NULL;
 972   if(!(h = dynamic_cast<TH3S*>(fContainer->At(kChi2)))) {
 973     AliWarning("No Histogram defined.");
 974     return NULL;
 975   }
 976   Int_t n = fkTrack->GetNumberOfTracklets();
 977   if(!n) return NULL;
 978
 979   h->Fill(n, fkTrack->GetChi2()/n, fCentralityClass);
 980   return h;
 981 }
 982
 983
 984 //_______________________________________________________
 985 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
 986   //
 987   // Plot the average pulse height
 988   //
 989   if(track) fkTrack = track;
 990   if(!fkTrack){
 991     AliDebug(4, "No Track defined.");
 992     return NULL;
 993   }
 994   TH3F *h = NULL;
 995   if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
 996     AliWarning("No Histogram defined.");
 997     return NULL;
 998   }
 999   AliTRDseedV1 *tracklet = NULL;
1000   AliTRDcluster *c = NULL;
1001   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1002     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1003     Int_t crossing = Int_t(tracklet->IsRowCross());
1004     Int_t detector = tracklet->GetDetector();
1005     tracklet->ResetClusterIter();
1006     while((c = tracklet->NextCluster())){
1007       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1008       Int_t localtime        = c->GetLocalTimeBin();
1009       Double_t absoluteCharge = TMath::Abs(c->GetQ());
1010       h->Fill(localtime, absoluteCharge, fCentralityClass);
1011       if(DebugLevel() > 3){
1012         Int_t inChamber = c->IsInChamber() ? 1 : 0;
1013         Double_t distance[2];
1014         GetDistanceToTracklet(distance, tracklet, c);
1015         Float_t theta = TMath::ATan(tracklet->GetZref(1));
1016         Float_t phi = TMath::ATan(tracklet->GetYref(1));
1017         AliExternalTrackParam *trdPar = fkTrack->GetTrackIn();
1018         Float_t momentumMC = 0, momentumRec = trdPar ? trdPar->P() : fkTrack->P(); // prefer Track Low
1019         Int_t pdg = 0;
1020         Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1021         UShort_t tpcNCLS = fkESD ? fkESD->GetTPCncls() : 0;
1022         if(fkMC){
1023           if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
1024           pdg = fkMC->GetPDG();
1025         }
1026         (*DebugStream()) << "PHt"
1027           << "Detector="        << detector
1028           << "Centrality="<< fCentralityClass
1029           << "crossing="        << crossing
1030           << "inChamber=" << inChamber
1031           << "Timebin="         << localtime
1032           << "Charge="          << absoluteCharge
1033           << "momentumMC="      << momentumMC
1034           << "momentumRec="     << momentumRec
1035           << "pdg="                             << pdg
1036           << "theta="                   << theta
1037           << "phi="                             << phi
1038           << "kinkIndex="       << kinkIndex
1039           << "TPCncls="         << tpcNCLS
1040           << "dy="        << distance[0]
1041           << "dz="        << distance[1]
1042           << "c.="        << c
1043           << "\n";
1044       }
1045     }
1046   }
1047   return h;
1048 }
1049
1050 //_______________________________________________________
1051 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1052   //
1053   // Plots the average pulse height vs the distance from the anode wire
1054   // (plus const anode wire offset)
1055   //
1056   if(track) fkTrack = track;
1057   if(!fkTrack){
1058      AliDebug(4, "No Track defined.");
1059      return NULL;
1060   }
1061   TH3 *h = NULL;
1062   if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
1063     AliWarning("No Histogram defined.");
1064     return NULL;
1065   }
1066   AliTRDseedV1 *tracklet(NULL);
1067   AliTRDcluster *c(NULL);
1068   Double_t xd(0.), dqdl(0.);
1069   TVectorD vq(AliTRDseedV1::kNtb), vxd(AliTRDseedV1::kNtb), vdqdl(AliTRDseedV1::kNtb);
1070   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1071     if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1072     Int_t det(tracklet->GetDetector());
1073     Bool_t rc(tracklet->IsRowCross());
1074     for(Int_t ic(0); ic<AliTRDseedV1::kNtb; ic++){
1075       Bool_t kFIRST(kFALSE);
1076       if(!(c = tracklet->GetClusters(ic))){
1077          if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1078       } else kFIRST=kTRUE;
1079       if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1080       xd = tracklet->GetX0() - c->GetX(); vxd[ic] = xd;
1081       dqdl=tracklet->GetdQdl(ic); vdqdl[ic] = dqdl;
1082       vq[ic]=c->GetQ();
1083       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) vq[ic]+=c->GetQ();
1084       h->Fill(xd, dqdl, fCentralityClass);
1085     }
1086     if(DebugLevel() > 3){
1087       (*DebugStream()) << "PHx"
1088         << "det="  << det
1089         << "rc="   << rc
1090         << "xd="   << &vxd
1091         << "q="    << &vq
1092         << "dqdl=" << &vdqdl
1093         << "\n";
1094     }
1095   }
1096   return h;
1097 }
1098
1099 //_______________________________________________________
1100 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1101   //
1102   // Plot the cluster charge
1103   //
1104   if(track) fkTrack = track;
1105   if(!fkTrack){
1106     AliDebug(4, "No Track defined.");
1107     return NULL;
1108   }
1109   TH2 *h = NULL;
1110   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeCluster)))){
1111     AliWarning("No Histogram defined.");
1112     return NULL;
1113   }
1114   AliTRDseedV1 *tracklet = NULL;
1115   AliTRDcluster *c = NULL;
1116   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1117     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1118     for(Int_t ic(0); ic < AliTRDseedV1::kNtb; ic++){
1119       Bool_t kFIRST(kFALSE);
1120       if(!(c = tracklet->GetClusters(ic))) {
1121         if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1122       } else kFIRST = kTRUE;
1123       Float_t q(c->GetQ());
1124       if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) q+=c->GetQ();
1125       h->Fill(q, fCentralityClass);
1126     }
1127   }
1128   return h;
1129 }
1130
1131 //_______________________________________________________
1132 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1133   //
1134   // Plot the charge deposit per chamber
1135   //
1136   if(track) fkTrack = track;
1137   if(!fkTrack){
1138     AliDebug(4, "No Track defined.");
1139     return NULL;
1140   }
1141   TH2 *h = NULL;
1142   if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeTracklet)))){
1143     AliWarning("No Histogram defined.");
1144     return NULL;
1145   }
1146   AliTRDseedV1 *tracklet = NULL;
1147   AliTRDcluster *c = NULL;
1148   Double_t qTot = 0;
1149   Int_t nTracklets =fkTrack->GetNumberOfTracklets();
1150   for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1151     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1152     qTot = 0.;
1153     for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1154       if(!(c = tracklet->GetClusters(ic))) continue;
1155       qTot += TMath::Abs(c->GetQ());
1156     }
1157     h->Fill(qTot, fCentralityClass);
1158     if(DebugLevel() > 3){
1159       Int_t crossing = (Int_t)tracklet->IsRowCross();
1160       Int_t detector = tracklet->GetDetector();
1161       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1162       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1163       Float_t momentum = 0.;
1164       Int_t pdg = 0;
1165       Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1166       UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1167       if(fkMC){
1168               if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1169         pdg = fkMC->GetPDG();
1170       }
1171       (*DebugStream()) << "ChargeTracklet"
1172         << "Detector="  << detector
1173         << "Centrality="<< fCentralityClass
1174         << "crossing="  << crossing
1175         << "momentum="  << momentum
1176         << "nTracklets="<< nTracklets
1177         << "pdg="                               << pdg
1178         << "theta="                     << theta
1179         << "phi="                               << phi
1180         << "kinkIndex=" << kinkIndex
1181         << "TPCncls="           << nclsTPC
1182         << "QT="        << qTot
1183         << "\n";
1184     }
1185   }
1186   return h;
1187 }
1188
1189 //_______________________________________________________
1190 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1191   //
1192   // Plot the number of tracks per Sector
1193   //
1194   if(track) fkTrack = track;
1195   if(!fkTrack){
1196     AliDebug(4, "No Track defined.");
1197     return NULL;
1198   }
1199   TH1 *h = NULL;
1200   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1201     AliWarning("No Histogram defined.");
1202     return NULL;
1203   }
1204
1205   // TODO we should compare with
1206   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1207
1208   AliTRDseedV1 *tracklet = NULL;
1209   Int_t sector = -1;
1210   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1211     if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1212     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1213     break;
1214   }
1215   h->Fill(sector);
1216   return h;
1217 }
1218
1219
1220 //________________________________________________________
1221 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1222 {
1223   Float_t x = c->GetX();
1224   dist[0] = c->GetY() - tracklet->GetYat(x);
1225   dist[1] = c->GetZ() - tracklet->GetZat(x);
1226 }
1227
1228 //________________________________________________________
1229 void AliTRDcheckDET::GetEtaPhiAt(const AliExternalTrackParam *track, Double_t x, Double_t &eta, Double_t &phi){
1230   //
1231   // Get phi and eta at a given radial position
1232   //
1233   if(!fEvent) return;
1234   AliExternalTrackParam workpar(*track);
1235
1236   Double_t posLocal[3];
1237   Bool_t sucPos = workpar.GetXYZAt(x, fEvent->GetRunInfo()->GetMagneticField(), posLocal);
1238   Double_t sagPhi = sucPos ? TMath::ATan2(posLocal[1], posLocal[0]) : 0.;
1239   phi = sagPhi;
1240   eta = workpar.Eta();
1241 }
1242
1243
1244 //_______________________________________________________
1245 TH1* AliTRDcheckDET::MakePlotChi2() const
1246 {
1247 // Plot chi2/track normalized to number of degree of freedom
1248 // (tracklets) and compare with the theoretical distribution.
1249 //
1250 // Alex Bercuci <A.Bercuci@gsi.de>
1251
1252   return NULL;
1253
1254 /*  TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1255   TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1256   f.SetParLimits(1,1, 1e100);
1257   TLegend *leg = new TLegend(.7,.7,.95,.95);
1258   leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
1259   TH1D *h1 = NULL;
1260   Bool_t kFIRST = kTRUE;
1261   for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1262     h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1263     if(h1->Integral()<50) continue;
1264     h1->Scale(1./h1->Integral());
1265     h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1266     h1->SetLineColor(il);h1->SetLineStyle(2);
1267     f.SetParameter(1, .5*il);f.SetLineColor(il);
1268     h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1269     leg->AddEntry(h1, Form("%d", il), "l");
1270     if(kFIRST){
1271       h1->GetXaxis()->SetRangeUser(0., 25.);
1272     }
1273     kFIRST = kFALSE;
1274   }
1275   leg->Draw();
1276   gPad->SetLogy();
1277   return h1;*/
1278 }
1279
1280
1281 //________________________________________________________
1282 TH1* AliTRDcheckDET::MakePlotNTracklets(){
1283         //
1284         // Make nice bar plot of the number of tracklets in each method
1285         //
1286         TH3 *hTracklets3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
1287         if(!hTracklets3D){
1288                 AliError("Tracklet Histogram not found");
1289                 return NULL;
1290         }
1291         TH1 *hBAR = hTracklets3D->Project3D("y");
1292         hBAR->SetName("hBAR");
1293         hBAR->SetTitle("Number of Tracklets");
1294         if(hBAR->Integral()) hBAR->Scale(1./hBAR->Integral());
1295         hBAR->SetLineColor(kBlack);
1296         hBAR->SetLineWidth(2);
1297         hBAR->Draw();
1298
1299         // Draw also centrality-dependent plots
1300         TH1 *hBarCent[AliTRDeventInfo::kCentralityClasses];
1301         Int_t nHistsCentrality = 0;
1302         for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1303                 hTracklets3D->GetZaxis()->SetRange(icent+2, icent+2);
1304                 hBarCent[icent] = hTracklets3D->Project3D("y");
1305                 if(!(hBarCent[icent] && hBarCent[icent]->GetEntries())){
1306                         delete hBarCent[icent];
1307                         hBarCent[icent] = NULL;
1308                         continue;
1309                 }
1310                 hBarCent[icent]->SetName(Form("hBarCent_%d", icent));
1311                 hBarCent[icent]->SetTitle("Number of Tracklets");
1312                 hBarCent[icent]->Scale(1./hBarCent[icent]->Integral());
1313                 hBarCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1314                 hBarCent[icent]->Draw("same");
1315                 nHistsCentrality++;
1316         }
1317         hTracklets3D->GetZaxis()->SetRange(0, hTracklets3D->GetNbinsZ());
1318         AliInfo(Form("Number of Centrality Classes: %d", nHistsCentrality));
1319         if(nHistsCentrality){
1320                 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1321                 leg->SetBorderSize(0);
1322                 leg->SetFillStyle(0);
1323                 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++)
1324                         if(hBarCent[icent]) leg->AddEntry(hBarCent[icent], Form("Centrality Class %d", icent), "l");
1325                 leg->Draw();
1326                 gPad->Update();
1327         }
1328         return hBAR;
1329
1330         /*
1331           TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1332           TH1D *hBAR = tmp->ProjectionY();
1333           TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1334           TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1335           TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1336           leg->SetBorderSize(1);
1337           leg->SetFillColor(0);
1338
1339           Float_t scale = hCON->Integral();
1340           if(scale) hCON->Scale(100./scale);
1341           hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1342           hCON->SetBarWidth(0.2);
1343           hCON->SetBarOffset(0.6);
1344           hCON->SetStats(kFALSE);
1345           hCON->GetYaxis()->SetRangeUser(0.,40.);
1346           hCON->GetYaxis()->SetTitleOffset(1.2);
1347           hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1348           hCON->SetMaximum(55.);
1349
1350           if(scale) hBAR->Scale(100./scale);
1351           hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1352           hBAR->SetBarWidth(0.2);
1353           hBAR->SetBarOffset(0.2);
1354           hBAR->SetTitle("");
1355           hBAR->SetStats(kFALSE);
1356           hBAR->GetYaxis()->SetRangeUser(0.,40.);
1357           hBAR->GetYaxis()->SetTitleOffset(1.2);
1358           hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1359
1360           if(scale) hSTA->Scale(100./scale);
1361           hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1362           hSTA->SetBarWidth(0.2);
1363           hSTA->SetBarOffset(0.4);
1364           hSTA->SetTitle("");
1365           hSTA->SetStats(kFALSE);
1366           hSTA->GetYaxis()->SetRangeUser(0.,40.);
1367           hSTA->GetYaxis()->SetTitleOffset(1.2);
1368           hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1369           leg->Draw();
1370           gPad->Update();
1371           return hCON;
1372          */
1373 }
1374
1375 //________________________________________________________
1376 void AliTRDcheckDET::MakePlotNclustersTrack(){
1377         //
1378         // Plot number of clusters
1379         // Put histos from all centrality classes into one pad
1380         //
1381         TH2 *hClusters = dynamic_cast<TH2 *>(fContainer->FindObject("hNcls"));
1382         if(!hClusters){
1383                 AliError("Cluster histogram not found in the output");
1384         }
1385         TH1 *hAllCentrality = hClusters->ProjectionX("hNcls_all");
1386         hAllCentrality->SetTitle("Number of clusters/track");
1387         if(hAllCentrality->Integral()>0.) hAllCentrality->Scale(1./hAllCentrality->Integral());
1388         hAllCentrality->SetLineColor(kBlack);
1389         hAllCentrality->SetLineWidth(2);
1390         hAllCentrality->GetYaxis()->SetRangeUser(0, 0.02);
1391         hAllCentrality->Draw();
1392         PutTrendValue("NClustersTrack", hAllCentrality->GetMean());
1393         PutTrendValue("NClustersTrackRMS", hAllCentrality->GetRMS());
1394
1395         // Now look at single centrality classes
1396         TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
1397         Int_t nHistsCentrality = 0;
1398         for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1399                 hProjCentral[icent] = hClusters->ProjectionX(Form("hNcls_%d", icent), icent+1, icent+1);
1400                 if(!hProjCentral[icent]->GetEntries()){
1401                         delete hProjCentral[icent];
1402                         hProjCentral[icent] = NULL;
1403                         continue;
1404                 }
1405                 hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
1406                 hProjCentral[icent]->SetTitle("Number of clusters/track");
1407                 hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
1408                 hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.03);
1409                 hProjCentral[icent]->Draw("same");
1410                 nHistsCentrality++;
1411         }
1412         AliInfo(Form("%d centrality classes found", nHistsCentrality));
1413         if(nHistsCentrality){
1414                 // Draw nice legend
1415                 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1416                 leg->SetBorderSize(0);
1417                 leg->SetFillStyle(0);
1418                 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1419                         if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
1420                 }
1421                 leg->Draw();
1422                 gPad->Update();
1423         }
1424 }
1425
1426 //________________________________________________________
1427 void AliTRDcheckDET::MakePlotNclustersTracklet(){
1428         //
1429         // Plot number of clusters for different centrality classes
1430         //
1431         TH2 *hClusters = dynamic_cast<TH2*>(fContainer->FindObject("hNclTls"));
1432         if(!hClusters){
1433                 AliError("Histogram for the number of clusters per tracklet not available");
1434                 return;
1435         }
1436         TH1 *hAllCentrality = hClusters->ProjectionX("hNclsTls_all");
1437         hAllCentrality->SetTitle("Number of clusters/track");
1438         if(hAllCentrality->Integral()>0.) hAllCentrality->Scale(1./hAllCentrality->Integral());
1439         hAllCentrality->SetLineColor(kBlack);
1440         hAllCentrality->SetLineWidth(2);
1441         hAllCentrality->GetYaxis()->SetRangeUser(0, 0.3);
1442         hAllCentrality->Draw("pc");
1443         PutTrendValue("NClustersTracklet", hAllCentrality->GetMean());
1444         PutTrendValue("NClustersTrackletRMS", hAllCentrality->GetRMS());
1445
1446         // Now look at single centrality classes
1447         TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
1448         Int_t nHistsCentrality = 0;
1449         for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1450                 hProjCentral[icent] = hClusters->ProjectionX(Form("hNclsTls_%d", icent), icent+1, icent+1);
1451                 if(!hProjCentral[icent]->GetEntries()){
1452                         delete hProjCentral[icent];
1453                         hProjCentral[icent] = NULL;
1454                         continue;
1455                 }
1456                 if(hProjCentral[icent]->Integral()>0.) hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
1457                 hProjCentral[icent]->SetTitle("Number of clusters/track");
1458                 hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
1459                 hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.3);
1460                 hProjCentral[icent]->Draw("pcsame");
1461                 nHistsCentrality++;
1462         }
1463         AliInfo(Form("%d centrality classes found", nHistsCentrality));
1464         if(nHistsCentrality){
1465                 // Draw nice legend
1466                 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1467                 leg->SetBorderSize(0);
1468                 leg->SetFillStyle(0);
1469                 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1470                         if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
1471                 }
1472                 leg->Draw();
1473                 gPad->Update();
1474         }
1475 }
1476
1477 //________________________________________________________
1478 void AliTRDcheckDET::MakePlotTrackletCharge(){
1479         //
1480         // Draw tracklet charge for different centrality classes
1481         //
1482         TH2 *hQt = dynamic_cast<TH2*>(fContainer->FindObject("hQtrklt"));
1483         if(!hQt){
1484                 AliError("Histogram for tracklet charges not found");
1485                 return;
1486         }
1487         // First Project all charhes
1488         TH1 *hQtAll = hQt->ProjectionX("hQtAll");
1489         hQtAll->SetTitle("Tracklet Charge");
1490         Double_t scalefactor = hQtAll->GetMaximum()>0.?(0.7 * hQtAll->Integral() / hQtAll->GetMaximum()):1;
1491         if(hQtAll->Integral()>0.) hQtAll->Scale(scalefactor/hQtAll->Integral());
1492         hQtAll->GetYaxis()->SetRangeUser(0., 1.);
1493         hQtAll->SetLineColor(kBlack);
1494         hQtAll->SetLineWidth(2);
1495         hQtAll->Draw("c");
1496         PutTrendValue("ChargeTracklet", hQtAll->GetMaximumBin());
1497         PutTrendValue("ChargeTrackletRMS", hQtAll->GetRMS());
1498
1499         TH1 *hQtCent[AliTRDeventInfo::kCentralityClasses];
1500         Int_t nHistsCentrality = 0;
1501         for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1502                 hQtCent[icent] = hQt->ProjectionX(Form("hQt_%d", icent), icent+1, icent+1);
1503                 if(!hQtCent[icent]->GetEntries()){
1504                         delete hQtCent[icent]; hQtCent[icent]=NULL;
1505                         continue;
1506                 }
1507                 hQtCent[icent]->SetTitle("Tracklet Charge");
1508                 if(hQtCent[icent]->Integral()>0.) hQtCent[icent]->Scale(scalefactor/hQtCent[icent]->Integral());
1509                 hQtCent[icent]->GetYaxis()->SetRangeUser(0., 1.);
1510                 hQtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1511                 hQtCent[icent]->Draw("csame");
1512                 nHistsCentrality++;
1513         }
1514         AliInfo(Form("%d centrality classes found", nHistsCentrality));
1515         if(nHistsCentrality){
1516                 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1517                 leg->SetBorderSize(0);
1518                 leg->SetFillStyle(0);
1519                 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1520                         if(hQtCent[icent]) leg->AddEntry(hQtCent[icent], Form("Centrality Class %d", icent), "l");
1521                 }
1522                 leg->Draw();
1523                 gPad->Update();
1524         }
1525 }
1526
1527 //________________________________________________________
1528 void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1529   //
1530   // Plot abundance of tracks with number of tracklets as function of momentum
1531   //
1532
1533   Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1534   TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1535   for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1536     g[itl] = new TGraphErrors();
1537     g[itl]->SetLineColor(color[itl]);
1538     g[itl]->SetMarkerColor(color[itl]);
1539     g[itl]->SetMarkerStyle(20 + itl);
1540   }
1541
1542   TH3 *hBar3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
1543   if(!hBar3D){
1544           AliError("Histogram for the number of tracklets vs p not available");
1545           return;
1546   }
1547   TH2 *hBar = (TH2 *)hBar3D->Project3D("yx");
1548   TAxis *ax(hBar->GetXaxis());
1549   Int_t np(ax->GetNbins());
1550   for(Int_t ipBin(1); ipBin<np; ipBin++){
1551     h = hBar->ProjectionY("npBin", ipBin, ipBin);
1552     if(!Int_t(h->Integral())) continue;
1553     h->Scale(100./h->Integral());
1554     Float_t p(ax->GetBinCenter(ipBin));
1555     Float_t dp(ax->GetBinWidth(ipBin));
1556     Int_t ip(g[0]->GetN());
1557     for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1558       g[itl]->SetPoint(ip, p, h->GetBinContent(itl+1));
1559       g[itl]->SetPointError(ip, dp/3.46, h->GetBinError(itl+1));
1560     }
1561   }
1562
1563   TLegend *leg = new TLegend(0.76, 0.6, 1., 0.9);
1564   leg->SetBorderSize(0);
1565   leg->SetHeader("Tracklet/Track");
1566   leg->SetFillStyle(0);
1567   h = hBar->ProjectionX("npxBin"); h->Reset();
1568   h->SetTitle("");
1569   h->GetYaxis()->SetRangeUser(1., 99.);
1570   h->GetYaxis()->SetMoreLogLabels();
1571   h->GetYaxis()->CenterTitle();
1572   h->GetYaxis()->SetTitleOffset(1.2);
1573   h->SetYTitle("Prob. [%]");
1574   h->GetXaxis()->SetRangeUser(0.4, 12.);
1575   h->GetXaxis()->SetMoreLogLabels();
1576   h->GetXaxis()->CenterTitle();
1577   h->Draw("p");
1578   for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1579     g[itl]->Draw("pc");
1580     leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
1581   }
1582
1583   leg->Draw();
1584   gPad->SetLogx();gPad->SetLogy();
1585 }
1586
1587 //________________________________________________________
1588 Bool_t AliTRDcheckDET::MakePlotPulseHeight(){
1589         //
1590         // Create Plot of the Pluse Height Spectrum
1591         //
1592         TCanvas *output = gPad->GetCanvas();
1593         output->Divide(2);
1594         output->cd(1);
1595
1596         TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1597         //TH3F *hPhx = dynamic_cast<TH3F *>(arr->At(0)),
1598         TH3F *hPht = dynamic_cast<TH3F *>(arr->At(1));
1599         if(!hPht) return kFALSE;
1600         // Project centrality of the 2 histograms
1601         //TH2 *hProjCentX = dynamic_cast<TH2 *>(hPhx->Project3D("yx")),
1602         TH2 *hProjCentT = dynamic_cast<TH2 *>(hPht->Project3D("yx"));
1603   if(!hProjCentT) return kFALSE;
1604         //hProjCentX->SetName("hProjCentX");
1605         hProjCentT->SetName("hProjCentT");
1606         // Draw 2D histogram versus time on pad 2
1607         output->cd(2);
1608         hProjCentT->SetTitle("2D Pulse Height vs. Time");
1609         hProjCentT->GetYaxis()->SetTitleOffset(1.8);
1610   hProjCentT->GetYaxis()->SetRangeUser(0., 600.);
1611         hProjCentT->SetStats(kFALSE);
1612         hProjCentT->Draw("colz");
1613
1614   // Fit each slice of the pulse height spectrum with a Landau function
1615   TGraphErrors *landaufit = new TGraphErrors;
1616   landaufit->SetTitle();
1617   landaufit->SetMarkerColor(kWhite);
1618   landaufit->SetLineColor(kWhite);
1619   landaufit->SetMarkerStyle(25);
1620   TH1 *projection;
1621   Int_t ntime = 0;
1622   for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetLast(); it++){
1623     projection = hProjCentT->ProjectionY("proj", it, it);
1624     if(!projection->GetEntries()){
1625       delete projection;
1626       continue;
1627     }
1628     TF1 fitfun("fitfun", "landau", 0, 1000);
1629     fitfun.SetParLimits(0, 1e-1, 1e9); fitfun.SetParameter(0, 1000);
1630     fitfun.SetParLimits(1, 10, 200); fitfun.SetParameter(1, 50);
1631     fitfun.SetParLimits(2, 1e-1, 200); fitfun.SetParameter(2, 10);
1632     projection->Fit(&fitfun, "QN", "", 0, 1000);
1633     landaufit->SetPoint(ntime, hProjCentT->GetXaxis()->GetBinCenter(it), fitfun.GetParameter(1));
1634     landaufit->SetPointError(ntime, hProjCentT->GetXaxis()->GetBinWidth(it)/2., fitfun.GetParError(1));
1635     ntime++;
1636     delete projection;
1637   }
1638   landaufit->Draw("lpesame");
1639
1640         // Fill 1D PHS as function of time resp. radius (same binning of the 2 histograms)
1641         //TH1 *hPhsX = new TH1F("hPhsX", "Average PH vs X", hProjCentT->GetXaxis()->GetNbins(), hProjCentT->GetXaxis()->GetXbins()->GetArray()),
1642   Int_t nbinsT = hProjCentT->GetXaxis()->GetNbins();
1643         TH1     *hPhsT = new TH1F("hPhsT", "Average Ph vs Time; Time (100 ns); Average Pulse Height (a.u.)", nbinsT, -0.5, nbinsT - 0.5),
1644                 *htmp = NULL;
1645         /*for(Int_t irad = 1; irad <= hProjCentX->GetXaxis()->GetNbins(); irad++){
1646                 htmp = hProjCentX->ProjectionY("tmp", irad, irad);
1647                 hPhsX->SetBinContent(irad, htmp->GetMean());
1648                 hPhsX->SetBinError(irad, htmp->GetMeanError());
1649                 delete htmp;
1650         }*/
1651         for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetNbins(); it++){
1652                 htmp = hProjCentT->ProjectionY("tmp", it, it);
1653                 hPhsT->SetBinContent(it, htmp->GetMean());
1654                 hPhsT->SetBinError(it, htmp->GetMeanError());
1655                 delete htmp;
1656         }
1657         output->cd(1);
1658         // Draw 1D histograms
1659         if(hPhsT->GetEntries()){
1660                 hPhsT->SetMarkerStyle(24);
1661                 hPhsT->SetMarkerColor(kBlack);
1662                 hPhsT->SetLineColor(kBlack);
1663                 hPhsT->GetYaxis()->SetTitleOffset(1.5);
1664                 hPhsT->Draw("e1");
1665                 // Now fit the PHS with respect to time to get plateau and slope
1666                 // Trending for the pulse height: plateau value, slope and timebin of the maximum
1667                 TLinearFitter fit(1,"pol1");
1668                 Double_t time = 0.;
1669                 for(Int_t itime = 10; itime <= 20; itime++){
1670                         time = Double_t(itime);
1671                         Double_t err(hPhsT->GetBinError(itime + 1));
1672                         if(err>1.e-10) fit.AddPoint(&time, hPhsT->GetBinContent(itime + 1), err);
1673                 }
1674                 if(!fit.Eval()){
1675                         Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1676                         Double_t slope = fit.GetParameter(1);
1677                         PutTrendValue("PHplateau", plateau);
1678                         PutTrendValue("PHslope", slope);
1679                         PutTrendValue("PHamplificationPeak", static_cast<Double_t>(hPhsT->GetMaximumBin()-1));
1680                         AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(hPhsT->GetMaximumBin()-1)));
1681                 }
1682         }
1683         /*if(hPhsX->GetEntries()){
1684                 hPhsX->SetMarkerStyle(22);
1685                 hPhsX->SetMarkerColor(kBlue);
1686                 hPhsX->SetLineColor(kBlue);
1687                 hPhsX->GetYaxis()->SetTitleOffset(1.5);
1688                 hPhsX->Draw("e1same");
1689                 // create axis according to the histogram dimensions of the original second histogram
1690                 TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1691                                 gPad->GetUymax(),
1692                                 gPad->GetUxmax(),
1693                                 gPad->GetUymax(),
1694                                 -0.08, 4.88, 510,"-L");
1695                 axis->SetLineColor(kBlue);
1696                 axis->SetLabelColor(kBlue);
1697                 axis->SetTextColor(kBlue);
1698                 axis->SetTitle("x_{0}-x_{c} [cm]");
1699                 axis->Draw();
1700                 gPad->Update();
1701         }*/
1702
1703         // Centrality-dependent Pulse-Height Spectrum
1704         TH1 *hPhtCent[AliTRDeventInfo::kCentralityClasses];
1705         TH2 *hPtmp;
1706         Int_t nHistsCentrality = 0;
1707   TArrayD ptbinning(hPht->GetNbinsX()+1);
1708   for(Int_t ibin = 1; ibin <=hPht->GetNbinsX(); ibin++) ptbinning[ibin-1] = hPht->GetXaxis()->GetBinLowEdge(ibin);
1709   ptbinning[hPht->GetNbinsX()] = hPht->GetXaxis()->GetBinUpEdge(hPht->GetNbinsX());
1710         for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1711                 hPht->GetZaxis()->SetRange(icent+2,icent+2);
1712                 hPtmp = dynamic_cast<TH2*>(hPht->Project3D("yx"));
1713                 if(!(hPtmp && hPtmp->GetEntries())){
1714                         hPhtCent[icent] = NULL;
1715                         continue;
1716                 }
1717                 hPhtCent[icent] = new TH1F(Form("hPhtCent_%d", icent), "Average Ph vs Time", hPtmp->GetNbinsX(), ptbinning.GetArray());
1718                 for(Int_t it = 1; it <= hPtmp->GetNbinsX(); it++){
1719                         htmp = hPtmp->ProjectionY("htmp", it, it);
1720                         hPhtCent[icent]->SetBinContent(it, htmp->GetMean());
1721                         hPhtCent[icent]->SetBinError(it, htmp->GetMeanError());
1722                         delete htmp;
1723                 }
1724                 delete hPtmp;
1725                 hPhtCent[icent]->SetMarkerStyle(24);
1726                 hPhtCent[icent]->SetMarkerColor(fkColorsCentrality[icent]);
1727                 hPhtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1728                 hPhtCent[icent]->Draw("e1same");
1729                 nHistsCentrality++;
1730         }
1731         hPht->GetZaxis()->SetRange(0, hPht->GetNbinsZ());
1732         if(nHistsCentrality){
1733                 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1734                 leg->SetBorderSize(0);
1735                 leg->SetFillStyle(0);
1736                 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1737                         if(hPhtCent[icent]) leg->AddEntry(hPhtCent[icent], Form("Centrality Class %d", icent), "p");
1738                 }
1739                 leg->Draw();
1740                 gPad->Update();
1741         }
1742
1743         // delete 2D Projection of the PHS vs x since it is only used to calculate the 1D projection
1744         //delete hProjCentX;
1745
1746         return kTRUE;
1747 }
1748
1749 //________________________________________________________
1750 void AliTRDcheckDET::MakePlotMeanClustersLayer(){
1751   //
1752   // Create Summary plot for the mean number of clusters per layer
1753   //
1754   TCanvas *output = gPad->GetCanvas();
1755   output->Divide(3,2);
1756   TObjArray *histos = (TObjArray *)fContainer->At(kNclustersLayer);
1757   if(!histos){
1758     AliWarning("Histos for each layer not found");
1759     return;
1760   }
1761   TProfile2D *hlayer = NULL;
1762   for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
1763     if(!(hlayer = dynamic_cast<TProfile2D *>(histos->At(ily)))) continue;
1764     output->cd(ily + 1);
1765     gPad->SetGrid(0,0);
1766     hlayer->Draw("colz");
1767   }
1768 }
1769
1770 //________________________________________________________
1771 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1772   //
1773   // Draw nice bar plots
1774   //
1775   if(!histo->GetEntries()) return kFALSE;
1776   if(histo->Integral()>0.) histo->Scale(100./histo->Integral());
1777   histo->SetFillColor(color);
1778   histo->SetBarOffset(.2);
1779   histo->SetBarWidth(.6);
1780   histo->Draw("bar1");
1781   return kTRUE;
1782 }