void genHistos(TString inFilename = "mc_hitcount.root", TString geantversion = "undefined") {
	bool savePictures = true;
	TString macroBasename = "genHistos.C";
	bool normalize = true;

	TFile * inFile = new TFile(inFilename);

	TTree * hitCounts = (TTree*) (inFile->Get("hitCount"));

	int colorChargePos(kBlue), colorChargeNeg(kRed);

	int nSimulatedPiPlus = hitCounts->GetEntries("PdgId == 211");
	int nSimulatedPiMinus = hitCounts->GetEntries("PdgId == -211");
	std::cout << "Simulated Pi+: " <<  nSimulatedPiPlus << std::endl;
	std::cout << "Simulated Pi-: " <<  nSimulatedPiMinus << std::endl;

	TH1D * hPipLayer1 = new TH1D("hPipLayer1", "#pi^{+}: #Hits in MVD Layer 1;# Hits In Layer 1;counts", 5, -0.5, 4.5);
	TH1D * hPimLayer1 = new TH1D("hPimLayer1", "#pi^{-}: #Hits in MVD Layer 1;# Hits In Layer 1;counts", 5, -0.5, 4.5);
	TH1D * hPipLayer2 = new TH1D("hPipLayer2", "#pi^{+}: #Hits in MVD Layer 2;# Hits In Layer 2;counts", 5, -0.5, 4.5);
	TH1D * hPimLayer2 = new TH1D("hPimLayer2", "#pi^{-}: #Hits in MVD Layer 2;# Hits In Layer 2;counts", 5, -0.5, 4.5);
	TH1D * hPipLayer3 = new TH1D("hPipLayer3", "#pi^{+}: #Hits in MVD Layer 3;# Hits In Layer 3;counts", 5, -0.5, 4.5);
	TH1D * hPimLayer3 = new TH1D("hPimLayer3", "#pi^{-}: #Hits in MVD Layer 3;# Hits In Layer 3;counts", 5, -0.5, 4.5);
	TH1D * hPipLayer4 = new TH1D("hPipLayer4", "#pi^{+}: #Hits in MVD Layer 4;# Hits In Layer 4;counts", 5, -0.5, 4.5);
	TH1D * hPimLayer4 = new TH1D("hPimLayer4", "#pi^{-}: #Hits in MVD Layer 4;# Hits In Layer 4;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk1 = new TH1D("hPipLayerDisk1", "#pi^{+}: #Hits in MVD Disk Layer 1;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk2 = new TH1D("hPipLayerDisk2", "#pi^{+}: #Hits in MVD Disk Layer 2;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk3 = new TH1D("hPipLayerDisk3", "#pi^{+}: #Hits in MVD Disk Layer 3;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk4 = new TH1D("hPipLayerDisk4", "#pi^{+}: #Hits in MVD Disk Layer 4;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk5 = new TH1D("hPipLayerDisk5", "#pi^{+}: #Hits in MVD Disk Layer 5;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk6 = new TH1D("hPipLayerDisk6", "#pi^{+}: #Hits in MVD Disk Layer 6;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk7 = new TH1D("hPipLayerDisk7", "#pi^{+}: #Hits in MVD Disk Layer 7;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPipLayerDisk8 = new TH1D("hPipLayerDisk8", "#pi^{+}: #Hits in MVD Disk Layer 8;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk1 = new TH1D("hPimLayerDisk1", "#pi^{-}: #Hits in MVD Disk Layer 1;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk2 = new TH1D("hPimLayerDisk2", "#pi^{-}: #Hits in MVD Disk Layer 2;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk3 = new TH1D("hPimLayerDisk3", "#pi^{-}: #Hits in MVD Disk Layer 3;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk4 = new TH1D("hPimLayerDisk4", "#pi^{-}: #Hits in MVD Disk Layer 4;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk5 = new TH1D("hPimLayerDisk5", "#pi^{-}: #Hits in MVD Disk Layer 5;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk6 = new TH1D("hPimLayerDisk6", "#pi^{-}: #Hits in MVD Disk Layer 6;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk7 = new TH1D("hPimLayerDisk7", "#pi^{-}: #Hits in MVD Disk Layer 7;# Hits In Layer;counts", 5, -0.5, 4.5);
	TH1D * hPimLayerDisk8 = new TH1D("hPimLayerDisk8", "#pi^{-}: #Hits in MVD Disk Layer 8;# Hits In Layer;counts", 5, -0.5, 4.5);

	hitCounts->Project("hPipLayer1", "HitInLayer1","PdgId == 211 && HitInLayer1 > 0");
	hitCounts->Project("hPimLayer1", "HitInLayer1","PdgId == -211 && HitInLayer1 > 0");
	hitCounts->Project("hPipLayer2", "HitInLayer2","PdgId == 211 && HitInLayer2 > 0");
	hitCounts->Project("hPimLayer2", "HitInLayer2","PdgId == -211 && HitInLayer2 > 0");
	hitCounts->Project("hPipLayer3", "HitInLayer3","PdgId == 211 && HitInLayer3 > 0");
	hitCounts->Project("hPimLayer3", "HitInLayer3","PdgId == -211 && HitInLayer3 > 0");
	hitCounts->Project("hPipLayer4", "HitInLayer4","PdgId == 211 && HitInLayer4 > 0");
	hitCounts->Project("hPimLayer4", "HitInLayer4","PdgId == -211 && HitInLayer4 > 0");
	hitCounts->Project("hPipLayerDisk1", "HitInDisk1","PdgId == 211 && HitInDisk1 > 0");
	hitCounts->Project("hPipLayerDisk2", "HitInDisk2","PdgId == 211 && HitInDisk2 > 0");
	hitCounts->Project("hPipLayerDisk3", "HitInDisk3","PdgId == 211 && HitInDisk3 > 0");
	hitCounts->Project("hPipLayerDisk4", "HitInDisk4","PdgId == 211 && HitInDisk4 > 0");
	hitCounts->Project("hPipLayerDisk5", "HitInDisk5","PdgId == 211 && HitInDisk5 > 0");
	hitCounts->Project("hPipLayerDisk6", "HitInDisk6","PdgId == 211 && HitInDisk6 > 0");
	hitCounts->Project("hPipLayerDisk7", "HitInDisk7","PdgId == 211 && HitInDisk7 > 0");
	hitCounts->Project("hPipLayerDisk8", "HitInDisk8","PdgId == 211 && HitInDisk8 > 0");
	hitCounts->Project("hPimLayerDisk1", "HitInDisk1","PdgId == -211 && HitInDisk1 > 0");
	hitCounts->Project("hPimLayerDisk2", "HitInDisk2","PdgId == -211 && HitInDisk2 > 0");
	hitCounts->Project("hPimLayerDisk3", "HitInDisk3","PdgId == -211 && HitInDisk3 > 0");
	hitCounts->Project("hPimLayerDisk4", "HitInDisk4","PdgId == -211 && HitInDisk4 > 0");
	hitCounts->Project("hPimLayerDisk5", "HitInDisk5","PdgId == -211 && HitInDisk5 > 0");
	hitCounts->Project("hPimLayerDisk6", "HitInDisk6","PdgId == -211 && HitInDisk6 > 0");
	hitCounts->Project("hPimLayerDisk7", "HitInDisk7","PdgId == -211 && HitInDisk7 > 0");
	hitCounts->Project("hPimLayerDisk8", "HitInDisk8","PdgId == -211 && HitInDisk8 > 0");

	if (normalize) {
		hPipLayer1->Scale(1/(double)nSimulatedPiPlus);
		hPimLayer1->Scale(1/(double)nSimulatedPiMinus);
		hPipLayer2->Scale(1/(double)nSimulatedPiPlus);
		hPimLayer2->Scale(1/(double)nSimulatedPiMinus);
		hPipLayer3->Scale(1/(double)nSimulatedPiPlus);
		hPimLayer3->Scale(1/(double)nSimulatedPiMinus);
		hPipLayer4->Scale(1/(double)nSimulatedPiPlus);
		hPimLayer4->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk1->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk1->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk2->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk2->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk3->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk3->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk4->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk4->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk5->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk5->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk6->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk6->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk7->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk7->Scale(1/(double)nSimulatedPiMinus);
		hPipLayerDisk8->Scale(1/(double)nSimulatedPiPlus);
		hPimLayerDisk8->Scale(1/(double)nSimulatedPiMinus);
		hPipLayer1->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayer1->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayer2->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayer2->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayer3->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayer3->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayer4->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayer4->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk1->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk1->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk2->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk2->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk3->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk3->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk4->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk4->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk5->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk5->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk6->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk6->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk7->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk7->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPipLayerDisk8->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayerDisk8->GetYaxis()->SetTitle("Ratio (a.u.)");
	}


	TH1D * hPipLayers = new TH1D("hPipLayers", "#pi^{+}: Hits in Layers (#beq 1 hit in layer);MVD Layer Name;Counts", 12, 0.5, 12.5);
	TH1D * hPimLayers = new TH1D("hPimLayers", "#pi^{-}: Hits in Layers (#beq 1 hit in layer);MVD Layer Name;Counts", 12, 0.5, 12.5);
	hPipLayers->SetLineColor(colorChargePos);
	hPimLayers->SetLineColor(colorChargeNeg);
	hPipLayers->SetBinContent(1, hPipLayer1->GetEntries());
	hPipLayers->SetBinContent(2, hPipLayer2->GetEntries());
	hPipLayers->SetBinContent(3, hPipLayer3->GetEntries());
	hPipLayers->SetBinContent(4, hPipLayer4->GetEntries());
	hPipLayers->SetBinContent(5, hPipLayerDisk1->GetEntries());
	hPipLayers->SetBinContent(6, hPipLayerDisk2->GetEntries());
	hPipLayers->SetBinContent(7, hPipLayerDisk3->GetEntries());
	hPipLayers->SetBinContent(8, hPipLayerDisk4->GetEntries());
	hPipLayers->SetBinContent(9, hPipLayerDisk5->GetEntries());
	hPipLayers->SetBinContent(10, hPipLayerDisk7->GetEntries());
	hPipLayers->SetBinContent(11, hPipLayerDisk6->GetEntries());
	hPipLayers->SetBinContent(12, hPipLayerDisk8->GetEntries());
	
	hPimLayers->SetBinContent(1, hPimLayer1->GetEntries());
	hPimLayers->SetBinContent(2, hPimLayer2->GetEntries());
	hPimLayers->SetBinContent(3, hPimLayer3->GetEntries());
	hPimLayers->SetBinContent(4, hPimLayer4->GetEntries());
	hPimLayers->SetBinContent(5, hPimLayerDisk1->GetEntries());
	hPimLayers->SetBinContent(6, hPimLayerDisk2->GetEntries());
	hPimLayers->SetBinContent(7, hPimLayerDisk3->GetEntries());
	hPimLayers->SetBinContent(8, hPimLayerDisk4->GetEntries());
	hPimLayers->SetBinContent(9, hPimLayerDisk5->GetEntries());
	hPimLayers->SetBinContent(10, hPimLayerDisk7->GetEntries());
	hPimLayers->SetBinContent(11, hPimLayerDisk6->GetEntries());
	hPimLayers->SetBinContent(12, hPimLayerDisk8->GetEntries());

	if (normalize) {
		hPipLayers->Scale(1/(double)nSimulatedPiPlus);
		hPimLayers->Scale(1/(double)nSimulatedPiMinus);
		hPipLayers->GetYaxis()->SetTitle("Ratio (a.u.)");
		hPimLayers->GetYaxis()->SetTitle("Ratio (a.u.)");
	}

	hPipLayers->SetTitle("#pi^{+}");
	hPimLayers->SetTitle("#pi^{-}");

	TCanvas * cPiLayers = new TCanvas("cPiLayers", "Hit Count Asymmetry in Layers", 0, 0, 800, 500);
	THStack * stPiLayers = new THStack();
	stPiLayers->SetName("stPiLayers");
	stPiLayers->SetTitle("#pi^{#pm}: #geq 1 Hit in Different MVD Layers");
	stPiLayers->Add(hPipLayers);
	stPiLayers->Add(hPimLayers);
	stPiLayers->Draw("NOSTACK");
	stPiLayers->GetXaxis()->SetTitle("MVD Layer Name");
	stPiLayers->GetXaxis()->SetTitleOffset(1.02);
	stPiLayers->GetYaxis()->SetTitle("Counts");
	if (normalize) stPiLayers->GetYaxis()->SetTitle("Ratio (a.u.)");

	stPiLayers->GetXaxis()->SetBinLabel(1, "Barrel 1");
	stPiLayers->GetXaxis()->SetBinLabel(2, "Barrel 2");
	stPiLayers->GetXaxis()->SetBinLabel(3, "Barrel 3");  // S
	stPiLayers->GetXaxis()->SetBinLabel(4, "Barrel 4");  // S
	stPiLayers->GetXaxis()->SetBinLabel(5, "Disk 1");
	stPiLayers->GetXaxis()->SetBinLabel(6, "Disk 2");
	stPiLayers->GetXaxis()->SetBinLabel(7, "Disk 3");
	stPiLayers->GetXaxis()->SetBinLabel(8, "Disk 4");
	stPiLayers->GetXaxis()->SetBinLabel(9, "Disk 5P");
	stPiLayers->GetXaxis()->SetBinLabel(10, "Disk 5S");
	stPiLayers->GetXaxis()->SetBinLabel(11, "Disk 6P");
	stPiLayers->GetXaxis()->SetBinLabel(12, "Disk 6S");

	TLegend * legPiLayers = andi::stackLegend(stPiLayers, geantversion);
	legPiLayers->Draw();

	cPiLayers->SaveAs("pion_asymmetry_" + geantversion + ".pdf");

}

TLegend * plainLegend(float x1 = 0.7, float y1 = 0.62, float x2 = 0.96, float y2 = 0.695) {
	TLegend * legend = new TLegend(x1, y1, x2, y2);
	legend->SetFillColor(kWhite);
	legend->SetTextSize(0.04);
	legend->SetEntrySeparation(0.1);
	legend->SetBorderSize(1);
}

TLegend * stackLegend(THStack * stack, TString header = "Histograms", double height = 0.75, double width = 0.7, double y2 = 0.93) {
	TLegend * tempLegend = plainLegend(width, height, 0.96, y2);
	tempLegend->SetHeader(header);
	int nHists = stack->GetHists()->GetEntries();
	for (int i = 0; i < nHists; i++) {
		tempLegend->AddEntry(stack->GetHists()->At(i), stack->GetHists()->At(i)->GetTitle(), "LPF");
	}
	tempLegend->SetTextFont(62);
	tempLegend->SetTextSize(0.038);  // empirically found

	return tempLegend;
}