EMCAL/AliEMCALCalibTimeDep.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
  16 /* $Id: AliEMCALCalibTimeDep.cxx $ */
  17
  18 //_________________________________________________________________________
  19 ///*-- Author:
  20 ///////////////////////////////////////////////////////////////////////////////
  21 //                                                                           //
  22 // class for EMCAL time-dep calibration
  23 // - supposed to run in preprocessor
  24 // we use input from the following sources:
  25 // AliEMCALCalibTempCoeff (APD temperature coefficients),
  26 // AliCaloCalibSignal (LED DA), AliEMCALSensorTempArray (ELMB DCS)
  27 // AliEMCALCalibReference: LED amplitude and temperature info at reference time
  28 //
  29 // output/result is in AliEMCALCalibTimeDepCorrection
  30 //                                                                           //
  31 ///////////////////////////////////////////////////////////////////////////////
  32
  33 #include <iostream>
  34 #include <TGraphSmooth.h>
  35 #include <TMath.h>
  36 #include "AliLog.h"
  37 #include "AliCDBEntry.h"
  38 #include "AliCDBManager.h"
  39 #include "AliEMCALSensorTempArray.h"
  40 #include "AliCaloCalibSignal.h"
  41 #include "AliEMCALCalibTempCoeff.h"
  42 #include "AliEMCALCalibReference.h"
  43 #include "AliEMCALCalibTimeDepCorrection.h"
  44 #include "AliEMCALCalibTimeDep.h"
  45
  46 /* first a bunch of constants.. */
  47 const double kSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
  48
  49 // some global variables for APD handling; values from Catania studies, best fit
  50 // TempCoeff = p0+p1*M (M=gain), where p0 and and p1 are functions of the dark current
  51 const double kTempCoeffP0Const = -0.903; //
  52 const double kTempCoeffP0Factor = -1.381e7; //
  53 const double kTempCoeffP1Const = -0.023; //
  54 const double kTempCoeffP1Factor = -4.966e5; //
  55
  56 const double kTempMaxDiffMedian = 2; // Temperature values should not be further away from median value within SM when considered in the average calc.
  57
  58 const double kErrorCode = -999; // to indicate that something went wrong
  59
  60 using namespace std;
  61
  62 ClassImp(AliEMCALCalibTimeDep)
  63
  64 //________________________________________________________________
  65 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
  66   fRun(0),
  67   fStartTime(0),
  68   fEndTime(0),
  69   fMinTemp(0),
  70   fMaxTemp(0),
  71   fMinTempVariation(0),
  72   fMaxTempVariation(0),
  73   fMinTempValid(15),
  74   fMaxTempValid(35),
  75   fMinTime(0),
  76   fMaxTime(0),
  77   fTemperatureResolution(0.1), // 0.1 deg C is default
  78   fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference
  79   fTimeBinsPerHour(2), // 2 30-min bins per hour is default
  80   fHighLowGainFactor(16), // factor ~16 between High gain and low gain
  81   fTempArray(NULL),
  82   fCalibSignal(NULL),
  83   fCalibTempCoeff(NULL),
  84   fCalibReference(NULL),
  85   fCalibTimeDepCorrection(NULL),
  86   fVerbosity(0)
  87 {
  88   // Constructor
  89 }
  90
  91 //________________________________________________________________
  92 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
  93   TObject(calibt),
  94   fRun(calibt.GetRunNumber()),
  95   fStartTime(calibt.GetStartTime()),
  96   fEndTime(calibt.GetEndTime()),
  97   fMinTemp(calibt.GetMinTemp()),
  98   fMaxTemp(calibt.GetMaxTemp()),
  99   fMinTempVariation(calibt.GetMinTempVariation()),
 100   fMaxTempVariation(calibt.GetMaxTempVariation()),
 101   fMinTempValid(calibt.GetMinTempValid()),
 102   fMaxTempValid(calibt.GetMaxTempValid()),
 103   fMinTime(calibt.GetMinTime()),
 104   fMaxTime(calibt.GetMaxTime()),
 105   fTemperatureResolution(calibt.GetTemperatureResolution()),
 106   fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()),
 107   fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
 108   fHighLowGainFactor(calibt.GetHighLowGainFactor()),
 109   fTempArray(calibt.GetTempArray()),
 110   fCalibSignal(calibt.GetCalibSignal()),
 111   fCalibTempCoeff(calibt.GetCalibTempCoeff()),
 112   fCalibReference(calibt.GetCalibReference()),
 113   fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
 114   fVerbosity(calibt.GetVerbosity())
 115 {
 116   // copy constructor
 117 }
 118
 119
 120 //________________________________________________________________
 121 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
 122 {
 123   // assignment operator; use copy ctor
 124   if (&calibt == this) return *this;
 125
 126   new (this) AliEMCALCalibTimeDep(calibt);
 127   return *this;
 128 }
 129
 130 //________________________________________________________________
 131 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
 132 {
 133   // Destructor
 134 }
 135
 136 //________________________________________________________________
 137 void  AliEMCALCalibTimeDep::Reset()
 138 {
 139   // clear variables to default
 140   fRun = 0;
 141   fStartTime = 0;
 142   fEndTime = 0;
 143   fMinTemp = 0;
 144   fMaxTemp = 0;
 145   fMinTempVariation = 0;
 146   fMaxTempVariation = 0;
 147   fMinTempValid = 15;
 148   fMaxTempValid = 35;
 149   fMinTime = 0;
 150   fMaxTime = 0;
 151   fTemperatureResolution = 0.1; // 0.1 deg C is default
 152   fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference
 153   fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
 154   fTempArray = NULL;
 155   fCalibSignal = NULL;
 156   fCalibTempCoeff = NULL;
 157   fCalibReference = NULL;
 158   fCalibTimeDepCorrection = NULL;
 159   fVerbosity = 0;
 160   return;
 161 }
 162
 163 //________________________________________________________________
 164 void  AliEMCALCalibTimeDep::PrintInfo() const
 165 {
 166   // print some info
 167   cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
 168   // basic variables, all 'publicly available' also
 169   cout << " VARIABLE DUMP: " << endl
 170        << " GetStartTime() " << GetStartTime() << endl
 171        << " GetEndTime() " << GetEndTime() << endl
 172        << " GetMinTime() " << GetMinTime() << endl
 173        << " GetMaxTime() " << GetMaxTime() << endl
 174        << " GetMinTemp() " << GetMinTemp() << endl
 175        << " GetMaxTemp() " << GetMaxTemp() << endl
 176        << " GetMinTempVariation() " << GetMinTempVariation() << endl
 177        << " GetMaxTempVariation() " << GetMaxTempVariation() << endl
 178        << " GetTemperatureResolution() " << GetTemperatureResolution() << endl;
 179   // run ranges
 180   cout << " RUN INFO: " << endl
 181        << " runnumber " << GetRunNumber() << endl
 182        << " length (in hours) " << GetLengthOfRunInHours() << endl
 183        << " length (in bins) " << GetLengthOfRunInBins() << endl
 184        << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
 185        << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
 186        << endl;
 187
 188   return;
 189 }
 190
 191 //________________________________________________________________
 192 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
 193 {
 194   return (fEndTime - fStartTime)*kSecToHour;
 195 }
 196
 197 //________________________________________________________________
 198 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
 199 {
 200   return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
 201 }
 202
 203 //________________________________________________________________
 204 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
 205 {
 206   return (fMaxTime - fMinTime)*kSecToHour;
 207 }
 208
 209 //________________________________________________________________
 210 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
 211 {
 212   return (fMaxTemp - fMinTemp);
 213 }
 214
 215 //________________________________________________________________
 216 void AliEMCALCalibTimeDep::Initialize(Int_t run,
 217                                       UInt_t startTime, UInt_t endTime)
 218 { // setup, and get temperature info
 219   Reset(); // start fresh
 220
 221   fRun = run;
 222   fStartTime = startTime;
 223   fEndTime = endTime;
 224
 225   // collect the needed information
 226   GetTemperatureInfo(); // temperature readings during the run
 227   ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
 228
 229   return;
 230 }
 231
 232 //________________________________________________________________
 233 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
 234 {// return estimate for this one SuperModule, if it had data
 235
 236   // first convert from seconds to hours..
 237   Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
 238
 239   int n = 0;
 240   Double_t valArr[8]={0}; // 8 sensors per SM
 241
 242   for (int i=0; i<fTempArray->NumSensors(); i++) {
 243
 244     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 245     int module = st->GetSector()*2 + st->GetSide();
 246     if ( module == imod ) { // right module
 247       // check if we had valid data for the time that is being asked for
 248       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
 249         AliSplineFit *f = st->GetFit();
 250         if (f) { // ok, looks like we have valid data/info
 251           // let's check what the expected value at the time appears to be
 252           Double_t val = f->Eval(timeHour);
 253           if ( fVerbosity > 0 ) {
 254             cout << " sensor i " << i << " val " << val << endl;
 255           }
 256           if (val>fMinTempValid && val<fMaxTempValid && n<8) {
 257             valArr[n] = val;
 258             n++;
 259           }
 260         }
 261       } // time
 262     }
 263
 264   } // loop over fTempArray
 265
 266   if (n>0) { // some valid data was found
 267     Double_t median = TMath::Median(n, valArr);
 268     Double_t average = 0;
 269     Int_t nval = 0;
 270     for (int is=0; is<n; is++) {
 271       if (TMath::Abs(valArr[is] - median) < kTempMaxDiffMedian) {
 272         average += valArr[is];
 273         nval++;
 274       }
 275     }
 276     if (nval >  0) {
 277       average /= nval;
 278       return average;
 279     }
 280     else { // this case should not happen, but kept for completeness (coverity etc)
 281       return median;
 282     }
 283   }
 284   else { // no good data
 285     return kErrorCode;
 286   }
 287
 288 }
 289
 290 //________________________________________________________________
 291 Int_t AliEMCALCalibTimeDep::CalcCorrection()
 292 { // OK, this is where the real action takes place - the heart of this class..
 293   /* The philosophy is as follows:
 294      0. Init corrections to 1.0 values, and see how many correction bins we need
 295      1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
 296      2. try to use temperature info + APD temperature coefficient info, to estimate correction.
 297      For now (from Dec 2009), we do not use LED info.
 298    */
 299
 300   // 0: Init
 301   // how many SuperModules do we have?
 302   Int_t nSM = fCalibReference->GetNSuperModule();
 303   // how many time-bins should we have for this run?
 304   Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
 305   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
 306
 307   // 1: get info on how much individual sensors might have changed during
 308   // the run (compare max-min for each sensor separately)
 309   if (fMaxTempVariation < fTemperatureResolution) {
 310     nBins = 1; // just one bin needed..
 311   }
 312   if (nBins == 1) {
 313     binSize = fEndTime - fStartTime;
 314   }
 315   if (fVerbosity > 0) {
 316     cout << " nBins " << nBins << " binSize " << binSize << endl;
 317   }
 318
 319   // set up a reasonable default (correction = 1.0)
 320   fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
 321   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
 322   fCalibTimeDepCorrection->SetStartTime(fStartTime);
 323   fCalibTimeDepCorrection->SetNTimeBins(nBins);
 324   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
 325
 326   // 2: try with Temperature correction
 327   Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
 328
 329   return nRemaining;
 330 }
 331
 332
 333 //________________________________________________________________
 334 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
 335 { // estimate the Temperature Coefficient, based on the dark current (IDark)
 336   // and the gain (M), based on Catania parameterizations
 337
 338   Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
 339   Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
 340
 341   Double_t dTC = dP0 + dP1*M;
 342   // from % numbers to regular ones..:
 343   dTC *= 0.01;
 344
 345   return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
 346 }
 347
 348 /* Next come the methods that do the work in picking up all the needed info..*/
 349 //________________________________________________________________
 350 void AliEMCALCalibTimeDep::GetTemperatureInfo()
 351 {
 352   // pick up Preprocessor output, based on fRun (most recent version)
 353   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
 354   if (entry) {
 355     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
 356   }
 357
 358   if (fTempArray) {
 359     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
 360                   fTempArray->NumSensors(),
 361                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
 362   }
 363   else {
 364     AliWarning( Form("AliEMCALSensorTempArray not found!") );
 365   }
 366
 367   return;
 368 }
 369
 370 //________________________________________________________________
 371 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
 372 {// assign max/min time and temperature values
 373
 374   fMinTemp = 999; // init to some large value (999 deg C)
 375   fMaxTemp = 0;
 376   fMinTempVariation = 999; // init to some large value (999 deg C)
 377   fMaxTempVariation = 0;
 378   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
 379   fMaxTime = 0;
 380
 381   Int_t n = 0; // number of valid readings
 382
 383   for (int i=0; i<fTempArray->NumSensors(); i++) {
 384
 385     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 386     if ( st->GetStartTime() == 0 ) { // no valid data
 387       continue;
 388     }
 389
 390     // check time ranges
 391     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
 392     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
 393
 394     // check temperature ranges
 395     AliSplineFit *f = st->GetFit();
 396
 397     if (f) { // ok, looks like we have valid data/info
 398       int np = f->GetKnots();
 399       Double_t *y0 = f->GetY0();
 400       // min and max values within the single sensor
 401       Double_t min = 999;
 402       Double_t max = 0;
 403       int nval = 0;
 404       for (int ip=0; ip<np; ip++) {
 405         if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
 406           if (min > y0[ip]) { min = y0[ip]; }
 407           if (max < y0[ip]) { max = y0[ip]; }
 408           nval++;
 409         }
 410       }
 411       if (nval>0) {
 412         if (fMinTemp > min) { fMinTemp = min; }
 413         if (fMaxTemp < max) { fMaxTemp = max; }
 414         Double_t variation = max - min;
 415         if (fMinTempVariation > variation) { fMinTempVariation = variation; }
 416         if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
 417
 418         n++;
 419       }
 420     }
 421   } // loop over fTempArray
 422
 423   if (n>0) { // some valid data was found
 424     return n;
 425   }
 426   else { // no good data
 427     return (Int_t) kErrorCode;
 428   }
 429
 430 }
 431
 432 //________________________________________________________________
 433 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
 434 {
 435   // pick up Preprocessor output, based on fRun (most recent version)
 436   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
 437   if (entry) {
 438     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
 439   }
 440
 441   if (fCalibSignal) {
 442     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
 443                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
 444                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
 445                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
 446                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
 447                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
 448   }
 449   else {
 450     AliWarning( Form("AliCaloCalibSignal not found!") );
 451   }
 452
 453   return;
 454 }
 455
 456 //________________________________________________________________
 457 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo()
 458 {
 459   // pick up Preprocessor output, based on fRun (most recent version)
 460   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
 461   // stored object should be a TTree; read the info
 462   if (entry) {
 463     fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
 464   }
 465
 466   if (fCalibTempCoeff) {
 467     AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
 468   }
 469   else {
 470     AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
 471   }
 472
 473   return;
 474 }
 475
 476 //________________________________________________________________
 477 void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
 478 {
 479   // pick up Preprocessor output, based on fRun (most recent version)
 480   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
 481   if (entry) {
 482     fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
 483   }
 484
 485   if (fCalibReference) {
 486     AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
 487   }
 488   else {
 489     AliWarning( Form("AliEMCALCalibReference not found!") );
 490   }
 491
 492   return;
 493 }
 494
 495 //________________________________________________________________
 496 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
 497 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
 498   // The correction factor we keep is c(T) = R(t0)/R(T)
 499   // T info from fCalibSignal, t0 info from fCalibReference
 500
 501   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
 502   // but one could upgrade this in the future
 503   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 504
 505   // sanity check; same SuperModule indices for corrections as for regular calibrations
 506   for (int i = 0; i < nSM; i++) {
 507     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
 508     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 509
 510     int iSMRef = dataCalibReference->GetSuperModuleNum();
 511     int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
 512     if (iSMRef != iSMCorr) {
 513       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
 514       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 515       return nRemaining;
 516     }
 517   }
 518
 519   int iSM = 0;
 520   int iCol = 0;
 521   int iRow = 0;
 522   int iStrip = 0;
 523   int iGain = 0;
 524
 525   // The fCalibSignal info is stored in TTrees
 526   // Note that the time-bins for the TTree's may not exactly match with our correction time bins
 527   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
 528   // fCalibSignal time info in seconds: Hour/kSecToHour
 529   // corrected for startTime difference: Hour/kSecToHour + timeDiff
 530   // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
 531
 532   // values for R(T), size of TArray = nBins
 533   // the [2] dimension below is for the low or high gain
 534   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 535   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 536   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 537   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 538
 539   // set up TArray's first
 540   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 541     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 542       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 543         for (iGain = 0; iGain < 2; iGain++) {
 544           // length of arrays
 545           ampT[iSM][iCol][iRow][iGain].Set(nBins);
 546           nT[iSM][iCol][iRow][iGain].Set(nBins);
 547           // content of arrys
 548           for (int j = 0; j < nBins; j++) {
 549             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
 550             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
 551           }
 552         }
 553       }
 554     }//iCol
 555     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 556       for (iGain = 0; iGain < 2; iGain++) {
 557         // length of arrays
 558         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
 559         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
 560         // content of arrys
 561         for (int j = 0; j < nBins; j++) {
 562           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 563           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 564         }
 565       }
 566     }//iStrip
 567   }
 568
 569   // OK, now loop over the TTrees and fill the arrays needed for R(T)
 570   TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 571   TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 572
 573   int iChannelNum = 0; // for regular towers
 574   int iRefNum = 0; // for LED
 575   double dHour = 0;
 576   double dAvgAmp = 0;
 577
 578   treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
 579   treeAvg->SetBranchAddress("fHour", &dHour);
 580   treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
 581
 582   int iBin = 0;
 583   // counters for how many values were seen per SuperModule
 584   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
 585   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
 586
 587   for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
 588     treeAvg->GetEntry(ient);
 589     // figure out where this info comes from
 590     fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
 591     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 592     // add value in the arrays
 593     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
 594     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
 595     nCount[iSM]++;
 596   }
 597
 598   treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
 599   treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
 600   treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
 601
 602   for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
 603     treeLEDRefAvg->GetEntry(ient);
 604     // figure out where this info comes from
 605     fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
 606     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 607     // add value in the arrays
 608     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
 609     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
 610     nCountLEDRef[iSM]++;
 611   }
 612
 613   // Normalize TArray values, and calculate average also
 614   Float_t norm = 0; // extra var, for readability
 615
 616   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 617     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
 618       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 619         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 620         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 621         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 622           for (iGain = 0; iGain < 2; iGain++) {
 623             // content of arrys
 624             for (int j = 0; j < nBins; j++) {
 625               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
 626                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
 627                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
 628               }
 629             }
 630           }
 631         }
 632       }//iCol
 633       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 634         for (iGain = 0; iGain < 2; iGain++) {
 635           for (int j = 0; j < nBins; j++) {
 636             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
 637               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
 638               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
 639             }
 640           }
 641         }
 642       }//iStrip
 643     }
 644   } // iSM
 645
 646
 647   // Calculate correction values, and store them
 648   // set kErrorCode values for those that could not be set
 649
 650   Float_t ratiot0 = 0;
 651   Float_t ratioT = 0;
 652   Float_t correction = 0; // c(T) = R(t0)/R(T)
 653   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
 654
 655   for (int i = 0; i < nSM; i++) {
 656     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
 657     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 658     iSM = dataCalibReference->GetSuperModuleNum();
 659
 660     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 661       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 662       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 663         refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
 664
 665       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 666
 667         // Calc. R(t0):
 668         AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
 669         iGain = refVal->GetHighLow(); // gain value used for reference calibration
 670         // valid amplitude values should be larger than 0
 671         if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
 672           ratiot0 =  refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
 673         }
 674         else {
 675           ratiot0 = kErrorCode;
 676         }
 677
 678         // Calc. R(T)
 679         for (int j = 0; j < nBins; j++) {
 680
 681           // calculate R(T) also; first try with individual tower:
 682           // same gain as for reference calibration is the default
 683           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
 684             // looks like valid data with the right gain combination
 685             ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
 686
 687             // if data appears to be saturated, and we are in high gain, then try with low gain instead
 688             int newGain = iGain;
 689             int newRefGain = refGain;
 690             if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
 691               newGain = 0;
 692             }
 693             if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
 694               newRefGain = 0;
 695             }
 696
 697             if (newGain!=iGain || newRefGain!=refGain) {
 698               // compensate for using different gain than in the reference calibration
 699               // we may need to have a custom H/L ratio value for each tower
 700               // later, but for now just use a common value, as the rest of the code does..
 701               ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
 702
 703               if (newGain<iGain) {
 704                 ratioT *= fHighLowGainFactor;
 705               }
 706               else if (newRefGain<refGain) {
 707                 ratioT /= fHighLowGainFactor;
 708               }
 709             }
 710           }
 711           else {
 712             ratioT = kErrorCode;
 713           }
 714
 715           // Calc. correction factor
 716           if (ratiot0>0 && ratioT>0) {
 717             correction = ratiot0/ratioT;
 718           }
 719           else {
 720             correction = kErrorCode;
 721             nRemaining++;
 722           }
 723
 724           // Store the value
 725           dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 726           /* Check that
 727           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
 728           also works OK */
 729         } // nBins
 730       }
 731     }
 732   }
 733
 734   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
 735   return nRemaining;
 736 }
 737
 738 //________________________________________________________________
 739 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
 740 { // use averages for each strip if no good values exist for some single tower
 741
 742   // go over fTimeDepCorrection info
 743   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 744
 745   // for calculating StripAverage info
 746   int nValidTower = 0;
 747   Float_t stripAverage = 0;
 748   Float_t val = 0;
 749
 750   int iSM = 0;
 751   int iCol = 0;
 752   int iRow = 0;
 753   int iStrip = 0;
 754   int firstCol = 0;
 755   int lastCol = 0;
 756
 757   for (int i = 0; i < nSM; i++) {
 758     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 759     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
 760
 761     for (int j = 0; j < nBins; j++) {
 762
 763       nValidTower = 0;
 764       stripAverage = 0;
 765
 766       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 767         firstCol = iStrip*2;
 768         if ((iSM%2)==1) { // C side
 769           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
 770         }
 771         lastCol = firstCol+1;
 772
 773         for (iCol = firstCol; iCol <= lastCol; iCol++) {
 774           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 775             val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
 776             if (val>0) { // valid value; error code is negative
 777               stripAverage += val;
 778               nValidTower++;
 779             }
 780           }
 781         }
 782
 783         // calc average over strip
 784         if (nValidTower>0) {
 785           stripAverage /= nValidTower;
 786           for (iCol = firstCol; iCol <= lastCol; iCol++) {
 787             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 788               val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
 789               if (val<0) { // invalid value; error code is negative
 790                 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
 791               }
 792             }
 793           }
 794         }
 795         else { // could not fill in unvalid entries
 796           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
 797         }
 798
 799       } // iStrip
 800     } // j, bins
 801   } // iSM
 802
 803   return nRemaining;
 804 }
 805
 806 //________________________________________________________________
 807 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize)
 808 { // OK, so we didn't have valid LED data that allowed us to do the correction only
 809   // with that info.
 810   // So, instead we'll rely on the temperature info and try to do the correction
 811   // based on that instead.
 812   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
 813   Int_t nRemaining = 0;
 814
 815   int iSM = 0;
 816   int iCol = 0;
 817   int iRow = 0;
 818
 819   Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
 820   memset(dTempCoeff, 0, sizeof(dTempCoeff));
 821   Double_t correction = 0;
 822   Double_t secondsPerBin = (Double_t) binSize;
 823
 824   for (int i = 0; i < nSM; i++) {
 825     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 826     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
 827
 828     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
 829     AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
 830
 831     // first get CalibTempCoeff info
 832     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 833       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 834
 835         dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
 836         if (fVerbosity > 1) {
 837           cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
 838                << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
 839         }
 840       }
 841     }
 842
 843     // figure out what the reference temperature is, from fCalibReference
 844     Double_t referenceTemperature = 0;
 845     int nVal = 0;
 846     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
 847       if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
 848         referenceTemperature += dataCalibReference->GetTemperature(iSensor);
 849         nVal++;
 850       }
 851     }
 852
 853     if (nVal>0) {
 854       referenceTemperature /= nVal; // valid values exist, we can look into corrections
 855
 856       for (int j = 0; j < nBins; j++) {
 857         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
 858         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
 859       // get the temperature at this time; use average over whole SM for now (TO BE CHECKED LATER - if we can do better with finer grained info)
 860         Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp);
 861
 862         Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
 863         if (fVerbosity > 0) {
 864           cout << " referenceTemperature " << referenceTemperature
 865                << " dSMTemperature " << dSMTemperature
 866                << " temperatureDiff " << temperatureDiff
 867                << endl;
 868         }
 869         // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down
 870         // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
 871         // dTempCoeff is a (unsigned) factor describing how many % the gain
 872         // changes with a degree change.
 873         // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
 874         // The correction we want to keep is what we should multiply our ADC value with as a function
 875         // of time, i.e. the inverse of the gain change..
 876         if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
 877              && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
 878           // significant enough difference that we need to consider it, and also not unreasonably large
 879
 880           // loop over all towers; effect of temperature change will depend on gain for this tower
 881           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 882             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 883
 884               // the correction should be inverse of modification in gain: (see discussion above)
 885               // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
 886               // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
 887               correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]);
 888               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 889               if (fVerbosity > 1) {
 890                 cout << " iSM " << iSM
 891                      << " iCol  " << iCol
 892                      << " iRow  " << iRow
 893                      << " j  " << j
 894                      << " correction  " << correction
 895                      << endl;
 896               }
 897             }
 898           }
 899
 900         } // if noteworthy temperature change
 901         else { // just set correction values to 1.0
 902           correction = 1;
 903           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 904             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 905               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 906             }
 907           }
 908         } // else
 909       } // j, Bins
 910
 911     } // if reference temperature exist
 912     else { // could not do the needed check.. signal that in the return code
 913       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 914     }
 915   } // iSM
 916
 917   return nRemaining;
 918 }
 919