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