[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPb / macros / production / DrawProduction.C

#include "TCanvas.h"
#include <TH1.h>
#include <TF1.h>
#include <TFile.h>
#include <Rtypes.h>
#include <TLegend.h>
#include <TString.h>

namespace RawProduction {
  class Output;
}

TF1* GetEfficency(const char* trigger, int cent, const char* pid)
{
  if( -1 == cent )
    cent = 0;
  char fileName[256]; sprintf(fileName, "eff_pi0_%s.root", trigger);
  char funcName[256]; sprintf(funcName, "eff_pi0_%s_cen%i", pid, cent);
  
  TFile* file = TFile::Open(fileName, "READ");
  TF1* func = (TF1*)file->Get(funcName);
  return func;
}

// void DrawAbs(const RawProduction::Output& output, const char* trigger, int cent, TStringToken& names) 
// {
//   names.NextToken();
// 
//   const char* pid = "All";
//   char canvName[256] = Form("PSBS_Abs_%s_c%03i_%s_%s", trigger, cent, pid, names.Data());
//   TCanvas* canv = new TCanvas(canvName, canvName);  
//   TLegend* leg = new TLegend(0.6,0.8,0.95,0.95);
// 
//   TH1* hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, pid, names.Data()));
//   if(names.Contains("mr")) hist->SetTitle(Form("Peak Position, %s, %s, %s", trigger, RawProduction::GetCentString(cent), pid));
//   if(names.Contains("sr")) hist->SetTitle(Form("Peak Width, %s, %s, %s", trigger, RawProduction::GetCentString(cent), pid));
//   if(names.Contains("mr")) hist->GetYaxis()->SetRangeUser(0.12, 0.15);
//   if(names.Contains("mr")) hist->GetYaxis()->SetTitle("Peak #mu");
//   if(names.Contains("sr")) hist->GetYaxis()->SetRangeUser(0., 0.012);
//   if(names.Contains("sr")) hist->GetYaxis()->SetTitle("Peak #sigma");
//   hist->GetXaxis()->SetTitle("p_{T}");
//   //Printf(hist->GetTitle());
//   hist->SetMarkerStyle(21);
//   //hist->SetMarkerSize(1.5);
//   hist->SetMarkerColor(kBlack);
//   hist->SetLineColor(kBlack);
//   hist->Draw();
//   leg->AddEntry(hist, "Pol1, Ratio", "lep");
//   
//   int marker = 21;
//   Color_t color[3] = {kRed, kBlue, kMagenta};
//   while( names.NextToken() ) {
//     hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "All", names.Data()));
//     //Printf(hist->GetName());
//     hist->SetMarkerStyle(++marker);
//     hist->SetMarkerColor(color[marker-22]);
//     hist->SetLineColor(color[marker-22]);
//     hist->Draw("same");
//     char legName[256] = "";
//     if(names.Contains("1")) sprintf(legName, "Pol1");
//     if(names.Contains("2")) sprintf(legName, "Pol2");
//     if( marker <23 ) sprintf(legName, "%s, Ratio", legName);
//     leg->AddEntry(hist, legName, "lep");
//   }
//   
//   hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, pid, names.Data()));
//   hist->Draw("same");
// 
//   leg->Draw();
//   
//   canv->SaveAs(Form("imgs/%s.png", canvName));
//   canv->SaveAs(Form("imgs/%s.pdf", canvName));
// }

void DrawRatios(const RawProduction::Output& output, const char* trigger, int cent)
{
  const int nProdTypes = 2;
  enum ProdType { RAW=0, EffCor=1 };
  
  for(int prodType=0; prodType < nProdTypes; ++prodType) {  
    TStringToken graphs("yr1 yr1int yr2 yr2int", " ");
    while ( graphs.NextToken() ) {
      char graphName[32]; sprintf(graphName, "%s", graphs.Data());

      // All
      TH1* hAll = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "All", graphName));
      char newName[256] = Form("E*%s_All", graphName);
      hAll = (TH1*) hAll->Clone(newName);
      hAll->SetTitle(Form("%s, %s, %s, %s", graphName, trigger, RawProduction::GetCentString(cent), graphName));
      hAll->GetXaxis()->SetTitle("IM_{#gamma #gamma}");
      hAll->SetLineColor(kBlack);
      hAll->SetMarkerColor(kBlack);
      hAll->SetMarkerStyle(20);
      hAll->GetXaxis()->SetLabelFont(63); //font in pixels
      hAll->GetXaxis()->SetLabelSize(16); //in pixels
      hAll->GetYaxis()->SetLabelFont(63); //font in pixels
      hAll->GetYaxis()->SetLabelSize(16); //in pixels

      // Allcore
      TH1* hAllcore = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "Allcore", graphName));
      sprintf(newName, "E*%s_Allcore", graphName);
      hAllcore = (TH1*) hAllcore->Clone(newName);
      hAllcore->SetLineColor(kCyan);
      hAllcore->SetMarkerColor(kCyan);
      hAllcore->SetMarkerStyle(25);
  
      // CPV
      TH1* hCPV = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "CPV", graphName));
      sprintf(newName, "E*%s_CPV", graphName);
      hCPV = (TH1*) hCPV->Clone(newName);
      hCPV->SetLineColor(kBlue);
      hCPV->SetMarkerColor(kBlue);
      hCPV->SetMarkerStyle(22);

      // CPVcore
      TH1* hCPVcore = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "CPVcore", graphName));
      sprintf(newName, "E*%s_CPVcore", graphName);
      hCPVcore = (TH1*) hCPVcore->Clone(newName);
      hCPVcore->SetLineColor(kGreen);
      hCPVcore->SetMarkerColor(kGreen);
      hCPVcore->SetMarkerStyle(26);

      // Both
      TH1* hBoth = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "Both", graphName));
      sprintf(newName, "E*%s_Both", graphName);
      hBoth = (TH1*) hBoth->Clone(newName);
      hBoth->SetLineColor(kRed);
      hBoth->SetMarkerColor(kRed);
      hBoth->SetMarkerStyle(33);

      // Efficiancy correction
      if( EffCor == prodType ) {
	hAll->GetYaxis()->SetTitle("#frac{d^{2}N_{#pi^{0}}}{p_{T}dp_{T}dy N_{ev}}");
	hAll->Divide(GetEfficency(trigger, cent, "All"));
	hAllcore->Divide(GetEfficency(trigger, cent, "Allcore"));
	hCPV->Divide(GetEfficency(trigger, cent, "CPV"));
	hCPVcore->Divide(GetEfficency(trigger, cent, "CPVcore"));
	hBoth->Divide(GetEfficency(trigger, cent, "Both"));
      }
      
      // hAll / hAllcore
      TH1* hAll_hAll = hAll->Clone(Form("%s_div_hAll", hAll->GetName()));
      hAll_hAll->Divide(hAll);
      // hAllcore / hAll
      TH1* hAllcore_hAll = hAllcore->Clone(Form("%s_div_hAll", hAllcore->GetName()));
      hAllcore_hAll->Divide(hAll);
      // hCPV / hAll
      TH1* hCPV_hAll = hCPV->Clone(Form("%s_div_hAll", hCPV->GetName()));
      hCPV_hAll->Divide(hAll);
      // hAllcore / hCPVcore
      TH1* hCPVcore_hAllcore = hCPVcore->Clone(Form("%s_div_Allcore", hCPVcore->GetName()));
      hCPVcore_hAllcore->Divide(hAllcore);
      // hBoth / hAll
      TH1* hBoth_hAll = hBoth->Clone(Form("%s_div_hAll", hBoth->GetName()));
      hBoth_hAll->Divide(hAll);

  
      TCanvas* canv = new TCanvas(Form("%s_c%03i_%s_ratio", trigger, cent, graphName), Form("%s_c%03i_%s_ratio", trigger, cent, graphName));
  
      canv->cd();
      TPad *pad1 = new TPad("pad1","pad1",0,0.4,1,1);
      pad1->SetBottomMargin(0);
      pad1->SetLogy();
      pad1->Draw();
      pad1->cd();
      hAll->GetYaxis()->SetTitleSize(0.06);
      hAll->GetYaxis()->SetTitleOffset(0.6);
      hAll->GetYaxis()->SetRangeUser(2.0e-7, 90.);
      hAll->DrawCopy();
      hAllcore->DrawCopy("same");
      hCPV->DrawCopy("same");
      hCPVcore->DrawCopy("same");
      hBoth->DrawCopy("same");
      TLegend* leg1 = new TLegend(0.7,0.6,0.85,0.88);
      leg1->AddEntry(hAll, "All", "lep");
      leg1->AddEntry(hAllcore, "Allcore", "lep");
      leg1->AddEntry(hCPV, "CPV", "lep");
      leg1->AddEntry(hCPVcore, "CPVcore", "lep");
      leg1->AddEntry(hBoth, "Both", "lep");
      leg1->Draw();
  
      canv->cd();
      TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.4);
      pad2->SetTopMargin(0);
      pad2->Draw();
      pad2->cd();
      pad2->SetGridy();
      hAll_hAll->GetYaxis()->SetRangeUser(0.2, 1.4);
      hAll_hAll->GetYaxis()->SetTitle("");
      hAll_hAll->SetTitle("");
      hAll_hAll->DrawCopy("AXIS");
      hAll_hAll->DrawCopy("AXIGsame");
      hAllcore_hAll->DrawCopy("same");
      hCPV_hAll->DrawCopy("same");
      hCPVcore_hAllcore->DrawCopy("same");
      hBoth_hAll->DrawCopy("same");
      TLegend* leg2 = new TLegend(0.7,0.63,0.85,0.98);
      //leg2->AddEntry(hAll_hAll, "All/All", "lep");
      leg2->AddEntry(hAllcore_hAll, "Allcore/All", "lep");
      leg2->AddEntry(hCPV_hAll, "CPV/All", "lep");
      leg2->AddEntry(hCPVcore_hAllcore, "CPVcore/Allcore", "lep");
      leg2->AddEntry(hBoth_hAll, "Both/All", "lep");
      leg2->Draw();
  
      if(RAW == prodType) {
	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s_RAW.pdf", trigger, cent, graphName));      
	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s_RAW.png", trigger, cent, graphName));      
      } else if (EffCor == prodType ) {
	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s.pdf", trigger, cent, graphName));      
	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s.png", trigger, cent, graphName));      
      }
      delete canv;
    }
  }
    
}


