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   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
 323   // set up a reasonable default (correction = 1.0)
 324   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
 325   fCalibTimeDepCorrection->SetStartTime(fStartTime);
 326   fCalibTimeDepCorrection->SetNTimeBins(nBins);
 327   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
 328
 329   // 1+2: try with LED corrections
 330   Int_t nRemaining = CalcLEDCorrection(nSM, nBins);
 331
 332   // 3: try with Temperature, if needed
 333   if (nRemaining>0) {
 334     nRemaining = CalcTemperatureCorrection(nSM, nBins);
 335   }
 336
 337   return nRemaining;
 338 }
 339
 340
 341 //________________________________________________________________
 342 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
 343 { // estimate the Temperature Coefficient, based on the dark current (IDark)
 344   // and the gain (M), based on Catania parameterizations
 345
 346   Double_t P0 = fkTempCoeffP0Const + fkTempCoeffP0Factor * IDark;
 347   Double_t P1 = fkTempCoeffP1Const + fkTempCoeffP1Factor * IDark;
 348
 349   Double_t TC = P0 + P1*M;
 350
 351   return TC;
 352 }
 353
 354 /* Next come the methods that do the work in picking up all the needed info..*/
 355 //________________________________________________________________
 356 void AliEMCALCalibTimeDep::GetTemperatureInfo()
 357 {
 358   // pick up Preprocessor output, based on fRun (most recent version)
 359   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
 360   if (entry) {
 361     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
 362   }
 363
 364   if (fTempArray) {
 365     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
 366                   fTempArray->NumSensors(),
 367                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
 368   }
 369   else {
 370     AliWarning( Form("AliEMCALSensorTempArray not found!") );
 371   }
 372
 373   return;
 374 }
 375
 376 //________________________________________________________________
 377 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
 378 {// assign max/min time and temperature values
 379
 380   fMinTemp = 999; // init to some large value (999 deg C)
 381   fMaxTemp = 0;
 382   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
 383   fMaxTime = 0;
 384
 385   Int_t n = 0; // number of valid readings
 386
 387   for (int i=0; i<fTempArray->NumSensors(); i++) {
 388
 389     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
 390
 391     // check time ranges
 392     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
 393     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
 394
 395     // check temperature ranges
 396     TGraph *g = st->GetGraph();
 397     if (g) { // ok, looks like we have valid data/info
 398       // let's check what the expected value at the time appears to be
 399       if (fMinTemp > g->GetMinimum()) { fMinTemp = g->GetMinimum(); }
 400       if (fMaxTemp < g->GetMaximum()) { fMaxTemp = g->GetMaximum(); }
 401       n++;
 402     }
 403   } // loop over fTempArray
 404
 405   if (n>0) { // some valid data was found
 406     return n;
 407   }
 408   else { // no good data
 409     return (Int_t) fkErrorCode;
 410   }
 411
 412 }
 413
 414 //________________________________________________________________
 415 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
 416 {
 417   // pick up Preprocessor output, based on fRun (most recent version)
 418   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
 419   if (entry) {
 420     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
 421   }
 422
 423   if (fCalibSignal) {
 424     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %d AvgEntries LEDRefEntries %d LEDRefAvgEntries %d",
 425                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
 426                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
 427                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
 428                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
 429                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
 430   }
 431   else {
 432     AliWarning( Form("AliCaloCalibSignal not found!") );
 433   }
 434
 435   return;
 436 }
 437
 438 //________________________________________________________________
 439 void AliEMCALCalibTimeDep::GetBiasAPDInfo()
 440 {
 441   // pick up Preprocessor output, based on fRun (most recent version)
 442   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/BiasAPD", fRun);
 443   if (entry) {
 444     fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
 445   }
 446
 447   if (fBiasAPD) {
 448     AliInfo( Form("BiasAPD: NSuperModule %d ", fBiasAPD->GetNSuperModule() ) );
 449   }
 450   else {
 451     AliWarning( Form("AliEMCALBiasAPD not found!") );
 452   }
 453
 454   return;
 455 }
 456
 457 //________________________________________________________________
 458 void AliEMCALCalibTimeDep::GetCalibMapAPDInfo()
 459 {
 460   // pick up Preprocessor output, based on fRun (most recent version)
 461   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
 462   if (entry) {
 463     fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
 464   }
 465
 466   if (fCalibMapAPD) {
 467     AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
 468   }
 469   else {
 470     AliWarning( Form("AliEMCALCalibMapAPD not found!") );
 471   }
 472
 473   return;
 474 }
 475
 476 //________________________________________________________________
 477 void AliEMCALCalibTimeDep::GetCalibAbsInfo()
 478 {
 479   // pick up Preprocessor output, based on fRun (most recent version)
 480   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
 481   if (entry) {
 482     fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject();
 483   }
 484
 485   if (fCalibAbs) {
 486     AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
 487   }
 488   else {
 489     AliWarning( Form("AliEMCALCalibAbs not found!") );
 490   }
 491
 492   return;
 493 }
 494
 495 //________________________________________________________________
 496 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
 497 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
 498   // The correction factor we keep is c(T) = R(t0)/R(T)
 499   // T info from fCalibSignal, t0 info from fCalibAbs
 500
 501   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
 502   // but one could upgrade this in the future
 503   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 504
 505   // sanity check; same SuperModule indices for corrections as for regular calibrations
 506   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 507   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 508
 509   for (int i = 0; i < nSM; i++) {
 510     int iSMAbs = CalibAbsData[i].fSuperModuleNum;
 511     int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 512     if (iSMAbs != iSMCorr) {
 513       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
 514       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 515       return nRemaining;
 516     }
 517   }
 518
 519   int iSM = 0;
 520   int iCol = 0;
 521   int iRow = 0;
 522   int iStrip = 0;
 523   int iGain = 0;
 524
 525   // The fCalibSignal info is stored in TTrees
 526   // Note that the time-bins for the TTree's may not exactly match with our correction time bins
 527   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
 528   // fCalibSignal time info in seconds: Hour/fkSecToHour
 529   // corrected for startTime difference: Hour/fkSecToHour + timeDiff
 530   // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
 531
 532   // values for R(T), size of TArray = nBins
 533   // the [2] dimension below is for the low or high gain
 534   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 535   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 536   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 537   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 538
 539   // set up TArray's first
 540   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 541     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 542       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 543         for (iGain = 0; iGain < 2; iGain++) {
 544           // length of arrays
 545           ampT[iSM][iCol][iRow][iGain].Set(nBins);
 546           nT[iSM][iCol][iRow][iGain].Set(nBins);
 547           // content of arrys
 548           for (int j = 0; j < nBins; j++) {
 549             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
 550             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
 551           }
 552         }
 553       }
 554     }//iCol
 555     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 556       for (iGain = 0; iGain < 2; iGain++) {
 557         // length of arrays
 558         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
 559         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
 560         // content of arrys
 561         for (int j = 0; j < nBins; j++) {
 562           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 563           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 564         }
 565       }
 566     }//iStrip
 567   }
 568
 569   // OK, now loop over the TTrees and fill the arrays needed for R(T)
 570   TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 571   TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 572
 573   int ChannelNum; // for regular towers
 574   int RefNum; // for LED
 575   double Hour;
 576   double AvgAmp;
 577
 578   TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
 579   TAvg->SetBranchAddress("fHour", &Hour);
 580   TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 581
 582   int iBin = 0;
 583   // counters for how many values were seen per SuperModule
 584   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
 585   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
 586
 587   for (int ient=0; ient<TAvg->GetEntries(); ient++) {
 588     TAvg->GetEntry(ient);
 589     // figure out where this info comes from
 590     fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
 591     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 592     // add value in the arrays
 593     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
 594     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
 595     nCount[iSM]++;
 596   }
 597
 598   TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
 599   TLEDRefAvg->SetBranchAddress("fHour", &Hour);
 600   TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 601
 602   for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
 603     TLEDRefAvg->GetEntry(ient);
 604     // figure out where this info comes from
 605     fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
 606     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 607     // add value in the arrays
 608     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
 609     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
 610     nCountLEDRef[iSM]++;
 611   }
 612
 613   // Normalize TArray values, and calculate average also
 614   Float_t norm = 0; // extra var, for readability
 615
 616   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 617     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
 618       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 619         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 620         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 621         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 622           for (iGain = 0; iGain < 2; iGain++) {
 623             // content of arrys
 624             for (int j = 0; j < nBins; j++) {
 625               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
 626                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
 627                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
 628               }
 629             }
 630           }
 631         }
 632       }//iCol
 633       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 634         for (iGain = 0; iGain < 2; iGain++) {
 635           for (int j = 0; j < nBins; j++) {
 636             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
 637               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
 638               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
 639             }
 640           }
 641         }
 642       }//iStrip
 643     }
 644   } // iSM
 645
 646
 647   // Calculate correction values, and store them
 648   // set fkErrorCode values for those that could not be set
 649
 650   Float_t Rt0 = 0;
 651   Float_t RT = 0;
 652   Float_t correction = 0; // c(T) = R(t0)/R(T)
 653   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
 654
 655   for (int i = 0; i < nSM; i++) {
 656     iSM = CalibAbsData[i].fSuperModuleNum;
 657     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 658       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 659       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 660       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 661
 662         // Calc. R(t0):
 663         AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
 664         iGain = v.fHighLow; // gain value used for abs. calibration
 665         refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration
 666
 667         // valid amplitude values should be larger than 0
 668         if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
 669           Rt0 =  v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
 670         }
 671         else {
 672           Rt0 = fkErrorCode;
 673         }
 674
 675         // Cal R(T)
 676         for (int j = 0; j < nBins; j++) {
 677
 678           // calculate R(T) also; first try with individual tower:
 679           // same gain as for abs. calibration is the default
 680           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
 681             // looks like valid data with the right gain combination
 682             RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
 683
 684             // if data appears to be saturated, and we are in high gain, then try with low gain instead
 685             if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
 686                  (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
 687               RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j);
 688               RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
 689             }
 690           }
 691           else {
 692             RT = fkErrorCode;
 693           }
 694
 695           // Calc. correction factor
 696           if (Rt0>0 && RT>0) {
 697             correction = Rt0/RT;
 698           }
 699           else {
 700             correction = fkErrorCode;
 701             nRemaining++;
 702           }
 703
 704           // Store the value
 705           CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 706           /* Check that
 707           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
 708           also works OK */
 709         } // nBins
 710       }
 711     }
 712   }
 713
 714   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
 715   return nRemaining;
 716 }
 717
 718 //________________________________________________________________
 719 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
 720 { // use averages for each strip if no good values exist for some single tower
 721
 722   // go over fTimeDepCorrection info
 723   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 724
 725   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 726
 727   // for calculating StripAverage info
 728   int nValidTower = 0;
 729   Float_t StripAverage = 0;
 730   Float_t val = 0;
 731
 732   int iSM = 0;
 733   int iCol = 0;
 734   int iRow = 0;
 735   int iStrip = 0;
 736   int firstCol = 0;
 737   int lastCol = 0;
 738
 739   for (int i = 0; i < nSM; i++) {
 740     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 741     for (int j = 0; j < nBins; j++) {
 742
 743       nValidTower = 0;
 744       StripAverage = 0;
 745
 746       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 747         firstCol = iStrip*2;
 748         if ((iSM%2)==1) { // C side
 749           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
 750         }
 751         lastCol = firstCol+1;
 752
 753         for (iCol = firstCol; iCol <= lastCol; iCol++) {
 754           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 755             val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 756             if (val>0) { // valid value; error code is negative
 757               StripAverage += val;
 758               nValidTower++;
 759             }
 760           }
 761         }
 762
 763         // calc average over strip
 764         if (nValidTower>0) {
 765           StripAverage /= nValidTower;
 766           for (iCol = firstCol; iCol <= lastCol; iCol++) {
 767             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 768               val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 769               if (val<0) { // invalid value; error code is negative
 770                 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
 771               }
 772             }
 773           }
 774         }
 775         else { // could not fill in unvalid entries
 776           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
 777         }
 778
 779       } // iStrip
 780     } // j, bins
 781   } // iSM
 782
 783   return nRemaining;
 784 }
 785
 786 //________________________________________________________________
 787 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins)
 788 { // OK, so we didn't have valid LED data that allowed us to do the correction only
 789   // with that info.
 790   // So, instead we'll rely on the temperature info and try to do the correction
 791   // based on that instead.
 792   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
 793   Int_t nRemaining = 0;
 794
 795   // info containers
 796   AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
 797   AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
 798   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 799   // correction container
 800   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 801
 802   int iSM = 0;
 803   int iCol = 0;
 804   int iRow = 0;
 805
 806   Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
 807   memset(TempCoeff, 0, sizeof(TempCoeff));
 808   Float_t MGain = 0;
 809   Double_t correction = 0;
 810   Double_t secondsPerBin = (3600/fTimeBinsPerHour);
 811
 812   for (int i = 0; i < nSM; i++) {
 813     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 814
 815     // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
 816     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 817       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 818         AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
 819         MGain =  fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2],
 820                                        BiasAPDData[i].fVoltage[iCol][iRow]);
 821         TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
 822       }
 823     }
 824
 825     // figure out what the reference temperature is, from fCalibAbs
 826     Double_t ReferenceTemperature = 0;
 827     int nVal = 0;
 828     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
 829       if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
 830         ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
 831         nVal++;
 832       }
 833     }
 834
 835     if (nVal>0) {
 836       ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
 837
 838       for (int j = 0; j < nBins; j++) {
 839
 840         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
 841         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
 842       // 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)
 843         Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp);
 844
 845         Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
 846         // if the new temperature is higher than the old/reference one, then the gain has gone down
 847         if (fabs(TemperatureDiff)>fTemperatureResolution) {
 848           // significant enough difference that we need to consider it
 849
 850           // loop over all towers; effect of temperature change will depend on gain for this tower
 851           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 852             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 853
 854               correction = TemperatureDiff * TempCoeff[iCol][iRow];
 855               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 856             }
 857           }
 858
 859         } // if noteworthy temperature change
 860         else { // just set correction values to 1.0
 861           correction = 1;
 862           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 863             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 864               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 865             }
 866           }
 867         } // else
 868       } // j, Bins
 869
 870     } // if reference temperature exist
 871     else { // could not do the needed check.. signal that in the return code
 872       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 873     }
 874   } // iSM
 875
 876   return nRemaining;
 877 }
 878