X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALCalibTimeDep.cxx;h=556f31809a4ca2efd2b19416c755499051f0cb55;hb=bb5d026fe994274064d8359c4bf43457108286c0;hp=2dce6cf924dddfa564f4cc8da7b518b00677fab0;hpb=bee992bce2382cf29e66d63476896285a3fd7211;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALCalibTimeDep.cxx b/EMCAL/AliEMCALCalibTimeDep.cxx index 2dce6cf924d..556f31809a4 100644 --- a/EMCAL/AliEMCALCalibTimeDep.cxx +++ b/EMCAL/AliEMCALCalibTimeDep.cxx @@ -19,28 +19,43 @@ ///*-- 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: +// AliEMCALCalibTempCoeff (APD temperature coefficients), +// AliCaloCalibSignal (LED DA), AliEMCALSensorTempArray (ELMB DCS) +// AliEMCALCalibReference: LED amplitude and temperature info at reference time +// +// output/result is in AliEMCALCalibTimeDepCorrection // // /////////////////////////////////////////////////////////////////////////////// #include #include +#include #include "AliLog.h" #include "AliCDBEntry.h" #include "AliCDBManager.h" #include "AliEMCALSensorTempArray.h" +#include "AliCaloCalibSignal.h" +#include "AliEMCALCalibTempCoeff.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; assume the same for all, at least to start with -const double fkTempSlope = 0.017; // 1.7% per deg. C, seems about right for all APDs, from studies at Catania, and by Rachid. -// Note that this is only valid at gain M~30; i.e. voltages at V30. -// At V50, it's more like 2.25%/C, and at V20 about 1.35%/C (also info from Catania) -const double fkNormTemp = 20; // let's say that 20 degrees C is normal as default +// 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 kTempCoeffP0Const = -0.903; // +const double kTempCoeffP0Factor = -1.381e7; // +const double kTempCoeffP1Const = -0.023; // +const double kTempCoeffP1Factor = -4.966e5; // + +const double kTempMaxDiffMedian = 2; // Temperature values should not be further away from median value within SM when considered in the average calc. -const double fkErrorCode = -999; // to indicate that something went wrong +const double kErrorCode = -999; // to indicate that something went wrong using namespace std; @@ -53,10 +68,22 @@ AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() : fEndTime(0), fMinTemp(0), fMaxTemp(0), + fMinTempVariation(0), + fMaxTempVariation(0), + fMinTempValid(15), + fMaxTempValid(35), fMinTime(0), fMaxTime(0), - fRefTemp(fkNormTemp), - fTempArray(NULL) + fTemperatureResolution(0.1), // 0.1 deg C is default + fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference + fTimeBinsPerHour(2), // 2 30-min bins per hour is default + fHighLowGainFactor(16), // factor ~16 between High gain and low gain + fTempArray(NULL), + fCalibSignal(NULL), + fCalibTempCoeff(NULL), + fCalibReference(NULL), + fCalibTimeDepCorrection(NULL), + fVerbosity(0) { // Constructor } @@ -69,10 +96,22 @@ AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) : fEndTime(calibt.GetEndTime()), fMinTemp(calibt.GetMinTemp()), fMaxTemp(calibt.GetMaxTemp()), + fMinTempVariation(calibt.GetMinTempVariation()), + fMaxTempVariation(calibt.GetMaxTempVariation()), + fMinTempValid(calibt.GetMinTempValid()), + fMaxTempValid(calibt.GetMaxTempValid()), fMinTime(calibt.GetMinTime()), fMaxTime(calibt.GetMaxTime()), - fRefTemp(calibt.GetRefTemp()), - fTempArray(calibt.GetTempArray()) + fTemperatureResolution(calibt.GetTemperatureResolution()), + fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()), + fTimeBinsPerHour(calibt.GetTimeBinsPerHour()), + fHighLowGainFactor(calibt.GetHighLowGainFactor()), + fTempArray(calibt.GetTempArray()), + fCalibSignal(calibt.GetCalibSignal()), + fCalibTempCoeff(calibt.GetCalibTempCoeff()), + fCalibReference(calibt.GetCalibReference()), + fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()), + fVerbosity(calibt.GetVerbosity()) { // copy constructor } @@ -103,10 +142,21 @@ void AliEMCALCalibTimeDep::Reset() fEndTime = 0; fMinTemp = 0; fMaxTemp = 0; + fMinTempVariation = 0; + fMaxTempVariation = 0; + fMinTempValid = 15; + fMaxTempValid = 35; fMinTime = 0; fMaxTime = 0; - fRefTemp = fkNormTemp; + fTemperatureResolution = 0.1; // 0.1 deg C is default + fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference + fTimeBinsPerHour = 2; // 2 30-min bins per hour is default fTempArray = NULL; + fCalibSignal = NULL; + fCalibTempCoeff = NULL; + fCalibReference = NULL; + fCalibTimeDepCorrection = NULL; + fVerbosity = 0; return; } @@ -114,44 +164,50 @@ void AliEMCALCalibTimeDep::Reset() void AliEMCALCalibTimeDep::PrintInfo() const { // print some info - cout << endl << " AliEMCALCalibTimeDep::Print() " << endl; + cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl; // basic variables, all 'publicly available' also cout << " VARIABLE DUMP: " << endl << " GetStartTime() " << GetStartTime() << endl << " GetEndTime() " << GetEndTime() << endl + << " GetMinTime() " << GetMinTime() << endl + << " GetMaxTime() " << GetMaxTime() << endl << " GetMinTemp() " << GetMinTemp() << endl << " GetMaxTemp() " << GetMaxTemp() << endl - << " GetRefTemp() " << GetRefTemp() << endl; + << " GetMinTempVariation() " << GetMinTempVariation() << endl + << " GetMaxTempVariation() " << GetMaxTempVariation() << endl + << " GetTemperatureResolution() " << GetTemperatureResolution() << endl; // run ranges cout << " RUN INFO: " << endl + << " runnumber " << GetRunNumber() << endl << " length (in hours) " << GetLengthOfRunInHours() << endl + << " length (in bins) " << GetLengthOfRunInBins() << endl << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours() << " (in deg. C) " << GetRangeOfTempMeasureInDegrees() << endl; - // range in correction values - double corrAtMinTemp = GetCorrection( fMinTemp ); - double corrAtMaxTemp = GetCorrection( fMaxTemp ); - double corrMaxMinDiff = 100*(corrAtMinTemp - corrAtMaxTemp); - cout << " CORRECTION INFO : " << endl - << " corrAtMinTemp " << corrAtMinTemp << endl - << " corrAtMaxTemp " << corrAtMaxTemp << endl - << " corrMaxMinDiff (~%) [=(corrAtMin - corrAtMax)*100] " - << corrMaxMinDiff << endl; return; } + +//________________________________________________________________ +Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const +{ + return (fEndTime - fStartTime)*kSecToHour; +} + //________________________________________________________________ -double AliEMCALCalibTimeDep::GetLengthOfRunInHours() const +Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const { - return (fEndTime - fStartTime)*fkSecToHour; + return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour; } + //________________________________________________________________ -double AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const +Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const { - return (fMaxTime - fMinTime)*fkSecToHour; + return (fMaxTime - fMinTime)*kSecToHour; } + //________________________________________________________________ -double AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const +Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const { return (fMaxTemp - fMinTemp); } @@ -159,7 +215,7 @@ double 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; @@ -168,55 +224,20 @@ void AliEMCALCalibTimeDep::Initialize(Int_t run, // collect the needed information GetTemperatureInfo(); // temperature readings during the run + ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp) return; } //________________________________________________________________ -double AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const -{// return estimate for all SuperModules and sensors, that had data - - // first convert from seconds to hours.. - double timeHour = (timeStamp - fStartTime) * fkSecToHour; - - double average = 0; - int n = 0; - - for (int i=0; iNumSensors(); i++) { - - AliEMCALSensorTemp *st = fTempArray->GetSensor(i); - - // check if we had valid data for the time that is being asked for - if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) { - AliSplineFit *f = st->GetFit(); - if (f) { // ok, looks like we have valid data/info - // let's check what the expected value at the time appears to be - double val = f->Eval(timeHour); - average += val; - n++; - } - } // time - } // loop over fTempArray - - if (n>0) { // some valid data was found - average /= n; - return average; - } - else { // no good data - return fkErrorCode; - } - -} - -//________________________________________________________________ -double AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const +Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const {// return estimate for this one SuperModule, if it had data // first convert from seconds to hours.. - double timeHour = (timeStamp - fStartTime) * fkSecToHour; + Double_t timeHour = (timeStamp - fStartTime) * kSecToHour; - double average = 0; int n = 0; + Double_t valArr[8]={0}; // 8 sensors per SM for (int i=0; iNumSensors(); i++) { @@ -228,10 +249,14 @@ double AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const AliSplineFit *f = st->GetFit(); if (f) { // ok, looks like we have valid data/info // let's check what the expected value at the time appears to be - double val = f->Eval(timeHour); - cout << " i " << i << " val " << val << endl; - average += val; - n++; + Double_t val = f->Eval(timeHour); + if ( fVerbosity > 0 ) { + cout << " sensor i " << i << " val " << val << endl; + } + if (val>fMinTempValid && val0) { // some valid data was found - average /= n; - return average; + Double_t median = TMath::Median(n, valArr); + Double_t average = 0; + Int_t nval = 0; + for (int is=0; is 0) { + average /= nval; + return average; + } + else { // this case should not happen, but kept for completeness (coverity etc) + return median; + } } else { // no good data - return fkErrorCode; + return kErrorCode; } } //________________________________________________________________ -double AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_t timeStamp) const -{// return estimate for this one SuperModule and sensor, if it had data +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, 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 temperature coefficient info, to estimate correction. + For now (from Dec 2009), we do not use LED info. + */ + + // 0: Init + // how many SuperModules do we have? + Int_t nSM = fCalibReference->GetNSuperModule(); + // how many time-bins should we have for this run? + Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1) + 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.. + } + if (nBins == 1) { + binSize = fEndTime - fStartTime; + } + if (fVerbosity > 0) { + cout << " nBins " << nBins << " binSize " << binSize << endl; + } - // first convert from seconds to hours.. - double timeHour = (timeStamp - fStartTime) * fkSecToHour; + // 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); - for (int i=0; iNumSensors(); i++) { - - AliEMCALSensorTemp *st = fTempArray->GetSensor(i); - int module = st->GetSector()*2 + st->GetSide(); - if ( module == imod && st->GetNum()==isens ) { // right module, and sensor - // check if we had valid data for the time that is being asked for - if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) { - AliSplineFit *f = st->GetFit(); - if (f) { // ok, looks like we have valid data/info - // let's check what the expected value at the time appears to be - double val = f->Eval(timeHour); - - return val; // no point to move further in for loop, we have found the sensor we were looking for - } - } // time - } - - } // loop over fTempArray - - // 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; + // 2: try with Temperature correction + Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize); + return nRemaining; } + //________________________________________________________________ -double AliEMCALCalibTimeDep::GetCorrection(double temperature) const -{ // gain decreases with increasing temperature - double gain = (1.0 - (temperature-fRefTemp)*fkTempSlope); - // i.e. multiplicative correction to ADC increases with increasing temperature - return 1.0/gain; +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 dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark; + Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark; + + Double_t dTC = dP0 + dP1*M; + // from % numbers to regular ones..: + dTC *= 0.01; + + return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion } -/* Next comes the method that does the work in picking up all the needed info..*/ +/* Next come the methods that do the work in picking up all the needed info..*/ //________________________________________________________________ void AliEMCALCalibTimeDep::GetTemperatureInfo() { @@ -309,3 +366,554 @@ void AliEMCALCalibTimeDep::GetTemperatureInfo() return; } + +//________________________________________________________________ +Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() +{// assign max/min time and temperature values + + 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; + + Int_t n = 0; // number of valid readings + + 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 + 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; + int nval = 0; + for (int ip=0; ipfMinTempValid && y0[ip] y0[ip]) { min = y0[ip]; } + if (max < y0[ip]) { max = y0[ip]; } + nval++; + } + } + if (nval>0) { + 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 + + if (n>0) { // some valid data was found + return n; + } + else { // no good data + return (Int_t) kErrorCode; + } + +} + +//________________________________________________________________ +void AliEMCALCalibTimeDep::GetCalibSignalInfo() +{ + // pick up Preprocessor output, based on fRun (most recent version) + AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun); + if (entry) { + fCalibSignal = (AliCaloCalibSignal *) entry->GetObject(); + } + + if (fCalibSignal) { + AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld", + fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(), + fCalibSignal->GetTreeAmpVsTime()->GetEntries(), + fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(), + fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(), + fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) ); + } + else { + AliWarning( Form("AliCaloCalibSignal not found!") ); + } + + return; +} + +//________________________________________________________________ +void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo() +{ + // pick up Preprocessor output, based on fRun (most recent version) + AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun); + // stored object should be a TTree; read the info + if (entry) { + fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject(); + } + + if (fCalibTempCoeff) { + AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) ); + } + else { + AliWarning( Form("AliEMCALCalibTempCoeff not found!") ); + } + + return; +} + +//________________________________________________________________ +void AliEMCALCalibTimeDep::GetCalibReferenceInfo() +{ + // pick up Preprocessor output, based on fRun (most recent version) + AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun); + if (entry) { + fCalibReference = (AliEMCALCalibReference *) entry->GetObject(); + } + + if (fCalibReference) { + AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) ); + } + else { + AliWarning( Form("AliEMCALCalibReference not found!") ); + } + + return; +} + +//________________________________________________________________ +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 fCalibReference + + // 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 + for (int i = 0; i < nSM; i++) { + 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; + } + } + + int iSM = 0; + int iCol = 0; + int iRow = 0; + int iStrip = 0; + int iGain = 0; + + // 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/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 + TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; + TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; + TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; + TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; + + // set up TArray's first + for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) { + for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + for (iGain = 0; iGain < 2; iGain++) { + // length of arrays + ampT[iSM][iCol][iRow][iGain].Set(nBins); + nT[iSM][iCol][iRow][iGain].Set(nBins); + // content of arrys + for (int j = 0; j < nBins; j++) { + ampT[iSM][iCol][iRow][iGain].AddAt(0, j); + nT[iSM][iCol][iRow][iGain].AddAt(0, j); + } + } + } + }//iCol + for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) { + for (iGain = 0; iGain < 2; iGain++) { + // length of arrays + ampLEDRefT[iSM][iStrip][iGain].Set(nBins); + nLEDRefT[iSM][iStrip][iGain].Set(nBins); + // content of arrys + for (int j = 0; j < nBins; j++) { + ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j); + nLEDRefT[iSM][iStrip][iGain].AddAt(0, j); + } + } + }//iStrip + } + + // OK, now loop over the TTrees and fill the arrays needed for R(T) + TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime(); + TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime(); + + int iChannelNum = 0; // for regular towers + int iRefNum = 0; // for LED + double dHour = 0; + double dAvgAmp = 0; + + 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++) { + treeAvg->GetEntry(ient); + // figure out where this info comes from + 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)+dAvgAmp, iBin ); + nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin ); + nCount[iSM]++; + } + + treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum); + treeLEDRefAvg->SetBranchAddress("fHour", &dHour); + treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp); + + for (int ient=0; ientGetEntries(); ient++) { + treeLEDRefAvg->GetEntry(ient); + // figure out where this info comes from + 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)+dAvgAmp, iBin ); + nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin ); + nCountLEDRef[iSM]++; + } + + // Normalize TArray values, and calculate average also + Float_t norm = 0; // extra var, for readability + + for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) { + if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data.. + 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 + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + for (iGain = 0; iGain < 2; iGain++) { + // content of arrys + for (int j = 0; j < nBins; j++) { + if (nT[iSM][iCol][iRow][iGain].At(j) > 0) { + norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j); + ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt + } + } + } + } + }//iCol + for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) { + for (iGain = 0; iGain < 2; iGain++) { + for (int j = 0; j < nBins; j++) { + if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) { + norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j); + ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt + } + } + } + }//iStrip + } + } // iSM + + + // Calculate correction values, and store them + // set kErrorCode values for those that could not be set + + 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++) { + 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): + AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow); + iGain = refVal->GetHighLow(); // gain value used for reference calibration + // valid amplitude values should be larger than 0 + if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) { + ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip); + } + else { + ratiot0 = kErrorCode; + } + + // Calc. R(T) + for (int j = 0; j < nBins; j++) { + + // calculate R(T) also; first try with individual tower: + // 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 + 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 + 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 && ratioT>0) { + correction = ratiot0/ratioT; + } + else { + correction = kErrorCode; + nRemaining++; + } + + // Store the value + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j); + /* Check that + fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction); + also works OK */ + } // nBins + } + } + } + + nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins); + return nRemaining; +} + +//________________________________________________________________ +Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) +{ // use averages for each strip if no good values exist for some single tower + + // go over fTimeDepCorrection info + Int_t nRemaining = 0; // we count the towers for which we could not get valid data + + // for calculating StripAverage info + int nValidTower = 0; + Float_t stripAverage = 0; + Float_t val = 0; + + int iSM = 0; + int iCol = 0; + int iRow = 0; + int iStrip = 0; + int firstCol = 0; + int lastCol = 0; + + for (int i = 0; i < nSM; i++) { + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + iSM = dataCalibTimeDepCorrection->GetSuperModuleNum(); + + for (int j = 0; j < nBins; j++) { + + nValidTower = 0; + stripAverage = 0; + + for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) { + firstCol = iStrip*2; + if ((iSM%2)==1) { // C side + firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2; + } + lastCol = firstCol+1; + + for (iCol = firstCol; iCol <= lastCol; iCol++) { + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j); + if (val>0) { // valid value; error code is negative + stripAverage += val; + nValidTower++; + } + } + } + + // calc average over strip + if (nValidTower>0) { + stripAverage /= nValidTower; + for (iCol = firstCol; iCol <= lastCol; iCol++) { + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j); + if (val<0) { // invalid value; error code is negative + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j); + } + } + } + } + else { // could not fill in unvalid entries + nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows; + } + + } // iStrip + } // j, bins + } // iSM + + return nRemaining; +} + +//________________________________________________________________ +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 + // based on that instead. + // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class + Int_t nRemaining = 0; + + int iSM = 0; + int iCol = 0; + int iRow = 0; + + Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows]; + memset(dTempCoeff, 0, sizeof(dTempCoeff)); + Double_t correction = 0; + Double_t secondsPerBin = (Double_t) binSize; + + for (int i = 0; i < nSM; i++) { + AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i); + iSM = dataCalibTimeDepCorrection->GetSuperModuleNum(); + + AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM); + AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM); + + // first get CalibTempCoeff info + for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + + dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow); + if (fVerbosity > 1) { + cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow + << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl; + } + } + } + + // figure out what the reference temperature is, from fCalibReference + Double_t referenceTemperature = 0; + int nVal = 0; + for (int iSensor = 0; iSensorGetTemperature(iSensor)>0) { // hopefully OK value + referenceTemperature += dataCalibReference->GetTemperature(iSensor); + nVal++; + } + } + + if (nVal>0) { + 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 dSMTemperature = GetTemperatureSM(iSM, timeStamp); + + Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new + if (fVerbosity > 0) { + cout << " referenceTemperature " << referenceTemperature + << " dSMTemperature " << dSMTemperature + << " temperatureDiff " << temperatureDiff + << endl; + } + // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down + // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up + // dTempCoeff is a (unsigned) factor describing how many % the gain + // changes with a 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 ( (TMath::Abs(temperatureDiff)>fTemperatureResolution) + && (TMath::Abs(temperatureDiff)GetCorrection(iCol,iRow)->AddAt(correction, j); + if (fVerbosity > 1) { + cout << " iSM " << iSM + << " iCol " << iCol + << " iRow " << iRow + << " j " << j + << " correction " << correction + << endl; + } + } + } + + } // if noteworthy temperature change + else { // just set correction values to 1.0 + correction = 1; + for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) { + for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) { + dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j); + } + } + } // else + } // j, Bins + + } // if reference temperature exist + else { // could not do the needed check.. signal that in the return code + nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins; + } + } // iSM + + return nRemaining; +} +