11 #include <TObjArray.h>
13 #include <TObjString.h>
17 #include <TProfile2D.h>
21 #include "AliTRDcluster.h"
22 #include "AliESDHeader.h"
23 #include "AliESDRun.h"
24 #include "AliESDtrack.h"
25 #include "AliTRDgeometry.h"
26 #include "AliTRDpadPlane.h"
27 #include "AliTRDSimParam.h"
28 #include "AliTRDseedV1.h"
29 #include "AliTRDtrackV1.h"
30 #include "AliTRDtrackerV1.h"
31 #include "AliTRDReconstructor.h"
32 #include "AliTrackReference.h"
33 #include "AliTrackPointArray.h"
34 #include "AliTracker.h"
35 #include "TTreeStream.h"
37 #include "info/AliTRDtrackInfo.h"
38 #include "info/AliTRDeventInfo.h"
39 #include "AliTRDcheckDET.h"
44 ////////////////////////////////////////////////////////////////////////////
46 // Reconstruction QA //
48 // Task doing basic checks for tracking and detector performance //
51 // Anton Andronic <A.Andronic@gsi.de> //
52 // Alexandru Bercuci <A.Bercuci@gsi.de> //
53 // Markus Fasel <M.Fasel@gsi.de> //
55 ////////////////////////////////////////////////////////////////////////////
57 //_______________________________________________________
58 AliTRDcheckDET::AliTRDcheckDET():
59 AliTRDrecoTask("checkDET", "Basic TRD data checker")
66 // Default constructor
68 DefineInput(1,AliTRDeventInfo::Class());
69 fReconstructor = new AliTRDReconstructor;
70 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
71 fGeo = new AliTRDgeometry;
75 //_______________________________________________________
76 AliTRDcheckDET::~AliTRDcheckDET(){
80 if(fTriggerNames) delete fTriggerNames;
81 delete fReconstructor;
85 //_______________________________________________________
86 void AliTRDcheckDET::ConnectInputData(Option_t *opt){
88 // Connect the Input data with the task
90 AliTRDrecoTask::ConnectInputData(opt);
91 fEventInfo = dynamic_cast<AliTRDeventInfo *>(GetInputData(1));
94 //_______________________________________________________
95 void AliTRDcheckDET::CreateOutputObjects(){
97 // Create Output Objects
99 OpenFile(0,"RECREATE");
100 fContainer = Histos();
101 if(!fTriggerNames) fTriggerNames = new TMap();
104 //_______________________________________________________
105 void AliTRDcheckDET::Exec(Option_t *opt){
107 // Execution function
108 // Filling TRD quality histos
110 if(!HasMCdata() && fEventInfo->GetEventHeader()->GetEventType() != 7) return; // For real data we select only physical events
111 AliTRDrecoTask::Exec(opt);
112 Int_t nTracks = 0; // Count the number of tracks per event
113 Int_t triggermask = fEventInfo->GetEventHeader()->GetTriggerMask();
114 TString triggername = fEventInfo->GetRunInfo()->GetFiredTriggerClasses(triggermask);
115 if(fDebugLevel > 6)printf("Trigger cluster: %d, Trigger class: %s\n", triggermask, triggername.Data());
116 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger))->Fill(triggermask);
117 for(Int_t iti = 0; iti < fTracks->GetEntriesFast(); iti++){
118 if(!fTracks->UncheckedAt(iti)) continue;
119 AliTRDtrackInfo *fTrackInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iti));
120 if(!fTrackInfo->GetTrack()) continue;
124 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks))->Fill(triggermask);
125 dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNtracksEvent))->Fill(nTracks);
127 if(triggermask <= 20 && !fTriggerNames->FindObject(Form("%d", triggermask))){
128 fTriggerNames->Add(new TObjString(Form("%d", triggermask)), new TObjString(triggername));
129 // also set the label for both histograms
130 TH1 *histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTriggerTracks));
131 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
132 histo = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kNeventsTrigger));
133 histo->GetXaxis()->SetBinLabel(histo->FindBin(triggermask), triggername);
135 PostData(0, fContainer);
138 //_______________________________________________________
139 void AliTRDcheckDET::Terminate(Option_t *){
141 // Terminate function
145 //_______________________________________________________
146 Bool_t AliTRDcheckDET::PostProcess(){
148 // Do Postprocessing (for the moment set the number of Reference histograms)
153 // Calculate of the trigger clusters purity
154 h = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTrigger"));
155 TH1F *h1 = dynamic_cast<TH1F *>(fContainer->FindObject("hEventsTriggerTracks"));
157 Float_t purities[20], val = 0;
158 TString triggernames[20];
159 Int_t nTriggerClasses = 0;
160 for(Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++){
161 if((val = h1->GetBinContent(ibin))){
162 purities[nTriggerClasses] = val;
163 triggernames[nTriggerClasses] = h1->GetXaxis()->GetBinLabel(ibin);
167 h = dynamic_cast<TH1F *>(fContainer->UncheckedAt(kTriggerPurity));
168 TAxis *ax = h->GetXaxis();
169 for(Int_t itrg = 0; itrg < nTriggerClasses; itrg++){
170 h->Fill(itrg, purities[itrg]);
171 ax->SetBinLabel(itrg+1, triggernames[itrg].Data());
173 ax->SetRangeUser(-0.5, nTriggerClasses+.5);
174 h->GetYaxis()->SetRangeUser(0,1);
181 //_______________________________________________________
182 Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){
184 // Setting Reference Figures
189 case kNclustersTrack:
190 ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl");
192 case kNclustersTracklet:
193 ((TH1F*)fContainer->FindObject("hNclTls"))->Draw("pc");
195 case kNtrackletsTrack:
196 MakePlotNTracklets();
198 case kNtrackletsCross:
199 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsCross"), kRed)) break;
201 case kNtrackletsFindable:
202 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtlsFindable"), kGreen)) break;
205 ((TH1F*)fContainer->FindObject("hNtrks"))->Draw("pl");
208 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hNtrksSector"), kGreen)) break;
214 MakePlotPulseHeight();
217 ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c");
220 case kChargeTracklet:
221 ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c");
223 case kNeventsTrigger:
224 ((TH1F*)fContainer->FindObject("hEventsTrigger"))->Draw("");
226 case kNeventsTriggerTracks:
227 ((TH1F*)fContainer->FindObject("hEventsTriggerTracks"))->Draw("");
230 if(!MakeBarPlot((TH1F*)fContainer->FindObject("hTriggerPurity"), kGreen)) break;
235 AliInfo(Form("Reference plot [%d] missing result", ifig));
239 //_______________________________________________________
240 TObjArray *AliTRDcheckDET::Histos(){
242 // Create QA histograms
244 if(fContainer) return fContainer;
246 fContainer = new TObjArray(20);
247 //fContainer->SetOwner(kTRUE);
249 // Register Histograms
251 if(!(h = (TH1F *)gROOT->FindObject("hNcls"))){
252 h = new TH1F("hNcls", "N_{clusters} / track", 181, -0.5, 180.5);
253 h->GetXaxis()->SetTitle("N_{clusters}");
254 h->GetYaxis()->SetTitle("Entries");
256 fContainer->AddAt(h, kNclustersTrack);
258 if(!(h = (TH1F *)gROOT->FindObject("hNclTls"))){
259 h = new TH1F("hNclTls","N_{clusters} / tracklet", 51, -0.5, 50.5);
260 h->GetXaxis()->SetTitle("N_{clusters}");
261 h->GetYaxis()->SetTitle("Entries");
263 fContainer->AddAt(h, kNclustersTracklet);
265 if(!(h = (TH1F *)gROOT->FindObject("hNtls"))){
266 h = new TH1F("hNtls", "N_{tracklets} / track", AliTRDgeometry::kNlayer, 0.5, 6.5);
267 h->GetXaxis()->SetTitle("N^{tracklet}");
268 h->GetYaxis()->SetTitle("freq. [%]");
270 fContainer->AddAt(h, kNtrackletsTrack);
272 if(!(h = (TH1F *)gROOT->FindObject("htlsSTA"))){
273 h = new TH1F("hNtlsSTA", "N_{tracklets} / track (Stand Alone)", AliTRDgeometry::kNlayer, 0.5, 6.5);
274 h->GetXaxis()->SetTitle("N^{tracklet}");
275 h->GetYaxis()->SetTitle("freq. [%]");
277 fContainer->AddAt(h, kNtrackletsSTA);
279 if(!(h = (TH1F *)gROOT->FindObject("htlsBAR"))){
280 h = new TH1F("hNtlsBAR", "N_{tracklets} / track (Barrel)", AliTRDgeometry::kNlayer, 0.5, 6.5);
281 h->GetXaxis()->SetTitle("N^{tracklet}");
282 h->GetYaxis()->SetTitle("freq. [%]");
284 fContainer->AddAt(h, kNtrackletsBAR);
287 if(!(h = (TH1F *)gROOT->FindObject("hNtlsCross"))){
288 h = new TH1F("hNtlsCross", "N_{tracklets}^{cross} / track", 7, -0.5, 6.5);
289 h->GetXaxis()->SetTitle("n_{row cross}");
290 h->GetYaxis()->SetTitle("freq. [%]");
292 fContainer->AddAt(h, kNtrackletsCross);
294 if(!(h = (TH1F *)gROOT->FindObject("hNtlsFindable"))){
295 h = new TH1F("hNtlsFindable", "Found/Findable Tracklets" , 101, -0.005, 1.005);
296 h->GetXaxis()->SetTitle("r [a.u]");
297 h->GetYaxis()->SetTitle("Entries");
299 fContainer->AddAt(h, kNtrackletsFindable);
301 if(!(h = (TH1F *)gROOT->FindObject("hNtrks"))){
302 h = new TH1F("hNtrks", "N_{tracks} / event", 100, 0, 100);
303 h->GetXaxis()->SetTitle("N_{tracks}");
304 h->GetYaxis()->SetTitle("Entries");
306 fContainer->AddAt(h, kNtracksEvent);
308 if(!(h = (TH1F *)gROOT->FindObject("hNtrksSector"))){
309 h = new TH1F("hNtrksSector", "N_{tracks} / sector", AliTRDgeometry::kNsector, -0.5, 17.5);
310 h->GetXaxis()->SetTitle("sector");
311 h->GetYaxis()->SetTitle("freq. [%]");
313 fContainer->AddAt(h, kNtracksSector);
316 TObjArray *arr = new TObjArray(2);
317 arr->SetOwner(kTRUE); arr->SetName("<PH>");
318 fContainer->AddAt(arr, kPH);
319 if(!(h = (TH1F *)gROOT->FindObject("hPHt"))){
320 h = new TProfile("hPHt", "<PH>", 31, -0.5, 30.5);
321 h->GetXaxis()->SetTitle("Time / 100ns");
322 h->GetYaxis()->SetTitle("<PH> [a.u]");
325 if(!(h = (TH1F *)gROOT->FindObject("hPHx")))
326 h = new TProfile("hPHx", "<PH>", 31, -0.08, 4.88);
331 if(!(h = (TH2S*)gROOT->FindObject("hChi2"))){
332 h = new TH2S("hChi2", "#chi^{2} per track", AliTRDgeometry::kNlayer, .5, AliTRDgeometry::kNlayer+.5, 100, 0, 50);
334 h->SetYTitle("#chi^{2}/ndf");
335 h->SetZTitle("entries");
337 fContainer->AddAt(h, kChi2);
339 if(!(h = (TH1F *)gROOT->FindObject("hQcl"))){
340 h = new TH1F("hQcl", "Q_{cluster}", 200, 0, 1200);
341 h->GetXaxis()->SetTitle("Q_{cluster} [a.u.]");
342 h->GetYaxis()->SetTitle("Entries");
344 fContainer->AddAt(h, kChargeCluster);
346 if(!(h = (TH1F *)gROOT->FindObject("hQtrklt"))){
347 h = new TH1F("hQtrklt", "Q_{tracklet}", 6000, 0, 6000);
348 h->GetXaxis()->SetTitle("Q_{tracklet} [a.u.]");
349 h->GetYaxis()->SetTitle("Entries");
351 fContainer->AddAt(h, kChargeTracklet);
354 if(!(h = (TH1F *)gROOT->FindObject("hEventsTrigger")))
355 h = new TH1F("hEventsTrigger", "Trigger Class", 100, 0, 100);
357 fContainer->AddAt(h, kNeventsTrigger);
359 if(!(h = (TH1F *)gROOT->FindObject("hEventsTriggerTracks")))
360 h = new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100);
362 fContainer->AddAt(h, kNeventsTriggerTracks);
364 if(!(h = (TH1F *)gROOT->FindObject("hTriggerPurity"))){
365 h = new TH1F("hTriggerPurity", "Trigger Purity", 10, -0.5, 9.5);
366 h->GetXaxis()->SetTitle("Trigger Cluster");
367 h->GetYaxis()->SetTitle("freq.");
369 fContainer->AddAt(h, kTriggerPurity);
378 //_______________________________________________________
379 TH1 *AliTRDcheckDET::PlotNClustersTracklet(const AliTRDtrackV1 *track){
381 // Plot the mean number of clusters per tracklet
383 if(track) fTrack = track;
385 AliWarning("No Track defined.");
389 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTracklet)))){
390 AliWarning("No Histogram defined.");
393 AliTRDseedV1 *tracklet = 0x0;
394 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
395 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
396 h->Fill(tracklet->GetN2());
401 //_______________________________________________________
402 TH1 *AliTRDcheckDET::PlotNClustersTrack(const AliTRDtrackV1 *track){
404 // Plot the number of clusters in one track
406 if(track) fTrack = track;
408 AliWarning("No Track defined.");
412 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNclustersTrack)))){
413 AliWarning("No Histogram defined.");
418 AliTRDseedV1 *tracklet = 0x0;
419 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
420 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
421 nclusters += tracklet->GetN();
423 Int_t crossing = Int_t(tracklet->IsRowCross());
424 Int_t detector = tracklet->GetDetector();
425 Float_t theta = TMath::ATan(tracklet->GetZref(1));
426 Float_t phi = TMath::ATan(tracklet->GetYref(1));
427 Float_t momentum = 0.;
429 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
430 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
432 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
435 (*fDebugStream) << "NClustersTrack"
436 << "Detector=" << detector
437 << "crossing=" << crossing
438 << "momentum=" << momentum
442 << "kinkIndex=" << kinkIndex
443 << "TPCncls=" << TPCncls
444 << "nclusters=" << nclusters
453 //_______________________________________________________
454 TH1 *AliTRDcheckDET::PlotNTrackletsTrack(const AliTRDtrackV1 *track){
456 // Plot the number of tracklets
458 if(track) fTrack = track;
460 AliWarning("No Track defined.");
463 TH1 *h = 0x0, *hMethod = 0x0;
464 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsTrack)))){
465 AliWarning("No Histogram defined.");
468 Int_t status = fESD->GetStatus();
469 /* 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);*/
470 if((status & AliESDtrack::kTRDin) != 0){
472 if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsBAR))))
473 AliWarning("Method: Barrel. Histogram not processed!");
476 if(!(hMethod = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsSTA))))
477 AliWarning("Method: StandAlone. Histogram not processed!");
479 Int_t nTracklets = fTrack->GetNumberOfTracklets();
481 hMethod->Fill(nTracklets);
484 // If we have one Tracklet, check in which layer this happens
486 AliTRDseedV1 *tracklet = 0x0;
487 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
488 if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){layer = il; break;}
490 (*fDebugStream) << "NTrackletsTrack"
499 //_______________________________________________________
500 TH1 *AliTRDcheckDET::PlotNTrackletsRowCross(const AliTRDtrackV1 *track){
502 // Plot the number of tracklets
504 if(track) fTrack = track;
506 AliWarning("No Track defined.");
510 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsCross)))){
511 AliWarning("No Histogram defined.");
516 AliTRDseedV1 *tracklet = 0x0;
517 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
518 if(!(tracklet = fTrack->GetTracklet(il)) || !tracklet->IsOK()) continue;
520 if(tracklet->IsRowCross()) ncross++;
526 //_______________________________________________________
527 TH1 *AliTRDcheckDET::PlotFindableTracklets(const AliTRDtrackV1 *track){
529 // Plots the ratio of number of tracklets vs.
530 // number of findable tracklets
532 // Findable tracklets are defined as track prolongation
533 // to layer i does not hit the dead area +- epsilon
535 // In order to check whether tracklet hist active area in Layer i,
536 // the track is refitted and the fitted position + an uncertainty
537 // range is compared to the chamber border (also with a different
540 // For the track fit two cases are distinguished:
541 // If the track is a stand alone track (defined by the status bit
542 // encoding, then the track is fitted with the tilted Rieman model
543 // Otherwise the track is fitted with the Kalman fitter in two steps:
544 // Since the track parameters are give at the outer point, we first
545 // fit in direction inwards. Afterwards we fit again in direction outwards
546 // to extrapolate the track to layers which are not reached by the track
547 // For the Kalman model, the radial track points have to be shifted by
548 // a distance epsilon in the direction that we want to fit
550 const Float_t epsilon = 0.01; // dead area tolerance
551 const Float_t epsilon_R = 1; // shift in radial direction of the anode wire position (Kalman filter only)
552 const Float_t delta_y = 0.7; // Tolerance in the track position in y-direction
553 const Float_t delta_z = 7.0; // Tolerance in the track position in z-direction (Padlength)
554 Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
556 if(track) fTrack = track;
558 AliWarning("No Track defined.");
562 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtrackletsFindable)))){
563 AliWarning("No Histogram defined.");
566 Int_t nFound = 0, nFindable = 0;
568 Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
569 Double_t y = 0., z = 0.;
570 AliTRDseedV1 *tracklet = 0x0;
572 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
573 if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
574 tracklet->SetReconstructor(fReconstructor);
578 // 2 Different cases:
579 // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
580 // 2nd barrel track: here we propagate the track to the layers
581 AliTrackPoint points[6];
583 memset(xyz, 0, sizeof(Float_t) * 3);
584 if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
586 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
587 xyz[0] = x_anode[il];
588 points[il].SetXYZ(xyz);
590 AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
596 // -> Kalman outwards
597 AliTRDtrackV1 copy_track(*fTrack); // Do Kalman on a (non-constant) copy of the track
598 AliTrackPoint points_inward[6], points_outward[6];
599 for(Int_t il = AliTRDgeometry::kNlayer; il--;){
600 // In order to avoid complications in the Kalman filter if the track points have the same radial
601 // position like the tracklets, we have to shift the radial postion of the anode wire by epsilon
602 // in the direction we want to go
603 // The track points have to be in reverse order for the Kalman Filter inwards
604 xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
605 points_inward[il].SetXYZ(xyz);
606 xyz[0] = x_anode[il] + epsilon_R;
607 points_outward[il].SetXYZ(xyz);
609 /*for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
610 printf("%d. X = %f\n", ipt, points[ipt].GetX());*/
612 AliTRDtrackerV1::FitKalman(©_track, 0x0, kFALSE, 6, points_inward);
613 memcpy(points, points_inward, sizeof(AliTrackPoint) * 6); // Preliminary store the inward results in the Array points
615 AliTRDtrackerV1::FitKalman(©_track, 0x0, kTRUE, 6, points_inward);
616 memcpy(points, points_outward, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
618 for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
619 y = points[il].GetY();
620 z = points[il].GetZ();
621 if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
622 pp = fGeo->GetPadPlane(il, stack);
623 ymin = pp->GetCol0() + epsilon;
624 ymax = pp->GetColEnd() - epsilon;
625 zmin = pp->GetRowEnd() + epsilon;
626 zmax = pp->GetRow0() - epsilon;
627 // ignore y-crossing (material)
628 if((z + delta_z > zmin && z - delta_z < zmax) && (y + delta_y > ymin && y - delta_y < ymax)) nFindable++;
630 Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetZfit(0) : 0};
631 Int_t hasTracklet = tracklet ? 1 : 0;
632 (*fDebugStream) << "FindableTracklets"
634 << "ytracklet=" << pos_tracklet[0]
636 << "ztracklet=" << pos_tracklet[1]
638 << "tracklet=" << hasTracklet
643 h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
644 if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
649 //_______________________________________________________
650 TH1 *AliTRDcheckDET::PlotChi2(const AliTRDtrackV1 *track){
652 // Plot the chi2 of the track
654 if(track) fTrack = track;
656 AliWarning("No Track defined.");
660 if(!(h = dynamic_cast<TH2S*>(fContainer->At(kChi2)))) {
661 AliWarning("No Histogram defined.");
664 Int_t n = fTrack->GetNumberOfTracklets();
667 h->Fill(n, fTrack->GetChi2()/n);
672 //_______________________________________________________
673 TH1 *AliTRDcheckDET::PlotPHt(const AliTRDtrackV1 *track){
675 // Plot the average pulse height
677 if(track) fTrack = track;
679 AliWarning("No Track defined.");
683 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(0)))){
684 AliWarning("No Histogram defined.");
687 AliTRDseedV1 *tracklet = 0x0;
688 AliTRDcluster *c = 0x0;
689 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
690 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
691 Int_t crossing = Int_t(tracklet->IsRowCross());
692 Int_t detector = tracklet->GetDetector();
693 tracklet->ResetClusterIter();
694 while((c = tracklet->NextCluster())){
695 if(!c->IsInChamber()) continue;
696 Int_t localtime = c->GetLocalTimeBin();
697 Double_t absolute_charge = TMath::Abs(c->GetQ());
698 h->Fill(localtime, absolute_charge);
700 Double_t distance[2];
701 GetDistanceToTracklet(distance, tracklet, c);
702 Float_t theta = TMath::ATan(tracklet->GetZref(1));
703 Float_t phi = TMath::ATan(tracklet->GetYref(1));
704 Float_t momentum = 0.;
706 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
707 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
709 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
712 (*fDebugStream) << "PHt"
713 << "Detector=" << detector
714 << "crossing=" << crossing
715 << "Timebin=" << localtime
716 << "Charge=" << absolute_charge
717 << "momentum=" << momentum
721 << "kinkIndex=" << kinkIndex
722 << "TPCncls=" << TPCncls
723 << "dy=" << distance[0]
724 << "dz=" << distance[1]
733 //_______________________________________________________
734 TH1 *AliTRDcheckDET::PlotPHx(const AliTRDtrackV1 *track){
736 // Plots the average pulse height vs the distance from the anode wire
737 // (plus const anode wire offset)
739 if(track) fTrack = track;
741 AliWarning("No Track defined.");
745 if(!(h = dynamic_cast<TProfile *>(((TObjArray*)(fContainer->At(kPH)))->At(1)))){
746 AliWarning("No Histogram defined.");
749 Float_t offset = .5*AliTRDgeometry::CamHght();
750 AliTRDseedV1 *tracklet = 0x0;
751 AliTRDcluster *c = 0x0;
752 Double_t distance = 0;
754 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
755 if(!(tracklet = fTrack->GetTracklet(itl)) || !(tracklet->IsOK())) continue;
756 tracklet->ResetClusterIter();
757 while((c = tracklet->NextCluster())){
758 if(!c->IsInChamber()) continue;
759 x = c->GetX()-AliTRDcluster::GetXcorr(c->GetLocalTimeBin());
760 y = c->GetY()-AliTRDcluster::GetYcorr(AliTRDgeometry::GetLayer(c->GetDetector()), c->GetCenter());
762 distance = tracklet->GetX0() - (c->GetX() + 0.3) + offset;
763 h->Fill(distance, TMath::Abs(c->GetQ()));
769 //_______________________________________________________
770 TH1 *AliTRDcheckDET::PlotChargeCluster(const AliTRDtrackV1 *track){
772 // Plot the cluster charge
774 if(track) fTrack = track;
776 AliWarning("No Track defined.");
780 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeCluster)))){
781 AliWarning("No Histogram defined.");
784 AliTRDseedV1 *tracklet = 0x0;
785 AliTRDcluster *c = 0x0;
786 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
787 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
788 for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
789 if(!(c = tracklet->GetClusters(itime))) continue;
796 //_______________________________________________________
797 TH1 *AliTRDcheckDET::PlotChargeTracklet(const AliTRDtrackV1 *track){
799 // Plot the charge deposit per chamber
801 if(track) fTrack = track;
803 AliWarning("No Track defined.");
807 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kChargeTracklet)))){
808 AliWarning("No Histogram defined.");
811 AliTRDseedV1 *tracklet = 0x0;
812 AliTRDcluster *c = 0x0;
814 Int_t nTracklets =fTrack->GetNumberOfTracklets();
815 for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
816 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
818 for(Int_t ic = AliTRDseedV1::kNclusters; ic--;){
819 if(!(c = tracklet->GetClusters(ic))) continue;
820 Qtot += TMath::Abs(c->GetQ());
824 Int_t crossing = (Int_t)tracklet->IsRowCross();
825 Int_t detector = tracklet->GetDetector();
826 Float_t theta = TMath::ATan(tracklet->GetZfit(1));
827 Float_t phi = TMath::ATan(tracklet->GetYfit(1));
828 Float_t momentum = 0.;
830 Int_t kinkIndex = fESD ? fESD->GetKinkIndex() : 0;
831 UShort_t TPCncls = fESD ? fESD->GetTPCncls() : 0;
833 if(fMC->GetTrackRef()) momentum = fMC->GetTrackRef()->P();
836 (*fDebugStream) << "ChargeTracklet"
837 << "Detector=" << detector
838 << "crossing=" << crossing
839 << "momentum=" << momentum
840 << "nTracklets="<< nTracklets
844 << "kinkIndex=" << kinkIndex
845 << "TPCncls=" << TPCncls
853 //_______________________________________________________
854 TH1 *AliTRDcheckDET::PlotNTracksSector(const AliTRDtrackV1 *track){
856 // Plot the number of tracks per Sector
858 if(track) fTrack = track;
860 AliWarning("No Track defined.");
864 if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNtracksSector)))){
865 AliWarning("No Histogram defined.");
869 // TODO we should compare with
870 // sector = Int_t(track->GetAlpha() / AliTRDgeometry::GetAlpha());
872 AliTRDseedV1 *tracklet = 0x0;
874 for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
875 if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
876 sector = static_cast<Int_t>(tracklet->GetDetector()/AliTRDgeometry::kNdets);
884 //________________________________________________________
885 void AliTRDcheckDET::SetRecoParam(AliTRDrecoParam *r)
888 fReconstructor->SetRecoParam(r);
891 //________________________________________________________
892 void AliTRDcheckDET::GetDistanceToTracklet(Double_t *dist, AliTRDseedV1 *tracklet, AliTRDcluster *c)
894 Float_t x = c->GetX();
895 dist[0] = c->GetY() - tracklet->GetYat(x);
896 dist[1] = c->GetZ() - tracklet->GetZat(x);
900 //_______________________________________________________
901 void AliTRDcheckDET::MakePlotChi2()
903 // Plot chi2/track normalized to number of degree of freedom
904 // (tracklets) and compare with the theoretical distribution.
906 // Alex Bercuci <A.Bercuci@gsi.de>
908 TH2S *h2 = (TH2S*)fContainer->At(kChi2);
909 TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.);
910 TLegend *leg = new TLegend(.7,.7,.95,.95);
911 leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track");
912 Bool_t kFIRST = kTRUE;
913 for(Int_t il=1; il<=h2->GetNbinsX(); il++){
914 TH1D *h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il);
915 if(h1->Integral()<50) continue;
916 h1->Scale(1./h1->Integral());
917 h1->SetMarkerStyle(7);h1->SetMarkerColor(il);
918 h1->SetLineColor(il);h1->SetLineStyle(2);
919 f.SetParameter(1, .5*il);f.SetLineColor(il);
920 h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame");
921 leg->AddEntry(h1, Form("%d", il), "l");
923 h1->GetXaxis()->SetRangeUser(0., 25.);
932 //________________________________________________________
933 void AliTRDcheckDET::MakePlotNTracklets(){
935 // Make nice bar plot of the number of tracklets in each method
937 TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR");
938 TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA");
939 TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls");
940 TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89);
942 Float_t scale = hCON->Integral();
943 hCON->Scale(100./scale);
944 hCON->SetFillColor(kRed);hCON->SetLineColor(kRed);
945 hCON->SetBarWidth(0.2);
946 hCON->SetBarOffset(0.6);
947 hCON->SetStats(kFALSE);
948 hCON->GetYaxis()->SetRangeUser(0.,40.);
949 hCON->GetYaxis()->SetTitleOffset(1.2);
950 hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f");
952 hBAR->Scale(100./scale);
953 hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen);
954 hBAR->SetBarWidth(0.2);
955 hBAR->SetBarOffset(0.2);
957 hBAR->SetStats(kFALSE);
958 hBAR->GetYaxis()->SetRangeUser(0.,40.);
959 hBAR->GetYaxis()->SetTitleOffset(1.2);
960 hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f");
962 hSTA->Scale(100./scale);
963 hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue);
964 hSTA->SetBarWidth(0.2);
965 hSTA->SetBarOffset(0.4);
967 hSTA->SetStats(kFALSE);
968 hSTA->GetYaxis()->SetRangeUser(0.,40.);
969 hSTA->GetYaxis()->SetTitleOffset(1.2);
970 hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f");
975 //________________________________________________________
976 void AliTRDcheckDET::MakePlotPulseHeight(){
978 // Create Plot of the Pluse Height Spectrum
981 TObjArray *arr = (TObjArray*)fContainer->FindObject("<PH>");
982 h = (TH1F*)arr->At(0);
983 h->SetMarkerStyle(24);
984 h->SetMarkerColor(kBlack);
985 h->SetLineColor(kBlack);
987 // copy the second histogram in a new one with the same x-dimension as the phs with respect to time
988 h1 = (TH1F *)arr->At(1);
989 h2 = new TH1F("hphs1","Average PH", 31, -0.5, 30.5);
990 for(Int_t ibin = h1->GetXaxis()->GetFirst(); ibin < h1->GetNbinsX(); ibin++)
991 h2->SetBinContent(ibin, h1->GetBinContent(ibin));
992 h2->SetMarkerStyle(22);
993 h2->SetMarkerColor(kBlue);
994 h2->SetLineColor(kBlue);
997 // create axis according to the histogram dimensions of the original second histogram
998 TGaxis *axis = new TGaxis(gPad->GetUxmin(),
1002 -0.08, 4.88, 510,"-L");
1003 axis->SetLineColor(kBlue);
1004 axis->SetLabelColor(kBlue);
1005 axis->SetTextColor(kBlue);
1006 axis->SetTitle("x_{0}-x_{c} [cm]");
1010 //________________________________________________________
1011 Bool_t AliTRDcheckDET::MakeBarPlot(TH1 *histo, Int_t color){
1013 // Draw nice bar plots
1015 if(!histo->GetEntries()) return kFALSE;
1016 histo->Scale(100./histo->Integral());
1017 histo->SetFillColor(color);
1018 histo->SetBarOffset(.2);
1019 histo->SetBarWidth(.6);
1020 histo->Draw("bar1");