put back previous default value, although decission needs to be made to which one...
[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     //cout << " n " << n << " nval " << nval << " median " << median << endl;
277     if (nval >  0) {
278       average /= nval;
279       //cout << " average " << average << endl; 
280       return average;
281     }
282     else { // this case should not happen, but kept for completeness (coverity etc)
283       return median;
284     }
285   }
286   else { // no good data
287     return kErrorCode;
288   }
289
290 }
291
292 //________________________________________________________________
293 Int_t AliEMCALCalibTimeDep::CalcCorrection()
294 { // OK, this is where the real action takes place - the heart of this class..
295   /* The philosophy is as follows:
296      0. Init corrections to 1.0 values, and see how many correction bins we need
297      1. Check how large temperature variations we have through the run - do we really need all the correction bias (otherwise adjust to single bin)
298      2. try to use temperature info + APD temperature coefficient info, to estimate correction.
299      For now (from Dec 2009), we do not use LED info.
300    */
301
302   // 0: Init
303   // how many SuperModules do we have?
304   Int_t nSM = fCalibReference->GetNSuperModule();
305   // how many time-bins should we have for this run?
306   Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
307   Int_t binSize = (Int_t) (3600/fTimeBinsPerHour); // in seconds
308
309   // 1: get info on how much individual sensors might have changed during
310   // the run (compare max-min for each sensor separately)
311   if (fMaxTempVariation < fTemperatureResolution) {
312     nBins = 1; // just one bin needed..
313   }
314   if (nBins == 1) {    
315     binSize = fEndTime - fStartTime;
316   }
317   if (fVerbosity > 0) {
318     cout << " nBins " << nBins << " binSize " << binSize << endl;
319   }
320
321   // set up a reasonable default (correction = 1.0)
322   fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
323   fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
324   fCalibTimeDepCorrection->SetStartTime(fStartTime);
325   fCalibTimeDepCorrection->SetNTimeBins(nBins);
326   fCalibTimeDepCorrection->SetTimeBinSize(binSize);
327
328   // 2: try with Temperature correction 
329   Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
330
331   return nRemaining;
332 }
333
334
335 //________________________________________________________________
336 Double_t AliEMCALCalibTimeDep::GetTempCoeff(Double_t IDark, Double_t M) const
337 { // estimate the Temperature Coefficient, based on the dark current (IDark)
338   // and the gain (M), based on Catania parameterizations
339
340   Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
341   Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
342  
343   Double_t dTC = dP0 + dP1*M;
344   // from % numbers to regular ones..:
345   dTC *= 0.01;
346
347   return TMath::Abs(dTC); // return the absolute value, to avoid any sign confusion
348 }
349
350 /* Next come the methods that do the work in picking up all the needed info..*/
351 //________________________________________________________________
352 void AliEMCALCalibTimeDep::GetTemperatureInfo() 
353 {
354   // pick up Preprocessor output, based on fRun (most recent version)
355   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
356   if (entry) {
357     fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
358   }
359
360   if (fTempArray) { 
361     AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
362                   fTempArray->NumSensors(),
363                   fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
364   }
365   else {
366     AliWarning( Form("AliEMCALSensorTempArray not found!") );
367   }
368   
369   return;
370 }
371
372 //________________________________________________________________
373 Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo() 
374 {// assign max/min time and temperature values
375
376   fMinTemp = 999; // init to some large value (999 deg C)
377   fMaxTemp = 0;
378   fMinTempVariation = 999; // init to some large value (999 deg C)
379   fMaxTempVariation = 0;
380   fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
381   fMaxTime = 0;
382
383   Int_t n = 0; // number of valid readings
384
385   for (int i=0; i<fTempArray->NumSensors(); i++) {
386     
387     AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
388     if ( st->GetStartTime() == 0 ) { // no valid data
389       continue;
390     }
391
392     // check time ranges
393     if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
394     if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
395  
396     // check temperature ranges
397     AliSplineFit *f = st->GetFit();
398
399     if (f) { // ok, looks like we have valid data/info
400       int np = f->GetKnots();
401       Double_t *y0 = f->GetY0();
402       // min and max values within the single sensor
403       Double_t min = 999;
404       Double_t max = 0;
405       int nval = 0;
406       for (int ip=0; ip<np; ip++) { 
407         if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
408           if (min > y0[ip]) { min = y0[ip]; }
409           if (max < y0[ip]) { max = y0[ip]; }
410           nval++;
411         }
412       }
413       if (nval>0) {
414         if (fMinTemp > min) { fMinTemp = min; }
415         if (fMaxTemp < max) { fMaxTemp = max; }
416         Double_t variation = max - min;
417         if (fMinTempVariation > variation) { fMinTempVariation = variation; }
418         if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
419         
420         n++;
421       }
422     }
423   } // loop over fTempArray
424   
425   if (n>0) { // some valid data was found
426     return n;
427   }
428   else { // no good data
429     return (Int_t) kErrorCode;
430   }
431
432 }
433
434 //________________________________________________________________
435 void AliEMCALCalibTimeDep::GetCalibSignalInfo() 
436 {
437   // pick up Preprocessor output, based on fRun (most recent version)
438   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
439   if (entry) {
440     fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
441   }
442
443   if (fCalibSignal) { 
444     AliInfo( Form("CalibSignal: NEvents %d NAcceptedEvents %d Entries %lld AvgEntries LEDRefEntries %lld LEDRefEntries %lld, LEDRefAvgEntries %lld",
445                   fCalibSignal->GetNEvents(), fCalibSignal->GetNAcceptedEvents(),
446                   fCalibSignal->GetTreeAmpVsTime()->GetEntries(),
447                   fCalibSignal->GetTreeAvgAmpVsTime()->GetEntries(),
448                   fCalibSignal->GetTreeLEDAmpVsTime()->GetEntries(),
449                   fCalibSignal->GetTreeLEDAvgAmpVsTime()->GetEntries() ) );               
450   }
451   else {
452     AliWarning( Form("AliCaloCalibSignal not found!") );
453   }
454   
455   return;
456 }
457
458 //________________________________________________________________
459 void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo() 
460 {
461   // pick up Preprocessor output, based on fRun (most recent version)
462   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
463   // stored object should be a TTree; read the info
464   if (entry) {
465     fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
466   }
467
468   if (fCalibTempCoeff) { 
469     AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
470   }
471   else {
472     AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
473   }
474   
475   return;
476 }
477
478 //________________________________________________________________
479 void AliEMCALCalibTimeDep::GetCalibReferenceInfo() 
480 {
481   // pick up Preprocessor output, based on fRun (most recent version)
482   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
483   if (entry) {
484     fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
485   }
486
487   if (fCalibReference) { 
488     AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
489   }
490   else {
491     AliWarning( Form("AliEMCALCalibReference not found!") );
492   }
493   
494   return;
495 }
496
497 //________________________________________________________________
498 Int_t AliEMCALCalibTimeDep::CalcLEDCorrection(Int_t nSM, Int_t nBins) 
499 {// Construct normalized ratios R(t)=LED(t)/LEDRef(t), for current time T and calibration time t0 
500   // The correction factor we keep is c(T) = R(t0)/R(T)
501   // T info from fCalibSignal, t0 info from fCalibReference
502
503   // NOTE: for now we don't use the RMS info either from fCalibSignal or fCalibReference
504   // but one could upgrade this in the future
505   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
506
507   // sanity check; same SuperModule indices for corrections as for regular calibrations
508   for (int i = 0; i < nSM; i++) {
509     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
510     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
511
512     int iSMRef = dataCalibReference->GetSuperModuleNum();
513     int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
514     if (iSMRef != iSMCorr) {
515       AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
516       nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
517       return nRemaining;
518     }
519   }
520  
521   int iSM = 0;
522   int iCol = 0;
523   int iRow = 0;
524   int iStrip = 0;
525   int iGain = 0;
526
527   // The fCalibSignal info is stored in TTrees
528   // Note that the time-bins for the TTree's may not exactly match with our correction time bins 
529   int timeDiff = fCalibSignal->GetStartTime() - fStartTime; // in seconds
530   // fCalibSignal time info in seconds: Hour/kSecToHour
531   // corrected for startTime difference: Hour/kSecToHour + timeDiff
532   // converted into a time-bin we use: (Hour + timeDiff*kSecToHour) * fTimeBinsPerHour
533
534   // values for R(T), size of TArray = nBins
535   // the [2] dimension below is for the low or high gain 
536   TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
537   TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2]; 
538   TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
539   TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2]; 
540
541   // set up TArray's first
542   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
543     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
544       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
545         for (iGain = 0; iGain < 2; iGain++) {
546           // length of arrays
547           ampT[iSM][iCol][iRow][iGain].Set(nBins);
548           nT[iSM][iCol][iRow][iGain].Set(nBins);
549           // content of arrys
550           for (int j = 0; j < nBins; j++) {
551             ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
552             nT[iSM][iCol][iRow][iGain].AddAt(0, j);
553           }
554         }
555       }
556     }//iCol
557     for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
558       for (iGain = 0; iGain < 2; iGain++) {
559         // length of arrays
560         ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
561         nLEDRefT[iSM][iStrip][iGain].Set(nBins);
562         // content of arrys
563         for (int j = 0; j < nBins; j++) {
564           ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
565           nLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
566         }
567       }
568     }//iStrip
569   }
570
571   // OK, now loop over the TTrees and fill the arrays needed for R(T)
572   TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
573   TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
574
575   int iChannelNum = 0; // for regular towers
576   int iRefNum = 0; // for LED
577   double dHour = 0;
578   double dAvgAmp = 0;
579
580   treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
581   treeAvg->SetBranchAddress("fHour", &dHour);
582   treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
583
584   int iBin = 0;
585   // counters for how many values were seen per SuperModule
586   int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
587   int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
588
589   for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
590     treeAvg->GetEntry(ient);
591     // figure out where this info comes from
592     fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
593     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
594     // add value in the arrays
595     ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
596     nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
597     nCount[iSM]++;
598   }
599
600   treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
601   treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
602   treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
603
604   for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
605     treeLEDRefAvg->GetEntry(ient);
606     // figure out where this info comes from
607     fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
608     iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour  ); // CHECK!!!
609     // add value in the arrays
610     ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
611     nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
612     nCountLEDRef[iSM]++;
613   }
614
615   // Normalize TArray values, and calculate average also
616   Float_t norm = 0; // extra var, for readability
617
618   for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
619     if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
620       for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
621         //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
622         iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
623         for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
624           for (iGain = 0; iGain < 2; iGain++) {
625             // content of arrys
626             for (int j = 0; j < nBins; j++) {
627               if (nT[iSM][iCol][iRow][iGain].At(j) > 0) { 
628                 norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j); 
629                 ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
630               } 
631             }
632           }
633         }
634       }//iCol
635       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
636         for (iGain = 0; iGain < 2; iGain++) {
637           for (int j = 0; j < nBins; j++) {
638             if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
639               norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j); 
640               ampLEDRefT[iSM][iStrip][iGain].AddAt(norm, j); // AddAt = SetAt
641             }
642           }
643         }
644       }//iStrip
645     }
646   } // iSM
647
648
649   // Calculate correction values, and store them
650   // set kErrorCode values for those that could not be set
651
652   Float_t ratiot0 = 0;
653   Float_t ratioT = 0;
654   Float_t correction = 0; // c(T) = R(t0)/R(T)
655   Int_t refGain = 0; // typically use low gain for LED reference amplitude (high gain typically well beyond saturation)
656
657   for (int i = 0; i < nSM; i++) {
658     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
659     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
660     iSM = dataCalibReference->GetSuperModuleNum();
661
662     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
663       //      iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
664       iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
665         refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration      
666
667       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
668
669         // Calc. R(t0):
670         AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
671         iGain = refVal->GetHighLow(); // gain value used for reference calibration      
672         // valid amplitude values should be larger than 0
673         if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
674           ratiot0 =  refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
675         }
676         else {
677           ratiot0 = kErrorCode;
678         }
679
680         // Calc. R(T)
681         for (int j = 0; j < nBins; j++) {
682
683           // calculate R(T) also; first try with individual tower:
684           // same gain as for reference calibration is the default
685           if (ampT[iSM][iCol][iRow][iGain].At(j)>0 && ampLEDRefT[iSM][iStrip][refGain].At(j)>0) {
686             // looks like valid data with the right gain combination
687             ratioT = ampT[iSM][iCol][iRow][iGain].At(j) / ampLEDRefT[iSM][iStrip][refGain].At(j); 
688
689             // if data appears to be saturated, and we are in high gain, then try with low gain instead
690             int newGain = iGain;
691             int newRefGain = refGain;
692             if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
693               newGain = 0;
694             }
695             if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) { 
696               newRefGain = 0;
697             }
698
699             if (newGain!=iGain || newRefGain!=refGain) {
700               // compensate for using different gain than in the reference calibration
701               // we may need to have a custom H/L ratio value for each tower
702               // later, but for now just use a common value, as the rest of the code does.. 
703               ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j); 
704
705               if (newGain<iGain) {
706                 ratioT *= fHighLowGainFactor;
707               }
708               else if (newRefGain<refGain) {
709                 ratioT /= fHighLowGainFactor;
710               }
711             }
712           }
713           else {
714             ratioT = kErrorCode;
715           }
716
717           // Calc. correction factor
718           if (ratiot0>0 && ratioT>0) { 
719             correction = ratiot0/ratioT;
720           }
721           else {
722             correction = kErrorCode;
723             nRemaining++;
724           }
725
726           // Store the value
727           dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
728           /* Check that
729           fTimeDepCorrection->SetCorrection(i, iCol, iRow, j, correction);
730           also works OK */
731         } // nBins
732       }
733     }
734   }
735
736   nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
737   return nRemaining;
738 }
739
740 //________________________________________________________________
741 Int_t AliEMCALCalibTimeDep::CalcLEDCorrectionStripBasis(Int_t nSM, Int_t nBins) 
742 { // use averages for each strip if no good values exist for some single tower 
743
744   // go over fTimeDepCorrection info
745   Int_t nRemaining = 0; // we count the towers for which we could not get valid data
746
747   // for calculating StripAverage info
748   int nValidTower = 0;
749   Float_t stripAverage = 0;
750   Float_t val = 0;
751
752   int iSM = 0;
753   int iCol = 0;
754   int iRow = 0;
755   int iStrip = 0;
756   int firstCol = 0;
757   int lastCol = 0;
758
759   for (int i = 0; i < nSM; i++) {
760     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
761     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
762
763     for (int j = 0; j < nBins; j++) {
764
765       nValidTower = 0;
766       stripAverage = 0;
767
768       for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
769         firstCol = iStrip*2;
770         if ((iSM%2)==1) { // C side
771           firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
772         }
773         lastCol = firstCol+1;
774
775         for (iCol = firstCol; iCol <= lastCol; iCol++) {
776           for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
777             val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
778             if (val>0) { // valid value; error code is negative
779               stripAverage += val;
780               nValidTower++;
781             } 
782           }
783         }
784
785         // calc average over strip
786         if (nValidTower>0) {
787           stripAverage /= nValidTower;
788           for (iCol = firstCol; iCol <= lastCol; iCol++) {
789             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
790               val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
791               if (val<0) { // invalid value; error code is negative
792                 dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
793               }
794             }
795           }
796         }
797         else { // could not fill in unvalid entries
798           nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
799         }
800
801       } // iStrip
802     } // j, bins
803   } // iSM
804
805   return nRemaining;
806 }
807
808 //________________________________________________________________
809 Int_t AliEMCALCalibTimeDep::CalcTemperatureCorrection(Int_t nSM, Int_t nBins, Int_t binSize) 
810 { // OK, so we didn't have valid LED data that allowed us to do the correction only 
811   // with that info.
812   // So, instead we'll rely on the temperature info and try to do the correction
813   // based on that instead.
814   // For this, we'll need the info from 3 classes (+temperature array), and output the values in a 4th class 
815   Int_t nRemaining = 0;
816
817   int iSM = 0;
818   int iCol = 0;
819   int iRow = 0;
820
821   Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
822   memset(dTempCoeff, 0, sizeof(dTempCoeff));
823   Double_t correction = 0;
824   Double_t secondsPerBin = (Double_t) binSize;
825
826   for (int i = 0; i < nSM; i++) {
827     AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
828     iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
829
830     AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
831     AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
832
833     // first get CalibTempCoeff info
834     for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
835       for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
836
837         dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
838         if (fVerbosity > 1) {
839           cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow 
840                << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl; 
841         }
842       }
843     }
844     
845     // figure out what the reference temperature is, from fCalibReference
846     Double_t referenceTemperature = 0;
847     int nVal = 0;
848     for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
849       if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
850         referenceTemperature += dataCalibReference->GetTemperature(iSensor);
851         nVal++;
852       }
853     }
854     
855     if (nVal>0) {
856       referenceTemperature /= nVal; // valid values exist, we can look into corrections
857       
858       Double_t dSMTemperature = 0;
859       for (int j = 0; j < nBins; j++) {
860         // what is the timestamp in the middle of this bin? (0.5 is for middle of bin)
861         UInt_t timeStamp = fStartTime + (UInt_t)((j+0.5)*secondsPerBin);
862       // 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)
863         Double_t oldSMTemperature = dSMTemperature;
864         dSMTemperature = GetTemperatureSM(iSM, timeStamp); 
865         if (j>0 && (dSMTemperature==kErrorCode)) { 
866           // if we have previous values, and retrieval of values failed - use that instead (hopefully good)
867           dSMTemperature = oldSMTemperature; 
868         }
869  
870         Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
871         if (fVerbosity > 0) {
872           cout << " referenceTemperature " << referenceTemperature
873                << " dSMTemperature " << dSMTemperature
874                << " temperatureDiff " << temperatureDiff
875                << endl;
876         }
877         // if the new temperature is higher than the old/reference one (diff<0), then the gain has gone down 
878         // if the new temperature is lower than the old/reference one (diff>0), then the gain has gone up
879         // dTempCoeff is a (unsigned) factor describing how many % the gain 
880         // changes with a degree change.  
881         // i.e. the product temperatureDiff * dTempCoeff increase when the gain goes up
882         // The correction we want to keep is what we should multiply our ADC value with as a function
883         // of time, i.e. the inverse of the gain change..
884         if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
885              && (TMath::Abs(temperatureDiff)<fMaxTemperatureDiff) ) {
886           // significant enough difference that we need to consider it, and also not unreasonably large
887
888           // loop over all towers; effect of temperature change will depend on gain for this tower
889           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
890             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
891
892               // the correction should be inverse of modification in gain: (see discussion above)
893               // modification in gain: 1.0 + (temperatureDiff * dTempCoeff[iCol][iRow])*0.01; 
894               // 1/(1+x) ~= 1 - x for small x, i.e. we arrive at:
895               correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]); 
896               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
897               if (fVerbosity > 1) {
898                 cout << " iSM " << iSM
899                      << " iCol  " << iCol 
900                      << " iRow  " << iRow 
901                      << " j  " << j 
902                      << " correction  " << correction 
903                      << endl;
904               }
905             }
906           }
907
908         } // if noteworthy temperature change
909         else { // just set correction values to 1.0
910           correction = 1;
911           for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
912             for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
913               dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
914             }
915           }
916         } // else
917       } // j, Bins
918
919     } // if reference temperature exist 
920     else { // could not do the needed check.. signal that in the return code
921       nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
922     }
923   } // iSM
924
925   return nRemaining;
926 }
927