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