EMCAL/AliEMCALQADataMakerRec.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 Based on the QA code for PHOS written by Yves Schutz July 2007
  17
  18 Authors:  J.Klay (Cal Poly) May 2008
  19           S. Salur LBL April 2008
  20
  21 */
  22
  23 // --- ROOT system ---
  24 #include <TClonesArray.h>
  25 #include <TFile.h>
  26 #include <TH1F.h>
  27 #include <TH1I.h>
  28 #include <TH2F.h>
  29 #include <TProfile.h>
  30
  31 // --- Standard library ---
  32
  33
  34 // --- AliRoot header files ---
  35 #include "AliESDCaloCluster.h"
  36 #include "AliESDCaloCells.h"
  37 #include "AliESDEvent.h"
  38 #include "AliLog.h"
  39 #include "AliEMCALQADataMakerRec.h"
  40 #include "AliQAChecker.h"
  41 #include "AliEMCALRecPoint.h"
  42 #include "AliEMCALRawUtils.h"
  43 #include "AliEMCALReconstructor.h"
  44 #include "AliEMCALRecParam.h"
  45 #include "AliRawReader.h"
  46 #include "AliCaloRawStream.h"
  47
  48 ClassImp(AliEMCALQADataMakerRec)
  49
  50 //____________________________________________________________________________
  51   AliEMCALQADataMakerRec::AliEMCALQADataMakerRec() :
  52     AliQADataMakerRec(AliQAv1::GetDetName(AliQAv1::kEMCAL), "EMCAL Quality Assurance Data Maker"),
  53     fSuperModules(4) // FIXME!!! number of SuperModules; 4 for 2009; update default to 12 for later runs..
  54 {
  55   // ctor
  56 }
  57
  58 //____________________________________________________________________________
  59 AliEMCALQADataMakerRec::AliEMCALQADataMakerRec(const AliEMCALQADataMakerRec& qadm) :
  60   AliQADataMakerRec()
  61 {
  62   //copy ctor
  63   SetName((const char*)qadm.GetName()) ;
  64   SetTitle((const char*)qadm.GetTitle());
  65   fSuperModules = qadm.GetSuperModules();
  66 }
  67
  68 //__________________________________________________________________
  69 AliEMCALQADataMakerRec& AliEMCALQADataMakerRec::operator = (const AliEMCALQADataMakerRec& qadm )
  70 {
  71   // Equal operator.
  72   this->~AliEMCALQADataMakerRec();
  73   new(this) AliEMCALQADataMakerRec(qadm);
  74   return *this;
  75 }
  76
  77 //____________________________________________________________________________
  78 void AliEMCALQADataMakerRec::EndOfDetectorCycle(AliQAv1::TASKINDEX_t task, TObjArray ** list)
  79 {
  80   //Detector specific actions at end of cycle
  81   // do the QA checking
  82   AliQAChecker::Instance()->Run(AliQAv1::kEMCAL, task, list) ;
  83 }
  84
  85 //____________________________________________________________________________
  86 void AliEMCALQADataMakerRec::InitESDs()
  87 {
  88   //Create histograms to controll ESD
  89   const Bool_t expert   = kTRUE ;
  90   const Bool_t image    = kTRUE ;
  91
  92   TH1F * h1 = new TH1F("hESDCaloClusterE",  "ESDs CaloCluster energy in EMCAL",    200, 0., 20.) ;
  93   h1->Sumw2() ;
  94   Add2ESDsList(h1, kESDCaloClusE, !expert, image)  ;
  95
  96   TH1I * h2 = new TH1I("hESDCaloClusterM", "ESDs CaloCluster multiplicity in EMCAL", 100, 0,  100) ;
  97   h2->Sumw2() ;
  98   Add2ESDsList(h2, kESDCaloClusM, !expert, image)  ;
  99
 100   TH1F * h3 = new TH1F("hESDCaloCellA",  "ESDs CaloCell amplitude in EMCAL",    500, 0., 250.) ;
 101   h3->Sumw2() ;
 102   Add2ESDsList(h3, kESDCaloCellA, !expert, image)  ;
 103
 104   TH1I * h4 = new TH1I("hESDCaloCellM", "ESDs CaloCell multiplicity in EMCAL", 200, 0,  1000) ;
 105   h4->Sumw2() ;
 106   Add2ESDsList(h4, kESDCaloCellM, !expert, image) ;
 107
 108 }
 109
 110 //____________________________________________________________________________
 111 void AliEMCALQADataMakerRec::InitRecPoints()
 112 {
 113   // create Reconstructed Points histograms in RecPoints subdir
 114   const Bool_t expert   = kTRUE ;
 115   const Bool_t image    = kTRUE ;
 116
 117   TH1F* h0 = new TH1F("hEMCALRpE","EMCAL RecPoint energies",200, 0.,20.); //GeV
 118   h0->Sumw2();
 119   Add2RecPointsList(h0,kRecPE, !expert, image);
 120
 121   TH1I* h1 = new TH1I("hEMCALRpM","EMCAL RecPoint multiplicities",100,0,100);
 122   h1->Sumw2();
 123   Add2RecPointsList(h1,kRecPM, !expert, image);
 124
 125   TH1I* h2 = new TH1I("hEMCALRpDigM","EMCAL RecPoint Digit Multiplicities",20,0,20);
 126   h2->Sumw2();
 127   Add2RecPointsList(h2,kRecPDigM, !expert, image);
 128
 129 }
 130
 131 //____________________________________________________________________________
 132 void AliEMCALQADataMakerRec::InitRaws()
 133 {
 134   // create Raws histograms in Raws subdir
 135    const Bool_t expert   = kTRUE ;
 136    const Bool_t saveCorr = kTRUE ;
 137    const Bool_t image    = kTRUE ;
 138
 139   int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
 140   int nTot = fSuperModules * nTowersPerSM; // max number of towers in all SuperModules
 141
 142   // counter info: number of channels per event (bins are SM index)
 143   TProfile * h0 = new TProfile("hLowEmcalSupermodules", "Low Gain EMC: # of towers vs SuperMod",
 144                                fSuperModules, -0.5, fSuperModules-0.5) ;
 145   Add2RawsList(h0, kNsmodLG, !expert, image, !saveCorr) ;
 146   TProfile * h1 = new TProfile("hHighEmcalSupermodules", "High Gain EMC: # of towers vs SuperMod",
 147                                fSuperModules, -0.5, fSuperModules-0.5) ;
 148   Add2RawsList(h1, kNsmodHG, !expert, image, !saveCorr) ;
 149
 150   // where did max sample occur? (bins are towers)
 151   TProfile * h2 = new TProfile("hLowEmcalRawtime", "Low Gain EMC: Time at Max vs towerId",
 152                                nTot, -0.5, nTot-0.5) ;
 153   Add2RawsList(h2, kTimeLG, !expert, image, !saveCorr) ;
 154   TProfile * h3 = new TProfile("hHighEmcalRawtime", "High Gain EMC: Time at Max vs towerId",
 155                                nTot, -0.5, nTot-0.5) ;
 156   Add2RawsList(h3, kTimeHG, !expert, image, !saveCorr) ;
 157
 158   // how much above pedestal was the max sample?  (bins are towers)
 159   TProfile * h4 = new TProfile("hLowEmcalRawMaxMinusMin", "Low Gain EMC: Max - Min vs towerId",
 160                                nTot, -0.5, nTot-0.5) ;
 161   Add2RawsList(h4, kSigLG, !expert, image, !saveCorr) ;
 162   TProfile * h5 = new TProfile("hHighEmcalRawMaxMinusMin", "High Gain EMC: Max - Min vs towerId",
 163                                nTot, -0.5, nTot-0.5) ;
 164   Add2RawsList(h5, kSigHG, !expert, image, !saveCorr) ;
 165
 166   // total counter: channels per event
 167   TH1I * h6 = new TH1I("hLowNtot", "Low Gain EMC: Total Number of found towers", 200, 0, nTot) ;
 168   h6->Sumw2() ;
 169   Add2RawsList(h6, kNtotLG, !expert, image, !saveCorr) ;
 170   TH1I * h7 = new TH1I("hHighNtot", "High Gain EMC: Total Number of found towers", 200,0, nTot) ;
 171   h7->Sumw2() ;
 172   Add2RawsList(h7, kNtotHG, !expert, image, !saveCorr) ;
 173
 174   // pedestal (bins are towers)
 175   TProfile * h8 = new TProfile("hLowEmcalRawPed", "Low Gain EMC: Pedestal vs towerId",
 176                                nTot, -0.5, nTot-0.5) ;
 177   Add2RawsList(h8, kPedLG, !expert, image, !saveCorr) ;
 178   TProfile * h9 = new TProfile("hHighEmcalRawPed", "High Gain EMC: Pedestal vs towerId",
 179                                nTot, -0.5, nTot-0.5) ;
 180   Add2RawsList(h9, kPedHG, !expert, image, !saveCorr) ;
 181
 182   // pedestal rms (standard dev = sqrt of variance estimator for pedestal) (bins are towers)
 183   TProfile * h10 = new TProfile("hLowEmcalRawPedRMS", "Low Gain EMC: Pedestal RMS vs towerId",
 184                                 nTot, -0.5, nTot-0.5) ;
 185   Add2RawsList(h10, kPedRMSLG, !expert, image, !saveCorr) ;
 186   TProfile * h11 = new TProfile("hHighEmcalRawPedRMS", "High Gain EMC: Pedestal RMS vs towerId",
 187                                 nTot, -0.5, nTot-0.5) ;
 188   Add2RawsList(h11, kPedRMSHG, !expert, image, !saveCorr) ;
 189
 190 }
 191
 192 //____________________________________________________________________________
 193 void AliEMCALQADataMakerRec::MakeESDs(AliESDEvent * esd)
 194 {
 195   // make QA data from ESDs
 196
 197   Int_t nTot = 0 ;
 198   for ( Int_t index = 0; index < esd->GetNumberOfCaloClusters() ; index++ ) {
 199     AliESDCaloCluster * clu = esd->GetCaloCluster(index) ;
 200     if( clu->IsEMCAL() ) {
 201       GetESDsData(kESDCaloClusE)->Fill(clu->E()) ;
 202       nTot++ ;
 203     }
 204   }
 205   GetESDsData(kESDCaloClusM)->Fill(nTot) ;
 206
 207   //fill calo cells
 208   AliESDCaloCells* cells = esd->GetEMCALCells();
 209   GetESDsData(kESDCaloCellM)->Fill(cells->GetNumberOfCells()) ;
 210
 211   for ( Int_t index = 0; index < cells->GetNumberOfCells() ; index++ ) {
 212     GetESDsData(kESDCaloCellA)->Fill(cells->GetAmplitude(index)) ;
 213   }
 214
 215 }
 216
 217 //____________________________________________________________________________
 218 void AliEMCALQADataMakerRec::MakeRaws(AliRawReader* rawReader)
 219 {
 220   //Fill prepared histograms with Raw digit properties
 221
 222   //Raw histogram filling not yet implemented
 223   //
 224   //Need to figure out how to get the info we want without having to
 225   //actually run Raw2Digits twice.
 226   //I suspect what we actually want is a raw digits method, not a true
 227   //emcal raw data method, but this doesn't seem to be allowed in
 228   //AliQADataMakerRec.h
 229
 230   // For now, to avoid redoing the expensive signal fits we just
 231   // look at max vs min of the signal spextra, a la online usage in
 232   // AliCaloCalibPedestal
 233
 234   rawReader->Reset() ;
 235   AliCaloRawStream in(rawReader,"EMCAL");
 236
 237   // setup
 238   int nTowersPerSM = 1152; // number of towers in a SuperModule; 24x48
 239   int nRows = 24; // number of rows per SuperModule
 240   int sampleMin = 0;
 241   int sampleMax = 0x3ff; // 1023 = 10-bit range
 242
 243   // for the pedestal calculation
 244   Bool_t selectPedestalSamples = kTRUE;
 245   int firstPedestalSample = 0;
 246   int lastPedestalSample = 15;
 247
 248   // SM counters; decl. should be safe, assuming we don't get more than 12 SuperModules..
 249   int nTotalSMLG[12] = {0};
 250   int nTotalSMHG[12] = {0};
 251
 252   // indices for the reading
 253   int iSM = 0;
 254   int sample = 0;
 255   int gain = 0;
 256   int time = 0;
 257   // counters, on sample level
 258   int i = 0; // the sample number in current event.
 259   int max = sampleMin, min = sampleMax;//Use these for picking the pedestal
 260   int maxTime = 0;
 261
 262   // for the pedestal calculation
 263   int sampleSum = 0; // sum of samples
 264   int squaredSampleSum = 0; // sum of samples squared
 265   int nSum = 0; // number of samples in sum
 266   // calc. quantities
 267   double meanPed = 0, squaredMean = 0, rmsPed = 0;
 268
 269   while (in.Next()) { // loop over input stream
 270     sample = in.GetSignal(); //Get the adc signal
 271     time = in.GetTime();
 272     if (sample < min) { min = sample; }
 273     if (sample > max) {
 274       max = sample;
 275       maxTime = time;
 276     }
 277     i++;
 278
 279     // should we add it for the pedestal calculation?
 280     if ( (firstPedestalSample<=time && time<=lastPedestalSample) || // sample time in range
 281          !selectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
 282       sampleSum += sample;
 283       squaredSampleSum += sample*sample;
 284       nSum++;
 285     }
 286
 287     if ( i >= in.GetTimeLength()) {
 288       // calculate pedesstal estimate: mean of possibly selected samples
 289       if (nSum > 0) {
 290         meanPed = sampleSum / (1.0 * nSum);
 291         squaredMean = squaredSampleSum / (1.0 * nSum);
 292         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
 293         rmsPed = sqrt(squaredMean - meanPed*meanPed);
 294       }
 295       else {
 296         meanPed = 0;
 297         squaredMean = 0;
 298         rmsPed  = 0;
 299       }
 300
 301       //If we're here then we're done with this tower
 302       gain = -1; // init to not valid value
 303       if ( in.IsLowGain() ) {
 304         gain = 0;
 305       }
 306       else if ( in.IsHighGain() ) {
 307         gain = 1;
 308       }
 309
 310       iSM = in.GetModule(); //The modules are numbered starting from 0
 311
 312       if (iSM>=0 && iSM<fSuperModules) { // valid module reading, can go on with filling
 313
 314         int towerId = iSM*nTowersPerSM + in.GetColumn()*nRows + in.GetRow();
 315
 316         if (gain == 0) {
 317           //fill the low gain histograms, and counters
 318           nTotalSMLG[iSM]++; // one more channel found
 319           GetRawsData(kSigLG)->Fill(towerId, max - min);
 320           GetRawsData(kTimeLG)->Fill(towerId, maxTime);
 321           if (nSum>0) { // only fill pedestal info in case it could be calculated
 322             GetRawsData(kPedLG)->Fill(towerId, meanPed);
 323             GetRawsData(kPedRMSLG)->Fill(towerId, rmsPed);
 324           }
 325         } // gain==0
 326         else if (gain == 1) {
 327           //fill the high gain ones
 328           nTotalSMHG[iSM]++; // one more channel found
 329           GetRawsData(kSigHG)->Fill(towerId, max - min);
 330           GetRawsData(kTimeHG)->Fill(towerId, maxTime);
 331           if (nSum>0) { // only fill pedestal info in case it could be calculated
 332             GetRawsData(kPedHG)->Fill(towerId, meanPed);
 333             GetRawsData(kPedRMSHG)->Fill(towerId, rmsPed);
 334           }
 335         }
 336       } // SM index OK
 337
 338       // reset counters
 339       max = sampleMin; min = sampleMax;
 340       maxTime = 0;
 341       i = 0;
 342       // also pedestal calc counters
 343       sampleSum = 0; // sum of samples
 344       squaredSampleSum = 0; // sum of samples squared
 345       nSum = 0; // number of samples in sum
 346
 347     }//End if, of channel
 348
 349   }//end while, of stream
 350
 351   // let's also fill the SM and event counter histograms
 352   int nTotalHG = 0;
 353   int nTotalLG = 0;
 354   for (iSM=0; iSM<fSuperModules; iSM++) {
 355     nTotalLG += nTotalSMLG[iSM];
 356     nTotalHG += nTotalSMHG[iSM];
 357     GetRawsData(kNsmodLG)->Fill(iSM, nTotalSMLG[iSM]);
 358     GetRawsData(kNsmodHG)->Fill(iSM, nTotalSMHG[iSM]);
 359   }
 360   GetRawsData(kNtotLG)->Fill(nTotalLG);
 361   GetRawsData(kNtotHG)->Fill(nTotalHG);
 362
 363   // just in case the next rawreader consumer forgets to reset; let's do it here again..
 364   rawReader->Reset() ;
 365
 366   return;
 367 }
 368
 369 //____________________________________________________________________________
 370 void AliEMCALQADataMakerRec::MakeRecPoints(TTree * clustersTree)
 371 {
 372   // makes data from RecPoints
 373   TBranch *emcbranch = clustersTree->GetBranch("EMCALECARP");
 374   if (!emcbranch) {
 375     AliError("can't get the branch with the EMCAL clusters !");
 376     return;
 377   }
 378   TObjArray * emcrecpoints = new TObjArray(100) ;
 379   emcbranch->SetAddress(&emcrecpoints);
 380   emcbranch->GetEntry(0);
 381
 382   GetRecPointsData(kRecPM)->Fill(emcrecpoints->GetEntriesFast()) ;
 383   TIter next(emcrecpoints) ;
 384   AliEMCALRecPoint * rp ;
 385   while ( (rp = dynamic_cast<AliEMCALRecPoint *>(next())) ) {
 386     GetRecPointsData(kRecPE)->Fill( rp->GetEnergy()) ;
 387     GetRecPointsData(kRecPDigM)->Fill(rp->GetMultiplicity());
 388   }
 389   emcrecpoints->Delete();
 390   delete emcrecpoints;
 391
 392 }
 393
 394 //____________________________________________________________________________
 395 void AliEMCALQADataMakerRec::StartOfDetectorCycle()
 396 {
 397   //Detector specific actions at start of cycle
 398
 399 }
 400