EMCAL/AliEMCALCalibTimeDep.cxx

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