[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / ITSsa / MakeCorrectedITSsaSpectraMultiBin.C

/////////////////////////////////////////////////////////
//Macro to read the results of the 3gausfit for spectra
//E. Biolcati, 13-feb-2010
/////////////////////////////////////////////////////////
#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TLatex.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TF1.h>
#include <TMath.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TSystem.h>
#include <TFile.h>
#include <TStyle.h>
#include <TPad.h>
#endif

void Labella(Int_t i);
void Legenda(TH1 *h2, TH1 *h3, TH1 *h4);

//_______________________________________________________
void MakeCorrectedITSsaSpectraMultiBin(Int_t multibin=0){

	TString dirNameSIM, dirNameDATA;
	TFile *fiSIM  = new TFile(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d/SpectraReco.root",multibin));
	TFile *fiDATA = new TFile(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d/SpectraReco.root",multibin));
	TFile *fout   = new TFile(Form("ITSsaSpectraCorr_MultiBin%d.root",multibin),"recreate");

	//dca correction
	TFile *fPanosPosP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_PosP_expostrange_MultBin%d.root",multibin));  
	TFile *fPanosNegP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_NegP_expostrange_MultBin%d.root",multibin)); 
	TH1F *hPanosPosP= (TH1F*)fPanosPosP->Get("fHistSecTOTCorrDATAMC"); 
	TH1F *hPanosNegP= (TH1F*)fPanosNegP->Get("fHistSecTOTCorrDATAMC"); 
	for(Int_t pbin=0;pbin<=hPanosPosP->GetNbinsX();pbin++)hPanosPosP->SetBinError(pbin,0); 
	for(Int_t pbin=0;pbin<=hPanosNegP->GetNbinsX();pbin++)hPanosNegP->SetBinError(pbin,0);

	//nevts
	TH1F *hstat = (TH1F*)fiDATA->Get("fHistNEvents");
	Double_t NEvts = hstat->GetBinContent(hstat->FindBin(1.));
	cout<<"Event number used for the normalization "<<NEvts<<endl;

	//canvas
	gStyle->SetOptStat(0);
	TCanvas *cs=new TCanvas("cs","cs",1200,700);
	cs->Divide(2,1);
	TCanvas *cmix=new TCanvas("cmix","cmix",800,900);
	cmix->Divide(1,3,0,0);

	TH1F *hPieff[2];
	TH1F *hKeff[2];
	TH1F *hPeff[2];
	TH1F *hPiDATA[2];
	TH1F *hKDATA[2];
	TH1F *hPDATA[2];
	TH1F *hPiDATANorm[2];
	TH1F *hKDATANorm[2];
	TH1F *hPDATANorm[2];
	TH1F *hKsuPi[2]; 
	TH1F *hPsuPi[2]; 
	TH1F *hPsuK[2]; 

	hPieff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg0");
	hKeff[0]  = (TH1F*)fiSIM->Get("hCorrFacNeg1");
	hPeff[0]  = (TH1F*)fiSIM->Get("hCorrFacNeg2");
	hPieff[1] = (TH1F*)fiSIM->Get("hCorrFacPos0");
	hKeff[1]  = (TH1F*)fiSIM->Get("hCorrFacPos1");
	hPeff[1]  = (TH1F*)fiSIM->Get("hCorrFacPos2");
	TCanvas *ceff=new TCanvas("ceff","ceff",800,500);
	ceff->Divide(2,1);
	ceff->cd(1);
	gPad->SetGridy();
	hPieff[0]->Draw();
	hKeff[0]->Draw("same");
	hPeff[0]->Draw("same");
	ceff->cd(2);
	gPad->SetGridy();
	hPieff[1]->Draw();
	hKeff[1]->Draw("same");
	hPeff[1]->Draw("same");

	hPiDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg0");
	hKDATA[0]  = (TH1F*)fiDATA->Get("hSpectraNeg1");
	hPDATA[0]  = (TH1F*)fiDATA->Get("hSpectraNeg2");
	hPiDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos0");
	hKDATA[1]  = (TH1F*)fiDATA->Get("hSpectraPos1");
	hPDATA[1]  = (TH1F*)fiDATA->Get("hSpectraPos2");

	hPiDATANorm[0] = new TH1F(*hPiDATA[0]);
	hPiDATANorm[0]->SetName("hSpectraPiPlusN");
	hPiDATANorm[1] = new TH1F(*hPiDATA[1]);
	hPiDATANorm[1]->SetName("hSpectraPiMinusN");
	hKDATANorm[0]  = new TH1F(*hKDATA[0]);
	hKDATANorm[0]->SetName("hSpectraKPlusN");
	hKDATANorm[1]  = new TH1F(*hKDATA[1]);
	hKDATANorm[1]->SetName("hSpectraKMinusN");
	hPDATANorm[0]  = new TH1F(*hPDATA[0]);
	hPDATANorm[0]->SetName("hSpectraPPlusN");
	hPDATANorm[1]  = new TH1F(*hPDATA[1]);
	hPDATANorm[1]->SetName("hSpectraPMinusN");

	for(Int_t i=0;i<2;i++){ //0==> negative, 1==>positive
		//line colors
		hPiDATA[i] -> SetLineColor(2);
		hKDATA[i]  -> SetLineColor(3);
		hPDATA[i]  -> SetLineColor(4);
		hPiDATA[i] -> SetMarkerStyle(27);
		hKDATA[i]  -> SetMarkerStyle(23);
		hPDATA[i]  -> SetMarkerStyle(24);
		hPiDATA[i] -> SetMarkerColor(2);
		hKDATA[i]  -> SetMarkerColor(3);
		hPDATA[i]  -> SetMarkerColor(4);
		hPiDATA[i] -> SetMarkerStyle(23);
		hKDATA[i]  -> SetMarkerStyle(23);
		hPDATA[i]  -> SetMarkerStyle(23);

		hPiDATANorm[i] -> SetLineColor(2);
		hKDATANorm[i]  -> SetLineColor(3);
		hPDATANorm[i]  -> SetLineColor(4);
		hPiDATANorm[i] -> SetMarkerStyle(27);
		hKDATANorm[i]  -> SetMarkerStyle(23);
		hPDATANorm[i]  -> SetMarkerStyle(24);
		hPiDATANorm[i] -> SetMarkerColor(2);
		hKDATANorm[i]  -> SetMarkerColor(3);
		hPDATANorm[i]  -> SetMarkerColor(4);
		hPiDATANorm[i] -> SetMarkerStyle(23);
		hKDATANorm[i]  -> SetMarkerStyle(23);
		hPDATANorm[i]  -> SetMarkerStyle(23);

		//division for efficiency
		hPiDATA[i] -> Divide(hPieff[i]);
		hKDATA[i]  -> Divide(hKeff[i]);
		hPDATA[i]  -> Divide(hPeff[i]);
		hPiDATANorm[i] -> Divide(hPieff[i]);
		hKDATANorm[i]  -> Divide(hKeff[i]);
		hPDATANorm[i]  -> Divide(hPeff[i]);

		//normalization number of events
		hPiDATANorm[i] -> Scale(1./NEvts);
		hKDATANorm[i]  -> Scale(1./NEvts);
		hPDATANorm[i]  -> Scale(1./NEvts);

		//correction factor based on fit to DCA distr for P and Pibar 
		hPDATA[0]->Multiply(hPanosPosP); 
		hPDATA[1]->Multiply(hPanosNegP);

		//drawing
		cs->cd(i+1);
		gPad->SetLogy();
		gPad->SetGridy();
		gPad->SetGridx();
		Float_t minim=0.01;
		Float_t maxim=10.;
		hPiDATANorm[i]->Draw("e");
		hPiDATANorm[i]->GetYaxis()->SetRangeUser(minim,maxim);
		hPiDATANorm[i]->GetYaxis()->SetTitle("#frac{1}{N}#frac{d^{2}N}{dp_{t}dy}");
		hPiDATANorm[i]->SetTitle("DATA from ITSsa - corrected");
		hKDATANorm[i]->Draw("esames");
		hPDATANorm[i]->Draw("esames");
		Labella(i);
		Legenda(hPiDATA[i],hKDATA[i],hPDATA[i]);
		cs->Update();
	}

	//fluka correction for pi
	TFile *fGeanFlukaPi= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root"));
	TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
	TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
	for(Int_t binPi=0;binPi<=hPiDATA[0]->GetNbinsX();binPi++){
		Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hPiDATA[0]->GetBinCenter(binPi)));
		Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hPiDATA[1]->GetBinCenter(binPi)));
		//cout<<"PiPos  "<<FlukaCorrPiPos<<"  "<<hPiDATA[0]->GetBinCenter(binPi)<<"  " <<binPi <<endl;
		//cout<<"PiNeg  "<<FlukaCorrPiNeg<<"  "<<hPiDATA[0]->GetBinCenter(binPi)<<"  " <<binPi <<endl;
		hPiDATA[0]->SetBinContent(binPi,hPiDATA[0]->GetBinContent(binPi)*FlukaCorrPiPos);
		hPiDATA[1]->SetBinContent(binPi,hPiDATA[1]->GetBinContent(binPi)*FlukaCorrPiNeg);
	}
	//fluka correction for pi
	TFile *fGeanFlukaK= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root"));
	TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
	TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
	for(Int_t binK=0;binK<=hKDATA[0]->GetNbinsX();binK++){
		Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hKDATA[0]->GetBinCenter(binK)));
		Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hKDATA[1]->GetBinCenter(binK)));
		//cout<<"KPos :"<<FlukaCorrKPos<<"  "<<hKDATA[0]->GetBinCenter(binK)<<"  " <<binK <<endl;
		//cout<<"KNeg :"<<FlukaCorrKNeg<<"  "<<hKDATA[0]->GetBinCenter(binK)<<"  " <<binK <<endl;
		hKDATA[0]->SetBinContent(binK,hKDATA[0]->GetBinContent(binK)*FlukaCorrKPos);
		hKDATA[1]->SetBinContent(binK,hKDATA[1]->GetBinContent(binK)*FlukaCorrKNeg);
	}
	//fluka correction for P
	//ITS specific file for protons/antiprotons
	Int_t kPos=0;
	Int_t kNeg=1;
	TFile* fITS = new TFile ("RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
	TH2D * hCorrFlukaITS[2];
	hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
	hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
	for(Int_t icharge = 0; icharge < 2; icharge++){
		Int_t nbins = hPDATA[0]->GetNbinsX();
		Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
		for(Int_t ibin = 0; ibin < nbins; ibin++){
			Float_t pt = hPDATA[0]->GetBinCenter(ibin);
			Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
			Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
			if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
			if (pt > maxPtCorrection) pt = maxPtCorrection;
			Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
			if(icharge==0){
				if (correction != 0) {// If the bin is empty this is a  0
					hPDATA[0]->SetBinContent(ibin,hPDATA[0]->GetBinContent(ibin)*correction);
					hPDATA[0]->SetBinError(ibin,hPDATA[0]->GetBinError  (ibin)*correction);
				}else if (hPDATA[0]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
					cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
					cout << " Bin content: " << hPDATA[0]->GetBinContent(ibin) << endl;
				}
			}
			if(icharge==1){
				if (correction != 0) {// If the bin is empty this is a  0
					hPDATA[1]->SetBinContent(ibin,hPDATA[1]->GetBinContent(ibin)*correction);
					hPDATA[1]->SetBinError(ibin,hPDATA[1]->GetBinError  (ibin)*correction);
				}else if (hPDATA[1]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
					cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
					cout << " Bin content: " << hPDATA[1]->GetBinContent(ibin) << endl;
				}
			}
		}
	}


	//mixed particle ratios
	for(Int_t i=0;i<2;i++){
		hKsuPi[i] = (TH1F*)hKDATA[i]->Clone("KsuPi");
		hPsuPi[i] = (TH1F*)hPDATA[i]->Clone("PsuPi");
		hPsuK[i]  = (TH1F*)hPDATA[i]->Clone("PsuK");
		hKsuPi[i] -> Divide(hPiDATA[i]);
		hPsuPi[i] -> Divide(hPiDATA[i]);
		hPsuK[i]  -> Divide(hKDATA[i]);
	}

	//positive/negative ratios
	TH1F *hPiratio = (TH1F*)hPiDATA[1]->Clone("PionsRatio");
	TH1F *hKratio  = (TH1F*)hKDATA[1]->Clone("KaonsRatio");
	TH1F *	hPratio  = (TH1F*)hPDATA[1]->Clone("ProtonsRatio");
	hPiratio -> Divide(hPiDATA[0]);
	hKratio  -> Divide(hKDATA[0]);
	hPratio  -> Divide(hPDATA[0]);


	//drawing positive/negative ratios
	TCanvas *cratio=new TCanvas("cratio","",980,600);
	cratio->Divide(1,3,0,0);
	cratio->SetBottomMargin(0.08);
	cratio->cd(1);
	gPad->SetGridy();
	hPiratio->GetYaxis()->SetTitle("");
	hPiratio->GetYaxis()->SetRangeUser(0.7,1.3);
	hPiratio->SetTitle("");
	hPiratio->Draw("mp");
	TLatex *ll1=new TLatex(0.7,0.7,"#pi^{+}/#pi^{-}");
	ll1->SetNDC();
	ll1->SetTextSize(0.14);
	ll1->Draw();
	cratio->cd(2);
	gPad->SetGridy();
	hKratio->SetTitle("");
	hKratio->GetYaxis()->SetRangeUser(0.7,1.3);
	hKratio->Draw("mp");
	TLatex *ll2=new TLatex(0.7,0.7,"K^{+}/K^{-}");
	ll2->SetNDC();
	ll2->SetTextSize(0.14);
	ll2->Draw();
	cratio->cd(3);
	gPad->SetGridy();
	hPratio->SetTitle("");
	hPratio->GetYaxis()->SetRangeUser(0.7,1.3);
	hPratio->GetXaxis()->SetTitle("p_{t} [GeV/c]");
	hPratio->Draw("mp");
	TLatex *ll3=new TLatex(0.7,0.7,"p/#bar{p}");
	ll3->SetNDC();
	ll3->SetTextSize(0.144);
	ll3->Draw();

	//drawing mixed particle ratios
	gStyle->SetOptTitle(0);
	cmix->cd(1);
	hKsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
	hKsuPi[0]->Draw();
	hKsuPi[1]->Draw("same");
	hKsuPi[0]->SetMinimum(0);
	hKsuPi[0]->SetMarkerStyle(23);
	hKsuPi[1]->SetMarkerStyle(24);
	TLegend *legm1=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
	legm1->AddEntry(hKsuPi[0],"K^{-}/#pi^{-}","p");
	legm1->AddEntry(hKsuPi[1],"K^{+}/#pi^{+}","p");
	legm1->SetFillColor(0);
	legm1->SetBorderSize(0);
	legm1->Draw();
	hKsuPi[0]->SetMarkerColor(2);
	hKsuPi[1]->SetMarkerColor(4);
	hKsuPi[0]->SetLineColor(2);
	hKsuPi[1]->SetLineColor(4);

	cmix->cd(2);
	hPsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
	hPsuPi[0]->Draw();
	hPsuPi[1]->Draw("same");
	hPsuPi[0]->SetMinimum(0);
	hPsuPi[0]->SetMarkerStyle(23);
	hPsuPi[1]->SetMarkerStyle(24);
	hPsuPi[0]->SetMarkerColor(2);
	hPsuPi[1]->SetMarkerColor(4);
	hPsuPi[0]->SetLineColor(2);
	hPsuPi[1]->SetLineColor(4);
	TLegend *legm2=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
	legm2->AddEntry(hPsuPi[0],"#bar{p}/#pi^{-}","p");
	legm2->AddEntry(hPsuPi[1],"p/#pi^{+}","p");
	legm2->SetFillColor(0);
	legm2->SetBorderSize(0);
	legm2->Draw();

	cmix->cd(3);
	hPsuK[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
	hPsuK[0]->Draw();
	hPsuK[1]->Draw("same");
	hPsuK[0]->SetMinimum(0);
	hPsuK[0]->SetMarkerStyle(23);
	hPsuK[1]->SetMarkerStyle(24);
	hPsuK[0]->SetMarkerColor(2);
	hPsuK[1]->SetMarkerColor(4);
	hPsuK[0]->SetLineColor(2);
	hPsuK[1]->SetLineColor(4);
	TLegend *legm3=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
	legm3->AddEntry(hPsuPi[0],"#bar{p}/K^{-}","p");
	legm3->AddEntry(hPsuPi[1],"p/K^{+}","p");
	legm3->SetFillColor(0);
	legm3->SetBorderSize(0);
	legm3->Draw();

  //save histograms in the root files
	fout->cd();
	hPiDATA[0]->Write();
	hKDATA[0] ->Write();
	hPDATA[0] ->Write();
	hPiDATA[1]->Write();
	hKDATA[1] ->Write();
	hPDATA[1] ->Write();
	hPiDATANorm[0]->Write();
	hKDATANorm[0] ->Write();
	hPDATANorm[0] ->Write();
	hPiDATANorm[1]->Write();
	hKDATANorm[1] ->Write();
	hPDATANorm[1] ->Write();
	//fout->Close();
	return;

}//end of the main


//_______________________________________________________
void Labella(Int_t i){
	Char_t txt[50];
	if(i==0) sprintf(txt,"negative particles");
	else sprintf(txt,"positive particles");
	TLatex *ltx=new TLatex(0.4,0.3,txt);
	ltx->SetNDC();
	ltx->SetTextColor(6);
	ltx->SetTextFont(22);
	ltx->Draw();
	return;
}

//_______________________________________________________
void Legenda(TH1 *h2, TH1 *h3, TH1 *h4){
	TLegend *leg=new TLegend(0.51,0.11,0.84,0.25,NULL,"brNDC");
	leg->SetFillColor(0);
	leg->SetBorderSize(0);
	TLegendEntry *entry2=leg->AddEntry(h2,"pions","p");
	entry2->SetTextColor(2);
	TLegendEntry *entry3=leg->AddEntry(h3,"kaons","p");
	entry3->SetTextColor(3);
	TLegendEntry *entry4=leg->AddEntry(h4,"protons","p");
	entry4->SetTextColor(4);
	leg->Draw("same");
}

//EOF
Commit	Line	Data
8bcdc72e	1	/////////////////////////////////////////////////////////
	2	//Macro to read the results of the 3gausfit for spectra
	3	//E. Biolcati, 13-feb-2010
	4	/////////////////////////////////////////////////////////
	5	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	6	#include <Riostream.h>
	7	#include <TLatex.h>
	8	#include <TH1F.h>
	9	#include <TH2D.h>
	10	#include <TF1.h>
	11	#include <TMath.h>
	12	#include <TLegend.h>
	13	#include <TLegendEntry.h>
	14	#include <TGraphErrors.h>
	15	#include <TGraphAsymmErrors.h>
	16	#include <TTree.h>
	17	#include <TCanvas.h>
	18	#include <TSystem.h>
	19	#include <TFile.h>
	20	#include <TStyle.h>
	21	#include <TPad.h>
	22	#endif
	23
	24	void Labella(Int_t i);
	25	void Legenda(TH1 h2, TH1 h3, TH1 *h4);
	26
	27	//_______________________________________________________
	28	void MakeCorrectedITSsaSpectraMultiBin(Int_t multibin=0){
	29
	30	TString dirNameSIM, dirNameDATA;
	31	TFile *fiSIM = new TFile(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d/SpectraReco.root",multibin));
	32	TFile *fiDATA = new TFile(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d/SpectraReco.root",multibin));
	33	TFile *fout = new TFile(Form("ITSsaSpectraCorr_MultiBin%d.root",multibin),"recreate");
	34
	35	//dca correction
	36	TFile *fPanosPosP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_PosP_expostrange_MultBin%d.root",multibin));
	37	TFile *fPanosNegP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_NegP_expostrange_MultBin%d.root",multibin));
	38	TH1F hPanosPosP= (TH1F)fPanosPosP->Get("fHistSecTOTCorrDATAMC");
	39	TH1F hPanosNegP= (TH1F)fPanosNegP->Get("fHistSecTOTCorrDATAMC");
	40	for(Int_t pbin=0;pbin<=hPanosPosP->GetNbinsX();pbin++)hPanosPosP->SetBinError(pbin,0);
	41	for(Int_t pbin=0;pbin<=hPanosNegP->GetNbinsX();pbin++)hPanosNegP->SetBinError(pbin,0);
	42
	43	//nevts
	44	TH1F hstat = (TH1F)fiDATA->Get("fHistNEvents");
	45	Double_t NEvts = hstat->GetBinContent(hstat->FindBin(1.));
	46	cout<<"Event number used for the normalization "<<NEvts<<endl;
	47
	48	//canvas
	49	gStyle->SetOptStat(0);
	50	TCanvas *cs=new TCanvas("cs","cs",1200,700);
	51	cs->Divide(2,1);
	52	TCanvas *cmix=new TCanvas("cmix","cmix",800,900);
	53	cmix->Divide(1,3,0,0);
	54
	55	TH1F *hPieff[2];
	56	TH1F *hKeff[2];
	57	TH1F *hPeff[2];
	58	TH1F *hPiDATA[2];
	59	TH1F *hKDATA[2];
	60	TH1F *hPDATA[2];
	61	TH1F *hPiDATANorm[2];
	62	TH1F *hKDATANorm[2];
	63	TH1F *hPDATANorm[2];
	64	TH1F *hKsuPi[2];
65	TH1F *hPsuPi[2];
66	TH1F *hPsuK[2];
67
68	hPieff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg0");
69	hKeff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg1");
70	hPeff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg2");
71	hPieff[1] = (TH1F*)fiSIM->Get("hCorrFacPos0");
72	hKeff[1] = (TH1F*)fiSIM->Get("hCorrFacPos1");
73	hPeff[1] = (TH1F*)fiSIM->Get("hCorrFacPos2");
74	TCanvas *ceff=new TCanvas("ceff","ceff",800,500);
75	ceff->Divide(2,1);
76	ceff->cd(1);
77	gPad->SetGridy();
78	hPieff[0]->Draw();
79	hKeff[0]->Draw("same");
80	hPeff[0]->Draw("same");
81	ceff->cd(2);
82	gPad->SetGridy();
83	hPieff[1]->Draw();
84	hKeff[1]->Draw("same");
85	hPeff[1]->Draw("same");
86
87	hPiDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg0");
88	hKDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg1");
89	hPDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg2");
90	hPiDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos0");
91	hKDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos1");
92	hPDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos2");
93
94	hPiDATANorm[0] = new TH1F(*hPiDATA[0]);
95	hPiDATANorm[0]->SetName("hSpectraPiPlusN");
96	hPiDATANorm[1] = new TH1F(*hPiDATA[1]);
97	hPiDATANorm[1]->SetName("hSpectraPiMinusN");
98	hKDATANorm[0] = new TH1F(*hKDATA[0]);
99	hKDATANorm[0]->SetName("hSpectraKPlusN");
100	hKDATANorm[1] = new TH1F(*hKDATA[1]);
101	hKDATANorm[1]->SetName("hSpectraKMinusN");
102	hPDATANorm[0] = new TH1F(*hPDATA[0]);
103	hPDATANorm[0]->SetName("hSpectraPPlusN");
104	hPDATANorm[1] = new TH1F(*hPDATA[1]);
105	hPDATANorm[1]->SetName("hSpectraPMinusN");
106
107	for(Int_t i=0;i<2;i++){ //0==> negative, 1==>positive
108	//line colors
109	hPiDATA[i] -> SetLineColor(2);
110	hKDATA[i] -> SetLineColor(3);
111	hPDATA[i] -> SetLineColor(4);
112	hPiDATA[i] -> SetMarkerStyle(27);
113	hKDATA[i] -> SetMarkerStyle(23);
114	hPDATA[i] -> SetMarkerStyle(24);
115	hPiDATA[i] -> SetMarkerColor(2);
116	hKDATA[i] -> SetMarkerColor(3);
117	hPDATA[i] -> SetMarkerColor(4);
118	hPiDATA[i] -> SetMarkerStyle(23);
119	hKDATA[i] -> SetMarkerStyle(23);
120	hPDATA[i] -> SetMarkerStyle(23);
121
122	hPiDATANorm[i] -> SetLineColor(2);
123	hKDATANorm[i] -> SetLineColor(3);
124	hPDATANorm[i] -> SetLineColor(4);
125	hPiDATANorm[i] -> SetMarkerStyle(27);
126	hKDATANorm[i] -> SetMarkerStyle(23);
127	hPDATANorm[i] -> SetMarkerStyle(24);
128	hPiDATANorm[i] -> SetMarkerColor(2);
129	hKDATANorm[i] -> SetMarkerColor(3);
130	hPDATANorm[i] -> SetMarkerColor(4);
131	hPiDATANorm[i] -> SetMarkerStyle(23);
132	hKDATANorm[i] -> SetMarkerStyle(23);
133	hPDATANorm[i] -> SetMarkerStyle(23);
134
135	//division for efficiency
136	hPiDATA[i] -> Divide(hPieff[i]);
137	hKDATA[i] -> Divide(hKeff[i]);
138	hPDATA[i] -> Divide(hPeff[i]);
139	hPiDATANorm[i] -> Divide(hPieff[i]);
140	hKDATANorm[i] -> Divide(hKeff[i]);
141	hPDATANorm[i] -> Divide(hPeff[i]);
142
143	//normalization number of events
144	hPiDATANorm[i] -> Scale(1./NEvts);
145	hKDATANorm[i] -> Scale(1./NEvts);
146	hPDATANorm[i] -> Scale(1./NEvts);
147
148	//correction factor based on fit to DCA distr for P and Pibar
149	hPDATA[0]->Multiply(hPanosPosP);
150	hPDATA[1]->Multiply(hPanosNegP);
151
152	//drawing
153	cs->cd(i+1);
154	gPad->SetLogy();
155	gPad->SetGridy();
156	gPad->SetGridx();
157	Float_t minim=0.01;
158	Float_t maxim=10.;
159	hPiDATANorm[i]->Draw("e");
160	hPiDATANorm[i]->GetYaxis()->SetRangeUser(minim,maxim);
161	hPiDATANorm[i]->GetYaxis()->SetTitle("#frac{1}{N}#frac{d^{2}N}{dp_{t}dy}");
162	hPiDATANorm[i]->SetTitle("DATA from ITSsa - corrected");
163	hKDATANorm[i]->Draw("esames");
164	hPDATANorm[i]->Draw("esames");
165	Labella(i);
166	Legenda(hPiDATA[i],hKDATA[i],hPDATA[i]);
167	cs->Update();
168	}
169
170	//fluka correction for pi
171	TFile *fGeanFlukaPi= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root"));
172	TH1F hGeantFlukaPiPos=(TH1F)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
173	TH1F hGeantFlukaPiNeg=(TH1F)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
174	for(Int_t binPi=0;binPi<=hPiDATA[0]->GetNbinsX();binPi++){
175	Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hPiDATA[0]->GetBinCenter(binPi)));
176	Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hPiDATA[1]->GetBinCenter(binPi)));
177	//cout<<"PiPos "<<FlukaCorrPiPos<<" "<<hPiDATA[0]->GetBinCenter(binPi)<<" " <<binPi <<endl;
178	//cout<<"PiNeg "<<FlukaCorrPiNeg<<" "<<hPiDATA[0]->GetBinCenter(binPi)<<" " <<binPi <<endl;
179	hPiDATA[0]->SetBinContent(binPi,hPiDATA[0]->GetBinContent(binPi)*FlukaCorrPiPos);
180	hPiDATA[1]->SetBinContent(binPi,hPiDATA[1]->GetBinContent(binPi)*FlukaCorrPiNeg);
181	}
182	//fluka correction for pi
183	TFile *fGeanFlukaK= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root"));
184	TH1F hGeantFlukaKPos=(TH1F)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
185	TH1F hGeantFlukaKNeg=(TH1F)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
186	for(Int_t binK=0;binK<=hKDATA[0]->GetNbinsX();binK++){
187	Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hKDATA[0]->GetBinCenter(binK)));
188	Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hKDATA[1]->GetBinCenter(binK)));
189	//cout<<"KPos :"<<FlukaCorrKPos<<" "<<hKDATA[0]->GetBinCenter(binK)<<" " <<binK <<endl;
190	//cout<<"KNeg :"<<FlukaCorrKNeg<<" "<<hKDATA[0]->GetBinCenter(binK)<<" " <<binK <<endl;
191	hKDATA[0]->SetBinContent(binK,hKDATA[0]->GetBinContent(binK)*FlukaCorrKPos);
192	hKDATA[1]->SetBinContent(binK,hKDATA[1]->GetBinContent(binK)*FlukaCorrKNeg);
193	}
194	//fluka correction for P
195	//ITS specific file for protons/antiprotons
196	Int_t kPos=0;
197	Int_t kNeg=1;
198	TFile* fITS = new TFile ("RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
199	TH2D * hCorrFlukaITS[2];
200	hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
201	hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
202	for(Int_t icharge = 0; icharge < 2; icharge++){
203	Int_t nbins = hPDATA[0]->GetNbinsX();
204	Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
205	for(Int_t ibin = 0; ibin < nbins; ibin++){
206	Float_t pt = hPDATA[0]->GetBinCenter(ibin);
207	Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
208	Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
209	if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
210	if (pt > maxPtCorrection) pt = maxPtCorrection;
211	Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
212	if(icharge==0){
213	if (correction != 0) {// If the bin is empty this is a 0
214	hPDATA[0]->SetBinContent(ibin,hPDATA[0]->GetBinContent(ibin)*correction);
215	hPDATA[0]->SetBinError(ibin,hPDATA[0]->GetBinError (ibin)*correction);
216	}else if (hPDATA[0]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
217	cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
218	cout << " Bin content: " << hPDATA[0]->GetBinContent(ibin) << endl;
219	}
220	}
221	if(icharge==1){
222	if (correction != 0) {// If the bin is empty this is a 0
223	hPDATA[1]->SetBinContent(ibin,hPDATA[1]->GetBinContent(ibin)*correction);
224	hPDATA[1]->SetBinError(ibin,hPDATA[1]->GetBinError (ibin)*correction);
225	}else if (hPDATA[1]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
226	cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
227	cout << " Bin content: " << hPDATA[1]->GetBinContent(ibin) << endl;
228	}
229	}
230	}
231	}
232
233
234	//mixed particle ratios
235	for(Int_t i=0;i<2;i++){
236	hKsuPi[i] = (TH1F*)hKDATA[i]->Clone("KsuPi");
237	hPsuPi[i] = (TH1F*)hPDATA[i]->Clone("PsuPi");
238	hPsuK[i] = (TH1F*)hPDATA[i]->Clone("PsuK");
239	hKsuPi[i] -> Divide(hPiDATA[i]);
240	hPsuPi[i] -> Divide(hPiDATA[i]);
241	hPsuK[i] -> Divide(hKDATA[i]);
242	}
243
244	//positive/negative ratios
245	TH1F hPiratio = (TH1F)hPiDATA[1]->Clone("PionsRatio");
246	TH1F hKratio = (TH1F)hKDATA[1]->Clone("KaonsRatio");
247	TH1F * hPratio = (TH1F*)hPDATA[1]->Clone("ProtonsRatio");
248	hPiratio -> Divide(hPiDATA[0]);
249	hKratio -> Divide(hKDATA[0]);
250	hPratio -> Divide(hPDATA[0]);
251
252
253	//drawing positive/negative ratios
254	TCanvas *cratio=new TCanvas("cratio","",980,600);
255	cratio->Divide(1,3,0,0);
256	cratio->SetBottomMargin(0.08);
257	cratio->cd(1);
258	gPad->SetGridy();
259	hPiratio->GetYaxis()->SetTitle("");
260	hPiratio->GetYaxis()->SetRangeUser(0.7,1.3);
261	hPiratio->SetTitle("");
262	hPiratio->Draw("mp");
263	TLatex *ll1=new TLatex(0.7,0.7,"#pi^{+}/#pi^{-}");
264	ll1->SetNDC();
265	ll1->SetTextSize(0.14);
266	ll1->Draw();
267	cratio->cd(2);
268	gPad->SetGridy();
269	hKratio->SetTitle("");
270	hKratio->GetYaxis()->SetRangeUser(0.7,1.3);
271	hKratio->Draw("mp");
272	TLatex *ll2=new TLatex(0.7,0.7,"K^{+}/K^{-}");
273	ll2->SetNDC();
274	ll2->SetTextSize(0.14);
275	ll2->Draw();
276	cratio->cd(3);
277	gPad->SetGridy();
278	hPratio->SetTitle("");
279	hPratio->GetYaxis()->SetRangeUser(0.7,1.3);
280	hPratio->GetXaxis()->SetTitle("p_{t} [GeV/c]");
281	hPratio->Draw("mp");
282	TLatex *ll3=new TLatex(0.7,0.7,"p/#bar{p}");
283	ll3->SetNDC();
284	ll3->SetTextSize(0.144);
285	ll3->Draw();
286
287	//drawing mixed particle ratios
288	gStyle->SetOptTitle(0);
289	cmix->cd(1);
290	hKsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
291	hKsuPi[0]->Draw();
292	hKsuPi[1]->Draw("same");
293	hKsuPi[0]->SetMinimum(0);
294	hKsuPi[0]->SetMarkerStyle(23);
295	hKsuPi[1]->SetMarkerStyle(24);
296	TLegend *legm1=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
297	legm1->AddEntry(hKsuPi[0],"K^{-}/#pi^{-}","p");
298	legm1->AddEntry(hKsuPi[1],"K^{+}/#pi^{+}","p");
299	legm1->SetFillColor(0);
300	legm1->SetBorderSize(0);
301	legm1->Draw();
302	hKsuPi[0]->SetMarkerColor(2);
303	hKsuPi[1]->SetMarkerColor(4);
304	hKsuPi[0]->SetLineColor(2);
305	hKsuPi[1]->SetLineColor(4);
306
307	cmix->cd(2);
308	hPsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
309	hPsuPi[0]->Draw();
310	hPsuPi[1]->Draw("same");
311	hPsuPi[0]->SetMinimum(0);
312	hPsuPi[0]->SetMarkerStyle(23);
313	hPsuPi[1]->SetMarkerStyle(24);
314	hPsuPi[0]->SetMarkerColor(2);
315	hPsuPi[1]->SetMarkerColor(4);
316	hPsuPi[0]->SetLineColor(2);
317	hPsuPi[1]->SetLineColor(4);
318	TLegend *legm2=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
319	legm2->AddEntry(hPsuPi[0],"#bar{p}/#pi^{-}","p");
320	legm2->AddEntry(hPsuPi[1],"p/#pi^{+}","p");
321	legm2->SetFillColor(0);
322	legm2->SetBorderSize(0);
323	legm2->Draw();
324
325	cmix->cd(3);
326	hPsuK[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
327	hPsuK[0]->Draw();
328	hPsuK[1]->Draw("same");
329	hPsuK[0]->SetMinimum(0);
330	hPsuK[0]->SetMarkerStyle(23);
331	hPsuK[1]->SetMarkerStyle(24);
332	hPsuK[0]->SetMarkerColor(2);
333	hPsuK[1]->SetMarkerColor(4);
334	hPsuK[0]->SetLineColor(2);
335	hPsuK[1]->SetLineColor(4);
336	TLegend *legm3=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
337	legm3->AddEntry(hPsuPi[0],"#bar{p}/K^{-}","p");
338	legm3->AddEntry(hPsuPi[1],"p/K^{+}","p");
339	legm3->SetFillColor(0);
340	legm3->SetBorderSize(0);
341	legm3->Draw();
342
343	//save histograms in the root files
344	fout->cd();
345	hPiDATA[0]->Write();
346	hKDATA[0] ->Write();
347	hPDATA[0] ->Write();
348	hPiDATA[1]->Write();
349	hKDATA[1] ->Write();
350	hPDATA[1] ->Write();
351	hPiDATANorm[0]->Write();
352	hKDATANorm[0] ->Write();
353	hPDATANorm[0] ->Write();
354	hPiDATANorm[1]->Write();
355	hKDATANorm[1] ->Write();
356	hPDATANorm[1] ->Write();
357	//fout->Close();
358	return;
359
360	}//end of the main
361
362
363	//_______________________________________________________
364	void Labella(Int_t i){
365	Char_t txt[50];
366	if(i==0) sprintf(txt,"negative particles");
367	else sprintf(txt,"positive particles");
368	TLatex *ltx=new TLatex(0.4,0.3,txt);
369	ltx->SetNDC();
370	ltx->SetTextColor(6);
371	ltx->SetTextFont(22);
372	ltx->Draw();
373	return;
374	}
375
376	//_______________________________________________________
377	void Legenda(TH1 h2, TH1 h3, TH1 *h4){
378	TLegend *leg=new TLegend(0.51,0.11,0.84,0.25,NULL,"brNDC");
379	leg->SetFillColor(0);
380	leg->SetBorderSize(0);
381	TLegendEntry *entry2=leg->AddEntry(h2,"pions","p");
382	entry2->SetTextColor(2);
383	TLegendEntry *entry3=leg->AddEntry(h3,"kaons","p");
384	entry3->SetTextColor(3);
385	TLegendEntry *entry4=leg->AddEntry(h4,"protons","p");
386	entry4->SetTextColor(4);
387	leg->Draw("same");
388	}
389
390	//EOF
391