added a few more variables to AliEMCALCalibTimeDepCorrection to make the info there...
[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   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
323   // set up a reasonable default (correction = 1.0)
324   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
325   fCalibTimeDepCorrection->SetStartTime(fStartTime);
326   fCalibTimeDepCorrection->SetNTimeBins(nBins);
327   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
328
329   // 1+2: try with LED corrections
330   Int_t nRemaining = CalcLEDCorrection(nSM, nBins);
331
332   // 3: try with Temperature, if needed
333   if (nRemaining>0) {
334     nRemaining = CalcTemperatureCorrection(nSM, nBins);
335   }
336
337   return nRemaining;
338 }
339
340
341 //________________________________________________________________
342 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
343 { // estimate the Temperature Coefficient, based on the dark current (IDark)
344   // and the gain (M), based on Catania parameterizations
345
346   Double_t P0 = fkTempCoeffP0Const + fkTempCoeffP0Factor * IDark;
347   Double_t P1 = fkTempCoeffP1Const + fkTempCoeffP1Factor * IDark;
348  
349   Double_t TC = P0 + P1*M;
350
351   return TC;
352 }
353
354 /* Next come the methods that do the work in picking up all the needed info..*/
355 //________________________________________________________________
356 void AliEMCALCalibTimeDep::GetTemperatureInfo() 
357 {
358   // pick up Preprocessor output, based on fRun (most recent version)
359   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
360   if (entry) {
361     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
362   }
363
364   if (fTempArray) { 
365     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
366                   fTempArray->NumSensors(),
367                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
368   }
369   else {
370     AliWarning( Form("AliEMCALSensorTempArray not found!") );
371   }
372   
373   return;
374 }
375
376 //________________________________________________________________
377 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() 
378 {// assign max/min time and temperature values
379
380   fMinTemp = 999; // init to some large value (999 deg C)
381   fMaxTemp = 0;
382   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
383   fMaxTime = 0;
384
385   Int_t n = 0; // number of valid readings
386
387   for (int i=0; i<fTempArray->NumSensors(); i++) {
388     
389     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
390
391     // check time ranges
392     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
393     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
394
395     // check temperature ranges
396     TGraph *g = st->GetGraph();
397     if (g) { // ok, looks like we have valid data/info
398       // let's check what the expected value at the time appears to be
399       if (fMinTemp > g->GetMinimum()) { fMinTemp = g->GetMinimum(); }
400       if (fMaxTemp < g->GetMaximum()) { fMaxTemp = g->GetMaximum(); }
401       n++;
402     }
403   } // loop over fTempArray
404   
405   if (n>0) { // some valid data was found
406     return n;
407   }
408   else { // no good data
409     return (Int_t) fkErrorCode;
410   }
411
412 }
413
414 //________________________________________________________________
415 void AliEMCALCalibTimeDep::GetCalibSignalInfo() 
416 {
417   // pick up Preprocessor output, based on fRun (most recent version)
418   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
419   if (entry) {
420     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
421   }
422
423   if (fCalibSignal) { 
424     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %d AvgEntries LEDRefEntries %d LEDRefAvgEntries %d",
425                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
426                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
427                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
428                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
429                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );               
430   }
431   else {
432     AliWarning( Form("AliCaloCalibSignal not found!") );
433   }
434   
435   return;
436 }
437
438 //________________________________________________________________
439 void AliEMCALCalibTimeDep::GetBiasAPDInfo() 
440 {
441   // pick up Preprocessor output, based on fRun (most recent version)
442   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/BiasAPD", fRun);
443   if (entry) {
444     fBiasAPD = (AliEMCALBiasAPD *) entry->GetObject();
445   }
446
447   if (fBiasAPD) { 
448     AliInfo( Form("BiasAPD: NSuperModule %d ", fBiasAPD->GetNSuperModule() ) );
449   }
450   else {
451     AliWarning( Form("AliEMCALBiasAPD not found!") );
452   }
453   
454   return;
455 }
456
457 //________________________________________________________________
458 void AliEMCALCalibTimeDep::GetCalibMapAPDInfo() 
459 {
460   // pick up Preprocessor output, based on fRun (most recent version)
461   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
462   if (entry) {
463     fCalibMapAPD = (AliEMCALCalibMapAPD *) entry->GetObject();
464   }
465
466   if (fCalibMapAPD) { 
467     AliInfo( Form("CalibMapAPD: NSuperModule %d ", fCalibMapAPD->GetNSuperModule() ) );
468   }
469   else {
470     AliWarning( Form("AliEMCALCalibMapAPD not found!") );
471   }
472   
473   return;
474 }
475
476 //________________________________________________________________
477 void AliEMCALCalibTimeDep::GetCalibAbsInfo() 
478 {
479   // pick up Preprocessor output, based on fRun (most recent version)
480   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/MapAPD", fRun);
481   if (entry) {
482     fCalibAbs = (AliEMCALCalibAbs *) entry->GetObject();
483   }
484
485   if (fCalibAbs) { 
486     AliInfo( Form("CalibAbs: NSuperModule %d ", fCalibAbs->GetNSuperModule() ) );
487   }
488   else {
489     AliWarning( Form("AliEMCALCalibAbs not found!") );
490   }
491   
492   return;
493 }
494
495 //________________________________________________________________
496 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) 
497 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0 
498   // The correction factor we keep is c(T) = R(t0)/R(T)
499   // T info from fCalibSignal, t0 info from fCalibAbs
500
501   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibAbs
502   // but one could upgrade this in the future
503   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
504
505   // sanity check; same SuperModule indices for corrections as for regular calibrations
506   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
507   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
508
509   for (int i = 0; i < nSM; i++) {
510     int iSMAbs = CalibAbsData[i].fSuperModuleNum;
511     int iSMCorr = CalibTimeDepCorrectionData[i].fSuperModuleNum;
512     if (iSMAbs != iSMCorr) {
513       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMAbs, iSMCorr) );
514       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
515       return nRemaining;
516     }
517   }
518  
519   int iSM = 0;
520   int iCol = 0;
521   int iRow = 0;
522   int iStrip = 0;
523   int iGain = 0;
524
525   // The fCalibSignal info is stored in TTrees
526   // Note that the time-bins for the TTree's may not exactly match with our correction time bins 
527   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
528   // fCalibSignal time info in seconds: Hour/fkSecToHour
529   // corrected for startTime difference: Hour/fkSecToHour + timeDiff
530   // converted into a time-bin we use: (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour
531
532   // values for R(T), size of TArray = nBins
533   // the [2] dimension below is for the low or high gain 
534   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
535   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
536   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
537   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
538
539   // set up TArray's first
540   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
541     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
542       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
543         for (iGain = 0; iGain < 2; iGain++) {
544           // length of arrays
545           ampT[iSM][iCol][iRow][iGain].Set(nBins);
546           nT[iSM][iCol][iRow][iGain].Set(nBins);
547           // content of arrys
548           for (int j = 0; j < nBins; j++) {
549             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
550             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
551           }
552         }
553       }
554     }//iCol
555     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
556       for (iGain = 0; iGain < 2; iGain++) {
557         // length of arrays
558         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
559         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
560         // content of arrys
561         for (int j = 0; j < nBins; j++) {
562           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
563           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
564         }
565       }
566     }//iStrip
567   }
568
569   // OK, now loop over the TTrees and fill the arrays needed for R(T)
570   TTree *TAvg = fCalibSignal->GetTreeAvgAmpVsTime();
571   TTree *TLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
572
573   int ChannelNum; // for regular towers
574   int RefNum; // for LED
575   double Hour;
576   double AvgAmp;
577
578   TAvg->SetBranchAddress("fChannelNum", &ChannelNum);
579   TAvg->SetBranchAddress("fHour", &Hour);
580   TAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
581
582   int iBin = 0;
583   // counters for how many values were seen per SuperModule
584   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
585   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
586
587   for (int ient=0; ient<TAvg->GetEntries(); ient++) {
588     TAvg->GetEntry(ient);
589     // figure out where this info comes from
590     fCalibSignal->DecodeChannelNum(ChannelNum, &iSM, &iCol, &iRow, &iGain);
591     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
592     // add value in the arrays
593     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+AvgAmp, iBin );
594     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
595     nCount[iSM]++;
596   }
597
598   TLEDRefAvg->SetBranchAddress("fRefNum", &RefNum);
599   TLEDRefAvg->SetBranchAddress("fHour", &Hour);
600   TLEDRefAvg->SetBranchAddress("fAvgAmp",&AvgAmp);
601
602   for (int ient=0; ient<TLEDRefAvg->GetEntries(); ient++) {
603     TLEDRefAvg->GetEntry(ient);
604     // figure out where this info comes from
605     fCalibSignal->DecodeRefNum(RefNum, &iSM, &iStrip, &iGain);
606     iBin = (int) ( (Hour + timeDiff*fkSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
607     // add value in the arrays
608     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+AvgAmp, iBin );
609     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
610     nCountLEDRef[iSM]++;
611   }
612
613   // Normalize TArray values, and calculate average also
614   Float_t norm = 0; // extra var, for readability
615
616   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
617     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
618       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
619         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
620         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
621         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
622           for (iGain = 0; iGain < 2; iGain++) {
623             // content of arrys
624             for (int j = 0; j < nBins; j++) {
625               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) { 
626                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j); 
627                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
628               } 
629             }
630           }
631         }
632       }//iCol
633       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
634         for (iGain = 0; iGain < 2; iGain++) {
635           for (int j = 0; j < nBins; j++) {
636             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
637               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j); 
638               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
639             }
640           }
641         }
642       }//iStrip
643     }
644   } // iSM
645
646
647   // Calculate correction values, and store them
648   // set fkErrorCode values for those that could not be set
649
650   Float_t Rt0 = 0;
651   Float_t RT = 0;
652   Float_t correction = 0; // c(T) = R(t0)/R(T)
653   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
654
655   for (int i = 0; i < nSM; i++) {
656     iSM = CalibAbsData[i].fSuperModuleNum;
657     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
658       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
659       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
660       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
661
662         // Calc. R(t0):
663         AliEMCALCalibAbs::AliEMCALCalibAbsVal &v = CalibAbsData[i].fAPDVal[iCol][iRow];
664         iGain = v.fHighLow; // gain value used for abs. calibration     
665         refGain = CalibAbsData[i].fLEDRefHighLow[iStrip]; // LED reference gain value used for abs. calibration 
666
667         // valid amplitude values should be larger than 0
668         if (v.fLEDAmp>0 && CalibAbsData[i].fLEDRefAmp[iStrip]>0) {
669           Rt0 =  v.fLEDAmp / CalibAbsData[i].fLEDRefAmp[iStrip];
670         }
671         else {
672           Rt0 = fkErrorCode;
673         }
674
675         // Cal R(T)
676         for (int j = 0; j < nBins; j++) {
677
678           // calculate R(T) also; first try with individual tower:
679           // same gain as for abs. calibration is the default
680           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
681             // looks like valid data with the right gain combination
682             RT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); 
683
684             // if data appears to be saturated, and we are in high gain, then try with low gain instead
685             if ( (ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1) ||
686                  (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
687               RT = ampT[iSM][iCol][iRow][0].At(j) / ampLEDRefT[iSM][iStrip][0].At(j); 
688               RT *= v.fHighLowRatio/CalibAbsData[i].fLEDRefHighLowRatio[iStrip]; // compensate for using different gain than in the absolute calibration
689             }
690           }
691           else {
692             RT = fkErrorCode;
693           }
694
695           // Calc. correction factor
696           if (Rt0>0 && RT>0) { 
697             correction = Rt0/RT;
698           }
699           else {
700             correction = fkErrorCode;
701             nRemaining++;
702           }
703
704           // Store the value
705           CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
706           /* Check that
707           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
708           also works OK */
709         } // nBins
710       }
711     }
712   }
713
714   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
715   return nRemaining;
716 }
717
718 //________________________________________________________________
719 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) 
720 { // use averages for each strip if no good values exist for some single tower 
721
722   // go over fTimeDepCorrection info
723   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
724
725   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
726
727   // for calculating StripAverage info
728   int nValidTower = 0;
729   Float_t StripAverage = 0;
730   Float_t val = 0;
731
732   int iSM = 0;
733   int iCol = 0;
734   int iRow = 0;
735   int iStrip = 0;
736   int firstCol = 0;
737   int lastCol = 0;
738
739   for (int i = 0; i < nSM; i++) {
740     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
741     for (int j = 0; j < nBins; j++) {
742
743       nValidTower = 0;
744       StripAverage = 0;
745
746       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
747         firstCol = iStrip*2;
748         if ((iSM%2)==1) { // C side
749           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
750         }
751         lastCol = firstCol+1;
752
753         for (iCol = firstCol; iCol <= lastCol; iCol++) {
754           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
755             val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
756             if (val>0) { // valid value; error code is negative
757               StripAverage += val;
758               nValidTower++;
759             } 
760           }
761         }
762
763         // calc average over strip
764         if (nValidTower>0) {
765           StripAverage /= nValidTower;
766           for (iCol = firstCol; iCol <= lastCol; iCol++) {
767             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
768               val = CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].At(j);
769               if (val<0) { // invalid value; error code is negative
770                 CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(val, j);
771               }
772             }
773           }
774         }
775         else { // could not fill in unvalid entries
776           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
777         }
778
779       } // iStrip
780     } // j, bins
781   } // iSM
782
783   return nRemaining;
784 }
785
786 //________________________________________________________________
787 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins) 
788 { // OK, so we didn't have valid LED data that allowed us to do the correction only 
789   // with that info.
790   // So, instead we'll rely on the temperature info and try to do the correction
791   // based on that instead.
792   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class 
793   Int_t nRemaining = 0;
794
795   // info containers
796   AliEMCALBiasAPD::AliEMCALSuperModuleBiasAPD * BiasAPDData = fBiasAPD->GetSuperModuleData();
797   AliEMCALCalibMapAPD::AliEMCALSuperModuleCalibMapAPD * CalibMapAPDData = fCalibMapAPD->GetSuperModuleData();
798   AliEMCALCalibAbs::AliEMCALSuperModuleCalibAbs * CalibAbsData = fCalibAbs->GetSuperModuleData();
799   // correction container
800   AliEMCALCalibTimeDepCorrection::AliEMCALSuperModuleCalibTimeDepCorrection * CalibTimeDepCorrectionData = fCalibTimeDepCorrection->GetSuperModuleData();
801
802   int iSM = 0;
803   int iCol = 0;
804   int iRow = 0;
805
806   Double_t TempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
807   memset(TempCoeff, 0, sizeof(TempCoeff));
808   Float_t MGain = 0; 
809   Double_t correction = 0;
810   Double_t secondsPerBin = (3600/fTimeBinsPerHour);
811
812   for (int i = 0; i < nSM; i++) {
813     iSM = CalibTimeDepCorrectionData[i].fSuperModuleNum;
814     
815     // first calculate the M=Gain values, and TemperatureCoeff, for all towers in this SuperModule, from BiasAPD and CalibMapAPD info
816     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
817       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
818         AliEMCALCalibMapAPD::AliEMCALCalibMapAPDVal &mapAPD = CalibMapAPDData[i].fAPDVal[iCol][iRow];
819         MGain =  fCalibMapAPD->GetGain(mapAPD.fPar[0], mapAPD.fPar[1], mapAPD.fPar[2], 
820                                        BiasAPDData[i].fVoltage[iCol][iRow]);
821         TempCoeff[iCol][iRow] = GetTempCoeff(mapAPD.fDarkCurrent, MGain);
822       }
823     }
824     
825     // figure out what the reference temperature is, from fCalibAbs
826     Double_t ReferenceTemperature = 0;
827     int nVal = 0;
828     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
829       if (CalibAbsData[i].fTemperature[iSensor]>0) { // hopefully OK value
830         ReferenceTemperature += CalibAbsData[i].fTemperature[iSensor];
831         nVal++;
832       }
833     }
834     
835     if (nVal>0) {
836       ReferenceTemperature /= nVal; // valid values exist, we can look into corrections
837       
838       for (int j = 0; j < nBins; j++) {
839
840         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
841         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
842       // 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)
843         Double_t SMTemperature = GetTemperatureSM(iSM, timeStamp); 
844
845         Double_t TemperatureDiff = ReferenceTemperature - SMTemperature; // old vs new
846         // if the new temperature is higher than the old/reference one, then the gain has gone down 
847         if (fabs(TemperatureDiff)>fTemperatureResolution) { 
848           // significant enough difference that we need to consider it
849
850           // loop over all towers; effect of temperature change will depend on gain for this tower
851           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
852             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
853
854               correction = TemperatureDiff * TempCoeff[iCol][iRow];
855               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
856             }
857           }
858
859         } // if noteworthy temperature change
860         else { // just set correction values to 1.0
861           correction = 1;
862           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
863             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
864               CalibTimeDepCorrectionData[i].fCorrection[iCol][iRow].AddAt(correction, j);
865             }
866           }
867         } // else
868       } // j, Bins
869
870     } // if reference temperature exist 
871     else { // could not do the needed check.. signal that in the return code
872       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
873     }
874   } // iSM
875
876   return nRemaining;
877 }
878