f132523bec030f88ae63b4ef72f567fa6779ad0e
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibPedestal.cxx
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
10 ClassImp(AliCaloCalibPedestal)
11
12 using namespace std;
13
14 // ctor; initialize everything in order to avoid compiler warnings
15 AliCaloCalibPedestal::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 //_____________________________________________________________________
120 AliCaloCalibPedestal::~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 //_____________________________________________________________________
128 AliCaloCalibPedestal::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 //_____________________________________________________________________
174 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
175 {
176   if (&source == this) return *this;
177
178   new (this) AliCaloCalibPedestal(source);
179   return *this;
180 }
181
182 //_____________________________________________________________________
183 void 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 //_____________________________________________________________________
214 Bool_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 //_____________________________________________________________________
266 Bool_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 //_____________________________________________________________________
299 Bool_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 //_____________________________________________________________________
331 void 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 //_____________________________________________________________________
416 void 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 //_____________________________________________________________________
462 void 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