File indexing completed on 2024-04-06 12:30:18
0001 #include <iostream>
0002 #include <vector>
0003 #include <string>
0004 #include <map>
0005
0006 int mstyle[15] = { 29, 29, 20, 21, 20, 21, 22, 21, 22, 20, 21, 22, 20, 20, 21};
0007 int mcolor[15] = { 1, 1, 3, 3, 4, 4, 4, 7, 7, 7, 2, 2, 2, 28, 28};
0008 float msiz[15] = {1.4,1.4,1.1,1.1,1.1,1.1,1.4,1.1,1.4,1.1,1.1,1.4,1.1,1.1,1.1};
0009 int lcolor[15] = { 1, 1, 3, 3, 4, 4, 4, 7, 7, 7, 2, 2, 2, 28, 28};
0010 int lstyle[15] = { 1, 1, 1, 2, 1, 2, 3, 2, 3, 1, 3, 2, 1, 1, 2};
0011 int lwidth[15] = { 1, 1, 1, 2, 1, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2};
0012
0013 void plotMomentum(char target[6], char list[20], char ene[6], char part[4],
0014 char dir[12]="histo", char g4ver[20]="G4.9.1.p01") {
0015
0016 setStyle();
0017 gStyle->SetOptLogy(1);
0018
0019 std::vector<std::string> types = typesOld();
0020 int energy = atoi(ene);
0021
0022 char ofile[100];
0023 sprintf (ofile, "%s/histo_%s%s_%s_%sGeV.root", dir, target, list, part, ene);
0024 std::cout << "Input file " << ofile << "\n";
0025 TFile *fout = TFile::Open(ofile);
0026 fout->cd();
0027
0028 char name[160], title[160], ctype[20], ytitle[20], cname[160];
0029 TH1F *hiParticle[5][20];
0030 for (unsigned int ii=0; ii<=(types.size()); ii++) {
0031 if (ii == 0) sprintf (ctype, "All Particles");
0032 else sprintf (ctype, "%s", types[ii-1].c_str());
0033 for (unsigned int jj=0; jj<5; jj++) {
0034 sprintf (name, "Particle%i_KE%s%s%sGeV(%s)",jj,target, list, ene, ctype);
0035 hiParticle[jj][ii] = (TH1F*)fout->FindObjectAny(name);
0036 hiParticle[jj][ii]->SetFillColor(30);
0037 }
0038 }
0039
0040 TCanvas *c[5];
0041 for (unsigned int jj=0; jj<2; jj++) {
0042 hiParticle[jj][11]->Rebin(25);
0043 hiParticle[jj][12]->Rebin(25);
0044 hiParticle[jj][13]->Rebin(25);
0045
0046 sprintf(cname, "c_%s%s%sGeV_nucleons_particle%i", target, list, ene, jj);
0047 c[jj] = new TCanvas(cname, cname, 800, 800);
0048 c[jj]->Divide(2,2);
0049 c[jj]->cd(1); hiParticle[jj][11]->Draw();
0050 c[jj]->cd(2); hiParticle[jj][12]->Draw();
0051 c[jj]->cd(3); hiParticle[jj][13]->Draw();
0052 c[jj]->cd(4); hiParticle[jj][14]->Draw();
0053 }
0054 }
0055
0056 void plotParticles(char target[6], char list[20], char ene[6], char part[4],
0057 char dir[12]="histo", char g4ver[20]="G4.9.1.p01") {
0058
0059 gStyle->SetOptLogy(1);
0060 gStyle->SetTitleX(.1);
0061 gStyle->SetTitleY(.9);
0062
0063
0064 std::vector<std::string> types = types();
0065 char ofile[100];
0066 sprintf (ofile, "%s/histo_%s%s_%s_%sGeV.root", dir, target, list, part, ene);
0067 std::cout << "Input file " << ofile << "\n";
0068 TFile *fout = TFile::Open(ofile);
0069 fout->cd();
0070
0071 char ctype[20], title[160], cname[160];
0072 sprintf(ctype,"Ions");
0073 TH1F *hProton[2], *hNeutron[2], *hHeavy[2], *hIon[2];
0074 TH1F *hProton_1[2], *hNeutron_1[2], *hHeavy_1[2], *hIon_1[2];
0075 for (int i=0; i<2; i++) {
0076 sprintf(title, "proton%i_%s%s%sGeV(%s)", i, target, list, ene, ctype);
0077 hProton[i] = (TH1F*)fout->FindObjectAny(title);
0078 sprintf(title, "proton%i_%s%s%sGeV", i, target, list, ene);
0079 hProton[i]->SetName(title);
0080 hProton[i]->SetTitle(title);
0081
0082 hProton_1[i] = (TH1F*) hProton[i]->Clone();
0083 sprintf(title, "new_%s", title);
0084 hProton_1[i]->SetName(title);
0085
0086 sprintf(title, "neutron%i_%s%s%sGeV(%s)", i, target, list, ene, ctype);
0087 hNeutron[i] = (TH1F*)fout->FindObjectAny(title);
0088 sprintf(title, "neutron%i_%s%s%sGeV", i, target, list, ene);
0089 hNeutron[i]->SetName(title);
0090 hNeutron[i]->SetTitle(title);
0091
0092 hNeutron_1[i] = (TH1F*) hNeutron[i]->Clone();
0093 sprintf(title, "new_%s", title);
0094 hNeutron_1[i]->SetName(title);
0095
0096 sprintf(title, "heavy%i_%s%s%sGeV(%s)", i, target, list, ene, ctype);
0097 hHeavy[i] = (TH1F*)fout->FindObjectAny(title);
0098 sprintf(title, "heavy%i_%s%s%sGeV", i, target, list, ene);
0099 hHeavy[i]->SetName(title);
0100 hHeavy[i]->SetTitle(title);
0101
0102 hHeavy_1[i] = (TH1F*) hHeavy[i]->Clone();
0103 sprintf(title, "new_%s", title);
0104 hHeavy_1[i]->SetName(title);
0105
0106 sprintf(title, "ion%i_%s%s%sGeV(%s)", i, target, list, ene, ctype);
0107 hIon[i] = (TH1F*)fout->FindObjectAny(title);
0108 sprintf(title, "ion%i_%s%s%sGeV", i, target, list, ene);
0109 hIon[i]->SetName(title);
0110 hIon[i]->SetTitle(title);
0111
0112 hIon_1[i] = (TH1F*) hIon[i]->Clone();
0113 sprintf(title, "new_%s", title);
0114 hIon_1[i]->SetName(title);
0115 }
0116
0117 int energy = atoi(ene);
0118 std::cout << "Energy " << energy << "\n";
0119 if (energy>=10) {
0120 hProton[0]->Rebin(5); hProton[1]->Rebin(5);
0121 hProton_1[0]->GetXaxis()->SetRangeUser(0.0, 5.0);
0122 hProton_1[1]->GetXaxis()->SetRangeUser(0.0, 2.5);
0123
0124 hNeutron[0]->Rebin(5); hNeutron[1]->Rebin(5);
0125 hNeutron_1[0]->GetXaxis()->SetRangeUser(0.0, 5.0);
0126 hNeutron_1[1]->GetXaxis()->SetRangeUser(0.0, 2.5);
0127
0128 hHeavy[0]->Rebin(5); hHeavy[1]->Rebin(5);
0129 hHeavy_1[0]->GetXaxis()->SetRangeUser(0.0, 5.0);
0130 hHeavy_1[1]->GetXaxis()->SetRangeUser(0.0, 5.0);
0131
0132 hIon[0]->GetXaxis()->SetRangeUser(0.0, 0.03);
0133 hIon[1]->GetXaxis()->SetRangeUser(0.0, 0.03);
0134 } else {
0135
0136 hProton_1[0]->GetXaxis()->SetRangeUser(0.0, 1.0);
0137 hProton_1[1]->GetXaxis()->SetRangeUser(0.0, 0.5);
0138
0139 hNeutron_1[0]->GetXaxis()->SetRangeUser(0.0, 1.0);
0140 hNeutron_1[1]->GetXaxis()->SetRangeUser(0.0, 0.5);
0141
0142 hHeavy_1[0]->GetXaxis()->SetRangeUser(0.0, 5.0);
0143 hHeavy_1[1]->GetXaxis()->SetRangeUser(0.0, 5.0);
0144
0145 hIon[0]->Rebin(5); hIon[1]->Rebin(5);
0146 hIon[0]->GetXaxis()->SetRangeUser(0.0, 1.5);
0147 hIon[1]->GetXaxis()->SetRangeUser(0.0, 1.5);
0148 }
0149
0150 for (int i=0; i<2; i++) {
0151 hProton[i]->SetFillColor(30);
0152 hNeutron[i]->SetFillColor(30);
0153 hHeavy[i]->SetFillColor(30);
0154 hIon[i]->SetFillColor(30);
0155
0156 hProton_1[i]->SetFillColor(30);
0157 hNeutron_1[i]->SetFillColor(30);
0158 hHeavy_1[i]->SetFillColor(30);
0159 hIon_1[i]->SetFillColor(30);
0160
0161 hProton[i]->GetXaxis()->SetRangeUser(0.0, energy);
0162 hNeutron[i]->GetXaxis()->SetRangeUser(0.0, energy);
0163 hHeavy[i]->GetXaxis()->SetRangeUser(0.0, energy);
0164 }
0165
0166 sprintf(cname, "c_%s%s%sGeV_protons", target, list, ene);
0167 TCanvas *cc4 = new TCanvas(cname, cname, 800, 800);
0168 cc4->Divide(2,2);
0169 cc4->cd(1); hProton[0]->Draw();
0170 cc4->cd(2); hProton_1[0]->Draw();
0171 cc4->cd(3); hProton[1]->Draw();
0172 cc4->cd(4); hProton_1[1]->Draw();
0173
0174 sprintf(cname, "c_%s%s%sGeV_neutrons", target, list, ene);
0175 TCanvas *cc5 = new TCanvas(cname, cname, 800, 800);
0176 cc5->Divide(2,2);
0177 cc5->cd(1); hNeutron[0]->Draw();
0178 cc5->cd(2); hNeutron_1[0]->Draw();
0179 cc5->cd(3); hNeutron[1]->Draw();
0180 cc5->cd(4); hNeutron_1[1]->Draw();
0181
0182 sprintf(cname, "c_%s%s%sGeV_Heavy", target, list, ene);
0183 TCanvas *cc6 = new TCanvas(cname, cname, 800, 800);
0184 cc6->Divide(2,2);
0185 cc6->cd(1); hHeavy[0]->Draw();
0186 cc6->cd(2); hHeavy_1[0]->Draw();
0187 cc6->cd(3); hHeavy[1]->Draw();
0188 cc6->cd(4); hHeavy_1[1]->Draw();
0189
0190 sprintf(cname, "c_%s%s%sGeV_Ion", target, list, ene);
0191 TCanvas *cc7 = new TCanvas(cname, cname, 800, 500);
0192 cc7->Divide(2,1);
0193 cc7->cd(1); hIon[0]->Draw();
0194 cc7->cd(2); hIon[1]->Draw();
0195 }
0196
0197
0198 void plotMultiplicity(char target[6], char list[20], char part[4], int ymax=25,
0199 char dir[12]="histo", char g4ver[20]="G4.9.1.p01",
0200 bool flag=true) {
0201
0202 setStyle();
0203 gStyle->SetOptTitle(0);
0204
0205 char name[1024], sym[10];
0206 if (part=="pim") sprintf(sym, "#pi^{-}");
0207 else if (part=="pip") sprintf(sym, "#pi^{+}");
0208 else sprintf(sym, "p");
0209
0210 std::map<string, double> means_300=getMean(target,list,part,"300.0","Multi",dir);
0211 std::map<string, double> means_200=getMean(target,list,part,"200.0","Multi",dir);
0212 std::map<string, double> means_150=getMean(target,list,part,"150.0","Multi",dir);
0213 std::map<string, double> means_100=getMean(target,list,part,"100.0","Multi",dir);
0214 std::map<string, double> means_50 =getMean(target,list,part,"50.0", "Multi",dir);
0215 std::map<string, double> means_30 =getMean(target,list,part,"30.0", "Multi",dir);
0216 std::map<string, double> means_20 =getMean(target,list,part,"20.0", "Multi",dir);
0217 std::map<string, double> means_15 =getMean(target,list,part,"15.0", "Multi",dir);
0218 std::map<string, double> means_9 =getMean(target,list,part,"9.0", "Multi",dir);
0219 std::map<string, double> means_7 =getMean(target,list,part,"7.0", "Multi",dir);
0220 std::map<string, double> means_5 =getMean(target,list,part,"5.0", "Multi",dir);
0221 std::map<string, double> means_3 =getMean(target,list,part,"3.0", "Multi",dir);
0222 std::map<string, double> means_2 =getMean(target,list,part,"2.0", "Multi",dir);
0223 std::map<string, double> means_1 =getMean(target,list,part,"1.0", "Multi",dir);
0224 if (flag) {
0225 std::map<string, double> means_10 =getMean(target,list,part,"10.0", "Multi",dir);
0226 std::map<string, double> means_8 =getMean(target,list,part,"8.0", "Multi",dir);
0227 std::map<string, double> means_6 =getMean(target,list,part,"6.0", "Multi",dir);
0228 std::map<string, double> means_4 =getMean(target,list,part,"4.0", "Multi",dir);
0229 }
0230
0231 char ctype[20];
0232 std::vector<std::string> types = types();
0233 std::vector<std::string> typeOld = typesOld();
0234
0235
0236 TGraph *gr[20];
0237 TLegend *leg = new TLegend(0.45, 0.53, 0.90, 0.90);
0238 char hdr[160];
0239 sprintf(hdr, "%s+%s (%s-%s)", sym, target, g4ver, list);
0240 leg->SetHeader(hdr); leg->SetFillColor(10); leg->SetMargin(0.45);
0241 leg->SetTextSize(.027);
0242 sprintf(name, "c_%s_%sMultiplicity_%s", part,target,list);
0243 TCanvas *cc = new TCanvas(name, name, 700, 700);
0244
0245 for (unsigned int ii=0; ii<=(types.size()); ii++) {
0246 if (ii == 0) sprintf (ctype, "All Particles");
0247 else sprintf (ctype, "%s", typeOld[ii-1].c_str());
0248
0249
0250
0251 string a(ctype);
0252 double vx[18], vy[18];
0253 int np=0;
0254 vx[np] = 300.0; vy[np] = means_300[a]; np++;
0255 vx[np] = 200.0; vy[np] = means_200[a]; np++;
0256 vx[np] = 150.0; vy[np] = means_150[a]; np++;
0257 vx[np] = 100.0; vy[np] = means_100[a]; np++;
0258 vx[np] = 50.0; vy[np] = means_50[a]; np++;
0259 vx[np] = 30.0; vy[np] = means_30[a]; np++;
0260 vx[np] = 20.0; vy[np] = means_20[a]; np++;
0261 vx[np] = 15.0; vy[np] = means_15[a]; np++;
0262 if (flag) { vx[np] = 10.0; vy[np] = means_10[a]; np++;}
0263 vx[np] = 9.0; vy[np] = means_9[a]; np++;
0264 if (flag) { vx[np] = 8.0; vy[np] = means_8[a]; np++;}
0265 vx[np] = 7.0; vy[np] = means_7[a]; np++;
0266 if (flag) { vx[np] = 6.0; vy[np] = means_6[a]; np++;}
0267 vx[np] = 5.0; vy[np] = means_5[a]; np++;
0268 if (flag && part != "pro") { vx[np] = 4.0; vy[np] = means_4[a]; np++;}
0269 vx[np] = 3.0; vy[np] = means_3[a]; np++;
0270 vx[np] = 2.0; vy[np] = means_2[a]; np++;
0271 vx[np] = 1.0; vy[np] = means_1[a]; np++;
0272
0273 if (ii > 20 ) {
0274 std::cout << ctype;
0275 for (int ix=0; ix<np; ix++) std::cout << " " << vx[ix] << " " << vy[ix];
0276 std::cout << "\n";
0277 }
0278
0279 gPad->SetLogx(1);
0280 gPad->SetGridx(1); gPad->SetGridy(1);
0281 gr[ii] = new TGraph(np, vx,vy);
0282 sprintf(name, "Multiplicity of secondaries %s-%s (%s %s)", sym, target, g4ver, list);
0283 gr[ii]->SetTitle(name);
0284 gr[ii]->GetXaxis()->SetTitle("Beam Momentum (GeV)");
0285 gr[ii]->GetYaxis()->SetTitle("Average Multiplicity");
0286
0287 gr[ii]->SetMarkerStyle(mstyle[ii]);
0288 gr[ii]->SetMarkerSize(msiz[ii]);
0289 gr[ii]->SetMarkerColor(mcolor[ii]);
0290 gr[ii]->SetLineColor(lcolor[ii]);
0291 gr[ii]->SetLineStyle(lstyle[ii]);
0292 gr[ii]->SetLineWidth(lwidth[ii]);
0293
0294 gr[ii]->GetYaxis()->SetRangeUser(-0.2, ymax);
0295 if (ii>1) {
0296 sprintf (ctype, "%s", types[ii-1].c_str());
0297 leg->AddEntry(gr[ii], ctype, "lP");
0298 }
0299 if (ii==2) gr[ii]->Draw("APl");
0300 else if (ii>2) gr[ii]->Draw("Pl");
0301 }
0302 leg->Draw("same");
0303 }
0304
0305 void plotMultiplicity(char target[6], char list[20], char ene[6], char part[4],
0306 char dir[12]="histo", char g4ver[20]="G4.9.1.p01") {
0307
0308 setStyle();
0309 gStyle->SetOptTitle(0);
0310
0311 char name[1024], sym[10];
0312 if (part=="pim") sprintf(sym, "#pi^{-}");
0313 else if (part=="pip") sprintf(sym, "#pi^{+}");
0314 else sprintf(sym, "p");
0315
0316 std::vector<std::string> typeOld = typesOld();
0317 int energy = atoi(ene);
0318
0319 char ofile[100];
0320 sprintf (ofile, "%s/histo_%s%s_%s_%sGeV.root", dir, target, list, part, ene);
0321 std::cout << "Input file " << ofile << "\n";
0322 TFile *fout = TFile::Open(ofile);
0323 fout->cd();
0324
0325 char name[160], title[160], ctype[20], ytitle[20], cname[160];
0326 TH1I *hiMulti[20];
0327 for (unsigned int ii=0; ii<=(typeOld.size()); ii++) {
0328 if (ii == 0) sprintf (ctype, "All Particles");
0329 else sprintf (ctype, "%s", typeOld[ii-1].c_str());
0330 sprintf (name, "Multi%s%s%sGeV(%s)", target, list, ene, ctype);
0331 hiMulti[ii] = (TH1I*)fout->FindObjectAny(name);
0332
0333 }
0334
0335 TCanvas *c[20];
0336 std::vector<std::string> types = types();
0337 for (unsigned int ii=0; ii<types.size(); ii++) {
0338 if (ii == 0) sprintf (ctype, "All Particles");
0339 else sprintf (ctype, "%s", types[ii-1].c_str());
0340 sprintf (cname, "Multiplicity (%s)", ctype);
0341 hiMulti[ii]->GetXaxis()->SetTitle(cname);
0342 hiMulti[ii]->SetMarkerStyle(mstyle[ii]);
0343 hiMulti[ii]->SetMarkerSize(msiz[ii]);
0344 hiMulti[ii]->SetMarkerColor(mcolor[ii]);
0345 hiMulti[ii]->SetLineColor(lcolor[ii]);
0346 hiMulti[ii]->SetLineStyle(lstyle[ii]);
0347 hiMulti[ii]->SetLineWidth(lwidth[ii]);
0348
0349 sprintf(cname, "c_%s%s_%s_%sGeV_Multiplicity(%s)", target, list, part,
0350 ene, ctype);
0351 c[ii] = new TCanvas(cname, cname, 800, 500);
0352 hiMulti[ii]->Draw();
0353
0354 TLegend *leg = new TLegend(0.35, 0.80, 0.8, 0.87);
0355 char hdr[160];
0356 sprintf(hdr, "%s+%s at %s GeV (%s-%s)", sym, target, ene, g4ver, list);
0357 leg->SetHeader(hdr); leg->SetFillColor(10); leg->SetMargin(0.45);
0358 leg->SetTextSize(.036); leg->Draw("same");
0359 }
0360 }
0361
0362 void plotTotalKE(char target[6], char list[20], char part[4],
0363 char dir[12]="histo", char g4ver[20]="G4.9.1.p01",
0364 bool flag=true) {
0365
0366 setStyle();
0367 gStyle->SetOptTitle(0);
0368
0369 char name[1024];
0370 char sym[10];
0371 if (part=="pim") sprintf(sym, "#pi^{-}");
0372 else if (part=="pip") sprintf(sym, "#pi^{+}");
0373 else sprintf(sym, "p");
0374
0375 std::map<string, double> means_300=getMean(target,list,part,"300.0","TotalKE",dir);
0376 std::map<string, double> means_200=getMean(target,list,part,"200.0","TotalKE",dir);
0377 std::map<string, double> means_150=getMean(target,list,part,"150.0","TotalKE",dir);
0378 std::map<string, double> means_100=getMean(target,list,part,"100.0","TotalKE",dir);
0379 std::map<string, double> means_50 =getMean(target,list,part,"50.0", "TotalKE",dir);
0380 std::map<string, double> means_30 =getMean(target,list,part,"30.0", "TotalKE",dir);
0381 std::map<string, double> means_20 =getMean(target,list,part,"20.0", "TotalKE",dir);
0382 std::map<string, double> means_15 =getMean(target,list,part,"15.0", "TotalKE",dir);
0383 std::map<string, double> means_9 =getMean(target,list,part,"9.0", "TotalKE",dir);
0384 std::map<string, double> means_7 =getMean(target,list,part,"7.0", "TotalKE",dir);
0385 std::map<string, double> means_5 =getMean(target,list,part,"5.0", "TotalKE",dir);
0386 std::map<string, double> means_3 =getMean(target,list,part,"3.0", "TotalKE",dir);
0387 std::map<string, double> means_2 =getMean(target,list,part,"2.0", "TotalKE",dir);
0388 std::map<string, double> means_1 =getMean(target,list,part,"1.0", "TotalKE",dir);
0389 if (flag) {
0390 std::map<string, double> means_10 =getMean(target,list,part,"10.0", "TotalKE",dir);
0391 std::map<string, double> means_8 =getMean(target,list,part,"8.0", "TotalKE",dir);
0392 std::map<string, double> means_6 =getMean(target,list,part,"6.0", "TotalKE",dir);
0393 std::map<string, double> means_4 =getMean(target,list,part,"4.0", "TotalKE",dir);
0394 }
0395
0396 char ctype[20];
0397 std::vector<std::string> types = types();
0398 std::vector<std::string> typeOld = typesOld();
0399
0400
0401 TGraph *gr[20];
0402 TLegend *leg = new TLegend(0.55, 0.45, 0.9, 0.80);
0403 char hdr[160];
0404 sprintf(hdr, "%s+%s (%s-%s)", sym, target, g4ver, list);
0405 leg->SetHeader(hdr);
0406 leg->SetFillColor(10);
0407 leg->SetMargin(0.45);
0408 leg->SetTextSize(.02);
0409 sprintf(name, "c_%s_%s_totalKE_%s", part,target,list);
0410 TCanvas *cc = new TCanvas(name, name, 700, 700);
0411
0412 for (unsigned int ii=0; ii<=(types.size()); ii++) {
0413 if (ii == 0) sprintf (ctype, "All Particles");
0414 else sprintf (ctype, "%s", typeOld[ii-1].c_str());
0415
0416 string a(ctype);
0417
0418 double vx[18], vy[18];
0419 int np=0;
0420 vx[np] = 300.0; vy[np] = means_300[a]; np++;
0421 vx[np] = 200.0; vy[np] = means_200[a]; np++;
0422 vx[np] = 150.0; vy[np] = means_150[a]; np++;
0423 vx[np] = 100.0; vy[np] = means_100[a]; np++;
0424 vx[np] = 50.0; vy[np] = means_50[a]; np++;
0425 vx[np] = 30.0; vy[np] = means_30[a]; np++;
0426 vx[np] = 20.0; vy[np] = means_20[a]; np++;
0427 vx[np] = 15.0; vy[np] = means_15[a]; np++;
0428 if (flag) { vx[np] = 10.0; vy[np] = means_10[a]; np++;}
0429 vx[np] = 9.0; vy[np] = means_9[a]; np++;
0430 if (flag) { vx[np] = 8.0; vy[np] = means_8[a]; np++;}
0431 vx[np] = 7.0; vy[np] = means_7[a]; np++;
0432 if (flag) { vx[np] = 6.0; vy[np] = means_6[a]; np++;}
0433 vx[np] = 5.0; vy[np] = means_5[a]; np++;
0434 if (flag && part != "pro") { vx[np] = 4.0; vy[np] = means_4[a]; np++;}
0435 vx[np] = 3.0; vy[np] = means_3[a]; np++;
0436 vx[np] = 2.0; vy[np] = means_2[a]; np++;
0437 vx[np] = 1.0; vy[np] = means_1[a]; np++;
0438
0439 for (int i=0; i<np; i++) vy[i] = vy[i]/vx[i];
0440
0441 gPad->SetLogx(1);
0442 gPad->SetGridx(1);
0443 gPad->SetGridy(1);
0444 gr[ii] = new TGraph(np, vx,vy);
0445 sprintf(name, "KE carried by secondaries in %s-%s (%s)", sym, target, list);
0446 gr[ii]->SetTitle(name);
0447 gr[ii]->GetXaxis()->SetTitle("Beam Momentum (GeV)");
0448 gr[ii]->GetYaxis()->SetTitle("Mean Total KE/Beam Momentum");
0449
0450 gr[ii]->SetMarkerStyle(mstyle[ii]);
0451 gr[ii]->SetMarkerSize(msiz[ii]);
0452 gr[ii]->SetMarkerColor(mcolor[ii]);
0453 gr[ii]->SetLineColor(lcolor[ii]);
0454 gr[ii]->SetLineStyle(lstyle[ii]);
0455 gr[ii]->SetLineWidth(lwidth[ii]);
0456
0457 gr[ii]->GetYaxis()->SetRangeUser(-0.02, 1.0);
0458 if (ii!= 0) sprintf (ctype, "%s", types[ii-1].c_str());
0459 if (ii!= 1) leg->AddEntry(gr[ii], ctype, "lP");
0460 if (ii==0) gr[ii]->Draw("APl");
0461 else if (ii>1) gr[ii]->Draw("Pl");
0462 }
0463 leg->Draw("same");
0464 }
0465
0466 void plotKE(char target[6], char list[20], char ene[6], char part[4],
0467 int typ=0, char dir[12]="histo", char g4ver[20]="G4.9.1.p01") {
0468
0469 setStyle();
0470 gStyle->SetOptTitle(0);
0471 gStyle->SetOptLogy(1);
0472
0473 char name[1024];
0474 char sym[10];
0475 if (part=="pim") sprintf(sym, "#pi^{-}");
0476 else if (part=="pip") sprintf(sym, "#pi^{+}");
0477 else sprintf(sym, "p");
0478
0479 std::vector<std::string> typeOld = typesOld();
0480 int energy = atoi(ene);
0481 int bins=energy/4;
0482 float ener = energy;
0483 std::cout << "Energy " << ener << "\n";
0484
0485 char ofile[100];
0486 sprintf (ofile, "%s/histo_%s%s_%s_%sGeV.root", dir, target, list, part, ene);
0487 std::cout << "Input file " << ofile << "\n";
0488 TFile *fout = TFile::Open(ofile);
0489 fout->cd();
0490
0491 char name[160], title[160], ctype[20], ytitle[20], cname[160], pre[10];
0492 TH1F *hiKE[20];
0493 if (typ == 0) sprintf (pre, "KE2");
0494 else sprintf (pre, "TotalKE");
0495 for (unsigned int ii=0; ii<=(typeOld.size()); ii++) {
0496 if (ii == 0) sprintf (ctype, "All Particles");
0497 else sprintf (ctype, "%s", typeOld[ii-1].c_str());
0498 sprintf (name, "%s%s%s%sGeV(%s)", pre, target, list, ene, ctype);
0499 hiKE[ii] = (TH1F*)fout->FindObjectAny(name);
0500
0501 }
0502
0503 TCanvas *c[25];
0504 std::vector<std::string> types = types();
0505 for (unsigned int ii=0; ii<types.size(); ii++) {
0506 if (ii == 0) sprintf (ctype, "All Particles");
0507 else sprintf (ctype, "%s", types[ii-1].c_str());
0508 if (typ == 0) sprintf (cname, "Kinetic Energy of %s (GeV)", ctype);
0509 else sprintf (cname, "Total Kinetic Energy of %s (GeV)", ctype);
0510 hiKE[ii]->GetXaxis()->SetTitle(cname);
0511 hiKE[ii]->SetMarkerStyle(mstyle[ii]);
0512 hiKE[ii]->SetMarkerSize(msiz[ii]);
0513 hiKE[ii]->SetMarkerColor(mcolor[ii]);
0514 hiKE[ii]->SetLineColor(lcolor[ii]);
0515 hiKE[ii]->SetLineStyle(lstyle[ii]);
0516 hiKE[ii]->SetLineWidth(lwidth[ii]);
0517 if (bins > 0) hiKE[ii]->Rebin(bins);
0518 hiKE[ii]->GetXaxis()->SetRangeUser(0.0, ener);
0519
0520 sprintf(cname, "c_%s%s_%s_%sGeV_%s(%s)", target,list,part,ene,pre,ctype);
0521 c[ii] = new TCanvas(cname, cname, 800, 500);
0522 hiKE[ii]->Draw();
0523
0524 TLegend *leg = new TLegend(0.35, 0.80, 0.8, 0.87);
0525 char hdr[160];
0526 sprintf(hdr, "%s+%s at %s GeV (%s-%s)", sym, target, ene, g4ver, list);
0527 leg->SetHeader(hdr); leg->SetFillColor(10); leg->SetMargin(0.45);
0528 leg->SetTextSize(.036); leg->Draw("same");
0529 }
0530
0531 TLegend *leg1 = new TLegend(0.50, 0.75, 0.90, 0.90);
0532 if (typ == 0) sprintf (cname, "Kinetic Energy (GeV)");
0533 else sprintf (cname, "Total Kinetic Energy (GeV)");
0534 hiKE[6]->GetXaxis()->SetTitle(cname);
0535 char hdr[160];
0536 sprintf(hdr, "%s+%s at %s GeV (%s-%s)", sym, target, ene, g4ver, list);
0537 leg1->SetHeader(hdr); leg1->SetFillColor(10); leg1->SetMargin(0.45);
0538 sprintf(cname, "c_%s%s_%s_%sGeV_%s(Pion)", target,list,part,ene,pre);
0539 leg1->SetTextSize(.030);
0540 c[19] = new TCanvas(cname, cname, 800, 500);
0541 hiKE[6]->Draw(); sprintf (ctype, "%s", types[5].c_str()); leg1->AddEntry(hiKE[6], ctype, "l");
0542 hiKE[5]->Draw("same"); sprintf (ctype, "%s", types[4].c_str()); leg1->AddEntry(hiKE[5], ctype, "l");
0543 hiKE[4]->Draw("same"); sprintf (ctype, "%s", types[3].c_str()); leg1->AddEntry(hiKE[4], ctype, "l"); leg->Draw("same");
0544
0545 TLegend *leg2 = new TLegend(0.50, 0.75, 0.90, 0.90);
0546 if (typ == 0) sprintf (cname, "Kinetic Energy (GeV)");
0547 else sprintf (cname, "Total Kinetic Energy (GeV)");
0548 hiKE[7]->GetXaxis()->SetTitle(cname);
0549 sprintf(hdr, "%s+%s at %s GeV (%s-%s)", sym, target, ene, g4ver, list);
0550 leg2->SetHeader(hdr); leg2->SetFillColor(10); leg2->SetMargin(0.45);
0551 sprintf(cname, "c_%s%s_%s_%sGeV_%s(Kaon)", target,list,part,ene,pre);
0552 leg2->SetTextSize(.030);
0553 c[20] = new TCanvas(cname, cname, 800, 500);
0554 hiKE[7]->Draw(); sprintf (ctype, "%s", types[6].c_str()); leg2->AddEntry(hiKE[7], ctype, "l");
0555 hiKE[8]->Draw("same"); sprintf (ctype, "%s", types[7].c_str()); leg2->AddEntry(hiKE[8], ctype, "l");
0556 hiKE[9]->Draw("same"); sprintf (ctype, "%s", types[8].c_str()); leg2->AddEntry(hiKE[9], ctype, "l"); leg2->Draw("same");
0557
0558 TLegend *leg3 = new TLegend(0.50, 0.75, 0.90, 0.90);
0559 if (typ == 0) sprintf (cname, "Kinetic Energy (GeV)");
0560 else sprintf (cname, "Total Kinetic Energy (GeV)");
0561 hiKE[12]->GetXaxis()->SetTitle(cname);
0562 sprintf(hdr, "%s+%s at %s GeV (%s-%s)", sym, target, ene, g4ver, list);
0563 leg3->SetHeader(hdr); leg3->SetFillColor(10); leg3->SetMargin(0.45);
0564 sprintf(cname, "c_%s%s_%s_%sGeV_%s(Nucleon)", target,list,part,ene,pre);
0565 leg3->SetTextSize(.030);
0566 c[21] = new TCanvas(cname, cname, 800, 500);
0567 hiKE[12]->Draw(); sprintf (ctype, "%s", types[11].c_str()); leg3->AddEntry(hiKE[12], ctype, "l");
0568 hiKE[11]->Draw("same"); sprintf (ctype, "%s", types[10].c_str()); leg3->AddEntry(hiKE[11], ctype, "l"); leg3->Draw("same");
0569 }
0570
0571 void printMeans(std::map<string, double> means) {
0572
0573 std::map<string, double>::iterator iter;
0574 for( iter = means.begin(); iter != means.end(); iter++ ) {
0575 std::cout << (*iter).first << " " << (*iter).second << "\n";
0576 }
0577 }
0578
0579 std::map<string, double> getMean(char target[6], char list[20], char part[5],
0580 char ene[6], char ctyp0[10]="Multi",
0581 char dir[12]="histo") {
0582
0583 std::vector<std::string> types = typesOld();
0584 std::map<string, double> means;
0585
0586 char ofile[100];
0587 sprintf (ofile, "%s/histo_%s%s_%s_%sGeV.root", dir, target, list, part, ene);
0588 std::cout << "Input File: " << ofile << "\n";
0589 TFile *fout = TFile::Open(ofile);
0590 fout->cd();
0591
0592 TH1I *hi[20];
0593 char name[160], title[160], ctype[20];
0594
0595 for (unsigned int ii=0; ii<=(types.size()); ii++) {
0596 if (ii == 0) sprintf (ctype, "All Particles");
0597 else sprintf (ctype, "%s", types[ii-1].c_str());
0598
0599 sprintf (name, "%s%s%s%sGeV(%s)", ctyp0, target, list, ene, ctype);
0600 hi[ii] = (TH1I*)fout->FindObjectAny(name);
0601
0602
0603 string a(ctype);
0604 means[a] = hi[ii]->GetMean();
0605 }
0606
0607
0608
0609 return means;
0610 }
0611
0612 std::vector<std::string> types() {
0613
0614 std::vector<string> tmp;
0615 tmp.push_back("Photon/Neutrino");
0616 tmp.push_back("e^{-}");
0617 tmp.push_back("e^{+}");
0618 tmp.push_back("#pi^{0}");
0619 tmp.push_back("#pi^{-}");
0620 tmp.push_back("#pi^{+}");
0621 tmp.push_back("K^{-}");
0622 tmp.push_back("K^{+}");
0623 tmp.push_back("K^{0}");
0624 tmp.push_back("AntiProton");
0625 tmp.push_back("p");
0626 tmp.push_back("n");
0627 tmp.push_back("Heavy Hadrons");
0628 tmp.push_back("Ions");
0629
0630 return tmp;
0631 }
0632
0633 std::vector<std::string> typesOld() {
0634
0635 std::vector<string> tmp;
0636 tmp.push_back("Photon/Neutrino");
0637 tmp.push_back("Electron");
0638 tmp.push_back("Positron");
0639 tmp.push_back("Pizero");
0640 tmp.push_back("Piminus");
0641 tmp.push_back("Piplus");
0642 tmp.push_back("Kminus");
0643 tmp.push_back("Kiplus");
0644 tmp.push_back("Kzero");
0645 tmp.push_back("AntiProton");
0646 tmp.push_back("Proton");
0647 tmp.push_back("Neutron/AntiNeutron");
0648 tmp.push_back("Heavy Hadrons");
0649 tmp.push_back("Ions");
0650
0651 return tmp;
0652 }
0653
0654 std::vector<double> massScan() {
0655
0656 std::vector<double> tmp;
0657 tmp.push_back(0.01);
0658 tmp.push_back(1.00);
0659 tmp.push_back(135.0);
0660 tmp.push_back(140.0);
0661 tmp.push_back(495.0);
0662 tmp.push_back(500.0);
0663 tmp.push_back(938.5);
0664 tmp.push_back(940.0);
0665 tmp.push_back(1850.0);
0666 std::cout << tmp.size() << " Mass regions for prtaicles: ";
0667 for (unsigned int i=0; i<tmp.size(); i++) {
0668 std::cout << tmp[i];
0669 if (i == tmp.size()-1) std::cout << " MeV\n";
0670 else std::cout << ", ";
0671 }
0672 return tmp;
0673 }
0674
0675 void setStyle() {
0676
0677 gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0678 gStyle->SetPadColor(kWhite); gStyle->SetFrameBorderMode(0);
0679 gStyle->SetFrameBorderSize(1); gStyle->SetFrameFillColor(0);
0680 gStyle->SetFrameFillStyle(0); gStyle->SetFrameLineColor(1);
0681 gStyle->SetFrameLineStyle(1); gStyle->SetFrameLineWidth(1);
0682 gStyle->SetTitleOffset(1.2,"Y"); gStyle->SetOptStat(0);
0683 gStyle->SetLegendBorderSize(1);
0684
0685 }
