Missing class from previous commit
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
CommitLineData
bd83f412 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15/* $Id: AliCaloCalibSignal.cxx $ */
16
17//________________________________________________________________________
18//
19// A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20// It can be created and used a la (ctor):
21/*
22 //Create the object for making the histograms
23 fSignals = new AliCaloCalibSignal( fDetType );
24 // AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
25 fNumModules = fSignals->GetModules();
26*/
27// fed an event:
28// fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
8bcca84c 29// get some info:
30// fSignals->GetXXX..()
bd83f412 31// etc.
32//________________________________________________________________________
33
8bcca84c 34#include "TProfile.h"
bd83f412 35#include "TFile.h"
36
f4fc542c 37#include "AliRawReader.h"
bd83f412 38#include "AliRawEventHeaderBase.h"
39#include "AliCaloRawStream.h"
40
41//The include file
42#include "AliCaloCalibSignal.h"
43
44ClassImp(AliCaloCalibSignal)
45
46using namespace std;
47
8bcca84c 48// variables for TTree filling; not sure if they should be static or not
5e99faca 49static int fChannelNum; // for regular towers
50static int fRefNum; // for LED
8bcca84c 51static double fAmp;
52static double fAvgAmp;
53static double fRMS;
54
bd83f412 55// ctor; initialize everything in order to avoid compiler warnings
56// put some reasonable defaults
57AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
58 TObject(),
59 fDetType(kNone),
60 fColumns(0),
61 fRows(0),
5e99faca 62 fLEDRefs(0),
bd83f412 63 fModules(0),
f4fc542c 64 fCaloString(),
65 fMapping(NULL),
bd83f412 66 fRunNumber(-1),
67 fStartTime(0),
68 fAmpCut(50),
69 fReqFractionAboveAmpCutVal(0.8),
70 fReqFractionAboveAmp(kTRUE),
71 fHour(0),
72 fLatestHour(0),
73 fUseAverage(kTRUE),
74 fSecInAverage(1800),
75 fNEvents(0),
8bcca84c 76 fNAcceptedEvents(0),
77 fTreeAmpVsTime(NULL),
5e99faca 78 fTreeAvgAmpVsTime(NULL),
79 fTreeLEDAmpVsTime(NULL),
80 fTreeLEDAvgAmpVsTime(NULL)
bd83f412 81{
82 //Default constructor. First we set the detector-type related constants.
83 if (detectorType == kPhos) {
84 fColumns = fgkPhosCols;
85 fRows = fgkPhosRows;
5e99faca 86 fLEDRefs = fgkPhosLEDRefs;
bd83f412 87 fModules = fgkPhosModules;
f4fc542c 88 fCaloString = "PHOS";
bd83f412 89 }
90 else {
91 //We'll just trust the enum to keep everything in line, so that if detectorType
92 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
93 //case, if someone intentionally gives another number
94 fColumns = fgkEmCalCols;
95 fRows = fgkEmCalRows;
5e99faca 96 fLEDRefs = fgkEmCalLEDRefs;
bd83f412 97 fModules = fgkEmCalModules;
f4fc542c 98 fCaloString = "EMCAL";
bd83f412 99 }
100
101 fDetType = detectorType;
102
8bcca84c 103 ResetInfo(); // trees and counters
bd83f412 104}
105
106// dtor
107//_____________________________________________________________________
108AliCaloCalibSignal::~AliCaloCalibSignal()
109{
8bcca84c 110 DeleteTrees();
bd83f412 111}
112
113//_____________________________________________________________________
8bcca84c 114void AliCaloCalibSignal::DeleteTrees()
bd83f412 115{
8bcca84c 116 // delete what was created in the ctor (TTrees)
117 if (fTreeAmpVsTime) delete fTreeAmpVsTime;
118 if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
5e99faca 119 if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
120 if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
8bcca84c 121 // and reset pointers
122 fTreeAmpVsTime = NULL;
123 fTreeAvgAmpVsTime = NULL;
5e99faca 124 fTreeLEDAmpVsTime = NULL;
125 fTreeLEDAvgAmpVsTime = NULL;
bd83f412 126
127 return;
128}
129
130// copy ctor
131//_____________________________________________________________________
132AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
133 TObject(sig),
134 fDetType(sig.GetDetectorType()),
135 fColumns(sig.GetColumns()),
136 fRows(sig.GetRows()),
5e99faca 137 fLEDRefs(sig.GetLEDRefs()),
bd83f412 138 fModules(sig.GetModules()),
f4fc542c 139 fCaloString(sig.GetCaloString()),
140 fMapping(NULL), //! note that we are not copying the map info
bd83f412 141 fRunNumber(sig.GetRunNumber()),
142 fStartTime(sig.GetStartTime()),
143 fAmpCut(sig.GetAmpCut()),
144 fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
145 fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
146 fHour(sig.GetHour()),
147 fLatestHour(sig.GetLatestHour()),
148 fUseAverage(sig.GetUseAverage()),
149 fSecInAverage(sig.GetSecInAverage()),
150 fNEvents(sig.GetNEvents()),
8bcca84c 151 fNAcceptedEvents(sig.GetNAcceptedEvents()),
152 fTreeAmpVsTime(NULL),
5e99faca 153 fTreeAvgAmpVsTime(NULL),
154 fTreeLEDAmpVsTime(NULL),
155 fTreeLEDAvgAmpVsTime(NULL)
bd83f412 156{
8bcca84c 157 // also the TTree contents
bd83f412 158 AddInfo(&sig);
159}
160
161// assignment operator; use copy ctor to make life easy..
162//_____________________________________________________________________
163AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
164{
165 // assignment operator; use copy ctor
166 if (&source == this) return *this;
167
168 new (this) AliCaloCalibSignal(source);
169 return *this;
170}
171
172//_____________________________________________________________________
8bcca84c 173void AliCaloCalibSignal::CreateTrees()
bd83f412 174{
8bcca84c 175 // initialize trees
176 // first, regular version
177 fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
178
179 fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
180 fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
181 fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
bd83f412 182
8bcca84c 183 // then, average version
184 fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
185
186 fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
187 fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
188 fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
189 fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
190
5e99faca 191 // then same for LED..
192 fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
193 fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
194 fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
195 fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
196
197 fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
198 fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
199 fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
200 fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
201 fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
202
8bcca84c 203 return;
bd83f412 204}
205
206//_____________________________________________________________________
8bcca84c 207void AliCaloCalibSignal::ResetInfo()
bd83f412 208{
209 Zero(); // set all counters to 0
8bcca84c 210 DeleteTrees(); // delete previous stuff
211 CreateTrees(); // and create some new ones
bd83f412 212 return;
213}
214
215//_____________________________________________________________________
216void AliCaloCalibSignal::Zero()
217{
8bcca84c 218 // set all counters to 0; not cuts etc. though
bd83f412 219 fHour = 0;
220 fLatestHour = 0;
221 fNEvents = 0;
222 fNAcceptedEvents = 0;
8bcca84c 223
5e99faca 224 // Set the number of points for each tower: Amp vs. Time
8bcca84c 225 memset(fNHighGain, 0, sizeof(fNHighGain));
226 memset(fNLowGain, 0, sizeof(fNLowGain));
5e99faca 227 // and LED reference
228 memset(fNRef, 0, sizeof(fNRef));
8bcca84c 229
bd83f412 230 return;
231}
232
233//_____________________________________________________________________
234Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
235{
236 int nAbove = 0;
237
238 int TowerNum = 0;
239 for (int i = 0; i<fModules; i++) {
240 for (int j = 0; j<fColumns; j++) {
241 for (int k = 0; k<fRows; k++) {
242 TowerNum = GetTowerNum(i,j,k);
243 if (AmpVal[TowerNum] > fAmpCut) {
244 nAbove++;
245 }
246 }
247 }
248 }
249
250 double fraction = (1.0*nAbove) / nTotChan;
251
252 if (fraction > fReqFractionAboveAmpCutVal) {
253 return true;
254 }
255 else return false;
256}
257
258//_____________________________________________________________________
259Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
260{
5e99faca 261 // 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.
262
8bcca84c 263 // add info from sig's TTrees to ours..
264 TTree *sigAmp = sig->GetTreeAmpVsTime();
265 TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
266
267 // we could try some merging via TList or what also as a more elegant approach
268 // but I wanted with the stupid/simple and hopefully safe approach of looping
269 // over what we want to add..
270
271 // associate variables for sigAmp and sigAvgAmp:
272 sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
273 sigAmp->SetBranchAddress("fHour",&fHour);
274 sigAmp->SetBranchAddress("fAmp",&fAmp);
275
276 // loop over the trees.. note that since we use the same variables we should not need
277 // to do any assignments between the getting and filling
278 for (int i=0; i<sigAmp->GetEntries(); i++) {
279 sigAmp->GetEntry(i);
280 fTreeAmpVsTime->Fill();
281 }
4fb34e09 282
8bcca84c 283 sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
284 sigAvgAmp->SetBranchAddress("fHour",&fHour);
285 sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
286 sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
bd83f412 287
8bcca84c 288 for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
289 sigAvgAmp->GetEntry(i);
290 fTreeAvgAmpVsTime->Fill();
291 }
bd83f412 292
5e99faca 293 // also LED..
294 TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
295 TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
296
297 // associate variables for sigAmp and sigAvgAmp:
298 sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
299 sigLEDAmp->SetBranchAddress("fHour",&fHour);
300 sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
301
302 // loop over the trees.. note that since we use the same variables we should not need
303 // to do any assignments between the getting and filling
304 for (int i=0; i<sigLEDAmp->GetEntries(); i++) {
305 sigLEDAmp->GetEntry(i);
306 fTreeLEDAmpVsTime->Fill();
307 }
308
309 sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
310 sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
311 sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
312 sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
313
314 for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
315 sigLEDAvgAmp->GetEntry(i);
316 fTreeLEDAvgAmpVsTime->Fill();
317 }
318
319
8bcca84c 320 return kTRUE;//We hopefully succesfully added info from the supplied object
bd83f412 321}
322
f4fc542c 323//_____________________________________________________________________
324Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
325{
326 // if fMapping is NULL the rawstream will crate its own mapping
327 AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
328
329 return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
330}
bd83f412 331
332//_____________________________________________________________________
333Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderBase *aliHeader)
334{
335 // Method to process=analyze one event in the data stream
336 if (!in) return kFALSE; //Return right away if there's a null pointer
337
338 fNEvents++; // one more event
339
5e99faca 340 // use maximum numbers to set array sizes
bd83f412 341 int AmpValHighGain[fgkMaxTowers];
342 int AmpValLowGain[fgkMaxTowers];
bd83f412 343 memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
344 memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
345
5e99faca 346 // also for LED reference
347 int LEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
348 memset(LEDAmpVal, 0, sizeof(LEDAmpVal));
349
8cdb1ddb 350 int sample, isample = 0; //The sample temp, and the sample number in current event.
bd83f412 351 int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
5e99faca 352 int gain = 0; // high or low gain
bd83f412 353
354 // Number of Low and High gain channels for this event:
355 int nLowChan = 0;
356 int nHighChan = 0;
357
5e99faca 358 int TowerNum = 0; // array index for regular towers
359 int RefNum = 0; // array index for LED references
bd83f412 360
361 // loop first to get the fraction of channels with amplitudes above cut
362 while (in->Next()) {
363 sample = in->GetSignal(); //Get the adc signal
364 if (sample < min) min = sample;
365 if (sample > max) max = sample;
8cdb1ddb 366 isample++;
367 if ( isample >= in->GetTimeLength()) {
bd83f412 368 //If we're here then we're done with this tower
5e99faca 369 if ( in->IsLowGain() ) {
370 gain = 0;
371 }
372 else if ( in->IsHighGain() ) {
373 gain = 1;
374 }
375 else if ( in->IsLEDMonData() ) {
376 gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
377 }
378
bd83f412 379 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
bd83f412 380 //Debug
381 if (arrayPos < 0 || arrayPos >= fModules) {
8bcca84c 382 printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
5e99faca 383 return kFALSE;
bd83f412 384 }
385
5e99faca 386 if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
387 // get tower number for AmpVal array
388 TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
389
390 if (gain == 0) {
391 // fill amplitude into the array
392 AmpValLowGain[TowerNum] = max - min;
393 nLowChan++;
394 }
395 else if (gain==1) {//fill the high gain ones
396 // fill amplitude into the array
397 AmpValHighGain[TowerNum] = max - min;
398 nHighChan++;
399 }//end if gain
400 } // regular tower
401 else if ( in->IsLEDMonData() ) { // LED ref.
402 RefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
403 LEDAmpVal[RefNum] = max - min;
404 } // end of LED ref
bd83f412 405
406 max = fgkSampleMin; min = fgkSampleMax;
8cdb1ddb 407 isample = 0;
bd83f412 408
409 }//End if end of tower
410
411 }//end while, of stream
412
413 // now check if it was a led event, only use high gain (that should be sufficient)
414 if (fReqFractionAboveAmp) {
452729ca 415 bool ok = false;
416 if (nHighChan > 0) {
417 ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
418 }
bd83f412 419 if (!ok) return false; // skip event
420 }
421
422 fNAcceptedEvents++; // one more event accepted
423
424 if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
425 fStartTime = aliHeader->Get("Timestamp");
426 }
427
428 fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
429 if (fLatestHour < fHour) {
430 fLatestHour = fHour;
431 }
432
8bcca84c 433 // it is a led event, now fill TTree
434 for(int i=0; i<fModules; i++){
435 for(int j=0; j<fColumns; j++){
436 for(int k=0; k<fRows; k++){
bd83f412 437
438 TowerNum = GetTowerNum(i, j, k);
439
440 if(AmpValHighGain[TowerNum]) {
8bcca84c 441 fAmp = AmpValHighGain[TowerNum];
442 fChannelNum = GetChannelNum(i,j,k,1);
443 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
444 fNHighGain[TowerNum]++;
bd83f412 445 }
446 if(AmpValLowGain[TowerNum]) {
8bcca84c 447 fAmp = AmpValLowGain[TowerNum];
448 fChannelNum = GetChannelNum(i,j,k,0);
449 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
450 fNLowGain[TowerNum]++;
bd83f412 451 }
5e99faca 452 } // rows
453 } // columns
454
455 // also LED refs
456 for(int j=0; j<fLEDRefs; j++){
457 for (gain=0; gain<2; gain++) {
458 fRefNum = GetRefNum(i, j, gain);
1bfe2e7a 459 if (LEDAmpVal[fRefNum]) {
460 fAmp = LEDAmpVal[fRefNum];
5e99faca 461 fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
462 fNRef[fRefNum]++;
463 }
bd83f412 464 }
465 }
5e99faca 466
467 } // modules
bd83f412 468
469 return kTRUE;
470}
471
472//_____________________________________________________________________
8bcca84c 473Bool_t AliCaloCalibSignal::Save(TString fileName)
bd83f412 474{
8bcca84c 475 //Saves all the TTrees to the designated file
bd83f412 476
477 TFile destFile(fileName, "recreate");
478
479 if (destFile.IsZombie()) {
480 return kFALSE;
481 }
482
483 destFile.cd();
484
8bcca84c 485 // save the trees
486 fTreeAmpVsTime->Write();
1bfe2e7a 487 fTreeLEDAmpVsTime->Write();
8bcca84c 488 if (fUseAverage) {
489 Analyze(); // get the latest and greatest averages
490 fTreeAvgAmpVsTime->Write();
1bfe2e7a 491 fTreeLEDAvgAmpVsTime->Write();
8bcca84c 492 }
493
494 destFile.Close();
495
496 return kTRUE;
497}
498
499//_____________________________________________________________________
500Bool_t AliCaloCalibSignal::Analyze()
501{
502 // Fill the tree holding the average values
503 if (!fUseAverage) { return kFALSE; }
504
505 // Reset the average TTree if Analyze has already been called earlier,
506 // meaning that the TTree could have been partially filled
507 if (fTreeAvgAmpVsTime->GetEntries() > 0) {
508 fTreeAvgAmpVsTime->Reset();
509 }
510
511 //0: setup variables for the TProfile plots that we'll use to do the averages
bd83f412 512 int numProfBins = 0;
513 double timeMin = 0;
514 double timeMax = 0;
8bcca84c 515 if (fSecInAverage > 0) {
516 numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
517 }
518 numProfBins += 2; // add extra buffer : first and last
519 double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
520 timeMin = - binSize;
521 timeMax = timeMin + numProfBins*binSize;
522
523 //1: set up TProfiles for the towers that had data
524 TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
525 memset(profile, 0, sizeof(profile));
526
527 char name[200]; // for profile id and title
528 int TowerNum = 0;
529
530 for (int i = 0; i<fModules; i++) {
531 for (int ic=0; ic<fColumns; ic++){
532 for (int ir=0; ir<fRows; ir++) {
533
534 TowerNum = GetTowerNum(i, ic, ir);
535 // high gain
536 if (fNHighGain[TowerNum] > 0) {
537 fChannelNum = GetChannelNum(i, ic, ir, 1);
538 sprintf(name, "profileChan%d", fChannelNum);
539 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
540 }
541
542 // same for low gain
543 if (fNLowGain[TowerNum] > 0) {
544 fChannelNum = GetChannelNum(i, ic, ir, 0);
545 sprintf(name, "profileChan%d", fChannelNum);
546 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
547 }
548
549 } // rows
550 } // columns
551 } // modules
552
553 //2: fill profiles by looping over tree
554 // Set addresses for tree-readback also
555 fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
556 fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
557 fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
558
559 for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
560 fTreeAmpVsTime->GetEntry(ient);
561 if (profile[fChannelNum]) {
562 // profile should always have been created above, for active channels
563 profile[fChannelNum]->Fill(fHour, fAmp);
bd83f412 564 }
bd83f412 565 }
566
8bcca84c 567 // re-associating the branch addresses here seems to be needed for OK 'average' storage
568 fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
569 fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
570 fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
571 fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
572
573 //3: fill avg tree by looping over the profiles
574 for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
575 if (profile[fChannelNum]) { // profile was created
576 if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
577 for(int it=0; it<numProfBins; it++) {
578 if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
579 fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
580 fHour = profile[fChannelNum]->GetBinCenter(it+1);
581 fRMS = profile[fChannelNum]->GetBinError(it+1);
582 fTreeAvgAmpVsTime->Fill();
583 } // some entries for this bin
584 } // loop over bins
585 } // some entries for this profile
586 } // profile exists
587 } // loop over all possible channels
588
5e99faca 589
590 // and finally, go through same exercise for LED also..
591
592 //1: set up TProfiles for the towers that had data
593 TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
594 memset(profileLED, 0, sizeof(profileLED));
595
596 for (int i = 0; i<fModules; i++) {
597 for(int j=0; j<fLEDRefs; j++){
598 for (int gain=0; gain<2; gain++) {
599 fRefNum = GetRefNum(i, j, gain);
600 if (fNRef[fRefNum] > 0) {
601 sprintf(name, "profileLEDRef%d", fRefNum);
602 profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
603 }
604 }// gain
605 }
606 } // modules
607
608 //2: fill profiles by looping over tree
609 // Set addresses for tree-readback also
610 fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
611 fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
612 fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
613
614 for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
615 fTreeLEDAmpVsTime->GetEntry(ient);
616 if (profileLED[fRefNum]) {
617 // profile should always have been created above, for active channels
618 profileLED[fRefNum]->Fill(fHour, fAmp);
619 }
620 }
621
622 // re-associating the branch addresses here seems to be needed for OK 'average' storage
623 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
624 fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
625 fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
626 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
627
628 //3: fill avg tree by looping over the profiles
629 for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
630 if (profileLED[fRefNum]) { // profile was created
631 if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
632 for(int it=0; it<numProfBins; it++) {
633 if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
634 fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
635 fHour = profileLED[fRefNum]->GetBinCenter(it+1);
636 fRMS = profileLED[fRefNum]->GetBinError(it+1);
637 fTreeLEDAvgAmpVsTime->Fill();
638 } // some entries for this bin
639 } // loop over bins
640 } // some entries for this profile
641 } // profile exists
642 } // loop over all possible channels
643
644 // OK, we're done..
645
bd83f412 646 return kTRUE;
647}