]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TRD/AliTRDcheckDET.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
macro to run standalone ESD filtering
[u/mrichter/AliRoot.git] / PWGPP / TRD / AliTRDcheckDET.cxx
CommitLineData
1ee39b3a 1////////////////////////////////////////////////////////////////////////////
2// //
3// //
4// Basic checks for tracking and detector performance //
5//
6// For the moment (15.10.2009) the following checks are implemented //
7// - Number of clusters/track
8// - Number of clusters/tracklet
9// - Number of tracklets/track from different sources (Barrel, StandAlone)
10// - Number of findable tracklets
11// - Number of tracks per event and TRD sector
12// - <PH>
13// - Chi2 distribution for tracks
14// - Charge distribution per cluster and tracklet
15// - Number of events and tracks per trigger source
16// - Trigger purity
17// - Track and Tracklet propagation status
18//
19// Authors: //
20// Anton Andronic <A.Andronic@gsi.de> //
21// Alexandru Bercuci <A.Bercuci@gsi.de> //
22// Markus Fasel <M.Fasel@gsi.de> //
23// //
24////////////////////////////////////////////////////////////////////////////
25
8a3dd7a0 26#include <TArrayD.h>
1ee39b3a 27#include <TAxis.h>
28#include <TCanvas.h>
29#include <TFile.h>
30#include <TH1F.h>
31#include <TH1I.h>
a5d9fe6f 32#include <TH2F.h>
2e8eb67c 33#include <TH3S.h>
34#include <TH3F.h>
1ee39b3a 35#include <TF1.h>
36#include <TGaxis.h>
37#include <TGraph.h>
147c3968 38#include <TGraphErrors.h>
1ee39b3a 39#include <TLegend.h>
3907f080 40#include <TLinearFitter.h>
1ee39b3a 41#include <TMath.h>
42#include <TMap.h>
44fd53ca 43#include <TProfile2D.h>
1ee39b3a 44#include <TObjArray.h>
45#include <TObject.h>
46#include <TObjString.h>
47
48#include <TPad.h>
49#include <TProfile.h>
50#include <TProfile2D.h>
51#include <TROOT.h>
f8f46e4d 52#include <TChain.h>
1ee39b3a 53
54#include "AliLog.h"
55#include "AliTRDcluster.h"
56#include "AliESDHeader.h"
57#include "AliESDRun.h"
58#include "AliESDtrack.h"
a5d9fe6f 59#include "AliExternalTrackParam.h"
1ee39b3a 60#include "AliTRDgeometry.h"
61#include "AliTRDpadPlane.h"
62#include "AliTRDSimParam.h"
63#include "AliTRDseedV1.h"
64#include "AliTRDtrackV1.h"
65#include "AliTRDtrackerV1.h"
66#include "AliTRDReconstructor.h"
67#include "AliTrackReference.h"
68#include "AliTrackPointArray.h"
69#include "AliTracker.h"
70#include "TTreeStream.h"
71
72#include "info/AliTRDtrackInfo.h"
73#include "info/AliTRDeventInfo.h"
fd7ffd88 74#include "AliTRDinfoGen.h"
1ee39b3a 75#include "AliTRDcheckDET.h"
2bfe5463 76#include "AliTRDpwgppHelper.h"
1ee39b3a 77
78#include <cstdio>
79#include <iostream>
80
81ClassImp(AliTRDcheckDET)
82
8a3dd7a0 83const Color_t AliTRDcheckDET::fkColorsCentrality[AliTRDeventInfo::kCentralityClasses] = {kTeal, kOrange, kGreen, kBlue ,kRed};
84
1ee39b3a 85//_______________________________________________________
86AliTRDcheckDET::AliTRDcheckDET():
f2e89a4c 87 AliTRDrecoTask()
2e8eb67c 88 ,fCentralityClass(-1)
db99a57a 89 ,fTriggerNames(NULL)
c732f879 90 ,fFlags(0)
1ee39b3a 91{
92 //
93 // Default constructor
94 //
4fa7d600 95 SetNameTitle("TRDcheckDET", "Basic TRD data checker");
f8f46e4d 96}
97
705f8b0a 98//_______________________________________________________
f8f46e4d 99AliTRDcheckDET::AliTRDcheckDET(char* name):
100 AliTRDrecoTask(name, "Basic TRD data checker")
2e8eb67c 101 ,fCentralityClass(-1)
db99a57a 102 ,fTriggerNames(NULL)
f8f46e4d 103 ,fFlags(0)
104{
105 //
106 // Default constructor
107 //
1ee39b3a 108 InitFunctorList();
109}
110
f8f46e4d 111
1ee39b3a 112//_______________________________________________________
113AliTRDcheckDET::~AliTRDcheckDET(){
114 //
115 // Destructor
116 //
117 if(fTriggerNames) delete fTriggerNames;
1ee39b3a 118}
119
1ee39b3a 120//_______________________________________________________
f8f46e4d 121void AliTRDcheckDET::UserCreateOutputObjects(){
1ee39b3a 122 //
123 // Create Output Objects
124 //
068e2c00 125 AliTRDrecoTask::UserCreateOutputObjects();
1ee39b3a 126 if(!fTriggerNames) fTriggerNames = new TMap();
127}
128
129//_______________________________________________________
f8f46e4d 130void AliTRDcheckDET::UserExec(Option_t *opt){
1ee39b3a 131 //
132 // Execution function
133 // Filling TRD quality histos
134 //
2e8eb67c 135 AliDebug(2, Form("EventInfo[%p] Header[%p]", (void*)fEvent, (void*)(fEvent?fEvent->GetEventHeader():NULL)));
136 if(fEvent) fCentralityClass = fEvent->GetCentrality();
137 else fCentralityClass = -1; // Assume pp
138
b4414720 139 AliTRDrecoTask::UserExec(opt);
b2f4ab8d 140
7fe4e88b 141 TH1F *histo(NULL); AliTRDtrackInfo *fTrackInfo(NULL); Int_t nTracks(0); // Count the number of tracks per event
1ee39b3a 142 for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
143 if(!fTracks->UncheckedAt(iti)) continue;
7fe4e88b 144 if(!(fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti)))) continue;
1ee39b3a 145 if(!fTrackInfo->GetTrack()) continue;
146 nTracks++;
147 }
b2f4ab8d 148 if(nTracks)
3ecd3327 149 if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent)))) histo->Fill(nTracks);
f8f46e4d 150
34b8a3ce 151 if(!fEvent->GetEventHeader()) return; // For trigger statistics event header is essential
2e8eb67c 152 Int_t triggermask = fEvent->GetEventHeader()->GetTriggerMask();
153 TString triggername = fEvent->GetRunInfo()->GetFiredTriggerClasses(triggermask);
b2f4ab8d 154 AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
7fe4e88b 155 if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger)))) histo->Fill(triggermask);
b2f4ab8d 156
7fe4e88b 157 if(nTracks){
158 if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks)))) histo->Fill(triggermask);
159 }
1ee39b3a 160 if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
161 fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
162 // also set the label for both histograms
7fe4e88b 163 if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))))
164 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
165 if((histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))))
166 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
1ee39b3a 167 }
1ee39b3a 168}
169
170
8a3dd7a0 171//_______________________________________________________
1ee39b3a 172//_______________________________________________________
173Bool_t AliTRDcheckDET::PostProcess(){
174 //
175 // Do Postprocessing (for the moment set the number of Reference histograms)
176 //
177
7fe4e88b 178 TH1F *h(NULL), *h1(NULL);
179
1ee39b3a 180 // Calculate of the trigger clusters purity
7fe4e88b 181 if((h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"))) &&
182 (h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks")))) {
183 h1->Divide(h);
184 Float_t purities[20], val = 0; memset(purities, 0, 20*sizeof(Float_t));
185 TString triggernames[20];
186 Int_t nTriggerClasses = 0;
187 for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
188 if((val = h1->GetBinContent(ibin))){
189 purities[nTriggerClasses] = val;
190 triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
191 nTriggerClasses++;
192 }
193 }
194
195 if((h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity)))){
196 TAxis *ax = h->GetXaxis();
197 for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
198 h->Fill(itrg, purities[itrg]);
199 ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
200 }
201 ax->SetRangeUser(-0.5, nTriggerClasses+.5);
202 h->GetYaxis()->SetRangeUser(0,1);
1ee39b3a 203 }
204 }
1ee39b3a 205
206 // track status
7fe4e88b 207 if((h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus)))){
208 Float_t ok = h->GetBinContent(1);
209 Int_t nerr = h->GetNbinsX();
210 for(Int_t ierr=nerr; ierr--;){
211 h->SetBinContent(ierr+1, ok>0.?1.e2*h->GetBinContent(ierr+1)/ok:0.);
212 }
213 h->SetBinContent(1, 0.);
1ee39b3a 214 }
1ee39b3a 215 // tracklet status
f8f46e4d 216
7fe4e88b 217 TObjArray *arr(NULL);
218 if(( arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus)))){
219 for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
220 if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
221 Float_t okB = h->GetBinContent(1);
222 Int_t nerrB = h->GetNbinsX();
223 for(Int_t ierr=nerrB; ierr--;){
224 h->SetBinContent(ierr+1, okB>0.?1.e2*h->GetBinContent(ierr+1)/okB:0.);
225 }
226 h->SetBinContent(1, 0.);
1ee39b3a 227 }
1ee39b3a 228 }
0c76cfa4 229 fNRefFigures = 17;
1ee39b3a 230
231 return kTRUE;
232}
233
b795b480 234void AliTRDcheckDET::MakeSummary(){
235 //
236 // Create summary plots for TRD check DET
237 // This function creates 2 summary plots:
238 // - General Quantities
239 // - PHS
240 // The function will reuse GetRefFigure
241 //
242
243 TCanvas *cOut = new TCanvas(Form("summary%s1", GetName()), Form("Summary 1 for task %s", GetName()), 1024, 768);
244 cOut->Divide(3,3);
245
246 // Create figures using GetRefFigure
247 cOut->cd(1); GetRefFigure(kFigNtracksEvent);
248 cOut->cd(2); GetRefFigure(kFigNtracksSector);
249 cOut->cd(3); GetRefFigure(kFigNclustersTrack);
250 cOut->cd(4); GetRefFigure(kFigNclustersTracklet);
251 cOut->cd(5); GetRefFigure(kFigNtrackletsTrack);
252 cOut->cd(6); GetRefFigure(kFigNTrackletsP);
253 cOut->cd(7); GetRefFigure(kFigChargeCluster);
254 cOut->cd(8); GetRefFigure(kFigChargeTracklet);
e77999ac 255 cOut->SaveAs("TRD_TrackQuality.gif");
b795b480 256 delete cOut;
257
258 // Second Plot: PHS
259 cOut = new TCanvas(Form("summary%s2", GetName()), Form("Summary 2 for task %s", GetName()), 1024, 512);
260 cOut->cd(); GetRefFigure(kFigPH);
e77999ac 261 cOut->SaveAs("TRD_PH.gif");
b795b480 262 delete cOut;
44fd53ca 263
264 // Third Plot: Mean Number of clusters as function of eta, phi and layer
265 cOut = new TCanvas(Form("summary%s3", GetName()), Form("Summary 3 for task %s", GetName()), 1024, 768);
266 cOut->cd(); MakePlotMeanClustersLayer();
e77999ac 267 cOut->SaveAs("TRD_MeanNclusters.gif");
44fd53ca 268 delete cOut;
269
b795b480 270}
271
1ee39b3a 272//_______________________________________________________
273Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
8a3dd7a0 274 //
275 // Setting Reference Figures
276 //
277 gPad->SetLogy(0);
278 gPad->SetLogx(0);
279 TH1 *h = NULL; TObjArray *arr=NULL;
280 TLegend *leg = NULL;
281 Bool_t kFIRST(1);
282 switch(ifig){
283 case kFigNclustersTrack:
284 MakePlotNclustersTrack();
285 return kTRUE;
286 case kFigNclustersTracklet:
287 MakePlotNclustersTracklet();
288 return kTRUE;
289 case kFigNtrackletsTrack:
290 h=MakePlotNTracklets();
291 if(h){
292 PutTrendValue("NTrackletsTrack", h->GetMean());
293 PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
294 }
295 return kTRUE;
296 case kFigNTrackletsP:
297 MakePlotnTrackletsVsP();
298 return kTRUE;
299 case kFigNtrackletsCross:
300 h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsCross"), -1);
301 if(!MakeBarPlot(h, kRed)) break;
302 PutTrendValue("NTrackletsCross", h->GetMean());
303 PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
304 return kTRUE;
305 case kFigNtrackletsFindable:
306 h = ProjectCentrality((TH2*)fContainer->FindObject("hNtlsFindable"), -1);
307 if(!MakeBarPlot(h, kGreen)) break;
308 PutTrendValue("NTrackletsFindable", h->GetMean());
309 PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
310 return kTRUE;
311 case kFigNtracksEvent:
312 (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
313 PutTrendValue("NTracksEvent", h->GetMean());
314 PutTrendValue("NTracksEventRMS", h->GetRMS());
315 return kTRUE;
316 case kFigNtracksSector:
317 h = (TH1F*)fContainer->FindObject("hNtrksSector");
318 if(!MakeBarPlot(h, kGreen)) break;
319 PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
320 return kTRUE;
321 case kFigTrackStatus:
322 if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
323 h->GetXaxis()->SetRangeUser(0.5, -1);
324 h->GetYaxis()->CenterTitle();
325 h->Draw("c");
326 PutTrendValue("TrackStatus", h->Integral());
327 gPad->SetLogy(0);
328 return kTRUE;
329 case kFigTrackletStatus:
330 if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
331 leg = new TLegend(.68, .7, .97, .97);
332 leg->SetBorderSize(0);leg->SetFillStyle(0);
333 leg->SetHeader("TRD layer");
334 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
335 if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
336 if(kFIRST){
337 h->Draw("pl");
338 h->GetXaxis()->SetRangeUser(0.5, -1);
339 h->GetYaxis()->CenterTitle();
340 kFIRST = kFALSE;
341 } else h->Draw("samepl");
342 leg->AddEntry(h, Form("ly = %d", ily), "l");
343 PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
344 }
345 leg->Draw();
346 gPad->SetLogy(0);
347 return kTRUE;
348 case kFigChi2:
349 MakePlotChi2();
350 return kTRUE;
351 case kFigPH:
352 gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
353 MakePlotPulseHeight();
354 gPad->SetLogy(0);
355 return kTRUE;
356 case kFigChargeCluster:
357 h = ProjectCentrality((TH2*)fContainer->FindObject("hQcl"), -1);
358 h->Draw("c");
359 gPad->SetLogy(1);
360 PutTrendValue("ChargeCluster", h->GetMaximumBin());
361 PutTrendValue("ChargeClusterRMS", h->GetRMS());
362 return kTRUE;
363 case kFigChargeTracklet:
364 MakePlotTrackletCharge();
365 return kTRUE;
366 case kFigNeventsTrigger:
367 ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
368 return kTRUE;
369 case kFigNeventsTriggerTracks:
370 ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
371 return kTRUE;
372 case kFigTriggerPurity:
373 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
374 break;
375 default:
376 break;
377 }
378 AliInfo(Form("Reference plot [%d] missing result", ifig));
379 return kFALSE;
1ee39b3a 380}
381
382//_______________________________________________________
383TObjArray *AliTRDcheckDET::Histos(){
384 //
385 // Create QA histograms
386 //
f8f46e4d 387
1ee39b3a 388 if(fContainer) return fContainer;
389
390 fContainer = new TObjArray(20);
246fe7f7 391 fContainer->SetOwner(kTRUE);
1ee39b3a 392
393 // Register Histograms
394 TH1 * h = NULL;
2e8eb67c 395 TH2 * h2 = NULL; // Pointer for two dimensional histograms
396 TH3 * h3 = NULL; // Pointer for tree dimensional histograms
1ee39b3a 397 TAxis *ax = NULL;
2e8eb67c 398 if(!(h2 = (TH2F *)gROOT->FindObject("hNcls"))){
399 h2 = new TH2F("hNcls", "N_{clusters}/track;N_{clusters};Centrality;Entries", 181, -0.5, 180.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
400 } else h2->Reset();
401 fContainer->AddAt(h2, kNclustersTrack);
402
44fd53ca 403 TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
404 arr->SetOwner(kTRUE); arr->SetName("clusters");
405 fContainer->AddAt(arr, kNclustersLayer);
406 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
407 if(!(h = (TProfile2D *)gROOT->FindObject(Form("hNcl%d", ily)))){
408 h = new TProfile2D(Form("hNcl%d", ily), Form("Mean Number of clusters in Layer %d", ily), 100, -1.0, 1.0, 50, -1.1*TMath::Pi(), 1.1*TMath::Pi());
409 h->GetXaxis()->SetTitle("#eta");
410 h->GetYaxis()->SetTitle("#phi");
411 } else h->Reset();
412 arr->AddAt(h, ily);
413 }
414
2e8eb67c 415 if(!(h2 = (TH2F *)gROOT->FindObject("hNclTls"))){
416 h2 = new TH2F("hNclTls","N_{clusters}/tracklet;N_{clusters};Entries", 51, -0.5, 50.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
417 } else h2->Reset();
418 fContainer->AddAt(h2, kNclustersTracklet);
1ee39b3a 419
2e8eb67c 420 if(!(h2 = (TH2F *)gROOT->FindObject("hNtls"))){
421 h2 = new TH2F("hNtls", "N_{tracklets}/track;N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
422 } else h2->Reset();
423 fContainer->AddAt(h2, kNtrackletsTrack);
1ee39b3a 424
2e8eb67c 425 if(!(h = (TH2F *)gROOT->FindObject("htlsSTA"))){
426 h = new TH2F("hNtlsSTA", "#splitline{N_{tracklets}/track}{Stand Alone};N^{tracklet};Centrality;freq.[%]", AliTRDgeometry::kNlayer, 0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
1ee39b3a 427 }
428 fContainer->AddAt(h, kNtrackletsSTA);
429
a5d9fe6f 430 // Binning for momentum dependent tracklet Plots
de791ef5 431 const Int_t kNp(30);
efb9a759 432 Float_t p=0.2;
2e8eb67c 433 Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1], binsCent[AliTRDeventInfo::kCentralityClasses+2];
efb9a759 434 for(Int_t i=0;i<kNp+1; i++,p+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=p;
147c3968 435 for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
2e8eb67c 436 for(Int_t icent = -1; icent < AliTRDeventInfo::kCentralityClasses + 1; icent++) binsCent[icent+1] = icent - 0.5;
437 if(!(h3 = (TH3F *)gROOT->FindObject("htlsBAR"))){
a5d9fe6f 438 // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
2e8eb67c 439 h3 = new TH3F("hNtlsBAR",
440 "N_{tracklets}/track;p [GeV/c];N^{tracklet};freq. [%]",
441 kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt, AliTRDeventInfo::kCentralityClasses + 1, binsCent);
1ee39b3a 442 }
2e8eb67c 443 fContainer->AddAt(h3, kNtrackletsBAR);
1ee39b3a 444
445 //
2e8eb67c 446 if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsCross"))){
447 h2 = new TH2F("hNtlsCross", "N_{tracklets}^{cross}/track;n_{row cross};Centrality;freq.[%]", 7, -0.5, 6.5, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
448 } else h2->Reset();
449 fContainer->AddAt(h2, kNtrackletsCross);
1ee39b3a 450
2e8eb67c 451 if(!(h2 = (TH2F *)gROOT->FindObject("hNtlsFindable"))){
452 h2 = new TH2F("hNtlsFindable", "Found/Findable Tracklets;r[a.u];Centrality;Entries" , 101, -0.005, 1.005, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
453 } else h2->Reset();
454 fContainer->AddAt(h2, kNtrackletsFindable);
1ee39b3a 455
456 if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
457 h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
458 h->GetXaxis()->SetTitle("N_{tracks}");
459 h->GetYaxis()->SetTitle("Entries");
460 } else h->Reset();
461 fContainer->AddAt(h, kNtracksEvent);
462
463 if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
464 h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
465 h->GetXaxis()->SetTitle("sector");
466 h->GetYaxis()->SetTitle("freq. [%]");
467 } else h->Reset();
468 fContainer->AddAt(h, kNtracksSector);
469
470 if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
471 const Int_t nerr = 7;
472 h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
473 const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
474 ax = h->GetXaxis();
475 for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
476 h->SetYTitle("Relative Error to Good [%]");
477 }
478 fContainer->AddAt(h, kTrackStatus);
479
44fd53ca 480 arr = new TObjArray(AliTRDgeometry::kNlayer);
1ee39b3a 481 arr->SetOwner(kTRUE); arr->SetName("TrackletStatus");
482 fContainer->AddAt(arr, kTrackletStatus);
483 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
484 if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
485 const Int_t nerr = 8;
486 h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
487 h->SetLineColor(ily+1);
488 const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
489 ax = h->GetXaxis();
490 for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
491 h->SetYTitle("Relative Error to Good [%]");
492 } else h->Reset();
493 arr->AddAt(h, ily);
494 }
495
496 // <PH> histos
c0a6ac06 497 arr = new TObjArray(3);
1ee39b3a 498 arr->SetOwner(kTRUE); arr->SetName("<PH>");
499 fContainer->AddAt(arr, kPH);
2e8eb67c 500 if(!(h3 = (TH3F *)gROOT->FindObject("hPHx"))){
40600e64 501 h3 = new TH3F("hPHx", "<PH>(x);x[mm];Charge[a.u.];Centrality", 31, -0.08, 4.88, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
2e8eb67c 502 } else h3->Reset();
503 arr->AddAt(h3, 0);
504 if(!(h3 = (TH3F *)gROOT->FindObject("hPHt"))){
40600e64 505 h3 = new TH3F("hPHt", "<PH>(t);time[100ns];Charge[a.u.];Centrality", 31, -0.5, 30.5, 100, 0, 1024, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
2e8eb67c 506 } else h3->Reset();
507 arr->AddAt(h3, 1);
1ee39b3a 508
509 // Chi2 histos
2e8eb67c 510 if(!(h3 = (TH3S*)gROOT->FindObject("hChi2"))){
511 h3 = new TH3S("hChi2", "#chi^{2}/track;ndf;#chi^{2}/ndf;Centrality", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
512 } else h3->Reset();
513 fContainer->AddAt(h3, kChi2);
1ee39b3a 514
2e8eb67c 515 if(!(h2 = (TH2F *)gROOT->FindObject("hQcl"))){
516 h2 = new TH2F("hQcl", "Q_{cluster};Charge[a.u.];Centrality;Entries", 200, 0, 1200, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
517 }else h2->Reset();
518 fContainer->AddAt(h2, kChargeCluster);
1ee39b3a 519
2e8eb67c 520 if(!(h2 = (TH2F *)gROOT->FindObject("hQtrklt"))){
521 h2 = new TH2F("hQtrklt", "Q_{tracklet};Charge[a.u.];Centrality;Entries", 6000, 0, 6000, AliTRDeventInfo::kCentralityClasses + 1, -1.5, AliTRDeventInfo::kCentralityClasses - 0.5);
522 }else h2->Reset();
523 fContainer->AddAt(h2, kChargeTracklet);
1ee39b3a 524
525
526 if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
527 h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
528 else h->Reset();
f8f46e4d 529
1ee39b3a 530 fContainer->AddAt(h, kNeventsTrigger);
531
532 if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
533 h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
534 else h->Reset();
535 fContainer->AddAt(h, kNeventsTriggerTracks);
536
537 if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
538 h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
539 h->GetXaxis()->SetTitle("Trigger Cluster");
540 h->GetYaxis()->SetTitle("freq.");
541 } else h->Reset();
542 fContainer->AddAt(h, kTriggerPurity);
543
544 return fContainer;
545}
546
8a3dd7a0 547//_______________________________________________________
548TH1 *AliTRDcheckDET::ProjectCentrality(TH2 *hIn, Int_t centralityBin){
549 //
550 // Project histogram to a given centrality Bin
551 //
552 if(!hIn) return NULL;
553 if(centralityBin >= AliTRDeventInfo::kCentralityClasses) centralityBin = -1;
554 Int_t binMin = centralityBin > -1 ? centralityBin + 2 : 0,
555 binMax = centralityBin > -1 ? centralityBin + 2 : -1;
556 return hIn->ProjectionX(Form("%s_%d", hIn->GetName(), centralityBin), binMin, binMax);
557}
558
1ee39b3a 559/*
560* Plotting Functions
561*/
562
563//_______________________________________________________
564TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
565{
566//
567// Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
568// PROL - track prolongation failure
569// PROP - track propagation failure
570// AJST - crossing sectors failure
571// SNP - too large bending
572// TINI - tracklet initialization failure
573// UPDT - track position/covariance update failure
574//
575// Performance plot looks as below:
576//Begin_Html
577//<img src="TRD/trackStatus.gif">
578//End_Html
579//
580 if(track) fkTrack = track;
581 if(!fkTrack){
4e7a049c 582 AliDebug(4, "No Track defined.");
db99a57a 583 return NULL;
1ee39b3a 584 }
db99a57a 585 TH1 *h = NULL;
1ee39b3a 586 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
587 AliWarning("No Histogram defined.");
db99a57a 588 return NULL;
1ee39b3a 589 }
590 h->Fill(fkTrack->GetStatusTRD());
591 return h;
592}
593
594//_______________________________________________________
595TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
596{
597//
598// Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
599// Geom -
600// Bound - tracklet too close to chamber walls
601// NoCl - no clusters in the track roads
602// NoAttach - fail to attach clusters
603// NoClTr - fail to use clusters for fit
604// NoFit - tracklet fit failled
605// Chi2 - chi2 tracklet-track over threshold
606//
607// Performance plot looks as below:
608//Begin_Html
609//<img src="TRD/trackletStatus.gif">
610//End_Html
611//
612 if(track) fkTrack = track;
613 if(!fkTrack){
4e7a049c 614 AliDebug(4, "No Track defined.");
db99a57a 615 return NULL;
1ee39b3a 616 }
db99a57a 617 TObjArray *arr =NULL;
1ee39b3a 618 if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
619 AliWarning("Histograms not defined.");
db99a57a 620 return NULL;
1ee39b3a 621 }
622
db99a57a 623 TH1 *h = NULL;
1ee39b3a 624 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
625 if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
626 AliWarning(Form("Missing histo for layer %d.", ily));
627 continue;
628 }
629 h->Fill(fkTrack->GetStatusTRD(ily));
630 }
631 return h;
632}
633
634//_______________________________________________________
635TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
636 //
637 // Plot the mean number of clusters per tracklet
638 //
639 if(track) fkTrack = track;
640 if(!fkTrack){
4e7a049c 641 AliDebug(4, "No Track defined.");
db99a57a 642 return NULL;
1ee39b3a 643 }
44fd53ca 644 AliExternalTrackParam *par = fkTrack->GetTrackIn() ? fkTrack->GetTrackIn() : fkTrack->GetTrackOut();
2e8eb67c 645 TH2 *h = NULL;
44fd53ca 646 TProfile2D *hlayer = NULL;
647 Double_t eta = 0., phi = 0.;
2e8eb67c 648 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTracklet)))){
1ee39b3a 649 AliWarning("No Histogram defined.");
db99a57a 650 return NULL;
1ee39b3a 651 }
db99a57a 652 AliTRDseedV1 *tracklet = NULL;
44fd53ca 653 TObjArray *histosLayer = dynamic_cast<TObjArray *>(fContainer->At(kNclustersLayer));
654 if(!histosLayer){
655 AliWarning("No Histograms for single layer defined");
656 }
1ee39b3a 657 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
658 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
2e8eb67c 659 h->Fill(tracklet->GetN2(), fCentralityClass);
44fd53ca 660 if(histosLayer && par){
661 if((hlayer = dynamic_cast<TProfile2D *>(histosLayer->At(itl)))){
662 GetEtaPhiAt(par, tracklet->GetX0(), eta, phi);
663 hlayer->Fill(eta, phi, tracklet->GetN2());
664 }
665 }
1ee39b3a 666 }
667 return h;
668}
669
670//_______________________________________________________
671TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
672 //
673 // Plot the number of clusters in one track
674 //
675 if(track) fkTrack = track;
676 if(!fkTrack){
4e7a049c 677 AliDebug(4, "No Track defined.");
db99a57a 678 return NULL;
1ee39b3a 679 }
2e8eb67c 680 TH2 *h = NULL;
681 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNclustersTrack)))){
1ee39b3a 682 AliWarning("No Histogram defined.");
db99a57a 683 return NULL;
1ee39b3a 684 }
685
686 Int_t nclusters = 0;
db99a57a 687 AliTRDseedV1 *tracklet = NULL;
a310e49b 688 AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
9653acd3 689 if(!par) return NULL;
0c76cfa4 690 Double_t momentumRec = par->P();
1ee39b3a 691 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
692 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
fc2ec42c 693 Int_t n(tracklet->GetN());
694 nclusters += n;
1ee39b3a 695 if(DebugLevel() > 2){
696 Int_t crossing = Int_t(tracklet->IsRowCross());
697 Int_t detector = tracklet->GetDetector();
698 Float_t theta = TMath::ATan(tracklet->GetZref(1));
699 Float_t phi = TMath::ATan(tracklet->GetYref(1));
0c76cfa4 700 Float_t momentumMC = 0.;
1ee39b3a 701 Int_t pdg = 0;
702 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
703 UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
704 if(fkMC){
0c76cfa4 705 if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
1ee39b3a 706 pdg = fkMC->GetPDG();
707 }
708 (*DebugStream()) << "NClustersTrack"
709 << "Detector=" << detector
2e8eb67c 710 << "Centrality="<< fCentralityClass
1ee39b3a 711 << "crossing=" << crossing
0c76cfa4 712 << "momentumMC="<< momentumMC
2e8eb67c 713 << "momentumRec="<< momentumRec
1ee39b3a 714 << "pdg=" << pdg
715 << "theta=" << theta
716 << "phi=" << phi
717 << "kinkIndex=" << kinkIndex
718 << "TPCncls=" << nclsTPC
fc2ec42c 719 << "TRDncls=" << n
1ee39b3a 720 << "\n";
721 }
722 }
2e8eb67c 723 h->Fill(nclusters, fCentralityClass);
1ee39b3a 724 return h;
725}
726
727
728//_______________________________________________________
729TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
730 //
731 // Plot the number of tracklets
732 //
733 if(track) fkTrack = track;
734 if(!fkTrack){
4e7a049c 735 AliDebug(4, "No Track defined.");
db99a57a 736 return NULL;
1ee39b3a 737 }
2e8eb67c 738 TH2 *h = NULL, *hSta = NULL; TH3 *hBarrel = NULL;
739 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsTrack)))){
1ee39b3a 740 AliWarning("No Histogram defined.");
db99a57a 741 return NULL;
1ee39b3a 742 }
743 Int_t nTracklets = fkTrack->GetNumberOfTracklets();
2e8eb67c 744 h->Fill(nTracklets, fCentralityClass);
1ee39b3a 745 if(!fkESD) return h;
746 Int_t status = fkESD->GetStatus();
a5d9fe6f 747
1ee39b3a 748/* printf("in/out/refit/pid: TRD[%d|%d|%d|%d]\n", status &AliESDtrack::kTRDin ? 1 : 0, status &AliESDtrack::kTRDout ? 1 : 0, status &AliESDtrack::kTRDrefit ? 1 : 0, status &AliESDtrack::kTRDpid ? 1 : 0);*/
a5d9fe6f 749 Double_t p = 0.;
0c76cfa4 750 Int_t method = -1; // to distinguish between stand alone and full barrel tracks in the debugging
1ee39b3a 751 if((status & AliESDtrack::kTRDin) != 0){
0c76cfa4 752 method = 1;
a5d9fe6f 753 // Full Barrel Track: Save momentum dependence
2e8eb67c 754 if(!(hBarrel = dynamic_cast<TH3F *>(fContainer->At(kNtrackletsBAR)))){
1ee39b3a 755 AliWarning("Method: Barrel. Histogram not processed!");
147c3968 756 return NULL;
a5d9fe6f 757 }
147c3968 758 AliExternalTrackParam *par(fkTrack->GetTrackIn());
759 if(!par){
760 AliError("Input track params missing");
761 return NULL;
762 }
b795b480 763 p = par->P(); // p needed later in the debug streaming
2e8eb67c 764 hBarrel->Fill(p, nTracklets, fCentralityClass);
1ee39b3a 765 } else {
a5d9fe6f 766 // Stand alone Track: momentum dependence not usefull
0c76cfa4 767 method = 0;
2e8eb67c 768 if(!(hSta = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsSTA)))) {
1ee39b3a 769 AliWarning("Method: StandAlone. Histogram not processed!");
147c3968 770 return NULL;
a5d9fe6f 771 }
2e8eb67c 772 hSta->Fill(nTracklets, fCentralityClass);
1ee39b3a 773 }
1ee39b3a 774
0c76cfa4 775 if(DebugLevel() > 2){
776 AliTRDseedV1 *tracklet = NULL;
fd7ffd88 777 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
b795b480 778 Int_t sector = -1, stack = -1, detector;
0c76cfa4 779 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
780 if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
b795b480 781 detector = tracklet->GetDetector();
fd7ffd88 782 sector = geo->GetSector(detector);
783 stack = geo->GetStack(detector);
0c76cfa4 784 break;
785 }
786 (*DebugStream()) << "NTrackletsTrack"
787 << "Sector=" << sector
2e8eb67c 788 << "Stack=" << stack
789 << "Centrality=" << fCentralityClass
0c76cfa4 790 << "NTracklets=" << nTracklets
791 << "Method=" << method
792 << "p=" << p
793 << "\n";
794 }
1ee39b3a 795 if(DebugLevel() > 3){
b795b480 796 AliTRDseedV1 *tracklet = NULL;
797 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
798 if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
0c76cfa4 799 (*DebugStream()) << "NTrackletsLayer"
b795b480 800 << "Layer=" << il
801 << "p=" << p
802 << "\n";
a5d9fe6f 803 }
1ee39b3a 804 }
805 }
806 return h;
807}
808
809
810//_______________________________________________________
811TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
812 //
813 // Plot the number of tracklets
814 //
815 if(track) fkTrack = track;
816 if(!fkTrack){
4e7a049c 817 AliDebug(4, "No Track defined.");
db99a57a 818 return NULL;
1ee39b3a 819 }
2e8eb67c 820 TH2 *h = NULL;
821 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsCross)))){
1ee39b3a 822 AliWarning("No Histogram defined.");
db99a57a 823 return NULL;
1ee39b3a 824 }
825
826 Int_t ncross = 0;
db99a57a 827 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 828 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
829 if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
830
831 if(tracklet->IsRowCross()) ncross++;
832 }
2e8eb67c 833 h->Fill(ncross, fCentralityClass);
1ee39b3a 834 return h;
835}
836
837//_______________________________________________________
838TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
839 //
840 // Plots the ratio of number of tracklets vs.
841 // number of findable tracklets
842 //
843 // Findable tracklets are defined as track prolongation
844 // to layer i does not hit the dead area +- epsilon
845 //
846 // In order to check whether tracklet hist active area in Layer i,
847 // the track is refitted and the fitted position + an uncertainty
848 // range is compared to the chamber border (also with a different
849 // uncertainty)
850 //
851 // For the track fit two cases are distinguished:
852 // If the track is a stand alone track (defined by the status bit
853 // encoding, then the track is fitted with the tilted Rieman model
854 // Otherwise the track is fitted with the Kalman fitter in two steps:
855 // Since the track parameters are give at the outer point, we first
856 // fit in direction inwards. Afterwards we fit again in direction outwards
857 // to extrapolate the track to layers which are not reached by the track
858 // For the Kalman model, the radial track points have to be shifted by
859 // a distance epsilon in the direction that we want to fit
860 //
861 const Float_t epsilon = 0.01; // dead area tolerance
862 const Float_t epsilonR = 1; // shift in radial direction of the anode wire position (Kalman filter only)
863 const Float_t deltaY = 0.7; // Tolerance in the track position in y-direction
864 const Float_t deltaZ = 7.0; // Tolerance in the track position in z-direction (Padlength)
865 Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
866
867 if(track) fkTrack = track;
868 if(!fkTrack){
4e7a049c 869 AliDebug(4, "No Track defined.");
db99a57a 870 return NULL;
1ee39b3a 871 }
2e8eb67c 872 TH2 *h = NULL;
873 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsFindable)))){
1ee39b3a 874 AliWarning("No Histogram defined.");
db99a57a 875 return NULL;
1ee39b3a 876 }
877 Int_t nFound = 0, nFindable = 0;
878 Int_t stack = -1;
879 Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
880 Double_t y = 0., z = 0.;
db99a57a 881 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 882 AliTRDpadPlane *pp;
883 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
884 if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
fd7ffd88 885 tracklet->SetReconstructor(AliTRDinfoGen::Reconstructor());
1ee39b3a 886 nFound++;
887 }
888 }
889 // 2 Different cases:
890 // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
891 // 2nd barrel track: here we propagate the track to the layers
892 AliTrackPoint points[6];
893 Float_t xyz[3];
894 memset(xyz, 0, sizeof(Float_t) * 3);
2cef552b 895 if(fkESD && ((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
1ee39b3a 896 // stand alone track
897 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
898 xyz[0] = xAnode[il];
899 points[il].SetXYZ(xyz);
900 }
db99a57a 901 AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
1ee39b3a 902 } else {
903 // barrel track
904 //
905 // 2 Steps:
906 // -> Kalman inwards
907 // -> Kalman outwards
908 AliTRDtrackV1 copyTrack(*fkTrack); // Do Kalman on a (non-constant) copy of the track
909 AliTrackPoint pointsInward[6], pointsOutward[6];
910 for(Int_t il = AliTRDgeometry::kNlayer; il--;){
911 // In order to avoid complications in the Kalman filter if the track points have the same radial
912 // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
913 // in the direction we want to go
914 // The track points have to be in reverse order for the Kalman Filter inwards
915 xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
916 pointsInward[il].SetXYZ(xyz);
917 xyz[0] = xAnode[il] + epsilonR;
918 pointsOutward[il].SetXYZ(xyz);
919 }
920 /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
921 printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
922 // Kalman inwards
db99a57a 923 AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
1ee39b3a 924 memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
925 // Kalman outwards
db99a57a 926 AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
1ee39b3a 927 memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
928 }
fd7ffd88 929 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
1ee39b3a 930 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
931 y = points[il].GetY();
932 z = points[il].GetZ();
fd7ffd88 933 if((stack = geo->GetStack(z, il)) < 0) continue; // Not findable
934 pp = geo->GetPadPlane(il, stack);
1ee39b3a 935 ymin = pp->GetCol0() + epsilon;
936 ymax = pp->GetColEnd() - epsilon;
937 zmin = pp->GetRowEnd() + epsilon;
938 zmax = pp->GetRow0() - epsilon;
939 // ignore y-crossing (material)
940 if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
941 if(DebugLevel() > 3){
942 Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
943 Int_t hasTracklet = tracklet ? 1 : 0;
944 (*DebugStream()) << "FindableTracklets"
945 << "layer=" << il
946 << "ytracklet=" << posTracklet[0]
947 << "ytrack=" << y
948 << "ztracklet=" << posTracklet[1]
949 << "ztrack=" << z
950 << "tracklet=" << hasTracklet
951 << "\n";
952 }
953 }
954
2e8eb67c 955 h->Fill((nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1), fCentralityClass);
1ee39b3a 956 AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
957 return h;
958}
959
960
961//_______________________________________________________
962TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
963 //
964 // Plot the chi2 of the track
965 //
966 if(track) fkTrack = track;
967 if(!fkTrack){
4e7a049c 968 AliDebug(4, "No Track defined.");
db99a57a 969 return NULL;
1ee39b3a 970 }
2e8eb67c 971 TH3 *h = NULL;
972 if(!(h = dynamic_cast<TH3S*>(fContainer->At(kChi2)))) {
1ee39b3a 973 AliWarning("No Histogram defined.");
db99a57a 974 return NULL;
1ee39b3a 975 }
976 Int_t n = fkTrack->GetNumberOfTracklets();
db99a57a 977 if(!n) return NULL;
1ee39b3a 978
2e8eb67c 979 h->Fill(n, fkTrack->GetChi2()/n, fCentralityClass);
1ee39b3a 980 return h;
981}
982
983
984//_______________________________________________________
985TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
986 //
987 // Plot the average pulse height
988 //
4e7a049c 989 if(track) fkTrack = track;
990 if(!fkTrack){
991 AliDebug(4, "No Track defined.");
db99a57a 992 return NULL;
1ee39b3a 993 }
2e8eb67c 994 TH3F *h = NULL;
995 if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
1ee39b3a 996 AliWarning("No Histogram defined.");
db99a57a 997 return NULL;
1ee39b3a 998 }
db99a57a 999 AliTRDseedV1 *tracklet = NULL;
1000 AliTRDcluster *c = NULL;
1ee39b3a 1001 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1002 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1003 Int_t crossing = Int_t(tracklet->IsRowCross());
1004 Int_t detector = tracklet->GetDetector();
1005 tracklet->ResetClusterIter();
1006 while((c = tracklet->NextCluster())){
c732f879 1007 if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1ee39b3a 1008 Int_t localtime = c->GetLocalTimeBin();
1009 Double_t absoluteCharge = TMath::Abs(c->GetQ());
2e8eb67c 1010 h->Fill(localtime, absoluteCharge, fCentralityClass);
1ee39b3a 1011 if(DebugLevel() > 3){
b795b480 1012 Int_t inChamber = c->IsInChamber() ? 1 : 0;
1ee39b3a 1013 Double_t distance[2];
1014 GetDistanceToTracklet(distance, tracklet, c);
1015 Float_t theta = TMath::ATan(tracklet->GetZref(1));
1016 Float_t phi = TMath::ATan(tracklet->GetYref(1));
b795b480 1017 AliExternalTrackParam *trdPar = fkTrack->GetTrackIn();
7fe4e88b 1018 Float_t momentumMC = 0, momentumRec = trdPar ? trdPar->P() : fkTrack->P(); // prefer Track Low
1ee39b3a 1019 Int_t pdg = 0;
1020 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
efb9a759 1021 UShort_t tpcNCLS = fkESD ? fkESD->GetTPCncls() : 0;
1ee39b3a 1022 if(fkMC){
b795b480 1023 if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
1ee39b3a 1024 pdg = fkMC->GetPDG();
1025 }
1026 (*DebugStream()) << "PHt"
1027 << "Detector=" << detector
2e8eb67c 1028 << "Centrality="<< fCentralityClass
1ee39b3a 1029 << "crossing=" << crossing
b795b480 1030 << "inChamber=" << inChamber
1ee39b3a 1031 << "Timebin=" << localtime
1032 << "Charge=" << absoluteCharge
b795b480 1033 << "momentumMC=" << momentumMC
1034 << "momentumRec=" << momentumRec
1ee39b3a 1035 << "pdg=" << pdg
1036 << "theta=" << theta
1037 << "phi=" << phi
1038 << "kinkIndex=" << kinkIndex
efb9a759 1039 << "TPCncls=" << tpcNCLS
1ee39b3a 1040 << "dy=" << distance[0]
1041 << "dz=" << distance[1]
1042 << "c.=" << c
1043 << "\n";
1044 }
1045 }
1046 }
1047 return h;
1048}
1049
1050//_______________________________________________________
1051TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1052 //
1053 // Plots the average pulse height vs the distance from the anode wire
1054 // (plus const anode wire offset)
1055 //
4e7a049c 1056 if(track) fkTrack = track;
1057 if(!fkTrack){
1058 AliDebug(4, "No Track defined.");
1059 return NULL;
1ee39b3a 1060 }
2e8eb67c 1061 TH3 *h = NULL;
1062 if(!(h = dynamic_cast<TH3F *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
1ee39b3a 1063 AliWarning("No Histogram defined.");
db99a57a 1064 return NULL;
1ee39b3a 1065 }
4e7a049c 1066 AliTRDseedV1 *tracklet(NULL);
1067 AliTRDcluster *c(NULL);
749e9d16 1068 Double_t xd(0.), dqdl(0.);
1069 TVectorD vq(AliTRDseedV1::kNtb), vxd(AliTRDseedV1::kNtb), vdqdl(AliTRDseedV1::kNtb);
1ee39b3a 1070 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1071 if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
749e9d16 1072 Int_t det(tracklet->GetDetector());
1073 Bool_t rc(tracklet->IsRowCross());
4e7a049c 1074 for(Int_t ic(0); ic<AliTRDseedV1::kNtb; ic++){
1075 Bool_t kFIRST(kFALSE);
1076 if(!(c = tracklet->GetClusters(ic))){
1077 if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1078 } else kFIRST=kTRUE;
c732f879 1079 if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
4e7a049c 1080 xd = tracklet->GetX0() - c->GetX(); vxd[ic] = xd;
1081 dqdl=tracklet->GetdQdl(ic); vdqdl[ic] = dqdl;
1082 vq[ic]=c->GetQ();
1083 if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) vq[ic]+=c->GetQ();
2e8eb67c 1084 h->Fill(xd, dqdl, fCentralityClass);
749e9d16 1085 }
1086 if(DebugLevel() > 3){
749e9d16 1087 (*DebugStream()) << "PHx"
1088 << "det=" << det
1089 << "rc=" << rc
1090 << "xd=" << &vxd
1091 << "q=" << &vq
1092 << "dqdl=" << &vdqdl
1093 << "\n";
1ee39b3a 1094 }
1095 }
1096 return h;
1097}
1098
1099//_______________________________________________________
1100TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1101 //
1102 // Plot the cluster charge
1103 //
1104 if(track) fkTrack = track;
1105 if(!fkTrack){
4e7a049c 1106 AliDebug(4, "No Track defined.");
db99a57a 1107 return NULL;
1ee39b3a 1108 }
2e8eb67c 1109 TH2 *h = NULL;
1110 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeCluster)))){
1ee39b3a 1111 AliWarning("No Histogram defined.");
db99a57a 1112 return NULL;
1ee39b3a 1113 }
db99a57a 1114 AliTRDseedV1 *tracklet = NULL;
1115 AliTRDcluster *c = NULL;
1ee39b3a 1116 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1117 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
4e7a049c 1118 for(Int_t ic(0); ic < AliTRDseedV1::kNtb; ic++){
1119 Bool_t kFIRST(kFALSE);
1120 if(!(c = tracklet->GetClusters(ic))) {
1121 if(!(c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) continue;
1122 } else kFIRST = kTRUE;
1123 Float_t q(c->GetQ());
1124 if(kFIRST && (c = tracklet->GetClusters(AliTRDseedV1::kNtb+ic))) q+=c->GetQ();
2e8eb67c 1125 h->Fill(q, fCentralityClass);
1ee39b3a 1126 }
1127 }
1128 return h;
1129}
1130
1131//_______________________________________________________
1132TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1133 //
1134 // Plot the charge deposit per chamber
1135 //
1136 if(track) fkTrack = track;
1137 if(!fkTrack){
4e7a049c 1138 AliDebug(4, "No Track defined.");
db99a57a 1139 return NULL;
1ee39b3a 1140 }
2e8eb67c 1141 TH2 *h = NULL;
1142 if(!(h = dynamic_cast<TH2F *>(fContainer->At(kChargeTracklet)))){
1ee39b3a 1143 AliWarning("No Histogram defined.");
db99a57a 1144 return NULL;
1ee39b3a 1145 }
db99a57a 1146 AliTRDseedV1 *tracklet = NULL;
1147 AliTRDcluster *c = NULL;
1ee39b3a 1148 Double_t qTot = 0;
1149 Int_t nTracklets =fkTrack->GetNumberOfTracklets();
8428f55d 1150 for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1ee39b3a 1151 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1152 qTot = 0.;
1153 for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1154 if(!(c = tracklet->GetClusters(ic))) continue;
1155 qTot += TMath::Abs(c->GetQ());
1156 }
2e8eb67c 1157 h->Fill(qTot, fCentralityClass);
1ee39b3a 1158 if(DebugLevel() > 3){
1159 Int_t crossing = (Int_t)tracklet->IsRowCross();
1160 Int_t detector = tracklet->GetDetector();
1161 Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1162 Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1163 Float_t momentum = 0.;
1164 Int_t pdg = 0;
1165 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1166 UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1167 if(fkMC){
1168 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1169 pdg = fkMC->GetPDG();
1170 }
1171 (*DebugStream()) << "ChargeTracklet"
1172 << "Detector=" << detector
2e8eb67c 1173 << "Centrality="<< fCentralityClass
1ee39b3a 1174 << "crossing=" << crossing
1175 << "momentum=" << momentum
1176 << "nTracklets="<< nTracklets
1177 << "pdg=" << pdg
1178 << "theta=" << theta
1179 << "phi=" << phi
1180 << "kinkIndex=" << kinkIndex
1181 << "TPCncls=" << nclsTPC
1182 << "QT=" << qTot
1183 << "\n";
1184 }
1185 }
1186 return h;
1187}
1188
1189//_______________________________________________________
1190TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1191 //
1192 // Plot the number of tracks per Sector
1193 //
1194 if(track) fkTrack = track;
1195 if(!fkTrack){
4e7a049c 1196 AliDebug(4, "No Track defined.");
db99a57a 1197 return NULL;
1ee39b3a 1198 }
db99a57a 1199 TH1 *h = NULL;
1ee39b3a 1200 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1201 AliWarning("No Histogram defined.");
db99a57a 1202 return NULL;
1ee39b3a 1203 }
1204
1205 // TODO we should compare with
1206 // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1207
db99a57a 1208 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 1209 Int_t sector = -1;
1210 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1211 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1212 sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1213 break;
1214 }
1215 h->Fill(sector);
1216 return h;
1217}
1218
1219
1ee39b3a 1220//________________________________________________________
1221void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1222{
1223 Float_t x = c->GetX();
1224 dist[0] = c->GetY() - tracklet->GetYat(x);
1225 dist[1] = c->GetZ() - tracklet->GetZat(x);
1226}
1227
44fd53ca 1228//________________________________________________________
efb9a759 1229void AliTRDcheckDET::GetEtaPhiAt(const AliExternalTrackParam *track, Double_t x, Double_t &eta, Double_t &phi){
44fd53ca 1230 //
1231 // Get phi and eta at a given radial position
1232 //
2cef552b 1233 if(!fEvent) return;
44fd53ca 1234 AliExternalTrackParam workpar(*track);
1235
1236 Double_t posLocal[3];
2e8eb67c 1237 Bool_t sucPos = workpar.GetXYZAt(x, fEvent->GetRunInfo()->GetMagneticField(), posLocal);
44fd53ca 1238 Double_t sagPhi = sucPos ? TMath::ATan2(posLocal[1], posLocal[0]) : 0.;
1239 phi = sagPhi;
1240 eta = workpar.Eta();
1241}
1242
1ee39b3a 1243
1244//_______________________________________________________
d48822f2 1245TH1* AliTRDcheckDET::MakePlotChi2() const
1ee39b3a 1246{
1247// Plot chi2/track normalized to number of degree of freedom
1248// (tracklets) and compare with the theoretical distribution.
1249//
1250// Alex Bercuci <A.Bercuci@gsi.de>
1251
c0a6ac06 1252 return NULL;
1253
7fe4e88b 1254/* TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1ee39b3a 1255 TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
c0a6ac06 1256 f.SetParLimits(1,1, 1e100);
1ee39b3a 1257 TLegend *leg = new TLegend(.7,.7,.95,.95);
1258 leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
db99a57a 1259 TH1D *h1 = NULL;
1ee39b3a 1260 Bool_t kFIRST = kTRUE;
1261 for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1262 h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1263 if(h1->Integral()<50) continue;
1264 h1->Scale(1./h1->Integral());
1265 h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1266 h1->SetLineColor(il);h1->SetLineStyle(2);
1267 f.SetParameter(1, .5*il);f.SetLineColor(il);
1268 h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1269 leg->AddEntry(h1, Form("%d", il), "l");
1270 if(kFIRST){
1271 h1->GetXaxis()->SetRangeUser(0., 25.);
1272 }
1273 kFIRST = kFALSE;
1274 }
1275 leg->Draw();
1276 gPad->SetLogy();
7fe4e88b 1277 return h1;*/
1ee39b3a 1278}
1279
1280
1281//________________________________________________________
1282TH1* AliTRDcheckDET::MakePlotNTracklets(){
8a3dd7a0 1283 //
1284 // Make nice bar plot of the number of tracklets in each method
1285 //
1286 TH3 *hTracklets3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
1287 if(!hTracklets3D){
1288 AliError("Tracklet Histogram not found");
1289 return NULL;
1290 }
1291 TH1 *hBAR = hTracklets3D->Project3D("y");
1292 hBAR->SetName("hBAR");
1293 hBAR->SetTitle("Number of Tracklets");
cc98ff07 1294 if(hBAR->Integral()) hBAR->Scale(1./hBAR->Integral());
8a3dd7a0 1295 hBAR->SetLineColor(kBlack);
1296 hBAR->SetLineWidth(2);
1297 hBAR->Draw();
1298
1299 // Draw also centrality-dependent plots
1300 TH1 *hBarCent[AliTRDeventInfo::kCentralityClasses];
1301 Int_t nHistsCentrality = 0;
1302 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1303 hTracklets3D->GetZaxis()->SetRange(icent+2, icent+2);
1304 hBarCent[icent] = hTracklets3D->Project3D("y");
1305 if(!(hBarCent[icent] && hBarCent[icent]->GetEntries())){
1306 delete hBarCent[icent];
1307 hBarCent[icent] = NULL;
1308 continue;
1309 }
1310 hBarCent[icent]->SetName(Form("hBarCent_%d", icent));
1311 hBarCent[icent]->SetTitle("Number of Tracklets");
1312 hBarCent[icent]->Scale(1./hBarCent[icent]->Integral());
1313 hBarCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1314 hBarCent[icent]->Draw("same");
1315 nHistsCentrality++;
1316 }
1317 hTracklets3D->GetZaxis()->SetRange(0, hTracklets3D->GetNbinsZ());
1318 AliInfo(Form("Number of Centrality Classes: %d", nHistsCentrality));
1319 if(nHistsCentrality){
1320 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1321 leg->SetBorderSize(0);
1322 leg->SetFillStyle(0);
1323 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++)
1324 if(hBarCent[icent]) leg->AddEntry(hBarCent[icent], Form("Centrality Class %d", icent), "l");
1325 leg->Draw();
1326 gPad->Update();
1327 }
1328 return hBAR;
1329
1330 /*
1331 TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1332 TH1D *hBAR = tmp->ProjectionY();
1333 TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1334 TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1335 TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1336 leg->SetBorderSize(1);
1337 leg->SetFillColor(0);
1338
1339 Float_t scale = hCON->Integral();
1340 if(scale) hCON->Scale(100./scale);
1341 hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1342 hCON->SetBarWidth(0.2);
1343 hCON->SetBarOffset(0.6);
1344 hCON->SetStats(kFALSE);
1345 hCON->GetYaxis()->SetRangeUser(0.,40.);
1346 hCON->GetYaxis()->SetTitleOffset(1.2);
1347 hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1348 hCON->SetMaximum(55.);
1349
1350 if(scale) hBAR->Scale(100./scale);
1351 hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1352 hBAR->SetBarWidth(0.2);
1353 hBAR->SetBarOffset(0.2);
1354 hBAR->SetTitle("");
1355 hBAR->SetStats(kFALSE);
1356 hBAR->GetYaxis()->SetRangeUser(0.,40.);
1357 hBAR->GetYaxis()->SetTitleOffset(1.2);
1358 hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1359
1360 if(scale) hSTA->Scale(100./scale);
1361 hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1362 hSTA->SetBarWidth(0.2);
1363 hSTA->SetBarOffset(0.4);
1364 hSTA->SetTitle("");
1365 hSTA->SetStats(kFALSE);
1366 hSTA->GetYaxis()->SetRangeUser(0.,40.);
1367 hSTA->GetYaxis()->SetTitleOffset(1.2);
1368 hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1369 leg->Draw();
1370 gPad->Update();
1371 return hCON;
1372 */
1373}
1374
1375//________________________________________________________
1376void AliTRDcheckDET::MakePlotNclustersTrack(){
1377 //
1378 // Plot number of clusters
1379 // Put histos from all centrality classes into one pad
1380 //
1381 TH2 *hClusters = dynamic_cast<TH2 *>(fContainer->FindObject("hNcls"));
1382 if(!hClusters){
1383 AliError("Cluster histogram not found in the output");
1384 }
1385 TH1 *hAllCentrality = hClusters->ProjectionX("hNcls_all");
1386 hAllCentrality->SetTitle("Number of clusters/track");
cc98ff07 1387 if(hAllCentrality->Integral()>0.) hAllCentrality->Scale(1./hAllCentrality->Integral());
8a3dd7a0 1388 hAllCentrality->SetLineColor(kBlack);
1389 hAllCentrality->SetLineWidth(2);
1390 hAllCentrality->GetYaxis()->SetRangeUser(0, 0.02);
1391 hAllCentrality->Draw();
1392 PutTrendValue("NClustersTrack", hAllCentrality->GetMean());
1393 PutTrendValue("NClustersTrackRMS", hAllCentrality->GetRMS());
1394
1395 // Now look at single centrality classes
1396 TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
1397 Int_t nHistsCentrality = 0;
1398 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1399 hProjCentral[icent] = hClusters->ProjectionX(Form("hNcls_%d", icent), icent+1, icent+1);
1400 if(!hProjCentral[icent]->GetEntries()){
1401 delete hProjCentral[icent];
1402 hProjCentral[icent] = NULL;
1403 continue;
1404 }
1405 hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
1406 hProjCentral[icent]->SetTitle("Number of clusters/track");
1407 hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
1408 hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.03);
1409 hProjCentral[icent]->Draw("same");
1410 nHistsCentrality++;
1411 }
1412 AliInfo(Form("%d centrality classes found", nHistsCentrality));
1413 if(nHistsCentrality){
1414 // Draw nice legend
1415 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1416 leg->SetBorderSize(0);
1417 leg->SetFillStyle(0);
1418 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1419 if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
1420 }
1421 leg->Draw();
1422 gPad->Update();
1423 }
1424}
1425
1426//________________________________________________________
1427void AliTRDcheckDET::MakePlotNclustersTracklet(){
1428 //
1429 // Plot number of clusters for different centrality classes
1430 //
1431 TH2 *hClusters = dynamic_cast<TH2*>(fContainer->FindObject("hNclTls"));
1432 if(!hClusters){
1433 AliError("Histogram for the number of clusters per tracklet not available");
1434 return;
1435 }
1436 TH1 *hAllCentrality = hClusters->ProjectionX("hNclsTls_all");
1437 hAllCentrality->SetTitle("Number of clusters/track");
cc98ff07 1438 if(hAllCentrality->Integral()>0.) hAllCentrality->Scale(1./hAllCentrality->Integral());
8a3dd7a0 1439 hAllCentrality->SetLineColor(kBlack);
1440 hAllCentrality->SetLineWidth(2);
1441 hAllCentrality->GetYaxis()->SetRangeUser(0, 0.3);
1442 hAllCentrality->Draw("pc");
1443 PutTrendValue("NClustersTracklet", hAllCentrality->GetMean());
1444 PutTrendValue("NClustersTrackletRMS", hAllCentrality->GetRMS());
1445
1446 // Now look at single centrality classes
1447 TH1 *hProjCentral[AliTRDeventInfo::kCentralityClasses];
1448 Int_t nHistsCentrality = 0;
1449 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1450 hProjCentral[icent] = hClusters->ProjectionX(Form("hNclsTls_%d", icent), icent+1, icent+1);
1451 if(!hProjCentral[icent]->GetEntries()){
1452 delete hProjCentral[icent];
1453 hProjCentral[icent] = NULL;
1454 continue;
1455 }
cc98ff07 1456 if(hProjCentral[icent]->Integral()>0.) hProjCentral[icent]->Scale(1./hProjCentral[icent]->Integral());
8a3dd7a0 1457 hProjCentral[icent]->SetTitle("Number of clusters/track");
1458 hProjCentral[icent]->SetLineColor(fkColorsCentrality[icent]);
1459 hProjCentral[icent]->GetYaxis()->SetRangeUser(0, 0.3);
1460 hProjCentral[icent]->Draw("pcsame");
1461 nHistsCentrality++;
1462 }
1463 AliInfo(Form("%d centrality classes found", nHistsCentrality));
1464 if(nHistsCentrality){
1465 // Draw nice legend
1466 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1467 leg->SetBorderSize(0);
1468 leg->SetFillStyle(0);
1469 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1470 if(hProjCentral[icent]) leg->AddEntry(hProjCentral[icent], Form("Centrality Class %d", icent), "l");
1471 }
1472 leg->Draw();
1473 gPad->Update();
1474 }
1475}
1476
1477//________________________________________________________
1478void AliTRDcheckDET::MakePlotTrackletCharge(){
1479 //
1480 // Draw tracklet charge for different centrality classes
1481 //
1482 TH2 *hQt = dynamic_cast<TH2*>(fContainer->FindObject("hQtrklt"));
1483 if(!hQt){
1484 AliError("Histogram for tracklet charges not found");
1485 return;
1486 }
1487 // First Project all charhes
1488 TH1 *hQtAll = hQt->ProjectionX("hQtAll");
1489 hQtAll->SetTitle("Tracklet Charge");
cc98ff07 1490 Double_t scalefactor = hQtAll->GetMaximum()>0.?(0.7 * hQtAll->Integral() / hQtAll->GetMaximum()):1;
1491 if(hQtAll->Integral()>0.) hQtAll->Scale(scalefactor/hQtAll->Integral());
8a3dd7a0 1492 hQtAll->GetYaxis()->SetRangeUser(0., 1.);
1493 hQtAll->SetLineColor(kBlack);
1494 hQtAll->SetLineWidth(2);
1495 hQtAll->Draw("c");
1496 PutTrendValue("ChargeTracklet", hQtAll->GetMaximumBin());
1497 PutTrendValue("ChargeTrackletRMS", hQtAll->GetRMS());
1498
1499 TH1 *hQtCent[AliTRDeventInfo::kCentralityClasses];
1500 Int_t nHistsCentrality = 0;
1501 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1502 hQtCent[icent] = hQt->ProjectionX(Form("hQt_%d", icent), icent+1, icent+1);
1503 if(!hQtCent[icent]->GetEntries()){
d2e80884 1504 delete hQtCent[icent]; hQtCent[icent]=NULL;
8a3dd7a0 1505 continue;
1506 }
1507 hQtCent[icent]->SetTitle("Tracklet Charge");
cc98ff07 1508 if(hQtCent[icent]->Integral()>0.) hQtCent[icent]->Scale(scalefactor/hQtCent[icent]->Integral());
8a3dd7a0 1509 hQtCent[icent]->GetYaxis()->SetRangeUser(0., 1.);
1510 hQtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1511 hQtCent[icent]->Draw("csame");
1512 nHistsCentrality++;
1513 }
1514 AliInfo(Form("%d centrality classes found", nHistsCentrality));
1515 if(nHistsCentrality){
1516 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1517 leg->SetBorderSize(0);
1518 leg->SetFillStyle(0);
1519 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1520 if(hQtCent[icent]) leg->AddEntry(hQtCent[icent], Form("Centrality Class %d", icent), "l");
1521 }
1522 leg->Draw();
1523 gPad->Update();
1524 }
1ee39b3a 1525}
1526
0c76cfa4 1527//________________________________________________________
1528void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1529 //
1530 // Plot abundance of tracks with number of tracklets as function of momentum
1531 //
147c3968 1532
147c3968 1533 Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1534 TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1535 for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1536 g[itl] = new TGraphErrors();
1537 g[itl]->SetLineColor(color[itl]);
1538 g[itl]->SetMarkerColor(color[itl]);
1539 g[itl]->SetMarkerStyle(20 + itl);
1540 }
1541
8a3dd7a0 1542 TH3 *hBar3D = dynamic_cast<TH3F *>(fContainer->FindObject("hNtlsBAR"));
1543 if(!hBar3D){
1544 AliError("Histogram for the number of tracklets vs p not available");
1545 return;
1546 }
1547 TH2 *hBar = (TH2 *)hBar3D->Project3D("yx");
147c3968 1548 TAxis *ax(hBar->GetXaxis());
1549 Int_t np(ax->GetNbins());
1550 for(Int_t ipBin(1); ipBin<np; ipBin++){
1551 h = hBar->ProjectionY("npBin", ipBin, ipBin);
81979445 1552 if(!Int_t(h->Integral())) continue;
147c3968 1553 h->Scale(100./h->Integral());
1554 Float_t p(ax->GetBinCenter(ipBin));
1555 Float_t dp(ax->GetBinWidth(ipBin));
81979445 1556 Int_t ip(g[0]->GetN());
147c3968 1557 for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
81979445 1558 g[itl]->SetPoint(ip, p, h->GetBinContent(itl+1));
1559 g[itl]->SetPointError(ip, dp/3.46, h->GetBinError(itl+1));
147c3968 1560 }
1561 }
1562
81979445 1563 TLegend *leg = new TLegend(0.76, 0.6, 1., 0.9);
94f7dff7 1564 leg->SetBorderSize(0);
1565 leg->SetHeader("Tracklet/Track");
1566 leg->SetFillStyle(0);
147c3968 1567 h = hBar->ProjectionX("npxBin"); h->Reset();
81979445 1568 h->SetTitle("");
1569 h->GetYaxis()->SetRangeUser(1., 99.);
1570 h->GetYaxis()->SetMoreLogLabels();
1571 h->GetYaxis()->CenterTitle();
1572 h->GetYaxis()->SetTitleOffset(1.2);
147c3968 1573 h->SetYTitle("Prob. [%]");
81979445 1574 h->GetXaxis()->SetRangeUser(0.4, 12.);
1575 h->GetXaxis()->SetMoreLogLabels();
1576 h->GetXaxis()->CenterTitle();
147c3968 1577 h->Draw("p");
1578 for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1579 g[itl]->Draw("pc");
1580 leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
0c76cfa4 1581 }
147c3968 1582
0c76cfa4 1583 leg->Draw();
81979445 1584 gPad->SetLogx();gPad->SetLogy();
0c76cfa4 1585}
1586
1ee39b3a 1587//________________________________________________________
c0a6ac06 1588Bool_t AliTRDcheckDET::MakePlotPulseHeight(){
8a3dd7a0 1589 //
1590 // Create Plot of the Pluse Height Spectrum
1591 //
1592 TCanvas *output = gPad->GetCanvas();
1593 output->Divide(2);
1594 output->cd(1);
1595
1596 TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1597 //TH3F *hPhx = dynamic_cast<TH3F *>(arr->At(0)),
1598 TH3F *hPht = dynamic_cast<TH3F *>(arr->At(1));
1599 if(!hPht) return kFALSE;
1600 // Project centrality of the 2 histograms
1601 //TH2 *hProjCentX = dynamic_cast<TH2 *>(hPhx->Project3D("yx")),
1602 TH2 *hProjCentT = dynamic_cast<TH2 *>(hPht->Project3D("yx"));
6475ec36 1603 if(!hProjCentT) return kFALSE;
8a3dd7a0 1604 //hProjCentX->SetName("hProjCentX");
1605 hProjCentT->SetName("hProjCentT");
1606 // Draw 2D histogram versus time on pad 2
1607 output->cd(2);
1608 hProjCentT->SetTitle("2D Pulse Height vs. Time");
1609 hProjCentT->GetYaxis()->SetTitleOffset(1.8);
6abdb9ed 1610 hProjCentT->GetYaxis()->SetRangeUser(0., 600.);
8a3dd7a0 1611 hProjCentT->SetStats(kFALSE);
1612 hProjCentT->Draw("colz");
1613
6abdb9ed 1614 // Fit each slice of the pulse height spectrum with a Landau function
1615 TGraphErrors *landaufit = new TGraphErrors;
1616 landaufit->SetTitle();
1617 landaufit->SetMarkerColor(kWhite);
1618 landaufit->SetLineColor(kWhite);
1619 landaufit->SetMarkerStyle(25);
1620 TH1 *projection;
1621 Int_t ntime = 0;
1622 for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetLast(); it++){
1623 projection = hProjCentT->ProjectionY("proj", it, it);
1624 if(!projection->GetEntries()){
1625 delete projection;
1626 continue;
1627 }
1628 TF1 fitfun("fitfun", "landau", 0, 1000);
1629 fitfun.SetParLimits(0, 1e-1, 1e9); fitfun.SetParameter(0, 1000);
1630 fitfun.SetParLimits(1, 10, 200); fitfun.SetParameter(1, 50);
1631 fitfun.SetParLimits(2, 1e-1, 200); fitfun.SetParameter(2, 10);
1632 projection->Fit(&fitfun, "QN", "", 0, 1000);
1633 landaufit->SetPoint(ntime, hProjCentT->GetXaxis()->GetBinCenter(it), fitfun.GetParameter(1));
1634 landaufit->SetPointError(ntime, hProjCentT->GetXaxis()->GetBinWidth(it)/2., fitfun.GetParError(1));
1635 ntime++;
1636 delete projection;
1637 }
1638 landaufit->Draw("lpesame");
1639
8a3dd7a0 1640 // Fill 1D PHS as function of time resp. radius (same binning of the 2 histograms)
1641 //TH1 *hPhsX = new TH1F("hPhsX", "Average PH vs X", hProjCentT->GetXaxis()->GetNbins(), hProjCentT->GetXaxis()->GetXbins()->GetArray()),
6abdb9ed 1642 Int_t nbinsT = hProjCentT->GetXaxis()->GetNbins();
1643 TH1 *hPhsT = new TH1F("hPhsT", "Average Ph vs Time; Time (100 ns); Average Pulse Height (a.u.)", nbinsT, -0.5, nbinsT - 0.5),
8a3dd7a0 1644 *htmp = NULL;
1645 /*for(Int_t irad = 1; irad <= hProjCentX->GetXaxis()->GetNbins(); irad++){
1646 htmp = hProjCentX->ProjectionY("tmp", irad, irad);
1647 hPhsX->SetBinContent(irad, htmp->GetMean());
1648 hPhsX->SetBinError(irad, htmp->GetMeanError());
1649 delete htmp;
1650 }*/
1651 for(Int_t it = 1; it <= hProjCentT->GetXaxis()->GetNbins(); it++){
1652 htmp = hProjCentT->ProjectionY("tmp", it, it);
1653 hPhsT->SetBinContent(it, htmp->GetMean());
1654 hPhsT->SetBinError(it, htmp->GetMeanError());
1655 delete htmp;
1656 }
1657 output->cd(1);
1658 // Draw 1D histograms
1659 if(hPhsT->GetEntries()){
1660 hPhsT->SetMarkerStyle(24);
1661 hPhsT->SetMarkerColor(kBlack);
1662 hPhsT->SetLineColor(kBlack);
1663 hPhsT->GetYaxis()->SetTitleOffset(1.5);
1664 hPhsT->Draw("e1");
1665 // Now fit the PHS with respect to time to get plateau and slope
1666 // Trending for the pulse height: plateau value, slope and timebin of the maximum
1667 TLinearFitter fit(1,"pol1");
1668 Double_t time = 0.;
1669 for(Int_t itime = 10; itime <= 20; itime++){
1670 time = Double_t(itime);
1671 Double_t err(hPhsT->GetBinError(itime + 1));
1672 if(err>1.e-10) fit.AddPoint(&time, hPhsT->GetBinContent(itime + 1), err);
1673 }
1674 if(!fit.Eval()){
1675 Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1676 Double_t slope = fit.GetParameter(1);
1677 PutTrendValue("PHplateau", plateau);
1678 PutTrendValue("PHslope", slope);
1679 PutTrendValue("PHamplificationPeak", static_cast<Double_t>(hPhsT->GetMaximumBin()-1));
1680 AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(hPhsT->GetMaximumBin()-1)));
1681 }
1682 }
1683 /*if(hPhsX->GetEntries()){
1684 hPhsX->SetMarkerStyle(22);
1685 hPhsX->SetMarkerColor(kBlue);
1686 hPhsX->SetLineColor(kBlue);
1687 hPhsX->GetYaxis()->SetTitleOffset(1.5);
1688 hPhsX->Draw("e1same");
1689 // create axis according to the histogram dimensions of the original second histogram
1690 TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1691 gPad->GetUymax(),
1692 gPad->GetUxmax(),
1693 gPad->GetUymax(),
1694 -0.08, 4.88, 510,"-L");
1695 axis->SetLineColor(kBlue);
1696 axis->SetLabelColor(kBlue);
1697 axis->SetTextColor(kBlue);
1698 axis->SetTitle("x_{0}-x_{c} [cm]");
1699 axis->Draw();
1700 gPad->Update();
1701 }*/
1702
1703 // Centrality-dependent Pulse-Height Spectrum
1704 TH1 *hPhtCent[AliTRDeventInfo::kCentralityClasses];
1705 TH2 *hPtmp;
1706 Int_t nHistsCentrality = 0;
40600e64 1707 TArrayD ptbinning(hPht->GetNbinsX()+1);
1708 for(Int_t ibin = 1; ibin <=hPht->GetNbinsX(); ibin++) ptbinning[ibin-1] = hPht->GetXaxis()->GetBinLowEdge(ibin);
1709 ptbinning[hPht->GetNbinsX()] = hPht->GetXaxis()->GetBinUpEdge(hPht->GetNbinsX());
8a3dd7a0 1710 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1711 hPht->GetZaxis()->SetRange(icent+2,icent+2);
1712 hPtmp = dynamic_cast<TH2*>(hPht->Project3D("yx"));
1713 if(!(hPtmp && hPtmp->GetEntries())){
1714 hPhtCent[icent] = NULL;
1715 continue;
1716 }
40600e64 1717 hPhtCent[icent] = new TH1F(Form("hPhtCent_%d", icent), "Average Ph vs Time", hPtmp->GetNbinsX(), ptbinning.GetArray());
8a3dd7a0 1718 for(Int_t it = 1; it <= hPtmp->GetNbinsX(); it++){
1719 htmp = hPtmp->ProjectionY("htmp", it, it);
1720 hPhtCent[icent]->SetBinContent(it, htmp->GetMean());
1721 hPhtCent[icent]->SetBinError(it, htmp->GetMeanError());
1722 delete htmp;
1723 }
1724 delete hPtmp;
1725 hPhtCent[icent]->SetMarkerStyle(24);
1726 hPhtCent[icent]->SetMarkerColor(fkColorsCentrality[icent]);
1727 hPhtCent[icent]->SetLineColor(fkColorsCentrality[icent]);
1728 hPhtCent[icent]->Draw("e1same");
1729 nHistsCentrality++;
1730 }
1731 hPht->GetZaxis()->SetRange(0, hPht->GetNbinsZ());
1732 if(nHistsCentrality){
1733 TLegend *leg = new TLegend(0.5, 0.6, 0.89, 0.89);
1734 leg->SetBorderSize(0);
1735 leg->SetFillStyle(0);
1736 for(Int_t icent = 0; icent < AliTRDeventInfo::kCentralityClasses; icent++){
1737 if(hPhtCent[icent]) leg->AddEntry(hPhtCent[icent], Form("Centrality Class %d", icent), "p");
1738 }
1739 leg->Draw();
1740 gPad->Update();
1741 }
1742
1743 // delete 2D Projection of the PHS vs x since it is only used to calculate the 1D projection
1744 //delete hProjCentX;
1745
1746 return kTRUE;
1ee39b3a 1747}
1748
44fd53ca 1749//________________________________________________________
1750void AliTRDcheckDET::MakePlotMeanClustersLayer(){
1751 //
1752 // Create Summary plot for the mean number of clusters per layer
1753 //
1754 TCanvas *output = gPad->GetCanvas();
1755 output->Divide(3,2);
1756 TObjArray *histos = (TObjArray *)fContainer->At(kNclustersLayer);
1757 if(!histos){
1758 AliWarning("Histos for each layer not found");
1759 return;
1760 }
1761 TProfile2D *hlayer = NULL;
1762 for(Int_t ily = 0; ily < AliTRDgeometry::kNlayer; ily++){
7fe4e88b 1763 if(!(hlayer = dynamic_cast<TProfile2D *>(histos->At(ily)))) continue;
44fd53ca 1764 output->cd(ily + 1);
1765 gPad->SetGrid(0,0);
1766 hlayer->Draw("colz");
1767 }
1768}
1769
1ee39b3a 1770//________________________________________________________
1771Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1772 //
1773 // Draw nice bar plots
1774 //
1775 if(!histo->GetEntries()) return kFALSE;
cc98ff07 1776 if(histo->Integral()>0.) histo->Scale(100./histo->Integral());
1ee39b3a 1777 histo->SetFillColor(color);
1778 histo->SetBarOffset(.2);
1779 histo->SetBarWidth(.6);
1780 histo->Draw("bar1");
1781 return kTRUE;
960a59e0 1782}