void a(){ gStyle->SetTitleOffset(1.3,"y"); Float_t beam_en = 11.2; Int_t cnt; float yield[103]; float eyield[103]; float lumi[103]; float elumi[103]; double ratio[103]; double eratio[103]; double x[103]; double ex[103]; // --------- float yield_h[275]; float eyield_h[275]; float lumi_h[275]; float elumi_h[275]; double ratio_h[275]; double eratio_h[275]; double x_h[275]; double ex_h[275] ; double x_m[275]; double ex_m[275]; memset(ex,0,sizeof(ex)); memset(ex_h,0,sizeof(ex_h)); memset(ex_m,0,sizeof(ex_m)); Double_t tagh_en[275]; Double_t tagm_en[103]; Int_t ind[300]; Float_t emin, emax; // FILE *file_yield = fopen("comp_61914_tagm.txt","r"); // FILE *file_yield = fopen("rrr_tagm.txt","r"); FILE *file_yield = fopen("rrr_tagm_new.txt","r"); FILE *file_lumi = fopen("61914_tagm_ps_acc_cor.txt","r"); // FILE *file_yield = fopen("comp_61915_tagm.txt","r"); // FILE *file_lumi = fopen("61915_tagm_ps_acc_cor.txt","r"); // FILE *file_yield = fopen("comp_61340_tagm.txt","r"); // FILE *file_lumi = fopen("61340_tagm_ps_acc_cor.txt","r"); // FILE *file_yield = fopen("comp_61950_tagm.txt","r"); // FILE *file_lumi = fopen("61950_tagm_ps_acc_cor.txt","r"); // FILE *file_yield = fopen("comp_61868_tagm.txt","r"); // FILE *file_lumi = fopen("61868_tagm_ps_acc_cor.txt","r"); // FILE *file_yield_h = fopen("comp_61914_tagh.txt","r"); // FILE *file_yield_h = fopen("rrr_tagh.txt","r"); // OLD SASCHA // FILE *file_yield_h = fopen("comp_61914_tagh_n.txt","r"); // NEW SASCHA // FILE *file_yield_h = fopen("comp_61914_tagh_NEW.txt","r"); // FILE *file_lumi_h = fopen("tagh_from_ccdb.txt","r"); // OLD SASCHA // FILE *file_lumi_h = fopen("61914_from_ccdb_acc_cor.txt","r"); // NEW SASCHA // FILE *file_lumi_h = fopen("61914_tagh_ps_acc_cor.txt","r"); // FILE *file_yield_h = fopen("comp_61915_tagh.txt","r"); // FILE *file_lumi_h = fopen("61915_tagh_ps_acc_cor.txt","r"); // FILE *file_yield_h = fopen("comp_61340_tagh.txt","r"); // FILE *file_lumi_h = fopen("61340_tagh_ps_acc_cor.txt","r"); FILE *file_yield_h = fopen("comp_61950_tagh.txt","r"); FILE *file_lumi_h = fopen("61950_tagh_ps_acc_cor.txt","r"); // FILE *file_yield_h = fopen("comp_61868_tagh.txt","r"); // FILE *file_lumi_h = fopen("61868_tagh_ps_acc_cor.txt","r"); FILE *file_tagm = fopen("tagm_en.txt","r"); FILE *file_tagh = fopen("tagh_en.txt","r"); int ii = 0; while (!feof(file_tagm)){ fscanf(file_tagm,"%d %f %f \n",&ind[ii], &emin, &emax); tagm_en[ii] = beam_en*(emin + emax)/2.; cout << ind[ii] << " x = " << emin << " " << emax << " TAGM energy = " << tagm_en[ii] << endl; ii++; } fclose(file_tagm); #if 1 ii = 0; while (!feof(file_tagh)){ fscanf(file_tagh,"%d %f %f \n",&ind[ii], &emin, &emax); tagh_en[ii] = beam_en*(emin + emax)/2.; cout << ind[ii] << " x = " << emin << " " << emax << " TAGH energy = " << tagh_en[ii] << endl; ii++; } fclose(file_tagh); #endif ii = 1; while (!feof(file_yield)){ fscanf(file_yield,"%d %f %f \n",&cnt, &yield[ii], &eyield[ii]); x[ii] = ii; ii++; } cout << " I am here 0" << endl; ii = 1; while (!feof(file_lumi)){ fscanf(file_lumi,"%d %f %f \n",&cnt, &lumi[ii], &elumi[ii]); ii++; } cout << " I am here 7" << endl; // --------- #if 1 ii = 1; while (!feof(file_yield_h)){ fscanf(file_yield_h,"%d %f %f \n",&cnt, &yield_h[ii], &eyield_h[ii]); x_h[ii] = ii; ii++; } cout << " I am here 0" << endl; ii = 1; while (!feof(file_lumi_h)){ fscanf(file_lumi_h,"%d %f %f \n",&cnt, &lumi_h[ii], &elumi_h[ii]); ii++; } #endif cout << " I am here " << endl; memset(ratio,0,sizeof(ratio)); memset(eratio,0,sizeof(eratio)); memset(ratio_h,0,sizeof(ratio_h)); memset(eratio_h,0,sizeof(eratio_h)); for(ii = 1; ii < 103; ii++){ if(lumi[ii] != 0){ ratio[ii] = yield[ii]/lumi[ii]; if(yield[ii] > 0) eratio[ii] = ratio[ii]*eyield[ii]/yield[ii]; else eratio[ii] = -10; } cout << " Ratio TAGM = " << ratio[ii] << " " << eratio[ii] << endl; } #if 1 for(ii = 1; ii < 275; ii++){ if(lumi_h[ii] != 0){ ratio_h[ii] = yield_h[ii]/lumi_h[ii]; if(yield_h[ii] > 0) eratio_h[ii] = ratio_h[ii]*eyield_h[ii]/yield_h[ii]; else eratio_h[ii] = -10; } cout << " Ratio TAGH = " << ratio_h[ii] << " " << eratio_h[ii] << " yield = " << yield_h[ii] << " " << eyield_h[ii] << " " << lumi_h[ii] << endl; } #endif TCanvas *c1 = new TCanvas("c1","A Simple Graph with error bars", 200, 10, 700, 500); TGraphErrors *gr = new TGraphErrors(274, tagh_en, ratio_h, ex_h, eratio_h); TGraphErrors *gr1 = new TGraphErrors(103, tagm_en, ratio, ex, eratio); // gr->SetMaximum(0.35); gr->SetMinimum(5e-6); gr->SetMaximum(40e-6); gr->SetMarkerStyle(20); gr->SetMarkerColor(4); gr->Draw("AP"); gr1->SetMarkerStyle(21); gr1->SetMarkerColor(2); gr1->Draw("Psame"); gr->SetTitle("Flux normalized yield"); gr->GetXaxis()->SetTitle("Tagger Energy (GeV) "); gr->GetYaxis()->SetTitle("Yield/Flux (A.U.)"); // Cross sections TCanvas *c2 = new TCanvas("c2","c2", 200, 10, 700, 500); gPad->SetGrid(); Double_t cross_tagm[103]; Double_t err_cross_tagm[103]; Double_t cross_tagh[275]; Double_t err_cross_tagh[275]; for(ii = 1; ii < 275; ii++){ Double_t p0 = -0.412474; Double_t p1 = 0.195072; Double_t p2 = -0.0231419; Double_t p3 = 0.000836805; Double_t eff = p0 + p1*tagh_en[ii] + p2*pow(tagh_en[ii],2) + p3*pow(tagh_en[ii],3); if(eff <= 0) eff = 100000; cross_tagh[ii] = 10.*1000.*ratio_h[ii]/5.64/eff; err_cross_tagh[ii] = 10.*1000.*eratio_h[ii]/5.64/eff; // cross_tagh[ii] = 10.*1000.*ratio_h[ii]/5.49/eff; // err_cross_tagh[ii] = 10.*1000.*eratio_h[ii]/5.49/eff; // Be // cross_tagh[ii] = 10.*1000.*ratio_h[ii]/2.194/eff; // err_cross_tagh[ii] = 10.*1000.*eratio_h[ii]/2.194/eff; cout << tagh_en[ii] << " " << eff << " " << ratio_h[ii] << " " << cross_tagh[ii] << " " << eff << endl; } for(ii = 1; ii < 103; ii++){ Double_t p0 = -0.412474; Double_t p1 = 0.195072; Double_t p2 = -0.0231419; Double_t p3 = 0.000836805; Double_t eff = p0 + p1*tagm_en[ii] + p2*pow(tagm_en[ii],2) + p3*pow(tagm_en[ii],3); if(eff <= 0) eff = 100000; cross_tagm[ii] = 10.*1000.*ratio[ii]/5.64/eff; err_cross_tagm[ii] = 10.*1000.*eratio[ii]/5.64/eff; // cross_tagm[ii] = 10.*1000.*ratio[ii]/5.49/eff; // err_cross_tagm[ii] = 10.*1000.*eratio[ii]/5.49/eff; // cross_tagm[ii] = 10.*1000.*ratio[ii]/2.194/eff; // err_cross_tagm[ii] = 10.*1000.*eratio[ii]/2.194/eff; // cout << tagh_en[ii] << " " << eff << " " << ratio_h[ii] << // " " << cross_tagh[ii] << " " << eff << endl; } TGraphErrors *gr11 = new TGraphErrors(274, tagm_en, cross_tagm, ex_m, err_cross_tagm); TGraphErrors *gr10 = new TGraphErrors(274, tagh_en, cross_tagh, ex_h, err_cross_tagh); Double_t x_nist[5] = {5, 6, 8, 10, 15}; // Helium Double_t cross_nist[5] = {5.542E-01, 4.699E-01, 3.621E-01, 2.956E-01, 2.043E-01}; // Berilium // Double_t cross_nist[5] = {11.08E-01, 9.399E-01, 7.241E-01, // 5.912E-01, 4.086E-01}; Double_t ex_nist[5]; memset(ex_nist,0,sizeof(ex_nist)); TGraphErrors *gr20 = new TGraphErrors(5, x_nist, cross_nist, ex_nist, ex_nist); auto axis = gr1->GetXaxis(); axis->SetLimits(5.5,11.); // gr10->SetMaximum(1.2); gr10->SetMaximum(0.6); gr10->SetMinimum(0.1); gr10->SetMarkerStyle(20); gr10->SetMarkerColor(4); gr10->Draw("AP"); gr10->GetXaxis()->SetRangeUser(5.5,11); gr10->Draw("AP"); // gr10->SetTitle("He target, Run 61914, 50 nA"); gr10->SetTitle("He target, Run 61950, 100 nA"); // gr10->SetTitle("Be target, Run 61340, 200 nA"); gr10->GetXaxis()->SetTitle("Beam Energy (GeV)"); gr10->GetYaxis()->SetTitle("Cross section (mb / atom)"); gr11->SetMarkerStyle(21); gr11->SetMarkerColor(8); gr11->Draw("Psame"); gr20->SetMarkerColor(2); gr20->SetLineColor(2); gr20->SetLineWidth(3); gr20->Draw("Csame"); // gr1->SetMarkerStyle(21); // gr1->SetMarkerColor(2); // gr1->Draw("Psame"); // gr->SetTitle("Flux normalized yield"); // gr->GetXaxis()->SetTitle("Tagger Energy (GeV) "); // gr->GetYaxis()->SetTitle("Yield/Flux (A.U.)"); leg = new TLegend(0.45,0.7,0.72,0.9); leg->SetFillColor(0); leg->SetTextSize(0.04); leg->AddEntry(gr20,"NIST ","L"); leg->Draw(); }