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