Correction suggested by Valgrind (Yves)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALQAChecker.cxx
CommitLineData
94594e5d 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
9e47432c 22
a2655076 23 For the moment we only implement the checking of raw data QA.
24 The checked for ESD and RecPoints will be implemented later.
9e47432c 25
26
9e47432c 27 */
28
94594e5d 29
30// --- ROOT system ---
31#include <TClass.h>
9e47432c 32#include <TH1.h>
33#include <TF1.h>
94594e5d 34#include <TH1I.h>
35#include <TIterator.h>
36#include <TKey.h>
37#include <TFile.h>
9e47432c 38#include <TLine.h>
39#include <TPaveText.h>
40#include <TMath.h>
94594e5d 41
42// --- Standard library ---
43
44// --- AliRoot header files ---
45#include "AliLog.h"
4e25ac79 46#include "AliQAv1.h"
94594e5d 47#include "AliQAChecker.h"
48#include "AliEMCALQAChecker.h"
49
50ClassImp(AliEMCALQAChecker)
51
9e47432c 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 *
486788fc 103AliEMCALQAChecker::Check(AliQAv1::ALITASK_t index, TObjArray ** list, const AliDetectorRecoParam * /*recoParam*/)
9e47432c 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{
a2655076 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.
6a754398 166
a2655076 167// -- Yaxian Mao, CCNU/CERN/LPSC
168
9e47432c 169 Float_t kThreshold = 80. ;
a2655076 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) ;
ffeb5d9c 180 hdata->GetListOfFunctions()->RemoveAll();
a2655076 181 if(hdata->GetEntries()==0)
ffeb5d9c 182 continue;
183 hdata->GetListOfFunctions()->Remove(fText);
9e47432c 184 Int_t nbin = hdata->GetNbinsX() ;
a2655076 185 Int_t nTower = nbin/4 ; // ! number of channels for each SM
9e47432c 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));
a2655076 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
9e47432c 197 TString projname = Form("proj_%d",iSM);
a2655076 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
ffeb5d9c 207 for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++) //channel loop to calculate the average
a2655076 208 aver += hdata->GetBinContent(iTower);
9e47432c 209
a2655076 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]) {
9e47432c 218 fLine[iSM] = new TLine((iSM+1)*nTower,ymin,(iSM+1)*nTower,ymax);
219 fLine[iSM]->SetLineColor(1);
220 fLine[iSM]->SetLineWidth(2);
ffeb5d9c 221 hdata->GetListOfFunctions()->Remove(fLine[iSM]);
9e47432c 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 }
a2655076 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));
6a754398 239
240 if(aver) init=TMath::Abs(mean/aver); //calculate the divation from the average
a2655076 241
242 //if mean or sigma is too huge or the sigma is too small, the fitting failed
6a754398 243 if((aver+mean) < 1. || width/aver >2. || width < 1.e-3)
244 flag = -1 ;
a2655076 245 else { //recalculate the average if the fitting didn't give the initial average
6a754398 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
a2655076 268
6a754398 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));
a2655076 271 } // end of recalculation
6a754398 272
273 rv+=init ; //define the deviation from the final average
274
275 //define the average line on each SM
a2655076 276 fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fHref[iSM]));
277 if(!fHref[iSM]) {
6a754398 278 fHref[iSM] = new TLine(iSM*nTower,aver,(iSM+1)*nTower,aver);
ffeb5d9c 279 hdata->GetListOfFunctions()->Remove(fHref[iSM]);
280 hdata->GetListOfFunctions()->Add(fHref[iSM]);
281 list[specie]->AddAt(hdata, kTowerHG) ;
a2655076 282 }
283 else {
6a754398 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) ;
a2655076 289 }
6a754398 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() ;
a2655076 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() ;
6a754398 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 }
a2655076 329
6a754398 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 }
a2655076 334
6a754398 335 if ( rv <=0.8 && rv >0.5 ) {
8976d7bd 336 AliError(Form("Got a large deviation of %2.2f from average, some error on the detector !!!",rv));
6a754398 337 test[specie] = 0.7;
338 }
339
340 if ( rv >0.8 ) {
8976d7bd 341 AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
6a754398 342 test[specie] = 0.9;
343 }
344
a2655076 345 } //end of checking raw
346
347 } //species loop
348 return test ;
9e47432c 349}
350
351//______________________________________________________________________________
352void 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 ;
9e47432c 360 lowValue[AliQAv1::kWARNING] = 0.1 ;
6a754398 361 hiValue[AliQAv1::kWARNING] = 0.5 ;
9e47432c 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}
94594e5d 368