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