Line data Source code
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 42 : ClassImp(AliEMCALCalibTimeDep)
63 :
64 : //________________________________________________________________
65 0 : AliEMCALCalibTimeDep::AliEMCALCalibTimeDep() :
66 0 : fRun(0),
67 0 : fStartTime(0),
68 0 : fEndTime(0),
69 0 : fMinTemp(0),
70 0 : fMaxTemp(0),
71 0 : fMinTempVariation(0),
72 0 : fMaxTempVariation(0),
73 0 : fMinTempValid(15),
74 0 : fMaxTempValid(35),
75 0 : fMinTime(0),
76 0 : fMaxTime(0),
77 0 : fTemperatureResolution(0.1), // 0.1 deg C is default
78 0 : fMaxTemperatureDiff(5), // 5 deg C is default max diff relative to reference
79 0 : fTimeBinsPerHour(2), // 2 30-min bins per hour is default
80 0 : fHighLowGainFactor(16), // factor ~16 between High gain and low gain
81 0 : fTempArray(NULL),
82 0 : fCalibSignal(NULL),
83 0 : fCalibTempCoeff(NULL),
84 0 : fCalibReference(NULL),
85 0 : fCalibTimeDepCorrection(NULL),
86 0 : fVerbosity(0)
87 0 : {
88 : // Constructor
89 0 : }
90 :
91 : //________________________________________________________________
92 : AliEMCALCalibTimeDep::AliEMCALCalibTimeDep(const AliEMCALCalibTimeDep& calibt) :
93 0 : TObject(calibt),
94 0 : fRun(calibt.GetRunNumber()),
95 0 : fStartTime(calibt.GetStartTime()),
96 0 : fEndTime(calibt.GetEndTime()),
97 0 : fMinTemp(calibt.GetMinTemp()),
98 0 : fMaxTemp(calibt.GetMaxTemp()),
99 0 : fMinTempVariation(calibt.GetMinTempVariation()),
100 0 : fMaxTempVariation(calibt.GetMaxTempVariation()),
101 0 : fMinTempValid(calibt.GetMinTempValid()),
102 0 : fMaxTempValid(calibt.GetMaxTempValid()),
103 0 : fMinTime(calibt.GetMinTime()),
104 0 : fMaxTime(calibt.GetMaxTime()),
105 0 : fTemperatureResolution(calibt.GetTemperatureResolution()),
106 0 : fMaxTemperatureDiff(calibt.GetMaxTemperatureDiff()),
107 0 : fTimeBinsPerHour(calibt.GetTimeBinsPerHour()),
108 0 : fHighLowGainFactor(calibt.GetHighLowGainFactor()),
109 0 : fTempArray(calibt.GetTempArray()),
110 0 : fCalibSignal(calibt.GetCalibSignal()),
111 0 : fCalibTempCoeff(calibt.GetCalibTempCoeff()),
112 0 : fCalibReference(calibt.GetCalibReference()),
113 0 : fCalibTimeDepCorrection(calibt.GetCalibTimeDepCorrection()),
114 0 : fVerbosity(calibt.GetVerbosity())
115 0 : {
116 : // copy constructor
117 0 : }
118 :
119 :
120 : //________________________________________________________________
121 : AliEMCALCalibTimeDep &AliEMCALCalibTimeDep::operator =(const AliEMCALCalibTimeDep& calibt)
122 : {
123 : // assignment operator; use copy ctor
124 0 : if (&calibt == this) return *this;
125 :
126 0 : new (this) AliEMCALCalibTimeDep(calibt);
127 0 : return *this;
128 0 : }
129 :
130 : //________________________________________________________________
131 : AliEMCALCalibTimeDep::~AliEMCALCalibTimeDep()
132 0 : {
133 : // Destructor
134 0 : }
135 :
136 : //________________________________________________________________
137 : void AliEMCALCalibTimeDep::Reset()
138 : {
139 : // clear variables to default
140 0 : fRun = 0;
141 0 : fStartTime = 0;
142 0 : fEndTime = 0;
143 0 : fMinTemp = 0;
144 0 : fMaxTemp = 0;
145 0 : fMinTempVariation = 0;
146 0 : fMaxTempVariation = 0;
147 0 : fMinTempValid = 15;
148 0 : fMaxTempValid = 35;
149 0 : fMinTime = 0;
150 0 : fMaxTime = 0;
151 0 : fTemperatureResolution = 0.1; // 0.1 deg C is default
152 0 : fMaxTemperatureDiff = 5; // 5 deg C is default max diff relative to reference
153 0 : fTimeBinsPerHour = 2; // 2 30-min bins per hour is default
154 0 : fTempArray = NULL;
155 0 : fCalibSignal = NULL;
156 0 : fCalibTempCoeff = NULL;
157 0 : fCalibReference = NULL;
158 0 : fCalibTimeDepCorrection = NULL;
159 0 : fVerbosity = 0;
160 0 : return;
161 : }
162 :
163 : //________________________________________________________________
164 : void AliEMCALCalibTimeDep::PrintInfo() const
165 : {
166 : // print some info
167 0 : cout << endl << " AliEMCALCalibTimeDep::PrintInfo() " << endl;
168 : // basic variables, all 'publicly available' also
169 0 : cout << " VARIABLE DUMP: " << endl
170 0 : << " GetStartTime() " << GetStartTime() << endl
171 0 : << " GetEndTime() " << GetEndTime() << endl
172 0 : << " GetMinTime() " << GetMinTime() << endl
173 0 : << " GetMaxTime() " << GetMaxTime() << endl
174 0 : << " GetMinTemp() " << GetMinTemp() << endl
175 0 : << " GetMaxTemp() " << GetMaxTemp() << endl
176 0 : << " GetMinTempVariation() " << GetMinTempVariation() << endl
177 0 : << " GetMaxTempVariation() " << GetMaxTempVariation() << endl
178 0 : << " GetTemperatureResolution() " << GetTemperatureResolution() << endl;
179 : // run ranges
180 0 : cout << " RUN INFO: " << endl
181 0 : << " runnumber " << GetRunNumber() << endl
182 0 : << " length (in hours) " << GetLengthOfRunInHours() << endl
183 0 : << " length (in bins) " << GetLengthOfRunInBins() << endl
184 0 : << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
185 0 : << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
186 0 : << endl;
187 :
188 0 : return;
189 : }
190 :
191 : //________________________________________________________________
192 : Double_t AliEMCALCalibTimeDep::GetLengthOfRunInHours() const
193 : {
194 0 : return (fEndTime - fStartTime)*kSecToHour;
195 : }
196 :
197 : //________________________________________________________________
198 : Double_t AliEMCALCalibTimeDep::GetLengthOfRunInBins() const
199 : {
200 0 : return (fEndTime - fStartTime)*kSecToHour*fTimeBinsPerHour;
201 : }
202 :
203 : //________________________________________________________________
204 : Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInHours() const
205 : {
206 0 : return (fMaxTime - fMinTime)*kSecToHour;
207 : }
208 :
209 : //________________________________________________________________
210 : Double_t AliEMCALCalibTimeDep::GetRangeOfTempMeasureInDegrees() const
211 : {
212 0 : 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 0 : Reset(); // start fresh
220 :
221 0 : fRun = run;
222 0 : fStartTime = startTime;
223 0 : fEndTime = endTime;
224 :
225 : // collect the needed information
226 0 : GetTemperatureInfo(); // temperature readings during the run
227 0 : ScanTemperatureInfo(); // see what the boundaries are (Min/Max Time/Temp)
228 :
229 0 : 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 0 : Double_t timeHour = (timeStamp - fStartTime) * kSecToHour;
238 :
239 : int n = 0;
240 0 : Double_t valArr[8]={0}; // 8 sensors per SM
241 :
242 0 : for (int i=0; i<fTempArray->NumSensors(); i++) {
243 :
244 0 : AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
245 0 : int module = st->GetSector()*2 + st->GetSide();
246 0 : if ( module == imod ) { // right module
247 : // check if we had valid data for the time that is being asked for
248 0 : if ( timeStamp>=st->GetStartTime() && timeStamp<=st->GetEndTime() ) {
249 0 : AliSplineFit *f = st->GetFit();
250 0 : 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 0 : Double_t val = f->Eval(timeHour);
253 0 : if ( fVerbosity > 0 ) {
254 0 : cout << " sensor i " << i << " val " << val << endl;
255 0 : }
256 0 : if (val>fMinTempValid && val<fMaxTempValid && n<8) {
257 0 : valArr[n] = val;
258 0 : n++;
259 0 : }
260 0 : }
261 0 : } // time
262 : }
263 :
264 : } // loop over fTempArray
265 :
266 0 : if (n>0) { // some valid data was found
267 0 : Double_t median = TMath::Median(n, valArr);
268 : Double_t average = 0;
269 : Int_t nval = 0;
270 0 : for (int is=0; is<n; is++) {
271 0 : if (TMath::Abs(valArr[is] - median) < kTempMaxDiffMedian) {
272 0 : average += valArr[is];
273 0 : nval++;
274 0 : }
275 : }
276 : //cout << " n " << n << " nval " << nval << " median " << median << endl;
277 0 : if (nval > 0) {
278 0 : average /= nval;
279 : //cout << " average " << average << endl;
280 0 : return average;
281 : }
282 : else { // this case should not happen, but kept for completeness (coverity etc)
283 0 : return median;
284 : }
285 : }
286 : else { // no good data
287 0 : return kErrorCode;
288 : }
289 :
290 0 : }
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 0 : Int_t nSM = fCalibReference->GetNSuperModule();
305 : // how many time-bins should we have for this run?
306 0 : Int_t nBins = (Int_t) (GetLengthOfRunInBins() + 1); // round-up (from double to int; always at least 1)
307 0 : 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 0 : if (fMaxTempVariation < fTemperatureResolution) {
312 : nBins = 1; // just one bin needed..
313 0 : }
314 0 : if (nBins == 1) {
315 0 : binSize = fEndTime - fStartTime;
316 0 : }
317 0 : if (fVerbosity > 0) {
318 0 : cout << " nBins " << nBins << " binSize " << binSize << endl;
319 0 : }
320 :
321 : // set up a reasonable default (correction = 1.0)
322 0 : fCalibTimeDepCorrection = new AliEMCALCalibTimeDepCorrection(nSM);
323 0 : fCalibTimeDepCorrection->InitCorrection(nSM, nBins, 1.0);
324 0 : fCalibTimeDepCorrection->SetStartTime(fStartTime);
325 0 : fCalibTimeDepCorrection->SetNTimeBins(nBins);
326 0 : fCalibTimeDepCorrection->SetTimeBinSize(binSize);
327 :
328 : // 2: try with Temperature correction
329 0 : Int_t nRemaining = CalcTemperatureCorrection(nSM, nBins, binSize);
330 :
331 0 : return nRemaining;
332 0 : }
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 0 : Double_t dP0 = kTempCoeffP0Const + kTempCoeffP0Factor * IDark;
341 0 : Double_t dP1 = kTempCoeffP1Const + kTempCoeffP1Factor * IDark;
342 :
343 0 : Double_t dTC = dP0 + dP1*M;
344 : // from % numbers to regular ones..:
345 0 : dTC *= 0.01;
346 :
347 0 : 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 0 : AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Temperature", fRun);
356 0 : if (entry) {
357 0 : fTempArray = (AliEMCALSensorTempArray *) entry->GetObject();
358 0 : }
359 :
360 0 : if (fTempArray) {
361 0 : AliInfo( Form("NumSensors %d - IdDCS: first %d last %d",
362 : fTempArray->NumSensors(),
363 : fTempArray->GetFirstIdDCS(), fTempArray->GetLastIdDCS() ) );
364 0 : }
365 : else {
366 0 : AliWarning( Form("AliEMCALSensorTempArray not found!") );
367 : }
368 :
369 : return;
370 0 : }
371 :
372 : //________________________________________________________________
373 : Int_t AliEMCALCalibTimeDep::ScanTemperatureInfo()
374 : {// assign max/min time and temperature values
375 :
376 0 : fMinTemp = 999; // init to some large value (999 deg C)
377 0 : fMaxTemp = 0;
378 0 : fMinTempVariation = 999; // init to some large value (999 deg C)
379 0 : fMaxTempVariation = 0;
380 0 : fMinTime = 2147483647; // init to a large value in the far future (0x7fffffff), year 2038 times..
381 0 : fMaxTime = 0;
382 :
383 : Int_t n = 0; // number of valid readings
384 :
385 0 : for (int i=0; i<fTempArray->NumSensors(); i++) {
386 :
387 0 : AliEMCALSensorTemp *st = fTempArray->GetSensor(i);
388 0 : if ( st->GetStartTime() == 0 ) { // no valid data
389 0 : continue;
390 : }
391 :
392 : // check time ranges
393 0 : if (fMinTime > st->GetStartTime()) { fMinTime = st->GetStartTime(); }
394 0 : if (fMaxTime < st->GetEndTime()) { fMaxTime = st->GetEndTime(); }
395 :
396 : // check temperature ranges
397 0 : AliSplineFit *f = st->GetFit();
398 :
399 0 : if (f) { // ok, looks like we have valid data/info
400 0 : int np = f->GetKnots();
401 0 : 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 0 : for (int ip=0; ip<np; ip++) {
407 0 : if (y0[ip]>fMinTempValid && y0[ip]<fMaxTempValid) {
408 0 : if (min > y0[ip]) { min = y0[ip]; }
409 0 : if (max < y0[ip]) { max = y0[ip]; }
410 0 : nval++;
411 0 : }
412 : }
413 0 : if (nval>0) {
414 0 : if (fMinTemp > min) { fMinTemp = min; }
415 0 : if (fMaxTemp < max) { fMaxTemp = max; }
416 0 : Double_t variation = max - min;
417 0 : if (fMinTempVariation > variation) { fMinTempVariation = variation; }
418 0 : if (fMaxTempVariation < variation) { fMaxTempVariation = variation; }
419 :
420 0 : n++;
421 0 : }
422 0 : }
423 0 : } // loop over fTempArray
424 :
425 0 : if (n>0) { // some valid data was found
426 0 : return n;
427 : }
428 : else { // no good data
429 0 : return (Int_t) kErrorCode;
430 : }
431 :
432 0 : }
433 :
434 : //________________________________________________________________
435 : void AliEMCALCalibTimeDep::GetCalibSignalInfo()
436 : {
437 : // pick up Preprocessor output, based on fRun (most recent version)
438 0 : AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/LED", fRun);
439 0 : if (entry) {
440 0 : fCalibSignal = (AliCaloCalibSignal *) entry->GetObject();
441 0 : }
442 :
443 0 : if (fCalibSignal) {
444 0 : 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 0 : }
451 : else {
452 0 : AliWarning( Form("AliCaloCalibSignal not found!") );
453 : }
454 :
455 : return;
456 0 : }
457 :
458 : //________________________________________________________________
459 : void AliEMCALCalibTimeDep::GetCalibTempCoeffInfo()
460 : {
461 : // pick up Preprocessor output, based on fRun (most recent version)
462 0 : AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/TempCoeff", fRun);
463 : // stored object should be a TTree; read the info
464 0 : if (entry) {
465 0 : fCalibTempCoeff = (AliEMCALCalibTempCoeff *) entry->GetObject();
466 0 : }
467 :
468 0 : if (fCalibTempCoeff) {
469 0 : AliInfo( Form("CalibTempCoeff: NSuperModule %d ", fCalibTempCoeff->GetNSuperModule() ) );
470 0 : }
471 : else {
472 0 : AliWarning( Form("AliEMCALCalibTempCoeff not found!") );
473 : }
474 :
475 : return;
476 0 : }
477 :
478 : //________________________________________________________________
479 : void AliEMCALCalibTimeDep::GetCalibReferenceInfo()
480 : {
481 : // pick up Preprocessor output, based on fRun (most recent version)
482 0 : AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Reference", fRun);
483 0 : if (entry) {
484 0 : fCalibReference = (AliEMCALCalibReference *) entry->GetObject();
485 0 : }
486 :
487 0 : if (fCalibReference) {
488 0 : AliInfo( Form("CalibReference: NSuperModule %d ", fCalibReference->GetNSuperModule() ) );
489 0 : }
490 : else {
491 0 : AliWarning( Form("AliEMCALCalibReference not found!") );
492 : }
493 :
494 : return;
495 0 : }
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 0 : for (int i = 0; i < nSM; i++) {
509 0 : AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
510 0 : AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
511 :
512 0 : int iSMRef = dataCalibReference->GetSuperModuleNum();
513 0 : int iSMCorr = dataCalibTimeDepCorrection->GetSuperModuleNum();
514 0 : if (iSMRef != iSMCorr) {
515 0 : AliWarning( Form("AliEMCALCalibTimeDep - SuperModule index mismatch: %d != %d", iSMRef, iSMCorr) );
516 0 : nRemaining = nSM * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
517 0 : return nRemaining;
518 : }
519 0 : }
520 :
521 0 : int iSM = 0;
522 0 : int iCol = 0;
523 0 : int iRow = 0;
524 0 : int iStrip = 0;
525 0 : 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 0 : 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 0 : TArrayF ampT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
537 0 : TArrayF nT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows][2];
538 0 : TArrayF ampLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
539 0 : TArrayF nLEDRefT[AliEMCALGeoParams::fgkEMCALModules][AliEMCALGeoParams::fgkEMCALLEDRefs][2];
540 :
541 : // set up TArray's first
542 0 : for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
543 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
544 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
545 0 : for (iGain = 0; iGain < 2; iGain++) {
546 : // length of arrays
547 0 : ampT[iSM][iCol][iRow][iGain].Set(nBins);
548 0 : nT[iSM][iCol][iRow][iGain].Set(nBins);
549 : // content of arrys
550 0 : for (int j = 0; j < nBins; j++) {
551 0 : ampT[iSM][iCol][iRow][iGain].AddAt(0, j);
552 0 : nT[iSM][iCol][iRow][iGain].AddAt(0, j);
553 : }
554 : }
555 : }
556 : }//iCol
557 0 : for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
558 0 : for (iGain = 0; iGain < 2; iGain++) {
559 : // length of arrays
560 0 : ampLEDRefT[iSM][iStrip][iGain].Set(nBins);
561 0 : nLEDRefT[iSM][iStrip][iGain].Set(nBins);
562 : // content of arrys
563 0 : for (int j = 0; j < nBins; j++) {
564 0 : ampLEDRefT[iSM][iStrip][iGain].AddAt(0, j);
565 0 : 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 0 : TTree *treeAvg = fCalibSignal->GetTreeAvgAmpVsTime();
573 0 : TTree *treeLEDRefAvg = fCalibSignal->GetTreeAvgAmpVsTime();
574 :
575 0 : int iChannelNum = 0; // for regular towers
576 0 : int iRefNum = 0; // for LED
577 0 : double dHour = 0;
578 0 : double dAvgAmp = 0;
579 :
580 0 : treeAvg->SetBranchAddress("fChannelNum", &iChannelNum);
581 0 : treeAvg->SetBranchAddress("fHour", &dHour);
582 0 : treeAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
583 :
584 : int iBin = 0;
585 : // counters for how many values were seen per SuperModule
586 0 : int nCount[AliEMCALGeoParams::fgkEMCALModules] = {0};
587 0 : int nCountLEDRef[AliEMCALGeoParams::fgkEMCALModules] = {0};
588 :
589 0 : for (int ient=0; ient<treeAvg->GetEntries(); ient++) {
590 0 : treeAvg->GetEntry(ient);
591 : // figure out where this info comes from
592 0 : fCalibSignal->DecodeChannelNum(iChannelNum, &iSM, &iCol, &iRow, &iGain);
593 0 : iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
594 : // add value in the arrays
595 0 : ampT[iSM][iCol][iRow][iGain].AddAt( ampT[iSM][iCol][iRow][iGain].At(iBin)+dAvgAmp, iBin );
596 0 : nT[iSM][iCol][iRow][iGain].AddAt( nT[iSM][iCol][iRow][iGain].At(iBin)+1, iBin );
597 0 : nCount[iSM]++;
598 : }
599 :
600 0 : treeLEDRefAvg->SetBranchAddress("fRefNum", &iRefNum);
601 0 : treeLEDRefAvg->SetBranchAddress("fHour", &dHour);
602 0 : treeLEDRefAvg->SetBranchAddress("fAvgAmp",&dAvgAmp);
603 :
604 0 : for (int ient=0; ient<treeLEDRefAvg->GetEntries(); ient++) {
605 0 : treeLEDRefAvg->GetEntry(ient);
606 : // figure out where this info comes from
607 0 : fCalibSignal->DecodeRefNum(iRefNum, &iSM, &iStrip, &iGain);
608 0 : iBin = (int) ( (dHour + timeDiff*kSecToHour) * fTimeBinsPerHour ); // CHECK!!!
609 : // add value in the arrays
610 0 : ampLEDRefT[iSM][iStrip][iGain].AddAt( ampLEDRefT[iSM][iStrip][iGain].At(iBin)+dAvgAmp, iBin );
611 0 : nLEDRefT[iSM][iStrip][iGain].AddAt( nLEDRefT[iSM][iStrip][iGain].At(iBin)+1, iBin );
612 0 : nCountLEDRef[iSM]++;
613 : }
614 :
615 : // Normalize TArray values, and calculate average also
616 : Float_t norm = 0; // extra var, for readability
617 :
618 0 : for (iSM = 0; iSM < AliEMCALGeoParams::fgkEMCALModules; iSM++) {
619 0 : if (nCount[iSM]>0 && nCountLEDRef[iSM]>0) { // avoid SuperModules with no data..
620 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
621 : // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
622 0 : iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
623 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
624 0 : for (iGain = 0; iGain < 2; iGain++) {
625 : // content of arrys
626 0 : for (int j = 0; j < nBins; j++) {
627 0 : if (nT[iSM][iCol][iRow][iGain].At(j) > 0) {
628 0 : norm = ampT[iSM][iCol][iRow][iGain].At(j) / nT[iSM][iCol][iRow][iGain].At(j);
629 0 : ampT[iSM][iCol][iRow][iGain].AddAt(norm, j); // AddAt = SetAt
630 : }
631 : }
632 : }
633 : }
634 : }//iCol
635 0 : for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
636 0 : for (iGain = 0; iGain < 2; iGain++) {
637 0 : for (int j = 0; j < nBins; j++) {
638 0 : if (nLEDRefT[iSM][iStrip][iGain].At(j) > 0) {
639 0 : norm = ampLEDRefT[iSM][iStrip][iGain].At(j) / nLEDRefT[iSM][iStrip][iGain].At(j);
640 0 : 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 0 : for (int i = 0; i < nSM; i++) {
658 0 : AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(i);
659 0 : AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
660 0 : iSM = dataCalibReference->GetSuperModuleNum();
661 :
662 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
663 : // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
664 0 : iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2; //TMP, FIXME
665 0 : refGain = dataCalibReference->GetLEDRefHighLow(iStrip); // LED reference gain value used for reference calibration
666 :
667 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
668 :
669 : // Calc. R(t0):
670 0 : AliEMCALCalibReferenceVal * refVal = dataCalibReference->GetAPDVal(iCol, iRow);
671 0 : iGain = refVal->GetHighLow(); // gain value used for reference calibration
672 : // valid amplitude values should be larger than 0
673 0 : if (refVal->GetLEDAmp()>0 && dataCalibReference->GetLEDRefAmp(iStrip)>0) {
674 0 : ratiot0 = refVal->GetLEDAmp() / dataCalibReference->GetLEDRefAmp(iStrip);
675 0 : }
676 : else {
677 : ratiot0 = kErrorCode;
678 : }
679 :
680 : // Calc. R(T)
681 0 : 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 0 : 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 0 : 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 0 : int newGain = iGain;
691 : int newRefGain = refGain;
692 0 : if ( ampT[iSM][iCol][iRow][iGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && iGain==1 ) {
693 : newGain = 0;
694 0 : }
695 0 : if ( ampLEDRefT[iSM][iStrip][refGain].At(j)>AliEMCALGeoParams::fgkOverflowCut && refGain==1 ) {
696 : newRefGain = 0;
697 0 : }
698 :
699 0 : 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 0 : ratioT = ampT[iSM][iCol][iRow][newGain].At(j) / ampLEDRefT[iSM][iStrip][newRefGain].At(j);
704 :
705 0 : if (newGain<iGain) {
706 0 : ratioT *= fHighLowGainFactor;
707 0 : }
708 0 : else if (newRefGain<refGain) {
709 0 : ratioT /= fHighLowGainFactor;
710 0 : }
711 : }
712 0 : }
713 : else {
714 : ratioT = kErrorCode;
715 : }
716 :
717 : // Calc. correction factor
718 0 : if (ratiot0>0 && ratioT>0) {
719 0 : correction = ratiot0/ratioT;
720 0 : }
721 : else {
722 : correction = kErrorCode;
723 0 : nRemaining++;
724 : }
725 :
726 : // Store the value
727 0 : 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 0 : nRemaining = CalcLEDCorrectionStripBasis(nSM, nBins);
737 : return nRemaining;
738 0 : }
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 0 : for (int i = 0; i < nSM; i++) {
760 0 : AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
761 0 : iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
762 :
763 0 : for (int j = 0; j < nBins; j++) {
764 :
765 : nValidTower = 0;
766 : stripAverage = 0;
767 :
768 0 : for (iStrip = 0; iStrip < AliEMCALGeoParams::fgkEMCALLEDRefs; iStrip++) {
769 0 : firstCol = iStrip*2;
770 0 : if ((iSM%2)==1) { // C side
771 0 : firstCol = (AliEMCALGeoParams::fgkEMCALLEDRefs-1 - iStrip)*2;
772 0 : }
773 0 : lastCol = firstCol+1;
774 :
775 0 : for (iCol = firstCol; iCol <= lastCol; iCol++) {
776 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
777 0 : val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
778 0 : if (val>0) { // valid value; error code is negative
779 0 : stripAverage += val;
780 0 : nValidTower++;
781 0 : }
782 : }
783 : }
784 :
785 : // calc average over strip
786 0 : if (nValidTower>0) {
787 0 : stripAverage /= nValidTower;
788 0 : for (iCol = firstCol; iCol <= lastCol; iCol++) {
789 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
790 0 : val = dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->At(j);
791 0 : if (val<0) { // invalid value; error code is negative
792 0 : dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(val, j);
793 0 : }
794 : }
795 : }
796 : }
797 : else { // could not fill in unvalid entries
798 0 : nRemaining += 2*AliEMCALGeoParams::fgkEMCALRows;
799 : }
800 :
801 : } // iStrip
802 : } // j, bins
803 : } // iSM
804 :
805 0 : 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 0 : Double_t dTempCoeff[AliEMCALGeoParams::fgkEMCALCols][AliEMCALGeoParams::fgkEMCALRows];
822 0 : memset(dTempCoeff, 0, sizeof(dTempCoeff));
823 : Double_t correction = 0;
824 0 : Double_t secondsPerBin = (Double_t) binSize;
825 :
826 0 : for (int i = 0; i < nSM; i++) {
827 0 : AliEMCALSuperModuleCalibTimeDepCorrection * dataCalibTimeDepCorrection = fCalibTimeDepCorrection->GetSuperModuleCalibTimeDepCorrectionNum(i);
828 0 : iSM = dataCalibTimeDepCorrection->GetSuperModuleNum();
829 :
830 0 : AliEMCALSuperModuleCalibReference * dataCalibReference = fCalibReference->GetSuperModuleCalibReferenceNum(iSM);
831 0 : AliEMCALSuperModuleCalibTempCoeff * dataCalibTempCoeff = fCalibTempCoeff->GetSuperModuleCalibTempCoeffNum(iSM);
832 :
833 : // first get CalibTempCoeff info
834 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
835 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
836 :
837 0 : dTempCoeff[iCol][iRow] = dataCalibTempCoeff->GetTC(iCol, iRow);
838 0 : if (fVerbosity > 1) {
839 0 : cout << " iSM " << iSM << " iCol " << iCol << " iRow " << iRow
840 0 : << " dTempCoeff " << dTempCoeff[iCol][iRow] << endl;
841 0 : }
842 : }
843 : }
844 :
845 : // figure out what the reference temperature is, from fCalibReference
846 : Double_t referenceTemperature = 0;
847 : int nVal = 0;
848 0 : for (int iSensor = 0; iSensor<AliEMCALGeoParams::fgkEMCALTempSensors ; iSensor++) {
849 0 : if (dataCalibReference->GetTemperature(iSensor)>0) { // hopefully OK value
850 0 : referenceTemperature += dataCalibReference->GetTemperature(iSensor);
851 0 : nVal++;
852 0 : }
853 : }
854 :
855 0 : if (nVal>0) {
856 0 : referenceTemperature /= nVal; // valid values exist, we can look into corrections
857 :
858 : Double_t dSMTemperature = 0;
859 0 : 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 0 : 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 0 : dSMTemperature = GetTemperatureSM(iSM, timeStamp);
865 0 : 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 0 : }
869 :
870 0 : Double_t temperatureDiff = referenceTemperature - dSMTemperature; // ref - new
871 0 : if (fVerbosity > 0) {
872 0 : cout << " referenceTemperature " << referenceTemperature
873 0 : << " dSMTemperature " << dSMTemperature
874 0 : << " temperatureDiff " << temperatureDiff
875 0 : << endl;
876 0 : }
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 0 : if ( (TMath::Abs(temperatureDiff)>fTemperatureResolution)
885 0 : && (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 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
890 0 : 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 0 : correction = 1.0 - (temperatureDiff * dTempCoeff[iCol][iRow]);
896 0 : dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
897 0 : if (fVerbosity > 1) {
898 0 : cout << " iSM " << iSM
899 0 : << " iCol " << iCol
900 0 : << " iRow " << iRow
901 0 : << " j " << j
902 0 : << " correction " << correction
903 0 : << endl;
904 0 : }
905 : }
906 : }
907 :
908 : } // if noteworthy temperature change
909 : else { // just set correction values to 1.0
910 : correction = 1;
911 0 : for (iCol = 0; iCol < AliEMCALGeoParams::fgkEMCALCols; iCol++) {
912 0 : for (iRow = 0; iRow < AliEMCALGeoParams::fgkEMCALRows; iRow++) {
913 0 : dataCalibTimeDepCorrection->GetCorrection(iCol,iRow)->AddAt(correction, j);
914 : }
915 : }
916 : } // else
917 : } // j, Bins
918 :
919 0 : } // if reference temperature exist
920 : else { // could not do the needed check.. signal that in the return code
921 0 : nRemaining += AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRows * nBins;
922 : }
923 : } // iSM
924 :
925 0 : return nRemaining;
926 0 : }
927 :
|
