TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TAxis.h>
   2 #include <TCanvas.h>
   3 #include <TFile.h>
   4 #include <TH1F.h>
   5 #include <TGaxis.h>
   6 #include <TGraph.h>
   7 #include <TMath.h>
   8 #include <TMap.h>
   9 #include <TObjArray.h>
  10 #include <TObject.h>
  11 #include <TObjString.h>
  12 #include <TPad.h>
  13 #include <TProfile.h>
  14 #include <TProfile2D.h>
  15 #include <TROOT.h>
  16
  17 #include "AliLog.h"
  18 #include "AliTRDcluster.h"
  19 #include "AliESDHeader.h"
  20 #include "AliESDRun.h"
  21 #include "AliESDtrack.h"
  22 #include "AliTRDgeometry.h"
  23 #include "AliTRDpadPlane.h"
  24 #include "AliTRDSimParam.h"
  25 #include "AliTRDseedV1.h"
  26 #include "AliTRDtrackV1.h"
  27 #include "AliTRDtrackerV1.h"
  28 #include "AliTRDReconstructor.h"
  29 #include "AliTrackReference.h"
  30 #include "AliTrackPointArray.h"
  31 #include "AliTracker.h"
  32 #include "TTreeStream.h"
  33
  34 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  35 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  36 #include "AliTRDcheckDetector.h"
  37
  38 #include <cstdio>
  39 #include <iostream>
  40
  41 ////////////////////////////////////////////////////////////////////////////
  42 //                                                                        //
  43 //  Reconstruction QA                                                     //
  44 //                                                                        //
  45 //  Task doing basic checks for tracking and detector performance         //
  46 //                                                                        //
  47 //  Authors:                                                              //
  48 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  49 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  50 //                                                                        //
  51 ////////////////////////////////////////////////////////////////////////////
  52
  53 //_______________________________________________________
  54 AliTRDcheckDetector::AliTRDcheckDetector():
  55   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  56   ,fEventInfo(0x0)
  57   ,fTriggerNames(0x0)
  58   ,fReconstructor(0x0)
  59   ,fGeo(0x0)
  60 {
  61   //
  62   // Default constructor
  63   //
  64   DefineInput(1,AliTRDeventInfo::Class());
  65   fReconstructor = new AliTRDReconstructor;
  66   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  67   fGeo = new AliTRDgeometry;
  68   InitFunctorList();
  69 }
  70
  71 //_______________________________________________________
  72 AliTRDcheckDetector::~AliTRDcheckDetector(){
  73   //
  74   // Destructor
  75   //
  76   if(fTriggerNames) delete fTriggerNames;
  77   delete fReconstructor;
  78   delete fGeo;
  79 }
  80
  81 //_______________________________________________________
  82 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
  83   //
  84   // Connect the Input data with the task
  85   //
  86   AliTRDrecoTask::ConnectInputData(opt);
  87   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
  88 }
  89
  90 //_______________________________________________________
  91 void AliTRDcheckDetector::CreateOutputObjects(){
  92   //
  93   // Create Output Objects
  94   //
  95   OpenFile(0,"RECREATE");
  96   fContainer = Histos();
  97   if(!fTriggerNames) fTriggerNames = new TMap();
  98 }
  99
 100 //_______________________________________________________
 101 void AliTRDcheckDetector::Exec(Option_t *opt){
 102   //
 103   // Execution function
 104   // Filling TRD quality histos
 105   //
 106   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 107   AliTRDrecoTask::Exec(opt);
 108   Int_t nTracks = 0;            // Count the number of tracks per event
 109   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 110   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 111   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
 112   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
 113   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 114     if(!fTracks->UncheckedAt(iti)) continue;
 115     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 116     if(!fTrackInfo->GetTrack()) continue;
 117     nTracks++;
 118   }
 119   if(nTracks){
 120     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
 121     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
 122   }
 123   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 124     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 125     // also set the label for both histograms
 126     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
 127     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 128     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
 129     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 130   }
 131   PostData(0, fContainer);
 132 }
 133
 134 //_______________________________________________________
 135 void AliTRDcheckDetector::Terminate(Option_t *){
 136   //
 137   // Terminate function
 138   //
 139 }
 140
 141 //_______________________________________________________
 142 Bool_t AliTRDcheckDetector::PostProcess(){
 143   //
 144   // Do Postprocessing (for the moment set the number of Reference histograms)
 145   //
 146
 147   TH1 * h = 0x0;
 148   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsTrack));
 149   if(h->GetEntries()) h->Scale(100./h->Integral());
 150
 151   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsCross));
 152   if(h->GetEntries()) h->Scale(100./h->Integral());
 153
 154   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksSector));
 155   if(h->GetEntries()) h->Scale(100./h->Integral());
 156
 157   // Calculate of the trigger clusters purity
 158   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
 159   TH1F *h1 = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
 160   h1->Divide(h);
 161   Float_t purities[20], val = 0;
 162   TString triggernames[20];
 163   Int_t nTriggerClasses = 0;
 164   for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
 165     if((val = h1->GetBinContent(ibin))){
 166       purities[nTriggerClasses] = val;
 167       triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
 168       nTriggerClasses++;
 169     }
 170   }
 171   h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
 172   TAxis *ax = h->GetXaxis();
 173   for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
 174     h->Fill(itrg, purities[itrg]);
 175     ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
 176   }
 177   ax->SetRangeUser(-0.5, nTriggerClasses+.5);
 178   h->GetYaxis()->SetRangeUser(0,1);
 179
 180   fNRefFigures = 11;
 181
 182   return kTRUE;
 183 }
 184
 185 //_______________________________________________________
 186 Bool_t AliTRDcheckDetector::GetRefFigure(Int_t ifig){
 187   //
 188   // Setting Reference Figures
 189   //
 190   TObjArray *arr = 0x0;
 191   TH1 *h = 0x0, *h1 = 0x0, *h2 = 0x0;
 192   TGaxis *axis = 0x0;
 193   switch(ifig){
 194   case kNclustersTrack:
 195     ((TH1F*)fContainer->At(kNclustersTrack))->Draw("pl");
 196     return kTRUE;
 197   case kNclustersTracklet:
 198     ((TH1F*)fContainer->At(kNclustersTracklet))->Draw("pc");
 199     return kTRUE;
 200   case kNtrackletsTrack:
 201     h = (TH1F*)fContainer->At(kNtrackletsTrack);
 202     if(!h->GetEntries()) break;
 203     h->SetFillColor(kGreen);
 204     h->SetBarOffset(.2);
 205     h->SetBarWidth(.6);
 206     h->Draw("bar1");
 207     return kTRUE;
 208   case kNtrackletsCross:
 209     h = (TH1F*)fContainer->At(kNtrackletsCross);
 210     if(!h->GetEntries()) break;
 211     h->SetFillColor(kRed);
 212     h->SetBarOffset(.2);
 213     h->SetBarWidth(.6);
 214     h->Draw("bar1");
 215     return kTRUE;
 216   case kNtrackletsFindable:
 217     h = (TH1F*)fContainer->At(kNtrackletsFindable);
 218     if(!h->GetEntries()) break;
 219     h->Scale(100./h->Integral());
 220     h->SetFillColor(kGreen);
 221     h->SetBarOffset(.2);
 222     h->SetBarWidth(.6);
 223     h->Draw("bar1");
 224     return kTRUE;
 225   case kNtracksEvent:
 226     ((TH1F*)fContainer->At(kNtracksEvent))->Draw("pl");
 227     return kTRUE;
 228   case kNtracksSector:
 229     h = (TH1F*)fContainer->At(kNtracksSector);
 230     if(!h->GetEntries()) break;
 231     h->SetFillColor(kGreen);
 232     h->SetBarOffset(.2);
 233     h->SetBarWidth(.6);
 234     h->Draw("bar1");
 235     return kTRUE;
 236   case kChi2:
 237     ((TH1F*)((TObjArray*)fContainer->At(kChi2))->At(0))->Draw("");
 238     return kTRUE;
 239   case kPH:
 240     arr = (TObjArray*)fContainer->At(kPH);
 241     h = (TH1F*)arr->At(0);
 242     h->SetMarkerStyle(24);
 243     h->SetMarkerColor(kBlack);
 244     h->SetLineColor(kBlack);
 245     h->Draw("e1");
 246     // copy the second histogram in a new one with the same x-dimension as the phs with respect to time
 247     h1 = (TH1F *)arr->At(1);
 248     h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
 249     for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
 250       h2->SetBinContent(ibin, h1->GetBinContent(ibin));
 251     h2->SetMarkerStyle(22);
 252     h2->SetMarkerColor(kBlue);
 253     h2->SetLineColor(kBlue);
 254     h2->Draw("e1same");
 255     gPad->Update();
 256     // create axis according to the histogram dimensions of the original second histogram
 257     axis = new TGaxis(gPad->GetUxmin(),
 258                       gPad->GetUymax(),
 259                       gPad->GetUxmax(),
 260                       gPad->GetUymax(),
 261                       -0.08, 4.88, 510,"-L");
 262     axis->SetLineColor(kBlue);
 263     axis->SetLabelColor(kBlue);
 264     axis->SetTextColor(kBlue);
 265     axis->SetTitle("x_{0}-x_{c} [cm]");
 266     axis->Draw();
 267     return kTRUE;
 268   case kChargeCluster:
 269     ((TH1F*)fContainer->At(kChargeCluster))->Draw("c");
 270     return kTRUE;
 271   case kChargeTracklet:
 272     ((TH1F*)fContainer->At(kChargeTracklet))->Draw("c");
 273     return kTRUE;
 274   case kNeventsTrigger:
 275     ((TH1F*)fContainer->At(kNeventsTrigger))->Draw("");
 276     return kTRUE;
 277   case kNeventsTriggerTracks:
 278     ((TH1F*)fContainer->At(kNeventsTriggerTracks))->Draw("");
 279     return kTRUE;
 280   case kTriggerPurity:
 281     h=(TH1F*)fContainer->At(kTriggerPurity);
 282     h->SetBarOffset(.2);
 283     h->SetBarWidth(.6);
 284     h->SetFillColor(kGreen);
 285     h->Draw("bar1");
 286     break;
 287   default:
 288     break;
 289   }
 290   AliInfo(Form("Reference plot [%d] missing result", ifig));
 291   return kFALSE;
 292 }
 293
 294 //_______________________________________________________
 295 TObjArray *AliTRDcheckDetector::Histos(){
 296   //
 297   // Create QA histograms
 298   //
 299   if(fContainer) return fContainer;
 300
 301   fContainer = new TObjArray(14);
 302   //fContainer->SetOwner(kTRUE);
 303
 304   // Register Histograms
 305   TH1 * h = 0x0;
 306   if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
 307     h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
 308     h->GetXaxis()->SetTitle("N_{clusters}");
 309     h->GetYaxis()->SetTitle("Entries");
 310   } else h->Reset();
 311   fContainer->AddAt(h, kNclustersTrack);
 312
 313   if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
 314     h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
 315     h->GetXaxis()->SetTitle("N_{clusters}");
 316     h->GetYaxis()->SetTitle("Entries");
 317   } else h->Reset();
 318   fContainer->AddAt(h, kNclustersTracklet);
 319
 320   if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
 321     h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
 322     h->GetXaxis()->SetTitle("N^{tracklet}");
 323     h->GetYaxis()->SetTitle("freq. [%]");
 324   } else h->Reset();
 325   fContainer->AddAt(h, kNtrackletsTrack);
 326
 327   //
 328   if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
 329     h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
 330     h->GetXaxis()->SetTitle("n_{row cross}");
 331     h->GetYaxis()->SetTitle("freq. [%]");
 332   } else h->Reset();
 333   fContainer->AddAt(h, kNtrackletsCross);
 334
 335   if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
 336     h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
 337     h->GetXaxis()->SetTitle("r [a.u]");
 338     h->GetYaxis()->SetTitle("Entries");
 339   } else h->Reset();
 340   fContainer->AddAt(h, kNtrackletsFindable);
 341
 342   if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
 343     h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
 344     h->GetXaxis()->SetTitle("N_{tracks}");
 345     h->GetYaxis()->SetTitle("Entries");
 346   } else h->Reset();
 347   fContainer->AddAt(h, kNtracksEvent);
 348
 349   if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
 350     h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
 351     h->GetXaxis()->SetTitle("sector");
 352     h->GetYaxis()->SetTitle("freq. [%]");
 353   } else h->Reset();
 354   fContainer->AddAt(h, kNtracksSector);
 355
 356   // <PH> histos
 357   TObjArray *arr = new TObjArray(2);
 358   arr->SetOwner(kTRUE);  arr->SetName("<PH>");
 359   fContainer->AddAt(arr, kPH);
 360   if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
 361     h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
 362     h->GetXaxis()->SetTitle("Time / 100ns");
 363     h->GetYaxis()->SetTitle("<PH> [a.u]");
 364   } else h->Reset();
 365   arr->AddAt(h, 0);
 366   if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
 367     h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
 368   else h->Reset();
 369   arr->AddAt(h, 1);
 370
 371   // Chi2 histos
 372   arr = new TObjArray(2);
 373   arr->SetOwner(kTRUE); arr->SetName("Chi2");
 374   fContainer->AddAt(arr, kChi2);
 375   if(!(h = (TH1F *)gROOT->FindObject("hChi2")))
 376     h = new TH1F("hChi2", "#Chi2", 200, 0, 20);
 377   else h->Reset();
 378   arr->AddAt(h, 0);
 379   if(!(h = (TH1F *)gROOT->FindObject("hChi2n")))
 380     h = new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5);
 381   else h->Reset();
 382   arr->AddAt(h, 1);
 383
 384
 385   if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
 386     h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
 387     h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
 388     h->GetYaxis()->SetTitle("Entries");
 389   }else h->Reset();
 390   fContainer->AddAt(h, kChargeCluster);
 391
 392   if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
 393     h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
 394     h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
 395     h->GetYaxis()->SetTitle("Entries");
 396   }else h->Reset();
 397   fContainer->AddAt(h, kChargeTracklet);
 398
 399
 400   if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
 401     h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
 402   else h->Reset();
 403   fContainer->AddAt(h, kNeventsTrigger);
 404
 405   if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
 406     h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
 407   else h->Reset();
 408   fContainer->AddAt(h, kNeventsTriggerTracks);
 409
 410   if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
 411     h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
 412     h->GetXaxis()->SetTitle("Trigger Cluster");
 413     h->GetYaxis()->SetTitle("freq.");
 414   } else h->Reset();
 415   fContainer->AddAt(h, kTriggerPurity);
 416
 417   return fContainer;
 418 }
 419
 420 /*
 421 * Plotting Functions
 422 */
 423
 424 //_______________________________________________________
 425 TH1 *AliTRDcheckDetector::PlotNClustersTracklet(const AliTRDtrackV1 *track){
 426   //
 427   // Plot the mean number of clusters per tracklet
 428   //
 429   if(track) fTrack = track;
 430   if(!fTrack){
 431     AliWarning("No Track defined.");
 432     return 0x0;
 433   }
 434   TH1 *h = 0x0;
 435   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
 436     AliWarning("No Histogram defined.");
 437     return 0x0;
 438   }
 439   AliTRDseedV1 *tracklet = 0x0;
 440   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 441     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 442     h->Fill(tracklet->GetN2());
 443   }
 444   return h;
 445 }
 446
 447 //_______________________________________________________
 448 TH1 *AliTRDcheckDetector::PlotNClustersTrack(const AliTRDtrackV1 *track){
 449   //
 450   // Plot the number of clusters in one track
 451   //
 452   if(track) fTrack = track;
 453   if(!fTrack){
 454     AliWarning("No Track defined.");
 455     return 0x0;
 456   }
 457   TH1 *h = 0x0;
 458   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
 459     AliWarning("No Histogram defined.");
 460     return 0x0;
 461   }
 462
 463   Int_t nclusters = 0;
 464   AliTRDseedV1 *tracklet = 0x0;
 465   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 466     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 467     nclusters += tracklet->GetN();
 468     if(fDebugLevel > 2){
 469       Int_t crossing = Int_t(tracklet->IsRowCross());
 470       Int_t detector = tracklet->GetDetector();
 471       Float_t theta = TMath::ATan(tracklet->GetZref(1));
 472       Float_t phi = TMath::ATan(tracklet->GetYref(1));
 473       Float_t momentum = 0.;
 474       Int_t pdg = 0;
 475       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 476       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 477       if(fMC){
 478         if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 479         pdg = fMC->GetPDG();
 480       }
 481       (*fDebugStream) << "NClustersTrack"
 482         << "Detector="  << detector
 483         << "crossing="  << crossing
 484         << "momentum="  << momentum
 485         << "pdg="                               << pdg
 486         << "theta="                     << theta
 487         << "phi="                               << phi
 488         << "kinkIndex=" << kinkIndex
 489         << "TPCncls="           << TPCncls
 490         << "nclusters=" << nclusters
 491         << "\n";
 492     }
 493   }
 494   h->Fill(nclusters);
 495   return h;
 496 }
 497
 498
 499 //_______________________________________________________
 500 TH1 *AliTRDcheckDetector::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
 501   //
 502   // Plot the number of tracklets
 503   //
 504   if(track) fTrack = track;
 505   if(!fTrack){
 506     AliWarning("No Track defined.");
 507     return 0x0;
 508   }
 509   TH1 *h = 0x0;
 510   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
 511     AliWarning("No Histogram defined.");
 512     return 0x0;
 513   }
 514   Int_t nTracklets = fTrack->GetNumberOfTracklets();
 515   h->Fill(nTracklets);
 516   if(fDebugLevel > 3){
 517     if(nTracklets == 1){
 518       // If we have one Tracklet, check in which layer this happens
 519       Int_t layer = -1;
 520       AliTRDseedV1 *tracklet = 0x0;
 521       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 522         if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
 523       }
 524       (*fDebugStream) << "NTrackletsTrack"
 525         << "Layer=" << layer
 526         << "\n";
 527     }
 528   }
 529   return h;
 530 }
 531
 532
 533 //_______________________________________________________
 534 TH1 *AliTRDcheckDetector::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
 535   //
 536   // Plot the number of tracklets
 537   //
 538   if(track) fTrack = track;
 539   if(!fTrack){
 540     AliWarning("No Track defined.");
 541     return 0x0;
 542   }
 543   TH1 *h = 0x0;
 544   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
 545     AliWarning("No Histogram defined.");
 546     return 0x0;
 547   }
 548
 549   Int_t ncross = 0;
 550   AliTRDseedV1 *tracklet = 0x0;
 551   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 552     if(!(tracklet = fTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
 553
 554     if(tracklet->IsRowCross()) ncross++;
 555   }
 556   h->Fill(ncross);
 557   return h;
 558 }
 559
 560 //_______________________________________________________
 561 TH1 *AliTRDcheckDetector::PlotFindableTracklets(const AliTRDtrackV1 *track){
 562   //
 563   // Plots the ratio of number of tracklets vs.
 564   // number of findable tracklets
 565   //
 566   // Findable tracklets are defined as track prolongation
 567   // to layer i does not hit the dead area +- epsilon
 568   //
 569   // In order to check whether tracklet hist active area in Layer i,
 570   // the track is refitted and the fitted position + an uncertainty
 571   // range is compared to the chamber border (also with a different
 572   // uncertainty)
 573   //
 574   // For the track fit two cases are distinguished:
 575   // If the track is a stand alone track (defined by the status bit
 576   // encoding, then the track is fitted with the tilted Rieman model
 577   // Otherwise the track is fitted with the Kalman fitter in two steps:
 578   // Since the track parameters are give at the outer point, we first
 579   // fit in direction inwards. Afterwards we fit again in direction outwards
 580   // to extrapolate the track to layers which are not reached by the track
 581   // For the Kalman model, the radial track points have to be shifted by
 582   // a distance epsilon in the direction that we want to fit
 583   //
 584   const Float_t epsilon = 0.01;   // dead area tolerance
 585   const Float_t epsilon_R = 1;    // shift in radial direction of the anode wire position (Kalman filter only)
 586   const Float_t delta_y = 0.7;    // Tolerance in the track position in y-direction
 587   const Float_t delta_z = 7.0;    // Tolerance in the track position in z-direction (Padlength)
 588   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 589
 590   if(track) fTrack = track;
 591   if(!fTrack){
 592     AliWarning("No Track defined.");
 593     return 0x0;
 594   }
 595   TH1 *h = 0x0;
 596   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
 597     AliWarning("No Histogram defined.");
 598     return 0x0;
 599   }
 600   Int_t nFound = 0, nFindable = 0;
 601   Int_t stack = -1;
 602   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 603   Double_t y = 0., z = 0.;
 604   AliTRDseedV1 *tracklet = 0x0;
 605   AliTRDpadPlane *pp;
 606   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 607     if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
 608       tracklet->SetReconstructor(fReconstructor);
 609       nFound++;
 610     }
 611   }
 612   // 2 Different cases:
 613   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 614   // 2nd barrel track: here we propagate the track to the layers
 615   AliTrackPoint points[6];
 616   Float_t xyz[3];
 617   memset(xyz, 0, sizeof(Float_t) * 3);
 618   if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 619     // stand alone track
 620     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 621       xyz[0] = x_anode[il];
 622       points[il].SetXYZ(xyz);
 623     }
 624     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
 625   } else {
 626     // barrel track
 627     //
 628     // 2 Steps:
 629     // -> Kalman inwards
 630     // -> Kalman outwards
 631     AliTRDtrackV1 copy_track(*fTrack);  // Do Kalman on a (non-constant) copy of the track
 632     AliTrackPoint points_inward[6], points_outward[6];
 633     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 634       // In order to avoid complications in the Kalman filter if the track points have the same radial
 635       // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
 636       // in the direction we want to go
 637       // The track points have to be in reverse order for the Kalman Filter inwards
 638       xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
 639       points_inward[il].SetXYZ(xyz);
 640       xyz[0] = x_anode[il] + epsilon_R;
 641       points_outward[il].SetXYZ(xyz);
 642     }
 643     /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 644       printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
 645     // Kalman inwards
 646     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kFALSE, 6, points_inward);
 647     memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
 648     // Kalman outwards
 649     AliTRDtrackerV1::FitKalman(&copy_track, 0x0, kTRUE, 6, points_inward);
 650     memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 651   }
 652   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 653     y = points[il].GetY();
 654     z = points[il].GetZ();
 655     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
 656     pp = fGeo->GetPadPlane(il, stack);
 657     ymin = pp->GetCol0() + epsilon;
 658     ymax = pp->GetColEnd() - epsilon;
 659     zmin = pp->GetRowEnd() + epsilon;
 660     zmax = pp->GetRow0() - epsilon;
 661     // ignore y-crossing (material)
 662     if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
 663       if(fDebugLevel > 3){
 664         Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetMeanz() : 0};
 665         Int_t hasTracklet = tracklet ? 1 : 0;
 666         (*fDebugStream)   << "FindableTracklets"
 667           << "layer="     << il
 668           << "ytracklet=" << pos_tracklet[0]
 669           << "ytrack="    << y
 670           << "ztracklet=" << pos_tracklet[1]
 671           << "ztrack="    << z
 672           << "tracklet="  << hasTracklet
 673           << "\n";
 674       }
 675   }
 676
 677   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
 678   if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 679   return h;
 680 }
 681
 682
 683 //_______________________________________________________
 684 TH1 *AliTRDcheckDetector::PlotChi2(const AliTRDtrackV1 *track){
 685   //
 686   // Plot the chi2 of the track
 687   //
 688   if(track) fTrack = track;
 689   if(!fTrack){
 690     AliWarning("No Track defined.");
 691     return 0x0;
 692   }
 693   TH1 *h = 0x0;
 694   if(!(h = dynamic_cast<TH1F *>(((TObjArray*)(fContainer->At(kChi2)))->At(0)))){
 695     AliWarning("No Histogram defined.");
 696     return 0x0;
 697   }
 698   h->Fill(fTrack->GetChi2());
 699   return h;
 700 }
 701
 702 //_______________________________________________________
 703 TH1 *AliTRDcheckDetector::PlotChi2Norm(const AliTRDtrackV1 *track){
 704   //
 705   // Plot the chi2 of the track
 706   //
 707   if(track) fTrack = track;
 708   if(!fTrack){
 709     AliWarning("No Track defined.");
 710     return 0x0;
 711   }
 712   TH1 *h = 0x0;
 713   if(!(h = dynamic_cast<TH1F *>(((TObjArray*)(fContainer->At(kChi2)))->At(1)))){
 714     AliWarning("No Histogram defined.");
 715     return 0x0;
 716   }
 717   Int_t nTracklets = 0;
 718   AliTRDseedV1 *tracklet = 0x0;
 719   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 720     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 721     nTracklets++;
 722   }
 723   h->Fill(fTrack->GetChi2()/nTracklets);
 724   return h;
 725 }
 726
 727
 728 //_______________________________________________________
 729 TH1 *AliTRDcheckDetector::PlotPHt(const AliTRDtrackV1 *track){
 730   //
 731   // Plot the average pulse height
 732   //
 733   if(track) fTrack = track;
 734   if(!fTrack){
 735     AliWarning("No Track defined.");
 736     return 0x0;
 737   }
 738   TProfile *h = 0x0;
 739   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
 740     AliWarning("No Histogram defined.");
 741     return 0x0;
 742   }
 743   AliTRDseedV1 *tracklet = 0x0;
 744   AliTRDcluster *c = 0x0;
 745   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 746     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 747     Int_t crossing = Int_t(tracklet->IsRowCross());
 748     Int_t detector = tracklet->GetDetector();
 749     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 750       if(!(c = tracklet->GetClusters(itime))) continue;
 751       Int_t localtime        = c->GetLocalTimeBin();
 752       Double_t absolute_charge = TMath::Abs(c->GetQ());
 753       h->Fill(localtime, absolute_charge);
 754       if(fDebugLevel > 3){
 755         Double_t distance[2];
 756         GetDistanceToTracklet(distance, tracklet, c);
 757         Float_t theta = TMath::ATan(tracklet->GetZref(1));
 758         Float_t phi = TMath::ATan(tracklet->GetYref(1));
 759         Float_t momentum = 0.;
 760         Int_t pdg = 0;
 761         Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 762         UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 763         if(fMC){
 764           if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 765           pdg = fMC->GetPDG();
 766         }
 767         (*fDebugStream) << "PHt"
 768           << "Detector="        << detector
 769           << "crossing="        << crossing
 770           << "Timebin="         << localtime
 771           << "Charge="          << absolute_charge
 772           << "momentum="        << momentum
 773           << "pdg="                             << pdg
 774           << "theta="                   << theta
 775           << "phi="                             << phi
 776           << "kinkIndex="       << kinkIndex
 777           << "TPCncls="         << TPCncls
 778           << "dy="        << distance[0]
 779           << "dz="        << distance[1]
 780           << "c.="        << c
 781           << "\n";
 782       }
 783     }
 784   }
 785   return h;
 786 }
 787
 788 //_______________________________________________________
 789 TH1 *AliTRDcheckDetector::PlotPHx(const AliTRDtrackV1 *track){
 790   //
 791   // Plots the average pulse height vs the distance from the anode wire
 792   // (plus const anode wire offset)
 793   //
 794   if(track) fTrack = track;
 795   if(!fTrack){
 796     AliWarning("No Track defined.");
 797     return 0x0;
 798   }
 799   TProfile *h = 0x0;
 800   if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
 801     AliWarning("No Histogram defined.");
 802     return 0x0;
 803   }
 804   Float_t offset = .5*AliTRDgeometry::CamHght();
 805   AliTRDseedV1 *tracklet = 0x0;
 806   AliTRDcluster *c = 0x0;
 807   Double_t distance = 0;
 808   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 809     if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
 810     tracklet->ResetClusterIter();
 811     while((c = tracklet->NextCluster())){
 812       distance = tracklet->GetX0() - c->GetX() + offset;
 813       h->Fill(distance, TMath::Abs(c->GetQ()));
 814     }
 815   }
 816   return h;
 817 }
 818
 819 //_______________________________________________________
 820 TH1 *AliTRDcheckDetector::PlotChargeCluster(const AliTRDtrackV1 *track){
 821   //
 822   // Plot the cluster charge
 823   //
 824   if(track) fTrack = track;
 825   if(!fTrack){
 826     AliWarning("No Track defined.");
 827     return 0x0;
 828   }
 829   TH1 *h = 0x0;
 830   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
 831     AliWarning("No Histogram defined.");
 832     return 0x0;
 833   }
 834   AliTRDseedV1 *tracklet = 0x0;
 835   AliTRDcluster *c = 0x0;
 836   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 837     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 838     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 839       if(!(c = tracklet->GetClusters(itime))) continue;
 840       h->Fill(c->GetQ());
 841     }
 842   }
 843   return h;
 844 }
 845
 846 //_______________________________________________________
 847 TH1 *AliTRDcheckDetector::PlotChargeTracklet(const AliTRDtrackV1 *track){
 848   //
 849   // Plot the charge deposit per chamber
 850   //
 851   if(track) fTrack = track;
 852   if(!fTrack){
 853     AliWarning("No Track defined.");
 854     return 0x0;
 855   }
 856   TH1 *h = 0x0;
 857   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
 858     AliWarning("No Histogram defined.");
 859     return 0x0;
 860   }
 861   AliTRDseedV1 *tracklet = 0x0;
 862   AliTRDcluster *c = 0x0;
 863   Double_t Qtot = 0;
 864   Int_t nTracklets =fTrack->GetNumberOfTracklets();
 865   for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
 866     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 867     Qtot = 0.;
 868     for(Int_t ic = AliTRDseed::knTimebins; ic--;){
 869       if(!(c = tracklet->GetClusters(ic))) continue;
 870       Qtot += TMath::Abs(c->GetQ());
 871     }
 872     h->Fill(Qtot);
 873     if(fDebugLevel > 3){
 874       Int_t crossing = (Int_t)tracklet->IsRowCross();
 875       Int_t detector = tracklet->GetDetector();
 876       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 877       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 878       Float_t momentum = 0.;
 879       Int_t pdg = 0;
 880       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 881       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 882       if(fMC){
 883               if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 884         pdg = fMC->GetPDG();
 885       }
 886       (*fDebugStream) << "ChargeTracklet"
 887         << "Detector="  << detector
 888         << "crossing="  << crossing
 889         << "momentum="  << momentum
 890         << "nTracklets="<< nTracklets
 891         << "pdg="                               << pdg
 892         << "theta="                     << theta
 893         << "phi="                               << phi
 894         << "kinkIndex=" << kinkIndex
 895         << "TPCncls="           << TPCncls
 896         << "QT="        << Qtot
 897         << "\n";
 898     }
 899   }
 900   return h;
 901 }
 902
 903 //_______________________________________________________
 904 TH1 *AliTRDcheckDetector::PlotNTracksSector(const AliTRDtrackV1 *track){
 905   //
 906   // Plot the number of tracks per Sector
 907   //
 908   if(track) fTrack = track;
 909   if(!fTrack){
 910     AliWarning("No Track defined.");
 911     return 0x0;
 912   }
 913   TH1 *h = 0x0;
 914   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
 915     AliWarning("No Histogram defined.");
 916     return 0x0;
 917   }
 918
 919   // TODO we should compare with
 920   // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
 921
 922   AliTRDseedV1 *tracklet = 0x0;
 923   Int_t sector = -1;
 924   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 925     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 926     sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
 927     break;
 928   }
 929   h->Fill(sector);
 930   return h;
 931 }
 932
 933
 934 //________________________________________________________
 935 void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
 936 {
 937
 938   fReconstructor->SetRecoParam(r);
 939 }
 940
 941 //________________________________________________________
 942 void AliTRDcheckDetector::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c)
 943 {
 944   Float_t x = c->GetX();
 945   dist[0] = c->GetY() - tracklet->GetYat(x);
 946   dist[1] = c->GetZ() - tracklet->GetZat(x);
 947 }