00001 #include <cmath>
00002
00003 #include "TCanvas.h"
00004 #include "TF1.h"
00005 #include "TH1F.h"
00006
00007 #include "MessageService/MsgService.h"
00008
00009 #include "ScintCal/ScintCalInfo.h"
00010
00011 ClassImp(ScintCalInfo)
00012
00013 CVSID("$Id: ScintCalInfo.cxx,v 1.6 2013/09/27 17:08:53 evansj Exp $");
00014
00015 using std::string;
00016
00017
00018 ScintCalInfo::ScintCalInfo() : TObject()
00019 {
00020 MSG("ScintCalInfo", Msg::kDebug)
00021 << "Default Info constructor"
00022 << endl;
00023 fepeak = 0.0;
00024 flength = 0.0;
00025 fname = "InfoObject";
00026 fnumhits = 0;
00027 fnummuons = 0;
00028 fnumtemps = 0;
00029 ftimeperiod = 0;
00030 ftotph = 0.0;
00031 ftotphsqrd = 0.0;
00032 ftottemp = 0.0;
00033 ftottempsqrd = 0.0;
00034 fpeak = 0.0;
00035
00036 fnumbins = 20000;
00037 flowmed = 0;
00038 fhighmed = 2000;
00039 string histname = "fhPH" + fname + Form("%d",ftimeperiod);
00040 fhPH = new TH1F(histname.c_str(),
00041 "Every PH in this time period",
00042 fnumbins, flowmed, fhighmed);
00043
00044 fhistogramFitted = false;
00045 }
00046
00047
00048 ScintCalInfo::ScintCalInfo(string name) : TObject()
00049 {
00050 MSG("ScintCalInfo", Msg::kDebug)
00051 << "Default Info constructor"
00052 << endl;
00053 fepeak = 0.0;
00054 flength = 0.0;
00055 fname = name;
00056 fnumhits = 0;
00057 fnummuons = 0;
00058 fnumtemps = 0;
00059 ftimeperiod = 0;
00060 ftotph = 0.0;
00061 ftotphsqrd = 0.0;
00062 ftottemp = 0.0;
00063 ftottempsqrd = 0.0;
00064 fpeak = 0.0;
00065
00066 fnumbins = 20000;
00067 flowmed = 0;
00068 fhighmed = 2000;
00069 string histname = "fhPH" + fname + Form("%d",ftimeperiod);
00070 fhPH = new TH1F(histname.c_str(),
00071 "Every PH in this time period",
00072 fnumbins, flowmed, fhighmed);
00073
00074 fhistogramFitted = false;
00075 }
00076
00077
00078 ScintCalInfo::ScintCalInfo(const UInt_t timeperiod,
00079 const Double_t length,
00080 const string name) : TObject()
00081 {
00082
00083 MSG("ScintCalInfo", Msg::kDebug)
00084 << "Creating info object for timeperiod "
00085 << timeperiod
00086 << endl;
00087 fepeak = 0.0;
00088 flength = length;
00089 fname = name;
00090 fnumhits = 0;
00091 fnummuons = 0;
00092 fnumtemps = 0;
00093 ftimeperiod = timeperiod;
00094 ftotph = 0.0;
00095 ftotphsqrd = 0.0;
00096 ftottemp = 0.0;
00097 ftottempsqrd = 0.0;
00098 fpeak = 0.0;
00099
00100 fnumbins = 20000;
00101 flowmed = 0;
00102 fhighmed = 2000;
00103 string histname = "fhPH" + fname + Form("%d",ftimeperiod);
00104 fhPH = new TH1F(histname.c_str(),
00105 "Every PH in this time period",
00106 fnumbins, flowmed, fhighmed);
00107
00108 fhistogramFitted = false;
00109 }
00110
00111
00112 ScintCalInfo::ScintCalInfo(const ScintCalInfo& source) : TObject()
00113 {
00114 MSG("ScintCalInfo", Msg::kDebug)
00115 << "Copy constructing info object for timeperiod "
00116 << source.ftimeperiod
00117 << endl;
00118 fepeak = source.fepeak;
00119 flength = source.flength;
00120 fname = source.fname;
00121 fnumhits = source.fnumhits;
00122 fnummuons = source.fnummuons;
00123 fnumtemps = source.fnumtemps;
00124 ftimeperiod = source.ftimeperiod;
00125 ftotph = source.ftotph;
00126 ftotphsqrd = source.ftotphsqrd;
00127 ftottemp = source.ftottemp;
00128 ftottempsqrd = source.ftottempsqrd;
00129 fpeak = source.fpeak;
00130
00131
00132
00133
00134
00135 fnumbins = source.fnumbins;
00136 flowmed = source.flowmed;
00137 fhighmed = source.fhighmed;
00138
00139
00140
00141
00142
00143 if (source.fhPH){fhPH = new TH1F(*(source.fhPH));}
00144
00145
00146
00147
00148
00149 fhistogramFitted = source.fhistogramFitted;
00150 }
00151
00152
00153 ScintCalInfo::~ScintCalInfo()
00154 {
00155 if (fhPH){delete fhPH; fhPH = 0;}
00156 }
00157
00158
00159 ScintCalInfo& ScintCalInfo::operator = (const ScintCalInfo& calInfo)
00160 {
00161 MSG("ScintCalInfo", Msg::kDebug)
00162 << "Overload =ing info object for timeperiod "
00163 << calInfo.ftimeperiod
00164 << endl;
00165 if (this != &calInfo){
00166 fepeak = calInfo.fepeak;
00167 flength = calInfo.flength;
00168 fname = calInfo.fname;
00169 fnumhits = calInfo.fnumhits;
00170 fnummuons = calInfo.fnummuons;
00171 fnumtemps = calInfo.fnumtemps;
00172 ftimeperiod = calInfo.ftimeperiod;
00173 ftotph = calInfo.ftotph;
00174 ftotphsqrd = calInfo.ftotphsqrd;
00175 ftottemp = calInfo.ftottemp;
00176 ftottempsqrd = calInfo.ftottempsqrd;
00177 fpeak = calInfo.fpeak;
00178
00179 fnumbins = calInfo.fnumbins;
00180 flowmed = calInfo.flowmed;
00181 fhighmed = calInfo.fhighmed;
00182
00183 if (fhPH) {delete fhPH; fhPH = 0;}
00184 if (calInfo.fhPH) {fhPH = new TH1F(*(calInfo.fhPH));}
00185
00186 fhistogramFitted = calInfo.fhistogramFitted;
00187 }
00188 return *this;
00189 }
00190
00191
00192 ScintCalInfo ScintCalInfo::operator+ (const ScintCalInfo& calInfo)
00193 {
00194 MSG("ScintCalInfo", Msg::kDebug)
00195 << "+operator"
00196 << endl;
00197 ScintCalInfo temp;
00198
00199 if (ftimeperiod != calInfo.ftimeperiod){
00200 MSG("ScintCalInfo", Msg::kWarning)
00201 << "Trying to amalgamate objects from different time periods"
00202 << endl;
00203 return temp;
00204 }
00205 if (flength != calInfo.flength){
00206 MSG("ScintCalInfo", Msg::kWarning)
00207 << "Trying to amalgamate objects with different lengths"
00208 << endl;
00209 return temp;
00210 }
00211 if (fnumbins != calInfo.fnumbins
00212 || flowmed != calInfo.flowmed
00213 || fhighmed != calInfo.fhighmed){
00214 MSG("ScintCalInfo", Msg::kWarning)
00215 << "Trying to amalgamate objects with incompatible histograms"
00216 << endl;
00217 return temp;
00218 }
00219
00220 temp.fepeak = 0.0;
00221 temp.flength = calInfo.flength;
00222 temp.fname = fname;
00223 temp.fnumhits = fnumhits + calInfo.fnumhits;
00224 temp.fnummuons = fnummuons + calInfo.fnummuons;
00225 temp.fnumtemps = fnumtemps + calInfo.fnumtemps;
00226 temp.ftimeperiod = calInfo.ftimeperiod;
00227 temp.ftotph = ftotph + calInfo.ftotph;
00228 temp.ftotphsqrd = ftotphsqrd + calInfo.ftotphsqrd;
00229 temp.ftottemp = ftottemp + calInfo.ftottemp;
00230 temp.ftottempsqrd = ftottempsqrd + calInfo.ftottempsqrd;
00231 temp.fpeak = 0.0;
00232
00233 temp.fnumbins = fnumbins;
00234 temp.flowmed = flowmed;
00235 temp.fhighmed = fhighmed;
00236
00237 if(temp.fhPH) {delete temp.fhPH; temp.fhPH = 0;}
00238 temp.fhPH = new TH1F(*fhPH);
00239 temp.fhPH->Add(calInfo.fhPH);
00240
00241 temp.fhistogramFitted = false;
00242 return temp;
00243 }
00244
00245
00246 void ScintCalInfo::DrawHist() const
00247 {
00248 TCanvas cPH("cPH",
00249 Form("PH for timeperiod %d",ftimeperiod),
00250 200,10,700,500);
00251 MSG("ScintCalDatabase",Msg::kInfo)
00252 << "Canvas made"
00253 << endl;
00254 fhPH->Draw();
00255 MSG("ScintCalDatabase",Msg::kInfo)
00256 << "Histogram on canvas"
00257 << endl;
00258 cPH.Print("BadBadBad.C");
00259 return;
00260 }
00261
00262
00263 void ScintCalInfo::DisassociateHist() const
00264 {
00265 fhPH->SetDirectory(0);
00266 return;
00267 }
00268
00269
00270 Float_t ScintCalInfo::Length() const
00271 {
00272 return flength;
00273 }
00274
00275
00276 Double_t ScintCalInfo::EMeanTemperature() const
00277 {
00278
00279 if (fnumtemps){
00280 Double_t tmean = this->MeanTemperature();
00281 Double_t tRMS = ftottempsqrd/((Double_t) fnumtemps);
00282 if (tRMS - tmean*tmean < 0){
00283 MSG("ScintCalInfo", Msg::kWarning)
00284 << "Temperatures in bin " << ftimeperiod
00285 << ": RMS is " << tRMS
00286 << " whilst mean^2 is " << tmean*tmean
00287 << ". This makes sigma imaginary."
00288 << endl;
00289 return 0.0;
00290 }
00291 Double_t tsigma = sqrt(tRMS - tmean*tmean);
00292 return tsigma/sqrt((Double_t) fnumtemps);
00293 }
00294 else {return 0.0;}
00295 }
00296
00297
00298 Double_t ScintCalInfo::EMean() const
00299 {
00300 return SigOvRtN();
00301 }
00302
00303
00304 Double_t ScintCalInfo::EMedian() const
00305 {
00306 return SigOvRtN();
00307 }
00308
00309
00310 Double_t ScintCalInfo::EPeak()
00311 {
00312 if (!fhistogramFitted){this->FitHistogram();}
00313 return fepeak;
00314 }
00315
00316
00317 void ScintCalInfo::FitHistogram()
00318 {
00319 MSG("ScintCalInfo", Msg::kInfo)
00320 << "Fitting histogram for bin "
00321 << ftimeperiod
00322 << endl;
00323 if (!fnumhits){
00324 fepeak = 0.0;
00325 fpeak = 0.0;
00326 fhistogramFitted = true;
00327 return;
00328 }
00329
00330 TCanvas cinfo("cinfo",
00331 Form("Info for timeperiod %d",ftimeperiod),
00332 200,10,700,500);
00333
00334 TH1* hFit = fhPH->Rebin(50,"hFit");
00335 hFit->Fit("gaus","","",300.0,500.0);
00336 TF1* fit = hFit->GetFunction("gaus");
00337 fpeak = fit->GetParameter(1);
00338
00339 fepeak = fit->GetParError(2);
00340 if (fit) {delete fit; fit = 0;}
00341 if (hFit) {delete hFit; hFit = 0;}
00342
00343 fhistogramFitted = true;
00344 return;
00345 }
00346
00347
00348 Double_t ScintCalInfo::Mean() const
00349 {
00350 if (!fnumhits){return 0.0;}
00351 return ftotph/((Double_t) fnumhits);
00352 }
00353
00354
00355 Double_t ScintCalInfo::MeanTemperature() const
00356 {
00357 if (fnumtemps){return ftottemp/((Double_t) fnumtemps);}
00358 else {return -1.0;}
00359 }
00360
00361
00362 Double_t ScintCalInfo::Median() const
00363 {
00364 Double_t median = 0.0;
00365 Double_t runningcount = 0.0;
00366 Bool_t underflowbin = true;
00367 Bool_t overflowbin = false;
00368 for (Int_t i=0; i <= fnumbins+1; ++i){
00369 if (i > 0){underflowbin = false;}
00370 runningcount += fhPH->GetBinContent(i);
00371 if (runningcount >= ((Double_t) fnumhits)/2.0 ){
00372 median = flowmed +
00373 ( ((Double_t) i)-0.5 )*(fhighmed-flowmed)/( (Double_t) fnumbins);
00374 break;
00375 }
00376 if (fnumbins+1 == i){overflowbin = true;}
00377 }
00378
00379 if(underflowbin){
00380 MAXMSG("ScintCalInfo", Msg::kWarning,1000)
00381 << "Median in underflow bin."
00382 << endl;
00383 median = -1.0;
00384 }
00385 if (overflowbin){
00386 MSG("ScintCalInfo", Msg::kWarning)
00387 << "Median in overflow bin."
00388 << endl;
00389 median = -1.0;
00390 }
00391 return median;
00392 }
00393
00394
00395 void ScintCalInfo::NewMuon()
00396 {
00397 ++fnummuons;
00398 return;
00399 }
00400
00401
00402 UInt_t ScintCalInfo::NumPlanes() const
00403 {
00404 return fnumhits;
00405 }
00406
00407
00408 UInt_t ScintCalInfo::NumMuons() const
00409 {
00410 return fnummuons;
00411 }
00412
00413
00414 Double_t ScintCalInfo::Peak()
00415 {
00416 if (!fhistogramFitted){this->FitHistogram();}
00417 return fpeak;
00418 }
00419
00420
00421 Double_t ScintCalInfo::Quantile50()
00422 {
00423 Double_t lowQuant = 0.0;
00424 Double_t highQuant = 0.0;
00425 Double_t runningcount = 0.0;
00426 Bool_t underflowbin = true;
00427 Bool_t overflowbin = false;
00428 for (Int_t i=0; i <= fnumbins+1; ++i){
00429 runningcount += fhPH->GetBinContent(i);
00430
00431 if (runningcount >= ((Double_t) fnumhits)/4.0 ){
00432 lowQuant = flowmed +
00433 ( ((Double_t) i)-0.5 )*(fhighmed-flowmed)/( (Double_t) fnumbins);
00434 if (i > 0){underflowbin = false;}
00435 }
00436 if (fnumbins+1 == i){overflowbin = true;}
00437
00438 if (runningcount >= 3.0*((Double_t) fnumhits)/4.0 ){
00439 highQuant = flowmed +
00440 ( ((Double_t) i)-0.5 )*(fhighmed-flowmed)/( (Double_t) fnumbins);
00441 break;
00442 }
00443 }
00444
00445 Double_t quantile50 = highQuant-lowQuant;
00446
00447
00448 if(underflowbin){
00449 MAXMSG("ScintCalInfo", Msg::kWarning,1000)
00450 << "Quantile in underflow bin."
00451 << endl;
00452 quantile50 = -1.0;
00453 }
00454 if (overflowbin){
00455 MSG("ScintCalInfo", Msg::kWarning)
00456 << "Median in overflow bin."
00457 << endl;
00458 quantile50 = -1.0;
00459 }
00460
00461 return quantile50;
00462 }
00463
00464
00465 Double_t ScintCalInfo::RMS() const
00466 {
00467 if (!fnumhits){return 0.0;}
00468 return ftotphsqrd/((Double_t) fnumhits);
00469 }
00470
00471
00472 void ScintCalInfo::SetBinLength(const Double_t binLength)
00473 {
00474 flength = binLength;
00475 return;
00476 }
00477
00478
00479 void ScintCalInfo::SetTimePeriod(const UInt_t timeperiod)
00480 {
00481 ftimeperiod = timeperiod;
00482 return;
00483 }
00484
00485
00486 Double_t ScintCalInfo::Sigma() const
00487 {
00488 if (!fnumhits){return 0.0;}
00489
00490 Double_t mean = this->Mean();
00491 Double_t RMS = this->RMS();
00492 if (RMS - mean*mean < 0){
00493 MSG("ScintCalInfo", Msg::kWarning)
00494 << "Pulse heights in bin " << ftimeperiod
00495 << ": RMS is " << RMS
00496 << " whilst mean^2 is " << mean*mean
00497 << ". This makes sigma imaginary."
00498 << endl;
00499 return 0.0;
00500 }
00501
00502 return sqrt(RMS - mean*mean);
00503 }
00504
00505
00506 Double_t ScintCalInfo::SigOvRtN() const
00507 {
00508 if (!fnumhits){return 0.0;}
00509
00510 Double_t sigma = this->Sigma();
00511 return sigma/sqrt((Double_t) fnumhits);
00512 }
00513
00514
00515
00516 void ScintCalInfo::SingleHitUpdate(const Float_t ph)
00517 {
00518 Double_t phd = (Double_t) ph;
00519
00520 fnumhits++;
00521 ftotph += phd;
00522 ftotphsqrd += (phd*phd);
00523 fhPH->Fill(ph);
00524
00525 fhistogramFitted = false;
00526 }
00527
00528
00529 void ScintCalInfo::SingleTempUpdate(const Float_t temp)
00530 {
00531 Double_t tempd = (Double_t) temp;
00532 fnumtemps++;
00533 ftottemp += tempd;
00534 ftottempsqrd += (tempd*tempd);
00535 }
00536
00537
00538 UInt_t ScintCalInfo::TimePeriod() const
00539 {
00540 return ftimeperiod;
00541 }