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