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