[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / macros / Plot_muons.C

#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <Riostream.h>
#include <assert.h>

#include "TVector2.h"
#include "TFile.h"
#include "TString.h"
#include "TF1.h"
#include "TF3.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TMath.h"
#include "TText.h"
#include "TRandom3.h"
#include "TArray.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TPaveStats.h"
#include "TASImage.h"

#define BohrR 1963.6885 // Bohr Radius for pions
#define FmToGeV 0.19733 // conversion to fm
#define PI 3.1415926
#define masspiC 0.1395702 // pi+ mass (GeV/c^2)
#define kappa3 0.16
#define kappa4 0.40

using namespace std;


void Plot_muons(){
  
 
  gStyle->SetOptStat(0);
  gStyle->SetOptDate(0);
  gStyle->SetOptFit(0111);
  
  int ChComb=1;// 0=Same-charge, 1=Mixed-charge
  int kTbin_L=5, kTbin_H=6;
  int binKT3=1;//1-2

  //TFile *_file0= new TFile("Results/PDC_12a17a_muons_R6.root","READ");
  //TFile *_file0= new TFile("Results/PDC_13b2_efix_p1_muons_R2.root","READ");
  TFile *_file0= new TFile("Results/PDC_12a17e_muons_R4.root","READ");
  
  TList *MyList=(TList*)_file0->Get("MyList");
  _file0->Close();

  TH1D *PionCandidates=(TH1D*)MyList->FindObject("fPionCandidates");
  PionCandidates->GetXaxis()->SetTitle("PDG code");
  //PionCandidates->Draw();
  //
  TH1D *MuonParentsPrimary=(TH1D*)MyList->FindObject("fMuonParents");
  MuonParentsPrimary->GetXaxis()->SetTitle("PDG code");
  MuonParentsPrimary->SetFillColor(1);
  //MuonParentsPrimary->Draw();
  //
  TH1D *MuonParentsSecondary=(TH1D*)MyList->FindObject("fSecondaryMuonParents");
  MuonParentsSecondary->GetXaxis()->SetTitle("PDG code");
  MuonParentsSecondary->SetFillColor(1);
  //MuonParentsSecondary->Draw();
  //
  // M0 R10-R6, M6 for R4, M17 for R2
  TH3D *PurityNum_3D = (TH3D*)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityNum");
  TH2D *PurityDen_2D = (TH2D*)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityDen");
  TH1D *PurityNum=PurityNum_3D->ProjectionX("PurityNum",kTbin_L,kTbin_H,1,20);
  double PurityNorm=PurityDen_2D->Integral(kTbin_L,kTbin_H,1,20);
  PurityNum->Scale(1/PurityNorm);
  char *namesAxis[15]={"e-e","e-mu","e-pi","e-k","e-p","mu-mu","mu-pi","mu-k","mu-p","pi-pi","pi-k","pi-p","k-k","k-p","p-p"};
  for(int i=1; i<=15; i++) PurityNum->GetXaxis()->SetBinLabel(i, namesAxis[i-1]);
  PurityNum->GetXaxis()->SetRange(1,15);
  PurityNum->GetYaxis()->SetTitle("Probability");
  PurityNum->Draw();
  //
  //
  TCanvas *can = new TCanvas("can", "can",800,0,800,800);// 11,53,700,500
  can->SetHighLightColor(2);
  gStyle->SetOptFit(0111);
  can->SetFillColor(10);//10
  can->SetBorderMode(0);
  can->SetBorderSize(2);
  can->SetFrameFillColor(0);
  can->SetFrameBorderMode(0);
  can->SetFrameBorderMode(0);
  can->cd();
  TPad *pad = new TPad("pad","pad",0.,0.,1.,1.);
  gPad->SetTickx();
  gPad->SetTicky();
  pad->SetGridx();
  pad->SetGridy();
  pad->SetTopMargin(0.02);//0.05
  pad->SetRightMargin(0.02);//3e-2
  pad->SetBottomMargin(0.1);//0.12
  pad->SetLeftMargin(0.1);
  pad->Draw();
  pad->cd();
  TLegend *legend = new TLegend(.5,.65, .9,.95,NULL,"brNDC");//.45 or .4 for x1
  legend->SetBorderSize(0);
  legend->SetFillColor(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.03);
  //
  TH3D *MuonSmearedNum2_3=(TH3D*)MyList->FindObject("fMuonContamSmearedNum2");
  TH3D *MuonSmearedDen2_3=(TH3D*)MyList->FindObject("fMuonContamSmearedDen2");
  TH3D *PionNum2_3=(TH3D*)MyList->FindObject("fMuonContamIdealNum2");
  TH3D *PionDen2_3=(TH3D*)MyList->FindObject("fMuonContamIdealDen2");
  TH3D *PionPionK2_3=(TH3D*)MyList->FindObject("fPionPionK2");
  //
  TH3D *MuonSmearedNum3_3=(TH3D*)MyList->FindObject("fMuonContamSmearedNum3");
  TH3D *MuonSmearedDen3_3=(TH3D*)MyList->FindObject("fMuonContamSmearedDen3");
  TH3D *PionNum3_3=(TH3D*)MyList->FindObject("fMuonContamIdealNum3");
  TH3D *PionDen3_3=(TH3D*)MyList->FindObject("fMuonContamIdealDen3");
  TH3D *PionPionK3_3=(TH3D*)MyList->FindObject("fPionPionK3");
  TH3D *MuonPionK3_3=(TH3D*)MyList->FindObject("fMuonPionK3");
  
  
  TH1D *MuonSmearedNum2[2];// SC/MC  
  TH1D *MuonSmearedDen2[2];
  TH1D *PionNum2[2];
  TH1D *PionDen2[2];
  TH1D *PionPionK2[2];
  //
  TH1D *MuonSmearedNum3[2];// SC/MC  
  TH1D *MuonSmearedDen3[2];
  TH1D *PionNum3[2];
  TH1D *PionDen3[2];
  TH1D *PionPionK3[2];
  //
  TH1D *C2muonSmeared[2];
  TH1D *C2pion[2];
  TH1D *C3muonSmeared[2];
  TH1D *C3pion[2];
  //
  for(int chtype=0; chtype<2; chtype++){
    TString *names[10];
    for(int i=0; i<10; i++) {names[i]=new TString("name_"); *names[i] += i; *names[i] += chtype;}

    MuonSmearedNum2[chtype]=(TH1D*)MuonSmearedNum2_3->ProjectionZ(names[0]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
    MuonSmearedDen2[chtype]=(TH1D*)MuonSmearedDen2_3->ProjectionZ(names[1]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
    PionNum2[chtype]=(TH1D*)PionNum2_3->ProjectionZ(names[2]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
    PionDen2[chtype]=(TH1D*)PionDen2_3->ProjectionZ(names[3]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
    PionPionK2[chtype]=(TH1D*)PionPionK2_3->ProjectionZ(names[4]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
    PionPionK2[chtype]->Divide(PionDen2[chtype]);
    ////////////////
    MuonSmearedNum3[chtype]=(TH1D*)MuonSmearedNum3_3->ProjectionZ(names[5]->Data(),chtype+1,chtype+1,binKT3,binKT3);
    MuonSmearedDen3[chtype]=(TH1D*)MuonSmearedDen3_3->ProjectionZ(names[6]->Data(),chtype+1,chtype+1,binKT3,binKT3);
    PionNum3[chtype]=(TH1D*)PionNum3_3->ProjectionZ(names[7]->Data(),chtype+1,chtype+1,binKT3,binKT3);
    PionDen3[chtype]=(TH1D*)PionDen3_3->ProjectionZ(names[8]->Data(),chtype+1,chtype+1,binKT3,binKT3);
    PionPionK3[chtype]=(TH1D*)PionPionK3_3->ProjectionZ(names[9]->Data(),chtype+1,chtype+1,binKT3,binKT3);
    PionPionK3[chtype]->Divide(PionDen3[chtype]);
    //
    C2muonSmeared[chtype]=(TH1D*)MuonSmearedNum2[chtype]->Clone();
    C2pion[chtype]=(TH1D*)PionNum2[chtype]->Clone();
    C2muonSmeared[chtype]->Divide(MuonSmearedDen2[chtype]);
    C2pion[chtype]->Divide(PionDen2[chtype]);
    //
    C3muonSmeared[chtype]=(TH1D*)MuonSmearedNum3[chtype]->Clone();
    C3pion[chtype]=(TH1D*)PionNum3[chtype]->Clone();
    C3muonSmeared[chtype]->Divide(MuonSmearedDen3[chtype]);
    C3pion[chtype]->Divide(PionDen3[chtype]);
    //
    //
    C2pion[chtype]->SetLineColor(4);
    C2muonSmeared[chtype]->SetLineColor(2);
    C2pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
    C2pion[chtype]->SetMinimum(0.98); C2pion[chtype]->SetMaximum(1.35);
    C2pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
    C2pion[chtype]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
    C2pion[chtype]->GetYaxis()->SetTitle("C_{2}K_{2}");
    //
    C3pion[chtype]->SetLineColor(4);
    C3muonSmeared[chtype]->SetLineColor(2);
    C3pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
    C3pion[chtype]->SetMinimum(0.90); C3pion[chtype]->SetMaximum(2.0);
    C3pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
    C3pion[chtype]->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
    C3pion[chtype]->GetYaxis()->SetTitle("C_{3}K_{3}");
  }
  //
  //
  
  C2pion[ChComb]->Draw();
  C2muonSmeared[ChComb]->Draw("same");
  legend->AddEntry(C2pion[ChComb],"Input Pion-Pion, C_{2}^{#pi-#pi,QS}K_{2}^{#pi-#pi,QS}","l");
  legend->AddEntry(C2muonSmeared[ChComb],"Pion-Muon residual, C_{2}^{#mu-#pi,QS}K_{2}^{#mu-#pi,QS}","l");
  legend->Draw("same");
  //PionPionK2->Draw("same");
  //
  //
  
  //
  //C3pion[ChComb]->Draw();
  //C3muonSmeared[ChComb]->Draw("same");
  //PionPionK3[0]->Draw("same");
  //legend->AddEntry(C3pion[ChComb],"Input Pion-Pion-Pion, C_{3}^{#pi-#pi-#pi,QS}K_{3}^{#pi-#pi-#pi,QS}","l");
  //legend->AddEntry(C3muonSmeared[ChComb],"Muon-Pion-Pion residual, C_{3}^{#mu-#pi-#pi,QS}K_{3}^{#mu-#pi-#pi,QS}","l");
  //legend->Draw("same");
    
  
  // corrections
  TFile *fout=new TFile("MuonCorrection_temp.root","RECREATE");
  TH1D *C2muonCorrection[2]; 
  C2muonCorrection[0] = new TH1D("C2muonCorrection_SC","",100,0,0.5);
  C2muonCorrection[1] = new TH1D("C2muonCorrection_MC","",100,0,0.5);
  TH1D *WeightmuonCorrection = new TH1D("WeightmuonCorrection","",100,0,0.5);
  TH1D *C3muonCorrection[2]; 
  C3muonCorrection[0] = new TH1D("C3muonCorrection_SC","",50,0,0.5);
  C3muonCorrection[1] = new TH1D("C3muonCorrection_MC","",50,0,0.5);
  //
  C2muonCorrection[0]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[0]->GetYaxis()->SetTitle("x_{2}"); 
  C2muonCorrection[1]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[1]->GetYaxis()->SetTitle("x_{2}"); 
  C3muonCorrection[0]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[0]->GetYaxis()->SetTitle("x_{3}"); 
  C3muonCorrection[1]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[1]->GetYaxis()->SetTitle("x_{3}"); 
  WeightmuonCorrection->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); WeightmuonCorrection->GetYaxis()->SetTitle("x_{2}^{w}"); 
  
  // 0.944 and 0.959
  float GoodPairFraction=1-0.93*(1-PurityNum->GetBinContent(10));
  cout<<"Pion Pair Purity = "<<PurityNum->GetBinContent(10)<<endl;
  cout<<"Effective Pion Pair Purity = "<<GoodPairFraction<<endl;
  float pionPurity=pow(GoodPairFraction,0.5);
  float muonPurity=1-pionPurity;
  for(int chtype=0; chtype<2; chtype++){
    for(int bin=1; bin<=100; bin++){
      
      bool emptybin2=kFALSE, emptybin3=kFALSE;
      if(PionPionK2[chtype]->GetBinContent(bin)==0) {PionPionK2[chtype]->SetBinContent(bin, 1.00001);}
      if(PionPionK3[chtype]->GetBinContent(bin)==0) {PionPionK3[chtype]->SetBinContent(bin, 1.00001);}
      if(bin > C2pion[chtype]->GetNbinsX()) emptybin2=kTRUE;
      if(bin > C3pion[chtype]->GetNbinsX()) emptybin3=kTRUE;
 
      double value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin);
      double den = (GoodPairFraction*C2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)*C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin);
      if(den > 0 && !emptybin2) C2muonCorrection[chtype]->SetBinContent(bin,value/den);
      else C2muonCorrection[chtype]->SetBinContent(bin, 1);
      //
      if(chtype==0){
	value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin) - 1.0;
	den = ((GoodPairFraction*C2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)*C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin)) - 1.0;
	if(den > 0 && !emptybin2) WeightmuonCorrection->SetBinContent(bin,value/den);
      }
      //
      value = C3pion[chtype]->GetBinContent(bin)/PionPionK3[chtype]->GetBinContent(bin);
      den = (pow(pionPurity,3)*C3pion[chtype]->GetBinContent(bin) + (3*pow(pionPurity,2)*muonPurity)*C3muonSmeared[chtype]->GetBinContent(bin))/PionPionK3[chtype]->GetBinContent(bin);
      if(den > 0 && !emptybin3) C3muonCorrection[chtype]->SetBinContent(bin,value/den);
      else C3muonCorrection[chtype]->SetBinContent(bin, 1);
    }
    //C2muonCorrection[chtype]->SetBinContent(1, C2muonCorrection[chtype]->GetBinContent(2));
    //C3muonCorrection[chtype]->SetBinContent(1, C3muonCorrection[chtype]->GetBinContent(2));
    C2muonCorrection[chtype]->Write();
    C3muonCorrection[chtype]->Write();
    if(chtype==0) {
      //WeightmuonCorrection->SetBinContent(1, WeightmuonCorrection->GetBinContent(2));
      WeightmuonCorrection->Write();
    }
  }
  
  //
  //C3muonCorrection[0]->SetMinimum(0.99); 
  //C3muonCorrection[0]->SetMaximum(1.05); 
  //C3muonCorrection[0]->GetYaxis()->SetTitleOffset(1.3);
  //C3muonCorrection[0]->Draw();
  //WeightmuonCorrection->GetYaxis()->SetTitleOffset(1.3);
  //WeightmuonCorrection->Draw();
  fout->Close();
  
  
}
Commit	Line	Data
3e748846	1	#include <math.h>
	2	#include <time.h>
	3	#include <stdio.h>
	4	#include <stdlib.h>
	5	#include <Riostream.h>
	6	#include <assert.h>
	7
	8	#include "TVector2.h"
	9	#include "TFile.h"
	10	#include "TString.h"
	11	#include "TF1.h"
	12	#include "TF3.h"
	13	#include "TH1.h"
	14	#include "TH2.h"
	15	#include "TH3.h"
	16	#include "TProfile.h"
	17	#include "TProfile2D.h"
	18	#include "TMath.h"
	19	#include "TText.h"
	20	#include "TRandom3.h"
	21	#include "TArray.h"
	22	#include "TLegend.h"
	23	#include "TStyle.h"
	24	#include "TMinuit.h"
	25	#include "TCanvas.h"
	26	#include "TLatex.h"
	27	#include "TPaveStats.h"
	28	#include "TASImage.h"
	29
	30	#define BohrR 1963.6885 // Bohr Radius for pions
	31	#define FmToGeV 0.19733 // conversion to fm
	32	#define PI 3.1415926
	33	#define masspiC 0.1395702 // pi+ mass (GeV/c^2)
	34	#define kappa3 0.16
	35	#define kappa4 0.40
	36
	37	using namespace std;
	38
	39
	40	void Plot_muons(){
	41
	42
	43	gStyle->SetOptStat(0);
	44	gStyle->SetOptDate(0);
	45	gStyle->SetOptFit(0111);
	46
	47	int ChComb=1;// 0=Same-charge, 1=Mixed-charge
	48	int kTbin_L=5, kTbin_H=6;
	49	int binKT3=1;//1-2
	50
	51	//TFile *_file0= new TFile("Results/PDC_12a17a_muons_R6.root","READ");
	52	//TFile *_file0= new TFile("Results/PDC_13b2_efix_p1_muons_R2.root","READ");
	53	TFile *_file0= new TFile("Results/PDC_12a17e_muons_R4.root","READ");
	54
	55	TList MyList=(TList)_file0->Get("MyList");
	56	_file0->Close();
	57
	58	TH1D PionCandidates=(TH1D)MyList->FindObject("fPionCandidates");
	59	PionCandidates->GetXaxis()->SetTitle("PDG code");
	60	//PionCandidates->Draw();
	61	//
	62	TH1D MuonParentsPrimary=(TH1D)MyList->FindObject("fMuonParents");
	63	MuonParentsPrimary->GetXaxis()->SetTitle("PDG code");
	64	MuonParentsPrimary->SetFillColor(1);
65	//MuonParentsPrimary->Draw();
66	//
67	TH1D MuonParentsSecondary=(TH1D)MyList->FindObject("fSecondaryMuonParents");
68	MuonParentsSecondary->GetXaxis()->SetTitle("PDG code");
69	MuonParentsSecondary->SetFillColor(1);
70	//MuonParentsSecondary->Draw();
71	//
72	// M0 R10-R6, M6 for R4, M17 for R2
73	TH3D PurityNum_3D = (TH3D)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityNum");
74	TH2D PurityDen_2D = (TH2D)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityDen");
75	TH1D *PurityNum=PurityNum_3D->ProjectionX("PurityNum",kTbin_L,kTbin_H,1,20);
76	double PurityNorm=PurityDen_2D->Integral(kTbin_L,kTbin_H,1,20);
77	PurityNum->Scale(1/PurityNorm);
78	char *namesAxis[15]={"e-e","e-mu","e-pi","e-k","e-p","mu-mu","mu-pi","mu-k","mu-p","pi-pi","pi-k","pi-p","k-k","k-p","p-p"};
79	for(int i=1; i<=15; i++) PurityNum->GetXaxis()->SetBinLabel(i, namesAxis[i-1]);
80	PurityNum->GetXaxis()->SetRange(1,15);
81	PurityNum->GetYaxis()->SetTitle("Probability");
82	PurityNum->Draw();
83	//
84	//
85	TCanvas *can = new TCanvas("can", "can",800,0,800,800);// 11,53,700,500
86	can->SetHighLightColor(2);
87	gStyle->SetOptFit(0111);
88	can->SetFillColor(10);//10
89	can->SetBorderMode(0);
90	can->SetBorderSize(2);
91	can->SetFrameFillColor(0);
92	can->SetFrameBorderMode(0);
93	can->SetFrameBorderMode(0);
94	can->cd();
95	TPad *pad = new TPad("pad","pad",0.,0.,1.,1.);
96	gPad->SetTickx();
97	gPad->SetTicky();
98	pad->SetGridx();
99	pad->SetGridy();
100	pad->SetTopMargin(0.02);//0.05
101	pad->SetRightMargin(0.02);//3e-2
102	pad->SetBottomMargin(0.1);//0.12
103	pad->SetLeftMargin(0.1);
104	pad->Draw();
105	pad->cd();
106	TLegend *legend = new TLegend(.5,.65, .9,.95,NULL,"brNDC");//.45 or .4 for x1
107	legend->SetBorderSize(0);
108	legend->SetFillColor(0);
109	legend->SetTextFont(42);
110	legend->SetTextSize(0.03);
111	//
112	TH3D MuonSmearedNum2_3=(TH3D)MyList->FindObject("fMuonContamSmearedNum2");
113	TH3D MuonSmearedDen2_3=(TH3D)MyList->FindObject("fMuonContamSmearedDen2");
114	TH3D PionNum2_3=(TH3D)MyList->FindObject("fMuonContamIdealNum2");
115	TH3D PionDen2_3=(TH3D)MyList->FindObject("fMuonContamIdealDen2");
116	TH3D PionPionK2_3=(TH3D)MyList->FindObject("fPionPionK2");
117	//
118	TH3D MuonSmearedNum3_3=(TH3D)MyList->FindObject("fMuonContamSmearedNum3");
119	TH3D MuonSmearedDen3_3=(TH3D)MyList->FindObject("fMuonContamSmearedDen3");
120	TH3D PionNum3_3=(TH3D)MyList->FindObject("fMuonContamIdealNum3");
121	TH3D PionDen3_3=(TH3D)MyList->FindObject("fMuonContamIdealDen3");
122	TH3D PionPionK3_3=(TH3D)MyList->FindObject("fPionPionK3");
123	TH3D MuonPionK3_3=(TH3D)MyList->FindObject("fMuonPionK3");
124
125
126	TH1D *MuonSmearedNum2[2];// SC/MC
127	TH1D *MuonSmearedDen2[2];
128	TH1D *PionNum2[2];
129	TH1D *PionDen2[2];
130	TH1D *PionPionK2[2];
131	//
132	TH1D *MuonSmearedNum3[2];// SC/MC
133	TH1D *MuonSmearedDen3[2];
134	TH1D *PionNum3[2];
135	TH1D *PionDen3[2];
136	TH1D *PionPionK3[2];
137	//
138	TH1D *C2muonSmeared[2];
139	TH1D *C2pion[2];
140	TH1D *C3muonSmeared[2];
141	TH1D *C3pion[2];
142	//
143	for(int chtype=0; chtype<2; chtype++){
144	TString *names[10];
145	for(int i=0; i<10; i++) {names[i]=new TString("name_"); names[i] += i; names[i] += chtype;}
146
147	MuonSmearedNum2[chtype]=(TH1D*)MuonSmearedNum2_3->ProjectionZ(names[0]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
148	MuonSmearedDen2[chtype]=(TH1D*)MuonSmearedDen2_3->ProjectionZ(names[1]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
149	PionNum2[chtype]=(TH1D*)PionNum2_3->ProjectionZ(names[2]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
150	PionDen2[chtype]=(TH1D*)PionDen2_3->ProjectionZ(names[3]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
151	PionPionK2[chtype]=(TH1D*)PionPionK2_3->ProjectionZ(names[4]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
152	PionPionK2[chtype]->Divide(PionDen2[chtype]);
153	////////////////
154	MuonSmearedNum3[chtype]=(TH1D*)MuonSmearedNum3_3->ProjectionZ(names[5]->Data(),chtype+1,chtype+1,binKT3,binKT3);
155	MuonSmearedDen3[chtype]=(TH1D*)MuonSmearedDen3_3->ProjectionZ(names[6]->Data(),chtype+1,chtype+1,binKT3,binKT3);
156	PionNum3[chtype]=(TH1D*)PionNum3_3->ProjectionZ(names[7]->Data(),chtype+1,chtype+1,binKT3,binKT3);
157	PionDen3[chtype]=(TH1D*)PionDen3_3->ProjectionZ(names[8]->Data(),chtype+1,chtype+1,binKT3,binKT3);
158	PionPionK3[chtype]=(TH1D*)PionPionK3_3->ProjectionZ(names[9]->Data(),chtype+1,chtype+1,binKT3,binKT3);
159	PionPionK3[chtype]->Divide(PionDen3[chtype]);
160	//
161	C2muonSmeared[chtype]=(TH1D*)MuonSmearedNum2[chtype]->Clone();
162	C2pion[chtype]=(TH1D*)PionNum2[chtype]->Clone();
163	C2muonSmeared[chtype]->Divide(MuonSmearedDen2[chtype]);
164	C2pion[chtype]->Divide(PionDen2[chtype]);
165	//
166	C3muonSmeared[chtype]=(TH1D*)MuonSmearedNum3[chtype]->Clone();
167	C3pion[chtype]=(TH1D*)PionNum3[chtype]->Clone();
168	C3muonSmeared[chtype]->Divide(MuonSmearedDen3[chtype]);
169	C3pion[chtype]->Divide(PionDen3[chtype]);
170	//
171	//
172	C2pion[chtype]->SetLineColor(4);
173	C2muonSmeared[chtype]->SetLineColor(2);
174	C2pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
175	C2pion[chtype]->SetMinimum(0.98); C2pion[chtype]->SetMaximum(1.35);
176	C2pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
177	C2pion[chtype]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
178	C2pion[chtype]->GetYaxis()->SetTitle("C_{2}K_{2}");
179	//
180	C3pion[chtype]->SetLineColor(4);
181	C3muonSmeared[chtype]->SetLineColor(2);
182	C3pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
183	C3pion[chtype]->SetMinimum(0.90); C3pion[chtype]->SetMaximum(2.0);
184	C3pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
185	C3pion[chtype]->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
186	C3pion[chtype]->GetYaxis()->SetTitle("C_{3}K_{3}");
187	}
188	//
189	//
190
191	C2pion[ChComb]->Draw();
192	C2muonSmeared[ChComb]->Draw("same");
193	legend->AddEntry(C2pion[ChComb],"Input Pion-Pion, C_{2}^{#pi-#pi,QS}K_{2}^{#pi-#pi,QS}","l");
194	legend->AddEntry(C2muonSmeared[ChComb],"Pion-Muon residual, C_{2}^{#mu-#pi,QS}K_{2}^{#mu-#pi,QS}","l");
195	legend->Draw("same");
196	//PionPionK2->Draw("same");
197	//
198	//
199
200	//
201	//C3pion[ChComb]->Draw();
202	//C3muonSmeared[ChComb]->Draw("same");
203	//PionPionK3[0]->Draw("same");
204	//legend->AddEntry(C3pion[ChComb],"Input Pion-Pion-Pion, C_{3}^{#pi-#pi-#pi,QS}K_{3}^{#pi-#pi-#pi,QS}","l");
205	//legend->AddEntry(C3muonSmeared[ChComb],"Muon-Pion-Pion residual, C_{3}^{#mu-#pi-#pi,QS}K_{3}^{#mu-#pi-#pi,QS}","l");
206	//legend->Draw("same");
207
208
209	// corrections
210	TFile *fout=new TFile("MuonCorrection_temp.root","RECREATE");
211	TH1D *C2muonCorrection[2];
212	C2muonCorrection[0] = new TH1D("C2muonCorrection_SC","",100,0,0.5);
213	C2muonCorrection[1] = new TH1D("C2muonCorrection_MC","",100,0,0.5);
214	TH1D *WeightmuonCorrection = new TH1D("WeightmuonCorrection","",100,0,0.5);
215	TH1D *C3muonCorrection[2];
216	C3muonCorrection[0] = new TH1D("C3muonCorrection_SC","",50,0,0.5);
217	C3muonCorrection[1] = new TH1D("C3muonCorrection_MC","",50,0,0.5);
218	//
219	C2muonCorrection[0]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[0]->GetYaxis()->SetTitle("x_{2}");
220	C2muonCorrection[1]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[1]->GetYaxis()->SetTitle("x_{2}");
221	C3muonCorrection[0]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[0]->GetYaxis()->SetTitle("x_{3}");
222	C3muonCorrection[1]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[1]->GetYaxis()->SetTitle("x_{3}");
223	WeightmuonCorrection->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); WeightmuonCorrection->GetYaxis()->SetTitle("x_{2}^{w}");
224
225	// 0.944 and 0.959
226	float GoodPairFraction=1-0.93*(1-PurityNum->GetBinContent(10));
227	cout<<"Pion Pair Purity = "<<PurityNum->GetBinContent(10)<<endl;
228	cout<<"Effective Pion Pair Purity = "<<GoodPairFraction<<endl;
229	float pionPurity=pow(GoodPairFraction,0.5);
230	float muonPurity=1-pionPurity;
231	for(int chtype=0; chtype<2; chtype++){
232	for(int bin=1; bin<=100; bin++){
233
234	bool emptybin2=kFALSE, emptybin3=kFALSE;
235	if(PionPionK2[chtype]->GetBinContent(bin)==0) {PionPionK2[chtype]->SetBinContent(bin, 1.00001);}
236	if(PionPionK3[chtype]->GetBinContent(bin)==0) {PionPionK3[chtype]->SetBinContent(bin, 1.00001);}
237	if(bin > C2pion[chtype]->GetNbinsX()) emptybin2=kTRUE;
238	if(bin > C3pion[chtype]->GetNbinsX()) emptybin3=kTRUE;
239
240	double value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin);
241	double den = (GoodPairFractionC2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin);
242	if(den > 0 && !emptybin2) C2muonCorrection[chtype]->SetBinContent(bin,value/den);
243	else C2muonCorrection[chtype]->SetBinContent(bin, 1);
244	//
245	if(chtype==0){
246	value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin) - 1.0;
247	den = ((GoodPairFractionC2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin)) - 1.0;
248	if(den > 0 && !emptybin2) WeightmuonCorrection->SetBinContent(bin,value/den);
249	}
250	//
251	value = C3pion[chtype]->GetBinContent(bin)/PionPionK3[chtype]->GetBinContent(bin);
252	den = (pow(pionPurity,3)C3pion[chtype]->GetBinContent(bin) + (3pow(pionPurity,2)muonPurity)C3muonSmeared[chtype]->GetBinContent(bin))/PionPionK3[chtype]->GetBinContent(bin);
253	if(den > 0 && !emptybin3) C3muonCorrection[chtype]->SetBinContent(bin,value/den);
254	else C3muonCorrection[chtype]->SetBinContent(bin, 1);
255	}
256	//C2muonCorrection[chtype]->SetBinContent(1, C2muonCorrection[chtype]->GetBinContent(2));
257	//C3muonCorrection[chtype]->SetBinContent(1, C3muonCorrection[chtype]->GetBinContent(2));
258	C2muonCorrection[chtype]->Write();
259	C3muonCorrection[chtype]->Write();
260	if(chtype==0) {
261	//WeightmuonCorrection->SetBinContent(1, WeightmuonCorrection->GetBinContent(2));
262	WeightmuonCorrection->Write();
263	}
264	}
265
266	//
267	//C3muonCorrection[0]->SetMinimum(0.99);
268	//C3muonCorrection[0]->SetMaximum(1.05);
269	//C3muonCorrection[0]->GetYaxis()->SetTitleOffset(1.3);
270	//C3muonCorrection[0]->Draw();
271	//WeightmuonCorrection->GetYaxis()->SetTitleOffset(1.3);
272	//WeightmuonCorrection->Draw();
273	fout->Close();
274
275
276	}