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 <TH2F.h>
  36 #include <TIterator.h>
  37 #include <TKey.h>
  38 #include <TFile.h>
  39 #include <TLine.h>
  40 #include <TText.h>
  41 #include <TPaveText.h>
  42 #include <TMath.h>
  43
  44 // --- Standard library ---
  45
  46 // --- AliRoot header files ---
  47 #include "AliLog.h"
  48 #include "AliCDBManager.h"
  49 #include "AliCDBEntry.h"
  50 #include "AliQAv1.h"
  51 #include "AliQAChecker.h"
  52 #include "AliQAThresholds.h"
  53 #include "AliEMCALQAChecker.h"
  54
  55 ClassImp(AliEMCALQAChecker)
  56
  57 //__________________________________________________________________
  58 AliEMCALQAChecker::AliEMCALQAChecker() :
  59 AliQACheckerBase("EMCAL","EMCAL Quality Assurance Data Maker"),
  60 fTextSM(new TText*[fgknSM]),
  61 fLineCol(new TLine(47.5,-0.5,47.5,119.5)),
  62 fText(new TPaveText(0.2,0.7,0.8,0.9,"NDC"))
  63 {
  64   // ctor
  65   fLineCol->SetLineColor(1);
  66   fLineCol->SetLineWidth(2);
  67
  68   fTextSM[0]= new TText(20, 12, "SM A0");
  69   fTextSM[1]= new TText(20, 36, "SM A1");
  70   fTextSM[2]= new TText(20, 60, "SM A2");
  71   fTextSM[3]= new TText(20, 84, "SM A3");
  72   fTextSM[4]= new TText(20, 108,"SM A4");
  73   fTextSM[5]= new TText(20, 132,"SM A5");
  74
  75   fTextSM[6]= new TText(64, 12, "SM C0");
  76   fTextSM[7]= new TText(64, 36, "SM C1");
  77   fTextSM[8]= new TText(64, 60, "SM C2");
  78   fTextSM[9]= new TText(64, 84, "SM C3");
  79   fTextSM[10]= new TText(64, 108,"SM C4");
  80   fTextSM[11]= new TText(64, 132,"SM C5");
  81
  82         for(int i = 0; i < 5; i++) {
  83                 fLineRow[i] = new TLine(-0.5,23.5+(24*i),95.5,23.5+(24*i));
  84                 fLineRow[i]->SetLineColor(1);
  85                 fLineRow[i]->SetLineWidth(2);
  86         }
  87
  88         for(int i = 0; i < 3; i++) {
  89                 fTextL1[i] = new TPaveText(0.2,0.8,0.8,0.9,"NDC");
  90         }
  91
  92
  93 }
  94
  95 //__________________________________________________________________
  96 AliEMCALQAChecker::~AliEMCALQAChecker()
  97 {
  98         /// dtor
  99   delete [] fTextSM ;
 100   delete fLineCol ;
 101   for (Int_t i=0; i<5; ++i) delete fLineRow[i] ;
 102   delete fText  ;
 103 }
 104
 105 //______________________________________________________________________________
 106 void AliEMCALQAChecker::Check(Double_t * test, AliQAv1::ALITASK_t index, TObjArray ** list, const AliDetectorRecoParam * /*recoParam*/)
 107 {
 108         /// Check objects in list
 109
 110         if ( index == AliQAv1::kRAW )
 111         {
 112     CheckRaws(test, list);
 113                 printf ("checkers for task %d \n", index) ;
 114         }
 115
 116         if ( index == AliQAv1::kREC)
 117         {
 118     CheckRecPoints(test, list);
 119         }
 120
 121         if ( index == AliQAv1::kESD )
 122         {
 123     CheckESD(test, list);
 124         }
 125         AliWarning(Form("Checker for task %d not implement for the moment",index));
 126 }
 127
 128 //______________________________________________________________________________
 129 TH1*
 130 AliEMCALQAChecker::GetHisto(TObjArray* list, const char* hname, Int_t specie) const
 131 {
 132         /// Get a given histo from the list
 133         TH1* h = static_cast<TH1*>(list->FindObject(Form("%s_%s",AliRecoParam::GetEventSpecieName(specie),hname)));
 134         if (!h)
 135         {
 136                 AliError(Form("Did not find expected histo %s",hname));
 137         }
 138         return h;
 139 }
 140
 141 //______________________________________________________________________________
 142 Double_t
 143 AliEMCALQAChecker::MarkHisto(TH1& histo, Double_t value) const
 144 {
 145         /// Mark histo as originator of some QA error/warning
 146
 147         if ( value != 1.0 )
 148         {
 149                 histo.SetBit(AliQAv1::GetQABit());
 150         }
 151
 152         return value;
 153 }
 154
 155
 156 //______________________________________________________________________________
 157 void AliEMCALQAChecker::CheckRaws(Double_t * test, TObjArray ** list)
 158 {
 159   //  Check RAW QA histograms
 160   //  -- Yaxian Mao, CCNU/CERN/LPSC
 161   //adding new checking method: 25/04/2010, Yaxian Mao
 162   //Comparing the amplitude from current run to the reference run, if the ratio in the range [0.8, .12], count as a good tower.
 163   //If more than 90% towers are good, EMCAL works fine, otherwise experts should be contacted.
 164
 165
 166   // Retrieve the thresholds
 167   Float_t ratioThresh = 0.9;   // threshold for calibration ratio = good towers/all towers (default 0.9)
 168   Float_t ThreshG     = 0.5;   // threshold for L1 Gamma triggers (default 0.5)
 169   Float_t ThreshJ     = 0.5;   // threshold for L1 Jet triggers (default 0.5)
 170   Int_t badLinkThresh = 1;     // threshold for bad links (default 1)
 171
 172   AliCDBManager* man = AliCDBManager::Instance();
 173   if(man){
 174     AliCDBEntry* entry = man->Get("GRP/Calib/QAThresholds");
 175     if(entry){
 176       TObjArray* branch = (TObjArray*) entry->GetObject();
 177       if(branch){
 178         AliQAThresholds* thresholds = (AliQAThresholds*) branch->FindObject("EMC");
 179         if(thresholds){
 180           TParameter<float>* ParamR  = (TParameter<float>*) thresholds->GetThreshold(0);
 181           TParameter<float>* ParamG  = (TParameter<float>*) thresholds->GetThreshold(1);
 182           TParameter<float>* ParamJ  = (TParameter<float>*) thresholds->GetThreshold(2);
 183           TParameter<int>* ParamL  = (TParameter<int>*) thresholds->GetThreshold(3);
 184           if(ParamR)
 185             ratioThresh = ParamR->GetVal();
 186           if(ParamG)
 187             ThreshG = ParamG->GetVal();
 188           if(ParamJ)
 189             ThreshJ = ParamJ->GetVal();
 190           if(ParamL)
 191             badLinkThresh = ParamL->GetVal();
 192         }
 193       }
 194     }
 195   }
 196
 197   Int_t nTowersPerSM = 24*48; // number of towers in a SuperModule; 24x48
 198   Double_t nTot = fgknSM * nTowersPerSM ;
 199   TList *lstF = 0;
 200   Int_t calibSpecieId = (Int_t)TMath::Log2( AliRecoParam::kCalib );
 201   for (Int_t specie = 0 ; specie < AliRecoParam::kNSpecies ; specie++) {
 202     test[specie] = 0.0 ;
 203     if ( !AliQAv1::Instance()->IsEventSpecieSet(specie)) continue ;
 204     if (list[specie]->GetEntries() == 0) test[specie] = 0. ; // nothing to check
 205     else {
 206                         //get calib histos
 207                         TH2F * hdata  = (TH2F*)list[specie]->At(k2DRatioAmp) ;
 208                         TH1F * ratio = (TH1F*)list[specie]->At(kRatioDist) ;
 209
 210                         //get L1 histos
 211                         TH2F *hL1GammaPatch = (TH2F*)list[specie]->At(kGL1);
 212                         TH2F *hL1JetPatch = (TH2F*)list[specie]->At(kJL1);
 213                         TH1I *hFrameR = (TH1I*)list[specie]->At(kSTUTRU);
 214
 215                         //calib histo checker first:
 216
 217       if(hdata->GetEntries()!=0 && ratio->GetEntries()!=0) {
 218                                 //adding the lines to distinguish different SMs
 219                                 lstF = hdata->GetListOfFunctions();
 220                                 CleanListOfFunctions(lstF);
 221                                 lstF->Add(fLineCol->Clone());
 222                                 for(Int_t iLine = 0; iLine < 4; iLine++) {
 223                                         lstF->Add(fLineRow[iLine]->Clone());
 224                                 }
 225                                 //Now adding the text to for each SM
 226                                 for(Int_t iSM = 0 ; iSM < fgknSM ; iSM++){  //number of SMs loop start
 227                                         lstF->Add(fTextSM[iSM]->Clone());
 228                                 }
 229                                 //
 230                                 lstF = ratio->GetListOfFunctions();
 231                                 CleanListOfFunctions(lstF);
 232                                 //
 233                                 //now check the ratio histogram
 234                                 Double_t binContent = 0. ;
 235                                 Int_t nGoodTower = 0 ;
 236                                 Double_t rv = 0. ;
 237                                 for(Int_t ix = 1; ix <= hdata->GetNbinsX(); ix++) {
 238                                         for(Int_t iy = 1; iy <= hdata->GetNbinsY(); iy++) {
 239                                                 binContent = hdata->GetBinContent(ix, iy) ;
 240                                                 if (binContent < 1.2 && binContent > 0.8) nGoodTower++ ;
 241                                         }
 242                                 }
 243                                 rv = nGoodTower/nTot ;
 244                                 printf("%2.2f %% towers out of range [0.8, 1.2]\n", (1-rv)*100);
 245                                 if(fText){
 246                                         lstF->Add(fText->Clone()) ;
 247                                         fText->Clear() ;
 248
 249                                         fText->AddText(Form("%2.2f %% towers out of range [0.8, 1.2]", (1-rv)*100));
 250                                         if (rv < ratioThresh) {
 251                                                 test[specie] = ratioThresh;
 252                                                 // 2 lines text info for quality
 253                                                 fText->SetFillColor(2) ;
 254                                                 fText->AddText(Form("EMCAL = NOK, CALL EXPERTS!!!"));
 255                                         }
 256                                         else {
 257                                                 test[specie] = 1 - ratioThresh;
 258                                                 fText->SetFillColor(3) ;
 259                                                 fText->AddText(Form("EMCAL = OK, ENJOY..."));
 260                                         }
 261                                 }//fText
 262                         } // calib histo checking done
 263
 264                         //now L1 checks:
 265       if (specie != calibSpecieId) {
 266                         if(hL1GammaPatch->GetEntries() !=0 ) {
 267                                 lstF = hL1GammaPatch->GetListOfFunctions();
 268                                 CleanListOfFunctions(lstF);
 269
 270
 271                                 // Checker for L1GammaPatch
 272                                 //Double_t dL1GmeanTrig    = 1./2961.;
 273                                 //Double_t dL1GmeanTrigTRU = 1./32.;
 274                                 //Int_t sigmaG    = 100; // deviation from mean value (increased to 100)
 275                                 //Int_t sigmaGTRU = 5; // deviation from mean value for TRUs
 276                                 Double_t dL1GEntries = hL1GammaPatch->GetEntries();
 277                                 Int_t badL1G[48][64]   = {{0}} ;
 278                                 Int_t badL1GTRU[2][16] = {{0}} ;
 279                                 Int_t nBadL1G    = 0;
 280                                 Int_t nBadL1GTRU = 0;
 281                                 Double_t binContentTRU[2][16] = {{0.}};
 282                                 for(Int_t ix = 1; ix <=  hL1GammaPatch->GetNbinsX(); ix++) {
 283                                   for(Int_t iy = 1; iy <=  hL1GammaPatch->GetNbinsY(); iy++) {
 284                                     Double_t binContent = hL1GammaPatch->GetBinContent(ix, iy) ;
 285                                     if (binContent != 0) {
 286
 287                                       // fill counter for TRUs
 288                                       binContentTRU[(Int_t)((ix-1)/24)][(Int_t)((iy-1)/4)] += binContent;
 289
 290                                       // OLD METHOD (if a patch triggers > sigmaG * mean value (1/#patch positions total) says "hot spot !")
 291                                       // if ((double)binContent/(double)dL1GEntries > sigmaG*dL1GmeanTrig) {
 292                                       //        badL1G[ix-1][iy-1] += 1;
 293                                       //        nBadL1G += 1;
 294                                       // }
 295
 296                                       // NEW METHOD (if Rate > Threshold * ( (Number of towers or TRUs * Average rate) - Rate ) --> "hot spot !")
 297                                       // Thresold = how much does the noisy tower/TRU contribute to the rate
 298                                       //            1.0 --> Rate of noisy tower/TRU = Rate of all other towers/TRUs
 299                                       if (binContent/dL1GEntries > ThreshG / ( 1 + ThreshG )) {
 300                                         badL1G[ix-1][iy-1] += 1;
 301                                         nBadL1G += 1;
 302                                       }
 303                                     }
 304                                   }
 305                                 }
 306
 307                                 // check TRUs
 308                                 for(Int_t ix = 1; ix <=  2; ix++) {
 309                                   for(Int_t iy = 1; iy <=  16; iy++) {
 310                                     if(binContentTRU[ix-1][iy-1]/dL1GEntries >  ThreshG / ( 1 + ThreshG )) {
 311                                       badL1GTRU[ix-1][iy-1] += 1;
 312                                       nBadL1GTRU += 1;
 313                                     }
 314                                   }
 315                                 }
 316
 317                                 if(fTextL1[0]){
 318                                         lstF->Add(fTextL1[0]->Clone()) ;
 319                                         fTextL1[0]->Clear() ;
 320
 321                                         if (nBadL1G == 0 && nBadL1GTRU == 0 ) {
 322                                                 fTextL1[0]->SetFillColor(3) ;
 323                                                 fTextL1[0]->AddText(Form("L1 GAMMA TRIGGER = OK, ENJOY..."));
 324                                         }
 325                                         else if (nBadL1G == 0){
 326                                                 fTextL1[0]->SetFillColor(2) ;
 327                                                 fTextL1[0]->AddText(Form("HOT SPOT IN L1 GAMMA TRIGGER (TRU) = CALL EXPERT!!"));
 328
 329                                         }
 330                                         else{
 331                                           fTextL1[0]->SetFillColor(2) ;
 332                                           fTextL1[0]->AddText(Form("HOT SPOT IN L1 GAMMA TRIGGER = CALL EXPERT!!"));
 333                                           /*
 334                                             for(Int_t ix = 1; ix <=  hL1GammaPatch->GetNbinsX(); ix++) {
 335                                             for(Int_t iy = 1; iy <=  hL1GammaPatch->GetNbinsY(); iy++) {
 336                                                                 if(badL1G[ix-1][iy-1] != 0) printf("L1 Gamma patch with position x = %d, y = %d is out of range\n",ix,iy);
 337                                                                 }
 338                                                                 }
 339                                           */
 340                                         }
 341                                 }//fTextL1[0]
 342                         }// L1 gamma patch checking done
 343
 344                         if(hL1JetPatch->GetEntries() !=0) {
 345                                 lstF = hL1JetPatch->GetListOfFunctions();
 346                                 CleanListOfFunctions(lstF);
 347
 348                                 // Checker for L1JetPatch
 349                                 //Double_t dL1JmeanTrig = 1/126.;
 350                                 //Int_t sigmaJ = 5; // deviation from  mean value
 351                                 Double_t dL1JEntries = hL1JetPatch->GetEntries();
 352                                 Int_t badL1J[12][16] = {{0}} ;
 353                                 Int_t nBadL1J = 0;
 354                                 for(Int_t ix = 1; ix <=  hL1JetPatch->GetNbinsX(); ix++) {
 355                                   for(Int_t iy = 1; iy <=  hL1JetPatch->GetNbinsY(); iy++) {
 356                                     Double_t binContent = hL1JetPatch->GetBinContent(ix, iy) ;
 357                                     if (binContent != 0) {
 358
 359                                       // OLD METHOD  (if a patch triggers > sigmaJ * mean value (1/#patch positions total) says "hot spot !")
 360                                       // if ((double)binContent/(double)dL1JEntries > sigmaJ*dL1JmeanTrig) {
 361                                       //        badL1J[ix-1][iy-1] += 1 ;
 362                                       //        nBadL1J += 1;
 363                                       // }
 364
 365                                       // NEW METHOD (if Rate > Threshold * ( (Number of towers or TRUs * Average rate) - Rate ) --> "hot spot !")
 366                                       // Threshold: same definitionas for Gamma
 367                                        if ((double)binContent/(double)dL1JEntries > ThreshJ / ( 1 + ThreshJ )) {
 368                                         badL1J[ix-1][iy-1] += 1 ;
 369                                         nBadL1J += 1;
 370                                        }
 371                                     }
 372                                   }
 373                                 }
 374
 375                                 if(fTextL1[1]){
 376                                         lstF->Add(fTextL1[1]->Clone()) ;
 377                                         fTextL1[1]->Clear() ;
 378
 379                                         if (nBadL1J == 0) {
 380                                                 fTextL1[1]->SetFillColor(3) ;
 381                                                 fTextL1[1]->AddText(Form("L1 JET TRIGGER = OK, ENJOY..."));
 382                                         }
 383                                         else {
 384                                                 fTextL1[1]->SetFillColor(2) ;
 385                                                 fTextL1[1]->AddText(Form("HOT SPOT IN L1 JET TRIGGER = CALL EXPERT!!"));
 386 /*
 387                                                 for(Int_t ix = 1; ix <=  hL1JetPatch->GetNbinsX(); ix++) {
 388                                                         for(Int_t iy = 1; iy <=  hL1JetPatch->GetNbinsY(); iy++) {
 389                                                                 if(badL1J[ix-1][iy-1] != 0) printf("L1 Jet patch with position x = %d, y = %d is out of range\n",(4*ix-4),(4*iy-4));
 390                                                         }
 391                                                 }
 392 */
 393
 394                                         }
 395                                 }//fTextL1[1]
 396                         } // L1 Jet patch checking done
 397       } // if (specie != calibSpecieId) ..
 398
 399                         if(hFrameR->GetEntries() !=0) {
 400                                 lstF = hFrameR->GetListOfFunctions();
 401                                 CleanListOfFunctions(lstF);
 402
 403                                 Int_t badLink[32] = {0};
 404                                 Int_t nBadLink = 0;
 405                                 for(Int_t ix = 1; ix <= hFrameR->GetNbinsX(); ix++) {
 406                                         Double_t binContent = hFrameR->GetBinContent(ix) ;
 407                                         if (binContent == 0) {
 408                                                 badLink[ix-1] += 1;
 409                                                 nBadLink += 1;
 410                                         }
 411                                 }
 412                                 if(fTextL1[2]){
 413                                         lstF->Add(fTextL1[2]->Clone()) ;
 414                                         fTextL1[2]->Clear() ;
 415
 416                                 if (nBadLink < badLinkThresh) {
 417                                                 fTextL1[2]->SetFillColor(3) ;
 418                                                 fTextL1[2]->AddText(Form("LINK TRU-STU = OK, ENJOY..."));
 419                                         }
 420                                         else {
 421                                                 fTextL1[2]->SetFillColor(2) ;
 422                                                 fTextL1[2]->AddText(Form("PROBLEM WITH TRU-STU LINK = CALL EXPERT!!"));
 423 /*
 424                                                 for(Int_t ix = 0; ix <= hFrameR->GetNbinsX(); ix++) {
 425                                                         if(badLink[ix] != 0) printf("STU link with TRU %d is out\n",ix);
 426                                                 }
 427 */
 428                                         }
 429                                 }//fTextL1[2]
 430                         } // Checker for link TRU-STU done
 431     }
 432
 433   } // specie
 434 }
 435 //______________________________________________________________________________
 436
 437 void AliEMCALQAChecker::Init(const AliQAv1::DETECTORINDEX_t det)
 438 {
 439         /// intialises QA and QA checker settings
 440         AliQAv1::Instance(det) ;
 441         Float_t hiValue[AliQAv1::kNBIT] ;
 442         Float_t lowValue[AliQAv1::kNBIT] ;
 443         lowValue[AliQAv1::kINFO]      = 0.0   ;
 444         hiValue[AliQAv1::kINFO]       = 0.1 ;
 445         lowValue[AliQAv1::kWARNING]   = 0.1 ;
 446   hiValue[AliQAv1::kWARNING]    = 0.5 ;
 447         lowValue[AliQAv1::kERROR]     = 0.5   ;
 448         hiValue[AliQAv1::kERROR]      = 0.8 ;
 449         lowValue[AliQAv1::kFATAL]     = 0.8   ;
 450         hiValue[AliQAv1::kFATAL]      = 1.0 ;
 451         SetHiLo(&hiValue[0], &lowValue[0]) ;
 452 }
 453
 454 //______________________________________________________________________________
 455
 456 void AliEMCALQAChecker::CleanListOfFunctions(TList *list)
 457 { // clean up
 458
 459         if (list) {
 460                 TObject *stats = list->FindObject("stats"); list->Remove(stats);
 461                 TObject *obj;
 462                 while ((obj = list->First())) { while(list->Remove(obj)) { } delete obj; }
 463                         if (stats) list->Add(stats);
 464                 }
 465         else {
 466                 AliWarning(Form("Checker : empty list of data functions; returning"));
 467                 return;
 468         }
 469
 470 }
 471
 472