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 = AliEMCALGeoParams::fgkEMCALCols;
  95     fRows = AliEMCALGeoParams::fgkEMCALRows;
  96     fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
  97     fModules = AliEMCALGeoParams::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 = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::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         // strip # is coded is 'column' in the channel maps
 417         RefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
 418         LEDAmpVal[RefNum] = max - min;
 419       } // end of LED ref
 420
 421     } // end while over channel
 422
 423   }//end while over DDL's, of input stream
 424
 425   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 426
 427   // now check if it was a led event, only use high gain (that should be sufficient)
 428   if (fReqFractionAboveAmp) {
 429     bool ok = false;
 430     if (nHighChan > 0) {
 431       ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
 432     }
 433     if (!ok) return false; // skip event
 434   }
 435
 436   fNAcceptedEvents++; // one more event accepted
 437
 438   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
 439     fStartTime = aliHeader->Get("Timestamp");
 440   }
 441
 442   fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
 443   if (fLatestHour < fHour) {
 444     fLatestHour = fHour;
 445   }
 446
 447   // it is a led event, now fill TTree
 448   for(int i=0; i<fModules; i++){
 449     for(int j=0; j<fColumns; j++){
 450       for(int k=0; k<fRows; k++){
 451
 452         TowerNum = GetTowerNum(i, j, k);
 453
 454         if(AmpValHighGain[TowerNum]) {
 455           fAmp = AmpValHighGain[TowerNum];
 456           fChannelNum = GetChannelNum(i,j,k,1);
 457           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
 458           fNHighGain[TowerNum]++;
 459         }
 460         if(AmpValLowGain[TowerNum]) {
 461           fAmp = AmpValLowGain[TowerNum];
 462           fChannelNum = GetChannelNum(i,j,k,0);
 463           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
 464           fNLowGain[TowerNum]++;
 465         }
 466       } // rows
 467     } // columns
 468
 469     // also LED refs
 470     for(int j=0; j<fLEDRefs; j++){
 471       for (gain=0; gain<2; gain++) {
 472         fRefNum = GetRefNum(i, j, gain);
 473         if (LEDAmpVal[fRefNum]) {
 474           fAmp = LEDAmpVal[fRefNum];
 475           fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
 476           fNRef[fRefNum]++;
 477         }
 478       }
 479     }
 480
 481   } // modules
 482
 483   return kTRUE;
 484 }
 485
 486 //_____________________________________________________________________
 487 Bool_t AliCaloCalibSignal::Save(TString fileName)
 488 {
 489   //Saves all the TTrees to the designated file
 490
 491   TFile destFile(fileName, "recreate");
 492
 493   if (destFile.IsZombie()) {
 494     return kFALSE;
 495   }
 496
 497   destFile.cd();
 498
 499   // save the trees
 500   fTreeAmpVsTime->Write();
 501   fTreeLEDAmpVsTime->Write();
 502   if (fUseAverage) {
 503     Analyze(); // get the latest and greatest averages
 504     fTreeAvgAmpVsTime->Write();
 505     fTreeLEDAvgAmpVsTime->Write();
 506   }
 507
 508   destFile.Close();
 509
 510   return kTRUE;
 511 }
 512
 513 //_____________________________________________________________________
 514 Bool_t AliCaloCalibSignal::Analyze()
 515 {
 516   // Fill the tree holding the average values
 517   if (!fUseAverage) { return kFALSE; }
 518
 519   // Reset the average TTree if Analyze has already been called earlier,
 520   // meaning that the TTree could have been partially filled
 521   if (fTreeAvgAmpVsTime->GetEntries() > 0) {
 522     fTreeAvgAmpVsTime->Reset();
 523   }
 524
 525   //0: setup variables for the TProfile plots that we'll use to do the averages
 526   int numProfBins = 0;
 527   double timeMin = 0;
 528   double timeMax = 0;
 529   if (fSecInAverage > 0) {
 530     numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
 531   }
 532   numProfBins += 2; // add extra buffer : first and last
 533   double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
 534   timeMin = - binSize;
 535   timeMax = timeMin + numProfBins*binSize;
 536
 537   //1: set up TProfiles for the towers that had data
 538   TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
 539   memset(profile, 0, sizeof(profile));
 540
 541   char name[200]; // for profile id and title
 542   int TowerNum = 0;
 543
 544   for (int i = 0; i<fModules; i++) {
 545     for (int ic=0; ic<fColumns; ic++){
 546       for (int ir=0; ir<fRows; ir++) {
 547
 548         TowerNum = GetTowerNum(i, ic, ir);
 549         // high gain
 550         if (fNHighGain[TowerNum] > 0) {
 551           fChannelNum = GetChannelNum(i, ic, ir, 1);
 552           sprintf(name, "profileChan%d", fChannelNum);
 553           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 554         }
 555
 556         // same for low gain
 557         if (fNLowGain[TowerNum] > 0) {
 558           fChannelNum = GetChannelNum(i, ic, ir, 0);
 559           sprintf(name, "profileChan%d", fChannelNum);
 560           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 561         }
 562
 563       } // rows
 564     } // columns
 565   } // modules
 566
 567   //2: fill profiles by looping over tree
 568   // Set addresses for tree-readback also
 569   fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
 570   fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
 571   fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
 572
 573   for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
 574     fTreeAmpVsTime->GetEntry(ient);
 575     if (profile[fChannelNum]) {
 576       // profile should always have been created above, for active channels
 577       profile[fChannelNum]->Fill(fHour, fAmp);
 578     }
 579   }
 580
 581   // re-associating the branch addresses here seems to be needed for OK 'average' storage
 582   fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
 583   fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
 584   fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
 585   fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
 586
 587   //3: fill avg tree by looping over the profiles
 588   for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
 589     if (profile[fChannelNum]) { // profile was created
 590       if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
 591         for(int it=0; it<numProfBins; it++) {
 592           if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
 593             fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
 594             fHour = profile[fChannelNum]->GetBinCenter(it+1);
 595             fRMS = profile[fChannelNum]->GetBinError(it+1);
 596             fTreeAvgAmpVsTime->Fill();
 597           } // some entries for this bin
 598         } // loop over bins
 599       } // some entries for this profile
 600     } // profile exists
 601   } // loop over all possible channels
 602
 603
 604   // and finally, go through same exercise for LED also..
 605
 606   //1: set up TProfiles for the towers that had data
 607   TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
 608   memset(profileLED, 0, sizeof(profileLED));
 609
 610   for (int i = 0; i<fModules; i++) {
 611     for(int j=0; j<fLEDRefs; j++){
 612       for (int gain=0; gain<2; gain++) {
 613         fRefNum = GetRefNum(i, j, gain);
 614         if (fNRef[fRefNum] > 0) {
 615           sprintf(name, "profileLEDRef%d", fRefNum);
 616           profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 617         }
 618       }// gain
 619     }
 620   } // modules
 621
 622   //2: fill profiles by looping over tree
 623   // Set addresses for tree-readback also
 624   fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
 625   fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
 626   fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
 627
 628   for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
 629     fTreeLEDAmpVsTime->GetEntry(ient);
 630     if (profileLED[fRefNum]) {
 631       // profile should always have been created above, for active channels
 632       profileLED[fRefNum]->Fill(fHour, fAmp);
 633     }
 634   }
 635
 636   // re-associating the branch addresses here seems to be needed for OK 'average' storage
 637   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
 638   fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
 639   fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
 640   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
 641
 642   //3: fill avg tree by looping over the profiles
 643   for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
 644     if (profileLED[fRefNum]) { // profile was created
 645       if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
 646         for(int it=0; it<numProfBins; it++) {
 647           if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
 648             fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
 649             fHour = profileLED[fRefNum]->GetBinCenter(it+1);
 650             fRMS = profileLED[fRefNum]->GetBinError(it+1);
 651             fTreeLEDAvgAmpVsTime->Fill();
 652           } // some entries for this bin
 653         } // loop over bins
 654       } // some entries for this profile
 655     } // profile exists
 656   } // loop over all possible channels
 657
 658   // OK, we're done..
 659
 660   return kTRUE;
 661 }