Fixing bug #57328
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALQAChecker.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  * 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 /*
17   Checks the quality assurance. 
18   By comparing with reference data
19
20   Based on PHOS code written by
21   Y. Schutz CERN July 2007
22
23  For the moment we only implement the checking of raw data QA.
24  The checked for ESD and RecPoints will be implemented later.
25  
26
27  */
28
29
30 // --- ROOT system ---
31 #include <TClass.h>
32 #include <TH1.h> 
33 #include <TF1.h> 
34 #include <TH1I.h> 
35 #include <TIterator.h> 
36 #include <TKey.h> 
37 #include <TFile.h> 
38 #include <TLine.h>
39 #include <TPaveText.h>
40 #include <TMath.h>
41
42 // --- Standard library ---
43
44 // --- AliRoot header files ---
45 #include "AliLog.h"
46 #include "AliQAv1.h"
47 #include "AliQAChecker.h"
48 #include "AliEMCALQAChecker.h"
49 #include "AliEMCALQADataMakerRec.h"
50
51 ClassImp(AliEMCALQAChecker)
52
53 //__________________________________________________________________
54 AliEMCALQAChecker::AliEMCALQAChecker() : 
55 AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
56 fLine(new TLine*[fknSM]),
57 fHref(new TLine*[fknSM]),
58 fText(NULL)
59 {
60         /// ctor
61   for (Int_t sm = 0 ; sm < fknSM ; sm++){
62     fLine[sm] = NULL ; 
63     fHref[sm] = NULL ; 
64   }
65 }          
66
67 //__________________________________________________________________
68 AliEMCALQAChecker::~AliEMCALQAChecker() 
69 {
70         /// dtor
71   delete [] fLine ; 
72   delete [] fHref ;
73   if (fText) 
74     delete fText ; 
75 }
76
77 //__________________________________________________________________
78 AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) : 
79 AliQACheckerBase(qac.GetName(), qac.GetTitle()), 
80 fLine(new TLine*[fknSM]),
81 fHref(new TLine*[fknSM]),
82 fText(qac.fText)
83 {
84         /// copy ctor 
85   for (Int_t sm = 0 ; sm < fknSM ; sm++){
86     fLine[sm] = qac.fLine[sm] ; 
87     fHref[sm] = qac.fHref[sm] ; 
88   }
89 }   
90 //__________________________________________________________________
91 AliEMCALQAChecker& AliEMCALQAChecker::operator = (const AliEMCALQAChecker &qac)
92 {
93   fText = qac.fText;
94
95   for (Int_t sm = 0 ; sm < fknSM ; sm++){
96     fLine[sm] = qac.fLine[sm] ; 
97     fHref[sm] = qac.fHref[sm] ; 
98   }
99   return *this ;
100 }
101
102 //______________________________________________________________________________
103 Double_t *
104 AliEMCALQAChecker::Check(AliQAv1::ALITASK_t index, TObjArray ** list, const AliDetectorRecoParam * /*recoParam*/)
105 {
106         /// Check objects in list
107         
108         if ( index == AliQAv1::kRAW ) 
109         {
110                 return CheckRaws(list);
111                 printf ("checkers for task %d \n", index) ;             
112         }
113         
114         if ( index == AliQAv1::kREC)
115         {
116                 return CheckRecPoints(list);
117         }
118         
119         if ( index == AliQAv1::kESD )
120         {
121                 return CheckESD(list);
122         }
123         AliWarning(Form("Checker for task %d not implement for the moment",index));
124         return NULL;
125 }
126
127 //______________________________________________________________________________
128 TH1* 
129 AliEMCALQAChecker::GetHisto(TObjArray* list, const char* hname, Int_t specie) const
130 {
131         /// Get a given histo from the list
132         TH1* h = static_cast<TH1*>(list->FindObject(Form("%s_%s",AliRecoParam::GetEventSpecieName(specie),hname)));
133         if (!h)
134         {
135                 AliError(Form("Did not find expected histo %s",hname));
136         }
137         return h;
138 }
139
140 //______________________________________________________________________________
141 Double_t 
142 AliEMCALQAChecker::MarkHisto(TH1& histo, Double_t value) const
143 {
144         /// Mark histo as originator of some QA error/warning
145         
146         if ( value != 1.0 )
147         {
148                 histo.SetBit(AliQAv1::GetQABit());
149         }
150         
151         return value;
152 }
153
154
155 //______________________________________________________________________________
156 Double_t *
157 AliEMCALQAChecker::CheckRaws(TObjArray ** list)
158 {
159 //  Check RAW QA histograms     
160 //  We count the times of the response for each tower, the propability for all towers should be the same (average is given value).
161 //  We skip the first few cycles since the statistics is not enough, the average should be always larger than 1 at least.
162 //  By calculating the difference between the counts for each tower and the average, we decide whether we should recalculate 
163 //  the average depending the the gaus fitting on the divation distribution. 
164 //  During the recalutation of the average, we count how many towers are used for the calculation.
165 //  From the fraction of towers used, we decide whether each SM works fine or not
166 //  From the divation of average, we set the QA flag for the full detetcor as INFO, WARNING, ERROR or FATAL.
167   
168 //  -- Yaxian Mao, CCNU/CERN/LPSC
169                                         
170   Float_t kThreshold = 80. ; 
171   
172   Double_t * test = new Double_t[AliRecoParam::kNSpecies] ;
173   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
174     test[specie] = 1.0 ; 
175     if ( !AliQAv1::Instance()->IsEventSpecieSet(specie)) 
176       continue ; 
177     if (list[specie]->GetEntries() == 0)  
178       test[specie] = 0. ; // nothing to check
179     else {
180       TH1 * hdata = (TH1*)list[specie]->At(kTowerHG) ; 
181       if(hdata->GetEntries()==0)
182         continue;
183       Int_t  nbin   = hdata->GetNbinsX() ;
184       Int_t  nTower = nbin/4 ;  // ! number of channels for each SM
185       if(fText) {
186         fText->DeleteText() ; 
187         fText->Clear() ;
188       }
189       else {
190         fText = new TPaveText(0.3,0.6,0.7,0.9,"NDC") ;
191       }
192       fText->AddText(Form("OK if more than %2.2f %% inside aver-sigma < HG counts < aver+sigma", kThreshold));
193       //TPaveText * fText[fknSM] = {0};
194       Double_t rv    = 0. ;   //value to define the flag for the full detector
195       for(Int_t iSM = 0 ; iSM < fknSM ; iSM++){  //number of SMs loop start
196         TString projname = Form("proj_%d",iSM);
197         Double_t aver  = 0. ;  //average calculated by the input histogram
198         Double_t recal = 0 ;   // recalculate the average if the divation is large
199         Double_t init  = 0. ;   // value to decide whether we should recalculate the average
200         Double_t mean  = 0. ;  //the divation from the average, from the gaus fitting peak 
201         Double_t width = 0. ;   // the sigma of the devation, from the gaus fitting
202         Int_t flag     = 0 ;  //counts for channels used for recalculation of average
203         Double_t ratio = 0. ;  //value to decide whether each SM works
204         TH1F * proj    = NULL  ;  //a temp histogram store the difference from the average for the fitting procedure       
205        
206        for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++)  //channel loop to calculate the average
207           aver += hdata->GetBinContent(iTower);
208         
209         aver /=nTower;
210         AliInfo(Form("SM: %d has average = %f\n",iSM, aver));
211        Double_t ymin = hdata->GetBinContent(hdata->GetMinimumBin());
212        Double_t ymax = hdata->GetBinContent(hdata->GetMaximumBin());
213         proj = new TH1F(projname,projname,nbin,-aver,aver);
214        fLine[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fLine[iSM]));
215         //initialize the lines separate different SM
216        if (!fLine[iSM]) {
217           fLine[iSM] = new TLine((iSM+1)*nTower,ymin,(iSM+1)*nTower,ymax);
218           fLine[iSM]->SetLineColor(1);
219           fLine[iSM]->SetLineWidth(2);
220           hdata->GetListOfFunctions()->Add(fLine[iSM]);  
221           list[specie]->AddAt(hdata, kTowerHG) ; 
222         }          
223         else {
224           fLine[iSM]->SetX1((iSM+1)*nTower) ; 
225           fLine[iSM]->SetY1(ymin) ; 
226           fLine[iSM]->SetX2((iSM+1)*nTower) ; 
227           fLine[iSM]->SetY2(ymax) ; 
228         }
229         //filling the histogram with the difference and fitting with gaus function to obtain mean and sigma
230        for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++){
231          proj->Fill(hdata->GetBinContent(iTower)-aver);
232        }        
233        proj->Fit("gaus","","QNO");
234        mean=proj->GetFunction("gaus")->GetParameter(1); // ! mean should be peaked at 0 in principal
235        width=proj->GetFunction("gaus")->GetParameter(2);
236        AliInfo(Form("aver = %f, mean = %f, sigma =%f\n",aver,mean, width));
237         
238        if(aver) init=TMath::Abs(mean/aver);  //calculate the divation from the average 
239        
240        //if mean or sigma is too huge or the sigma is too small, the fitting failed 
241         if((aver+mean) < 1. || width/aver >2. || width < 1.e-3) 
242           flag = -1 ;
243        else {   //recalculate the average if the fitting didn't give the initial average  
244                aver+=mean;  //average corrected after fitting is fitting works
245              //if(aver <= 1. ) break ;
246             // if divation is too large, we conside to recalate the average by excluding channels too far from the average
247          while(init>0.01 && aver >= 1.){
248           if(flag) flag = 0 ;    //for each time recalculation, reset flag 
249            for(Int_t iTower = iSM*nTower ; iTower < (iSM+1)*nTower ; iTower++){
250              if(hdata->GetBinContent(iTower)>(aver-width) && hdata->GetBinContent(iTower)<(aver+width)){
251                flag++ ;
252                      recal += hdata->GetBinContent(iTower);
253              }  
254            }  //end of channels loop 
255            
256            if(flag == 0 || recal < 1.) {
257              flag = -1 ;
258              break ;
259            }
260            else {
261              recal/=flag ;
262                    init=(aver-recal)/(aver+recal) ; //difference between the new average and the pervious one
263             aver =(aver+recal)/2 ; //recalculate the average by the last two values  
264            }  
265          } //out of while condition
266           
267          ratio=100.0*flag/nTower ;  //counting how many towers used for average recalculation
268                AliInfo(Form("SM %d has %2.2f %% chanhel works \n", iSM, ratio));
269        }  // end of recalculation
270         
271         rv+=init ; //define the deviation from the final average
272     
273         //define the average line on each SM
274        fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fHref[iSM]));
275        if(!fHref[iSM]) {
276          fHref[iSM] = new TLine(iSM*nTower,aver,(iSM+1)*nTower,aver);
277                hdata->GetListOfFunctions()->Add(fHref[iSM]);  
278                list[specie]->AddAt(hdata, kTowerHG) ; 
279        }          
280        else {
281          fHref[iSM]->Clear() ; 
282                fHref[iSM]->SetX1(iSM*nTower) ; 
283          fHref[iSM]->SetY1(aver) ; 
284          fHref[iSM]->SetX2((iSM+1)*nTower) ; 
285          fHref[iSM]->SetY2(aver) ; 
286        }
287         
288         hdata->Paint() ; 
289  
290         // if channels used for average recalculation smaller than the threshold,
291         // then too much noise channels or channels didn't work 
292         if(ratio<= kThreshold && flag >0){                             
293           //fText->SetFillColor(2);
294           fHref[iSM]->SetLineColor(2);
295                 fHref[iSM]->SetLineWidth(2);
296                 fText->SetFillColor(2);
297                 fText->AddText(Form("SM %d: NOK = %2.0f %% channels OK!!!",iSM,ratio));
298                 //fText[iSM]->AddText(Form("SM %d NOT OK, only %5.2f %% works!!!",iSM,flag/nTower));
299         }         
300         else if (flag == -1 || flag == 0 ) {  //fitting failed or recalculation didn't call
301           fText->SetFillColor(2);
302           fHref[iSM]->SetLineColor(2);
303           fHref[iSM]->SetLineWidth(2);
304           fText->SetFillColor(2);
305           fText->AddText(Form("SM %d: NOK Fitting failed",iSM,ratio));
306                 //fText[iSM]->AddText(Form("SM %d has %5.2f %% towers wok normally",iSM,flag/nTower));
307         }  
308         else {
309           fText->SetFillColor(3);
310                 fHref[iSM]->SetLineColor(3);
311                 fHref[iSM]->SetLineWidth(2);
312                 fText->AddText(Form("SM %d: OK = %2.0f %% channels OK",iSM,ratio));
313                 //fText[iSM]->AddText(Form("SM %d has %5.2f %% towers wok normally",iSM,flag/nTower));
314         }
315         hdata->Paint() ; 
316        //hdata->GetListOfFunctions()->Add(fText[iSM]);
317        delete proj ; 
318       }  // end of SM loop
319       rv/=fknSM;
320       hdata->GetListOfFunctions()->Add(fText);
321       hdata->Paint() ;  
322       if ( rv <=0.1 ) {
323         AliInfo(Form("The deviation rv = %2.2f is small compared to average, detector works fine",rv));
324         test[specie] =  0.05;
325       }
326       
327       if ( rv <=0.5 && rv >0.1 )  {
328         AliWarning(Form("The deviation rv = %2.2f is acceptable from average,  the detector works not perfect!",rv));
329         test[specie] =  0.3;
330       }
331       
332       if ( rv <=0.8 && rv >0.5 ) {
333         AliError(Form("Got a large deviation of %2.2f from average, some error on the detector !!!",rv));
334         test[specie] =  0.7;
335       }
336       
337       if ( rv >0.8 ) {
338         AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
339         test[specie] =  0.9; 
340       }
341     
342     } //end of checking raw
343     
344   }  //species loop 
345   return test ;
346 }
347
348 //______________________________________________________________________________
349 void AliEMCALQAChecker::Init(const AliQAv1::DETECTORINDEX_t det) 
350 {
351         /// intialises QA and QA checker settings
352         AliQAv1::Instance(det) ; 
353         Float_t hiValue[AliQAv1::kNBIT] ; 
354         Float_t lowValue[AliQAv1::kNBIT] ;
355         lowValue[AliQAv1::kINFO]      = 0.0   ; 
356         hiValue[AliQAv1::kINFO]       = 0.1 ; 
357         lowValue[AliQAv1::kWARNING]   = 0.1 ; 
358   hiValue[AliQAv1::kWARNING]    = 0.5 ; 
359         lowValue[AliQAv1::kERROR]     = 0.5   ; 
360         hiValue[AliQAv1::kERROR]      = 0.8 ; 
361         lowValue[AliQAv1::kFATAL]     = 0.8   ; 
362         hiValue[AliQAv1::kFATAL]      = 1.0 ; 
363         SetHiLo(&hiValue[0], &lowValue[0]) ; 
364 }
365