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