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 "AliCaloRawStream.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 = fgkEmCalCols;
 105     fRows = fgkEmCalRows;
 106     fModules = 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   AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 332   return ProcessEvent(&rawStream);
 333 }
 334
 335 //_____________________________________________________________________
 336 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *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   // counters
 347   int i = 0; // the sample number in current event.
 348   int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
 349
 350   // for the pedestal calculation
 351   int sampleSum = 0; // sum of samples
 352   int squaredSampleSum = 0; // sum of samples squared
 353   int nSum = 0; // number of samples in sum
 354   // calc. quantities
 355   double mean = 0, squaredMean = 0, rms = 0;
 356
 357   while (in->Next()) {
 358     sample = in->GetSignal(); //Get the adc signal
 359     if (sample < min) min = sample;
 360     if (sample > max) max = sample;
 361     i++;
 362
 363     // should we add it for the pedestal calculation?
 364     time = in->GetTime();
 365     if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
 366          !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
 367       sampleSum += sample;
 368       squaredSampleSum += sample*sample;
 369       nSum++;
 370     }
 371
 372     if ( i >= in->GetTimeLength()) {
 373       // calculate pedesstal estimate: mean of possibly selected samples
 374       if (nSum > 0) {
 375         mean = sampleSum / (1.0 * nSum);
 376         squaredMean = squaredSampleSum / (1.0 * nSum);
 377         // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
 378         rms = sqrt(squaredMean - mean*mean);
 379       }
 380       else {
 381         mean = 0;
 382         squaredMean = 0;
 383         rms  = 0;
 384       }
 385
 386       //If we're here then we're done with this tower
 387       gain = -1; // init to not valid value
 388       if ( in->IsLowGain() ) {
 389         gain = 0;
 390       }
 391       else if ( in->IsHighGain() ) {
 392         gain = 1;
 393       }
 394
 395       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 396       if (arrayPos >= fModules) {
 397         //TODO: return an error message, if appopriate (perhaps if debug>0?)
 398         return kFALSE;
 399       }
 400
 401       //Debug
 402       if (arrayPos < 0 || arrayPos >= fModules) {
 403         printf("Oh no: arrayPos = %i.\n", arrayPos);
 404       }
 405
 406       fNChanFills++; // one more channel found, and profile to be filled
 407       //NOTE: coordinates are (column, row) for the profiles
 408       if (gain == 0) {
 409         //fill the low gain histograms
 410         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 411         if (nSum>0) { // only fill pedestal info in case it could be calculated
 412           ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 413           ((TProfile2D*)fPedestalRMSLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 414         }
 415       }
 416       else if (gain == 1) {
 417         //fill the high gain ones
 418         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
 419         if (nSum>0) { // only fill pedestal info in case it could be calculated
 420           ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), mean);
 421           ((TProfile2D*)fPedestalRMSHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), rms);
 422         }
 423       }//end if gain
 424
 425       // reset counters
 426       max = fgkSampleMin; min = fgkSampleMax;
 427       i = 0;
 428       // also pedestal calc counters
 429       sampleSum = 0; // sum of samples
 430       squaredSampleSum = 0; // sum of samples squared
 431       nSum = 0; // number of samples in sum
 432
 433     }//End if end of tower
 434
 435   }//end while, of stream
 436
 437   return kTRUE;
 438 }
 439
 440 //_____________________________________________________________________
 441 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
 442 {
 443   //Saves all the histograms (or profiles, to be accurate) to the designated file
 444
 445   TFile destFile(fileName, "recreate");
 446
 447   if (destFile.IsZombie()) {
 448     return kFALSE;
 449   }
 450
 451   destFile.cd();
 452
 453   for (int i = 0; i < fModules; i++) {
 454     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
 455       fPeakMinusPedLowGain[i]->Write();
 456     }
 457     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
 458       fPeakMinusPedHighGain[i]->Write();
 459     }
 460     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
 461       fPedestalLowGain[i]->Write();
 462     }
 463     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
 464       fPedestalHighGain[i]->Write();
 465       Printf("save %d", i);
 466     }
 467     if( ((TProfile2D *)fPedestalRMSLowGain[i])->GetEntries() || saveEmptyHistos) {
 468       fPedestalRMSLowGain[i]->Write();
 469     }
 470     if( ((TProfile2D *)fPedestalRMSHighGain[i])->GetEntries() || saveEmptyHistos) {
 471       fPedestalRMSHighGain[i]->Write();
 472     }
 473   }
 474
 475   destFile.Close();
 476
 477   return kTRUE;
 478 }
 479
 480 //_____________________________________________________________________
 481 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
 482 {
 483
 484   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
 485   TH1::AddDirectory(kFALSE);
 486
 487   TFile *sourceFile = new TFile(fileName);
 488   if (sourceFile->IsZombie()) {
 489     return kFALSE;//We couldn't load the reference
 490   }
 491
 492   if (fReference) delete fReference;//Delete the reference object, if it already exists
 493   fReference = 0;
 494
 495   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
 496
 497   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
 498     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
 499     fReference = 0;
 500     return kFALSE;
 501   }
 502
 503   delete sourceFile;
 504
 505   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
 506   //if we are called by someone who has set it to false...
 507   TH1::AddDirectory(kTRUE);
 508
 509   return kTRUE;//We succesfully loaded the object
 510 }
 511
 512 //_____________________________________________________________________
 513 void AliCaloCalibPedestal::ValidateComparisonProfiles()
 514 {
 515   //Make sure the comparison histos exist
 516   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
 517   //the same time
 518
 519
 520   //Then, loop for the requested number of modules
 521   TString title, name;
 522   for (int i = 0; i < fModules; i++) {
 523     //Pedestals, low gain
 524     name = "hPedlowgainDiff";
 525     name += i;
 526     title = "Pedestals difference, low gain, module ";
 527     title += i;
 528     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
 529                                             fColumns, 0.0, fColumns,
 530                                             fRows, fRowMin, fRowMax,"s"));
 531
 532     //Pedestals, high gain
 533     name = "hPedhighgainDiff";
 534     name += i;
 535     title = "Pedestals difference, high gain, module ";
 536     title += i;
 537     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
 538                                              fColumns, 0.0, fColumns,
 539                                              fRows, fRowMin, fRowMax,"s"));
 540
 541     //Peak-Pedestals, high gain
 542     name = "hPeakMinusPedhighgainDiff";
 543     name += i;
 544     title = "Peak-Pedestal difference, high gain, module ";
 545     title += i;
 546     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
 547                                                  fColumns, 0.0, fColumns,
 548                                                  fRows, fRowMin, fRowMax,"s"));
 549
 550     //Pedestals, low gain
 551     name = "hPedlowgainRatio";
 552     name += i;
 553     title = "Pedestals ratio, low gain, module ";
 554     title += i;
 555     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
 556                                              fColumns, 0.0, fColumns,
 557                                              fRows, fRowMin, fRowMax,"s"));
 558
 559     //Pedestals, high gain
 560     name = "hPedhighgainRatio";
 561     name += i;
 562     title = "Pedestals ratio, high gain, module ";
 563     title += i;
 564     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
 565                                               fColumns, 0.0, fColumns,
 566                                               fRows, fRowMin, fRowMax,"s"));
 567
 568     //Peak-Pedestals, low gain
 569     name = "hPeakMinusPedlowgainRatio";
 570     name += i;
 571     title = "Peak-Pedestal ratio, low gain, module ";
 572     title += i;
 573     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
 574                                                  fColumns, 0.0, fColumns,
 575                                                  fRows, fRowMin, fRowMax,"s"));
 576
 577     //Peak-Pedestals, high gain
 578     name = "hPeakMinusPedhighgainRatio";
 579     name += i;
 580     title = "Peak-Pedestal ratio, high gain, module ";
 581     title += i;
 582     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
 583                                                   fColumns, 0.0, fColumns,
 584                                                   fRows, fRowMin, fRowMax,"s"));
 585
 586   }//end for nModules create the histograms
 587 }
 588
 589 //_____________________________________________________________________
 590 void AliCaloCalibPedestal::ComputeDiffAndRatio()
 591 {
 592   // calculate differences and ratios relative to a reference
 593   ValidateComparisonProfiles();//Make sure the comparison histos exist
 594
 595   if (!fReference) {
 596     return;//Return if the reference object isn't loaded
 597   }
 598
 599   for (int i = 0; i < fModules; i++) {
 600     //Compute the ratio of the histograms
 601
 602     ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
 603     ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
 604     ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
 605     ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
 606
 607     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
 608     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
 609     for (int j = 0; j <= fColumns; j++) {
 610       for (int k = 0; k <= fRows; k++) {
 611         int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
 612         double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
 613         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 614         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
 615
 616         diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
 617         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 618         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
 619
 620         diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
 621         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 622         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
 623
 624         diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
 625         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 626         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
 627
 628       } // rows
 629     } // columns
 630
 631   } // modules
 632
 633 }
 634
 635 //_____________________________________________________________________
 636 void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
 637 {
 638   //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
 639   int countTot = 0;
 640   int countNew = 0;
 641   int countRes = 0;
 642   ofstream * fout = 0;
 643   ofstream * diff = 0;
 644   char name[512];//Quite a long temp buffer, just in case the filename includes a path
 645
 646   if (deadMapFile) {
 647     snprintf(name, 512, "%s.txt", deadMapFile);
 648     fout = new ofstream(name);
 649     snprintf(name, 512, "%sdiff.txt", deadMapFile);
 650     diff = new ofstream(name);
 651     if (!fout->is_open()) {
 652       delete fout;
 653       fout = 0;//Set the pointer to empty if the file was not opened
 654     }
 655     if (!diff->is_open()) {
 656       delete diff;
 657       fout = 0;//Set the pointer to empty if the file was not opened
 658     }
 659   }
 660
 661   for (int i = 0; i < fModules; i++) {
 662     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
 663       for (int j = 1; j <= fColumns; j++) {
 664         for (int k = 1; k <= fRows; k++) {
 665
 666           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {//It's dead
 667             countTot++;//One more dead total
 668             if (fout) {
 669               (*fout) << i << " "
 670                       << (fRows - k) << " "
 671                       << (j-1) << " "
 672                       << "1" << " "
 673                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
 674             }
 675
 676             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= threshold) {
 677               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
 678               countNew++;//This tower wasn't dead before!
 679               if (diff) {
 680                 ( *diff) << i << " "
 681                          << (fRows - k) << " "
 682                          << (j - 1) << " "
 683                          << "1" << " "
 684                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
 685               }
 686             }
 687             else {
 688               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
 689             }
 690           }
 691           else { //It's ALIVE!!
 692             //Don't bother with writing the live ones.
 693             //if (fout)
 694             //  (*fout) << i << " "
 695             //     << (fRows - k) << " "
 696             //     << (j - 1) << " "
 697             //     << "1" << " "
 698             //     << "1" << endl;//Write the status to the deadmap file, if the file is open.
 699             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {
 700               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
 701               countRes++; //This tower was dead before => it's a miracle! :P
 702               if (diff) {
 703                 (*diff) << i << " "
 704                         << (fRows - k) << " "
 705                         << (j - 1) << " "
 706                         << "1" << " "
 707                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
 708               }
 709             }
 710             else {
 711               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
 712             }
 713           }
 714
 715         }//end for k/rows
 716       }//end for j/columns
 717     }//end if GetEntries >= 0
 718
 719   }//end for modules
 720
 721  if (fout) {
 722    fout->close();
 723    delete fout;
 724  }
 725
 726  fDeadTowers = countTot;
 727  fNewDeadTowers = countNew;
 728  fResurrectedTowers = countRes;
 729 }
 730