void DrawProductions(const RawProduction::Output& output, const char* trigger, int cent)
{
  // TStringToken mrst("mr1r;mr2r;mr1;mr2", ";");
  // DrawAbs(output, trigger, cent, mrst);
    
  DrawRatios(output, trigger, cent);
}

void DrawProduction()
{
  gROOT->LoadMacro("MakeRawProduction.C+g");
  RawProduction::Output output;
  gStyle->SetOptStat(0);
  
  
//   DrawProductions(output, "kMB", -10);
//   DrawProductions(output, "kPHOSPb", -10);

  DrawProductions(output, "kCentral", -1);
//   DrawProductions(output, "kMB", -1);
//   DrawProductions(output, "kPHOSPb", -1);
// 
//   DrawProductions(output, "kSemiCentral", -11);
//   DrawProductions(output, "kMB", -11);
//   DrawProductions(output, "kPHOSPb", -11);
// 
//   DrawProductions(output, "kPHOSPb", -6);
//   DrawProductions(output, "kMB", -6);
}
Commit	Line	Data
f4b109c4	1	#include "TCanvas.h"
	2	#include <TH1.h>
	3	#include <TF1.h>
	4	#include <TFile.h>
	5	#include <Rtypes.h>
	6	#include <TLegend.h>
	7	#include <TString.h>
	8
	9	namespace RawProduction {
	10	class Output;
	11	}
	12
	13	TF1* GetEfficency(const char* trigger, int cent, const char* pid)
	14	{
	15	if( -1 == cent )
	16	cent = 0;
	17	char fileName[256]; sprintf(fileName, "eff_pi0_%s.root", trigger);
	18	char funcName[256]; sprintf(funcName, "eff_pi0_%s_cen%i", pid, cent);
	19
	20	TFile* file = TFile::Open(fileName, "READ");
	21	TF1* func = (TF1*)file->Get(funcName);
	22	return func;
	23	}
	24
	25	// void DrawAbs(const RawProduction::Output& output, const char* trigger, int cent, TStringToken& names)
	26	// {
	27	// names.NextToken();
	28	//
	29	// const char* pid = "All";
	30	// char canvName[256] = Form("PSBS_Abs_%s_c%03i_%s_%s", trigger, cent, pid, names.Data());
	31	// TCanvas* canv = new TCanvas(canvName, canvName);
	32	// TLegend* leg = new TLegend(0.6,0.8,0.95,0.95);
	33	//
	34	// TH1* hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, pid, names.Data()));
	35	// if(names.Contains("mr")) hist->SetTitle(Form("Peak Position, %s, %s, %s", trigger, RawProduction::GetCentString(cent), pid));
	36	// if(names.Contains("sr")) hist->SetTitle(Form("Peak Width, %s, %s, %s", trigger, RawProduction::GetCentString(cent), pid));
	37	// if(names.Contains("mr")) hist->GetYaxis()->SetRangeUser(0.12, 0.15);
	38	// if(names.Contains("mr")) hist->GetYaxis()->SetTitle("Peak #mu");
	39	// if(names.Contains("sr")) hist->GetYaxis()->SetRangeUser(0., 0.012);
	40	// if(names.Contains("sr")) hist->GetYaxis()->SetTitle("Peak #sigma");
	41	// hist->GetXaxis()->SetTitle("p_{T}");
	42	// //Printf(hist->GetTitle());
	43	// hist->SetMarkerStyle(21);
	44	// //hist->SetMarkerSize(1.5);
	45	// hist->SetMarkerColor(kBlack);
	46	// hist->SetLineColor(kBlack);
	47	// hist->Draw();
	48	// leg->AddEntry(hist, "Pol1, Ratio", "lep");
	49	//
	50	// int marker = 21;
	51	// Color_t color[3] = {kRed, kBlue, kMagenta};
	52	// while( names.NextToken() ) {
	53	// hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "All", names.Data()));
	54	// //Printf(hist->GetName());
	55	// hist->SetMarkerStyle(++marker);
	56	// hist->SetMarkerColor(color[marker-22]);
	57	// hist->SetLineColor(color[marker-22]);
	58	// hist->Draw("same");
	59	// char legName[256] = "";
	60	// if(names.Contains("1")) sprintf(legName, "Pol1");
	61	// if(names.Contains("2")) sprintf(legName, "Pol2");
	62	// if( marker <23 ) sprintf(legName, "%s, Ratio", legName);
	63	// leg->AddEntry(hist, legName, "lep");
	64	// }
