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