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