65	//
66	// hist = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, pid, names.Data()));
67	// hist->Draw("same");
68	//
69	// leg->Draw();
70	//
71	// canv->SaveAs(Form("imgs/%s.png", canvName));
72	// canv->SaveAs(Form("imgs/%s.pdf", canvName));
73	// }
74
75	void DrawRatios(const RawProduction::Output& output, const char* trigger, int cent)
76	{
77	const int nProdTypes = 2;
78	enum ProdType { RAW=0, EffCor=1 };
79
80	for(int prodType=0; prodType < nProdTypes; ++prodType) {
81	TStringToken graphs("yr1 yr1int yr2 yr2int", " ");
82	while ( graphs.NextToken() ) {
83	char graphName[32]; sprintf(graphName, "%s", graphs.Data());
84
85	// All
86	TH1* hAll = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "All", graphName));
87	char newName[256] = Form("E*%s_All", graphName);
88	hAll = (TH1*) hAll->Clone(newName);
89	hAll->SetTitle(Form("%s, %s, %s, %s", graphName, trigger, RawProduction::GetCentString(cent), graphName));
90	hAll->GetXaxis()->SetTitle("IM_{#gamma #gamma}");
91	hAll->SetLineColor(kBlack);
92	hAll->SetMarkerColor(kBlack);
93	hAll->SetMarkerStyle(20);
94	hAll->GetXaxis()->SetLabelFont(63); //font in pixels
95	hAll->GetXaxis()->SetLabelSize(16); //in pixels
96	hAll->GetYaxis()->SetLabelFont(63); //font in pixels
97	hAll->GetYaxis()->SetLabelSize(16); //in pixels
98
99	// Allcore
100	TH1* hAllcore = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "Allcore", graphName));
101	sprintf(newName, "E*%s_Allcore", graphName);
102	hAllcore = (TH1*) hAllcore->Clone(newName);
103	hAllcore->SetLineColor(kCyan);
104	hAllcore->SetMarkerColor(kCyan);
105	hAllcore->SetMarkerStyle(25);
106
107	// CPV
108	TH1* hCPV = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "CPV", graphName));
109	sprintf(newName, "E*%s_CPV", graphName);
110	hCPV = (TH1*) hCPV->Clone(newName);
111	hCPV->SetLineColor(kBlue);
112	hCPV->SetMarkerColor(kBlue);
113	hCPV->SetMarkerStyle(22);
114
115	// CPVcore
116	TH1* hCPVcore = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "CPVcore", graphName));
117	sprintf(newName, "E*%s_CPVcore", graphName);
118	hCPVcore = (TH1*) hCPVcore->Clone(newName);
119	hCPVcore->SetLineColor(kGreen);
120	hCPVcore->SetMarkerColor(kGreen);
121	hCPVcore->SetMarkerStyle(26);
122
123	// Both
124	TH1* hBoth = output.GetHistogram(Form("%s/c%03i/%s/%s", trigger, cent, "Both", graphName));
125	sprintf(newName, "E*%s_Both", graphName);
126	hBoth = (TH1*) hBoth->Clone(newName);
127	hBoth->SetLineColor(kRed);
128	hBoth->SetMarkerColor(kRed);
129	hBoth->SetMarkerStyle(33);
130
131	// Efficiancy correction
132	if( EffCor == prodType ) {
133	hAll->GetYaxis()->SetTitle("#frac{d^{2}N_{#pi^{0}}}{p_{T}dp_{T}dy N_{ev}}");
134	hAll->Divide(GetEfficency(trigger, cent, "All"));
135	hAllcore->Divide(GetEfficency(trigger, cent, "Allcore"));
136	hCPV->Divide(GetEfficency(trigger, cent, "CPV"));
137	hCPVcore->Divide(GetEfficency(trigger, cent, "CPVcore"));
138	hBoth->Divide(GetEfficency(trigger, cent, "Both"));
139	}
140
141	// hAll / hAllcore
142	TH1* hAll_hAll = hAll->Clone(Form("%s_div_hAll", hAll->GetName()));
143	hAll_hAll->Divide(hAll);
144	// hAllcore / hAll
145	TH1* hAllcore_hAll = hAllcore->Clone(Form("%s_div_hAll", hAllcore->GetName()));
146	hAllcore_hAll->Divide(hAll);
147	// hCPV / hAll
148	TH1* hCPV_hAll = hCPV->Clone(Form("%s_div_hAll", hCPV->GetName()));
149	hCPV_hAll->Divide(hAll);
150	// hAllcore / hCPVcore
151	TH1* hCPVcore_hAllcore = hCPVcore->Clone(Form("%s_div_Allcore", hCPVcore->GetName()));
152	hCPVcore_hAllcore->Divide(hAllcore);
153	// hBoth / hAll
154	TH1* hBoth_hAll = hBoth->Clone(Form("%s_div_hAll", hBoth->GetName()));
155	hBoth_hAll->Divide(hAll);
156
157
158	TCanvas* canv = new TCanvas(Form("%s_c%03i_%s_ratio", trigger, cent, graphName), Form("%s_c%03i_%s_ratio", trigger, cent, graphName));
159
160	canv->cd();
161	TPad *pad1 = new TPad("pad1","pad1",0,0.4,1,1);
162	pad1->SetBottomMargin(0);
163	pad1->SetLogy();
164	pad1->Draw();
165	pad1->cd();
166	hAll->GetYaxis()->SetTitleSize(0.06);
167	hAll->GetYaxis()->SetTitleOffset(0.6);
168	hAll->GetYaxis()->SetRangeUser(2.0e-7, 90.);
169	hAll->DrawCopy();
170	hAllcore->DrawCopy("same");
171	hCPV->DrawCopy("same");
172	hCPVcore->DrawCopy("same");
173	hBoth->DrawCopy("same");
174	TLegend* leg1 = new TLegend(0.7,0.6,0.85,0.88);
175	leg1->AddEntry(hAll, "All", "lep");
176	leg1->AddEntry(hAllcore, "Allcore", "lep");
177	leg1->AddEntry(hCPV, "CPV", "lep");
178	leg1->AddEntry(hCPVcore, "CPVcore", "lep");
179	leg1->AddEntry(hBoth, "Both", "lep");
180	leg1->Draw();
181
182	canv->cd();
183	TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.4);
184	pad2->SetTopMargin(0);
185	pad2->Draw();
186	pad2->cd();
187	pad2->SetGridy();
188	hAll_hAll->GetYaxis()->SetRangeUser(0.2, 1.4);
189	hAll_hAll->GetYaxis()->SetTitle("");
190	hAll_hAll->SetTitle("");
191	hAll_hAll->DrawCopy("AXIS");
192	hAll_hAll->DrawCopy("AXIGsame");
193	hAllcore_hAll->DrawCopy("same");
194	hCPV_hAll->DrawCopy("same");
195	hCPVcore_hAllcore->DrawCopy("same");
196	hBoth_hAll->DrawCopy("same");
197	TLegend* leg2 = new TLegend(0.7,0.63,0.85,0.98);
198	//leg2->AddEntry(hAll_hAll, "All/All", "lep");
199	leg2->AddEntry(hAllcore_hAll, "Allcore/All", "lep");
200	leg2->AddEntry(hCPV_hAll, "CPV/All", "lep");
201	leg2->AddEntry(hCPVcore_hAllcore, "CPVcore/Allcore", "lep");
202	leg2->AddEntry(hBoth_hAll, "Both/All", "lep");
203	leg2->Draw();
204
205	if(RAW == prodType) {
206	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s_RAW.pdf", trigger, cent, graphName));
207	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s_RAW.png", trigger, cent, graphName));
208	} else if (EffCor == prodType ) {
209	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s.pdf", trigger, cent, graphName));
210	canv->SaveAs(Form("imgs/prod_ratios_%s_%i_%s.png", trigger, cent, graphName));
211	}
212	delete canv;
213	}
214	}
215
216	}
217
218
219	void DrawProductions(const RawProduction::Output& output, const char* trigger, int cent)
220	{
221	// TStringToken mrst("mr1r;mr2r;mr1;mr2", ";");
222	// DrawAbs(output, trigger, cent, mrst);
223
224	DrawRatios(output, trigger, cent);
225	}
226
227	void DrawProduction()
228	{
229	gROOT->LoadMacro("MakeRawProduction.C+g");
230	RawProduction::Output output;
231	gStyle->SetOptStat(0);
232
233
234	// DrawProductions(output, "kMB", -10);
235	// DrawProductions(output, "kPHOSPb", -10);
236
237	DrawProductions(output, "kCentral", -1);
238	// DrawProductions(output, "kMB", -1);
239	// DrawProductions(output, "kPHOSPb", -1);
240	//
241	// DrawProductions(output, "kSemiCentral", -11);
242	// DrawProductions(output, "kMB", -11);
243	// DrawProductions(output, "kPHOSPb", -11);
244	//
245	// DrawProductions(output, "kPHOSPb", -6);
246	// DrawProductions(output, "kMB", -6);
247	}