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