Updated version of ITS QA Checker and related modifications (Melinda)
[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"
9e47432c 49#include "AliEMCALQADataMakerRec.h"
94594e5d 50
51ClassImp(AliEMCALQAChecker)
52
9e47432c 53//__________________________________________________________________
54AliEMCALQAChecker::AliEMCALQAChecker() :
55AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
56fLine(new TLine*[fknSM]),
57fHref(new TLine*[fknSM]),
58fText(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//__________________________________________________________________
68AliEMCALQAChecker::~AliEMCALQAChecker()
69{
70 /// dtor
71 delete [] fLine ;
72 delete [] fHref ;
73 if (fText)
74 delete fText ;
75}
76
77//__________________________________________________________________
78AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) :
79AliQACheckerBase(qac.GetName(), qac.GetTitle()),
80fLine(new TLine*[fknSM]),
81fHref(new TLine*[fknSM]),
82fText(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//__________________________________________________________________
91AliEMCALQAChecker& 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//______________________________________________________________________________
103Double_t *
104AliEMCALQAChecker::Check(AliQAv1::ALITASK_t index, TObjArray ** list, 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//______________________________________________________________________________
128TH1*
129AliEMCALQAChecker::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//______________________________________________________________________________
141Double_t
142AliEMCALQAChecker::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//______________________________________________________________________________
156Double_t *
157AliEMCALQAChecker::CheckRaws(TObjArray ** list)
158{
a2655076 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.
6a754398 167
a2655076 168// -- Yaxian Mao, CCNU/CERN/LPSC
169
9e47432c 170 Float_t kThreshold = 80. ;
a2655076 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;
9e47432c 183 Int_t nbin = hdata->GetNbinsX() ;
a2655076 184 Int_t nTower = nbin/4 ; // ! number of channels for each SM
9e47432c 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));
a2655076 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
9e47432c 196 TString projname = Form("proj_%d",iSM);
a2655076 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);
9e47432c 208
a2655076 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]) {
9e47432c 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 }
a2655076 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));
6a754398 237
238 if(aver) init=TMath::Abs(mean/aver); //calculate the divation from the average
a2655076 239
240 //if mean or sigma is too huge or the sigma is too small, the fitting failed
6a754398 241 if((aver+mean) < 1. || width/aver >2. || width < 1.e-3)
242 flag = -1 ;
a2655076 243 else { //recalculate the average if the fitting didn't give the initial average
6a754398 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
a2655076 266
6a754398 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));
a2655076 269 } // end of recalculation
6a754398 270
271 rv+=init ; //define the deviation from the final average
272
273 //define the average line on each SM
a2655076 274 fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fHref[iSM]));
275 if(!fHref[iSM]) {
6a754398 276 fHref[iSM] = new TLine(iSM*nTower,aver,(iSM+1)*nTower,aver);
277 hdata->GetListOfFunctions()->Add(fHref[iSM]);
278 list[specie]->AddAt(hdata, kTowerHG) ;
a2655076 279 }
280 else {
6a754398 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) ;
a2655076 286 }
6a754398 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() ;
a2655076 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() ;
6a754398 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 }
a2655076 326
6a754398 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 }
a2655076 331
6a754398 332 if ( rv <=0.8 && rv >0.5 ) {
8976d7bd 333 AliError(Form("Got a large deviation of %2.2f from average, some error on the detector !!!",rv));
6a754398 334 test[specie] = 0.7;
335 }
336
337 if ( rv >0.8 ) {
8976d7bd 338 AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
6a754398 339 test[specie] = 0.9;
340 }
341
a2655076 342 } //end of checking raw
343
344 } //species loop
345 return test ;
9e47432c 346}
347
348//______________________________________________________________________________
349void 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 ;
9e47432c 357 lowValue[AliQAv1::kWARNING] = 0.1 ;
6a754398 358 hiValue[AliQAv1::kWARNING] = 0.5 ;
9e47432c 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}
94594e5d 365