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