replacing fabs with TMath::Abs
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALCalibTimeDep.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliEMCALCalibTimeDep.cxx $ */
17
18 //_________________________________________________________________________
19 ///*-- Author: 
20 ///////////////////////////////////////////////////////////////////////////////
21 //                                                                           //
22 // class for EMCAL time-dep calibration   
23 // - 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
28 //
29 // output/result is in AliEMCALCalibTimeDepCorrection 
30 //                                                                           //
31 ///////////////////////////////////////////////////////////////////////////////
32
33 #include <iostream>
34 #include <TGraphSmooth.h>
35 #include <TMath.h>
36 #include "AliLog.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"
45
46 /* first a bunch of constants.. */
47 const double kSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
48
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; //
55  
56 const double kErrorCode = -999; // to indicate that something went wrong
57
58 using namespace std;
59
60 ClassImp(AliEMCALCalibTimeDep)
61
62 //________________________________________________________________
63 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
64   fRun(0),
65   fStartTime(0),
66   fEndTime(0),
67   fMinTemp(0),
68   fMaxTemp(0),
69   fMinTempVariation(0),
70   fMaxTempVariation(0),
71   fMinTime(0),
72   fMaxTime(0),
73   fTemperatureResolution(0.1), // 0.1 deg C is default
74   fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference 
75   fTimeBinsPerHour(2), // 2 30-min bins per hour is default
76   fHighLowGainFactor(16), // factor ~16 between High gain and low gain
77   fTempArray(NULL),
78   fCalibSignal(NULL),
79   fCalibTempCoeff(NULL),
80   fCalibReference(NULL),
81   fCalibTimeDepCorrection(NULL),
82   fVerbosity(0)
83 {
84   // Constructor
85 }
86
87 //________________________________________________________________
88 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
89   TObject(calibt),
90   fRun(calibt.GetRunNumber()),
91   fStartTime(calibt.GetStartTime()),
92   fEndTime(calibt.GetEndTime()),
93   fMinTemp(calibt.GetMinTemp()),
94   fMaxTemp(calibt.GetMaxTemp()),
95   fMinTempVariation(calibt.GetMinTempVariation()),
96   fMaxTempVariation(calibt.GetMaxTempVariation()),
97   fMinTime(calibt.GetMinTime()),
98   fMaxTime(calibt.GetMaxTime()),
99   fTemperatureResolution(calibt.GetTemperatureResolution()),
100   fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()),
101   fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
102   fHighLowGainFactor(calibt.GetHighLowGainFactor()),
103   fTempArray(calibt.GetTempArray()),
104   fCalibSignal(calibt.GetCalibSignal()),
105   fCalibTempCoeff(calibt.GetCalibTempCoeff()),
106   fCalibReference(calibt.GetCalibReference()),
107   fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
108   fVerbosity(calibt.GetVerbosity())
109 {
110   // copy constructor
111 }
112
113
114 //________________________________________________________________
115 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
116 {
117   // assignment operator; use copy ctor
118   if (&calibt == this) return *this;
119
120   new (this) AliEMCALCalibTimeDep(calibt);
121   return *this;
122 }
123
124 //________________________________________________________________
125 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
126 {
127   // Destructor
128 }
129
130 //________________________________________________________________
131 void  AliEMCALCalibTimeDep::Reset() 
132 {
133   // clear variables to default
134   fRun = 0;
135   fStartTime = 0;
136   fEndTime = 0;
137   fMinTemp = 0;
138   fMaxTemp = 0;
139   fMinTempVariation = 0;
140   fMaxTempVariation = 0;
141   fMinTime = 0;
142   fMaxTime = 0;
143   fTemperatureResolution = 0.1; // 0.1 deg C is default
144   fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference 
145   fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
146   fTempArray = NULL;
147   fCalibSignal = NULL;
148   fCalibTempCoeff = NULL;
149   fCalibReference = NULL;
150   fCalibTimeDepCorrection = NULL;
151   fVerbosity = 0;
152   return;
153 }
154
155 //________________________________________________________________
156 void  AliEMCALCalibTimeDep::PrintInfo() const
157 {
158   // print some info
159   cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
160   // basic variables, all 'publicly available' also
161   cout << " VARIABLE DUMP: " << endl
162        << " GetStartTime() " << GetStartTime() << endl
163        << " GetEndTime() " << GetEndTime() << endl
164        << " GetMinTime() " << GetMinTime() << endl
165        << " GetMaxTime() " << GetMaxTime() << endl
166        << " GetMinTemp() " << GetMinTemp() << endl
167        << " GetMaxTemp() " << GetMaxTemp() << endl
168        << " GetMinTempVariation() " << GetMinTempVariation() << endl
169        << " GetMaxTempVariation() " << GetMaxTempVariation() << endl;
170   // run ranges
171   cout << " RUN INFO: " << endl
172        << " runnumber " << GetRunNumber() << endl
173        << " length (in hours) " << GetLengthOfRunInHours() << endl
174        << " length (in bins) " << GetLengthOfRunInBins() << endl
175        << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
176        << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
177        << endl;
178
179   return;
180 }
181
182 //________________________________________________________________ 
183 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
184 {
185   return (fEndTime - fStartTime)*kSecToHour;
186 }
187
188 //________________________________________________________________ 
189 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
190 {
191   return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
192 }
193
194 //________________________________________________________________ 
195 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
196 {
197   return (fMaxTime - fMinTime)*kSecToHour;
198 }
199
200 //________________________________________________________________ 
201 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
202 {
203   return (fMaxTemp - fMinTemp);
204 }
205
206 //________________________________________________________________
207 void AliEMCALCalibTimeDep::Initialize(Int_t run, 
208                                       UInt_t startTime, UInt_t endTime)
209 { // setup, and get temperature info
210   Reset(); // start fresh
211
212   fRun = run;
213   fStartTime = startTime;
214   fEndTime = endTime;
215   
216   // collect the needed information
217   GetTemperatureInfo(); // temperature readings during the run
218   ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
219
220   return;
221 }
222
223 //________________________________________________________________
224 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
225 {// return estimate for this one SuperModule, if it had data 
226
227   // first convert from seconds to hours..
228   Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
229
230   Double_t average = 0;
231   int n = 0;
232
233   for (int i=0; i<fTempArray->NumSensors(); i++) {
234     
235     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
236     int module = st->GetSector()*2 + st->GetSide();
237     if ( module == imod ) { // right module
238       // check if we had valid data for the time that is being asked for
239       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
240         AliSplineFit *f = st->GetFit();
241         if (f) { // ok, looks like we have valid data/info
242           // let's check what the expected value at the time appears to be
243           Double_t val = f->Eval(timeHour);
244           if ( fVerbosity > 0 ) {
245             cout << " sensor i " << i << " val " << val << endl;
246           }
247           average += val;
248           n++;
249         }
250       } // time
251     }
252     
253   } // loop over fTempArray
254   
255   if (n>0) { // some valid data was found
256     average /= n;
257     return average;
258   }
259   else { // no good data
260     return kErrorCode;
261   }
262
263 }
264
265 //________________________________________________________________
266 Int_t AliEMCALCalibTimeDep::CalcCorrection()
267 { // OK, this is where the real action takes place - the heart of this class..
268   /* The philosophy is as follows:
269      0. Init corrections to 1.0 values, and see how many correction bins we need
270      1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
271      2. try to use temperature info + APD temperature coefficient info, to estimate correction.
272      For now (from Dec 2009), we do not use LED info.
273    */
274
275   // 0: Init
276   // how many SuperModules do we have?
277   Int_t nSM = fCalibReference->GetNSuperModule();
278   // how many time-bins should we have for this run?
279   Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
280   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
281
282   // 1: get info on how much individual sensors might have changed during
283   // the run (compare max-min for each sensor separately)
284   if (fMaxTempVariation < fTemperatureResolution) {
285     nBins = 1; // just one bin needed..
286     binSize = fEndTime - fStartTime;
287   }
288
289   // set up a reasonable default (correction = 1.0)
290   fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
291   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
292   fCalibTimeDepCorrection->SetStartTime(fStartTime);
293   fCalibTimeDepCorrection->SetNTimeBins(nBins);
294   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
295
296   // 2: try with Temperature correction 
297   Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
298
299   return nRemaining;
300 }
301
302
303 //________________________________________________________________
304 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
305 { // estimate the Temperature Coefficient, based on the dark current (IDark)
306   // and the gain (M), based on Catania parameterizations
307
308   Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
309   Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
310  
311   Double_t dTC = dP0 + dP1*M;
312   // from % numbers to regular ones..:
313   dTC *= 0.01;
314
315   return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
316 }
317
318 /* Next come the methods that do the work in picking up all the needed info..*/
319 //________________________________________________________________
320 void AliEMCALCalibTimeDep::GetTemperatureInfo() 
321 {
322   // pick up Preprocessor output, based on fRun (most recent version)
323   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
324   if (entry) {
325     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
326   }
327
328   if (fTempArray) { 
329     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
330                   fTempArray->NumSensors(),
331                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
332   }
333   else {
334     AliWarning( Form("AliEMCALSensorTempArray not found!") );
335   }
336   
337   return;
338 }
339
340 //________________________________________________________________
341 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() 
342 {// assign max/min time and temperature values
343
344   fMinTemp = 999; // init to some large value (999 deg C)
345   fMaxTemp = 0;
346   fMinTempVariation = 999; // init to some large value (999 deg C)
347   fMaxTempVariation = 0;
348   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
349   fMaxTime = 0;
350
351   Int_t n = 0; // number of valid readings
352
353   for (int i=0; i<fTempArray->NumSensors(); i++) {
354     
355     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
356     if ( st->GetStartTime() == 0 ) { // no valid data
357       continue;
358     }
359
360     // check time ranges
361     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
362     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
363  
364     // check temperature ranges
365     AliSplineFit *f = st->GetFit();
366
367     if (f) { // ok, looks like we have valid data/info
368       int np = f->GetKnots();
369       Double_t *y0 = f->GetY0();
370       // min and max values within the single sensor
371       Double_t min = 999;
372       Double_t max = 0;
373       for (int ip=0; ip<np; ip++) { 
374         if (min > y0[ip]) { min = y0[ip]; }
375         if (max < y0[ip]) { max = y0[ip]; }
376       }
377       if (fMinTemp > min) { fMinTemp = min; }
378       if (fMaxTemp < max) { fMaxTemp = max; }
379       Double_t variation = max - min;
380       if (fMinTempVariation > variation) { fMinTempVariation = variation; }
381       if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
382
383       n++;
384     }
385   } // loop over fTempArray
386   
387   if (n>0) { // some valid data was found
388     return n;
389   }
390   else { // no good data
391     return (Int_t) kErrorCode;
392   }
393
394 }
395
396 //________________________________________________________________
397 void AliEMCALCalibTimeDep::GetCalibSignalInfo() 
398 {
399   // pick up Preprocessor output, based on fRun (most recent version)
400   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
401   if (entry) {
402     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
403   }
404
405   if (fCalibSignal) { 
406     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
407                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
408                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
409                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
410                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
411                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );               
412   }
413   else {
414     AliWarning( Form("AliCaloCalibSignal not found!") );
415   }
416   
417   return;
418 }
419
420 //________________________________________________________________
421 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo() 
422 {
423   // pick up Preprocessor output, based on fRun (most recent version)
424   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
425   // stored object should be a TTree; read the info
426   if (entry) {
427     fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
428   }
429
430   if (fCalibTempCoeff) { 
431     AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
432   }
433   else {
434     AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
435   }
436   
437   return;
438 }
439
440 //________________________________________________________________
441 void AliEMCALCalibTimeDep::GetCalibReferenceInfo() 
442 {
443   // pick up Preprocessor output, based on fRun (most recent version)
444   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
445   if (entry) {
446     fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
447   }
448
449   if (fCalibReference) { 
450     AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
451   }
452   else {
453     AliWarning( Form("AliEMCALCalibReference not found!") );
454   }
455   
456   return;
457 }
458
459 //________________________________________________________________
460 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) 
461 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0 
462   // The correction factor we keep is c(T) = R(t0)/R(T)
463   // T info from fCalibSignal, t0 info from fCalibReference
464
465   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
466   // but one could upgrade this in the future
467   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
468
469   // sanity check; same SuperModule indices for corrections as for regular calibrations
470   for (int i = 0; i < nSM; i++) {
471     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
472     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
473
474     int iSMRef = dataCalibReference->GetSuperModuleNum();
475     int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
476     if (iSMRef != iSMCorr) {
477       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
478       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
479       return nRemaining;
480     }
481   }
482  
483   int iSM = 0;
484   int iCol = 0;
485   int iRow = 0;
486   int iStrip = 0;
487   int iGain = 0;
488
489   // The fCalibSignal info is stored in TTrees
490   // Note that the time-bins for the TTree's may not exactly match with our correction time bins 
491   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
492   // fCalibSignal time info in seconds: Hour/kSecToHour
493   // corrected for startTime difference: Hour/kSecToHour + timeDiff
494   // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
495
496   // values for R(T), size of TArray = nBins
497   // the [2] dimension below is for the low or high gain 
498   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
499   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
500   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
501   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
502
503   // set up TArray's first
504   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
505     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
506       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
507         for (iGain = 0; iGain < 2; iGain++) {
508           // length of arrays
509           ampT[iSM][iCol][iRow][iGain].Set(nBins);
510           nT[iSM][iCol][iRow][iGain].Set(nBins);
511           // content of arrys
512           for (int j = 0; j < nBins; j++) {
513             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
514             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
515           }
516         }
517       }
518     }//iCol
519     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
520       for (iGain = 0; iGain < 2; iGain++) {
521         // length of arrays
522         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
523         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
524         // content of arrys
525         for (int j = 0; j < nBins; j++) {
526           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
527           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
528         }
529       }
530     }//iStrip
531   }
532
533   // OK, now loop over the TTrees and fill the arrays needed for R(T)
534   TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
535   TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
536
537   int iChannelNum = 0; // for regular towers
538   int iRefNum = 0; // for LED
539   double dHour = 0;
540   double dAvgAmp = 0;
541
542   treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
543   treeAvg->SetBranchAddress("fHour", &dHour);
544   treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
545
546   int iBin = 0;
547   // counters for how many values were seen per SuperModule
548   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
549   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
550
551   for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
552     treeAvg->GetEntry(ient);
553     // figure out where this info comes from
554     fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
555     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
556     // add value in the arrays
557     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
558     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
559     nCount[iSM]++;
560   }
561
562   treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
563   treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
564   treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
565
566   for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
567     treeLEDRefAvg->GetEntry(ient);
568     // figure out where this info comes from
569     fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
570     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
571     // add value in the arrays
572     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
573     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
574     nCountLEDRef[iSM]++;
575   }
576
577   // Normalize TArray values, and calculate average also
578   Float_t norm = 0; // extra var, for readability
579
580   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
581     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
582       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
583         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
584         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
585         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
586           for (iGain = 0; iGain < 2; iGain++) {
587             // content of arrys
588             for (int j = 0; j < nBins; j++) {
589               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) { 
590                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j); 
591                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
592               } 
593             }
594           }
595         }
596       }//iCol
597       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
598         for (iGain = 0; iGain < 2; iGain++) {
599           for (int j = 0; j < nBins; j++) {
600             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
601               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j); 
602               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
603             }
604           }
605         }
606       }//iStrip
607     }
608   } // iSM
609
610
611   // Calculate correction values, and store them
612   // set kErrorCode values for those that could not be set
613
614   Float_t ratiot0 = 0;
615   Float_t ratioT = 0;
616   Float_t correction = 0; // c(T) = R(t0)/R(T)
617   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
618
619   for (int i = 0; i < nSM; i++) {
620     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
621     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
622     iSM = dataCalibReference->GetSuperModuleNum();
623
624     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
625       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
626       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
627         refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration      
628
629       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
630
631         // Calc. R(t0):
632         AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
633         iGain = refVal->GetHighLow(); // gain value used for reference calibration      
634         // valid amplitude values should be larger than 0
635         if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
636           ratiot0 =  refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
637         }
638         else {
639           ratiot0 = kErrorCode;
640         }
641
642         // Calc. R(T)
643         for (int j = 0; j < nBins; j++) {
644
645           // calculate R(T) also; first try with individual tower:
646           // same gain as for reference calibration is the default
647           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
648             // looks like valid data with the right gain combination
649             ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); 
650
651             // if data appears to be saturated, and we are in high gain, then try with low gain instead
652             int newGain = iGain;
653             int newRefGain = refGain;
654             if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
655               newGain = 0;
656             }
657             if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) { 
658               newRefGain = 0;
659             }
660
661             if (newGain!=iGain || newRefGain!=refGain) {
662               // compensate for using different gain than in the reference calibration
663               // we may need to have a custom H/L ratio value for each tower
664               // later, but for now just use a common value, as the rest of the code does.. 
665               ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j); 
666
667               if (newGain<iGain) {
668                 ratioT *= fHighLowGainFactor;
669               }
670               else if (newRefGain<refGain) {
671                 ratioT /= fHighLowGainFactor;
672               }
673             }
674           }
675           else {
676             ratioT = kErrorCode;
677           }
678
679           // Calc. correction factor
680           if (ratiot0>0 && ratioT>0) { 
681             correction = ratiot0/ratioT;
682           }
683           else {
684             correction = kErrorCode;
685             nRemaining++;
686           }
687
688           // Store the value
689           dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
690           /* Check that
691           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
692           also works OK */
693         } // nBins
694       }
695     }
696   }
697
698   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
699   return nRemaining;
700 }
701
702 //________________________________________________________________
703 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) 
704 { // use averages for each strip if no good values exist for some single tower 
705
706   // go over fTimeDepCorrection info
707   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
708
709   // for calculating StripAverage info
710   int nValidTower = 0;
711   Float_t stripAverage = 0;
712   Float_t val = 0;
713
714   int iSM = 0;
715   int iCol = 0;
716   int iRow = 0;
717   int iStrip = 0;
718   int firstCol = 0;
719   int lastCol = 0;
720
721   for (int i = 0; i < nSM; i++) {
722     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
723     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
724
725     for (int j = 0; j < nBins; j++) {
726
727       nValidTower = 0;
728       stripAverage = 0;
729
730       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
731         firstCol = iStrip*2;
732         if ((iSM%2)==1) { // C side
733           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
734         }
735         lastCol = firstCol+1;
736
737         for (iCol = firstCol; iCol <= lastCol; iCol++) {
738           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
739             val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
740             if (val>0) { // valid value; error code is negative
741               stripAverage += val;
742               nValidTower++;
743             } 
744           }
745         }
746
747         // calc average over strip
748         if (nValidTower>0) {
749           stripAverage /= nValidTower;
750           for (iCol = firstCol; iCol <= lastCol; iCol++) {
751             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
752               val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
753               if (val<0) { // invalid value; error code is negative
754                 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
755               }
756             }
757           }
758         }
759         else { // could not fill in unvalid entries
760           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
761         }
762
763       } // iStrip
764     } // j, bins
765   } // iSM
766
767   return nRemaining;
768 }
769
770 //________________________________________________________________
771 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize) 
772 { // OK, so we didn't have valid LED data that allowed us to do the correction only 
773   // with that info.
774   // So, instead we'll rely on the temperature info and try to do the correction
775   // based on that instead.
776   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class 
777   Int_t nRemaining = 0;
778
779   int iSM = 0;
780   int iCol = 0;
781   int iRow = 0;
782
783   Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
784   memset(dTempCoeff, 0, sizeof(dTempCoeff));
785   Double_t correction = 0;
786   Double_t secondsPerBin = (Double_t) binSize;
787
788   for (int i = 0; i < nSM; i++) {
789     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
790     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
791
792     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
793     AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
794
795     // first get CalibTempCoeff info
796     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
797       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
798
799         dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
800         if (fVerbosity > 1) {
801           cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow 
802                << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl; 
803         }
804       }
805     }
806     
807     // figure out what the reference temperature is, from fCalibReference
808     Double_t referenceTemperature = 0;
809     int nVal = 0;
810     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
811       if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
812         referenceTemperature += dataCalibReference->GetTemperature(iSensor);
813         nVal++;
814       }
815     }
816     
817     if (nVal>0) {
818       referenceTemperature /= nVal; // valid values exist, we can look into corrections
819       
820       for (int j = 0; j < nBins; j++) {
821         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
822         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
823       // 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)
824         Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp); 
825
826         Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
827         if (fVerbosity > 0) {
828           cout << " referenceTemperature " << referenceTemperature
829                << " dSMTemperature " << dSMTemperature
830                << " temperatureDiff " << temperatureDiff
831                << endl;
832         }
833         // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down 
834         // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
835         // dTempCoeff is a (unsigned) factor describing how many % the gain 
836         // changes with a degree change.  
837         // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
838         // The correction we want to keep is what we should multiply our ADC value with as a function
839         // of time, i.e. the inverse of the gain change..
840         if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
841              && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
842           // significant enough difference that we need to consider it, and also not unreasonably large
843
844           // loop over all towers; effect of temperature change will depend on gain for this tower
845           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
846             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
847
848               // the correction should be inverse of modification in gain: (see discussion above)
849               // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01; 
850               // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
851               correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]); 
852               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
853               if (fVerbosity > 1) {
854                 cout << " iSM " << iSM
855                      << " iCol  " << iCol 
856                      << " iRow  " << iRow 
857                      << " j  " << j 
858                      << " correction  " << correction 
859                      << endl;
860               }
861             }
862           }
863
864         } // if noteworthy temperature change
865         else { // just set correction values to 1.0
866           correction = 1;
867           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
868             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
869               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
870             }
871           }
872         } // else
873       } // j, Bins
874
875     } // if reference temperature exist 
876     else { // could not do the needed check.. signal that in the return code
877       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
878     }
879   } // iSM
880
881   return nRemaining;
882 }
883