]> 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
Updates needed for shuttle preprocessor
[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) {
444 fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
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);
462 if (entry) {
463 fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
464 }
465
466 if (fCalibMapAPD) {
467 AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
468 }
469 else {
470 AliWarning( Form("AliEMCALCalibMapAPD not found!") );
471 }
472
473 return;
474}
475
476//________________________________________________________________
477void AliEMCALCalibTimeDep::GetCalibAbsInfo()
478{
479 // pick up Preprocessor output, based on fRun (most recent version)
480 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
481 if (entry) {
482 fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject();
483 }
484
485 if (fCalibAbs) {
486 AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
487 }
488 else {
489 AliWarning( Form("AliEMCALCalibAbs not found!") );
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)
499 // T info from fCalibSignal, t0 info from fCalibAbs
500
501 // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
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
506 AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
507 AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
508
509 for (int i = 0; i < nSM; i++) {
510 int iSMAbs = CalibAbsData[i].fSuperModuleNum;
511 int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
512 if (iSMAbs != iSMCorr) {
513 AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
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
528 // fCalibSignal time info in seconds: Hour/fkSecToHour
529 // corrected for startTime difference: Hour/fkSecToHour + timeDiff
530 // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
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)
570 TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
571 TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
572
573 int ChannelNum; // for regular towers
574 int RefNum; // for LED
575 double Hour;
576 double AvgAmp;
577
578 TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
579 TAvg->SetBranchAddress("fHour", &Hour);
580 TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
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
587 for (int ient=0; ient<TAvg->GetEntries(); ient++) {
588 TAvg->GetEntry(ient);
589 // figure out where this info comes from
590 fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
591 iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour ); // CHECK!!!
592 // add value in the arrays
593 ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
594 nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
595 nCount[iSM]++;
596 }
597
598 TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
599 TLEDRefAvg->SetBranchAddress("fHour", &Hour);
600 TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
601
602 for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
603 TLEDRefAvg->GetEntry(ient);
604 // figure out where this info comes from
605 fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
606 iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour ); // CHECK!!!
607 // add value in the arrays
608 ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
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
648 // set fkErrorCode values for those that could not be set
649
650 Float_t Rt0 = 0;
651 Float_t RT = 0;
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++) {
656 iSM = CalibAbsData[i].fSuperModuleNum;
657 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
658 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
659 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
660 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
661
662 // Calc. R(t0):
663 AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
664 iGain = v.fHighLow; // gain value used for abs. calibration
665 refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration
666
667 // valid amplitude values should be larger than 0
668 if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
669 Rt0 = v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
670 }
671 else {
672 Rt0 = fkErrorCode;
673 }
674
675 // Cal R(T)
676 for (int j = 0; j < nBins; j++) {
677
678 // calculate R(T) also; first try with individual tower:
679 // same gain as for abs. calibration is the default
680 if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
681 // looks like valid data with the right gain combination
682 RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
683
684 // if data appears to be saturated, and we are in high gain, then try with low gain instead
685 if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
686 (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
687 RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j);
688 RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
689 }
690 }
691 else {
692 RT = fkErrorCode;
693 }
694
695 // Calc. correction factor
696 if (Rt0>0 && RT>0) {
697 correction = Rt0/RT;
698 }
699 else {
700 correction = fkErrorCode;
701 nRemaining++;
702 }
703
704 // Store the value
705 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
706 /* Check that
707 fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
708 also works OK */
709 } // nBins
710 }
711 }
712 }
713
714 nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
715 return nRemaining;
716}
717
718//________________________________________________________________
719Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
720{ // use averages for each strip if no good values exist for some single tower
721
722 // go over fTimeDepCorrection info
723 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
724
725 AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
726
727 // for calculating StripAverage info
728 int nValidTower = 0;
729 Float_t StripAverage = 0;
730 Float_t val = 0;
731
732 int iSM = 0;
733 int iCol = 0;
734 int iRow = 0;
735 int iStrip = 0;
736 int firstCol = 0;
737 int lastCol = 0;
738
739 for (int i = 0; i < nSM; i++) {
740 iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
741 for (int j = 0; j < nBins; j++) {
742
743 nValidTower = 0;
744 StripAverage = 0;
745
746 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
747 firstCol = iStrip*2;
748 if ((iSM%2)==1) { // C side
749 firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
750 }
751 lastCol = firstCol+1;
752
753 for (iCol = firstCol; iCol <= lastCol; iCol++) {
754 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
755 val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
756 if (val>0) { // valid value; error code is negative
757 StripAverage += val;
758 nValidTower++;
759 }
760 }
761 }
762
763 // calc average over strip
764 if (nValidTower>0) {
765 StripAverage /= nValidTower;
766 for (iCol = firstCol; iCol <= lastCol; iCol++) {
767 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
768 val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
769 if (val<0) { // invalid value; error code is negative
770 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
771 }
772 }
773 }
774 }
775 else { // could not fill in unvalid entries
776 nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
777 }
778
779 } // iStrip
780 } // j, bins
781 } // iSM
782
783 return nRemaining;
784}
785
786//________________________________________________________________
787Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins)
788{ // OK, so we didn't have valid LED data that allowed us to do the correction only
789 // with that info.
790 // So, instead we'll rely on the temperature info and try to do the correction
791 // based on that instead.
792 // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
793 Int_t nRemaining = 0;
794
795 // info containers
796 AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
797 AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
798 AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
799 // correction container
800 AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
801
802 int iSM = 0;
803 int iCol = 0;
804 int iRow = 0;
805
806 Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
807 memset(TempCoeff, 0, sizeof(TempCoeff));
808 Float_t MGain = 0;
809 Double_t correction = 0;
220ed45a 810 Double_t secondsPerBin = (3600/fTimeBinsPerHour);
d81e6423 811
812 for (int i = 0; i < nSM; i++) {
813 iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
814
815 // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
816 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
817 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
818 AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
819 MGain = fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2],
820 BiasAPDData[i].fVoltage[iCol][iRow]);
821 TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
822 }
823 }
824
825 // figure out what the reference temperature is, from fCalibAbs
826 Double_t ReferenceTemperature = 0;
827 int nVal = 0;
828 for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
829 if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
830 ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
831 nVal++;
832 }
833 }
834
835 if (nVal>0) {
836 ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
837
838 for (int j = 0; j < nBins; j++) {
839
840 // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
841 UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
842 // 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)
843 Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp);
844
845 Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
846 // if the new temperature is higher than the old/reference one, then the gain has gone down
847 if (fabs(TemperatureDiff)>fTemperatureResolution) {
848 // significant enough difference that we need to consider it
849
850 // loop over all towers; effect of temperature change will depend on gain for this tower
851 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
852 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
853
854 correction = TemperatureDiff * TempCoeff[iCol][iRow];
855 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
856 }
857 }
858
859 } // if noteworthy temperature change
860 else { // just set correction values to 1.0
861 correction = 1;
862 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
863 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
864 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
865 }
866 }
867 } // else
868 } // j, Bins
869
870 } // if reference temperature exist
871 else { // could not do the needed check.. signal that in the return code
872 nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
873 }
874 } // iSM
875
876 return nRemaining;
877}
878