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