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