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 //                                                                           //
  24 ///////////////////////////////////////////////////////////////////////////////
  25
  26 #include <iostream>
  27 #include <TGraphSmooth.h>
  28 #include "AliLog.h"
  29 #include "AliCDBEntry.h"
  30 #include "AliCDBManager.h"
  31 #include "AliEMCALSensorTempArray.h"
  32 #include "AliCaloCalibSignal.h"
  33 #include "AliEMCALBiasAPD.h"
  34 #include "AliEMCALCalibMapAPD.h"
  35 #include "AliEMCALCalibAbs.h"
  36 #include "AliEMCALCalibTimeDepCorrection.h"
  37 #include "AliEMCALCalibTimeDep.h"
  38
  39 /* first a bunch of constants.. */
  40 const double fkSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
  41
  42 // some global variables for APD handling; values from Catania studies, best fit
  43 // TempCoeff = p0+p1*M (M=gain), where p0 and and p1 are functions of the dark current
  44 const double fkTempCoeffP0Const = -0.903; //
  45 const double fkTempCoeffP0Factor = -1.381e7; //
  46 const double fkTempCoeffP1Const = -0.023; //
  47 const double fkTempCoeffP1Factor = -4.966e5; //
  48
  49 const double fkErrorCode = -999; // to indicate that something went wrong
  50
  51 using namespace std;
  52
  53 ClassImp(AliEMCALCalibTimeDep)
  54
  55 //________________________________________________________________
  56 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
  57   fRun(0),
  58   fStartTime(0),
  59   fEndTime(0),
  60   fMinTemp(0),
  61   fMaxTemp(0),
  62   fMinTime(0),
  63   fMaxTime(0),
  64   fTemperatureResolution(0.1), // 0.1 deg C is default
  65   fTimeBinsPerHour(2), // 2 30-min bins per hour is default
  66   fTempArray(NULL),
  67   fCalibSignal(NULL),
  68   fBiasAPD(NULL),
  69   fCalibMapAPD(NULL),
  70   fCalibAbs(NULL),
  71   fCalibTimeDepCorrection(NULL)
  72 {
  73   // Constructor
  74 }
  75
  76 //________________________________________________________________
  77 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
  78   TObject(calibt),
  79   fRun(calibt.GetRunNumber()),
  80   fStartTime(calibt.GetStartTime()),
  81   fEndTime(calibt.GetEndTime()),
  82   fMinTemp(calibt.GetMinTemp()),
  83   fMaxTemp(calibt.GetMaxTemp()),
  84   fMinTime(calibt.GetMinTime()),
  85   fMaxTime(calibt.GetMaxTime()),
  86   fTemperatureResolution(calibt.GetTemperatureResolution()),
  87   fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
  88   fTempArray(calibt.GetTempArray()),
  89   fCalibSignal(calibt.GetCalibSignal()),
  90   fBiasAPD(calibt.GetBiasAPD()),
  91   fCalibMapAPD(calibt.GetCalibMapAPD()),
  92   fCalibAbs(calibt.GetCalibAbs()),
  93   fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection())
  94 {
  95   // copy constructor
  96 }
  97
  98
  99 //________________________________________________________________
 100 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
 101 {
 102   // assignment operator; use copy ctor
 103   if (&calibt == this) return *this;
 104
 105   new (this) AliEMCALCalibTimeDep(calibt);
 106   return *this;
 107 }
 108
 109 //________________________________________________________________
 110 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
 111 {
 112   // Destructor
 113 }
 114
 115 //________________________________________________________________
 116 void  AliEMCALCalibTimeDep::Reset()
 117 {
 118   // clear variables to default
 119   fRun = 0;
 120   fStartTime = 0;
 121   fEndTime = 0;
 122   fMinTemp = 0;
 123   fMaxTemp = 0;
 124   fMinTime = 0;
 125   fMaxTime = 0;
 126   fTemperatureResolution = 0.1; // 0.1 deg C is default
 127   fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
 128   fTempArray = NULL;
 129   fCalibSignal = NULL;
 130   fBiasAPD = NULL;
 131   fCalibMapAPD = NULL;
 132   fCalibAbs = NULL;
 133   fCalibTimeDepCorrection = NULL;
 134   return;
 135 }
 136
 137 //________________________________________________________________
 138 void  AliEMCALCalibTimeDep::PrintInfo() const
 139 {
 140   // print some info
 141   cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
 142   // basic variables, all 'publicly available' also
 143   cout << " VARIABLE DUMP: " << endl
 144        << " GetStartTime() " << GetStartTime() << endl
 145        << " GetEndTime() " << GetEndTime() << endl
 146        << " GetMinTemp() " << GetMinTemp() << endl
 147        << " GetMaxTemp() " << GetMaxTemp() << endl;
 148   // run ranges
 149   cout << " RUN INFO: " << endl
 150        << " length (in hours) " << GetLengthOfRunInHours() << endl
 151        << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
 152        << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
 153        << endl;
 154
 155   return;
 156 }
 157
 158 //________________________________________________________________
 159 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
 160 {
 161   return (fEndTime - fStartTime)*fkSecToHour;
 162 }
 163
 164 //________________________________________________________________
 165 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
 166 {
 167   return (fEndTime - fStartTime)*fkSecToHour*fTimeBinsPerHour;
 168 }
 169
 170 //________________________________________________________________
 171 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
 172 {
 173   return (fMaxTime - fMinTime)*fkSecToHour;
 174 }
 175
 176 //________________________________________________________________
 177 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
 178 {
 179   return (fMaxTemp - fMinTemp);
 180 }
 181
 182 //________________________________________________________________
 183 void AliEMCALCalibTimeDep::Initialize(Int_t run,
 184                                       UInt_t startTime, UInt_t endTime)
 185 {
 186   Reset(); // start fresh
 187
 188   fRun = run;
 189   fStartTime = startTime;
 190   fEndTime = endTime;
 191
 192   // collect the needed information
 193   GetTemperatureInfo(); // temperature readings during the run
 194   ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
 195
 196   return;
 197 }
 198
 199 //________________________________________________________________
 200 Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const
 201 {// return estimate for all SuperModules and sensors, that had data
 202
 203   // first convert from seconds to hours..
 204   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
 205
 206   Double_t average = 0;
 207   int n = 0;
 208
 209   for (int i=0; i<fTempArray->NumSensors(); i++) {
 210
 211     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 212
 213     // check if we had valid data for the time that is being asked for
 214     if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
 215       AliSplineFit *f = st->GetFit();
 216       if (f) { // ok, looks like we have valid data/info
 217         // let's check what the expected value at the time appears to be
 218         Double_t val = f->Eval(timeHour);
 219         average += val;
 220         n++;
 221       }
 222     } // time
 223   } // loop over fTempArray
 224
 225   if (n>0) { // some valid data was found
 226     average /= n;
 227     return average;
 228   }
 229   else { // no good data
 230     return fkErrorCode;
 231   }
 232
 233 }
 234
 235 //________________________________________________________________
 236 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
 237 {// return estimate for this one SuperModule, if it had data
 238
 239   // first convert from seconds to hours..
 240   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
 241
 242   Double_t average = 0;
 243   int n = 0;
 244
 245   for (int i=0; i<fTempArray->NumSensors(); i++) {
 246
 247     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 248     int module = st->GetSector()*2 + st->GetSide();
 249     if ( module == imod ) { // right module
 250       // check if we had valid data for the time that is being asked for
 251       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
 252         AliSplineFit *f = st->GetFit();
 253         if (f) { // ok, looks like we have valid data/info
 254           // let's check what the expected value at the time appears to be
 255           Double_t val = f->Eval(timeHour);
 256           cout << " i " << i << " val " << val << endl;
 257           average += val;
 258           n++;
 259         }
 260       } // time
 261     }
 262
 263   } // loop over fTempArray
 264
 265   if (n>0) { // some valid data was found
 266     average /= n;
 267     return average;
 268   }
 269   else { // no good data
 270     return fkErrorCode;
 271   }
 272
 273 }
 274
 275 //________________________________________________________________
 276 Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_t timeStamp) const
 277 {// return estimate for this one SuperModule and sensor, if it had data
 278
 279   // first convert from seconds to hours..
 280   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
 281
 282   for (int i=0; i<fTempArray->NumSensors(); i++) {
 283
 284     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 285     int module = st->GetSector()*2 + st->GetSide();
 286     if ( module == imod && st->GetNum()==isens ) { // right module, and sensor
 287       // check if we had valid data for the time that is being asked for
 288       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
 289         AliSplineFit *f = st->GetFit();
 290         if (f) { // ok, looks like we have valid data/info
 291           // let's check what the expected value at the time appears to be
 292           Double_t val = f->Eval(timeHour);
 293
 294           return val; // no point to move further in for loop, we have found the sensor we were looking for
 295         }
 296       } // time
 297     }
 298
 299   } // loop over fTempArray
 300
 301   // if we made it all here, it means that we didn't find the sensor we were looking for
 302   // i.e. no good data
 303   return fkErrorCode;
 304
 305 }
 306
 307 //________________________________________________________________
 308 Int_t AliEMCALCalibTimeDep::CalcCorrection()
 309 { // OK, this is where the real action takes place - the heart of this class..
 310   /* The philosophy is as follows:
 311      0. Init corrections to 1.0 values
 312      1: if we have LED info for the tower, use it
 313      2. if not 1, we rely on LED info averaged over strip
 314      3. if not 2 either, we try to use temperature info + APD bias and calibration info
 315    */
 316
 317   // 0: Init
 318   // how many SuperModules do we have?
 319   Int_t nSM = fCalibAbs->GetNSuperModule();
 320   // how many time-bins should we have for this run?
 321   Int_t nBins = (Int_t) GetLengthOfRunInBins(); // round-down (from double to int)
 322   // set up a reasonable default (correction = 1.0)
 323   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
 324
 325   // 1+2: try with LED corrections
 326   Int_t nRemaining = CalcLEDCorrection(nSM, nBins);
 327
 328   // 3: try with Temperature, if needed
 329   if (nRemaining>0) {
 330     nRemaining = CalcTemperatureCorrection(nSM, nBins);
 331   }
 332
 333   return nRemaining;
 334 }
 335
 336
 337 //________________________________________________________________
 338 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
 339 { // estimate the Temperature Coefficient, based on the dark current (IDark)
 340   // and the gain (M), based on Catania parameterizations
 341
 342   Double_t P0 = fkTempCoeffP0Const + fkTempCoeffP0Factor * IDark;
 343   Double_t P1 = fkTempCoeffP1Const + fkTempCoeffP1Factor * IDark;
 344
 345   Double_t TC = P0 + P1*M;
 346
 347   return TC;
 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   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
 379   fMaxTime = 0;
 380
 381   Int_t n = 0; // number of valid readings
 382
 383   for (int i=0; i<fTempArray->NumSensors(); i++) {
 384
 385     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 386
 387     // check time ranges
 388     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
 389     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
 390
 391     // check temperature ranges
 392     TGraph *g = st->GetGraph();
 393     if (g) { // ok, looks like we have valid data/info
 394       // let's check what the expected value at the time appears to be
 395       if (fMinTemp > g->GetMinimum()) { fMinTemp = g->GetMinimum(); }
 396       if (fMaxTemp < g->GetMaximum()) { fMaxTemp = g->GetMaximum(); }
 397       n++;
 398     }
 399   } // loop over fTempArray
 400
 401   if (n>0) { // some valid data was found
 402     return n;
 403   }
 404   else { // no good data
 405     return (Int_t) fkErrorCode;
 406   }
 407
 408 }
 409
 410 //________________________________________________________________
 411 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
 412 {
 413   // pick up Preprocessor output, based on fRun (most recent version)
 414   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
 415   if (entry) {
 416     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
 417   }
 418
 419   if (fCalibSignal) {
 420     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %d AvgEntries LEDRefEntries %d LEDRefAvgEntries %d",
 421                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
 422                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
 423                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
 424                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
 425                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
 426   }
 427   else {
 428     AliWarning( Form("AliCaloCalibSignal not found!") );
 429   }
 430
 431   return;
 432 }
 433
 434 //________________________________________________________________
 435 void AliEMCALCalibTimeDep::GetBiasAPDInfo()
 436 {
 437   // pick up Preprocessor output, based on fRun (most recent version)
 438   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/BiasAPD", fRun);
 439   if (entry) {
 440     fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
 441   }
 442
 443   if (fBiasAPD) {
 444     AliInfo( Form("BiasAPD: NSuperModule %d ", fBiasAPD->GetNSuperModule() ) );
 445   }
 446   else {
 447     AliWarning( Form("AliEMCALBiasAPD not found!") );
 448   }
 449
 450   return;
 451 }
 452
 453 //________________________________________________________________
 454 void AliEMCALCalibTimeDep::GetCalibMapAPDInfo()
 455 {
 456   // pick up Preprocessor output, based on fRun (most recent version)
 457   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
 458   if (entry) {
 459     fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
 460   }
 461
 462   if (fCalibMapAPD) {
 463     AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
 464   }
 465   else {
 466     AliWarning( Form("AliEMCALCalibMapAPD not found!") );
 467   }
 468
 469   return;
 470 }
 471
 472 //________________________________________________________________
 473 void AliEMCALCalibTimeDep::GetCalibAbsInfo()
 474 {
 475   // pick up Preprocessor output, based on fRun (most recent version)
 476   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
 477   if (entry) {
 478     fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject();
 479   }
 480
 481   if (fCalibAbs) {
 482     AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
 483   }
 484   else {
 485     AliWarning( Form("AliEMCALCalibAbs not found!") );
 486   }
 487
 488   return;
 489 }
 490
 491 //________________________________________________________________
 492 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
 493 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
 494   // The correction factor we keep is c(T) = R(t0)/R(T)
 495   // T info from fCalibSignal, t0 info from fCalibAbs
 496
 497   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
 498   // but one could upgrade this in the future
 499   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 500
 501   // sanity check; same SuperModule indices for corrections as for regular calibrations
 502   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 503   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 504
 505   for (int i = 0; i < nSM; i++) {
 506     int iSMAbs = CalibAbsData[i].fSuperModuleNum;
 507     int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 508     if (iSMAbs != iSMCorr) {
 509       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
 510       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 511       return nRemaining;
 512     }
 513   }
 514
 515   int iSM = 0;
 516   int iCol = 0;
 517   int iRow = 0;
 518   int iStrip = 0;
 519   int iGain = 0;
 520
 521   // The fCalibSignal info is stored in TTrees
 522   // Note that the time-bins for the TTree's may not exactly match with our correction time bins
 523   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
 524   // fCalibSignal time info in seconds: Hour/fkSecToHour
 525   // corrected for startTime difference: Hour/fkSecToHour + timeDiff
 526   // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
 527
 528   // values for R(T), size of TArray = nBins
 529   // the [2] dimension below is for the low or high gain
 530   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 531   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 532   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 533   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 534
 535   // set up TArray's first
 536   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 537     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 538       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 539         for (iGain = 0; iGain < 2; iGain++) {
 540           // length of arrays
 541           ampT[iSM][iCol][iRow][iGain].Set(nBins);
 542           nT[iSM][iCol][iRow][iGain].Set(nBins);
 543           // content of arrys
 544           for (int j = 0; j < nBins; j++) {
 545             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
 546             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
 547           }
 548         }
 549       }
 550     }//iCol
 551     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 552       for (iGain = 0; iGain < 2; iGain++) {
 553         // length of arrays
 554         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
 555         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
 556         // content of arrys
 557         for (int j = 0; j < nBins; j++) {
 558           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 559           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 560         }
 561       }
 562     }//iStrip
 563   }
 564
 565   // OK, now loop over the TTrees and fill the arrays needed for R(T)
 566   TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 567   TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 568
 569   int ChannelNum; // for regular towers
 570   int RefNum; // for LED
 571   double Hour;
 572   double AvgAmp;
 573
 574   TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
 575   TAvg->SetBranchAddress("fHour", &Hour);
 576   TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 577
 578   int iBin = 0;
 579   // counters for how many values were seen per SuperModule
 580   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
 581   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
 582
 583   for (int ient=0; ient<TAvg->GetEntries(); ient++) {
 584     TAvg->GetEntry(ient);
 585     // figure out where this info comes from
 586     fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
 587     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 588     // add value in the arrays
 589     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
 590     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
 591     nCount[iSM]++;
 592   }
 593
 594   TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
 595   TLEDRefAvg->SetBranchAddress("fHour", &Hour);
 596   TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 597
 598   for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
 599     TLEDRefAvg->GetEntry(ient);
 600     // figure out where this info comes from
 601     fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
 602     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 603     // add value in the arrays
 604     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
 605     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
 606     nCountLEDRef[iSM]++;
 607   }
 608
 609   // Normalize TArray values, and calculate average also
 610   Float_t norm = 0; // extra var, for readability
 611
 612   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 613     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
 614       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 615         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 616         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 617         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 618           for (iGain = 0; iGain < 2; iGain++) {
 619             // content of arrys
 620             for (int j = 0; j < nBins; j++) {
 621               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
 622                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
 623                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
 624               }
 625             }
 626           }
 627         }
 628       }//iCol
 629       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 630         for (iGain = 0; iGain < 2; iGain++) {
 631           for (int j = 0; j < nBins; j++) {
 632             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
 633               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
 634               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
 635             }
 636           }
 637         }
 638       }//iStrip
 639     }
 640   } // iSM
 641
 642
 643   // Calculate correction values, and store them
 644   // set fkErrorCode values for those that could not be set
 645
 646   Float_t Rt0 = 0;
 647   Float_t RT = 0;
 648   Float_t correction = 0; // c(T) = R(t0)/R(T)
 649   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
 650
 651   for (int i = 0; i < nSM; i++) {
 652     iSM = CalibAbsData[i].fSuperModuleNum;
 653     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 654       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 655       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 656       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 657
 658         // Calc. R(t0):
 659         AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
 660         iGain = v.fHighLow; // gain value used for abs. calibration
 661         refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration
 662
 663         // valid amplitude values should be larger than 0
 664         if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
 665           Rt0 =  v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
 666         }
 667         else {
 668           Rt0 = fkErrorCode;
 669         }
 670
 671         // Cal R(T)
 672         for (int j = 0; j < nBins; j++) {
 673
 674           // calculate R(T) also; first try with individual tower:
 675           // same gain as for abs. calibration is the default
 676           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
 677             // looks like valid data with the right gain combination
 678             RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
 679
 680             // if data appears to be saturated, and we are in high gain, then try with low gain instead
 681             if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
 682                  (ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1) ) { // presumably the gains should then both be changed.. can look into possibly only changing one in the future
 683               RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j);
 684               RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
 685             }
 686           }
 687           else {
 688             RT = fkErrorCode;
 689           }
 690
 691           // Calc. correction factor
 692           if (Rt0>0 && RT>0) {
 693             correction = Rt0/RT;
 694           }
 695           else {
 696             correction = fkErrorCode;
 697             nRemaining++;
 698           }
 699
 700           // Store the value
 701           CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 702           /* Check that
 703           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
 704           also works OK */
 705         } // nBins
 706       }
 707     }
 708   }
 709
 710   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
 711   return nRemaining;
 712 }
 713
 714 //________________________________________________________________
 715 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
 716 { // use averages for each strip if no good values exist for some single tower
 717
 718   // go over fTimeDepCorrection info
 719   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 720
 721   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 722
 723   // for calculating StripAverage info
 724   int nValidTower = 0;
 725   Float_t StripAverage = 0;
 726   Float_t val = 0;
 727
 728   int iSM = 0;
 729   int iCol = 0;
 730   int iRow = 0;
 731   int iStrip = 0;
 732   int firstCol = 0;
 733   int lastCol = 0;
 734
 735   for (int i = 0; i < nSM; i++) {
 736     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 737     for (int j = 0; j < nBins; j++) {
 738
 739       nValidTower = 0;
 740       StripAverage = 0;
 741
 742       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 743         firstCol = iStrip*2;
 744         if ((iSM%2)==1) { // C side
 745           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
 746         }
 747         lastCol = firstCol+1;
 748
 749         for (iCol = firstCol; iCol <= lastCol; iCol++) {
 750           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 751             val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 752             if (val>0) { // valid value; error code is negative
 753               StripAverage += val;
 754               nValidTower++;
 755             }
 756           }
 757         }
 758
 759         // calc average over strip
 760         if (nValidTower>0) {
 761           StripAverage /= nValidTower;
 762           for (iCol = firstCol; iCol <= lastCol; iCol++) {
 763             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 764               val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 765               if (val<0) { // invalid value; error code is negative
 766                 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
 767               }
 768             }
 769           }
 770         }
 771         else { // could not fill in unvalid entries
 772           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
 773         }
 774
 775       } // iStrip
 776     } // j, bins
 777   } // iSM
 778
 779   return nRemaining;
 780 }
 781
 782 //________________________________________________________________
 783 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins)
 784 { // OK, so we didn't have valid LED data that allowed us to do the correction only
 785   // with that info.
 786   // So, instead we'll rely on the temperature info and try to do the correction
 787   // based on that instead.
 788   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
 789   Int_t nRemaining = 0;
 790
 791   // info containers
 792   AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
 793   AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
 794   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 795   // correction container
 796   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 797
 798   int iSM = 0;
 799   int iCol = 0;
 800   int iRow = 0;
 801
 802   Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
 803   memset(TempCoeff, 0, sizeof(TempCoeff));
 804   Float_t MGain = 0;
 805   Double_t correction = 0;
 806   Double_t secondsPerBin = (1/fkSecToHour)*(1/fTimeBinsPerHour);
 807
 808   for (int i = 0; i < nSM; i++) {
 809     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 810
 811     // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
 812     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 813       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 814         AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
 815         MGain =  fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2],
 816                                        BiasAPDData[i].fVoltage[iCol][iRow]);
 817         TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
 818       }
 819     }
 820
 821     // figure out what the reference temperature is, from fCalibAbs
 822     Double_t ReferenceTemperature = 0;
 823     int nVal = 0;
 824     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
 825       if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
 826         ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
 827         nVal++;
 828       }
 829     }
 830
 831     if (nVal>0) {
 832       ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
 833
 834       for (int j = 0; j < nBins; j++) {
 835
 836         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
 837         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
 838       // 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)
 839         Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp);
 840
 841         Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
 842         // if the new temperature is higher than the old/reference one, then the gain has gone down
 843         if (fabs(TemperatureDiff)>fTemperatureResolution) {
 844           // significant enough difference that we need to consider it
 845
 846           // loop over all towers; effect of temperature change will depend on gain for this tower
 847           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 848             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 849
 850               correction = TemperatureDiff * TempCoeff[iCol][iRow];
 851               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 852             }
 853           }
 854
 855         } // if noteworthy temperature change
 856         else { // just set correction values to 1.0
 857           correction = 1;
 858           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 859             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 860               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 861             }
 862           }
 863         } // else
 864       } // j, Bins
 865
 866     } // if reference temperature exist
 867     else { // could not do the needed check.. signal that in the return code
 868       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 869     }
 870   } // iSM
 871
 872   return nRemaining;
 873 }
 874