Move AliWarning and AliInfo to AliDebug to reduce prints during reconstruction
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALCalibTimeDep.cxx
CommitLineData
a42992b7 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// //
0ce5c45d 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
a42992b7 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"
d81e6423 39#include "AliCaloCalibSignal.h"
40#include "AliEMCALBiasAPD.h"
41#include "AliEMCALCalibMapAPD.h"
82d90a2f 42#include "AliEMCALCalibReference.h"
d81e6423 43#include "AliEMCALCalibTimeDepCorrection.h"
a42992b7 44#include "AliEMCALCalibTimeDep.h"
45
46/* first a bunch of constants.. */
0ce5c45d 47const double kSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
a42992b7 48
d81e6423 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
0ce5c45d 51const double kTempCoeffP0Const = -0.903; //
52const double kTempCoeffP0Factor = -1.381e7; //
53const double kTempCoeffP1Const = -0.023; //
54const double kTempCoeffP1Factor = -4.966e5; //
d81e6423 55
0ce5c45d 56const double kErrorCode = -999; // to indicate that something went wrong
621ff010 57
a42992b7 58using namespace std;
59
60ClassImp(AliEMCALCalibTimeDep)
61
62//________________________________________________________________
63AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
64 fRun(0),
65 fStartTime(0),
66 fEndTime(0),
67 fMinTemp(0),
68 fMaxTemp(0),
69 fMinTime(0),
70 fMaxTime(0),
d81e6423 71 fTemperatureResolution(0.1), // 0.1 deg C is default
72 fTimeBinsPerHour(2), // 2 30-min bins per hour is default
82d90a2f 73 fHighLowGainFactor(16), // factor ~16 between High gain and low gain
d81e6423 74 fTempArray(NULL),
75 fCalibSignal(NULL),
76 fBiasAPD(NULL),
77 fCalibMapAPD(NULL),
82d90a2f 78 fCalibReference(NULL),
d81e6423 79 fCalibTimeDepCorrection(NULL)
a42992b7 80{
81 // Constructor
82}
83
84//________________________________________________________________
85AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
621ff010 86 TObject(calibt),
a42992b7 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()),
d81e6423 94 fTemperatureResolution(calibt.GetTemperatureResolution()),
95 fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
82d90a2f 96 fHighLowGainFactor(calibt.GetHighLowGainFactor()),
d81e6423 97 fTempArray(calibt.GetTempArray()),
98 fCalibSignal(calibt.GetCalibSignal()),
99 fBiasAPD(calibt.GetBiasAPD()),
100 fCalibMapAPD(calibt.GetCalibMapAPD()),
82d90a2f 101 fCalibReference(calibt.GetCalibReference()),
d81e6423 102 fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection())
a42992b7 103{
104 // copy constructor
105}
106
107
108//________________________________________________________________
109AliEMCALCalibTimeDep &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//________________________________________________________________
119AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
120{
121 // Destructor
122}
123
124//________________________________________________________________
125void 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;
d81e6423 135 fTemperatureResolution = 0.1; // 0.1 deg C is default
136 fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
a42992b7 137 fTempArray = NULL;
d81e6423 138 fCalibSignal = NULL;
139 fBiasAPD = NULL;
140 fCalibMapAPD = NULL;
82d90a2f 141 fCalibReference = NULL;
d81e6423 142 fCalibTimeDepCorrection = NULL;
a42992b7 143 return;
144}
145
146//________________________________________________________________
621ff010 147void AliEMCALCalibTimeDep::PrintInfo() const
a42992b7 148{
149 // print some info
d81e6423 150 cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
a42992b7 151 // basic variables, all 'publicly available' also
152 cout << " VARIABLE DUMP: " << endl
153 << " GetStartTime() " << GetStartTime() << endl
154 << " GetEndTime() " << GetEndTime() << endl
155 << " GetMinTemp() " << GetMinTemp() << endl
d81e6423 156 << " GetMaxTemp() " << GetMaxTemp() << endl;
a42992b7 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;
a42992b7 163
164 return;
165}
d81e6423 166
a42992b7 167//________________________________________________________________
d81e6423 168Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
a42992b7 169{
0ce5c45d 170 return (fEndTime - fStartTime)*kSecToHour;
a42992b7 171}
d81e6423 172
173//________________________________________________________________
174Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
175{
0ce5c45d 176 return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
d81e6423 177}
178
a42992b7 179//________________________________________________________________
d81e6423 180Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
a42992b7 181{
0ce5c45d 182 return (fMaxTime - fMinTime)*kSecToHour;
a42992b7 183}
d81e6423 184
a42992b7 185//________________________________________________________________
d81e6423 186Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
a42992b7 187{
188 return (fMaxTemp - fMinTemp);
189}
190
191//________________________________________________________________
192void AliEMCALCalibTimeDep::Initialize(Int_t run,
193 UInt_t startTime, UInt_t endTime)
0ce5c45d 194{ // setup, and get temperature info
a42992b7 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
d81e6423 203 ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
a42992b7 204
205 return;
206}
207
208//________________________________________________________________
d81e6423 209Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const
a42992b7 210{// return estimate for all SuperModules and sensors, that had data
211
212 // first convert from seconds to hours..
0ce5c45d 213 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
621ff010 214
d81e6423 215 Double_t average = 0;
621ff010 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
d81e6423 227 Double_t val = f->Eval(timeHour);
621ff010 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
0ce5c45d 239 return kErrorCode;
621ff010 240 }
a42992b7 241
a42992b7 242}
243
244//________________________________________________________________
d81e6423 245Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
a42992b7 246{// return estimate for this one SuperModule, if it had data
247
248 // first convert from seconds to hours..
0ce5c45d 249 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
621ff010 250
d81e6423 251 Double_t average = 0;
621ff010 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
d81e6423 264 Double_t val = f->Eval(timeHour);
621ff010 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
0ce5c45d 279 return kErrorCode;
621ff010 280 }
a42992b7 281
a42992b7 282}
283
284//________________________________________________________________
d81e6423 285Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_t timeStamp) const
a42992b7 286{// return estimate for this one SuperModule and sensor, if it had data
621ff010 287
a42992b7 288 // first convert from seconds to hours..
0ce5c45d 289 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
621ff010 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
d81e6423 301 Double_t val = f->Eval(timeHour);
621ff010 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
a42992b7 309
621ff010 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
0ce5c45d 312 return kErrorCode;
621ff010 313
a42992b7 314}
315
316//________________________________________________________________
d81e6423 317Int_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?
82d90a2f 328 Int_t nSM = fCalibReference->GetNSuperModule();
d81e6423 329 // how many time-bins should we have for this run?
330 Int_t nBins = (Int_t) GetLengthOfRunInBins(); // round-down (from double to int)
220ed45a 331 Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
d81e6423 332 // set up a reasonable default (correction = 1.0)
333 fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
220ed45a 334 fCalibTimeDepCorrection->SetStartTime(fStartTime);
335 fCalibTimeDepCorrection->SetNTimeBins(nBins);
336 fCalibTimeDepCorrection->SetTimeBinSize(binSize);
d81e6423 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;
a42992b7 347}
348
d81e6423 349
350//________________________________________________________________
351Double_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
0ce5c45d 355 Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
356 Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
d81e6423 357
0ce5c45d 358 Double_t dTC = dP0 + dP1*M;
d81e6423 359
0ce5c45d 360 return dTC;
d81e6423 361}
362
363/* Next come the methods that do the work in picking up all the needed info..*/
a42992b7 364//________________________________________________________________
365void 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
621ff010 373 if (fTempArray) {
a42992b7 374 AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
375 fTempArray->NumSensors(),
376 fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
621ff010 377 }
378 else {
379 AliWarning( Form("AliEMCALSensorTempArray not found!") );
a42992b7 380 }
381
382 return;
383}
d81e6423 384
385//________________________________________________________________
386Int_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
0ce5c45d 418 return (Int_t) kErrorCode;
d81e6423 419 }
420
421}
422
423//________________________________________________________________
424void 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//________________________________________________________________
448void 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) {
2f17a269 453 fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
d81e6423 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//________________________________________________________________
467void AliEMCALCalibTimeDep::GetCalibMapAPDInfo()
468{
469 // pick up Preprocessor output, based on fRun (most recent version)
470 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
61917ab3 471 // stored object should be a TTree; read the info
d81e6423 472 if (entry) {
2f17a269 473 fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
d81e6423 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//________________________________________________________________
82d90a2f 487void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
d81e6423 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) {
82d90a2f 492 fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
d81e6423 493 }
494
82d90a2f 495 if (fCalibReference) {
496 AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
d81e6423 497 }
498 else {
82d90a2f 499 AliWarning( Form("AliEMCALCalibReference not found!") );
d81e6423 500 }
501
502 return;
503}
504
505//________________________________________________________________
506Int_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)
82d90a2f 509 // T info from fCalibSignal, t0 info from fCalibReference
d81e6423 510
82d90a2f 511 // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
d81e6423 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
d81e6423 516 for (int i = 0; i < nSM; i++) {
0ce5c45d 517 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
518 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
2f17a269 519
0ce5c45d 520 int iSMRef = dataCalibReference->GetSuperModuleNum();
521 int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
82d90a2f 522 if (iSMRef != iSMCorr) {
523 AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
d81e6423 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
0ce5c45d 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
d81e6423 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)
0ce5c45d 580 TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
581 TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
d81e6423 582
0ce5c45d 583 int iChannelNum; // for regular towers
584 int iRefNum; // for LED
585 double dHour;
586 double dAvgAmp;
d81e6423 587
0ce5c45d 588 treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
589 treeAvg->SetBranchAddress("fHour", &dHour);
590 treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
d81e6423 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
0ce5c45d 597 for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
598 treeAvg->GetEntry(ient);
d81e6423 599 // figure out where this info comes from
0ce5c45d 600 fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
601 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
d81e6423 602 // add value in the arrays
0ce5c45d 603 ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
d81e6423 604 nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
605 nCount[iSM]++;
606 }
607
0ce5c45d 608 treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
609 treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
610 treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
d81e6423 611
0ce5c45d 612 for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
613 treeLEDRefAvg->GetEntry(ient);
d81e6423 614 // figure out where this info comes from
0ce5c45d 615 fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
616 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
d81e6423 617 // add value in the arrays
0ce5c45d 618 ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
d81e6423 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
0ce5c45d 658 // set kErrorCode values for those that could not be set
d81e6423 659
0ce5c45d 660 Float_t ratiot0 = 0;
661 Float_t ratioT = 0;
d81e6423 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++) {
0ce5c45d 666 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
667 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
668 iSM = dataCalibReference->GetSuperModuleNum();
2f17a269 669
d81e6423 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
0ce5c45d 673 refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
82d90a2f 674
d81e6423 675 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
676
677 // Calc. R(t0):
0ce5c45d 678 AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
82d90a2f 679 iGain = refVal->GetHighLow(); // gain value used for reference calibration
d81e6423 680 // valid amplitude values should be larger than 0
0ce5c45d 681 if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
682 ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
d81e6423 683 }
684 else {
0ce5c45d 685 ratiot0 = kErrorCode;
d81e6423 686 }
687
82d90a2f 688 // Calc. R(T)
d81e6423 689 for (int j = 0; j < nBins; j++) {
690
691 // calculate R(T) also; first try with individual tower:
82d90a2f 692 // same gain as for reference calibration is the default
d81e6423 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
0ce5c45d 695 ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
d81e6423 696
697 // if data appears to be saturated, and we are in high gain, then try with low gain instead
82d90a2f 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..
0ce5c45d 711 ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
82d90a2f 712
713 if (newGain<iGain) {
0ce5c45d 714 ratioT *= fHighLowGainFactor;
82d90a2f 715 }
716 else if (newRefGain<refGain) {
0ce5c45d 717 ratioT /= fHighLowGainFactor;
82d90a2f 718 }
d81e6423 719 }
720 }
721 else {
0ce5c45d 722 ratioT = kErrorCode;
d81e6423 723 }
724
725 // Calc. correction factor
0ce5c45d 726 if (ratiot0>0 && ratioT>0) {
727 correction = ratiot0/ratioT;
d81e6423 728 }
729 else {
0ce5c45d 730 correction = kErrorCode;
d81e6423 731 nRemaining++;
732 }
733
734 // Store the value
0ce5c45d 735 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
d81e6423 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//________________________________________________________________
749Int_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
d81e6423 755 // for calculating StripAverage info
756 int nValidTower = 0;
0ce5c45d 757 Float_t stripAverage = 0;
d81e6423 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++) {
0ce5c45d 768 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
769 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
2f17a269 770
d81e6423 771 for (int j = 0; j < nBins; j++) {
772
773 nValidTower = 0;
0ce5c45d 774 stripAverage = 0;
d81e6423 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++) {
0ce5c45d 785 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
d81e6423 786 if (val>0) { // valid value; error code is negative
0ce5c45d 787 stripAverage += val;
d81e6423 788 nValidTower++;
789 }
790 }
791 }
792
793 // calc average over strip
794 if (nValidTower>0) {
0ce5c45d 795 stripAverage /= nValidTower;
d81e6423 796 for (iCol = firstCol; iCol <= lastCol; iCol++) {
797 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
0ce5c45d 798 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
d81e6423 799 if (val<0) { // invalid value; error code is negative
0ce5c45d 800 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
d81e6423 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//________________________________________________________________
817Int_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
d81e6423 825 int iSM = 0;
826 int iCol = 0;
827 int iRow = 0;
828
0ce5c45d 829 Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
830 memset(dTempCoeff, 0, sizeof(dTempCoeff));
831 Float_t gainM = 0;
d81e6423 832 Double_t correction = 0;
220ed45a 833 Double_t secondsPerBin = (3600/fTimeBinsPerHour);
d81e6423 834
835 for (int i = 0; i < nSM; i++) {
0ce5c45d 836 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(iSM);
837 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
2f17a269 838
0ce5c45d 839 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
840 AliEMCALSuperModuleCalibMapAPD * dataCalibMapAPD = fCalibMapAPD->GetSuperModuleCalibMapAPDNum(iSM);
841 AliEMCALSuperModuleBiasAPD * dataBiasAPD = fBiasAPD->GetSuperModuleBiasAPDNum(iSM);
d81e6423 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++) {
0ce5c45d 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);
d81e6423 850 }
851 }
852
82d90a2f 853 // figure out what the reference temperature is, from fCalibReference
0ce5c45d 854 Double_t referenceTemperature = 0;
d81e6423 855 int nVal = 0;
856 for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
0ce5c45d 857 if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
858 referenceTemperature += dataCalibReference->GetTemperature(iSensor);
d81e6423 859 nVal++;
860 }
861 }
862
863 if (nVal>0) {
0ce5c45d 864 referenceTemperature /= nVal; // valid values exist, we can look into corrections
d81e6423 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)
0ce5c45d 871 Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp);
d81e6423 872
0ce5c45d 873 Double_t temperatureDiff = referenceTemperature - dSMTemperature; // old vs new
d81e6423 874 // if the new temperature is higher than the old/reference one, then the gain has gone down
0ce5c45d 875 if (fabs(temperatureDiff)>fTemperatureResolution) {
d81e6423 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
0ce5c45d 882 correction = temperatureDiff * dTempCoeff[iCol][iRow];
883 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
d81e6423 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++) {
0ce5c45d 892 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
d81e6423 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