d8c655a18db7277f626ea22621eeacf855a37c30
[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 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include <iostream>
27 #include <TGraphSmooth.h>
28 #include "AliLog.h"
29 #include "AliCDBEntry.h"
30 #include "AliCDBManager.h"
31 #include "AliEMCALSensorTempArray.h"
32 #include "AliCaloCalibSignal.h"
33 #include "AliEMCALBiasAPD.h"
34 #include "AliEMCALCalibMapAPD.h"
35 #include "AliEMCALCalibAbs.h"
36 #include "AliEMCALCalibTimeDepCorrection.h" 
37 #include "AliEMCALCalibTimeDep.h"
38
39 /* first a bunch of constants.. */
40 const double fkSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
41
42 // some global variables for APD handling; values from Catania studies, best fit
43 // TempCoeff = p0+p1*M (M=gain), where p0 and and p1 are functions of the dark current
44 const double fkTempCoeffP0Const = -0.903; // 
45 const double fkTempCoeffP0Factor = -1.381e7; // 
46 const double fkTempCoeffP1Const = -0.023; // 
47 const double fkTempCoeffP1Factor = -4.966e5; //
48  
49 const double fkErrorCode = -999; // to indicate that something went wrong
50
51 using namespace std;
52
53 ClassImp(AliEMCALCalibTimeDep)
54
55 //________________________________________________________________
56 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
57   fRun(0),
58   fStartTime(0),
59   fEndTime(0),
60   fMinTemp(0),
61   fMaxTemp(0),
62   fMinTime(0),
63   fMaxTime(0),
64   fTemperatureResolution(0.1), // 0.1 deg C is default
65   fTimeBinsPerHour(2), // 2 30-min bins per hour is default
66   fTempArray(NULL),
67   fCalibSignal(NULL),
68   fBiasAPD(NULL),
69   fCalibMapAPD(NULL),
70   fCalibAbs(NULL),
71   fCalibTimeDepCorrection(NULL)
72 {
73   // Constructor
74 }
75
76 //________________________________________________________________
77 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
78   TObject(calibt),
79   fRun(calibt.GetRunNumber()),
80   fStartTime(calibt.GetStartTime()),
81   fEndTime(calibt.GetEndTime()),
82   fMinTemp(calibt.GetMinTemp()),
83   fMaxTemp(calibt.GetMaxTemp()),
84   fMinTime(calibt.GetMinTime()),
85   fMaxTime(calibt.GetMaxTime()),
86   fTemperatureResolution(calibt.GetTemperatureResolution()),
87   fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
88   fTempArray(calibt.GetTempArray()),
89   fCalibSignal(calibt.GetCalibSignal()),
90   fBiasAPD(calibt.GetBiasAPD()),
91   fCalibMapAPD(calibt.GetCalibMapAPD()),
92   fCalibAbs(calibt.GetCalibAbs()),
93   fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection())
94 {
95   // copy constructor
96 }
97
98
99 //________________________________________________________________
100 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
101 {
102   // assignment operator; use copy ctor
103   if (&calibt == this) return *this;
104
105   new (this) AliEMCALCalibTimeDep(calibt);
106   return *this;
107 }
108
109 //________________________________________________________________
110 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
111 {
112   // Destructor
113 }
114
115 //________________________________________________________________
116 void  AliEMCALCalibTimeDep::Reset() 
117 {
118   // clear variables to default
119   fRun = 0;
120   fStartTime = 0;
121   fEndTime = 0;
122   fMinTemp = 0;
123   fMaxTemp = 0;
124   fMinTime = 0;
125   fMaxTime = 0;
126   fTemperatureResolution = 0.1; // 0.1 deg C is default
127   fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
128   fTempArray = NULL;
129   fCalibSignal = NULL;
130   fBiasAPD = NULL;
131   fCalibMapAPD = NULL;
132   fCalibAbs = NULL;
133   fCalibTimeDepCorrection = NULL;
134   return;
135 }
136
137 //________________________________________________________________
138 void  AliEMCALCalibTimeDep::PrintInfo() const
139 {
140   // print some info
141   cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
142   // basic variables, all 'publicly available' also
143   cout << " VARIABLE DUMP: " << endl
144        << " GetStartTime() " << GetStartTime() << endl
145        << " GetEndTime() " << GetEndTime() << endl
146        << " GetMinTemp() " << GetMinTemp() << endl
147        << " GetMaxTemp() " << GetMaxTemp() << endl;
148   // run ranges
149   cout << " RUN INFO: " << endl
150        << " length (in hours) " << GetLengthOfRunInHours() << endl
151        << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
152        << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
153        << endl;
154
155   return;
156 }
157
158 //________________________________________________________________ 
159 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
160 {
161   return (fEndTime - fStartTime)*fkSecToHour;
162 }
163
164 //________________________________________________________________ 
165 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
166 {
167   return (fEndTime - fStartTime)*fkSecToHour*fTimeBinsPerHour;
168 }
169
170 //________________________________________________________________ 
171 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
172 {
173   return (fMaxTime - fMinTime)*fkSecToHour;
174 }
175
176 //________________________________________________________________ 
177 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
178 {
179   return (fMaxTemp - fMinTemp);
180 }
181
182 //________________________________________________________________
183 void AliEMCALCalibTimeDep::Initialize(Int_t run, 
184                                       UInt_t startTime, UInt_t endTime)
185 {
186   Reset(); // start fresh
187
188   fRun = run;
189   fStartTime = startTime;
190   fEndTime = endTime;
191   
192   // collect the needed information
193   GetTemperatureInfo(); // temperature readings during the run
194   ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
195
196   return;
197 }
198
199 //________________________________________________________________
200 Double_t AliEMCALCalibTimeDep::GetTemperature(UInt_t timeStamp) const
201 {// return estimate for all SuperModules and sensors, that had data 
202
203   // first convert from seconds to hours..
204   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
205
206   Double_t average = 0;
207   int n = 0;
208
209   for (int i=0; i<fTempArray->NumSensors(); i++) {
210     
211     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
212
213     // check if we had valid data for the time that is being asked for
214     if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
215       AliSplineFit *f = st->GetFit();
216       if (f) { // ok, looks like we have valid data/info
217         // let's check what the expected value at the time appears to be
218         Double_t val = f->Eval(timeHour);
219         average += val;
220         n++;
221       }
222     } // time
223   } // loop over fTempArray
224   
225   if (n>0) { // some valid data was found
226     average /= n;
227     return average;
228   }
229   else { // no good data
230     return fkErrorCode;
231   }
232
233 }
234
235 //________________________________________________________________
236 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
237 {// return estimate for this one SuperModule, if it had data 
238
239   // first convert from seconds to hours..
240   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
241
242   Double_t average = 0;
243   int n = 0;
244
245   for (int i=0; i<fTempArray->NumSensors(); i++) {
246     
247     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
248     int module = st->GetSector()*2 + st->GetSide();
249     if ( module == imod ) { // right module
250       // check if we had valid data for the time that is being asked for
251       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
252         AliSplineFit *f = st->GetFit();
253         if (f) { // ok, looks like we have valid data/info
254           // let's check what the expected value at the time appears to be
255           Double_t val = f->Eval(timeHour);
256           cout << " i " << i << " val " << val << endl;
257           average += val;
258           n++;
259         }
260       } // time
261     }
262     
263   } // loop over fTempArray
264   
265   if (n>0) { // some valid data was found
266     average /= n;
267     return average;
268   }
269   else { // no good data
270     return fkErrorCode;
271   }
272
273 }
274
275 //________________________________________________________________
276 Double_t AliEMCALCalibTimeDep::GetTemperatureSMSensor(int imod, int isens, UInt_t timeStamp) const
277 {// return estimate for this one SuperModule and sensor, if it had data 
278
279   // first convert from seconds to hours..
280   Double_t timeHour = (timeStamp - fStartTime) * fkSecToHour;
281
282   for (int i=0; i<fTempArray->NumSensors(); i++) {
283     
284     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
285     int module = st->GetSector()*2 + st->GetSide();
286     if ( module == imod && st->GetNum()==isens ) { // right module, and sensor
287       // check if we had valid data for the time that is being asked for
288       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
289         AliSplineFit *f = st->GetFit();
290         if (f) { // ok, looks like we have valid data/info
291           // let's check what the expected value at the time appears to be
292           Double_t val = f->Eval(timeHour);
293
294           return val; // no point to move further in for loop, we have found the sensor we were looking for
295         }
296       } // time
297     }
298     
299   } // loop over fTempArray
300   
301   // if we made it all here, it means that we didn't find the sensor we were looking for
302   // i.e. no good data
303   return fkErrorCode;
304
305 }
306
307 //________________________________________________________________
308 Int_t AliEMCALCalibTimeDep::CalcCorrection()
309 { // OK, this is where the real action takes place - the heart of this class..
310   /* The philosophy is as follows:
311      0. Init corrections to 1.0 values
312      1: if we have LED info for the tower, use it
313      2. if not 1, we rely on LED info averaged over strip
314      3. if not 2 either, we try to use temperature info + APD bias and calibration info
315    */
316
317   // 0: Init
318   // how many SuperModules do we have?
319   Int_t nSM = fCalibAbs->GetNSuperModule();
320   // how many time-bins should we have for this run?
321   Int_t nBins = (Int_t) GetLengthOfRunInBins(); // round-down (from double to int)
322   // set up a reasonable default (correction = 1.0)
323   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
324
325   // 1+2: try with LED corrections
326   Int_t nRemaining = CalcLEDCorrection(nSM, nBins);
327
328   // 3: try with Temperature, if needed
329   if (nRemaining>0) {
330     nRemaining = CalcTemperatureCorrection(nSM, nBins);
331   }
332
333   return nRemaining;
334 }
335
336
337 //________________________________________________________________
338 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
339 { // estimate the Temperature Coefficient, based on the dark current (IDark)
340   // and the gain (M), based on Catania parameterizations
341
342   Double_t P0 = fkTempCoeffP0Const + fkTempCoeffP0Factor * IDark;
343   Double_t P1 = fkTempCoeffP1Const + fkTempCoeffP1Factor * IDark;
344  
345   Double_t TC = P0 + P1*M;
346
347   return TC;
348 }
349
350 /* Next come the methods that do the work in picking up all the needed info..*/
351 //________________________________________________________________
352 void AliEMCALCalibTimeDep::GetTemperatureInfo() 
353 {
354   // pick up Preprocessor output, based on fRun (most recent version)
355   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
356   if (entry) {
357     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
358   }
359
360   if (fTempArray) { 
361     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
362                   fTempArray->NumSensors(),
363                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
364   }
365   else {
366     AliWarning( Form("AliEMCALSensorTempArray not found!") );
367   }
368   
369   return;
370 }
371
372 //________________________________________________________________
373 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() 
374 {// assign max/min time and temperature values
375
376   fMinTemp = 999; // init to some large value (999 deg C)
377   fMaxTemp = 0;
378   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
379   fMaxTime = 0;
380
381   Int_t n = 0; // number of valid readings
382
383   for (int i=0; i<fTempArray->NumSensors(); i++) {
384     
385     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
386
387     // check time ranges
388     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
389     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
390
391     // check temperature ranges
392     TGraph *g = st->GetGraph();
393     if (g) { // ok, looks like we have valid data/info
394       // let's check what the expected value at the time appears to be
395       if (fMinTemp > g->GetMinimum()) { fMinTemp = g->GetMinimum(); }
396       if (fMaxTemp < g->GetMaximum()) { fMaxTemp = g->GetMaximum(); }
397       n++;
398     }
399   } // loop over fTempArray
400   
401   if (n>0) { // some valid data was found
402     return n;
403   }
404   else { // no good data
405     return (Int_t) fkErrorCode;
406   }
407
408 }
409
410 //________________________________________________________________
411 void AliEMCALCalibTimeDep::GetCalibSignalInfo() 
412 {
413   // pick up Preprocessor output, based on fRun (most recent version)
414   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
415   if (entry) {
416     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
417   }
418
419   if (fCalibSignal) { 
420     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %d AvgEntries LEDRefEntries %d LEDRefAvgEntries %d",
421                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
422                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
423                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
424                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
425                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );               
426   }
427   else {
428     AliWarning( Form("AliCaloCalibSignal not found!") );
429   }
430   
431   return;
432 }
433
434 //________________________________________________________________
435 void AliEMCALCalibTimeDep::GetBiasAPDInfo() 
436 {
437   // pick up Preprocessor output, based on fRun (most recent version)
438   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/BiasAPD", fRun);
439   if (entry) {
440     fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
441   }
442
443   if (fBiasAPD) { 
444     AliInfo( Form("BiasAPD: NSuperModule %d ", fBiasAPD->GetNSuperModule() ) );
445   }
446   else {
447     AliWarning( Form("AliEMCALBiasAPD not found!") );
448   }
449   
450   return;
451 }
452
453 //________________________________________________________________
454 void AliEMCALCalibTimeDep::GetCalibMapAPDInfo() 
455 {
456   // pick up Preprocessor output, based on fRun (most recent version)
457   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
458   if (entry) {
459     fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
460   }
461
462   if (fCalibMapAPD) { 
463     AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
464   }
465   else {
466     AliWarning( Form("AliEMCALCalibMapAPD not found!") );
467   }
468   
469   return;
470 }
471
472 //________________________________________________________________
473 void AliEMCALCalibTimeDep::GetCalibAbsInfo() 
474 {
475   // pick up Preprocessor output, based on fRun (most recent version)
476   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
477   if (entry) {
478     fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject();
479   }
480
481   if (fCalibAbs) { 
482     AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
483   }
484   else {
485     AliWarning( Form("AliEMCALCalibAbs not found!") );
486   }
487   
488   return;
489 }
490
491 //________________________________________________________________
492 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) 
493 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0 
494   // The correction factor we keep is c(T) = R(t0)/R(T)
495   // T info from fCalibSignal, t0 info from fCalibAbs
496
497   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
498   // but one could upgrade this in the future
499   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
500
501   // sanity check; same SuperModule indices for corrections as for regular calibrations
502   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
503   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
504
505   for (int i = 0; i < nSM; i++) {
506     int iSMAbs = CalibAbsData[i].fSuperModuleNum;
507     int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
508     if (iSMAbs != iSMCorr) {
509       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
510       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
511       return nRemaining;
512     }
513   }
514  
515   int iSM = 0;
516   int iCol = 0;
517   int iRow = 0;
518   int iStrip = 0;
519   int iGain = 0;
520
521   // The fCalibSignal info is stored in TTrees
522   // Note that the time-bins for the TTree's may not exactly match with our correction time bins 
523   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
524   // fCalibSignal time info in seconds: Hour/fkSecToHour
525   // corrected for startTime difference: Hour/fkSecToHour + timeDiff
526   // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
527
528   // values for R(T), size of TArray = nBins
529   // the [2] dimension below is for the low or high gain 
530   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
531   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
532   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
533   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
534
535   // set up TArray's first
536   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
537     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
538       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
539         for (iGain = 0; iGain < 2; iGain++) {
540           // length of arrays
541           ampT[iSM][iCol][iRow][iGain].Set(nBins);
542           nT[iSM][iCol][iRow][iGain].Set(nBins);
543           // content of arrys
544           for (int j = 0; j < nBins; j++) {
545             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
546             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
547           }
548         }
549       }
550     }//iCol
551     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
552       for (iGain = 0; iGain < 2; iGain++) {
553         // length of arrays
554         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
555         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
556         // content of arrys
557         for (int j = 0; j < nBins; j++) {
558           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
559           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
560         }
561       }
562     }//iStrip
563   }
564
565   // OK, now loop over the TTrees and fill the arrays needed for R(T)
566   TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
567   TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
568
569   int ChannelNum; // for regular towers
570   int RefNum; // for LED
571   double Hour;
572   double AvgAmp;
573
574   TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
575   TAvg->SetBranchAddress("fHour", &Hour);
576   TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
577
578   int iBin = 0;
579   // counters for how many values were seen per SuperModule
580   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
581   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
582
583   for (int ient=0; ient<TAvg->GetEntries(); ient++) {
584     TAvg->GetEntry(ient);
585     // figure out where this info comes from
586     fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
587     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
588     // add value in the arrays
589     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
590     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
591     nCount[iSM]++;
592   }
593
594   TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
595   TLEDRefAvg->SetBranchAddress("fHour", &Hour);
596   TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
597
598   for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
599     TLEDRefAvg->GetEntry(ient);
600     // figure out where this info comes from
601     fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
602     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
603     // add value in the arrays
604     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
605     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
606     nCountLEDRef[iSM]++;
607   }
608
609   // Normalize TArray values, and calculate average also
610   Float_t norm = 0; // extra var, for readability
611
612   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
613     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
614       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
615         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
616         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
617         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
618           for (iGain = 0; iGain < 2; iGain++) {
619             // content of arrys
620             for (int j = 0; j < nBins; j++) {
621               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) { 
622                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j); 
623                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
624               } 
625             }
626           }
627         }
628       }//iCol
629       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
630         for (iGain = 0; iGain < 2; iGain++) {
631           for (int j = 0; j < nBins; j++) {
632             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
633               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j); 
634               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
635             }
636           }
637         }
638       }//iStrip
639     }
640   } // iSM
641
642
643   // Calculate correction values, and store them
644   // set fkErrorCode values for those that could not be set
645
646   Float_t Rt0 = 0;
647   Float_t RT = 0;
648   Float_t correction = 0; // c(T) = R(t0)/R(T)
649   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
650
651   for (int i = 0; i < nSM; i++) {
652     iSM = CalibAbsData[i].fSuperModuleNum;
653     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
654       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
655       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
656       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
657
658         // Calc. R(t0):
659         AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
660         iGain = v.fHighLow; // gain value used for abs. calibration     
661         refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration 
662
663         // valid amplitude values should be larger than 0
664         if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
665           Rt0 =  v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
666         }
667         else {
668           Rt0 = fkErrorCode;
669         }
670
671         // Cal R(T)
672         for (int j = 0; j < nBins; j++) {
673
674           // calculate R(T) also; first try with individual tower:
675           // same gain as for abs. calibration is the default
676           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
677             // looks like valid data with the right gain combination
678             RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); 
679
680             // if data appears to be saturated, and we are in high gain, then try with low gain instead
681             if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
682                  (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
683               RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j); 
684               RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
685             }
686           }
687           else {
688             RT = fkErrorCode;
689           }
690
691           // Calc. correction factor
692           if (Rt0>0 && RT>0) { 
693             correction = Rt0/RT;
694           }
695           else {
696             correction = fkErrorCode;
697             nRemaining++;
698           }
699
700           // Store the value
701           CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
702           /* Check that
703           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
704           also works OK */
705         } // nBins
706       }
707     }
708   }
709
710   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
711   return nRemaining;
712 }
713
714 //________________________________________________________________
715 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) 
716 { // use averages for each strip if no good values exist for some single tower 
717
718   // go over fTimeDepCorrection info
719   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
720
721   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
722
723   // for calculating StripAverage info
724   int nValidTower = 0;
725   Float_t StripAverage = 0;
726   Float_t val = 0;
727
728   int iSM = 0;
729   int iCol = 0;
730   int iRow = 0;
731   int iStrip = 0;
732   int firstCol = 0;
733   int lastCol = 0;
734
735   for (int i = 0; i < nSM; i++) {
736     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
737     for (int j = 0; j < nBins; j++) {
738
739       nValidTower = 0;
740       StripAverage = 0;
741
742       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
743         firstCol = iStrip*2;
744         if ((iSM%2)==1) { // C side
745           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
746         }
747         lastCol = firstCol+1;
748
749         for (iCol = firstCol; iCol <= lastCol; iCol++) {
750           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
751             val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
752             if (val>0) { // valid value; error code is negative
753               StripAverage += val;
754               nValidTower++;
755             } 
756           }
757         }
758
759         // calc average over strip
760         if (nValidTower>0) {
761           StripAverage /= nValidTower;
762           for (iCol = firstCol; iCol <= lastCol; iCol++) {
763             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
764               val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
765               if (val<0) { // invalid value; error code is negative
766                 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
767               }
768             }
769           }
770         }
771         else { // could not fill in unvalid entries
772           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
773         }
774
775       } // iStrip
776     } // j, bins
777   } // iSM
778
779   return nRemaining;
780 }
781
782 //________________________________________________________________
783 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins) 
784 { // OK, so we didn't have valid LED data that allowed us to do the correction only 
785   // with that info.
786   // So, instead we'll rely on the temperature info and try to do the correction
787   // based on that instead.
788   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class 
789   Int_t nRemaining = 0;
790
791   // info containers
792   AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
793   AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
794   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
795   // correction container
796   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
797
798   int iSM = 0;
799   int iCol = 0;
800   int iRow = 0;
801
802   Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
803   memset(TempCoeff, 0, sizeof(TempCoeff));
804   Float_t MGain = 0; 
805   Double_t correction = 0;
806   Double_t secondsPerBin = (1/fkSecToHour)*(1/fTimeBinsPerHour);
807
808   for (int i = 0; i < nSM; i++) {
809     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
810     
811     // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
812     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
813       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
814         AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
815         MGain =  fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2], 
816                                        BiasAPDData[i].fVoltage[iCol][iRow]);
817         TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
818       }
819     }
820     
821     // figure out what the reference temperature is, from fCalibAbs
822     Double_t ReferenceTemperature = 0;
823     int nVal = 0;
824     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
825       if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
826         ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
827         nVal++;
828       }
829     }
830     
831     if (nVal>0) {
832       ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
833       
834       for (int j = 0; j < nBins; j++) {
835
836         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
837         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
838       // 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)
839         Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp); 
840
841         Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
842         // if the new temperature is higher than the old/reference one, then the gain has gone down 
843         if (fabs(TemperatureDiff)>fTemperatureResolution) { 
844           // significant enough difference that we need to consider it
845
846           // loop over all towers; effect of temperature change will depend on gain for this tower
847           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
848             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
849
850               correction = TemperatureDiff * TempCoeff[iCol][iRow];
851               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
852             }
853           }
854
855         } // if noteworthy temperature change
856         else { // just set correction values to 1.0
857           correction = 1;
858           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
859             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
860               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
861             }
862           }
863         } // else
864       } // j, Bins
865
866     } // if reference temperature exist 
867     else { // could not do the needed check.. signal that in the return code
868       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
869     }
870   } // iSM
871
872   return nRemaining;
873 }
874