TCS and BH signal event generator

src/programs/Simulation/SConscript

@@ -3,7 +3,7 @@ import os,sbms
 
 Import('*')
 
-subdirs = ['genr8', 'GEN2HDDM', 'genr8_2_hddm', 'HDGeant', 'mcsmear', 'bggen', 'gen_2k', 'gen_2pi', 'gen_2pi_amp', 'gen_2pi_primakoff','gen_3pi', 'gen_pi0', 'gen_omega_3pi', 'gen_omega_radiative' , 'nullgen', 'gen_amp', 'gen_amp_V2', 'BGRate_calc', 'genEtaRegge', 'gen_ee', 'gen_ee_hb', 'genScalarRegge', 'gen_compton', 'gen_omegapi', 'gen_vec_ps', 'gen_compton_simple', 'gen_primex_eta_he4', 'gen_whizard', 'MC_GEN', 'bggen_jpsi', 'gen_2pi0_primakoff', 'gen_EtaPb', 'bggen_upd', 'gen_generic_root']
+subdirs = ['genr8', 'GEN2HDDM', 'genr8_2_hddm', 'HDGeant', 'mcsmear', 'bggen', 'gen_2k', 'gen_2pi', 'gen_2pi_amp', 'gen_2pi_primakoff','gen_3pi', 'gen_pi0', 'gen_omega_3pi', 'gen_omega_radiative' , 'nullgen', 'gen_amp', 'gen_amp_V2', 'BGRate_calc', 'genEtaRegge', 'gen_ee', 'gen_ee_hb', 'genScalarRegge', 'gen_compton', 'gen_omegapi', 'gen_tcs_bh', 'gen_vec_ps', 'gen_compton_simple', 'gen_primex_eta_he4', 'gen_whizard', 'MC_GEN', 'bggen_jpsi', 'gen_2pi0_primakoff', 'gen_EtaPb', 'bggen_upd', 'gen_generic_root']
 
 
 # only build if	    EvtGen is installed

src/programs/Simulation/gen_tcs_bh/CMakeLists.txt

