Coverity corrections
[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);
8bcca84c 29// get some info:
30// fSignals->GetXXX..()
bd83f412 31// etc.
32//________________________________________________________________________
b07ee441 33#include <string>
34#include <sstream>
35#include <fstream>
bd83f412 36
8bcca84c 37#include "TProfile.h"
bd83f412 38#include "TFile.h"
39
f4fc542c 40#include "AliRawReader.h"
32cd4c24 41#include "AliCaloRawStreamV3.h"
bd83f412 42
43//The include file
44#include "AliCaloCalibSignal.h"
45
46ClassImp(AliCaloCalibSignal)
47
48using namespace std;
49
8bcca84c 50// variables for TTree filling; not sure if they should be static or not
18831b6c 51static int fChannelNum = 0; // for regular towers
52static int fRefNum = 0; // for LED
53static double fAmp = 0;
54static double fAvgAmp = 0;
55static double fRMS = 0;
8bcca84c 56
bd83f412 57// ctor; initialize everything in order to avoid compiler warnings
58// put some reasonable defaults
59AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
60 TObject(),
61 fDetType(kNone),
62 fColumns(0),
63 fRows(0),
5e99faca 64 fLEDRefs(0),
bd83f412 65 fModules(0),
f4fc542c 66 fCaloString(),
67 fMapping(NULL),
bd83f412 68 fRunNumber(-1),
69 fStartTime(0),
5ea60b80 70 fAmpCut(40), // min. 40 ADC counts as default
71 fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
bd83f412 72 fReqFractionAboveAmp(kTRUE),
30aa89b0 73 fAmpCutLEDRef(100), // min. 100 ADC counts as default
74 fReqLEDRefAboveAmpCutVal(kTRUE),
bd83f412 75 fHour(0),
76 fLatestHour(0),
77 fUseAverage(kTRUE),
78 fSecInAverage(1800),
79 fNEvents(0),
8bcca84c 80 fNAcceptedEvents(0),
81 fTreeAmpVsTime(NULL),
5e99faca 82 fTreeAvgAmpVsTime(NULL),
83 fTreeLEDAmpVsTime(NULL),
84 fTreeLEDAvgAmpVsTime(NULL)
bd83f412 85{
86 //Default constructor. First we set the detector-type related constants.
87 if (detectorType == kPhos) {
88 fColumns = fgkPhosCols;
89 fRows = fgkPhosRows;
5e99faca 90 fLEDRefs = fgkPhosLEDRefs;
bd83f412 91 fModules = fgkPhosModules;
f4fc542c 92 fCaloString = "PHOS";
bd83f412 93 }
94 else {
95 //We'll just trust the enum to keep everything in line, so that if detectorType
96 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
97 //case, if someone intentionally gives another number
bd3cd9c2 98 fColumns = AliEMCALGeoParams::fgkEMCALCols;
99 fRows = AliEMCALGeoParams::fgkEMCALRows;
100 fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
101 fModules = AliEMCALGeoParams::fgkEMCALModules;
f4fc542c 102 fCaloString = "EMCAL";
bd83f412 103 }
104
105 fDetType = detectorType;
106
8bcca84c 107 ResetInfo(); // trees and counters
bd83f412 108}
109
110// dtor
111//_____________________________________________________________________
112AliCaloCalibSignal::~AliCaloCalibSignal()
113{
8bcca84c 114 DeleteTrees();
bd83f412 115}
116
117//_____________________________________________________________________
8bcca84c 118void AliCaloCalibSignal::DeleteTrees()
bd83f412 119{
8bcca84c 120 // delete what was created in the ctor (TTrees)
121 if (fTreeAmpVsTime) delete fTreeAmpVsTime;
122 if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
5e99faca 123 if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
124 if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
8bcca84c 125 // and reset pointers
126 fTreeAmpVsTime = NULL;
127 fTreeAvgAmpVsTime = NULL;
5e99faca 128 fTreeLEDAmpVsTime = NULL;
129 fTreeLEDAvgAmpVsTime = NULL;
bd83f412 130
131 return;
132}
133
134// copy ctor
135//_____________________________________________________________________
136AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
137 TObject(sig),
138 fDetType(sig.GetDetectorType()),
139 fColumns(sig.GetColumns()),
140 fRows(sig.GetRows()),
5e99faca 141 fLEDRefs(sig.GetLEDRefs()),
bd83f412 142 fModules(sig.GetModules()),
f4fc542c 143 fCaloString(sig.GetCaloString()),
a51e676d 144 fMapping(), //! note that we are not copying the map info
bd83f412 145 fRunNumber(sig.GetRunNumber()),
146 fStartTime(sig.GetStartTime()),
147 fAmpCut(sig.GetAmpCut()),
148 fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
149 fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
30aa89b0 150 fAmpCutLEDRef(sig.GetAmpCutLEDRef()),
151 fReqLEDRefAboveAmpCutVal(sig.GetReqLEDRefAboveAmpCutVal()),
bd83f412 152 fHour(sig.GetHour()),
153 fLatestHour(sig.GetLatestHour()),
154 fUseAverage(sig.GetUseAverage()),
155 fSecInAverage(sig.GetSecInAverage()),
156 fNEvents(sig.GetNEvents()),
8bcca84c 157 fNAcceptedEvents(sig.GetNAcceptedEvents()),
a51e676d 158 fTreeAmpVsTime(),
159 fTreeAvgAmpVsTime(),
160 fTreeLEDAmpVsTime(),
161 fTreeLEDAvgAmpVsTime()
bd83f412 162{
8bcca84c 163 // also the TTree contents
bd83f412 164 AddInfo(&sig);
a51e676d 165 for (Int_t i = 0; i<2*fgkMaxTowers; i++) {
166 if(i < fgkMaxTowers){
167 fNHighGain[i] = sig.fNHighGain[i];
168 fNLowGain[i] = sig.fNLowGain[i];
169 }
170 fNRef[i] = sig.fNRef[i];
171 }
172
173
174
bd83f412 175}
176
177// assignment operator; use copy ctor to make life easy..
178//_____________________________________________________________________
179AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
180{
181 // assignment operator; use copy ctor
182 if (&source == this) return *this;
183
184 new (this) AliCaloCalibSignal(source);
185 return *this;
186}
187
188//_____________________________________________________________________
8bcca84c 189void AliCaloCalibSignal::CreateTrees()
bd83f412 190{
8bcca84c 191 // initialize trees
192 // first, regular version
193 fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
194
195 fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
196 fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
197 fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
bd83f412 198
8bcca84c 199 // then, average version
200 fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
201
202 fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
203 fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
204 fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
205 fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
206
5e99faca 207 // then same for LED..
208 fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
209 fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
210 fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
211 fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
212
213 fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
214 fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
215 fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
216 fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
217 fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
218
8bcca84c 219 return;
bd83f412 220}
221
222//_____________________________________________________________________
8bcca84c 223void AliCaloCalibSignal::ResetInfo()
ab962f7b 224{ // reset trees and counters
bd83f412 225 Zero(); // set all counters to 0
8bcca84c 226 DeleteTrees(); // delete previous stuff
227 CreateTrees(); // and create some new ones
bd83f412 228 return;
229}
230
231//_____________________________________________________________________
232void AliCaloCalibSignal::Zero()
233{
8bcca84c 234 // set all counters to 0; not cuts etc. though
bd83f412 235 fHour = 0;
236 fLatestHour = 0;
237 fNEvents = 0;
238 fNAcceptedEvents = 0;
8bcca84c 239
5ea60b80 240 // Set the number of points for each tower: Amp vs. Time
8bcca84c 241 memset(fNHighGain, 0, sizeof(fNHighGain));
242 memset(fNLowGain, 0, sizeof(fNLowGain));
5e99faca 243 // and LED reference
244 memset(fNRef, 0, sizeof(fNRef));
8bcca84c 245
bd83f412 246 return;
247}
248
249//_____________________________________________________________________
ab962f7b 250Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal,
5ea60b80 251 int resultArray[])
ab962f7b 252{ // check fraction of towers, per column, that are above amplitude cut
5ea60b80 253 Bool_t returnCode = false;
bd83f412 254
ab962f7b 255 int iTowerNum = 0;
5ea60b80 256 double fraction = 0;
bd83f412 257 for (int i = 0; i<fModules; i++) {
258 for (int j = 0; j<fColumns; j++) {
5ea60b80 259 int nAbove = 0;
bd83f412 260 for (int k = 0; k<fRows; k++) {
ab962f7b 261 iTowerNum = GetTowerNum(i,j,k);
262 if (iAmpVal[iTowerNum] > fAmpCut) {
bd83f412 263 nAbove++;
264 }
265 }
5ea60b80 266 resultArray[i*fColumns +j] = 0; // init. to denied
267 if (nAbove > 0) {
268 fraction = (1.0*nAbove) / fRows;
269 /*
270 printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
271 i, j, nAbove, fraction);
272 */
273 if (fraction > fReqFractionAboveAmpCutVal) {
274 resultArray[i*fColumns + j] = nAbove;
275 returnCode = true;
276 }
277 }
bd83f412 278 }
5ea60b80 279 } // modules loop
bd83f412 280
5ea60b80 281 return returnCode;
bd83f412 282}
283
30aa89b0 284
285//_____________________________________________________________________
286Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal,
287 int resultArray[])
288{ // check which LEDRef/Mon strips are above amplitude cut
289 Bool_t returnCode = false;
290
291 int iRefNum = 0;
292 int gain = 1; // look at high gain; this should be rather saturated usually..
293 for (int i = 0; i<fModules; i++) {
294 for (int j = 0; j<fLEDRefs; j++) {
295 iRefNum = GetRefNum(i, j, gain);
296 if (iAmpVal[iRefNum] > fAmpCutLEDRef) {
297 resultArray[i*fLEDRefs +j] = 1; // enough signal
298 returnCode = true;
299 }
300 else {
301 resultArray[i*fLEDRefs +j] = 0; // not enough signal
302 }
303
304 /*
305 printf("DS mod %d LEDRef %d ampVal %d\n",
306 i, j, iAmpVal[iRefNum]);
307 */
308 } // LEDRefs
309 } // modules loop
310
311 return returnCode;
312}
313
b07ee441 314// Parameter/cut handling
315//_____________________________________________________________________
316void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
ab962f7b 317{ // set parameters from file
b07ee441 318 static const string delimitor("::");
319
320 // open, check input file
321 ifstream in( parameterFile );
322 if( !in ) {
323 printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
324 return;
325 }
326
327 // Note: this method is a bit more complicated than it really has to be
328 // - allowing for multiple entries per line, arbitrary order of the
329 // different variables etc. But I wanted to try and do this in as
330 // correct a C++ way as I could (as an exercise).
331
332 // read in
333 char readline[1024];
334 while ((in.rdstate() & ios::failbit) == 0 ) {
335
336 // Read into the raw char array and then construct a string
337 // to do the searching
338 in.getline(readline, 1024);
339 istringstream s(readline);
340
341 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
342
ab962f7b 343 string keyValue;
344 s >> keyValue;
b07ee441 345
346 // check stream status
1d059662 347 if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break;
b07ee441 348
349 // skip rest of line if comments found
ab962f7b 350 if( keyValue.substr( 0, 2 ) == "//" ) break;
b07ee441 351
ab962f7b 352 // look for "::" in keyValue pair
353 size_t position = keyValue.find( delimitor );
b07ee441 354 if( position == string::npos ) {
ab962f7b 355 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
b07ee441 356 }
357
ab962f7b 358 // split keyValue pair
359 string key( keyValue.substr( 0, position ) );
360 string value( keyValue.substr( position+delimitor.size(),
361 keyValue.size()-delimitor.size() ) );
b07ee441 362
363 // check value does not contain a new delimitor
364 if( value.find( delimitor ) != string::npos ) {
ab962f7b 365 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
b07ee441 366 }
367
368 // debug: check key value pair
369 // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
370
371 // if the key matches with something we expect, we assign the new value
30aa89b0 372 if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") ||
373 (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ) {
b07ee441 374 istringstream iss(value);
375 printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
376
377 if (key == "fAmpCut") {
378 iss >> fAmpCut;
379 }
380 else if (key == "fReqFractionAboveAmpCutVal") {
381 iss >> fReqFractionAboveAmpCutVal;
382 }
30aa89b0 383 else if (key == "fAmpCutLEDRef") {
384 iss >> fAmpCutLEDRef;
385 }
b07ee441 386 else if (key == "fSecInAverage") {
387 iss >> fSecInAverage;
388 }
389 } // some match found/expected
390
391 }
392 }
393
394 in.close();
30aa89b0 395 return;
b07ee441 396}
397
398//_____________________________________________________________________
399void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
ab962f7b 400{ // write parameters to file
b07ee441 401 static const string delimitor("::");
402 ofstream out( parameterFile );
403 out << "// " << parameterFile << endl;
404 out << "fAmpCut" << "::" << fAmpCut << endl;
405 out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
30aa89b0 406 out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl;
b07ee441 407 out << "fSecInAverage" << "::" << fSecInAverage << endl;
408
409 out.close();
410 return;
411}
412
bd83f412 413//_____________________________________________________________________
414Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
415{
5e99faca 416 // note/FIXME: we are not yet adding correctly the info for fN{HighGain,LowGain,Ref} here - but consider this a feature for now (20080905): we'll do Analyze() unless entries were found for a tower in this original object.
417
8bcca84c 418 // add info from sig's TTrees to ours..
419 TTree *sigAmp = sig->GetTreeAmpVsTime();
420 TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
421
422 // we could try some merging via TList or what also as a more elegant approach
423 // but I wanted with the stupid/simple and hopefully safe approach of looping
424 // over what we want to add..
425
426 // associate variables for sigAmp and sigAvgAmp:
427 sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
428 sigAmp->SetBranchAddress("fHour",&fHour);
429 sigAmp->SetBranchAddress("fAmp",&fAmp);
430
431 // loop over the trees.. note that since we use the same variables we should not need
432 // to do any assignments between the getting and filling
433 for (int i=0; i<sigAmp->GetEntries(); i++) {
434 sigAmp->GetEntry(i);
435 fTreeAmpVsTime->Fill();
436 }
4fb34e09 437
8bcca84c 438 sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
439 sigAvgAmp->SetBranchAddress("fHour",&fHour);
440 sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
441 sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
bd83f412 442
8bcca84c 443 for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
444 sigAvgAmp->GetEntry(i);
445 fTreeAvgAmpVsTime->Fill();
446 }
bd83f412 447
5e99faca 448 // also LED..
449 TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
450 TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
451
452 // associate variables for sigAmp and sigAvgAmp:
453 sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
454 sigLEDAmp->SetBranchAddress("fHour",&fHour);
455 sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
456
457 // loop over the trees.. note that since we use the same variables we should not need
458 // to do any assignments between the getting and filling
459 for (int i=0; i<sigLEDAmp->GetEntries(); i++) {
460 sigLEDAmp->GetEntry(i);
461 fTreeLEDAmpVsTime->Fill();
462 }
463
464 sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
465 sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
466 sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
467 sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
468
469 for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
470 sigLEDAvgAmp->GetEntry(i);
471 fTreeLEDAvgAmpVsTime->Fill();
472 }
473
474
8bcca84c 475 return kTRUE;//We hopefully succesfully added info from the supplied object
bd83f412 476}
477
f4fc542c 478//_____________________________________________________________________
479Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
480{
481 // if fMapping is NULL the rawstream will crate its own mapping
30aa89b0 482 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
483 if (fDetType == kEmCal) {
484 rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
485 }
b07ee441 486 return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
f4fc542c 487}
bd83f412 488
489//_____________________________________________________________________
b07ee441 490Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
bd83f412 491{
492 // Method to process=analyze one event in the data stream
493 if (!in) return kFALSE; //Return right away if there's a null pointer
494
495 fNEvents++; // one more event
496
5e99faca 497 // use maximum numbers to set array sizes
ab962f7b 498 int iAmpValHighGain[fgkMaxTowers];
499 int iAmpValLowGain[fgkMaxTowers];
500 memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
501 memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
bd83f412 502
5e99faca 503 // also for LED reference
ab962f7b 504 int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
505 memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
5e99faca 506
18831b6c 507 int sample = 0; // temporary value
5e99faca 508 int gain = 0; // high or low gain
bd83f412 509
30aa89b0 510 // Number of Low and High gain, and LED Ref, channels for this event:
bd83f412 511 int nLowChan = 0;
512 int nHighChan = 0;
30aa89b0 513 int nLEDRefChan = 0;
bd83f412 514
ab962f7b 515 int iTowerNum = 0; // array index for regular towers
516 int iRefNum = 0; // array index for LED references
bd83f412 517
518 // loop first to get the fraction of channels with amplitudes above cut
32cd4c24 519
520 while (in->NextDDL()) {
521 while (in->NextChannel()) {
522
523 // counters
bd3cd9c2 524 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
56613d93 525 int nsamples = 0;
526
32cd4c24 527 while (in->NextBunch()) {
528 const UShort_t *sig = in->GetSignals();
56613d93 529 nsamples += in->GetBunchLength();
32cd4c24 530 for (Int_t i = 0; i < in->GetBunchLength(); i++) {
531 sample = sig[i];
532
533 // check if it's a min or max value
534 if (sample < min) min = sample;
535 if (sample > max) max = sample;
536
537 } // loop over samples in bunch
538 } // loop over bunches
539
56613d93 540 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
541
32cd4c24 542 gain = -1; // init to not valid value
bd83f412 543 //If we're here then we're done with this tower
5e99faca 544 if ( in->IsLowGain() ) {
545 gain = 0;
546 }
547 else if ( in->IsHighGain() ) {
548 gain = 1;
549 }
550 else if ( in->IsLEDMonData() ) {
551 gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
552 }
553
32cd4c24 554 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
bd83f412 555 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
bd83f412 556 //Debug
557 if (arrayPos < 0 || arrayPos >= fModules) {
8bcca84c 558 printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
5e99faca 559 return kFALSE;
bd83f412 560 }
561
5e99faca 562 if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
563 // get tower number for AmpVal array
ab962f7b 564 iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
5e99faca 565
566 if (gain == 0) {
567 // fill amplitude into the array
ab962f7b 568 iAmpValLowGain[iTowerNum] = max - min;
5e99faca 569 nLowChan++;
570 }
571 else if (gain==1) {//fill the high gain ones
572 // fill amplitude into the array
ab962f7b 573 iAmpValHighGain[iTowerNum] = max - min;
5e99faca 574 nHighChan++;
575 }//end if gain
576 } // regular tower
6e9e8153 577 else if ( in->IsLEDMonData() ) { // LED ref.;
578 // strip # is coded is 'column' in the channel maps
ab962f7b 579 iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
580 iLEDAmpVal[iRefNum] = max - min;
30aa89b0 581 nLEDRefChan++;
5e99faca 582 } // end of LED ref
56613d93 583
584 } // nsamples>0 check, some data found for this channel; not only trailer/header
32cd4c24 585 } // end while over channel
bd83f412 586
32cd4c24 587 }//end while over DDL's, of input stream
bd83f412 588
32cd4c24 589 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
590
30aa89b0 591 // now check if it was an LED event, using the LED Reference info per strip
5ea60b80 592
593 // by default all columns are accepted (init check to > 0)
594 int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
595 for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) {
596 checkResultArray[ia] = 1;
597 }
bd83f412 598 if (fReqFractionAboveAmp) {
452729ca 599 bool ok = false;
600 if (nHighChan > 0) {
ab962f7b 601 ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray);
452729ca 602 }
bd83f412 603 if (!ok) return false; // skip event
604 }
605
30aa89b0 606 // by default all columns are accepted (init check to > 0)
607 int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs];
608 for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) {
609 checkResultArrayLEDRef[ia] = 1;
610 }
611 if (fReqLEDRefAboveAmpCutVal) {
612 bool ok = false;
613 if (nLEDRefChan > 0) {
614 ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef);
615 }
616 if (!ok) return false; // skip event
617 }
618
bd83f412 619 fNAcceptedEvents++; // one more event accepted
620
621 if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
b07ee441 622 fStartTime = Timestamp;
bd83f412 623 }
624
b07ee441 625 fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
bd83f412 626 if (fLatestHour < fHour) {
627 fLatestHour = fHour;
628 }
629
8bcca84c 630 // it is a led event, now fill TTree
30aa89b0 631 // We also do the activity check for LEDRefs/Strips, but need to translate between column
632 // and strip indices for that; based on these relations:
633 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
634 // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
635 // which leads to
636 // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
637
8bcca84c 638 for(int i=0; i<fModules; i++){
30aa89b0 639 for(int j=0; j<fColumns; j++) {
640 int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
641 if (checkResultArray[i*fColumns + j]>0 && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check
8bcca84c 642 for(int k=0; k<fRows; k++){
bd83f412 643
ab962f7b 644 iTowerNum = GetTowerNum(i, j, k);
bd83f412 645
ab962f7b 646 if(iAmpValHighGain[iTowerNum]) {
647 fAmp = iAmpValHighGain[iTowerNum];
8bcca84c 648 fChannelNum = GetChannelNum(i,j,k,1);
ab962f7b 649 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
650 fNHighGain[iTowerNum]++;
bd83f412 651 }
ab962f7b 652 if(iAmpValLowGain[iTowerNum]) {
653 fAmp = iAmpValLowGain[iTowerNum];
8bcca84c 654 fChannelNum = GetChannelNum(i,j,k,0);
ab962f7b 655 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
656 fNLowGain[iTowerNum]++;
bd83f412 657 }
5e99faca 658 } // rows
30aa89b0 659 } // column passed check, and LED Ref for strip passed check (if any)
5e99faca 660 } // columns
661
662 // also LED refs
663 for(int j=0; j<fLEDRefs; j++){
5ea60b80 664 int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
30aa89b0 665 if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check
666 (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks
667 for (gain=0; gain<2; gain++) {
668 fRefNum = GetRefNum(i, j, gain);
669 if (iLEDAmpVal[fRefNum]) {
670 fAmp = iLEDAmpVal[fRefNum];
671 fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
672 fNRef[fRefNum]++;
673 }
674 } // gain
675 } // at least one column in strip passed check, and LED Ref passed check (if any)
bd83f412 676 }
5e99faca 677
678 } // modules
bd83f412 679
680 return kTRUE;
681}
682
683//_____________________________________________________________________
8bcca84c 684Bool_t AliCaloCalibSignal::Save(TString fileName)
bd83f412 685{
8bcca84c 686 //Saves all the TTrees to the designated file
bd83f412 687
688 TFile destFile(fileName, "recreate");
689
690 if (destFile.IsZombie()) {
691 return kFALSE;
692 }
693
694 destFile.cd();
695
8bcca84c 696 // save the trees
697 fTreeAmpVsTime->Write();
1bfe2e7a 698 fTreeLEDAmpVsTime->Write();
8bcca84c 699 if (fUseAverage) {
700 Analyze(); // get the latest and greatest averages
701 fTreeAvgAmpVsTime->Write();
1bfe2e7a 702 fTreeLEDAvgAmpVsTime->Write();
8bcca84c 703 }
704
705 destFile.Close();
706
707 return kTRUE;
708}
709
710//_____________________________________________________________________
711Bool_t AliCaloCalibSignal::Analyze()
712{
713 // Fill the tree holding the average values
714 if (!fUseAverage) { return kFALSE; }
715
716 // Reset the average TTree if Analyze has already been called earlier,
717 // meaning that the TTree could have been partially filled
718 if (fTreeAvgAmpVsTime->GetEntries() > 0) {
719 fTreeAvgAmpVsTime->Reset();
720 }
721
722 //0: setup variables for the TProfile plots that we'll use to do the averages
bd83f412 723 int numProfBins = 0;
724 double timeMin = 0;
725 double timeMax = 0;
8bcca84c 726 if (fSecInAverage > 0) {
727 numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
728 }
729 numProfBins += 2; // add extra buffer : first and last
730 double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
731 timeMin = - binSize;
732 timeMax = timeMin + numProfBins*binSize;
733
734 //1: set up TProfiles for the towers that had data
735 TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
736 memset(profile, 0, sizeof(profile));
a51e676d 737 const Int_t buffersize = 200;
738 char name[buffersize]; // for profile id and title
ab962f7b 739 int iTowerNum = 0;
8bcca84c 740
741 for (int i = 0; i<fModules; i++) {
742 for (int ic=0; ic<fColumns; ic++){
743 for (int ir=0; ir<fRows; ir++) {
744
ab962f7b 745 iTowerNum = GetTowerNum(i, ic, ir);
8bcca84c 746 // high gain
ab962f7b 747 if (fNHighGain[iTowerNum] > 0) {
8bcca84c 748 fChannelNum = GetChannelNum(i, ic, ir, 1);
a51e676d 749 snprintf(name,buffersize,"profileChan%d", fChannelNum);
8bcca84c 750 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
751 }
752
753 // same for low gain
ab962f7b 754 if (fNLowGain[iTowerNum] > 0) {
8bcca84c 755 fChannelNum = GetChannelNum(i, ic, ir, 0);
a51e676d 756 snprintf(name,buffersize,"profileChan%d", fChannelNum);
8bcca84c 757 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
758 }
759
760 } // rows
761 } // columns
762 } // modules
763
764 //2: fill profiles by looping over tree
765 // Set addresses for tree-readback also
766 fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
767 fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
768 fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
769
770 for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
771 fTreeAmpVsTime->GetEntry(ient);
772 if (profile[fChannelNum]) {
773 // profile should always have been created above, for active channels
774 profile[fChannelNum]->Fill(fHour, fAmp);
bd83f412 775 }
bd83f412 776 }
777
8bcca84c 778 // re-associating the branch addresses here seems to be needed for OK 'average' storage
779 fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
780 fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
781 fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
782 fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
783
784 //3: fill avg tree by looping over the profiles
785 for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
786 if (profile[fChannelNum]) { // profile was created
787 if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
788 for(int it=0; it<numProfBins; it++) {
789 if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
790 fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
791 fHour = profile[fChannelNum]->GetBinCenter(it+1);
792 fRMS = profile[fChannelNum]->GetBinError(it+1);
793 fTreeAvgAmpVsTime->Fill();
794 } // some entries for this bin
795 } // loop over bins
796 } // some entries for this profile
797 } // profile exists
798 } // loop over all possible channels
799
5e99faca 800
801 // and finally, go through same exercise for LED also..
802
803 //1: set up TProfiles for the towers that had data
804 TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
805 memset(profileLED, 0, sizeof(profileLED));
806
807 for (int i = 0; i<fModules; i++) {
808 for(int j=0; j<fLEDRefs; j++){
809 for (int gain=0; gain<2; gain++) {
810 fRefNum = GetRefNum(i, j, gain);
811 if (fNRef[fRefNum] > 0) {
812 sprintf(name, "profileLEDRef%d", fRefNum);
813 profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
814 }
815 }// gain
816 }
817 } // modules
818
819 //2: fill profiles by looping over tree
820 // Set addresses for tree-readback also
821 fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
822 fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
823 fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
824
825 for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
826 fTreeLEDAmpVsTime->GetEntry(ient);
827 if (profileLED[fRefNum]) {
828 // profile should always have been created above, for active channels
829 profileLED[fRefNum]->Fill(fHour, fAmp);
830 }
831 }
832
833 // re-associating the branch addresses here seems to be needed for OK 'average' storage
834 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
835 fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
836 fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
837 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
838
839 //3: fill avg tree by looping over the profiles
840 for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
841 if (profileLED[fRefNum]) { // profile was created
842 if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
843 for(int it=0; it<numProfBins; it++) {
844 if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
845 fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
846 fHour = profileLED[fRefNum]->GetBinCenter(it+1);
847 fRMS = profileLED[fRefNum]->GetBinError(it+1);
848 fTreeLEDAvgAmpVsTime->Fill();
849 } // some entries for this bin
850 } // loop over bins
851 } // some entries for this profile
852 } // profile exists
853 } // loop over all possible channels
854
855 // OK, we're done..
856
bd83f412 857 return kTRUE;
858}
ab962f7b 859
860//_____________________________________________________________________
861Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId,
862 int *imod, int *icol, int *irow, int *igain) const
863{ // return the module, column, row, and gain for a given channel number
864 *igain = chanId/(fModules*fColumns*fRows);
865 *imod = (chanId/(fColumns*fRows)) % fModules;
866 *icol = (chanId/fRows) % fColumns;
867 *irow = chanId % fRows;
868 return kTRUE;
869}
870
871//_____________________________________________________________________
872Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId,
873 int *imod, int *istripMod, int *igain) const
874{ // return the module, stripModule, and gain for a given reference number
875 *igain = refId/(fModules*fLEDRefs);
876 *imod = (refId/(fLEDRefs)) % fModules;
877 *istripMod = refId % fLEDRefs;
878 return kTRUE;
879}