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