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