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   static const string delimitor("::");
 316
 317   // open, check input file
 318   ifstream in( parameterFile );
 319   if( !in ) {
 320     printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
 321     return;
 322   }
 323
 324   // Note: this method is a bit more complicated than it really has to be
 325   // - allowing for multiple entries per line, arbitrary order of the
 326   // different variables etc. But I wanted to try and do this in as
 327   // correct a C++ way as I could (as an exercise).
 328
 329   // read in
 330   char readline[1024];
 331   while ((in.rdstate() & ios::failbit) == 0 ) {
 332
 333     // Read into the raw char array and then construct a string
 334     // to do the searching
 335     in.getline(readline, 1024);
 336     istringstream s(readline);
 337
 338     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
 339
 340       string keyValue;
 341       s >> keyValue;
 342
 343       // check stream status
 344       if( s.rdstate() & ios::failbit ) break;
 345
 346       // skip rest of line if comments found
 347       if( keyValue.substr( 0, 2 ) == "//" ) break;
 348
 349       // look for "::" in keyValue pair
 350       size_t position = keyValue.find( delimitor );
 351       if( position == string::npos ) {
 352         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 353       }
 354
 355       // split keyValue pair
 356       string key( keyValue.substr( 0, position ) );
 357       string value( keyValue.substr( position+delimitor.size(),
 358                                       keyValue.size()-delimitor.size() ) );
 359
 360       // check value does not contain a new delimitor
 361       if( value.find( delimitor ) != string::npos ) {
 362         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 363       }
 364
 365       // debug: check key value pair
 366       // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 367
 368       // if the key matches with something we expect, we assign the new value
 369       istringstream iss(value);
 370       // the comparison strings defined at the beginning of this method
 371       if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") ) {
 372         printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 373
 374         if (key == "fFirstPedestalSample") {
 375           iss >> fFirstPedestalSample;
 376         }
 377         else if (key == "fLastPedestalSample") {
 378           iss >> fLastPedestalSample;
 379         }
 380       } // some match
 381
 382     }
 383   }
 384
 385   in.close();
 386   return;
 387
 388 }
 389
 390 //_____________________________________________________________________
 391 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
 392 {
 393   static const string delimitor("::");
 394   ofstream out( parameterFile );
 395   out << "// " << parameterFile << endl;
 396   out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
 397   out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
 398
 399   out.close();
 400   return;
 401 }
 402
 403 //_____________________________________________________________________
 404 Bool_t AliCaloCalibPedestal::AddInfo(const AliCaloCalibPedestal *ped)
 405 {
 406   // just do this for the basic histograms/profiles that get filled in ProcessEvent
 407   // may not have data for all modules, but let's just Add everything..
 408   for (int i = 0; i < fModules; i++) {
 409     GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
 410     GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
 411     GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
 412     GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
 413   }//end for nModules
 414
 415   // DeadMap; Diff profiles etc would need to be redone after this operation
 416
 417   return kTRUE;//We succesfully added info from the supplied object
 418 }
 419
 420 //_____________________________________________________________________
 421 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
 422 {
 423   // if fMapping is NULL the rawstream will crate its own mapping
 424   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 425   return ProcessEvent(&rawStream);
 426 }
 427
 428 //_____________________________________________________________________
 429 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
 430 {
 431   // Method to process=analyze one event in the data stream
 432   if (!in) return kFALSE; //Return right away if there's a null pointer
 433   fNEvents++; // one more event
 434
 435   // indices for the reading
 436   int sample = 0;
 437   int gain = 0;
 438   int time = 0;
 439   int i = 0; // sample counter
 440   int startBin = 0;
 441
 442   // start loop over input stream
 443   while (in->NextDDL()) {
 444     while (in->NextChannel()) {
 445
 446       // counters
 447       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
 448       int nsamples = 0;
 449
 450       // for the pedestal calculation
 451       int sampleSum = 0; // sum of samples
 452       int squaredSampleSum = 0; // sum of samples squared
 453       int nSum = 0; // number of samples in sum
 454       // calc. quantities
 455       double mean = 0, squaredMean = 0, rms = 0;
 456
 457       while (in->NextBunch()) {
 458         const UShort_t *sig = in->GetSignals();
 459         startBin = in->GetStartTimeBin();
 460         nsamples += in->GetBunchLength();
 461         for (i = 0; i < in->GetBunchLength(); i++) {
 462           sample = sig[i];
 463           time = startBin--;
 464
 465           // check if it's a min or max value
 466           if (sample < min) min = sample;
 467           if (sample > max) max = sample;
 468
 469           // should we add it for the pedestal calculation?
 470           if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
 471                !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
 472             sampleSum += sample;
 473             squaredSampleSum += sample*sample;
 474             nSum++;
 475           }
 476
 477         } // loop over samples in bunch
 478       } // loop over bunches
 479
 480       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
 481
 482       // calculate pedesstal estimate: mean of possibly selected samples
 483       if (nSum > 0) {
 484         mean = sampleSum / (1.0 * nSum);
 485         squaredMean = squaredSampleSum / (1.0 * nSum);
 486         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
 487         rms = sqrt(squaredMean - mean*mean);
 488       }
 489       else {
 490         mean = 0;
 491         squaredMean = 0;
 492         rms  = 0;
 493       }
 494
 495       // we're done with the calc. for this channel; let's prepare to fill histo
 496       gain = -1; // init to not valid value
 497       if ( in->IsLowGain() ) {
 498         gain = 0;
 499       }
 500       else if ( in->IsHighGain() ) {
 501         gain = 1;
 502       }
 503
 504       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
 505       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 506       if (arrayPos >= fModules) {
 507         //TODO: return an error message, if appopriate (perhaps if debug>0?)
 508         return kFALSE;
 509       }
 510       //Debug
 511       if (arrayPos < 0 || arrayPos >= fModules) {
 512         printf("Oh no: arrayPos = %i.\n", arrayPos);
 513       }
 514
 515       fNChanFills++; // one more channel found, and profile to be filled
 516       //NOTE: coordinates are (column, row) for the profiles
 517       if (gain == 0) {
 518         //fill the low gain histograms
 519         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 520         if (nSum>0) { // only fill pedestal info in case it could be calculated
 521           ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 522           ((TProfile2D*)fPedestalRMSLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 523         }
 524       }
 525       else if (gain == 1) {
 526         //fill the high gain ones
 527         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 528         if (nSum>0) { // only fill pedestal info in case it could be calculated
 529           ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 530           ((TProfile2D*)fPedestalRMSHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 531         }
 532       }//end if valid gain
 533
 534       } // nsamples>0 check, some data found for this channel; not only trailer/header
 535     }// end while over channel
 536   }//end while over DDL's, of input stream
 537
 538   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 539
 540   return kTRUE;
 541 }
 542
 543 //_____________________________________________________________________
 544 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
 545 {
 546   //Saves all the histograms (or profiles, to be accurate) to the designated file
 547
 548   TFile destFile(fileName, "recreate");
 549
 550   if (destFile.IsZombie()) {
 551     return kFALSE;
 552   }
 553
 554   destFile.cd();
 555
 556   for (int i = 0; i < fModules; i++) {
 557     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
 558       fPeakMinusPedLowGain[i]->Write();
 559     }
 560     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
 561       fPeakMinusPedHighGain[i]->Write();
 562     }
 563     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
 564       fPedestalLowGain[i]->Write();
 565     }
 566     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
 567       fPedestalHighGain[i]->Write();
 568       Printf("save %d", i);
 569     }
 570     if( ((TProfile2D *)fPedestalRMSLowGain[i])->GetEntries() || saveEmptyHistos) {
 571       fPedestalRMSLowGain[i]->Write();
 572     }
 573     if( ((TProfile2D *)fPedestalRMSHighGain[i])->GetEntries() || saveEmptyHistos) {
 574       fPedestalRMSHighGain[i]->Write();
 575     }
 576   }
 577
 578   destFile.Close();
 579
 580   return kTRUE;
 581 }
 582
 583 //_____________________________________________________________________
 584 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
 585 {
 586
 587   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
 588   TH1::AddDirectory(kFALSE);
 589
 590   TFile *sourceFile = new TFile(fileName);
 591   if (sourceFile->IsZombie()) {
 592     return kFALSE;//We couldn't load the reference
 593   }
 594
 595   if (fReference) delete fReference;//Delete the reference object, if it already exists
 596   fReference = 0;
 597
 598   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
 599
 600   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
 601     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
 602     fReference = 0;
 603     return kFALSE;
 604   }
 605
 606   delete sourceFile;
 607
 608   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
 609   //if we are called by someone who has set it to false...
 610   TH1::AddDirectory(kTRUE);
 611
 612   return kTRUE;//We succesfully loaded the object
 613 }
 614
 615 //_____________________________________________________________________
 616 void AliCaloCalibPedestal::ValidateComparisonProfiles()
 617 {
 618   //Make sure the comparison histos exist
 619   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
 620   //the same time
 621
 622
 623   //Then, loop for the requested number of modules
 624   TString title, name;
 625   for (int i = 0; i < fModules; i++) {
 626     //Pedestals, low gain
 627     name = "hPedlowgainDiff";
 628     name += i;
 629     title = "Pedestals difference, low gain, module ";
 630     title += i;
 631     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
 632                                             fColumns, 0.0, fColumns,
 633                                             fRows, fRowMin, fRowMax,"s"));
 634
 635     //Pedestals, high gain
 636     name = "hPedhighgainDiff";
 637     name += i;
 638     title = "Pedestals difference, high gain, module ";
 639     title += i;
 640     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
 641                                              fColumns, 0.0, fColumns,
 642                                              fRows, fRowMin, fRowMax,"s"));
 643
 644     //Peak-Pedestals, high gain
 645     name = "hPeakMinusPedhighgainDiff";
 646     name += i;
 647     title = "Peak-Pedestal difference, high gain, module ";
 648     title += i;
 649     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
 650                                                  fColumns, 0.0, fColumns,
 651                                                  fRows, fRowMin, fRowMax,"s"));
 652
 653     //Pedestals, low gain
 654     name = "hPedlowgainRatio";
 655     name += i;
 656     title = "Pedestals ratio, low gain, module ";
 657     title += i;
 658     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
 659                                              fColumns, 0.0, fColumns,
 660                                              fRows, fRowMin, fRowMax,"s"));
 661
 662     //Pedestals, high gain
 663     name = "hPedhighgainRatio";
 664     name += i;
 665     title = "Pedestals ratio, high gain, module ";
 666     title += i;
 667     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
 668                                               fColumns, 0.0, fColumns,
 669                                               fRows, fRowMin, fRowMax,"s"));
 670
 671     //Peak-Pedestals, low gain
 672     name = "hPeakMinusPedlowgainRatio";
 673     name += i;
 674     title = "Peak-Pedestal ratio, low gain, module ";
 675     title += i;
 676     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
 677                                                  fColumns, 0.0, fColumns,
 678                                                  fRows, fRowMin, fRowMax,"s"));
 679
 680     //Peak-Pedestals, high gain
 681     name = "hPeakMinusPedhighgainRatio";
 682     name += i;
 683     title = "Peak-Pedestal ratio, high gain, module ";
 684     title += i;
 685     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
 686                                                   fColumns, 0.0, fColumns,
 687                                                   fRows, fRowMin, fRowMax,"s"));
 688
 689   }//end for nModules create the histograms
 690 }
 691
 692 //_____________________________________________________________________
 693 void AliCaloCalibPedestal::ComputeDiffAndRatio()
 694 {
 695   // calculate differences and ratios relative to a reference
 696   ValidateComparisonProfiles();//Make sure the comparison histos exist
 697
 698   if (!fReference) {
 699     return;//Return if the reference object isn't loaded
 700   }
 701
 702   for (int i = 0; i < fModules; i++) {
 703     //Compute the ratio of the histograms
 704
 705     ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
 706     ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
 707     ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
 708     ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
 709
 710     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
 711     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
 712     for (int j = 0; j <= fColumns; j++) {
 713       for (int k = 0; k <= fRows; k++) {
 714         int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
 715         double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
 716         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 717         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
 718
 719         diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
 720         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 721         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
 722
 723         diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
 724         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 725         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
 726
 727         diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
 728         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 729         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
 730
 731       } // rows
 732     } // columns
 733
 734   } // modules
 735
 736 }
 737
 738 //_____________________________________________________________________
 739 void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
 740 {
 741   //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
 742   int countTot = 0;
 743   int countNew = 0;
 744   int countRes = 0;
 745   ofstream * fout = 0;
 746   ofstream * diff = 0;
 747   char name[512];//Quite a long temp buffer, just in case the filename includes a path
 748
 749   if (deadMapFile) {
 750     snprintf(name, 512, "%s.txt", deadMapFile);
 751     fout = new ofstream(name);
 752     snprintf(name, 512, "%sdiff.txt", deadMapFile);
 753     diff = new ofstream(name);
 754     if (!fout->is_open()) {
 755       delete fout;
 756       fout = 0;//Set the pointer to empty if the file was not opened
 757     }
 758     if (!diff->is_open()) {
 759       delete diff;
 760       fout = 0;//Set the pointer to empty if the file was not opened
 761     }
 762   }
 763
 764   for (int i = 0; i < fModules; i++) {
 765     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
 766       for (int j = 1; j <= fColumns; j++) {
 767         for (int k = 1; k <= fRows; k++) {
 768
 769           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {//It's dead
 770             countTot++;//One more dead total
 771             if (fout) {
 772               (*fout) << i << " "
 773                       << (fRows - k) << " "
 774                       << (j-1) << " "
 775                       << "1" << " "
 776                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
 777             }
 778
 779             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= threshold) {
 780               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
 781               countNew++;//This tower wasn't dead before!
 782               if (diff) {
 783                 ( *diff) << i << " "
 784                          << (fRows - k) << " "
 785                          << (j - 1) << " "
 786                          << "1" << " "
 787                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
 788               }
 789             }
 790             else {
 791               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
 792             }
 793           }
 794           else { //It's ALIVE!!
 795             //Don't bother with writing the live ones.
 796             //if (fout)
 797             //  (*fout) << i << " "
 798             //     << (fRows - k) << " "
 799             //     << (j - 1) << " "
 800             //     << "1" << " "
 801             //     << "1" << endl;//Write the status to the deadmap file, if the file is open.
 802             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {
 803               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
 804               countRes++; //This tower was dead before => it's a miracle! :P
 805               if (diff) {
 806                 (*diff) << i << " "
 807                         << (fRows - k) << " "
 808                         << (j - 1) << " "
 809                         << "1" << " "
 810                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
 811               }
 812             }
 813             else {
 814               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
 815             }
 816           }
 817
 818         }//end for k/rows
 819       }//end for j/columns
 820     }//end if GetEntries >= 0
 821
 822   }//end for modules
 823
 824  if (fout) {
 825    fout->close();
 826    delete fout;
 827  }
 828
 829  fDeadTowers = countTot;
 830  fNewDeadTowers = countNew;
 831  fResurrectedTowers = countRes;
 832 }
 833