Update by Sylwester
[u/mrichter/AliRoot.git] / TRD / AliTRDqaBlackEvents.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * withount fee, provided that the abov copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 /* $Id: AliTRDqaBlackEvents.cxx 23387 2008-01-17 17:25:16Z cblume $ */\r
17 \r
18 ////////////////////////////////////////////////////////////////////////////\r
19 //                                                                        //\r
20 //  QA of black events                                                    //\r
21 //                                                                        //\r
22 //  Author:                                                               //\r
23 //    Sylwester Radomski (radomski@physi.uni-heidelberg.de)               //\r
24 //                                                                        //\r
25 ////////////////////////////////////////////////////////////////////////////\r
26 \r
27 #include "TH1D.h"\r
28 #include "TH2D.h"\r
29 #include "TH2S.h"\r
30 #include "TH3F.h"\r
31 #include "TF1.h"\r
32 #include "TFile.h"\r
33 #include "TCanvas.h"\r
34 #include "TPad.h"\r
35 #include "TLatex.h"\r
36 #include "TStyle.h"\r
37 \r
38 #include "AliTRDgeometry.h"\r
39 #include "AliTRDrawStreamTB.h"\r
40 #include "AliTRDqaBlackEvents.h"\r
41 \r
42 ClassImp(AliTRDqaBlackEvents)\r
43 \r
44 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
45 \r
46 AliTRDqaBlackEvents::AliTRDqaBlackEvents() \r
47   :TObject() \r
48   ,fnEvents(0)\r
49   ,fCreateFull(0)\r
50   ,fThresh(0)\r
51   ,fCount(0)\r
52   ,fOccupancy(0)\r
53   ,fDetRob(0)\r
54   ,fTBEvent(0)\r
55   ,fFitType(0)\r
56   ,fMinNoise(0.5)\r
57   ,fMaxNoise(2) \r
58 {\r
59   //\r
60   // Constructor \r
61   // to create the histograms call Init()\r
62   //\r
63 }\r
64 \r
65 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
66 \r
67 AliTRDqaBlackEvents::AliTRDqaBlackEvents(const AliTRDqaBlackEvents &qa) \r
68   :TObject(qa) \r
69   ,fnEvents(0)\r
70   ,fCreateFull(0)\r
71   ,fThresh(0)\r
72   ,fCount(0)\r
73   ,fOccupancy(0)\r
74   ,fDetRob(0)\r
75   ,fTBEvent(0)\r
76   ,fFitType(0)\r
77   ,fMinNoise(0.5)\r
78   ,fMaxNoise(2) \r
79 {\r
80   //\r
81   // Copy constructor \r
82   // to create the histograms call Init()\r
83   //\r
84 }\r
85 \r
86 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
87 \r
88 void AliTRDqaBlackEvents::Init() \r
89 {\r
90   //\r
91   // creates histograms \r
92   // \r
93 \r
94   //TFile *file = new \r
95   //Info("Init", "Statring");\r
96 \r
97   fnEvents = 0;\r
98 \r
99   // histograms for chambers\r
100   for(Int_t i=0; i<kDET; i++) {\r
101     fNPoint[i]  = new TH2D(Form("entries_%d", i), "",  16, -0.5, 15.5, 144, -0.5, 143.5);\r
102     fData[i]    = new TH3F(Form("data_%d", i), "", 16, -0.5, 15.5, 144, -0.5, 143.5, 50, -0.5, 49.5);\r
103     fChPed[i]   = new TH2D(Form("ped_%d", i), "", 16, -0.5, 15.5, 144, -0.5, 143.5);\r
104     fChNoise[i] = new TH2D(Form("noise_%d", i), "", 16, -0.5, 15.5, 144, -0.5, 143.5);\r
105     fPed[i]     = new TH1D(Form("pedDist_%d", i), ";pedestals (ADC counts)", 100, 5, 15);\r
106 \r
107     fNoise[i]   = new TH1D(Form("noiseDist_%d", i), ";noise (ADC counts)", 100, 0, 5); \r
108     fSignal[i]  = new TH1D(Form("signal_%d", i), ";signal (ADC counts)", 100, -0.5, 99.5);\r
109 \r
110     fnEntriesRM[i] = new TH2D(Form("entriesRM_%d", i), ";ROB,MCM", 8, -0.5, 7.5, 16, -0.5, 15.5);\r
111   }\r
112 \r
113   // histogram for each MCM\r
114   for(Int_t i=0; i < kDET * kROB * kMCM; i++)\r
115     fFullCounter[i] = 0;\r
116 \r
117   // histograms from the whole detector\r
118   fOccupancy = new TH1D("occupancy", "", 20, -0.5, 19.5);\r
119   fDetRob    = new TH2D("DetRob", ";detector;ROB", kDET, -0.5, 539.5, 8, -0.5, 7.5);\r
120   fTBEvent   = new TH2D("tbEvent", ";event ID;time bin", 100, -0.5, 99.5, 30, -0.5, 29.5);\r
121 \r
122   //Info("Init", "Done");\r
123 }\r
124 \r
125 \r
126 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
127 \r
128 void AliTRDqaBlackEvents::Reset() \r
129 {\r
130   //\r
131   // Resets the histograms\r
132   //\r
133 \r
134   for(Int_t i=0; i<kDET; i++) {\r
135     fData[i]->Reset();\r
136     fChPed[i]->Reset();\r
137     fChNoise[i]->Reset();\r
138   }\r
139 }\r
140 \r
141 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
142 \r
143 Int_t AliTRDqaBlackEvents::AddEvent(AliTRDrawStreamTB *data) \r
144 {\r
145   //\r
146   // Add an event\r
147   //\r
148 \r
149   // structure to keep track if particular chanel is used\r
150   Char_t isUsed[kDET][kCOL][kPAD]; \r
151   for(Int_t i=0; i<kDET; i++)\r
152     for(Int_t j=0; j<kCOL; j++)\r
153       for(Int_t k=0; k<kPAD; k++)\r
154         isUsed[i][j][k] = 0;\r
155   \r
156   Int_t nb = 0;\r
157   Int_t rob_last = -1;\r
158   Int_t mcm_last = -1;\r
159   Int_t sm_01 = -1;\r
160 \r
161   while (data->Next()) {\r
162 \r
163     Int_t det = data->GetDet();\r
164 \r
165     Int_t row = data->GetRow();\r
166     Int_t col = data->GetCol();\r
167 \r
168     Int_t rob = data->GetROB();\r
169     Int_t mcm = data->GetMCM();\r
170     Int_t adc = data->GetADC();\r
171 \r
172     Int_t *sig = data->GetSignals();\r
173     nb++;\r
174 \r
175     // ugly hook\r
176     if (det == 0 && row == 0 && col == 0) {\r
177       if (sm_01 > -1) printf("\t\t\t!!! !!! second data set !!! !!!\n");\r
178       sm_01++; \r
179     }\r
180 \r
181     det += sm_01 * 30;   /// ugly\r
182 \r
183     // memory coruption protection\r
184     if (det<0 || det>=kDET) continue;\r
185     \r
186     // check the ROBs\r
187     fDetRob->Fill(det, rob, 1./(kMCM*18));\r
188 \r
189     isUsed[det][row][col]++;\r
190 \r
191     // check if mcm signal is continuus\r
192     if ((rob_last != rob) || (mcm_last != mcm)) {\r
193       rob_last = rob;\r
194       mcm_last = mcm;\r
195       fnEntriesRM[det]->Fill(rob,mcm);\r
196     }\r
197     \r
198     // number of entries for each channels\r
199     fNPoint[det]->Fill(row, col);\r
200     \r
201 \r
202     // create a structure for an MCM if needed\r
203     Int_t mcmIndex = det * (kMCM * kROB) + rob * kMCM + mcm;\r
204     if (fCreateFull && !fFullSignal[mcmIndex])\r
205       fFullSignal[mcmIndex] = new TH2S(Form("mcm_%d_%d_%d", det, rob, mcm), \r
206                                        Form("mcm-%d-%d-%d;ADC;time bin", det, rob,mcm),\r
207                                        21, -0.5, 20.5, 30, -0.5, 29.5);\r
208     \r
209 \r
210     // loop over Time Bins and fill histograms\r
211     for(Int_t k=0; k<kTB; k++) { /// to be corrected\r
212 \r
213       //if (sig[k] < 1) \r
214       //printf("det = %d rob = %d mcm = %d adc = %d k = %d S = %d\n", det, rob, mcm, adc, k, sig[k]);\r
215       \r
216       fSignal[det]->Fill(sig[k]);\r
217       fData[det]->Fill(row, col, sig[k]);\r
218       \r
219       // check if data strange enought\r
220       if (fCreateFull && fFullSignal[mcmIndex]) {\r
221         if (sig[k] > fThresh || sig[k] < 1) fFullCounter[mcmIndex]++;\r
222         fFullSignal[mcmIndex]->Fill(adc, k, sig[k]);\r
223       }\r
224       \r
225       // noisy chamber\r
226       if (det == 29) {\r
227         fTBEvent->Fill(fnEvents, k, sig[k]);\r
228       }\r
229 \r
230     }\r
231   }\r
232   \r
233   // is the dead-alive status changing during the run\r
234   for(Int_t i=0; i<kDET; i++) {\r
235     for(Int_t j=0; j<kCOL; j++)\r
236       for(Int_t k=0; k<kPAD; k++)\r
237         fOccupancy->Fill(isUsed[i][j][k]);\r
238   }\r
239 \r
240   fnEvents++;\r
241   return nb;\r
242 }\r
243 \r
244 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
245 \r
246 void AliTRDqaBlackEvents::Process(const char *filename) \r
247 {\r
248   //\r
249   // Process something\r
250   //\r
251 \r
252   Int_t map[kDET];\r
253   \r
254   TH1D *hist = new TH1D("fitSignal", "", 50, -0.5, 49.5);\r
255   TF1 *fit = new TF1("fit", "gaus(0)", 0, 20);\r
256   fit->SetParameters(1e3, 10, 1);\r
257     \r
258   for(Int_t i=0; i<kDET; i++) {\r
259     \r
260     map[i] = 0;\r
261     if (fData[i]->GetSum() < 10) continue;\r
262     map[i] = 1;\r
263 \r
264     Info("process", "processing chamber %d", i);\r
265 \r
266     for(Int_t j=0; j<fData[i]->GetXaxis()->GetNbins(); j++) {\r
267       for(Int_t k=0; k<fData[i]->GetYaxis()->GetNbins(); k++) {\r
268         \r
269         // project the histogramm\r
270         hist->Reset();\r
271         for(Int_t bb=0; bb<50; bb++) {\r
272           Int_t dataBin = fData[i]->FindBin(j, k, bb);\r
273           Double_t v = fData[i]->GetBinContent(dataBin);\r
274           hist->SetBinContent(bb+1, v);\r
275         }\r
276 \r
277         //TH1D *hist = fData[i]->ProjectionZ(Form("pad_%_%d_%d", i, j, k), j+1, j+1, k+1, k+1);\r
278         \r
279         Int_t bin = fChPed[i]->FindBin(j, k);\r
280 \r
281         if (hist->GetSum() > 1) {\r
282           \r
283           Double_t ped = 0, noise = 0;\r
284 \r
285           if (fFitType == 0) {\r
286             fit->SetParameters(1e3, 10, 1);\r
287             hist->Fit(fit, "q0", "goff", 0, 20);\r
288             TF1 *f = hist->GetFunction("fit");\r
289             ped = TMath::Abs(f->GetParameter(1));\r
290             noise = TMath::Abs(f->GetParameter(2));\r
291           } else {\r
292             ped = hist->GetMean();\r
293             noise = hist->GetRMS();\r
294           }\r
295 \r
296           fChPed[i]->SetBinContent(bin, ped);\r
297           fChNoise[i]->SetBinContent(bin, noise);\r
298           \r
299           fPed[i]->Fill(ped);\r
300           fNoise[i]->Fill(noise);\r
301 \r
302         } else {\r
303           fChPed[i]->SetBinContent(bin, 0);\r
304           fChNoise[i]->SetBinContent(bin, 0);\r
305         }\r
306         \r
307         //delete hist;\r
308       }\r
309     }\r
310   }\r
311 \r
312   Info("Process", "Number of events = %d", fnEvents);\r
313 \r
314   // normalize number of entries histos\r
315   Int_t max = 0;\r
316   for(Int_t i=0; i<kDET; i++) { \r
317     if (!map[i]) continue;\r
318     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {\r
319       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {\r
320         Int_t dataBin = fNPoint[i]->FindBin(j, k);\r
321         Double_t v = fNPoint[i]->GetBinContent(dataBin);\r
322         if (v > max) max = (Int_t)v;\r
323       }\r
324     }\r
325   }\r
326   \r
327   char entriesDistName[100];\r
328   \r
329   for(Int_t i=0; i<kDET; i++) {\r
330     \r
331     if (!map[i]) continue;\r
332     \r
333     sprintf(entriesDistName, "entriesDist_%d", i);\r
334     fNPointDist[i] = new TH1D(entriesDistName, ";number of events", max+2, -0.5, max+1.5);\r
335     \r
336     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {\r
337       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {\r
338         Int_t dataBin = fNPoint[i]->FindBin(j, k);\r
339         Double_t v = fNPoint[i]->GetBinContent(dataBin);\r
340         //if (v > fnEvents) printf("N = %d V = %lf\n", fnEvents, v);\r
341         fNPointDist[i]->Fill(v); \r
342       }\r
343     }\r
344     \r
345     fNPoint[i]->Scale(1./fnEvents);\r
346   }\r
347   \r
348 \r
349   for(Int_t i=0; i<kDET; i++) {\r
350     fnEntriesRM[i]->SetMaximum(fnEvents * 1.5);\r
351   }\r
352 \r
353   // save histograms\r
354 \r
355   TFile *file = new TFile(filename, "recreate");\r
356   for(Int_t i=0; i<kDET; i++) {\r
357     if (!map[i]) continue; \r
358     fChPed[i]->Write();\r
359     fChNoise[i]->Write();\r
360     fNPoint[i]->Write();\r
361     fNPointDist[i]->Write();\r
362     fPed[i]->Write();\r
363     fNoise[i]->Write();\r
364     fSignal[i]->Write();\r
365     fnEntriesRM[i]->Write();\r
366   }\r
367 \r
368   Int_t nMcm = 0;\r
369   for(Int_t i=0; i < kDET * kROB * kMCM; i++) {\r
370     if (fFullSignal[i] && fFullCounter[i] > fCount) {\r
371       fFullSignal[i]->Write();\r
372       nMcm++;\r
373     }\r
374   }\r
375 \r
376   printf("Number of saved MCMs = %d\n", nMcm);\r
377 \r
378   fOccupancy->Write();\r
379   fDetRob->Write();\r
380   fTBEvent->Write();\r
381   file->Close();\r
382   delete file;\r
383 }\r
384 \r
385 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
386 \r
387 void AliTRDqaBlackEvents::DrawChamber(const char *filename, Int_t det, Int_t w, Int_t h) \r
388 {\r
389   //\r
390   // Draw raport for one chamber: \r
391   // pedestal map, noise map, distribution of pedestal and noise\r
392   // \r
393   // input:\r
394   // name of the file with histograms (created with Process())\r
395   // detector Id (0 - 539)\r
396   // \r
397 \r
398   // setup global style\r
399   gStyle->SetPalette(1);\r
400   gStyle->SetOptStat(0);\r
401   gStyle->SetPadTopMargin(0.02);\r
402   gStyle->SetPadBottomMargin(0.05);\r
403 \r
404   TFile *file = new TFile(filename, "READ");\r
405 \r
406   TCanvas *c = new TCanvas("blackEvents",Form("blackEvents %d",det), w, h);\r
407   c->SetVertical(kFALSE);\r
408   c->Divide(3,1, 0.01, 0.01);\r
409   c->cd(3);\r
410   \r
411   TPad *mPad = (TPad*) gPad;\r
412   mPad->Divide(1,2,0.01,0.01);\r
413   \r
414   c->cd(1);\r
415   TH2D *h2 = (TH2D*)file->Get(Form("ped_%d",det));\r
416   h2->SetMinimum(5);\r
417   h2->SetMaximum(15);\r
418   h2->SetTitle(";Z direction;#phi direction");\r
419   h2->Draw("colz");\r
420   \r
421   c->cd(2);\r
422   h2 = (TH2D*)file->Get(Form("noise_%d",det));\r
423   h2->SetMinimum(fMinNoise);\r
424   h2->SetMaximum(fMaxNoise);\r
425   h2->SetTitle(";Z direction;#phi direction");\r
426   h2->Draw("colz");\r
427   \r
428   mPad->cd(1);\r
429   //gPad->SetLogy();\r
430   TH1D *h1 = (TH1D*)file->Get(Form("pedDist_%d", det));\r
431   h1->Draw();\r
432   \r
433   mPad->cd(2);\r
434   gPad->SetLogy();\r
435   h1 = (TH1D*)file->Get(Form("noiseDist_%d", det));\r
436   h1->Draw();                    \r
437   \r
438   h1->Fit("gaus");\r
439   TF1 *f = h1->GetFunction("gaus");\r
440   const char *tt = Form("#mu = %.2f #sigma = %0.2f ", f->GetParameter(1),f->GetParameter(2));\r
441   TLatex *ll = new TLatex(2, 100, tt);\r
442   ll->SetTextSize(0.06);\r
443   ll->Draw();\r
444 }\r
445 \r
446 ///////////////////////////////////////////////////////////////////////////////////////////////////\r
447 \r
448 void AliTRDqaBlackEvents::DrawSm(const char *filename, Int_t sm, Int_t w, Int_t h) \r
449 {\r
450   //\r
451   // ????????????\r
452   //\r
453   \r
454   gStyle->SetPalette(1);\r
455   gStyle->SetOptStat(0);\r
456   \r
457   gStyle->SetPadTopMargin(0.02);\r
458   //gStyle->SetPadBottomMargin(0.05);  \r
459   //gStyle->SetPadLeftMargin(0.02);  \r
460   //gStyle->SetPadRightMargin(0.02);\r
461 \r
462   TFile *file = new TFile(filename, "READ");\r
463 \r
464   TCanvas *c = new TCanvas("blackEventsSM",Form("blackEvents SM %d",sm), w, h);\r
465   c->SetVertical(kFALSE);\r
466   c->Divide(5, 6, 0.001, 0.01);\r
467   \r
468   for(Int_t i=0; i<30; i++) {\r
469     \r
470     TH2D *h2 = (TH2D*)file->Get(Form("noise_%d",i+30*sm));\r
471     if (!h2) continue;\r
472     h2->SetMinimum(fMinNoise);\r
473     h2->SetMaximum(fMaxNoise);\r
474 \r
475     // to be replaced by the official calculation\r
476     Int_t stack = i/6;\r
477     Int_t layer = i%6;\r
478     Int_t index = (5-layer)*5 + stack + 1;\r
479     //printf("%d %d %d %d\n", i, stack, layer, index);\r
480     c->cd(index);\r
481     gPad->SetBottomMargin(0.02);\r
482     gPad->SetTopMargin(0.02);\r
483 \r
484     h2->Draw("col");\r
485   }\r
486 }\r
487 \r
488 ///////////////////////////////////////////////////////////////////////////////////////////////////\r