Subversion Repositories f9daq

#include <stdlib.h>
#include <stdio.h>
#include <TString.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TMultiGraph.h>
#include <math.h>

int IUdraw(const char *datoteka =  "C:/home/data/meritev-iv.dat"){
auto c1 = new TCanvas("c1","A Simple Graph Example",200,10,1400,1000);
double sigmaGauss = 0.3;
double xError = 0.5;
double rangeX = 5;
//c1->SetLogy();
//auto mg = new TMultiGraph();
 
// printf("#%s#\n", datoteka);
 /*
string s = datoteka;
s.erase(s.find_last_of("."), string::npos);
fprintf(stderr,"%s\n", s.c_str());
std::string key ("/");
string ime;
std::size_t found = s.rfind(key);
if (found!=std::string::npos) {
  ime = s.substr(found+1,s.length());
}
*/
TString ime = gSystem->BaseName(datoteka);
TString path = gSystem->DirName(datoteka);
ime.Remove(ime.Index(".dat"));
fprintf(stderr,"Ime=%s\n", ime.Data());
TString formatted = path + TString("/") + ime; // + TString(".png");
//formatted.Form("%s.png", );
fprintf(stderr,"png=%s\n", formatted.Data());
TString formattedGraph1;
TString formattedGraph2;
TString formattedGraph3;
TString formattedGraph4;
TString formattedGraph5;
TString formatted1;
TString formatted2;
TString formatted3;
TString formatted4;
TString formatted5;

for (Int_t d = 0; d<1; d++) {
  FILE * fp=fopen(datoteka,"r");
  if (!fp) continue;
  float f[4];
  Int_t j=0;  
  Int_t ndim=400;
  char line[400];  
  double fpXX[250];
  double fpYY[250];
//  Int_t k = 0;
  while (fgets(line,ndim,fp)!=NULL) {
        if (line[0] == '#') continue;
    sscanf(line,"%f%f%f%f",&f[0],&f[1],&f[2],&f[3]);
//      printf("%f %f %f %f \n",f[0],f[1],f[2],f[3]);
        fpXX[j]=f[2];
        fpYY[j]=f[3];
        j++;
  }
  // float *fpx =    new float[j+1];
  // float *fpy =    new float[j+1];
  // for (Int_t i=0; i<j-3; i++) {
    // if (fpXX[i]!= 0) {
                // k++;
                // fpx[i]=fpXX[i+2];
                // fpy[i]=fpYY[i+2];
                // printf("%2.6f\t",fpx[i]);
                // printf("\n");
                // printf("%2.6f\t",fpy[i]);
                // printf("\n");
        // }
  // }
  fclose(fp);
 
  auto graf11 = new TGraph();
  graf11->SetMarkerStyle(20);
  graf11->SetDrawOption("AP");
  graf11->SetLineWidth(3);
  graf11->SetFillStyle(0);
 
  auto graf12 = new TGraph();
  graf12->SetMarkerStyle(20);
  graf12->SetDrawOption("AP");
  graf12->SetLineWidth(3);
  graf12->SetFillStyle(0);
 
  auto graf2 = new TGraph();
  graf2->SetMarkerStyle(20);
  graf2->SetDrawOption("AP");
  graf2->SetLineWidth(3);
  graf2->SetFillStyle(0);
 
  auto graf3 = new TGraph();
  graf3->SetMarkerStyle(20);
  graf3->SetDrawOption("AP");
  graf3->SetLineWidth(3);
  graf3->SetFillStyle(0);
 
  auto graf4 = new TGraph();
  graf4->SetMarkerStyle(20);
  graf4->SetDrawOption("AP");
  graf4->SetLineWidth(3);
  graf4->SetFillStyle(0);
 
  auto graf5 = new TGraph();
  graf5->SetMarkerStyle(20);
  graf5->SetDrawOption("AP");
  graf5->SetLineWidth(3);
  graf5->SetFillStyle(0);
 
  for (int i=0; i<j; i++) {
        if (fpYY[i]<=0) fpYY[i] = 1e-12;  
        graf11->SetPoint(i,fpXX[i],log(fpYY[i]));
        graf12->SetPoint(i,fpXX[i],log(fpYY[i]));
        graf5->SetPoint(i,log(fpYY[i]),fpXX[i]);
  }
 
  double max2 = 0;
  double max2x = 0;
  int k = 0;
  for (int i=0; i<j-1; i++) {
        if (fpYY[i]<=0) fpYY[i] = 1e-12;
        double x2 = (fpXX[i+1]+fpXX[i])/2;
        double y2 = (log(fpYY[i+1])-log(fpYY[i]))/(fpXX[i+1]-fpXX[i]);
        graf2->SetPoint(i,x2, y2);
        if (y2 > 1.1) {
                double y3 = 1 / y2;
                graf3->SetPoint(k,x2, y3);
                k++;
        }
        if (y2 > max2) {
                max2 = y2;
                max2x = x2;
        }
  }
 
  double max4 = 0;
  double max4x = 0;
  double x4stari = 0;
  double y4stari = 0;
  double x4novi = 0;
  double y4novi = 0;
  for (int i=0; i<j-2; i++) {
        if (fpYY[i]<=0) fpYY[i] = 1e-12;
        double x4 = (fpXX[i+2]+fpXX[i])/4+fpXX[i+1]/2;
        double y4 = ((log(fpYY[i+2])-log(fpYY[i+1]))/(fpXX[i+2]-fpXX[i+1]) - (log(fpYY[i+1])-log(fpYY[i]))/(fpXX[i+1]-fpXX[i]))/((fpXX[i+2]-fpXX[i])/2);
        if (i==0) {
                x4novi = x4;
                y4novi = y4;
        }
        else {
                x4novi = (x4 + x4stari)/2;
                y4novi = (y4 + y4stari)/2;
        }
        graf4->SetPoint(i,x4novi, y4novi);
        x4stari = x4;
        y4stari = y4;
        if (y4 > max4) {
                max4 = y4;
                max4x = x4;
        }
  }
}

TF1 *fa1 = new TF1("fa1","pol1",61,67);
fa1->SetParameters(-28.4467,0.0875386);
graf11->Fit("fa1", "R");
graf11->Draw("APL");
c1->Update();
TF1 *fa2 = new TF1("fa2","pol1",69.7,70.1);
graf11->Fit("fa2", "R+");
//graf12->Draw("APL");
printf("p0=%e, p1=%e\n", fa1->GetParameter(0), fa1->GetParameter(1));
double par0 = fa1->GetParameter(0);
double par1 = fa1->GetParameter(1);
printf("p2=%e, p3=%e\n", fa2->GetParameter(0), fa2->GetParameter(1));
double par2 = fa2->GetParameter(0);
double par3 = fa2->GetParameter(1);
double x1 = (par0 - par2)/(par3 - par1);
formattedGraph1.Form("I ( V ) -> %s, V(1) = %f ; V [V]; ln(I) [A]", ime.Data(), x1);
graf11->SetTitle(formattedGraph1);
graf12->SetTitle(formattedGraph1);
//graf13->Draw("APL");
formatted1.Form("%s-1.png", formatted.Data());
c1->SaveAs(formatted1);

auto c2 = new TCanvas("c2","A Simple Graph Example",200,10,1400,1000);
graf2->Draw("APL");
graf2->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
TF1 *fa2 = new TF1("fa2","gaus(0)+pol1(3)",max2x - sigmaGauss, max2x + sigmaGauss);
 fa2->SetParameters(60,x1,0.1339,0,0);
 fa2->SetParLimits(1, max2x - sigmaGauss, max2x + sigmaGauss);
// fa1->SetParLimits(1, 10**(-14), 10**(-10));
 graf2->Fit("fa2", "R");
// printf("p20=%e, p21=%e, p22=%e, max2x=%f\n", fa2->GetParameter(0), fa2->GetParameter(1), fa2->GetParameter(2), max2x);
// double par2 = fa2->GetParameter(2);
// double par3 = fa2->GetParameter(3);
double x2 = fa2->GetParameter(1);
formattedGraph2.Form("I ( V ) -> %s, V(2) = %f ; Voltage [V]; d(ln(I)/dV)", ime.Data(), x2);
graf2->SetTitle(formattedGraph2);
formatted2.Form("%s-2.png", formatted.Data());
gStyle->SetOptFit(1);
c2->SaveAs(formatted2);

auto c3 = new TCanvas("c3","A Simple Graph Example",200,10,1400,1000);
graf3->Draw("APL");
graf3->Fit("pol3");
graf3->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
TF1 *fa3 = new TF1("fa3", "pol1", x2, x2 + 0.5);
fa3->SetParameters(-30.6306,0.4398);
//fa3->SetParLimits(0, 10**(-14), 10**(-10));
//fa3->SetParLimits(1, 10**(-14), 10**(-10));
graf3->Fit("fa3", "R");
double n3 = fa3->GetParameter(0);
double k3 = fa3->GetParameter(1);
double x3 = - n3 / k3;
formattedGraph3.Form("I ( V ) -> %s, V(3) = %f ; Voltage [V]; 1/(d(ln(I)/dV))", ime.Data(), x3);
graf3->SetTitle(formattedGraph3);
//printf("p30=%e, p31=%e\n", fa3->GetParameter(0), fa3->GetParameter(1));
formatted3.Form("%s-3.png", formatted.Data());
c3->SaveAs(formatted3);

auto c4 = new TCanvas("c4","A Simple Graph Example",200,10,1400,1000);
graf4->Draw("APL");
graf4->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
TF1 *fa4 = new TF1("fa4","gaus(0)",max4x - sigmaGauss, max4x + sigmaGauss);
 fa4->SetParameters(60,x2,0.1339);
 fa4->SetParLimits(1, x2 - xError, x2 + xError);
 fa4->SetParLimits(0, 0, 1000);
 graf4->Fit("fa4", "R");
// printf("p40=%e, p41=%e, p42=%e, max4x=%f\n", fa4->GetParameter(0), fa4->GetParameter(1), fa4->GetParameter(2), max4x);
// double par2 = fa2->GetParameter(2);
// double par3 = fa2->GetParameter(3);
double x4 = fa4->GetParameter(1);
formatted4.Form("%s-4.png", formatted.Data());
// TString stringx4;
// stringx4.Form("V(4) = %f", x4);
   // TText *t = new TText(60,70,stringx4);
   // t->SetTextAlign(22);
   // t->SetTextColor(1);
   // t->SetTextFont(43);
   // t->SetTextSize(20);
   // t->Draw();
formattedGraph4.Form("I ( V ) -> %s, V(4) = %f ; Voltage [V]; d^2(ln(I)/dV^2)", ime.Data(), x4);
graf4->SetTitle(formattedGraph4);
gPad->Modified(); gPad->Update();
c4->SaveAs(formatted4);

auto c5 = new TCanvas("c5","A Simple Graph Example",200,10,1400,1000);
graf5->Draw("APL");
TF1 *fa5 = new TF1("fa5","-[0] - [1]*x - [2]*x**2",par0 + par1*x1 +4,par0 + par1*x1+11);
// fa5->SetParameters(28.4467,0.0875386,4);
graf5->Fit("fa5", "R+");
printf("p4=%e, p5=%e, p6=%e\n", fa5->GetParameter(0), fa5->GetParameter(1), fa5->GetParameter(2));
double par4 = fa5->GetParameter(0);
double par5 = fa5->GetParameter(1);
double par6 = fa5->GetParameter(2);
double xx5 = - (double)(par5/(2*par6));
double x5 = - par4 - par5 * xx5 - par6 * xx5**2;
// double f5(double p);
// double f5(double p) {
        // double par0, par1, par4, par5, par6;
    // double a=-(double)(par0/par1) + (double)(1/par1) * p + par4 + par5 * p + par6 * p**2;
    // return a;
// }
// cout.precision(4);        
// cout.setf(ios::fixed);
// double a5,b5,cc5,e5,faa5,fb5,fc5;
// a5 = x1-5;
// b5 = x1+5;
// e5 = 0.001;
    // if (f5(a5)*f5(b5)>0) print("No root exists between a and b.\n");
    // else {
    // while (fabs(a5-b5)>=e5) {
        // cc5=(a5+b5)/2.0;        //bisect the interval and find the value of c
        // faa5=f5(a5);
        // fb5=f5(b5);
        // fc5=f5(cc5);  
        // if (fc5==0) {
            // printf("The root of the equation is %f \n",cc5);
            // break;
        // }
        // if (fa5*fc5>0) a5=cc5;
        // else if (fa5*fc5<0) b5=cc5;
    // }
    // }            //The loop ends when the difference between a and b becomes less than the desired accuracy ie now the value stored in 'c' can be called the approximate root of the equation        
    // printf("The root of the equation is %f \n",cc5);  
formattedGraph5.Form("V ( I ) -> %s, V(5) = %f ; ln(I) [A] ; V [V]", ime.Data(), x5);
graf5->SetTitle(formattedGraph5);
formatted5.Form("%s-5.png", formatted.Data());
c5->SaveAs(formatted5);

printf("RESULT: V(1) = %f, V(2) = %f, V(3) = %f, V(4) = %f, V(5) = %f", x1, x2, x3, x4, x5);
return 0;
}