Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:39

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   //  std::cout << "Number of types: " << types.size() << "\n";
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     // std::cout<<"ii "<<ii<<"  ctype "<<ctype<<std::endl;
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     //    std::cout << ii << " (" << ctype << ") " << name << " " << hiMulti[ii] << "\n";
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   //  std::cout << "Number of types " << types.size() << "\n";
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     //    std::cout<<a<<" "<< means_300[a]<<std::endl;
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     //    std::cout << ii << " (" << ctype << ") " << name << " " << hiKE[ii] <<"\n";
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     //    std::cout << "Histo " << ii << " Name " << name << " " << hi[ii] << " " << hi[ii]->GetMean() << "\n";
0602     
0603     string a(ctype);
0604     means[a] = hi[ii]->GetMean();
0605   }
0606 
0607   //  printMeans(means);
0608 
0609   return means;
0610 }
0611 
0612 std::vector<std::string> types() {
0613 
0614   std::vector<string> tmp;
0615   tmp.push_back("Photon/Neutrino");     // 1
0616   tmp.push_back("e^{-}");               // 2
0617   tmp.push_back("e^{+}");               // 3 
0618   tmp.push_back("#pi^{0}");             // 4 
0619   tmp.push_back("#pi^{-}");             // 5
0620   tmp.push_back("#pi^{+}");             // 6
0621   tmp.push_back("K^{-}");               // 7
0622   tmp.push_back("K^{+}");               // 8
0623   tmp.push_back("K^{0}");               // 9
0624   tmp.push_back("AntiProton");          // 10
0625   tmp.push_back("p");                   // 11
0626   tmp.push_back("n");                   // 12
0627   tmp.push_back("Heavy Hadrons");       // 13
0628   tmp.push_back("Ions");                // 14
0629 
0630   return tmp;
0631 }
0632 
0633 std::vector<std::string> typesOld() {
0634 
0635   std::vector<string> tmp;
0636   tmp.push_back("Photon/Neutrino");     // 1
0637   tmp.push_back("Electron");            // 2
0638   tmp.push_back("Positron");            // 3 
0639   tmp.push_back("Pizero");              // 4 
0640   tmp.push_back("Piminus");             // 5
0641   tmp.push_back("Piplus");              // 6
0642   tmp.push_back("Kminus");              // 7
0643   tmp.push_back("Kiplus");              // 8
0644   tmp.push_back("Kzero");               // 9
0645   tmp.push_back("AntiProton");          // 10
0646   tmp.push_back("Proton");              // 11
0647   tmp.push_back("Neutron/AntiNeutron"); // 12
0648   tmp.push_back("Heavy Hadrons");       // 13
0649   tmp.push_back("Ions");                // 14
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 }