]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EMCAL/AliEMCALQAChecker.cxx
This is a backward incompatible change in AliRoot. The following methods have been...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALQAChecker.cxx
... / ...
CommitLineData
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
50ClassImp(AliEMCALQAChecker)
51
52//__________________________________________________________________
53AliEMCALQAChecker::AliEMCALQAChecker() :
54AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
55fLine(new TLine*[fknSM]),
56fHref(new TLine*[fknSM]),
57fText(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//__________________________________________________________________
67AliEMCALQAChecker::~AliEMCALQAChecker()
68{
69 /// dtor
70 delete [] fLine ;
71 delete [] fHref ;
72 if (fText)
73 delete fText ;
74}
75
76//__________________________________________________________________
77AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) :
78AliQACheckerBase(qac.GetName(), qac.GetTitle()),
79fLine(new TLine*[fknSM]),
80fHref(new TLine*[fknSM]),
81fText(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//__________________________________________________________________
90AliEMCALQAChecker& 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//______________________________________________________________________________
102Double_t *
103AliEMCALQAChecker::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//______________________________________________________________________________
127TH1*
128AliEMCALQAChecker::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//______________________________________________________________________________
140Double_t
141AliEMCALQAChecker::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//______________________________________________________________________________
155Double_t *
156AliEMCALQAChecker::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 if(hdata->GetEntries()==0)
181 continue;
182 Int_t nbin = hdata->GetNbinsX() ;
183 Int_t nTower = nbin/4 ; // ! number of channels for each SM
184 if(fText) {
185 fText->DeleteText() ;
186 fText->Clear() ;
187 }
188 else {
189 fText = new TPaveText(0.3,0.6,0.7,0.9,"NDC") ;
190 }
191 fText->AddText(Form("OK if more than %2.2f %% inside aver-sigma < HG counts < aver+sigma", kThreshold));
192 //TPaveText * fText[fknSM] = {0};
193 Double_t rv = 0. ; //value to define the flag for the full detector
194 for(Int_t iSM = 0 ; iSM < fknSM ; iSM++){ //number of SMs loop start
195 TString projname = Form("proj_%d",iSM);
196 Double_t aver = 0. ; //average calculated by the input histogram
197 Double_t recal = 0 ; // recalculate the average if the divation is large
198 Double_t init = 0. ; // value to decide whether we should recalculate the average
199 Double_t mean = 0. ; //the divation from the average, from the gaus fitting peak
200 Double_t width = 0. ; // the sigma of the devation, from the gaus fitting
201 Int_t flag = 0 ; //counts for channels used for recalculation of average
202 Double_t ratio = 0. ; //value to decide whether each SM works
203 TH1F * proj = NULL ; //a temp histogram store the difference from the average for the fitting procedure
204
205 for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++) //channel loop to calculate the average
206 aver += hdata->GetBinContent(iTower);
207
208 aver /=nTower;
209 AliInfo(Form("SM: %d has average = %f\n",iSM, aver));
210 Double_t ymin = hdata->GetBinContent(hdata->GetMinimumBin());
211 Double_t ymax = hdata->GetBinContent(hdata->GetMaximumBin());
212 proj = new TH1F(projname,projname,nbin,-aver,aver);
213 fLine[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fLine[iSM]));
214 //initialize the lines separate different SM
215 if (!fLine[iSM]) {
216 fLine[iSM] = new TLine((iSM+1)*nTower,ymin,(iSM+1)*nTower,ymax);
217 fLine[iSM]->SetLineColor(1);
218 fLine[iSM]->SetLineWidth(2);
219 hdata->GetListOfFunctions()->Add(fLine[iSM]);
220 list[specie]->AddAt(hdata, kTowerHG) ;
221 }
222 else {
223 fLine[iSM]->SetX1((iSM+1)*nTower) ;
224 fLine[iSM]->SetY1(ymin) ;
225 fLine[iSM]->SetX2((iSM+1)*nTower) ;
226 fLine[iSM]->SetY2(ymax) ;
227 }
228 //filling the histogram with the difference and fitting with gaus function to obtain mean and sigma
229 for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++){
230 proj->Fill(hdata->GetBinContent(iTower)-aver);
231 }
232 proj->Fit("gaus","","QNO");
233 mean=proj->GetFunction("gaus")->GetParameter(1); // ! mean should be peaked at 0 in principal
234 width=proj->GetFunction("gaus")->GetParameter(2);
235 AliInfo(Form("aver = %f, mean = %f, sigma =%f\n",aver,mean, width));
236
237 if(aver) init=TMath::Abs(mean/aver); //calculate the divation from the average
238
239 //if mean or sigma is too huge or the sigma is too small, the fitting failed
240 if((aver+mean) < 1. || width/aver >2. || width < 1.e-3)
241 flag = -1 ;
242 else { //recalculate the average if the fitting didn't give the initial average
243 aver+=mean; //average corrected after fitting is fitting works
244 //if(aver <= 1. ) break ;
245 // if divation is too large, we conside to recalate the average by excluding channels too far from the average
246 while(init>0.01 && aver >= 1.){
247 if(flag) flag = 0 ; //for each time recalculation, reset flag
248 for(Int_t iTower = iSM*nTower ; iTower < (iSM+1)*nTower ; iTower++){
249 if(hdata->GetBinContent(iTower)>(aver-width) && hdata->GetBinContent(iTower)<(aver+width)){
250 flag++ ;
251 recal += hdata->GetBinContent(iTower);
252 }
253 } //end of channels loop
254
255 if(flag == 0 || recal < 1.) {
256 flag = -1 ;
257 break ;
258 }
259 else {
260 recal/=flag ;
261 init=(aver-recal)/(aver+recal) ; //difference between the new average and the pervious one
262 aver =(aver+recal)/2 ; //recalculate the average by the last two values
263 }
264 } //out of while condition
265
266 ratio=100.0*flag/nTower ; //counting how many towers used for average recalculation
267 AliInfo(Form("SM %d has %2.2f %% chanhel works \n", iSM, ratio));
268 } // end of recalculation
269
270 rv+=init ; //define the deviation from the final average
271
272 //define the average line on each SM
273 fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fHref[iSM]));
274 if(!fHref[iSM]) {
275 fHref[iSM] = new TLine(iSM*nTower,aver,(iSM+1)*nTower,aver);
276 hdata->GetListOfFunctions()->Add(fHref[iSM]);
277 list[specie]->AddAt(hdata, kTowerHG) ;
278 }
279 else {
280 fHref[iSM]->Clear() ;
281 fHref[iSM]->SetX1(iSM*nTower) ;
282 fHref[iSM]->SetY1(aver) ;
283 fHref[iSM]->SetX2((iSM+1)*nTower) ;
284 fHref[iSM]->SetY2(aver) ;
285 }
286
287 hdata->Paint() ;
288
289 // if channels used for average recalculation smaller than the threshold,
290 // then too much noise channels or channels didn't work
291 if(ratio<= kThreshold && flag >0){
292 //fText->SetFillColor(2);
293 fHref[iSM]->SetLineColor(2);
294 fHref[iSM]->SetLineWidth(2);
295 fText->SetFillColor(2);
296 fText->AddText(Form("SM %d: NOK = %2.0f %% channels OK!!!",iSM,ratio));
297 //fText[iSM]->AddText(Form("SM %d NOT OK, only %5.2f %% works!!!",iSM,flag/nTower));
298 }
299 else if (flag == -1 || flag == 0 ) { //fitting failed or recalculation didn't call
300 fText->SetFillColor(2);
301 fHref[iSM]->SetLineColor(2);
302 fHref[iSM]->SetLineWidth(2);
303 fText->SetFillColor(2);
304 fText->AddText(Form("SM %d: NOK Fitting failed",iSM,ratio));
305 //fText[iSM]->AddText(Form("SM %d has %5.2f %% towers wok normally",iSM,flag/nTower));
306 }
307 else {
308 fText->SetFillColor(3);
309 fHref[iSM]->SetLineColor(3);
310 fHref[iSM]->SetLineWidth(2);
311 fText->AddText(Form("SM %d: OK = %2.0f %% channels OK",iSM,ratio));
312 //fText[iSM]->AddText(Form("SM %d has %5.2f %% towers wok normally",iSM,flag/nTower));
313 }
314 hdata->Paint() ;
315 //hdata->GetListOfFunctions()->Add(fText[iSM]);
316 delete proj ;
317 } // end of SM loop
318 rv/=fknSM;
319 hdata->GetListOfFunctions()->Add(fText);
320 hdata->Paint() ;
321 if ( rv <=0.1 ) {
322 AliInfo(Form("The deviation rv = %2.2f is small compared to average, detector works fine",rv));
323 test[specie] = 0.05;
324 }
325
326 if ( rv <=0.5 && rv >0.1 ) {
327 AliWarning(Form("The deviation rv = %2.2f is acceptable from average, the detector works not perfect!",rv));
328 test[specie] = 0.3;
329 }
330
331 if ( rv <=0.8 && rv >0.5 ) {
332 AliError(Form("Got a large deviation of %2.2f from average, some error on the detector !!!",rv));
333 test[specie] = 0.7;
334 }
335
336 if ( rv >0.8 ) {
337 AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
338 test[specie] = 0.9;
339 }
340
341 } //end of checking raw
342
343 } //species loop
344 return test ;
345}
346
347//______________________________________________________________________________
348void AliEMCALQAChecker::Init(const AliQAv1::DETECTORINDEX_t det)
349{
350 /// intialises QA and QA checker settings
351 AliQAv1::Instance(det) ;
352 Float_t hiValue[AliQAv1::kNBIT] ;
353 Float_t lowValue[AliQAv1::kNBIT] ;
354 lowValue[AliQAv1::kINFO] = 0.0 ;
355 hiValue[AliQAv1::kINFO] = 0.1 ;
356 lowValue[AliQAv1::kWARNING] = 0.1 ;
357 hiValue[AliQAv1::kWARNING] = 0.5 ;
358 lowValue[AliQAv1::kERROR] = 0.5 ;
359 hiValue[AliQAv1::kERROR] = 0.8 ;
360 lowValue[AliQAv1::kFATAL] = 0.8 ;
361 hiValue[AliQAv1::kFATAL] = 1.0 ;
362 SetHiLo(&hiValue[0], &lowValue[0]) ;
363}
364