TRD/qaRec/AliTRDcheckDetector.cxx

   1 #include <TFile.h>
   2 #include <TH1F.h>
   3 #include <TGraph.h>
   4 #include <TMath.h>
   5 #include <TMap.h>
   6 #include <TObjArray.h>
   7 #include <TObject.h>
   8 #include <TObjString.h>
   9 #include <TProfile.h>
  10 #include <TProfile2D.h>
  11 #include <TROOT.h>
  12
  13 #include "AliLog.h"
  14 #include "AliTRDcluster.h"
  15 #include "AliESDHeader.h"
  16 #include "AliESDRun.h"
  17 #include "AliESDtrack.h"
  18 #include "AliTRDgeometry.h"
  19 #include "AliTRDpadPlane.h"
  20 #include "AliTRDseedV1.h"
  21 #include "AliTRDtrackV1.h"
  22 #include "AliTRDtrackerV1.h"
  23 #include "AliTRDReconstructor.h"
  24 #include "AliTrackReference.h"
  25 #include "AliTrackPointArray.h"
  26 #include "AliTracker.h"
  27 #include "TTreeStream.h"
  28
  29 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  30 #include "AliTRDtrackInfo/AliTRDeventInfo.h"
  31 #include "AliTRDcheckDetector.h"
  32
  33 #include <cstdio>
  34 #include <iostream>
  35
  36 ////////////////////////////////////////////////////////////////////////////
  37 //                                                                        //
  38 //  Reconstruction QA                                                     //
  39 //                                                                        //
  40 //  Task doing basic checks for tracking and detector performance         //
  41 //                                                                        //
  42 //  Authors:                                                              //
  43 //    Anton Andronic <A.Andronic@gsi.de>                                  //
  44 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  45 //                                                                        //
  46 ////////////////////////////////////////////////////////////////////////////
  47
  48 //_______________________________________________________
  49 AliTRDcheckDetector::AliTRDcheckDetector():
  50   AliTRDrecoTask("DetChecker", "Basic Detector Checker")
  51   ,fEventInfo(0x0)
  52   ,fTriggerNames(0x0)
  53   ,fReconstructor(0x0)
  54 {
  55   //
  56   // Default constructor
  57   //
  58   DefineInput(1,AliTRDeventInfo::Class());
  59   fReconstructor = new AliTRDReconstructor;
  60   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  61   InitFunctorList();
  62 }
  63
  64 //_______________________________________________________
  65 AliTRDcheckDetector::~AliTRDcheckDetector(){
  66   //
  67   // Destructor
  68   //
  69   if(fEventInfo) delete fEventInfo;
  70   if(fTriggerNames) delete fTriggerNames;
  71   delete fReconstructor;
  72 }
  73
  74 //_______________________________________________________
  75 void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
  76   //
  77   // Connect the Input data with the task
  78   //
  79   AliTRDrecoTask::ConnectInputData(opt);
  80   fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
  81 }
  82
  83 //_______________________________________________________
  84 void AliTRDcheckDetector::CreateOutputObjects(){
  85   //
  86   // Create Output Objects
  87   //
  88   OpenFile(0,"RECREATE");
  89   fContainer = Histos();
  90   if(!fTriggerNames) fTriggerNames = new TMap();
  91 }
  92
  93 //_______________________________________________________
  94 void AliTRDcheckDetector::Exec(Option_t *opt){
  95   //
  96   // Execution function
  97   // Filling TRD quality histos
  98   //
  99   if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
 100   AliTRDrecoTask::Exec(opt);
 101   Int_t nTracks = 0;            // Count the number of tracks per event
 102   Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
 103   TString triggername =  fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
 104   if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
 105   dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger))->Fill(triggermask);
 106   for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
 107     if(!fTracks->UncheckedAt(iti)) continue;
 108     AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
 109     if(!fTrackInfo->GetTrack()) continue;
 110     nTracks++;
 111   }
 112   if(nTracks){
 113     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks))->Fill(triggermask);
 114     dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist))->Fill(nTracks);
 115   }
 116   if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
 117     fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
 118     // also set the label for both histograms
 119     TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 120     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 121     histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 122     histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
 123   }
 124   PostData(0, fContainer);
 125 }
 126
 127 //_______________________________________________________
 128 void AliTRDcheckDetector::Terminate(Option_t *){
 129   //
 130   // Terminate function
 131   //
 132 }
 133
 134 //_______________________________________________________
 135 Bool_t AliTRDcheckDetector::PostProcess(){
 136   //
 137   // Do Postprocessing (for the moment set the number of Reference histograms)
 138   //
 139
 140   TH1 * histo = 0x0;
 141   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksEventHist));
 142   histo->GetXaxis()->SetTitle("Number of Tracks");
 143   histo->GetYaxis()->SetTitle("Events");
 144   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclustersHist));
 145   histo->GetXaxis()->SetTitle("Number of Clusters");
 146   histo->GetYaxis()->SetTitle("Entries");
 147   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtrackletsHist));
 148   histo->GetXaxis()->SetTitle("Number of Tracklets");
 149   histo->GetYaxis()->SetTitle("Entries");
 150   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNclusterTrackletHist));
 151   histo->GetXaxis()->SetTitle("Number of Clusters");
 152   histo->GetYaxis()->SetTitle("Entries");
 153   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChi2));
 154   histo->GetXaxis()->SetTitle("#chi^2");
 155   histo->GetYaxis()->SetTitle("Entries");
 156   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNTracksSectorHist));
 157   histo->GetXaxis()->SetTitle("Sector");
 158   histo->GetYaxis()->SetTitle("Number of Tracks");
 159   histo = dynamic_cast<TProfile *>(fContainer->UncheckedAt(kPulseHeight));
 160   histo->GetXaxis()->SetTitle("Time / 100ns");
 161   histo->GetYaxis()->SetTitle("Average Pulse Height (a. u.)");
 162   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kClusterCharge));
 163   histo->GetXaxis()->SetTitle("Cluster Charge (a.u.)");
 164   histo->GetYaxis()->SetTitle("Entries");
 165   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kChargeDeposit));
 166   histo->GetXaxis()->SetTitle("Charge Deposit (a.u.)");
 167   histo->GetYaxis()->SetTitle("Entries");
 168
 169   // Calculate the purity of the trigger clusters
 170   histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTrigger));
 171   TH1F *histoTracks = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNEventsTriggerTracks));
 172   histoTracks->Divide(histo);
 173   Float_t purities[20], val = 0;
 174   TString triggernames[20];
 175   Int_t nTriggerClasses = 0;
 176   for(Int_t ibin = 1; ibin <= histo->GetNbinsX(); ibin++){
 177     if((val = histoTracks->GetBinContent(ibin))){
 178       purities[nTriggerClasses] = val;
 179       triggernames[nTriggerClasses] = histoTracks->GetXaxis()->GetBinLabel(ibin);
 180       nTriggerClasses++;
 181     }
 182   }
 183   TH1F *hTriggerInf = new TH1F("fTriggerInf", "Trigger Information", TMath::Max(nTriggerClasses, 1), 0, TMath::Max(nTriggerClasses, 1));
 184   for(Int_t ibin = 1; ibin <= nTriggerClasses; ibin++){
 185     hTriggerInf->SetBinContent(ibin, purities[ibin-1]);
 186     hTriggerInf->GetXaxis()->SetBinLabel(ibin, triggernames[ibin-1].Data());
 187   }
 188   hTriggerInf->GetXaxis()->SetTitle("Trigger Cluster");
 189   hTriggerInf->GetYaxis()->SetTitle("Ratio");
 190   hTriggerInf->GetYaxis()->SetRangeUser(0,1);
 191 //      hTriggerInf->SetMarkerColor(kBlue);
 192 //      hTriggerInf->SetMarkerStyle(22);
 193   fContainer->Add(hTriggerInf);
 194   fNRefFigures = 10;
 195   return kTRUE;
 196 }
 197
 198 //_______________________________________________________
 199 void AliTRDcheckDetector::GetRefFigure(Int_t ifig){
 200   //
 201   // Setting Reference Figures
 202   //
 203   TH1 *h = 0x0;
 204   switch(ifig){
 205   case 0:
 206     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 207     break;
 208   case 1:
 209     ((TH1F*)fContainer->At(kNclustersHist))->Draw("pl");
 210     break;
 211   case 2:
 212     h = (TH1F*)fContainer->At(kNtrackletsHist);
 213     if(!h->GetEntries()) break;
 214     h->Scale(100./h->Integral());
 215     h->GetXaxis()->SetRangeUser(.5, 6.5);
 216     h->SetFillColor(kGreen);
 217     h->SetBarOffset(.2);
 218     h->SetBarWidth(.6);
 219     h->Draw("bar1");
 220     break;
 221   case 3:
 222     ((TH1F*)fContainer->At(kNclusterTrackletHist))->Draw("pc");
 223     break;
 224   case 4:
 225     ((TH1F*)fContainer->At(kChi2))->Draw("");
 226     break;
 227   case 5:
 228     h = (TH1F*)fContainer->At(kNTracksSectorHist);
 229     if(!h->GetEntries()) break;
 230     h->Scale(100./h->Integral());
 231     h->SetFillColor(kGreen);
 232     h->SetBarOffset(.2);
 233     h->SetBarWidth(.6);
 234     h->Draw("bar1");
 235     break;
 236   case 6:
 237     h = (TH1F*)fContainer->At(kPulseHeight);
 238     h->SetMarkerStyle(24);
 239     h->Draw("e1");
 240     break;
 241   case 7:
 242     ((TH1F*)fContainer->At(kClusterCharge))->Draw("c");
 243     break;
 244   case 8:
 245     ((TH1F*)fContainer->At(kChargeDeposit))->Draw("c");
 246     break;
 247   case 9:
 248     h=(TH1F*)fContainer->At(kPurity);
 249     h->SetBarOffset(.2);
 250     h->SetBarWidth(.6);
 251     h->SetFillColor(kGreen);
 252     h->Draw("bar1");
 253     break;
 254   default:
 255     ((TH1F*)fContainer->At(kNTracksEventHist))->Draw("pl");
 256     break;
 257   }
 258 }
 259
 260 //_______________________________________________________
 261 TObjArray *AliTRDcheckDetector::Histos(){
 262   //
 263   // Create QA histograms
 264   //
 265   if(fContainer) return fContainer;
 266
 267   fContainer = new TObjArray(25);
 268   // Register Histograms
 269   fContainer->AddAt(new TH1F("hNtrks", "Number of Tracks per event", 100, 0, 100), kNTracksEventHist);
 270   fContainer->AddAt(new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100), kNEventsTriggerTracks);
 271   fContainer->AddAt(new TH1F("hNcls", "Nr. of clusters per track", 181, -0.5, 180.5), kNclustersHist);
 272   fContainer->AddAt(new TH1F("hNtls", "Nr. tracklets per track", 7, -0.5, 6.5), kNtrackletsHist);
 273   fContainer->AddAt(new TH1F("hNtlsFindable", "Ratio of found/findable Tracklets" , 11, -0.05, 1.05), kNTrackletsVsFindable);
 274   fContainer->AddAt(new TH1F("hNclTls","Mean Number of clusters per tracklet", 31, -0.5, 30.5), kNclusterTrackletHist);
 275   fContainer->AddAt(new TH1F("hChi2", "Chi2", 200, 0, 20), kChi2);
 276   fContainer->AddAt(new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5), kChi2Normalized);
 277   fContainer->AddAt(new TH1F("hSM", "Track Counts in Supermodule", 18, -0.5, 17.5), kNTracksSectorHist);
 278   // Detector signal on Detector-by-Detector basis
 279   fContainer->AddAt(new TProfile("hPHdetector", "Average PH", 31, -0.5, 30.5), kPulseHeight);
 280   fContainer->AddAt(new TH1F("hQclDetector", "Cluster charge", 200, 0, 1200), kClusterCharge);
 281   fContainer->AddAt(new TH1F("hQTdetector", "Total Charge Deposit", 6000, 0, 6000), kChargeDeposit);
 282   fContainer->AddAt(new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100), kNEventsTrigger);
 283
 284   return fContainer;
 285 }
 286
 287 /*
 288 * Plotting Functions
 289 */
 290
 291 //_______________________________________________________
 292 TH1 *AliTRDcheckDetector::PlotMeanNClusters(const AliTRDtrackV1 *track){
 293   //
 294   // Plot the mean number of clusters per tracklet
 295   //
 296   if(track) fTrack = track;
 297   if(!fTrack){
 298     AliWarning("No Track defined.");
 299     return 0x0;
 300   }
 301   TH1 *h = 0x0;
 302   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclusterTrackletHist)))){
 303     AliWarning("No Histogram defined.");
 304     return 0x0;
 305   }
 306   AliTRDseedV1 *tracklet = 0x0;
 307   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 308     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 309     h->Fill(tracklet->GetN());
 310   }
 311   return h;
 312 }
 313
 314 //_______________________________________________________
 315 TH1 *AliTRDcheckDetector::PlotNClusters(const AliTRDtrackV1 *track){
 316   //
 317   // Plot the number of clusters in one track
 318   //
 319   if(track) fTrack = track;
 320   if(!fTrack){
 321     AliWarning("No Track defined.");
 322     return 0x0;
 323   }
 324   TH1 *h = 0x0;
 325   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersHist)))){
 326     AliWarning("No Histogram defined.");
 327     return 0x0;
 328   }
 329
 330   Int_t nclusters = 0;
 331   AliTRDseedV1 *tracklet = 0x0;
 332   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 333     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 334     nclusters += tracklet->GetN();
 335     if(fDebugLevel > 2){
 336       Int_t crossing = tracklet->GetNChange();
 337       AliTRDcluster *c = 0x0;
 338       for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 339         if(!(c = tracklet->GetClusters(itime))) continue;
 340         break;
 341       }
 342       Int_t detector = c->GetDetector();
 343       Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 344       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 345       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 346       Float_t momentum = 0.;
 347       Int_t pdg = 0;
 348       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 349       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 350       if(fMC){
 351         if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 352         pdg = fMC->GetPDG();
 353       }
 354       (*fDebugStream) << "NClusters"
 355         << "Detector="  << detector
 356         << "Sector="    << sector
 357         << "crossing="  << crossing
 358         << "momentum="  << momentum
 359         << "pdg="                               << pdg
 360         << "theta="                     << theta
 361         << "phi="                               << phi
 362         << "kinkIndex=" << kinkIndex
 363         << "TPCncls="           << TPCncls
 364         << "nclusters=" << nclusters
 365         << "\n";
 366     }
 367   }
 368   h->Fill(nclusters);
 369   return h;
 370 }
 371
 372 //_______________________________________________________
 373 TH1 *AliTRDcheckDetector::PlotChi2(const AliTRDtrackV1 *track){
 374   //
 375   // Plot the chi2 of the track
 376   //
 377   if(track) fTrack = track;
 378   if(!fTrack){
 379     AliWarning("No Track defined.");
 380     return 0x0;
 381   }
 382   TH1 *h = 0x0;
 383   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2)))){
 384     AliWarning("No Histogram defined.");
 385     return 0x0;
 386   }
 387   h->Fill(fTrack->GetChi2());
 388   return h;
 389 }
 390
 391 //_______________________________________________________
 392 TH1 *AliTRDcheckDetector::PlotNormalizedChi2(const AliTRDtrackV1 *track){
 393   //
 394   // Plot the chi2 of the track
 395   //
 396   if(track) fTrack = track;
 397   if(!fTrack){
 398     AliWarning("No Track defined.");
 399     return 0x0;
 400   }
 401   TH1 *h = 0x0;
 402   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChi2Normalized)))){
 403     AliWarning("No Histogram defined.");
 404     return 0x0;
 405   }
 406   Int_t nTracklets = 0;
 407   AliTRDseedV1 *tracklet = 0x0;
 408   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 409     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 410     nTracklets++;
 411   }
 412   h->Fill(fTrack->GetChi2()/nTracklets);
 413   return h;
 414 }
 415
 416
 417 //_______________________________________________________
 418 TH1 *AliTRDcheckDetector::PlotNTracklets(const AliTRDtrackV1 *track){
 419   //
 420   // Plot the number of tracklets
 421   //
 422   if(track) fTrack = track;
 423   if(!fTrack){
 424     AliWarning("No Track defined.");
 425     return 0x0;
 426   }
 427   TH1 *h = 0x0;
 428   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsHist)))){
 429     AliWarning("No Histogram defined.");
 430     return 0x0;
 431   }
 432   Int_t nTracklets = GetNTracklets(track);
 433   h->Fill(nTracklets);
 434   if(fDebugLevel > 3){
 435     if(nTracklets == 1){
 436       // If we have one Tracklet, check in which layer this happens
 437       Int_t layer = -1;
 438       AliTRDseedV1 *tracklet = 0x0;
 439       for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 440         if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()){layer =  il; break;}
 441       }
 442       (*fDebugStream) << "PlotNTracklets"
 443         << "Layer=" << layer
 444         << "\n";
 445     }
 446   }
 447   return h;
 448 }
 449
 450 //_______________________________________________________
 451 TH1 *AliTRDcheckDetector::PlotTrackletsVsFindable(const AliTRDtrackV1 *track){
 452   //
 453   // Plots the ratio of number of tracklets vs.
 454   // number of findable tracklets
 455   //
 456   // Findable tracklets are defined as track prolongation
 457   // to layer i does not hit the dead area +- epsilon
 458   //
 459   const Float_t epsilon = 0.01;   // dead area tolerance
 460   const Float_t epsilon_R = 1;
 461   Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
 462
 463   if(track) fTrack = track;
 464   if(!fTrack){
 465     AliWarning("No Track defined.");
 466     return 0x0;
 467   }
 468   TH1 *h = 0x0;
 469   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTrackletsVsFindable)))){
 470     AliWarning("No Histogram defined.");
 471     return 0x0;
 472   }
 473   Int_t nFound = 0, nFindable = 0;
 474   Int_t stack = -1;
 475   Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
 476   Double_t y = 0., z = 0.;
 477   AliTRDseedV1 *tracklet = 0x0;
 478   AliTRDpadPlane *pp;
 479   AliTRDgeometry *fGeo = new AliTRDgeometry;
 480   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 481     if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
 482       tracklet->SetReconstructor(fReconstructor);
 483       nFound++;
 484       // printout tracklet
 485     }
 486   }
 487   // 2 Different cases:
 488   // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
 489   // 2nd barrel track: here we propagate the track to the layers
 490   AliTrackPoint points[6];
 491   Float_t xyz[3];
 492   memset(xyz, 0, sizeof(Float_t) * 3);
 493   if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
 494     // stand alone track
 495     for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 496       xyz[0] = x_anode[il];
 497       points[il].SetXYZ(xyz);
 498     }
 499     AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
 500   } else {
 501     // barrel track
 502     for(Int_t il = AliTRDgeometry::kNlayer; il--;){
 503       // The track points have to be in reverse order for the Kalman Filter (since we move back)
 504       xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
 505       points[il].SetXYZ(xyz);
 506     }
 507     for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
 508       printf("%d. X = %f\n", ipt, points[ipt].GetX());
 509     AliTRDtrackerV1::FitKalman(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kFALSE, 6, points);
 510     AliTrackPoint tempPointCont[AliTRDgeometry::kNlayer];
 511     memcpy(tempPointCont, points, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
 512     for(Int_t il = AliTRDgeometry::kNlayer; il--; ){
 513       tempPointCont[il].GetXYZ(xyz);
 514       points[AliTRDgeometry::kNlayer - il -1].SetXYZ(xyz);
 515     }
 516   }
 517   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 518     y = points[il].GetY();
 519     z = points[il].GetZ();
 520     if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
 521     pp = fGeo->GetPadPlane(il, stack);
 522     ymin = pp->GetCol0() + epsilon;
 523     ymax = pp->GetColEnd() - epsilon;
 524     zmin = pp->GetRowEnd() + epsilon;
 525     zmax = pp->GetRow0() - epsilon;
 526     // ignore y-crossing (material)
 527     if((z > zmin && z < zmax) && (y > ymin && y < ymax)) nFindable++;
 528 /*
 529       if(fDebugLevel > 3){
 530         Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetMeanz() : 0};
 531         Int_t hasTracklet = tracklet ? 1 : 0;
 532         (*fDebugStream)   << "GetFindableTracklets"
 533           << "layer="     << il
 534           << "ytracklet=" << pos_tracklet[0]
 535           << "ytrack="    << y
 536           << "ztracklet=" << pos_tracklet[1]
 537           << "ztrack="    << z
 538           << "tracklet="  << hasTracklet
 539           << "\n";
 540       }
 541 */
 542   }
 543
 544   delete fGeo;
 545   h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
 546   if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
 547   return h;
 548 }
 549
 550 //_______________________________________________________
 551 TH1 *AliTRDcheckDetector::PlotPulseHeight(const AliTRDtrackV1 *track){
 552   //
 553   // Plot the average pulse height
 554   //
 555   if(track) fTrack = track;
 556   if(!fTrack){
 557     AliWarning("No Track defined.");
 558     return 0x0;
 559   }
 560   TProfile *h = 0x0;
 561   if(!(h = dynamic_cast<TProfile *>(fContainer->At(kPulseHeight)))){
 562     AliWarning("No Histogram defined.");
 563     return 0x0;
 564   }
 565   AliTRDseedV1 *tracklet = 0x0;
 566   AliTRDcluster *c = 0x0;
 567   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 568     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 569     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 570       if(!(c = tracklet->GetClusters(itime))) continue;
 571       Int_t localtime        = c->GetLocalTimeBin();
 572       Double_t absolute_charge = TMath::Abs(c->GetQ());
 573       h->Fill(localtime, absolute_charge);
 574       if(fDebugLevel > 3){
 575         Int_t crossing = tracklet->GetNChange();
 576         Int_t detector = c->GetDetector();
 577         Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 578         Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 579         Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 580         Float_t momentum = 0.;
 581         Int_t pdg = 0;
 582         Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 583         UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 584         if(fMC){
 585           if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 586           pdg = fMC->GetPDG();
 587         }
 588         (*fDebugStream) << "PulseHeight"
 589           << "Detector="        << detector
 590           << "Sector="          << sector
 591           << "crossing="        << crossing
 592           << "Timebin="         << localtime
 593           << "Charge="          << absolute_charge
 594           << "momentum="        << momentum
 595           << "pdg="                             << pdg
 596           << "theta="                   << theta
 597           << "phi="                             << phi
 598           << "kinkIndex="       << kinkIndex
 599           << "TPCncls="         << TPCncls
 600           << "\n";
 601       }
 602     }
 603   }
 604   return h;
 605 }
 606
 607 //_______________________________________________________
 608 TH1 *AliTRDcheckDetector::PlotClusterCharge(const AliTRDtrackV1 *track){
 609   //
 610   // Plot the cluster charge
 611   //
 612   if(track) fTrack = track;
 613   if(!fTrack){
 614     AliWarning("No Track defined.");
 615     return 0x0;
 616   }
 617   TH1 *h = 0x0;
 618   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kClusterCharge)))){
 619     AliWarning("No Histogram defined.");
 620     return 0x0;
 621   }
 622   AliTRDseedV1 *tracklet = 0x0;
 623   AliTRDcluster *c = 0x0;
 624   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 625     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
 626     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 627       if(!(c = tracklet->GetClusters(itime))) continue;
 628       h->Fill(c->GetQ());
 629     }
 630   }
 631   return h;
 632 }
 633
 634 //_______________________________________________________
 635 TH1 *AliTRDcheckDetector::PlotChargeDeposit(const AliTRDtrackV1 *track){
 636   //
 637   // Plot the charge deposit per chamber
 638   //
 639   if(track) fTrack = track;
 640   if(!fTrack){
 641     AliWarning("No Track defined.");
 642     return 0x0;
 643   }
 644   TH1 *h = 0x0;
 645   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeDeposit)))){
 646     AliWarning("No Histogram defined.");
 647     return 0x0;
 648   }
 649   AliTRDseedV1 *tracklet = 0x0;
 650   AliTRDcluster *c = 0x0, *c1 = 0x0;    // c1 for the Debug Stream
 651   Double_t Qtot = 0;
 652   Int_t nTracklets = 0;
 653   if(fDebugLevel > 3)
 654     nTracklets = GetNTracklets(fTrack); // fill NTracklet to the Debug Stream
 655   for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
 656     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 657     Qtot = 0;
 658     c1 = 0x0;
 659     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 660       if(!(c = tracklet->GetClusters(itime))) continue;
 661       if(!c1) c1 = c;
 662       Qtot += TMath::Abs(c->GetQ());
 663     }
 664     h->Fill(Qtot);
 665     if(fDebugLevel > 3){
 666       Int_t crossing = tracklet->GetNChange();
 667       Int_t detector = c1->GetDetector();
 668       Float_t sector = static_cast<Int_t>(detector/AliTRDgeometry::kNdets);
 669       Float_t theta = TMath::ATan(tracklet->GetZfit(1));
 670       Float_t phi = TMath::ATan(tracklet->GetYfit(1));
 671       Float_t momentum = 0.;
 672       Int_t pdg = 0;
 673       Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
 674       UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
 675       if(fMC){
 676               if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
 677         pdg = fMC->GetPDG();
 678       }
 679       (*fDebugStream) << "ChargeDeposit"
 680         << "Detector="  << detector
 681         << "Sector="    << sector
 682         << "crossing="  << crossing
 683         << "momentum="  << momentum
 684         << "nTracklets="<< nTracklets
 685         << "pdg="                               << pdg
 686         << "theta="                     << theta
 687         << "phi="                               << phi
 688         << "kinkIndex=" << kinkIndex
 689         << "TPCncls="           << TPCncls
 690         << "QT="        << Qtot
 691         << "\n";
 692     }
 693   }
 694   return h;
 695 }
 696
 697 //_______________________________________________________
 698 TH1 *AliTRDcheckDetector::PlotTracksSector(const AliTRDtrackV1 *track){
 699   //
 700   // Plot the number of tracks per Sector
 701   //
 702   if(track) fTrack = track;
 703   if(!fTrack){
 704     AliWarning("No Track defined.");
 705     return 0x0;
 706   }
 707   TH1 *h = 0x0;
 708   if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTracksSectorHist)))){
 709     AliWarning("No Histogram defined.");
 710     return 0x0;
 711   }
 712   AliTRDseedV1 *tracklet = 0x0;
 713   AliTRDcluster *c = 0x0;
 714   Int_t sector = -1;
 715   for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
 716     if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
 717     for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
 718       if(!(c = tracklet->GetClusters(itime))) continue;
 719       sector = static_cast<Int_t>(c->GetDetector()/AliTRDgeometry::kNdets);
 720     }
 721     break;
 722   }
 723   h->Fill(sector);
 724   return h;
 725 }
 726
 727 //_______________________________________________________
 728 Int_t AliTRDcheckDetector::GetNTracklets(const AliTRDtrackV1 *track){
 729   //
 730   // Count the number of tracklets per track
 731   //
 732   if(!track) return 0;
 733   Int_t nTracklets = 0;
 734   AliTRDseedV1 *tracklet = 0x0;
 735   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 736     if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()) nTracklets++;
 737   }
 738   return nTracklets;
 739 }
 740
 741 //________________________________________________________
 742 void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
 743 {
 744
 745   fReconstructor->SetRecoParam(r);
 746 }