]>
Commit | Line | Data |
---|---|---|
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 | ||
52 | ClassImp(AliCaloCalibSignal) | |
53 | ||
54 | using namespace std; | |
55 | ||
56 | // variables for TTree filling; not sure if they should be static or not | |
57 | static int fChannelNum = 0; // for regular towers | |
58 | static int fRefNum = 0; // for LED | |
59 | static double fAmp = 0; | |
60 | static double fAvgAmp = 0; | |
61 | static double fRMS = 0; | |
62 | ||
63 | // ctor; initialize everything in order to avoid compiler warnings | |
64 | // put some reasonable defaults | |
65 | AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) : | |
66 | TObject(), | |
67 | fDetType(kNone), | |
68 | fColumns(0), | |
69 | fRows(0), | |
70 | fLEDRefs(0), | |
71 | fModules(0), | |
72 | fCaloString(), | |
73 | fMapping(NULL), | |
74 | fFittingAlgorithm(0), | |
75 | fRawAnalyzer(0), | |
76 | fRunNumber(-1), | |
77 | fStartTime(0), | |
78 | fAmpCut(40), // min. 40 ADC counts as default | |
79 | fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default | |
80 | fReqFractionAboveAmp(kTRUE), | |
81 | fAmpCutLEDRef(100), // min. 100 ADC counts as default | |
82 | fReqLEDRefAboveAmpCutVal(kTRUE), | |
83 | fHour(0), | |
84 | fLatestHour(0), | |
85 | fUseAverage(kTRUE), | |
86 | fSecInAverage(1800), | |
87 | fDownscale(10), | |
88 | fNEvents(0), | |
89 | fNAcceptedEvents(0), | |
90 | fTreeAmpVsTime(NULL), | |
91 | fTreeAvgAmpVsTime(NULL), | |
92 | fTreeLEDAmpVsTime(NULL), | |
93 | fTreeLEDAvgAmpVsTime(NULL) | |
94 | { | |
95 | //Default constructor. First we set the detector-type related constants. | |
96 | if (detectorType == kPhos) { | |
97 | fColumns = fgkPhosCols; | |
98 | fRows = fgkPhosRows; | |
99 | fLEDRefs = fgkPhosLEDRefs; | |
100 | fModules = fgkPhosModules; | |
101 | fCaloString = "PHOS"; | |
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 | |
107 | fColumns = AliEMCALGeoParams::fgkEMCALCols; | |
108 | fRows = AliEMCALGeoParams::fgkEMCALRows; | |
109 | fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs; | |
110 | fModules = AliEMCALGeoParams::fgkEMCALModules; | |
111 | fCaloString = "EMCAL"; | |
112 | } | |
113 | ||
114 | fDetType = detectorType; | |
115 | SetFittingAlgorithm(Algo::kStandard); | |
116 | ResetInfo(); // trees and counters | |
117 | } | |
118 | ||
119 | // dtor | |
120 | //_____________________________________________________________________ | |
121 | AliCaloCalibSignal::~AliCaloCalibSignal() | |
122 | { | |
123 | DeleteTrees(); | |
124 | } | |
125 | ||
126 | //_____________________________________________________________________ | |
127 | void AliCaloCalibSignal::DeleteTrees() | |
128 | { | |
129 | // delete what was created in the ctor (TTrees) | |
130 | if (fTreeAmpVsTime) delete fTreeAmpVsTime; | |
131 | if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime; | |
132 | if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime; | |
133 | if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime; | |
134 | // and reset pointers | |
135 | fTreeAmpVsTime = NULL; | |
136 | fTreeAvgAmpVsTime = NULL; | |
137 | fTreeLEDAmpVsTime = NULL; | |
138 | fTreeLEDAvgAmpVsTime = NULL; | |
139 | ||
140 | return; | |
141 | } | |
142 | ||
143 | // copy ctor | |
144 | //_____________________________________________________________________ | |
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()), | |
165 | // fDownscale(sig.GetDownscale()), | |
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 | // | |
186 | // assignment operator; use copy ctor to make life easy.. | |
187 | //_____________________________________________________________________ | |
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 | //} | |
196 | ||
197 | //_____________________________________________________________________ | |
198 | void AliCaloCalibSignal::CreateTrees() | |
199 | { | |
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"); | |
207 | ||
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 | ||
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 | ||
228 | return; | |
229 | } | |
230 | ||
231 | //_____________________________________________________________________ | |
232 | void AliCaloCalibSignal::ResetInfo() | |
233 | { // reset trees and counters | |
234 | Zero(); // set all counters to 0 | |
235 | DeleteTrees(); // delete previous stuff | |
236 | CreateTrees(); // and create some new ones | |
237 | return; | |
238 | } | |
239 | ||
240 | //_____________________________________________________________________ | |
241 | void AliCaloCalibSignal::Zero() | |
242 | { | |
243 | // set all counters to 0; not cuts etc. though | |
244 | fHour = 0; | |
245 | fLatestHour = 0; | |
246 | fNEvents = 0; | |
247 | fNAcceptedEvents = 0; | |
248 | ||
249 | // Set the number of points for each tower: Amp vs. Time | |
250 | memset(fNHighGain, 0, sizeof(fNHighGain)); | |
251 | memset(fNLowGain, 0, sizeof(fNLowGain)); | |
252 | // and LED reference | |
253 | memset(fNRef, 0, sizeof(fNRef)); | |
254 | ||
255 | return; | |
256 | } | |
257 | ||
258 | //_____________________________________________________________________ | |
259 | Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal, | |
260 | int resultArray[]) const | |
261 | { // check fraction of towers, per column, that are above amplitude cut | |
262 | Bool_t returnCode = false; | |
263 | ||
264 | int iTowerNum = 0; | |
265 | double fraction = 0; | |
266 | for (int i = 0; i<fModules; i++) { | |
267 | for (int j = 0; j<fColumns; j++) { | |
268 | int nAbove = 0; | |
269 | for (int k = 0; k<fRows; k++) { | |
270 | iTowerNum = GetTowerNum(i,j,k); | |
271 | if (iAmpVal[iTowerNum] > fAmpCut) { | |
272 | nAbove++; | |
273 | } | |
274 | } | |
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 | } | |
287 | } | |
288 | } // modules loop | |
289 | ||
290 | return returnCode; | |
291 | } | |
292 | ||
293 | ||
294 | //_____________________________________________________________________ | |
295 | Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal, | |
296 | int resultArray[]) const | |
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 | ||
323 | // Parameter/cut handling | |
324 | //_____________________________________________________________________ | |
325 | void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile) | |
326 | { // set parameters from file | |
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 | ||
352 | string keyValue; | |
353 | s >> keyValue; | |
354 | ||
355 | // check stream status | |
356 | if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break; | |
357 | ||
358 | // skip rest of line if comments found | |
359 | if( keyValue.substr( 0, 2 ) == "//" ) break; | |
360 | ||
361 | // look for "::" in keyValue pair | |
362 | size_t position = keyValue.find( delimitor ); | |
363 | if( position == string::npos ) { | |
364 | printf("wrong format for key::value pair: %s\n", keyValue.c_str()); | |
365 | } | |
366 | ||
367 | // split keyValue pair | |
368 | string key( keyValue.substr( 0, position ) ); | |
369 | string value( keyValue.substr( position+delimitor.size(), | |
370 | keyValue.size()-delimitor.size() ) ); | |
371 | ||
372 | // check value does not contain a new delimitor | |
373 | if( value.find( delimitor ) != string::npos ) { | |
374 | printf("wrong format for key::value pair: %s\n", keyValue.c_str()); | |
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 | |
381 | if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") || | |
382 | (key == "fAmpCutLEDRef") || (key == "fSecInAverage") || | |
383 | (key == "fFittingAlgorithm") || (key == "fDownscale") ) { | |
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 | } | |
393 | else if (key == "fAmpCutLEDRef") { | |
394 | iss >> fAmpCutLEDRef; | |
395 | } | |
396 | else if (key == "fSecInAverage") { | |
397 | iss >> fSecInAverage; | |
398 | } | |
399 | else if (key == "fFittingAlgorithm") { | |
400 | iss >> fFittingAlgorithm; | |
401 | SetFittingAlgorithm( fFittingAlgorithm ); | |
402 | } | |
403 | else if (key == "fDownscale") { | |
404 | iss >> fDownscale; | |
405 | } | |
406 | } // some match found/expected | |
407 | ||
408 | } | |
409 | } | |
410 | ||
411 | in.close(); | |
412 | return; | |
413 | } | |
414 | ||
415 | //_____________________________________________________________________ | |
416 | void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile) | |
417 | { // write parameters to file | |
418 | static const string delimitor("::"); | |
419 | ofstream out( parameterFile ); | |
420 | out << "// " << parameterFile << endl; | |
421 | out << "fAmpCut" << "::" << fAmpCut << endl; | |
422 | out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl; | |
423 | out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl; | |
424 | out << "fSecInAverage" << "::" << fSecInAverage << endl; | |
425 | out << "fFittingAlgorithm" << "::" << fFittingAlgorithm << endl; | |
426 | out << "fDownscale" << "::" << fDownscale << endl; | |
427 | ||
428 | out.close(); | |
429 | return; | |
430 | } | |
431 | ||
432 | //_____________________________________________________________________ | |
433 | void 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 | //_____________________________________________________________________ | |
442 | Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig) | |
443 | { | |
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 | ||
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 | } | |
465 | ||
466 | sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum); | |
467 | sigAvgAmp->SetBranchAddress("fHour",&fHour); | |
468 | sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp); | |
469 | sigAvgAmp->SetBranchAddress("fRMS",&fRMS); | |
470 | ||
471 | for (int i=0; i<sigAvgAmp->GetEntries(); i++) { | |
472 | sigAvgAmp->GetEntry(i); | |
473 | fTreeAvgAmpVsTime->Fill(); | |
474 | } | |
475 | ||
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 | ||
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(); | |
517 | fDownscale = sig->GetDownscale(); | |
518 | fNEvents = sig->GetNEvents(); | |
519 | fNAcceptedEvents = sig->GetNAcceptedEvents(); | |
520 | ||
521 | return kTRUE;//We hopefully succesfully added info from the supplied object | |
522 | } | |
523 | ||
524 | //_____________________________________________________________________ | |
525 | Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader) | |
526 | { | |
527 | // if fMapping is NULL the rawstream will crate its own mapping | |
528 | AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping); | |
529 | if (fDetType == kEmCal) { | |
530 | rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range | |
531 | } | |
532 | ||
533 | return ProcessEvent( &rawStream, rawReader->GetTimestamp() ); | |
534 | } | |
535 | ||
536 | //_____________________________________________________________________ | |
537 | Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp) | |
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 | ||
544 | if ( (fNEvents%fDownscale)!=0 ) return kFALSE; // mechanism to skip some of the input events, if we want | |
545 | ||
546 | // use maximum numbers to set array sizes | |
547 | int iAmpValHighGain[fgkMaxTowers]; | |
548 | int iAmpValLowGain[fgkMaxTowers]; | |
549 | memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain)); | |
550 | memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain)); | |
551 | ||
552 | // also for LED reference | |
553 | int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values | |
554 | memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal)); | |
555 | ||
556 | int gain = 0; // high or low gain | |
557 | ||
558 | // Number of Low and High gain, and LED Ref, channels for this event: | |
559 | int nLowChan = 0; | |
560 | int nHighChan = 0; | |
561 | int nLEDRefChan = 0; | |
562 | ||
563 | int iTowerNum = 0; // array index for regular towers | |
564 | int iRefNum = 0; // array index for LED references | |
565 | ||
566 | // loop first to get the fraction of channels with amplitudes above cut | |
567 | ||
568 | while (in->NextDDL()) { | |
569 | while (in->NextChannel()) { | |
570 | ||
571 | vector<AliCaloBunchInfo> bunchlist; | |
572 | while (in->NextBunch()) { | |
573 | bunchlist.push_back( AliCaloBunchInfo(in->GetStartTimeBin(), in->GetBunchLength(), in->GetSignals() ) ); | |
574 | } | |
575 | if (bunchlist.size() == 0) continue; | |
576 | ||
577 | gain = -1; // init to not valid value | |
578 | //If we're here then we're done with this tower | |
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 | } | |
588 | else { continue; } // don't try to fit TRU.. | |
589 | ||
590 | // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now | |
591 | int arrayPos = in->GetModule(); //The modules are numbered starting from 0 | |
592 | //Debug | |
593 | if (arrayPos < 0 || arrayPos >= fModules) { | |
594 | printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); | |
595 | return kFALSE; | |
596 | } | |
597 | ||
598 | AliCaloFitResults res = fRawAnalyzer->Evaluate( bunchlist, in->GetAltroCFG1(), in->GetAltroCFG2()); | |
599 | if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower | |
600 | // get tower number for AmpVal array | |
601 | iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); | |
602 | ||
603 | if (gain == 0) { | |
604 | // fill amplitude into the array | |
605 | iAmpValLowGain[iTowerNum] = (int) res.GetAmp(); | |
606 | nLowChan++; | |
607 | } | |
608 | else if (gain==1) {//fill the high gain ones | |
609 | // fill amplitude into the array | |
610 | iAmpValHighGain[iTowerNum] = (int) res.GetAmp(); | |
611 | nHighChan++; | |
612 | }//end if gain | |
613 | } // regular tower | |
614 | else if ( in->IsLEDMonData() ) { // LED ref.; | |
615 | // strip # is coded is 'column' in the channel maps | |
616 | iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain); | |
617 | iLEDAmpVal[iRefNum] = (int) res.GetAmp(); | |
618 | nLEDRefChan++; | |
619 | } // end of LED ref | |
620 | ||
621 | } // end while over channel | |
622 | ||
623 | }//end while over DDL's, of input stream | |
624 | ||
625 | in->Reset(); // just in case the next customer forgets to check if the stream was reset.. | |
626 | ||
627 | // now check if it was an LED event, using the LED Reference info per strip | |
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 | } | |
634 | if (fReqFractionAboveAmp) { | |
635 | bool ok = false; | |
636 | if (nHighChan > 0) { | |
637 | ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray); | |
638 | } | |
639 | if (!ok) return false; // skip event | |
640 | } | |
641 | ||
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 | ||
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 | |
658 | fStartTime = Timestamp; | |
659 | } | |
660 | ||
661 | fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr; | |
662 | if (fLatestHour < fHour) { | |
663 | fLatestHour = fHour; | |
664 | } | |
665 | ||
666 | // it is a led event, now fill TTree | |
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 | ||
674 | for(int i=0; i<fModules; i++){ | |
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 | |
678 | for(int k=0; k<fRows; k++){ | |
679 | ||
680 | iTowerNum = GetTowerNum(i, j, k); | |
681 | ||
682 | if(iAmpValHighGain[iTowerNum]) { | |
683 | fAmp = iAmpValHighGain[iTowerNum]; | |
684 | fChannelNum = GetChannelNum(i,j,k,1); | |
685 | fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]); | |
686 | fNHighGain[iTowerNum]++; | |
687 | } | |
688 | if(iAmpValLowGain[iTowerNum]) { | |
689 | fAmp = iAmpValLowGain[iTowerNum]; | |
690 | fChannelNum = GetChannelNum(i,j,k,0); | |
691 | fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]); | |
692 | fNLowGain[iTowerNum]++; | |
693 | } | |
694 | } // rows | |
695 | } // column passed check, and LED Ref for strip passed check (if any) | |
696 | } // columns | |
697 | ||
698 | // also LED refs | |
699 | for(int j=0; j<fLEDRefs; j++){ | |
700 | int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!! | |
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) | |
712 | } | |
713 | ||
714 | } // modules | |
715 | ||
716 | return kTRUE; | |
717 | } | |
718 | ||
719 | //_____________________________________________________________________ | |
720 | Bool_t AliCaloCalibSignal::Save(TString fileName) | |
721 | { | |
722 | //Saves all the TTrees to the designated file | |
723 | ||
724 | TFile destFile(fileName, "recreate"); | |
725 | ||
726 | if (destFile.IsZombie()) { | |
727 | return kFALSE; | |
728 | } | |
729 | ||
730 | destFile.cd(); | |
731 | ||
732 | // save the trees | |
733 | fTreeAmpVsTime->Write(); | |
734 | fTreeLEDAmpVsTime->Write(); | |
735 | if (fUseAverage) { | |
736 | Analyze(); // get the latest and greatest averages | |
737 | fTreeAvgAmpVsTime->Write(); | |
738 | fTreeLEDAvgAmpVsTime->Write(); | |
739 | } | |
740 | ||
741 | destFile.Close(); | |
742 | ||
743 | return kTRUE; | |
744 | } | |
745 | ||
746 | //_____________________________________________________________________ | |
747 | Bool_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 | |
759 | int numProfBins = 0; | |
760 | double timeMin = 0; | |
761 | double timeMax = 0; | |
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)); | |
773 | const Int_t buffersize = 200; | |
774 | char name[buffersize]; // for profile id and title | |
775 | int iTowerNum = 0; | |
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 | ||
781 | iTowerNum = GetTowerNum(i, ic, ir); | |
782 | // high gain | |
783 | if (fNHighGain[iTowerNum] > 0) { | |
784 | fChannelNum = GetChannelNum(i, ic, ir, 1); | |
785 | snprintf(name,buffersize,"profileChan%d", fChannelNum); | |
786 | profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s"); | |
787 | } | |
788 | ||
789 | // same for low gain | |
790 | if (fNLowGain[iTowerNum] > 0) { | |
791 | fChannelNum = GetChannelNum(i, ic, ir, 0); | |
792 | snprintf(name,buffersize,"profileChan%d", fChannelNum); | |
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); | |
811 | } | |
812 | } | |
813 | ||
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 | ||
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) { | |
848 | snprintf(name, buffersize, "profileLEDRef%d", fRefNum); | |
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 | ||
893 | return kTRUE; | |
894 | } | |
895 | ||
896 | //_____________________________________________________________________ | |
897 | Bool_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 | //_____________________________________________________________________ | |
908 | Bool_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 | } |