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     //cout << " n " << n << " nval " << nval << " median " << median << endl;
 277     if (nval >  0) {
 278       average /= nval;
 279       //cout << " average " << average << endl;
 280       return average;
 281     }
 282     else { // this case should not happen, but kept for completeness (coverity etc)
 283       return median;
 284     }
 285   }
 286   else { // no good data
 287     return kErrorCode;
 288   }
 289
 290 }
 291
 292 //________________________________________________________________
 293 Int_t AliEMCALCalibTimeDep::CalcCorrection()
 294 { // OK, this is where the real action takes place - the heart of this class..
 295   /* The philosophy is as follows:
 296      0. Init corrections to 1.0 values, and see how many correction bins we need
 297      1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
 298      2. try to use temperature info + APD temperature coefficient info, to estimate correction.
 299      For now (from Dec 2009), we do not use LED info.
 300    */
 301
 302   // 0: Init
 303   // how many SuperModules do we have?
 304   Int_t nSM = fCalibReference->GetNSuperModule();
 305   // how many time-bins should we have for this run?
 306   Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
 307   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
 308
 309   // 1: get info on how much individual sensors might have changed during
 310   // the run (compare max-min for each sensor separately)
 311   if (fMaxTempVariation < fTemperatureResolution) {
 312     nBins = 1; // just one bin needed..
 313   }
 314   if (nBins == 1) {
 315     binSize = fEndTime - fStartTime;
 316   }
 317   if (fVerbosity > 0) {
 318     cout << " nBins " << nBins << " binSize " << binSize << endl;
 319   }
 320
 321   // set up a reasonable default (correction = 1.0)
 322   fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
 323   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
 324   fCalibTimeDepCorrection->SetStartTime(fStartTime);
 325   fCalibTimeDepCorrection->SetNTimeBins(nBins);
 326   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
 327
 328   // 2: try with Temperature correction
 329   Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
 330
 331   return nRemaining;
 332 }
 333
 334
 335 //________________________________________________________________
 336 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
 337 { // estimate the Temperature Coefficient, based on the dark current (IDark)
 338   // and the gain (M), based on Catania parameterizations
 339
 340   Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
 341   Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
 342
 343   Double_t dTC = dP0 + dP1*M;
 344   // from % numbers to regular ones..:
 345   dTC *= 0.01;
 346
 347   return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
 348 }
 349
 350 /* Next come the methods that do the work in picking up all the needed info..*/
 351 //________________________________________________________________
 352 void AliEMCALCalibTimeDep::GetTemperatureInfo()
 353 {
 354   // pick up Preprocessor output, based on fRun (most recent version)
 355   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
 356   if (entry) {
 357     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
 358   }
 359
 360   if (fTempArray) {
 361     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
 362                   fTempArray->NumSensors(),
 363                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
 364   }
 365   else {
 366     AliWarning( Form("AliEMCALSensorTempArray not found!") );
 367   }
 368
 369   return;
 370 }
 371
 372 //________________________________________________________________
 373 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
 374 {// assign max/min time and temperature values
 375
 376   fMinTemp = 999; // init to some large value (999 deg C)
 377   fMaxTemp = 0;
 378   fMinTempVariation = 999; // init to some large value (999 deg C)
 379   fMaxTempVariation = 0;
 380   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
 381   fMaxTime = 0;
 382
 383   Int_t n = 0; // number of valid readings
 384
 385   for (int i=0; i<fTempArray->NumSensors(); i++) {
 386
 387     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 388     if ( st->GetStartTime() == 0 ) { // no valid data
 389       continue;
 390     }
 391
 392     // check time ranges
 393     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
 394     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
 395
 396     // check temperature ranges
 397     AliSplineFit *f = st->GetFit();
 398
 399     if (f) { // ok, looks like we have valid data/info
 400       int np = f->GetKnots();
 401       Double_t *y0 = f->GetY0();
 402       // min and max values within the single sensor
 403       Double_t min = 999;
 404       Double_t max = 0;
 405       int nval = 0;
 406       for (int ip=0; ip<np; ip++) {
 407         if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
 408           if (min > y0[ip]) { min = y0[ip]; }
 409           if (max < y0[ip]) { max = y0[ip]; }
 410           nval++;
 411         }
 412       }
 413       if (nval>0) {
 414         if (fMinTemp > min) { fMinTemp = min; }
 415         if (fMaxTemp < max) { fMaxTemp = max; }
 416         Double_t variation = max - min;
 417         if (fMinTempVariation > variation) { fMinTempVariation = variation; }
 418         if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
 419
 420         n++;
 421       }
 422     }
 423   } // loop over fTempArray
 424
 425   if (n>0) { // some valid data was found
 426     return n;
 427   }
 428   else { // no good data
 429     return (Int_t) kErrorCode;
 430   }
 431
 432 }
 433
 434 //________________________________________________________________
 435 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
 436 {
 437   // pick up Preprocessor output, based on fRun (most recent version)
 438   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
 439   if (entry) {
 440     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
 441   }
 442
 443   if (fCalibSignal) {
 444     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
 445                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
 446                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
 447                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
 448                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
 449                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
 450   }
 451   else {
 452     AliWarning( Form("AliCaloCalibSignal not found!") );
 453   }
 454
 455   return;
 456 }
 457
 458 //________________________________________________________________
 459 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo()
 460 {
 461   // pick up Preprocessor output, based on fRun (most recent version)
 462   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
 463   // stored object should be a TTree; read the info
 464   if (entry) {
 465     fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
 466   }
 467
 468   if (fCalibTempCoeff) {
 469     AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
 470   }
 471   else {
 472     AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
 473   }
 474
 475   return;
 476 }
 477
 478 //________________________________________________________________
 479 void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
 480 {
 481   // pick up Preprocessor output, based on fRun (most recent version)
 482   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
 483   if (entry) {
 484     fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
 485   }
 486
 487   if (fCalibReference) {
 488     AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
 489   }
 490   else {
 491     AliWarning( Form("AliEMCALCalibReference not found!") );
 492   }
 493
 494   return;
 495 }
 496
 497 //________________________________________________________________
 498 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
 499 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
 500   // The correction factor we keep is c(T) = R(t0)/R(T)
 501   // T info from fCalibSignal, t0 info from fCalibReference
 502
 503   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
 504   // but one could upgrade this in the future
 505   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 506
 507   // sanity check; same SuperModule indices for corrections as for regular calibrations
 508   for (int i = 0; i < nSM; i++) {
 509     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
 510     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 511
 512     int iSMRef = dataCalibReference->GetSuperModuleNum();
 513     int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
 514     if (iSMRef != iSMCorr) {
 515       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
 516       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 517       return nRemaining;
 518     }
 519   }
 520
 521   int iSM = 0;
 522   int iCol = 0;
 523   int iRow = 0;
 524   int iStrip = 0;
 525   int iGain = 0;
 526
 527   // The fCalibSignal info is stored in TTrees
 528   // Note that the time-bins for the TTree's may not exactly match with our correction time bins
 529   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
 530   // fCalibSignal time info in seconds: Hour/kSecToHour
 531   // corrected for startTime difference: Hour/kSecToHour + timeDiff
 532   // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
 533
 534   // values for R(T), size of TArray = nBins
 535   // the [2] dimension below is for the low or high gain
 536   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 537   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 538   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 539   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 540
 541   // set up TArray's first
 542   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 543     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 544       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 545         for (iGain = 0; iGain < 2; iGain++) {
 546           // length of arrays
 547           ampT[iSM][iCol][iRow][iGain].Set(nBins);
 548           nT[iSM][iCol][iRow][iGain].Set(nBins);
 549           // content of arrys
 550           for (int j = 0; j < nBins; j++) {
 551             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
 552             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
 553           }
 554         }
 555       }
 556     }//iCol
 557     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 558       for (iGain = 0; iGain < 2; iGain++) {
 559         // length of arrays
 560         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
 561         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
 562         // content of arrys
 563         for (int j = 0; j < nBins; j++) {
 564           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 565           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 566         }
 567       }
 568     }//iStrip
 569   }
 570
 571   // OK, now loop over the TTrees and fill the arrays needed for R(T)
 572   TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 573   TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 574
 575   int iChannelNum = 0; // for regular towers
 576   int iRefNum = 0; // for LED
 577   double dHour = 0;
 578   double dAvgAmp = 0;
 579
 580   treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
 581   treeAvg->SetBranchAddress("fHour", &dHour);
 582   treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
 583
 584   int iBin = 0;
 585   // counters for how many values were seen per SuperModule
 586   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
 587   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
 588
 589   for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
 590     treeAvg->GetEntry(ient);
 591     // figure out where this info comes from
 592     fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
 593     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 594     // add value in the arrays
 595     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
 596     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
 597     nCount[iSM]++;
 598   }
 599
 600   treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
 601   treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
 602   treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
 603
 604   for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
 605     treeLEDRefAvg->GetEntry(ient);
 606     // figure out where this info comes from
 607     fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
 608     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 609     // add value in the arrays
 610     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
 611     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
 612     nCountLEDRef[iSM]++;
 613   }
 614
 615   // Normalize TArray values, and calculate average also
 616   Float_t norm = 0; // extra var, for readability
 617
 618   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 619     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
 620       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 621         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 622         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 623         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 624           for (iGain = 0; iGain < 2; iGain++) {
 625             // content of arrys
 626             for (int j = 0; j < nBins; j++) {
 627               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
 628                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
 629                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
 630               }
 631             }
 632           }
 633         }
 634       }//iCol
 635       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 636         for (iGain = 0; iGain < 2; iGain++) {
 637           for (int j = 0; j < nBins; j++) {
 638             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
 639               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
 640               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
 641             }
 642           }
 643         }
 644       }//iStrip
 645     }
 646   } // iSM
 647
 648
 649   // Calculate correction values, and store them
 650   // set kErrorCode values for those that could not be set
 651
 652   Float_t ratiot0 = 0;
 653   Float_t ratioT = 0;
 654   Float_t correction = 0; // c(T) = R(t0)/R(T)
 655   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
 656
 657   for (int i = 0; i < nSM; i++) {
 658     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
 659     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 660     iSM = dataCalibReference->GetSuperModuleNum();
 661
 662     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 663       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 664       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 665         refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
 666
 667       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 668
 669         // Calc. R(t0):
 670         AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
 671         iGain = refVal->GetHighLow(); // gain value used for reference calibration
 672         // valid amplitude values should be larger than 0
 673         if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
 674           ratiot0 =  refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
 675         }
 676         else {
 677           ratiot0 = kErrorCode;
 678         }
 679
 680         // Calc. R(T)
 681         for (int j = 0; j < nBins; j++) {
 682
 683           // calculate R(T) also; first try with individual tower:
 684           // same gain as for reference calibration is the default
 685           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
 686             // looks like valid data with the right gain combination
 687             ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
 688
 689             // if data appears to be saturated, and we are in high gain, then try with low gain instead
 690             int newGain = iGain;
 691             int newRefGain = refGain;
 692             if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
 693               newGain = 0;
 694             }
 695             if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
 696               newRefGain = 0;
 697             }
 698
 699             if (newGain!=iGain || newRefGain!=refGain) {
 700               // compensate for using different gain than in the reference calibration
 701               // we may need to have a custom H/L ratio value for each tower
 702               // later, but for now just use a common value, as the rest of the code does..
 703               ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
 704
 705               if (newGain<iGain) {
 706                 ratioT *= fHighLowGainFactor;
 707               }
 708               else if (newRefGain<refGain) {
 709                 ratioT /= fHighLowGainFactor;
 710               }
 711             }
 712           }
 713           else {
 714             ratioT = kErrorCode;
 715           }
 716
 717           // Calc. correction factor
 718           if (ratiot0>0 && ratioT>0) {
 719             correction = ratiot0/ratioT;
 720           }
 721           else {
 722             correction = kErrorCode;
 723             nRemaining++;
 724           }
 725
 726           // Store the value
 727           dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 728           /* Check that
 729           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
 730           also works OK */
 731         } // nBins
 732       }
 733     }
 734   }
 735
 736   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
 737   return nRemaining;
 738 }
 739
 740 //________________________________________________________________
 741 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
 742 { // use averages for each strip if no good values exist for some single tower
 743
 744   // go over fTimeDepCorrection info
 745   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 746
 747   // for calculating StripAverage info
 748   int nValidTower = 0;
 749   Float_t stripAverage = 0;
 750   Float_t val = 0;
 751
 752   int iSM = 0;
 753   int iCol = 0;
 754   int iRow = 0;
 755   int iStrip = 0;
 756   int firstCol = 0;
 757   int lastCol = 0;
 758
 759   for (int i = 0; i < nSM; i++) {
 760     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 761     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
 762
 763     for (int j = 0; j < nBins; j++) {
 764
 765       nValidTower = 0;
 766       stripAverage = 0;
 767
 768       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 769         firstCol = iStrip*2;
 770         if ((iSM%2)==1) { // C side
 771           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
 772         }
 773         lastCol = firstCol+1;
 774
 775         for (iCol = firstCol; iCol <= lastCol; iCol++) {
 776           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 777             val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
 778             if (val>0) { // valid value; error code is negative
 779               stripAverage += val;
 780               nValidTower++;
 781             }
 782           }
 783         }
 784
 785         // calc average over strip
 786         if (nValidTower>0) {
 787           stripAverage /= nValidTower;
 788           for (iCol = firstCol; iCol <= lastCol; iCol++) {
 789             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 790               val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
 791               if (val<0) { // invalid value; error code is negative
 792                 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
 793               }
 794             }
 795           }
 796         }
 797         else { // could not fill in unvalid entries
 798           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
 799         }
 800
 801       } // iStrip
 802     } // j, bins
 803   } // iSM
 804
 805   return nRemaining;
 806 }
 807
 808 //________________________________________________________________
 809 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize)
 810 { // OK, so we didn't have valid LED data that allowed us to do the correction only
 811   // with that info.
 812   // So, instead we'll rely on the temperature info and try to do the correction
 813   // based on that instead.
 814   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
 815   Int_t nRemaining = 0;
 816
 817   int iSM = 0;
 818   int iCol = 0;
 819   int iRow = 0;
 820
 821   Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
 822   memset(dTempCoeff, 0, sizeof(dTempCoeff));
 823   Double_t correction = 0;
 824   Double_t secondsPerBin = (Double_t) binSize;
 825
 826   for (int i = 0; i < nSM; i++) {
 827     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
 828     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
 829
 830     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
 831     AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
 832
 833     // first get CalibTempCoeff info
 834     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 835       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 836
 837         dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
 838         if (fVerbosity > 1) {
 839           cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
 840                << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
 841         }
 842       }
 843     }
 844
 845     // figure out what the reference temperature is, from fCalibReference
 846     Double_t referenceTemperature = 0;
 847     int nVal = 0;
 848     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
 849       if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
 850         referenceTemperature += dataCalibReference->GetTemperature(iSensor);
 851         nVal++;
 852       }
 853     }
 854
 855     if (nVal>0) {
 856       referenceTemperature /= nVal; // valid values exist, we can look into corrections
 857
 858       Double_t dSMTemperature = 0;
 859       for (int j = 0; j < nBins; j++) {
 860         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
 861         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
 862       // 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)
 863         Double_t oldSMTemperature = dSMTemperature;
 864         dSMTemperature = GetTemperatureSM(iSM, timeStamp);
 865         if (j>0 && (dSMTemperature==kErrorCode)) {
 866           // if we have previous values, and retrieval of values failed - use that instead (hopefully good)
 867           dSMTemperature = oldSMTemperature;
 868         }
 869
 870         Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
 871         if (fVerbosity > 0) {
 872           cout << " referenceTemperature " << referenceTemperature
 873                << " dSMTemperature " << dSMTemperature
 874                << " temperatureDiff " << temperatureDiff
 875                << endl;
 876         }
 877         // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down
 878         // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
 879         // dTempCoeff is a (unsigned) factor describing how many % the gain
 880         // changes with a degree change.
 881         // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
 882         // The correction we want to keep is what we should multiply our ADC value with as a function
 883         // of time, i.e. the inverse of the gain change..
 884         if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
 885              && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
 886           // significant enough difference that we need to consider it, and also not unreasonably large
 887
 888           // loop over all towers; effect of temperature change will depend on gain for this tower
 889           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 890             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 891
 892               // the correction should be inverse of modification in gain: (see discussion above)
 893               // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
 894               // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
 895               correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]);
 896               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 897               if (fVerbosity > 1) {
 898                 cout << " iSM " << iSM
 899                      << " iCol  " << iCol
 900                      << " iRow  " << iRow
 901                      << " j  " << j
 902                      << " correction  " << correction
 903                      << endl;
 904               }
 905             }
 906           }
 907
 908         } // if noteworthy temperature change
 909         else { // just set correction values to 1.0
 910           correction = 1;
 911           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 912             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 913               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
 914             }
 915           }
 916         } // else
 917       } // j, Bins
 918
 919     } // if reference temperature exist
 920     else { // could not do the needed check.. signal that in the return code
 921       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 922     }
 923   } // iSM
 924
 925   return nRemaining;
 926 }
 927