EMCAL/AliCaloCalibPedestal.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 //* $Id$ */
  16
  17 //________________________________________________________________________
  18 //
  19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
  20 // It can be created and used a la (ctor):
  21 /*
  22   //Create the object for making the histograms
  23   fPedestals = new AliCaloCalibPedestal( fDetType );
  24   // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL
  25   fNumModules = fPedestals->GetModules();
  26 */
  27 // fed an event:
  28 //  fPedestals->ProcessEvent(fCaloRawStream);
  29 // asked to draw histograms:
  30 //  fPedestals->GetDeadMap(i)->Draw("col");
  31 // or
  32 //  fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz");
  33 // etc.
  34 // The pseudo-code examples above were from the first implementation in MOOD (summer 2007).
  35 //________________________________________________________________________
  36
  37 //#include "TCanvas.h"
  38 #include "TH1.h"
  39 #include "TF1.h"
  40 #include "TFile.h"
  41 #include <fstream>
  42 #include <sstream>
  43 #include <iostream>
  44 #include <cmath>
  45
  46 #include "AliRawReader.h"
  47 #include "AliCaloRawStreamV3.h"
  48
  49 //The include file
  50 #include "AliCaloCalibPedestal.h"
  51
  52 ClassImp(AliCaloCalibPedestal)
  53
  54 using namespace std;
  55
  56 // ctor; initialize everything in order to avoid compiler warnings
  57 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
  58   TObject(),
  59   fPedestalLowGain(),
  60   fPedestalHighGain(),
  61   fPedestalRMSLowGain(),
  62   fPedestalRMSHighGain(),
  63   fPeakMinusPedLowGain(),
  64   fPeakMinusPedHighGain(),
  65   fPeakMinusPedHighGainHisto(),
  66   fPedestalLowGainDiff(),
  67   fPedestalHighGainDiff(),
  68   fPeakMinusPedLowGainDiff(),
  69   fPeakMinusPedHighGainDiff(),
  70   fPedestalLowGainRatio(),
  71   fPedestalHighGainRatio(),
  72   fPeakMinusPedLowGainRatio(),
  73   fPeakMinusPedHighGainRatio(),
  74   fDeadMap(),
  75   fNEvents(0),
  76   fNChanFills(0),
  77   fDeadTowers(0),
  78   fNewDeadTowers(0),
  79   fResurrectedTowers(0),
  80   fReference(0),
  81   fDetType(kNone),
  82   fColumns(0),
  83   fRows(0),
  84   fModules(0),
  85   fRowMin(0),
  86   fRowMax(0),
  87   fRowMultiplier(0),
  88   fCaloString(),
  89   fMapping(NULL),
  90   fRunNumber(-1),
  91   fSelectPedestalSamples(kTRUE),
  92   fFirstPedestalSample(0),
  93   fLastPedestalSample(15),
  94   fDeadThreshold(5),
  95   fWarningThreshold(50),
  96   fWarningFraction(0.002),
  97   fHotSigma(5)
  98 {
  99   //Default constructor. First we set the detector-type related constants.
 100   if (detectorType == kPhos) {
 101     fColumns = fgkPhosCols;
 102     fRows = fgkPhosRows;
 103     fModules = fgkPhosModules;
 104     fCaloString = "PHOS";
 105     fRowMin = -1*fRows;
 106     fRowMax = 0;
 107     fRowMultiplier = -1;
 108   }
 109   else {
 110     //We'll just trust the enum to keep everything in line, so that if detectorType
 111     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
 112     //case, if someone intentionally gives another number
 113     fColumns = AliEMCALGeoParams::fgkEMCALCols;
 114     fRows = AliEMCALGeoParams::fgkEMCALRows;
 115     fModules = AliEMCALGeoParams::fgkEMCALModules;
 116     fCaloString = "EMCAL";
 117     fRowMin = 0;
 118     fRowMax = fRows;
 119     fRowMultiplier = 1;
 120   }
 121   fDetType = detectorType;
 122
 123   //Then, loop for the requested number of modules
 124   TString title, name;
 125   for (int i = 0; i < fModules; i++) {
 126     //Pedestals, low gain
 127     name = "hPedlowgain";
 128     name += i;
 129     title = "Pedestals, low gain, module ";
 130     title += i;
 131     fPedestalLowGain.Add(new TProfile2D(name, title,
 132                                         fColumns, 0.0, fColumns,
 133                                         fRows, fRowMin, fRowMax,"s"));
 134
 135     //Pedestals, high gain
 136     name = "hPedhighgain";
 137     name += i;
 138     title = "Pedestals, high gain, module ";
 139     title += i;
 140     fPedestalHighGain.Add(new TProfile2D(name, title,
 141                                          fColumns, 0.0, fColumns,
 142                                          fRows, fRowMin, fRowMax,"s"));
 143     //All Samples, low gain
 144     name = "hPedestalRMSlowgain";
 145     name += i;
 146     title = "Pedestal RMS, low gain, module ";
 147     title += i;
 148     fPedestalRMSLowGain.Add(new TProfile2D(name, title,
 149                                         fColumns, 0.0, fColumns,
 150                                         fRows, fRowMin, fRowMax,"s"));
 151
 152     //All Samples, high gain
 153     name = "hPedestalRMShighgain";
 154     name += i;
 155     title = "Pedestal RMS, high gain, module ";
 156     title += i;
 157     fPedestalRMSHighGain.Add(new TProfile2D(name, title,
 158                                          fColumns, 0.0, fColumns,
 159                                          fRows, fRowMin, fRowMax,"s"));
 160
 161     //Peak-Pedestals, low gain
 162     name = "hPeakMinusPedlowgain";
 163     name += i;
 164     title = "Peak-Pedestal, low gain, module ";
 165     title += i;
 166     fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
 167                                             fColumns, 0.0, fColumns,
 168                                             fRows, fRowMin, fRowMax,"s"));
 169
 170     //Peak-Pedestals, high gain
 171     name = "hPeakMinusPedhighgain";
 172     name += i;
 173     title = "Peak-Pedestal, high gain, module ";
 174     title += i;
 175     fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
 176                                              fColumns, 0.0, fColumns,
 177                                              fRows, fRowMin, fRowMax,"s"));
 178
 179     //Peak-Pedestals, high gain - TH2F histo
 180     name = "hPeakMinusPedhighgainHisto";
 181     name += i;
 182     title = "Peak-Pedestal, high gain, module ";
 183     title += i;
 184     fPeakMinusPedHighGainHisto.Add(new TH2F(name, title,
 185                                             fColumns*fRows, 0.0, fColumns*fRows,
 186                                             100, 0, 1000));
 187
 188     name = "hDeadMap";
 189     name += i;
 190     title = "Dead map, module ";
 191     title += i;
 192     fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
 193                           fRows, fRowMin, fRowMax));
 194
 195   }//end for nModules create the histograms
 196
 197   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
 198   fPedestalLowGain.Compress();
 199   fPedestalHighGain.Compress();
 200   fPedestalRMSLowGain.Compress();
 201   fPedestalRMSHighGain.Compress();
 202   fPeakMinusPedLowGain.Compress();
 203   fPeakMinusPedHighGain.Compress();
 204   fPeakMinusPedHighGainHisto.Compress();
 205   fDeadMap.Compress();
 206   //Make them the owners of the profiles, so we don't need to care about deleting them
 207   //fPedestalLowGain.SetOwner();
 208   //fPedestalHighGain.SetOwner();
 209   //fPeakMinusPedLowGain.SetOwner();
 210   //fPeakMinusPedHighGain.SetOwner();
 211
 212 }
 213
 214 // dtor
 215 //_____________________________________________________________________
 216 AliCaloCalibPedestal::~AliCaloCalibPedestal()
 217 {
 218   if (fReference) delete fReference;//Delete the reference object, if it has been loaded
 219   //TObjArray will delete the histos/profiles when it is deleted.
 220 }
 221
 222 // copy ctor
 223 //_____________________________________________________________________
 224 AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
 225   TObject(ped),
 226   fPedestalLowGain(),
 227   fPedestalHighGain(),
 228   fPedestalRMSLowGain(),
 229   fPedestalRMSHighGain(),
 230   fPeakMinusPedLowGain(),
 231   fPeakMinusPedHighGain(),
 232   fPeakMinusPedHighGainHisto(),
 233   fPedestalLowGainDiff(),
 234   fPedestalHighGainDiff(),
 235   fPeakMinusPedLowGainDiff(),
 236   fPeakMinusPedHighGainDiff(),
 237   fPedestalLowGainRatio(),
 238   fPedestalHighGainRatio(),
 239   fPeakMinusPedLowGainRatio(),
 240   fPeakMinusPedHighGainRatio(),
 241   fDeadMap(),
 242   fNEvents(ped.GetNEvents()),
 243   fNChanFills(ped.GetNChanFills()),
 244   fDeadTowers(ped.GetDeadTowerCount()),
 245   fNewDeadTowers(ped.GetDeadTowerNew()),
 246   fResurrectedTowers(ped.GetDeadTowerResurrected()),
 247   fReference( 0 ), //! note that we do not try to copy the reference info here
 248   fDetType(ped.GetDetectorType()),
 249   fColumns(ped.GetColumns()),
 250   fRows(ped.GetRows()),
 251   fModules(ped.GetModules()),
 252   fRowMin(ped.GetRowMin()),
 253   fRowMax(ped.GetRowMax()),
 254   fRowMultiplier(ped.GetRowMultiplier()),
 255   fCaloString(ped.GetCaloString()),
 256   fMapping(NULL), //! note that we are not copying the map info
 257   fRunNumber(ped.GetRunNumber()),
 258   fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
 259   fFirstPedestalSample(ped.GetFirstPedestalSample()),
 260   fLastPedestalSample(ped.GetLastPedestalSample()),
 261   fDeadThreshold(ped.GetDeadThreshold()),
 262   fWarningThreshold(ped.GetWarningThreshold()),
 263   fWarningFraction(ped.GetWarningFraction()),
 264   fHotSigma(ped.GetHotSigma())
 265 {
 266   // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
 267   //DS: this has not really been tested yet..
 268   for (int i = 0; i < fModules; i++) {
 269     fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
 270     fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
 271     fPedestalRMSLowGain.Add( ped.GetPedRMSProfileLowGain(i) );
 272     fPedestalRMSHighGain.Add( ped.GetPedRMSProfileHighGain(i) );
 273     fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
 274     fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
 275     fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) );
 276
 277     fDeadMap.Add( ped.GetDeadMap(i) );
 278   }//end for nModules
 279
 280   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
 281   fPedestalLowGain.Compress();
 282   fPedestalHighGain.Compress();
 283   fPedestalRMSLowGain.Compress();
 284   fPedestalRMSHighGain.Compress();
 285   fPeakMinusPedLowGain.Compress();
 286   fPeakMinusPedHighGain.Compress();
 287   fPeakMinusPedHighGainHisto.Compress();
 288
 289   fDeadMap.Compress();
 290 }
 291
 292 // assignment operator; use copy ctor to make life easy..
 293 //_____________________________________________________________________
 294 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
 295 {
 296   // assignment operator; use copy ctor
 297   if (&source == this) return *this;
 298
 299   new (this) AliCaloCalibPedestal(source);
 300   return *this;
 301 }
 302
 303 //_____________________________________________________________________
 304 void AliCaloCalibPedestal::Reset()
 305 {
 306   // Reset all arrays/histograms
 307   for (int i = 0; i < fModules; i++) {
 308     GetPedProfileLowGain(i)->Reset();
 309     GetPedProfileHighGain(i)->Reset();
 310     GetPeakProfileLowGain(i)->Reset();
 311     GetPeakProfileHighGain(i)->Reset();
 312     GetPeakHighGainHisto(i)->Reset();
 313     GetDeadMap(i)->Reset();
 314
 315     if (!fPedestalLowGainDiff.IsEmpty()) {
 316       //This means that the comparison profiles have been created.
 317
 318       GetPedProfileLowGainDiff(i)->Reset();
 319       GetPedProfileHighGainDiff(i)->Reset();
 320       GetPeakProfileLowGainDiff(i)->Reset();
 321       GetPeakProfileHighGainDiff(i)->Reset();
 322
 323       GetPedProfileLowGainRatio(i)->Reset();
 324       GetPedProfileHighGainRatio(i)->Reset();
 325       GetPeakProfileLowGainRatio(i)->Reset();
 326       GetPeakProfileHighGainRatio(i)->Reset();
 327     }
 328   }
 329   fNEvents = 0;
 330   fNChanFills = 0;
 331   fDeadTowers = 0;
 332   fNewDeadTowers = 0;
 333   fResurrectedTowers = 0;
 334
 335   //To think about: should fReference be deleted too?... let's not do it this time, at least...
 336 }
 337
 338 // Parameter/cut handling
 339 //_____________________________________________________________________
 340 void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
 341 {
 342   // Note: this method is a bit more complicated than it really has to be
 343   // - allowing for multiple entries per line, arbitrary order of the
 344   // different variables etc. But I wanted to try and do this in as
 345   // correct a C++ way as I could (as an exercise).
 346
 347   static const string delimitor("::");
 348
 349   // open, check input file
 350   ifstream in( parameterFile );
 351   if( !in ) {
 352     printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
 353     return;
 354   }
 355
 356
 357   // read in
 358   char readline[1024];
 359   while ((in.rdstate() & ios::failbit) == 0 ) {
 360
 361     // Read into the raw char array and then construct a string
 362     // to do the searching
 363     in.getline(readline, 1024);
 364     istringstream s(readline);
 365
 366     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
 367
 368       string keyValue;
 369       s >> keyValue;
 370
 371       // check stream status
 372       if( s.rdstate() & ios::failbit ) break;
 373
 374       // skip rest of line if comments found
 375       if( keyValue.substr( 0, 2 ) == "//" ) break;
 376
 377       // look for "::" in keyValue pair
 378       size_t position = keyValue.find( delimitor );
 379       if( position == string::npos ) {
 380         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 381       }
 382
 383       // split keyValue pair
 384       string key( keyValue.substr( 0, position ) );
 385       string value( keyValue.substr( position+delimitor.size(),
 386                                       keyValue.size()-delimitor.size() ) );
 387
 388       // check value does not contain a new delimitor
 389       if( value.find( delimitor ) != string::npos ) {
 390         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 391       }
 392
 393       // debug: check key value pair
 394       // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 395
 396       // if the key matches with something we expect, we assign the new value
 397       istringstream iss(value);
 398       // the comparison strings defined at the beginning of this method
 399       if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) {
 400         printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 401
 402         if (key == "fFirstPedestalSample") {
 403           iss >> fFirstPedestalSample;
 404         }
 405         else if (key == "fLastPedestalSample") {
 406           iss >> fLastPedestalSample;
 407         }
 408         else if (key == "fDeadThreshold") {
 409           iss >> fDeadThreshold;
 410         }
 411         else if (key == "fWarningThreshold") {
 412           iss >> fWarningThreshold;
 413         }
 414         else if (key == "fWarningFraction") {
 415           iss >> fWarningFraction;
 416         }
 417         else if (key == "fHotSigma") {
 418           iss >> fHotSigma;
 419         }
 420
 421       } // some match
 422
 423     }
 424   }
 425
 426   in.close();
 427   return;
 428
 429 }
 430
 431 //_____________________________________________________________________
 432 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
 433 {
 434   //Write parameters in file.
 435
 436   static const string delimitor("::");
 437   ofstream out( parameterFile );
 438   out << "// " << parameterFile << endl;
 439   out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
 440   out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
 441   out << "fDeadThreshold" << "::" << fDeadThreshold << endl;
 442   out << "fWarningThreshold" << "::" << fWarningThreshold << endl;
 443   out << "fWarningFraction" << "::" << fWarningFraction << endl;
 444   out << "fHotSigma" << "::" << fHotSigma << endl;
 445
 446   out.close();
 447   return;
 448 }
 449
 450 //_____________________________________________________________________
 451 Bool_t AliCaloCalibPedestal::AddInfo(const AliCaloCalibPedestal *ped)
 452 {
 453   // just do this for the basic histograms/profiles that get filled in ProcessEvent
 454   // may not have data for all modules, but let's just Add everything..
 455   for (int i = 0; i < fModules; i++) {
 456     GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
 457     GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
 458     GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
 459     GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
 460     GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) );
 461
 462   }//end for nModules
 463
 464   // DeadMap; Diff profiles etc would need to be redone after this operation
 465
 466   return kTRUE;//We succesfully added info from the supplied object
 467 }
 468
 469 //_____________________________________________________________________
 470 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
 471 {
 472   // if fMapping is NULL the rawstream will crate its own mapping
 473   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 474   if (fDetType == kEmCal) {
 475     rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
 476   }
 477   return ProcessEvent(&rawStream);
 478 }
 479
 480 //_____________________________________________________________________
 481 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
 482 {
 483   // Method to process=analyze one event in the data stream
 484   if (!in) return kFALSE; //Return right away if there's a null pointer
 485   fNEvents++; // one more event
 486
 487   // indices for the reading
 488   int sample = 0;
 489   int gain = 0;
 490   int time = 0;
 491   int i = 0; // sample counter
 492   int startBin = 0;
 493
 494   // start loop over input stream
 495   while (in->NextDDL()) {
 496     while (in->NextChannel()) {
 497
 498       // counters
 499       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
 500       int nsamples = 0;
 501
 502       // for the pedestal calculation
 503       int sampleSum = 0; // sum of samples
 504       int squaredSampleSum = 0; // sum of samples squared
 505       int nSum = 0; // number of samples in sum
 506       // calc. quantities
 507       double mean = 0, squaredMean = 0, rms = 0;
 508
 509       while (in->NextBunch()) {
 510         const UShort_t *sig = in->GetSignals();
 511         startBin = in->GetStartTimeBin();
 512         nsamples += in->GetBunchLength();
 513         for (i = 0; i < in->GetBunchLength(); i++) {
 514           sample = sig[i];
 515           time = startBin--;
 516
 517           // check if it's a min or max value
 518           if (sample < min) min = sample;
 519           if (sample > max) max = sample;
 520
 521           // should we add it for the pedestal calculation?
 522           if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
 523                !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
 524             sampleSum += sample;
 525             squaredSampleSum += sample*sample;
 526             nSum++;
 527           }
 528
 529         } // loop over samples in bunch
 530       } // loop over bunches
 531
 532       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
 533
 534       // calculate pedesstal estimate: mean of possibly selected samples
 535       if (nSum > 0) {
 536         mean = sampleSum / (1.0 * nSum);
 537         squaredMean = squaredSampleSum / (1.0 * nSum);
 538         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
 539         rms = sqrt(squaredMean - mean*mean);
 540       }
 541       else {
 542         mean = 0;
 543         squaredMean = 0;
 544         rms  = 0;
 545       }
 546
 547       // we're done with the calc. for this channel; let's prepare to fill histo
 548       gain = -1; // init to not valid value
 549       if ( in->IsLowGain() ) {
 550         gain = 0;
 551       }
 552       else if ( in->IsHighGain() ) {
 553         gain = 1;
 554       }
 555
 556       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
 557       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 558       if (arrayPos >= fModules) {
 559         //TODO: return an error message, if appopriate (perhaps if debug>0?)
 560         return kFALSE;
 561       }
 562       //Debug
 563       if (arrayPos < 0 || arrayPos >= fModules) {
 564         printf("Oh no: arrayPos = %i.\n", arrayPos);
 565       }
 566
 567       fNChanFills++; // one more channel found, and profile to be filled
 568       //NOTE: coordinates are (column, row) for the profiles
 569       if (gain == 0) {
 570         //fill the low gain histograms
 571         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 572         if (nSum>0) { // only fill pedestal info in case it could be calculated
 573           ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 574           ((TProfile2D*)fPedestalRMSLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 575         }
 576       }
 577       else if (gain == 1) {
 578         //fill the high gain ones
 579         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 580         if (nSum>0) { // only fill pedestal info in case it could be calculated
 581           ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 582           ((TProfile2D*)fPedestalRMSHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 583         }
 584         // for warning checks
 585         int idx = in->GetRow() + fRows * in->GetColumn();
 586         ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min);
 587       }//end if valid gain
 588
 589       } // nsamples>0 check, some data found for this channel; not only trailer/header
 590     }// end while over channel
 591   }//end while over DDL's, of input stream
 592
 593   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 594
 595   return kTRUE;
 596 }
 597
 598 //_____________________________________________________________________
 599 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
 600 {
 601   //Saves all the histograms (or profiles, to be accurate) to the designated file
 602
 603   TFile destFile(fileName, "recreate");
 604
 605   if (destFile.IsZombie()) {
 606     return kFALSE;
 607   }
 608
 609   destFile.cd();
 610
 611   for (int i = 0; i < fModules; i++) {
 612     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
 613       fPeakMinusPedLowGain[i]->Write();
 614     }
 615     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
 616       fPeakMinusPedHighGain[i]->Write();
 617     }
 618     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
 619       fPedestalLowGain[i]->Write();
 620     }
 621     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
 622       fPedestalHighGain[i]->Write();
 623     }
 624     if( ((TProfile2D *)fPedestalRMSLowGain[i])->GetEntries() || saveEmptyHistos) {
 625       fPedestalRMSLowGain[i]->Write();
 626     }
 627     if( ((TProfile2D *)fPedestalRMSHighGain[i])->GetEntries() || saveEmptyHistos) {
 628       fPedestalRMSHighGain[i]->Write();
 629     }
 630     if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) {
 631       fPeakMinusPedHighGainHisto[i]->Write();
 632     }
 633
 634   }
 635
 636   destFile.Close();
 637
 638   return kTRUE;
 639 }
 640
 641 //_____________________________________________________________________
 642 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
 643 {
 644
 645   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
 646   TH1::AddDirectory(kFALSE);
 647
 648   TFile *sourceFile = new TFile(fileName);
 649   if (sourceFile->IsZombie()) {
 650     return kFALSE;//We couldn't load the reference
 651   }
 652
 653   if (fReference) delete fReference;//Delete the reference object, if it already exists
 654   fReference = 0;
 655
 656   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
 657
 658   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
 659     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
 660     fReference = 0;
 661     return kFALSE;
 662   }
 663
 664   delete sourceFile;
 665
 666   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
 667   //if we are called by someone who has set it to false...
 668   TH1::AddDirectory(kTRUE);
 669
 670   return kTRUE;//We succesfully loaded the object
 671 }
 672
 673 //_____________________________________________________________________
 674 void AliCaloCalibPedestal::ValidateComparisonProfiles()
 675 {
 676   //Make sure the comparison histos exist
 677   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
 678   //the same time
 679
 680
 681   //Then, loop for the requested number of modules
 682   TString title, name;
 683   for (int i = 0; i < fModules; i++) {
 684     //Pedestals, low gain
 685     name = "hPedlowgainDiff";
 686     name += i;
 687     title = "Pedestals difference, low gain, module ";
 688     title += i;
 689     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
 690                                             fColumns, 0.0, fColumns,
 691                                             fRows, fRowMin, fRowMax,"s"));
 692
 693     //Pedestals, high gain
 694     name = "hPedhighgainDiff";
 695     name += i;
 696     title = "Pedestals difference, high gain, module ";
 697     title += i;
 698     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
 699                                              fColumns, 0.0, fColumns,
 700                                              fRows, fRowMin, fRowMax,"s"));
 701
 702     //Peak-Pedestals, high gain
 703     name = "hPeakMinusPedhighgainDiff";
 704     name += i;
 705     title = "Peak-Pedestal difference, high gain, module ";
 706     title += i;
 707     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
 708                                                  fColumns, 0.0, fColumns,
 709                                                  fRows, fRowMin, fRowMax,"s"));
 710
 711     //Pedestals, low gain
 712     name = "hPedlowgainRatio";
 713     name += i;
 714     title = "Pedestals ratio, low gain, module ";
 715     title += i;
 716     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
 717                                              fColumns, 0.0, fColumns,
 718                                              fRows, fRowMin, fRowMax,"s"));
 719
 720     //Pedestals, high gain
 721     name = "hPedhighgainRatio";
 722     name += i;
 723     title = "Pedestals ratio, high gain, module ";
 724     title += i;
 725     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
 726                                               fColumns, 0.0, fColumns,
 727                                               fRows, fRowMin, fRowMax,"s"));
 728
 729     //Peak-Pedestals, low gain
 730     name = "hPeakMinusPedlowgainRatio";
 731     name += i;
 732     title = "Peak-Pedestal ratio, low gain, module ";
 733     title += i;
 734     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
 735                                                  fColumns, 0.0, fColumns,
 736                                                  fRows, fRowMin, fRowMax,"s"));
 737
 738     //Peak-Pedestals, high gain
 739     name = "hPeakMinusPedhighgainRatio";
 740     name += i;
 741     title = "Peak-Pedestal ratio, high gain, module ";
 742     title += i;
 743     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
 744                                                   fColumns, 0.0, fColumns,
 745                                                   fRows, fRowMin, fRowMax,"s"));
 746
 747   }//end for nModules create the histograms
 748 }
 749
 750 //_____________________________________________________________________
 751 void AliCaloCalibPedestal::ComputeDiffAndRatio()
 752 {
 753   // calculate differences and ratios relative to a reference
 754   ValidateComparisonProfiles();//Make sure the comparison histos exist
 755
 756   if (!fReference) {
 757     return;//Return if the reference object isn't loaded
 758   }
 759
 760   for (int i = 0; i < fModules; i++) {
 761     //Compute the ratio of the histograms
 762
 763     ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
 764     ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
 765     ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
 766     ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
 767
 768     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
 769     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
 770     for (int j = 0; j <= fColumns; j++) {
 771       for (int k = 0; k <= fRows; k++) {
 772         int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
 773         double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
 774         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 775         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
 776
 777         diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
 778         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 779         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
 780
 781         diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
 782         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 783         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
 784
 785         diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
 786         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 787         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
 788
 789       } // rows
 790     } // columns
 791
 792   } // modules
 793
 794 }
 795
 796 //_____________________________________________________________________
 797 void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile)
 798 { // look for hot/noisy towers
 799   ofstream * fout = 0;
 800   char name[512];//Quite a long temp buffer, just in case the filename includes a path
 801
 802   if (hotMapFile) {
 803     snprintf(name, 512, "%s.txt", hotMapFile);
 804     fout = new ofstream(name);
 805     if (!fout->is_open()) {
 806       delete fout;
 807       fout = 0;//Set the pointer to empty if the file was not opened
 808     }
 809   }
 810
 811   for(int i = 0; i < fModules; i++){
 812
 813     //first we compute the peak-pedestal distribution for each supermodule...
 814     if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) {
 815       double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin());
 816       double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin());
 817       TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1);
 818
 819       for (int j = 1; j <= fColumns; j++) {
 820         for (int k = 1; k <= fRows; k++) {
 821           hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k));
 822         }
 823       }
 824
 825       //...and then we fit it
 826       hPeakFit->Fit("gaus", "OQ", "",  hPeakFit->GetMean() - 3*hPeakFit->GetRMS(),
 827                     hPeakFit->GetMean() + 3*hPeakFit->GetRMS());
 828       //hPeakFit->Draw();
 829       double mean  = hPeakFit->GetFunction("gaus")->GetParameter(1);
 830       double sigma = hPeakFit->GetFunction("gaus")->GetParameter(2);
 831
 832       delete hPeakFit;
 833
 834       //Then we look for warm/hot towers
 835       TH2F * hPeak2D = GetPeakHighGainHisto(i);
 836       hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) );
 837
 838       int idx = 0 ;
 839       int warnCounter = 0;
 840       for (int j = 1; j <= fColumns; j++) {
 841         for (int k = 1; k <= fRows; k++) {
 842           //we start looking for warm/warning towers...
 843           // histogram x-axis index
 844           idx = k-1 + fRows*(j-1); // this is what is used in the Fill call
 845           hPeak2D->GetXaxis()->SetRangeUser(idx, idx);
 846           warnCounter = (int) hPeak2D->Integral();
 847           if(warnCounter > fNEvents * fWarningFraction) {
 848             ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning);
 849             /* printf("mod %d col %d row %d warnCounter %d - status %d\n",
 850                i, j-1, k-1, warnCounter, (int) (kWarning)); */
 851           }
 852           //...then we look for hot ones (towers whose values are greater than mean + X*sigma)
 853           if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) {
 854             ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot);
 855             /* printf("mod %d col %d row %d  binc %d - status %d\n",
 856                i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */
 857           }
 858
 859           //Write the status to the hot/warm map file, if the file is open.
 860           // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum)
 861           if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) {
 862
 863             (*fout) << i << " "
 864                     << (j - 1) << " "
 865                     << (k - 1) << " "
 866                     << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl;               }
 867
 868         }
 869       }
 870
 871     }
 872   }
 873   return;
 874 }
 875
 876 //_____________________________________________________________________
 877 void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile)
 878 {
 879   //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
 880   int countTot = 0;
 881   int countNew = 0;
 882   int countRes = 0;
 883   ofstream * fout = 0;
 884   ofstream * diff = 0;
 885   char name[512];//Quite a long temp buffer, just in case the filename includes a path
 886
 887   if (deadMapFile) {
 888     snprintf(name, 512, "%s.txt", deadMapFile);
 889     fout = new ofstream(name);
 890     snprintf(name, 512, "%sdiff.txt", deadMapFile);
 891     diff = new ofstream(name);
 892     if (!fout->is_open()) {
 893       delete fout;
 894       fout = 0;//Set the pointer to empty if the file was not opened
 895     }
 896     if (!diff->is_open()) {
 897       delete diff;
 898       fout = 0;//Set the pointer to empty if the file was not opened
 899     }
 900   }
 901
 902   for (int i = 0; i < fModules; i++) {
 903     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
 904       for (int j = 1; j <= fColumns; j++) {
 905         for (int k = 1; k <= fRows; k++) {
 906
 907           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead
 908             countTot++;//One more dead total
 909             if (fout) {
 910               (*fout) << i << " "
 911                       << (j - 1) << " "
 912                       << (k - 1) << " "
 913                       << "1" << " "
 914                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
 915             }
 916
 917             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) {
 918               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
 919               countNew++;//This tower wasn't dead before!
 920               if (diff) {
 921                 ( *diff) << i << " "
 922                          << (j - 1) << " "
 923                          << (k - 1) << " "
 924                          << "1" << " "
 925                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
 926               }
 927             }
 928             else {
 929               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
 930             }
 931           }
 932           else { //It's ALIVE!!
 933             //Don't bother with writing the live ones.
 934             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {
 935               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
 936               countRes++; //This tower was dead before => it's a miracle! :P
 937               if (diff) {
 938                 (*diff) << i << " "
 939                         << (j - 1) << " "
 940                         << (k - 1) << " "
 941                         << "1" << " "
 942                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
 943               }
 944             }
 945             else {
 946               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
 947             }
 948           }
 949
 950         }//end for k/rows
 951       }//end for j/columns
 952     }//end if GetEntries >= 0
 953
 954   }//end for modules
 955
 956  if (fout) {
 957    fout->close();
 958    delete fout;
 959  }
 960
 961  fDeadTowers = countTot;
 962  fNewDeadTowers = countNew;
 963  fResurrectedTowers = countRes;
 964 }
 965
 966 //_____________________________________________________________________
 967 Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const
 968 {
 969   //Check if channel is dead or hot.
 970   Int_t status =  (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) );
 971   if(status == kAlive)
 972     return kFALSE;
 973   else
 974     return kTRUE;
 975
 976 }
 977
 978 //_____________________________________________________________________
 979 void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status)
 980 {
 981   //Set status of channel dead, hot, alive ...
 982   ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status);
 983 }