@@ -0,0 +1,60 @@
+# Minimum CMake version required
+cmake_minimum_required(VERSION 3.10)
+
+# Project name and C++ standard
+project(gen_tcs_bh CXX)
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Compiler flags
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -fPIC -ansi -Wall -Wno-overloaded-virtual -Wno-long-long -fno-common -O2")
+
+# Find ROOT
+find_package(ROOT REQUIRED COMPONENTS RIO Net Hist Tree)
+include_directories(${ROOT_INCLUDE_DIRS})
+link_directories(${ROOT_LIBRARY_DIRS})
+add_definitions(${ROOT_DEFINITIONS})
+
+# Source files
+set(ANA_C
+    FillTables.cc
+    RadProcess.cc
+    TableProvider.cc
+    TabUtils.cc
+    interpol.cc
+    Utils.cc
+    Options_tcs.cc
+    FormFactors.cc
+    TreatOptions.cc
+    PDFparam.cc
+    Constants.cc
+    Polar.cc
+    PDGinfo.cc
+    Kinematics.cc
+    LeptonPair.cc
+    PartUtils.cc
+    ReactionKinTwo.cc
+)
+
+# Main executable
+add_executable(gen_tcs_bh main.cc ${ANA_C})
+
+# Link ROOT libraries
+target_link_libraries(gen_tcs_bh ${ROOT_LIBRARIES})
+
+
+# Custom target for cleaning files, copying scripts, and executing the script
+add_custom_target(custom_clean_and_setup
+    # Cleaning up generated files
+    COMMAND ${CMAKE_COMMAND} -E remove -f *.o *~ main includes/*.o *.o Data/*
+    COMMENT "Cleaning up generated files"
+
+    # Copying files
+    COMMAND ${CMAKE_COMMAND} -E copy ../set.csh ${CMAKE_BINARY_DIR}/set.csh
+#    COMMAND ${CMAKE_COMMAND} -E copy ../TCS.input ${CMAKE_BINARY_DIR}/TCS.input
+#    COMMAND ${CMAKE_COMMAND} -E copy ../PS_EEPHOTO_FIX.input ${CMAKE_BINARY_DIR}/PS_EEPHOTO_FIX.input
+    # Executing the script
+    COMMAND ./set.csh
+    WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
+    COMMENT "Copying set.csh, then executing set.csh"
+)

src/programs/Simulation/gen_tcs_bh/Constants.cc

@@ -0,0 +1,2 @@
+#define Constants_cxx
+#include "Constants.h"

src/programs/Simulation/gen_tcs_bh/Constants.h

@@ -0,0 +1,18 @@
+#ifndef Constants_h
+#define Constants_h
+#include <iostream>
+#include <fstream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+const double alphaEM= 0.007297352; //1./137.036, 
+const double m_el=510.998910e-6;
+const double PI=3.141592, M_Neutron=0.939566, M_Proton=0.9382723 , M_Muon=0.1056583;
+const double M_pizero=0.1349766, M_piplus=0.13957, M_eta=0.547853, M_rho=0.77549, M_omega=0.78265, M_etac=2.9803, M_jpsi=3.09916;
+const double hbarsquare = 0.0003894; // convert GeV/c2 to barn // check unit !!! factor 100
+
+#endif
+
+
+

src/programs/Simulation/gen_tcs_bh/FillTables.cc

@@ -0,0 +1,202 @@
+#define FillTables_cxx
+#include "FillTables.h"
+#include <limits>
+#include <unistd.h>
+using namespace std;
+
+int FillTables::FillTable(int reaction, TableProvider& table_provider){
+
+	genSettings_t genSettings0;
+	TString Tablepath0 = genSettings0.xsecTablepath;
+        printf("Start to fill all the tables ...");
+
+        ifstream infile1x,infile1z;
+        string ff1x, ff1z, spoub;
+        char proc[24];
+        double em=0.,tm=0.,Qpm=0.,Thm=0.,phim=0.,psim=0.; //Qm=0.,PLm=0.,
+        int ii=0, jj=0, kk=0, ll=0, mm=0,nn=0,countercheck=0; //oo=0,
+        int nn_xx=-1,nn_tt=-1,nn_qp=-1,nn_th=-1,nn_pc=-1,nn_ps=-1; //,nn_qq=-1,nn_pl=-1
+        double res[15]={0.};
+        float CosThetaMax=1.;
+
+        if (reaction==1) strcpy(proc, "tcs");
+
+        cout<<"reaction index: "<<reaction<<endl;
+
+	if (strcmp(proc, "tcs")==0){
+
+                if (protonorneutron==1){
+                        if (targetpoldir==1 || targetpoldir==2){
+                                cout<<"Read table for circularly polarized photon and transversely polarized target"<<endl;
+                                ff1x=Form(Tablepath0+"/grid_table_tcs_circperp.dat");
+                                //ff1x=Form("Data/grid_table_tcs_circperp.dat");
+                        } else {
+                                if (beampolartype==1){
+                                        cout<<"Read table for linearly polarized photon and longitudinally polarized target"<<endl;
+                                        ff1x=Form(Tablepath0+"/grid_table_tcs_linlong.dat");
+                                        //ff1x=Form("Data/grid_table_tcs_linlong.dat");
+                                } else {
+                                        cout<<"Read table for circularly polarized photon and longitudinally polarized target"<<endl;
+                                        ff1x=Form(Tablepath0+"/grid_table_tcs_circlong.dat");
+                                        //ff1x=Form("Data/grid_table_tcs_circlong.dat");
+                                }
+                        }
+                } else if (protonorneutron==2){
+                        cout<<"Read cross section off neutron"<<endl;
+                        if (targetpoldir==1 || targetpoldir==2){
+                                cout<<"Read table for circularly polarized photon and transversely polarized target"<<endl;
+                                ff1x=Form(Tablepath0+"/grid_table_tcs_neutron_circperp.dat");
+                                //ff1x=Form("Data/grid_table_tcs_neutron_circperp.dat");
+                        } else {
+                                if (beampolartype==1){
+                                        cout<<"Read table for linearly polarized photon and longitudinally polarized target"<<endl;
+                                        ff1x=Form(Tablepath0+"/grid_table_tcs_neutron_linlong.dat");
+                                        //ff1x=Form("Data/grid_table_tcs_neutron_linlong.dat");
+                                } else {
+                                        cout<<"Read table for circularly polarized photon and longitudinally polarized target"<<endl;
+                                        ff1x=Form(Tablepath0+"/grid_table_tcs_neutron_circlong.dat");
+                                        //ff1x=Form("Data/grid_table_tcs_neutron_circlong.dat");
+                                }
+                        }
+
+                } else {
+                        cout<<"ERROR: target = proton or neutron"<<endl;
+                        return 5;
+                }
+                infile1x.open(ff1x.c_str());
+                if (!infile1x) {
+                        cout<<"no input for TCS cross section"<<endl;
+                        cout<<"EXIT program with ERROR"<<endl;
+                        return 0;
+                }
+                if (targetpoldir==1 || targetpoldir==2){
+                        infile1x>>nn_xx>>nn_tt>>nn_qp>>nn_th>>nn_pc>>nn_ps>>CosThetaMax;
+                } else {
+                        if (beampolartype==1){
+                                infile1x>>nn_xx>>nn_tt>>nn_qp>>nn_th>>nn_pc>>nn_ps>>CosThetaMax;
+                        } else {
+                                infile1x>>nn_xx>>nn_tt>>nn_qp>>nn_th>>nn_pc>>CosThetaMax;
+                        }
+                }
+
+                if  (nn_xx!=NEB || nn_tt!=NT || nn_qp!=NQp2
+                        || nn_th!=NTh || nn_pc!=NPhi){
+                        cout<<"WARNING: not the same number of bins in table and in generator"<<endl;
+                        cout<<"->check TreatOptions.h"<<endl;
+                        cout<<"Eb: "<<nn_xx<<" "<<NEB<<endl;
+                        cout<<"-t: "<<nn_tt<<" "<<NT<<endl;
+                        cout<<"Qp: "<<nn_qp<<" "<<NQp2<<endl;
+                        cout<<"Th: "<<nn_th<<" "<<NTh<<endl;
+                        cout<<"Ph: "<<nn_pc<<" "<<NPhi<<endl;
+                } else cout<<"initialization, N bins OK"<<endl;
+                if (targetpoldir==1 || targetpoldir==2){
+                        if (nn_ps != NPhis){
+                                cout<<"WARNING: not the same number of bins in table and in generator"<<endl;
+                                cout<<"phi_s: "<<nn_ps<<" "<<NPhis<<endl;
+                        }
+                }
+                if (beampolartype==1){
+                        if (nn_ps != NPsis){
+                                cout<<"WARNING: not the same number of bins in table and in generator"<<endl;
+                                cout<<"psi_s: "<<nn_ps<<" "<<NPsis<<endl;
+                        }
+                }
+                if (beampolartype==1 && (targetpoldir==1 || targetpoldir==2)){
+                        cout<<"WARNING: transversaly polarized target + linearly polarized beam is not included in TCS event generator"<<endl;
+                        cout<<"photon beam is switch to circular polarization"<<endl;
+                        beampolartype=0;
+                }
+                std::cout<<"Max theta_gamma-gamma = "<<CosThetaMax<<std::endl;
+                float Egamma_bintable[NEB+1], tt_bintable[NT+1], Qp2_bintable[NQp2+1],phi_bintable[NPhi+1],theta_bintable[NTh+1], Phis_bintable[NPhis+1], Psis_bintable[NPsis+1];
+                for (int i=0; i<=NEB; i++){
+                      Egamma_bintable[i] = (float) LinearBins(NEB,EMINI,EMAXI,i);
+                }
+                for (int i=0; i<=NT; i++){
+                      tt_bintable[i] = (float) LinearBins(NT,TMINI,TMAXI,i);
+                }
+                for (int i=0; i<=NQp2; i++){
+                        Qp2_bintable[i] = (float) LinearBins(NQp2,QP2MINI,QP2MAXI,i);
+                }
+                for (int i=0; i<=NPhi; i++){
+                        phi_bintable[i] = (float) LinearBins(NPhi,PHIMINI,PHIMAXI,i);
+                }
+                for (int i=0; i<=NTh; i++){
+                        theta_bintable[i] = (float) LinearBins(NTh,THMINI,THMAXI,i);
+                }
+                if (targetpoldir==1 || targetpoldir==2){
+                     for (int i=0; i<=NPhis; i++){
+                        Phis_bintable[i] = (float) LinearBins(NPhis,PHISMINI,PHISMAXI,i);
+                     }
+                }
+                if (beampolartype==1){
+                     for (int i=0; i<=NPsis; i++){
+                        Psis_bintable[i] = (float) LinearBins(NPsis,PSISMINI,PSISMAXI,i);
+                     }
+                }
+
+                cout<<"\n start: insert values associating cross sections to kinematics"<<endl;
+                cout<<">>> this step can take few minutes"<<endl;
+                ii=0;jj=0;kk=0;ll=0;mm=0;nn=0;
+                while (infile1x.good()){
+                        if (!(infile1x>>em)) break;
+                        if (targetpoldir==1 || targetpoldir==2){
+                                infile1x>> tm>>Qpm>>Thm>>phim>>psim>> res[0]>>res[1]>>res[2]>>res[3]>>res[4]>>res[5];
+                        } else {
+                                if (beampolartype==1){
+                                        infile1x>> tm>>Qpm>>Thm>>phim>>psim>> res[1]>>res[2]>>res[3]>>res[4]>>res[5]>>res[6]>>res[7];
+
+                                } else {
+                                        infile1x>> tm>>Qpm>>Thm>>phim>> res[0]>>res[1]>>res[2]>>res[3]>>res[4]>>res[5];
+                                }
+                        }
+
+                        ii=SetBins(em+0.001,Egamma_bintable,NEB);
+                        jj=SetBins(tm+0.001,tt_bintable,NT);
+                        kk=SetBins(Qpm+0.001,Qp2_bintable,NQp2);
+                        ll=SetBins(Thm+0.001,theta_bintable,NTh);
+                        mm=SetBins(phim+0.001,phi_bintable,NPhi);
+                        if (beampolartype==1) nn=SetBins(psim+0.001,Psis_bintable,NPsis);
+                        if (targetpoldir==1 || targetpoldir==2) nn=SetBins(psim+0.001,Phis_bintable,NPhis);
+
+                        if (targetpoldir==1 || targetpoldir==2 || beampolartype==1){
+
+                            if(!table_provider.insert_tcs_pol( {ii,jj,kk,ll,mm,nn}, {res[0],res[1],res[2],res[3],res[4],res[5], res[6], res[7] })) {
+                                                                std::cout << "[Fill Table TCS pol] FAILED TO INSERT " << Vars_Cin_tcs_pol({ii,jj,kk,ll,mm,nn})
+                                                                << " -> " << Vars_Sef_tcs_pol({res[0],res[1],res[2],res[3],res[4],res[5]}) << std::endl;
+return 5;
+                            } else {
+                                countercheck++;
+                            }
+                        } else {
+                            if(!table_provider.insert_tcs( {ii,jj,kk,ll,mm}, {res[0],res[1],res[2],res[3],res[4],res[5]})) {
+                                                                std::cout << "[Fill Table] FAILED TO INSERT " << Vars_Cin_tcs({ii,jj,kk,ll,mm})
+                                                                << " -> " << Vars_Sef_tcs({res[0],res[1],res[2],res[3],res[4],res[5]}) << std::endl;
+return 5;
+                            } else {
+                                countercheck++;
+                            }
+                        }
+
+                }
+                cout<<"size of table: "<<countercheck<<endl;
+                infile1x.close();
+                return 1;
+
+	} // end tcs
+	 else {
+                printf ("\n process= use lower case call help");
+                cout<<"EXIT with ERROR"<<endl;
+                return 0;
+        }
+
+} // end fill table
+
+
+
+
+
+
+
+
+
+

src/programs/Simulation/gen_tcs_bh/FillTables.h

@@ -0,0 +1,33 @@
+#ifndef FillTables_h
+#define FillTables_h
+#include <iostream>
+#include <fstream>
+#include "Riostream.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <cmath>
+#include "TROOT.h"
+#include "TFile.h"
+#include "TableProvider.h"
+#include "TabUtils.h"
+#include "Options_tcs.h"
+#include "TreatOptions.h"
+#include "Utils.h"
+#include <unistd.h>
+#include "Constants.h"
+#include "genSettings_t.h"
+
+class FillTables {
+
+public:
+int FillTable(int reaction, TableProvider& table_provider);
+
+private:
+
+};
+
+#endif
+
+
+

src/programs/Simulation/gen_tcs_bh/FormFactors.cc

@@ -0,0 +1,61 @@
+#define FormFactors_cxx
+#include "FormFactors.h"
+#include <limits>
+#include <unistd.h>
+using namespace std;
+
+//----------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------
+
+float G_E(int nucleon, float QQ){
+        float result=0;
+        float MagM = 2.7928;
+        if (nucleon==1) result = G_M(nucleon,QQ)/MagM * (1. - 0.1306 * QQ + 0.004174 * pow(QQ,2.) - 0.000752 * pow(QQ,3.));
+        return result;
+
+}
+
+//----------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------
+
+float G_M (int nucleon, float QQ){
+        float result=0;
+        float MagM = 2.7928;
+        if (nucleon==1) result = MagM /( 1. + 0.116 * sqrt(QQ) + 2.874 * QQ
+                        + 0.241 * pow(QQ,1.5) + 1.006 * pow(QQ,2.) + 0.345 * pow(QQ,2.5));
+        return result;
+}
+
+//----------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------
+
+float FormFactors(int FF,int nucleon, float QQ){
+
+        float result=0;
+        float G_E = 0, G_M=0;
+        float MagM = 2.7928;
+
+        if (nucleon==1){
+                G_M = MagM /( 1. + 0.116 * sqrt(QQ) + 2.874 * QQ
+                        + 0.241 * pow(QQ,1.5) + 1.006 * pow(QQ,2.) + 0.345 * pow(QQ,2.5));
+
+                G_E = G_M/MagM * (1. - 0.1306 * QQ + 0.004174 * pow(QQ,2.) - 0.000752 * pow(QQ,3.));
+
+        } else if (nucleon==2){
+                cout<<"not implemented yet"<<endl;
+                G_E=0;
+                G_M=0;
+
+        } else {
+                G_E=0;
+                G_M=0;
+        }
+
+        if (FF==1){
+                result = ( QQ / pow(2. * M_Nucleon, 2.) * G_M + G_E ) /( 1. + QQ / pow(2. * M_Nucleon, 2.) );
+        } else if (FF ==2){
+                result =  ( G_M - G_E ) / ( 1. + QQ / pow(2. * M_Nucleon, 2.) );
+        }
+        return result;
+}
+