diff --git a/f17do14an/figs/det_arrange_f17do14an_sim001.png b/f17do14an/figs/det_arrange_f17do14an_sim001.png new file mode 100644 index 0000000000000000000000000000000000000000..48f57a0b290b95e9340668d43fbcf439760a81eb Binary files /dev/null and b/f17do14an/figs/det_arrange_f17do14an_sim001.png differ diff --git a/f17do14an/figs/det_arrange_f17do14an_sim001.root b/f17do14an/figs/det_arrange_f17do14an_sim001.root new file mode 100644 index 0000000000000000000000000000000000000000..eeb928b65af23a854a841b999ce2f68c8978e35e Binary files /dev/null and b/f17do14an/figs/det_arrange_f17do14an_sim001.root differ diff --git a/f17do14an/figs/f17do14an_sim001_cm.png b/f17do14an/figs/f17do14an_sim001_cm.png new file mode 100644 index 0000000000000000000000000000000000000000..fedfad4b457f835c384c99bc8b797aa4b6e9af5a Binary files /dev/null and b/f17do14an/figs/f17do14an_sim001_cm.png differ diff --git a/f17do14an/figs/f17do14an_sim001_kin.png b/f17do14an/figs/f17do14an_sim001_kin.png new file mode 100644 index 0000000000000000000000000000000000000000..cee0d49c96df8366475ce6ae8a4df1b59eb67626 Binary files /dev/null and b/f17do14an/figs/f17do14an_sim001_kin.png differ diff --git a/f17do14an/figs/f17do14an_sim001_lab.png b/f17do14an/figs/f17do14an_sim001_lab.png new file mode 100644 index 0000000000000000000000000000000000000000..280c9c0a646e5a6a80cc48ef47bab70014b6e814 Binary files /dev/null and b/f17do14an/figs/f17do14an_sim001_lab.png differ diff --git a/f17do14an/figs/f17do14an_sim001_pos.png b/f17do14an/figs/f17do14an_sim001_pos.png new file mode 100644 index 0000000000000000000000000000000000000000..70546abbdfbe809032dd7c711af263b5e613333e Binary files /dev/null and b/f17do14an/figs/f17do14an_sim001_pos.png differ diff --git a/f17do14an/figs/f17do14an_sim001_rel.png b/f17do14an/figs/f17do14an_sim001_rel.png new file mode 100644 index 0000000000000000000000000000000000000000..43303a21b90bbaa90db263e40c658a0cad2d6dff Binary files /dev/null and b/f17do14an/figs/f17do14an_sim001_rel.png differ diff --git a/f17do14an/figs/sa_coin_f17do14an_sim001.pdf b/f17do14an/figs/sa_coin_f17do14an_sim001.pdf new file mode 100644 index 0000000000000000000000000000000000000000..f664ecddf251ed11d59d2127f535b10f31ae2c95 Binary files /dev/null and b/f17do14an/figs/sa_coin_f17do14an_sim001.pdf differ diff --git a/f17do14an/figs/sa_coin_f17do14an_sim001.png b/f17do14an/figs/sa_coin_f17do14an_sim001.png new file mode 100644 index 0000000000000000000000000000000000000000..f0900909eb46e33b709bb336dffae7771a8705cd Binary files /dev/null and b/f17do14an/figs/sa_coin_f17do14an_sim001.png differ diff --git a/f17do14an/macro/Makefile b/f17do14an/macro/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..2c9043a1385e13a39b74b1e51e5610e2363c3303 --- /dev/null +++ b/f17do14an/macro/Makefile @@ -0,0 +1,87 @@ +# name of .C program +#TARGET = plot_yield_coin +#TARGET = check_sa_lab +TARGET = check_sa_cm + +DEP = $(TARGET).d +SRCS = $(TARGET).C +OBJS = $(addsuffix .o, $(basename $(SRCS))) + +ROOTLIBDIR = $(shell root-config --libdir) +ROOTCFLAGS = $(shell root-config --cflags) +CXXFLAGS = $(ROOTCFLAGS) -Wall +F77 = gfortran + +# 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 + +# 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) + +.f.o: + ${F77} -c $< + +dep: + $(CC) -MM $(SRCS) > $(DEP) $(CXXFLAGS) +-include $(DEP) + +clean: + 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/f17do14an/macro/check_sa_cm b/f17do14an/macro/check_sa_cm new file mode 100755 index 0000000000000000000000000000000000000000..106d3953e5e4db40f497cf9878aedb9886559418 Binary files /dev/null and b/f17do14an/macro/check_sa_cm differ diff --git a/f17do14an/macro/check_sa_cm.C b/f17do14an/macro/check_sa_cm.C new file mode 100644 index 0000000000000000000000000000000000000000..e65ba4f31113bee0cdea198470b7826e8003e517 --- /dev/null +++ b/f17do14an/macro/check_sa_cm.C @@ -0,0 +1,348 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* #include */ +/* #include */ +/* #include */ +#include +#include +#include + +using namespace std; + + +//#include "../prm_reac.C" // parameter file <-- +#include "../prm/prm_reac.f17do14an_sim001.C" // parameter file <-- + +//Bool_t inv_thlab = true; +Bool_t inv_thlab = false; + +TFile *f = new TFile(Form("../root/%s_%s.root",label,simnum)); +TTree *tr = (TTree*)f->Get("sim"); + +TF1 *fn1 = new TF1("fn1","gaus",-100,100); + +int main(){ + const Double_t pi = TMath::Pi(); + + // No. of injected beam particles + const Double_t Ibeam = 8.e5; // pps + const Double_t Dbeam = 9.; // days + //const Double_t nbeam = 3.e10; // <=> 5e4 pps, 7 days + //const Double_t nbeam = Ibeam*Dbeam*24.*60.*60.; + const Double_t nbeam = Ibeam*60.*60.; // 1 hour + + // No. of injected beam particles + const Double_t Nav = 6.0221e23; + //const Double_t C_cm2_mb_mg = Nav/mtrg*1.e-30; + const Double_t C_cm2_mb_mg = 2.*Nav/(12. + 2.*mtrg)*1.e-30; + //cout << C_cm2_mb_mg << endl; + //const Double_t cs = 19.2; //mb + const Double_t cs = 0.13; //mb + const Double_t qfeff = 0.01; + + // Bin deffinition + //const Int_t necm = 140; // Ebin = 50 keV + const Int_t necm = 200; // Ebin = 100 keV + const Double_t ecmmin = 0.; + const Double_t ecmmax = 11.; + const Double_t decm = (ecmmax - ecmmin)/double(necm); + const Int_t nelab1 = 120; + const Double_t elab1min = 50.; + const Double_t elab1max = 110.; + const Double_t delab1 = (elab1max - elab1min)/double(nelab1); + + + + //const Int_t nthlab = 18; // theta bin = 5 deg. + const Int_t nthlab = 90; // theta bin = 1 deg. + //const Int_t nthlab = 10; // theta bin = 9 deg. + const Double_t thlabmin = 0.; + // const Double_t thlabmax = 180.; + const Double_t thlabmax = 90.; + const Double_t dthlab = (thlabmax - thlabmin)/double(nthlab); + const Double_t d2r = pi/180.; + + //Int_t ibinmin = 0; // including under flow + Int_t ibinmin = 1; // fisrt bin + Int_t ibinmax = 1; + Double_t dsdO_max = 800.; + + // Cut condition + const char* cut = + //"hit.idet==0&&iex==0"; + //"hit.idet>=0&&hit.idet<=2&&iex==0" + //"hit.any&&iex==0"; + "hit.idet[0]==0&&hit.idet[1]==3&&iex==0"; + //"&&pos.ztar>-236."; // to avoid events near the entrance + //"&&pos.ztar>-283."; // to avoid events near the entrance + // "&&(E.lab[0] < 32. || E.lab[0]>=32.&&(E.cm > 0.19*E.lab[0]+0.22))"; + + const Double_t lbsize = 0.05; + //const Double_t txsize = 0.06; + + TApplication theApp("App", NULL, NULL); + gStyle->SetTextFont(62); + gStyle->SetTextSize(lbsize); + + gStyle->SetLabelFont(62,"XYZt"); + gStyle->SetLabelSize(lbsize,"XYZt"); + + //gStyle->SetTitleFont(62); + gStyle->SetTitleFont(62,"XYZt"); + gStyle->SetTitleSize(lbsize,"XYZt"); + gStyle->SetTitleFontSize(lbsize); + gStyle->SetTitleYOffset(1.2); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetPadBottomMargin(0.1); + + // gStyle->SetStatW(0.3); + // gStyle->SetStatH(0.2); + // gStyle->SetStatX(1); + // gStyle->SetStatY(0.4); + gStyle->SetOptStat(0); + + + // Canvas + //TCanvas *c1=new TCanvas("c1","c1",0,0,880,800); + TCanvas *c1=new TCanvas("c1","c1",1100,0,1400,300); + //c1->Divide(3,3); + c1->Divide(5,1); + c1->cd(1); + + TCanvas *c2=new TCanvas("c2","c2",1100,350,1400,300); + //c2->Divide(3,3); + c2->Divide(5,1); + c2->cd(1); + + + //Char_t hname[100]; + TString hname; + + // TH2 for Nhit (Ecm vs thlab) + //sprintf(hname,"hhit"); + hname = "hhit1"; + TH2D *hhit1 = new TH2D(hname.Data(),"Simulated yield hitting detectors;" + "E_{c.m.} (MeV);#theta_{lab. 1} (deg.)", + //necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + if(inv_thlab) + //tr->Draw(Form("180. - th.lab[0]:E.cm>>%s",hname.Data()),cut); + tr->Draw(Form("180. - th.lab[0]:E.cm>>%s",hname.Data()),cut); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),cut); + else + //tr->Draw(Form("th.lab[0]:E.cm>>%s",hname.Data()),cut); + tr->Draw(Form("th.lab[0]:E.cm>>%s",hname.Data()),cut); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),cut); + + hname = "hhit2"; + TH2D *hhit2 = new TH2D(hname.Data(),"Simulated yield hitting detectors;" + "E_{c.m.} (MeV);#theta_{lab. 2} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[1]:E.cm>>%s",hname.Data()),cut); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),cut); + else + tr->Draw(Form("th.lab[1]:E.cm>>%s",hname.Data()),cut); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),cut); + + + // TH2 for Ntot (Ecm vs thlab) + //sprintf(hname,"htot"); + hname = "htot1"; + TH2D *htot1 = new TH2D(hname.Data(),"Total simulated yield;" + "E_{c.m.} (MeV);#theta_{lab. 1} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[0]:E.cm>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + else + tr->Draw(Form("th.lab[0]:E.cm>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + + hname = "htot2"; + TH2D *htot2 = new TH2D(hname.Data(),"Total simulated yield;" + "E_{c.m.} (MeV);#theta_{lab. 2} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[1]:E.cm>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + else + tr->Draw(Form("th.lab[1]:E.cm>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + + + //sprintf(hname,"hsin_int"); + hname = "hsin_int"; + TH2D *hsin_int = new TH2D(hname.Data(), + "Integrated sin(#theta) in #theta bin size;" + "E_{c.m.} (MeV);#theta_{lab. 1} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + Double_t ecm = ecmmin + decm/2.; + for(Int_t i=0; iFill( ecm,thlab, - ( cos(d2r*(thlab+dthlab/2.)) + - cos(d2r*(thlab-dthlab/2.)) ) ); + // hsin_int->Fill( ecm,thlab, - ( cos(d2r*((180.-2.*thlab)+dthlab)) + // - cos(d2r*((180.-2.*thlab)-dthlab)) ) ); + thlab += dthlab; + } + ecm += decm; + } + + // dO: TH2 Divide (Ecm vs thlab) + TH2D *hdO_1 = new TH2D("hdO_1","Solid angle (msr);" + "E_{c.m.} (MeV);#theta_{lab. 1} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + hdO_1->Divide(hhit1,htot1); //,4.*pi); + /* hdO->Draw("colz"); */ + /* c1->Update(); gSystem->ProcessEvents(); getchar(); */ + + /* /\* // Reading yield data *\/ */ + /* /\* //sprintf(hname,"hyld"); *\/ */ + /* /\* hname = "hyld"; *\/ */ + /* /\* TH2D *hyld = new TH2D(hname.Data(),"Experimental yield;" *\/ */ + /* /\* "E_{c.m.} (MeV);#theta_{lab.} (deg.)", *\/ */ + /* /\* necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); *\/ */ + /* /\* ifstream fdat("yield_si26aa_ecm_thlab_tel1.dat"); *\/ */ + /* /\* Double_t xdat, ydat; *\/ */ + /* /\* while( !fdat.eof() ){ *\/ */ + /* /\* fdat >> xdat >> ydat; *\/ */ + /* /\* hyld->Fill(xdat,ydat); *\/ */ + /* /\* } *\/ */ + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thlab) + TH2D *hdsdO_1 = new TH2D("hdsdO_1","Cross section (mb/sr);" + "E_{c.m.} (MeV);#theta_{lab. 1} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + + hdO_1->Multiply(hdO_1,hsin_int,2.*pi*1000.); + + // dO: TH2 Divide (Ecm vs thlab) + TH2D *hdO_2 = new TH2D("hdO_2","Solid angle (msr);" + "E_{c.m.} (MeV);#theta_{lab. 2} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + hdO_2->Divide(hhit2,htot2); //,4.*pi); + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thlab) + TH2D *hdsdO_2 = new TH2D("hdsdO_2","Cross section (mb/sr);" + "E_{c.m.} (MeV);#theta_{lab. 2} (deg.)", + necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + + hdO_2->Multiply(hdO_2,hsin_int,2.*pi*1000.); + + + c1->cd(1); + hhit1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(2); + htot1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(3); + hsin_int->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(4); + gPad->SetLogz(0); + hdO_1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(5); + TH1D * hdO_1_px = hdO_1->ProjectionX(); + hdO_1_px->Draw("hist"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + ////// + + c2->cd(1); + hhit2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(2); + htot2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(3); + hsin_int->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(4); + gPad->SetLogz(0); + hdO_2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(5); + TH1D * hdO_2_px =hdO_2->ProjectionX(); + hdO_2_px->Draw("hist"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + /* c1->Print(Form("../figs/sa_%s_%s_telall.png", */ + /* label,simnum)); */ + /* c1->Print(Form("../figs/sa_%s_%s_telall.pdf", */ + /* label,simnum)); */ + + TCanvas *c3=new TCanvas("c3","c3",0,700,1400,300); + c3->Divide(4,1); + + c3->cd(1); + /* hname = "hhit_cm_lab"; */ + /* TH2D *hhit_cm_lab = new TH2D(hname.Data(),"Simulated yield hitting detectors;" */ + /* "E_{c.m.} (MeV);E_{c.m.} (MeV)", */ + /* necm,ecmmin,ecmmax,necm,ecmmin,ecmmax); */ + /* tr->Draw(Form("E.cm:E.cm>>%s",hname.Data()),cut); */ + /* hhit_cm_lab->Draw("colz"); */ + hname = "hhit"; + TH1D *hhit = new TH1D(hname.Data(),"Simulated yield hitting detectors;" + "E_{c.m.} (MeV);Counts/bin", + necm,ecmmin,ecmmax); + tr->Draw(Form("E.cm>>%s",hname.Data()),cut); + hhit->Draw(); + + hname = "hhit30"; + TH1D *hhit30 = new TH1D(hname.Data(),"Simulated yield hitting detectors;" + "E_{c.m.} (MeV);Counts/bin", + necm,ecmmin,ecmmax); + tr->Draw(Form("E.cm>>%s",hname.Data()),Form("abs(p.spc)<30&&%s",cut)); + hhit->Draw(); + hhit30->SetLineColor(kRed); + hhit30->Draw("same"); + + + + c3->Update(); gSystem->ProcessEvents(); + + c3->cd(2); + hname = "hyld"; + TH1D *hyld = new TH1D(hname.Data(),"Yield;" + "E_{c.m.} (MeV);Counts/bin", + necm,ecmmin,ecmmax); + hyld->Multiply(hdO_1_px,hdO_2_px,nbeam*C_cm2_mb_mg*thk_tar*cs/16./pi/pi*qfeff); + //hyld->Multiply(hyld,hdO_2_px); + hyld->Draw("hist"); + + Double_t yld_int = hyld->Integral(0,necm); + cout << hyld->GetMaximum() << endl; + TLatex *tx1 = new TLatex(1,0.9*hyld->GetMaximum(), + Form("Yield = %7.2f/h",yld_int)); + tx1->Draw(); + + c3->Update(); gSystem->ProcessEvents(); + + theApp.Run(); +} diff --git a/f17do14an/macro/check_sa_lab b/f17do14an/macro/check_sa_lab new file mode 100755 index 0000000000000000000000000000000000000000..de2d0f83840634d4e85f35b483bd9cd717e25dd7 Binary files /dev/null and b/f17do14an/macro/check_sa_lab differ diff --git a/f17do14an/macro/check_sa_lab.C b/f17do14an/macro/check_sa_lab.C new file mode 100644 index 0000000000000000000000000000000000000000..6561fd9fbf6d82e6cdb1062a5a7df608fba68459 --- /dev/null +++ b/f17do14an/macro/check_sa_lab.C @@ -0,0 +1,340 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* #include */ +/* #include */ +/* #include */ +#include +#include +#include + +using namespace std; + + +//#include "../prm_reac.C" // parameter file <-- +#include "../prm/prm_reac.f17do14an_sim001.C" // parameter file <-- + +//Bool_t inv_thlab = true; +Bool_t inv_thlab = false; + +TFile *f = new TFile(Form("../root/%s_%s.root",label,simnum)); +TTree *tr = (TTree*)f->Get("sim"); + +TF1 *fn1 = new TF1("fn1","gaus",-100,100); + +int main(){ + const Double_t pi = TMath::Pi(); + + // No. of injected beam particles + const Double_t Ibeam = 5.e5; // pps + const Double_t Dbeam = 7.; // days + //const Double_t nbeam = 3.e10; // <=> 5e4 pps, 7 days + //const Double_t nbeam = Ibeam*Dbeam*24.*60.*60.; + const Double_t nbeam = Ibeam*60.; + + // No. of injected beam particles + const Double_t Nav = 6.0221e23; + //const Double_t C_cm2_mb_mg = Nav/mtrg*1.e-30; + const Double_t C_cm2_mb_mg = 2.*Nav/(12. + 2.*mtrg)*1.e-30; + //cout << C_cm2_mb_mg << endl; + const Double_t cs = 19.2; //mb + const Double_t qfeff = 0.05; + + // Bin deffinition + //const Int_t necm = 140; // Ebin = 50 keV + const Int_t necm = 110; // Ebin = 100 keV + const Double_t ecmmin = 0.; + const Double_t ecmmax = 11.; + const Double_t decm = (ecmmax - ecmmin)/double(necm); + const Int_t nelab1 = 120; + const Double_t elab1min = 50.; + const Double_t elab1max = 110.; + const Double_t delab1 = (elab1max - elab1min)/double(nelab1); + + + + //const Int_t nthlab = 18; // theta bin = 5 deg. + const Int_t nthlab = 90; // theta bin = 1 deg. + //const Int_t nthlab = 10; // theta bin = 9 deg. + const Double_t thlabmin = 0.; + // const Double_t thlabmax = 180.; + const Double_t thlabmax = 90.; + const Double_t dthlab = (thlabmax - thlabmin)/double(nthlab); + const Double_t d2r = pi/180.; + + //Int_t ibinmin = 0; // including under flow + Int_t ibinmin = 1; // fisrt bin + Int_t ibinmax = 1; + Double_t dsdO_max = 800.; + + // Cut condition + const char* cut = + //"hit.idet==0&&iex==0"; + //"hit.idet>=0&&hit.idet<=2&&iex==0" + //"hit.any&&iex==0"; + "hit.idet[0]==0&&hit.idet[1]==3&&iex==0"; + //"&&pos.ztar>-236."; // to avoid events near the entrance + //"&&pos.ztar>-283."; // to avoid events near the entrance + // "&&(E.lab[0] < 32. || E.lab[0]>=32.&&(E.cm > 0.19*E.lab[0]+0.22))"; + + const Double_t lbsize = 0.05; + //const Double_t txsize = 0.06; + + TApplication theApp("App", NULL, NULL); + gStyle->SetTextFont(62); + gStyle->SetTextSize(lbsize); + + gStyle->SetLabelFont(62,"XYZt"); + gStyle->SetLabelSize(lbsize,"XYZt"); + + //gStyle->SetTitleFont(62); + gStyle->SetTitleFont(62,"XYZt"); + gStyle->SetTitleSize(lbsize,"XYZt"); + gStyle->SetTitleFontSize(lbsize); + gStyle->SetTitleYOffset(1.2); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetPadBottomMargin(0.1); + + // gStyle->SetStatW(0.3); + // gStyle->SetStatH(0.2); + // gStyle->SetStatX(1); + // gStyle->SetStatY(0.4); + gStyle->SetOptStat(0); + + + // Canvas + //TCanvas *c1=new TCanvas("c1","c1",0,0,880,800); + TCanvas *c1=new TCanvas("c1","c1",1100,0,1400,300); + //c1->Divide(3,3); + c1->Divide(5,1); + c1->cd(1); + + TCanvas *c2=new TCanvas("c2","c2",1100,350,1400,300); + //c2->Divide(3,3); + c2->Divide(5,1); + c2->cd(1); + + + //Char_t hname[100]; + TString hname; + + // TH2 for Nhit (Ecm vs thlab) + //sprintf(hname,"hhit"); + hname = "hhit1"; + TH2D *hhit1 = new TH2D(hname.Data(),"Simulated yield hitting detectors;" + "E_{lab. 1} (MeV);#theta_{lab. 1} (deg.)", + //necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + if(inv_thlab) + //tr->Draw(Form("180. - th.lab[0]:E.cm>>%s",hname.Data()),cut); + tr->Draw(Form("180. - th.lab[0]:E.lab[0]>>%s",hname.Data()),cut); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),cut); + else + //tr->Draw(Form("th.lab[0]:E.cm>>%s",hname.Data()),cut); + tr->Draw(Form("th.lab[0]:E.lab[0]>>%s",hname.Data()),cut); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),cut); + + hname = "hhit2"; + TH2D *hhit2 = new TH2D(hname.Data(),"Simulated yield hitting detectors;" + "E_{lab. 1} (MeV);#theta_{lab. 2} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[1]:E.lab[0]>>%s",hname.Data()),cut); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),cut); + else + tr->Draw(Form("th.lab[1]:E.lab[0]>>%s",hname.Data()),cut); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),cut); + + + // TH2 for Ntot (Ecm vs thlab) + //sprintf(hname,"htot"); + hname = "htot1"; + TH2D *htot1 = new TH2D(hname.Data(),"Total simulated yield;" + "E_{lab. 1} (MeV);#theta_{lab. 1} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[0]:E.lab[0]>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + else + tr->Draw(Form("th.lab[0]:E.lab[0]>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + + hname = "htot2"; + TH2D *htot2 = new TH2D(hname.Data(),"Total simulated yield;" + "E_{lab. 1} (MeV);#theta_{lab. 2} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + if(inv_thlab) + tr->Draw(Form("180. - th.lab[1]:E.lab[0]>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("180. - th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + else + tr->Draw(Form("th.lab[1]:E.lab[0]>>%s",hname.Data()),"iex==0"); + //tr->Draw(Form("th.lab2_MC:E.cm>>%s",hname.Data()),"iex==0"); + + + //sprintf(hname,"hsin_int"); + hname = "hsin_int"; + TH2D *hsin_int = new TH2D(hname.Data(), + "Integrated sin(#theta) in #theta bin size;" + "E_{lab. 1} (MeV);#theta_{lab. 1} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + Double_t elab1 = elab1min + delab1/2.; + for(Int_t i=0; iFill( elab1,thlab, - ( cos(d2r*(thlab+dthlab/2.)) + - cos(d2r*(thlab-dthlab/2.)) ) ); + // hsin_int->Fill( ecm,thlab, - ( cos(d2r*((180.-2.*thlab)+dthlab)) + // - cos(d2r*((180.-2.*thlab)-dthlab)) ) ); + thlab += dthlab; + } + elab1 += delab1; + } + + // dO: TH2 Divide (Ecm vs thlab) + TH2D *hdO_1 = new TH2D("hdO_1","Solid angle (msr);" + "E_{lab. 1} (MeV);#theta_{lab. 1} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + hdO_1->Divide(hhit1,htot1); //,4.*pi); + /* hdO->Draw("colz"); */ + /* c1->Update(); gSystem->ProcessEvents(); getchar(); */ + + /* /\* // Reading yield data *\/ */ + /* /\* //sprintf(hname,"hyld"); *\/ */ + /* /\* hname = "hyld"; *\/ */ + /* /\* TH2D *hyld = new TH2D(hname.Data(),"Experimental yield;" *\/ */ + /* /\* "E_{c.m.} (MeV);#theta_{lab.} (deg.)", *\/ */ + /* /\* necm,ecmmin,ecmmax,nthlab,thlabmin,thlabmax); *\/ */ + /* /\* ifstream fdat("yield_si26aa_ecm_thlab_tel1.dat"); *\/ */ + /* /\* Double_t xdat, ydat; *\/ */ + /* /\* while( !fdat.eof() ){ *\/ */ + /* /\* fdat >> xdat >> ydat; *\/ */ + /* /\* hyld->Fill(xdat,ydat); *\/ */ + /* /\* } *\/ */ + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thlab) + TH2D *hdsdO_1 = new TH2D("hdsdO_1","Cross section (mb/sr);" + "E_{lab. 1} (MeV);#theta_{lab. 1} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + + hdO_1->Multiply(hdO_1,hsin_int,2.*pi*1000.); + + // dO: TH2 Divide (Ecm vs thlab) + TH2D *hdO_2 = new TH2D("hdO_2","Solid angle (msr);" + "E_{lab. 1} (MeV);#theta_{lab. 2} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + hdO_2->Divide(hhit2,htot2); //,4.*pi); + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thlab) + TH2D *hdsdO_2 = new TH2D("hdsdO_2","Cross section (mb/sr);" + "E_{lab. 1} (MeV);#theta_{lab. 2} (deg.)", + nelab1,elab1min,elab1max,nthlab,thlabmin,thlabmax); + + hdO_2->Multiply(hdO_2,hsin_int,2.*pi*1000.); + + + c1->cd(1); + hhit1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(2); + htot1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(3); + hsin_int->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(4); + gPad->SetLogz(0); + hdO_1->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(5); + TH1D * hdO_1_px = hdO_1->ProjectionX(); + hdO_1_px->Draw("hist"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + ////// + + c2->cd(1); + hhit2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(2); + htot2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(3); + hsin_int->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(4); + gPad->SetLogz(0); + hdO_2->Draw("colz"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + c2->cd(5); + TH1D * hdO_2_px =hdO_2->ProjectionX(); + hdO_2_px->Draw("hist"); + c2->Update(); gSystem->ProcessEvents(); //getchar(); + + /* c1->Print(Form("../figs/sa_%s_%s_telall.png", */ + /* label,simnum)); */ + /* c1->Print(Form("../figs/sa_%s_%s_telall.pdf", */ + /* label,simnum)); */ + + TCanvas *c3=new TCanvas("c3","c3",0,700,1400,300); + c3->Divide(4,1); + + c3->cd(1); + /* hname = "hhit_cm_lab"; */ + /* TH2D *hhit_cm_lab = new TH2D(hname.Data(),"Simulated yield hitting detectors;" */ + /* "E_{lab. 1} (MeV);E_{c.m.} (MeV)", */ + /* nelab1,elab1min,elab1max,necm,ecmmin,ecmmax); */ + /* tr->Draw(Form("E.cm:E.lab[0]>>%s",hname.Data()),cut); */ + /* hhit_cm_lab->Draw("colz"); */ + hname = "hhit"; + TH1D *hhit = new TH1D(hname.Data(),"Simulated yield hitting detectors;" + "E_{lab. 1} (MeV);Counts/bin", + nelab1,elab1min,elab1max); + tr->Draw(Form("E.lab[0]>>%s",hname.Data()),cut); + hhit->Draw(); + + hname = "hhit30"; + TH1D *hhit30 = new TH1D(hname.Data(),"Simulated yield hitting detectors;" + "E_{lab. 1} (MeV);Counts/bin", + nelab1,elab1min,elab1max); + tr->Draw(Form("E.lab[0]>>%s",hname.Data()),Form("abs(p.spc)<30&&%s",cut)); + hhit->Draw(); + hhit30->SetLineColor(kRed); + hhit30->Draw("same"); + + + + c3->Update(); gSystem->ProcessEvents(); + + c3->cd(2); + hname = "hyld"; + TH1D *hyld = new TH1D(hname.Data(),"Yield;" + "E_{lab. 1} (MeV);Counts/bin", + nelab1,elab1min,elab1max); + hyld->Multiply(hhit,hdO_1_px,nbeam*C_cm2_mb_mg*thk_tar*cs/16./pi/pi*qfeff); + hyld->Multiply(hyld,hdO_2_px); + hyld->Draw("hist"); + c3->Update(); gSystem->ProcessEvents(); + + theApp.Run(); +} diff --git a/f17do14an/macro/chk_sim.C b/f17do14an/macro/chk_sim.C new file mode 100644 index 0000000000000000000000000000000000000000..23cec6077d7b5acd39517b73883672f84808199a --- /dev/null +++ b/f17do14an/macro/chk_sim.C @@ -0,0 +1,371 @@ +{ + #include "prm/prm_reac.f17do14an_sim001.C" + + /* Char_t label[50] = "n14dc11an"; */ + /* Char_t simnum[10] = "sim001"; */ + //TFile *f = new TFile(Form("../root/%s_%s.root",label,sumnum)); + + //q2[0] = (m14o+ma - m17f-mh)*u*1000.; + q2[0] = 0.; + /* cout << q2[0] << endl; */ + /* getchar(); */ + + // { + gStyle->SetStatW(0.3); + gStyle->SetStatH(0.3); + gStyle->SetStatFont(62); + gStyle->SetLabelFont(62,"XY"); + gStyle->SetLabelSize(0.06,"XY"); + gStyle->SetTitleFont(62,"XY"); + gStyle->SetTitleSize(0.06,"XY"); + gStyle->SetTextFont(62); + gStyle->SetTextSize(0.06); + // } + + TCut pscut50 = "abs(p.spc)<50."; + TCut pscut30 = "abs(p.spc)<30."; + Int_t color50 = kOrange-3; + Int_t color30 = kRed; + Int_t color_coin = kCyan; + + Double_t Ecm_min = 0.; + Double_t Ecm_max = 12.; + + Double_t Elab0_min = 40.; + Double_t Elab0_max = 120.; + Double_t Elab1_min = 0.; + Double_t Elab1_max = 60.; + + Double_t thlab0_min = 0.; + Double_t thlab0_max = 15.; + Double_t thlab1_min = 0.; + Double_t thlab1_max = 40.; + + Double_t pspec_max = 250.; + + Int_t cn_x0 = 1080; + Int_t cn_y0 = 300; + Int_t cn_dx = 50; + Int_t cn_dy = 50; + + //Double_t Q3 = -5.14649; // Q value of 3-body channel + Double_t Q3 = (mprj+mtrg-mres_gs-mrec-mspc)*u*1000.; + //cout << Q3 << endl; + + // { + TCanvas *c1=new TCanvas("c1","Kinematics",cn_x0,cn_y0,1200,800); + c1->Divide(3,2); + + c1->cd(1); + sim->Draw(Form("E.rel[0]:E.rel[2]>>h11(40,0,%f,40,0,%f)", + // Ecm_max+q2[0],Ecm_max+q2[0]),"",""); + 10.,10.),"",""); + sim->SetMarkerColor(color50); + sim->Draw("E.rel[0]:E.rel[2]",pscut50,"same"); + sim->SetMarkerColor(color30); + sim->Draw("E.rel[0]:E.rel[2]",pscut30,"same"); + sim->SetMarkerColor(kBlack); + TLatex *tx1=new TLatex(5,8,"All"); + TLatex *tx2=new TLatex(5,7.2,"|p_{s}| < 50 MeV/c"); + TLatex *tx3=new TLatex(5,6.4,"|p_{s}| < 30 MeV/c"); + tx1->SetTextColor(kBlack); + tx2->SetTextColor(color50); + tx3->SetTextColor(color30); + tx1->Draw(); + tx2->Draw(); + tx3->Draw(); + + + c1->cd(2); + sim->Draw(Form("E.lab[1]:E.lab[0]>>h12(80,%f,%f,40,%f,%f)", + Elab0_min,Elab0_max,Elab1_min,Elab1_max), + "",""); + sim->SetMarkerColor(color50); + sim->Draw("E.lab[1]:E.lab[0]",pscut50,"same"); + sim->SetMarkerColor(color30); + sim->Draw("E.lab[1]:E.lab[0]",pscut30,"same"); + sim->SetMarkerColor(kBlack); + + c1->cd(3); + sim->Draw(Form("th.lab[1]:th.lab[0]>>h13(80,%f,%f,40,%f,%f)", + thlab0_min,thlab0_max,thlab1_min,thlab1_max), + "",""); + sim->SetMarkerColor(color50); + sim->Draw("th.lab[1]:th.lab[0]",pscut50,"same"); + sim->SetMarkerColor(color30); + sim->Draw("th.lab[1]:th.lab[0]",pscut30,"same"); + sim->SetMarkerColor(color_coin); + sim->SetMarkerStyle(6); + sim->Draw("th.lab[1]:th.lab[0]","hit.both","same"); + sim->SetMarkerColor(kBlack); + sim->SetMarkerStyle(1); + TLatex *tx4=new TLatex(8,36,"Coincidence"); + tx4->SetTextColor(color_coin); + tx4->Draw(); + + c1->cd(4); + sim->Draw(Form("E.cm-%f:abs(p.spc)>>h14(80,0,%f,40,%f,%f)", + q2[0],pspec_max,Ecm_min,Ecm_max), + "",""); + sim->SetMarkerColor(color50); + sim->Draw(Form("E.cm-%f:abs(p.spc)",q2[0]),pscut50,"same"); + sim->SetMarkerColor(color30); + sim->Draw(Form("E.cm-%f:abs(p.spc)",q2[0]),pscut30,"same"); + sim->SetMarkerColor(kBlack); + + c1->cd(5); + sim->SetMarkerColor(kBlack); + sim->SetLineColor(kBlack); + sim->Draw(Form("E.cm-%f>>h15(100,%f,%f)",q2[0],Ecm_min,Ecm_max),pscut50*"w.evt","E"); + sim->SetMarkerColor(color50); + sim->SetLineColor(color50); + sim->Draw(Form("E.cm-%f",q2[0]),pscut50*"hit.both&&w.evt","Esame"); + sim->SetMarkerColor(color30); + sim->SetLineColor(color30); + sim->Draw(Form("E.cm-%f",q2[0]),pscut30*"hit.both&&w.evt","Esame"); + gStyle->SetStatX(0.4); + sim->SetLineColor(kBlack); + h15->SetStats(0); + + c1->cd(6); + sim->SetMarkerColor(color50); + sim->Draw(Form("th.cm:E.cm-%f>>h16(100,%f,%f,100,0,180)",q2[0],Ecm_min,Ecm_max), + pscut50*"w.evt"); + sim->SetMarkerColor(color_coin); + sim->SetMarkerStyle(7); + sim->Draw(Form("th.cm:E.cm-%f>>h16(100,%f,%f,100,0,180)",q2[0],Ecm_min,Ecm_max), + pscut50*"w.evt*hit.both","same"); + sim->SetMarkerColor(kBlack); + sim->SetMarkerStyle(1); + gStyle->SetStatX(0.9); + + // } + + + /* // { */ + /* TCanvas *c2=new TCanvas("c2","E.lab & th.lab", */ + /* cn_x0+cn_dx,cn_y0+cn_dy,1600,800); */ + /* c2->Divide(4,2); */ + + /* // E.lab */ + /* c2->cd(1); */ + /* sim->Draw(Form("E.lab_det[0]:(E.lab_det_MC[0] - E.lab_det[0])" */ + /* ">>h21(80,-0.3,0.3,40,%f,%f)",Elab0_min,Elab0_max), */ + /* "hit.both",""); */ + /* c2->cd(5); */ + /* sim->Draw("E.lab_det_MC[0] - E.lab_det[0]" */ + /* ">>h23(100,-0.3,0.3)","hit.both"); */ + + /* c2->cd(2); */ + /* sim->Draw(Form("E.lab_det[1]:(E.lab_det_MC[1] - E.lab_det[1])" */ + /* ">>h22(80,-0.3,0.3,40,%f,%f)",Elab0_min,Elab0_max), */ + /* "hit.both",""); */ + /* c2->cd(6); */ + /* sim->Draw("E.lab_det_MC[1] - E.lab_det[1]" */ + /* ">>h24(100,-0.3,0.3)","hit.both"); */ + /* /\* c2->cd(1); *\/ */ + /* /\* sim->Draw("E.lab[0] - E.lab_det[0]" *\/ */ + /* /\* ">>h21(100,-1.,3.)","hit.both",""); *\/ */ + /* /\* c2->cd(5); *\/ */ + /* /\* sim->Draw("E.lab_MC[0] - E.lab_det_MC[0]" *\/ */ + /* /\* ">>h23(100,-1.,3.)","hit.both"); *\/ */ + + /* /\* c2->cd(2); *\/ */ + /* /\* sim->Draw("E.lab[1] - E.lab_det[1]" *\/ */ + /* /\* ">>h22(100,-0.3,0.3)","hit.both",""); *\/ */ + /* /\* c2->cd(6); *\/ */ + /* /\* sim->Draw("E.lab_MC[1] - E.lab_det_MC[1]" *\/ */ + /* /\* ">>h24(100,-0.3,0.3)","hit.both"); *\/ */ + /* /\* return; *\/ */ + + + /* // th.lab */ + /* c2->cd(3); */ + /* sim->Draw("th.lab[0]:(th.lab_MC[0] - th.lab[0])" */ + /* ">>h25(80,-3,3,40,0,30)","hit.both"); */ + /* c2->cd(7); */ + /* sim->Draw("th.lab_MC[0] - th.lab[0]" */ + /* ">>h27(100,-3,3)","hit.both"); */ + /* c2->cd(4); */ + /* sim->Draw("th.lab[1]:(th.lab_MC[1] - th.lab[1])" */ + /* ">>h26(80,-3,3,40,0,70)","hit.both"); */ + /* c2->cd(8); */ + /* sim->Draw("th.lab_MC[1] - th.lab[1]" */ + /* ">>h28(100,-3,3)","hit.both"); */ + /* // } //return; */ + + /* // { */ + /* TCanvas *c3=new TCanvas("c3","pos", */ + /* cn_x0+2*cn_dx,cn_y0+2*cn_dy,1600,800); */ + /* c3->Divide(4,2); */ + + /* // pos.xhit */ + /* c3->cd(1); */ + /* // sim->Draw("pos.xhit:(pos.xhit_MC - pos.xhit)" */ + /* // ">>h31(200,-4.,4.,200,-30,30)","hit.both","colz"); */ + /* sim->Draw("pos.xhit_MC:pos.xhit" */ + /* ">>h31(100,-30.,30.,100,-30,30)","hit.both","colz"); */ + /* gStyle->SetStatY(0.); */ + /* c3->cd(5); */ + /* sim->Draw("pos.xhit_MC - pos.xhit" */ + /* ">>h33(100,-3.,3.)","hit.both&&abs(pos.xhit)<30"); */ + + /* // pos.yhit */ + /* c3->cd(2); */ + /* sim->Draw("pos.yhit_MC:pos.yhit" */ + /* ">>h32(100,-30.,30.,100,-30,30)","hit.both","colz"); */ + /* c3->cd(6); */ + /* sim->Draw("pos.yhit_MC - pos.yhit" */ + /* ">>h34(100,-3.,3.)","hit.both&&abs(pos.yhit)<30"); */ + + /* // pos.*tar */ + /* sim->SetMarkerStyle(7); */ + /* c3->cd(3); */ + /* sim->Draw("pos.ytar:pos.xtar" */ + /* ">>h35(80,-40,40,40,-40,40)","","colz"); */ + /* gStyle->SetStatY(0.9); */ + /* c3->cd(4); */ + /* sim->Draw("pos.ytar:pos.ytar_MC-pos.ytar" */ + /* ">>h36(80,-5,5,40,-40,40)","","colz"); */ + /* c3->cd(7); */ + /* sim->Draw("pos.xtar_MC-pos.xtar:pos.xtar" */ + /* ">>h37(80,-40,40,40,-5,5)","","colz"); */ + /* sim->SetMarkerStyle(1); */ + /* // } //return; */ + + + /* // { */ + /* TCanvas *c4=new TCanvas("c4","rel",cn_x0+3*cn_dx,cn_y0+3*cn_dy,1600,800); */ + /* c4->Divide(4,2); */ + + /* // E.beam */ + /* c4->cd(1); */ + /* sim->Draw(Form("E.beam_dec:(E.beam_mid - E.beam_dec)" */ + /* ">>h41(80,-1.,1.,40,%f,%f)",Ebeam_ref_mid-2,Ebeam_ref_mid+2), */ + /* "",""); */ + /* c4->cd(5); */ + /* sim->Draw("E.beam_mid - E.beam_dec" */ + /* ">>h45(100,-1.,1.)",""); */ + + /* // E.rel[0] */ + /* c4->cd(2); */ + /* sim->Draw("E.rel[0]:(E.rel_MC[0] - E.rel[0])" */ + /* ">>h42(80,-0.6,0.6,40,0,6)","hit.both",""); */ + /* c4->cd(6); */ + /* sim->Draw("E.rel_MC[0] - E.rel[0]" */ + /* ">>h46(100,-0.6,0.6)","hit.both"); */ + + /* // E.rel[1] */ + /* c4->cd(3); */ + /* sim->Draw("E.rel[1]:(E.rel_MC[1] - E.rel[1])" */ + /* ">>h43(80,-0.6,0.6,40,0,6)","hit.both",""); */ + /* c4->cd(7); */ + /* sim->Draw("E.rel_MC[1] - E.rel[1]" */ + /* ">>h47(100,-0.6,0.6)","hit.both"); */ + + /* // E.rel[2] */ + /* c4->cd(4); */ + /* sim->Draw("E.rel[2]:(E.rel_MC[2] - E.rel[2])" */ + /* ">>h44(80,-0.6,0.6,40,0,6)","hit.both",""); */ + /* c4->cd(8); */ + /* sim->Draw("E.rel_MC[2] - E.rel[2]" */ + /* ">>h48(100,-0.6,0.6)","hit.both"); */ + + + /* // } //return; */ + + + /* // { */ + /* TCanvas *c5=new TCanvas("c5","cm",cn_x0+4*cn_dx,cn_y0+4*cn_dy,1600,800); */ + /* c5->Divide(4,2); */ + + /* // E.cm */ + /* c5->cd(1); */ + /* sim->Draw(Form("E.cm-%f:(E.cm_MC - E.cm)" */ + /* ">>h51(80,-0.3,0.3,40,%f,%f)",q2[0],Ecm_min,Ecm_max), */ + /* "hit.both",""); */ + /* c5->cd(2); */ + /* sim->SetLineColor(kBlack); */ + /* sim->Draw("E.cm_MC - E.cm" */ + /* ">>h55(100,-0.3,0.3)","hit.both"); */ + /* gStyle->SetOptFit(1); */ + /* /\* gStyle->SetStatX(0.5); *\/ */ + /* /\* gStyle->SetStatY(0.7); *\/ */ + /* h55->Fit("gaus"); */ + + /* c5->cd(3); */ + /* sim->SetLineColor(kBlack); */ + /* gStyle->SetOptFit(1); */ + /* /\* gStyle->SetStatX(0.5); *\/ */ + /* /\* gStyle->SetStatY(0.7); *\/ */ + /* sim->Draw("E.cm_MC - E.cm>>h54_3(100,-0.3,0.3)", */ + /* Form("hit.both&&E.cm-%f>0.9&&E.cm-%f<1.1",q2[0],q2[0]), */ + /* ""); */ + /* h54_3->Fit("gaus"); */ + + + + /* /\* sim->SetLineColor(kBlack); *\/ */ + /* /\* sim->Draw("E.cm_MC - E.cm" *\/ */ + /* /\* ">>h55(100,-0.3,0.3)","hit.both"); *\/ */ + /* //h55->SetStats(0); */ + /* /\* sim->SetLineColor(color50); *\/ */ + /* /\* sim->Draw("E.cm_MC - E.cm>>h54_2",pscut50*"hit.both","same"); *\/ */ + /* /\* h54_2->Fit("gaus"); *\/ */ + /* /\* gStyle->SetOptFit(1); *\/ */ + /* /\* h54_2->Draw("same"); *\/ */ + /* //h54_3->Draw(""); */ + + + + /* // th.cm */ + /* /\* c5->cd(2); *\/ */ + /* /\* sim->Draw("th.cm:(th.cm_MC - th.cm)" *\/ */ + /* /\* ">>h52(80,-10,10,40,0,180)","hit.both"); *\/ */ + /* /\* c5->cd(6); *\/ */ + /* /\* sim->Draw("th.cm_MC - th.cm" *\/ */ + /* /\* ">>h56(100,-10,10)","hit.both"); *\/ */ + /* /\* sim->SetLineColor(color50); *\/ */ + /* /\* sim->Draw("th.cm_MC - th.cm",pscut50*"hit.both","same"); *\/ */ + /* /\* sim->SetLineColor(kBlack); *\/ */ + + /* /\* // Q value *\/ */ + /* /\* c5->cd(3); *\/ */ + /* /\* sim->Draw(Form("E.Q3_mid>>h53(80,%f,%f)",Q3-1,Q3+1)); *\/ */ + /* /\* c5->cd(7); *\/ */ + /* /\* sim->Draw(Form("E.Q3_dec>>h57(80,%f,%f)",Q3-1,Q3+1)); *\/ */ + + /* /\* // Q value *\/ */ + /* /\* c5->cd(4); *\/ */ + /* /\* sim->Draw("p.spc:(p.spc_MC - p.spc)" *\/ */ + /* /\* ">>h54(80,-100,100,80,-20-,200)","hit.both"); *\/ */ + /* /\* c5->cd(8); *\/ */ + /* /\* sim->Draw("p.spc_MC - p.spc" *\/ */ + /* /\* ">>h58(100,-100,100)","hit.both"); *\/ */ + + + /* // } //return; */ + + /* TString fname; */ + + /* fname = Form("figs/%s_%s_kin.png",label,simnum); */ + /* c1->Update(); */ + /* c1->Print(fname.Data()); */ + + /* fname = Form("figs/%s_%s_lab.png",label,simnum); */ + /* c2->Update(); */ + /* c2->Print(fname.Data()); */ + + /* fname = Form("figs/%s_%s_pos.png",label,simnum); */ + /* c3->Update(); */ + /* c3->Print(fname.Data()); */ + + /* fname = Form("figs/%s_%s_rel.png",label,simnum); */ + /* c4->Update(); */ + /* c4->Print(fname.Data()); */ + + /* fname = Form("figs/%s_%s_cm.png",label,simnum); */ + /* c5->Update(); */ + /* c5->Print(fname.Data()); */ + +} diff --git a/f17do14an/macro/cs2yield_ecm_thcm.C b/f17do14an/macro/cs2yield_ecm_thcm.C new file mode 100644 index 0000000000000000000000000000000000000000..e01223a406ec5ed9dfa8178b724d0b2246d58a91 --- /dev/null +++ b/f17do14an/macro/cs2yield_ecm_thcm.C @@ -0,0 +1,303 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* #include */ +/* #include */ +/* #include */ +#include +#include +#include + +using namespace std; + +#include "Math/Interpolator.h" +//#include "Math/Polynomial.h" +#include "Math/InterpolationTypes.h" + +//#include "../prm_reac.C" // parameter file <-- +#include "../prm/prm_reac.si26ap_sim001.C" // parameter file <-- + +TFile *f = new TFile(Form("../root/%s_%s.root",label,simnum)); +TTree *tr = (TTree*)f->Get("sim"); + +TF1 *fn1 = new TF1("fn1","gaus",-100,100); + +extern "C" { + /* void enewzsub_(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_(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_(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); */ +} + +// effective thickness // +Double_t eff_thk(Double_t Ecm, Double_t dEcm){ // beam energy in GeV + for(Int_t ichar=0;ichar<33;ichar++){ + if(he[ichar]=='\0'){ + he[ichar]=' '; + he[ichar+1]='\0'; + } + } + float xpu; + float range1, range2, drange; + xpu = (Ecm + dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,he,&range1); + xpu = (Ecm - dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,he,&range2); + drange = fabs(range1 - range2); + + return(drange); +} + + + +int main(){ + /////////////////////////////////////////////////// + // Parameters + //Bool_t inv_thcm = true; + Bool_t inv_thcm = false; + + const Double_t pi = TMath::Pi(); + + // No. of injected beam particles + const Double_t Ibeam = 2.e4; // pps + //const Double_t Dbeam = 5.; // days + const Double_t Dbeam = 1.; // days + //const Double_t nbeam = 3.e10; // <=> 5e4 pps, 7 days + //const Double_t nbeam = Ibeam*Dbeam*24.*60.*60.; + const Double_t nbeam = 1231020000; + // No. of injected beam particles + const Double_t Nav = 6.0221e23; + const Double_t C_cm2_mb_mg = Nav/mtrg*1.e-30; + + // Bin deffinition + //const Int_t necm = 140; // Ebin = 50 keV + const Int_t necm = 70; // Ebin = 100 keV + const Double_t ecmmin = 0.; + const Double_t ecmmax = 7.; + const Double_t decm = (ecmmax - ecmmin)/double(necm); + + const Int_t nthcm = 90; // theta bin = 2 deg. + const Double_t thcmmin = 0.; + const Double_t thcmmax = 180.; + const Double_t dthcm = (thcmmax - thcmmin)/double(nthcm); + const Double_t d2r = pi/180.; + + // Cut condition + const char* cut = + //"hit.idet==0&&iex==0" + //"hit.idet>=0&&hit.idet<=2&&iex==0" + "hit.any&&iex==0" + //"&&pos.ztar>-236."; // to avoid events near the entrance + "&&pos.ztar>-266."; // to avoid events near the entrance + // "&&(E.lab[0] < 32. || E.lab[0]>=32.&&(E.cm > 0.19*E.lab[0]+0.22))"; + + // cross section data file + const Char_t fcs[100] = "../dat/nonsmoker_cs.txt"; // NON-SMOKER calculation + + const Double_t lbsize = 0.05; + //const Double_t txsize = 0.06; + /////////////////////////////////////////////////// + + + TApplication theApp("App", NULL, NULL); + gStyle->SetTextFont(62); + gStyle->SetTextSize(lbsize); + + gStyle->SetLabelFont(62,"XYZt"); + gStyle->SetLabelSize(lbsize,"XYZt"); + + //gStyle->SetTitleFont(62); + gStyle->SetTitleFont(62,"XYZt"); + gStyle->SetTitleSize(lbsize,"XYZt"); + gStyle->SetTitleFontSize(lbsize); + gStyle->SetTitleYOffset(1.2); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetPadBottomMargin(0.1); + + // gStyle->SetStatW(0.3); + // gStyle->SetStatH(0.2); + // gStyle->SetStatX(1); + // gStyle->SetStatY(0.4); + gStyle->SetOptStat(0); + + // Canvas + TCanvas *c1=new TCanvas("c1","c1",0,0,880,800); + c1->Divide(3,3); + //c1->Divide(2,2); + c1->cd(1); + + Char_t hname[100]; + + // TH2 for Nhit (Ecm vs thcm) + sprintf(hname,"hhit"); + TH2D *hhit = new TH2D(hname,"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + if(inv_thcm) + tr->Draw(Form("180. - th.cm:E.cm>>%s",hname),cut); + else + tr->Draw(Form("th.cm:E.cm>>%s",hname),cut); + + // TH2 for Ntot (Ecm vs thcm) + sprintf(hname,"htot"); + TH2D *htot = new TH2D(hname,"Total events;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + if(inv_thcm) + tr->Draw(Form("180. - th.cm:E.cm>>%s",hname),"iex==0"); + else + tr->Draw(Form("th.cm:E.cm>>%s",hname),"iex==0"); + + sprintf(hname,"hsin_int"); + TH2D *hsin_int = new TH2D(hname, + "Integrated sin(#theta) in #theta bin size;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // dO: TH2 Divide (Ecm vs thcm) + TH2D *hdO = new TH2D("hdO","Solid angle (msr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO->Divide(hhit,htot); //,4.*pi); + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thcm) + TH2D *hdsdO = new TH2D("hdsdO","Cross section (mb/sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // TH2 Fill: thk/dE*dE: Eff. thick. vs Ecm + TH2D *hthk = new TH2D("hthk",";" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // Reading cross section data + const Int_t nhf = 90; // no. of data + Double_t ehf[nhf], notused, cshf[nhf]; + ifstream datahf(fcs); + for(Int_t i=0; i> ehf[i] >> notused >> notused >> cshf[i]; + cshf[i] *= 1000./4./pi; // b -> mb/sr + //cout << i << " " << ehf[i] <<" "<< cshf[i] << endl; + } + datahf.close(); + ROOT::Math::Interpolator interhf(nhf, ROOT::Math::Interpolation::kLINEAR); + interhf.SetData(nhf, ehf, cshf); // MeV, mb + + Double_t ecm = ecmmin + decm/2.; + for(Int_t i=0; iFill(ecm,thcm,w_cshf); + hthk->Fill(ecm,thcm,w_efth); + hsin_int->Fill( ecm,thcm, - ( cos(d2r*(thcm+dthcm/2.)) + - cos(d2r*(thcm-dthcm/2.)) ) ); + thcm += dthcm; + } + ecm += decm; + } + + hdO->Multiply(hdO,hsin_int,2.*pi*1000.); + + /////////////////////////////////////////////////// + // Draw histograms + c1->cd(1); + hhit->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(2); + htot->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(3); + hsin_int->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(4); + gPad->SetLogz(0); + hdO->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(5); + gPad->SetLogz(1); + hdsdO->GetZaxis()->SetRangeUser(1.e-8,1.e2); + hdsdO->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(6); + gPad->SetLogz(0); + hthk->SetTitle(Form("Effective thickness / %5.3f MeV (mg/cm^2)",decm)); + hthk->SetTitleFont(62); + hthk->SetTitleSize(lbsize); + hthk->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(8); + /* + Yiled + = Nbeam * ds/dO(mb/sr) * dO(sr) * thk(mb/cm^2) * const(mg*mb/cm^2) + */ + TH2D *hyld = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + gPad->SetLogz(1); + hyld->Multiply(hdO,hdsdO,1/1000.); + hyld->Multiply(hyld,hthk,nbeam*C_cm2_mb_mg); + hyld->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); + hyld->SetTitleFont(62); + hyld->SetTitleSize(lbsize); + hyld->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(7); + TH1D *hdO_ecm = hdO->ProjectionX("hdO_ecm",0,nthcm); + gPad->SetLogy(1); + hdO_ecm->GetYaxis()->SetRangeUser(10.,2000.); + hdO_ecm->Draw("HIST"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(9); + TH1D *hyld_ecm = hyld->ProjectionX("hyld_ecm",0,nthcm); + gPad->SetLogy(1); + hyld_ecm->GetYaxis()->SetRangeUser(1.,1000.); + hyld_ecm->SetTitle(Form("Yield / %5.3f MeV",decm)); + hyld_ecm->SetTitleFont(62); + hyld_ecm->SetTitleSize(lbsize); + hyld_ecm->Draw("HIST"); + Double_t ntot_yld = hyld_ecm->Integral(0,necm); + //cout << hyld_ecm->Integral(0,necm) << endl; + Double_t rate_yld = ntot_yld/Dbeam/24.; // counts/hour + TLatex *t1=new TLatex(0.2,400,Form("Total: %.1f counts",ntot_yld)); + TLatex *t2=new TLatex(0.2,200,Form("Rate: %.1f counts/h",rate_yld)); + t1->Draw(); + t2->Draw(); + c1->Update(); gSystem->ProcessEvents(); + + gSystem->Exec("mkdir -vp ../figs");//Creates figs/ directory if there is not + c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_telall.png",label,simnum)); + c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_telall.pdf",label,simnum)); + + theApp.Run(); +} diff --git a/f17do14an/macro/cs2yield_ecm_thcm_tel.C b/f17do14an/macro/cs2yield_ecm_thcm_tel.C new file mode 100644 index 0000000000000000000000000000000000000000..738724993973fbb7697c22b60796e744a5a860aa --- /dev/null +++ b/f17do14an/macro/cs2yield_ecm_thcm_tel.C @@ -0,0 +1,312 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* #include */ +/* #include */ +/* #include */ +#include +#include +#include + +using namespace std; + +#include "Math/Interpolator.h" +//#include "Math/Polynomial.h" +#include "Math/InterpolationTypes.h" + +//#include "../prm_reac.C" // parameter file <-- +#include "../prm/prm_reac.si26ap_sim001.C" // parameter file <-- + +TFile *f = new TFile(Form("../root/%s_%s.root",label,simnum)); +TTree *tr = (TTree*)f->Get("sim"); + +TF1 *fn1 = new TF1("fn1","gaus",-100,100); + +extern "C" { + /* void enewzsub_(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_(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_(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); */ +} + +// effective thickness // +Double_t eff_thk(Double_t Ecm, Double_t dEcm){ // beam energy in GeV + for(Int_t ichar=0;ichar<33;ichar++){ + if(he[ichar]=='\0'){ + he[ichar]=' '; + he[ichar+1]='\0'; + } + } + float xpu; + float range1, range2, drange; + xpu = (Ecm + dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,he,&range1); + xpu = (Ecm - dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,he,&range2); + drange = fabs(range1 - range2); + + return(drange); +} + + + +int main(){ + /////////////////////////////////////////////////// + // Parameters + //Bool_t inv_thcm = true; + Bool_t inv_thcm = false; + + const Double_t pi = TMath::Pi(); + + // No. of injected beam particles + //const Double_t Ibeam = 2.e4; // pps + //const Double_t Dbeam = 5.; // days + const Double_t Dbeam = 1.; // days + //const Double_t nbeam = 3.e10; // <=> 5e4 pps, 7 days + //const Double_t nbeam = Ibeam*Dbeam*24.*60.*60.; + const Double_t nbeam = 1231020000; + // No. of injected beam particles + const Double_t Nav = 6.0221e23; + const Double_t C_cm2_mb_mg = Nav/mtrg*1.e-30; + + // Bin deffinition + //const Int_t necm = 140; // Ebin = 50 keV + const Int_t necm = 80; // Ebin = 100 keV + const Double_t ecmmin = 0.; + const Double_t ecmmax = 8.; + const Double_t decm = (ecmmax - ecmmin)/double(necm); + + const Int_t nthcm = 90; // theta bin = 2 deg. + const Double_t thcmmin = 0.; + const Double_t thcmmax = 180.; + const Double_t dthcm = (thcmmax - thcmmin)/double(nthcm); + const Double_t d2r = pi/180.; + + // Cut condition + const Int_t idet = 4; + const char* cut = Form( + "hit.idet==%d&&iex==0" + // "hit.any&&iex==0" + "&&pos.ztar>-319." + ,idet + ); // to avoid events near the entrance + // "&&(E.lab[0] < 32. || E.lab[0]>=32.&&(E.cm > 0.19*E.lab[0]+0.22))"; + + // cross section data file + const Char_t fcs[100] = "../dat/nonsmoker_cs.txt"; // NON-SMOKER calculation + + const Double_t lbsize = 0.05; + //const Double_t txsize = 0.06; + /////////////////////////////////////////////////// + + + TApplication theApp("App", NULL, NULL); + gStyle->SetTextFont(62); + gStyle->SetTextSize(lbsize); + + gStyle->SetLabelFont(62,"XYZt"); + gStyle->SetLabelSize(lbsize,"XYZt"); + + //gStyle->SetTitleFont(62); + gStyle->SetTitleFont(62,"XYZt"); + gStyle->SetTitleSize(lbsize,"XYZt"); + gStyle->SetTitleFontSize(lbsize); + gStyle->SetTitleYOffset(1.2); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetPadBottomMargin(0.1); + + // gStyle->SetStatW(0.3); + // gStyle->SetStatH(0.2); + // gStyle->SetStatX(1); + // gStyle->SetStatY(0.4); + gStyle->SetOptStat(0); + + // Canvas + TCanvas *c1=new TCanvas("c1","c1",0,0,880,800); + c1->Divide(3,3); + //c1->Divide(2,2); + c1->cd(1); + + Char_t hname[100]; + + // TH2 for Nhit (Ecm vs thcm) + sprintf(hname,"hhit"); + TH2D *hhit = new TH2D(hname,"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + if(inv_thcm) + tr->Draw(Form("180. - th.cm:E.cm>>%s",hname),cut); + else + tr->Draw(Form("th.cm:E.cm>>%s",hname),cut); + + // TH2 for Ntot (Ecm vs thcm) + sprintf(hname,"htot"); + TH2D *htot = new TH2D(hname,"Total events;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + if(inv_thcm) + tr->Draw(Form("180. - th.cm:E.cm>>%s",hname),"iex==0"); + else + tr->Draw(Form("th.cm:E.cm>>%s",hname),"iex==0"); + + sprintf(hname,"hsin_int"); + TH2D *hsin_int = new TH2D(hname, + "Integrated sin(#theta) in #theta bin size;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // dO: TH2 Divide (Ecm vs thcm) + TH2D *hdO = new TH2D("hdO","Solid angle (msr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO->Divide(hhit,htot); //,4.*pi); + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thcm) + TH2D *hdsdO = new TH2D("hdsdO","Cross section (mb/sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // TH2 Fill: thk/dE*dE: Eff. thick. vs Ecm + TH2D *hthk = new TH2D("hthk",";" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // Reading cross section data + const Int_t nhf = 90; // no. of data + Double_t ehf[nhf], notused, cshf[nhf]; + ifstream datahf(fcs); + for(Int_t i=0; i> ehf[i] >> notused >> notused >> cshf[i]; + cshf[i] *= 1000./4./pi; // b -> mb/sr + //cout << i << " " << ehf[i] <<" "<< cshf[i] << endl; + } + datahf.close(); + ROOT::Math::Interpolator interhf(nhf, ROOT::Math::Interpolation::kLINEAR); + interhf.SetData(nhf, ehf, cshf); // MeV, mb + + Double_t ecm = ecmmin + decm/2.; + for(Int_t i=0; iFill(ecm,thcm,w_cshf); + hthk->Fill(ecm,thcm,w_efth); + hsin_int->Fill( ecm,thcm, - ( cos(d2r*(thcm+dthcm/2.)) + - cos(d2r*(thcm-dthcm/2.)) ) ); + thcm += dthcm; + } + ecm += decm; + } + + hdO->Multiply(hdO,hsin_int,2.*pi*1000.); + + /////////////////////////////////////////////////// + // Draw histograms + c1->SetFillStyle(0); + + + c1->cd(1); + gPad->SetFillStyle(0); + hhit->SetFillStyle(0); + hhit->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(2); + htot->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(3); + hsin_int->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(4); + gPad->SetLogz(0); + hdO->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(5); + gPad->SetLogz(1); + hdsdO->GetZaxis()->SetRangeUser(1.e-8,1.e2); + hdsdO->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(6); + gPad->SetLogz(0); + hthk->SetTitle(Form("Effective thickness / %5.3f MeV (mg/cm^2)",decm)); + hthk->SetTitleFont(62); + hthk->SetTitleSize(lbsize); + hthk->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(8); + /* + Yiled + = Nbeam * ds/dO(mb/sr) * dO(sr) * thk(mb/cm^2) * const(mg*mb/cm^2) + */ + TH2D *hyld = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + gPad->SetLogz(1); + hyld->Multiply(hdO,hdsdO,1/1000.); + hyld->Multiply(hyld,hthk,nbeam*C_cm2_mb_mg); + hyld->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); + hyld->SetTitleFont(62); + hyld->SetTitleSize(lbsize); + hyld->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(7); + TH1D *hdO_ecm = hdO->ProjectionX("hdO_ecm",0,nthcm); + gPad->SetLogy(1); + hdO_ecm->GetYaxis()->SetRangeUser(10.,2000.); + hdO_ecm->Draw("HIST"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(9); + TH1D *hyld_ecm = hyld->ProjectionX("hyld_ecm",0,nthcm); + Int_t ymax9 = 100.; + gPad->SetLogy(1); + hyld_ecm->GetYaxis()->SetRangeUser(1.,ymax9); + hyld_ecm->SetTitle(Form("Yield / %5.3f MeV",decm)); + hyld_ecm->SetTitleFont(62); + hyld_ecm->SetTitleSize(lbsize); + hyld_ecm->Draw("HIST"); + Double_t ntot_yld = hyld_ecm->Integral(0,necm); + //cout << hyld_ecm->Integral(0,necm) << endl; + Double_t rate_yld = ntot_yld/Dbeam/24.; // counts/hour + TLatex *t1=new TLatex(0.2,ymax9*0.4,Form("Total: %.1f counts",ntot_yld)); + TLatex *t2=new TLatex(0.2,ymax9*0.2,Form("Rate: %.1f counts/h",rate_yld)); + t1->Draw(); + t2->Draw(); + c1->Update(); gSystem->ProcessEvents(); + + gSystem->Exec("mkdir -vp ../figs");//Creates figs/ directory if there is not + c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_24h_tel%d.png",label,simnum,idet+1)); + c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_24h_tel%d.pdf",label,simnum,idet+1)); + /* c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_24h_telall.png",label,simnum)); */ + /* c1->Print(Form("../figs/yield_ecm-thcm_%s_%s_24h_telall.pdf",label,simnum)); */ + + theApp.Run(); +} diff --git a/f17do14an/macro/plot_yield_coin b/f17do14an/macro/plot_yield_coin new file mode 100755 index 0000000000000000000000000000000000000000..d4494a898dfca57f12b6997553aa878205cf7d49 Binary files /dev/null and b/f17do14an/macro/plot_yield_coin differ diff --git a/f17do14an/macro/plot_yield_coin.C b/f17do14an/macro/plot_yield_coin.C new file mode 100644 index 0000000000000000000000000000000000000000..0c5ca9aaa09f97b27da2ca73a1525a5c2beb24f6 --- /dev/null +++ b/f17do14an/macro/plot_yield_coin.C @@ -0,0 +1,434 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* #include */ +/* #include */ +/* #include */ +#include +#include +#include + +using namespace std; + +#include "Math/Interpolator.h" +//#include "Math/Polynomial.h" +#include "Math/InterpolationTypes.h" + +#include "../prm/prm_reac.f17do14an_sim001.C" // parameter file <-- + +TFile *f = new TFile(Form("../root/%s_%s.root",label,simnum)); +TTree *tr = (TTree*)f->Get("sim"); +TF1 *fn1 = new TF1("fn1","gaus",-100,100); + +extern "C" { + /* void enewzsub_(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_(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_(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); */ +} + +// effective thickness // +Double_t eff_thk(Double_t Ecm, Double_t dEcm){ // beam energy in GeV + for(Int_t ichar=0;ichar<33;ichar++){ + if(cd2[ichar]=='\0'){ + cd2[ichar]=' '; + cd2[ichar+1]='\0'; + } + } + float xpu; + float range1, range2, drange; + xpu = (Ecm + dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,cd2,&range1); + xpu = (Ecm - dEcm/2.)*(mprj + mtrg)/mprj/mtrg; + rangezsub_(&zprj,&mprj,&xpu,cd2,&range2); + drange = fabs(range1 - range2); + + return(drange); +} + + + +int main(){ + const Double_t pi = TMath::Pi(); + const Double_t lbsize = 0.05; + //const Double_t txsize = 0.04; + + TApplication theApp("App", NULL, NULL); + gStyle->SetTextFont(62); + gStyle->SetTextSize(lbsize); + + gStyle->SetLabelFont(62,"XYZt"); + gStyle->SetLabelSize(lbsize,"XYZt"); + + //gStyle->SetTitleFont(62); + gStyle->SetTitleFont(62,"XYZt"); + gStyle->SetTitleSize(lbsize,"XYZt"); + gStyle->SetTitleFontSize(lbsize); + gStyle->SetTitleYOffset(1.2); + + gStyle->SetPadTopMargin(0.06); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetPadBottomMargin(0.1); + + // gStyle->SetStatW(0.3); + // gStyle->SetStatH(0.2); + // gStyle->SetStatX(1); + // gStyle->SetStatY(0.4); + gStyle->SetOptStat(0); + + // Reading cross section data + /* const Int_t ndata = 181; // no. of data */ + /* Double_t thdata[ndata], csdata[ndata]; */ + /* //ifstream datahf("../../../nonsmoker_cs.txt"); */ + /* ifstream fdata("cs_he6pt_green.txt"); */ + /* for(Int_t i=0; i> thdata[i] >> csdata[i]; */ + /* //cshdata[i] *= 1000./4./pi; // b -> mb/sr */ + /* //cout << i << " " << ehf[i] <<" "<< cshf[i] << endl; */ + /* } */ + /* fdata.close(); */ + /* ROOT::Math::Interpolator intercs(ndata, ROOT::Math::Interpolation::kLINEAR); */ + /* intercs.SetData(ndata, thdata, csdata); // MeV, mb */ + Double_t w_cs = 19.2/16./pi/pi; // constant cs. 19.2 mb @Ef17 = 120 MeV + + // No. of injected beam particles + const Double_t Ibeam = 5.5e5; // pps + const Double_t Dbeam = 5.; // days + //const Double_t nbeam = 3.e10; // <=> 5e4 pps, 7 days + const Double_t nbeam = Ibeam*Dbeam*24.*60.*60.; + + // No. of injected beam particles + const Double_t Nav = 6.0221e23; // mol^-1 + //const Double_t C_cm2_mb_mg = Nav/mtrg*1.e-30; // molecules/mg * cm^2/mb + const Double_t C_cm2_mb_mg = Nav/(12. + 2.*mtrg)*2.*1.e-30; // CD2 + //cout << C_cm2_mb_mg << endl; + + // Bin deffinition + const Int_t necm = 200; // Ebin = 50 keV + //const Int_t necm = 70; // Ebin = 100 keV + + const Double_t csmin = 0.01; + const Double_t csmax = 1.; + + const Double_t ecmmin = 0.; + const Double_t ecmmax = 10.; + const Double_t decm = (ecmmax - ecmmin)/double(necm); + + const Int_t nthcm = 45; + const Double_t thcmmin = 0.; + const Double_t thcmmax = 180.; + + const Double_t samin = 0.; + const Double_t samax = .001; + + const Double_t dthcm = (thcmmax - thcmmin)/double(nthcm); + const Double_t d2r = pi/180.; + + // Cut condition + const char* cut_a_all = + "hit.any[0]&&iex==0"; + const char* cut_o1a4 = + "hit.idet[0]==0&&hit.idet[1]==3&&iex==0&&abs(p.spc)<50*w.evt"; + const char* cut_t_all = + "hit.any[1]&&iex==0"; + const char* cut_o2a3 = + "hit.idet[0]==1&&hit.idet[1]==2&&iex==0&&abs(p.spc)<50*w.evt"; + //Char_t hname[100]; + TString hname; + + // TH2 for Nhit (Ecm vs thcm) + //sprintf(hname,"hhit_a_all"); + hname = "hhit_a_all"; + TH2D *hhit_a_all = new TH2D(hname.Data(),"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + tr->Draw(Form("th.cm:E.cm>>%s",hname.Data()),cut_a_all); + /* TH2D *hhit_a[ndet]; */ + /* for(Int_t i=0; iDraw(Form("th.cm:E.cm>>%s",hname),Form("hit.idet[0]==%d&&iex==0",i)); */ + /* } */ + //sprintf(hname,"hhit_o1a4"); + hname = "hhit_o1a4"; + TH2D *hhit_o1a4 = new TH2D(hname.Data(),"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + tr->Draw(Form("th.cm:E.cm>>%s",hname.Data()),cut_o1a4); + //sprintf(hname,"hhit_o2a3"); + hname = "hhit_o2a3"; + TH2D *hhit_o2a3 = new TH2D(hname.Data(),"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + tr->Draw(Form("th.cm:E.cm>>%s",hname.Data()),cut_o2a3); + + + //sprintf(hname,"hhit_t_all"); + hname = "hhit_t_all"; + TH2D *hhit_t_all = new TH2D(hname.Data(),"Events hitting detectors;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + tr->Draw(Form("th.cm:E.cm>>%s",hname.Data()),cut_t_all); + /* TH2D *hhit_t[ndet]; */ + /* for(Int_t i=0; iDraw(Form("th.cm:E.cm>>%s",hname),Form("hit.idet[1]==%d&&iex==0",i)); */ + /* } */ + + + // TH2 for Ntot (Ecm vs thcm) + //sprintf(hname,"htot"); + hname = "htot"; + TH2D *htot = new TH2D(hname.Data(),"Total events;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + tr->Draw(Form("th.cm:E.cm>>%s",hname.Data()),"iex==0"); + + //sprintf(hname,"hsin_int"); + hname = "hsin_int"; + TH2D *hsin_int = new TH2D(hname.Data(), + "Integrated sin(#theta) in #theta bin size;" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // dO: TH2 Divide (Ecm vs thcm) + TH2D *hdO_a_all = new TH2D("hdO_a_all","Solid angle #alpha (sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO_a_all->Divide(hhit_a_all,htot); //,4.*pi); + /* TH2D *hdO_a[ndet]; */ + /* for(Int_t i=0; iDivide(hhit_a[i],htot); */ + /* } */ + TH2D *hdO_o1a4 = new TH2D("hdO_o1a4","Solid angle #alpha (sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO_o1a4->Divide(hhit_o1a4,htot); + TH2D *hdO_o2a3 = new TH2D("hdO_o2a3","Solid angle #alpha (sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO_o2a3->Divide(hhit_o2a3,htot); + + + TH2D *hdO_t_all = new TH2D("hdO_t_all","Solid angle t (sr);" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO_t_all->Divide(hhit_t_all,htot); //,4.*pi); + /* TH2D *hdO_t[ndet]; */ + /* for(Int_t i=0; iDivide(hhit_t[i],htot); */ + /* } */ + + /* hdO->Draw("colz"); */ + /* c1->Update(); gSystem->ProcessEvents(); getchar(); */ + + // TH2 Fill: ds/dO: Diff. cross section vs Ecm (vs thcm) + TH2D *hdsdO = new TH2D("hdsdO","Cross section (mb/sr^{2});" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + // TH2 Fill: thk/dE*dE: Eff. thick. vs Ecm + TH2D *hthk = new TH2D("hthk",";" + "E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + + Double_t ecm = ecmmin + decm/2.; + for(Int_t i=0; iFill(ecm,thcm,w_cs); + hthk->Fill(ecm,thcm,w_efth); + hsin_int->Fill( ecm,thcm, - ( cos(d2r*(thcm+dthcm/2.)) + - cos(d2r*(thcm-dthcm/2.)) ) ); + thcm += dthcm; + } + ecm += decm; + } + + hdO_a_all->Multiply(hdO_a_all,hsin_int,2.*pi); + hdO_t_all->Multiply(hdO_t_all,hsin_int,2.*pi); + /* for(Int_t i=0; iMultiply(hdO_a[i],hsin_int,2.*pi); */ + /* hdO_t[i]->Multiply(hdO_t[i],hsin_int,2.*pi); */ + /* } */ + hdO_o1a4->Multiply(hdO_o1a4,hsin_int,2.*pi); + hdO_o2a3->Multiply(hdO_o2a3,hsin_int,2.*pi); + + // Canvas + TCanvas *c1=new TCanvas("c1","c1",1200,0,1320,600); + c1->Divide(4,2); + + c1->cd(1); + gPad->SetLogz(1); + hdsdO->GetXaxis()->SetNdivisions(502); + hdsdO->GetZaxis()->SetRangeUser(csmin,csmax); + hdsdO->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(2); + gPad->SetLogz(0); + hthk->GetXaxis()->SetNdivisions(502); + hthk->SetTitle(Form("Effective thickness / %5.3f MeV (mg/cm^2)",decm)); + hthk->SetTitleFont(62); + hthk->SetTitleSize(lbsize); + hthk->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(3); + TH2D *hdO_all_coin = new TH2D("hdO_all_coin",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + hdO_all_coin->Add(hdO_all_coin,hdO_o1a4); + hdO_all_coin->Add(hdO_all_coin,hdO_o2a3); + hdO_all_coin->Draw("colz"); + + c1->cd(4); + gPad->SetLogy(0); + TH1D *hdO_all_coin_thcm = hdO_all_coin->ProjectionY("hdO_all_coin_thcm",0,necm); + hdO_all_coin_thcm->Scale(1./necm); + hdO_all_coin_thcm->GetYaxis()->SetRangeUser(samin,samax); + hdO_all_coin_thcm->GetXaxis()->SetTitleOffset(0.96); + hdO_all_coin_thcm->SetLineColor(kBlack); + hdO_all_coin_thcm->SetLineStyle(2); + hdO_all_coin_thcm->Draw("HIST"); + TH1D *hdO_o1a4_thcm = hdO_o1a4->ProjectionY("hdO_o1a4_thcm",0,necm); + hdO_o1a4_thcm->Scale(1./necm); + hdO_o1a4_thcm->SetLineColor(kRed); + //hdO_o1a4_thcm->SetLineStyle(2); + hdO_o1a4_thcm->Draw("HISTsame"); + TH1D *hdO_o2a3_thcm = hdO_o2a3->ProjectionY("hdO_o2a3_thcm",0,necm); + hdO_o2a3_thcm->Scale(1./necm); + hdO_o2a3_thcm->SetLineColor(kBlue); + //hdO_o2a3_thcm->SetLineStyle(2); + hdO_o2a3_thcm->Draw("HISTsame"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + c1->cd(5); + TH2D *hyld_o1a4 = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + gPad->SetLogz(1); + hyld_o1a4->Multiply(hdO_o1a4,hdsdO); + hyld_o1a4->Multiply(hyld_o1a4,hthk,nbeam*C_cm2_mb_mg); + hyld_o1a4->GetXaxis()->SetNdivisions(502); + hyld_o1a4->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); + hyld_o1a4->SetTitleFont(62); + hyld_o1a4->SetTitleSize(lbsize); + hyld_o1a4->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + /* c1->cd(6); */ + TH2D *hyld_o2a3 = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", + necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); + gPad->SetLogz(1); + hyld_o2a3->Multiply(hdO_o2a3,hdsdO); + hyld_o2a3->Multiply(hyld_o2a3,hthk,nbeam*C_cm2_mb_mg); + hyld_o2a3->GetXaxis()->SetNdivisions(502); + hyld_o2a3->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); + hyld_o2a3->SetTitleFont(62); + hyld_o2a3->SetTitleSize(lbsize); + hyld_o2a3->Draw("colz"); + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + /* c1->cd(6); */ + /* TH2D *hyld_a2t1 = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", */ + /* necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); */ + /* gPad->SetLogz(1); */ + /* hyld_a2t1->Multiply(hdO_a2t1,hdsdO); */ + /* hyld_a2t1->Multiply(hyld_a2t1,hthk,nbeam*C_cm2_mb_mg); */ + /* hyld_a2t1->GetXaxis()->SetNdivisions(502); */ + /* hyld_a2t1->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); */ + /* hyld_a2t1->SetTitleFont(62); */ + /* hyld_a2t1->SetTitleSize(lbsize); */ + /* //hyld_a2t1->Draw("colsame"); */ + /* hyld_a2t1->Draw("colz"); */ + /* c1->Update(); gSystem->ProcessEvents(); //getchar(); */ + + /* /\* c1->cd(6); *\/ */ + /* TH2D *hyld_a2t3 = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", */ + /* necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); */ + /* /\* gPad->SetLogz(1); *\/ */ + /* /\* hyld_a2t3->Multiply(hdO_a2t3,hdsdO); *\/ */ + /* /\* hyld_a2t3->Multiply(hyld_a2t3,hthk,nbeam*C_cm2_mb_mg); *\/ */ + /* /\* hyld_a2t3->GetXaxis()->SetNdivisions(502); *\/ */ + /* /\* hyld_a2t3->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); *\/ */ + /* /\* hyld_a2t3->SetTitleFont(62); *\/ */ + /* /\* hyld_a2t3->SetTitleSize(lbsize); *\/ */ + /* /\* hyld_a2t3->Draw("colz"); *\/ */ + /* /\* c1->Update(); gSystem->ProcessEvents(); //getchar(); *\/ */ + + /* /\* c1->cd(7); *\/ */ + /* TH2D *hyld_a3t2 = new TH2D("hyld",";E_{c.m.} (MeV);#theta_{c.m.} (deg.)", */ + /* necm,ecmmin,ecmmax,nthcm,thcmmin,thcmmax); */ + /* /\* gPad->SetLogz(1); *\/ */ + /* /\* hyld_a3t2->Multiply(hdO_a3t2,hdsdO); *\/ */ + /* /\* hyld_a3t2->Multiply(hyld_a3t2,hthk,nbeam*C_cm2_mb_mg); *\/ */ + /* /\* hyld_a3t2->GetXaxis()->SetNdivisions(502); *\/ */ + /* /\* hyld_a3t2->SetTitle(Form("Yield / %5.3f MeV / %3.1f deg.",decm,dthcm)); *\/ */ + /* /\* hyld_a3t2->SetTitleFont(62); *\/ */ + /* /\* hyld_a3t2->SetTitleSize(lbsize); *\/ */ + /* /\* hyld_a3t2->Draw("colz"); *\/ */ + /* /\* c1->Update(); gSystem->ProcessEvents(); //getchar(); *\/ */ + + c1->cd(8); + TH1D *hyld_o1a4_thcm = hyld_o1a4->ProjectionY("hyld_o1a4_thcm",0,necm); + TH1D *hyld_o2a3_thcm = hyld_o2a3->ProjectionY("hyld_o2a3_thcm",0,necm); + TH1D *hyld_all_thcm = new TH1D("hyld_all_thcm","",nthcm,thcmmin,thcmmax); + hyld_all_thcm->Add(hyld_all_thcm, hyld_o1a4_thcm); + hyld_all_thcm->Add(hyld_all_thcm, hyld_o2a3_thcm); + hyld_all_thcm->SetLineColor(kBlack); + hyld_all_thcm->SetLineStyle(2); + hyld_all_thcm->Draw("HIST"); + + hyld_o1a4_thcm->SetLineColor(kRed); + hyld_o1a4_thcm->Draw("HISTsame"); + hyld_o2a3_thcm->SetLineColor(kBlue); + hyld_o2a3_thcm->Draw("HISTsame"); + Double_t ntot_yld = hyld_all_thcm->Integral(0,nthcm); + //cout << hyld_thcm->Integral(0,nthcm) << endl; + Double_t rate_yld = ntot_yld/Dbeam/24.; // counts/hour + Double_t ymax=hyld_all_thcm->GetMaximum(); + TLatex *t1=new TLatex(80,0.92*ymax,Form("Total: %.1f counts",ntot_yld)); + TLatex *t2=new TLatex(80,0.85*ymax,Form("Rate: %.1f counts/h",rate_yld)); + t1->Draw(); + t2->Draw(); + + c1->Update(); gSystem->ProcessEvents(); //getchar(); + + gSystem->Exec("mkdir -vp ../figs");//Creates figs/ directory if there is not + c1->Print(Form("../figs/sa_coin_%s_%s.png",label,simnum)); + c1->Print(Form("../figs/sa_coin_%s_%s.pdf",label,simnum)); + + + theApp.Run(); +} diff --git a/f17do14an/prm/prm_reac.f17do14an_sim001.C b/f17do14an/prm/prm_reac.f17do14an_sim001.C new file mode 100644 index 0000000000000000000000000000000000000000..5ee4a82f623fcfdb4d4725073cb4c7ef5157f8b2 --- /dev/null +++ b/f17do14an/prm/prm_reac.f17do14an_sim001.C @@ -0,0 +1,139 @@ +/* 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 */ +const Long64_t ntot = 2.e5; // number of events +Char_t dir[50] = "f17do14an"; // +Char_t label[50] = "f17do14an"; // +Char_t simnum[10] = "sim001"; // No beam spread, detector resolutions (in energy, space) +////////////////////////// + + +/* 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 m17f = 17.002095237; +Int_t z9 = 9; +float m14o = 14.008596706; +Int_t z8 = 8; +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 = z9; +float mprj = m17f; +Int_t ztrg = z1; +float mtrg = md; +float mprt = mh; +Int_t zres = z8; +float mres[Nexci], q2[Nexci]; +float mres_gs = m14o; +Int_t zrec = z2; +float mrec = ma; +Int_t zspc = 0; +float mspc = mn; +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 carbon[34] = "c"; +Char_t lif[34] = "lif"; +Char_t he[34] = "he"; +Char_t si[34] = "si"; +//Char_t mat_tar[34], mat_bak[34]; +Char_t mat_tar[34] = "cd2"; +Char_t mat_bak[34] = "c"; +/* strcpy(mat_tar, cd2); */ +/* strcpy(mat_bak, carbon); */ +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.15; // Target thickness in mg/cm^2 // sim001 +//float thk_tar = 0.; // Target thickness in mg/cm^2 // sim001 +float thk_bak = 0.; +float r_sphere = 260.; // 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 +////////////////////////// + + +/* 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 = 120.; // 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 +Double_t ytar0 = 0.; // mm +Double_t xtar0_sgm = 0.; // mm sim001 +Double_t ytar0_sgm = 0.; // mm sim001 +/* Double_t xtar0_sgm = 5.; // sim002 */ +/* Double_t ytar0_sgm = 10.; // sim002 */ +// +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 +Double_t sgm_ppac_b = 0.; // mm, PPAC b intrinsic resolution // sim001 +/* Double_t sgm_ppac_a = 0.5; // mm, PPAC a intrinsic resolution sim002 */ +/* Double_t sgm_ppac_b = 0.5; // mm, PPAC b intrinsic resolution sim002 */ +////////////////////////// + + +/* Detector arrangement */ +const Int_t ndet = 4; +Double_t thdet[ndet] = {0., 0., 0., 0.}; // in deg +Double_t ldet[ndet] = {250., 250.,150., 150.}; +Double_t xdet_off[ndet] = {23., -23., 50., -50.}; +Double_t zdet_off[ndet] = {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_xhalf[ndet] = {12., 12., 22.5, 22.5}; +Double_t psd_yhalf[ndet] = {13.5, 13.5, 22.5, 22.5}; +Int_t n_xstrip[ndet] = {48, 48, 22, 22}; // sim001 +Int_t n_ystrip[ndet] = {3, 3, 22, 22}; // sim001 +/* Int_t n_xstrip[ndet] = {8, 8, 8, 8}; // sim001 */ +/* Int_t n_ystrip[ndet] = {3, 3, 8, 8}; // sim001 */ +//Double_t strip_half[ndet]; +//Double_t dth[ndet]; // 03-04-2018 +////////////////////////// + +//Double_t dt = 1.; // TOF resolution in nsec. + + diff --git a/f17do14an/root/f17do14an_sim001.root b/f17do14an/root/f17do14an_sim001.root new file mode 100644 index 0000000000000000000000000000000000000000..6b570ccbd36182e9ac62aa3b5662015061c09963 Binary files /dev/null and b/f17do14an/root/f17do14an_sim001.root differ diff --git a/prm_reac.C b/prm_reac.C index b16c767f2443639740dd2672daead45e74aa7303..cf37ac686411bc3b8a213cc560ba2354b531f9b6 120000 --- a/prm_reac.C +++ b/prm_reac.C @@ -1 +1 @@ -n14dc11an/prm/prm_reac.n14dc11an_sim005.C \ No newline at end of file +f17do14an/prm/prm_reac.f17do14an_sim001.C \ No newline at end of file diff --git a/simthm.C b/simthm.C index a32a695008adaa03d205e8776c608c027d33f37b..8deba5c2ac9657e3a4ab0e51cfe8e4f8aca3b2e0 100644 --- a/simthm.C +++ b/simthm.C @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ #include #include // for the seed of the dandom number #include +#include //#include "Math/Interpolator.h" //#include "Math/Polynomial.h" @@ -182,12 +184,12 @@ int main() { // (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]; + TVector3 vdet[ndet]; //, voff[ndet]; // Detectors for(Int_t idet=0; idetMakeBox("PSD",Si,psd_half[idet],psd_half[idet],0.5); // in mm + psd[idet] = sc->MakeBox("PSD",Si,psd_xhalf[idet],psd_yhalf[idet],0.5); // in mm + psd[idet]->SetFillColor(kBlue); 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); @@ -199,23 +201,97 @@ int main() { // 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]); + //voff[idet].SetXYZ(xdet_off[idet],0.,zdet_off[idet]); + } + + // Detector strips + //TPolyLine3D *xstrip[100][ndet]; + const Int_t n_xstrip_max = *max_element(begin(n_xstrip), end(n_xstrip)); + const Int_t n_ystrip_max = *max_element(begin(n_ystrip), end(n_ystrip)); + TPolyLine3D *xstrip[n_xstrip_max][ndet]; + TPolyLine3D *ystrip[n_ystrip_max][ndet]; + for(Int_t idet=0; idetSetPoint(0,vdet[idet].X() + + double((-n_xstrip[idet]/2.+i+1)*2./n_xstrip[idet])*psd_xhalf[idet]*cos(d2r*thdet[idet]), + vdet[idet].Y()-psd_yhalf[idet], + vdet[idet].Z() + + double((-n_xstrip[idet]/2.+i+1)*2./n_xstrip[idet])*psd_xhalf[idet]*sin(d2r*thdet[idet])); + xstrip[i][idet]->SetPoint(1,vdet[idet].X() + + double((-n_xstrip[idet]/2.+i+1)*2./n_xstrip[idet])*psd_xhalf[idet]*cos(d2r*thdet[idet]), + vdet[idet].Y()+psd_yhalf[idet], + vdet[idet].Z() + + double((-n_xstrip[idet]/2.+i+1)*2./n_xstrip[idet])*psd_xhalf[idet]*sin(d2r*thdet[idet])); + xstrip[i][idet]->SetLineColor(kGreen); + } + for(Int_t i=0; iSetPoint(0,vdet[idet].X() - psd_yhalf[idet]*cos(d2r*thdet[idet]), + vdet[idet].Y()+double((-n_ystrip[idet]/2.+i+1)*2./n_ystrip[idet])*psd_yhalf[idet], + vdet[idet].Z() - psd_yhalf[idet]*sin(d2r*thdet[idet])); + ystrip[i][idet]->SetPoint(1,vdet[idet].X() + psd_yhalf[idet]*cos(d2r*thdet[idet]), + vdet[idet].Y()+double((-n_ystrip[idet]/2.+i+1)*2./n_ystrip[idet])*psd_yhalf[idet], + vdet[idet].Z() + psd_yhalf[idet]*sin(d2r*thdet[idet])); + ystrip[i][idet]->SetLineColor(kRed); + } + } + + // To show short axes at (0,0,0) + Double_t l_ax = 50; + TPolyLine3D *xaxis = new TPolyLine3D(2); + xaxis->SetPoint(0,-l_ax,0,0); + xaxis->SetPoint(1,l_ax,0,0); + xaxis->SetLineColor(kRed); + TPolyLine3D *yaxis = new TPolyLine3D(2); + yaxis->SetPoint(0,0,-l_ax,0); + yaxis->SetPoint(1,0,l_ax,0); + yaxis->SetLineColor(kGreen); + TPolyLine3D *zaxis = new TPolyLine3D(2); + zaxis->SetPoint(0,0,0,-l_ax); + zaxis->SetPoint(1,0,0,l_ax); + zaxis->SetLineColor(kBlue); + // Target shape at (0,0,0) + Double_t r_tar = 15; + TPolyLine3D *target = new TPolyLine3D(32); + for(Int_t i=0; i<32; i++){ + target->SetPoint(i,r_tar*cos(2.*pi/32.*double(i)), + r_tar*sin(2.*pi/32.*double(i)),0); } gGeoManager->CloseGeometry(); - top->SetLineColor(kRed); + top->SetLineColor(kGray); + top->SetLineWidth(1); 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(); + TView3D *vw1 = (TView3D*)TView::CreateView(1); + vw1->SetRange(-r_sphere,-r_sphere,-r_sphere,r_sphere,r_sphere,r_sphere); + vw1->ShowAxis(); + vw1->SetPerspective(); + vw1->RotateView(120,140); + vw1->ZoomView(nullptr,2.5); + xaxis->Draw(); + yaxis->Draw(); + zaxis->Draw(); + target->Draw(); + for(Int_t idet=0; idetDraw(); + for(Int_t i=0; iDraw(); + } 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)); + TFile *f0 = new TFile(Form("%s/figs/det_arrange_%s_%s.root", + dir,label,simnum),"RECREATE"); + f0->cd(); + c0->Write(); + f0->Close(); ////////////////////////////////////// // Output tree file @@ -473,7 +549,8 @@ int main() { 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 + if(istate[iprt]) + 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 } @@ -519,16 +596,16 @@ int main() { /* 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]] ); + n_xstrip[hit.idet[iprt]], + psd_xhalf[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] ] ); + n_ystrip[ hit.idet[iprt] ], + psd_yhalf[ 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 @@ -583,6 +660,10 @@ int main() { geve++; } // if istate2 + /* cout << hit.idet[0] << " " */ + /* << hit.idet[1] << " " */ + /* << hit.any << " " */ + /* << hit.both << endl; */ //if(rnum<=eff_thk(ps)){ sim->Fill(); @@ -597,7 +678,7 @@ int main() { cout << "Done!" << endl; //return 0; - //theApp.Run(); + theApp.Run(); getchar(); return 0;