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
  50 ClassImp(AliEMCALQAChecker)
  51
  52 //__________________________________________________________________
  53 AliEMCALQAChecker::AliEMCALQAChecker() :
  54 AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
  55 fLine(new TLine*[fknSM]),
  56 fHref(new TLine*[fknSM]),
  57 fText(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 //__________________________________________________________________
  67 AliEMCALQAChecker::~AliEMCALQAChecker()
  68 {
  69         /// dtor
  70   delete [] fLine ;
  71   delete [] fHref ;
  72   if (fText)
  73     delete fText ;
  74 }
  75
  76 //__________________________________________________________________
  77 AliEMCALQAChecker::AliEMCALQAChecker(const AliEMCALQAChecker& qac) :
  78 AliQACheckerBase(qac.GetName(), qac.GetTitle()),
  79 fLine(new TLine*[fknSM]),
  80 fHref(new TLine*[fknSM]),
  81 fText(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 //__________________________________________________________________
  90 AliEMCALQAChecker& 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 //______________________________________________________________________________
 102 Double_t *
 103 AliEMCALQAChecker::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 //______________________________________________________________________________
 127 TH1*
 128 AliEMCALQAChecker::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 //______________________________________________________________________________
 140 Double_t
 141 AliEMCALQAChecker::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 //______________________________________________________________________________
 155 Double_t *
 156 AliEMCALQAChecker::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 //______________________________________________________________________________
 348 void 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