Add more user oriented options to costumize the tender - Jiri
[u/mrichter/AliRoot.git] / EMCAL / beamtest07 / AliEMCALCalibTestBeam.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: 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.. */
34 const double fkSmoothBandwidth = 5.0; // 5 seems to be a good value..; see $ROOTSYS/tutorials/graphs/motorcycle.C for reference
35 const double fkSecToHour = 1.0/3600.0; // conversion factor from seconds to hours
36
37 // some global variables for APD handling
38 const double fkTempSlope = 0.017; // 1.7% per deg. C, seems about right for all APDs, from studies at Catania, and by Rachid
39 const 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
42 const 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..
44 const int fkMinValidTemp = 0;
45 const int fkMaxValidTemp = 100;
46
47 using namespace std;
48
49 ClassImp(AliEMCALCalibTestBeam)
50
51 //________________________________________________________________
52 AliEMCALCalibTestBeam::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 //________________________________________________________________
62 AliEMCALCalibTestBeam::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 //________________________________________________________________
77 AliEMCALCalibTestBeam &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 //________________________________________________________________
88 AliEMCALCalibTestBeam::~AliEMCALCalibTestBeam()
89 {
90   // Destructor
91   if (fTempGraph) delete fTempGraph;
92   if (fTempSpline) delete fTempSpline;
93 }
94
95 //________________________________________________________________
96 void  AliEMCALCalibTestBeam::Reset() 
97 { // clear variables
98   ResetRunInfo();
99   ResetTempInfo();
100   return;
101 }
102 //________________________________________________________________
103 void  AliEMCALCalibTestBeam::ResetRunInfo() 
104 {
105   // clear variables
106   fTimeStart = 0;
107   fTimeStop = 0;
108   fNEvents = 0;
109 }
110
111 //________________________________________________________________
112 void  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 //________________________________________________________________
130 void  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 //________________________________________________________________ 
161 double AliEMCALCalibTestBeam::GetLengthOfRunInHours() const
162 { // run length in hours
163   return (fTimeStop - fTimeStart)*fkSecToHour;
164 }
165 //________________________________________________________________ 
166 double AliEMCALCalibTestBeam::GetRangeOfTempMeasureInHours() const
167 { // interval with temperature measurements
168   return (fMaxTime - fMinTime)*fkSecToHour;
169 }
170 //________________________________________________________________ 
171 double AliEMCALCalibTestBeam::GetRangeOfTempMeasureInDegrees() const
172 { // range of temperatures
173   return (fMaxTemp - fMinTemp);
174 }
175
176 //________________________________________________________________
177 void AliEMCALCalibTestBeam::Init(const int runno)
178 { // init/clear info
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 //________________________________________________________________
197 double AliEMCALCalibTestBeam::GetCorrection(int secSinceRunStart) const
198 { // calculate correction
199   double temperature = GetTemperature(secSinceRunStart);
200   return GetCorrection(temperature);
201 }
202
203 //________________________________________________________________
204 double 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 //________________________________________________________________
212 double 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 //________________________________________________________________
221 void AliEMCALCalibTestBeam::GetRunTime(const int runno, 
222                                        const char *filename)
223 { // get run info
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 //________________________________________________________________
275 void AliEMCALCalibTestBeam::GetTemperatureInfo(const char* filename)
276 { // get temperature info
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 //________________________________________________________________
342 void AliEMCALCalibTestBeam::GetEMCALLogbookInfo(const int runno,                
343                                                 const char *filename)
344 { // get logbook info
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