Updated geometry including for the first time ACORDE
[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"
39
40#include <fstream>
41#include <iostream>
42
a235e2bc 43#include "AliCaloRawStream.h"
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(),
57 fPeakMinusPedLowGain(),
58 fPeakMinusPedHighGain(),
59 fPedestalLowGainDiff(),
60 fPedestalHighGainDiff(),
61 fPeakMinusPedLowGainDiff(),
62 fPeakMinusPedHighGainDiff(),
63 fPedestalLowGainRatio(),
64 fPedestalHighGainRatio(),
65 fPeakMinusPedLowGainRatio(),
66 fPeakMinusPedHighGainRatio(),
67 fDeadMap(),
68 fDeadTowers(0),
69 fNewDeadTowers(0),
70 fResurrectedTowers(0),
71 fReference(0),
72 fDetType(kNone),
73 fColumns(0),
74 fRows(0),
75 fModules(0),
76 fRunNumber(-1)
77{
78 //Default constructor. First we set the detector-type related constants.
79 if (detectorType == kPhos) {
80 fColumns = fgkPhosCols;
81 fRows = fgkPhosRows;
82 fModules = fgkPhosModules;
83 }
84 else {
85 //We'll just trust the enum to keep everything in line, so that if detectorType
86 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
87 //case, if someone intentionally gives another number
88 fColumns = fgkEmCalCols;
89 fRows = fgkEmCalRows;
90 fModules = fgkEmCalModules;
91 }
92 fDetType = detectorType;
93
94 //Then, loop for the requested number of modules
95 TString title, name;
96 for (int i = 0; i < fModules; i++) {
97 //Pedestals, low gain
98 name = "hPedlowgain";
99 name += i;
100 title = "Pedestals, low gain, module ";
101 title += i;
102 fPedestalLowGain.Add(new TProfile2D(name, title,
103 fColumns, 0.0, fColumns,
104 fRows, -fRows, 0.0));
105
106 //Pedestals, high gain
107 name = "hPedhighgain";
108 name += i;
109 title = "Pedestals, high gain, module ";
110 title += i;
111 fPedestalHighGain.Add(new TProfile2D(name, title,
112 fColumns, 0.0, fColumns,
113 fRows, -fRows, 0.0));
114
115 //Peak-Pedestals, low gain
116 name = "hPeakMinusPedlowgain";
117 name += i;
118 title = "Peak-Pedestal, low gain, module ";
119 title += i;
120 fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
121 fColumns, 0.0, fColumns,
122 fRows, -fRows, 0.0));
123
124 //Peak-Pedestals, high gain
125 name = "hPeakMinusPedhighgain";
126 name += i;
127 title = "Peak-Pedestal, high gain, module ";
128 title += i;
129 fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
130 fColumns, 0.0, fColumns,
131 fRows, -fRows, 0.0));
132
133 name = "hDeadMap";
134 name += i;
135 title = "Dead map, module ";
136 title += i;
137 fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
138 fRows, -fRows, 0.0));
139
140 }//end for nModules create the histograms
141
142 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
143 fPedestalLowGain.Compress();
144 fPedestalHighGain.Compress();
145 fPeakMinusPedLowGain.Compress();
146 fPeakMinusPedHighGain.Compress();
147 fDeadMap.Compress();
148 //Make them the owners of the profiles, so we don't need to care about deleting them
149 //fPedestalLowGain.SetOwner();
150 //fPedestalHighGain.SetOwner();
151 //fPeakMinusPedLowGain.SetOwner();
152 //fPeakMinusPedHighGain.SetOwner();
153
154}
155
156// dtor
157//_____________________________________________________________________
158AliCaloCalibPedestal::~AliCaloCalibPedestal()
159{
160 if (fReference) delete fReference;//Delete the reference object, if it has been loaded
161 //TObjArray will delete the histos/profiles when it is deleted.
162}
163
164// copy ctor
165//_____________________________________________________________________
166AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
167 TObject(ped),
168 fPedestalLowGain(),
169 fPedestalHighGain(),
170 fPeakMinusPedLowGain(),
171 fPeakMinusPedHighGain(),
172 fPedestalLowGainDiff(),
173 fPedestalHighGainDiff(),
174 fPeakMinusPedLowGainDiff(),
175 fPeakMinusPedHighGainDiff(),
176 fPedestalLowGainRatio(),
177 fPedestalHighGainRatio(),
178 fPeakMinusPedLowGainRatio(),
179 fPeakMinusPedHighGainRatio(),
180 fDeadMap(),
181 fDeadTowers(ped.GetDeadTowerCount()),
182 fNewDeadTowers(ped.GetDeadTowerNew()),
183 fResurrectedTowers(ped.GetDeadTowerResurrected()),
184 fReference( 0 ), //! note that we do not try to copy the reference info here
185 fDetType(ped.GetDetectorType()),
186 fColumns(ped.GetColumns()),
187 fRows(ped.GetRows()),
188 fModules(ped.GetModules()),
189 fRunNumber(ped.GetRunNumber())
190{
191 // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
192 //DS: this has not really been tested yet..
193 for (int i = 0; i < fModules; i++) {
194 fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
195 fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
196 fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
197 fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
198
199 fDeadMap.Add( ped.GetDeadMap(i) );
200 }//end for nModules
201
202 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
203 fPedestalLowGain.Compress();
204 fPedestalHighGain.Compress();
205 fPeakMinusPedLowGain.Compress();
206 fPeakMinusPedHighGain.Compress();
207 fDeadMap.Compress();
208}
209
210// assignment operator; use copy ctor to make life easy..
211//_____________________________________________________________________
212AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
213{
a235e2bc 214 // assignment operator; use copy ctor
356c3e0c 215 if (&source == this) return *this;
216
217 new (this) AliCaloCalibPedestal(source);
218 return *this;
219}
220
221//_____________________________________________________________________
222void AliCaloCalibPedestal::Reset()
223{
a235e2bc 224 // Reset all arrays/histograms
356c3e0c 225 for (int i = 0; i < fModules; i++) {
226 GetPedProfileLowGain(i)->Reset();
227 GetPedProfileHighGain(i)->Reset();
228 GetPeakProfileLowGain(i)->Reset();
229 GetPeakProfileHighGain(i)->Reset();
230 GetDeadMap(i)->Reset();
231
232 if (!fPedestalLowGainDiff.IsEmpty()) {
233 //This means that the comparison profiles have been created.
234
235 GetPedProfileLowGainDiff(i)->Reset();
236 GetPedProfileHighGainDiff(i)->Reset();
237 GetPeakProfileLowGainDiff(i)->Reset();
238 GetPeakProfileHighGainDiff(i)->Reset();
239
240 GetPedProfileLowGainRatio(i)->Reset();
241 GetPedProfileHighGainRatio(i)->Reset();
242 GetPeakProfileLowGainRatio(i)->Reset();
243 GetPeakProfileHighGainRatio(i)->Reset();
244 }
245 }
246 fDeadTowers = 0;
247 fNewDeadTowers = 0;
248 fResurrectedTowers = 0;
249
250 //To think about: should fReference be deleted too?... let's not do it this time, at least...
251}
252
253//_____________________________________________________________________
254Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
255{
a235e2bc 256 // Method to process=analyze one event in the data stream
356c3e0c 257 if (!in) return kFALSE; //Return right away if there's a null pointer
258
356c3e0c 259 int sample, i = 0; //The sample temp, and the sample number in current event.
260 int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
261 int gain = 0;
262
263 while (in->Next()) {
264 sample = in->GetSignal(); //Get the adc signal
265 if (sample < min) min = sample;
266 if (sample > max) max = sample;
267 i++;
268 if ( i >= in->GetTimeLength()) {
269 //If we're here then we're done with this tower
270 gain = 1 - in->IsLowGain();
271
272 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
273 if (arrayPos >= fModules) {
274 //TODO: return an error message, if appopriate (perhaps if debug>0?)
275 return kFALSE;
276 }
277
278 //Debug
279 if (arrayPos < 0 || arrayPos >= fModules) {
280 printf("Oh no: arrayPos = %i.\n", arrayPos);
281 }
282
283 //NOTE: coordinates are (column, row) for the profiles
284 if (gain == 0) {
285 //fill the low gain histograms
286 ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
287 ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
288 }
289 else {//fill the high gain ones
290 ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
291 ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
292 }//end if gain
293
294 max = fgkSampleMin; min = fgkSampleMax;
295 i = 0;
296
297 }//End if end of tower
298
299 }//end while, of stream
300
301 return kTRUE;
302}
303
304//_____________________________________________________________________
305Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
306{
307 //Saves all the histograms (or profiles, to be accurate) to the designated file
308
309 TFile destFile(fileName, "recreate");
310
311 if (destFile.IsZombie()) {
312 return kFALSE;
313 }
314
315 destFile.cd();
316
317 for (int i = 0; i < fModules; i++) {
318 if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
319 fPeakMinusPedLowGain[i]->Write();
320 }
321 if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
322 fPeakMinusPedHighGain[i]->Write();
323 }
324 if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
325 fPedestalLowGain[i]->Write();
326 }
327 if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
328 fPedestalHighGain[i]->Write();
329 }
330 }
331
332 destFile.Close();
333
334 return kTRUE;
335}
336
337//_____________________________________________________________________
338Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
339{
340
341 //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
342 TH1::AddDirectory(kFALSE);
343
344 TFile *sourceFile = new TFile(fileName);
345 if (sourceFile->IsZombie()) {
346 return kFALSE;//We couldn't load the reference
347 }
348
349 if (fReference) delete fReference;//Delete the reference object, if it already exists
350 fReference = 0;
351
352 fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
353
354 if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
355 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
356 fReference = 0;
357 return kFALSE;
358 }
359
360 delete sourceFile;
361
362 //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
363 //if we are called by someone who has set it to false...
364 TH1::AddDirectory(kTRUE);
365
366 return kTRUE;//We succesfully loaded the object
367}
368
369//_____________________________________________________________________
370void AliCaloCalibPedestal::ValidateComparisonProfiles()
371{
a235e2bc 372 //Make sure the comparison histos exist
356c3e0c 373 if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
374 //the same time
375
376
377 //Then, loop for the requested number of modules
378 TString title, name;
379 for (int i = 0; i < fModules; i++) {
380 //Pedestals, low gain
381 name = "hPedlowgainDiff";
382 name += i;
383 title = "Pedestals difference, low gain, module ";
384 title += i;
385 fPedestalLowGainDiff.Add(new TProfile2D(name, title,
386 fColumns, 0.0, fColumns,
387 fRows, -fRows, 0.0));
388
389 //Pedestals, high gain
390 name = "hPedhighgainDiff";
391 name += i;
392 title = "Pedestals difference, high gain, module ";
393 title += i;
394 fPedestalHighGainDiff.Add(new TProfile2D(name, title,
395 fColumns, 0.0, fColumns,
396 fRows, -fRows, 0.0));
397
398 //Peak-Pedestals, low gain
399 name = "hPeakMinusPedlowgainDiff";
400 name += i;
401 title = "Peak-Pedestal difference, low gain, module ";
402 title += i;
403 fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
404 fColumns, 0.0, fColumns,
405 fRows, -fRows, 0.0));
406
407 //Peak-Pedestals, high gain
408 name = "hPeakMinusPedhighgainDiff";
409 name += i;
410 title = "Peak-Pedestal difference, high gain, module ";
411 title += i;
412 fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
413 fColumns, 0.0, fColumns,
414 fRows, -fRows, 0.0));
415
416 //Pedestals, low gain
417 name = "hPedlowgainRatio";
418 name += i;
419 title = "Pedestals ratio, low gain, module ";
420 title += i;
421 fPedestalLowGainRatio.Add(new TProfile2D(name, title,
422 fColumns, 0.0, fColumns,
423 fRows, -fRows, 0.0));
424
425 //Pedestals, high gain
426 name = "hPedhighgainRatio";
427 name += i;
428 title = "Pedestals ratio, high gain, module ";
429 title += i;
430 fPedestalHighGainRatio.Add(new TProfile2D(name, title,
431 fColumns, 0.0, fColumns,
432 fRows, -fRows, 0.0));
433
434 //Peak-Pedestals, low gain
435 name = "hPeakMinusPedlowgainRatio";
436 name += i;
437 title = "Peak-Pedestal ratio, low gain, module ";
438 title += i;
439 fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
440 fColumns, 0.0, fColumns,
441 fRows, -fRows, 0.0));
442
443 //Peak-Pedestals, high gain
444 name = "hPeakMinusPedhighgainRatio";
445 name += i;
446 title = "Peak-Pedestal ratio, high gain, module ";
447 title += i;
448 fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
449 fColumns, 0.0, fColumns,
450 fRows, -fRows, 0.0));
451
452 }//end for nModules create the histograms
453}
454
455//_____________________________________________________________________
456void AliCaloCalibPedestal::ComputeDiffAndRatio()
457{
a235e2bc 458 // calculate differences and ratios relative to a reference
356c3e0c 459 ValidateComparisonProfiles();//Make sure the comparison histos exist
460
461 if (!fReference) {
462 return;//Return if the reference object isn't loaded
463 }
464
465 for (int i = 0; i < fModules; i++) {
466 //Compute the ratio of the histograms
467
468 ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
469 ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
470 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
471 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
472
473 //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
474 //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
475 for (int j = 0; j <= fColumns; j++) {
476 for (int k = 0; k <= fRows; k++) {
477 int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
478 double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
479 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
480 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
481
482 diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
483 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
484 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
485
486 diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
487 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
488 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
489
490 diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
491 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
492 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
493
494 } // rows
495 } // columns
496
497 } // modules
498
499}
500
501//_____________________________________________________________________
502void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
a235e2bc 503{
504 //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
356c3e0c 505 int countTot = 0;
506 int countNew = 0;
507 int countRes = 0;
508 ofstream * fout = 0;
509 ofstream * diff = 0;
510 char name[512];//Quite a long temp buffer, just in case the filename includes a path
511
512 if (deadMapFile) {
513 snprintf(name, 512, "%s.txt", deadMapFile);
514 fout = new ofstream(name);
515 snprintf(name, 512, "%sdiff.txt", deadMapFile);
516 diff = new ofstream(name);
517 if (!fout->is_open()) {
518 delete fout;
519 fout = 0;//Set the pointer to empty if the file was not opened
520 }
521 if (!diff->is_open()) {
522 delete diff;
523 fout = 0;//Set the pointer to empty if the file was not opened
524 }
525 }
526
527 for (int i = 0; i < fModules; i++) {
528 if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
529 for (int j = 1; j <= fColumns; j++) {
530 for (int k = 1; k <= fRows; k++) {
531
532 if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {//It's dead
533 countTot++;//One more dead total
534 if (fout) {
535 (*fout) << i << " "
536 << (fRows - k) << " "
537 << (j-1) << " "
538 << "1" << " "
539 << "0" << endl;//Write the status to the deadmap file, if the file is open.
540 }
541
542 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= threshold) {
543 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
544 countNew++;//This tower wasn't dead before!
545 if (diff) {
546 ( *diff) << i << " "
547 << (fRows - k) << " "
548 << (j - 1) << " "
549 << "1" << " "
550 << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
551 }
552 }
553 else {
554 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
555 }
556 }
557 else { //It's ALIVE!!
558 //Don't bother with writing the live ones.
559 //if (fout)
560 // (*fout) << i << " "
561 // << (fRows - k) << " "
562 // << (j - 1) << " "
563 // << "1" << " "
564 // << "1" << endl;//Write the status to the deadmap file, if the file is open.
565 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {
566 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
567 countRes++; //This tower was dead before => it's a miracle! :P
568 if (diff) {
569 (*diff) << i << " "
570 << (fRows - k) << " "
571 << (j - 1) << " "
572 << "1" << " "
573 << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
574 }
575 }
576 else {
577 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
578 }
579 }
580
581 }//end for k/rows
582 }//end for j/columns
583 }//end if GetEntries >= 0
584
585 }//end for modules
586
587 if (fout) {
588 fout->close();
589 delete fout;
590 }
591
592 fDeadTowers = countTot;
593 fNewDeadTowers = countNew;
594 fResurrectedTowers = countRes;
595}
596