12 #include <TObjArray.h>
14 #include <TObjString.h>
18 #include <TProfile2D.h>
22 #include "AliTRDcluster.h"
23 #include "AliESDHeader.h"
24 #include "AliESDRun.h"
25 #include "AliESDtrack.h"
26 #include "AliTRDgeometry.h"
27 #include "AliTRDpadPlane.h"
28 #include "AliTRDSimParam.h"
29 #include "AliTRDseedV1.h"
30 #include "AliTRDtrackV1.h"
31 #include "AliTRDtrackerV1.h"
32 #include "AliTRDReconstructor.h"
33 #include "AliTrackReference.h"
34 #include "AliTrackPointArray.h"
35 #include "AliTracker.h"
36 #include "TTreeStream.h"
38 #include "info/AliTRDtrackInfo.h"
39 #include "info/AliTRDeventInfo.h"
40 #include "AliTRDcheckDET.h"
45 ////////////////////////////////////////////////////////////////////////////
47 // Reconstruction QA //
49 // Task doing basic checks for tracking and detector performance //
52 // Anton Andronic <A.Andronic@gsi.de> //
53 // Alexandru Bercuci <A.Bercuci@gsi.de> //
54 // Markus Fasel <M.Fasel@gsi.de> //
56 ////////////////////////////////////////////////////////////////////////////
58 //_______________________________________________________
59 AliTRDcheckDET::AliTRDcheckDET():
60 AliTRDrecoTask("checkDET", "Basic TRD data checker")
67 // Default constructor
69 DefineInput(1,AliTRDeventInfo::Class());
70 fReconstructor = new AliTRDReconstructor;
71 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
72 fGeo = new AliTRDgeometry;
76 //_______________________________________________________
77 AliTRDcheckDET::~AliTRDcheckDET(){
81 if(fTriggerNames) delete fTriggerNames;
82 delete fReconstructor;
86 //_______________________________________________________
87 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
89 // Connect the Input data with the task
91 AliTRDrecoTask::ConnectInputData(opt);
92 fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
95 //_______________________________________________________
96 void AliTRDcheckDET::CreateOutputObjects(){
98 // Create Output Objects
100 OpenFile(0,"RECREATE");
101 fContainer = Histos();
102 if(!fTriggerNames) fTriggerNames = new TMap();
105 //_______________________________________________________
106 void AliTRDcheckDET::Exec(Option_t *opt){
108 // Execution function
109 // Filling TRD quality histos
111 if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
112 AliTRDrecoTask::Exec(opt);
113 Int_t nTracks = 0; // Count the number of tracks per event
114 Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
115 TString triggername = fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
116 if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
117 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
118 for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
119 if(!fTracks->UncheckedAt(iti)) continue;
120 AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
121 if(!fTrackInfo->GetTrack()) continue;
125 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
126 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
128 if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
129 fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
130 // also set the label for both histograms
131 TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
132 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
133 histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
134 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
136 PostData(0, fContainer);
139 //_______________________________________________________
140 void AliTRDcheckDET::Terminate(Option_t *){
142 // Terminate function
146 //_______________________________________________________
147 Bool_t AliTRDcheckDET::PostProcess(){
149 // Do Postprocessing (for the moment set the number of Reference histograms)
154 // Calculate of the trigger clusters purity
155 h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
156 TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
158 Float_t purities[20], val = 0;
159 TString triggernames[20];
160 Int_t nTriggerClasses = 0;
161 for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
162 if((val = h1->GetBinContent(ibin))){
163 purities[nTriggerClasses] = val;
164 triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
168 h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
169 TAxis *ax = h->GetXaxis();
170 for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
171 h->Fill(itrg, purities[itrg]);
172 ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
174 ax->SetRangeUser(-0.5, nTriggerClasses+.5);
175 h->GetYaxis()->SetRangeUser(0,1);
182 //_______________________________________________________
183 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
185 // Setting Reference Figures
190 case kNclustersTrack:
191 ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
193 case kNclustersTracklet:
194 ((TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
196 case kNtrackletsTrack:
197 MakePlotNTracklets();
199 case kNtrackletsCross:
200 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsCross"), kRed)) break;
202 case kNtrackletsFindable:
203 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsFindable"), kGreen)) break;
206 ((TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
209 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtrksSector"), kGreen)) break;
212 ((TH1I *)fContainer->FindObject("hTrackStatus"))->Draw("c");
214 case kTrackletStatus:
215 ((TH2S *)fContainer->FindObject("hTrackletStatus"))->Draw("colz");
221 MakePlotPulseHeight();
224 ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
227 case kChargeTracklet:
228 ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
230 case kNeventsTrigger:
231 ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
233 case kNeventsTriggerTracks:
234 ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
237 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
242 AliInfo(Form("Reference plot [%d] missing result", ifig));
246 //_______________________________________________________
247 TObjArray *AliTRDcheckDET::Histos(){
249 // Create QA histograms
251 if(fContainer) return fContainer;
253 fContainer = new TObjArray(20);
254 //fContainer->SetOwner(kTRUE);
256 // Register Histograms
259 if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
260 h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
261 h->GetXaxis()->SetTitle("N_{clusters}");
262 h->GetYaxis()->SetTitle("Entries");
264 fContainer->AddAt(h, kNclustersTrack);
266 if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
267 h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
268 h->GetXaxis()->SetTitle("N_{clusters}");
269 h->GetYaxis()->SetTitle("Entries");
271 fContainer->AddAt(h, kNclustersTracklet);
273 if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
274 h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
275 h->GetXaxis()->SetTitle("N^{tracklet}");
276 h->GetYaxis()->SetTitle("freq. [%]");
278 fContainer->AddAt(h, kNtrackletsTrack);
280 if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
281 h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
282 h->GetXaxis()->SetTitle("N^{tracklet}");
283 h->GetYaxis()->SetTitle("freq. [%]");
285 fContainer->AddAt(h, kNtrackletsSTA);
287 if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
288 h = new TH1F("hNtlsBAR", "N_{tracklets} / track (Barrel)", AliTRDgeometry::kNlayer, 0.5, 6.5);
289 h->GetXaxis()->SetTitle("N^{tracklet}");
290 h->GetYaxis()->SetTitle("freq. [%]");
292 fContainer->AddAt(h, kNtrackletsBAR);
295 if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
296 h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
297 h->GetXaxis()->SetTitle("n_{row cross}");
298 h->GetYaxis()->SetTitle("freq. [%]");
300 fContainer->AddAt(h, kNtrackletsCross);
302 if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
303 h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
304 h->GetXaxis()->SetTitle("r [a.u]");
305 h->GetYaxis()->SetTitle("Entries");
307 fContainer->AddAt(h, kNtrackletsFindable);
309 if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
310 h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
311 h->GetXaxis()->SetTitle("N_{tracks}");
312 h->GetYaxis()->SetTitle("Entries");
314 fContainer->AddAt(h, kNtracksEvent);
316 if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
317 h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
318 h->GetXaxis()->SetTitle("sector");
319 h->GetYaxis()->SetTitle("freq. [%]");
321 fContainer->AddAt(h, kNtracksSector);
323 if(!(h = (TH1I *)gROOT->FindObject("hTrackStatus"))){
324 h = new TH1I("hTrackStatus", "Track Status", 7,0,7);
325 axis = h->GetXaxis();
326 axis->SetBinLabel(axis->GetFirst() + 0, "OK");
327 axis->SetBinLabel(axis->GetFirst() + 1, "PROL");
328 axis->SetBinLabel(axis->GetFirst() + 2, "PROP");
329 axis->SetBinLabel(axis->GetFirst() + 3, "ADJ");
330 axis->SetBinLabel(axis->GetFirst() + 4, "SNP");
331 axis->SetBinLabel(axis->GetFirst() + 5, "TLIN");
332 axis->SetBinLabel(axis->GetFirst() + 6, "UP");
334 fContainer->AddAt(h, kTrackStatus);
336 if(!(h = (TH2S *)gROOT->FindObject("hTrackletStatus"))){
337 h = new TH2S("hTrackletStatus", "Tracklet status", 6,0,6,8,0,8);
338 axis = h->GetYaxis();
339 axis->SetBinLabel(axis->GetFirst() + 0, "OK");
340 axis->SetBinLabel(axis->GetFirst() + 1, "Geom");
341 axis->SetBinLabel(axis->GetFirst() + 2, "Boun");
342 axis->SetBinLabel(axis->GetFirst() + 3, "NoCl");
343 axis->SetBinLabel(axis->GetFirst() + 4, "NoAttach");
344 axis->SetBinLabel(axis->GetFirst() + 5, "NoClTr");
345 axis->SetBinLabel(axis->GetFirst() + 6, "NoFit");
346 axis->SetBinLabel(axis->GetFirst() + 7, "Chi2");
348 fContainer->AddAt(h, kTrackletStatus);
351 TObjArray *arr = new TObjArray(2);
352 arr->SetOwner(kTRUE); arr->SetName("<PH>");
353 fContainer->AddAt(arr, kPH);
354 if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
355 h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
356 h->GetXaxis()->SetTitle("Time / 100ns");
357 h->GetYaxis()->SetTitle("<PH> [a.u]");
360 if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
361 h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
366 if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
367 h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
369 h->SetYTitle("#chi^{2}/ndf");
370 h->SetZTitle("entries");
372 fContainer->AddAt(h, kChi2);
374 if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
375 h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
376 h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
377 h->GetYaxis()->SetTitle("Entries");
379 fContainer->AddAt(h, kChargeCluster);
381 if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
382 h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
383 h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
384 h->GetYaxis()->SetTitle("Entries");
386 fContainer->AddAt(h, kChargeTracklet);
389 if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
390 h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
392 fContainer->AddAt(h, kNeventsTrigger);
394 if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
395 h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
397 fContainer->AddAt(h, kNeventsTriggerTracks);
399 if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
400 h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
401 h->GetXaxis()->SetTitle("Trigger Cluster");
402 h->GetYaxis()->SetTitle("freq.");
404 fContainer->AddAt(h, kTriggerPurity);
413 //_______________________________________________________
414 TH1 *AliTRDcheckDET::PlotTrackStatus(const AliTRDtrackV1 *track){
416 // Plot the track status
418 if(track) fTrack = track;
420 AliWarning("No Track defined.");
424 if(!(h = dynamic_cast<TH1I *>(fContainer->At(kTrackStatus)))){
425 AliWarning("No Histogram defined.");
428 h->Fill(fTrack->GetStatusTRD());
432 //_______________________________________________________
433 TH1 *AliTRDcheckDET::PlotTrackletStatus(const AliTRDtrackV1 *track){
435 // Plot the track status
437 if(track) fTrack = track;
439 AliWarning("No Track defined.");
443 if(!(h = dynamic_cast<TH2S *>(fContainer->At(kTrackletStatus)))){
444 AliWarning("No Histogram defined.");
448 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
449 status = fTrack->GetStatusTRD(il);
455 //_______________________________________________________
456 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
458 // Plot the mean number of clusters per tracklet
460 if(track) fTrack = track;
462 AliWarning("No Track defined.");
466 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
467 AliWarning("No Histogram defined.");
470 AliTRDseedV1 *tracklet = 0x0;
471 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
472 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
473 h->Fill(tracklet->GetN2());
478 //_______________________________________________________
479 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
481 // Plot the number of clusters in one track
483 if(track) fTrack = track;
485 AliWarning("No Track defined.");
489 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
490 AliWarning("No Histogram defined.");
495 AliTRDseedV1 *tracklet = 0x0;
496 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
497 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
498 nclusters += tracklet->GetN();
500 Int_t crossing = Int_t(tracklet->IsRowCross());
501 Int_t detector = tracklet->GetDetector();
502 Float_t theta = TMath::ATan(tracklet->GetZref(1));
503 Float_t phi = TMath::ATan(tracklet->GetYref(1));
504 Float_t momentum = 0.;
506 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
507 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
509 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
512 (*fDebugStream) << "NClustersTrack"
513 << "Detector=" << detector
514 << "crossing=" << crossing
515 << "momentum=" << momentum
519 << "kinkIndex=" << kinkIndex
520 << "TPCncls=" << TPCncls
521 << "nclusters=" << nclusters
530 //_______________________________________________________
531 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
533 // Plot the number of tracklets
535 if(track) fTrack = track;
537 AliWarning("No Track defined.");
540 TH1 *h = 0x0, *hMethod = 0x0;
541 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
542 AliWarning("No Histogram defined.");
545 Int_t status = fESD->GetStatus();
546 /* 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);*/
547 if((status & AliESDtrack::kTRDin) != 0){
549 if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsBAR))))
550 AliWarning("Method: Barrel. Histogram not processed!");
553 if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA))))
554 AliWarning("Method: StandAlone. Histogram not processed!");
556 Int_t nTracklets = fTrack->GetNumberOfTracklets();
558 hMethod->Fill(nTracklets);
561 // If we have one Tracklet, check in which layer this happens
563 AliTRDseedV1 *tracklet = 0x0;
564 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
565 if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer = il; break;}
567 (*fDebugStream) << "NTrackletsTrack"
576 //_______________________________________________________
577 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
579 // Plot the number of tracklets
581 if(track) fTrack = track;
583 AliWarning("No Track defined.");
587 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
588 AliWarning("No Histogram defined.");
593 AliTRDseedV1 *tracklet = 0x0;
594 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
595 if(!(tracklet = fTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
597 if(tracklet->IsRowCross()) ncross++;
603 //_______________________________________________________
604 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
606 // Plots the ratio of number of tracklets vs.
607 // number of findable tracklets
609 // Findable tracklets are defined as track prolongation
610 // to layer i does not hit the dead area +- epsilon
612 // In order to check whether tracklet hist active area in Layer i,
613 // the track is refitted and the fitted position + an uncertainty
614 // range is compared to the chamber border (also with a different
617 // For the track fit two cases are distinguished:
618 // If the track is a stand alone track (defined by the status bit
619 // encoding, then the track is fitted with the tilted Rieman model
620 // Otherwise the track is fitted with the Kalman fitter in two steps:
621 // Since the track parameters are give at the outer point, we first
622 // fit in direction inwards. Afterwards we fit again in direction outwards
623 // to extrapolate the track to layers which are not reached by the track
624 // For the Kalman model, the radial track points have to be shifted by
625 // a distance epsilon in the direction that we want to fit
627 const Float_t epsilon = 0.01; // dead area tolerance
628 const Float_t epsilon_R = 1; // shift in radial direction of the anode wire position (Kalman filter only)
629 const Float_t delta_y = 0.7; // Tolerance in the track position in y-direction
630 const Float_t delta_z = 7.0; // Tolerance in the track position in z-direction (Padlength)
631 Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
633 if(track) fTrack = track;
635 AliWarning("No Track defined.");
639 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
640 AliWarning("No Histogram defined.");
643 Int_t nFound = 0, nFindable = 0;
645 Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
646 Double_t y = 0., z = 0.;
647 AliTRDseedV1 *tracklet = 0x0;
649 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
650 if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
651 tracklet->SetReconstructor(fReconstructor);
655 // 2 Different cases:
656 // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
657 // 2nd barrel track: here we propagate the track to the layers
658 AliTrackPoint points[6];
660 memset(xyz, 0, sizeof(Float_t) * 3);
661 if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
663 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
664 xyz[0] = x_anode[il];
665 points[il].SetXYZ(xyz);
667 AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
673 // -> Kalman outwards
674 AliTRDtrackV1 copy_track(*fTrack); // Do Kalman on a (non-constant) copy of the track
675 AliTrackPoint points_inward[6], points_outward[6];
676 for(Int_t il = AliTRDgeometry::kNlayer; il--;){
677 // In order to avoid complications in the Kalman filter if the track points have the same radial
678 // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
679 // in the direction we want to go
680 // The track points have to be in reverse order for the Kalman Filter inwards
681 xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
682 points_inward[il].SetXYZ(xyz);
683 xyz[0] = x_anode[il] + epsilon_R;
684 points_outward[il].SetXYZ(xyz);
686 /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
687 printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
689 AliTRDtrackerV1::FitKalman(©_track, 0x0, kFALSE, 6, points_inward);
690 memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
692 AliTRDtrackerV1::FitKalman(©_track, 0x0, kTRUE, 6, points_inward);
693 memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
695 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
696 y = points[il].GetY();
697 z = points[il].GetZ();
698 if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
699 pp = fGeo->GetPadPlane(il, stack);
700 ymin = pp->GetCol0() + epsilon;
701 ymax = pp->GetColEnd() - epsilon;
702 zmin = pp->GetRowEnd() + epsilon;
703 zmax = pp->GetRow0() - epsilon;
704 // ignore y-crossing (material)
705 if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
707 Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
708 Int_t hasTracklet = tracklet ? 1 : 0;
709 (*fDebugStream) << "FindableTracklets"
711 << "ytracklet=" << pos_tracklet[0]
713 << "ztracklet=" << pos_tracklet[1]
715 << "tracklet=" << hasTracklet
720 h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
721 if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
726 //_______________________________________________________
727 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
729 // Plot the chi2 of the track
731 if(track) fTrack = track;
733 AliWarning("No Track defined.");
737 if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
738 AliWarning("No Histogram defined.");
741 Int_t n = fTrack->GetNumberOfTracklets();
744 h->Fill(n, fTrack->GetChi2()/n);
749 //_______________________________________________________
750 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
752 // Plot the average pulse height
754 if(track) fTrack = track;
756 AliWarning("No Track defined.");
760 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
761 AliWarning("No Histogram defined.");
764 AliTRDseedV1 *tracklet = 0x0;
765 AliTRDcluster *c = 0x0;
766 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
767 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
768 Int_t crossing = Int_t(tracklet->IsRowCross());
769 Int_t detector = tracklet->GetDetector();
770 tracklet->ResetClusterIter();
771 while((c = tracklet->NextCluster())){
772 if(!c->IsInChamber()) continue;
773 Int_t localtime = c->GetLocalTimeBin();
774 Double_t absolute_charge = TMath::Abs(c->GetQ());
775 h->Fill(localtime, absolute_charge);
777 Double_t distance[2];
778 GetDistanceToTracklet(distance, tracklet, c);
779 Float_t theta = TMath::ATan(tracklet->GetZref(1));
780 Float_t phi = TMath::ATan(tracklet->GetYref(1));
781 Float_t momentum = 0.;
783 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
784 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
786 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
789 (*fDebugStream) << "PHt"
790 << "Detector=" << detector
791 << "crossing=" << crossing
792 << "Timebin=" << localtime
793 << "Charge=" << absolute_charge
794 << "momentum=" << momentum
798 << "kinkIndex=" << kinkIndex
799 << "TPCncls=" << TPCncls
800 << "dy=" << distance[0]
801 << "dz=" << distance[1]
810 //_______________________________________________________
811 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
813 // Plots the average pulse height vs the distance from the anode wire
814 // (plus const anode wire offset)
816 if(track) fTrack = track;
818 AliWarning("No Track defined.");
822 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
823 AliWarning("No Histogram defined.");
826 Float_t offset = .5*AliTRDgeometry::CamHght();
827 AliTRDseedV1 *tracklet = 0x0;
828 AliTRDcluster *c = 0x0;
829 Double_t distance = 0;
831 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
832 if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
833 tracklet->ResetClusterIter();
834 while((c = tracklet->NextCluster())){
835 if(!c->IsInChamber()) continue;
836 x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
837 y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
839 distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
840 h->Fill(distance, TMath::Abs(c->GetQ()));
846 //_______________________________________________________
847 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
849 // Plot the cluster charge
851 if(track) fTrack = track;
853 AliWarning("No Track defined.");
857 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
858 AliWarning("No Histogram defined.");
861 AliTRDseedV1 *tracklet = 0x0;
862 AliTRDcluster *c = 0x0;
863 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
864 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
865 for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
866 if(!(c = tracklet->GetClusters(itime))) continue;
873 //_______________________________________________________
874 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
876 // Plot the charge deposit per chamber
878 if(track) fTrack = track;
880 AliWarning("No Track defined.");
884 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
885 AliWarning("No Histogram defined.");
888 AliTRDseedV1 *tracklet = 0x0;
889 AliTRDcluster *c = 0x0;
891 Int_t nTracklets =fTrack->GetNumberOfTracklets();
892 for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
893 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
895 for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
896 if(!(c = tracklet->GetClusters(ic))) continue;
897 Qtot += TMath::Abs(c->GetQ());
901 Int_t crossing = (Int_t)tracklet->IsRowCross();
902 Int_t detector = tracklet->GetDetector();
903 Float_t theta = TMath::ATan(tracklet->GetZfit(1));
904 Float_t phi = TMath::ATan(tracklet->GetYfit(1));
905 Float_t momentum = 0.;
907 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
908 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
910 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
913 (*fDebugStream) << "ChargeTracklet"
914 << "Detector=" << detector
915 << "crossing=" << crossing
916 << "momentum=" << momentum
917 << "nTracklets="<< nTracklets
921 << "kinkIndex=" << kinkIndex
922 << "TPCncls=" << TPCncls
930 //_______________________________________________________
931 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
933 // Plot the number of tracks per Sector
935 if(track) fTrack = track;
937 AliWarning("No Track defined.");
941 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
942 AliWarning("No Histogram defined.");
946 // TODO we should compare with
947 // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
949 AliTRDseedV1 *tracklet = 0x0;
951 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
952 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
953 sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
961 //________________________________________________________
962 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
965 fReconstructor->SetRecoParam(r);
968 //________________________________________________________
969 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c)
971 Float_t x = c->GetX();
972 dist[0] = c->GetY() - tracklet->GetYat(x);
973 dist[1] = c->GetZ() - tracklet->GetZat(x);
977 //_______________________________________________________
978 void AliTRDcheckDET::MakePlotChi2()
980 // Plot chi2/track normalized to number of degree of freedom
981 // (tracklets) and compare with the theoretical distribution.
983 // Alex Bercuci <A.Bercuci@gsi.de>
985 TH2S *h2 = (TH2S*)fContainer->At(kChi2);
986 TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
987 TLegend *leg = new TLegend(.7,.7,.95,.95);
988 leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
989 Bool_t kFIRST = kTRUE;
990 for(Int_t il=1; il<=h2->GetNbinsX(); il++){
991 TH1D *h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
992 if(h1->Integral()<50) continue;
993 h1->Scale(1./h1->Integral());
994 h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
995 h1->SetLineColor(il);h1->SetLineStyle(2);
996 f.SetParameter(1, .5*il);f.SetLineColor(il);
997 h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
998 leg->AddEntry(h1, Form("%d", il), "l");
1000 h1->GetXaxis()->SetRangeUser(0., 25.);
1009 //________________________________________________________
1010 void AliTRDcheckDET::MakePlotNTracklets(){
1012 // Make nice bar plot of the number of tracklets in each method
1014 TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR");
1015 TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
1016 TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
1017 TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
1019 Float_t scale = hCON->Integral();
1020 hCON->Scale(100./scale);
1021 hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
1022 hCON->SetBarWidth(0.2);
1023 hCON->SetBarOffset(0.6);
1024 hCON->SetStats(kFALSE);
1025 hCON->GetYaxis()->SetRangeUser(0.,40.);
1026 hCON->GetYaxis()->SetTitleOffset(1.2);
1027 hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
1029 hBAR->Scale(100./scale);
1030 hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
1031 hBAR->SetBarWidth(0.2);
1032 hBAR->SetBarOffset(0.2);
1034 hBAR->SetStats(kFALSE);
1035 hBAR->GetYaxis()->SetRangeUser(0.,40.);
1036 hBAR->GetYaxis()->SetTitleOffset(1.2);
1037 hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
1039 hSTA->Scale(100./scale);
1040 hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
1041 hSTA->SetBarWidth(0.2);
1042 hSTA->SetBarOffset(0.4);
1044 hSTA->SetStats(kFALSE);
1045 hSTA->GetYaxis()->SetRangeUser(0.,40.);
1046 hSTA->GetYaxis()->SetTitleOffset(1.2);
1047 hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
1052 //________________________________________________________
1053 void AliTRDcheckDET::MakePlotPulseHeight(){
1055 // Create Plot of the Pluse Height Spectrum
1058 TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
1059 h = (TH1F*)arr->At(0);
1060 h->SetMarkerStyle(24);
1061 h->SetMarkerColor(kBlack);
1062 h->SetLineColor(kBlack);
1064 // copy the second histogram in a new one with the same x-dimension as the phs with respect to time
1065 h1 = (TH1F *)arr->At(1);
1066 h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
1067 for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
1068 h2->SetBinContent(ibin, h1->GetBinContent(ibin));
1069 h2->SetMarkerStyle(22);
1070 h2->SetMarkerColor(kBlue);
1071 h2->SetLineColor(kBlue);
1074 // create axis according to the histogram dimensions of the original second histogram
1075 TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1079 -0.08, 4.88, 510,"-L");
1080 axis->SetLineColor(kBlue);
1081 axis->SetLabelColor(kBlue);
1082 axis->SetTextColor(kBlue);
1083 axis->SetTitle("x_{0}-x_{c} [cm]");
1087 //________________________________________________________
1088 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1090 // Draw nice bar plots
1092 if(!histo->GetEntries()) return kFALSE;
1093 histo->Scale(100./histo->Integral());
1094 histo->SetFillColor(color);
1095 histo->SetBarOffset(.2);
1096 histo->SetBarWidth(.6);
1097 histo->Draw("bar1");