use the new AliCaloRawStreamV3 decoder (from Yuri), based on the new AliAltroRawStrea...
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
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
34 #include "TProfile.h"
35 #include "TFile.h"
36
37 #include "AliRawReader.h"
38 #include "AliRawEventHeaderBase.h"
39 #include "AliCaloRawStreamV3.h"
40
41 //The include file
42 #include "AliCaloCalibSignal.h"
43
44 ClassImp(AliCaloCalibSignal)
45
46 using namespace std;
47
48 // variables for TTree filling; not sure if they should be static or not
49 static int fChannelNum; // for regular towers
50 static int fRefNum; // for LED
51 static double fAmp;
52 static double fAvgAmp;
53 static double fRMS;
54
55 // ctor; initialize everything in order to avoid compiler warnings
56 // put some reasonable defaults
57 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :  
58   TObject(),
59   fDetType(kNone),
60   fColumns(0),
61   fRows(0),
62   fLEDRefs(0),
63   fModules(0),
64   fCaloString(),
65   fMapping(NULL),
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),
76   fNAcceptedEvents(0),
77   fTreeAmpVsTime(NULL),
78   fTreeAvgAmpVsTime(NULL),
79   fTreeLEDAmpVsTime(NULL),
80   fTreeLEDAvgAmpVsTime(NULL)
81 {
82   //Default constructor. First we set the detector-type related constants.
83   if (detectorType == kPhos) {
84     fColumns = fgkPhosCols;
85     fRows = fgkPhosRows;
86     fLEDRefs = fgkPhosLEDRefs;
87     fModules = fgkPhosModules;
88     fCaloString = "PHOS";
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;
96     fLEDRefs = fgkEmCalLEDRefs;
97     fModules = fgkEmCalModules;
98     fCaloString = "EMCAL";
99   }
100
101   fDetType = detectorType;
102
103   ResetInfo(); // trees and counters
104
105
106 // dtor
107 //_____________________________________________________________________
108 AliCaloCalibSignal::~AliCaloCalibSignal()
109 {
110   DeleteTrees();
111 }
112
113 //_____________________________________________________________________
114 void AliCaloCalibSignal::DeleteTrees()
115 {
116   // delete what was created in the ctor (TTrees)
117   if (fTreeAmpVsTime) delete fTreeAmpVsTime;
118   if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
119   if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
120   if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
121   // and reset pointers
122   fTreeAmpVsTime = NULL;
123   fTreeAvgAmpVsTime = NULL;
124   fTreeLEDAmpVsTime = NULL;
125   fTreeLEDAvgAmpVsTime = NULL;
126
127   return;
128 }
129
130 // copy ctor
131 //_____________________________________________________________________
132 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
133   TObject(sig),
134   fDetType(sig.GetDetectorType()),
135   fColumns(sig.GetColumns()),
136   fRows(sig.GetRows()),
137   fLEDRefs(sig.GetLEDRefs()),
138   fModules(sig.GetModules()),
139   fCaloString(sig.GetCaloString()),
140   fMapping(NULL), //! note that we are not copying the map info
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()),
151   fNAcceptedEvents(sig.GetNAcceptedEvents()),
152   fTreeAmpVsTime(NULL),
153   fTreeAvgAmpVsTime(NULL),
154   fTreeLEDAmpVsTime(NULL),
155   fTreeLEDAvgAmpVsTime(NULL)
156 {
157   // also the TTree contents
158   AddInfo(&sig);
159 }
160
161 // assignment operator; use copy ctor to make life easy..
162 //_____________________________________________________________________
163 AliCaloCalibSignal& 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 //_____________________________________________________________________
173 void AliCaloCalibSignal::CreateTrees()
174 {
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");
182
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
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
203   return;
204 }
205
206 //_____________________________________________________________________
207 void AliCaloCalibSignal::ResetInfo()
208 {
209   Zero(); // set all counters to 0
210   DeleteTrees(); // delete previous stuff
211   CreateTrees(); // and create some new ones
212   return;
213 }
214
215 //_____________________________________________________________________
216 void AliCaloCalibSignal::Zero()
217 {
218   // set all counters to 0; not cuts etc. though
219   fHour = 0;
220   fLatestHour = 0;
221   fNEvents = 0;
222   fNAcceptedEvents = 0;
223
224   // Set the number of points for each tower: Amp vs. Time 
225   memset(fNHighGain, 0, sizeof(fNHighGain));
226   memset(fNLowGain, 0, sizeof(fNLowGain));
227   // and LED reference
228   memset(fNRef, 0, sizeof(fNRef));
229
230   return;
231 }
232
233 //_____________________________________________________________________
234 Bool_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 //_____________________________________________________________________
259 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
260 {
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
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   }
282
283   sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
284   sigAvgAmp->SetBranchAddress("fHour",&fHour);
285   sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
286   sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
287
288   for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
289     sigAvgAmp->GetEntry(i);
290     fTreeAvgAmpVsTime->Fill();
291   }
292
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
320   return kTRUE;//We hopefully succesfully added info from the supplied object
321 }
322
323 //_____________________________________________________________________
324 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
325 {
326   // if fMapping is NULL the rawstream will crate its own mapping
327   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
328
329   return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
330 }
331
332 //_____________________________________________________________________
333 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *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
340   // use maximum numbers to set array sizes
341   int AmpValHighGain[fgkMaxTowers];
342   int AmpValLowGain[fgkMaxTowers];
343   memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
344   memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
345
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
350   int sample; // temporary value
351   int gain = 0; // high or low gain
352   
353   // Number of Low and High gain channels for this event:
354   int nLowChan = 0; 
355   int nHighChan = 0; 
356
357   int TowerNum = 0; // array index for regular towers
358   int RefNum = 0; // array index for LED references
359
360   // loop first to get the fraction of channels with amplitudes above cut
361
362   while (in->NextDDL()) {
363     while (in->NextChannel()) {
364
365       // counters
366       int max = fgkSampleMin, min = fgkSampleMax; // min and max sample values
367       
368       while (in->NextBunch()) {
369         const UShort_t *sig = in->GetSignals();
370         for (Int_t i = 0; i < in->GetBunchLength(); i++) {
371           sample = sig[i];
372
373           // check if it's a min or max value
374           if (sample < min) min = sample;
375           if (sample > max) max = sample;
376
377         } // loop over samples in bunch
378       } // loop over bunches
379
380       gain = -1; // init to not valid value
381       //If we're here then we're done with this tower
382       if ( in->IsLowGain() ) {
383         gain = 0;
384       }
385       else if ( in->IsHighGain() ) {
386         gain = 1;
387       }
388       else if ( in->IsLEDMonData() ) {
389         gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
390       }
391
392       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
393       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
394       //Debug
395       if (arrayPos < 0 || arrayPos >= fModules) {
396         printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); 
397         return kFALSE;
398       }
399
400       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
401         // get tower number for AmpVal array
402         TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
403
404         if (gain == 0) {
405           // fill amplitude into the array         
406           AmpValLowGain[TowerNum] = max - min;
407           nLowChan++;
408         } 
409         else if (gain==1) {//fill the high gain ones
410           // fill amplitude into the array
411           AmpValHighGain[TowerNum] = max - min;
412           nHighChan++;
413         }//end if gain
414       } // regular tower
415       else if ( in->IsLEDMonData() ) { // LED ref.
416         RefNum = GetRefNum(arrayPos, in->GetColumn(), gain); 
417         LEDAmpVal[RefNum] = max - min;
418       } // end of LED ref
419       
420     } // end while over channel 
421    
422   }//end while over DDL's, of input stream
423   
424   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
425
426   // now check if it was a led event, only use high gain (that should be sufficient)
427   if (fReqFractionAboveAmp) {
428     bool ok = false;
429     if (nHighChan > 0) { 
430       ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan); 
431     }
432     if (!ok) return false; // skip event
433   }
434
435   fNAcceptedEvents++; // one more event accepted
436
437   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
438     fStartTime = aliHeader->Get("Timestamp");
439   }
440
441   fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
442   if (fLatestHour < fHour) {
443     fLatestHour = fHour; 
444   }
445   
446   // it is a led event, now fill TTree
447   for(int i=0; i<fModules; i++){
448     for(int j=0; j<fColumns; j++){
449       for(int k=0; k<fRows; k++){
450         
451         TowerNum = GetTowerNum(i, j, k); 
452
453         if(AmpValHighGain[TowerNum]) {
454           fAmp = AmpValHighGain[TowerNum];
455           fChannelNum = GetChannelNum(i,j,k,1);
456           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
457           fNHighGain[TowerNum]++;
458         }
459         if(AmpValLowGain[TowerNum]) {
460           fAmp = AmpValLowGain[TowerNum];
461           fChannelNum = GetChannelNum(i,j,k,0);
462           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
463           fNLowGain[TowerNum]++;
464         }
465       } // rows
466     } // columns
467
468     // also LED refs
469     for(int j=0; j<fLEDRefs; j++){
470       for (gain=0; gain<2; gain++) {
471         fRefNum = GetRefNum(i, j, gain); 
472         if (LEDAmpVal[fRefNum]) {
473           fAmp = LEDAmpVal[fRefNum];
474           fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
475           fNRef[fRefNum]++;
476         }
477       }
478     }
479
480   } // modules
481   
482   return kTRUE;
483 }
484
485 //_____________________________________________________________________
486 Bool_t AliCaloCalibSignal::Save(TString fileName)
487 {
488   //Saves all the TTrees to the designated file
489   
490   TFile destFile(fileName, "recreate");
491   
492   if (destFile.IsZombie()) {
493     return kFALSE;
494   }
495   
496   destFile.cd();
497
498   // save the trees
499   fTreeAmpVsTime->Write();
500   fTreeLEDAmpVsTime->Write();
501   if (fUseAverage) { 
502     Analyze(); // get the latest and greatest averages
503     fTreeAvgAmpVsTime->Write();
504     fTreeLEDAvgAmpVsTime->Write();
505   }
506
507   destFile.Close();
508   
509   return kTRUE;
510 }
511
512 //_____________________________________________________________________
513 Bool_t AliCaloCalibSignal::Analyze()
514 {
515   // Fill the tree holding the average values
516   if (!fUseAverage) { return kFALSE; }
517
518   // Reset the average TTree if Analyze has already been called earlier,
519   // meaning that the TTree could have been partially filled
520   if (fTreeAvgAmpVsTime->GetEntries() > 0) {
521     fTreeAvgAmpVsTime->Reset();
522   }
523
524   //0: setup variables for the TProfile plots that we'll use to do the averages
525   int numProfBins = 0;
526   double timeMin = 0;
527   double timeMax = 0;
528   if (fSecInAverage > 0) {
529     numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
530   }
531   numProfBins += 2; // add extra buffer : first and last
532   double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
533   timeMin = - binSize;
534   timeMax = timeMin + numProfBins*binSize;
535
536   //1: set up TProfiles for the towers that had data
537   TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
538   memset(profile, 0, sizeof(profile));
539
540   char name[200]; // for profile id and title
541   int TowerNum = 0;
542
543   for (int i = 0; i<fModules; i++) {
544     for (int ic=0; ic<fColumns; ic++){
545       for (int ir=0; ir<fRows; ir++) {
546
547         TowerNum = GetTowerNum(i, ic, ir);
548         // high gain
549         if (fNHighGain[TowerNum] > 0) {
550           fChannelNum = GetChannelNum(i, ic, ir, 1); 
551           sprintf(name, "profileChan%d", fChannelNum);
552           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
553         }
554
555         // same for low gain
556         if (fNLowGain[TowerNum] > 0) {
557           fChannelNum = GetChannelNum(i, ic, ir, 0); 
558           sprintf(name, "profileChan%d", fChannelNum);
559           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
560         }
561
562       } // rows
563     } // columns
564   } // modules
565
566   //2: fill profiles by looping over tree
567   // Set addresses for tree-readback also
568   fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
569   fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
570   fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
571
572   for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
573     fTreeAmpVsTime->GetEntry(ient);
574     if (profile[fChannelNum]) { 
575       // profile should always have been created above, for active channels
576       profile[fChannelNum]->Fill(fHour, fAmp);
577     }
578   }
579
580   // re-associating the branch addresses here seems to be needed for OK 'average' storage           
581   fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
582   fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
583   fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
584   fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
585
586   //3: fill avg tree by looping over the profiles
587   for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
588     if (profile[fChannelNum]) { // profile was created
589       if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
590         for(int it=0; it<numProfBins; it++) {
591           if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
592             fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
593             fHour = profile[fChannelNum]->GetBinCenter(it+1);
594             fRMS = profile[fChannelNum]->GetBinError(it+1);
595             fTreeAvgAmpVsTime->Fill();
596           } // some entries for this bin
597         } // loop over bins
598       } // some entries for this profile
599     } // profile exists  
600   } // loop over all possible channels
601
602
603   // and finally, go through same exercise for LED also.. 
604
605   //1: set up TProfiles for the towers that had data
606   TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
607   memset(profileLED, 0, sizeof(profileLED));
608
609   for (int i = 0; i<fModules; i++) {
610     for(int j=0; j<fLEDRefs; j++){
611       for (int gain=0; gain<2; gain++) {
612         fRefNum = GetRefNum(i, j, gain);
613         if (fNRef[fRefNum] > 0) { 
614           sprintf(name, "profileLEDRef%d", fRefNum);
615           profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
616         } 
617       }// gain
618     } 
619   } // modules
620
621   //2: fill profiles by looping over tree
622   // Set addresses for tree-readback also
623   fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
624   fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
625   fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
626
627   for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
628     fTreeLEDAmpVsTime->GetEntry(ient);
629     if (profileLED[fRefNum]) { 
630       // profile should always have been created above, for active channels
631       profileLED[fRefNum]->Fill(fHour, fAmp);
632     }
633   }
634
635   // re-associating the branch addresses here seems to be needed for OK 'average' storage           
636   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
637   fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
638   fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
639   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
640
641   //3: fill avg tree by looping over the profiles
642   for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
643     if (profileLED[fRefNum]) { // profile was created
644       if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
645         for(int it=0; it<numProfBins; it++) {
646           if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
647             fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
648             fHour = profileLED[fRefNum]->GetBinCenter(it+1);
649             fRMS = profileLED[fRefNum]->GetBinError(it+1);
650             fTreeLEDAvgAmpVsTime->Fill();
651           } // some entries for this bin
652         } // loop over bins
653       } // some entries for this profile
654     } // profile exists  
655   } // loop over all possible channels
656
657   // OK, we're done..
658
659   return kTRUE;
660 }