Index: trunk/examples/createBSatBDependenceTable.cpp
===================================================================
--- trunk/examples/createBSatBDependenceTable.cpp (revision 472)
+++ trunk/examples/createBSatBDependenceTable.cpp (revision 473)
@@ -1,133 +1,295 @@
//==============================================================================
// createBSatBDependenceTable.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to create impact parameter lookup table.
// [Developer only]
-//==============================================================================
-#include
+//==============================================================================
+#include
#include
#include
#include
#include
#include "TableGeneratorNucleus.h"
#include "TableGeneratorSettings.h"
#include "TFile.h"
#include "TH2F.h"
#include "TVector3.h"
#include "Constants.h"
#include "DipoleModelParameters.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
-double overlapFunctionT(double);
-void createLookUpTableOnFile(unsigned int, unsigned int);
+double overlapFunctionT0(double);
+double overlapFunctionT1bSat(double);
+double overlapFunctionT1bNonSat(double);
+
+void createLookUpTableOnFile0(unsigned int, int, int, int);
+void createLookUpTableOnFile1(unsigned int, int, int, int, int);
DipoleModelParameters *gParameters;
int main(int argc, const char* argv[])
{
- if (argc != 3) {
- cout << "Usage: " << argv[0] << " A numberOfConfigurations" << endl;
+ if (argc != 7) {
+ cout << "Usage: " << argv[0] << " A numberOfConfigurations starting ending structure model" << endl;
+ cout << "structure:"<numConf) ending=numConf;
+
cout<<"Create table for A="<BG(); // GeV^-2
double arg = (b*b/hbarc2) / (2*BG);
return 1/(2*M_PI*BG) * exp(-arg);
-}
-
+}
-void createLookUpTableOnFile(unsigned int A, unsigned int numberOfConfigurations)
+double overlapFunctionT1bSat(double b)
{
+ // For bSat subnucleonic fluctuations, use a modified Gaussian (BE),
+ //
+ // Parameters for COM shifted STU, bSat:
+ // Bqc = 4.5
+ // Bq = 1.26
+ // Nq = 3
+ // Sg = 0.3
+ // Sigma = 0.5
+ double Bq = 1.26;
+ double Sg = 0.3;
+ double Nq=3; // GeV^-2
+ double arg = (b*b/hbarc2) / (2*Bq);
+ double bDep= 1/(2*M_PI*Bq*Nq) * 1/(exp(arg)-Sg);
+
+ return bDep;
+}
+
+double overlapFunctionT1bNonSat(double b)
+{
+ // For bSat subnucleonic fluctuations, use a modified Gaussian (BE),
+ //
+ // Parameters for COM shifted STU, bNonSat:
+ // Bqc = 4.5
+ // Bq = 1.26
+ // Nq = 3
+ // Sigma = 0.5 (For Qs fluctuations)
+
+ double Bq = 1.26;
+ double Nq=3; // GeV^-2
+ double arg = (b*b/hbarc2) / (2*Bq);
+ double bDep= 1/(2*M_PI*Bq*Nq) * exp(-arg);
+
+ return bDep;
+}
+
+void createLookUpTableOnFile0(unsigned int A, int numConf,
+ int starting, int ending)
+{
//
// A function to create a lookup table of nucleon configurations and save them on a file
// This function is supposed to be independent from the rest of the code
// and be called in the Main program with the nucleus number as an argument.
//
+ PR(starting);
+ PR(ending);
+ int numberOfConfigurations=ending-starting; //#TT Test
TableGeneratorNucleus myNucleus(A);
- // myNucleus.init(A);
- TableGeneratorSettings::instance()->setSeed(42); // sets also seed for ROOT
+ myNucleus.init(A);
+ TableGeneratorSettings::instance()->setSeed(0); // sets also seed for ROOT
vector hConfigurations(numberOfConfigurations);
int numBbins=1e3, numAngleBins=1e3;
double bRange=upperIntegrationLimit*myNucleus.radius()*1.1, angleRange=2*M_PI;
ostringstream histoName;
- for(unsigned int iNuclei=0; iNucleiGetXaxis()->GetBinCenter(ib);
for(int iphi=1; iphi <= numAngleBins; iphi++) {
double phi=hConfigurations[iNuclei]->GetYaxis()->GetBinCenter(iphi);
TVector3 bvector=TVector3(b*cos(phi), b*sin(phi), 0.);
double sumOfT=0.;
for(unsigned int iA=0; iA < A; iA++){
- sumOfT+=overlapFunctionT((bvector-myNucleus.configuration.at(iA).position()).Perp());
+ sumOfT+=overlapFunctionT0((bvector-myNucleus.configuration.at(iA).position()).Perp());
} //for iA
hConfigurations[iNuclei]->SetBinContent(ib, iphi, sumOfT);
}//for iphi
} //for ib
- if ((numberOfConfigurations/10) && iNuclei % (numberOfConfigurations/10) == 0) {
- cout<Write();
delete hConfigurations[i];
hConfigurations[i] = 0;
}
lufile->Close();
cout<<" done, bye!"<setSeed(0); // sets also seed for ROOT
+
+ vector hConfigurations(numConf);
+ int numBbins=1e3, numAngleBins=1e3, Nq=3;
+ double bRange=upperIntegrationLimit*myNucleus.radius()*1.1, angleRange=2*M_PI;
+ ostringstream histoName;
+
+ double sigma=0.5; //should also be changed in TableGeneratorNucleus
+ double satNorm=exp(sigma*sigma/2);
+
+ for(int iNuclei=starting; iNuclei positions;
+ positions.clear();
+ vector satScales;
+ satScales.clear();
+ //Generate a Nucleon configuration:
+ while(!myNucleus.generate()){}
+ for(unsigned int iA=0; iAGetXaxis()->GetBinCenter(ib);
+ for(int iphi=1; iphi <= numAngleBins; iphi++) {
+ double phi=hConfigurations[iNuclei]->GetYaxis()->GetBinCenter(iphi);
+ TVector3 bvector=TVector3(b*cos(phi), b*sin(phi), 0.);
+ double sumOfT_i=0, sumOfT_j;
+ for(unsigned int iA=0; iA < A; iA++){
+ sumOfT_j=0;
+ for ( int jQ=0; jQSetBinContent(ib, iphi, sumOfT_i);
+ }//for iphi
+ } //for ib
+ if ((numberOfConfigurations/10) && iNuclei-starting % (numberOfConfigurations/10) == 0) {
+ cout<Write();
+ delete hConfigurations[i];
+ hConfigurations[i] = 0;
+ }
+ lufile->Close();
+ cout<<" done, bye!"<