test/StatisticalErrors/MakePlot.C

0001 #include <iostream>
0002 #include <fstream>
0003 #include <sstream>
0004 #include "TH1F.h"
0005 #include "TCanvas.h"
0006 #include <vector>
0007 #include <string>
0008
0009 void makePlot(std::vector<double> vec, double& meanP, double& rmsP, double& sigmaP, std::string parNum);
0010 void skimValues(std::vector<double>& vec, double mean, double rms);
0011
0012 void MakePlot() {
0013   ifstream file("Values.txt");
0014   std::string valueString;
0015   std::vector<double> values;
0016   bool first=true;
0017   std::string numPar;
0018
0019   while( getline(file, valueString) ) {
0020     if(first) {
0021       numPar=valueString.substr(10,valueString.length());
0022       first=false;
0023     }
0024     else {
0025       std::stringstream ss(valueString);
0026       double value;
0027       ss >> value;
0028       values.push_back(value);
0029     }
0030   }
0031   if(file.is_open()) file.close();
0032
0033   double meanPlot=0, rmsPlot=0, sigmaPlot=0;
0034
0035   makePlot(values, meanPlot, rmsPlot, sigmaPlot, numPar);
0036
0037   // Uncomment if you don't want "adjustments"...
0038   //   std::cout << "sigma_final " << sigmaPlot << std::endl;
0039   //   return;
0040
0041   if( rmsPlot==0 || sigmaPlot==0 ) {
0042     std::cout << "sigma_final " << sigmaPlot << std::endl;
0043     return;
0044   }
0045
0046   if( sigmaPlot/rmsPlot<5. && rmsPlot/sigmaPlot<5. ) {
0047     std::cout << "sigma_final " << sigmaPlot << std::endl;
0048     return;
0049   }
0050
0051   std::cout << " Difference between RMS and sigma too large." << std::endl;
0052
0053   int prevDim=values.size();
0054   int iterN=1;
0055
0056   skimValues(values, meanPlot, rmsPlot);
0057   int lostHits=prevDim-values.size();
0058   int lostHitsTot=lostHits;
0059   prevDim=values.size();
0060
0061   while( (sigmaPlot/rmsPlot<5. || rmsPlot/sigmaPlot<5.) && lostHits!=0 ) {
0062     std::cout << " After iteration " << iterN << ", " << lostHits << " values rejected (" << lostHitsTot << " lost)." << std::endl;
0063     makePlot(values, meanPlot, rmsPlot, sigmaPlot, numPar);
0064     skimValues(values, meanPlot, rmsPlot);
0065     lostHits=prevDim-values.size();
0066     lostHitsTot+=lostHits;
0067     prevDim=values.size();
0068     iterN++;
0069   }
0070
0071   std::cout << "sigma_final " << sigmaPlot << std::endl;
0072   return;
0073 }
0074
0075 void makePlot(std::vector<double> vec, double& meanP, double& rmsP, double& sigmaP, std::string parNum) {
0076   int cnt=vec.size();
0077   double minV=99999., maxV=-99999.;
0078
0079   for(int kk=0; kk<cnt; ++kk) {
0080     if(vec[kk]<minV)
0081       minV=vec[kk];
0082     if(vec[kk]>maxV)
0083       maxV=vec[kk];
0084   }
0085
0086   double minH = minV - 0.1*(maxV-minV);
0087   double maxH = maxV + 0.1*(maxV-minV);
0088   TH1F *histo = new TH1F("value", ("Parameter "+parNum).c_str(), 100, minH, maxH);
0089   histo->GetXaxis()->SetTitle(("Value of parameter "+parNum).c_str());
0090   histo->GetYaxis()->SetTitle("Entries");
0091
0092   for(int i=0; i<cnt; ++i) {
0093     histo->Fill(vec[i]);
0094   }
0095   histo->Fit("gaus");
0096   TCanvas *cc = new TCanvas("cc","cc");
0097   histo->Draw();
0098   cc->SaveAs(("plot_param_"+parNum+".gif").c_str());
0099   meanP = histo->GetMean();
0100   rmsP = histo->GetRMS();
0101   sigmaP = ((TF1*)histo->GetListOfFunctions()->First())->GetParameter(2);
0102   if(histo) delete histo;
0103   if(cc) delete cc;
0104 }
0105
0106 void skimValues(std::vector<double>& vec, double mean, double rms) {
0107   int dim=vec.size();
0108   for(int jj=0; jj<dim;) {
0109     double it=vec[jj];
0110     if( fabs( (it-mean)/rms ) > 5. ) {
0111       vec.erase(vec.begin()+jj);
0112       dim--;
0113     }
0114     else {
0115       jj++;
0116     }
0117   }
0118 }