FitOverROI.C

00001 #include "TH1.h"
00002 #include "TF1.h"
00003 #include "TMath.h"
00004 #include "TFitResultPtr.h"
00005 #include "TFitResult.h"
00006 #include "TSpectrum.h"
00007 #include "TTree.h"
00008 #include "TROOT.h"
00009 #include "TCanvas.h"
00010 #include "TPad.h"
00011 #include "TLine.h"
00012 #include "TAxis.h"
00013 #include "TFile.h"
00014 #include "TList.h"
00015 #include "TCut.h"
00016 #include <algorithm>
00017 #include <numeric>
00018 #include <string>
00019 #include <iomanip>
00020 #include <iostream>
00021 #include <sstream>
00022 #include "TMatrixDSym.h"
00023 #include "TGraph.h"
00024 #include <fstream>
00025 #include "ChanFitSettings.hh"
00026 
00027 using namespace std;
00028 
00029 #define HISTOGRAMWIDTH 7
00030 #define NPEAKS 7
00031 
00032 enum PARAMETERS {CONSTANT, LAMBDA, MEAN, SIGMA, AMP_E, P_E , SHOTNOISE, PEDMEAN, SHOTNOISE2, PEDWEIGHT,  NPAR };
00033 
00034 #define mCON params[CONSTANT]
00035 #define mLAM params[LAMBDA]
00036 #define mMEAN params[MEAN]
00037 #define mSIG params[SIGMA]
00038 #define mAMP params[AMP_E]
00039 #define mPE params[P_E]
00040 #define mSHT params[SHOTNOISE]
00041 #define mPDM params[PEDMEAN]
00042 #define mSHT2 params[SHOTNOISE2]
00043 #define mPDWT params[PEDWEIGHT]
00044 
00045 
00046 const char* names[] = {"CONSTANT","LAMBDA", "SPE_MEAN", "SPE_SIGMA", "AMP_E", "P_E", "SHOTNOISE" ,"PEDMEAN","SHOTNOISE2","PEDWEIGHT","NPAR"};
00047 
00048 
00049 
00050 
00051 Double_t background_funca(Double_t* x, Double_t* params)
00052 {
00053         double y = x[0] - mPDM;
00054         /*double exp_term= mPE/mAMP * 
00055                 TMath::Power(TMath::E(), (mSHT*mSHT - 2* mAMP *x[0])/(2 * mAMP*mAMP)) *
00056                 (1+ TMath::Erf(-mSHT*mSHT/(sqrt(2) * mAMP * mSHT)  )  );
00057                 */
00058                 
00059                 //mathematica version (my own convolution)
00060         double exp_term=(TMath::Power(TMath::E(),(-2*mAMP*y + mSHT*mSHT)/
00061                                         (2.*mAMP*mAMP))*mPE*
00062                                         (1 + TMath::Erf(((mAMP*y)/mSHT - mSHT)/(sqrt(2)*mAMP))))/
00063                                         (2.*mAMP);
00064         return mCON*TMath::Poisson(1,mLAM)*exp_term;
00065 }
00066 
00067 Double_t background_funcb(Double_t* x, Double_t* params)
00068 {
00069         double y = x[0] - mPDM;
00070         /*double exp_term= mPE/mAMP * 
00071                 TMath::Power(TMath::E(), (mSHT*mSHT - 2* mAMP *x[0])/(2 * mAMP*mAMP)) *
00072                 (1+ TMath::Erf(-mSHT*mSHT/(sqrt(2) * mAMP * mSHT)  )  );
00073                 */
00074                 
00075                 //mathematica version (my own convolution)
00076         double exp_term=(TMath::Power(TMath::E(),(-2*mAMP*y + mSHT2*mSHT2)/
00077                                         (2.*mAMP*mAMP))*mPE*
00078                                         (1 + TMath::Erf(((mAMP*y)/mSHT2 - mSHT2)/(sqrt(2)*mAMP))))/
00079                                         (2.*mAMP);
00080         return mCON*TMath::Poisson(1,mLAM)*exp_term;
00081 }
00082 
00083 Double_t background_func(Double_t* x, Double_t* params) {
00084   return mPDWT*background_funca(x, params) + (1-mPDWT)*background_funcb(x, params);
00085 }
00086 
00087 Double_t response_0a(Double_t* x, Double_t* params) {
00088   double y = x[0] - mPDM;
00089   return mPDWT*mCON*TMath::Poisson(0,mLAM)*TMath::Gaus(y, 0, mSHT, true);
00090 }
00091 
00092 Double_t response_0b(Double_t* x, Double_t* params) {
00093   double y = x[0] - mPDM;
00094   return (1-mPDWT)*mCON*TMath::Poisson(0,mLAM)*TMath::Gaus(y, 0, mSHT2, true);
00095 } //NOTE: response_0a is modified to include mPDWT
00096 
00097 Double_t response_0(Double_t* x, Double_t* params) {
00098   return response_0a(x, params) + response_0b(x, params);
00099 }
00100 
00101 Double_t response_1a(Double_t* x, Double_t* params)
00102 {
00103         double y = x[0] - mPDM;
00104         double exp_term= background_func(x, params);
00105         double sigma_1=sqrt(mSHT*mSHT + mSIG*mSIG);     
00106         
00107         double gauss_term= (1-mPE)*(TMath::Power(TMath::E(), - (mMEAN-y)*(mMEAN-y)/(2*sigma_1*sigma_1)) *
00108                 (1 + TMath::Erf((mMEAN*mSHT*mSHT + mSIG*mSIG*y)/(sqrt(2) * mSIG*mSHT *sigma_1))) /
00109                 (sqrt(2*TMath::Pi())*sigma_1 * (1+ TMath::Erf(mMEAN/(sqrt(2)*mSIG)))));
00110         
00111         
00112         double response=mCON*TMath::Poisson(1,mLAM)*gauss_term+exp_term;
00113         return response;
00114 }
00115 
00116 Double_t response_1b(Double_t* x, Double_t* params)
00117 {
00118         double y = x[0] - mPDM;
00119         double exp_term= background_func(x, params);
00120         double sigma_1=sqrt(mSHT2*mSHT2 + mSIG*mSIG);   
00121         
00122         double gauss_term= (1-mPE)*(TMath::Power(TMath::E(), - (mMEAN-y)*(mMEAN-y)/(2*sigma_1*sigma_1)) *
00123                 (1 + TMath::Erf((mMEAN*mSHT2*mSHT2 + mSIG*mSIG*y)/(sqrt(2) * mSIG*mSHT2 *sigma_1))) /
00124                 (sqrt(2*TMath::Pi())*sigma_1 * (1+ TMath::Erf(mMEAN/(sqrt(2)*mSIG)))));
00125         
00126         
00127         double response=mCON*TMath::Poisson(1,mLAM)*gauss_term+exp_term;
00128         return response;
00129 }
00130 
00131 Double_t response_1(Double_t* x, Double_t* params) {
00132   return mPDWT*response_1a(x, params) + (1-mPDWT)*response_1b(x, params);
00133 }
00134 
00135 
00136 
00137 Double_t response_2a(Double_t* x, Double_t* params)
00138 {
00139         double y = x[0] - mPDM;
00140         double response=(TMath::Power(TMath::E(),(-2*mAMP*y + TMath::Power(mSHT,2))/
00141         (2.*TMath::Power(mAMP,2)))*mPE*mPE*
00142                 (mAMP*y - TMath::Power(mSHT,2)))/TMath::Power(mAMP,3) + 
00143                 (2*TMath::Power(-1 + mPE,2)*sqrt(2/TMath::Pi()))/
00144                 (TMath::Power(TMath::E(),TMath::Power(-2*mMEAN + y,2)/
00145         (2.*(TMath::Power(mSHT,2) + 
00146                 2*TMath::Power(mSIG,2))))*mSHT*
00147                 sqrt(2/TMath::Power(mSHT,2) + TMath::Power(mSIG,-2))*
00148                 mSIG*TMath::Power(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG)),2)
00149                 ) + (2*TMath::Power(TMath::E(),
00150                 (2*mAMP*mMEAN - 2*mAMP*y + TMath::Power(mSHT,2) + 
00151         TMath::Power(mSIG,2))/(2.*TMath::Power(mAMP,2)))*
00152                 (-1 + mPE)*mPE*
00153                 (-1 + TMath::Erf((mAMP*(mMEAN - y) + 
00154         TMath::Power(mSHT,2) + TMath::Power(mSIG,2))/
00155         (mAMP*sqrt(2*TMath::Power(mSHT,2) + 
00156         2*TMath::Power(mSIG,2))))
00157         ))/(mAMP*(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG))));
00158         return TMath::Poisson(2,mLAM)*mCON*response;
00159         /*double sigma_1=sqrt(mSHT*mSHT + mSIG*mSIG);   
00160 
00161         double response= TMath::Poisson(2,mLAM)*(mPE*mPE*y/(mAMP*mAMP) * TMath::Power(TMath::E(), -y/mAMP)+ 
00162                 2 * (1-mPE)*mPE/(sqrt(2*TMath::Pi())*sigma_1) * TMath::Power(TMath::E(), -0.5*(y-mMEAN-mAMP)*(y-mMEAN-mAMP)/(sigma_1*sigma_1)) + 
00163                 (1-mPE)*(1-mPE)/(2 *sqrt(TMath::Pi())*sigma_1) *TMath::Power(TMath::E(),-0.5 *(y- 2*mMEAN)*(y- 2*mMEAN)/(sigma_1 * sigma_1 *2)));
00164         return mCON*response;*/
00165 }
00166         
00167 Double_t response_2b(Double_t* x, Double_t* params)
00168 {
00169         double y = x[0] - mPDM;
00170         double response=(TMath::Power(TMath::E(),(-2*mAMP*y + TMath::Power(mSHT2,2))/
00171         (2.*TMath::Power(mAMP,2)))*mPE*mPE*
00172                 (mAMP*y - TMath::Power(mSHT2,2)))/TMath::Power(mAMP,3) + 
00173                 (2*TMath::Power(-1 + mPE,2)*sqrt(2/TMath::Pi()))/
00174                 (TMath::Power(TMath::E(),TMath::Power(-2*mMEAN + y,2)/
00175         (2.*(TMath::Power(mSHT2,2) + 
00176                 2*TMath::Power(mSIG,2))))*mSHT2*
00177                 sqrt(2/TMath::Power(mSHT2,2) + TMath::Power(mSIG,-2))*
00178                 mSIG*TMath::Power(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG)),2)
00179                 ) + (2*TMath::Power(TMath::E(),
00180                 (2*mAMP*mMEAN - 2*mAMP*y + TMath::Power(mSHT2,2) + 
00181         TMath::Power(mSIG,2))/(2.*TMath::Power(mAMP,2)))*
00182                 (-1 + mPE)*mPE*
00183                 (-1 + TMath::Erf((mAMP*(mMEAN - y) + 
00184         TMath::Power(mSHT2,2) + TMath::Power(mSIG,2))/
00185         (mAMP*sqrt(2*TMath::Power(mSHT2,2) + 
00186         2*TMath::Power(mSIG,2))))
00187         ))/(mAMP*(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG))));
00188         return TMath::Poisson(2,mLAM)*mCON*response;
00189         
00190 }
00191 
00192 Double_t response_2(Double_t* x, Double_t* params) {
00193   return mPDWT*response_2a(x, params) + (1-mPDWT)*response_2b(x, params);
00194 }
00195 
00196 Double_t m_n(Double_t*x, Double_t* params)
00197 {
00198 return mMEAN + mAMP*mPE - mMEAN*mPE + 
00199    (sqrt(2/TMath::Pi())*mSIG*
00200       (1/(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG))) + 
00201         mPE/(-2 + TMath::Erfc(mMEAN/(sqrt(2)*mSIG)))))/
00202     TMath::Power(TMath::E(),TMath::Power(mMEAN,2)/(2.*TMath::Power(mSIG,2)));
00203 }
00204 
00205 Double_t sigma_n(Double_t*x, Double_t* params){
00206 return sqrt(-(TMath::Power(mMEAN,2)*(-1 + mPE)) + 
00207     2*TMath::Power(mAMP,2)*mPE  - 
00208     (-1 + mPE)*TMath::Power(mSIG,2) - 
00209     TMath::Power(mMEAN + mAMP*mPE - mMEAN*mPE + 
00210       (sqrt(2/TMath::Pi())*mSIG*
00211          (1/(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG))) + 
00212            mPE/(-2 + TMath::Erfc(mMEAN/(sqrt(2)*mSIG)))
00213            ))/
00214        TMath::Power(TMath::E(),TMath::Power(mMEAN,2)/
00215          (2.*TMath::Power(mSIG,2))),2) + 
00216     (mMEAN*sqrt(2/TMath::Pi())*mSIG*
00217        (1/(1 + TMath::Erf(mMEAN/(sqrt(2)*mSIG))) + 
00218          mPE/(-2 + TMath::Erfc(mMEAN/(sqrt(2)*mSIG)))))
00219       /TMath::Power(TMath::E(),TMath::Power(mMEAN,2)/
00220        (2.*TMath::Power(mSIG,2))));
00221 }
00222            
00223 Double_t response_multia(Double_t* x, Double_t* params)
00224 {       
00225         double y = x[0] - mPDM;
00226         
00227                 
00228         double response=0;
00229         
00230         for(int i=3; i<=NPEAKS; i++)
00231         {
00232                 response += TMath::Poisson(i,mLAM) *TMath::Gaus(y,m_n(x,params)*i,sqrt(i*sigma_n(x,params)*sigma_n(x,params)+TMath::Power(mSHT,2)),true);
00233         }
00234         return mCON*response;
00235 }       
00236         
00237 Double_t response_multib(Double_t* x, Double_t* params)
00238 {       
00239         double y = x[0] - mPDM;
00240 
00241                 
00242         double response=0;
00243         
00244         for(int i=3; i<=NPEAKS; i++)
00245         {
00246                 response += TMath::Poisson(i,mLAM) *TMath::Gaus(y,m_n(x,params)*i,sqrt(i*sigma_n(x,params)*sigma_n(x,params)+TMath::Power(mSHT2,2)),true);
00247         }
00248         return mCON*response;
00249 }       
00250 
00251 Double_t response_multi(Double_t* x, Double_t* params) {
00252   return mPDWT*response_multia(x, params) + (1-mPDWT)*response_multib(x, params);
00253 }
00254 
00255 Double_t SPEFunc(Double_t* x, Double_t* params)
00256 {
00257         double sig = (response_0(x, params)+response_1(x,params)+response_2(x,params)+response_multi(x,params));
00258         /*double fixerparams[] = {mCON, 0.00143997, 83.9085, 29.0772, 149.998, 0.704112 , mSHT };
00259         double base= (response_1(x,fixerparams)+response_2(x,fixerparams)+response_multi(x,fixerparams));
00260         return sig*(TMath::Poisson(0,0.00143997))+base;*/
00261         return sig;
00262 
00263 }
00264 
00265 Double_t FitSPE(TH1* spe, ChanFitSettings& CFS, int ntriggers,
00266 bool allow_bg = true,
00267 bool force_old = false)
00268 {
00269 
00270 TF1* spefunc = (TF1*)gROOT->GetFunction("spefunc");
00271         int nEvtsInRange = spe->Integral(0,spe->GetNbinsX()+1,"width");
00272         
00273         spe->GetXaxis()->SetRangeUser(CFS.range_min, CFS.range_max);
00274         
00275         double fitmin, fitmax;
00276         fitmin = spe->GetBinLowEdge(spe->GetXaxis()->GetFirst());
00277         fitmax = 0;
00278         //find the last non-zero bin
00279         int bin = spe->GetXaxis()->GetLast();
00280         //bin = spe->GetNbinsX()+1;
00281         while(spe->GetBinContent(bin) == 0) {bin--;}
00282         fitmax = spe->GetBinLowEdge(bin);
00283         fitmax = fitmax*1; 
00284         cout<<"99999999999999999999                     "<<fitmax<<endl;
00285         //Pedestal centering
00286         
00287         Double_t pedmean = CFS.pedmean_min_bound;
00288         if(fitmin<0 && !force_old && (max(fitmin,CFS.pedrange_min)< min(CFS.pedrange_max,fitmax))){
00289                 spe->GetXaxis()->SetRangeUser(max(fitmin,CFS.pedrange_min), min(CFS.pedrange_max,fitmax));
00290                 TFitResultPtr pedfit = spe->Fit("gaus","MIS");
00291                 pedmean = pedfit->Value(1);
00292                 cout<<"!!!!!!!PED MEAN: " << pedmean<<endl;
00293         }
00294         if(force_old){
00295                         pedmean=spefunc->GetParameter(PEDMEAN);
00296                         cout<<"!!!!!! Force old. PEDMEAN: "<<pedmean<<endl;
00297                 }
00298         
00299         
00300         
00301         double params[NPAR];
00302         //find the likely place for the single p.e. peak
00303         bin = spe->FindBin(0);
00304         //params[SHOTNOISE] = FindExtremum(bin, spe, FORWARD, MINIMUM)/3.;
00305 
00306 
00307         
00308         if(!force_old || !spefunc){
00309                 params[SHOTNOISE] = CFS.shotnoise_start_value;
00310                 params[MEAN] = CFS.mean_start_value;
00311                 params[SIGMA] = CFS.sigma_start_value;
00312                 params[LAMBDA] = CFS.lambda_start_value;
00313                 if(CFS.constant_start_value){ params[CONSTANT] = nEvtsInRange;} else{  params[CONSTANT] = CFS.constant_start_value;}
00314                 params[AMP_E] = CFS.amp_E_start_value;
00315                 params[P_E] = CFS.p_E_start_value;  
00316                 
00317                 params[SHOTNOISE2] = CFS.shotnoise2_start_value;
00318                 params[PEDWEIGHT] = CFS.pedweight_start_value;
00319                 
00320                 spefunc = new TF1("spefunc", SPEFunc, fitmin, fitmax, NPAR);
00321 
00322                 spefunc->SetParameters(params);
00323                 for(int i=0; i<NPAR; i++)
00324                 spefunc->SetParName(i, names[i]);
00325                 //spefunc->SetParLimits(CONSTANT, 0 , params[CONSTANT]);
00326                 //DIVOT
00327                 spefunc->SetParLimits(LAMBDA, CFS.lambda_min_bound, CFS.lambda_max_bound);
00328                 spefunc->SetParLimits(MEAN, CFS.mean_min_bound, CFS.mean_max_bound);//spe_min, spe_max);
00329                 spefunc->SetParLimits(SHOTNOISE, CFS.shotnoise_min_bound, CFS.shotnoise_max_bound);//spe_min, spe_max);         
00330                 spefunc->SetParLimits(SIGMA, CFS.sigma_min_bound, CFS.sigma_max_bound);
00331                 spefunc->SetParLimits(P_E, CFS.p_E_min_bound, CFS.p_E_max_bound);
00332                 spefunc->SetParLimits(AMP_E, CFS.amp_E_min_bound, CFS.amp_E_max_bound);//Focus here
00333                 if(!allow_bg){
00334                         cout<<"Disallowing expoential background for this fit"<<endl;
00335                         spefunc->SetParameter(AMP_E,0);
00336                         spefunc->SetParameter(P_E,0);
00337                         spefunc->FixParameter(AMP_E,0);
00338                         spefunc->FixParameter(P_E,0);
00339                 }
00340                 
00341                 spefunc->SetParLimits(SHOTNOISE2, CFS.shotnoise2_min_bound, CFS.shotnoise2_max_bound);
00342                 spefunc->SetParLimits(PEDWEIGHT, CFS.pedweight_min_bound, CFS.pedweight_max_bound);
00343 
00344                 
00345                 spefunc->SetLineStyle(1);
00346                 spefunc->SetLineColor(kBlue);
00347         }
00348         else{
00349                 for(int i=0; i<NPAR; i++){
00350                         if(i == CONSTANT || i == LAMBDA)
00351                         continue;
00352                         double width=0;
00353                         spefunc->SetParLimits(i,
00354                                 std::max(0. ,spefunc->GetParameter(i)-width*spefunc->GetParError(i)),
00355                                 spefunc->GetParameter(i) + width*spefunc->GetParError(i));
00356                 }
00357         }
00358         spefunc->FixParameter(CONSTANT, nEvtsInRange);
00359         spefunc->FixParameter(PEDMEAN, pedmean);
00360         spe->GetXaxis()->SetRangeUser(fitmin, fitmax);
00361 
00362 
00364         spe->Draw();
00365         spefunc->Draw("same");
00366 
00368         //spe->Fit(spefunc,"MRES");
00369         //spe->Fit(spefunc,"MRL");
00370         TFitResultPtr fitResult = spe->Fit(spefunc,"MRS");
00371 
00372 
00373         //spefunc->DrawCopy("same");
00374         std::cout<<"Fit results: \n"<<
00375         "\t Chi2/NDF = "<<spefunc->GetChisquare()<<"/"<<spefunc->GetNDF()
00376         <<"\n\t Prob = "<<spefunc->GetProb()<<std::endl;
00377         for(int i=0; i<NPAR; i++)
00378                 params[i] = spefunc->GetParameter(i);
00379         TList* funclist = spe->GetListOfFunctions();
00380         static TF1* background = new TF1("background",background_func,fitmin,fitmax,NPAR);
00381         background->SetRange(fitmin, fitmax);
00382         background->SetLineColor(kRed);
00383         background->SetParameters(spefunc->GetParameters());
00384         funclist->Add(background->Clone());
00385         //TF1* bgcopy = background->DrawCopy("same");
00386         //funclist->Add(bgcopy);
00387 
00388         TF1* response_0_f = (TF1*)gROOT->GetFunction("response_0_f");   
00389         if(!response_0_f){response_0_f = new TF1("response_0_f",response_0,fitmin,fitmax,NPAR);}
00390         response_0_f->SetRange(fitmin, fitmax);
00391         response_0_f->SetLineColor(kGreen); 
00392         response_0_f->SetParameters(spefunc->GetParameters());
00393         response_0_f->SetRange(fitmin,70);
00394         if(response_0_f){funclist->Add(response_0_f->Clone());}
00395         
00396         TF1* response_1_f = (TF1*)gROOT->GetFunction("response_1_f");   
00397         if(!response_1_f){response_1_f = new TF1("response_1_f",response_1,fitmin,fitmax,NPAR);}
00398         response_1_f->SetRange(fitmin, fitmax);
00399         response_1_f->SetLineColor(kMagenta); 
00400         response_1_f->SetParameters(spefunc->GetParameters());
00401         if(response_1_f){funclist->Add(response_1_f->Clone());}
00402         
00403         TF1* response_2_f = (TF1*)gROOT->GetFunction("response_2_f");   
00404         if(!response_2_f){response_2_f = new TF1("response_2_f",response_2,fitmin,fitmax,NPAR);}
00405         response_2_f->SetRange(fitmin, fitmax);
00406         response_2_f->SetLineColor(kGreen); 
00407         response_2_f->SetParameters(spefunc->GetParameters());
00408         if(response_2_f){funclist->Add(response_2_f->Clone());}
00409         
00410         TF1* response_multi_f = (TF1*)gROOT->GetFunction("response_multi_f");   
00411         if(!response_multi_f){response_multi_f = new TF1("response_multi_f",response_multi,fitmin,fitmax,NPAR);}
00412         response_multi_f->SetRange(fitmin, fitmax);
00413         response_multi_f->SetLineColor(kMagenta); 
00414         response_multi_f->SetParameters(spefunc->GetParameters());
00415         if(response_multi_f){funclist->Add(response_multi_f->Clone());}
00416         
00417         TF1* response_0a_f = (TF1*)gROOT->GetFunction("response_0a_f");//new TF1("response_0a_f",response_0a,fitmin,fitmax,NPAR);
00418         if(!response_0a_f){response_0a_f = new TF1("response_0a_f",response_0a,fitmin,fitmax,NPAR);}
00419         response_0a_f->SetRange(fitmin, fitmax);
00420         response_0a_f->SetLineColor(kTeal);
00421         response_0a_f->SetParameters(spefunc->GetParameters());
00422         response_0a_f->SetRange(fitmin,70);
00423         if(response_0a_f){funclist->Add(response_0a_f->Clone());}
00424 
00425         TF1* response_0b_f = (TF1*)gROOT->GetFunction("response_0b_f"); //new TF1("response_0b_f",response_0b,fitmin,fitmax,NPAR);
00426         if(!response_0b_f){response_0b_f = new TF1("response_0b_f",response_0b,fitmin,fitmax,NPAR);}
00427         response_0b_f->SetRange(fitmin, fitmax);
00428         response_0b_f->SetLineColor(kTeal);
00429         response_0b_f->SetParameters(spefunc->GetParameters());
00430         response_0b_f->SetRange(fitmin,70);
00431         if(response_0b_f){funclist->Add(response_0b_f->Clone());}
00432         
00433 
00434         double thresh = spefunc->GetParameter(MEAN) - 2.*spefunc->GetParameter(SIGMA);
00435         //TLine* twosigma = new TLine(thresh, spe->GetYaxis()->GetBinLowEdge(1), thresh, 1.1*spe->GetMaximum());
00436         //twosigma->SetLineColor(kGreen);
00437         //twosigma->Draw();
00438 
00439         cout<<"Valid Events collected: "<<ntriggers<<std::endl;
00440         cout<<"Events Passing Cuts:" <<spe->GetEntries()<<std::endl;
00441         cout<<"Noise fraction: "
00442         <<spefunc->GetParameter(AMP_E)/spefunc->GetParameter(CONSTANT)
00443         <<std::endl;
00444         cout<<"Average photoelectrons per pulse: "<<spefunc->GetParameter(LAMBDA)
00445         <<std::endl;
00446         cout<<"Sigma width of pedestal: "<<spefunc->GetParameter(SHOTNOISE)<<endl;
00447         //double mean = spefunc->GetParameter(MEAN);
00448         
00449         double sigma_1=sqrt(mSHT*mSHT + mSIG*mSIG);     
00450         double m_n=(1-mPE)*(mMEAN+ sigma_1/(sqrt(2 * TMath::Pi()) * 0.5 *(1- TMath::Erf(-mMEAN/(sigma_1*sqrt(2))))) *TMath::Power(TMath::E(),-0.5 *mMEAN*mMEAN/(sigma_1 * sigma_1))) + mPE * mAMP;
00451         double sigma_n=sqrt((1-mPE)*(mMEAN*mMEAN+sigma_1*sigma_1+ sigma_1*mMEAN/(sqrt(2 * TMath::Pi()) * 0.5 *(1- TMath::Erf(-mMEAN/(sigma_1*sqrt(2))))) *TMath::Power(TMath::E(),-0.5 *mMEAN*mMEAN/(sigma_1 * sigma_1))) + 2* mPE * mAMP*mAMP-m_n*m_n);
00452         
00453         double conversion = 0.00049 /*V/ct*/ * 4E-9 /*ns/samp*/ * 1./25. /*ohm*/ *
00454         0.1 /*amplifier*/ * 1.E12 /*pC/C*/;
00455         
00456         cout<<"Single photoelectron peak: \n\t"<<m_n<<" +- " <<sigma_n <<" counts\n\t"
00457         <<m_n*conversion<<" pC\n\t"<< m_n*conversion / 1.6E-7<<" gain.\n";
00458         
00459         cout<<"Two-sigma threshold "<<thresh<<"counts.\n";
00460 
00461         //afan-----------
00462         TMatrixDSym cov = fitResult->GetCovarianceMatrix();
00463         double pdfmean_approx = mPE*mAMP+(1-mPE)*mMEAN;
00464         double pdfmean_error_uncorr = sqrt(TMath::Power((mAMP-mMEAN)*spefunc->GetParError(P_E),2)+TMath::Power((1-mPE)*spefunc->GetParError(MEAN),2)+TMath::Power(mPE*spefunc->GetParError(AMP_E),2));
00465         double pdfmean_error_corr = pdfmean_error_uncorr + 2*(mAMP-mMEAN)*(1-mPE)*cov[P_E][MEAN] + 2*(mAMP-mMEAN)*mPE*cov[P_E][AMP_E] + 2*(1-mPE)*mPE*cov[MEAN][AMP_E];
00466         cout<<"Approximated pdfmean (using corr): "<<pdfmean_approx<<" +- "<<pdfmean_error_corr<<endl;
00467         return spefunc->GetParameter(LAMBDA);
00468 }
00469 
00470 /* int ProcessSPEFile(const char* fname, Long_t roi = -1, int channel = -1, 
00471 double emax = -1, bool force_old = false)
00472 {
00473         if(!gPad) new TCanvas;
00474         gPad->SetLogy();
00475         gPad->SetTitle(fname);
00476         static bool loaded = false;
00477         if(!loaded){
00478                 gROOT->ProcessLine(".L lib/libDict.so");
00479                 loaded = true;
00480         }
00481 
00482         TFile* fin = new TFile(fname);
00483         if(!fin->IsOpen()){
00484                 std::cerr<<"Unable to open file "<<fname<<std::endl;
00485                 return 1;
00486         }
00487 
00488         TTree* Events = (TTree*)(fin->Get("Events"));
00489         if(!Events){
00490                 std::cerr<<"Unable to load Events tree from file "<<fname<<std::endl;
00491                 return 2;
00492         }
00493         TString data_source;
00494         if(roi == -1) data_source = "channels[].pulses.integral";
00495         else data_source = TString("-channels[].regions[") + roi
00496         + TString("].integral");
00497         if( emax < 0){
00498                 Events->Draw(data_source+" >> htemp",data_source+" > 0");
00499                 TH1* htemp = (TH1*)(gROOT->FindObject("htemp"));
00500                 emax = htemp->GetMean()*HISTOGRAMWIDTH;
00501         }
00502         
00503         TCut min_en = (data_source+" > 0").Data();
00504         char chstring[100];
00505         sprintf(chstring,data_source+" < %.0f",emax);
00506         TCut max_en = chstring;
00507         sprintf(chstring,"channels[].channel_id == %d",channel);
00508         TCut chan_cut = (channel == -1 ? "" : chstring);
00509 
00510         TCut time_cut = (roi == -1 ? get_time_cut(Events, chan_cut) : "" );
00511 
00512         TCut total_cut = min_en && max_en && time_cut && chan_cut;
00513         
00514         Events->Draw(data_source+" >> hspec",total_cut,"e");
00515         TH1* hspec = (TH1*)(gROOT->FindObject("hspec"));
00516 
00517 
00518         FitSPE(hspec, Events->GetEntries("channels[].baseline.found_baseline"),
00519         force_old);
00520 
00521 
00522         return 0;
00523 } */
00524