1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 Checks the quality assurance.
18 By comparing with reference data
20 Based on PHOS code written by
21 Y. Schutz CERN July 2007
23 For the moment we only implement the checking of raw data QA.
24 The checked for ESD and RecPoints will be implemented later.
30 // --- ROOT system ---
35 #include <TIterator.h>
39 #include <TPaveText.h>
42 // --- Standard library ---
44 // --- AliRoot header files ---
47 #include "AliQAChecker.h"
48 #include "AliEMCALQAChecker.h"
49 #include "AliEMCALQADataMakerRec.h"
51 ClassImp(AliEMCALQAChecker)
53 //__________________________________________________________________
54 AliEMCALQAChecker::AliEMCALQAChecker() :
55 AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
56 fLine(new TLine*[fknSM]),
57 fHref(new TLine*[fknSM]),
61 for (Int_t sm = 0 ; sm < fknSM ; sm++){
67 //__________________________________________________________________
68 AliEMCALQAChecker::~AliEMCALQAChecker()
77 //__________________________________________________________________
78 AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) :
79 AliQACheckerBase(qac.GetName(), qac.GetTitle()),
80 fLine(new TLine*[fknSM]),
81 fHref(new TLine*[fknSM]),
85 for (Int_t sm = 0 ; sm < fknSM ; sm++){
86 fLine[sm] = qac.fLine[sm] ;
87 fHref[sm] = qac.fHref[sm] ;
90 //__________________________________________________________________
91 AliEMCALQAChecker& AliEMCALQAChecker::operator = (const AliEMCALQAChecker &qac)
95 for (Int_t sm = 0 ; sm < fknSM ; sm++){
96 fLine[sm] = qac.fLine[sm] ;
97 fHref[sm] = qac.fHref[sm] ;
102 //______________________________________________________________________________
104 AliEMCALQAChecker::Check(AliQAv1::ALITASK_t index, TObjArray ** list, const AliDetectorRecoParam * /*recoParam*/)
106 /// Check objects in list
108 if ( index == AliQAv1::kRAW )
110 return CheckRaws(list);
111 printf ("checkers for task %d \n", index) ;
114 if ( index == AliQAv1::kREC)
116 return CheckRecPoints(list);
119 if ( index == AliQAv1::kESD )
121 return CheckESD(list);
123 AliWarning(Form("Checker for task %d not implement for the moment",index));
127 //______________________________________________________________________________
129 AliEMCALQAChecker::GetHisto(TObjArray* list, const char* hname, Int_t specie) const
131 /// Get a given histo from the list
132 TH1* h = static_cast<TH1*>(list->FindObject(Form("%s_%s",AliRecoParam::GetEventSpecieName(specie),hname)));
135 AliError(Form("Did not find expected histo %s",hname));
140 //______________________________________________________________________________
142 AliEMCALQAChecker::MarkHisto(TH1& histo, Double_t value) const
144 /// Mark histo as originator of some QA error/warning
148 histo.SetBit(AliQAv1::GetQABit());
155 //______________________________________________________________________________
157 AliEMCALQAChecker::CheckRaws(TObjArray ** list)
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.
168 // -- Yaxian Mao, CCNU/CERN/LPSC
170 Float_t kThreshold = 80. ;
172 Double_t * test = new Double_t[AliRecoParam::kNSpecies] ;
173 for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
175 if ( !AliQAv1::Instance()->IsEventSpecieSet(specie))
177 if (list[specie]->GetEntries() == 0)
178 test[specie] = 0. ; // nothing to check
180 TH1 * hdata = (TH1*)list[specie]->At(kTowerHG) ;
181 if(hdata->GetEntries()==0)
183 Int_t nbin = hdata->GetNbinsX() ;
184 Int_t nTower = nbin/4 ; // ! number of channels for each SM
186 fText->DeleteText() ;
190 fText = new TPaveText(0.3,0.6,0.7,0.9,"NDC") ;
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
206 for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++) //channel loop to calculate the average
207 aver += hdata->GetBinContent(iTower);
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
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) ;
224 fLine[iSM]->SetX1((iSM+1)*nTower) ;
225 fLine[iSM]->SetY1(ymin) ;
226 fLine[iSM]->SetX2((iSM+1)*nTower) ;
227 fLine[iSM]->SetY2(ymax) ;
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);
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));
238 if(aver) init=TMath::Abs(mean/aver); //calculate the divation from the average
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)
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)){
252 recal += hdata->GetBinContent(iTower);
254 } //end of channels loop
256 if(flag == 0 || recal < 1.) {
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
265 } //out of while condition
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
271 rv+=init ; //define the deviation from the final average
273 //define the average line on each SM
274 fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(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) ;
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) ;
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));
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));
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));
316 //hdata->GetListOfFunctions()->Add(fText[iSM]);
320 hdata->GetListOfFunctions()->Add(fText);
323 AliInfo(Form("The deviation rv = %2.2f is small compared to average, detector works fine",rv));
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));
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));
338 AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
342 } //end of checking raw
348 //______________________________________________________________________________
349 void AliEMCALQAChecker::Init(const AliQAv1::DETECTORINDEX_t det)
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]) ;