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