X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALCalibTimeDep.cxx;h=39d0d779eec0c1cc1bf4174631ffe1fedc0a49f9;hb=03358a32f3a81ff04986da7e9ad70a2d3d8ae208;hp=33895968710b7e770093cb775c2e4a1cbde6e965;hpb=220ed45a1b69a4c518291b7428c8821e7a724737;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALCalibTimeDep.cxx b/EMCAL/AliEMCALCalibTimeDep.cxx index 33895968710..39d0d779eec 100644 --- a/EMCAL/AliEMCALCalibTimeDep.cxx +++ b/EMCAL/AliEMCALCalibTimeDep.cxx @@ -19,7 +19,14 @@ ///*-- Author: /////////////////////////////////////////////////////////////////////////////// // // -// class for EMCAL time-dep calibration // +// class for EMCAL time-dep calibration +// - supposed to run in preprocessor +// we use input from the following sources: +// AliEMCALBiasAPD (bias values), AliEMCALCalibMapAPD (APD calibration and location info), +// AliCaloCalibSignal (LED DA), AliEMCALSensorTempArray (ELMB DCS) +// AliEMCALCalibReference: LED amplitude and temperature info at reference time +// +// output/result is in AliEMCALCalibTimeDepCorrection // // /////////////////////////////////////////////////////////////////////////////// @@ -32,21 +39,21 @@ #include "AliCaloCalibSignal.h" #include "AliEMCALBiasAPD.h" #include "AliEMCALCalibMapAPD.h" -#include "AliEMCALCalibAbs.h" +#include "AliEMCALCalibReference.h" #include "AliEMCALCalibTimeDepCorrection.h" #include "AliEMCALCalibTimeDep.h" /* first a bunch of constants.. */ -const double fkSecToHour = 1.0/3600.0; // conversion factor from seconds to hours +const double kSecToHour = 1.0/3600.0; // conversion factor from seconds to hours // some global variables for APD handling; values from Catania studies, best fit // TempCoeff = p0+p1*M (M=gain), where p0 and and p1 are functions of the dark current -const double fkTempCoeffP0Const = -0.903; // -const double fkTempCoeffP0Factor = -1.381e7; // -const double fkTempCoeffP1Const = -0.023; // -const double fkTempCoeffP1Factor = -4.966e5; // +const double kTempCoeffP0Const = -0.903; // +const double kTempCoeffP0Factor = -1.381e7; // +const double kTempCoeffP1Const = -0.023; // +const double kTempCoeffP1Factor = -4.966e5; // -const double fkErrorCode = -999; // to indicate that something went wrong +const double kErrorCode = -999; // to indicate that something went wrong using namespace std; @@ -59,16 +66,20 @@ AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() : fEndTime(0), fMinTemp(0), fMaxTemp(0), + fMinTempVariation(0), + fMaxTempVariation(0), fMinTime(0), fMaxTime(0), fTemperatureResolution(0.1), // 0.1 deg C is default fTimeBinsPerHour(2), // 2 30-min bins per hour is default + fHighLowGainFactor(16), // factor ~16 between High gain and low gain fTempArray(NULL), fCalibSignal(NULL), fBiasAPD(NULL), fCalibMapAPD(NULL), - fCalibAbs(NULL), - fCalibTimeDepCorrection(NULL) + fCalibReference(NULL), + fCalibTimeDepCorrection(NULL), + fVerbosity(0) { // Constructor } @@ -81,16 +92,20 @@ AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) : fEndTime(calibt.GetEndTime()), fMinTemp(calibt.GetMinTemp()), fMaxTemp(calibt.GetMaxTemp()), + fMinTempVariation(calibt.GetMinTempVariation()), + fMaxTempVariation(calibt.GetMaxTempVariation()), fMinTime(calibt.GetMinTime()), fMaxTime(calibt.GetMaxTime()), fTemperatureResolution(calibt.GetTemperatureResolution()), fTimeBinsPerHour(calibt.GetTimeBinsPerHour()), + fHighLowGainFactor(calibt.GetHighLowGainFactor()), fTempArray(calibt.GetTempArray()), fCalibSignal(calibt.GetCalibSignal()), fBiasAPD(calibt.GetBiasAPD()), fCalibMapAPD(calibt.GetCalibMapAPD()), - fCalibAbs(calibt.GetCalibAbs()), - fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()) + fCalibReference(calibt.GetCalibReference()), + fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()), + fVerbosity(calibt.GetVerbosity()) { // copy constructor } @@ -121,6 +136,8 @@ void AliEMCALCalibTimeDep::Reset() fEndTime = 0; fMinTemp = 0; fMaxTemp = 0; + fMinTempVariation = 0; + fMaxTempVariation = 0; fMinTime = 0; fMaxTime = 0; fTemperatureResolution = 0.1; // 0.1 deg C is default @@ -129,8 +146,9 @@ void AliEMCALCalibTimeDep::Reset() fCalibSignal = NULL; fBiasAPD = NULL; fCalibMapAPD = NULL; - fCalibAbs = NULL; + fCalibReference = NULL; fCalibTimeDepCorrection = NULL; + fVerbosity = 0; return; } @@ -143,10 +161,15 @@ void AliEMCALCalibTimeDep::PrintInfo() const cout << " VARIABLE DUMP: " << endl << " GetStartTime() " << GetStartTime() << endl << " GetEndTime() " << GetEndTime() << endl + << " GetMinTime() " << GetMinTime() << endl + << " GetMaxTime() " << GetMaxTime() << endl << " GetMinTemp() " << GetMinTemp() << endl - << " GetMaxTemp() " << GetMaxTemp() << endl; + << " GetMaxTemp() " << GetMaxTemp() << endl + << " GetMinTempVariation() " << GetMinTempVariation() << endl + << " GetMaxTempVariation() " << GetMaxTempVariation() << endl; // run ranges cout << " RUN INFO: " << endl + << " runnumber " << GetRunNumber() << endl << " length (in hours) " << GetLengthOfRunInHours() << endl << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours() << " (in deg. C) " << GetRangeOfTempMeasureInDegrees() @@ -158,19 +181,19 @@ void AliEMCALCalibTimeDep::PrintInfo() const //________________________________________________________________ Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const { - return (fEndTime - fStartTime)*fkSecToHour; + return (fEndTime - fStartTime)*kSecToHour; } //________________________________________________________________ Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const { - return (fEndTime - fStartTime)*fkSecToHour*fTimeBinsPerHour; + return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour; } //________________________________________________________________ Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const { - return (fMaxTime - fMinTime)*fkSecToHour; + return (fMaxTime - fMinTime)*kSecToHour; } //________________________________________________________________ @@ -182,7 +205,7 @@ Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const //________________________________________________________________ void AliEMCALCalibTimeDep::Initialize(Int_t run, UInt_t startTime, UInt_t endTime) -{ +{ // setup, and get temperature info Reset(); // start fresh fRun = run; @@ -201,7 +224,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const {// return estimate for all SuperModules and sensors, that had data // first convert from seconds to hours.. - Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour; + Double_t timeHour = (timeStamp - fStartTime) * kSecToHour; Double_t average = 0; int n = 0; @@ -227,7 +250,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const return average; } else { // no good data - return fkErrorCode; + return kErrorCode; } } @@ -237,7 +260,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) cons {// return estimate for this one SuperModule, if it had data // first convert from seconds to hours.. - Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour; + Double_t timeHour = (timeStamp - fStartTime) * kSecToHour; Double_t average = 0; int n = 0; @@ -253,7 +276,9 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) cons if (f) { // ok, looks like we have valid data/info // let's check what the expected value at the time appears to be Double_t val = f->Eval(timeHour); - cout << " i " << i << " val " << val << endl; + if ( fVerbosity > 0 ) { + cout << " sensor i " << i << " val " << val << endl; + } average += val; n++; } @@ -267,7 +292,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) cons return average; } else { // no good data - return fkErrorCode; + return kErrorCode; } } @@ -277,7 +302,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_ {// return estimate for this one SuperModule and sensor, if it had data // first convert from seconds to hours.. - Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour; + Double_t timeHour = (timeStamp - fStartTime) * kSecToHour; for (int i=0; iNumSensors(); i++) { @@ -300,7 +325,7 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_ // if we made it all here, it means that we didn't find the sensor we were looking for // i.e. no good data - return fkErrorCode; + return kErrorCode; } @@ -308,31 +333,35 @@ Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_ Int_t AliEMCALCalibTimeDep::CalcCorrection() { // OK, this is where the real action takes place - the heart of this class.. /* The philosophy is as follows: - 0. Init corrections to 1.0 values - 1: if we have LED info for the tower, use it - 2. if not 1, we rely on LED info averaged over strip - 3. if not 2 either, we try to use temperature info + APD bias and calibration info + 0. Init corrections to 1.0 values, and see how many correction bins we need + 1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin) + 2. try to use temperature info + APD bias and calibration info, to estimate correction. + For now (from Dec 2009), we do not use LED info, since we are not taking LED triggers during the run. */ // 0: Init // how many SuperModules do we have? - Int_t nSM = fCalibAbs->GetNSuperModule(); + Int_t nSM = fCalibReference->GetNSuperModule(); // how many time-bins should we have for this run? Int_t nBins = (Int_t) GetLengthOfRunInBins(); // round-down (from double to int) Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds + + // 1: get info on how much individual sensors might have changed during + // the run (compare max-min for each sensor separately) + if (fMaxTempVariation < fTemperatureResolution) { + nBins = 1; // just one bin needed.. + binSize = fEndTime - fStartTime; + } + // set up a reasonable default (correction = 1.0) + fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM); fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0); fCalibTimeDepCorrection->SetStartTime(fStartTime); fCalibTimeDepCorrection->SetNTimeBins(nBins); fCalibTimeDepCorrection->SetTimeBinSize(binSize); - // 1+2: try with LED corrections - Int_t nRemaining = CalcLEDCorrection(nSM, nBins); - - // 3: try with Temperature, if needed - if (nRemaining>0) { - nRemaining = CalcTemperatureCorrection(nSM, nBins); - } + // 2: try with Temperature correction + Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize); return nRemaining; } @@ -343,12 +372,12 @@ Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const { // estimate the Temperature Coefficient, based on the dark current (IDark) // and the gain (M), based on Catania parameterizations - Double_t P0 = fkTempCoeffP0Const + fkTempCoeffP0Factor * IDark; - Double_t P1 = fkTempCoeffP1Const + fkTempCoeffP1Factor * IDark; + Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark; + Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark; - Double_t TC = P0 + P1*M; + Double_t dTC = dP0 + dP1*M; - return TC; + return dTC; } /* Next come the methods that do the work in picking up all the needed info..*/ @@ -379,6 +408,8 @@ Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() fMinTemp = 999; // init to some large value (999 deg C) fMaxTemp = 0; + fMinTempVariation = 999; // init to some large value (999 deg C) + fMaxTempVariation = 0; fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times.. fMaxTime = 0; @@ -387,17 +418,33 @@ Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() for (int i=0; iNumSensors(); i++) { AliEMCALSensorTemp *st = fTempArray->GetSensor(i); + if ( st->GetStartTime() == 0 ) { // no valid data + continue; + } // check time ranges if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); } if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); } - + // check temperature ranges - TGraph *g = st->GetGraph(); - if (g) { // ok, looks like we have valid data/info - // let's check what the expected value at the time appears to be - if (fMinTemp > g->GetMinimum()) { fMinTemp = g->GetMinimum(); } - if (fMaxTemp < g->GetMaximum()) { fMaxTemp = g->GetMaximum(); } + AliSplineFit *f = st->GetFit(); + + if (f) { // ok, looks like we have valid data/info + int np = f->GetKnots(); + Double_t *y0 = f->GetY0(); + // min and max values within the single sensor + Double_t min = 999; + Double_t max = 0; + for (int ip=0; ip y0[ip]) { min = y0[ip]; } + if (max < y0[ip]) { max = y0[ip]; } + } + if (fMinTemp > min) { fMinTemp = min; } + if (fMaxTemp < max) { fMaxTemp = max; } + Double_t variation = max - min; + if (fMinTempVariation > variation) { fMinTempVariation = variation; } + if (fMaxTempVariation < variation) { fMaxTempVariation = variation; } + n++; } } // loop over fTempArray @@ -406,7 +453,7 @@ Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() return n; } else { // no good data - return (Int_t) fkErrorCode; + return (Int_t) kErrorCode; } } @@ -459,6 +506,7 @@ void AliEMCALCalibTimeDep::GetCalibMapAPDInfo() { // pick up Preprocessor output, based on fRun (most recent version) AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun); + // stored object should be a TTree; read the info if (entry) { fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject(); } @@ -474,19 +522,19 @@ void AliEMCALCalibTimeDep::GetCalibMapAPDInfo() } //________________________________________________________________ -void AliEMCALCalibTimeDep::GetCalibAbsInfo() +void AliEMCALCalibTimeDep::GetCalibReferenceInfo() { // pick up Preprocessor output, based on fRun (most recent version) - AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun); + AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun); if (entry) { - fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject(); + fCalibReference = (AliEMCALCalibReference *) entry->GetObject(); } - if (fCalibAbs) { - AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) ); + if (fCalibReference) { + AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) ); } else { - AliWarning( Form("AliEMCALCalibAbs not found!") ); + AliWarning( Form("AliEMCALCalibReference not found!") ); } return; @@ -496,21 +544,21 @@ void AliEMCALCalibTimeDep::GetCalibAbsInfo() Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0 // The correction factor we keep is c(T) = R(t0)/R(T) - // T info from fCalibSignal, t0 info from fCalibAbs + // T info from fCalibSignal, t0 info from fCalibReference - // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs + // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference // but one could upgrade this in the future Int_t nRemaining = 0; // we count the towers for which we could not get valid data // sanity check; same SuperModule indices for corrections as for regular calibrations - AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData(); - AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData(); - for (int i = 0; i < nSM; i++) { - int iSMAbs = CalibAbsData[i].fSuperModuleNum; - int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum; - if (iSMAbs != iSMCorr) { - AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) ); + AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i); + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + + int iSMRef = dataCalibReference->GetSuperModuleNum(); + int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum(); + if (iSMRef != iSMCorr) { + AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) ); nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins; return nRemaining; } @@ -525,9 +573,9 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) // The fCalibSignal info is stored in TTrees // Note that the time-bins for the TTree's may not exactly match with our correction time bins int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds - // fCalibSignal time info in seconds: Hour/fkSecToHour - // corrected for startTime difference: Hour/fkSecToHour + timeDiff - // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour + // fCalibSignal time info in seconds: Hour/kSecToHour + // corrected for startTime difference: Hour/kSecToHour + timeDiff + // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour // values for R(T), size of TArray = nBins // the [2] dimension below is for the low or high gain @@ -567,45 +615,45 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) } // OK, now loop over the TTrees and fill the arrays needed for R(T) - TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime(); - TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime(); + TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime(); + TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime(); - int ChannelNum; // for regular towers - int RefNum; // for LED - double Hour; - double AvgAmp; + int iChannelNum; // for regular towers + int iRefNum; // for LED + double dHour; + double dAvgAmp; - TAvg->SetBranchAddress("fChannelNum", &ChannelNum); - TAvg->SetBranchAddress("fHour", &Hour); - TAvg->SetBranchAddress("fAvgAmp",&AvgAmp); + treeAvg->SetBranchAddress("fChannelNum", &iChannelNum); + treeAvg->SetBranchAddress("fHour", &dHour); + treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp); int iBin = 0; // counters for how many values were seen per SuperModule int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0}; int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0}; - for (int ient=0; ientGetEntries(); ient++) { - TAvg->GetEntry(ient); + for (int ient=0; ientGetEntries(); ient++) { + treeAvg->GetEntry(ient); // figure out where this info comes from - fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain); - iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour ); // CHECK!!! + fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain); + iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!! // add value in the arrays - ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin ); + ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin ); nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin ); nCount[iSM]++; } - TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum); - TLEDRefAvg->SetBranchAddress("fHour", &Hour); - TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp); + treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum); + treeLEDRefAvg->SetBranchAddress("fHour", &dHour); + treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp); - for (int ient=0; ientGetEntries(); ient++) { - TLEDRefAvg->GetEntry(ient); + for (int ient=0; ientGetEntries(); ient++) { + treeLEDRefAvg->GetEntry(ient); // figure out where this info comes from - fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain); - iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour ); // CHECK!!! + fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain); + iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!! // add value in the arrays - ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin ); + ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin ); nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin ); nCountLEDRef[iSM]++; } @@ -645,64 +693,84 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) // Calculate correction values, and store them - // set fkErrorCode values for those that could not be set + // set kErrorCode values for those that could not be set - Float_t Rt0 = 0; - Float_t RT = 0; + Float_t ratiot0 = 0; + Float_t ratioT = 0; Float_t correction = 0; // c(T) = R(t0)/R(T) Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation) for (int i = 0; i < nSM; i++) { - iSM = CalibAbsData[i].fSuperModuleNum; + AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i); + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + iSM = dataCalibReference->GetSuperModuleNum(); + for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol); iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME + refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { // Calc. R(t0): - AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow]; - iGain = v.fHighLow; // gain value used for abs. calibration - refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration - + AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow); + iGain = refVal->GetHighLow(); // gain value used for reference calibration // valid amplitude values should be larger than 0 - if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) { - Rt0 = v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip]; + if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) { + ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip); } else { - Rt0 = fkErrorCode; + ratiot0 = kErrorCode; } - // Cal R(T) + // Calc. R(T) for (int j = 0; j < nBins; j++) { // calculate R(T) also; first try with individual tower: - // same gain as for abs. calibration is the default + // same gain as for reference calibration is the default if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) { // looks like valid data with the right gain combination - RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); + ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); // if data appears to be saturated, and we are in high gain, then try with low gain instead - if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) || - (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 - RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j); - RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration + int newGain = iGain; + int newRefGain = refGain; + if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) { + newGain = 0; + } + if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) { + newRefGain = 0; + } + + if (newGain!=iGain || newRefGain!=refGain) { + // compensate for using different gain than in the reference calibration + // we may need to have a custom H/L ratio value for each tower + // later, but for now just use a common value, as the rest of the code does.. + ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j); + + if (newGain0 && RT>0) { - correction = Rt0/RT; + if (ratiot0>0 && ratioT>0) { + correction = ratiot0/ratioT; } else { - correction = fkErrorCode; + correction = kErrorCode; nRemaining++; } // Store the value - CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j); + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j); /* Check that fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction); also works OK */ @@ -722,11 +790,9 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) // go over fTimeDepCorrection info Int_t nRemaining = 0; // we count the towers for which we could not get valid data - AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData(); - // for calculating StripAverage info int nValidTower = 0; - Float_t StripAverage = 0; + Float_t stripAverage = 0; Float_t val = 0; int iSM = 0; @@ -737,11 +803,13 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) int lastCol = 0; for (int i = 0; i < nSM; i++) { - iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum; + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + iSM = dataCalibTimeDepCorrection->GetSuperModuleNum(); + for (int j = 0; j < nBins; j++) { nValidTower = 0; - StripAverage = 0; + stripAverage = 0; for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) { firstCol = iStrip*2; @@ -752,9 +820,9 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) for (iCol = firstCol; iCol <= lastCol; iCol++) { for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { - val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j); + val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j); if (val>0) { // valid value; error code is negative - StripAverage += val; + stripAverage += val; nValidTower++; } } @@ -762,12 +830,12 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) // calc average over strip if (nValidTower>0) { - StripAverage /= nValidTower; + stripAverage /= nValidTower; for (iCol = firstCol; iCol <= lastCol; iCol++) { for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { - val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j); + val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j); if (val<0) { // invalid value; error code is negative - CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j); + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j); } } } @@ -784,7 +852,7 @@ Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) } //________________________________________________________________ -Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins) +Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize) { // OK, so we didn't have valid LED data that allowed us to do the correction only // with that info. // So, instead we'll rely on the temperature info and try to do the correction @@ -792,67 +860,95 @@ Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins) // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class Int_t nRemaining = 0; - // info containers - AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData(); - AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData(); - AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData(); - // correction container - AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData(); - int iSM = 0; int iCol = 0; int iRow = 0; - Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows]; - memset(TempCoeff, 0, sizeof(TempCoeff)); - Float_t MGain = 0; + Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows]; + memset(dTempCoeff, 0, sizeof(dTempCoeff)); + Float_t gainM = 0; Double_t correction = 0; - Double_t secondsPerBin = (3600/fTimeBinsPerHour); + Double_t secondsPerBin = (Double_t) binSize; for (int i = 0; i < nSM; i++) { - iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum; - + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + iSM = dataCalibTimeDepCorrection->GetSuperModuleNum(); + + AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM); + AliEMCALSuperModuleCalibMapAPD * dataCalibMapAPD = fCalibMapAPD->GetSuperModuleCalibMapAPDNum(iSM); + AliEMCALSuperModuleBiasAPD * dataBiasAPD = fBiasAPD->GetSuperModuleBiasAPDNum(iSM); + // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { - AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow]; - MGain = fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2], - BiasAPDData[i].fVoltage[iCol][iRow]); - TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain); + AliEMCALCalibMapAPDVal * mapAPD = dataCalibMapAPD->GetAPDVal(iCol, iRow); + gainM = fCalibMapAPD->GetGain(mapAPD->GetPar(0), mapAPD->GetPar(1), mapAPD->GetPar(2), + dataBiasAPD->GetVoltage(iCol, iRow)); + dTempCoeff[iCol][iRow] = GetTempCoeff(mapAPD->GetDarkCurrent(), gainM); + if (fVerbosity > 1) { + cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow + << " par0 " << mapAPD->GetPar(0) + << " par1 " << mapAPD->GetPar(1) + << " par2 " << mapAPD->GetPar(2) + << " bias " << dataBiasAPD->GetVoltage(iCol, iRow) + << " gainM " << gainM << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl; + } } } - // figure out what the reference temperature is, from fCalibAbs - Double_t ReferenceTemperature = 0; + // figure out what the reference temperature is, from fCalibReference + Double_t referenceTemperature = 0; int nVal = 0; for (int iSensor = 0; iSensor0) { // hopefully OK value - ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor]; + if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value + referenceTemperature += dataCalibReference->GetTemperature(iSensor); nVal++; } } if (nVal>0) { - ReferenceTemperature /= nVal; // valid values exist, we can look into corrections + referenceTemperature /= nVal; // valid values exist, we can look into corrections for (int j = 0; j < nBins; j++) { - // what is the timestamp in the middle of this bin? (0.5 is for middle of bin) UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin); // 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) - Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp); - - Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new + Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp); + + Double_t temperatureDiff = dSMTemperature - referenceTemperature ; // new - old + if (fVerbosity > 0) { + cout << " referenceTemperature " << referenceTemperature + << " dSMTemperature " << dSMTemperature + << " temperatureDiff " << temperatureDiff + << endl; + } // if the new temperature is higher than the old/reference one, then the gain has gone down - if (fabs(TemperatureDiff)>fTemperatureResolution) { + // if the new temperature is lower than the old/reference one, then the gain has gone up + // dTempCoeff is a negative number describing how many % (hence 0.01 factor below) the gain + // changes with a positive degree change. + // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up + // The correction we want to keep is what we should multiply our ADC value with as a function + // of time, i.e. the inverse of the gain change.. + if (fabs(temperatureDiff)>fTemperatureResolution) { // significant enough difference that we need to consider it // loop over all towers; effect of temperature change will depend on gain for this tower for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { - correction = TemperatureDiff * TempCoeff[iCol][iRow]; - CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j); + // the correction should be inverse of modification in gain: (see discussion above) + // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01; + // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at: + correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow])*0.01; + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j); + if (fVerbosity > 1) { + cout << " iSM " << iSM + << " iCol " << iCol + << " iRow " << iRow + << " j " << j + << " correction " << correction + << endl; + } } } @@ -861,7 +957,7 @@ Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins) correction = 1; for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { - CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j); + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j); } } } // else