[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / corrusingfitsigma.C

//This macro calculates the correction factor for the difference of the n sigma TOF distributions in data and mc  


class AliSpectraBothHistoManager;
AliSpectraBothHistoManager* managerdata=0x0;
AliSpectraBothHistoManager* managermc=0x0;

Double_t Gaus_plus_tail(Double_t* x, Double_t* par )
{

	Double_t mean=par[0];
	Double_t rms=par[1];
	Double_t c=par[2];
	Double_t slope=par[3]/par[1];//df/dx continues
	Double_t cut=par[3]; 


	Double_t one_over_sqrt_2pi=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
	Double_t returnvalue=0.0;

	Double_t n=0.5*(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/(TMath::Abs(rms)*slope);
	if (x[0]<mean+cut*rms)
		returnvalue=TMath::Exp(-1.0*(x[0]-mean)*(x[0]-mean)/(2.0*rms*rms))*one_over_sqrt_2pi/TMath::Abs(rms);
	else
		returnvalue=TMath::Exp(slope*(mean+cut*rms-x[0]))*TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/TMath::Abs(rms);
                                                                                                         
	return c*returnvalue/n;	
}

//Double_t background(Double_t x, Double_t scale, Double_t slope, Double_t startpoint)
Double_t background(Double_t* x,Double_t* par)
{
	//at start ponit it is eqaul to scale
	Double_t scale=par[0];
	Double_t slope=par[1];
	Double_t startpoint=par[2];
	return scale*TMath::Exp((startpoint-x[0])*slope);
}


Double_t Gaus_on_background(Double_t* x,Double_t* par)
{
	//par[0] peak position 
	//par[1] sigma gaus
	//par[2] Normalization of signal
	//par[3] cut
	//par[4] scale of bg
	//par[5] slope of bg
	//par[6] start point of bg const 
	//cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<par[6]<<par[7]<<" "<<x[0]<< endl;
	Double_t returnvalue=0.0;
	 returnvalue+=Gaus_plus_tail(x,par);
	//returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);
	returnvalue+=background(x,&(par[4]));
	return returnvalue;
}

Int_t   OpenFile(TString nameFile,TString outputname,Int_t mg1)
{
	

	TFile *file = TFile::Open(nameFile.Data());
	if(!file)
	{
		cout<<"no file"<<endl;
		return -1;
	}	
	TDirectoryFile *dir=(TDirectoryFile*)file->Get(outputname.Data());
	if(!dir)
	{
		cout<<"no dir "<<outputname.Data()<<endl;
		return -1;
	}
	if(!dir->Get("SpectraHistos"))
		return -1;
	if(mg1==1)
	{	
		managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
		return managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();

	}
	else 
	{
		managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
		return managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();
	}
	return -1;
}


void corrusingfitsigma(TString nameFile1,TString outputname1,TString nameFile2,TString outputname2, Float_t nsigmacut=3.0, Float_t minpt=0.6,Float_t maxpt=1.4, Bool_t saveoption=kFALSE)
{	
	ofstream oufiletxt("TOFcorrPID.txt");
	Float_t nsigmacut=TMath::Abs(nsigmacut);
	TString pidmethods[3]={"TPC","TOF","TPCTOF"};
	TString Particle[]={"Pion","Kaon","Proton"};
	gStyle->SetOptStat(0);	
	TH1::AddDirectory(kFALSE);
	gSystem->Load("libCore");
	gSystem->Load("libPhysics");
	gSystem->Load("libTree");
	gSystem->Load("libMatrix");
	gSystem->Load("libSTEERBase");
	gSystem->Load("libESD");
	gSystem->Load("libAOD");
	gSystem->Load("libANALYSIS");
	gSystem->Load("libOADB");
	gSystem->Load("libANALYSISalice");
	gSystem->Load("libTender");
	gSystem->Load("libCORRFW");
	gSystem->Load("libPWGTools");
	gSystem->Load("libPWGLFspectra");
	Int_t ne1= OpenFile(nameFile1,outputname1,1);
	Int_t ne2= OpenFile(nameFile2,outputname2,2);
	Int_t imethod=1;
		oufiletxt<<"file 1 "<<nameFile1.Data()<<endl;
	oufiletxt<<"outputname 1 "<<outputname1.Data()<<endl;
	oufiletxt<<"file 2 "<<nameFile2.Data()<<endl;
	oufiletxt<<"outputname 2 "<<outputname2.Data()<<endl;
	oufiletxt<<"nsigma "<<nsigmacut<<endl;
	oufiletxt<<"minpt "<<minpt<<endl;
	oufiletxt<<"maxpt "<<maxpt<<endl;
	
	TList* lout=new TList();

	for(Int_t ipart=0;ipart<3;ipart++)
	{
		if(!managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
			continue;		
		TH2F *nsig_data = (TH2F*)((TH2F*)managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
		if(!nsig_data)
			continue;

		if(!managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
			continue;
		TH2F *nsig_mc = (TH2F*)((TH2F*)managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
		if(!nsig_mc)
			continue;
		Int_t firstbin=nsig_data->GetXaxis()->FindBin(minpt);
		Int_t lastbin=nsig_data->GetXaxis()->FindBin(maxpt);
		Float_t ptlow=nsig_data->GetXaxis()->GetBinLowEdge(firstbin);
		Float_t ptup=nsig_data->GetXaxis()->GetBinUpEdge(lastbin);

		TH1D* hist1data= new TH1D(Form("datahist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
		TH1D* hist2data= new TH1D(Form("datahist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
		TH1D* sigmadata= new TH1D(Form("sigmadata%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	
		hist1data->SetLineColor(kBlack);
		hist2data->SetLineColor(kBlack);
		sigmadata->SetMarkerColor(kBlack);
		sigmadata->SetMarkerStyle(22);
		hist1data->SetLineStyle(1);
		hist2data->SetLineStyle(2);
		hist1data->SetLineWidth(2);
		hist2data->SetLineWidth(2);

	
		TH1D* hist1mc= new TH1D(Form("mchist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
		TH1D* hist2mc= new TH1D(Form("mchist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
		TH1D* sigmamc= new TH1D(Form("sigmamc%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);

		hist1mc->SetLineColor(kRed);
		hist2mc->SetLineColor(kRed);
		sigmamc->SetMarkerColor(kRed);
		sigmamc->SetMarkerStyle(23);
		sigmamc->SetLineColor(kRed);
		hist1mc->SetLineStyle(1);
		hist2mc->SetLineStyle(2);
		hist1mc->SetLineWidth(2);
		hist2mc->SetLineWidth(2);


		for(int ibin=firstbin;ibin<lastbin;ibin++)
		{
				
			TH1F *nsig_data_Proj1=(TH1F*)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
			nsig_data_Proj1->Sumw2();
			nsig_data_Proj1->GetXaxis()->SetRangeUser(-6,6);
			nsig_data_Proj1->SetLineColor(kBlack);		
			TH1F *nsig_mc_Proj1=(TH1F*)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
                	nsig_mc_Proj1->Sumw2();
                	nsig_mc_Proj1->GetXaxis()->SetRangeUser(-6,6);
			nsig_mc_Proj1->SetLineColor(kRed);

			TF1* gfundata=new TF1("fun_for_fitdata",Gaus_on_background,-5,5,7);
			gfundata->SetParameter(0,0.0);
			gfundata->SetParameter(1,1.0);
			gfundata->SetParameter(2,nsig_data_Proj1->GetBinContent(nsig_data_Proj1->GetXaxis()->FindBin(0.0)));
			gfundata->FixParameter(3,1.0);
			gfundata->FixParameter(4,0.0);
			gfundata->FixParameter(5,1.0);
			gfundata->FixParameter(6,-5);
			gfundata->SetLineColor(kBlack);
			nsig_data_Proj1->Fit(gfundata,"SR+","E");
			Float_t sigmadatav=gfundata->GetParameter(1);
			Float_t meandatav=gfundata->GetParameter(0);
			sigmadata->SetBinContent(ibin-firstbin+1,gfundata->GetParameter(1));	
			sigmadata->SetBinError(ibin-firstbin+1,gfundata->GetParError(1));


			TF1* gfunmc=new TF1("fun_for_fitmc",Gaus_on_background,-5,5,7);
			gfunmc->SetParameter(0,0.0);
			gfunmc->SetParameter(1,1.0);
			gfunmc->SetParameter(2,nsig_mc_Proj1->GetBinContent(nsig_mc_Proj1->GetXaxis()->FindBin(0.0)));
			gfunmc->FixParameter(3,1.0);
			gfunmc->FixParameter(4,0.0);
			gfunmc->FixParameter(5,1.0);
			gfunmc->FixParameter(6,-5);
			gfunmc->SetLineColor(kRed);
			nsig_mc_Proj1->Fit(gfunmc,"SR+","E");
			sigmamc->SetBinContent(ibin-firstbin+1,gfunmc->GetParameter(1));	
			sigmamc->SetBinError(ibin-firstbin+1,gfunmc->GetParError(1));
			Float_t sigmamcv=gfunmc->GetParameter(1);
			Float_t meanmcv=gfundata->GetParameter(0);


			hist1data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));			
			Float_t width=nsig_data_Proj1->GetXaxis()->GetBinWidth(-1.0*nsigmacut);
			Int_t lowbind=nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav);
			Int_t upbind=nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav);
			Float_t lowbinpartd=TMath::Abs(nsig_data_Proj1->GetXaxis()->GetBinLowEdge(lowbind)+1.0*nsigmacut*sigmadatav-meandatav);
			Float_t upbinpartd=TMath::Abs(nsigmacut*sigmadatav+meandatav-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbind));
			

			Float_t corrd=(lowbinpartd*nsig_data_Proj1->GetBinContent(lowbind)+upbinpartd*nsig_data_Proj1->GetBinContent(upbind))/width;

			hist2data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))-corrd);
			cout<<corrd/nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))<<endl;

			hist1mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));
			Int_t lowbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv);
			Int_t upbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv);
			Float_t lowbinpartmc=TMath::Abs(nsig_mc_Proj1->GetXaxis()->GetBinLowEdge(lowbinmc)+1.0*nsigmacut*sigmamcv-meanmcv);
			Float_t upbinpartmc=TMath::Abs(nsigmacut*sigmamcv+meanmcv-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbinmc));
			Float_t corrmc=(lowbinpartmc*nsig_mc_Proj1->GetBinContent(lowbinmc)+upbinpartmc*nsig_mc_Proj1->GetBinContent(upbinmc))/width;
			cout<<corrmc/nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))<<" aaa "<<endl;

			hist2mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))-corrmc);
			TCanvas* ctmp=new TCanvas(Form("c%d%s",ibin,Particle[ipart].Data()),Form("c%d%s",ibin,Particle[ipart].Data()),800,800);
			ctmp->cd()->SetLogy();

			nsig_data_Proj1->GetXaxis()->SetTitle("#sigma_{TOF}");
			nsig_data_Proj1->GetYaxis()->SetTitle("N_{entries}");
			nsig_data_Proj1->SetTitle(Form("PionTOF[%.2f,%.2f]",nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)));
			nsig_data_Proj1->Draw("E1");
			nsig_mc_Proj1->Draw("E1same");
			if(saveoption)
				ctmp->SaveAs(Form("c%d%s.png",ibin,Particle[ipart].Data()));
			delete ctmp;
			delete gfundata;
			delete gfunmc;
			delete nsig_data_Proj1;
			delete nsig_mc_Proj1;
			
		}
	
			hist1data->Sumw2();
			hist1data->Print("all");
			hist1data->Sumw2();
			hist1mc->Sumw2();
			hist2mc->Sumw2();

		hist2data->GetXaxis()->SetTitle("p_{T} (GeV/c)");
		hist2data->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
		hist2mc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
		hist2mc->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
		hist2data->Divide(hist2data,hist1data,1,1,"B");
		hist2mc->Divide(hist2mc,hist1mc,1,1,"B");
		lout->Add(hist2data);	
		lout->Add(hist2mc);	

		sigmadata->GetXaxis()->SetTitle("p_{T} (GeV/c)");
		sigmadata->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
		sigmamc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
		sigmamc->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
		

		lout->Add(sigmadata);
		lout->Add(sigmamc);

		TH1D* corr=(TH1D*)hist2data->Clone(Form("corr%s",Particle[ipart].Data()));	
		corr->Divide(hist2mc);
		corr->SetLineColor(kBlue);
		corr->GetYaxis()->SetTitle("correction");
		corr->Fit("pol0","R0","",minpt,maxpt);
		
		oufiletxt<<Particle[ipart].Data()<<" "<<((TF1*)corr->GetListOfFunctions()->At(0))->GetParameter(0)<<endl;
		lout->Add(corr);
		TCanvas* corrcanvas=new TCanvas(Form("corrcanvas%s",Particle[ipart].Data()), Form("corrcanvas%s",Particle[ipart].Data()),600,600);
		corrcanvas->SetMargin(0.12,0.05,0.15,0.1);
		corr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
		corr->GetYaxis()->SetTitle("TOF PID mismatch corr.");
		corr->GetYaxis()->SetTitleOffset(1.5*corr->GetYaxis()->GetTitleOffset());
		corr->GetYaxis()->SetRangeUser(0.98,1.04);
	
		corrcanvas->cd();
		corr->Draw("E1");
		corrcanvas->SaveAs(Form("Corr%s.eps",Particle[ipart].Data()));
		corrcanvas->SaveAs(Form("Corr%s.png",Particle[ipart].Data()));

	}
	TFile* fileout=TFile::Open("TOFcorrPID.root","Recreate");
	
	lout->Write("output",TObject::kSingleKey);
	fileout->Close();
}
Commit	Line	Data
f887e5c3	1	//This macro calculates the correction factor for the difference of the n sigma TOF distributions in data and mc
	2
	3
	4	class AliSpectraBothHistoManager;
	5	AliSpectraBothHistoManager* managerdata=0x0;
	6	AliSpectraBothHistoManager* managermc=0x0;
	7
	8	Double_t Gaus_plus_tail(Double_t* x, Double_t* par )
	9	{
	10
	11	Double_t mean=par[0];
	12	Double_t rms=par[1];
	13	Double_t c=par[2];
	14	Double_t slope=par[3]/par[1];//df/dx continues
	15	Double_t cut=par[3];
	16
	17
	18
	19	Double_t one_over_sqrt_2pi=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
	20	Double_t returnvalue=0.0;
	21
	22	Double_t n=0.5(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cutcut0.5)one_over_sqrt_2pi/(TMath::Abs(rms)*slope);
	23	if (x[0]<mean+cut*rms)
	24	returnvalue=TMath::Exp(-1.0(x[0]-mean)(x[0]-mean)/(2.0rmsrms))*one_over_sqrt_2pi/TMath::Abs(rms);
	25	else
	26	returnvalue=TMath::Exp(slope(mean+cutrms-x[0]))TMath::Exp(-cutcut0.5)one_over_sqrt_2pi/TMath::Abs(rms);
	27
	28	return c*returnvalue/n;
	29	}
	30
	31	//Double_t background(Double_t x, Double_t scale, Double_t slope, Double_t startpoint)
	32	Double_t background(Double_t* x,Double_t* par)
	33	{
	34	//at start ponit it is eqaul to scale
	35	Double_t scale=par[0];
	36	Double_t slope=par[1];
	37	Double_t startpoint=par[2];
	38	return scaleTMath::Exp((startpoint-x[0])slope);
	39	}
	40
	41
	42	Double_t Gaus_on_background(Double_t* x,Double_t* par)
	43	{
	44	//par[0] peak position
	45	//par[1] sigma gaus
	46	//par[2] Normalization of signal
	47	//par[3] cut
	48	//par[4] scale of bg
	49	//par[5] slope of bg
	50	//par[6] start point of bg const
	51	//cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<par[6]<<par[7]<<" "<<x[0]<< endl;
	52	Double_t returnvalue=0.0;
	53	returnvalue+=Gaus_plus_tail(x,par);
	54	//returnvalue+=par[3]TMath::Exp((par[5]-x[0])par[4]);
	55	returnvalue+=background(x,&(par[4]));
	56	return returnvalue;
	57	}
	58
	59	Int_t OpenFile(TString nameFile,TString outputname,Int_t mg1)
	60	{
	61
	62
	63	TFile *file = TFile::Open(nameFile.Data());
	64	if(!file)
65	{
66	cout<<"no file"<<endl;
67	return -1;
68	}
69	TDirectoryFile dir=(TDirectoryFile)file->Get(outputname.Data());
70	if(!dir)
71	{
72	cout<<"no dir "<<outputname.Data()<<endl;
73	return -1;
74	}
75	if(!dir->Get("SpectraHistos"))
76	return -1;
77	if(mg1==1)
78	{
79	managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
80	return managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();
81
82	}
83	else
84	{
85	managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
86	return managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();
87	}
88	return -1;
89	}
90
91
92	void corrusingfitsigma(TString nameFile1,TString outputname1,TString nameFile2,TString outputname2, Float_t nsigmacut=3.0, Float_t minpt=0.6,Float_t maxpt=1.4, Bool_t saveoption=kFALSE)
93	{
94	ofstream oufiletxt("TOFcorrPID.txt");
95	Float_t nsigmacut=TMath::Abs(nsigmacut);
96	TString pidmethods[3]={"TPC","TOF","TPCTOF"};
97	TString Particle[]={"Pion","Kaon","Proton"};
98	gStyle->SetOptStat(0);
99	TH1::AddDirectory(kFALSE);
4070f709	100	gSystem->Load("libCore");
4070f709	101	gSystem->Load("libPhysics");
f887e5c3	102	gSystem->Load("libTree");
	103	gSystem->Load("libMatrix");
	104	gSystem->Load("libSTEERBase");
	105	gSystem->Load("libESD");
	106	gSystem->Load("libAOD");
	107	gSystem->Load("libANALYSIS");
	108	gSystem->Load("libOADB");
	109	gSystem->Load("libANALYSISalice");
af472fff	110	gSystem->Load("libTender");
f887e5c3	111	gSystem->Load("libCORRFW");
	112	gSystem->Load("libPWGTools");
	113	gSystem->Load("libPWGLFspectra");
	114	Int_t ne1= OpenFile(nameFile1,outputname1,1);
	115	Int_t ne2= OpenFile(nameFile2,outputname2,2);
	116	Int_t imethod=1;
	117	oufiletxt<<"file 1 "<<nameFile1.Data()<<endl;
	118	oufiletxt<<"outputname 1 "<<outputname1.Data()<<endl;
	119	oufiletxt<<"file 2 "<<nameFile2.Data()<<endl;
	120	oufiletxt<<"outputname 2 "<<outputname2.Data()<<endl;
	121	oufiletxt<<"nsigma "<<nsigmacut<<endl;
	122	oufiletxt<<"minpt "<<minpt<<endl;
	123	oufiletxt<<"maxpt "<<maxpt<<endl;
	124
	125	TList* lout=new TList();
	126
	127	for(Int_t ipart=0;ipart<3;ipart++)
	128	{
	129	if(!managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
	130	continue;
	131	TH2F nsig_data = (TH2F)((TH2F*)managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
	132	if(!nsig_data)
	133	continue;
	134
	135	if(!managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
	136	continue;
	137	TH2F nsig_mc = (TH2F)((TH2F*)managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
	138	if(!nsig_mc)
	139	continue;
	140	Int_t firstbin=nsig_data->GetXaxis()->FindBin(minpt);
	141	Int_t lastbin=nsig_data->GetXaxis()->FindBin(maxpt);
	142	Float_t ptlow=nsig_data->GetXaxis()->GetBinLowEdge(firstbin);
	143	Float_t ptup=nsig_data->GetXaxis()->GetBinUpEdge(lastbin);
	144
	145	TH1D* hist1data= new TH1D(Form("datahist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	146	TH1D* hist2data= new TH1D(Form("datahist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	147	TH1D* sigmadata= new TH1D(Form("sigmadata%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	148
	149	hist1data->SetLineColor(kBlack);
	150	hist2data->SetLineColor(kBlack);
	151	sigmadata->SetMarkerColor(kBlack);
	152	sigmadata->SetMarkerStyle(22);
	153	hist1data->SetLineStyle(1);
	154	hist2data->SetLineStyle(2);
	155	hist1data->SetLineWidth(2);
	156	hist2data->SetLineWidth(2);
	157
	158
	159	TH1D* hist1mc= new TH1D(Form("mchist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	160	TH1D* hist2mc= new TH1D(Form("mchist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	161	TH1D* sigmamc= new TH1D(Form("sigmamc%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
	162
	163	hist1mc->SetLineColor(kRed);
	164	hist2mc->SetLineColor(kRed);
	165	sigmamc->SetMarkerColor(kRed);
	166	sigmamc->SetMarkerStyle(23);
	167	sigmamc->SetLineColor(kRed);
	168	hist1mc->SetLineStyle(1);
	169	hist2mc->SetLineStyle(2);
	170	hist1mc->SetLineWidth(2);
	171	hist2mc->SetLineWidth(2);
	172
	173
	174
175	for(int ibin=firstbin;ibin<lastbin;ibin++)
176	{
177
178	TH1F nsig_data_Proj1=(TH1F)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
179	nsig_data_Proj1->Sumw2();
180	nsig_data_Proj1->GetXaxis()->SetRangeUser(-6,6);
181	nsig_data_Proj1->SetLineColor(kBlack);
182	TH1F nsig_mc_Proj1=(TH1F)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
183	nsig_mc_Proj1->Sumw2();
184	nsig_mc_Proj1->GetXaxis()->SetRangeUser(-6,6);
185	nsig_mc_Proj1->SetLineColor(kRed);
186
187	TF1* gfundata=new TF1("fun_for_fitdata",Gaus_on_background,-5,5,7);
188	gfundata->SetParameter(0,0.0);
189	gfundata->SetParameter(1,1.0);
190	gfundata->SetParameter(2,nsig_data_Proj1->GetBinContent(nsig_data_Proj1->GetXaxis()->FindBin(0.0)));
191	gfundata->FixParameter(3,1.0);
192	gfundata->FixParameter(4,0.0);
193	gfundata->FixParameter(5,1.0);
194	gfundata->FixParameter(6,-5);
195	gfundata->SetLineColor(kBlack);
196	nsig_data_Proj1->Fit(gfundata,"SR+","E");
197	Float_t sigmadatav=gfundata->GetParameter(1);
198	Float_t meandatav=gfundata->GetParameter(0);
199	sigmadata->SetBinContent(ibin-firstbin+1,gfundata->GetParameter(1));
200	sigmadata->SetBinError(ibin-firstbin+1,gfundata->GetParError(1));
201
202
203	TF1* gfunmc=new TF1("fun_for_fitmc",Gaus_on_background,-5,5,7);
204	gfunmc->SetParameter(0,0.0);
205	gfunmc->SetParameter(1,1.0);
206	gfunmc->SetParameter(2,nsig_mc_Proj1->GetBinContent(nsig_mc_Proj1->GetXaxis()->FindBin(0.0)));
207	gfunmc->FixParameter(3,1.0);
208	gfunmc->FixParameter(4,0.0);
209	gfunmc->FixParameter(5,1.0);
210	gfunmc->FixParameter(6,-5);
211	gfunmc->SetLineColor(kRed);
212	nsig_mc_Proj1->Fit(gfunmc,"SR+","E");
213	sigmamc->SetBinContent(ibin-firstbin+1,gfunmc->GetParameter(1));
214	sigmamc->SetBinError(ibin-firstbin+1,gfunmc->GetParError(1));
215	Float_t sigmamcv=gfunmc->GetParameter(1);
216	Float_t meanmcv=gfundata->GetParameter(0);
217
218
176a37e3	219	hist1data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));
	220	Float_t width=nsig_data_Proj1->GetXaxis()->GetBinWidth(-1.0*nsigmacut);
	221	Int_t lowbind=nsig_data_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmadatav+meandatav);
	222	Int_t upbind=nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav);
	223	Float_t lowbinpartd=TMath::Abs(nsig_data_Proj1->GetXaxis()->GetBinLowEdge(lowbind)+1.0nsigmacutsigmadatav-meandatav);
	224	Float_t upbinpartd=TMath::Abs(nsigmacut*sigmadatav+meandatav-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbind));
	225
	226
	227
	228	Float_t corrd=(lowbinpartdnsig_data_Proj1->GetBinContent(lowbind)+upbinpartdnsig_data_Proj1->GetBinContent(upbind))/width;
	229
	230	hist2data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))-corrd);
	231	cout<<corrd/nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))<<endl;
	232
f887e5c3	233	hist1mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));
176a37e3	234	Int_t lowbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmamcv+meanmcv);
	235	Int_t upbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv);
	236	Float_t lowbinpartmc=TMath::Abs(nsig_mc_Proj1->GetXaxis()->GetBinLowEdge(lowbinmc)+1.0nsigmacutsigmamcv-meanmcv);
	237	Float_t upbinpartmc=TMath::Abs(nsigmacut*sigmamcv+meanmcv-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbinmc));
	238	Float_t corrmc=(lowbinpartmcnsig_mc_Proj1->GetBinContent(lowbinmc)+upbinpartmcnsig_mc_Proj1->GetBinContent(upbinmc))/width;
	239	cout<<corrmc/nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))<<" aaa "<<endl;
	240
	241	hist2mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0nsigmacutsigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))-corrmc);
