]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/qaRec/AliTRDcheckDetector.cxx
! MAJOR ! update in the TRD tracking
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDcheckDetector.cxx
CommitLineData
2cdfacc5 1#include <TAxis.h>
2#include <TCanvas.h>
95cda7cf 3#include <TFile.h>
23abf4db 4#include <TH1F.h>
2cdfacc5 5#include <TGaxis.h>
abc70aaf 6#include <TGraph.h>
23abf4db 7#include <TMath.h>
2a4a428a 8#include <TMap.h>
23abf4db 9#include <TObjArray.h>
2a4a428a 10#include <TObject.h>
11#include <TObjString.h>
2cdfacc5 12#include <TPad.h>
23abf4db 13#include <TProfile.h>
9e548ae7 14#include <TProfile2D.h>
95cda7cf 15#include <TROOT.h>
23abf4db 16
107fde80 17#include "AliLog.h"
23abf4db 18#include "AliTRDcluster.h"
a24151d1 19#include "AliESDHeader.h"
20#include "AliESDRun.h"
3cfaffa4 21#include "AliESDtrack.h"
2374fb85 22#include "AliTRDgeometry.h"
3cfaffa4 23#include "AliTRDpadPlane.h"
2cdfacc5 24#include "AliTRDSimParam.h"
23abf4db 25#include "AliTRDseedV1.h"
26#include "AliTRDtrackV1.h"
3cfaffa4 27#include "AliTRDtrackerV1.h"
28#include "AliTRDReconstructor.h"
4e488f26 29#include "AliTrackReference.h"
3cfaffa4 30#include "AliTrackPointArray.h"
31#include "AliTracker.h"
23abf4db 32#include "TTreeStream.h"
33
34#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
a24151d1 35#include "AliTRDtrackInfo/AliTRDeventInfo.h"
23abf4db 36#include "AliTRDcheckDetector.h"
37
95cda7cf 38#include <cstdio>
56f8e2b1 39#include <iostream>
95cda7cf 40
23abf4db 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//_______________________________________________________
54AliTRDcheckDetector::AliTRDcheckDetector():
93e41bce 55 AliTRDrecoTask("DetChecker", "Basic Detector Checker")
a24151d1 56 ,fEventInfo(0x0)
2a4a428a 57 ,fTriggerNames(0x0)
3cfaffa4 58 ,fReconstructor(0x0)
7dc3c50c 59 ,fGeo(0x0)
23abf4db 60{
2b468513 61 //
62 // Default constructor
63 //
a24151d1 64 DefineInput(1,AliTRDeventInfo::Class());
3cfaffa4 65 fReconstructor = new AliTRDReconstructor;
66 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
7dc3c50c 67 fGeo = new AliTRDgeometry;
107fde80 68 InitFunctorList();
23abf4db 69}
70
71//_______________________________________________________
72AliTRDcheckDetector::~AliTRDcheckDetector(){
2b468513 73 //
74 // Destructor
75 //
2a4a428a 76 if(fTriggerNames) delete fTriggerNames;
3cfaffa4 77 delete fReconstructor;
7dc3c50c 78 delete fGeo;
23abf4db 79}
80
a24151d1 81//_______________________________________________________
82void AliTRDcheckDetector::ConnectInputData(Option_t *opt){
a391a274 83 //
84 // Connect the Input data with the task
85 //
86 AliTRDrecoTask::ConnectInputData(opt);
87 fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
a24151d1 88}
89
23abf4db 90//_______________________________________________________
91void AliTRDcheckDetector::CreateOutputObjects(){
2b468513 92 //
93 // Create Output Objects
94 //
2374fb85 95 OpenFile(0,"RECREATE");
107fde80 96 fContainer = Histos();
2374fb85 97 if(!fTriggerNames) fTriggerNames = new TMap();
23abf4db 98}
99
100//_______________________________________________________
107fde80 101void AliTRDcheckDetector::Exec(Option_t *opt){
2b468513 102 //
103 // Execution function
104 // Filling TRD quality histos
105 //
9e548ae7 106 if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
107fde80 107 AliTRDrecoTask::Exec(opt);
2b468513 108 Int_t nTracks = 0; // Count the number of tracks per event
a24151d1 109 Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
afefec95 110 TString triggername = fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
9e548ae7 111 if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
b1957d3c 112 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
2b468513 113 for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
107fde80 114 if(!fTracks->UncheckedAt(iti)) continue;
115 AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
116 if(!fTrackInfo->GetTrack()) continue;
2b468513 117 nTracks++;
118 }
959f550d 119 if(nTracks){
b1957d3c 120 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
121 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
959f550d 122 }
9e548ae7 123 if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
2a4a428a 124 fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
9e548ae7 125 // also set the label for both histograms
b1957d3c 126 TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
9e548ae7 127 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
b1957d3c 128 histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
9e548ae7 129 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
130 }
2b468513 131 PostData(0, fContainer);
23abf4db 132}
133
134//_______________________________________________________
135void AliTRDcheckDetector::Terminate(Option_t *){
2b468513 136 //
137 // Terminate function
138 //
23abf4db 139}
140
95cda7cf 141//_______________________________________________________
142Bool_t AliTRDcheckDetector::PostProcess(){
a391a274 143 //
144 // Do Postprocessing (for the moment set the number of Reference histograms)
145 //
146
b1957d3c 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());
a391a274 156
b1957d3c 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);
a391a274 161 Float_t purities[20], val = 0;
162 TString triggernames[20];
163 Int_t nTriggerClasses = 0;
b1957d3c 164 for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
165 if((val = h1->GetBinContent(ibin))){
a391a274 166 purities[nTriggerClasses] = val;
b1957d3c 167 triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
a391a274 168 nTriggerClasses++;
169 }
170 }
b1957d3c 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
7dc3c50c 180 fNRefFigures = 11;
b1957d3c 181
a391a274 182 return kTRUE;
95cda7cf 183}
184
185//_______________________________________________________
e15179be 186Bool_t AliTRDcheckDetector::GetRefFigure(Int_t ifig){
a391a274 187 //
188 // Setting Reference Figures
189 //
b1957d3c 190 TObjArray *arr = 0x0;
2cdfacc5 191 TH1 *h = 0x0, *h1 = 0x0, *h2 = 0x0;
192 TGaxis *axis = 0x0;
a391a274 193 switch(ifig){
b1957d3c 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);
a391a274 202 if(!h->GetEntries()) break;
a391a274 203 h->SetFillColor(kGreen);
204 h->SetBarOffset(.2);
205 h->SetBarWidth(.6);
206 h->Draw("bar1");
b1957d3c 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);
7dc3c50c 218 if(!h->GetEntries()) break;
219 h->Scale(100./h->Integral());
7dc3c50c 220 h->SetFillColor(kGreen);
221 h->SetBarOffset(.2);
222 h->SetBarWidth(.6);
223 h->Draw("bar1");
b1957d3c 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);
a391a274 230 if(!h->GetEntries()) break;
a391a274 231 h->SetFillColor(kGreen);
232 h->SetBarOffset(.2);
233 h->SetBarWidth(.6);
234 h->Draw("bar1");
b1957d3c 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);
a391a274 242 h->SetMarkerStyle(24);
2cdfacc5 243 h->SetMarkerColor(kBlack);
244 h->SetLineColor(kBlack);
a391a274 245 h->Draw("e1");
2cdfacc5 246 // copy the second histogram in a new one with the same x-dimension as the phs with respect to time
b1957d3c 247 h1 = (TH1F *)arr->At(1);
2cdfacc5 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);
b1957d3c 265 axis->SetTitle("x_{0}-x_{c} [cm]");
2cdfacc5 266 axis->Draw();
b1957d3c 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);
a391a274 282 h->SetBarOffset(.2);
283 h->SetBarWidth(.6);
284 h->SetFillColor(kGreen);
285 h->Draw("bar1");
286 break;
287 default:
a391a274 288 break;
289 }
b1957d3c 290 AliInfo(Form("Reference plot [%d] missing result", ifig));
291 return kFALSE;
95cda7cf 292}
293
107fde80 294//_______________________________________________________
295TObjArray *AliTRDcheckDetector::Histos(){
296 //
297 // Create QA histograms
298 //
299 if(fContainer) return fContainer;
300
b1957d3c 301 fContainer = new TObjArray(14);
94f34623 302 //fContainer->SetOwner(kTRUE);
303
107fde80 304 // Register Histograms
b1957d3c 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);
107fde80 416
417 return fContainer;
418}
419
420/*
421* Plotting Functions
422*/
423
424//_______________________________________________________
52e4836c 425TH1 *AliTRDcheckDetector::PlotNClustersTracklet(const AliTRDtrackV1 *track){
107fde80 426 //
427 // Plot the mean number of clusters per tracklet
428 //
74b2e03d 429 if(track) fTrack = track;
107fde80 430 if(!fTrack){
74b2e03d 431 AliWarning("No Track defined.");
432 return 0x0;
107fde80 433 }
434 TH1 *h = 0x0;
b1957d3c 435 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
107fde80 436 AliWarning("No Histogram defined.");
437 return 0x0;
438 }
439 AliTRDseedV1 *tracklet = 0x0;
2374fb85 440 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
107fde80 441 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
52e4836c 442 h->Fill(tracklet->GetN2());
107fde80 443 }
444 return h;
445}
446
447//_______________________________________________________
52e4836c 448TH1 *AliTRDcheckDetector::PlotNClustersTrack(const AliTRDtrackV1 *track){
107fde80 449 //
450 // Plot the number of clusters in one track
451 //
74b2e03d 452 if(track) fTrack = track;
107fde80 453 if(!fTrack){
74b2e03d 454 AliWarning("No Track defined.");
455 return 0x0;
107fde80 456 }
457 TH1 *h = 0x0;
b1957d3c 458 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
107fde80 459 AliWarning("No Histogram defined.");
460 return 0x0;
461 }
2374fb85 462
107fde80 463 Int_t nclusters = 0;
464 AliTRDseedV1 *tracklet = 0x0;
2374fb85 465 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
107fde80 466 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
467 nclusters += tracklet->GetN();
468 if(fDebugLevel > 2){
52e4836c 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));
107fde80 473 Float_t momentum = 0.;
474 Int_t pdg = 0;
52a79f8d 475 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
476 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
107fde80 477 if(fMC){
52a79f8d 478 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
75d00a6d 479 pdg = fMC->GetPDG();
107fde80 480 }
b1957d3c 481 (*fDebugStream) << "NClustersTrack"
107fde80 482 << "Detector=" << detector
107fde80 483 << "crossing=" << crossing
484 << "momentum=" << momentum
485 << "pdg=" << pdg
486 << "theta=" << theta
487 << "phi=" << phi
52a79f8d 488 << "kinkIndex=" << kinkIndex
489 << "TPCncls=" << TPCncls
107fde80 490 << "nclusters=" << nclusters
491 << "\n";
492 }
493 }
494 h->Fill(nclusters);
495 return h;
496}
497
107fde80 498
499//_______________________________________________________
b1957d3c 500TH1 *AliTRDcheckDetector::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
107fde80 501 //
b1957d3c 502 // Plot the number of tracklets
107fde80 503 //
74b2e03d 504 if(track) fTrack = track;
505 if(!fTrack){
506 AliWarning("No Track defined.");
507 return 0x0;
107fde80 508 }
509 TH1 *h = 0x0;
b1957d3c 510 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
107fde80 511 AliWarning("No Histogram defined.");
512 return 0x0;
513 }
b1957d3c 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 }
2374fb85 528 }
107fde80 529 return h;
530}
531
532
533//_______________________________________________________
b1957d3c 534TH1 *AliTRDcheckDetector::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
107fde80 535 //
536 // Plot the number of tracklets
537 //
74b2e03d 538 if(track) fTrack = track;
107fde80 539 if(!fTrack){
74b2e03d 540 AliWarning("No Track defined.");
541 return 0x0;
107fde80 542 }
543 TH1 *h = 0x0;
b1957d3c 544 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
107fde80 545 AliWarning("No Histogram defined.");
546 return 0x0;
547 }
b1957d3c 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++;
3cfaffa4 555 }
b1957d3c 556 h->Fill(ncross);
3cfaffa4 557 return h;
558}
559
560//_______________________________________________________
b1957d3c 561TH1 *AliTRDcheckDetector::PlotFindableTracklets(const AliTRDtrackV1 *track){
3cfaffa4 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 //
7dc3c50c 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 //
3cfaffa4 584 const Float_t epsilon = 0.01; // dead area tolerance
7dc3c50c 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)
3cfaffa4 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;
b1957d3c 596 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
3cfaffa4 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.;
107fde80 604 AliTRDseedV1 *tracklet = 0x0;
3cfaffa4 605 AliTRDpadPlane *pp;
3cfaffa4 606 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
607 if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
608 tracklet->SetReconstructor(fReconstructor);
609 nFound++;
3cfaffa4 610 }
107fde80 611 }
3cfaffa4 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
7dc3c50c 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];
3cfaffa4 633 for(Int_t il = AliTRDgeometry::kNlayer; il--;){
7dc3c50c 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
3cfaffa4 638 xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
7dc3c50c 639 points_inward[il].SetXYZ(xyz);
640 xyz[0] = x_anode[il] + epsilon_R;
641 points_outward[il].SetXYZ(xyz);
3cfaffa4 642 }
7dc3c50c 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);
3cfaffa4 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)
7dc3c50c 662 if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
3cfaffa4 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;
b1957d3c 666 (*fDebugStream) << "FindableTracklets"
3cfaffa4 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 }
3cfaffa4 675 }
676
3cfaffa4 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)));
107fde80 679 return h;
680}
681
b1957d3c 682
683//_______________________________________________________
684TH1 *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//_______________________________________________________
703TH1 *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
107fde80 728//_______________________________________________________
b1957d3c 729TH1 *AliTRDcheckDetector::PlotPHt(const AliTRDtrackV1 *track){
107fde80 730 //
731 // Plot the average pulse height
732 //
74b2e03d 733 if(track) fTrack = track;
107fde80 734 if(!fTrack){
74b2e03d 735 AliWarning("No Track defined.");
736 return 0x0;
107fde80 737 }
738 TProfile *h = 0x0;
b1957d3c 739 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
107fde80 740 AliWarning("No Histogram defined.");
741 return 0x0;
742 }
107fde80 743 AliTRDseedV1 *tracklet = 0x0;
744 AliTRDcluster *c = 0x0;
2374fb85 745 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
107fde80 746 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
52e4836c 747 Int_t crossing = Int_t(tracklet->IsRowCross());
748 Int_t detector = tracklet->GetDetector();
3cfaffa4 749 for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
107fde80 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){
2cdfacc5 755 Double_t distance[2];
756 GetDistanceToTracklet(distance, tracklet, c);
52e4836c 757 Float_t theta = TMath::ATan(tracklet->GetZref(1));
758 Float_t phi = TMath::ATan(tracklet->GetYref(1));
107fde80 759 Float_t momentum = 0.;
760 Int_t pdg = 0;
52a79f8d 761 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
762 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
107fde80 763 if(fMC){
75d00a6d 764 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
765 pdg = fMC->GetPDG();
107fde80 766 }
b1957d3c 767 (*fDebugStream) << "PHt"
107fde80 768 << "Detector=" << detector
2374fb85 769 << "crossing=" << crossing
107fde80 770 << "Timebin=" << localtime
771 << "Charge=" << absolute_charge
772 << "momentum=" << momentum
773 << "pdg=" << pdg
774 << "theta=" << theta
775 << "phi=" << phi
52a79f8d 776 << "kinkIndex=" << kinkIndex
777 << "TPCncls=" << TPCncls
2cdfacc5 778 << "dy=" << distance[0]
779 << "dz=" << distance[1]
780 << "c.=" << c
107fde80 781 << "\n";
782 }
783 }
784 }
785 return h;
786}
787
2cdfacc5 788//_______________________________________________________
b1957d3c 789TH1 *AliTRDcheckDetector::PlotPHx(const AliTRDtrackV1 *track){
2cdfacc5 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;
b1957d3c 800 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
2cdfacc5 801 AliWarning("No Histogram defined.");
802 return 0x0;
803 }
a589a812 804 Float_t offset = .5*AliTRDgeometry::CamHght();
2cdfacc5 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
107fde80 819//_______________________________________________________
b1957d3c 820TH1 *AliTRDcheckDetector::PlotChargeCluster(const AliTRDtrackV1 *track){
107fde80 821 //
822 // Plot the cluster charge
823 //
74b2e03d 824 if(track) fTrack = track;
107fde80 825 if(!fTrack){
74b2e03d 826 AliWarning("No Track defined.");
827 return 0x0;
107fde80 828 }
829 TH1 *h = 0x0;
b1957d3c 830 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
107fde80 831 AliWarning("No Histogram defined.");
832 return 0x0;
833 }
107fde80 834 AliTRDseedV1 *tracklet = 0x0;
835 AliTRDcluster *c = 0x0;
2374fb85 836 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
107fde80 837 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
3cfaffa4 838 for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
107fde80 839 if(!(c = tracklet->GetClusters(itime))) continue;
840 h->Fill(c->GetQ());
841 }
842 }
843 return h;
844}
845
846//_______________________________________________________
b1957d3c 847TH1 *AliTRDcheckDetector::PlotChargeTracklet(const AliTRDtrackV1 *track){
107fde80 848 //
849 // Plot the charge deposit per chamber
850 //
74b2e03d 851 if(track) fTrack = track;
107fde80 852 if(!fTrack){
74b2e03d 853 AliWarning("No Track defined.");
854 return 0x0;
107fde80 855 }
856 TH1 *h = 0x0;
b1957d3c 857 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
107fde80 858 AliWarning("No Histogram defined.");
859 return 0x0;
860 }
107fde80 861 AliTRDseedV1 *tracklet = 0x0;
52e4836c 862 AliTRDcluster *c = 0x0;
107fde80 863 Double_t Qtot = 0;
b1957d3c 864 Int_t nTracklets =fTrack->GetNumberOfTracklets();
2374fb85 865 for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
107fde80 866 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
52e4836c 867 Qtot = 0.;
868 for(Int_t ic = AliTRDseed::knTimebins; ic--;){
869 if(!(c = tracklet->GetClusters(ic))) continue;
107fde80 870 Qtot += TMath::Abs(c->GetQ());
871 }
872 h->Fill(Qtot);
873 if(fDebugLevel > 3){
52e4836c 874 Int_t crossing = (Int_t)tracklet->IsRowCross();
875 Int_t detector = tracklet->GetDetector();
107fde80 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;
52a79f8d 880 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
881 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
107fde80 882 if(fMC){
75d00a6d 883 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
884 pdg = fMC->GetPDG();
107fde80 885 }
b1957d3c 886 (*fDebugStream) << "ChargeTracklet"
107fde80 887 << "Detector=" << detector
107fde80 888 << "crossing=" << crossing
889 << "momentum=" << momentum
3cfaffa4 890 << "nTracklets="<< nTracklets
107fde80 891 << "pdg=" << pdg
892 << "theta=" << theta
893 << "phi=" << phi
52a79f8d 894 << "kinkIndex=" << kinkIndex
895 << "TPCncls=" << TPCncls
107fde80 896 << "QT=" << Qtot
897 << "\n";
898 }
899 }
900 return h;
901}
902
903//_______________________________________________________
b1957d3c 904TH1 *AliTRDcheckDetector::PlotNTracksSector(const AliTRDtrackV1 *track){
107fde80 905 //
906 // Plot the number of tracks per Sector
907 //
74b2e03d 908 if(track) fTrack = track;
107fde80 909 if(!fTrack){
74b2e03d 910 AliWarning("No Track defined.");
911 return 0x0;
107fde80 912 }
913 TH1 *h = 0x0;
b1957d3c 914 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
107fde80 915 AliWarning("No Histogram defined.");
916 return 0x0;
917 }
52e4836c 918
919 // TODO we should compare with
920 // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
921
107fde80 922 AliTRDseedV1 *tracklet = 0x0;
107fde80 923 Int_t sector = -1;
2374fb85 924 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
925 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
52e4836c 926 sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
107fde80 927 break;
928 }
929 h->Fill(sector);
930 return h;
931}
932
3cfaffa4 933
934//________________________________________________________
935void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
936{
937
938 fReconstructor->SetRecoParam(r);
939}
2cdfacc5 940
941//________________________________________________________
b1957d3c 942void 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);
2cdfacc5 947}