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