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