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 #include "AliEMCALQADataMakerRec.h"
  50
  51 ClassImp(AliEMCALQAChecker)
  52
  53 //__________________________________________________________________
  54 AliEMCALQAChecker::AliEMCALQAChecker() :
  55 AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
  56 fLine(new TLine*[fknSM]),
  57 fHref(new TLine*[fknSM]),
  58 fText(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 //__________________________________________________________________
  68 AliEMCALQAChecker::~AliEMCALQAChecker()
  69 {
  70         /// dtor
  71   delete [] fLine ;
  72   delete [] fHref ;
  73   if (fText)
  74     delete fText ;
  75 }
  76
  77 //__________________________________________________________________
  78 AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) :
  79 AliQACheckerBase(qac.GetName(), qac.GetTitle()),
  80 fLine(new TLine*[fknSM]),
  81 fHref(new TLine*[fknSM]),
  82 fText(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 //__________________________________________________________________
  91 AliEMCALQAChecker& 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 //______________________________________________________________________________
 103 Double_t *
 104 AliEMCALQAChecker::Check(AliQAv1::ALITASK_t index, TObjArray ** list, const 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 //______________________________________________________________________________
 128 TH1*
 129 AliEMCALQAChecker::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 //______________________________________________________________________________
 141 Double_t
 142 AliEMCALQAChecker::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 //______________________________________________________________________________
 156 Double_t *
 157 AliEMCALQAChecker::CheckRaws(TObjArray ** list)
 158 {
 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.
 167
 168 //  -- Yaxian Mao, CCNU/CERN/LPSC
 169
 170   Float_t kThreshold = 80. ;
 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;
 183       Int_t  nbin   = hdata->GetNbinsX() ;
 184       Int_t  nTower = nbin/4 ;  // ! number of channels for each SM
 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));
 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
 205
 206        for(Int_t iTower = iSM*nTower ; iTower<(iSM+1)*nTower ; iTower++)  //channel loop to calculate the average
 207           aver += hdata->GetBinContent(iTower);
 208
 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]) {
 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         }
 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));
 237
 238        if(aver) init=TMath::Abs(mean/aver);  //calculate the divation from the average
 239
 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)
 242           flag = -1 ;
 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)){
 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
 266
 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
 270
 271         rv+=init ; //define the deviation from the final average
 272
 273         //define the average line on each SM
 274        fHref[iSM] = dynamic_cast<TLine*>(hdata->GetListOfFunctions()->FindObject(fHref[iSM]));
 275        if(!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) ;
 279        }
 280        else {
 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) ;
 286        }
 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() ;
 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() ;
 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       }
 326
 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       }
 331
 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));
 334         test[specie] =  0.7;
 335       }
 336
 337       if ( rv >0.8 ) {
 338         AliError(Form("Got too large deviation of %2.2f from average, detector out of control ???",rv));
 339         test[specie] =  0.9;
 340       }
 341
 342     } //end of checking raw
 343
 344   }  //species loop
 345   return test ;
 346 }
 347
 348 //______________________________________________________________________________
 349 void 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 ;
 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]) ;
 364 }
 365