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