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