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