silvermy@ornl.gov - adding AliCaloCalibSignal, from Josh H (UT) and David S (ORNL)
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
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: AliCaloCalibSignal.cxx $ */
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   fSignals = new AliCaloCalibSignal( fDetType );
24   // AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
25   fNumModules = fSignals->GetModules();
26 */
27 // fed an event:
28 //  fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
29 // asked to draw graphs or profiles:
30 //  fSignals->GetGraphAmpVsTimeHighGain(module,column,row)->Draw("ap");
31 // or
32 //  fSignals->GetProfAmpVsTimeHighGain(module,column,row)->Draw();
33 // etc.
34 //________________________________________________________________________
35
36 #include "TFile.h"
37
38 #include "AliRawEventHeaderBase.h"
39 #include "AliCaloRawStream.h"
40
41 //The include file
42 #include "AliCaloCalibSignal.h"
43
44 ClassImp(AliCaloCalibSignal)
45
46 using namespace std;
47
48 // ctor; initialize everything in order to avoid compiler warnings
49 // put some reasonable defaults
50 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :  
51   TObject(),
52   fDetType(kNone),
53   fColumns(0),
54   fRows(0),
55   fModules(0),
56   fRunNumber(-1),
57   fStartTime(0),
58   fAmpCut(50),
59   fReqFractionAboveAmpCutVal(0.8),
60   fReqFractionAboveAmp(kTRUE),
61   fHour(0),
62   fLatestHour(0),
63   fUseAverage(kTRUE),
64   fSecInAverage(1800), 
65   fNEvents(0),
66   fNAcceptedEvents(0)
67 {
68   //Default constructor. First we set the detector-type related constants.
69   if (detectorType == kPhos) {
70     fColumns = fgkPhosCols;
71     fRows = fgkPhosRows;
72     fModules = fgkPhosModules;
73   } 
74   else {
75     //We'll just trust the enum to keep everything in line, so that if detectorType
76     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
77     //case, if someone intentionally gives another number
78     fColumns = fgkEmCalCols;
79     fRows = fgkEmCalRows;
80     fModules = fgkEmCalModules;
81   }
82
83   fDetType = detectorType;
84
85   // Set the number of points for each Amp vs. Time graph to 0
86   memset(fNHighGain, 0, sizeof(fNHighGain));
87   memset(fNLowGain, 0, sizeof(fNLowGain));
88
89   CreateGraphs(); // set up the TGraphs
90
91   // init TProfiles to NULL=0 also
92   memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
93   memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
94
95
96 // dtor
97 //_____________________________________________________________________
98 AliCaloCalibSignal::~AliCaloCalibSignal()
99 {
100   ClearObjects();
101 }
102
103 //_____________________________________________________________________
104 void AliCaloCalibSignal::ClearObjects()
105 {
106   // delete what was created in the ctor (TGraphs), and possible later (TProfiles)
107   for (int i=0; i<fgkMaxTowers; i++) {
108     if ( fGraphAmpVsTimeHighGain[i] ) { delete fGraphAmpVsTimeHighGain[i]; }
109     if ( fGraphAmpVsTimeLowGain[i] ) { delete fGraphAmpVsTimeLowGain[i]; }
110     if ( fProfAmpVsTimeHighGain[i] ) { delete fProfAmpVsTimeHighGain[i]; }
111     if ( fProfAmpVsTimeLowGain[i] ) { delete fProfAmpVsTimeLowGain[i]; }
112   }
113   // set pointers
114   memset(fGraphAmpVsTimeHighGain, 0, sizeof(fGraphAmpVsTimeHighGain));
115   memset(fGraphAmpVsTimeLowGain, 0, sizeof(fGraphAmpVsTimeLowGain));
116   memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
117   memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
118
119   return;
120 }
121
122 // copy ctor
123 //_____________________________________________________________________
124 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
125   TObject(sig),
126   fDetType(sig.GetDetectorType()),
127   fColumns(sig.GetColumns()),
128   fRows(sig.GetRows()),
129   fModules(sig.GetModules()),
130   fRunNumber(sig.GetRunNumber()),
131   fStartTime(sig.GetStartTime()),
132   fAmpCut(sig.GetAmpCut()),
133   fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
134   fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
135   fHour(sig.GetHour()),
136   fLatestHour(sig.GetLatestHour()),
137   fUseAverage(sig.GetUseAverage()),
138   fSecInAverage(sig.GetSecInAverage()),
139   fNEvents(sig.GetNEvents()),
140   fNAcceptedEvents(sig.GetNAcceptedEvents())
141 {
142   // also the TGraph contents
143   AddInfo(&sig);
144 }
145
146 // assignment operator; use copy ctor to make life easy..
147 //_____________________________________________________________________
148 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
149 {
150   // assignment operator; use copy ctor
151   if (&source == this) return *this;
152
153   new (this) AliCaloCalibSignal(source);
154   return *this;
155 }
156
157 //_____________________________________________________________________
158 void AliCaloCalibSignal::CreateGraphs()
159 {
160   //Then, loop for the requested number of modules
161   TString title, name;
162   for (int i = 0; i < fModules; i++) {
163     
164     // Amplitude vs. Time graph for each channel
165     for(int ic=0;ic < fColumns;ic++){
166       for(int ir=0;ir < fRows;ir++){
167         
168         int id = GetTowerNum(i, ic, ir);
169         
170         // high gain graph
171         name = "fGraphAmpVsTimeHighGain_"; name += i;
172         name += "_"; name += ic;
173         name += "_"; name += ir;
174         title = "Amp vs. Time High Gain Mod "; title += i;
175         title += " Col "; title += ic;
176         title += " Row "; title += ir;
177         
178         fGraphAmpVsTimeHighGain[id] = new TGraph();
179         fGraphAmpVsTimeHighGain[id]->SetName(name);
180         fGraphAmpVsTimeHighGain[id]->SetTitle(title);
181         fGraphAmpVsTimeHighGain[id]->SetMarkerStyle(20);
182         
183         // Low Gain
184         name = "fGraphAmpVsTimeLowGain_"; name += i;
185         name += "_"; name += ic;
186         name += "_"; name += ir;
187         title = "Amp vs. Time Low Gain Mod "; title += i;
188         title += " Col "; title += ic;
189         title += " Row "; title += ir;
190         
191         fGraphAmpVsTimeLowGain[id] = new TGraph();
192         fGraphAmpVsTimeLowGain[id]->SetName(name);
193         fGraphAmpVsTimeLowGain[id]->SetTitle(title);
194         fGraphAmpVsTimeLowGain[id]->SetMarkerStyle(20);
195         
196       }
197     }
198
199   }//end for nModules 
200 }
201
202 //_____________________________________________________________________
203 void AliCaloCalibSignal::Reset()
204 {
205   Zero(); // set all counters to 0
206   ClearObjects(); // delete previous TGraphs and TProfiles
207   CreateGraphs(); // and create some new ones
208   return;
209 }
210
211 //_____________________________________________________________________
212 void AliCaloCalibSignal::Zero()
213 {
214   // set all counters to 0; not cuts etc.though
215   fHour = 0;
216   fLatestHour = 0;
217   fNEvents = 0;
218   fNAcceptedEvents = 0;
219   return;
220 }
221
222 //_____________________________________________________________________
223 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
224 {
225   int nAbove = 0;
226     
227   int TowerNum = 0;
228   for (int i = 0; i<fModules; i++) {
229     for (int j = 0; j<fColumns; j++) {
230       for (int k = 0; k<fRows; k++) {
231         TowerNum = GetTowerNum(i,j,k);
232         if (AmpVal[TowerNum] > fAmpCut) { 
233           nAbove++;
234         }
235       }
236     }
237   }
238   
239   double fraction = (1.0*nAbove) / nTotChan;
240   
241   if (fraction > fReqFractionAboveAmpCutVal) {  
242     return true;
243   }
244   else return false;
245 }
246
247 //_____________________________________________________________________
248 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
249 {
250   // just do this for the basic graphs/profiles that get filled in ProcessEvent
251   // may not have data for all channels, but let's just Add everything..
252   for (int i = 0; i < fModules; i++) {
253     for (int j = 0; j < fColumns; j++) {
254       for (int k = 0; k < fRows; k++) {
255         
256         int id = GetTowerNum(i,j,k);
257         if(fUseAverage){
258           GetProfAmpVsTimeHighGain(id)->Add(sig->GetProfAmpVsTimeHighGain(id));
259           GetProfAmpVsTimeLowGain(id)->Add(sig->GetProfAmpVsTimeLowGain(id));
260         }
261         else{     
262           //DS
263 //        sig->GetGraphAmpVsTimeHighGain(i,j,k);
264 //        sig->GetGraphAmpVsTimeLowGain(i,j,k);
265         }
266
267       }//end for nModules 
268     }//end for nColumns
269   }//end for nRows
270
271   return kTRUE;//We succesfully added info from the supplied object
272 }
273
274
275 //_____________________________________________________________________
276 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderBase *aliHeader)
277
278   // Method to process=analyze one event in the data stream
279   if (!in) return kFALSE; //Return right away if there's a null pointer
280   
281   fNEvents++; // one more event
282
283   // PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes
284   int AmpValHighGain[fgkMaxTowers];
285   int AmpValLowGain[fgkMaxTowers];
286
287   memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
288   memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
289
290   int sample, i = 0; //The sample temp, and the sample number in current event.
291   int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
292   int gain = 0;
293   
294   // Number of Low and High gain channels for this event:
295   int nLowChan = 0; 
296   int nHighChan = 0; 
297
298   int TowerNum = 0; // array index for TGraphs etc.
299
300   // loop first to get the fraction of channels with amplitudes above cut
301   while (in->Next()) {
302     sample = in->GetSignal(); //Get the adc signal
303     if (sample < min) min = sample;
304     if (sample > max) max = sample;
305     i++;
306     if ( i >= in->GetTimeLength()) {
307       //If we're here then we're done with this tower
308       gain = 1 - in->IsLowGain();
309       
310       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
311       if (arrayPos >= fModules) {
312         //TODO: return an error message, if appopriate (perhaps if debug>0?)
313         return kFALSE;
314       } 
315       
316       //Debug
317       if (arrayPos < 0 || arrayPos >= fModules) {
318         printf("Oh no: arrayPos = %i.\n", arrayPos); 
319       }
320
321       // get tower number for AmpVal array
322       TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
323
324       if (gain == 0) {
325         // fill amplitude into the array           
326         AmpValLowGain[TowerNum] = max - min;
327         nLowChan++;
328       } 
329       else if (gain==1) {//fill the high gain ones
330         // fill amplitude into the array
331         AmpValHighGain[TowerNum] = max - min;
332         nHighChan++;
333       }//end if gain
334
335       
336       max = fgkSampleMin; min = fgkSampleMax;
337       i = 0;
338       
339     }//End if end of tower
340    
341   }//end while, of stream
342   
343   // now check if it was a led event, only use high gain (that should be sufficient)
344   if (fReqFractionAboveAmp) {
345     bool ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
346     if (!ok) return false; // skip event
347   }
348
349   fNAcceptedEvents++; // one more event accepted
350
351   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
352     fStartTime = aliHeader->Get("Timestamp");
353   }
354
355   fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
356   if (fLatestHour < fHour) {
357     fLatestHour = fHour; 
358   }
359   
360   // it is a led event, now fill graphs (maybe profiles later)
361   for(int i=0;i<fModules;i++){
362     for(int j=0;j<fColumns;j++){
363       for(int k=0;k<fRows;k++){
364         
365         TowerNum = GetTowerNum(i, j, k); 
366
367         if(AmpValHighGain[TowerNum]) {
368           fGraphAmpVsTimeHighGain[TowerNum]->SetPoint(fNHighGain[TowerNum]++,fHour,AmpValHighGain[TowerNum]);
369         }
370         if(AmpValLowGain[TowerNum]) {
371           fGraphAmpVsTimeLowGain[TowerNum]->SetPoint(fNLowGain[TowerNum]++,fHour,AmpValLowGain[TowerNum]);
372         }
373       }
374     }
375   }
376   
377   return kTRUE;
378 }
379
380 //_____________________________________________________________________
381 void AliCaloCalibSignal::CreateProfile(int imod, int ic, int ir, int towerId, int gain,
382                                        int nbins, double min, double max)
383 { //! create/setup a TProfile
384   TString title, name;   
385   if (gain == 0) { 
386     name = "fProfAmpVsTimeLowGain_";   
387     title = "Amp vs. Time Low Gain Mod "; 
388   } 
389   else if (gain == 1) { 
390     name = "fProfAmpVsTimeHighGain_"; 
391     title = "Amp vs. Time High Gain Mod "; 
392   } 
393   name += imod;
394   name += "_"; name += ic;
395   name += "_"; name += ir;
396   title += imod;
397   title += " Col "; title += ic;
398   title += " Row "; title += ir;
399             
400   // use "s" option for RMS
401   if (gain==0) { 
402     fProfAmpVsTimeLowGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
403   }
404   else if (gain==1) {
405     fProfAmpVsTimeHighGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
406   }
407
408   return;
409 }
410 //_____________________________________________________________________
411 Bool_t AliCaloCalibSignal::Save(TString fileName, Bool_t saveEmptyGraphs)
412 {
413   //Saves all the histograms (or profiles, to be accurate) to the designated file
414   
415   TFile destFile(fileName, "recreate");
416   
417   if (destFile.IsZombie()) {
418     return kFALSE;
419   }
420   
421   destFile.cd();
422
423   // setup variables for the TProfile plot
424   int numProfBins = 0;
425   double timeMin = 0;
426   double timeMax = 0;
427   if (fUseAverage) {
428     if (fSecInAverage > 0) {
429       numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
430     }
431     numProfBins += 2; // add extra buffer : first and last
432     double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
433     timeMin = - binSize;
434     timeMax = timeMin + numProfBins*binSize;
435   }
436
437   int numGraphPoints= 0;
438   int TowerNum = 0;    
439   for (int i = 0; i < fModules; i++) {
440     
441     for(int ic=0;ic < fColumns;ic++){
442       for(int ir=0;ir < fRows;ir++){
443
444         TowerNum = GetTowerNum(i, ic, ir); 
445
446         // 1st: high gain
447         numGraphPoints= fGraphAmpVsTimeHighGain[TowerNum]->GetN();
448         if( numGraphPoints>0 || saveEmptyGraphs) {
449           
450           // average the graphs points over time if requested and put them in a profile plot
451           if(fUseAverage && numGraphPoints>0) {
452             
453             // get the values
454             double *graphX = fGraphAmpVsTimeHighGain[TowerNum]->GetX();
455             double *graphY = fGraphAmpVsTimeHighGain[TowerNum]->GetY();
456
457             // create the TProfile: 1 is for High gain              
458             CreateProfile(i, ic, ir, TowerNum, 1,
459                           numProfBins, timeMin, timeMax);
460
461             // loop over graph points and fill profile
462             for(int ip=0; ip < numGraphPoints; ip++){
463               fProfAmpVsTimeHighGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
464             }
465             
466             fProfAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
467             fProfAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
468             fProfAmpVsTimeHighGain[TowerNum]->Write();
469
470           }
471            else{
472              //otherwise, just save the graphs and forget the profiling
473              fGraphAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
474              fGraphAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
475              fGraphAmpVsTimeHighGain[TowerNum]->Write();
476            }
477           
478         } // low gain graph info should be saved in one form or another
479         
480         // 2nd: now go to the low gain case
481         numGraphPoints= fGraphAmpVsTimeLowGain[TowerNum]->GetN();
482         if( numGraphPoints>0 || saveEmptyGraphs) {
483           
484           // average the graphs points over time if requested and put them in a profile plot
485           if(fUseAverage && numGraphPoints>0) {
486             
487             double *graphX = fGraphAmpVsTimeLowGain[TowerNum]->GetX();
488             double *graphY = fGraphAmpVsTimeLowGain[TowerNum]->GetY();
489             
490             // create the TProfile: 0 is for Low gain       
491             CreateProfile(i, ic, ir, TowerNum, 0,
492                           numProfBins, timeMin, timeMax);
493
494             // loop over graph points and fill profile
495             for(int ip=0; ip < numGraphPoints; ip++){
496               fProfAmpVsTimeLowGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
497             }
498             
499             fProfAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
500             fProfAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
501             fProfAmpVsTimeLowGain[TowerNum]->Write();
502
503           }
504           else{
505              //otherwise, just save the graphs and forget the profiling
506             fGraphAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
507             fGraphAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
508             fGraphAmpVsTimeLowGain[TowerNum]->Write();
509           }
510           
511         } // low gain graph info should be saved in one form or another
512
513       } // fRows
514     } // fColumns
515
516   } // fModules
517   destFile.Close();
518   
519   return kTRUE;
520 }