coverity
[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,
afae9650 250 int resultArray[]) const
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,
afae9650 286 int resultArray[]) const
30aa89b0 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)
d50632ab 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
d50632ab 473 // We should also copy other pieces of info: counters and parameters
474 // (not number of columns and rows etc which should be the same)
475 // note that I just assign them here rather than Add them, but we
476 // normally just Add (e.g. in Preprocessor) one object so this should be fine.
477 fRunNumber = sig->GetRunNumber();
478 fStartTime = sig->GetStartTime();
479 fAmpCut = sig->GetAmpCut();
480 fReqFractionAboveAmpCutVal = sig->GetReqFractionAboveAmpCutVal();
481 fReqFractionAboveAmp = sig->GetReqFractionAboveAmp();
482 fAmpCutLEDRef = sig->GetAmpCutLEDRef();
483 fReqLEDRefAboveAmpCutVal = sig->GetReqLEDRefAboveAmpCutVal();
484 fHour = sig->GetHour();
485 fLatestHour = sig->GetLatestHour();
486 fUseAverage = sig->GetUseAverage();
487 fSecInAverage = sig->GetSecInAverage();
488 fNEvents = sig->GetNEvents();
489 fNAcceptedEvents = sig->GetNAcceptedEvents();
5e99faca 490
8bcca84c 491 return kTRUE;//We hopefully succesfully added info from the supplied object
bd83f412 492}
493
f4fc542c 494//_____________________________________________________________________
495Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
496{
497 // if fMapping is NULL the rawstream will crate its own mapping
30aa89b0 498 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
499 if (fDetType == kEmCal) {
500 rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
501 }
b07ee441 502 return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
f4fc542c 503}
bd83f412 504
505//_____________________________________________________________________
b07ee441 506Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
bd83f412 507{
508 // Method to process=analyze one event in the data stream
509 if (!in) return kFALSE; //Return right away if there's a null pointer
510
511 fNEvents++; // one more event
512
5e99faca 513 // use maximum numbers to set array sizes
ab962f7b 514 int iAmpValHighGain[fgkMaxTowers];
515 int iAmpValLowGain[fgkMaxTowers];
516 memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
517 memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
bd83f412 518
5e99faca 519 // also for LED reference
ab962f7b 520 int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
521 memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
5e99faca 522
18831b6c 523 int sample = 0; // temporary value
5e99faca 524 int gain = 0; // high or low gain
bd83f412 525
30aa89b0 526 // Number of Low and High gain, and LED Ref, channels for this event:
bd83f412 527 int nLowChan = 0;
528 int nHighChan = 0;
30aa89b0 529 int nLEDRefChan = 0;
bd83f412 530
ab962f7b 531 int iTowerNum = 0; // array index for regular towers
532 int iRefNum = 0; // array index for LED references
bd83f412 533
534 // loop first to get the fraction of channels with amplitudes above cut
32cd4c24 535
536 while (in->NextDDL()) {
537 while (in->NextChannel()) {
538
539 // counters
bd3cd9c2 540 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
56613d93 541 int nsamples = 0;
542
32cd4c24 543 while (in->NextBunch()) {
544 const UShort_t *sig = in->GetSignals();
56613d93 545 nsamples += in->GetBunchLength();
32cd4c24 546 for (Int_t i = 0; i < in->GetBunchLength(); i++) {
547 sample = sig[i];
548
549 // check if it's a min or max value
550 if (sample < min) min = sample;
551 if (sample > max) max = sample;
552
553 } // loop over samples in bunch
554 } // loop over bunches
555
56613d93 556 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
557
32cd4c24 558 gain = -1; // init to not valid value
bd83f412 559 //If we're here then we're done with this tower
5e99faca 560 if ( in->IsLowGain() ) {
561 gain = 0;
562 }
563 else if ( in->IsHighGain() ) {
564 gain = 1;
565 }
566 else if ( in->IsLEDMonData() ) {
567 gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
568 }
569
32cd4c24 570 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
bd83f412 571 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
bd83f412 572 //Debug
573 if (arrayPos < 0 || arrayPos >= fModules) {
8bcca84c 574 printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
5e99faca 575 return kFALSE;
bd83f412 576 }
577
5e99faca 578 if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
579 // get tower number for AmpVal array
ab962f7b 580 iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
5e99faca 581
582 if (gain == 0) {
583 // fill amplitude into the array
ab962f7b 584 iAmpValLowGain[iTowerNum] = max - min;
5e99faca 585 nLowChan++;
586 }
587 else if (gain==1) {//fill the high gain ones
588 // fill amplitude into the array
ab962f7b 589 iAmpValHighGain[iTowerNum] = max - min;
5e99faca 590 nHighChan++;
591 }//end if gain
592 } // regular tower
6e9e8153 593 else if ( in->IsLEDMonData() ) { // LED ref.;
594 // strip # is coded is 'column' in the channel maps
ab962f7b 595 iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
596 iLEDAmpVal[iRefNum] = max - min;
30aa89b0 597 nLEDRefChan++;
5e99faca 598 } // end of LED ref
56613d93 599
600 } // nsamples>0 check, some data found for this channel; not only trailer/header
32cd4c24 601 } // end while over channel
bd83f412 602
32cd4c24 603 }//end while over DDL's, of input stream
bd83f412 604
32cd4c24 605 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
606
30aa89b0 607 // now check if it was an LED event, using the LED Reference info per strip
5ea60b80 608
609 // by default all columns are accepted (init check to > 0)
610 int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
611 for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) {
612 checkResultArray[ia] = 1;
613 }
bd83f412 614 if (fReqFractionAboveAmp) {
452729ca 615 bool ok = false;
616 if (nHighChan > 0) {
ab962f7b 617 ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray);
452729ca 618 }
bd83f412 619 if (!ok) return false; // skip event
620 }
621
30aa89b0 622 // by default all columns are accepted (init check to > 0)
623 int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs];
624 for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) {
625 checkResultArrayLEDRef[ia] = 1;
626 }
627 if (fReqLEDRefAboveAmpCutVal) {
628 bool ok = false;
629 if (nLEDRefChan > 0) {
630 ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef);
631 }
632 if (!ok) return false; // skip event
633 }
634
bd83f412 635 fNAcceptedEvents++; // one more event accepted
636
637 if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
b07ee441 638 fStartTime = Timestamp;
bd83f412 639 }
640
b07ee441 641 fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
bd83f412 642 if (fLatestHour < fHour) {
643 fLatestHour = fHour;
644 }
645
8bcca84c 646 // it is a led event, now fill TTree
30aa89b0 647 // We also do the activity check for LEDRefs/Strips, but need to translate between column
648 // and strip indices for that; based on these relations:
649 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
650 // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
651 // which leads to
652 // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
653
8bcca84c 654 for(int i=0; i<fModules; i++){
30aa89b0 655 for(int j=0; j<fColumns; j++) {
656 int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
657 if (checkResultArray[i*fColumns + j]>0 && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check
8bcca84c 658 for(int k=0; k<fRows; k++){
bd83f412 659
ab962f7b 660 iTowerNum = GetTowerNum(i, j, k);
bd83f412 661
ab962f7b 662 if(iAmpValHighGain[iTowerNum]) {
663 fAmp = iAmpValHighGain[iTowerNum];
8bcca84c 664 fChannelNum = GetChannelNum(i,j,k,1);
ab962f7b 665 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
666 fNHighGain[iTowerNum]++;
bd83f412 667 }
ab962f7b 668 if(iAmpValLowGain[iTowerNum]) {
669 fAmp = iAmpValLowGain[iTowerNum];
8bcca84c 670 fChannelNum = GetChannelNum(i,j,k,0);
ab962f7b 671 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
672 fNLowGain[iTowerNum]++;
bd83f412 673 }
5e99faca 674 } // rows
30aa89b0 675 } // column passed check, and LED Ref for strip passed check (if any)
5e99faca 676 } // columns
677
678 // also LED refs
679 for(int j=0; j<fLEDRefs; j++){
5ea60b80 680 int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
30aa89b0 681 if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check
682 (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks
683 for (gain=0; gain<2; gain++) {
684 fRefNum = GetRefNum(i, j, gain);
685 if (iLEDAmpVal[fRefNum]) {
686 fAmp = iLEDAmpVal[fRefNum];
687 fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
688 fNRef[fRefNum]++;
689 }
690 } // gain
691 } // at least one column in strip passed check, and LED Ref passed check (if any)
bd83f412 692 }
5e99faca 693
694 } // modules
bd83f412 695
696 return kTRUE;
697}
698
699//_____________________________________________________________________
8bcca84c 700Bool_t AliCaloCalibSignal::Save(TString fileName)
bd83f412 701{
8bcca84c 702 //Saves all the TTrees to the designated file
bd83f412 703
704 TFile destFile(fileName, "recreate");
705
706 if (destFile.IsZombie()) {
707 return kFALSE;
708 }
709
710 destFile.cd();
711
8bcca84c 712 // save the trees
713 fTreeAmpVsTime->Write();
1bfe2e7a 714 fTreeLEDAmpVsTime->Write();
8bcca84c 715 if (fUseAverage) {
716 Analyze(); // get the latest and greatest averages
717 fTreeAvgAmpVsTime->Write();
1bfe2e7a 718 fTreeLEDAvgAmpVsTime->Write();
8bcca84c 719 }
720
721 destFile.Close();
722
723 return kTRUE;
724}
725
726//_____________________________________________________________________
727Bool_t AliCaloCalibSignal::Analyze()
728{
729 // Fill the tree holding the average values
730 if (!fUseAverage) { return kFALSE; }
731
732 // Reset the average TTree if Analyze has already been called earlier,
733 // meaning that the TTree could have been partially filled
734 if (fTreeAvgAmpVsTime->GetEntries() > 0) {
735 fTreeAvgAmpVsTime->Reset();
736 }
737
738 //0: setup variables for the TProfile plots that we'll use to do the averages
bd83f412 739 int numProfBins = 0;
740 double timeMin = 0;
741 double timeMax = 0;
8bcca84c 742 if (fSecInAverage > 0) {
743 numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
744 }
745 numProfBins += 2; // add extra buffer : first and last
746 double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
747 timeMin = - binSize;
748 timeMax = timeMin + numProfBins*binSize;
749
750 //1: set up TProfiles for the towers that had data
751 TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
752 memset(profile, 0, sizeof(profile));
a51e676d 753 const Int_t buffersize = 200;
754 char name[buffersize]; // for profile id and title
ab962f7b 755 int iTowerNum = 0;
8bcca84c 756
757 for (int i = 0; i<fModules; i++) {
758 for (int ic=0; ic<fColumns; ic++){
759 for (int ir=0; ir<fRows; ir++) {
760
ab962f7b 761 iTowerNum = GetTowerNum(i, ic, ir);
8bcca84c 762 // high gain
ab962f7b 763 if (fNHighGain[iTowerNum] > 0) {
8bcca84c 764 fChannelNum = GetChannelNum(i, ic, ir, 1);
a51e676d 765 snprintf(name,buffersize,"profileChan%d", fChannelNum);
8bcca84c 766 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
767 }
768
769 // same for low gain
ab962f7b 770 if (fNLowGain[iTowerNum] > 0) {
8bcca84c 771 fChannelNum = GetChannelNum(i, ic, ir, 0);
a51e676d 772 snprintf(name,buffersize,"profileChan%d", fChannelNum);
8bcca84c 773 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
774 }
775
776 } // rows
777 } // columns
778 } // modules
779
780 //2: fill profiles by looping over tree
781 // Set addresses for tree-readback also
782 fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
783 fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
784 fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
785
786 for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
787 fTreeAmpVsTime->GetEntry(ient);
788 if (profile[fChannelNum]) {
789 // profile should always have been created above, for active channels
790 profile[fChannelNum]->Fill(fHour, fAmp);
bd83f412 791 }
bd83f412 792 }
793
8bcca84c 794 // re-associating the branch addresses here seems to be needed for OK 'average' storage
795 fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
796 fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
797 fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
798 fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
799
800 //3: fill avg tree by looping over the profiles
801 for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
802 if (profile[fChannelNum]) { // profile was created
803 if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
804 for(int it=0; it<numProfBins; it++) {
805 if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
806 fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
807 fHour = profile[fChannelNum]->GetBinCenter(it+1);
808 fRMS = profile[fChannelNum]->GetBinError(it+1);
809 fTreeAvgAmpVsTime->Fill();
810 } // some entries for this bin
811 } // loop over bins
812 } // some entries for this profile
813 } // profile exists
814 } // loop over all possible channels
815
5e99faca 816
817 // and finally, go through same exercise for LED also..
818
819 //1: set up TProfiles for the towers that had data
820 TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
821 memset(profileLED, 0, sizeof(profileLED));
822
823 for (int i = 0; i<fModules; i++) {
824 for(int j=0; j<fLEDRefs; j++){
825 for (int gain=0; gain<2; gain++) {
826 fRefNum = GetRefNum(i, j, gain);
827 if (fNRef[fRefNum] > 0) {
b48d1356 828 snprintf(name, buffersize, "profileLEDRef%d", fRefNum);
5e99faca 829 profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
830 }
831 }// gain
832 }
833 } // modules
834
835 //2: fill profiles by looping over tree
836 // Set addresses for tree-readback also
837 fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
838 fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
839 fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
840
841 for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
842 fTreeLEDAmpVsTime->GetEntry(ient);
843 if (profileLED[fRefNum]) {
844 // profile should always have been created above, for active channels
845 profileLED[fRefNum]->Fill(fHour, fAmp);
846 }
847 }
848
849 // re-associating the branch addresses here seems to be needed for OK 'average' storage
850 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
851 fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
852 fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
853 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
854
855 //3: fill avg tree by looping over the profiles
856 for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
857 if (profileLED[fRefNum]) { // profile was created
858 if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
859 for(int it=0; it<numProfBins; it++) {
860 if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
861 fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
862 fHour = profileLED[fRefNum]->GetBinCenter(it+1);
863 fRMS = profileLED[fRefNum]->GetBinError(it+1);
864 fTreeLEDAvgAmpVsTime->Fill();
865 } // some entries for this bin
866 } // loop over bins
867 } // some entries for this profile
868 } // profile exists
869 } // loop over all possible channels
870
871 // OK, we're done..
872
bd83f412 873 return kTRUE;
874}
ab962f7b 875
876//_____________________________________________________________________
877Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId,
878 int *imod, int *icol, int *irow, int *igain) const
879{ // return the module, column, row, and gain for a given channel number
880 *igain = chanId/(fModules*fColumns*fRows);
881 *imod = (chanId/(fColumns*fRows)) % fModules;
882 *icol = (chanId/fRows) % fColumns;
883 *irow = chanId % fRows;
884 return kTRUE;
885}
886
887//_____________________________________________________________________
888Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId,
889 int *imod, int *istripMod, int *igain) const
890{ // return the module, stripModule, and gain for a given reference number
891 *igain = refId/(fModules*fLEDRefs);
892 *imod = (refId/(fLEDRefs)) % fModules;
893 *istripMod = refId % fLEDRefs;
894 return kTRUE;
895}