diff --git a/prm/prm_reac.c11li6n14pd_sim001.C b/prm/prm_reac.c11li6n14pd_sim001.C index 29baf9c0d3db0cdfff421b37a265134b54f40086..ccd9abb31b5349350a026db45ebba3de0bed14eb 100644 --- a/prm/prm_reac.c11li6n14pd_sim001.C +++ b/prm/prm_reac.c11li6n14pd_sim001.C @@ -92,9 +92,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance -//Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle (MeV) +Double_t dEbeam = 0.; // Beam energy spread (MeV) // sim001 +//Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 Double_t depe2 = 0.; // Relative proton energy resolution DE/E // sim001, sim002 //Double_t depe2 = 0.01; // Relative proton energy resolution DE/E Double_t xtar0 = 0.; // mm diff --git a/prm/prm_reac.c11li6n14pd_sim002.C b/prm/prm_reac.c11li6n14pd_sim002.C index 29baf9c0d3db0cdfff421b37a265134b54f40086..0581dcd5ce94d68b4e61acece074690a9982c945 100644 --- a/prm/prm_reac.c11li6n14pd_sim002.C +++ b/prm/prm_reac.c11li6n14pd_sim002.C @@ -92,9 +92,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle //Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 +Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 Double_t depe2 = 0.; // Relative proton energy resolution DE/E // sim001, sim002 //Double_t depe2 = 0.01; // Relative proton energy resolution DE/E Double_t xtar0 = 0.; // mm diff --git a/prm/prm_reac.c11li6n14pd_sim003.C b/prm/prm_reac.c11li6n14pd_sim003.C index 025e71103b8094e4cf0b2007dc8d26f06cc89555..7f835bc6d6edc683ac88c17885d1c30f192ecbcd 100644 --- a/prm/prm_reac.c11li6n14pd_sim003.C +++ b/prm/prm_reac.c11li6n14pd_sim003.C @@ -92,9 +92,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance -Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 -//Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle +Double_t dEbeam = 0.; // Beam energy spread (MeV) // sim001 sim003 +//Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 //Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 Double_t xtar0 = 0.; // mm diff --git a/prm/prm_reac.c11li6n14pd_sim004.C b/prm/prm_reac.c11li6n14pd_sim004.C index 9464b88d1a08f2c20e9ca5faa941cb270eb692be..f7b15f6520895cc052fa16b4ad8c7de935a2c38e 100644 --- a/prm/prm_reac.c11li6n14pd_sim004.C +++ b/prm/prm_reac.c11li6n14pd_sim004.C @@ -92,8 +92,8 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance -Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle +Double_t dEbeam = 0.; // Beam energy spread (MeV) // sim001 sim003 sim004 //Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 //Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 diff --git a/prm/prm_reac.c11li6n14pd_sim005.C b/prm/prm_reac.c11li6n14pd_sim005.C index b33f3407356563f83082d6f2742ee5379a2fd8a0..397f7762fecd72b453cc1f892c0bbbdae1fece24 100644 --- a/prm/prm_reac.c11li6n14pd_sim005.C +++ b/prm/prm_reac.c11li6n14pd_sim005.C @@ -92,9 +92,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle //Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 sim005 +Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 sim005 //Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 sim005 Double_t xtar0 = 0.; // mm diff --git a/prm/prm_reac.c11li6n14pd_sim006.C b/prm/prm_reac.c11li6n14pd_sim006.C index 1068036ff4accb04825932dc11abd65d95e01316..eae783b0ce7f8fd722f233c9277e01c28ec21030 100644 --- a/prm/prm_reac.c11li6n14pd_sim006.C +++ b/prm/prm_reac.c11li6n14pd_sim006.C @@ -92,8 +92,8 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance -Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle +Double_t dEbeam = 0.; // Beam energy spread (MeV) // sim001 sim003 sim004 sim006 //Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 //Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 sim005 sim006 diff --git a/prm/prm_reac.c11li6n14pd_sim007.C b/prm/prm_reac.c11li6n14pd_sim007.C index 27f26e1cc0d8f8b84380ba80e5a428539bae1ebf..12a0273439c3e2cb0682f30d8174267c5ed066fa 100644 --- a/prm/prm_reac.c11li6n14pd_sim007.C +++ b/prm/prm_reac.c11li6n14pd_sim007.C @@ -94,9 +94,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle //Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 sim005 +Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 sim005 //Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 sim005 sim006 Double_t xtar0 = 0.; // mm diff --git a/prm/prm_reac.c11li6n14pd_sim008.C b/prm/prm_reac.c11li6n14pd_sim008.C index ab9d2b581806920706873ac5480db6f14a3c4712..716c5d1fc31a58a36b4d5602c4c218044147db46 100644 --- a/prm/prm_reac.c11li6n14pd_sim008.C +++ b/prm/prm_reac.c11li6n14pd_sim008.C @@ -81,8 +81,6 @@ 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.6; // Target thickness in mg/cm^2 // sim007 -//float r_sphere = 250.; // radius of the TGeoVolume in mm -float r_sphere = 320.; // radius of the TGeoVolume in mm //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 @@ -96,9 +94,9 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance +Double_t Ebeam_ref_mid = 11.; // Mean beam energy at middle //Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 sim005 +Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 sim005 //Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 Double_t depe2 = 0.01; // Relative energy resolution DE/E @5MeV sim003 sim005 sim006 Double_t xtar0 = 0.; // mm @@ -120,7 +118,8 @@ Double_t sgm_ppac_b = 0.5; // mm, PPAC b intrinsic resolution sim002 sim005 sim0 /* 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 +Double_t thdet[ndet] = {-60., -40., -12., 12., 40., 60.}; // in deg sim001-007 +//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 xdet_off[ndet] = {0., 0., 0., 0., 0., 0.}; @@ -134,6 +133,8 @@ Double_t psd_half[ndet] = {25., 25., 20., 20., 25., 25.}; //sim008 Int_t n_strip[ndet] = {16, 16, 40, 40, 16, 16}; // sim004 sim005 sim006 sim007 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/prm_reac.c11li6n14pd_sim009.C b/prm/prm_reac.c11li6n14pd_sim009.C new file mode 100644 index 0000000000000000000000000000000000000000..b015e89751e7ca58e075addcb994651736ee008f --- /dev/null +++ b/prm/prm_reac.c11li6n14pd_sim009.C @@ -0,0 +1,144 @@ +/* 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 different forward detector setting +////////////////////////// + + +/* 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"; +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.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 +//Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 +Double_t dEbeam = 0.21; // Beam energy spread (MeV) // sim002 sim005 +//Double_t depe2 = 0.; // Relative energy resolution DE/E @5MeV sim001, sim002 sim004 +Double_t depe2 = 0.01; // 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., -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 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/prm_reac.c11li6n14pd_sim010.C b/prm/prm_reac.c11li6n14pd_sim010.C new file mode 100644 index 0000000000000000000000000000000000000000..9ccca434c1876ef52bf188f25c4d0ab8f45fee4a --- /dev/null +++ b/prm/prm_reac.c11li6n14pd_sim010.C @@ -0,0 +1,151 @@ +/* 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 index ab9d2b581806920706873ac5480db6f14a3c4712..9ccca434c1876ef52bf188f25c4d0ab8f45fee4a 100644 --- a/prm_reac.C +++ b/prm_reac.C @@ -1,6 +1,6 @@ /* Common constants */ const Double_t u = 931.49432/1000.; // GeV -const Double_t c = 299.7924; // speed of light, mm/ns +//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.; @@ -18,7 +18,9 @@ Char_t label[50] = "c11li6n14pd"; // /* 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] = "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 ////////////////////////// @@ -62,8 +64,9 @@ float Eaft, Eaft_beam; //Char_t ch2[34] = "POLYETHYLENE"; //Char_t cd2[34] = "cd2"; Char_t lif[34] = "lif"; -Char_t he[34] = "he"; +//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 @@ -79,10 +82,10 @@ float temperature = 299.; 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.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 r_sphere = 250.; // radius of the TGeoVolume in mm -float r_sphere = 320.; // radius of the TGeoVolume in mm //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 @@ -96,11 +99,12 @@ float Ebfr, Ebfr_beam; // Beam energy E in MeV after energy loss in matter // |---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 eprj_mid = 11./1000.; // Mean beam energy at entrance -//Double_t deprj = 0./1000.; // Beam energy spread (GeV) // sim001 sim003 sim004 -Double_t deprj = 0.21/1000.; // Beam energy spread (GeV) // sim002 sim005 +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.01; // Relative energy resolution DE/E @5MeV sim003 sim005 sim006 +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 */ @@ -120,20 +124,26 @@ Double_t sgm_ppac_b = 0.5; // mm, PPAC b intrinsic resolution sim002 sim005 sim0 /* 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 -//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 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 +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}; // 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/simthm_c11li6n14pd.C b/simthm_c11li6n14pd.C index 5af182341a74330c02e78f3cc4fc51a66c6e8fae..09c950c6922de45e41069f90bfac7d7351a66f28 100644 --- a/simthm_c11li6n14pd.C +++ b/simthm_c11li6n14pd.C @@ -112,10 +112,10 @@ int main() { // For enewz for(Int_t ichar=0;ichar<33;ichar++){ - if(he[ichar]=='\0'){ - he[ichar]=' '; - he[ichar+1]='\0'; - } + /* if(he[ichar]=='\0'){ */ + /* he[ichar]=' '; */ + /* he[ichar+1]='\0'; */ + /* } */ if(si[ichar]=='\0'){ si[ichar]=' '; si[ichar+1]='\0'; @@ -124,11 +124,70 @@ int main() { lif[ichar]=' '; lif[ichar+1]='\0'; } + if(c[ichar]=='\0'){ + c[ichar]=' '; + c[ichar+1]='\0'; + } } // to store the parameter file gROOT->ProcessLine(Form(".! cp -pv prm_reac.C prm/prm_reac.%s_%s.C", label,simnum)); + + ////////////////////////////////////// + // 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("mkdir -vp figs"); // Creates figs/ directory if there is not */ + /* c0->Print(Form("figs/det_arrange_%s_%s.png",label,simnum)); */ + + ////////////////////////////////////// + // Output tree file + gSystem->Exec("mkdir -vp root"); // Creates root/ directory if there is not Char_t fname[200]; sprintf(fname,"root/%s_%s.root",label,simnum); cout << "Output file name: " << fname << endl; @@ -201,50 +260,10 @@ int main() { Int_t iex; sim->Branch("iex",&iex,"iex/I"); for(Int_t iexci=0; iexciSetPalette(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(); - /////////////////////////////////////////////////////////////// TVector3 vhit[2][ndet]; TVector3 vhit_MC[2]; @@ -256,11 +275,6 @@ int main() { TLorentzVector LVsum[3], LVsum_MC[3]; // #1+#2, #2+#3, #3+#1 TVector3 vel[3], vel_MC[3]; - TCanvas *c0 = new TCanvas("c0","c0",1450,0,800,800); - top->Draw(); - c0->Update(); - gSystem->ProcessEvents(); - TLorentzVector LV_target(0.0, 0.0, 0.0, u*mtrg); // px, py, pz, e TGenPhaseSpace event; @@ -321,8 +335,10 @@ int main() { // |---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 = 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, lif, &unit_pressure, &pressure, @@ -375,7 +391,7 @@ int main() { 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]*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; @@ -439,9 +455,13 @@ int main() { // 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]; + float fmass = mass[iprt]/u; enewzsub_etot_(&atom_num[iprt], &fmass, &Ebfr, lif, &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, c, &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] @@ -473,8 +493,13 @@ int main() { // 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, c, &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_(&zrec, &mrec, &Eaft, he, &unit_pressure, &pressure, + eoldzsub_etot_(&atom_num[iprt], &fmass, &Eaft, lif, &unit_pressure, &pressure, &temperature, &unit_mgpcm2, &fthk_mid_end_MC, &Ebfr); // -> E.lab_MC[0,1] E.lab_MC[iprt] = Ebfr; // MeV @@ -503,7 +528,7 @@ int main() { 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]*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; @@ -524,12 +549,6 @@ int main() { cout << "Done!" << endl; //return 0; - TView3D *view = (TView3D*)TView::CreateView(1); - view->SetRange(-250,-250,-250,250,250,250); - view->Front(); - view->ShowAxis(); - //c0->Update();gSystem->ProcessEvents();getchar(); - //c0->Print("detector_array_sim01.png"); theApp.Run();