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