loon: ScintCalInfo.cxx Source File

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 //____________________________________________________________________72
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 //____________________________________________________________________72
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 //____________________________________________________________________72
00078 ScintCalInfo::ScintCalInfo(const UInt_t timeperiod,
00079                            const Double_t length,
00080                            const string name) : TObject()
00081 {
00082   //Accepts length in seconds.
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 //___________________ _________________________________________________72
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 //   MSG("ScintCalInfo", Msg::kInfo) 
00132 //     << "Info::Copy1 "
00133 //     << source.ftimeperiod
00134 //     << endl;
00135   fnumbins = source.fnumbins;
00136   flowmed = source.flowmed;
00137   fhighmed = source.fhighmed;
00138 //   source.DrawHist();
00139 //   MSG("ScintCalInfo", Msg::kInfo) 
00140 //     << "Info::Copy2 "
00141 //     << source.ftimeperiod
00142 //     << endl;
00143   if (source.fhPH){fhPH = new TH1F(*(source.fhPH));}
00144 //   MSG("ScintCalInfo", Msg::kInfo) 
00145 //     << "Info::Copy3 "
00146 //     << source.ftimeperiod
00147 //     << endl;
00148 
00149   fhistogramFitted = source.fhistogramFitted;
00150 }
00151 
00152 //____________________________________________________________________72
00153 ScintCalInfo::~ScintCalInfo()
00154 {
00155   if (fhPH){delete fhPH; fhPH = 0;}
00156 }
00157 
00158 //____________________________________________________________________72
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 //____________________________________________________________________72
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  //____________________________________________________________________72
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 //____________________________________________________________________72
00263 void ScintCalInfo::DisassociateHist() const
00264 {
00265   fhPH->SetDirectory(0);
00266   return;
00267 }
00268 
00269 //____________________________________________________________________72
00270 Float_t ScintCalInfo::Length() const
00271 {
00272   return flength;
00273 }
00274 
00275 //____________________________________________________________________72
00276 Double_t ScintCalInfo::EMeanTemperature() const
00277 {
00278   //Returns sigma/sqrt(n) for the temperature entries.
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 //____________________________________________________________________72
00298 Double_t ScintCalInfo::EMean() const
00299 {
00300   return SigOvRtN();
00301 }
00302 
00303 //____________________________________________________________________72
00304 Double_t ScintCalInfo::EMedian() const
00305 {
00306   return SigOvRtN();
00307 }
00308 
00309 //____________________________________________________________________72
00310 Double_t ScintCalInfo::EPeak()
00311 {
00312   if (!fhistogramFitted){this->FitHistogram();}
00313   return fepeak;
00314 }
00315 
00316 //____________________________________________________________________72
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   //Now calculate the truncated mean with a gaussian fit of the peak
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   //  fepeak = fit->GetParameter(2)/sqrt((Double_t) fnumhits);
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 //____________________________________________________________________72
00348 Double_t ScintCalInfo::Mean() const
00349 {
00350   if (!fnumhits){return 0.0;}
00351   return ftotph/((Double_t) fnumhits);
00352 }
00353 
00354 //____________________________________________________________________72
00355 Double_t ScintCalInfo::MeanTemperature() const
00356 {
00357   if (fnumtemps){return ftottemp/((Double_t) fnumtemps);}
00358   else {return -1.0;}
00359 }
00360 
00361 //____________________________________________________________________72
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   //Deal with median in underflow & overflow bins
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 //____________________________________________________________________72
00395 void ScintCalInfo::NewMuon()
00396 {
00397   ++fnummuons;
00398   return;
00399 }
00400 
00401 //____________________________________________________________________72
00402 UInt_t ScintCalInfo::NumPlanes() const
00403 {
00404   return fnumhits;
00405 }
00406 
00407 //____________________________________________________________________72
00408 UInt_t ScintCalInfo::NumMuons() const
00409 {
00410   return fnummuons;
00411 }
00412 
00413 //____________________________________________________________________72
00414 Double_t ScintCalInfo::Peak()
00415 {
00416   if (!fhistogramFitted){this->FitHistogram();}
00417   return fpeak;
00418 }
00419 
00420 //____________________________________________________________________72
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   //Deal with quantiles in underflow & overflow bins
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 //____________________________________________________________________72
00465 Double_t ScintCalInfo::RMS() const
00466 {
00467   if (!fnumhits){return 0.0;}
00468   return ftotphsqrd/((Double_t) fnumhits);
00469 }
00470 
00471 //____________________________________________________________________72
00472 void ScintCalInfo::SetBinLength(const Double_t binLength)
00473 {
00474   flength = binLength;
00475   return;
00476 }
00477 
00478 //____________________________________________________________________72
00479 void ScintCalInfo::SetTimePeriod(const UInt_t timeperiod)
00480 {
00481   ftimeperiod = timeperiod;
00482   return;
00483 }
00484 
00485 //____________________________________________________________________72
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 //____________________________________________________________________72
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 //____________________________________________________________________72
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 //____________________________________________________________________72
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 //____________________________________________________________________72
00538 UInt_t ScintCalInfo::TimePeriod() const
00539 {
00540   return ftimeperiod;
00541 }