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 #include <string>
  34 #include <sstream>
  35 #include <fstream>
  36
  37 #include "TProfile.h"
  38 #include "TFile.h"
  39
  40 #include "AliRawReader.h"
  41 #include "AliCaloRawStreamV3.h"
  42
  43 //The include file
  44 #include "AliCaloCalibSignal.h"
  45
  46 ClassImp(AliCaloCalibSignal)
  47
  48 using namespace std;
  49
  50 // variables for TTree filling; not sure if they should be static or not
  51 static int fChannelNum; // for regular towers
  52 static int fRefNum; // for LED
  53 static double fAmp;
  54 static double fAvgAmp;
  55 static double fRMS;
  56
  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),
  64   fLEDRefs(0),
  65   fModules(0),
  66   fCaloString(),
  67   fMapping(NULL),
  68   fRunNumber(-1),
  69   fStartTime(0),
  70   fAmpCut(40), // min. 40 ADC counts as default
  71   fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
  72   fReqFractionAboveAmp(kTRUE),
  73   fHour(0),
  74   fLatestHour(0),
  75   fUseAverage(kTRUE),
  76   fSecInAverage(1800),
  77   fNEvents(0),
  78   fNAcceptedEvents(0),
  79   fTreeAmpVsTime(NULL),
  80   fTreeAvgAmpVsTime(NULL),
  81   fTreeLEDAmpVsTime(NULL),
  82   fTreeLEDAvgAmpVsTime(NULL)
  83 {
  84   //Default constructor. First we set the detector-type related constants.
  85   if (detectorType == kPhos) {
  86     fColumns = fgkPhosCols;
  87     fRows = fgkPhosRows;
  88     fLEDRefs = fgkPhosLEDRefs;
  89     fModules = fgkPhosModules;
  90     fCaloString = "PHOS";
  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
  96     fColumns = AliEMCALGeoParams::fgkEMCALCols;
  97     fRows = AliEMCALGeoParams::fgkEMCALRows;
  98     fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
  99     fModules = AliEMCALGeoParams::fgkEMCALModules;
 100     fCaloString = "EMCAL";
 101   }
 102
 103   fDetType = detectorType;
 104
 105   ResetInfo(); // trees and counters
 106 }
 107
 108 // dtor
 109 //_____________________________________________________________________
 110 AliCaloCalibSignal::~AliCaloCalibSignal()
 111 {
 112   DeleteTrees();
 113 }
 114
 115 //_____________________________________________________________________
 116 void AliCaloCalibSignal::DeleteTrees()
 117 {
 118   // delete what was created in the ctor (TTrees)
 119   if (fTreeAmpVsTime) delete fTreeAmpVsTime;
 120   if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
 121   if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
 122   if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
 123   // and reset pointers
 124   fTreeAmpVsTime = NULL;
 125   fTreeAvgAmpVsTime = NULL;
 126   fTreeLEDAmpVsTime = NULL;
 127   fTreeLEDAvgAmpVsTime = NULL;
 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()),
 139   fLEDRefs(sig.GetLEDRefs()),
 140   fModules(sig.GetModules()),
 141   fCaloString(sig.GetCaloString()),
 142   fMapping(NULL), //! note that we are not copying the map info
 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()),
 153   fNAcceptedEvents(sig.GetNAcceptedEvents()),
 154   fTreeAmpVsTime(NULL),
 155   fTreeAvgAmpVsTime(NULL),
 156   fTreeLEDAmpVsTime(NULL),
 157   fTreeLEDAvgAmpVsTime(NULL)
 158 {
 159   // also the TTree contents
 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 //_____________________________________________________________________
 175 void AliCaloCalibSignal::CreateTrees()
 176 {
 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");
 184
 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
 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
 205   return;
 206 }
 207
 208 //_____________________________________________________________________
 209 void AliCaloCalibSignal::ResetInfo()
 210 {
 211   Zero(); // set all counters to 0
 212   DeleteTrees(); // delete previous stuff
 213   CreateTrees(); // and create some new ones
 214   return;
 215 }
 216
 217 //_____________________________________________________________________
 218 void AliCaloCalibSignal::Zero()
 219 {
 220   // set all counters to 0; not cuts etc. though
 221   fHour = 0;
 222   fLatestHour = 0;
 223   fNEvents = 0;
 224   fNAcceptedEvents = 0;
 225
 226   // Set the number of points for each tower: Amp vs. Time
 227   memset(fNHighGain, 0, sizeof(fNHighGain));
 228   memset(fNLowGain, 0, sizeof(fNLowGain));
 229   // and LED reference
 230   memset(fNRef, 0, sizeof(fNRef));
 231
 232   return;
 233 }
 234
 235 //_____________________________________________________________________
 236 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal,
 237                                                  int resultArray[])
 238 {
 239   Bool_t returnCode = false;
 240
 241   int TowerNum = 0;
 242   double fraction = 0;
 243   for (int i = 0; i<fModules; i++) {
 244     for (int j = 0; j<fColumns; j++) {
 245       int nAbove = 0;
 246       for (int k = 0; k<fRows; k++) {
 247         TowerNum = GetTowerNum(i,j,k);
 248         if (AmpVal[TowerNum] > fAmpCut) {
 249           nAbove++;
 250         }
 251       }
 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       }
 264     }
 265   } // modules loop
 266
 267   return returnCode;
 268 }
 269
 270 // Parameter/cut handling
 271 //_____________________________________________________________________
 272 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
 273 {
 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
 299       string key_value;
 300       s >> key_value;
 301
 302       // check stream status
 303       if( s.rdstate() & ios::failbit ) break;
 304
 305       // skip rest of line if comments found
 306       if( key_value.substr( 0, 2 ) == "//" ) break;
 307
 308       // look for "::" in key_value pair
 309       size_t position = key_value.find( delimitor );
 310       if( position == string::npos ) {
 311         printf("wrong format for key::value pair: %s\n", key_value.c_str());
 312       }
 313
 314       // split key_value pair
 315       string key( key_value.substr( 0, position ) );
 316       string value( key_value.substr( position+delimitor.size(),
 317                                       key_value.size()-delimitor.size() ) );
 318
 319       // check value does not contain a new delimitor
 320       if( value.find( delimitor ) != string::npos ) {
 321         printf("wrong format for key::value pair: %s\n", key_value.c_str());
 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)
 353 {
 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
 365 //_____________________________________________________________________
 366 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
 367 {
 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
 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   }
 389
 390   sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
 391   sigAvgAmp->SetBranchAddress("fHour",&fHour);
 392   sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
 393   sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
 394
 395   for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
 396     sigAvgAmp->GetEntry(i);
 397     fTreeAvgAmpVsTime->Fill();
 398   }
 399
 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
 427   return kTRUE;//We hopefully succesfully added info from the supplied object
 428 }
 429
 430 //_____________________________________________________________________
 431 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
 432 {
 433   // if fMapping is NULL the rawstream will crate its own mapping
 434   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 435
 436   return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
 437 }
 438
 439 //_____________________________________________________________________
 440 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
 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
 447   // use maximum numbers to set array sizes
 448   int AmpValHighGain[fgkMaxTowers];
 449   int AmpValLowGain[fgkMaxTowers];
 450   memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
 451   memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
 452
 453   // also for LED reference
 454   int LEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
 455   memset(LEDAmpVal, 0, sizeof(LEDAmpVal));
 456
 457   int sample; // temporary value
 458   int gain = 0; // high or low gain
 459
 460   // Number of Low and High gain channels for this event:
 461   int nLowChan = 0;
 462   int nHighChan = 0;
 463
 464   int TowerNum = 0; // array index for regular towers
 465   int RefNum = 0; // array index for LED references
 466
 467   // loop first to get the fraction of channels with amplitudes above cut
 468
 469   while (in->NextDDL()) {
 470     while (in->NextChannel()) {
 471
 472       // counters
 473       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
 474       int nsamples = 0;
 475
 476       while (in->NextBunch()) {
 477         const UShort_t *sig = in->GetSignals();
 478         nsamples += in->GetBunchLength();
 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
 489       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
 490
 491       gain = -1; // init to not valid value
 492       //If we're here then we're done with this tower
 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
 503       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
 504       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 505       //Debug
 506       if (arrayPos < 0 || arrayPos >= fModules) {
 507         printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
 508         return kFALSE;
 509       }
 510
 511       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
 512         // get tower number for AmpVal array
 513         TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
 514
 515         if (gain == 0) {
 516           // fill amplitude into the array
 517           AmpValLowGain[TowerNum] = max - min;
 518           nLowChan++;
 519         }
 520         else if (gain==1) {//fill the high gain ones
 521           // fill amplitude into the array
 522           AmpValHighGain[TowerNum] = max - min;
 523           nHighChan++;
 524         }//end if gain
 525       } // regular tower
 526       else if ( in->IsLEDMonData() ) { // LED ref.;
 527         // strip # is coded is 'column' in the channel maps
 528         RefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
 529         LEDAmpVal[RefNum] = max - min;
 530       } // end of LED ref
 531
 532       } // nsamples>0 check, some data found for this channel; not only trailer/header
 533     } // end while over channel
 534
 535   }//end while over DDL's, of input stream
 536
 537   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 538
 539   // now check if it was a led event, only use high gain (that should be sufficient)
 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
 547   if (fReqFractionAboveAmp) {
 548     bool ok = false;
 549     if (nHighChan > 0) {
 550       ok = CheckFractionAboveAmp(AmpValHighGain, checkResultArray);
 551     }
 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
 558     fStartTime = Timestamp;
 559   }
 560
 561   fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
 562   if (fLatestHour < fHour) {
 563     fLatestHour = fHour;
 564   }
 565
 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++){
 569       if (checkResultArray[i*fColumns + j] > 0) { // column passed check
 570       for(int k=0; k<fRows; k++){
 571
 572         TowerNum = GetTowerNum(i, j, k);
 573
 574         if(AmpValHighGain[TowerNum]) {
 575           fAmp = AmpValHighGain[TowerNum];
 576           fChannelNum = GetChannelNum(i,j,k,1);
 577           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
 578           fNHighGain[TowerNum]++;
 579         }
 580         if(AmpValLowGain[TowerNum]) {
 581           fAmp = AmpValLowGain[TowerNum];
 582           fChannelNum = GetChannelNum(i,j,k,0);
 583           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
 584           fNLowGain[TowerNum]++;
 585         }
 586       } // rows
 587       } // column passed check
 588     } // columns
 589
 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
 597     // also LED refs
 598     for(int j=0; j<fLEDRefs; j++){
 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
 601       for (gain=0; gain<2; gain++) {
 602         fRefNum = GetRefNum(i, j, gain);
 603         if (LEDAmpVal[fRefNum]) {
 604           fAmp = LEDAmpVal[fRefNum];
 605           fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
 606           fNRef[fRefNum]++;
 607         }
 608       } // gain
 609       } // at least one column in strip passed check
 610     }
 611
 612   } // modules
 613
 614   return kTRUE;
 615 }
 616
 617 //_____________________________________________________________________
 618 Bool_t AliCaloCalibSignal::Save(TString fileName)
 619 {
 620   //Saves all the TTrees to the designated file
 621
 622   TFile destFile(fileName, "recreate");
 623
 624   if (destFile.IsZombie()) {
 625     return kFALSE;
 626   }
 627
 628   destFile.cd();
 629
 630   // save the trees
 631   fTreeAmpVsTime->Write();
 632   fTreeLEDAmpVsTime->Write();
 633   if (fUseAverage) {
 634     Analyze(); // get the latest and greatest averages
 635     fTreeAvgAmpVsTime->Write();
 636     fTreeLEDAvgAmpVsTime->Write();
 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
 657   int numProfBins = 0;
 658   double timeMin = 0;
 659   double timeMax = 0;
 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
 673   int TowerNum = 0;
 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
 679         TowerNum = GetTowerNum(i, ic, ir);
 680         // high gain
 681         if (fNHighGain[TowerNum] > 0) {
 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
 688         if (fNLowGain[TowerNum] > 0) {
 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);
 709     }
 710   }
 711
 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
 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
 791   return kTRUE;
 792 }