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