Subversion Repositories f9daq

Rev

Rev 281 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <TString.h>
  4. #include <TGraph.h>
  5. #include <TCanvas.h>
  6. #include <TMultiGraph.h>
  7. #include <math.h>
  8.  
  9. int IUdraw(const char *datoteka =  "C:/home/data/meritev-iv.dat"){
  10. auto c1 = new TCanvas("c1","A Simple Graph Example",200,10,1400,1000);
  11. double sigmaGauss = 0.3;
  12. double xError = 0.5;
  13. double rangeX = 5;
  14. //c1->SetLogy();
  15. //auto mg = new TMultiGraph();
  16.  
  17. // printf("#%s#\n", datoteka);
  18.  /*
  19. string s = datoteka;
  20. s.erase(s.find_last_of("."), string::npos);
  21. fprintf(stderr,"%s\n", s.c_str());
  22. std::string key ("/");
  23. string ime;
  24. std::size_t found = s.rfind(key);
  25. if (found!=std::string::npos) {
  26.   ime = s.substr(found+1,s.length());
  27. }
  28. */
  29. TString ime = gSystem->BaseName(datoteka);
  30. TString path = gSystem->DirName(datoteka);
  31. ime.Remove(ime.Index(".dat"));
  32. fprintf(stderr,"Ime=%s\n", ime.Data());
  33. TString formatted = path + TString("/") + ime; // + TString(".png");
  34. //formatted.Form("%s.png", );
  35. fprintf(stderr,"png=%s\n", formatted.Data());
  36. TString formattedGraph1;
  37. TString formattedGraph2;
  38. TString formattedGraph3;
  39. TString formattedGraph4;
  40. TString formattedGraph5;
  41. TString formatted1;
  42. TString formatted2;
  43. TString formatted3;
  44. TString formatted4;
  45. TString formatted5;
  46.  
  47. for (Int_t d = 0; d<1; d++) {
  48.   FILE * fp=fopen(datoteka,"r");
  49.   if (!fp) continue;
  50.   float f[4];
  51.   Int_t j=0;  
  52.   Int_t ndim=400;
  53.   char line[400];  
  54.   double fpXX[250];
  55.   double fpYY[250];
  56. //  Int_t k = 0;
  57.   while (fgets(line,ndim,fp)!=NULL) {
  58.         if (line[0] == '#') continue;
  59.     sscanf(line,"%f%f%f%f",&f[0],&f[1],&f[2],&f[3]);
  60. //      printf("%f %f %f %f \n",f[0],f[1],f[2],f[3]);
  61.         fpXX[j]=f[2];
  62.         fpYY[j]=f[3];
  63.         j++;
  64.   }
  65.   // float *fpx =    new float[j+1];
  66.   // float *fpy =    new float[j+1];
  67.   // for (Int_t i=0; i<j-3; i++) {
  68.     // if (fpXX[i]!= 0) {
  69.                 // k++;
  70.                 // fpx[i]=fpXX[i+2];
  71.                 // fpy[i]=fpYY[i+2];
  72.                 // printf("%2.6f\t",fpx[i]);
  73.                 // printf("\n");
  74.                 // printf("%2.6f\t",fpy[i]);
  75.                 // printf("\n");
  76.         // }
  77.   // }
  78.   fclose(fp);
  79.  
  80.   auto graf11 = new TGraph();
  81.   graf11->SetMarkerStyle(20);
  82.   graf11->SetDrawOption("AP");
  83.   graf11->SetLineWidth(3);
  84.   graf11->SetFillStyle(0);
  85.  
  86.   auto graf12 = new TGraph();
  87.   graf12->SetMarkerStyle(20);
  88.   graf12->SetDrawOption("AP");
  89.   graf12->SetLineWidth(3);
  90.   graf12->SetFillStyle(0);
  91.  
  92.   auto graf2 = new TGraph();
  93.   graf2->SetMarkerStyle(20);
  94.   graf2->SetDrawOption("AP");
  95.   graf2->SetLineWidth(3);
  96.   graf2->SetFillStyle(0);
  97.  
  98.   auto graf3 = new TGraph();
  99.   graf3->SetMarkerStyle(20);
  100.   graf3->SetDrawOption("AP");
  101.   graf3->SetLineWidth(3);
  102.   graf3->SetFillStyle(0);
  103.  
  104.   auto graf4 = new TGraph();
  105.   graf4->SetMarkerStyle(20);
  106.   graf4->SetDrawOption("AP");
  107.   graf4->SetLineWidth(3);
  108.   graf4->SetFillStyle(0);
  109.  
  110.   auto graf5 = new TGraph();
  111.   graf5->SetMarkerStyle(20);
  112.   graf5->SetDrawOption("AP");
  113.   graf5->SetLineWidth(3);
  114.   graf5->SetFillStyle(0);
  115.  
  116.   for (int i=0; i<j; i++) {
  117.         if (fpYY[i]<=0) fpYY[i] = 1e-12;  
  118.         graf11->SetPoint(i,fpXX[i],log(fpYY[i]));
  119.         graf12->SetPoint(i,fpXX[i],log(fpYY[i]));
  120.         graf5->SetPoint(i,log(fpYY[i]),fpXX[i]);
  121.   }
  122.  
  123.   double max2 = 0;
  124.   double max2x = 0;
  125.   int k = 0;
  126.   for (int i=0; i<j-1; i++) {
  127.         if (fpYY[i]<=0) fpYY[i] = 1e-12;
  128.         double x2 = (fpXX[i+1]+fpXX[i])/2;
  129.         double y2 = (log(fpYY[i+1])-log(fpYY[i]))/(fpXX[i+1]-fpXX[i]);
  130.         graf2->SetPoint(i,x2, y2);
  131.         if (y2 > 1.1) {
  132.                 double y3 = 1 / y2;
  133.                 graf3->SetPoint(k,x2, y3);
  134.                 k++;
  135.         }
  136.         if (y2 > max2) {
  137.                 max2 = y2;
  138.                 max2x = x2;
  139.         }
  140.   }
  141.  
  142.   double max4 = 0;
  143.   double max4x = 0;
  144.   double x4stari = 0;
  145.   double y4stari = 0;
  146.   double x4novi = 0;
  147.   double y4novi = 0;
  148.   for (int i=0; i<j-2; i++) {
  149.         if (fpYY[i]<=0) fpYY[i] = 1e-12;
  150.         double x4 = (fpXX[i+2]+fpXX[i])/4+fpXX[i+1]/2;
  151.         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);
  152.         if (i==0) {
  153.                 x4novi = x4;
  154.                 y4novi = y4;
  155.         }
  156.         else {
  157.                 x4novi = (x4 + x4stari)/2;
  158.                 y4novi = (y4 + y4stari)/2;
  159.         }
  160.         graf4->SetPoint(i,x4novi, y4novi);
  161.         x4stari = x4;
  162.         y4stari = y4;
  163.         if (y4 > max4) {
  164.                 max4 = y4;
  165.                 max4x = x4;
  166.         }
  167.   }
  168. }
  169.  
  170. TF1 *fa1 = new TF1("fa1","pol1",61,67);
  171. fa1->SetParameters(-28.4467,0.0875386);
  172. graf11->Fit("fa1", "R");
  173. graf11->Draw("APL");
  174. c1->Update();
  175. TF1 *fa2 = new TF1("fa2","pol1",69.7,70.1);
  176. graf11->Fit("fa2", "R+");
  177. //graf12->Draw("APL");
  178. printf("p0=%e, p1=%e\n", fa1->GetParameter(0), fa1->GetParameter(1));
  179. double par0 = fa1->GetParameter(0);
  180. double par1 = fa1->GetParameter(1);
  181. printf("p2=%e, p3=%e\n", fa2->GetParameter(0), fa2->GetParameter(1));
  182. double par2 = fa2->GetParameter(0);
  183. double par3 = fa2->GetParameter(1);
  184. double x1 = (par0 - par2)/(par3 - par1);
  185. formattedGraph1.Form("I ( V ) -> %s, V(1) = %f ; V [V]; ln(I) [A]", ime.Data(), x1);
  186. graf11->SetTitle(formattedGraph1);
  187. graf12->SetTitle(formattedGraph1);
  188. //graf13->Draw("APL");
  189. formatted1.Form("%s-1.png", formatted.Data());
  190. c1->SaveAs(formatted1);
  191.  
  192. auto c2 = new TCanvas("c2","A Simple Graph Example",200,10,1400,1000);
  193. graf2->Draw("APL");
  194. graf2->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
  195. TF1 *fa2 = new TF1("fa2","gaus(0)+pol1(3)",max2x - sigmaGauss, max2x + sigmaGauss);
  196.  fa2->SetParameters(60,x1,0.1339,0,0);
  197.  fa2->SetParLimits(1, max2x - sigmaGauss, max2x + sigmaGauss);
  198. // fa1->SetParLimits(1, 10**(-14), 10**(-10));
  199.  graf2->Fit("fa2", "R");
  200. // printf("p20=%e, p21=%e, p22=%e, max2x=%f\n", fa2->GetParameter(0), fa2->GetParameter(1), fa2->GetParameter(2), max2x);
  201. // double par2 = fa2->GetParameter(2);
  202. // double par3 = fa2->GetParameter(3);
  203. double x2 = fa2->GetParameter(1);
  204. formattedGraph2.Form("I ( V ) -> %s, V(2) = %f ; Voltage [V]; d(ln(I)/dV)", ime.Data(), x2);
  205. graf2->SetTitle(formattedGraph2);
  206. formatted2.Form("%s-2.png", formatted.Data());
  207. gStyle->SetOptFit(1);
  208. c2->SaveAs(formatted2);
  209.  
  210. auto c3 = new TCanvas("c3","A Simple Graph Example",200,10,1400,1000);
  211. graf3->Draw("APL");
  212. graf3->Fit("pol3");
  213. graf3->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
  214. TF1 *fa3 = new TF1("fa3", "pol1", x2, x2 + 0.5);
  215. fa3->SetParameters(-30.6306,0.4398);
  216. //fa3->SetParLimits(0, 10**(-14), 10**(-10));
  217. //fa3->SetParLimits(1, 10**(-14), 10**(-10));
  218. graf3->Fit("fa3", "R");
  219. double n3 = fa3->GetParameter(0);
  220. double k3 = fa3->GetParameter(1);
  221. double x3 = - n3 / k3;
  222. formattedGraph3.Form("I ( V ) -> %s, V(3) = %f ; Voltage [V]; 1/(d(ln(I)/dV))", ime.Data(), x3);
  223. graf3->SetTitle(formattedGraph3);
  224. //printf("p30=%e, p31=%e\n", fa3->GetParameter(0), fa3->GetParameter(1));
  225. formatted3.Form("%s-3.png", formatted.Data());
  226. c3->SaveAs(formatted3);
  227.  
  228. auto c4 = new TCanvas("c4","A Simple Graph Example",200,10,1400,1000);
  229. graf4->Draw("APL");
  230. graf4->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);
  231. TF1 *fa4 = new TF1("fa4","gaus(0)",max4x - sigmaGauss, max4x + sigmaGauss);
  232.  fa4->SetParameters(60,x2,0.1339);
  233.  fa4->SetParLimits(1, x2 - xError, x2 + xError);
  234.  fa4->SetParLimits(0, 0, 1000);
  235.  graf4->Fit("fa4", "R");
  236. // printf("p40=%e, p41=%e, p42=%e, max4x=%f\n", fa4->GetParameter(0), fa4->GetParameter(1), fa4->GetParameter(2), max4x);
  237. // double par2 = fa2->GetParameter(2);
  238. // double par3 = fa2->GetParameter(3);
  239. double x4 = fa4->GetParameter(1);
  240. formatted4.Form("%s-4.png", formatted.Data());
  241. // TString stringx4;
  242. // stringx4.Form("V(4) = %f", x4);
  243.    // TText *t = new TText(60,70,stringx4);
  244.    // t->SetTextAlign(22);
  245.    // t->SetTextColor(1);
  246.    // t->SetTextFont(43);
  247.    // t->SetTextSize(20);
  248.    // t->Draw();
  249. formattedGraph4.Form("I ( V ) -> %s, V(4) = %f ; Voltage [V]; d^2(ln(I)/dV^2)", ime.Data(), x4);
  250. graf4->SetTitle(formattedGraph4);
  251. gPad->Modified(); gPad->Update();
  252. c4->SaveAs(formatted4);
  253.  
  254. auto c5 = new TCanvas("c5","A Simple Graph Example",200,10,1400,1000);
  255. graf5->Draw("APL");
  256. TF1 *fa5 = new TF1("fa5","-[0] - [1]*x - [2]*x**2",par0 + par1*x1 +4,par0 + par1*x1+11);
  257. // fa5->SetParameters(28.4467,0.0875386,4);
  258. graf5->Fit("fa5", "R+");
  259. printf("p4=%e, p5=%e, p6=%e\n", fa5->GetParameter(0), fa5->GetParameter(1), fa5->GetParameter(2));
  260. double par4 = fa5->GetParameter(0);
  261. double par5 = fa5->GetParameter(1);
  262. double par6 = fa5->GetParameter(2);
  263. double xx5 = - (double)(par5/(2*par6));
  264. double x5 = - par4 - par5 * xx5 - par6 * xx5**2;
  265. // double f5(double p);
  266. // double f5(double p) {
  267.         // double par0, par1, par4, par5, par6;
  268.     // double a=-(double)(par0/par1) + (double)(1/par1) * p + par4 + par5 * p + par6 * p**2;
  269.     // return a;
  270. // }
  271. // cout.precision(4);        
  272. // cout.setf(ios::fixed);
  273. // double a5,b5,cc5,e5,faa5,fb5,fc5;
  274. // a5 = x1-5;
  275. // b5 = x1+5;
  276. // e5 = 0.001;
  277.     // if (f5(a5)*f5(b5)>0) print("No root exists between a and b.\n");
  278.     // else {
  279.     // while (fabs(a5-b5)>=e5) {
  280.         // cc5=(a5+b5)/2.0;        //bisect the interval and find the value of c
  281.         // faa5=f5(a5);
  282.         // fb5=f5(b5);
  283.         // fc5=f5(cc5);  
  284.         // if (fc5==0) {
  285.             // printf("The root of the equation is %f \n",cc5);
  286.             // break;
  287.         // }
  288.         // if (fa5*fc5>0) a5=cc5;
  289.         // else if (fa5*fc5<0) b5=cc5;
  290.     // }
  291.     // }            //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        
  292.     // printf("The root of the equation is %f \n",cc5);  
  293. formattedGraph5.Form("V ( I ) -> %s, V(5) = %f ; ln(I) [A] ; V [V]", ime.Data(), x5);
  294. graf5->SetTitle(formattedGraph5);
  295. formatted5.Form("%s-5.png", formatted.Data());
  296. c5->SaveAs(formatted5);
  297.  
  298. printf("RESULT: V(1) = %f, V(2) = %f, V(3) = %f, V(4) = %f, V(5) = %f", x1, x2, x3, x4, x5);
  299. return 0;
  300. }
  301.