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->ReadTreeBiasAPDInfo( (TTree *) 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   // stored object should be a TTree; read the info
 463   if (entry) {
 464     fCalibMapAPD->ReadTreeCalibMapAPDInfo( (TTree *) entry->GetObject() );
 465   }
 466
 467   if (fCalibMapAPD) {
 468     AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
 469   }
 470   else {
 471     AliWarning( Form("AliEMCALCalibMapAPD not found!") );
 472   }
 473
 474   return;
 475 }
 476
 477 //________________________________________________________________
 478 void AliEMCALCalibTimeDep::GetCalibAbsInfo()
 479 {
 480   // pick up Preprocessor output, based on fRun (most recent version)
 481   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
 482   if (entry) {
 483     fCalibAbs->ReadTreeCalibAbsInfo( (TTree *) entry->GetObject() );
 484   }
 485
 486   if (fCalibAbs) {
 487     AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
 488   }
 489   else {
 490     AliWarning( Form("AliEMCALCalibAbs not found!") );
 491   }
 492
 493   return;
 494 }
 495
 496 //________________________________________________________________
 497 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
 498 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
 499   // The correction factor we keep is c(T) = R(t0)/R(T)
 500   // T info from fCalibSignal, t0 info from fCalibAbs
 501
 502   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
 503   // but one could upgrade this in the future
 504   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 505
 506   // sanity check; same SuperModule indices for corrections as for regular calibrations
 507   AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 508   AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 509
 510   for (int i = 0; i < nSM; i++) {
 511     int iSMAbs = CalibAbsData[i].fSuperModuleNum;
 512     int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 513     if (iSMAbs != iSMCorr) {
 514       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
 515       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 516       return nRemaining;
 517     }
 518   }
 519
 520   int iSM = 0;
 521   int iCol = 0;
 522   int iRow = 0;
 523   int iStrip = 0;
 524   int iGain = 0;
 525
 526   // The fCalibSignal info is stored in TTrees
 527   // Note that the time-bins for the TTree's may not exactly match with our correction time bins
 528   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
 529   // fCalibSignal time info in seconds: Hour/fkSecToHour
 530   // corrected for startTime difference: Hour/fkSecToHour + timeDiff
 531   // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
 532
 533   // values for R(T), size of TArray = nBins
 534   // the [2] dimension below is for the low or high gain
 535   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 536   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
 537   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 538   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
 539
 540   // set up TArray's first
 541   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 542     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 543       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 544         for (iGain = 0; iGain < 2; iGain++) {
 545           // length of arrays
 546           ampT[iSM][iCol][iRow][iGain].Set(nBins);
 547           nT[iSM][iCol][iRow][iGain].Set(nBins);
 548           // content of arrys
 549           for (int j = 0; j < nBins; j++) {
 550             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
 551             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
 552           }
 553         }
 554       }
 555     }//iCol
 556     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 557       for (iGain = 0; iGain < 2; iGain++) {
 558         // length of arrays
 559         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
 560         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
 561         // content of arrys
 562         for (int j = 0; j < nBins; j++) {
 563           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 564           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
 565         }
 566       }
 567     }//iStrip
 568   }
 569
 570   // OK, now loop over the TTrees and fill the arrays needed for R(T)
 571   TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 572   TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
 573
 574   int ChannelNum; // for regular towers
 575   int RefNum; // for LED
 576   double Hour;
 577   double AvgAmp;
 578
 579   TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
 580   TAvg->SetBranchAddress("fHour", &Hour);
 581   TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 582
 583   int iBin = 0;
 584   // counters for how many values were seen per SuperModule
 585   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
 586   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
 587
 588   for (int ient=0; ient<TAvg->GetEntries(); ient++) {
 589     TAvg->GetEntry(ient);
 590     // figure out where this info comes from
 591     fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
 592     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 593     // add value in the arrays
 594     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
 595     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
 596     nCount[iSM]++;
 597   }
 598
 599   TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
 600   TLEDRefAvg->SetBranchAddress("fHour", &Hour);
 601   TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
 602
 603   for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
 604     TLEDRefAvg->GetEntry(ient);
 605     // figure out where this info comes from
 606     fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
 607     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
 608     // add value in the arrays
 609     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
 610     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
 611     nCountLEDRef[iSM]++;
 612   }
 613
 614   // Normalize TArray values, and calculate average also
 615   Float_t norm = 0; // extra var, for readability
 616
 617   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
 618     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
 619       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 620         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 621         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 622         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 623           for (iGain = 0; iGain < 2; iGain++) {
 624             // content of arrys
 625             for (int j = 0; j < nBins; j++) {
 626               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
 627                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
 628                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
 629               }
 630             }
 631           }
 632         }
 633       }//iCol
 634       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 635         for (iGain = 0; iGain < 2; iGain++) {
 636           for (int j = 0; j < nBins; j++) {
 637             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
 638               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
 639               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
 640             }
 641           }
 642         }
 643       }//iStrip
 644     }
 645   } // iSM
 646
 647
 648   // Calculate correction values, and store them
 649   // set fkErrorCode values for those that could not be set
 650
 651   Float_t Rt0 = 0;
 652   Float_t RT = 0;
 653   Float_t correction = 0; // c(T) = R(t0)/R(T)
 654   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
 655
 656   for (int i = 0; i < nSM; i++) {
 657     iSM = CalibAbsData[i].fSuperModuleNum;
 658     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 659       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 660       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
 661       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 662
 663         // Calc. R(t0):
 664         AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
 665         iGain = v.fHighLow; // gain value used for abs. calibration
 666         refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration
 667
 668         // valid amplitude values should be larger than 0
 669         if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
 670           Rt0 =  v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
 671         }
 672         else {
 673           Rt0 = fkErrorCode;
 674         }
 675
 676         // Cal R(T)
 677         for (int j = 0; j < nBins; j++) {
 678
 679           // calculate R(T) also; first try with individual tower:
 680           // same gain as for abs. calibration is the default
 681           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
 682             // looks like valid data with the right gain combination
 683             RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
 684
 685             // if data appears to be saturated, and we are in high gain, then try with low gain instead
 686             if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
 687                  (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
 688               RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j);
 689               RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
 690             }
 691           }
 692           else {
 693             RT = fkErrorCode;
 694           }
 695
 696           // Calc. correction factor
 697           if (Rt0>0 && RT>0) {
 698             correction = Rt0/RT;
 699           }
 700           else {
 701             correction = fkErrorCode;
 702             nRemaining++;
 703           }
 704
 705           // Store the value
 706           CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 707           /* Check that
 708           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
 709           also works OK */
 710         } // nBins
 711       }
 712     }
 713   }
 714
 715   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
 716   return nRemaining;
 717 }
 718
 719 //________________________________________________________________
 720 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
 721 { // use averages for each strip if no good values exist for some single tower
 722
 723   // go over fTimeDepCorrection info
 724   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
 725
 726   AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 727
 728   // for calculating StripAverage info
 729   int nValidTower = 0;
 730   Float_t StripAverage = 0;
 731   Float_t val = 0;
 732
 733   int iSM = 0;
 734   int iCol = 0;
 735   int iRow = 0;
 736   int iStrip = 0;
 737   int firstCol = 0;
 738   int lastCol = 0;
 739
 740   for (int i = 0; i < nSM; i++) {
 741     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 742     for (int j = 0; j < nBins; j++) {
 743
 744       nValidTower = 0;
 745       StripAverage = 0;
 746
 747       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
 748         firstCol = iStrip*2;
 749         if ((iSM%2)==1) { // C side
 750           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
 751         }
 752         lastCol = firstCol+1;
 753
 754         for (iCol = firstCol; iCol <= lastCol; iCol++) {
 755           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 756             val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 757             if (val>0) { // valid value; error code is negative
 758               StripAverage += val;
 759               nValidTower++;
 760             }
 761           }
 762         }
 763
 764         // calc average over strip
 765         if (nValidTower>0) {
 766           StripAverage /= nValidTower;
 767           for (iCol = firstCol; iCol <= lastCol; iCol++) {
 768             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 769               val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
 770               if (val<0) { // invalid value; error code is negative
 771                 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
 772               }
 773             }
 774           }
 775         }
 776         else { // could not fill in unvalid entries
 777           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
 778         }
 779
 780       } // iStrip
 781     } // j, bins
 782   } // iSM
 783
 784   return nRemaining;
 785 }
 786
 787 //________________________________________________________________
 788 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins)
 789 { // OK, so we didn't have valid LED data that allowed us to do the correction only
 790   // with that info.
 791   // So, instead we'll rely on the temperature info and try to do the correction
 792   // based on that instead.
 793   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
 794   Int_t nRemaining = 0;
 795
 796   // info containers
 797   AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
 798   AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
 799   AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
 800   // correction container
 801   AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
 802
 803   int iSM = 0;
 804   int iCol = 0;
 805   int iRow = 0;
 806
 807   Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
 808   memset(TempCoeff, 0, sizeof(TempCoeff));
 809   Float_t MGain = 0;
 810   Double_t correction = 0;
 811   Double_t secondsPerBin = (3600/fTimeBinsPerHour);
 812
 813   for (int i = 0; i < nSM; i++) {
 814     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
 815
 816     // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
 817     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 818       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 819         AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
 820         MGain =  fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2],
 821                                        BiasAPDData[i].fVoltage[iCol][iRow]);
 822         TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
 823       }
 824     }
 825
 826     // figure out what the reference temperature is, from fCalibAbs
 827     Double_t ReferenceTemperature = 0;
 828     int nVal = 0;
 829     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
 830       if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
 831         ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
 832         nVal++;
 833       }
 834     }
 835
 836     if (nVal>0) {
 837       ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
 838
 839       for (int j = 0; j < nBins; j++) {
 840
 841         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
 842         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
 843       // 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)
 844         Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp);
 845
 846         Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
 847         // if the new temperature is higher than the old/reference one, then the gain has gone down
 848         if (fabs(TemperatureDiff)>fTemperatureResolution) {
 849           // significant enough difference that we need to consider it
 850
 851           // loop over all towers; effect of temperature change will depend on gain for this tower
 852           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 853             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 854
 855               correction = TemperatureDiff * TempCoeff[iCol][iRow];
 856               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 857             }
 858           }
 859
 860         } // if noteworthy temperature change
 861         else { // just set correction values to 1.0
 862           correction = 1;
 863           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
 864             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
 865               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
 866             }
 867           }
 868         } // else
 869       } // j, Bins
 870
 871     } // if reference temperature exist
 872     else { // could not do the needed check.. signal that in the return code
 873       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
 874     }
 875   } // iSM
 876
 877   return nRemaining;
 878 }
 879