Rename AliTRDrawStream to AliTRDrawStreamOld
[u/mrichter/AliRoot.git] / TRD / AliTRDqaBlackEvents.cxx
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  * withount fee, provided thats the abov 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: AliTRDqaBlackEvents.cxx 23387 2008-01-17 17:25:16Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  QA of black events                                                    //
21 //                                                                        //
22 //  Author:                                                               //
23 //    Sylwester Radomski (radomski@physi.uni-heidelberg.de)               //
24 //                                                                        //
25 ////////////////////////////////////////////////////////////////////////////
26
27 #include "TH1D.h"
28 #include "TH2D.h"
29 #include "TH2S.h"
30 //#include "TH3F.h"
31 #include "TF1.h"
32 #include "TFile.h"
33 #include "TCanvas.h"
34 #include "TPad.h"
35 #include "TLatex.h"
36 #include "TStyle.h"
37 #include "TGraph.h"
38
39 #include "AliRawReader.h"
40
41 #include "AliTRDrawStreamOld.h"
42 #include "AliTRDqaBlackEvents.h"
43
44 ClassImp(AliTRDqaBlackEvents)
45
46 ///////////////////////////////////////////////////////////////////////////////////////////////////
47
48 AliTRDqaBlackEvents::AliTRDqaBlackEvents() 
49   :TObject() 
50   ,fnEvents(0)
51   ,fCreateFull(0)
52   ,fThresh(0)
53   ,fCount(0)
54   ,fRefEv(0)
55   ,fOccupancy(0)
56   ,fDetRob(0)
57   ,fTBEvent(0)
58   ,fRefHistPed(0)
59   ,fRefHistNoise(0)
60   ,fErrorHC(0)
61   ,fErrorMCM(0)
62   ,fErrorADC(0)
63   ,fErrorSMHC(0)
64   ,fErrorSMMCM(0)
65   ,fErrorSMADC(0)
66   ,fErrorGraphHC(0)
67   ,fErrorGraphMCM(0)
68   ,fErrorGraphADC(0)
69   ,fGraphMCM(0)
70   ,fMcmTracks(0)
71   ,fMapMCM(0)
72   ,fFracMCM(0)
73   ,fSMHCped(0)
74   ,fSMHCerr(0)
75   ,fNoiseTotal(0)
76   ,fPP(0)
77   ,fMinNoise(0.5)
78   ,fMaxNoise(2)
79   ,fFitType(0) 
80   //  ,fRefFileName("")
81 {
82   //
83   // Constructor 
84   // to create the histograms call Init()
85   //
86
87   strcpy(fRefFileName, "");
88 }
89
90 ///////////////////////////////////////////////////////////////////////////////////////////////////
91
92 AliTRDqaBlackEvents::AliTRDqaBlackEvents(const AliTRDqaBlackEvents &qa) 
93   :TObject(qa) 
94   ,fnEvents(0)
95   ,fCreateFull(0)
96   ,fThresh(0)
97   ,fCount(0)
98   ,fRefEv(0)
99   ,fOccupancy(0)
100   ,fDetRob(0)
101   ,fTBEvent(0)
102   ,fRefHistPed(0)
103   ,fRefHistNoise(0)
104   ,fErrorHC(0)
105   ,fErrorMCM(0)
106   ,fErrorADC(0)
107   ,fErrorSMHC(0)
108   ,fErrorSMMCM(0)
109   ,fErrorSMADC(0)
110   ,fErrorGraphHC(0)
111   ,fErrorGraphMCM(0)
112   ,fErrorGraphADC(0)
113   ,fGraphMCM(0)
114   ,fMcmTracks(0)
115   ,fMapMCM(0)
116   ,fFracMCM(0)
117   ,fSMHCped(0)
118   ,fSMHCerr(0)
119   ,fNoiseTotal(0)
120   ,fPP(0)
121   ,fMinNoise(0.5)
122   ,fMaxNoise(2) 
123   ,fFitType(0)
124    //,fRefFileName("")
125 {
126   //
127   // Copy constructor 
128   // to create the histograms call Init()
129   //
130   
131   strcpy(fRefFileName, "");
132 }
133
134 ///////////////////////////////////////////////////////////////////////////////////////////////////
135
136 void AliTRDqaBlackEvents::Init() 
137 {
138   //
139   // creates histograms 
140   // 
141
142   //TFile *file = new 
143   //Info("Init", "Statring");
144
145   fnEvents = 0;
146
147   // histograms for chambers
148   for(Int_t det=0; det<kDET; det++) {
149
150     fNPoint[det]  = new TH2D(Form("entries_%d", det), "",  16, -0.5, 15.5, 144, -0.5, 143.5);
151     //fData[det]    = new TH3F(Form("data_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5, 50, -0.5, 49.5);
152
153     // pedestal noise maps using RMS and Fit
154     fChPed[det]   = new TH2D(Form("ped_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
155     fChNoise[det] = new TH2D(Form("noise_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
156     
157     //fChPed[det]   = new TH2D(Form("ped_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
158     //fChNoise[det] = new TH2D(Form("noise_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);    
159
160     // distribution per detector
161     fPed[det]     = new TH1D(Form("pedDist_%d", det), ";pedestals (ADC counts)", 100, 5, 15);
162     fNoise[det]   = new TH1D(Form("noiseDist_%d", det), ";noise (ADC counts)", 100, 0, 5); 
163     fSignal[det]  = new TH1D(Form("signal_%d", det), ";signal (ADC counts)", 100, -0.5, 99.5);
164     fChPP[det]    = new TH1D(Form("pp_%d", det), ";pp (ADC)", 200, -0.5, 199.5);
165
166     fnEntriesRM[det] = new TH2D(Form("entriesRM_%d", det), ";ROB,MCM", 8, -0.5, 7.5, 16, -0.5, 15.5);
167     
168     // histograms after reference subtraction
169     fChPedRes[det]    = new TH2D(Form("pedRef_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
170     fChNoiseRes[det]  = new TH2D(Form("noiseRef_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
171
172
173     // error codes
174     fErrorLocMCM[det] = new TH2D(Form("errorLocMCM_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
175     fErrorLocADC[det] = new TH2D(Form("errorLocADC_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);    
176     fErrorLocHC[det]  = new TH2D(Form("errorLocHC_%d", det), "", 16, -0.5, 15.5, 144, -0.5, 143.5);
177   }
178
179   // histogram for each MCM
180   for(Int_t i=0; i < kDET * kROB * kMCM; i++)
181     fFullCounter[i] = 0;
182
183   // histograms from the whole detector
184   fOccupancy = new TH1D("occupancy", "", 20, -0.5, 19.5);
185   fDetRob    = new TH2D("DetRob", ";detector;ROB", kDET, -0.5, 539.5, 8, -0.5, 7.5);
186   fTBEvent   = new TH2D("tbEvent", ";event ID;time bin", 100, -0.5, 99.5, 30, -0.5, 29.5);
187
188   // errors statistics and location
189   fErrorHC  = new TH1D("errorHC", ";error ID;", 18, -3.5, 14.5);
190   fErrorMCM = new TH1D("errorMCM", ";error ID;", 18, -3.5, 14.5);
191   fErrorADC = new TH1D("errorADC", ";error ID;", 18, -3.5, 14.5);
192   
193   fErrorSMHC  = new TH1D("errorSM_HC", ";SM id", 18, -0.5, 17.5);
194   fErrorSMMCM = new TH1D("errorSM_MCM", ";SM id", 18, -0.5, 17.5);
195   fErrorSMADC = new TH1D("errorSM_ADC", ";SM id", 18, -0.5, 17.5);
196
197
198   fErrorGraphHC  = new TGraph();
199   fErrorGraphMCM = new TGraph();
200   fErrorGraphADC = new TGraph();
201
202   fGraphMCM = new TGraph();
203
204   for(Int_t i=0; i<3; i++) {
205     fGraphPP[i] = new TGraph();
206   }
207
208
209   fMapMCM = new TH2D("mapMCM", ";det;mcm", 540, -0.5, 539.5, kROB*kMCM, -0.5, kROB*kMCM-0.5);
210   fFracMCM = new TH1D("fracMCM", ";frequency", 100, 0, 1);
211   
212
213   fErrorGraphHC->GetHistogram()->SetTitle("fraction of events with HC error;event number");
214   fErrorGraphMCM->GetHistogram()->SetTitle("fraction of events with MCM error;event number;");
215   fErrorGraphADC->GetHistogram()->SetTitle("fraction of events with ADC error;event number;"); 
216
217
218   fSMHCped = new TH2D("smHcPed", ";super module;half chamber", 18, -0.5, 17.5, 60, -0.5, 59.5);
219   
220   // link monitor
221   const char *linkName[3] = {"smLink", "smBeaf", "smData"};
222   const char *linkGrName[3] = {"grSmLink", "grSmBeaf", "grSmData"};
223   for(Int_t i=0; i<3; i++) {
224     fSMLink[i] = new TH2D(linkName[i], ";super module;link", 18, -0.5, 17.5, 60, -0.5, 59.5);
225     fGrLink[i] = new TGraph();
226     fGrLink[i]->SetName(linkGrName[i]);
227   }
228
229   //fZSsize = new TH1D("zssizeSingle", ";threshold;nADC", 40, -0.5, 39.5);
230   
231
232   //Info("Init", "Done");
233
234   // number of ADC channels fired per SM and in total
235   for(Int_t sm=0; sm<kSM+1; sm++)
236     fNumberADC[sm] = new TGraph();
237
238   //
239   fNoiseTotal = new TH1D("noiseTotal", "noise (ADC)", 250, 0, 10);
240   fPP = new TH1D("peakPeak", "p-p (ADC)", 200, -0.5, 199.5);
241
242   for(Int_t sm=0; sm<kSM; sm++) {
243     fSmNoiseRms[sm] = new TH1D(Form("noiseRms_sm%d", sm), ";noise from RMS (ADC)", 100, 0, 10);
244     fSmNoiseFit[sm] = new TH1D(Form("noiseFit_sm%d", sm), ";noise frim Fit (ADC)", 100, 0, 10);
245     fSmPP[sm] = new TH1D(Form("peakPeak_sm%d", sm), ";peak-peak (ADC)", 200, -0.5, 199.5); 
246   }
247
248   // event number consistancy
249   for(Int_t i=0; i<1000; i++) {
250     fEvNoDist[i] = new TH1D(Form("mcmEvDist_%d", i), ";#Delta Events", 201, -100.5, 100.5);
251   }
252   fRefEv = -1;
253
254   fMcmTracks = new TObjArray();
255
256   // clean data direct
257   for(Int_t det=0; det<kDET; det++) 
258     for(Int_t row=0; row<kROW; row++)
259       for(Int_t pad=0; pad<kPAD; pad++)
260         for(Int_t ch=0; ch<kCH; ch++) {
261           fDataDirect[det][row][pad][ch] = 0;
262           fSignalDirect[det][ch] = 0;  // overdone
263         }
264 }
265
266
267 ///////////////////////////////////////////////////////////////////////////////////////////////////
268
269 void AliTRDqaBlackEvents::Reset() 
270 {
271   //
272   // Resets the histograms
273   //
274
275   for(Int_t i=0; i<kDET; i++) {
276     //fData[i]->Reset();
277     fChPed[i]->Reset();
278     fChNoise[i]->Reset();
279   }
280 }
281
282
283 ///////////////////////////////////////////////////////////////////////////////////////////////////
284
285 void AliTRDqaBlackEvents::SetRefFile(const char *filename) {
286   
287   strcpy(fRefFileName, filename);
288 }
289
290 ///////////////////////////////////////////////////////////////////////////////////////////////////
291
292 void AliTRDqaBlackEvents::ReadRefHists(Int_t det) {
293   //  
294   // Read the reference histograms
295   //
296
297   fRefHistPed = 0;
298   fRefHistNoise = 0;
299   
300   TFile *file = 0;
301   if (fRefFileName) TFile::Open(fRefFileName);
302   if (!file) return;
303
304   fRefHistPed   = (TH2D*)file->Get(Form("ped_%d",det));
305   fRefHistNoise = (TH2D*)file->Get(Form("noise_%d", det));
306
307   if (file) file->Close();
308 }
309
310 ///////////////////////////////////////////////////////////////////////////////////////////////////
311
312 void AliTRDqaBlackEvents::StartEvent()
313 {
314   //
315   // start an event
316   //
317
318   // clear the mcm data
319   for(Int_t i=0; i < kDET * kROB * kMCM; i++) {
320     if (fFullSignal[i]) fFullSignal[i]->Reset();
321     fFullCounter[i] = 0;
322   }
323
324   for(Int_t i=0; i<2; i++) {
325     fnErrorHC[i] = 0;
326     fnErrorMCM[i] = 0;
327     fnErrorADC[i] = 0;
328   }
329
330  
331   Int_t ppThresh[3] = {10, 20, 40};
332   for(Int_t i=0; i<3; i++) {
333     fppThresh[i] = ppThresh[i];
334     fnPP[i] = 0;
335     fnLink[i] = 0;
336   }
337
338   for(Int_t sm=0; sm<kSM+1; sm++) fnADCinSM[sm] = 0;
339
340   if (fRefEv > 0) fRefEv++;
341   fEvNoDist[999]->Reset();  // keep only the last event
342
343 }
344
345 ///////////////////////////////////////////////////////////////////////////////////////////////////
346
347 void AliTRDqaBlackEvents::AddBuffer(AliTRDrawStreamOld *data, AliRawReader * const reader) 
348 {
349   
350   //printf ("try to read data\n");
351   Int_t nextBuff  = data->NextBuffer();
352   //printf("done ...\n");
353
354   if (nextBuff == 0) return;
355  
356   Int_t sm = reader->GetEquipmentId() - 1024;
357   //printf("reading SM %d\n", sm);
358   AliInfo(Form("reading SM %d", sm));
359   
360   if (sm < 0 || sm > 18) return;
361
362   // lopp over stacks, links ...
363
364   for (Int_t istack = 0; istack < 5; istack++) {        
365     for (Int_t ilink = 0; ilink < 12; ilink++) {
366       
367       //printf("HC = %d %d\n", istack, ilink);
368
369       Int_t det = sm * 30 + istack * 6 + ilink/2;       
370       
371       // check if data delivered
372       if (!(data->IsLinkActiveInStack(istack, ilink))) continue;
373       fSMLink[0]->Fill(sm, istack * 12 + ilink);
374       fnLink[0]++;
375
376       // check if beaf-beaf
377       if (data->GetLinkMonitorError(istack, ilink)) {
378         fSMLink[1]->Fill(sm, istack * 12 + ilink);
379         fnLink[1]++;
380         continue;
381       }
382         
383       // fill histogram with HC header errors
384       Int_t nErrHc = 0;
385       Int_t nErrHcTot = 0;
386       
387       nErrHc = FillBits(fErrorHC, data->GetH0ErrorCode(istack, ilink), 0);
388       if (!nErrHc) fErrorHC->Fill(-3);
389       nErrHcTot += nErrHc;
390       
391       nErrHc = FillBits(fErrorHC, data->GetH1ErrorCode(istack, ilink), 2);
392       if (!nErrHc) fErrorHC->Fill(-2);
393       nErrHcTot += nErrHc;
394       
395       nErrHc = FillBits(fErrorHC, data->GetHCErrorCode(istack, ilink), 4);
396       if (!nErrHc) fErrorHC->Fill(-1);
397       nErrHcTot += nErrHc;
398       
399       // trending
400       fnErrorHC[0]++;
401       if (nErrHcTot > 0) { 
402         fnErrorHC[1]++;
403         fErrorSMHC->Fill(sm);
404       }
405       
406       // data integrity protection 
407
408       //if (data->GetHCErrorCode(istack, ilink) > 0) continue;
409       if (data->GetH0ErrorCode(istack, ilink) > 0) continue;
410       if (data->GetH1ErrorCode(istack, ilink) > 0) continue;
411       
412       fSMLink[2]->Fill(sm, istack * 12 + ilink);        
413       fnLink[2]++;
414              
415       
416       for (Int_t imcm = 0; imcm < data->GetHCMCMmax(istack, ilink); imcm++ ){
417           
418         //printf("mcm = %d %d %d\n", istack, ilink, imcm);
419         
420         // fill MCM error code
421         
422         Int_t nErrMcm = 0;
423         Int_t nErrMcmTot = 0;
424         
425         nErrMcm = FillBits(fErrorMCM, data->GetMCMhdErrorCode(istack, ilink, imcm), 0);
426         if (!nErrMcm) fErrorMCM->Fill(-3);
427         nErrMcmTot += nErrMcm;
428         
429         nErrMcm = FillBits(fErrorMCM, data->GetMCMADCMaskErrorCode(istack, ilink, imcm), 5);
430         if (!nErrMcm) fErrorMCM->Fill(-2);
431         nErrMcmTot += nErrMcm;
432         
433         nErrMcm = FillBits(fErrorMCM, data->GetMCMErrorCode(istack, ilink, imcm), 10);
434         if (!nErrMcm) fErrorMCM->Fill(-1);
435         nErrMcmTot += nErrMcm;                    
436         
437         // trending
438         fnErrorMCM[0]++;
439         if (nErrMcmTot > 0) { 
440           fnErrorMCM[1]++;
441           fErrorSMMCM->Fill(sm);
442         }
443         
444         // MCM protection
445         if ( (data->GetMCMhdErrorCode(istack,ilink,imcm)) & 2 ) continue;
446         //if ((data->GetMCMADCMaskErrorCode(istack,ilink,imcm))) continue;
447         //if ((data->GetMCMErrorCode(istack,ilink,imcm))) continue;
448         
449         Int_t mcmEvent = data->GetEventNumber(istack, ilink, imcm);
450         
451         // set the reference event number 
452         if (fRefEv < 0) {
453           fRefEv = mcmEvent;
454           printf("Reference Event Number = %d (%d %d %d)\n", fRefEv, istack, ilink, imcm);
455         }
456         
457         // fill event distribution
458         if (!(fnEvents%10)) {
459           fEvNoDist[fnEvents/10]->Fill(mcmEvent - fRefEv);
460         }
461         
462         fEvNoDist[999]->Fill(mcmEvent - fRefEv);
463         
464         Int_t mcm = data->GetMCM(istack, ilink, imcm);
465         Int_t rob = data->GetROB(istack, ilink, imcm);
466
467         // create a structure for an MCM if needed
468         Int_t mcmIndex = det * (kMCM * kROB) + rob * kMCM + mcm;
469         if (fCreateFull && !fFullSignal[mcmIndex])
470           fFullSignal[mcmIndex] = 
471             new TH2S(Form("mcm_%d_%d_%d_%d_%d", sm, istack, ilink/2, rob, mcm), 
472                      Form("mcm-%d-%d-%d-%d-%d;ADC;time bin", sm, istack, ilink/2, rob, mcm),
473                      21, -0.5, 20.5, 30, -0.5, 29.5);
474         
475         
476         //Int_t zsADC[21][40];
477         /*
478           for(Int_t ina=0; ina<21; ina++) 
479           for(Int_t th=0; th<40; th++)
480           zsADC[ina][th] = 0;
481         */
482         
483         // first loop over ADC chanels 
484         
485         for (Int_t iadc=0; iadc < data->GetADCcount(istack, ilink, imcm); iadc++) {
486           
487           //printf("ADC = %d\n", iadc);
488           
489
490           // fill ADC error bits
491           Int_t nErrAdc = FillBits(fErrorADC, data->GetADCErrorCode(), 0);
492           if (!nErrAdc) fErrorADC->Fill(-1);
493           
494           fnErrorADC[0]++;
495           if (nErrAdc > 0) {
496             fnErrorADC[1]++;
497             fErrorSMADC->Fill(sm);
498           }
499           
500           // ADC protection
501           if ((data->GetADCErrorCode(istack,ilink,imcm,iadc))) continue;
502
503           Int_t minV = 1024;
504           Int_t maxV = 0;
505             
506           Int_t *sig = data->GetSignalDirect(istack, ilink, imcm, iadc);
507
508           //Int_t adc = data->GetADCnumber(istack, ilink, imcm, iadc);
509           Int_t row = data->GetRow(istack, ilink, imcm);
510           Int_t col = data->GetCol(istack, ilink, imcm, iadc);              
511             
512           // loop over Time Bins and fill histograms
513           for(Int_t k=0; k < data->GetNumberOfTimeBins(istack, ilink); k++) { 
514                       
515             //fSignal[det]->Fill(sig[k]);
516             //fData[det]->Fill(row, col, sig[k]); // slow
517
518             if (sig[k] < kCH) {
519               fSignalDirect[det][sig[k]]++;
520               fDataDirect[det][row][col][sig[k]]++; // direct data
521             }
522
523             // peak-peak
524             minV = (minV < sig[k]) ? minV : sig[k];
525             maxV = (maxV > sig[k]) ? maxV : sig[k];
526             
527             // check for active MCMs
528             if (fCreateFull && fFullSignal[mcmIndex]) {
529               if (sig[k] > fThresh || sig[k] < 0) fFullCounter[mcmIndex]++;
530               //if (sm == 0 && istack == 0 && ilink/2 == 1 && rob == 1 && mcm == 15) fFullCounter[mcmIndex]++; // special
531               //fFullSignal[mcmIndex]->Fill(adc, k, sig[k]); // slow
532             }
533             
534             // zero suppresion tests
535             /*
536               for(Int_t th=0; th<40; th++) {
537               if (sig[k] > th) {
538               zsADC[iadc][th] = 1;
539               if (iadc > 0) zsADC[iadc-1][th] = 1;
540               if (iadc < 10) zsADC[iadc+1][th] = 1;
541               }
542               }
543             */
544             
545           } // tb
546           
547           if (maxV > 0) { 
548             fnADCinSM[sm]++;
549             fnADCinSM[kSM]++;
550           }
551           
552           Int_t adcPP = maxV - minV;
553           //if (adcPP == 100) fFullCounter[mcmIndex] += 10;
554
555           fPP->Fill(adcPP);
556           fChPP[det]->Fill(adcPP);
557           fSmPP[sm]->Fill(adcPP);
558           
559           for(Int_t i=0; i<3; i++) {
560             if ((adcPP) > fppThresh[i]) fnPP[i]++;
561           }
562           
563             
564         } // adc
565
566           // fill ZS histos
567           /*
568             for(Int_t th=0; th<40; th++) {
569             Int_t nnADC = 0;
570             for(Int_t ins=0; ins<21; ins++) 
571             nnADC += zsADC[ins][th];
572             fZSsize->Fill(th, nnADC);
573             }
574           */
575
576         // fill active MCMs
577
578         if (fCreateFull && fFullSignal[mcmIndex] && (fFullCounter[mcmIndex] > fCount)) {
579           
580           for (Int_t iadc=0; iadc < data->GetADCcount(istack, ilink, imcm); iadc++) {
581             
582             // ADC protection
583             if ((data->GetADCErrorCode(istack,ilink,imcm,iadc))) continue;        
584             
585             //Int_t row = data->GetRow(istack, ilink, imcm);
586             //Int_t col = data->GetCol(istack, ilink, imcm, iadc);                  
587             Int_t adc = data->GetADCnumber(istack, ilink, imcm, iadc);
588             
589             Int_t *sig = data->GetSignalDirect(istack, ilink, imcm, iadc);
590             
591             // loop over Time Bins and fill histograms
592             for(Int_t k=0; k < data->GetNumberOfTimeBins(istack, ilink); k++) { 
593               fFullSignal[mcmIndex]->Fill(adc, k, sig[k]); // slow
594             }
595           } // tb
596         }
597         
598       } // mcm 
599     } // link
600   } // stack  
601   
602   // printf("end of loops\n");
603
604
605
606 ///////////////////////////////////////////////////////////////////////////////////////////////////
607
608 void AliTRDqaBlackEvents::FinishEvent()
609 {
610   //
611   // Processing at the end of the current event
612   //
613
614   for(Int_t i=0; i<3; i++) {
615     fGraphPP[i]->SetPoint(fnEvents, fnEvents, fnPP[i]);
616   }
617
618   // trend of the number of links
619   for(Int_t i=0; i<3; i++) {
620     fGrLink[i]->SetPoint(fnEvents, fnEvents, fnLink[i]);
621   }
622
623   // save interesting histos
624   Int_t mcmTrackCandidate = 0;
625   for(Int_t i = 0; i < kDET * kROB * kMCM; i++) { 
626     if ((fFullCounter[i] > fCount) && fFullSignal[i] && CheckMCM(i) )  {
627       
628       fMcmTracks->AddLast(fFullSignal[i]->Clone(Form("event_%d_%s", fnEvents, fFullSignal[i]->GetName())));
629       mcmTrackCandidate++;
630       
631       Int_t mcmTrackletDet = i/(kROB * kMCM);   
632       Int_t mcmTrackletMcm = i%(kROB * kMCM);
633       fMapMCM->Fill(mcmTrackletDet, mcmTrackletMcm);
634     }
635   }
636   
637   fGraphMCM->SetPoint(fnEvents, fnEvents, mcmTrackCandidate);
638   AliInfo(Form("Number of MCM track candidates = %d\n", mcmTrackCandidate));
639   
640   
641   // update fraction of error graphs
642   Double_t err;
643   
644   err = (fnErrorHC[0] > 0)? 100.*fnErrorHC[1]/fnErrorHC[0] : -1;
645   fErrorGraphHC->SetPoint(fnEvents, fnEvents, err);
646   
647   err = (fnErrorMCM[0] > 0)? 100.*fnErrorMCM[1]/fnErrorMCM[0] : -1;
648   fErrorGraphMCM->SetPoint(fnEvents, fnEvents, err);
649   
650   err = (fnErrorADC[0] > 0)? 100.*fnErrorADC[1]/fnErrorADC[0] : -1;
651   fErrorGraphADC->SetPoint(fnEvents, fnEvents, err);
652
653   // number of fired ADC per SM
654   for(Int_t sm=0; sm<kSM+1; sm++) 
655     fNumberADC[sm]->SetPoint(fnEvents, fnEvents, fnADCinSM[sm]);
656
657   fnEvents++;
658 }
659
660 ///////////////////////////////////////////////////////////////////////////////////////////////////
661
662 void AliTRDqaBlackEvents::Process(const char *filename) 
663 {
664   //
665   // Process something
666   //
667   
668   char fn[256];
669   strcpy(fn, filename);
670   
671   //AliInfo(Form("FILENAME = %s (%s)\n", filename, fn));
672
673   Int_t map[kDET];
674   
675   TH1D *hist = new TH1D("fitSignal", "", 50, -0.5, 49.5);
676   TF1 *fit = new TF1("fit", "gaus(0)", 0, 20);
677   fit->SetParameters(1e3, 10, 1);
678     
679   for(Int_t det=0; det<kDET; det++) {
680
681     //AliInfo(Form("processing chamber %d\n", det));   
682
683     map[det] = 0;
684     //if (fData[det]->GetSum() < 10) continue;
685     //if (fDataDirect[det][10][10][10] < 20) continue;
686     //map[det] = 1;
687
688
689     // rewrite signal-direct
690     for(Int_t ch=0; ch<kCH; ch++) {
691       fSignal[det]->Fill(ch, fSignalDirect[det][ch]);
692     }
693
694     // read reference distributions
695     ReadRefHists(det);
696
697     //for(Int_t row=0; row<fData[det]->GetXaxis()->GetNbins(); row++) {
698     //for(Int_t pad=0; pad<fData[det]->GetYaxis()->GetNbins(); pad++) {
699         
700     for(Int_t row=0; row<kROW; row++) {
701       for(Int_t pad=0; pad<kPAD; pad++) {
702
703         // project the histogramm
704         hist->Reset();
705         //for(Int_t bb=0; bb<50; bb++) {
706         for(Int_t bb=0; bb<kCH; bb++) {
707           //Int_t dataBin = fData[det]->FindBin(row, pad, bb);
708           //Double_t v = fData[det]->GetBinContent(dataBin);
709           hist->SetBinContent(bb+1, fDataDirect[det][row][pad][bb]);
710         }
711
712         Int_t bin = fChPed[det]->FindBin(row, pad);
713
714         if (hist->GetSum() > 1) {
715           
716           map[det] = 1;
717           Double_t ped = 0, noise = 0;
718
719           if (fFitType == 0) {
720             fit->SetParameters(1e3, 10, 1);
721             hist->Fit(fit, "q0", "goff", 0, 20);
722             TF1 *f = hist->GetFunction("fit");
723             ped = TMath::Abs(f->GetParameter(1));
724             noise = TMath::Abs(f->GetParameter(2));
725             fSmNoiseFit[det/30]->Fill(noise);
726           } else {
727             ped = hist->GetMean();
728             noise = hist->GetRMS();
729             fSmNoiseRms[det/30]->Fill(noise);
730             //if (pad == 0)
731             //  AliInfo(Form("data %f %f %f\n", hist->GetSum(), ped, noise));
732           }
733
734           fChPed[det]->SetBinContent(bin, ped);
735           fChNoise[det]->SetBinContent(bin, noise);
736           fNoiseTotal->Fill(noise);
737
738           // subtract reference values
739           Double_t refped = 0;
740           Double_t refnoise = 0;
741           
742           if (fRefHistPed)   refped   = fRefHistPed->GetBinContent(bin);
743           if (fRefHistNoise) refnoise = fRefHistPed->GetBinContent(bin);
744
745           fChPedRes[det]->SetBinContent(bin, ped-refped);
746           fChNoiseRes[det]->SetBinContent(bin, noise-refnoise);
747           
748           fPed[det]->Fill(ped);
749           fNoise[det]->Fill(noise);
750
751           // fill SM-HC plot
752           Int_t sm = det / 30;
753           Int_t hc = (pad < kPAD/2) ? 2* (det % 30) : 2* (det % 30) + 1;
754           if (ped > 9. && ped < 11) fSMHCped->Fill(sm, hc, 1./1152.); // number of pads in HC
755
756         } else {
757           
758           // not enought data found 
759           fChPed[det]->SetBinContent(bin, 0);
760           fChNoise[det]->SetBinContent(bin, 0);
761           fChPedRes[det]->SetBinContent(bin, 0);
762           fChNoiseRes[det]->SetBinContent(bin, 0);
763         }
764         
765         //delete hist;
766       }
767     }
768   }
769
770
771   //AliInfo(Form("Number of events = %d\n", fnEvents));
772
773   // normalize number of entries histos
774   Int_t max = 0;
775   for(Int_t i=0; i<kDET; i++) { 
776     if (!map[i]) continue;
777     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {
778       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {
779         Int_t dataBin = fNPoint[i]->FindBin(j, k);
780         Double_t v = fNPoint[i]->GetBinContent(dataBin);
781         if (v > max) max = (Int_t)v;
782       }
783     }
784   }
785   
786   char entriesDistName[100];
787
788   for(Int_t i=0; i<kDET; i++) {
789     
790     if (!map[i]) continue;
791     
792     sprintf(entriesDistName, "entriesDist_%d", i);
793     fNPointDist[i] = new TH1D(entriesDistName, ";number of events", max+2, -0.5, max+1.5);
794     
795     for(Int_t j=0; j<fNPoint[i]->GetXaxis()->GetNbins(); j++) {
796       for(Int_t k=0; k<fNPoint[i]->GetYaxis()->GetNbins(); k++) {
797         Int_t dataBin = fNPoint[i]->FindBin(j, k);
798         Double_t v = fNPoint[i]->GetBinContent(dataBin);
799         //if (v > fnEvents) AliInfo(Form("N = %d V = %lf\n", fnEvents, v));
800         fNPointDist[i]->Fill(v); 
801       }
802     }
803     
804     fNPoint[i]->Scale(1./fnEvents);
805   }
806   
807
808   for(Int_t i=0; i<kDET; i++) {
809     fnEntriesRM[i]->SetMaximum(fnEvents * 1.5);
810   }
811
812   // save histograms
813
814   //AliInfo(Form("FILENAME 2 = %s (%d)\n", fn, fn));
815   TFile *file = new TFile(fn, "recreate");
816   for(Int_t det = 0; det < kDET; det++) {
817     if (!map[det]) continue; 
818     fChPed[det]->Write();
819     fChNoise[det]->Write();
820     fNPoint[det]->Write();
821     fNPointDist[det]->Write();
822     fPed[det]->Write();
823     fNoise[det]->Write();
824     fSignal[det]->Write();
825     fnEntriesRM[det]->Write();
826     fChPP[det]->Write();
827
828     fChPedRes[det]->Write();
829     fChNoiseRes[det]->Write();
830
831     // save error hists
832     fErrorLocMCM[det]->SetMinimum(0);
833     fErrorLocMCM[det]->SetMaximum(fnEvents);
834     fErrorLocMCM[det]->Write();
835
836     fErrorLocADC[det]->SetMinimum(0);
837     fErrorLocADC[det]->SetMaximum(fnEvents);
838     fErrorLocADC[det]->Write();
839   }
840
841   for(Int_t sm=0; sm<kSM; sm++) {
842     fSmNoiseRms[sm]->Write();
843     fSmNoiseFit[sm]->Write();
844     fSmPP[sm]->Write();
845   }
846
847
848
849   Int_t nMcm = 0;
850   for(Int_t i=0; i < kDET * kROB * kMCM; i++) {
851     if (fFullSignal[i] && fFullCounter[i] > fCount) {
852       fFullSignal[i]->Write();
853       nMcm++;
854     }
855   }
856
857   AliInfo(Form("Number of saved MCMs = %d\n", nMcm));
858   
859   fMcmTracks->Write();
860   AliInfo(Form("Number of tracks = %d\n", fMcmTracks->GetEntries()));
861   
862   // permanently problematic MCMs
863   for(Int_t det=0; det<kDET; det++) {
864     for(Int_t mcm=0; mcm<kROB*kMCM; mcm++) {
865       
866       Int_t mRob = mcm / kMCM;
867       Int_t mMcm = mcm % kMCM;
868       Int_t bin = fMapMCM->FindBin(det, mcm);
869       Double_t frac = 1. * fMapMCM->GetBinContent(bin) / fnEvents;      
870       fFracMCM->Fill(frac);
871       
872       if (frac > 0.7) {
873         AliInfo(Form("{%d, %d, %d, %f}, \n", det, mRob, mMcm, frac));
874       }      
875     }
876   }
877
878
879   fOccupancy->Write();
880   fDetRob->Write();
881   fTBEvent->Write();
882   
883   // error hists
884   fErrorHC->Write();
885   fErrorMCM->Write();
886   fErrorADC->Write();
887
888   fErrorSMHC->Write();
889   fErrorSMMCM->Write();
890   fErrorSMADC->Write();  
891   
892   // write graphs
893   fErrorGraphHC->Write("trendErrorHC");
894   fErrorGraphMCM->Write("trendErrorMCM");
895   fErrorGraphADC->Write("trendErrorADC");
896   
897   fGraphMCM->Write("trendMCM");
898   
899   for(Int_t i=0; i<3; i++) {
900     fGraphPP[i]->Write(Form("fracPP_%d", i));
901   }
902   
903   //fZSsize->Scale(1./fnEvents);
904   //fZSsize->Write();
905
906   fMapMCM->SetMaximum(fnEvents);
907   fMapMCM->Write();
908   fFracMCM->Write();
909   
910   fSMHCped->Write();
911
912   for(Int_t i=0; i<3; i++ ) {
913     fSMLink[i]->Write();
914     fGrLink[i]->Write();
915   }
916
917   for(Int_t sm=0; sm<kSM; sm++)
918     fNumberADC[sm]->Write(Form("nADCinSM%d",sm));
919   
920   fNumberADC[kSM]->Write("nADCinEvent");
921
922   fNoiseTotal->Write();
923   fPP->Write();
924
925   for(Int_t i=0; i<1000; i++) {
926     if (fEvNoDist[i]->GetSum() > 0) fEvNoDist[i]->Write();
927   }
928
929
930   file->Close();
931   delete file;
932 }
933
934 ///////////////////////////////////////////////////////////////////////////////////////////////////
935
936 Int_t AliTRDqaBlackEvents::CheckMCM(Int_t index) const {
937   //
938   // Checks a single MCM
939   //  
940
941   return 1;
942   
943   static Int_t data[21][3] = {
944     {1, 0, 1}, 
945     {242, 0, 0}, 
946     {242, 0, 1}, 
947     {242, 0, 2}, 
948     {242, 0, 4}, 
949     {242, 0, 5}, 
950     {242, 0, 6}, 
951     {242, 0, 8}, 
952     {242, 0, 12}, 
953     {251, 7, 7}, 
954     {254, 3, 11}, 
955     {259, 3, 14}, 
956     {260, 1, 9}, 
957     {260, 3, 15}, 
958     {273, 1, 7}, 
959     {273, 1, 15}, 
960     {276, 5, 11}, 
961     {280, 6, 2}, 
962     {299, 6, 4}, 
963     {511, 2, 9}, 
964     {517, 7, 15}
965   };
966   
967   for(Int_t i=0; i<21; i++) {
968     Int_t wIndex = data[i][0] * kROB*kMCM + data[i][1] * kMCM + data[i][2];
969     if (index == wIndex) return 0;
970   }
971
972   return 1;
973 }
974
975 ///////////////////////////////////////////////////////////////////////////////////////////////////
976
977
978
979
980
981
982
983 void AliTRDqaBlackEvents::DrawChamber(const char *filename, Int_t det, Int_t w, Int_t h) 
984 {
985   //
986   // Draw raport for one chamber: 
987   // pedestal map, noise map, distribution of pedestal and noise
988   // 
989   // input:
990   // name of the file with histograms (created with Process())
991   // detector Id (0 - 539)
992   // 
993
994   // setup global style
995   gStyle->SetPalette(1);
996   gStyle->SetOptStat(0);
997   gStyle->SetPadTopMargin(0.02);
998   gStyle->SetPadBottomMargin(0.05);
999
1000   TFile *file = new TFile(filename, "READ");
1001
1002   TCanvas *c = new TCanvas("blackEvents",Form("blackEvents %d",det), w, h);
1003   c->SetVertical(kFALSE);
1004   c->Divide(3,1, 0.01, 0.01);
1005   c->cd(3);
1006   
1007   TPad *mPad = (TPad*) gPad;
1008   mPad->Divide(1,2,0.01,0.01);
1009   
1010   c->cd(1);
1011   TH2D *h2 = (TH2D*)file->Get(Form("ped_%d",det));
1012   h2->SetMinimum(5);
1013   h2->SetMaximum(15);
1014   h2->SetTitle(";Z direction;#phi direction");
1015   h2->Draw("colz");
1016   
1017   c->cd(2);
1018   h2 = (TH2D*)file->Get(Form("noise_%d",det));
1019   h2->SetMinimum(fMinNoise);
1020   h2->SetMaximum(fMaxNoise);
1021   h2->SetTitle(";Z direction;#phi direction");
1022   h2->Draw("colz");
1023   
1024   mPad->cd(1);
1025   //gPad->SetLogy();
1026   TH1D *h1 = (TH1D*)file->Get(Form("pedDist_%d", det));
1027   h1->Draw();
1028   
1029   mPad->cd(2);
1030   gPad->SetLogy();
1031   h1 = (TH1D*)file->Get(Form("noiseDist_%d", det));
1032   h1->Draw();                    
1033   
1034   h1->Fit("gaus");
1035   TF1 *f = h1->GetFunction("gaus");
1036   const char *tt = Form("#mu = %.2f #sigma = %0.2f ", f->GetParameter(1),f->GetParameter(2));
1037   TLatex *ll = new TLatex(2, 100, tt);
1038   ll->SetTextSize(0.06);
1039   ll->Draw();
1040 }
1041
1042 ///////////////////////////////////////////////////////////////////////////////////////////////////
1043
1044 void AliTRDqaBlackEvents::DrawSm(const char *filename, Int_t sm, Int_t w, Int_t h) 
1045 {
1046   //
1047   // ????????????
1048   //
1049   
1050   gStyle->SetPalette(1);
1051   gStyle->SetOptStat(0);
1052   
1053   gStyle->SetPadTopMargin(0.02);
1054   //gStyle->SetPadBottomMargin(0.05);  
1055   //gStyle->SetPadLeftMargin(0.02);  
1056   //gStyle->SetPadRightMargin(0.02);
1057
1058   TFile *file = new TFile(filename, "READ");
1059
1060   TCanvas *c = new TCanvas("blackEventsSM",Form("blackEvents SM %d",sm), w, h);
1061   c->SetVertical(kFALSE);
1062   c->Divide(5, 6, 0.001, 0.01);
1063   
1064   for(Int_t i=0; i<30; i++) {
1065     
1066     TH2D *h2 = (TH2D*)file->Get(Form("noise_%d",i+30*sm));
1067     if (!h2) continue;
1068     h2->SetMinimum(fMinNoise);
1069     h2->SetMaximum(fMaxNoise);
1070
1071     // to be replaced by the official calculation
1072     Int_t stack = i/6;
1073     Int_t layer = i%6;
1074     Int_t index = (5-layer)*5 + stack + 1;
1075     //AliInfo(Form("%d %d %d %d\n", i, stack, layer, index));
1076     c->cd(index);
1077     gPad->SetBottomMargin(0.02);
1078     gPad->SetTopMargin(0.02);
1079
1080     h2->Draw("col");
1081   }
1082 }
1083
1084 ///////////////////////////////////////////////////////////////////////////////////////////////////
1085
1086 Int_t AliTRDqaBlackEvents::FillBits(TH1D *hist, Int_t code, Int_t offset) {
1087   //
1088   // Fill bits
1089   //
1090
1091   Int_t nb = 0;
1092   UInt_t test = 1;
1093   for(Int_t i=0; i<8; i++) {
1094     if (code & test) {
1095       hist->Fill(i+offset);
1096       nb++;
1097     }
1098     test *= 2;       
1099   }
1100   
1101   return nb;
1102 }
1103
1104 ///////////////////////////////////////////////////////////////////////////////////////////////////