Bug fix in the constructor (thanks to A.Marin)
[u/mrichter/AliRoot.git] / 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 "TCanvas.h"
38 #include "TH1.h"
39 #include "TF1.h"
40 #include "TFile.h"
41 #include <fstream>
42 #include <sstream>
43 #include <iostream>
44 #include <stdexcept>
45 #include <cmath>
46
47 #include "AliRawReader.h"
48 #include "AliCaloRawStreamV3.h"
49
50 //The include file
51 #include "AliCaloCalibPedestal.h"
52
53 ClassImp(AliCaloCalibPedestal)
54
55 using namespace std;
56
57 // ctor; initialize everything in order to avoid compiler warnings
58 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :  
59   TObject(),
60   fPedestalLowGain(),
61   fPedestalHighGain(),
62   fPedestalLEDRefLowGain(),
63   fPedestalLEDRefHighGain(),
64   fPeakMinusPedLowGain(),
65   fPeakMinusPedHighGain(),
66   fPeakMinusPedHighGainHisto(),
67   fPedestalLowGainDiff(),
68   fPedestalHighGainDiff(),
69   fPedestalLEDRefLowGainDiff(),
70   fPedestalLEDRefHighGainDiff(),
71   fPeakMinusPedLowGainDiff(),
72   fPeakMinusPedHighGainDiff(),
73   fPedestalLowGainRatio(),
74   fPedestalHighGainRatio(),
75   fPedestalLEDRefLowGainRatio(),
76   fPedestalLEDRefHighGainRatio(),
77   fPeakMinusPedLowGainRatio(),
78   fPeakMinusPedHighGainRatio(),
79   fDeadMap(),
80   fNEvents(0),
81   fNChanFills(0),
82   fDeadTowers(0),
83   fNewDeadTowers(0),
84   fResurrectedTowers(0),
85   fReference(0),
86   fDetType(kNone),
87   fColumns(0),
88   fRows(0),
89   fLEDRefs(0),
90   fModules(0),
91   fRowMin(0),
92   fRowMax(0),
93   fRowMultiplier(0),
94   fCaloString(),
95   fMapping(NULL),
96   fRunNumber(-1),
97   fSelectPedestalSamples(kTRUE), 
98   fFirstPedestalSample(0),
99   fLastPedestalSample(15),
100   fDeadThreshold(5),
101   fWarningThreshold(50),
102   fWarningFraction(0.002),
103   fHotSigma(5)
104 {
105   //Default constructor. First we set the detector-type related constants.
106   if (detectorType == kPhos) {
107     fColumns = fgkPhosCols;
108     fRows = fgkPhosRows;
109     fLEDRefs = fgkPhosLEDRefs;
110     fModules = fgkPhosModules;
111     fCaloString = "PHOS";
112     fRowMin = -1*fRows;
113     fRowMax = 0;
114     fRowMultiplier = -1;
115   } 
116   else {
117     //We'll just trust the enum to keep everything in line, so that if detectorType
118     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
119     //case, if someone intentionally gives another number
120     fColumns = AliEMCALGeoParams::fgkEMCALCols;
121     fRows = AliEMCALGeoParams::fgkEMCALRows;
122     fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
123     fModules = AliEMCALGeoParams::fgkEMCALModules;
124     fCaloString = "EMCAL";
125     fRowMin = 0;
126     fRowMax = fRows;
127     fRowMultiplier = 1;
128   } 
129   fDetType = detectorType;
130
131   // ValidateProfiles(); // not to be done in ctor; info from Axel N. 
132 }
133
134 //_____________________________________________________________________
135 void AliCaloCalibPedestal::ValidateProfiles()
136 {
137   //Make sure the basic histos exist
138   if (!fPedestalLowGain.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
139   //the same time
140
141   //Then, loop for the requested number of modules
142   TString title, name;
143   for (int i = 0; i < fModules; i++) {
144     //Pedestals, low gain
145     name = "hPedlowgain";
146     name += i;
147     title = "Pedestals, low gain, module ";
148     title += i; 
149     fPedestalLowGain.Add(new TProfile2D(name, title,
150                                         fColumns, 0.0, fColumns, 
151                                         fRows, fRowMin, fRowMax,"s"));
152   
153     //Pedestals, high gain
154     name = "hPedhighgain";
155     name += i;
156     title = "Pedestals, high gain, module ";
157     title += i; 
158     fPedestalHighGain.Add(new TProfile2D(name, title,
159                                          fColumns, 0.0, fColumns, 
160                                          fRows, fRowMin, fRowMax,"s"));
161
162     //LED Ref/Mon pedestals, low gain
163     name = "hPedestalLEDReflowgain";
164     name += i;
165     title = "Pedestal LEDRef, low gain, module ";
166     title += i; 
167     fPedestalLEDRefLowGain.Add(new TProfile(name, title,
168                                             fLEDRefs, 0.0, fLEDRefs, "s"));
169     
170     //LED Ref/Mon pedestals, high gain
171     name = "hPedestalLEDRefhighgain";
172     name += i;
173     title = "Pedestal LEDRef, high gain, module ";
174     title += i; 
175     fPedestalLEDRefHighGain.Add(new TProfile(name, title,
176                                              fLEDRefs, 0.0, fLEDRefs, "s"));
177   
178     //Peak-Pedestals, low gain
179     name = "hPeakMinusPedlowgain";
180     name += i;
181     title = "Peak-Pedestal, low gain, module ";
182     title += i; 
183     fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
184                                             fColumns, 0.0, fColumns, 
185                                             fRows, fRowMin, fRowMax,"s"));
186   
187     //Peak-Pedestals, high gain
188     name = "hPeakMinusPedhighgain";
189     name += i;
190     title = "Peak-Pedestal, high gain, module ";
191     title += i; 
192     fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
193                                              fColumns, 0.0, fColumns, 
194                                              fRows, fRowMin, fRowMax,"s"));
195
196     //Peak-Pedestals, high gain - TH2F histo
197     name = "hPeakMinusPedhighgainHisto";
198     name += i;
199     title = "Peak-Pedestal, high gain, module ";
200     title += i; 
201     fPeakMinusPedHighGainHisto.Add(new TH2F(name, title,
202                                             fColumns*fRows, 0.0, fColumns*fRows, 
203                                             100, 0, 1000));
204  
205     name = "hDeadMap";
206     name += i;
207     title = "Dead map, module ";
208     title += i;
209     fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns, 
210                           fRows, fRowMin, fRowMax));
211   
212   }//end for nModules create the histograms
213
214   CompressAndSetOwner();
215 }
216
217 //_____________________________________________________________________
218 void AliCaloCalibPedestal::CompressAndSetOwner()
219
220   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
221   fPedestalLowGain.Compress();
222   fPedestalHighGain.Compress();
223   fPedestalLEDRefLowGain.Compress();
224   fPedestalLEDRefHighGain.Compress();
225   fPeakMinusPedLowGain.Compress();
226   fPeakMinusPedHighGain.Compress();
227   fPeakMinusPedHighGainHisto.Compress();
228   fDeadMap.Compress();
229
230   // set owner ship for everyone
231   fPedestalLowGain.SetOwner(kTRUE);
232   fPedestalHighGain.SetOwner(kTRUE);
233   fPedestalLEDRefLowGain.SetOwner(kTRUE);
234   fPedestalLEDRefHighGain.SetOwner(kTRUE);
235   fPeakMinusPedLowGain.SetOwner(kTRUE);
236   fPeakMinusPedHighGain.SetOwner(kTRUE);
237   fPeakMinusPedHighGainHisto.SetOwner(kTRUE);
238   fPedestalLowGainDiff.SetOwner(kTRUE);
239   fPedestalHighGainDiff.SetOwner(kTRUE);
240   fPedestalLEDRefLowGainDiff.SetOwner(kTRUE);
241   fPedestalLEDRefHighGainDiff.SetOwner(kTRUE);
242   fPeakMinusPedLowGainDiff.SetOwner(kTRUE);
243   fPeakMinusPedHighGainDiff.SetOwner(kTRUE);
244   fPedestalLowGainRatio.SetOwner(kTRUE);
245   fPedestalHighGainRatio.SetOwner(kTRUE);
246   fPedestalLEDRefLowGainRatio.SetOwner(kTRUE);
247   fPedestalLEDRefHighGainRatio.SetOwner(kTRUE);
248   fPeakMinusPedLowGainRatio.SetOwner(kTRUE);
249   fPeakMinusPedHighGainRatio.SetOwner(kTRUE);
250   fDeadMap.SetOwner(kTRUE);
251 }
252
253 // dtor
254 //_____________________________________________________________________
255 AliCaloCalibPedestal::~AliCaloCalibPedestal()
256 {
257   //dtor
258   
259   if (fReference) delete fReference;//Delete the reference object, if it has been loaded
260   
261   // delete also TObjArray's 
262   fPedestalLowGain.Delete(); 
263   fPedestalHighGain.Delete();
264   fPedestalLEDRefLowGain.Delete();
265   fPedestalLEDRefHighGain.Delete();
266   fPeakMinusPedLowGain.Delete();
267   fPeakMinusPedHighGain.Delete();
268   fPeakMinusPedHighGainHisto.Delete();
269   fPedestalLowGainDiff.Delete();
270   fPedestalHighGainDiff.Delete();
271   fPedestalLEDRefLowGainDiff.Delete();
272   fPedestalLEDRefHighGainDiff.Delete();
273   fPeakMinusPedLowGainDiff.Delete();
274   fPeakMinusPedHighGainDiff.Delete();
275   fPedestalLowGainRatio.Delete();
276   fPedestalHighGainRatio.Delete();
277   fPedestalLEDRefLowGainRatio.Delete();
278   fPedestalLEDRefHighGainRatio.Delete();
279   fPeakMinusPedLowGainRatio.Delete();
280   fPeakMinusPedHighGainRatio.Delete();
281   fDeadMap.Delete();
282   
283 }
284
285 // copy ctor
286 //_____________________________________________________________________
287 AliCaloCalibPedestal::AliCaloCalibPedestal(AliCaloCalibPedestal &ped) :
288   TObject(ped),
289   fPedestalLowGain(),
290   fPedestalHighGain(),
291   fPedestalLEDRefLowGain(),
292   fPedestalLEDRefHighGain(),
293   fPeakMinusPedLowGain(),
294   fPeakMinusPedHighGain(),
295   fPeakMinusPedHighGainHisto(),
296   fPedestalLowGainDiff(),
297   fPedestalHighGainDiff(),
298   fPedestalLEDRefLowGainDiff(),
299   fPedestalLEDRefHighGainDiff(),
300   fPeakMinusPedLowGainDiff(),
301   fPeakMinusPedHighGainDiff(),
302   fPedestalLowGainRatio(),
303   fPedestalHighGainRatio(),
304   fPedestalLEDRefLowGainRatio(),
305   fPedestalLEDRefHighGainRatio(),
306   fPeakMinusPedLowGainRatio(),
307   fPeakMinusPedHighGainRatio(),
308   fDeadMap(),
309   fNEvents(ped.GetNEvents()),
310   fNChanFills(ped.GetNChanFills()),
311   fDeadTowers(ped.GetDeadTowerCount()),
312   fNewDeadTowers(ped.GetDeadTowerNew()),
313   fResurrectedTowers(ped.GetDeadTowerResurrected()),
314   fReference( 0 ), //! note that we do not try to copy the reference info here
315   fDetType(ped.GetDetectorType()),
316   fColumns(ped.GetColumns()),
317   fRows(ped.GetRows()),
318   fLEDRefs(ped.GetLEDRefs()),
319   fModules(ped.GetModules()),
320   fRowMin(ped.GetRowMin()),
321   fRowMax(ped.GetRowMax()),
322   fRowMultiplier(ped.GetRowMultiplier()),
323   fCaloString(ped.GetCaloString()),
324   fMapping(NULL), //! note that we are not copying the map info
325   fRunNumber(ped.GetRunNumber()),
326   fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
327   fFirstPedestalSample(ped.GetFirstPedestalSample()),
328   fLastPedestalSample(ped.GetLastPedestalSample()),
329   fDeadThreshold(ped.GetDeadThreshold()),
330   fWarningThreshold(ped.GetWarningThreshold()),
331   fWarningFraction(ped.GetWarningFraction()),
332   fHotSigma(ped.GetHotSigma())
333 {
334   // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
335   //DS: this has not really been tested yet..
336   for (int i = 0; i < fModules; i++) {
337     fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
338     fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
339     fPedestalLEDRefLowGain.Add( ped.GetPedLEDRefProfileLowGain(i) );
340     fPedestalLEDRefHighGain.Add( ped.GetPedLEDRefProfileHighGain(i) );
341     fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
342     fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
343     fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) );
344
345     fDeadMap.Add( ped.GetDeadMap(i) );  
346   }//end for nModules 
347  
348   CompressAndSetOwner();
349 }
350
351 // assignment operator; use copy ctor to make life easy..
352 //_____________________________________________________________________
353 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (AliCaloCalibPedestal &source)
354 {
355   // assignment operator; use copy ctor
356   if (&source == this) return *this;
357
358   new (this) AliCaloCalibPedestal(source);
359   return *this;
360 }
361
362 //_____________________________________________________________________
363 void AliCaloCalibPedestal::Reset()
364 {
365   ValidateProfiles(); // make sure histos/profiles exist
366   // Reset all arrays/histograms
367   for (int i = 0; i < fModules; i++) {
368     GetPedProfileLowGain(i)->Reset();
369     GetPedProfileHighGain(i)->Reset();
370     GetPedLEDRefProfileLowGain(i)->Reset();
371     GetPedLEDRefProfileHighGain(i)->Reset();
372     GetPeakProfileLowGain(i)->Reset();
373     GetPeakProfileHighGain(i)->Reset();
374     GetPeakHighGainHisto(i)->Reset();
375     GetDeadMap(i)->Reset();
376     
377     if (!fPedestalLowGainDiff.IsEmpty()) {
378       //This means that the comparison profiles have been created.
379   
380       GetPedProfileLowGainDiff(i)->Reset();
381       GetPedProfileHighGainDiff(i)->Reset();
382       GetPedLEDRefProfileLowGainDiff(i)->Reset();
383       GetPedLEDRefProfileHighGainDiff(i)->Reset();
384       GetPeakProfileLowGainDiff(i)->Reset();
385       GetPeakProfileHighGainDiff(i)->Reset();
386       
387       GetPedProfileLowGainRatio(i)->Reset();
388       GetPedProfileHighGainRatio(i)->Reset();
389       GetPedLEDRefProfileLowGainRatio(i)->Reset();
390       GetPedLEDRefProfileHighGainRatio(i)->Reset();
391       GetPeakProfileLowGainRatio(i)->Reset();
392       GetPeakProfileHighGainRatio(i)->Reset();
393     }
394   }
395   fNEvents = 0;
396   fNChanFills = 0;
397   fDeadTowers = 0;
398   fNewDeadTowers = 0;
399   fResurrectedTowers = 0;
400  
401   //To think about: should fReference be deleted too?... let's not do it this time, at least...
402 }
403
404 // Parameter/cut handling
405 //_____________________________________________________________________
406 void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
407 {  
408   // Note: this method is a bit more complicated than it really has to be
409   // - allowing for multiple entries per line, arbitrary order of the
410   // different variables etc. But I wanted to try and do this in as
411   // correct a C++ way as I could (as an exercise).
412
413   static const string delimitor("::");
414         
415   // open, check input file
416   ifstream in( parameterFile );
417   if( !in ) {
418     printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
419     return;
420   } 
421
422
423   // read in
424   char readline[1024];
425   while ((in.rdstate() & ios::failbit) == 0 ) {
426     
427     // Read into the raw char array and then construct a string
428     // to do the searching
429     in.getline(readline, 1024);
430     istringstream s(readline);          
431                 
432     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
433                         
434       string keyValue; 
435       s >> keyValue;
436       
437       // check stream status
438       if( ( s.rdstate() & ios::failbit ) == ios::failbit) break;
439                         
440       // skip rest of line if comments found
441       if( keyValue.substr( 0, 2 ) == "//" ) break;
442                         
443       // look for "::" in keyValue pair
444       size_t position = keyValue.find( delimitor );
445       if( position == string::npos ) {
446         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
447       }
448                                 
449       // split keyValue pair
450       string key( keyValue.substr( 0, position ) );
451       string value( keyValue.substr( position+delimitor.size(), 
452                                       keyValue.size()-delimitor.size() ) );
453                         
454       // check value does not contain a new delimitor
455       if( value.find( delimitor ) != string::npos ) {
456         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
457       }
458       
459       // debug: check key value pair
460       // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
461
462       // if the key matches with something we expect, we assign the new value
463       istringstream iss(value);
464       // the comparison strings defined at the beginning of this method
465       if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) {
466         printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
467
468         if (key == "fFirstPedestalSample") { 
469           iss >> fFirstPedestalSample; 
470         }
471         else if (key == "fLastPedestalSample") { 
472           iss >> fLastPedestalSample; 
473         }
474         else if (key == "fDeadThreshold") { 
475           iss >> fDeadThreshold; 
476         }
477         else if (key == "fWarningThreshold") { 
478           iss >> fWarningThreshold; 
479         }
480         else if (key == "fWarningFraction") { 
481           iss >> fWarningFraction; 
482         }
483         else if (key == "fHotSigma") { 
484           iss >> fHotSigma; 
485         }
486
487       } // some match
488
489     }           
490   }
491
492   in.close();
493   return;
494         
495 }
496
497 //_____________________________________________________________________
498 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
499 {
500   //Write parameters in file.
501         
502   static const string delimitor("::");
503   ofstream out( parameterFile );
504   out << "// " << parameterFile << endl;
505   out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
506   out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
507   out << "fDeadThreshold" << "::" << fDeadThreshold << endl;
508   out << "fWarningThreshold" << "::" << fWarningThreshold << endl;
509   out << "fWarningFraction" << "::" << fWarningFraction << endl;
510   out << "fHotSigma" << "::" << fHotSigma << endl;
511
512   out.close();
513   return;
514 }
515
516 //_____________________________________________________________________
517 Bool_t AliCaloCalibPedestal::AddInfo(AliCaloCalibPedestal *ped)
518 {
519   ValidateProfiles(); // make sure histos/profiles exist
520   // just do this for the basic histograms/profiles that get filled in ProcessEvent
521   // may not have data for all modules, but let's just Add everything..
522   for (int i = 0; i < fModules; i++) {
523     GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
524     GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
525     GetPedLEDRefProfileLowGain(i)->Add( ped->GetPedLEDRefProfileLowGain(i) );
526     GetPedLEDRefProfileHighGain(i)->Add( ped->GetPedLEDRefProfileHighGain(i) );
527     GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
528     GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
529     GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) );
530
531   }//end for nModules 
532
533   // We should also copy other pieces of info: counters and parameters 
534   // (not number of columns and rows etc which should be the same)
535   // note that I just assign them here rather than Add them, but we
536   // normally just Add (e.g. in Preprocessor) one object so this should be fine.
537   fNEvents = ped->GetNEvents();
538   fNChanFills = ped->GetNChanFills();
539   fDeadTowers = ped->GetDeadTowerCount();
540   fNewDeadTowers = ped->GetDeadTowerNew();
541   fResurrectedTowers = ped->GetDeadTowerResurrected();
542   fRunNumber = ped->GetRunNumber();
543   fSelectPedestalSamples = ped->GetSelectPedestalSamples();
544   fFirstPedestalSample = ped->GetFirstPedestalSample();
545   fLastPedestalSample = ped->GetLastPedestalSample();
546   fDeadThreshold = ped->GetDeadThreshold();
547   fWarningThreshold = ped->GetWarningThreshold();
548   fWarningFraction = ped->GetWarningFraction();
549   fHotSigma = ped->GetHotSigma();
550
551   // DeadMap; Diff profiles etc would need to be redone after this operation
552
553   return kTRUE;//We succesfully added info from the supplied object
554 }
555
556 //_____________________________________________________________________
557 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
558
559   // if fMapping is NULL the rawstream will crate its own mapping
560   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
561   if (fDetType == kEmCal) {
562     rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range 
563   }
564   return ProcessEvent(&rawStream);
565 }
566
567 //_____________________________________________________________________
568 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
569
570   // Method to process=analyze one event in the data stream
571   if (!in) return kFALSE; //Return right away if there's a null pointer
572   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
573
574   fNEvents++; // one more event
575
576   if (fNEvents==1) ValidateProfiles(); // 1st event, make sure histos/profiles exist
577   
578   // indices for the reading
579   int sample = 0;
580   int time = 0;
581   int i = 0; // sample counter
582   int startBin = 0;
583
584   // start loop over input stream 
585   while (in->NextDDL()) {
586     while (in->NextChannel()) {
587
588       // counters
589       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
590       int nsamples = 0;
591
592       // pedestal samples
593       int nPed = 0;
594       vector<int> pedSamples; 
595
596       while (in->NextBunch()) {
597         const UShort_t *sig = in->GetSignals();
598         startBin = in->GetStartTimeBin();
599         nsamples += in->GetBunchLength();
600         for (i = 0; i < in->GetBunchLength(); i++) {
601           sample = sig[i];
602           time = startBin--;
603
604           // check if it's a min or max value
605           if (sample < min) min = sample;
606           if (sample > max) max = sample;
607           
608           // should we add it for the pedestal calculation?
609           if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
610                !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all 
611             pedSamples.push_back( sig[i] );
612             nPed++;
613           }
614           
615         } // loop over samples in bunch
616       } // loop over bunches
617
618       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
619
620       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
621       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
622       if (arrayPos >= fModules) {
623         //TODO: return an error message, if appopriate (perhaps if debug>0?)
624         return kFALSE;
625       }     
626       //Debug
627       if (arrayPos < 0 || arrayPos >= fModules) {
628         printf("Oh no: arrayPos = %i.\n", arrayPos); 
629       }
630       
631       fNChanFills++; // one more channel found, and profile to be filled
632       //NOTE: coordinates are (column, row) for the profiles
633       if ( in->IsLowGain() ) {
634         //fill the low gain histograms
635         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
636         if (nPed>0) { // only fill pedestal info in case it could be calculated
637           for ( i=0; i<nPed; i++) {
638             ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]); 
639           }
640         }
641       } 
642       else if ( in->IsHighGain() ) {    
643         //fill the high gain ones
644         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
645         if (nPed>0) { // only fill pedestal info in case it could be calculated
646           for ( i=0; i<nPed; i++) {
647             ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]); 
648           }       
649         }
650         // for warning checks
651         int idx = in->GetRow() + fRows * in->GetColumn();
652         ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min);
653       } 
654       else if ( in->IsLEDMonData() ) {
655         // for LED Mon data, the mapping class holds the gain info in the Row variable
656         // and the Strip number in the Column..
657         int gain = in->GetRow(); 
658         int stripId = in->GetColumn();
659         if (nPed>0 && stripId<fLEDRefs) {
660           if (gain == 0) {
661             for ( i=0; i<nPed; i++) {
662               ((TProfile*)fPedestalLEDRefLowGain[arrayPos])->Fill(stripId, pedSamples[i]);
663             }
664           }
665           else {
666             for ( i=0; i<nPed; i++) {
667               ((TProfile*)fPedestalLEDRefHighGain[arrayPos])->Fill(stripId, pedSamples[i]);
668             }
669           }
670         }
671       }
672
673       } // nsamples>0 check, some data found for this channel; not only trailer/header
674     }// end while over channel   
675   }//end while over DDL's, of input stream 
676
677  
678   return kTRUE;
679 }
680
681 //_____________________________________________________________________
682 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
683 {
684   //Saves all the histograms (or profiles, to be accurate) to the designated file
685   ValidateProfiles(); // make sure histos/profiles exist
686   TFile destFile(fileName, "recreate");
687   
688   if (destFile.IsZombie()) {
689     return kFALSE;
690   }
691   
692   destFile.cd();
693   
694   for (int i = 0; i < fModules; i++) {
695     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
696       fPeakMinusPedLowGain[i]->Write();
697     }
698     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) { 
699       fPeakMinusPedHighGain[i]->Write();
700     }
701     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
702       fPedestalLowGain[i]->Write();
703     }
704     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
705       fPedestalHighGain[i]->Write();
706     }
707     if( ((TProfile *)fPedestalLEDRefLowGain[i])->GetEntries() || saveEmptyHistos) {
708       fPedestalLEDRefLowGain[i]->Write();
709     }
710     if( ((TProfile *)fPedestalLEDRefHighGain[i])->GetEntries() || saveEmptyHistos) {
711       fPedestalLEDRefHighGain[i]->Write();
712     }
713     if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) { 
714       fPeakMinusPedHighGainHisto[i]->Write();
715     }
716
717   } 
718   
719   destFile.Close();
720   
721   return kTRUE;
722 }
723
724 //_____________________________________________________________________
725 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
726 {
727   
728   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
729   TH1::AddDirectory(kFALSE);
730   
731   TFile *sourceFile = new TFile(fileName);
732   if (sourceFile->IsZombie()) {
733     return kFALSE;//We couldn't load the reference
734   }
735
736   if (fReference) delete fReference;//Delete the reference object, if it already exists
737   fReference = 0;
738   
739   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
740  
741   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
742     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
743     fReference = 0;
744     return kFALSE;
745   }
746         
747   delete sourceFile;
748
749   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
750   //if we are called by someone who has set it to false...
751   TH1::AddDirectory(kTRUE);
752  
753   return kTRUE;//We succesfully loaded the object
754 }
755
756
757 //_____________________________________________________________________
758 Bool_t AliCaloCalibPedestal::SetReference(AliCaloCalibPedestal *ref)
759 {
760   if (fReference) delete fReference;//Delete the reference object, if it already exists
761   fReference = 0;
762   
763   fReference = ref;
764  
765   if (!fReference || (fReference->GetDetectorType() != fDetType)) {
766     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
767     fReference = 0;
768     return kFALSE;
769   }
770
771   return kTRUE;//We succesfully loaded the object
772 }
773
774 //_____________________________________________________________________
775 void AliCaloCalibPedestal::ValidateComparisonProfiles()
776 {
777   //Make sure the comparison histos exist
778   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
779   //the same time
780                                                 
781                                                 
782   //Then, loop for the requested number of modules
783   TString title, name;
784   for (int i = 0; i < fModules; i++) {
785     //Pedestals, low gain
786     name = "hPedlowgainDiff";
787     name += i;
788     title = "Pedestals difference, low gain, module ";
789     title += i; 
790     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
791                                             fColumns, 0.0, fColumns, 
792                                             fRows, fRowMin, fRowMax,"s"));
793   
794     //Pedestals, high gain
795     name = "hPedhighgainDiff";
796     name += i;
797     title = "Pedestals difference, high gain, module ";
798     title += i; 
799     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
800                                              fColumns, 0.0, fColumns, 
801                                              fRows, fRowMin, fRowMax,"s"));
802
803     //LED Ref/Mon pedestals, low gain
804     name = "hPedestalLEDReflowgainDiff";
805     name += i;
806     title = "Pedestal difference LEDRef, low gain, module ";
807     title += i; 
808     fPedestalLEDRefLowGainDiff.Add(new TProfile(name, title,
809                                                 fLEDRefs, 0.0, fLEDRefs, "s"));
810     
811     //LED Ref/Mon pedestals, high gain
812     name = "hPedestalLEDRefhighgainDiff";
813     name += i;
814     title = "Pedestal difference LEDRef, high gain, module ";
815     title += i; 
816     fPedestalLEDRefHighGainDiff.Add(new TProfile(name, title,
817                                                  fLEDRefs, 0.0, fLEDRefs, "s"));
818
819     //Peak-Pedestals, high gain
820     name = "hPeakMinusPedhighgainDiff";
821     name += i;
822     title = "Peak-Pedestal difference, high gain, module ";
823     title += i; 
824     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
825                                                  fColumns, 0.0, fColumns, 
826                                                  fRows, fRowMin, fRowMax,"s"));
827
828     //Peak-Pedestals, low gain
829     name = "hPeakMinusPedlowgainDiff";
830     name += i;
831     title = "Peak-Pedestal difference, low gain, module ";
832     title += i; 
833     fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
834                                                 fColumns, 0.0, fColumns, 
835                                                 fRows, fRowMin, fRowMax,"s"));
836   
837     //Pedestals, low gain
838     name = "hPedlowgainRatio";
839     name += i;
840     title = "Pedestals ratio, low gain, module ";
841     title += i; 
842     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
843                                              fColumns, 0.0, fColumns, 
844                                              fRows, fRowMin, fRowMax,"s"));
845   
846     //Pedestals, high gain
847     name = "hPedhighgainRatio";
848     name += i;
849     title = "Pedestals ratio, high gain, module ";
850     title += i; 
851     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
852                                               fColumns, 0.0, fColumns, 
853                                               fRows, fRowMin, fRowMax,"s"));
854
855     //LED Ref/Mon pedestals, low gain
856     name = "hPedestalLEDReflowgainRatio";
857     name += i;
858     title = "Pedestal ratio LEDRef, low gain, module ";
859     title += i; 
860     fPedestalLEDRefLowGainRatio.Add(new TProfile(name, title,
861                                                  fLEDRefs, 0.0, fLEDRefs, "s"));
862     
863     //LED Ref/Mon pedestals, high gain
864     name = "hPedestalLEDRefhighgainRatio";
865     name += i;
866     title = "Pedestal ratio LEDRef, high gain, module ";
867     title += i; 
868     fPedestalLEDRefHighGainRatio.Add(new TProfile(name, title,
869                                                   fLEDRefs, 0.0, fLEDRefs, "s"));
870   
871     //Peak-Pedestals, low gain
872     name = "hPeakMinusPedlowgainRatio";
873     name += i;
874     title = "Peak-Pedestal ratio, low gain, module ";
875     title += i; 
876     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
877                                                  fColumns, 0.0, fColumns, 
878                                                  fRows, fRowMin, fRowMax,"s"));
879   
880     //Peak-Pedestals, high gain
881     name = "hPeakMinusPedhighgainRatio";
882     name += i;
883     title = "Peak-Pedestal ratio, high gain, module ";
884     title += i; 
885     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
886                                                   fColumns, 0.0, fColumns, 
887                                                   fRows, fRowMin, fRowMax,"s"));
888     
889   }//end for nModules create the histograms
890 }
891
892 //_____________________________________________________________________
893 void AliCaloCalibPedestal::ComputeDiffAndRatio()
894 {
895   // calculate differences and ratios relative to a reference
896   ValidateProfiles(); // make sure histos/profiles exist
897   ValidateComparisonProfiles();//Make sure the comparison histos exist
898  
899   if (!fReference) {
900     return;//Return if the reference object isn't loaded
901   }
902
903   int bin = 0;
904   double diff = 0;
905   double ratio = 1;
906   for (int i = 0; i < fModules; i++) {
907     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
908     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
909     for (int j = 0; j < fColumns; j++) {
910       for (int k = 0; k < fRows; k++) {
911         bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
912
913         if (fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) > 0) {
914           diff = GetPeakProfileHighGain(i)->GetBinContent(bin) - fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
915           ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(bin, diff);
916           ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
917           ratio = GetPeakProfileHighGain(i)->GetBinContent(bin) / fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);  
918           ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinContent(bin, ratio);
919           ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinEntries(bin, 1);
920         }
921
922         if (fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) > 0) {
923           diff = GetPeakProfileLowGain(i)->GetBinContent(bin) - fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
924           ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(bin, diff);
925           ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
926           ratio = GetPeakProfileLowGain(i)->GetBinContent(bin) / fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);  
927           ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinContent(bin, ratio);
928           ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinEntries(bin, 1);
929         }
930
931         if (fReference->GetPedProfileHighGain(i)->GetBinContent(bin) > 0) {
932           diff = GetPedProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
933           ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(bin, diff);
934           ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
935           ratio = GetPedProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedProfileHighGain(i)->GetBinContent(bin);  
936           ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinContent(bin, ratio);
937           ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinEntries(bin, 1);
938         }
939
940         if (fReference->GetPedProfileLowGain(i)->GetBinContent(bin) > 0) {
941           diff = GetPedProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
942           ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(bin, diff);
943           ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
944           ratio = GetPedProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedProfileLowGain(i)->GetBinContent(bin);  
945           ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinContent(bin, ratio);
946           ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinEntries(bin, 1);
947         }
948
949       } // rows
950     } // columns
951
952     // same for LED Ref/Mon channels
953     for (int j = 0; j <= fLEDRefs; j++) {    
954       bin = j+1;//Note that we assume here that all histos have the same structure...
955
956       if (fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) > 0) {
957         diff = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
958         ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinContent(bin, diff);
959         ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinEntries(bin, 1);
960         ratio = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);  
961         ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinContent(bin, ratio);
962         ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinEntries(bin, 1);
963       }
964
965       if (fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) > 0) {
966         diff = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
967         ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinContent(bin, diff);
968         ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinEntries(bin, 1);
969         ratio = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);  
970         ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinContent(bin, ratio);
971         ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinEntries(bin, 1);
972       } 
973      
974     }
975
976   } // modules
977  
978 }
979
980 //_____________________________________________________________________
981 void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile)
982 { // look for hot/noisy towers
983   ValidateProfiles(); // make sure histos/profiles exist
984   ofstream * fout = 0;
985   char name[512];//Quite a long temp buffer, just in case the filename includes a path
986
987   if (hotMapFile) {
988     snprintf(name, 512, "%s.txt", hotMapFile);
989     fout = new ofstream(name);
990     if (!fout->is_open()) {
991       delete fout;
992       fout = 0;//Set the pointer to empty if the file was not opened
993     }
994   }
995  
996   for(int i = 0; i < fModules; i++){
997                 
998     //first we compute the peak-pedestal distribution for each supermodule...
999     if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) {
1000       double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin());
1001       double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin());
1002       TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1);
1003
1004       for (int j = 1; j <= fColumns; j++) {
1005         for (int k = 1; k <= fRows; k++) {        
1006           hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k));
1007         }
1008       }
1009
1010       //...and then we try to fit it
1011       double mean  = hPeakFit->GetMean();
1012       double sigma = hPeakFit->GetRMS();
1013       try {
1014         hPeakFit->Fit("gaus", "OQ", "",  mean - 3*sigma, mean + 3*sigma);
1015         mean  = hPeakFit->GetFunction("gaus")->GetParameter(1);
1016         sigma = hPeakFit->GetFunction("gaus")->GetParameter(2);
1017       }
1018       catch (const std::exception & e) {
1019         printf("AliCaloCalibPedestal: TH1D PeakFit exception %s", e.what()); 
1020       }      
1021       //hPeakFit->Draw();
1022
1023       delete hPeakFit;
1024
1025       //Then we look for warm/hot towers
1026       TH2F * hPeak2D = GetPeakHighGainHisto(i);
1027       hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) );
1028
1029       int idx = 0 ;
1030       int warnCounter = 0;
1031       for (int j = 1; j <= fColumns; j++) {
1032         for (int k = 1; k <= fRows; k++) {                              
1033           //we start looking for warm/warning towers...
1034           // histogram x-axis index
1035           idx = k-1 + fRows*(j-1); // this is what is used in the Fill call
1036           hPeak2D->GetXaxis()->SetRangeUser(idx, idx);
1037           warnCounter = (int) hPeak2D->Integral();
1038           if(warnCounter > fNEvents * fWarningFraction) {
1039             ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning);            
1040             /* printf("mod %d col %d row %d warnCounter %d - status %d\n", 
1041                i, j-1, k-1, warnCounter, (int) (kWarning)); */  
1042           }
1043           //...then we look for hot ones (towers whose values are greater than mean + X*sigma)  
1044           if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) {
1045             ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot); 
1046             /* printf("mod %d col %d row %d  binc %d - status %d\n", 
1047                i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */      
1048           }
1049
1050           //Write the status to the hot/warm map file, if the file is open.
1051           // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum)
1052           if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) {
1053             
1054             (*fout) << i << " " 
1055                     << (j - 1) << " " 
1056                     << (k - 1) << " " 
1057                     << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl;               }
1058           
1059         }
1060       }
1061
1062     } 
1063   }
1064   return;
1065 }
1066
1067 //_____________________________________________________________________
1068 void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile)
1069 {
1070   ValidateProfiles(); // make sure histos/profiles exist
1071   //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
1072   int countTot = 0;
1073   int countNew = 0;
1074   int countRes = 0;
1075   ofstream * fout = 0;
1076   ofstream * diff = 0;
1077   char name[512];//Quite a long temp buffer, just in case the filename includes a path
1078   
1079   if (deadMapFile) {
1080     snprintf(name, 512, "%s.txt", deadMapFile);
1081     fout = new ofstream(name);
1082     snprintf(name, 512, "%sdiff.txt", deadMapFile);
1083     diff = new ofstream(name);
1084     if (!fout->is_open()) {
1085       delete fout;
1086       fout = 0;//Set the pointer to empty if the file was not opened
1087     }
1088     if (!diff->is_open()) {
1089       delete diff;
1090       diff = 0;//Set the pointer to empty if the file was not opened
1091     }
1092   }
1093  
1094   for (int i = 0; i < fModules; i++) {
1095     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
1096       for (int j = 1; j <= fColumns; j++) {
1097         for (int k = 1; k <= fRows; k++) {
1098
1099           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead
1100             countTot++;//One more dead total
1101             if (fout) {
1102               (*fout) << i << " " 
1103                       << (j - 1) << " " 
1104                       << (k - 1) << " " 
1105                       << "1" << " " 
1106                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
1107             }
1108             
1109             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) {
1110               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased); 
1111               countNew++;//This tower wasn't dead before!
1112               if (diff) {
1113                 ( *diff) << i << " " 
1114                          << (j - 1) << " " 
1115                          << (k - 1) << " " 
1116                          << "1" << " " 
1117                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
1118               }
1119             } 
1120             else {
1121               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new         
1122             }
1123           } 
1124           else { //It's ALIVE!!
1125             //Don't bother with writing the live ones.
1126             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {
1127               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
1128               countRes++; //This tower was dead before => it's a miracle! :P
1129               if (diff) {
1130                 (*diff) << i << " " 
1131                         << (j - 1) << " " 
1132                         << (k - 1) << " " 
1133                         << "1" << " " 
1134                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
1135               }
1136             } 
1137             else {
1138               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
1139             }
1140           }
1141             
1142         }//end for k/rows
1143       }//end for j/columns
1144     }//end if GetEntries >= 0
1145   
1146   }//end for modules
1147  
1148  if (fout) {
1149    fout->close();
1150    delete fout;
1151  }
1152  
1153  fDeadTowers = countTot;
1154  fNewDeadTowers = countNew;
1155  fResurrectedTowers = countRes;
1156 }
1157
1158 //_____________________________________________________________________
1159 Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const
1160 {
1161   // Check if status info histo exists
1162   if (!fDeadMap[imod]) { 
1163     return kFALSE;
1164   }
1165
1166   //Check if channel is dead or hot.  
1167   Int_t status =  (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) );
1168   if(status == kAlive)
1169     return kFALSE;
1170   else 
1171     return kTRUE;
1172   
1173 }
1174
1175 //_____________________________________________________________________
1176 void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status)
1177 {
1178   ValidateProfiles(); // make sure histos/profiles exist
1179   //Set status of channel dead, hot, alive ...  
1180   ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status);   
1181 }