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