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