rtang/SOLARIS_Analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
SOLARIS_Analysis/Cleopatra/ExtractXSec.h

547 lines
16 KiB
C++
Raw Blame History

/***********************************************************************
*
* This is ExtractXSec for *.out for Ptolemy
*
* This program will extract cross section distribution from Ptolmey output.
* save as *.Xsec.txt and *.root for distribution
*
* save ExPtolemy.txt for excitation energies and total X-section
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ ExtractXSec.C -o ExtractXSec `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <stdlib.h>
#include <vector>
#include <TROOT.h>
#include <TFile.h>
#include <TString.h>
#include <TMath.h>
#include <TGraph.h>
#include <TMacro.h>
#include <TF1.h>
#include <TObjArray.h>
#include "../Armory/AnalysisLib.h"
using namespace std;
TObjArray * gList = NULL;
double distfunct(double *x, double *par){
TGraph * gTemp = (TGraph *) gList->At(par[0]);
return (gTemp->Eval(x[0]))* sin(x[0] * TMath::DegToRad());
}
bool isFloat( string str ) {
int length = str.length();
for( int i = 0; i < length; i++){
string it = str.substr(i,1);
if( it == " " || it == "."|| it == "1"|| it == "2"||
it == "3"|| it == "4"|| it == "5"|| it == "6"|| it == "7"|| it == "8"|| it == "9"|| it == "0"){
continue;
}else{
return false;
}
}
return true;
}
int ExtractXSec (string readFile, int indexForElastic=1) {
//indexForElastic = 1 ; for Ratio
//indexForElastic = 2 ; for Total
//indexForElastic = 3 ; for Rutherford
//--- open file.out
ifstream file_in;
file_in.open(readFile.c_str(), ios::in);
if( !file_in){
printf("Unable to open %s \n", readFile.c_str());
exit(1);
}
//---- variable for Xsec
vector<string> title;
vector<vector<double>> dataMatrix;
vector<double> angle;
vector<double> Ex;
vector<double> totalXsec;
vector<string> reaction;
double angleMin = 0;
double angleMax = 0;
double angleStep = -1;
//================================== read file.out
string line;
int lineNum = 0;
vector<double> dataXsec;
bool startExtract = false;
bool angleFilled = false;
int numCal = 0;
int reactionFlag = 0;
bool preFindForElastic = false;
printf("======================================================\n");
while(getline(file_in, line)){
lineNum ++;
//printf("%d , %s\n", lineNum, line.c_str());
size_t pos;
string findStr;
int findLen;
//---- find Title
findStr = "$============================================";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos){
title.push_back(line.substr(pos + findLen+1));
pos = title.back().find("(");
string ex = title.back().substr(3, pos-3);
Ex.push_back( atof(ex.c_str()) );
//clear dataXsex for new data set
dataXsec.clear();
numCal ++;
continue;
}
//---- find Reaction
findStr = "0INPUT... CHANNEL";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
reaction.push_back( line.substr(pos + findLen + 1) );
reactionFlag = 1; // 1 for (d,d) or (p,p)
///printf("%d |----------- (d,d), %d\n", lineNum, reactionFlag);
continue; // nextline
}
findStr = "REACTION:";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
reaction.push_back( line.substr(pos + findLen + 1) );
reactionFlag = 2; // 2 for (d,p) or (p,d)
///printf("%d |----------- (d,p), %d\n", lineNum, reactionFlag);
continue; // nextline
}
//----- pre find for Elastic scattering
findStr = "0 OPTICAL MODEL SCATTERING FOR THE OUTGOING CHANNEL";
pos = line.find(findStr);
if( pos != string::npos ) {
preFindForElastic = true;
///printf("%d |----------- pre Elastics Flag : %d\n", lineNum, preFindForElastic);
continue;
}
//----- find angle stetting when not known
if( angleStep == -1 ){
findStr = "anglemin=";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos){
size_t pos1 = line.find("anglemax=");
angleMin = atof( line.substr(pos+findLen, pos1 - pos - findLen).c_str() );
pos = pos1;
pos1 = line.find("anglestep=");
angleMax = atof( line.substr(pos+findLen, pos1 - pos - findLen).c_str() );
angleStep = atof( line.substr(pos1 + findLen+1).c_str() );
//printf("%d |---- angle range found.\n", lineNum);
}
continue; //nextline
}
//----- check if start extracting Xsec or not
findStr = "dumpdumpdump";
if( reactionFlag == 1 && preFindForElastic ){
findStr = "C.M. LAB RUTHERFORD";
}
if( reactionFlag == 2 ){
findStr = "0 C.M. REACTION REACTION LOW L HIGH L % FROM";
}
pos = line.find(findStr);
if( pos != string::npos ) {
startExtract = true;
//printf("%d |----- start extraction \n", lineNum);
continue;
}
//----- start extracting Xsec
if( startExtract ){
if( line.length() < 20 ) continue;
///printf("%d |%s \n", line.length(), line.c_str());
double numAngle, numXsec;
if( reactionFlag == 1 ){ // Elastics (d,d) or (p,p)
numAngle = atof( line.substr(0, 7).c_str());
if( indexForElastic == 1 ){
numXsec = atof( line.substr(15, 15).c_str());
}else if( indexForElastic == 2 ){
if ( line.length() > 60 ) {
numXsec = atof( line.substr(30, 13).c_str());
}else{
numXsec = -404;
}
}else{
if ( line.length() > 60 ) {
numXsec = atof( line.substr(57, 14).c_str());
}else{
numXsec = -404;
}
}
}
if( reactionFlag == 2 ){ // InElastics (d,p) or (p,d)
if( isFloat( line.substr(0, 7) ) ){
numAngle = atof( line.substr(0, 7).c_str());
numXsec = atof( line.substr(7, 19).c_str());
}else{
numAngle = -1.0;
numXsec = -1.0;
}
}
if( numAngle != numXsec && numXsec > 0. ){
if( !angleFilled ) angle.push_back(numAngle);
dataXsec.push_back(numXsec);
///printf(" %f , %f \n", angle.back(), dataXsec.back());
}
}
//------ find total Xsec, if found stop extraction
findStr = "dumpdumpdump";
if( reactionFlag == 1 && preFindForElastic){
findStr = "0TOTAL REACTION CROSS SECTION =";
}
if( reactionFlag == 2){
findStr = "0TOTAL:";
}
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
totalXsec.push_back( atof(line.substr(pos + findLen).c_str()) );
printf("%2d | %s | total Xsec(4pi): %f mb \n", numCal, title.back().c_str(), totalXsec.back());
//push back dataXsec to dataMatrix
dataMatrix.push_back(dataXsec);
//printf("%d |----- end extraction \n", lineNum);
angleFilled = true;
startExtract = false;
reactionFlag = 0;
preFindForElastic = false;
continue;
}
}
file_in.close();
//================================== summary
printf("====================== Total number of line read : %d \n", lineNum);
printf("Angle : %5.2f, %5.2f | step : %5.2f \n", angleMin, angleMax, angleStep);
printf("Number of Angle read : %lu \n", angle.size());
printf("Number of data read : %lu \n", dataXsec.size());
printf("Number of Reaction : %lu \n", reaction.size());
printf("----------------------------- list of Calculation \n");
//... find suitable lenght for displaying reaction string
int reactionStrLen = 0;
for( int i = 0; i < numCal ; i++){
int len = reaction[i].length();
if( len > reactionStrLen ) reactionStrLen = len;
}
vector<double> partialXsec;
partialXsec.clear();
for( int i = 0; i < numCal ; i++){
double partialSumXsec = 0.0;
for( int j = 0; j < (dataMatrix[i]).size() ; j++ ){
//double theta = (angle[j] + angleStep/2.) * TMath::DegToRad();
double theta = (angle[j]) * TMath::DegToRad();
double dTheta = angleStep * TMath::DegToRad();
double phi = TMath::TwoPi();
partialSumXsec += dataMatrix[i][j] * sin( theta ) * dTheta * phi ;
}
partialXsec.push_back(partialSumXsec);
size_t pos = title[i].find(")");
printf("%*s| %s | Xsec(%3.0f-%3.0f deg) : %f mb\n", reactionStrLen + 3, reaction[i].c_str(), title[i].substr(pos+1).c_str(), angleMin, angleMax, partialSumXsec);
}
printf("---------------------------------------------------\n");
//================================== save file.Xsec.txt
string saveFileName = readFile;
int len = saveFileName.length();
saveFileName = saveFileName.substr(0, len - 4);
saveFileName += ".Xsec.txt";
printf("Output : %s \n", saveFileName.c_str());
FILE * file_out;
file_out = fopen(saveFileName.c_str(), "w+");
for( int i = 0; i < numCal ; i++){
fprintf(file_out, "#%14s\n", reaction[i].c_str());
}
int space = 19;
fprintf(file_out, "%8s\t", "Angle");
for( int i = 0; i < numCal ; i++){
fprintf(file_out, "%*s", space, title[i].c_str());
}
fprintf(file_out, "\n");
for( int i = 0; i < angle.size() ; i ++){
fprintf(file_out, "%8.3f\t", angle[i]);
for( int j = 0; j < numCal ; j++){
fprintf(file_out, "%*f", space, dataMatrix[j][i]);
}
fprintf(file_out, "\n");
}
fclose(file_out);
//================================== Make TMacro for ExList
TMacro ExList;
TMacro ReactList;
ExList.AddLine("#---Ex relative_xsec SF sigma_in_MeV");
for( int i = 0; i < numCal ; i++){
ExList.AddLine(Form("%9.5f %9.5f 1.0 0.000", Ex[i], partialXsec[i]));
ReactList.AddLine(reaction[i].c_str());
}
//================================== Save in ROOT
len = saveFileName.length();
saveFileName = saveFileName.substr(0, len - 9);
TString fileName = saveFileName;
fileName += ".root";
printf("Output : %s \n", fileName.Data());
TFile * fileOut = new TFile(fileName, "RECREATE" );
gList = new TObjArray(); ///no SetTitle() method for TObjArray
gList->SetName("TGraph of d.s.c");
TObjArray * fList = new TObjArray();
fList->SetName("TF1 of distributions = d.s.c. * sin()");
TGraph ** gGraph = new TGraph *[numCal];
TF1 ** dist = new TF1*[numCal];
for( int i = 0; i < numCal ; i++){
gGraph[i] = new TGraph();
TString name = reaction[i];
name += "|";
name += title[i];
gGraph[i]->SetName(name);
for( int j = 0; j < angle.size() ; j++){
gGraph[i]->SetPoint(j, angle[j], dataMatrix[i][j]);
}
gList->Add(gGraph[i]);
name.Form("dist%d", i);
dist[i] = new TF1(name, distfunct, angleMin, angleMax, 1 );
dist[i]->SetParameter(0, i);
fList->Add(dist[i]);
}
gList->Write("thetaCM_TGraph", 1);
fList->Write("thetaCM_TF1", 1);
ExList.Write("ExList");
ReactList.Write("ReactionList");
fileOut->Write();
fileOut->Close();
printf("---------------------------------------------------\n");
return 0;
}
int ExtractXSecFromText(string readFile){
//read file
ifstream file_in;
file_in.open(readFile.c_str(), ios::in);
if( !file_in){
printf("Unable to open %s \n", readFile.c_str());
exit(1);
}
//---- variable for Xsec
vector<vector<double>> dataMatrix;
vector<double> angle;
vector<double> Ex;
//================================== read file.out
string line;
int lineNum = 0;
vector<double> dataXsec;
int numCal = 0;
printf("======================================================\n");
bool headerDone = false;
while(getline(file_in, line)){
lineNum ++;
if( line.substr(0,1) == "#" ) continue;
//printf("%d | %s\n", lineNum, line.c_str());
//after the comment line, the next line must be column name
vector<string> header= AnalysisLib::SplitStr(line, " ");
//printf("---%lu #", header.size());
//for( int i = 0; i < header.size(); i++){
// printf("%s|", header[i].c_str());
//}
//printf("\n");
if ( !headerDone ) {
for( int i = 1 ; i < header.size(); i++){
string energy = header[i].substr(3, header[i].length());
Ex.push_back(atof(energy.c_str()));
//printf("%f---\n", Ex.back());
}
headerDone = true;
}else{
vector<double> temp;
for( int i = 0; i < header.size(); i++){
if( i == 0 ) angle.push_back(atof(header[i].c_str()));
if( i > 0 ) temp.push_back(atof(header[i].c_str()));
}
dataMatrix.push_back(temp);
}
}
file_in.close();
double angleMin = angle.front();
double angleMax = angle.back();
double angleStep = (angleMax - angleMin)/(angle.size()-1);
//------ print out read data
printf("====================== data read as:\n");
printf("%5s, ", "Ex");
for(int i = 0; i < Ex.size(); i++){
printf("%6.3f|", Ex[i]);
}
printf("\n");
for(int i = 0; i < dataMatrix.size(); i++){
printf("%5.2f, ", angle[i]);
for( int j = 0; j < dataMatrix[i].size(); j++) printf("%6.3f|", dataMatrix[i][j]);
printf("\n");
}
printf("---------------------------------\n");
printf("angle min :%f\n", angleMin);
printf("angle max :%f\n", angleMax);
printf("angle step:%f\n", angleStep);
printf("---------------------------------\n");
//------- calculate summed cross section
vector<double> partialXsec(Ex.size());
for( int i = 0; i < dataMatrix.size() ; i++){
for( int j = 0; j < (dataMatrix[i]).size() ; j++ ){
//double theta = (angle[j] + angleStep/2.) * TMath::DegToRad();
double theta = (angle[i]) * TMath::DegToRad();
double dTheta = angleStep * TMath::DegToRad();
double phi = TMath::TwoPi();
partialXsec[j] += dataMatrix[i][j] * sin( theta ) * dTheta * phi ;
}
}
for(int i = 0; i < Ex.size(); i++){
printf("Ex=%6.3f| Xsec(%3.0f-%3.0f deg) : %f mb\n", Ex[i] , angleMin, angleMax, partialXsec[i]);
}
printf("---------------------------------------------------\n");
//================================== save *.Ex.txt
string saveExName = readFile;
int len = saveExName.length();
saveExName = saveExName.substr(0, len - 4);
saveExName += ".Ex.txt";
printf("Output : %s \n", saveExName.c_str());
FILE * file_Ex;
file_Ex = fopen(saveExName.c_str(), "w+");
fprintf(file_Ex, "//Generated_by_ExtractXsec.h___Ex___Xsec___SF____sigma\n");
for( int i = 0; i < Ex.size() ; i++){
fprintf(file_Ex, "%9.5f %9.5f 1.0 0.000\n", Ex[i], partialXsec[i]);
}
fprintf(file_Ex, "#=====End_of_File\n");
fclose(file_Ex);
//================================== Save in ROOT
string saveFileName = readFile;
len = saveFileName.length();
saveFileName = saveFileName.substr(0, len - 4);
TString fileName = saveFileName;
fileName += ".root";
printf("Output : %s \n", fileName.Data());
TFile * fileOut = new TFile(fileName, "RECREATE" );
gList = new TObjArray();
gList->SetName("TGraph of d.s.c");
TObjArray * fList = new TObjArray();
fList->SetName("TF1 of distributions = d.s.c. * sin()");
TGraph ** gGraph = new TGraph *[numCal];
TF1 ** dist = new TF1*[numCal];
for( int i = 0; i < Ex.size() ; i++){
gGraph[i] = new TGraph();
TString name ;
name.Form("Ex=%6.3f MeV", Ex[i]);
gGraph[i]->SetName(name);
for( int j = 0; j < angle.size() ; j++){
gGraph[i]->SetPoint(j, angle[j], dataMatrix[j][i]);
}
gList->Add(gGraph[i]);
name.Form("dist%d", i);
dist[i] = new TF1(name, distfunct, angleMin, angleMax, 1 );
dist[i]->SetParameter(0, i);
fList->Add(dist[i]);
//delete tempFunc;
}
gList->Write("qList", 1);
fList->Write("pList", 1);
fileOut->Write();
fileOut->Close();
printf("---------------------------------------------------\n");
return 0;
}
Reference in New Issue View Git Blame Copy Permalink