f887e5c3	242	TCanvas* ctmp=new TCanvas(Form("c%d%s",ibin,Particle[ipart].Data()),Form("c%d%s",ibin,Particle[ipart].Data()),800,800);
	243	ctmp->cd()->SetLogy();
	244
	245	nsig_data_Proj1->GetXaxis()->SetTitle("#sigma_{TOF}");
	246	nsig_data_Proj1->GetYaxis()->SetTitle("N_{entries}");
	247	nsig_data_Proj1->SetTitle(Form("PionTOF[%.2f,%.2f]",nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)));
	248	nsig_data_Proj1->Draw("E1");
	249	nsig_mc_Proj1->Draw("E1same");
	250	if(saveoption)
	251	ctmp->SaveAs(Form("c%d%s.png",ibin,Particle[ipart].Data()));
	252	delete ctmp;
	253	delete gfundata;
	254	delete gfunmc;
	255	delete nsig_data_Proj1;
	256	delete nsig_mc_Proj1;
	257
	258	}
	259
176a37e3	260	hist1data->Sumw2();
	261	hist1data->Print("all");
	262	hist1data->Sumw2();
	263	hist1mc->Sumw2();
	264	hist2mc->Sumw2();
f887e5c3	265
	266	hist2data->GetXaxis()->SetTitle("p_{T} (GeV/c)");
	267	hist2data->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
	268	hist2mc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
	269	hist2mc->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
	270	hist2data->Divide(hist2data,hist1data,1,1,"B");
	271	hist2mc->Divide(hist2mc,hist1mc,1,1,"B");
	272	lout->Add(hist2data);
	273	lout->Add(hist2mc);
	274
	275	sigmadata->GetXaxis()->SetTitle("p_{T} (GeV/c)");
	276	sigmadata->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
	277	sigmamc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
	278	sigmamc->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
	279
	280
	281
	282	lout->Add(sigmadata);
	283	lout->Add(sigmamc);
	284
	285	TH1D* corr=(TH1D*)hist2data->Clone(Form("corr%s",Particle[ipart].Data()));
	286	corr->Divide(hist2mc);
	287	corr->SetLineColor(kBlue);
	288	corr->GetYaxis()->SetTitle("correction");
