1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliEMCALCalibTimeDep.cxx $ */
18 //_________________________________________________________________________
20 ///////////////////////////////////////////////////////////////////////////////
22 // class for EMCAL time-dep calibration
23 // - supposed to run in preprocessor
24 // we use input from the following sources:
25 // AliEMCALBiasAPD (bias values), AliEMCALCalibMapAPD (APD calibration and location info),
26 // AliCaloCalibSignal (LED DA), AliEMCALSensorTempArray (ELMB DCS)
27 // AliEMCALCalibReference: LED amplitude and temperature info at reference time
29 // output/result is in AliEMCALCalibTimeDepCorrection
31 ///////////////////////////////////////////////////////////////////////////////
34 #include <TGraphSmooth.h>
36 #include "AliCDBEntry.h"
37 #include "AliCDBManager.h"
38 #include "AliEMCALSensorTempArray.h"
39 #include "AliCaloCalibSignal.h"
40 #include "AliEMCALBiasAPD.h"
41 #include "AliEMCALCalibMapAPD.h"
42 #include "AliEMCALCalibReference.h"
43 #include "AliEMCALCalibTimeDepCorrection.h"
44 #include "AliEMCALCalibTimeDep.h"
46 /* first a bunch of constants.. */
47 const double kSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
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
51 const double kTempCoeffP0Const = -0.903; //
52 const double kTempCoeffP0Factor = -1.381e7; //
53 const double kTempCoeffP1Const = -0.023; //
54 const double kTempCoeffP1Factor = -4.966e5; //
56 const double kErrorCode = -999; // to indicate that something went wrong
60 ClassImp(AliEMCALCalibTimeDep)
62 //________________________________________________________________
63 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
73 fTemperatureResolution(0.1), // 0.1 deg C is default
74 fTimeBinsPerHour(2), // 2 30-min bins per hour is default
75 fHighLowGainFactor(16), // factor ~16 between High gain and low gain
80 fCalibReference(NULL),
81 fCalibTimeDepCorrection(NULL),
87 //________________________________________________________________
88 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
90 fRun(calibt.GetRunNumber()),
91 fStartTime(calibt.GetStartTime()),
92 fEndTime(calibt.GetEndTime()),
93 fMinTemp(calibt.GetMinTemp()),
94 fMaxTemp(calibt.GetMaxTemp()),
95 fMinTempVariation(calibt.GetMinTempVariation()),
96 fMaxTempVariation(calibt.GetMaxTempVariation()),
97 fMinTime(calibt.GetMinTime()),
98 fMaxTime(calibt.GetMaxTime()),
99 fTemperatureResolution(calibt.GetTemperatureResolution()),
100 fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
101 fHighLowGainFactor(calibt.GetHighLowGainFactor()),
102 fTempArray(calibt.GetTempArray()),
103 fCalibSignal(calibt.GetCalibSignal()),
104 fBiasAPD(calibt.GetBiasAPD()),
105 fCalibMapAPD(calibt.GetCalibMapAPD()),
106 fCalibReference(calibt.GetCalibReference()),
107 fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
108 fVerbosity(calibt.GetVerbosity())
114 //________________________________________________________________
115 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
117 // assignment operator; use copy ctor
118 if (&calibt == this) return *this;
120 new (this) AliEMCALCalibTimeDep(calibt);
124 //________________________________________________________________
125 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
130 //________________________________________________________________
131 void AliEMCALCalibTimeDep::Reset()
133 // clear variables to default
139 fMinTempVariation = 0;
140 fMaxTempVariation = 0;
143 fTemperatureResolution = 0.1; // 0.1 deg C is default
144 fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
149 fCalibReference = NULL;
150 fCalibTimeDepCorrection = NULL;
155 //________________________________________________________________
156 void AliEMCALCalibTimeDep::PrintInfo() const
159 cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
160 // basic variables, all 'publicly available' also
161 cout << " VARIABLE DUMP: " << endl
162 << " GetStartTime() " << GetStartTime() << endl
163 << " GetEndTime() " << GetEndTime() << endl
164 << " GetMinTime() " << GetMinTime() << endl
165 << " GetMaxTime() " << GetMaxTime() << endl
166 << " GetMinTemp() " << GetMinTemp() << endl
167 << " GetMaxTemp() " << GetMaxTemp() << endl
168 << " GetMinTempVariation() " << GetMinTempVariation() << endl
169 << " GetMaxTempVariation() " << GetMaxTempVariation() << endl;
171 cout << " RUN INFO: " << endl
172 << " runnumber " << GetRunNumber() << endl
173 << " length (in hours) " << GetLengthOfRunInHours() << endl
174 << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
175 << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
181 //________________________________________________________________
182 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
184 return (fEndTime - fStartTime)*kSecToHour;
187 //________________________________________________________________
188 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
190 return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
193 //________________________________________________________________
194 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
196 return (fMaxTime - fMinTime)*kSecToHour;
199 //________________________________________________________________
200 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
202 return (fMaxTemp - fMinTemp);
205 //________________________________________________________________
206 void AliEMCALCalibTimeDep::Initialize(Int_t run,
207 UInt_t startTime, UInt_t endTime)
208 { // setup, and get temperature info
209 Reset(); // start fresh
212 fStartTime = startTime;
215 // collect the needed information
216 GetTemperatureInfo(); // temperature readings during the run
217 ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
222 //________________________________________________________________
223 Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const
224 {// return estimate for all SuperModules and sensors, that had data
226 // first convert from seconds to hours..
227 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
229 Double_t average = 0;
232 for (int i=0; i<fTempArray->NumSensors(); i++) {
234 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
236 // check if we had valid data for the time that is being asked for
237 if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
238 AliSplineFit *f = st->GetFit();
239 if (f) { // ok, looks like we have valid data/info
240 // let's check what the expected value at the time appears to be
241 Double_t val = f->Eval(timeHour);
246 } // loop over fTempArray
248 if (n>0) { // some valid data was found
252 else { // no good data
258 //________________________________________________________________
259 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
260 {// return estimate for this one SuperModule, if it had data
262 // first convert from seconds to hours..
263 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
265 Double_t average = 0;
268 for (int i=0; i<fTempArray->NumSensors(); i++) {
270 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
271 int module = st->GetSector()*2 + st->GetSide();
272 if ( module == imod ) { // right module
273 // check if we had valid data for the time that is being asked for
274 if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
275 AliSplineFit *f = st->GetFit();
276 if (f) { // ok, looks like we have valid data/info
277 // let's check what the expected value at the time appears to be
278 Double_t val = f->Eval(timeHour);
279 if ( fVerbosity > 0 ) {
280 cout << " sensor i " << i << " val " << val << endl;
288 } // loop over fTempArray
290 if (n>0) { // some valid data was found
294 else { // no good data
300 //________________________________________________________________
301 Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_t timeStamp) const
302 {// return estimate for this one SuperModule and sensor, if it had data
304 // first convert from seconds to hours..
305 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
307 for (int i=0; i<fTempArray->NumSensors(); i++) {
309 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
310 int module = st->GetSector()*2 + st->GetSide();
311 if ( module == imod && st->GetNum()==isens ) { // right module, and sensor
312 // check if we had valid data for the time that is being asked for
313 if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
314 AliSplineFit *f = st->GetFit();
315 if (f) { // ok, looks like we have valid data/info
316 // let's check what the expected value at the time appears to be
317 Double_t val = f->Eval(timeHour);
319 return val; // no point to move further in for loop, we have found the sensor we were looking for
324 } // loop over fTempArray
326 // if we made it all here, it means that we didn't find the sensor we were looking for
332 //________________________________________________________________
333 Int_t AliEMCALCalibTimeDep::CalcCorrection()
334 { // OK, this is where the real action takes place - the heart of this class..
335 /* The philosophy is as follows:
336 0. Init corrections to 1.0 values, and see how many correction bins we need
337 1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
338 2. try to use temperature info + APD bias and calibration info, to estimate correction.
339 For now (from Dec 2009), we do not use LED info, since we are not taking LED triggers during the run.
343 // how many SuperModules do we have?
344 Int_t nSM = fCalibReference->GetNSuperModule();
345 // how many time-bins should we have for this run?
346 Int_t nBins = (Int_t) GetLengthOfRunInBins(); // round-down (from double to int)
347 Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
349 // 1: get info on how much individual sensors might have changed during
350 // the run (compare max-min for each sensor separately)
351 if (fMaxTempVariation < fTemperatureResolution) {
352 nBins = 1; // just one bin needed..
353 binSize = fEndTime - fStartTime;
356 // set up a reasonable default (correction = 1.0)
357 fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
358 fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
359 fCalibTimeDepCorrection->SetStartTime(fStartTime);
360 fCalibTimeDepCorrection->SetNTimeBins(nBins);
361 fCalibTimeDepCorrection->SetTimeBinSize(binSize);
363 // 2: try with Temperature correction
364 Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
370 //________________________________________________________________
371 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
372 { // estimate the Temperature Coefficient, based on the dark current (IDark)
373 // and the gain (M), based on Catania parameterizations
375 Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
376 Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
378 Double_t dTC = dP0 + dP1*M;
383 /* Next come the methods that do the work in picking up all the needed info..*/
384 //________________________________________________________________
385 void AliEMCALCalibTimeDep::GetTemperatureInfo()
387 // pick up Preprocessor output, based on fRun (most recent version)
388 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
390 fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
394 AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
395 fTempArray->NumSensors(),
396 fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
399 AliWarning( Form("AliEMCALSensorTempArray not found!") );
405 //________________________________________________________________
406 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
407 {// assign max/min time and temperature values
409 fMinTemp = 999; // init to some large value (999 deg C)
411 fMinTempVariation = 999; // init to some large value (999 deg C)
412 fMaxTempVariation = 0;
413 fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
416 Int_t n = 0; // number of valid readings
418 for (int i=0; i<fTempArray->NumSensors(); i++) {
420 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
421 if ( st->GetStartTime() == 0 ) { // no valid data
426 if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
427 if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
429 // check temperature ranges
430 AliSplineFit *f = st->GetFit();
432 if (f) { // ok, looks like we have valid data/info
433 int np = f->GetKnots();
434 Double_t *y0 = f->GetY0();
435 // min and max values within the single sensor
438 for (int ip=0; ip<np; ip++) {
439 if (min > y0[ip]) { min = y0[ip]; }
440 if (max < y0[ip]) { max = y0[ip]; }
442 if (fMinTemp > min) { fMinTemp = min; }
443 if (fMaxTemp < max) { fMaxTemp = max; }
444 Double_t variation = max - min;
445 if (fMinTempVariation > variation) { fMinTempVariation = variation; }
446 if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
450 } // loop over fTempArray
452 if (n>0) { // some valid data was found
455 else { // no good data
456 return (Int_t) kErrorCode;
461 //________________________________________________________________
462 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
464 // pick up Preprocessor output, based on fRun (most recent version)
465 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
467 fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
471 AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
472 fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
473 fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
474 fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
475 fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
476 fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
479 AliWarning( Form("AliCaloCalibSignal not found!") );
485 //________________________________________________________________
486 void AliEMCALCalibTimeDep::GetBiasAPDInfo()
488 // pick up Preprocessor output, based on fRun (most recent version)
489 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/BiasAPD", fRun);
491 fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
495 AliInfo( Form("BiasAPD: NSuperModule %d ", fBiasAPD->GetNSuperModule() ) );
498 AliWarning( Form("AliEMCALBiasAPD not found!") );
504 //________________________________________________________________
505 void AliEMCALCalibTimeDep::GetCalibMapAPDInfo()
507 // pick up Preprocessor output, based on fRun (most recent version)
508 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
509 // stored object should be a TTree; read the info
511 fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
515 AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
518 AliWarning( Form("AliEMCALCalibMapAPD not found!") );
524 //________________________________________________________________
525 void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
527 // pick up Preprocessor output, based on fRun (most recent version)
528 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
530 fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
533 if (fCalibReference) {
534 AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
537 AliWarning( Form("AliEMCALCalibReference not found!") );
543 //________________________________________________________________
544 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
545 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
546 // The correction factor we keep is c(T) = R(t0)/R(T)
547 // T info from fCalibSignal, t0 info from fCalibReference
549 // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
550 // but one could upgrade this in the future
551 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
553 // sanity check; same SuperModule indices for corrections as for regular calibrations
554 for (int i = 0; i < nSM; i++) {
555 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
556 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
558 int iSMRef = dataCalibReference->GetSuperModuleNum();
559 int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
560 if (iSMRef != iSMCorr) {
561 AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
562 nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
573 // The fCalibSignal info is stored in TTrees
574 // Note that the time-bins for the TTree's may not exactly match with our correction time bins
575 int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
576 // fCalibSignal time info in seconds: Hour/kSecToHour
577 // corrected for startTime difference: Hour/kSecToHour + timeDiff
578 // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
580 // values for R(T), size of TArray = nBins
581 // the [2] dimension below is for the low or high gain
582 TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
583 TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
584 TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
585 TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
587 // set up TArray's first
588 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
589 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
590 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
591 for (iGain = 0; iGain < 2; iGain++) {
593 ampT[iSM][iCol][iRow][iGain].Set(nBins);
594 nT[iSM][iCol][iRow][iGain].Set(nBins);
596 for (int j = 0; j < nBins; j++) {
597 ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
598 nT[iSM][iCol][iRow][iGain].AddAt(0, j);
603 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
604 for (iGain = 0; iGain < 2; iGain++) {
606 ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
607 nLEDRefT[iSM][iStrip][iGain].Set(nBins);
609 for (int j = 0; j < nBins; j++) {
610 ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
611 nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
617 // OK, now loop over the TTrees and fill the arrays needed for R(T)
618 TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
619 TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
621 int iChannelNum = 0; // for regular towers
622 int iRefNum = 0; // for LED
626 treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
627 treeAvg->SetBranchAddress("fHour", &dHour);
628 treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
631 // counters for how many values were seen per SuperModule
632 int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
633 int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
635 for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
636 treeAvg->GetEntry(ient);
637 // figure out where this info comes from
638 fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
639 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
640 // add value in the arrays
641 ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
642 nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
646 treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
647 treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
648 treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
650 for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
651 treeLEDRefAvg->GetEntry(ient);
652 // figure out where this info comes from
653 fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
654 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
655 // add value in the arrays
656 ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
657 nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
661 // Normalize TArray values, and calculate average also
662 Float_t norm = 0; // extra var, for readability
664 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
665 if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
666 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
667 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
668 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
669 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
670 for (iGain = 0; iGain < 2; iGain++) {
672 for (int j = 0; j < nBins; j++) {
673 if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
674 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
675 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
681 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
682 for (iGain = 0; iGain < 2; iGain++) {
683 for (int j = 0; j < nBins; j++) {
684 if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
685 norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
686 ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
695 // Calculate correction values, and store them
696 // set kErrorCode values for those that could not be set
700 Float_t correction = 0; // c(T) = R(t0)/R(T)
701 Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
703 for (int i = 0; i < nSM; i++) {
704 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
705 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
706 iSM = dataCalibReference->GetSuperModuleNum();
708 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
709 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
710 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
711 refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
713 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
716 AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
717 iGain = refVal->GetHighLow(); // gain value used for reference calibration
718 // valid amplitude values should be larger than 0
719 if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
720 ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
723 ratiot0 = kErrorCode;
727 for (int j = 0; j < nBins; j++) {
729 // calculate R(T) also; first try with individual tower:
730 // same gain as for reference calibration is the default
731 if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
732 // looks like valid data with the right gain combination
733 ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
735 // if data appears to be saturated, and we are in high gain, then try with low gain instead
737 int newRefGain = refGain;
738 if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
741 if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
745 if (newGain!=iGain || newRefGain!=refGain) {
746 // compensate for using different gain than in the reference calibration
747 // we may need to have a custom H/L ratio value for each tower
748 // later, but for now just use a common value, as the rest of the code does..
749 ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
752 ratioT *= fHighLowGainFactor;
754 else if (newRefGain<refGain) {
755 ratioT /= fHighLowGainFactor;
763 // Calc. correction factor
764 if (ratiot0>0 && ratioT>0) {
765 correction = ratiot0/ratioT;
768 correction = kErrorCode;
773 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
775 fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
782 nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
786 //________________________________________________________________
787 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
788 { // use averages for each strip if no good values exist for some single tower
790 // go over fTimeDepCorrection info
791 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
793 // for calculating StripAverage info
795 Float_t stripAverage = 0;
805 for (int i = 0; i < nSM; i++) {
806 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
807 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
809 for (int j = 0; j < nBins; j++) {
814 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
816 if ((iSM%2)==1) { // C side
817 firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
819 lastCol = firstCol+1;
821 for (iCol = firstCol; iCol <= lastCol; iCol++) {
822 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
823 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
824 if (val>0) { // valid value; error code is negative
831 // calc average over strip
833 stripAverage /= nValidTower;
834 for (iCol = firstCol; iCol <= lastCol; iCol++) {
835 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
836 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
837 if (val<0) { // invalid value; error code is negative
838 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
843 else { // could not fill in unvalid entries
844 nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
854 //________________________________________________________________
855 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize)
856 { // OK, so we didn't have valid LED data that allowed us to do the correction only
858 // So, instead we'll rely on the temperature info and try to do the correction
859 // based on that instead.
860 // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
861 Int_t nRemaining = 0;
867 Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
868 memset(dTempCoeff, 0, sizeof(dTempCoeff));
870 Double_t correction = 0;
871 Double_t secondsPerBin = (Double_t) binSize;
873 for (int i = 0; i < nSM; i++) {
874 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
875 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
877 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
878 AliEMCALSuperModuleCalibMapAPD * dataCalibMapAPD = fCalibMapAPD->GetSuperModuleCalibMapAPDNum(iSM);
879 AliEMCALSuperModuleBiasAPD * dataBiasAPD = fBiasAPD->GetSuperModuleBiasAPDNum(iSM);
881 // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
882 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
883 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
884 AliEMCALCalibMapAPDVal * mapAPD = dataCalibMapAPD->GetAPDVal(iCol, iRow);
885 gainM = fCalibMapAPD->GetGain(mapAPD->GetPar(0), mapAPD->GetPar(1), mapAPD->GetPar(2),
886 dataBiasAPD->GetVoltage(iCol, iRow));
887 dTempCoeff[iCol][iRow] = GetTempCoeff(mapAPD->GetDarkCurrent(), gainM);
888 if (fVerbosity > 1) {
889 cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
890 << " par0 " << mapAPD->GetPar(0)
891 << " par1 " << mapAPD->GetPar(1)
892 << " par2 " << mapAPD->GetPar(2)
893 << " bias " << dataBiasAPD->GetVoltage(iCol, iRow)
894 << " gainM " << gainM << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
899 // figure out what the reference temperature is, from fCalibReference
900 Double_t referenceTemperature = 0;
902 for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
903 if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
904 referenceTemperature += dataCalibReference->GetTemperature(iSensor);
910 referenceTemperature /= nVal; // valid values exist, we can look into corrections
912 for (int j = 0; j < nBins; j++) {
913 // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
914 UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
915 // 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)
916 Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp);
918 Double_t temperatureDiff = dSMTemperature - referenceTemperature ; // new - old
919 if (fVerbosity > 0) {
920 cout << " referenceTemperature " << referenceTemperature
921 << " dSMTemperature " << dSMTemperature
922 << " temperatureDiff " << temperatureDiff
925 // if the new temperature is higher than the old/reference one, then the gain has gone down
926 // if the new temperature is lower than the old/reference one, then the gain has gone up
927 // dTempCoeff is a negative number describing how many % (hence 0.01 factor below) the gain
928 // changes with a positive degree change.
929 // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
930 // The correction we want to keep is what we should multiply our ADC value with as a function
931 // of time, i.e. the inverse of the gain change..
932 if (fabs(temperatureDiff)>fTemperatureResolution) {
933 // significant enough difference that we need to consider it
935 // loop over all towers; effect of temperature change will depend on gain for this tower
936 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
937 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
939 // the correction should be inverse of modification in gain: (see discussion above)
940 // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
941 // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
942 correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
943 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
944 if (fVerbosity > 1) {
945 cout << " iSM " << iSM
949 << " correction " << correction
955 } // if noteworthy temperature change
956 else { // just set correction values to 1.0
958 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
959 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
960 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
966 } // if reference temperature exist
967 else { // could not do the needed check.. signal that in the return code
968 nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;