PWG1/TRD/AliTRDcheckDET.cxx

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