added low and high flux parameter options
[u/mrichter/AliRoot.git] / EMCAL / beamtest07 / AliEMCALCalibTestBeam.cxx
CommitLineData
82a2d831 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: AliEMCALCalibTestBeam.cxx $ */
17
18//_________________________________________________________________________
19///*-- Author:
20///////////////////////////////////////////////////////////////////////////////
21// //
22// class for EMCAL testbeam calibration //
23// //
24///////////////////////////////////////////////////////////////////////////////
25
26#include <iostream>
27#include <fstream>
28#include <TGraphSmooth.h>
29#include <TFile.h>
30#include <TTree.h>
31#include "AliEMCALCalibTestBeam.h"
32
33/* first a bunch of constants.. */
34const double fkSmoothBandwidth = 5.0; // 5 seems to be a good value..; see $ROOTSYS/tutorials/graphs/motorcycle.C for reference
35const double fkSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
36
37// some global variables for APD handling
38const double fkTempSlope = 0.017; // 1.7% per deg. C, seems about right for all APDs, from studies at Catania, and by Rachid
39const double fkNormTemp = 20; // let's say that 20 degrees is normal, i.e. we'll do the corrections relative to that reference temperature
40
41// some global constants and ROOT objects for temperature handling
42const int fkSensor = 1; // let's just look at temperature data from one sensor
43// let's also have some ranges for valid readings of the temperatures..
44const int fkMinValidTemp = 0;
45const int fkMaxValidTemp = 100;
46
47using namespace std;
48
49ClassImp(AliEMCALCalibTestBeam)
50
51//________________________________________________________________
52AliEMCALCalibTestBeam::AliEMCALCalibTestBeam(const int runNumber)
53{
54 // Constructor
55 // set pointers even though they will be re-assigned in Init()/Reset() later..
56 fTempGraph = NULL;
57 fTempSpline = NULL;
58 Init(runNumber); // collect and setup info; including the spline
59}
60
61//________________________________________________________________
62AliEMCALCalibTestBeam::AliEMCALCalibTestBeam(const AliEMCALCalibTestBeam& calibtb) :
63 fTimeStart(calibtb.GetTimeStart()),
64 fTimeStop(calibtb.GetTimeStop()),
65 fMinTemp(calibtb.GetMinTemp()),
66 fMaxTemp(calibtb.GetMaxTemp()),
67 fMinTime(calibtb.GetMinTime()),
68 fMaxTime(calibtb.GetMaxTime())
69{
70 // copy constructor
71 fTempGraph = calibtb.GetTemperatureGraph();
72 fTempSpline = calibtb.GetTemperatureSpline();
73}
74
75
76//________________________________________________________________
77AliEMCALCalibTestBeam &AliEMCALCalibTestBeam::operator =(const AliEMCALCalibTestBeam& calibtb)
78{
79 // assignment operator; use copy ctor
80 if (&calibtb == this) return *this;
81
82 new (this) AliEMCALCalibTestBeam(calibtb);
83 return *this;
84
85}
86
87//________________________________________________________________
88AliEMCALCalibTestBeam::~AliEMCALCalibTestBeam()
89{
90 // Destructor
91 if (fTempGraph) delete fTempGraph;
92 if (fTempSpline) delete fTempSpline;
93}
94
95//________________________________________________________________
96void AliEMCALCalibTestBeam::Reset()
97{
98 ResetRunInfo();
99 ResetTempInfo();
100 return;
101}
102//________________________________________________________________
103void AliEMCALCalibTestBeam::ResetRunInfo()
104{
105 // clear variables
106 fTimeStart = 0;
107 fTimeStop = 0;
108 fNEvents = 0;
109}
110
111//________________________________________________________________
112void AliEMCALCalibTestBeam::ResetTempInfo()
113{
114 // clear variables
115 fMinTemp = 0;
116 fMaxTemp = 0;
117 fMinTime = 0;
118 fMaxTime = 0;
119
120 if (fTempGraph) delete fTempGraph;
121 fTempGraph = new TGraph();
122
123 if (fTempSpline) delete fTempSpline;
124 fTempSpline = NULL;
125
126 return;
127}
128
129//________________________________________________________________
130void AliEMCALCalibTestBeam::Print() const
131{
132 // print some info
133 cout << endl << " AliEMCALCalibTestBeam::Print() " << endl;
134 // basic variables, all 'publicly available' also
135 cout << " VARIABLE DUMP: " << endl
136 << " GetTimeStart() " << GetTimeStart() << endl
137 << " GetTimeStop() " << GetTimeStop() << endl
138 << " GetNTempVal() " << GetNTempVal() << endl
139 << " GetMinTemp() " << GetMinTemp() << endl
140 << " GetMaxTemp() " << GetMaxTemp() << endl;
141 // run ranges
142 cout << " RUN INFO: " << endl
143 << " NEvents " << GetNEvents() << endl
144 << " length (in hours) " << GetLengthOfRunInHours() << endl
145 << " range of temperature measurements (in hours) " << GetRangeOfTempMeasureInHours()
146 << " (in deg. C) " << GetRangeOfTempMeasureInDegrees()
147 << endl;
148 // range in correction values
149 double corrAtMinTemp = GetCorrection( fMinTemp );
150 double corrAtMaxTemp = GetCorrection( fMaxTemp );
151 double corrMaxMinDiff = 100*(corrAtMinTemp - corrAtMaxTemp);
152 cout << " CORRECTION INFO : " << endl
153 << " corrAtMinTemp " << corrAtMinTemp << endl
154 << " corrAtMaxTemp " << corrAtMaxTemp << endl
155 << " corrMaxMinDiff (~%) [=(corrAtMin - corrAtMax)*100] "
156 << corrMaxMinDiff << endl;
157
158 return;
159}
160//________________________________________________________________
161double AliEMCALCalibTestBeam::GetLengthOfRunInHours() const
162{
163 return (fTimeStop - fTimeStart)*fkSecToHour;
164}
165//________________________________________________________________
166double AliEMCALCalibTestBeam::GetRangeOfTempMeasureInHours() const
167{
168 return (fMaxTime - fMinTime)*fkSecToHour;
169}
170//________________________________________________________________
171double AliEMCALCalibTestBeam::GetRangeOfTempMeasureInDegrees() const
172{
173 return (fMaxTemp - fMinTemp);
174}
175
176//________________________________________________________________
177void AliEMCALCalibTestBeam::Init(const int runno)
178{
179 Reset(); // start fresh
180
181 // collect the needed information
182 GetRunTime(runno); // when was the run exactly..
183 GetTemperatureInfo(); // temperature readings during the run
184 GetEMCALLogbookInfo(runno); // logbook notes during the run
185
186 // make a smooth graph and spline, if we can
187 if (fTempGraph->GetN() > 0) { // need to have some points
188 TGraphSmooth smooth;
189 TGraph * graphSmoothed = smooth.SmoothKern(fTempGraph,"normal", fkSmoothBandwidth);
190 // use a cubic spline on the smoothed graph..
191 fTempSpline = new TSpline3("fTempSpline",graphSmoothed);
192 }
193 return;
194}
195
196//________________________________________________________________
197double AliEMCALCalibTestBeam::GetCorrection(int secSinceRunStart) const
198{
199 double temperature = GetTemperature(secSinceRunStart);
200 return GetCorrection(temperature);
201}
202
203//________________________________________________________________
204double AliEMCALCalibTestBeam::GetTemperature(int secSinceRunStart) const
205{
206 // first convert from seconds to hours..
207 double timeHour = secSinceRunStart * fkSecToHour;
208 return fTempSpline->Eval(timeHour);
209}
210
211//________________________________________________________________
212double AliEMCALCalibTestBeam::GetCorrection(double temperature) const
213{ // gain decreases with increasing temperature
214 double gain = (1.0 - (temperature-fkNormTemp)*fkTempSlope);
215 // i.e. multiplicative correction to ADC increases with increasing temperature
216 return 1.0/gain;
217}
218
219/* Next comes the method that does the work in picking up all the needed info..*/
220//________________________________________________________________
221void AliEMCALCalibTestBeam::GetRunTime(const int runno,
222 const char *filename)
223{
224 TFile *fR = new TFile(filename);
225 if (!fR) {
226 printf("file %s could not be found. Perhaps you don't have afs working?\n",filename);
227 printf("Then you can try to get it with:\n wget http://cern.ch/dsilverm/testbeam07/calib/daqLogbook.root");
228 }
229
230 TTree *daqLogbook = (TTree*)gDirectory->Get("daqLogbook");
231
232 //Declaration of leaves types
233 Int_t runNumber;
234 Int_t timeStart;
235 Int_t runDuration;
236 Int_t totalEvents;
237 Float_t totalData;
238 Float_t averageDataRate;
239 Float_t averageEventsPerSecond;
240
241 // Set branch addresses.
242 daqLogbook->SetBranchAddress("runNumber",&runNumber);
243 daqLogbook->SetBranchAddress("timeStart",&timeStart);
244 daqLogbook->SetBranchAddress("runDuration",&runDuration);
245 daqLogbook->SetBranchAddress("totalEvents",&totalEvents);
246 daqLogbook->SetBranchAddress("totalData",&totalData);
247 daqLogbook->SetBranchAddress("averageDataRate",&averageDataRate);
248 daqLogbook->SetBranchAddress("averageEventsPerSecond",&averageEventsPerSecond);
249
250 Long64_t nentries = daqLogbook->GetEntries();
251 Long64_t nbytes = 0;
252
253 for (Long64_t i=0; i<nentries;i++) {
254 nbytes += daqLogbook->GetEntry(i);
255 if (runNumber == runno) {
256 cout << " getRunTime: runNumber " << runNumber
257 << " timeStart " << timeStart
258 << " runDuration (seconds) " << runDuration
259 << " totalEvents " << totalEvents
260 << " totalData (kByte) " << totalData
261 << endl;
262
263 // assign the start and stop unix-timestamps
264 fTimeStart = timeStart;
265 fTimeStop = timeStart + runDuration;
266 fNEvents = totalEvents;
267 return; // then we are done..; should only pick up info for one run maximum..
268 }
269 }
270
271 return;
272}
273
274//________________________________________________________________
275void AliEMCALCalibTestBeam::GetTemperatureInfo(const char* filename)
276{
277 fTempGraph->Clear();
278
279 TFile *fT = new TFile(filename);
280 if (!fT) {
281 printf("file %s could not be found. Perhaps you don't have afs working?\n",filename);
282 printf("Then you can try to get it with:\n wget http://cern.ch/dsilverm/testbeam07/calib/temperature-merged-alldays.root");
283 }
284 TTree *tree = (TTree*)gDirectory->Get("tree");
285
286 //Declaration of leaves types
287 Int_t unixTime;
288 Float_t temperature;
289 Int_t SensorNumber;
290
291 // Set branch addresses.
292 tree->SetBranchAddress("unixTime",&unixTime);
293 tree->SetBranchAddress("temperature",&temperature);
294 tree->SetBranchAddress("SensorNumber",&SensorNumber);
295
296 Long64_t nentries = tree->GetEntries();
297 Long64_t nbytes = 0;
298
299 // init counter and range info
300 int np = 0; // number of points
301 fMinTemp = fkMaxValidTemp;
302 fMaxTemp = fkMinValidTemp;
303 fMinTime = 2000000000; // some timestamp far in the future..
304 fMaxTime = 0;
305
306 for (Long64_t i=0; i<nentries;i++) {
307 nbytes += tree->GetEntry(i);
308 if (unixTime>=fTimeStart && unixTime<=fTimeStop && SensorNumber==fkSensor) {
309 if (temperature>=fkMinValidTemp && temperature<=fkMaxValidTemp) {
310 /*
311 cout << " i " << i
312 << " unixTime " << unixTime
313 << " temperature " << temperature
314 << " SensorNumber " << SensorNumber
315 << endl;
316 */
317 if (fMinTemp > temperature) fMinTemp = temperature;
318 if (fMaxTemp < temperature) fMaxTemp = temperature;
319 if (fMinTime > unixTime) fMinTime = unixTime;
320 if (fMaxTime < unixTime) fMaxTime = unixTime;
321
322 fTempGraph->SetPoint(np, (unixTime-fTimeStart)*fkSecToHour, temperature);
323 np++;
324 } // valid temperature reading
325 } // valid time reading
326 }
327
328 // what if no values were found? Then we reset the max/min values etc. to the defaults
329 if (np==0) {
330 ResetTempInfo(); // partial Reset()
331 }
332
333 cout << " Number of temperature points found: np= " << np
334 << " for sensor " << fkSensor
335 << " : maxtemp " << fMaxTemp
336 << " mintemp " << fMinTemp
337 << endl;
338 return;
339}
340
341//________________________________________________________________
342void AliEMCALCalibTestBeam::GetEMCALLogbookInfo(const int runno,
343 const char *filename)
344{
345 Int_t Run,DAQEvt,CCUSB,Counter,XPos,YPos,Momentum;
346 char TriggerAndComment[1000];
347
348 ifstream fin;
349 fin.open(filename);
350 if(!fin.good()) {
351 printf("file %s could not be found. Perhaps you don't have afs working?\n",filename);
352 printf("Then you can try to get it with:\n wget http://cern.ch/dsilverm/testbeam07/calib/EMCAL_Logbook_SPS_and_PS.csv");
353 }
354
355 char line[1000];
356 int nlines=0;
357 while (!fin.getline(line,1000).eof()) {
358 if(!fin.good()) break;
359
360 sscanf(line,"%i,%i,%i,%i,%i,%i,%i,%s",
361 &Run,&DAQEvt,&CCUSB,&Counter,&XPos,&YPos,&Momentum,
362 TriggerAndComment);
363
364 if (Run == runno) { // we found the info
365 cout << " EMCAL Logbook info line: " << line << endl;
366 }
367
368 nlines++;
369 }
370 return;
371}
372
373