]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALCalibTimeDep.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated SNM Glauber fit
[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
a42992b7 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 }
4cf7bc35 276 //cout << " n " << n << " nval " << nval << " median " << median << endl;
1740b16f 277 if (nval > 0) {
278 average /= nval;
4cf7bc35 279 //cout << " average " << average << endl;
1740b16f 280 return average;
281 }
282 else { // this case should not happen, but kept for completeness (coverity etc)
283 return median;
284 }
621ff010 285 }
286 else { // no good data
0ce5c45d 287 return kErrorCode;
621ff010 288 }
a42992b7 289
a42992b7 290}
291
a42992b7 292//________________________________________________________________
d81e6423 293Int_t AliEMCALCalibTimeDep::CalcCorrection()
294{ // OK, this is where the real action takes place - the heart of this class..
295 /* The philosophy is as follows:
716fca62 296 0. Init corrections to 1.0 values, and see how many correction bins we need
297 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 298 2. try to use temperature info + APD temperature coefficient info, to estimate correction.
299 For now (from Dec 2009), we do not use LED info.
d81e6423 300 */
301
302 // 0: Init
303 // how many SuperModules do we have?
82d90a2f 304 Int_t nSM = fCalibReference->GetNSuperModule();
d81e6423 305 // how many time-bins should we have for this run?
ed3db319 306 Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
220ed45a 307 Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
716fca62 308
309 // 1: get info on how much individual sensors might have changed during
310 // the run (compare max-min for each sensor separately)
311 if (fMaxTempVariation < fTemperatureResolution) {
312 nBins = 1; // just one bin needed..
2c62d6a3 313 }
314 if (nBins == 1) {
716fca62 315 binSize = fEndTime - fStartTime;
316 }
1740b16f 317 if (fVerbosity > 0) {
2c62d6a3 318 cout << " nBins " << nBins << " binSize " << binSize << endl;
319 }
716fca62 320
d81e6423 321 // set up a reasonable default (correction = 1.0)
716fca62 322 fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
d81e6423 323 fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
220ed45a 324 fCalibTimeDepCorrection->SetStartTime(fStartTime);
325 fCalibTimeDepCorrection->SetNTimeBins(nBins);
326 fCalibTimeDepCorrection->SetTimeBinSize(binSize);
d81e6423 327
716fca62 328 // 2: try with Temperature correction
329 Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
d81e6423 330
331 return nRemaining;
a42992b7 332}
333
d81e6423 334
335//________________________________________________________________
336Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
337{ // estimate the Temperature Coefficient, based on the dark current (IDark)
338 // and the gain (M), based on Catania parameterizations
339
0ce5c45d 340 Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
341 Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
d81e6423 342
0ce5c45d 343 Double_t dTC = dP0 + dP1*M;
ed3db319 344 // from % numbers to regular ones..:
345 dTC *= 0.01;
d81e6423 346
e1a60af4 347 return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
d81e6423 348}
349
350/* Next come the methods that do the work in picking up all the needed info..*/
a42992b7 351//________________________________________________________________
352void AliEMCALCalibTimeDep::GetTemperatureInfo()
353{
354 // pick up Preprocessor output, based on fRun (most recent version)
355 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
356 if (entry) {
357 fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
358 }
359
621ff010 360 if (fTempArray) {
a42992b7 361 AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
362 fTempArray->NumSensors(),
363 fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
621ff010 364 }
365 else {
366 AliWarning( Form("AliEMCALSensorTempArray not found!") );
a42992b7 367 }
368
369 return;
370}
d81e6423 371
372//________________________________________________________________
373Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
374{// assign max/min time and temperature values
375
376 fMinTemp = 999; // init to some large value (999 deg C)
377 fMaxTemp = 0;
716fca62 378 fMinTempVariation = 999; // init to some large value (999 deg C)
379 fMaxTempVariation = 0;
d81e6423 380 fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
381 fMaxTime = 0;
382
383 Int_t n = 0; // number of valid readings
384
385 for (int i=0; i<fTempArray->NumSensors(); i++) {
386
387 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
716fca62 388 if ( st->GetStartTime() == 0 ) { // no valid data
389 continue;
390 }
d81e6423 391
392 // check time ranges
393 if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
394 if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
716fca62 395
d81e6423 396 // check temperature ranges
716fca62 397 AliSplineFit *f = st->GetFit();
398
399 if (f) { // ok, looks like we have valid data/info
400 int np = f->GetKnots();
401 Double_t *y0 = f->GetY0();
402 // min and max values within the single sensor
403 Double_t min = 999;
404 Double_t max = 0;
1740b16f 405 int nval = 0;
716fca62 406 for (int ip=0; ip<np; ip++) {
1740b16f 407 if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
408 if (min > y0[ip]) { min = y0[ip]; }
409 if (max < y0[ip]) { max = y0[ip]; }
410 nval++;
411 }
412 }
413 if (nval>0) {
414 if (fMinTemp > min) { fMinTemp = min; }
415 if (fMaxTemp < max) { fMaxTemp = max; }
416 Double_t variation = max - min;
417 if (fMinTempVariation > variation) { fMinTempVariation = variation; }
418 if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
419
420 n++;
716fca62 421 }
d81e6423 422 }
423 } // loop over fTempArray
424
425 if (n>0) { // some valid data was found
426 return n;
427 }
428 else { // no good data
0ce5c45d 429 return (Int_t) kErrorCode;
d81e6423 430 }
431
432}
433
434//________________________________________________________________
435void AliEMCALCalibTimeDep::GetCalibSignalInfo()
436{
437 // pick up Preprocessor output, based on fRun (most recent version)
438 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
439 if (entry) {
440 fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
441 }
442
443 if (fCalibSignal) {
29b7e56e 444 AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
d81e6423 445 fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
446 fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
447 fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
29b7e56e 448 fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
d81e6423 449 fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
450 }
451 else {
452 AliWarning( Form("AliCaloCalibSignal not found!") );
453 }
454
455 return;
456}
457
458//________________________________________________________________
ed3db319 459void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo()
d81e6423 460{
461 // pick up Preprocessor output, based on fRun (most recent version)
ed3db319 462 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
61917ab3 463 // stored object should be a TTree; read the info
d81e6423 464 if (entry) {
ed3db319 465 fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
d81e6423 466 }
467
ed3db319 468 if (fCalibTempCoeff) {
469 AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
d81e6423 470 }
471 else {
ed3db319 472 AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
d81e6423 473 }
474
475 return;
476}
477
478//________________________________________________________________
82d90a2f 479void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
d81e6423 480{
481 // pick up Preprocessor output, based on fRun (most recent version)
716fca62 482 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
d81e6423 483 if (entry) {
82d90a2f 484 fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
d81e6423 485 }
486
82d90a2f 487 if (fCalibReference) {
488 AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
d81e6423 489 }
490 else {
82d90a2f 491 AliWarning( Form("AliEMCALCalibReference not found!") );
d81e6423 492 }
493
494 return;
495}
496
497//________________________________________________________________
498Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
499{// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
500 // The correction factor we keep is c(T) = R(t0)/R(T)
82d90a2f 501 // T info from fCalibSignal, t0 info from fCalibReference
d81e6423 502
82d90a2f 503 // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
d81e6423 504 // but one could upgrade this in the future
505 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
506
507 // sanity check; same SuperModule indices for corrections as for regular calibrations
d81e6423 508 for (int i = 0; i < nSM; i++) {
0ce5c45d 509 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
510 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
2f17a269 511
0ce5c45d 512 int iSMRef = dataCalibReference->GetSuperModuleNum();
513 int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
82d90a2f 514 if (iSMRef != iSMCorr) {
515 AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
d81e6423 516 nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
517 return nRemaining;
518 }
519 }
520
521 int iSM = 0;
522 int iCol = 0;
523 int iRow = 0;
524 int iStrip = 0;
525 int iGain = 0;
526
527 // The fCalibSignal info is stored in TTrees
528 // Note that the time-bins for the TTree's may not exactly match with our correction time bins
529 int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
0ce5c45d 530 // fCalibSignal time info in seconds: Hour/kSecToHour
531 // corrected for startTime difference: Hour/kSecToHour + timeDiff
532 // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
d81e6423 533
534 // values for R(T), size of TArray = nBins
535 // the [2] dimension below is for the low or high gain
536 TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
537 TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
538 TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
539 TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
540
541 // set up TArray's first
542 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
543 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
544 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
545 for (iGain = 0; iGain < 2; iGain++) {
546 // length of arrays
547 ampT[iSM][iCol][iRow][iGain].Set(nBins);
548 nT[iSM][iCol][iRow][iGain].Set(nBins);
549 // content of arrys
550 for (int j = 0; j < nBins; j++) {
551 ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
552 nT[iSM][iCol][iRow][iGain].AddAt(0, j);
553 }
554 }
555 }
556 }//iCol
557 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
558 for (iGain = 0; iGain < 2; iGain++) {
559 // length of arrays
560 ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
561 nLEDRefT[iSM][iStrip][iGain].Set(nBins);
562 // content of arrys
563 for (int j = 0; j < nBins; j++) {
564 ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
565 nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
566 }
567 }
568 }//iStrip
569 }
570
571 // OK, now loop over the TTrees and fill the arrays needed for R(T)
0ce5c45d 572 TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
573 TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
d81e6423 574
18831b6c 575 int iChannelNum = 0; // for regular towers
576 int iRefNum = 0; // for LED
577 double dHour = 0;
578 double dAvgAmp = 0;
d81e6423 579
0ce5c45d 580 treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
581 treeAvg->SetBranchAddress("fHour", &dHour);
582 treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
d81e6423 583
584 int iBin = 0;
585 // counters for how many values were seen per SuperModule
586 int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
587 int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
588
0ce5c45d 589 for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
590 treeAvg->GetEntry(ient);
d81e6423 591 // figure out where this info comes from
0ce5c45d 592 fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
593 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
d81e6423 594 // add value in the arrays
0ce5c45d 595 ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
d81e6423 596 nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
597 nCount[iSM]++;
598 }
599
0ce5c45d 600 treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
601 treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
602 treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
d81e6423 603
0ce5c45d 604 for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
605 treeLEDRefAvg->GetEntry(ient);
d81e6423 606 // figure out where this info comes from
0ce5c45d 607 fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
608 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
d81e6423 609 // add value in the arrays
0ce5c45d 610 ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
d81e6423 611 nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
612 nCountLEDRef[iSM]++;
613 }
614
615 // Normalize TArray values, and calculate average also
616 Float_t norm = 0; // extra var, for readability
617
618 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
619 if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
620 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
621 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
622 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
623 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
624 for (iGain = 0; iGain < 2; iGain++) {
625 // content of arrys
626 for (int j = 0; j < nBins; j++) {
627 if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
628 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
629 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
630 }
631 }
632 }
633 }
634 }//iCol
635 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
636 for (iGain = 0; iGain < 2; iGain++) {
637 for (int j = 0; j < nBins; j++) {
638 if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
639 norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
640 ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
641 }
642 }
643 }
644 }//iStrip
645 }
646 } // iSM
647
648
649 // Calculate correction values, and store them
0ce5c45d 650 // set kErrorCode values for those that could not be set
d81e6423 651
0ce5c45d 652 Float_t ratiot0 = 0;
653 Float_t ratioT = 0;
d81e6423 654 Float_t correction = 0; // c(T) = R(t0)/R(T)
655 Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
656
657 for (int i = 0; i < nSM; i++) {
0ce5c45d 658 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
659 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
660 iSM = dataCalibReference->GetSuperModuleNum();
2f17a269 661
d81e6423 662 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
663 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
664 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
0ce5c45d 665 refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
82d90a2f 666
d81e6423 667 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
668
669 // Calc. R(t0):
0ce5c45d 670 AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
82d90a2f 671 iGain = refVal->GetHighLow(); // gain value used for reference calibration
d81e6423 672 // valid amplitude values should be larger than 0
0ce5c45d 673 if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
674 ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
d81e6423 675 }
676 else {
0ce5c45d 677 ratiot0 = kErrorCode;
d81e6423 678 }
679
82d90a2f 680 // Calc. R(T)
d81e6423 681 for (int j = 0; j < nBins; j++) {
682
683 // calculate R(T) also; first try with individual tower:
82d90a2f 684 // same gain as for reference calibration is the default
d81e6423 685 if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
686 // looks like valid data with the right gain combination
0ce5c45d 687 ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
d81e6423 688
689 // if data appears to be saturated, and we are in high gain, then try with low gain instead
82d90a2f 690 int newGain = iGain;
691 int newRefGain = refGain;
692 if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
693 newGain = 0;
694 }
695 if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
696 newRefGain = 0;
697 }
698
699 if (newGain!=iGain || newRefGain!=refGain) {
700 // compensate for using different gain than in the reference calibration
701 // we may need to have a custom H/L ratio value for each tower
702 // later, but for now just use a common value, as the rest of the code does..
0ce5c45d 703 ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
82d90a2f 704
705 if (newGain<iGain) {
0ce5c45d 706 ratioT *= fHighLowGainFactor;
82d90a2f 707 }
708 else if (newRefGain<refGain) {
0ce5c45d 709 ratioT /= fHighLowGainFactor;
82d90a2f 710 }
d81e6423 711 }
712 }
713 else {
0ce5c45d 714 ratioT = kErrorCode;
d81e6423 715 }
716
717 // Calc. correction factor
0ce5c45d 718 if (ratiot0>0 && ratioT>0) {
719 correction = ratiot0/ratioT;
d81e6423 720 }
721 else {
0ce5c45d 722 correction = kErrorCode;
d81e6423 723 nRemaining++;
724 }
725
726 // Store the value
0ce5c45d 727 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
d81e6423 728 /* Check that
729 fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
730 also works OK */
731 } // nBins
732 }
733 }
734 }
735
736 nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
737 return nRemaining;
738}
739
740//________________________________________________________________
741Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
742{ // use averages for each strip if no good values exist for some single tower
743
744 // go over fTimeDepCorrection info
745 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
746
d81e6423 747 // for calculating StripAverage info
748 int nValidTower = 0;
0ce5c45d 749 Float_t stripAverage = 0;
d81e6423 750 Float_t val = 0;
751
752 int iSM = 0;
753 int iCol = 0;
754 int iRow = 0;
755 int iStrip = 0;
756 int firstCol = 0;
757 int lastCol = 0;
758
759 for (int i = 0; i < nSM; i++) {
0ce5c45d 760 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
761 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
2f17a269 762
d81e6423 763 for (int j = 0; j < nBins; j++) {
764
765 nValidTower = 0;
0ce5c45d 766 stripAverage = 0;
d81e6423 767
768 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
769 firstCol = iStrip*2;
770 if ((iSM%2)==1) { // C side
771 firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
772 }
773 lastCol = firstCol+1;
774
775 for (iCol = firstCol; iCol <= lastCol; iCol++) {
776 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
0ce5c45d 777 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
d81e6423 778 if (val>0) { // valid value; error code is negative
0ce5c45d 779 stripAverage += val;
d81e6423 780 nValidTower++;
781 }
782 }
783 }
784
785 // calc average over strip
786 if (nValidTower>0) {
0ce5c45d 787 stripAverage /= nValidTower;
d81e6423 788 for (iCol = firstCol; iCol <= lastCol; iCol++) {
789 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
0ce5c45d 790 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
d81e6423 791 if (val<0) { // invalid value; error code is negative
0ce5c45d 792 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
d81e6423 793 }
794 }
795 }
796 }
797 else { // could not fill in unvalid entries
798 nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
799 }
800
801 } // iStrip
802 } // j, bins
803 } // iSM
804
805 return nRemaining;
806}
807
808//________________________________________________________________
716fca62 809Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize)
d81e6423 810{ // OK, so we didn't have valid LED data that allowed us to do the correction only
811 // with that info.
812 // So, instead we'll rely on the temperature info and try to do the correction
813 // based on that instead.
814 // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
815 Int_t nRemaining = 0;
816
d81e6423 817 int iSM = 0;
818 int iCol = 0;
819 int iRow = 0;
820
0ce5c45d 821 Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
822 memset(dTempCoeff, 0, sizeof(dTempCoeff));
d81e6423 823 Double_t correction = 0;
716fca62 824 Double_t secondsPerBin = (Double_t) binSize;
d81e6423 825
826 for (int i = 0; i < nSM; i++) {
716fca62 827 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
0ce5c45d 828 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
2f17a269 829
0ce5c45d 830 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
ed3db319 831 AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
716fca62 832
ed3db319 833 // first get CalibTempCoeff info
d81e6423 834 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
835 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
ed3db319 836
837 dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
716fca62 838 if (fVerbosity > 1) {
839 cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
ed3db319 840 << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
716fca62 841 }
d81e6423 842 }
843 }
844
82d90a2f 845 // figure out what the reference temperature is, from fCalibReference
0ce5c45d 846 Double_t referenceTemperature = 0;
d81e6423 847 int nVal = 0;
848 for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
0ce5c45d 849 if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
850 referenceTemperature += dataCalibReference->GetTemperature(iSensor);
d81e6423 851 nVal++;
852 }
853 }
854
855 if (nVal>0) {
0ce5c45d 856 referenceTemperature /= nVal; // valid values exist, we can look into corrections
d81e6423 857
4cf7bc35 858 Double_t dSMTemperature = 0;
d81e6423 859 for (int j = 0; j < nBins; j++) {
d81e6423 860 // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
861 UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
862 // 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)
4cf7bc35 863 Double_t oldSMTemperature = dSMTemperature;
864 dSMTemperature = GetTemperatureSM(iSM, timeStamp);
865 if (j>0 && (dSMTemperature==kErrorCode)) {
866 // if we have previous values, and retrieval of values failed - use that instead (hopefully good)
867 dSMTemperature = oldSMTemperature;
868 }
869
ed3db319 870 Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
716fca62 871 if (fVerbosity > 0) {
872 cout << " referenceTemperature " << referenceTemperature
873 << " dSMTemperature " << dSMTemperature
874 << " temperatureDiff " << temperatureDiff
875 << endl;
876 }
ed3db319 877 // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down
878 // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
879 // dTempCoeff is a (unsigned) factor describing how many % the gain
880 // changes with a degree change.
716fca62 881 // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
882 // The correction we want to keep is what we should multiply our ADC value with as a function
883 // of time, i.e. the inverse of the gain change..
e1a60af4 884 if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
885 && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
ed3db319 886 // significant enough difference that we need to consider it, and also not unreasonably large
d81e6423 887
888 // loop over all towers; effect of temperature change will depend on gain for this tower
889 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
890 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
891
716fca62 892 // the correction should be inverse of modification in gain: (see discussion above)
893 // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
894 // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
ed3db319 895 correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]);
0ce5c45d 896 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
716fca62 897 if (fVerbosity > 1) {
898 cout << " iSM " << iSM
899 << " iCol " << iCol
900 << " iRow " << iRow
901 << " j " << j
902 << " correction " << correction
903 << endl;
904 }
d81e6423 905 }
906 }
907
908 } // if noteworthy temperature change
909 else { // just set correction values to 1.0
910 correction = 1;
911 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
912 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
0ce5c45d 913 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
d81e6423 914 }
915 }
916 } // else
917 } // j, Bins
918
919 } // if reference temperature exist
920 else { // could not do the needed check.. signal that in the return code
921 nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
922 }
923 } // iSM
924
925 return nRemaining;
926}
927