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