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