remove exec
[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 kTempMaxDiffMedian = 2; // Temperature values should not be further away from median value within SM when considered in the average calc.
57
58 const double kErrorCode = -999; // to indicate that something went wrong
59
60 using namespace std;
61
62 ClassImp(AliEMCALCalibTimeDep)
63
64 //________________________________________________________________
65 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
66   fRun(0),
67   fStartTime(0),
68   fEndTime(0),
69   fMinTemp(0),
70   fMaxTemp(0),
71   fMinTempVariation(0),
72   fMaxTempVariation(0),
73   fMinTempValid(15),
74   fMaxTempValid(35),
75   fMinTime(0),
76   fMaxTime(0),
77   fTemperatureResolution(0.1), // 0.1 deg C is default
78   fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference 
79   fTimeBinsPerHour(2), // 2 30-min bins per hour is default
80   fHighLowGainFactor(16), // factor ~16 between High gain and low gain
81   fTempArray(NULL),
82   fCalibSignal(NULL),
83   fCalibTempCoeff(NULL),
84   fCalibReference(NULL),
85   fCalibTimeDepCorrection(NULL),
86   fVerbosity(0)
87 {
88   // Constructor
89 }
90
91 //________________________________________________________________
92 AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
93   TObject(calibt),
94   fRun(calibt.GetRunNumber()),
95   fStartTime(calibt.GetStartTime()),
96   fEndTime(calibt.GetEndTime()),
97   fMinTemp(calibt.GetMinTemp()),
98   fMaxTemp(calibt.GetMaxTemp()),
99   fMinTempVariation(calibt.GetMinTempVariation()),
100   fMaxTempVariation(calibt.GetMaxTempVariation()),
101   fMinTempValid(calibt.GetMinTempValid()),
102   fMaxTempValid(calibt.GetMaxTempValid()),
103   fMinTime(calibt.GetMinTime()),
104   fMaxTime(calibt.GetMaxTime()),
105   fTemperatureResolution(calibt.GetTemperatureResolution()),
106   fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()),
107   fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
108   fHighLowGainFactor(calibt.GetHighLowGainFactor()),
109   fTempArray(calibt.GetTempArray()),
110   fCalibSignal(calibt.GetCalibSignal()),
111   fCalibTempCoeff(calibt.GetCalibTempCoeff()),
112   fCalibReference(calibt.GetCalibReference()),
113   fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
114   fVerbosity(calibt.GetVerbosity())
115 {
116   // copy constructor
117 }
118
119
120 //________________________________________________________________
121 AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
122 {
123   // assignment operator; use copy ctor
124   if (&calibt == this) return *this;
125
126   new (this) AliEMCALCalibTimeDep(calibt);
127   return *this;
128 }
129
130 //________________________________________________________________
131 AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
132 {
133   // Destructor
134 }
135
136 //________________________________________________________________
137 void  AliEMCALCalibTimeDep::Reset() 
138 {
139   // clear variables to default
140   fRun = 0;
141   fStartTime = 0;
142   fEndTime = 0;
143   fMinTemp = 0;
144   fMaxTemp = 0;
145   fMinTempVariation = 0;
146   fMaxTempVariation = 0;
147   fMinTempValid = 15; 
148   fMaxTempValid = 35;
149   fMinTime = 0;
150   fMaxTime = 0;
151   fTemperatureResolution = 0.1; // 0.1 deg C is default
152   fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference 
153   fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
154   fTempArray = NULL;
155   fCalibSignal = NULL;
156   fCalibTempCoeff = NULL;
157   fCalibReference = NULL;
158   fCalibTimeDepCorrection = NULL;
159   fVerbosity = 0;
160   return;
161 }
162
163 //________________________________________________________________
164 void  AliEMCALCalibTimeDep::PrintInfo() const
165 {
166   // print some info
167   cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
168   // basic variables, all 'publicly available' also
169   cout << " VARIABLE DUMP: " << endl
170        << " GetStartTime() " << GetStartTime() << endl
171        << " GetEndTime() " << GetEndTime() << endl
172        << " GetMinTime() " << GetMinTime() << endl
173        << " GetMaxTime() " << GetMaxTime() << endl
174        << " GetMinTemp() " << GetMinTemp() << endl
175        << " GetMaxTemp() " << GetMaxTemp() << endl
176        << " GetMinTempVariation() " << GetMinTempVariation() << endl
177        << " GetMaxTempVariation() " << GetMaxTempVariation() << endl
178        << " GetTemperatureResolution() " << GetTemperatureResolution() << endl;
179   // run ranges
180   cout << " RUN INFO: " << endl
181        << " runnumber " << GetRunNumber() << endl
182        << " length (in hours) " << GetLengthOfRunInHours() << endl
183        << " length (in bins) " << GetLengthOfRunInBins() << endl
184        << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
185        << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
186        << endl;
187
188   return;
189 }
190
191 //________________________________________________________________ 
192 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
193 {
194   return (fEndTime - fStartTime)*kSecToHour;
195 }
196
197 //________________________________________________________________ 
198 Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
199 {
200   return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
201 }
202
203 //________________________________________________________________ 
204 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
205 {
206   return (fMaxTime - fMinTime)*kSecToHour;
207 }
208
209 //________________________________________________________________ 
210 Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
211 {
212   return (fMaxTemp - fMinTemp);
213 }
214
215 //________________________________________________________________
216 void AliEMCALCalibTimeDep::Initialize(Int_t run, 
217                                       UInt_t startTime, UInt_t endTime)
218 { // setup, and get temperature info
219   Reset(); // start fresh
220
221   fRun = run;
222   fStartTime = startTime;
223   fEndTime = endTime;
224   
225   // collect the needed information
226   GetTemperatureInfo(); // temperature readings during the run
227   ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
228
229   return;
230 }
231
232 //________________________________________________________________
233 Double_t AliEMCALCalibTimeDep::GetTemperatureSM(int imod, UInt_t timeStamp) const
234 {// return estimate for this one SuperModule, if it had data 
235
236   // first convert from seconds to hours..
237   Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
238
239   int n = 0;
240   Double_t valArr[8]={0}; // 8 sensors per SM
241
242   for (int i=0; i<fTempArray->NumSensors(); i++) {
243     
244     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
245     int module = st->GetSector()*2 + st->GetSide();
246     if ( module == imod ) { // right module
247       // check if we had valid data for the time that is being asked for
248       if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
249         AliSplineFit *f = st->GetFit();
250         if (f) { // ok, looks like we have valid data/info
251           // let's check what the expected value at the time appears to be
252           Double_t val = f->Eval(timeHour);
253           if ( fVerbosity > 0 ) {
254             cout << " sensor i " << i << " val " << val << endl;
255           }
256           if (val>fMinTempValid && val<fMaxTempValid && n<8) { 
257             valArr[n] = val;
258             n++;
259           }
260         }
261       } // time
262     }
263     
264   } // loop over fTempArray
265   
266   if (n>0) { // some valid data was found
267     Double_t median = TMath::Median(n, valArr);
268     Double_t average = 0;
269     Int_t nval = 0;
270     for (int is=0; is<n; is++) {
271       if (TMath::Abs(valArr[is] - median) < kTempMaxDiffMedian) {
272         average += valArr[is];
273         nval++;
274       }
275     }
276     if (nval >  0) {
277       average /= nval;
278       return average;
279     }
280     else { // this case should not happen, but kept for completeness (coverity etc)
281       return median;
282     }
283   }
284   else { // no good data
285     return kErrorCode;
286   }
287
288 }
289
290 //________________________________________________________________
291 Int_t AliEMCALCalibTimeDep::CalcCorrection()
292 { // OK, this is where the real action takes place - the heart of this class..
293   /* The philosophy is as follows:
294      0. Init corrections to 1.0 values, and see how many correction bins we need
295      1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
296      2. try to use temperature info + APD temperature coefficient info, to estimate correction.
297      For now (from Dec 2009), we do not use LED info.
298    */
299
300   // 0: Init
301   // how many SuperModules do we have?
302   Int_t nSM = fCalibReference->GetNSuperModule();
303   // how many time-bins should we have for this run?
304   Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
305   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
306
307   // 1: get info on how much individual sensors might have changed during
308   // the run (compare max-min for each sensor separately)
309   if (fMaxTempVariation < fTemperatureResolution) {
310     nBins = 1; // just one bin needed..
311   }
312   if (nBins == 1) {    
313     binSize = fEndTime - fStartTime;
314   }
315   if (fVerbosity > 0) {
316     cout << " nBins " << nBins << " binSize " << binSize << endl;
317   }
318
319   // set up a reasonable default (correction = 1.0)
320   fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
321   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
322   fCalibTimeDepCorrection->SetStartTime(fStartTime);
323   fCalibTimeDepCorrection->SetNTimeBins(nBins);
324   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
325
326   // 2: try with Temperature correction 
327   Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
328
329   return nRemaining;
330 }
331
332
333 //________________________________________________________________
334 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
335 { // estimate the Temperature Coefficient, based on the dark current (IDark)
336   // and the gain (M), based on Catania parameterizations
337
338   Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
339   Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
340  
341   Double_t dTC = dP0 + dP1*M;
342   // from % numbers to regular ones..:
343   dTC *= 0.01;
344
345   return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
346 }
347
348 /* Next come the methods that do the work in picking up all the needed info..*/
349 //________________________________________________________________
350 void AliEMCALCalibTimeDep::GetTemperatureInfo() 
351 {
352   // pick up Preprocessor output, based on fRun (most recent version)
353   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
354   if (entry) {
355     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
356   }
357
358   if (fTempArray) { 
359     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
360                   fTempArray->NumSensors(),
361                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
362   }
363   else {
364     AliWarning( Form("AliEMCALSensorTempArray not found!") );
365   }
366   
367   return;
368 }
369
370 //________________________________________________________________
371 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() 
372 {// assign max/min time and temperature values
373
374   fMinTemp = 999; // init to some large value (999 deg C)
375   fMaxTemp = 0;
376   fMinTempVariation = 999; // init to some large value (999 deg C)
377   fMaxTempVariation = 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     if ( st->GetStartTime() == 0 ) { // no valid data
387       continue;
388     }
389
390     // check time ranges
391     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
392     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
393  
394     // check temperature ranges
395     AliSplineFit *f = st->GetFit();
396
397     if (f) { // ok, looks like we have valid data/info
398       int np = f->GetKnots();
399       Double_t *y0 = f->GetY0();
400       // min and max values within the single sensor
401       Double_t min = 999;
402       Double_t max = 0;
403       int nval = 0;
404       for (int ip=0; ip<np; ip++) { 
405         if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
406           if (min > y0[ip]) { min = y0[ip]; }
407           if (max < y0[ip]) { max = y0[ip]; }
408           nval++;
409         }
410       }
411       if (nval>0) {
412         if (fMinTemp > min) { fMinTemp = min; }
413         if (fMaxTemp < max) { fMaxTemp = max; }
414         Double_t variation = max - min;
415         if (fMinTempVariation > variation) { fMinTempVariation = variation; }
416         if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
417         
418         n++;
419       }
420     }
421   } // loop over fTempArray
422   
423   if (n>0) { // some valid data was found
424     return n;
425   }
426   else { // no good data
427     return (Int_t) kErrorCode;
428   }
429
430 }
431
432 //________________________________________________________________
433 void AliEMCALCalibTimeDep::GetCalibSignalInfo() 
434 {
435   // pick up Preprocessor output, based on fRun (most recent version)
436   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
437   if (entry) {
438     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
439   }
440
441   if (fCalibSignal) { 
442     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
443                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
444                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
445                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
446                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
447                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );               
448   }
449   else {
450     AliWarning( Form("AliCaloCalibSignal not found!") );
451   }
452   
453   return;
454 }
455
456 //________________________________________________________________
457 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo() 
458 {
459   // pick up Preprocessor output, based on fRun (most recent version)
460   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
461   // stored object should be a TTree; read the info
462   if (entry) {
463     fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
464   }
465
466   if (fCalibTempCoeff) { 
467     AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
468   }
469   else {
470     AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
471   }
472   
473   return;
474 }
475
476 //________________________________________________________________
477 void AliEMCALCalibTimeDep::GetCalibReferenceInfo() 
478 {
479   // pick up Preprocessor output, based on fRun (most recent version)
480   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
481   if (entry) {
482     fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
483   }
484
485   if (fCalibReference) { 
486     AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
487   }
488   else {
489     AliWarning( Form("AliEMCALCalibReference 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 fCalibReference
500
501   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
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   for (int i = 0; i < nSM; i++) {
507     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
508     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
509
510     int iSMRef = dataCalibReference->GetSuperModuleNum();
511     int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
512     if (iSMRef != iSMCorr) {
513       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, 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/kSecToHour
529   // corrected for startTime difference: Hour/kSecToHour + timeDiff
530   // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * 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 *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
571   TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
572
573   int iChannelNum = 0; // for regular towers
574   int iRefNum = 0; // for LED
575   double dHour = 0;
576   double dAvgAmp = 0;
577
578   treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
579   treeAvg->SetBranchAddress("fHour", &dHour);
580   treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
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<treeAvg->GetEntries(); ient++) {
588     treeAvg->GetEntry(ient);
589     // figure out where this info comes from
590     fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
591     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
592     // add value in the arrays
593     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
594     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
595     nCount[iSM]++;
596   }
597
598   treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
599   treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
600   treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
601
602   for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
603     treeLEDRefAvg->GetEntry(ient);
604     // figure out where this info comes from
605     fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
606     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
607     // add value in the arrays
608     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, 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 kErrorCode values for those that could not be set
649
650   Float_t ratiot0 = 0;
651   Float_t ratioT = 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     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
657     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
658     iSM = dataCalibReference->GetSuperModuleNum();
659
660     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
661       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
662       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
663         refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration      
664
665       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
666
667         // Calc. R(t0):
668         AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
669         iGain = refVal->GetHighLow(); // gain value used for reference calibration      
670         // valid amplitude values should be larger than 0
671         if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
672           ratiot0 =  refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
673         }
674         else {
675           ratiot0 = kErrorCode;
676         }
677
678         // Calc. R(T)
679         for (int j = 0; j < nBins; j++) {
680
681           // calculate R(T) also; first try with individual tower:
682           // same gain as for reference calibration is the default
683           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
684             // looks like valid data with the right gain combination
685             ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); 
686
687             // if data appears to be saturated, and we are in high gain, then try with low gain instead
688             int newGain = iGain;
689             int newRefGain = refGain;
690             if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
691               newGain = 0;
692             }
693             if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) { 
694               newRefGain = 0;
695             }
696
697             if (newGain!=iGain || newRefGain!=refGain) {
698               // compensate for using different gain than in the reference calibration
699               // we may need to have a custom H/L ratio value for each tower
700               // later, but for now just use a common value, as the rest of the code does.. 
701               ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j); 
702
703               if (newGain<iGain) {
704                 ratioT *= fHighLowGainFactor;
705               }
706               else if (newRefGain<refGain) {
707                 ratioT /= fHighLowGainFactor;
708               }
709             }
710           }
711           else {
712             ratioT = kErrorCode;
713           }
714
715           // Calc. correction factor
716           if (ratiot0>0 && ratioT>0) { 
717             correction = ratiot0/ratioT;
718           }
719           else {
720             correction = kErrorCode;
721             nRemaining++;
722           }
723
724           // Store the value
725           dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
726           /* Check that
727           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
728           also works OK */
729         } // nBins
730       }
731     }
732   }
733
734   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
735   return nRemaining;
736 }
737
738 //________________________________________________________________
739 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) 
740 { // use averages for each strip if no good values exist for some single tower 
741
742   // go over fTimeDepCorrection info
743   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
744
745   // for calculating StripAverage info
746   int nValidTower = 0;
747   Float_t stripAverage = 0;
748   Float_t val = 0;
749
750   int iSM = 0;
751   int iCol = 0;
752   int iRow = 0;
753   int iStrip = 0;
754   int firstCol = 0;
755   int lastCol = 0;
756
757   for (int i = 0; i < nSM; i++) {
758     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
759     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
760
761     for (int j = 0; j < nBins; j++) {
762
763       nValidTower = 0;
764       stripAverage = 0;
765
766       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
767         firstCol = iStrip*2;
768         if ((iSM%2)==1) { // C side
769           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
770         }
771         lastCol = firstCol+1;
772
773         for (iCol = firstCol; iCol <= lastCol; iCol++) {
774           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
775             val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
776             if (val>0) { // valid value; error code is negative
777               stripAverage += val;
778               nValidTower++;
779             } 
780           }
781         }
782
783         // calc average over strip
784         if (nValidTower>0) {
785           stripAverage /= nValidTower;
786           for (iCol = firstCol; iCol <= lastCol; iCol++) {
787             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
788               val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
789               if (val<0) { // invalid value; error code is negative
790                 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
791               }
792             }
793           }
794         }
795         else { // could not fill in unvalid entries
796           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
797         }
798
799       } // iStrip
800     } // j, bins
801   } // iSM
802
803   return nRemaining;
804 }
805
806 //________________________________________________________________
807 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize) 
808 { // OK, so we didn't have valid LED data that allowed us to do the correction only 
809   // with that info.
810   // So, instead we'll rely on the temperature info and try to do the correction
811   // based on that instead.
812   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class 
813   Int_t nRemaining = 0;
814
815   int iSM = 0;
816   int iCol = 0;
817   int iRow = 0;
818
819   Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
820   memset(dTempCoeff, 0, sizeof(dTempCoeff));
821   Double_t correction = 0;
822   Double_t secondsPerBin = (Double_t) binSize;
823
824   for (int i = 0; i < nSM; i++) {
825     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
826     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
827
828     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
829     AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
830
831     // first get CalibTempCoeff info
832     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
833       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
834
835         dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
836         if (fVerbosity > 1) {
837           cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow 
838                << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl; 
839         }
840       }
841     }
842     
843     // figure out what the reference temperature is, from fCalibReference
844     Double_t referenceTemperature = 0;
845     int nVal = 0;
846     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
847       if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
848         referenceTemperature += dataCalibReference->GetTemperature(iSensor);
849         nVal++;
850       }
851     }
852     
853     if (nVal>0) {
854       referenceTemperature /= nVal; // valid values exist, we can look into corrections
855       
856       for (int j = 0; j < nBins; j++) {
857         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
858         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
859       // 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)
860         Double_t dSMTemperature = GetTemperatureSM(iSM, timeStamp); 
861
862         Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
863         if (fVerbosity > 0) {
864           cout << " referenceTemperature " << referenceTemperature
865                << " dSMTemperature " << dSMTemperature
866                << " temperatureDiff " << temperatureDiff
867                << endl;
868         }
869         // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down 
870         // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
871         // dTempCoeff is a (unsigned) factor describing how many % the gain 
872         // changes with a degree change.  
873         // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
874         // The correction we want to keep is what we should multiply our ADC value with as a function
875         // of time, i.e. the inverse of the gain change..
876         if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
877              && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
878           // significant enough difference that we need to consider it, and also not unreasonably large
879
880           // loop over all towers; effect of temperature change will depend on gain for this tower
881           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
882             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
883
884               // the correction should be inverse of modification in gain: (see discussion above)
885               // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01; 
886               // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
887               correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]); 
888               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
889               if (fVerbosity > 1) {
890                 cout << " iSM " << iSM
891                      << " iCol  " << iCol 
892                      << " iRow  " << iRow 
893                      << " j  " << j 
894                      << " correction  " << correction 
895                      << endl;
896               }
897             }
898           }
899
900         } // if noteworthy temperature change
901         else { // just set correction values to 1.0
902           correction = 1;
903           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
904             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
905               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
906             }
907           }
908         } // else
909       } // j, Bins
910
911     } // if reference temperature exist 
912     else { // could not do the needed check.. signal that in the return code
913       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
914     }
915   } // iSM
916
917   return nRemaining;
918 }
919