157 lines
3.8 KiB
C++
157 lines
3.8 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: Histogram_RAW_File.cc
|
|
|
|
Purpose: Produce a histogram of matrix elements from a RAW file.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
if (argc != 7) {
|
|
cout << "Argument needed: rawfile, iKmin, iKmax, omega, domega, logscale_int" << endl;
|
|
ABACUSerror("");
|
|
}
|
|
|
|
int n = 1;
|
|
const char* rawfilename = argv[n++];
|
|
int iKmin = atoi(argv[n++]);
|
|
int iKmax = atoi(argv[n++]);
|
|
DP om = atof(argv[n++]);
|
|
DP dom = atof(argv[n++]);
|
|
int logscale_int = atoi(argv[n++]); // see below for meaning
|
|
|
|
ifstream RAW_infile;
|
|
RAW_infile.open(rawfilename);
|
|
if (RAW_infile.fail()) {
|
|
cout << rawfilename << endl;
|
|
ABACUSerror("Could not open RAW_infile... ");
|
|
}
|
|
|
|
stringstream outfilename;
|
|
string outfilename_string;
|
|
outfilename << rawfilename << "_hist";
|
|
if (iKmin != iKmax) outfilename << "_iKmin_" << iKmin << "_iKmax_" << iKmax;
|
|
outfilename << "_om_" << om << "_dom_" << dom << "_li_" << logscale_int;
|
|
outfilename_string = outfilename.str();
|
|
const char* outfilename_c_str = outfilename_string.c_str();
|
|
|
|
ofstream outfile;
|
|
outfile.open(outfilename_c_str);
|
|
outfile.precision(16);
|
|
|
|
|
|
|
|
// Use the binning in current ABACUS version's threads,
|
|
// with logscale unit = 2^{1/64} \simeq 1.01
|
|
int npts = 6400/logscale_int;
|
|
DP logscale_used = logscale_int * (1.0/64) * log(2.0);
|
|
|
|
DP omega_min = om - 0.5* dom;
|
|
DP omega_max = om + 0.5* dom;
|
|
|
|
DP omega;
|
|
int iK;
|
|
DP FF;
|
|
DP dev;
|
|
string label;
|
|
|
|
Vect<int> nFF (0, npts);
|
|
|
|
int nread = 0;
|
|
int nfound = 0;
|
|
int il = 0;
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
|
|
nread++;
|
|
|
|
if (iK >= iKmin && iK <= iKmax && omega > omega_min && omega < omega_max) {
|
|
nfound++;
|
|
il = int(-log(FF*FF)/logscale_used);
|
|
if (il < 0) il=0;
|
|
if (il > npts) continue;
|
|
nFF[il] += 1;
|
|
}
|
|
}
|
|
RAW_infile.close();
|
|
|
|
if (nfound == 0) {
|
|
cout << "Didn't find any contributing state." << endl;
|
|
return(1);
|
|
}
|
|
|
|
for (int i = 0; i < npts; ++i) {
|
|
if (i > 0) outfile << endl;
|
|
outfile << i << "\t" << exp(-i * logscale_used) << "\t" << nFF[i] << "\t" << DP(nFF[i])/nfound;
|
|
}
|
|
outfile.close();
|
|
|
|
|
|
// Now that we know the total number of entries, produce animated data:
|
|
// the total number of entries is sliced into 100 pieces
|
|
// and binned in the same way as before
|
|
|
|
int stepsize = nfound/100;
|
|
|
|
outfilename << "_sliced";
|
|
outfilename_string = outfilename.str();
|
|
const char* sliced_c_str = outfilename_string.c_str();
|
|
|
|
ofstream sliced_file;
|
|
sliced_file.open(sliced_c_str);
|
|
sliced_file.precision(16);
|
|
|
|
RAW_infile.open(rawfilename);
|
|
if (RAW_infile.fail()) {
|
|
cout << rawfilename << endl;
|
|
ABACUSerror("Could not open RAW_infile... ");
|
|
}
|
|
|
|
nread = 0;
|
|
nfound = 0;
|
|
int nsteps = 0;
|
|
for (int i = 0; i < npts; ++i) nFF[i] = 0;
|
|
|
|
while (RAW_infile.peek() != EOF) {
|
|
RAW_infile >> omega >> iK >> FF >> dev >> label;
|
|
|
|
nread++;
|
|
|
|
if (iK >= iKmin && iK <= iKmax && omega > omega_min && omega < omega_max) {
|
|
nfound++;
|
|
il = int(-log(FF*FF)/logscale_used);
|
|
if (il < 0) il=0;
|
|
if (il > npts) continue;
|
|
nFF[il] += 1;
|
|
}
|
|
|
|
if (nfound == stepsize) {// || RAW_infile.peek() == EOF) {
|
|
if (nsteps++ > 0) sliced_file << endl;
|
|
sliced_file << DP(nFF[0])/stepsize;
|
|
for (int i = 1; i < npts; ++i) sliced_file << "\t" << DP(nFF[i])/stepsize;
|
|
for (int i = 0; i < npts; ++i) nFF[i] = 0;
|
|
nfound = 0;
|
|
}
|
|
}
|
|
|
|
RAW_infile.close();
|
|
|
|
sliced_file.close();
|
|
|
|
return(0);
|
|
}
|