silvermy@ornl.gov - adding AliCaloCalibSignal, from Josh H (UT) and David S (ORNL)
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
CommitLineData
bd83f412 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
44ClassImp(AliCaloCalibSignal)
45
46using namespace std;
47
48// ctor; initialize everything in order to avoid compiler warnings
49// put some reasonable defaults
50AliCaloCalibSignal::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//_____________________________________________________________________
98AliCaloCalibSignal::~AliCaloCalibSignal()
99{
100 ClearObjects();
101}
102
103//_____________________________________________________________________
104void 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//_____________________________________________________________________
124AliCaloCalibSignal::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//_____________________________________________________________________
148AliCaloCalibSignal& 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//_____________________________________________________________________
158void 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//_____________________________________________________________________
203void 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//_____________________________________________________________________
212void 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//_____________________________________________________________________
223Bool_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//_____________________________________________________________________
248Bool_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//_____________________________________________________________________
276Bool_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//_____________________________________________________________________
381void 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//_____________________________________________________________________
411Bool_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}