From 47a2460591d7f68cad09f5c7f3d5157c6dbb056c Mon Sep 17 00:00:00 2001 From: Seiya Hayakawa Date: Thu, 28 Sep 2023 18:33:21 +0900 Subject: [PATCH] Updated files to apply to generic reactions --- .gitignore | 12 +- Makefile | 68 ++++++- prm_reac.C | 152 +------------- simthm.C | 572 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 637 insertions(+), 167 deletions(-) mode change 100644 => 120000 prm_reac.C create mode 100644 simthm.C diff --git a/.gitignore b/.gitignore index df9533d..50fc17b 100644 --- a/.gitignore +++ b/.gitignore @@ -6,9 +6,9 @@ *.d .DS_Store *~ -simthm_c11li6n14pd -simthm_c11li6n14pd.?*?.C -root/*.root -figs/*.png -figs/*.pdf -figs/*.ps \ No newline at end of file +simthm +simthm.?*?.C +*/root/*.root +*/figs/*.png +*/figs/*.pdf +*/figs/*.ps \ No newline at end of file diff --git a/Makefile b/Makefile index 10ad437..fb59831 100644 --- a/Makefile +++ b/Makefile @@ -1,27 +1,40 @@ # name of .C program -TARGET = simthm_c11li6n14pd +TARGET = simthm DEP = $(TARGET).d SRCS = $(TARGET).C OBJS = $(addsuffix .o, $(basename $(SRCS))) +ROOTLIBDIR = $(shell root-config --libdir) ROOTCFLAGS = $(shell root-config --cflags) -ROOTLIBS = $(shell root-config --libs) CXXFLAGS = $(ROOTCFLAGS) -Wall -GSLLIBS = -lm -CXXLIBS = $(ROOTLIBS) $(GSLLIBS) -lGeom -lMathMore -CC = g++ -F77=gfortran +F77 = gfortran -ENEWZLIB = /usr/local/lib # where libenewzlib.a locates +# Edit below for your circumstances +ifeq ($(shell uname),Linux) + ROOTLIBS = $(shell root-config --libs) + CXXLIBS = $(ROOTLIBS) -lGeom $(GSLLIBS) -lMathMore + CC = clang++-7 + ENEWZLIB = /usr/local/lib #mac <- used + LFLAGSC = -L$(ENEWZLIB) -lenewzlib \ + /usr/lib/gcc/x86_64-linux-gnu/9/libstdc++.a +else ifeq ($(shell uname),Darwin) + ROOTLIBS = -L$(ROOTLIBDIR) -lCore -lImt -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lROOTVecOps -lTree -lTreePlayer -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -lMultiProc -lROOTDataFrame -Wl,-rpath,$(ROOTLIBDIR) -lpthread -lm -ldl + CXXLIBS = $(ROOTLIBS) -lGeom -lstdc++ -lc++ $(GSLLIBS) -lMathMore + CC = clang++ +# CC = g++ + ENEWZLIB = /usr/local/lib + LFLAGSC = -L$(ENEWZLIB) -lenewzlib #/usr/local/lib/gcc/11/libstdc++.a +endif -LFLAGSC = -L$(ENEWZLIB) -lenewzlib /usr/lib/gcc/x86_64-linux-gnu/9/libstdc++.a -# where libstdc++.a locates +# cp_prm: all: dep $(TARGET) $(TARGET): $(OBJS) +# cp -p enewz_prm_${ISOTOPE}.C enewz_prm.C $(F77) $(OBJS) -o $@ $(CXXLIBS) ${LFLAGSC} +# $(F77) $(OBJS) -o $@ $(ROOTLIBS) ${LFLAGSC} $(GSLLIBS) .cc.o: $(CC) -c $< $(CXXFLAGS) @@ -34,4 +47,39 @@ dep: -include $(DEP) clean: - rm -f $(TARGET) $(OBJS) + rm -f $(TARGET) $(OBJS) $(DEP) + +# DEP = $(TARGET).d +# SRCS = $(TARGET).C +# OBJS = $(addsuffix .o, $(basename $(SRCS))) + +# ROOTCFLAGS = $(shell root-config --cflags) +# ROOTLIBS = $(shell root-config --libs) +# CXXFLAGS = $(ROOTCFLAGS) -Wall +# GSLLIBS = -lm +# CXXLIBS = $(ROOTLIBS) $(GSLLIBS) -lGeom -lMathMore +# CC = g++ +# F77=gfortran + +# ENEWZLIB = /usr/local/lib # where libenewzlib.a locates + +# LFLAGSC = -L$(ENEWZLIB) -lenewzlib /usr/lib/gcc/x86_64-linux-gnu/9/libstdc++.a +# # where libstdc++.a locates + +# all: dep $(TARGET) + +# $(TARGET): $(OBJS) +# $(F77) $(OBJS) -o $@ $(CXXLIBS) ${LFLAGSC} + +# .cc.o: +# $(CC) -c $< $(CXXFLAGS) + +# .f.o: +# ${F77} -c $< + +# dep: +# $(CC) -MM $(SRCS) > $(DEP) $(CXXFLAGS) +# -include $(DEP) + +# clean: +# rm -f $(TARGET) $(OBJS) $(DEP) diff --git a/prm_reac.C b/prm_reac.C deleted file mode 100644 index 9ccca43..0000000 --- a/prm_reac.C +++ /dev/null @@ -1,151 +0,0 @@ -/* Common constants */ -const Double_t u = 931.49432/1000.; // GeV -//const Double_t c = 299.7924; // speed of light, mm/ns -const Double_t pi = 4.*atan(1.); -const Double_t d2r = pi/180.; -const Double_t dummy = -1000.; -////////////////////////// - -/* Event information */ -//Long64_t ntot = 2.e9; // number of events -//Long64_t ntot = 1.e8; // number of events -const Long64_t ntot = 2.e5; // number of events -Char_t label[50] = "c11li6n14pd"; // -//Char_t simnum[10] = "sim001"; // No beam spread, detector resolutions (in energy, space) -//Char_t simnum[10] = "sim002"; // only beam (energy, spacial, PPAC resolution) spread -//Char_t simnum[10] = "sim003"; // only SSD energy resolution */ -// Char_t simnum[10] = "sim004"; // only SSD strip width */ -/* Char_t simnum[10] = "sim005"; // including beam spread, SSD resolutions */ -//Char_t simnum[10] = "sim006"; // Only target thickness and beam energy spread are off -//Char_t simnum[10] = "sim007"; // same as sim005 but the target thickness is x4 -//Char_t simnum[10] = "sim008"; // same as sim005 but different forward detector setting -//Char_t simnum[10] = "sim009"; // same as sim005 but W no. of strips -Char_t simnum[10] = "sim010"; // better beam energy spread, W no. of strips -////////////////////////// - - -/* particles information */ -Int_t z1 = 1; -float mh= 1.00782503; -float mn= 1.00866492; -float md= 2.01410178; -//float mt= 3.01604928; -Int_t z2 = 2; -float ma= 4.00260325; -Int_t z3 = 3; -float m6li= 6.01512279; -Int_t z6 = 6; -float m11c = 11.01143361; -Int_t z7 = 7; -float m14n = 14.00307400; -const Int_t Nexci = 1; -Double_t exci[Nexci] = {0.}; // in MeV, 14N -//Double_t exci[Nexci] = {0., 2.312798, 3.94810, 4.9151}; // in MeV, 14N -Int_t zprj = z6; -float mprj = m11c; -Int_t ztrg = z3; -float mtrg = m6li; -Int_t zres = z7; -float mprt = ma; -float mres[Nexci], q2[Nexci]; -float mres_gs = m14n; -Int_t zrec = z1; -float mrec = mh; -Int_t zspc = z1; -float mspc = md; -Int_t atom_num[3] = {zres, zrec, zspc} ; // atomic numbers -////////////////////////// - - -/* parameters for eoldz */ -float Eaftpu, Eaftpu_beam; -float Eaft, Eaft_beam; -//Char_t al[34] = "al"; -//Char_t ch2[34] = "POLYETHYLENE"; -//Char_t cd2[34] = "cd2"; -Char_t lif[34] = "lif"; -//Char_t he[34] = "he"; -Char_t si[34] = "si"; -Char_t c[34] = "c"; -Int_t unit_pressure = 1; // Torr -// float pressure0 = 0.; //dummy -float pressure = 0.; //dummy -// float dpressure = 0.; // sigma -//float dpressure = 1.; // sigma -//float pressure, temperature; // sigma -Int_t unit_temperature = 1; // K, for srho_gas -//float temperature0 = 293.; -//float temperature0 = 299.; -float temperature = 299.; -//float dtemperature = 0.; // sigma -//float dtemperature = 1.; // sigma -Int_t unit_mm = 1; // in mm -Int_t unit_mgpcm2 = 2; // in mg/cm^2 -//float thk_tar = 0.; // Target thickness in mg/cm^2 // sim001 sim003 sim004 sim006 -//float thk_tar = 0.15; // Target thickness in mg/cm^2 // sim002 sim005 -float thk_tar = 0.1; // Target thickness in mg/cm^2 // sim010 -float thk_bak = 0.1; // Backing thickness in mg/cm^2 // sim010 -//float thk_tar = 0.6; // Target thickness in mg/cm^2 // sim007 -//float dEbin = 0.05; // energy step for effective thickness calculation -float Ebfrpu, Ebfrpu_beam; // Beam energy E in MeV/u after energy loss in matter -float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter -////////////////////////// - - -/* Beam condition */ -// 24.4 MeV --(2xPPACs 19 um Mylar)--(1/2 LiF target 75 ug/cm^2)-- 11 MeV -// 24.4 * 2 *10/16/100 * 0.68(1 sigma) = 0.21 MeV -// Monte Carlo reaction depth in target -// |---x---| Reference: middle -> E.beam_mid = (eprj_mid +/- deprj)*1000. -// |---|-x-| Practice: -1/2 thk_tar -- +1/2 thk_tar -> E.beam_MC -// -1------1 -Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle in MeV -//Double_t dEbeam = 0.; // Beam energy spread (GeV) // sim001 sim003 sim004 -//Double_t dEbeam = 0.21; // Beam energy spread (MeV) -Double_t dEbeam = 0.5; // Beam energy spread (MeV) // sim010 -//Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 -Double_t depe2 = 0.006; // Relative energy resolution DE/E @5MeV sim003 sim005 sim006 -Double_t xtar0 = 0.; // mm -Double_t ytar0 = 0.; // mm -/* Double_t xtar0_sgm = 0.; // mm sim001 sim003 sim004 */ -/* Double_t ytar0_sgm = 0.; // mm sim001 sim003 sim004 */ -Double_t xtar0_sgm = 5.; // sim002 sim005 sim006 -Double_t ytar0_sgm = 10.; // sim002 sim005 sim006 -// -Double_t za = - 600.; // mm, PPAC a z position -Double_t zb = za + 500.; // mm, PPAC b z position -/* Double_t sgm_ppac_a = 0.; // mm, PPAC a intrinsic resolution // sim001 sim003 sim004 */ -/* Double_t sgm_ppac_b = 0.; // mm, PPAC b intrinsic resolution // sim001 sim003 sim004 */ -Double_t sgm_ppac_a = 0.5; // mm, PPAC a intrinsic resolution sim002 sim005 sim006 -Double_t sgm_ppac_b = 0.5; // mm, PPAC b intrinsic resolution sim002 sim005 sim006 -/* ////////////////////////// */ - - -/* Detector arrangement */ -const Int_t ndet = 6; -//Double_t thdet[ndet] = {-56., -34., -12., 12., 34., 56.}; // in deg -//Double_t thdet[ndet] = {-60., -40., -12., 12., 40., 60.}; // in deg sim001-008 -Double_t thdet[ndet] = {-60., -40., -11., 11., 40., 60.}; // in deg sim010 -//Double_t thdet[ndet] = {-60., -40., -9., 9., 40., 60.}; // in deg sim008 -//Double_t ldet[ndet] = {200., 200., 200., 200., 200., 200.}; //sim001-007 -//Double_t ldet[ndet] = {200., 200., 300., 300., 200., 200.}; // sim008 -Double_t ldet[ndet] = {200., 200., 200., 200., 200., 200.}; //sim010 -Double_t xdet_off[ndet] = {0., 0., 0., 0., 0., 0.}; -Double_t zdet_off[ndet] = {0., 0., 0., 0., 0., 0.}; -// float int_det = 10.; // interval between DE and E layers (mm) -// float de_thick = 0.02; // DE thickness (mm) -Double_t psd_half[ndet] = {25., 25., 25., 25., 25., 25.}; //sim001-007 -//Double_t psd_half[ndet] = {25., 25., 20., 20., 25., 25.}; //sim008 -//Int_t n_strip[ndet] = {0}; // sim001 sim002 sim003 -//Int_t n_strip[ndet] = {16, 16, 16, 16, 16, 16}; // sim004 sim005 sim006 sim007 -//Int_t n_strip[ndet] = {16, 16, 40, 40, 16, 16}; // sim008 -Int_t n_strip[ndet] = {32, 32, 32, 32, 32, 32}; // sim009 -Double_t strip_half[ndet]; -//Double_t dth[ndet]; // 03-04-2018 -float r_sphere = 250.; // radius of the TGeoVolume in mm -//float r_sphere = 320.; // radius of the TGeoVolume in mm -////////////////////////// - -//Double_t dt = 1.; // TOF resolution in nsec. - - diff --git a/prm_reac.C b/prm_reac.C new file mode 120000 index 0000000..b16c767 --- /dev/null +++ b/prm_reac.C @@ -0,0 +1 @@ +n14dc11an/prm/prm_reac.n14dc11an_sim005.C \ No newline at end of file diff --git a/simthm.C b/simthm.C new file mode 100644 index 0000000..91693ac --- /dev/null +++ b/simthm.C @@ -0,0 +1,572 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TGeoPatternFinder.h" +#include "TSystem.h" +#include +#include +#include +#include +#include +#include "TApplication.h" +#include "TF1.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TH1.h" +#include "TLegend.h" +#include "TLegendEntry.h" +#include "TColor.h" +#include +#include +#include // for the seed of the dandom number +#include + +//#include "Math/Interpolator.h" +//#include "Math/Polynomial.h" +//#include "Math/InterpolationTypes.h" +//#include "Math/Types.h" + +#include "prm_reac.C" // parameter file + +using namespace std; + +extern "C" { + void enewzsub_etot_(int *z1, float *m1, float *bfr_ene, char matter1[33], + int *unit_pressure, float *pressure, float *temperature, + int *unit_thick, float *thick1, float *aft_ene); + void eoldzsub_etot_(int *z1, float *m1, float *aft_ene, char matter1[33], + int *unit_pressure, float *pressure, float *temperature, + int *unit_thick, float *thick1, float *bfr_ene); + /* void rangezsub_etot_(int *z1, float *m1, float *bfr_ene, char matter1[33], */ + /* float *range); */ + /* void srhogas_(float *pressure, float *temperature, int *unit_pressure, */ + /* int *unit_temperature, float *mol_weight, */ + /* float *density); */ +} + +Double_t sgm_ppac_ab(Double_t z){ + Double_t sgm = sqrt(pow((zb-z)/(zb-za),2)*sgm_ppac_a + + pow((za-z)/(zb-za),2)*sgm_ppac_b); + return sgm; +}; + +Double_t theta_cm(Double_t mprj, Double_t mtrg, TLorentzVector LV_beam, TLorentzVector LV3, TVector3 vel[3]){ + TVector3 vprj = pow(u*mprj,-1)/LV_beam.Gamma() *(LV_beam.Vect()); // velocity + TLorentzVector ptran = - LV3; // 4-momentum of the transferred particle + TVector3 vtran = 1./sqrt(pow(mtrg*u - LV3.Energy(),2) + LV3.Mag2()) + /ptran.Gamma() *ptran.Vect(); // velocity, See Slaus et al + TVector3 ki = vprj - vtran; // target breakup 13-11-2017 + TVector3 kf = vel[0] - vel[1]; // Final relative velocity (1<->2); vel[0] is beamlike (heavy ion) + return (ki.Angle(kf))/d2r; +} + +Double_t mom_spec(TVector3 mom){ + Double_t scp; + if( cos(mom.Phi())>0 ) // x>0 direction + scp=1; + else if ( cos(mom.Phi())<0 ) // x<0 direction + scp=-1.; + else + scp=0.; + return mom.Mag()*scp*1000.; // in MeV/c target breakup +} + +Double_t strip_position(Double_t x, Int_t n_strip, Double_t psd_half){ + Double_t x_strip = x; + Double_t w_strip = psd_half*2./double(n_strip); // strip width + Double_t x0 = - psd_half; + Double_t x1 = x0 + w_strip; + for(Int_t i=0; i x0 && x <= x1){ + x_strip = (x0 + x1)/2.; + break; + } + x0 = x1; + x1 += w_strip; + } + return x_strip; +} + +int main() { + TApplication theApp("App", NULL, NULL); + + // For enewz + for(Int_t ichar=0;ichar<33;ichar++){ + /* if(he[ichar]=='\0'){ */ + /* he[ichar]=' '; */ + /* he[ichar+1]='\0'; */ + /* } */ + /* if(si[ichar]=='\0'){ */ + /* si[ichar]=' '; */ + /* si[ichar+1]='\0'; */ + /* } */ + /* if(cd2[ichar]=='\0'){ */ + /* cd2[ichar]=' '; */ + /* cd2[ichar+1]='\0'; */ + /* } */ + /* if(lif[ichar]=='\0'){ */ + /* lif[ichar]=' '; */ + /* lif[ichar+1]='\0'; */ + /* } */ + if(mat_tar[ichar]=='\0'){ + mat_tar[ichar]=' '; + mat_tar[ichar+1]='\0'; + } + if(mat_bak[ichar]=='\0'){ + mat_bak[ichar]=' '; + mat_bak[ichar+1]='\0'; + } + /* if(c[ichar]=='\0'){ */ + /* c[ichar]=' '; */ + /* c[ichar+1]='\0'; */ + /* } */ + } + + ////////////////////////////////////// + // Geometrical definitions + gStyle->SetPalette(1); + gSystem->Load("libPhysics.so"); + gSystem->Load("libGeom.so"); + TGeoManager *sc = new TGeoManager("scattcham","Scattering Chamber"); + TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum",0,0,0); + TGeoMaterial *matSi = new TGeoMaterial("Si",28.086,14,2.321); + TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1,matVacuum); + TGeoMedium *Si = new TGeoMedium("Root Material",2,matSi); + TGeoVolume *top = sc->MakeSphere("TOP",Vacuum,0.,r_sphere,0.,180.,0.,360.); //mm + sc->SetTopVolume(top); + TGeoVolume *psd[ndet]; + // (X,Y,Z,rotation angle theta originating at (X,Y,Z) ) + // tan(theta) = X/Z + TGeoCombiTrans *posdet[ndet]; + //Char_t rotname[10]; + TVector3 vdet[ndet], voff[ndet]; + + // Detectors + for(Int_t idet=0; idetMakeBox("PSD",Si,psd_half[idet],psd_half[idet],0.5); // in mm + vdet[idet].SetXYZ( ldet[idet]*sin(d2r*thdet[idet])+xdet_off[idet], + 0,ldet[idet]*cos(d2r*thdet[idet])+zdet_off[idet] ); + //sprintf(rotname,"rot%d",idet); + posdet[idet] = new + // Rotate detectors + TGeoCombiTrans(vdet[idet].X(),vdet[idet].Y(),vdet[idet].Z(), + new TGeoRotation(Form("rot%d",idet),90,thdet[idet],0)); + // Put detectors at coordinates + // TGeoCombiTrans(xdet[idet],0.,zdet[idet], + // new TGeoRotation(Form("rot%d",idet),90,rotdet[idet],0)); + top->AddNode(psd[idet],idet,posdet[idet]); + voff[idet].SetXYZ(xdet_off[idet],0.,zdet_off[idet]); + } + + gGeoManager->CloseGeometry(); + top->SetLineColor(kRed); + gGeoManager->SetTopVisible(); + + TCanvas *c0 = new TCanvas("c0","c0",1450,0,800,800); + top->Draw(); + TView3D *view = (TView3D*)TView::CreateView(1); + view->SetRange(-r_sphere,-r_sphere,-r_sphere,r_sphere,r_sphere,r_sphere); + view->Front(); + view->ShowAxis(); + c0->Update(); + gSystem->ProcessEvents(); + gSystem->Exec(Form("mkdir -vp %s/figs",dir)); // Creates figs/ directory if there is not + c0->Print(Form("%s/figs/det_arrange_%s_%s.png",dir,label,simnum)); + + ////////////////////////////////////// + // Output tree file + gSystem->Exec(Form("mkdir -vp %s/root",dir)); // Creates root/ directory if there is not + //Char_t fname[200]; + TString fname; + //sprintf(fname,"%s/root/%s_%s.root",dir,label,simnum); + fname = Form("%s/root/%s_%s.root",dir,label,simnum); + cout << "Output file name: " << fname.Data() << endl; + cout << "OK to run?" << endl; getchar(); + + // to store the parameter file + gSystem->Exec(Form("mkdir -vp %s/prm",dir)); // Creates prm/ directory if there is not + gROOT->ProcessLine(Form(".! cp -pv prm_reac.C %s/prm/prm_reac.%s_%s.C", + dir,label,simnum)); + + TFile *fl2 = new TFile(fname.Data(),"RECREATE"); + TTree* sim = new TTree("sim",""); + // variables for tree + // E.beam_mid: fixed beam energy at the middle of the target + // E.beam_mid_MC: same as _mid but deviated according to energy spread + // E.beam_dec: beam energy at arbitrary target depth + typedef struct{ + Double_t beam_mid, beam_mid_MC, beam_dec; + Double_t lab[3], lab_det[2], lab_det_MC[2], lab_MC[3]; + Double_t rel[3], rel_MC[3], cm, cm_MC; + Double_t Q3_mid, Q3_dec; // 3-body Q-value + } ENERGY; + ENERGY E; + sim->Branch("E",&E,"beam_mid/D:beam_mid_MC:beam_dec:" + "lab[3]:lab_det[2]:lab_det_MC[2]:lab_MC[3]:" + "rel[3]:rel_MC[3]:cm:cm_MC:Q3_mid:Q3_dec"); + + typedef struct{ + Double_t lab[3], lab_MC[3], cm, cm_MC, cm_dev, cm_MC_nr; + } THETA; + THETA th; + sim->Branch("th",&th,"lab[3]/D:lab_MC[3]:cm:cm_MC:cm_dev:cm_MC_nr"); + + typedef struct { + Double_t lab[3], lab_MC[3]; + } PHI; + PHI phi; + sim->Branch("phi",&phi,"lab[3]/D:lab_MC[3]"); + + typedef struct { + Double_t spc, spc_MC; + } MOMENTUM; + MOMENTUM p; + sim->Branch("p",&p,"spc/D:spc_MC"); + + typedef struct { + Double_t xtar, ytar, rtar; + Double_t ztar; + Double_t xtar_MC, ytar_MC, rtar_MC; + Double_t xhit[2][ndet], yhit[2][ndet], xhit_MC[2][ndet], yhit_MC[2][ndet]; + } POSITION; + POSITION pos; + sim->Branch("pos",&pos, + Form("xtar/D:ytar:rtar:ztar:" + "xtar_MC:ytar_MC:rtar_MC:" + "xhit[2][%d]:yhit[2][%d]:" + "xhit_MC[2][%d]:yhit_MC[2][%d]", + ndet,ndet,ndet,ndet) + ); + + typedef struct{ + Int_t idet[2]; // for particle 1&2 + Bool_t any; + Bool_t both; + } HIT; + HIT hit; + sim->Branch("hit",&hit,"idet[2]/I:any/O:both/O"); + + typedef struct { + Double_t evt; //, efth, hf, yld; + } WEIGHT; + WEIGHT w; + sim->Branch("w",&w,"evt/D"); //:efth:hf:yld"); + // + Int_t iex; + sim->Branch("iex",&iex,"iex/I"); + for(Int_t iexci=0; iexciRndm()*Nexci; // 0-1 * Nexci + Double_t mass[3] = { u*mres[iex], u*mrec, u*mspc} ; // masses of 14N, p, d + + // Definition of beam energy + // Monte Carlo reaction depth in target + // |---x---| _mid Reference: middle -> E.beam_mid + // |---x---| _mid_MC with spread: middle -> E.beam_mid_MC + // |---|-x-| _dec Practice: -1/2 thk_tar -- +1/2 thk_tar -> E.beam_dec + // -1------1 + /* E.beam_mid = eprj_mid*1000.; // defined energy at middle */ + /* E.beam_mid_MC = ran->Gaus(eprj_mid,deprj)*1000.; // energy at middle with initial beam spread */ + E.beam_mid = Ebeam_ref_mid; // defined energy at middle + E.beam_mid_MC = ran->Gaus(Ebeam_ref_mid,dEbeam); // energy at middle with initial beam spread + float fthk_mid_dec = thk_tar/2.*ran->Uniform(-1.,1.); // distance from the middle + Ebfr = E.beam_mid_MC; + enewzsub_etot_(&zprj, &mprj, &Ebfr, mat_tar, &unit_pressure, &pressure, + &temperature, &unit_mgpcm2, &fthk_mid_dec, &Eaft); + E.beam_dec = Eaft; + + // Beam hitting position + pos.xtar = ran->Gaus(xtar0,xtar0_sgm); + pos.ytar = ran->Gaus(ytar0,ytar0_sgm); + pos.ztar = 0.; + pos.rtar = sqrt(pos.xtar*pos.xtar + pos.ytar*pos.ytar); + //TVector3 vtar(xtar,ytar,0); // vector from origin to target pos. + TVector3 vdec(pos.xtar,pos.ytar,pos.ztar); // vector to the reaction point. + // pos.xtar -> pos.xtar_MC + // pos.ytar -> pos.ytar_MC + // -> pos.rtar_MC + Double_t sgm_track = sgm_ppac_ab(pos.ztar); + pos.xtar_MC = ran->Gaus(pos.xtar,sgm_track); + pos.ytar_MC = ran->Gaus(pos.ytar,sgm_track); + pos.rtar_MC = sqrt(pos.xtar_MC*pos.xtar_MC + pos.ytar_MC*pos.ytar_MC); + TVector3 vmid_MC(pos.xtar_MC,pos.ytar_MC,pos.ztar); // ztar for solid target + + // Beam, compound, decay + // E^2 = m^2*c^4 + c^2*p^2 + // E = K + m*c^2 + // c*p = sqrt(K^2 + 2*m*c^2*K) + //TLorentzVector beam_dec(0.0, 0.0, sqrt(2.*u*mprj*E.beam_dec/1000.), u*mprj+E.beam_dec/1000.); // px, py, pz, e + //TLorentzVector beam_mid(0.0, 0.0, sqrt(2.*u*mprj*E.beam_mid/1000.), u*mprj+E.beam_mid/1000.); // px, py, pz, e + TLorentzVector LV_beam_dec(0.0, 0.0, + sqrt(pow(E.beam_dec/1000.,2) + 2.*E.beam_dec*u*mprj/1000.), + u*mprj+E.beam_dec/1000.); // px, py, pz, e + TLorentzVector LV_beam_mid(0.0, 0.0, + sqrt(pow(E.beam_mid/1000.,2) + 2.*E.beam_mid*u*mprj/1000.), + u*mprj+E.beam_mid/1000.); // px, py, pz, e + TLorentzVector LV_W_dec = LV_beam_dec + LV_target; // coumpound of beam + target + TLorentzVector LV_W_mid = LV_beam_mid + LV_target; // coumpound of beam + target + event.SetDecay(LV_W_dec, 3, mass); // W: beam + target -> 3-body decay: mres+mrec+m3 + w.evt = event.Generate(); + + // Decay kinematics + for(Int_t iprt=0; iprt<3; iprt++){ + LV[iprt] = *event.GetDecay(iprt); // lorentz vector + mom[iprt] = LV[iprt].Vect(); // momenta + vel[iprt] = pow(mass[iprt],-1)/LV[iprt].Gamma() *mom[iprt]; // velocity + th.lab[iprt] = (LV[iprt].Theta())/d2r; + phi.lab[iprt] = (LV[iprt].Phi())/d2r; + E.lab[iprt] = (LV[iprt].Energy()-mass[iprt])*1000.; // MeV + } + for(Int_t iprt=0; iprt<3; iprt++){ // After LV[iprt] are all defined... + LVsum[iprt] = LV[iprt] + LV[(iprt+1)%3]; // (iprt+1)%3 = 1, 2, 0 + E.rel[iprt] = (LVsum[iprt].M() - (mass[iprt] + mass[(iprt+1)%3]))*1000.; // MeV + } + E.cm = E.rel[0] - q2[iex]; // MeV + th.cm = theta_cm(mprj, mtrg, LV_beam_dec, LV[2], vel); + p.spc = mom_spec(mom[2]); // in MeV/c, tartget breakup + E.Q3_mid = E.lab[0] + E.lab[1] + E.lab[2] - E.beam_mid; + + // Tracking decaying particles + /// Particle 1: residual heavy ion + /// Particle 2: recoil light particle + /// Particle 3: spectator + Double_t xx, yy, zz; + Double_t snext; + Bool_t istate[2] = {false, false}; + Int_t idnode[2]; + const Char_t *path; + Char_t sdet[2]; + for(Int_t iprt = 0; iprt<2; iprt++){ + pdir[iprt] = mom[iprt].Unit(); // Unit vector of the particle motion + xx = pdir[iprt](0); // Each element of the vector + yy = pdir[iprt](1); + zz = pdir[iprt](2); + gGeoManager->SetCurrentPoint(pos.xtar,pos.ytar,pos.ztar); // Reaction point + gGeoManager->SetCurrentDirection(xx,yy,zz); // Reaction direction + gGeoManager->GetCurrentNode(); + gGeoManager->FindNode(); + gGeoManager->FindNextBoundary(r_sphere); // finds distance to next boundary -> fStep + // for solid target, r_sphere, for gas target, r_sphere*2. + snext = gGeoManager->GetStep(); // Dist. reaction point <-> detector hit point + gGeoManager->Step(); + istate[iprt] = gGeoManager->IsStepEntering(); + idnode[iprt] = gGeoManager->GetCurrentNodeId(); + //Double_t safety2 = gGeoManager->GetSafeDistance(); + path = gGeoManager->GetPath(); // /TOP_1/PSD_3 etc. + strncpy(sdet, path+11, 1); // extract the 12th character from pathq + sdet[1] = '\0'; + hit.idet[iprt] = atoi(sdet); // Number of the hitting + vdec_det[iprt].SetXYZ(snext*xx,snext*yy,snext*zz); // Vector reaction point <-> detector hit point + } + + // If particles hit any detectors + if(istate[0] || istate[1]){ + hit.any = true; + /* cout << "any" << endl; */ + } + + // If both particle 1 & 2 hit different detectors respectively + if(istate[0] && istate[1] && idnode[0]!=idnode[1]){ + hit.both = true; + ////////////////////// + // Resolution Study // + ////////////////////// + for(Int_t iprt = 0; iprt<2; iprt++){ // Particle 1 & 2 + // Integerization: no. of the hit detector // Warning: when no hit, it returns 0! + vhit[iprt][hit.idet[iprt]] + = vdec_det[iprt] + vdec - vdet[hit.idet[iprt]]; + vhit[iprt][hit.idet[iprt]].RotateY(-d2r*thdet[hit.idet[iprt]]); + pos.xhit[iprt][hit.idet[iprt]] = vhit[iprt][hit.idet[iprt]].X(); + pos.yhit[iprt][hit.idet[iprt]] = vhit[iprt][hit.idet[iprt]].Y(); + + //////////////////////////////////////// + // Energy loss by the target material // + //////////////////////////////////////// + // E.lab -> E.lab_det + float fthk_dec_end = (thk_tar/2. - fthk_mid_dec)/cos(th.lab[iprt]*d2r); // Dist. reaction point <-> detector hit + Ebfr = E.lab[iprt]; // MeV + float fmass = mass[iprt]/u; + enewzsub_etot_(&atom_num[iprt], &fmass, &Ebfr, mat_tar, &unit_pressure, &pressure, + &temperature, &unit_mgpcm2, &fthk_dec_end, &Eaft); + float fthk_bak = thk_bak/cos(th.lab[iprt]*d2r); + Ebfr = Eaft; + enewzsub_etot_(&atom_num[iprt], &fmass, &Ebfr, mat_bak, &unit_pressure, &pressure, + &temperature, &unit_mgpcm2, &fthk_bak, &Eaft); + E.lab_det[iprt] = Eaft; // MeV + + // E.lab_det[0,1] -> E.lab_det_MC[0,1] + // pos.xhit -> pos.xhit_MC + /* pos.xhit_MC[iprt][hit.idet[iprt]] */ + /* = ran->Uniform(pos.xhit[iprt][hit.idet[iprt]]-strip_half[hit.idet[iprt]], */ + /* pos.xhit[iprt][hit.idet[iprt]]+strip_half[hit.idet[iprt]]); // for strip */ + pos.xhit_MC[iprt][hit.idet[iprt]] + = strip_position( pos.xhit[iprt][hit.idet[iprt]], + n_strip[hit.idet[iprt]], + psd_half[hit.idet[iprt]] ); + // pos.yhit -> pos.yhit_MC + /* pos.yhit_MC[iprt][hit.idet[iprt]] */ + /* = ran->Uniform(pos.yhit[iprt][hit.idet[iprt]]-strip_half[hit.idet[iprt]], */ + /* pos.yhit[iprt][hit.idet[iprt]]+strip_half[hit.idet[iprt]]); // for strip */ + pos.yhit_MC[iprt][hit.idet[iprt]] + = strip_position( pos.yhit[iprt][ hit.idet[iprt] ], + n_strip[ hit.idet[iprt] ], + psd_half[ hit.idet[iprt] ] ); + vhit_MC[iprt].SetXYZ( pos.xhit_MC[iprt][hit.idet[iprt]], + pos.yhit_MC[iprt][hit.idet[iprt]], 0 ); //vector on det. surface + vhit_MC[iprt].RotateY(d2r*thdet[hit.idet[iprt]]); // Rotate from z axis to det. angle + E.lab_det_MC[iprt] = ran->Gaus(E.lab_det[iprt], // based on resolution at 5 MeV + sqrt(5.*E.lab_det[iprt])*depe2); // MeV + // vectors from the mid target to the detection point with resolution + vmid_det_MC[iprt] = vhit_MC[iprt] + vdet[hit.idet[iprt]] - vmid_MC; + // th.lab_MC + th.lab_MC[iprt] = (vmid_det_MC[iprt].Theta())/d2r; + // phi.lab_MC + phi.lab_MC[iprt] = (vmid_det_MC[iprt].Phi())/d2r; + Eaft = E.lab_det_MC[iprt]; + float fthk_bak_MC = thk_bak/cos(th.lab_MC[iprt]*d2r); + fmass = mass[iprt]/u; + eoldzsub_etot_(&atom_num[iprt], &fmass, &Eaft, mat_bak, &unit_pressure, &pressure, + &temperature, &unit_mgpcm2, &fthk_bak_MC, &Ebfr); + Eaft = Ebfr; + float fthk_mid_end_MC = (thk_tar/2.)/cos(th.lab_MC[iprt]*d2r); + eoldzsub_etot_(&atom_num[iprt], &fmass, &Eaft, mat_tar, &unit_pressure, &pressure, + &temperature, &unit_mgpcm2, &fthk_mid_end_MC, &Ebfr); + // -> E.lab_MC[0,1] + E.lab_MC[iprt] = Ebfr; // MeV + + // Definition of TLorentzVector // -> LV_MC[1,0] + LV_MC[iprt].SetE(mass[iprt]+E.lab_MC[iprt]/1000.); + //LV_MC[iprt].SetVect(vmid_det_MC[iprt].Unit()*sqrt(2.*mass[iprt]*E.lab_MC[iprt]/1000.)); + LV_MC[iprt].SetVect(vmid_det_MC[iprt].Unit() + *sqrt(pow(E.lab_MC[iprt]/1000.,2) + 2.*mass[iprt]*E.lab_MC[iprt]/1000.)); + mom_MC[iprt] = LV_MC[iprt].Vect(); + vel_MC[iprt] = pow(mass[iprt],-1)/LV_MC[iprt].Gamma() *mom_MC[iprt]; // velocity + } // for iprt < 2 + + // Particle 3 + LV_MC[2] = LV_W_mid - LV_MC[0] - LV_MC[1]; + mom_MC[2] = LV_MC[2].Vect(); + vel_MC[2] = pow(mass[2],-1)/LV[2].Gamma() *mom[2]; // velocity + // TVector3 pp3vec_MC = beam_mid.Vect() - mom_MC[0] - mom_MC[1]; // momentum conservation + // E.lab_MC[2] = pp3vec_MC.Mag2()/2./mass[2]*1000.; + E.lab_MC[2] = (LV_MC[2].E() - mass[2])*1000.; + th.lab_MC[2] = (LV_MC[2].Theta())/d2r; + phi.lab_MC[2] = (LV_MC[2].Phi())/d2r; + + // Relative energy, E_cm, theta_cm, Q value + for(Int_t iprt=0; iprt<3; iprt++){ // After LV[iprt] are all defined... + LVsum_MC[iprt] = LV_MC[iprt] + LV_MC[(iprt+1)%3]; // (iprt+1)%3 = 1, 2, 0 + E.rel_MC[iprt] = (LVsum_MC[iprt].M() - (mass[iprt] + mass[(iprt+1)%3]))*1000.; // MeV + } + E.cm_MC = E.rel_MC[0] - q2[iex]; // MeV + th.cm_MC = theta_cm(mprj, mtrg, LV_beam_mid, LV_MC[2], vel_MC); + p.spc_MC = mom_spec(mom_MC[2]); // in MeV/c, tartget breakup + E.Q3_dec = E.lab_MC[0] + E.lab_MC[1] + E.lab_MC[2] - E.beam_dec; + + geve++; + } // if istate2 + + //if(rnum<=eff_thk(ps)){ + sim->Fill(); + //geve++; + //} + + } // for ievt + fl2->Write(); + //fl2->Close(); + cout << "Number of coincidence events inside the detectors =" << geve << endl; + cout << "Total number of events =" << ntot << endl; + cout << "Done!" << endl; + //return 0; + + //theApp.Run(); + getchar(); + + return 0; + +} -- GitLab