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