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 // AliEMCALCalibTempCoeff (APD temperature coefficients),
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>
37 #include "AliCDBEntry.h"
38 #include "AliCDBManager.h"
39 #include "AliEMCALSensorTempArray.h"
40 #include "AliCaloCalibSignal.h"
41 #include "AliEMCALCalibTempCoeff.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 kTempMaxDiffMedian = 2; // Temperature values should not be further away from median value within SM when considered in the average calc.
58 const double kErrorCode = -999; // to indicate that something went wrong
62 ClassImp(AliEMCALCalibTimeDep)
64 //________________________________________________________________
65 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
77 fTemperatureResolution(0.1), // 0.1 deg C is default
78 fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference
79 fTimeBinsPerHour(2), // 2 30-min bins per hour is default
80 fHighLowGainFactor(16), // factor ~16 between High gain and low gain
83 fCalibTempCoeff(NULL),
84 fCalibReference(NULL),
85 fCalibTimeDepCorrection(NULL),
91 //________________________________________________________________
92 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
94 fRun(calibt.GetRunNumber()),
95 fStartTime(calibt.GetStartTime()),
96 fEndTime(calibt.GetEndTime()),
97 fMinTemp(calibt.GetMinTemp()),
98 fMaxTemp(calibt.GetMaxTemp()),
99 fMinTempVariation(calibt.GetMinTempVariation()),
100 fMaxTempVariation(calibt.GetMaxTempVariation()),
101 fMinTempValid(calibt.GetMinTempValid()),
102 fMaxTempValid(calibt.GetMaxTempValid()),
103 fMinTime(calibt.GetMinTime()),
104 fMaxTime(calibt.GetMaxTime()),
105 fTemperatureResolution(calibt.GetTemperatureResolution()),
106 fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()),
107 fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
108 fHighLowGainFactor(calibt.GetHighLowGainFactor()),
109 fTempArray(calibt.GetTempArray()),
110 fCalibSignal(calibt.GetCalibSignal()),
111 fCalibTempCoeff(calibt.GetCalibTempCoeff()),
112 fCalibReference(calibt.GetCalibReference()),
113 fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
114 fVerbosity(calibt.GetVerbosity())
120 //________________________________________________________________
121 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
123 // assignment operator; use copy ctor
124 if (&calibt == this) return *this;
126 new (this) AliEMCALCalibTimeDep(calibt);
130 //________________________________________________________________
131 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
136 //________________________________________________________________
137 void AliEMCALCalibTimeDep::Reset()
139 // clear variables to default
145 fMinTempVariation = 0;
146 fMaxTempVariation = 0;
151 fTemperatureResolution = 0.1; // 0.1 deg C is default
152 fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference
153 fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
156 fCalibTempCoeff = NULL;
157 fCalibReference = NULL;
158 fCalibTimeDepCorrection = NULL;
163 //________________________________________________________________
164 void AliEMCALCalibTimeDep::PrintInfo() const
167 cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
168 // basic variables, all 'publicly available' also
169 cout << " VARIABLE DUMP: " << endl
170 << " GetStartTime() " << GetStartTime() << endl
171 << " GetEndTime() " << GetEndTime() << endl
172 << " GetMinTime() " << GetMinTime() << endl
173 << " GetMaxTime() " << GetMaxTime() << endl
174 << " GetMinTemp() " << GetMinTemp() << endl
175 << " GetMaxTemp() " << GetMaxTemp() << endl
176 << " GetMinTempVariation() " << GetMinTempVariation() << endl
177 << " GetMaxTempVariation() " << GetMaxTempVariation() << endl
178 << " GetTemperatureResolution() " << GetTemperatureResolution() << endl;
180 cout << " RUN INFO: " << endl
181 << " runnumber " << GetRunNumber() << endl
182 << " length (in hours) " << GetLengthOfRunInHours() << endl
183 << " length (in bins) " << GetLengthOfRunInBins() << endl
184 << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
185 << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
191 //________________________________________________________________
192 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
194 return (fEndTime - fStartTime)*kSecToHour;
197 //________________________________________________________________
198 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
200 return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
203 //________________________________________________________________
204 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
206 return (fMaxTime - fMinTime)*kSecToHour;
209 //________________________________________________________________
210 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
212 return (fMaxTemp - fMinTemp);
215 //________________________________________________________________
216 void AliEMCALCalibTimeDep::Initialize(Int_t run,
217 UInt_t startTime, UInt_t endTime)
218 { // setup, and get temperature info
219 Reset(); // start fresh
222 fStartTime = startTime;
225 // collect the needed information
226 GetTemperatureInfo(); // temperature readings during the run
227 ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
232 //________________________________________________________________
233 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
234 {// return estimate for this one SuperModule, if it had data
236 // first convert from seconds to hours..
237 Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
240 Double_t valArr[8]={0}; // 8 sensors per SM
242 for (int i=0; i<fTempArray->NumSensors(); i++) {
244 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
245 int module = st->GetSector()*2 + st->GetSide();
246 if ( module == imod ) { // right module
247 // check if we had valid data for the time that is being asked for
248 if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
249 AliSplineFit *f = st->GetFit();
250 if (f) { // ok, looks like we have valid data/info
251 // let's check what the expected value at the time appears to be
252 Double_t val = f->Eval(timeHour);
253 if ( fVerbosity > 0 ) {
254 cout << " sensor i " << i << " val " << val << endl;
256 if (val>fMinTempValid && val<fMaxTempValid && n<8) {
264 } // loop over fTempArray
266 if (n>0) { // some valid data was found
267 Double_t median = TMath::Median(n, valArr);
268 Double_t average = 0;
270 for (int is=0; is<n; is++) {
271 if (TMath::Abs(valArr[is] - median) < kTempMaxDiffMedian) {
272 average += valArr[is];
276 //cout << " n " << n << " nval " << nval << " median " << median << endl;
279 //cout << " average " << average << endl;
282 else { // this case should not happen, but kept for completeness (coverity etc)
286 else { // no good data
292 //________________________________________________________________
293 Int_t AliEMCALCalibTimeDep::CalcCorrection()
294 { // OK, this is where the real action takes place - the heart of this class..
295 /* The philosophy is as follows:
296 0. Init corrections to 1.0 values, and see how many correction bins we need
297 1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
298 2. try to use temperature info + APD temperature coefficient info, to estimate correction.
299 For now (from Dec 2009), we do not use LED info.
303 // how many SuperModules do we have?
304 Int_t nSM = fCalibReference->GetNSuperModule();
305 // how many time-bins should we have for this run?
306 Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
307 Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
309 // 1: get info on how much individual sensors might have changed during
310 // the run (compare max-min for each sensor separately)
311 if (fMaxTempVariation < fTemperatureResolution) {
312 nBins = 1; // just one bin needed..
315 binSize = fEndTime - fStartTime;
317 if (fVerbosity > 0) {
318 cout << " nBins " << nBins << " binSize " << binSize << endl;
321 // set up a reasonable default (correction = 1.0)
322 fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
323 fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
324 fCalibTimeDepCorrection->SetStartTime(fStartTime);
325 fCalibTimeDepCorrection->SetNTimeBins(nBins);
326 fCalibTimeDepCorrection->SetTimeBinSize(binSize);
328 // 2: try with Temperature correction
329 Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
335 //________________________________________________________________
336 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
337 { // estimate the Temperature Coefficient, based on the dark current (IDark)
338 // and the gain (M), based on Catania parameterizations
340 Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
341 Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
343 Double_t dTC = dP0 + dP1*M;
344 // from % numbers to regular ones..:
347 return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
350 /* Next come the methods that do the work in picking up all the needed info..*/
351 //________________________________________________________________
352 void AliEMCALCalibTimeDep::GetTemperatureInfo()
354 // pick up Preprocessor output, based on fRun (most recent version)
355 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
357 fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
361 AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
362 fTempArray->NumSensors(),
363 fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
366 AliWarning( Form("AliEMCALSensorTempArray not found!") );
372 //________________________________________________________________
373 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
374 {// assign max/min time and temperature values
376 fMinTemp = 999; // init to some large value (999 deg C)
378 fMinTempVariation = 999; // init to some large value (999 deg C)
379 fMaxTempVariation = 0;
380 fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
383 Int_t n = 0; // number of valid readings
385 for (int i=0; i<fTempArray->NumSensors(); i++) {
387 AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
388 if ( st->GetStartTime() == 0 ) { // no valid data
393 if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
394 if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
396 // check temperature ranges
397 AliSplineFit *f = st->GetFit();
399 if (f) { // ok, looks like we have valid data/info
400 int np = f->GetKnots();
401 Double_t *y0 = f->GetY0();
402 // min and max values within the single sensor
406 for (int ip=0; ip<np; ip++) {
407 if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
408 if (min > y0[ip]) { min = y0[ip]; }
409 if (max < y0[ip]) { max = y0[ip]; }
414 if (fMinTemp > min) { fMinTemp = min; }
415 if (fMaxTemp < max) { fMaxTemp = max; }
416 Double_t variation = max - min;
417 if (fMinTempVariation > variation) { fMinTempVariation = variation; }
418 if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
423 } // loop over fTempArray
425 if (n>0) { // some valid data was found
428 else { // no good data
429 return (Int_t) kErrorCode;
434 //________________________________________________________________
435 void AliEMCALCalibTimeDep::GetCalibSignalInfo()
437 // pick up Preprocessor output, based on fRun (most recent version)
438 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
440 fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
444 AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
445 fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
446 fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
447 fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
448 fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
449 fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );
452 AliWarning( Form("AliCaloCalibSignal not found!") );
458 //________________________________________________________________
459 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo()
461 // pick up Preprocessor output, based on fRun (most recent version)
462 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
463 // stored object should be a TTree; read the info
465 fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
468 if (fCalibTempCoeff) {
469 AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
472 AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
478 //________________________________________________________________
479 void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
481 // pick up Preprocessor output, based on fRun (most recent version)
482 AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
484 fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
487 if (fCalibReference) {
488 AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
491 AliWarning( Form("AliEMCALCalibReference not found!") );
497 //________________________________________________________________
498 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins)
499 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0
500 // The correction factor we keep is c(T) = R(t0)/R(T)
501 // T info from fCalibSignal, t0 info from fCalibReference
503 // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
504 // but one could upgrade this in the future
505 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
507 // sanity check; same SuperModule indices for corrections as for regular calibrations
508 for (int i = 0; i < nSM; i++) {
509 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
510 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
512 int iSMRef = dataCalibReference->GetSuperModuleNum();
513 int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
514 if (iSMRef != iSMCorr) {
515 AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
516 nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
527 // The fCalibSignal info is stored in TTrees
528 // Note that the time-bins for the TTree's may not exactly match with our correction time bins
529 int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
530 // fCalibSignal time info in seconds: Hour/kSecToHour
531 // corrected for startTime difference: Hour/kSecToHour + timeDiff
532 // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
534 // values for R(T), size of TArray = nBins
535 // the [2] dimension below is for the low or high gain
536 TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
537 TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
538 TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
539 TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
541 // set up TArray's first
542 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
543 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
544 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
545 for (iGain = 0; iGain < 2; iGain++) {
547 ampT[iSM][iCol][iRow][iGain].Set(nBins);
548 nT[iSM][iCol][iRow][iGain].Set(nBins);
550 for (int j = 0; j < nBins; j++) {
551 ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
552 nT[iSM][iCol][iRow][iGain].AddAt(0, j);
557 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
558 for (iGain = 0; iGain < 2; iGain++) {
560 ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
561 nLEDRefT[iSM][iStrip][iGain].Set(nBins);
563 for (int j = 0; j < nBins; j++) {
564 ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
565 nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
571 // OK, now loop over the TTrees and fill the arrays needed for R(T)
572 TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
573 TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
575 int iChannelNum = 0; // for regular towers
576 int iRefNum = 0; // for LED
580 treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
581 treeAvg->SetBranchAddress("fHour", &dHour);
582 treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
585 // counters for how many values were seen per SuperModule
586 int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
587 int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
589 for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
590 treeAvg->GetEntry(ient);
591 // figure out where this info comes from
592 fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
593 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
594 // add value in the arrays
595 ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
596 nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
600 treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
601 treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
602 treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
604 for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
605 treeLEDRefAvg->GetEntry(ient);
606 // figure out where this info comes from
607 fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
608 iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
609 // add value in the arrays
610 ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
611 nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
615 // Normalize TArray values, and calculate average also
616 Float_t norm = 0; // extra var, for readability
618 for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
619 if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
620 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
621 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
622 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
623 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
624 for (iGain = 0; iGain < 2; iGain++) {
626 for (int j = 0; j < nBins; j++) {
627 if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
628 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
629 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
635 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
636 for (iGain = 0; iGain < 2; iGain++) {
637 for (int j = 0; j < nBins; j++) {
638 if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
639 norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
640 ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
649 // Calculate correction values, and store them
650 // set kErrorCode values for those that could not be set
654 Float_t correction = 0; // c(T) = R(t0)/R(T)
655 Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
657 for (int i = 0; i < nSM; i++) {
658 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
659 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
660 iSM = dataCalibReference->GetSuperModuleNum();
662 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
663 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
664 iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
665 refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
667 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
670 AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
671 iGain = refVal->GetHighLow(); // gain value used for reference calibration
672 // valid amplitude values should be larger than 0
673 if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
674 ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
677 ratiot0 = kErrorCode;
681 for (int j = 0; j < nBins; j++) {
683 // calculate R(T) also; first try with individual tower:
684 // same gain as for reference calibration is the default
685 if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
686 // looks like valid data with the right gain combination
687 ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j);
689 // if data appears to be saturated, and we are in high gain, then try with low gain instead
691 int newRefGain = refGain;
692 if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
695 if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
699 if (newGain!=iGain || newRefGain!=refGain) {
700 // compensate for using different gain than in the reference calibration
701 // we may need to have a custom H/L ratio value for each tower
702 // later, but for now just use a common value, as the rest of the code does..
703 ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
706 ratioT *= fHighLowGainFactor;
708 else if (newRefGain<refGain) {
709 ratioT /= fHighLowGainFactor;
717 // Calc. correction factor
718 if (ratiot0>0 && ratioT>0) {
719 correction = ratiot0/ratioT;
722 correction = kErrorCode;
727 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
729 fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
736 nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
740 //________________________________________________________________
741 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins)
742 { // use averages for each strip if no good values exist for some single tower
744 // go over fTimeDepCorrection info
745 Int_t nRemaining = 0; // we count the towers for which we could not get valid data
747 // for calculating StripAverage info
749 Float_t stripAverage = 0;
759 for (int i = 0; i < nSM; i++) {
760 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
761 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
763 for (int j = 0; j < nBins; j++) {
768 for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
770 if ((iSM%2)==1) { // C side
771 firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
773 lastCol = firstCol+1;
775 for (iCol = firstCol; iCol <= lastCol; iCol++) {
776 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
777 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
778 if (val>0) { // valid value; error code is negative
785 // calc average over strip
787 stripAverage /= nValidTower;
788 for (iCol = firstCol; iCol <= lastCol; iCol++) {
789 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
790 val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
791 if (val<0) { // invalid value; error code is negative
792 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
797 else { // could not fill in unvalid entries
798 nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
808 //________________________________________________________________
809 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize)
810 { // OK, so we didn't have valid LED data that allowed us to do the correction only
812 // So, instead we'll rely on the temperature info and try to do the correction
813 // based on that instead.
814 // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class
815 Int_t nRemaining = 0;
821 Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
822 memset(dTempCoeff, 0, sizeof(dTempCoeff));
823 Double_t correction = 0;
824 Double_t secondsPerBin = (Double_t) binSize;
826 for (int i = 0; i < nSM; i++) {
827 AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
828 iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
830 AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
831 AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
833 // first get CalibTempCoeff info
834 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
835 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
837 dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
838 if (fVerbosity > 1) {
839 cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
840 << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
845 // figure out what the reference temperature is, from fCalibReference
846 Double_t referenceTemperature = 0;
848 for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
849 if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
850 referenceTemperature += dataCalibReference->GetTemperature(iSensor);
856 referenceTemperature /= nVal; // valid values exist, we can look into corrections
858 Double_t dSMTemperature = 0;
859 for (int j = 0; j < nBins; j++) {
860 // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
861 UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
862 // get the temperature at this time; use average over whole SM for now (TO BE CHECKED LATER - if we can do better with finer grained info)
863 Double_t oldSMTemperature = dSMTemperature;
864 dSMTemperature = GetTemperatureSM(iSM, timeStamp);
865 if (j>0 && (dSMTemperature==kErrorCode)) {
866 // if we have previous values, and retrieval of values failed - use that instead (hopefully good)
867 dSMTemperature = oldSMTemperature;
870 Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
871 if (fVerbosity > 0) {
872 cout << " referenceTemperature " << referenceTemperature
873 << " dSMTemperature " << dSMTemperature
874 << " temperatureDiff " << temperatureDiff
877 // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down
878 // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
879 // dTempCoeff is a (unsigned) factor describing how many % the gain
880 // changes with a degree change.
881 // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
882 // The correction we want to keep is what we should multiply our ADC value with as a function
883 // of time, i.e. the inverse of the gain change..
884 if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
885 && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
886 // significant enough difference that we need to consider it, and also not unreasonably large
888 // loop over all towers; effect of temperature change will depend on gain for this tower
889 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
890 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
892 // the correction should be inverse of modification in gain: (see discussion above)
893 // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01;
894 // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
895 correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]);
896 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
897 if (fVerbosity > 1) {
898 cout << " iSM " << iSM
902 << " correction " << correction
908 } // if noteworthy temperature change
909 else { // just set correction values to 1.0
911 for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
912 for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
913 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
919 } // if reference temperature exist
920 else { // could not do the needed check.. signal that in the return code
921 nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;