]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDcheckDET.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add variable TRD segmentation levels. Possible options are
[u/mrichter/AliRoot.git] / PWG1 / 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
26#include <TAxis.h>
27#include <TCanvas.h>
28#include <TFile.h>
29#include <TH1F.h>
30#include <TH1I.h>
a5d9fe6f 31#include <TH2F.h>
1ee39b3a 32#include <TF1.h>
33#include <TGaxis.h>
34#include <TGraph.h>
147c3968 35#include <TGraphErrors.h>
1ee39b3a 36#include <TLegend.h>
3907f080 37#include <TLinearFitter.h>
1ee39b3a 38#include <TMath.h>
39#include <TMap.h>
40#include <TObjArray.h>
41#include <TObject.h>
42#include <TObjString.h>
43
44#include <TPad.h>
45#include <TProfile.h>
46#include <TProfile2D.h>
47#include <TROOT.h>
f8f46e4d 48#include <TChain.h>
1ee39b3a 49
50#include "AliLog.h"
51#include "AliTRDcluster.h"
52#include "AliESDHeader.h"
53#include "AliESDRun.h"
54#include "AliESDtrack.h"
a5d9fe6f 55#include "AliExternalTrackParam.h"
1ee39b3a 56#include "AliTRDgeometry.h"
57#include "AliTRDpadPlane.h"
58#include "AliTRDSimParam.h"
59#include "AliTRDseedV1.h"
60#include "AliTRDtrackV1.h"
61#include "AliTRDtrackerV1.h"
62#include "AliTRDReconstructor.h"
63#include "AliTrackReference.h"
64#include "AliTrackPointArray.h"
65#include "AliTracker.h"
66#include "TTreeStream.h"
67
68#include "info/AliTRDtrackInfo.h"
69#include "info/AliTRDeventInfo.h"
70#include "AliTRDcheckDET.h"
71
72#include <cstdio>
73#include <iostream>
74
75ClassImp(AliTRDcheckDET)
76
77//_______________________________________________________
78AliTRDcheckDET::AliTRDcheckDET():
f2e89a4c 79 AliTRDrecoTask()
db99a57a 80 ,fEventInfo(NULL)
81 ,fTriggerNames(NULL)
82 ,fReconstructor(NULL)
83 ,fGeo(NULL)
c732f879 84 ,fFlags(0)
1ee39b3a 85{
86 //
87 // Default constructor
88 //
f2e89a4c 89 SetNameTitle("checkDET", "Basic TRD data checker");
f8f46e4d 90}
91
705f8b0a 92//_______________________________________________________
f8f46e4d 93AliTRDcheckDET::AliTRDcheckDET(char* name):
94 AliTRDrecoTask(name, "Basic TRD data checker")
db99a57a 95 ,fEventInfo(NULL)
96 ,fTriggerNames(NULL)
97 ,fReconstructor(NULL)
98 ,fGeo(NULL)
f8f46e4d 99 ,fFlags(0)
100{
101 //
102 // Default constructor
103 //
104 DefineInput(2, AliTRDeventInfo::Class());
105
1ee39b3a 106 fReconstructor = new AliTRDReconstructor;
107 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
108 fGeo = new AliTRDgeometry;
109 InitFunctorList();
110}
111
f8f46e4d 112
1ee39b3a 113//_______________________________________________________
114AliTRDcheckDET::~AliTRDcheckDET(){
115 //
116 // Destructor
117 //
118 if(fTriggerNames) delete fTriggerNames;
119 delete fReconstructor;
120 delete fGeo;
121}
122
123//_______________________________________________________
124void AliTRDcheckDET::ConnectInputData(Option_t *opt){
125 //
126 // Connect the Input data with the task
127 //
128 AliTRDrecoTask::ConnectInputData(opt);
f8f46e4d 129 fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));
1ee39b3a 130}
131
132//_______________________________________________________
f8f46e4d 133void AliTRDcheckDET::UserCreateOutputObjects(){
1ee39b3a 134 //
135 // Create Output Objects
136 //
f8f46e4d 137 OpenFile(1,"RECREATE");
1ee39b3a 138 fContainer = Histos();
139 if(!fTriggerNames) fTriggerNames = new TMap();
140}
141
142//_______________________________________________________
f8f46e4d 143void AliTRDcheckDET::UserExec(Option_t *opt){
1ee39b3a 144 //
145 // Execution function
146 // Filling TRD quality histos
147 //
148 if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
b4414720 149 AliTRDrecoTask::UserExec(opt);
1ee39b3a 150 Int_t nTracks = 0; // Count the number of tracks per event
151 Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
152 TString triggername = fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
153 AliDebug(6, Form("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data()));
154 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
155 for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
156 if(!fTracks->UncheckedAt(iti)) continue;
157 AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
158 if(!fTrackInfo->GetTrack()) continue;
159 nTracks++;
160 }
f8f46e4d 161
1ee39b3a 162 if(nTracks){
163 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
164 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
165 }
166 if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
167 fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
168 // also set the label for both histograms
169 TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
170 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
171 histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
172 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
173 }
f8f46e4d 174 PostData(1, fContainer);
1ee39b3a 175}
176
177
178//_______________________________________________________
179Bool_t AliTRDcheckDET::PostProcess(){
180 //
181 // Do Postprocessing (for the moment set the number of Reference histograms)
182 //
183
db99a57a 184 TH1 * h = NULL;
1ee39b3a 185
186 // Calculate of the trigger clusters purity
187 h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
188 TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
189 h1->Divide(h);
190 Float_t purities[20], val = 0;
191 TString triggernames[20];
192 Int_t nTriggerClasses = 0;
193 for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
194 if((val = h1->GetBinContent(ibin))){
195 purities[nTriggerClasses] = val;
196 triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
197 nTriggerClasses++;
198 }
199 }
200 h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
201 TAxis *ax = h->GetXaxis();
202 for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
203 h->Fill(itrg, purities[itrg]);
204 ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
205 }
206 ax->SetRangeUser(-0.5, nTriggerClasses+.5);
207 h->GetYaxis()->SetRangeUser(0,1);
208
209 // track status
210 h=dynamic_cast<TH1F*>(fContainer->At(kTrackStatus));
211 Float_t ok = h->GetBinContent(1);
212 Int_t nerr = h->GetNbinsX();
213 for(Int_t ierr=nerr; ierr--;){
214 h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/ok);
215 }
216 h->SetBinContent(1, 0.);
217
218 // tracklet status
f8f46e4d 219
1ee39b3a 220 TObjArray *arr = dynamic_cast<TObjArray*>(fContainer->UncheckedAt(kTrackletStatus));
221 for(Int_t ily = AliTRDgeometry::kNlayer; ily--;){
222 h=dynamic_cast<TH1F*>(arr->At(ily));
223 Float_t okB = h->GetBinContent(1);
224 Int_t nerrB = h->GetNbinsX();
225 for(Int_t ierr=nerrB; ierr--;){
226 h->SetBinContent(ierr+1, 1.e2*h->GetBinContent(ierr+1)/okB);
227 }
228 h->SetBinContent(1, 0.);
229 }
230
0c76cfa4 231 fNRefFigures = 17;
1ee39b3a 232
233 return kTRUE;
234}
235
236//_______________________________________________________
237Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
238 //
239 // Setting Reference Figures
240 //
0c76cfa4 241 enum FigureType_t{
242 kFigNclustersTrack,
243 kFigNclustersTracklet,
244 kFigNtrackletsTrack,
245 kFigNTrackletsP,
246 kFigNtrackletsCross,
247 kFigNtrackletsFindable,
248 kFigNtracksEvent,
249 kFigNtracksSector,
250 kFigTrackStatus,
251 kFigTrackletStatus,
252 kFigChi2,
253 kFigPH,
254 kFigChargeCluster,
255 kFigChargeTracklet,
256 kFigNeventsTrigger,
257 kFigNeventsTriggerTracks,
258 kFigTriggerPurity
259 };
1ee39b3a 260 gPad->SetLogy(0);
261 gPad->SetLogx(0);
db99a57a 262 TH1 *h = NULL; TObjArray *arr=NULL;
263 TLegend *leg = NULL;
1ee39b3a 264 Bool_t kFIRST(1);
265 switch(ifig){
0c76cfa4 266 case kFigNclustersTrack:
1ee39b3a 267 (h=(TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
268 PutTrendValue("NClustersTrack", h->GetMean());
269 PutTrendValue("NClustersTrackRMS", h->GetRMS());
270 return kTRUE;
0c76cfa4 271 case kFigNclustersTracklet:
1ee39b3a 272 (h =(TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
273 PutTrendValue("NClustersTracklet", h->GetMean());
274 PutTrendValue("NClustersTrackletRMS", h->GetRMS());
275 return kTRUE;
0c76cfa4 276 case kFigNtrackletsTrack:
1ee39b3a 277 h=MakePlotNTracklets();
278 PutTrendValue("NTrackletsTrack", h->GetMean());
279 PutTrendValue("NTrackletsTrackRMS", h->GetRMS());
280 return kTRUE;
0c76cfa4 281 case kFigNTrackletsP:
282 MakePlotnTrackletsVsP();
283 return kTRUE;
284 case kFigNtrackletsCross:
1ee39b3a 285 h = (TH1F*)fContainer->FindObject("hNtlsCross");
286 if(!MakeBarPlot(h, kRed)) break;
287 PutTrendValue("NTrackletsCross", h->GetMean());
288 PutTrendValue("NTrackletsCrossRMS", h->GetRMS());
289 return kTRUE;
0c76cfa4 290 case kFigNtrackletsFindable:
1ee39b3a 291 h = (TH1F*)fContainer->FindObject("hNtlsFindable");
292 if(!MakeBarPlot(h, kGreen)) break;
293 PutTrendValue("NTrackletsFindable", h->GetMean());
294 PutTrendValue("NTrackletsFindableRMS", h->GetRMS());
295 return kTRUE;
0c76cfa4 296 case kFigNtracksEvent:
1ee39b3a 297 (h = (TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
298 PutTrendValue("NTracksEvent", h->GetMean());
299 PutTrendValue("NTracksEventRMS", h->GetRMS());
300 return kTRUE;
0c76cfa4 301 case kFigNtracksSector:
1ee39b3a 302 h = (TH1F*)fContainer->FindObject("hNtrksSector");
303 if(!MakeBarPlot(h, kGreen)) break;
304 PutTrendValue("NTracksSector", h->Integral()/h->GetNbinsX());
305 return kTRUE;
0c76cfa4 306 case kFigTrackStatus:
1ee39b3a 307 if(!(h=(TH1F *)fContainer->FindObject("hTrackStatus"))) break;
308 h->GetXaxis()->SetRangeUser(0.5, -1);
309 h->GetYaxis()->CenterTitle();
310 h->Draw("c");
311 PutTrendValue("TrackStatus", h->Integral());
312 gPad->SetLogy(0);
313 return kTRUE;
0c76cfa4 314 case kFigTrackletStatus:
1ee39b3a 315 if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))) break;
147c3968 316 leg = new TLegend(.68, .7, .97, .97);
317 leg->SetBorderSize(0);leg->SetFillStyle(0);
1ee39b3a 318 leg->SetHeader("TRD layer");
147c3968 319 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
1ee39b3a 320 if(!(h=dynamic_cast<TH1F*>(arr->At(ily)))) continue;
321 if(kFIRST){
147c3968 322 h->Draw("pl");
1ee39b3a 323 h->GetXaxis()->SetRangeUser(0.5, -1);
324 h->GetYaxis()->CenterTitle();
325 kFIRST = kFALSE;
147c3968 326 } else h->Draw("samepl");
327 leg->AddEntry(h, Form("ly = %d", ily), "l");
1ee39b3a 328 PutTrendValue(Form("TrackletStatus%d", ily), h->Integral());
329 }
330 leg->Draw();
331 gPad->SetLogy(0);
332 return kTRUE;
0c76cfa4 333 case kFigChi2:
1ee39b3a 334 MakePlotChi2();
335 return kTRUE;
0c76cfa4 336 case kFigPH:
94f7dff7 337 gPad->SetMargin(0.125, 0.015, 0.1, 0.1);
1ee39b3a 338 MakePlotPulseHeight();
339 gPad->SetLogy(0);
340 return kTRUE;
0c76cfa4 341 case kFigChargeCluster:
1ee39b3a 342 (h = (TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
343 gPad->SetLogy(1);
344 PutTrendValue("ChargeCluster", h->GetMaximumBin());
345 PutTrendValue("ChargeClusterRMS", h->GetRMS());
346 return kTRUE;
0c76cfa4 347 case kFigChargeTracklet:
1ee39b3a 348 (h=(TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
349 PutTrendValue("ChargeTracklet", h->GetMaximumBin());
350 PutTrendValue("ChargeTrackletRMS", h->GetRMS());
351 return kTRUE;
0c76cfa4 352 case kFigNeventsTrigger:
1ee39b3a 353 ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
354 return kTRUE;
0c76cfa4 355 case kFigNeventsTriggerTracks:
1ee39b3a 356 ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
357 return kTRUE;
0c76cfa4 358 case kFigTriggerPurity:
1ee39b3a 359 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
360 break;
361 default:
362 break;
363 }
364 AliInfo(Form("Reference plot [%d] missing result", ifig));
365 return kFALSE;
366}
367
368//_______________________________________________________
369TObjArray *AliTRDcheckDET::Histos(){
370 //
371 // Create QA histograms
372 //
f8f46e4d 373
1ee39b3a 374 if(fContainer) return fContainer;
375
376 fContainer = new TObjArray(20);
377 //fContainer->SetOwner(kTRUE);
378
379 // Register Histograms
380 TH1 * h = NULL;
381 TAxis *ax = NULL;
382 if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
383 h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
384 h->GetXaxis()->SetTitle("N_{clusters}");
385 h->GetYaxis()->SetTitle("Entries");
386 } else h->Reset();
387 fContainer->AddAt(h, kNclustersTrack);
388
389 if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
390 h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
391 h->GetXaxis()->SetTitle("N_{clusters}");
392 h->GetYaxis()->SetTitle("Entries");
393 } else h->Reset();
394 fContainer->AddAt(h, kNclustersTracklet);
395
396 if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
397 h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
398 h->GetXaxis()->SetTitle("N^{tracklet}");
399 h->GetYaxis()->SetTitle("freq. [%]");
400 } else h->Reset();
401 fContainer->AddAt(h, kNtrackletsTrack);
402
403 if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
404 h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
405 h->GetXaxis()->SetTitle("N^{tracklet}");
406 h->GetYaxis()->SetTitle("freq. [%]");
407 }
408 fContainer->AddAt(h, kNtrackletsSTA);
409
a5d9fe6f 410 // Binning for momentum dependent tracklet Plots
de791ef5 411 const Int_t kNp(30);
412 Float_t P=0.2;
147c3968 413 Float_t binsP[kNp+1], binsTrklt[AliTRDgeometry::kNlayer+1];
d239caab 414 for(Int_t i=0;i<kNp+1; i++,P+=(TMath::Exp(i*i*.001)-1.)) binsP[i]=P;
147c3968 415 for(Int_t il = 0; il <= AliTRDgeometry::kNlayer; il++) binsTrklt[il] = 0.5 + il;
1ee39b3a 416 if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
a5d9fe6f 417 // Make tracklets for barrel tracks momentum dependent (if we do not exceed min and max values)
147c3968 418 h = new TH2F("hNtlsBAR",
419 "N_{tracklets} / track;p [GeV/c];N^{tracklet};freq. [%]",
420 kNp, binsP, AliTRDgeometry::kNlayer, binsTrklt);
1ee39b3a 421 }
422 fContainer->AddAt(h, kNtrackletsBAR);
423
424 //
425 if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
426 h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
427 h->GetXaxis()->SetTitle("n_{row cross}");
428 h->GetYaxis()->SetTitle("freq. [%]");
429 } else h->Reset();
430 fContainer->AddAt(h, kNtrackletsCross);
431
432 if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
433 h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
434 h->GetXaxis()->SetTitle("r [a.u]");
435 h->GetYaxis()->SetTitle("Entries");
436 } else h->Reset();
437 fContainer->AddAt(h, kNtrackletsFindable);
438
439 if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
440 h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
441 h->GetXaxis()->SetTitle("N_{tracks}");
442 h->GetYaxis()->SetTitle("Entries");
443 } else h->Reset();
444 fContainer->AddAt(h, kNtracksEvent);
445
446 if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
447 h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
448 h->GetXaxis()->SetTitle("sector");
449 h->GetYaxis()->SetTitle("freq. [%]");
450 } else h->Reset();
451 fContainer->AddAt(h, kNtracksSector);
452
453 if(!(h = (TH1F*)gROOT->FindObject("hTrackStatus"))){
454 const Int_t nerr = 7;
455 h = new TH1F("hTrackStatus", "Track Status", nerr, -0.5, nerr-0.5);
456 const Char_t *label[nerr] = {"OK", "PROL", "PROP", "AJST", "SNP", "TINI", "UPDT"};
457 ax = h->GetXaxis();
458 for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
459 h->SetYTitle("Relative Error to Good [%]");
460 }
461 fContainer->AddAt(h, kTrackStatus);
462
463 TObjArray *arr = new TObjArray(AliTRDgeometry::kNlayer);
464 arr->SetOwner(kTRUE); arr->SetName("TrackletStatus");
465 fContainer->AddAt(arr, kTrackletStatus);
466 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
467 if(!(h = (TH1F *)gROOT->FindObject(Form("hTrackletStatus%d", ily)))){
468 const Int_t nerr = 8;
469 h = new TH1F(Form("hTrackletStatus%d", ily), "Tracklet status", nerr, -0.5, nerr-0.5);
470 h->SetLineColor(ily+1);
471 const Char_t *label[nerr] = {"OK", "Geom", "Bound", "NoCl", "NoAttach", "NoClTr", "NoFit", "Chi2"};
472 ax = h->GetXaxis();
473 for(Int_t ierr=nerr; ierr--;) ax->SetBinLabel(ierr+1, label[ierr]);
474 h->SetYTitle("Relative Error to Good [%]");
475 } else h->Reset();
476 arr->AddAt(h, ily);
477 }
478
479 // <PH> histos
480 arr = new TObjArray(2);
481 arr->SetOwner(kTRUE); arr->SetName("<PH>");
482 fContainer->AddAt(arr, kPH);
483 if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
484 h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
485 h->GetXaxis()->SetTitle("Time / 100ns");
486 h->GetYaxis()->SetTitle("<PH> [a.u]");
487 } else h->Reset();
488 arr->AddAt(h, 0);
489 if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
490 h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
491 else h->Reset();
492 arr->AddAt(h, 1);
493
494 // Chi2 histos
495 if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
496 h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
497 h->SetXTitle("ndf");
498 h->SetYTitle("#chi^{2}/ndf");
499 h->SetZTitle("entries");
500 } else h->Reset();
501 fContainer->AddAt(h, kChi2);
502
503 if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
504 h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
505 h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
506 h->GetYaxis()->SetTitle("Entries");
507 }else h->Reset();
508 fContainer->AddAt(h, kChargeCluster);
509
510 if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
511 h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
512 h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
513 h->GetYaxis()->SetTitle("Entries");
514 }else h->Reset();
515 fContainer->AddAt(h, kChargeTracklet);
516
517
518 if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
519 h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
520 else h->Reset();
f8f46e4d 521 printf("Histos Adding \n");
522
1ee39b3a 523 fContainer->AddAt(h, kNeventsTrigger);
524
525 if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
526 h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
527 else h->Reset();
528 fContainer->AddAt(h, kNeventsTriggerTracks);
529
530 if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
531 h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
532 h->GetXaxis()->SetTitle("Trigger Cluster");
533 h->GetYaxis()->SetTitle("freq.");
534 } else h->Reset();
535 fContainer->AddAt(h, kTriggerPurity);
536
537 return fContainer;
538}
539
540/*
541* Plotting Functions
542*/
543
544//_______________________________________________________
545TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track)
546{
547//
548// Plot the track propagation status. The following errors are defined (see AliTRDtrackV1::ETRDtrackError)
549// PROL - track prolongation failure
550// PROP - track propagation failure
551// AJST - crossing sectors failure
552// SNP - too large bending
553// TINI - tracklet initialization failure
554// UPDT - track position/covariance update failure
555//
556// Performance plot looks as below:
557//Begin_Html
558//<img src="TRD/trackStatus.gif">
559//End_Html
560//
561 if(track) fkTrack = track;
562 if(!fkTrack){
563 AliWarning("No Track defined.");
db99a57a 564 return NULL;
1ee39b3a 565 }
db99a57a 566 TH1 *h = NULL;
1ee39b3a 567 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kTrackStatus)))){
568 AliWarning("No Histogram defined.");
db99a57a 569 return NULL;
1ee39b3a 570 }
571 h->Fill(fkTrack->GetStatusTRD());
572 return h;
573}
574
575//_______________________________________________________
576TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track)
577{
578//
579// Plot the tracklet propagation status. The following errors are defined for tracklet (see AliTRDtrackV1::ETRDlayerError)
580// Geom -
581// Bound - tracklet too close to chamber walls
582// NoCl - no clusters in the track roads
583// NoAttach - fail to attach clusters
584// NoClTr - fail to use clusters for fit
585// NoFit - tracklet fit failled
586// Chi2 - chi2 tracklet-track over threshold
587//
588// Performance plot looks as below:
589//Begin_Html
590//<img src="TRD/trackletStatus.gif">
591//End_Html
592//
593 if(track) fkTrack = track;
594 if(!fkTrack){
595 AliWarning("No Track defined.");
db99a57a 596 return NULL;
1ee39b3a 597 }
db99a57a 598 TObjArray *arr =NULL;
1ee39b3a 599 if(!(arr = dynamic_cast<TObjArray*>(fContainer->At(kTrackletStatus)))){
600 AliWarning("Histograms not defined.");
db99a57a 601 return NULL;
1ee39b3a 602 }
603
db99a57a 604 TH1 *h = NULL;
1ee39b3a 605 for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){
606 if(!(h = dynamic_cast<TH1F*>(arr->At(ily)))){
607 AliWarning(Form("Missing histo for layer %d.", ily));
608 continue;
609 }
610 h->Fill(fkTrack->GetStatusTRD(ily));
611 }
612 return h;
613}
614
615//_______________________________________________________
616TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
617 //
618 // Plot the mean number of clusters per tracklet
619 //
620 if(track) fkTrack = track;
621 if(!fkTrack){
622 AliWarning("No Track defined.");
db99a57a 623 return NULL;
1ee39b3a 624 }
db99a57a 625 TH1 *h = NULL;
1ee39b3a 626 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
627 AliWarning("No Histogram defined.");
db99a57a 628 return NULL;
1ee39b3a 629 }
db99a57a 630 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 631 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
632 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
633 h->Fill(tracklet->GetN2());
634 }
635 return h;
636}
637
638//_______________________________________________________
639TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
640 //
641 // Plot the number of clusters in one track
642 //
643 if(track) fkTrack = track;
644 if(!fkTrack){
645 AliWarning("No Track defined.");
db99a57a 646 return NULL;
1ee39b3a 647 }
db99a57a 648 TH1 *h = NULL;
1ee39b3a 649 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
650 AliWarning("No Histogram defined.");
db99a57a 651 return NULL;
1ee39b3a 652 }
653
654 Int_t nclusters = 0;
db99a57a 655 AliTRDseedV1 *tracklet = NULL;
a310e49b 656 AliExternalTrackParam *par = fkTrack->GetTrackOut() ? fkTrack->GetTrackOut() : fkTrack->GetTrackIn();
0c76cfa4 657 Double_t momentumRec = par->P();
1ee39b3a 658 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
659 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
fc2ec42c 660 Int_t n(tracklet->GetN());
661 nclusters += n;
1ee39b3a 662 if(DebugLevel() > 2){
663 Int_t crossing = Int_t(tracklet->IsRowCross());
664 Int_t detector = tracklet->GetDetector();
665 Float_t theta = TMath::ATan(tracklet->GetZref(1));
666 Float_t phi = TMath::ATan(tracklet->GetYref(1));
0c76cfa4 667 Float_t momentumMC = 0.;
1ee39b3a 668 Int_t pdg = 0;
669 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
670 UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
671 if(fkMC){
0c76cfa4 672 if(fkMC->GetTrackRef()) momentumMC = fkMC->GetTrackRef()->P();
1ee39b3a 673 pdg = fkMC->GetPDG();
674 }
675 (*DebugStream()) << "NClustersTrack"
676 << "Detector=" << detector
677 << "crossing=" << crossing
0c76cfa4 678 << "momentumMC="<< momentumMC
679 << "momentumRec=" << momentumRec
1ee39b3a 680 << "pdg=" << pdg
681 << "theta=" << theta
682 << "phi=" << phi
683 << "kinkIndex=" << kinkIndex
684 << "TPCncls=" << nclsTPC
fc2ec42c 685 << "TRDncls=" << n
1ee39b3a 686 << "\n";
687 }
688 }
689 h->Fill(nclusters);
690 return h;
691}
692
693
694//_______________________________________________________
695TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
696 //
697 // Plot the number of tracklets
698 //
699 if(track) fkTrack = track;
700 if(!fkTrack){
701 AliWarning("No Track defined.");
db99a57a 702 return NULL;
1ee39b3a 703 }
a5d9fe6f 704 TH1 *h = NULL, *hSta = NULL; TH2 *hBarrel = NULL;
1ee39b3a 705 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
706 AliWarning("No Histogram defined.");
db99a57a 707 return NULL;
1ee39b3a 708 }
709 Int_t nTracklets = fkTrack->GetNumberOfTracklets();
710 h->Fill(nTracklets);
711 if(!fkESD) return h;
712 Int_t status = fkESD->GetStatus();
a5d9fe6f 713
1ee39b3a 714/* 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 715 Double_t p = 0.;
0c76cfa4 716 Int_t method = -1; // to distinguish between stand alone and full barrel tracks in the debugging
1ee39b3a 717 if((status & AliESDtrack::kTRDin) != 0){
0c76cfa4 718 method = 1;
a5d9fe6f 719 // Full Barrel Track: Save momentum dependence
720 if(!(hBarrel = dynamic_cast<TH2F *>(fContainer->At(kNtrackletsBAR)))){
1ee39b3a 721 AliWarning("Method: Barrel. Histogram not processed!");
147c3968 722 return NULL;
a5d9fe6f 723 }
147c3968 724 AliExternalTrackParam *par(fkTrack->GetTrackIn());
725 if(!par){
726 AliError("Input track params missing");
727 return NULL;
728 }
729 hBarrel->Fill(par->P(), nTracklets);
1ee39b3a 730 } else {
a5d9fe6f 731 // Stand alone Track: momentum dependence not usefull
0c76cfa4 732 method = 0;
a5d9fe6f 733 if(!(hSta = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA)))) {
1ee39b3a 734 AliWarning("Method: StandAlone. Histogram not processed!");
147c3968 735 return NULL;
a5d9fe6f 736 }
147c3968 737 hSta->Fill(nTracklets);
1ee39b3a 738 }
1ee39b3a 739
0c76cfa4 740 if(DebugLevel() > 2){
741 AliTRDseedV1 *tracklet = NULL;
742 Int_t sector = -1;
743 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
744 if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
745 sector = fGeo->GetSector(tracklet->GetDetector());
746 break;
747 }
748 (*DebugStream()) << "NTrackletsTrack"
749 << "Sector=" << sector
750 << "NTracklets=" << nTracklets
751 << "Method=" << method
752 << "p=" << p
753 << "\n";
754 }
1ee39b3a 755 if(DebugLevel() > 3){
756 if(nTracklets == 1){
757 // If we have one Tracklet, check in which layer this happens
758 Int_t layer = -1;
db99a57a 759 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 760 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
761 if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){layer = il; break;}
762 }
a5d9fe6f 763 if(layer >= 0){
0c76cfa4 764 (*DebugStream()) << "NTrackletsLayer"
a5d9fe6f 765 << "Layer=" << layer
766 << "p=" << p
767 << "\n";
768 }
1ee39b3a 769 }
770 }
771 return h;
772}
773
774
775//_______________________________________________________
776TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
777 //
778 // Plot the number of tracklets
779 //
780 if(track) fkTrack = track;
781 if(!fkTrack){
782 AliWarning("No Track defined.");
db99a57a 783 return NULL;
1ee39b3a 784 }
db99a57a 785 TH1 *h = NULL;
1ee39b3a 786 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
787 AliWarning("No Histogram defined.");
db99a57a 788 return NULL;
1ee39b3a 789 }
790
791 Int_t ncross = 0;
db99a57a 792 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 793 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
794 if(!(tracklet = fkTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
795
796 if(tracklet->IsRowCross()) ncross++;
797 }
798 h->Fill(ncross);
799 return h;
800}
801
802//_______________________________________________________
803TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
804 //
805 // Plots the ratio of number of tracklets vs.
806 // number of findable tracklets
807 //
808 // Findable tracklets are defined as track prolongation
809 // to layer i does not hit the dead area +- epsilon
810 //
811 // In order to check whether tracklet hist active area in Layer i,
812 // the track is refitted and the fitted position + an uncertainty
813 // range is compared to the chamber border (also with a different
814 // uncertainty)
815 //
816 // For the track fit two cases are distinguished:
817 // If the track is a stand alone track (defined by the status bit
818 // encoding, then the track is fitted with the tilted Rieman model
819 // Otherwise the track is fitted with the Kalman fitter in two steps:
820 // Since the track parameters are give at the outer point, we first
821 // fit in direction inwards. Afterwards we fit again in direction outwards
822 // to extrapolate the track to layers which are not reached by the track
823 // For the Kalman model, the radial track points have to be shifted by
824 // a distance epsilon in the direction that we want to fit
825 //
826 const Float_t epsilon = 0.01; // dead area tolerance
827 const Float_t epsilonR = 1; // shift in radial direction of the anode wire position (Kalman filter only)
828 const Float_t deltaY = 0.7; // Tolerance in the track position in y-direction
829 const Float_t deltaZ = 7.0; // Tolerance in the track position in z-direction (Padlength)
830 Double_t xAnode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
831
832 if(track) fkTrack = track;
833 if(!fkTrack){
834 AliWarning("No Track defined.");
db99a57a 835 return NULL;
1ee39b3a 836 }
db99a57a 837 TH1 *h = NULL;
1ee39b3a 838 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
839 AliWarning("No Histogram defined.");
db99a57a 840 return NULL;
1ee39b3a 841 }
842 Int_t nFound = 0, nFindable = 0;
843 Int_t stack = -1;
844 Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
845 Double_t y = 0., z = 0.;
db99a57a 846 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 847 AliTRDpadPlane *pp;
848 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
849 if((tracklet = fkTrack->GetTracklet(il)) && tracklet->IsOK()){
850 tracklet->SetReconstructor(fReconstructor);
851 nFound++;
852 }
853 }
854 // 2 Different cases:
855 // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
856 // 2nd barrel track: here we propagate the track to the layers
857 AliTrackPoint points[6];
858 Float_t xyz[3];
859 memset(xyz, 0, sizeof(Float_t) * 3);
860 if(((fkESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fkESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
861 // stand alone track
862 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
863 xyz[0] = xAnode[il];
864 points[il].SetXYZ(xyz);
865 }
db99a57a 866 AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fkTrack), NULL, kTRUE, 6, points);
1ee39b3a 867 } else {
868 // barrel track
869 //
870 // 2 Steps:
871 // -> Kalman inwards
872 // -> Kalman outwards
873 AliTRDtrackV1 copyTrack(*fkTrack); // Do Kalman on a (non-constant) copy of the track
874 AliTrackPoint pointsInward[6], pointsOutward[6];
875 for(Int_t il = AliTRDgeometry::kNlayer; il--;){
876 // In order to avoid complications in the Kalman filter if the track points have the same radial
877 // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
878 // in the direction we want to go
879 // The track points have to be in reverse order for the Kalman Filter inwards
880 xyz[0] = xAnode[AliTRDgeometry::kNlayer - il - 1] - epsilonR;
881 pointsInward[il].SetXYZ(xyz);
882 xyz[0] = xAnode[il] + epsilonR;
883 pointsOutward[il].SetXYZ(xyz);
884 }
885 /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
886 printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
887 // Kalman inwards
db99a57a 888 AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kFALSE, 6, pointsInward);
1ee39b3a 889 memcpy(points, pointsInward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
890 // Kalman outwards
db99a57a 891 AliTRDtrackerV1::FitKalman(&copyTrack, NULL, kTRUE, 6, pointsInward);
1ee39b3a 892 memcpy(points, pointsOutward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
893 }
894 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
895 y = points[il].GetY();
896 z = points[il].GetZ();
897 if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
898 pp = fGeo->GetPadPlane(il, stack);
899 ymin = pp->GetCol0() + epsilon;
900 ymax = pp->GetColEnd() - epsilon;
901 zmin = pp->GetRowEnd() + epsilon;
902 zmax = pp->GetRow0() - epsilon;
903 // ignore y-crossing (material)
904 if((z + deltaZ > zmin && z - deltaZ < zmax) && (y + deltaY > ymin && y - deltaY < ymax)) nFindable++;
905 if(DebugLevel() > 3){
906 Double_t posTracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
907 Int_t hasTracklet = tracklet ? 1 : 0;
908 (*DebugStream()) << "FindableTracklets"
909 << "layer=" << il
910 << "ytracklet=" << posTracklet[0]
911 << "ytrack=" << y
912 << "ztracklet=" << posTracklet[1]
913 << "ztrack=" << z
914 << "tracklet=" << hasTracklet
915 << "\n";
916 }
917 }
918
919 h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
920 AliDebug(2, Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
921 return h;
922}
923
924
925//_______________________________________________________
926TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
927 //
928 // Plot the chi2 of the track
929 //
930 if(track) fkTrack = track;
931 if(!fkTrack){
932 AliWarning("No Track defined.");
db99a57a 933 return NULL;
1ee39b3a 934 }
db99a57a 935 TH1 *h = NULL;
1ee39b3a 936 if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
937 AliWarning("No Histogram defined.");
db99a57a 938 return NULL;
1ee39b3a 939 }
940 Int_t n = fkTrack->GetNumberOfTracklets();
db99a57a 941 if(!n) return NULL;
1ee39b3a 942
943 h->Fill(n, fkTrack->GetChi2()/n);
944 return h;
945}
946
947
948//_______________________________________________________
949TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
950 //
951 // Plot the average pulse height
952 //
953 if(track) fkTrack = track;
954 if(!fkTrack){
955 AliWarning("No Track defined.");
db99a57a 956 return NULL;
1ee39b3a 957 }
db99a57a 958 TProfile *h = NULL;
1ee39b3a 959 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
960 AliWarning("No Histogram defined.");
db99a57a 961 return NULL;
1ee39b3a 962 }
db99a57a 963 AliTRDseedV1 *tracklet = NULL;
964 AliTRDcluster *c = NULL;
1ee39b3a 965 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
966 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
967 Int_t crossing = Int_t(tracklet->IsRowCross());
968 Int_t detector = tracklet->GetDetector();
969 tracklet->ResetClusterIter();
970 while((c = tracklet->NextCluster())){
c732f879 971 if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1ee39b3a 972 Int_t localtime = c->GetLocalTimeBin();
973 Double_t absoluteCharge = TMath::Abs(c->GetQ());
974 h->Fill(localtime, absoluteCharge);
975 if(DebugLevel() > 3){
976 Double_t distance[2];
977 GetDistanceToTracklet(distance, tracklet, c);
978 Float_t theta = TMath::ATan(tracklet->GetZref(1));
979 Float_t phi = TMath::ATan(tracklet->GetYref(1));
980 Float_t momentum = 0.;
981 Int_t pdg = 0;
982 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
983 UShort_t TPCncls = fkESD ? fkESD->GetTPCncls() : 0;
984 if(fkMC){
985 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
986 pdg = fkMC->GetPDG();
987 }
988 (*DebugStream()) << "PHt"
989 << "Detector=" << detector
990 << "crossing=" << crossing
991 << "Timebin=" << localtime
992 << "Charge=" << absoluteCharge
993 << "momentum=" << momentum
994 << "pdg=" << pdg
995 << "theta=" << theta
996 << "phi=" << phi
997 << "kinkIndex=" << kinkIndex
998 << "TPCncls=" << TPCncls
999 << "dy=" << distance[0]
1000 << "dz=" << distance[1]
1001 << "c.=" << c
1002 << "\n";
1003 }
1004 }
1005 }
1006 return h;
1007}
1008
1009//_______________________________________________________
1010TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
1011 //
1012 // Plots the average pulse height vs the distance from the anode wire
1013 // (plus const anode wire offset)
1014 //
1015 if(track) fkTrack = track;
1016 if(!fkTrack){
1017 AliWarning("No Track defined.");
db99a57a 1018 return NULL;
1ee39b3a 1019 }
db99a57a 1020 TProfile *h = NULL;
1ee39b3a 1021 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
1022 AliWarning("No Histogram defined.");
db99a57a 1023 return NULL;
1ee39b3a 1024 }
1025 Float_t offset = .5*AliTRDgeometry::CamHght();
db99a57a 1026 AliTRDseedV1 *tracklet = NULL;
1027 AliTRDcluster *c = NULL;
1ee39b3a 1028 Double_t distance = 0;
1029 Double_t x, y;
1030 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1031 if(!(tracklet = fkTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
1032 tracklet->ResetClusterIter();
1033 while((c = tracklet->NextCluster())){
c732f879 1034 if(!IsUsingClustersOutsideChamber() && !c->IsInChamber()) continue;
1ee39b3a 1035 x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
1036 y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
1037
1038 distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
1039 h->Fill(distance, TMath::Abs(c->GetQ()));
1040 }
1041 }
1042 return h;
1043}
1044
1045//_______________________________________________________
1046TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
1047 //
1048 // Plot the cluster charge
1049 //
1050 if(track) fkTrack = track;
1051 if(!fkTrack){
1052 AliWarning("No Track defined.");
db99a57a 1053 return NULL;
1ee39b3a 1054 }
db99a57a 1055 TH1 *h = NULL;
1ee39b3a 1056 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
1057 AliWarning("No Histogram defined.");
db99a57a 1058 return NULL;
1ee39b3a 1059 }
db99a57a 1060 AliTRDseedV1 *tracklet = NULL;
1061 AliTRDcluster *c = NULL;
1ee39b3a 1062 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1063 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
1064 for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
1065 if(!(c = tracklet->GetClusters(itime))) continue;
1066 h->Fill(c->GetQ());
1067 }
1068 }
1069 return h;
1070}
1071
1072//_______________________________________________________
1073TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
1074 //
1075 // Plot the charge deposit per chamber
1076 //
1077 if(track) fkTrack = track;
1078 if(!fkTrack){
1079 AliWarning("No Track defined.");
db99a57a 1080 return NULL;
1ee39b3a 1081 }
db99a57a 1082 TH1 *h = NULL;
1ee39b3a 1083 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
1084 AliWarning("No Histogram defined.");
db99a57a 1085 return NULL;
1ee39b3a 1086 }
db99a57a 1087 AliTRDseedV1 *tracklet = NULL;
1088 AliTRDcluster *c = NULL;
1ee39b3a 1089 Double_t qTot = 0;
1090 Int_t nTracklets =fkTrack->GetNumberOfTracklets();
8428f55d 1091 for(Int_t itl(0); itl < AliTRDgeometry::kNlayer; itl++){
1ee39b3a 1092 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1093 qTot = 0.;
1094 for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
1095 if(!(c = tracklet->GetClusters(ic))) continue;
1096 qTot += TMath::Abs(c->GetQ());
1097 }
1098 h->Fill(qTot);
1099 if(DebugLevel() > 3){
1100 Int_t crossing = (Int_t)tracklet->IsRowCross();
1101 Int_t detector = tracklet->GetDetector();
1102 Float_t theta = TMath::ATan(tracklet->GetZfit(1));
1103 Float_t phi = TMath::ATan(tracklet->GetYfit(1));
1104 Float_t momentum = 0.;
1105 Int_t pdg = 0;
1106 Int_t kinkIndex = fkESD ? fkESD->GetKinkIndex() : 0;
1107 UShort_t nclsTPC = fkESD ? fkESD->GetTPCncls() : 0;
1108 if(fkMC){
1109 if(fkMC->GetTrackRef()) momentum = fkMC->GetTrackRef()->P();
1110 pdg = fkMC->GetPDG();
1111 }
1112 (*DebugStream()) << "ChargeTracklet"
1113 << "Detector=" << detector
1114 << "crossing=" << crossing
1115 << "momentum=" << momentum
1116 << "nTracklets="<< nTracklets
1117 << "pdg=" << pdg
1118 << "theta=" << theta
1119 << "phi=" << phi
1120 << "kinkIndex=" << kinkIndex
1121 << "TPCncls=" << nclsTPC
1122 << "QT=" << qTot
1123 << "\n";
1124 }
1125 }
1126 return h;
1127}
1128
1129//_______________________________________________________
1130TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
1131 //
1132 // Plot the number of tracks per Sector
1133 //
1134 if(track) fkTrack = track;
1135 if(!fkTrack){
1136 AliWarning("No Track defined.");
db99a57a 1137 return NULL;
1ee39b3a 1138 }
db99a57a 1139 TH1 *h = NULL;
1ee39b3a 1140 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
1141 AliWarning("No Histogram defined.");
db99a57a 1142 return NULL;
1ee39b3a 1143 }
1144
1145 // TODO we should compare with
1146 // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
1147
db99a57a 1148 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 1149 Int_t sector = -1;
1150 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
1151 if(!(tracklet = fkTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
1152 sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
1153 break;
1154 }
1155 h->Fill(sector);
1156 return h;
1157}
1158
1159
1160//________________________________________________________
1161void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
1162{
1163
1164 fReconstructor->SetRecoParam(r);
1165}
1166
1167//________________________________________________________
1168void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 * const tracklet, AliTRDcluster * const c)
1169{
1170 Float_t x = c->GetX();
1171 dist[0] = c->GetY() - tracklet->GetYat(x);
1172 dist[1] = c->GetZ() - tracklet->GetZat(x);
1173}
1174
1175
1176//_______________________________________________________
1177TH1* AliTRDcheckDET::MakePlotChi2()
1178{
1179// Plot chi2/track normalized to number of degree of freedom
1180// (tracklets) and compare with the theoretical distribution.
1181//
1182// Alex Bercuci <A.Bercuci@gsi.de>
1183
1184 TH2S *h2 = (TH2S*)fContainer->At(kChi2);
1185 TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
1186 TLegend *leg = new TLegend(.7,.7,.95,.95);
1187 leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
db99a57a 1188 TH1D *h1 = NULL;
1ee39b3a 1189 Bool_t kFIRST = kTRUE;
1190 for(Int_t il=1; il<=h2->GetNbinsX(); il++){
1191 h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
1192 if(h1->Integral()<50) continue;
1193 h1->Scale(1./h1->Integral());
1194 h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
1195 h1->SetLineColor(il);h1->SetLineStyle(2);
1196 f.SetParameter(1, .5*il);f.SetLineColor(il);
1197 h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
1198 leg->AddEntry(h1, Form("%d", il), "l");
1199 if(kFIRST){
1200 h1->GetXaxis()->SetRangeUser(0., 25.);
1201 }
1202 kFIRST = kFALSE;
1203 }
1204 leg->Draw();
1205 gPad->SetLogy();
1206 return h1;
1207}
1208
1209
1210//________________________________________________________
1211TH1* AliTRDcheckDET::MakePlotNTracklets(){
1212 //
1213 // Make nice bar plot of the number of tracklets in each method
1214 //
a5d9fe6f 1215 TH2F *tmp = (TH2F *)fContainer->FindObject("hNtlsBAR");
1216 TH1D *hBAR = tmp->ProjectionY();
1ee39b3a 1217 TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1218 TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1219 TLegend *leg = new TLegend(0.13, 0.75, 0.39, 0.89);
1220 leg->SetBorderSize(1);
1221 leg->SetFillColor(0);
1222
1223 Float_t scale = hCON->Integral();
1224 hCON->Scale(100./scale);
1225 hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1226 hCON->SetBarWidth(0.2);
1227 hCON->SetBarOffset(0.6);
1228 hCON->SetStats(kFALSE);
1229 hCON->GetYaxis()->SetRangeUser(0.,40.);
1230 hCON->GetYaxis()->SetTitleOffset(1.2);
1231 hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1232 hCON->SetMaximum(55.);
1233
1234 hBAR->Scale(100./scale);
1235 hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1236 hBAR->SetBarWidth(0.2);
1237 hBAR->SetBarOffset(0.2);
1238 hBAR->SetTitle("");
1239 hBAR->SetStats(kFALSE);
1240 hBAR->GetYaxis()->SetRangeUser(0.,40.);
1241 hBAR->GetYaxis()->SetTitleOffset(1.2);
1242 hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1243
1244 hSTA->Scale(100./scale);
1245 hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1246 hSTA->SetBarWidth(0.2);
1247 hSTA->SetBarOffset(0.4);
1248 hSTA->SetTitle("");
1249 hSTA->SetStats(kFALSE);
1250 hSTA->GetYaxis()->SetRangeUser(0.,40.);
1251 hSTA->GetYaxis()->SetTitleOffset(1.2);
1252 hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1253 leg->Draw();
1254 gPad->Update();
1255 return hCON;
1256}
1257
0c76cfa4 1258//________________________________________________________
1259void AliTRDcheckDET::MakePlotnTrackletsVsP(){
1260 //
1261 // Plot abundance of tracks with number of tracklets as function of momentum
1262 //
147c3968 1263
1264
1265
1266
1267 Color_t color[AliTRDgeometry::kNlayer] = {kBlue, kOrange, kBlack, kGreen, kCyan, kRed};
1268 TH1 *h(NULL); TGraphErrors *g[AliTRDgeometry::kNlayer];
1269 for(Int_t itl(0); itl<AliTRDgeometry::kNlayer; itl++){
1270 g[itl] = new TGraphErrors();
1271 g[itl]->SetLineColor(color[itl]);
1272 g[itl]->SetMarkerColor(color[itl]);
1273 g[itl]->SetMarkerStyle(20 + itl);
1274 }
1275
0c76cfa4 1276 TH2 *hBar = (TH2F *)fContainer->FindObject("hNtlsBAR");
147c3968 1277 TAxis *ax(hBar->GetXaxis());
1278 Int_t np(ax->GetNbins());
1279 for(Int_t ipBin(1); ipBin<np; ipBin++){
1280 h = hBar->ProjectionY("npBin", ipBin, ipBin);
81979445 1281 if(!Int_t(h->Integral())) continue;
147c3968 1282 h->Scale(100./h->Integral());
1283 Float_t p(ax->GetBinCenter(ipBin));
1284 Float_t dp(ax->GetBinWidth(ipBin));
81979445 1285 Int_t ip(g[0]->GetN());
147c3968 1286 for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
81979445 1287 g[itl]->SetPoint(ip, p, h->GetBinContent(itl+1));
1288 g[itl]->SetPointError(ip, dp/3.46, h->GetBinError(itl+1));
147c3968 1289 }
1290 }
1291
81979445 1292 TLegend *leg = new TLegend(0.76, 0.6, 1., 0.9);
94f7dff7 1293 leg->SetBorderSize(0);
1294 leg->SetHeader("Tracklet/Track");
1295 leg->SetFillStyle(0);
147c3968 1296 h = hBar->ProjectionX("npxBin"); h->Reset();
81979445 1297 h->SetTitle("");
1298 h->GetYaxis()->SetRangeUser(1., 99.);
1299 h->GetYaxis()->SetMoreLogLabels();
1300 h->GetYaxis()->CenterTitle();
1301 h->GetYaxis()->SetTitleOffset(1.2);
147c3968 1302 h->SetYTitle("Prob. [%]");
81979445 1303 h->GetXaxis()->SetRangeUser(0.4, 12.);
1304 h->GetXaxis()->SetMoreLogLabels();
1305 h->GetXaxis()->CenterTitle();
147c3968 1306 h->Draw("p");
1307 for(Int_t itl(AliTRDgeometry::kNlayer); itl--;){
1308 g[itl]->Draw("pc");
1309 leg->AddEntry(g[itl], Form("n = %d", itl+1),"pl");
0c76cfa4 1310 }
147c3968 1311
0c76cfa4 1312 leg->Draw();
81979445 1313 gPad->SetLogx();gPad->SetLogy();
0c76cfa4 1314}
1315
1ee39b3a 1316//________________________________________________________
1317TH1* AliTRDcheckDET::MakePlotPulseHeight(){
1318 //
1319 // Create Plot of the Pluse Height Spectrum
1320 //
1321 TH1 *h, *h1, *h2;
1322 TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1323 h = (TH1F*)arr->At(0);
1324 h->SetMarkerStyle(24);
1325 h->SetMarkerColor(kBlack);
1326 h->SetLineColor(kBlack);
1327 h->Draw("e1");
3907f080 1328 // Trending for the pulse height: plateau value, slope and timebin of the maximum
1329 TLinearFitter fit(1,"pol1");
1330 Double_t time = 0.;
1331 for(Int_t itime = 10; itime <= 20; itime++){
1332 time = static_cast<Double_t>(itime);
1333 fit.AddPoint(&time, h->GetBinContent(itime + 1), h->GetBinError(itime + 1));
1334 }
1335 fit.Eval();
1336 Double_t plateau = fit.GetParameter(0) + 12 * fit.GetParameter(1);
1337 Double_t slope = fit.GetParameter(1);
1338 PutTrendValue("PHplateau", plateau);
1339 PutTrendValue("PHslope", slope);
1340 PutTrendValue("PHamplificationPeak", static_cast<Double_t>(h->GetMaximumBin()-1));
1341 AliDebug(1, Form("plateau %f, slope %f, MaxTime %f", plateau, slope, static_cast<Double_t>(h->GetMaximumBin()-1)));
1ee39b3a 1342// copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1343 h1 = (TH1F *)arr->At(1);
1344 h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1345 for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
1346 h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1347 h2->SetMarkerStyle(22);
1348 h2->SetMarkerColor(kBlue);
1349 h2->SetLineColor(kBlue);
1350 h2->Draw("e1same");
1351 gPad->Update();
1352// create axis according to the histogram dimensions of the original second histogram
1353 TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1354 gPad->GetUymax(),
1355 gPad->GetUxmax(),
1356 gPad->GetUymax(),
1357 -0.08, 4.88, 510,"-L");
1358 axis->SetLineColor(kBlue);
1359 axis->SetLabelColor(kBlue);
1360 axis->SetTextColor(kBlue);
1361 axis->SetTitle("x_{0}-x_{c} [cm]");
1362 axis->Draw();
1363 return h1;
1364}
1365
1366//________________________________________________________
1367Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1368 //
1369 // Draw nice bar plots
1370 //
1371 if(!histo->GetEntries()) return kFALSE;
1372 histo->Scale(100./histo->Integral());
1373 histo->SetFillColor(color);
1374 histo->SetBarOffset(.2);
1375 histo->SetBarWidth(.6);
1376 histo->Draw("bar1");
1377 return kTRUE;
960a59e0 1378}