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