{
#include <string.h>
#include <algorithm> // std::min
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <Riostream.h>
#include <stdio.h>
//#include <stdlib.h>
//#include <iostream.h>
#include <sstream>
#include <string>
#include <cmath>
#include "TROOT.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TString.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TLeaf.h"
#include "TApplication.h" //"TGClient.h"
#include "TGClient.h"
gROOT->Reset();
ofstream out, out2;
// Theis file records ii jj a1 a2 w[ii][jj] P1 P2 P3 P4 S1 S2 S3 S4 B[ii][jj] in columns
out2.open("output2.dat");
//Theis file records ii jj (Artificial spectrum) Artificial->GetBinContent(ii,jj) B[ii][jj] (background area)
out.open("output.dat");
int nrange;
int m,ii,jj;
int i, seed;
seed = time(NULL);
srand(seed);
// to control the number of bins and iteration
cout<<"Number of bins"<< endl;
cin >> nrange;
cout<<"Iterations"<< endl;
cin >> m;
float v[nrange][nrange];
double w[nrange][nrange];
double G[nrange][nrange], G1[nrange][nrange], G2[nrange][nrange], G3[nrange][nrange];
double B[nrange][nrange];
int C[nrange][nrange];
double rnd,P1,P2,P3,P4,S1,S2,S3,S4,L1,L2,L3,L4,M1,M2,M3,M4,a1,a2,K1,K2,K3,K4,Q1,Q2,Q3,Q4,D1,D2,D3,D4;
Int_t nlines = 0;
Int_t nbins=nrange, nbins2d=nrange;
int minE2d=0, maxE2d=nrange-1, coinmin=0, coinmax=nrange-1;
TRandom *r3 = new TRandom3();
TFile *f = new TFile("cuts.root","RECREATE");
// Artificial Gaussian Spectrum
TH2F *Artificial = new TH2F("Artif","",nbins2d,minE2d,maxE2d,nbins2d,minE2d,maxE2d);
// Spectrum after subtracting background
TH2F *Bg_Sub = new TH2F("Bg_Sub","",nbins2d,minE2d,maxE2d,nbins2d,minE2d,maxE2d);
// background
TH2F *Bg = new TH2F("Bg","",nbins2d,minE2d,maxE2d,nbins2d,minE2d,maxE2d);
// Artificial Gaussian Spectrum after the background subtraction
TH2F *Artificial1 = new TH2F("Artif1","",nbins2d,minE2d,maxE2d,nbins2d,minE2d,maxE2d);
// Genertaing the artficial spectrum
for(ii=0; ii<nrange; ii++)
{
for(jj=0; jj<nrange; jj++)
{
// Genertaing the artficial spectrum
//-------------------------------------RANDOM----------------------------------------------------------------------------
// Ridges
G[ii][jj]=((r3->Rndm())*5+50)*exp(-(pow(((double)ii-64),2))/20)+((r3->Rndm())*5+50)*exp(-(pow(((double)jj-128),2))/20);
// peak
out<<G1[ii][jj]<<endl;
G1[ii][jj]=300*exp(-(((pow(((double)ii-64),2))/20)+(pow(((double)jj-128),2))/20));
// Random backdround
G3[ii][jj]=(r3->Rndm())*5;
// Total spectrum before subtraction
G2[ii][jj]=G[ii][jj]+G1[ii][jj]+G3[ii][jj];
Artificial->SetBinContent(ii,jj,G2[ii][jj]);
Artificial1->SetBinContent(ii,jj,G1[ii][jj]);
if (ii==64&&jj==128) cout << G2[ii][jj] << endl;
}
}
//--------------------------------------------------------case3 (LLS)----------------------------------------------------
for (int ii=1;ii<coinmax;ii++)
{
for (int jj=1;jj<coinmax;jj++)
{
v[ii][jj] = log(log(sqrt(Artificial->GetBinContent(ii,jj)+1)+1)+1);
}
}
for(p=1; p<=m;p++)
{
cout << "Iteration number " << p << endl;
for(ii=p;ii<coinmax-p;ii++)
{
for( jj=p;jj<coinmax-p;jj++)
{
Q1=v[ii+p][jj+p];
Q2=v[ii-p][jj+p];
Q3=v[ii+p][jj-p];
Q4=v[ii-p][jj-p];
D1=v[ii+p][jj];
D2=v[ii][jj+p];
D3=v[ii][jj-p];
D4=v[ii-p][jj];
S1=((exp(exp(D1)-1)-1) * (exp(exp(D1)-1)-1))-1;
S2=((exp(exp(D2)-1)-1) * (exp(exp(D2)-1)-1))-1;
S3=((exp(exp(D3)-1)-1) * (exp(exp(D3)-1)-1))-1;
S4=((exp(exp(D4)-1)-1) * (exp(exp(D4)-1)-1))-1;
P1=((exp(exp(Q1)-1)-1) * (exp(exp(Q1)-1)-1))-1;
P2=((exp(exp(Q2)-1)-1) * (exp(exp(Q2)-1)-1))-1;
P3=((exp(exp(Q3)-1)-1) * (exp(exp(Q3)-1)-1))-1;
P4=((exp(exp(Q4)-1)-1) * (exp(exp(Q4)-1)-1))-1;
L1=(max(S1,(P1+P3)/2)); // -((P1+P3)/2);
L2=(max(S2,(P1+P2)/2)); // -((P1+P2)/2);
L3=(max(S3,(P3+P4)/2)); // -((P3+P4)/2);
L4=(max(S4,(P2+P4)/2)); // -((P2+P4)/2);
M1=L1-((P1+P3)/2);
M2=L2-((P1+P2)/2);
M3=L3-((P3+P4)/2);
M4=L4-((P2+P4)/2);
a1=((exp(exp(v[ii][jj])-1)-1) * (exp(exp(v[ii][jj])-1)-1))-1;
a2=((M1+M4)/2)+((M2+M3)/2)+((P1+P2+P3+P4)/4);
w[ii][jj]=min(a1,a2);
}
}
for(ii=p;ii<coinmax-p;ii++)
{
for( jj=p;jj<coinmax-p;jj++)
{
v[ii][jj]=w[ii][jj];
}
}
}
for (int ii=1;ii<coinmax;ii++)
{
for (int jj=1;jj<coinmax;jj++)
{
B[ii][jj]= w[ii][jj];
Bg_Sub->SetBinContent(ii,jj,((Artificial->GetBinContent(ii,jj))-B[ii][jj]));
}
}
//------------------------------------------------------------------------------------------------------------------------------
*/
TCanvas *c1 = new TCanvas("ccc","ccc",800,1000);
Artificial->Draw("lego");
gStyle->SetPalette(1);
gStyle->SetOptStat(0);
//Add12->SetMaximum(400.0);
//Add12->SetMinimum(0.0);
TCanvas *c2 = new TCanvas("ccc2","ccc2",800,1000);
Bg_Sub->Draw("lego");
//Bg_Sub->SetMaximum(400.0);
//Bg_Sub->SetMinimum(0.0);
//c1->cd(3);
//Bg->Draw("lego");
TCanvas *c3 = new TCanvas("ccc3","ccc3",800,1000);
Artificial1->Draw("lego");
TCutG *cut21 = new TCutG("cut21",5);
cut21->SetPoint(0,60.75,124.5);
cut21->SetPoint(1,60.75,129.5);
cut21->SetPoint(2,65.75,129.5);
cut21->SetPoint(3,65.75,124.5);
cut21->SetPoint(4,60.75,124.5);
int integral1 = cut21->IntegralHist(Artificial);
cout << "net area peak before random=" << " "<< integral1 <<endl;
int integral3 = cut21->IntegralHist(Artificial1);
cout << "net area peak before =" << " "<< integral3 <<endl;
int integral2 = cut21->IntegralHist(Bg_Sub);
cout << "net area peak after =" <<" "<< integral2 <<endl;
TCutG *cut1 = new TCutG("cut1",5);
cut1->SetPoint(0,120,40);
cut1->SetPoint(1,200,40);
cut1->SetPoint(2,200,100);
cut1->SetPoint(3,120,100);
cut1->SetPoint(4,120,40);
int integral4 = cut1->IntegralHist(Artificial);
cout << "net area peak before =" << " "<< integral4 <<endl;
//int integral5 = cut1->IntegralHist(Artificial1);
//cout << "net area peak before =" << " "<< integral5 <<endl;
int integral6 = cut1->IntegralHist(Bg_Sub);
cout << "net area peak after =" <<" "<< integral6 <<endl;
//in.close();
out.close();
out2.close();
b = new TBrowser();
f->Write();
}
Plz help I straggled here for one week
You need to find which call to the “exp” function in your code gets too big parameter (from the reported error, I assume it tried to calculate “exp(999.)”, while the maximum allowed value for double precision floating point numbers is “exp(709.)”).
Thank you for your reply.
I changed the double to float
same problem
what I have to do
You need to find the place where you try to calculate “exp(999.)” yourself and protect your source code against usage of “improper” parameters’ values.
C/C++ “double” usually means: Double-precision floating-point format
C/C++ “float” usually means: Single-precision floating-point format
BTW. Your “SomeMacro.cxx” file should look like:
#include <string.h>
// ... all includes "outside" of the "SomeMacro" function's source code
#include "TGClient.h"
void SomeMacro(void) { // use the same name as the file's basename
// note: do NOT execute gROOT->Reset(); unless you REALLY know what it does
ofstream out, out2;
// ... the remaining source code comes here
f->Write();
}
Thank you so much
I read what did you send me but unfortunately I did not understand I am still beginner.
I tried this but it did not work
I am sorry for this
Really I need your help
I am stuck here for one week
could you plz help
same problem I got
As Wile E suggested:
One way is to compile your file (.x MyCode.C+ - note the trailing ‘+’). For that you need to transform your code into proper C++ which should be worth it either way.
Then you can attach a debugger do set a breakpoint at the exp calls and see which one fails.
Or you type .T in ROOT and see which line is executed, one by one, and at which point CINT is complaining about the wrong input to exp.
Cheers, Axel.
Try to use “safe_exp” defined as follows:
#include <iostream>
#include <cmath>
double safe_exp(double x) {
#if 0 /* 0 or 1 */
if (x < -745) x = -745; else if (x > 709) x = 709; else if (std::isnan(x)) x = 0;
#endif /* 0 or 1 */
double v = std::exp(x);
if ((!std::isfinite(x)) || (!std::isfinite(v)))
std::cout << "Warning : exp( " << x << " ) = " << v << std::endl;
return v;
}
1 Like
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.