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