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