176a37e3	289	corr->Fit("pol0","R0","",minpt,maxpt);
f887e5c3	290
	291	oufiletxt<<Particle[ipart].Data()<<" "<<((TF1*)corr->GetListOfFunctions()->At(0))->GetParameter(0)<<endl;
	292	lout->Add(corr);
176a37e3	293	TCanvas* corrcanvas=new TCanvas(Form("corrcanvas%s",Particle[ipart].Data()), Form("corrcanvas%s",Particle[ipart].Data()),600,600);
	294	corrcanvas->SetMargin(0.12,0.05,0.15,0.1);
	295	corr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
	296	corr->GetYaxis()->SetTitle("TOF PID mismatch corr.");
	297	corr->GetYaxis()->SetTitleOffset(1.5*corr->GetYaxis()->GetTitleOffset());
	298	corr->GetYaxis()->SetRangeUser(0.98,1.04);
	299
	300	corrcanvas->cd();
	301	corr->Draw("E1");
	302	corrcanvas->SaveAs(Form("Corr%s.eps",Particle[ipart].Data()));
	303	corrcanvas->SaveAs(Form("Corr%s.png",Particle[ipart].Data()));
	304
f887e5c3	305	}
	306	TFile* fileout=TFile::Open("TOFcorrPID.root","Recreate");
	307
	308	lout->Write("output",TObject::kSingleKey);
	309	fileout->Close();
	310	}