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