Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /* $Id: AliTRDqaRecPoints.cxx 23387 2008-01-17 17:25:16Z cblume $ */
17 :
18 : ////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // Produces the data needed to calculate the quality assurance. //
21 : // All data must be mergeable objects. //
22 : // //
23 : // Author: //
24 : // Sylwester Radomski (radomski@physi.uni-heidelberg.de) //
25 : // //
26 : ////////////////////////////////////////////////////////////////////////////
27 :
28 : // --- ROOT system ---
29 : #include <TClonesArray.h>
30 : #include <TFile.h>
31 : #include <TH1D.h>
32 : #include <TH2D.h>
33 : #include <TH3D.h>
34 : #include <TProfile.h>
35 : #include <TF1.h>
36 : #include <TCanvas.h>
37 :
38 : // --- AliRoot header files ---
39 : #include "AliESDEvent.h"
40 : #include "AliLog.h"
41 : #include "AliTRDcluster.h"
42 : #include "AliTRDqaRecPoints.h"
43 : #include "AliTRDgeometry.h"
44 :
45 : #include "AliQAChecker.h"
46 :
47 16 : ClassImp(AliTRDqaRecPoints)
48 :
49 : //____________________________________________________________________________
50 :
51 : AliTRDqaRecPoints::AliTRDqaRecPoints() :
52 0 : TObject(),
53 0 : fnEvents(0),
54 0 : fHist(0),
55 0 : fnPad(0),
56 0 : fRef(0)
57 0 : {
58 : //
59 : // Default constructor
60 : //
61 :
62 0 : for (Int_t i = 0; i < 540; i++) {
63 0 : fRefHist[i] = 0x0;
64 : }
65 :
66 0 : }
67 :
68 : //____________________________________________________________________________
69 :
70 : AliTRDqaRecPoints::AliTRDqaRecPoints(const AliTRDqaRecPoints &/*qa*/) :
71 0 : TObject(),
72 0 : fnEvents(0),
73 0 : fHist(0),
74 0 : fnPad(0),
75 0 : fRef(0)
76 0 : {
77 : //
78 : // Copy constructor
79 : //
80 :
81 0 : for (Int_t i = 0; i < 540; i++) {
82 0 : fRefHist[i] = 0x0;
83 : }
84 :
85 0 : }
86 :
87 : //____________________________________________________________________________
88 : void AliTRDqaRecPoints::Process(const char* filename)
89 : {
90 : //
91 : // Detector specific actions at end of cycle
92 : //
93 : //TStopwatch watch;
94 : //watch.Start();
95 :
96 0 : AliInfo("End of TRD cycle");
97 :
98 : //if (task == AliQAv1::kRECPOINTS) {
99 :
100 0 : TH1D *hist = new TH1D("fitHist", "", 200, -0.5, 199.5);
101 : //fHist->Print();
102 :
103 : // fill detector map;
104 0 : for(int i=0; i<540; i++) {
105 0 : Double_t v = ((TH1D*)fHist->At(0))->GetBinContent(i+1);
106 0 : Int_t sm = i/30;
107 0 : Int_t det = i%30;
108 :
109 0 : TH2D *detMap = (TH2D*)fHist->At(87);
110 0 : Int_t bin = detMap->FindBin(sm, det);
111 0 : detMap->SetBinContent(bin, v);
112 : }
113 :
114 :
115 : // Rec points full chambers
116 0 : for (Int_t i=0; i<540; i++) {
117 :
118 : //AliInfo(Form("I = %d", i));
119 :
120 : //TH1D *h = ((TH2D*)fHist->At(1))->ProjectionY(Form("qaTRD_recPoints_amp_%d",i), i+1, i+1);
121 0 : hist->Reset();
122 0 : for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
123 0 : Double_t xvalue = hist->GetBinCenter(b);
124 0 : Int_t bin = ((TH2D*)fHist->At(1))->FindBin(i,xvalue);
125 0 : Double_t value = ((TH2D*)fHist->At(1))->GetBinContent(bin);
126 0 : hist->SetBinContent(b, value);
127 : }
128 :
129 : //printf("Sum = %d %f\n", i, hist->GetSum());
130 0 : if (hist->GetSum() < 100) continue; // chamber not present
131 :
132 0 : hist->Fit("landau", "q0", "goff", 10, 180);
133 0 : TF1 *fit = hist->GetFunction("landau");
134 0 : ((TH1D*)fHist->At(12))->Fill(fit->GetParameter(1));
135 0 : ((TH1D*)fHist->At(13))->Fill(fit->GetParameter(2));
136 0 : }
137 :
138 : // time-bin by time-bin sm by sm
139 0 : for(Int_t i=0; i<18; i++) { // loop over super-modules
140 :
141 0 : for(Int_t j=0; j<35; j++) { // loop over time bins
142 :
143 : //TH1D *h = ((TH3D*)fHist->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, j+1, j+1);
144 0 : hist->Reset();
145 0 : for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
146 0 : Double_t xvalue = hist->GetBinCenter(b);
147 : Double_t svalue = 0;
148 :
149 0 : for(Int_t det=i*30; det<(i+1)*30; det++) { // loop over detectors
150 0 : Int_t bin = ((TH3D*)fHist->At(10))->FindBin(det,j,xvalue);
151 0 : Double_t value = ((TH3D*)fHist->At(10))->GetBinContent(bin);
152 0 : svalue += value;
153 : }
154 : //printf("v = %f\n", value);
155 0 : hist->SetBinContent(b, svalue);
156 : }
157 :
158 0 : if (hist->GetSum() < 100) continue;
159 : //printf("fitting %d %d %f\n", i, j, hist->GetSum());
160 :
161 0 : hist->Fit("landau", "q0", "goff", 10, 180);
162 0 : TF1 *fit = hist->GetFunction("landau");
163 :
164 0 : TH1D *h1 = (TH1D*)fHist->At(14+18+i);
165 0 : Int_t bin = h1->FindBin(j);
166 : // printf("%d %d %d\n", det, j, bin);
167 0 : h1->SetBinContent(bin, TMath::Abs(fit->GetParameter(1)));
168 0 : }
169 : }
170 :
171 :
172 : // time-bin by time-bin chamber by chamber
173 :
174 0 : for (Int_t i=0; i<540; i++) {
175 :
176 : //TH1D *test = ((TH3D*)fHist->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, 0, 35);
177 : //if (test->GetSum() < 100) continue;
178 :
179 : //AliInfo(Form("fitting det = %d", i));
180 :
181 0 : for(Int_t j=0; j<35; j++) {
182 :
183 : //TH1D *h = ((TH3D*)fHist->At(10))->ProjectionZ(Form("ampTime_%d",i), i+1, i+1, j+1, j+1);
184 0 : hist->Reset();
185 0 : for(Int_t b=1; b<hist->GetXaxis()->GetNbins()-1; b++) {
186 0 : Double_t xvalue = hist->GetBinCenter(b);
187 0 : Int_t bin = ((TH3D*)fHist->At(10))->FindBin(i,j,xvalue);
188 0 : Double_t value = ((TH3D*)fHist->At(10))->GetBinContent(bin);
189 : //printf("v = %f\n", value);
190 0 : hist->SetBinContent(b, value);
191 : }
192 :
193 0 : if (hist->GetSum() < 100) continue;
194 : //printf("fitting %d %d %f\n", i, j, hist->GetSum());
195 :
196 0 : hist->Fit("landau", "q0", "goff", 10, 180);
197 0 : TF1 *fit = hist->GetFunction("landau");
198 :
199 0 : Int_t sm = i/30;
200 0 : Int_t det = i%30;
201 0 : TH2D *h2 = (TH2D*)fHist->At(14+sm);
202 0 : Int_t bin = h2->FindBin(det,j);
203 : // printf("%d %d %d\n", det, j, bin);
204 0 : h2->SetBinContent(bin, TMath::Abs(fit->GetParameter(1)));
205 0 : h2->SetBinError(bin,fit->GetParError(1));
206 0 : }
207 : }
208 :
209 0 : if (hist) delete hist;
210 :
211 :
212 0 : TFile *outFile = new TFile(filename, "RECREATE");
213 0 : outFile->mkdir("TRD");
214 0 : gDirectory->cd("TRD");
215 0 : gDirectory->mkdir("RecPoints");
216 0 : gDirectory->cd("RecPoints");
217 0 : fHist->Write();
218 :
219 0 : if (fRef) {
220 0 : for(Int_t i=0; i<540; i++) {
221 : //fRefHist[i]->Scale(1./fnEvents);
222 0 : fRefHist[i]->Write();
223 : }
224 0 : }
225 :
226 0 : outFile->Close();
227 :
228 0 : }
229 :
230 : //____________________________________________________________________________
231 :
232 : void AliTRDqaRecPoints::Init()
233 : {
234 : //
235 : // Create Reconstructed Points histograms in RecPoints subdir
236 : //
237 :
238 : //const Int_t kNhist = 14 + 4 * 18 + 2;
239 : const Int_t kNhist = 88+2*18+1;
240 0 : TH1 *hist[kNhist];
241 :
242 0 : hist[0] = new TH1D("qaTRD_recPoints_det", ";Detector ID of the cluster", 540, -0.5, 539.5);
243 0 : hist[1] = new TH2D("qaTRD_recPoints_amp", ";Amplitude", 540, -0.5, 539, 200, -0.5, 199.5);
244 0 : hist[2] = new TH1D("qaTRD_recPoints_npad", ";Number of Pads", 12, -0.5, 11.5);
245 :
246 0 : hist[3] = new TH1D("qaTRD_recPoints_dist2", ";residuals [2pad]", 100, -1, 1);
247 0 : hist[4] = new TH1D("qaTRD_recPoints_dist3", ";residuals [3pad]", 100, -1, 1);
248 0 : hist[5] = new TH1D("qaTRD_recPoints_dist4", ";residuals [4pad]", 100, -1, 1);
249 0 : hist[6] = new TH1D("qaTRD_recPoints_dist5", ";residuals [5pad]", 100, -1, 1);
250 :
251 0 : hist[7] = new TH2D("qaTRD_recPoints_rowCol", ";row;col", 16, -0.5, 15.5, 145, -0.5, 144.5);
252 0 : hist[8] = new TH1D("qaTRD_recPoints_time", ";time bin", 35, -0.5, 34.5);
253 0 : hist[9] = new TH1D("qaTRD_recPoints_nCls", ";number of clusters", 500, -0.5, 499.5);
254 :
255 0 : hist[10] = new TH3D("qaTRD_recPoints_sigTime", ";chamber;time bin;signal",
256 : 540, -0.5, 539.5, 35, -0.5, 34.5, 200, -0.5, 199.5);
257 0 : hist[11] = new TProfile("qaTRD_recPoints_prf", ";distance;center of gravity"
258 : , 120, -0.6, 0.6, -1.2, 1.2, "");
259 :
260 0 : hist[12] = new TH1D("qaTRD_recPoints_ampMPV", ";amplitude MPV", 150, 0, 150);
261 0 : hist[13] = new TH1D("qaTRD_recPoints_ampSigma", ";amplitude Sigma", 200, 0, 200);
262 :
263 : // chamber by chamber
264 0 : for(Int_t i=0; i<18; i++) {
265 0 : hist[14+i] = new TH2D(Form("qaTRD_recPoints_sigTime_sm%d",i), Form("sm%d;det;time bin",i),
266 : 30, -0.5, 29.5, 35, -0.5, 34.5);
267 0 : hist[14+i]->SetMinimum(20);
268 0 : hist[14+i]->SetMaximum(40);
269 : }
270 :
271 : // time bin by time bin sm-by-sm
272 0 : for(Int_t i=0; i<18; i++) {
273 0 : hist[14+18+i] = new TH1D(Form("qaTRD_recPoints_sigTimeShape_sm%d", i),
274 0 : Form("sm%d;time bin;signal",i),
275 : 35, -0.5, 34.5);
276 :
277 0 : hist[14+18+i]->SetMaximum(120);
278 : }
279 :
280 : // str = 50
281 0 : for(Int_t i=0; i<18; i++) {
282 0 : hist[50+i] = new TH1D(Form("qaTRD_recPoints_nCls_sm%d",i),
283 0 : Form("sm%d;time bin;number of clusters",i),
284 : 35, -0.5, 34.5);
285 : }
286 :
287 : // str = 68
288 0 : for(Int_t i=0; i<18; i++) {
289 0 : hist[68+i] = new TH1D(Form("qaTRD_recPoints_totalCharge_sm%d", i),
290 0 : Form("sm%d;time bin;total charge", i),
291 : 35, -0.5, 34.5);
292 : }
293 :
294 0 : hist[86] = new TH1D("qaTRD_recPoints_signal", ";amplitude", 200, -0.5, 199.5);
295 0 : hist[87] = new TH2D("qaTRD_recPoints_detMap", ";sm;chamber", 18, -0.5, 17.5, 30, -0.5, 29.5);
296 :
297 0 : for(Int_t i=0; i<18; i++)
298 0 : hist[88+i] = new TH2D(Form("qaTRD_recPoints_XY_sm%d", i),
299 0 : Form("SM%d;Y;X",i), 240, -60, 60, 200, 290, 370);
300 :
301 0 : for(Int_t i=0; i<18; i++)
302 0 : hist[106+i] = new TH2D(Form("qaTRD_recPoints_XPad_sm%d", i),
303 0 : Form("SM%d;Y;X",i), 144, -0.5, 143.5, 200, 290, 370);
304 :
305 :
306 0 : hist[124] = new TH1D("qRef", "ref", 100, 0, 1e4);
307 :
308 0 : fHist = new TObjArray(200);
309 0 : for(Int_t i=0; i<kNhist; i++) {
310 : //hist[i]->Sumw2();
311 0 : fHist->AddAt(hist[i], i);
312 : }
313 :
314 : // reference histograms
315 0 : if (fRef) {
316 0 : for(Int_t i=0; i<540; i++) {
317 0 : fRefHist[i] = new TH2D(Form("refRecPoints_sm%d", i), "",
318 : 16, -0.5, 15.5, 144, -0.5, 143.5); //, 30, -0.5, 29.5);
319 : }
320 0 : } else {
321 :
322 0 : TFile *refFile = new TFile("outRef.root");
323 0 : refFile->cd("TRD/RecPoints");
324 :
325 0 : for(Int_t i=0; i<540; i++) {
326 0 : fRefHist[i] = (TH2D*)gDirectory->Get(Form("refRecPoints_sm%d", i));
327 : }
328 : }
329 0 : }
330 :
331 : //____________________________________________________________________________
332 :
333 : void AliTRDqaRecPoints::AddEvent(TTree * clustersTree)
334 : {
335 : //
336 : // Makes data from RecPoints
337 : //
338 :
339 : // Info("MakeRecPoints", "making");
340 :
341 0 : Int_t nsize = Int_t(clustersTree->GetTotBytes() / (sizeof(AliTRDcluster)));
342 0 : TObjArray *clusterArray = new TObjArray(nsize+1000);
343 :
344 0 : TBranch *branch = clustersTree->GetBranch("TRDcluster");
345 0 : if (!branch) {
346 0 : AliError("Can't get the branch !");
347 0 : return;
348 : }
349 0 : branch->SetAddress(&clusterArray);
350 :
351 : // Loop through all entries in the tree
352 0 : Int_t nEntries = (Int_t) clustersTree->GetEntries();
353 : Int_t nbytes = 0;
354 : AliTRDcluster *c = 0;
355 0 : Int_t nDet[540];
356 0 : for (Int_t i=0; i<540; i++) nDet[i] = 0;
357 :
358 0 : for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
359 :
360 : //printf("Entry = %d\n", iEntry);
361 :
362 : // Import the tree
363 0 : nbytes += clustersTree->GetEvent(iEntry);
364 :
365 : // Get the number of points in the detector
366 0 : Int_t nCluster = clusterArray->GetEntriesFast();
367 :
368 : // Loop through all TRD digits
369 0 : for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
370 0 : c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
371 :
372 0 : Int_t iDet = c->GetDetector();
373 0 : if (iDet < 0 || iDet > 539) continue;
374 :
375 :
376 0 : Int_t iSM = iDet / 30;
377 : //Int_t iStack = iDet % 30;
378 0 : Int_t nPad = c->GetNPads();
379 :
380 0 : if (fnPad && nPad != fnPad) continue;
381 :
382 : //if (iSM == 0 && iStack == 29) continue;
383 : //if (iSM == 8 && iStack == 11) continue;
384 : //if (iSM == 8 && iStack == 7) continue;
385 :
386 0 : Int_t padRow = c->GetPadRow();
387 0 : Int_t padCol = c->GetPadCol();
388 : //Int_t timeBin = c->GetPadTime();
389 :
390 : Double_t refQ = 0;
391 :
392 0 : if (fRef) {
393 0 : fRefHist[iDet]->Fill(padRow, padCol, c->GetQ());
394 0 : } else {
395 0 : Int_t bin = fRefHist[iDet]->FindBin(padRow, padCol);
396 0 : refQ = fRefHist[iDet]->GetBinContent(bin);
397 : //printf("bin = %d\n", bin);
398 : }
399 :
400 0 : ((TH1D*)fHist->At(124))->Fill(refQ);
401 : //printf("ref Q = %lf\n", refQ);
402 :
403 0 : Double_t charge = c->GetQ() - (refQ / (490. * 30));
404 0 : if (charge < 0) continue;
405 :
406 0 : if (charge > 20) {
407 0 : ((TH2D*)fHist->At(88+iSM))->Fill(c->GetY(), c->GetX());
408 0 : ((TH2D*)fHist->At(106+iSM))->Fill(c->GetPadCol(), c->GetX());
409 0 : }
410 :
411 0 : nDet[iDet]++;
412 0 : ((TH1D*)fHist->At(0))->Fill(iDet);
413 0 : ((TH1D*)fHist->At(86))->Fill(charge);
414 0 : ((TH1D*)fHist->At(1))->Fill(iDet, charge);
415 0 : ((TH1D*)fHist->At(2))->Fill(c->GetNPads());
416 0 : if (c->GetNPads() < 6)
417 0 : ((TH1D*)fHist->At(1+c->GetNPads()))->Fill(c->GetCenter());
418 :
419 : //if (c->GetPadTime() < 5)
420 0 : ((TH2D*)fHist->At(7))->Fill(padRow, c->GetPadCol());
421 0 : ((TH1D*)fHist->At(8))->Fill(c->GetPadTime());
422 :
423 0 : ((TH3D*)fHist->At(10))->Fill(iDet, c->GetPadTime(), charge);
424 :
425 0 : ((TH1D*)fHist->At(50+iSM))->Fill(c->GetPadTime());
426 0 : ((TH1D*)fHist->At(68+iSM))->Fill(c->GetPadTime(), charge);
427 :
428 : // PRF for 2pad
429 : //if (c->GetNPads() == 2) {
430 0 : Short_t *sig = c->GetSignals();
431 : Double_t frac = -10;
432 :
433 0 : if (sig[0] == 0 && sig[1] == 0 && sig[2] == 0 && sig[5] == 0 && sig[6] == 0)
434 0 : frac = 1. * sig[4] / (sig[3] + sig[4]);
435 :
436 0 : if (sig[0] == 0 && sig[1] == 0 && sig[4] == 0 && sig[5] == 0 && sig[6] == 0)
437 0 : frac = -1. * sig[2] / (sig[2] + sig[3]);
438 :
439 0 : if (frac > -10) ((TProfile*)fHist->At(11))->Fill(c->GetCenter(), frac);
440 :
441 : //}
442 0 : }
443 : }
444 :
445 0 : for(Int_t i=0; i<540; i++)
446 0 : if (nDet[i] > 0) ((TH1D*)fHist->At(9))->Fill(nDet[i]);
447 :
448 0 : delete clusterArray;
449 :
450 : /*
451 : TFile *outFile = new TFile("outQA.root", "RECREATE");
452 : outFile->mkdir("TRD");
453 : gDirectory->cd("TRD");
454 : gDirectory->mkdir("RecPoints");
455 : gDirectory->cd("RecPoints");
456 : fHist->Write();
457 : outFile->Close();
458 : */
459 0 : }
460 :
461 : //____________________________________________________________________________
|
