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

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TFile.h>
#include <TF1.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLatex.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TFractionFitter.h>
#endif

enum species {kPion, kKaon, kProton};
enum chargesign {kPositive,kNegative};

Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax);

void FitDCADistributions(TString period="LHC10d",
			 TString MCname="LHC10f6a",  
			 Int_t iSpecies=2, 
			 Int_t cSign=1
			 )
{
  gROOT->SetStyle("Plain");
  gStyle->SetFillColor(0);
  gStyle->SetOptStat(0000);
  //binning
  const Int_t nptbins = 22;
  Double_t ptbins[nptbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
  //Final correction
  TH1F *hPrimAllDATA=new TH1F("hPrimAllDATA","hPrimAllDATA",nptbins,ptbins);
  TH1F *hPrimAllMC=new TH1F("hPrimAllMC","hPrimAllMC",nptbins,ptbins);
  TH1F *hPrimAllDATAMC=new TH1F("hPrimAllDATAMC","hPrimAllDATAMC",nptbins,ptbins);
  
  
  TString fnameMC=Form("%s/AnalysisResults.root",MCname.Data());
  TString lnameMC="clistITSsaMult-1to-1";
  TString lnameCutMC="DCACutMult-1to-1";
  TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
  TString lnameDATA="clistITSsaMult-1to-1";
  TString lnameCutDATA="DCACutMult-1to-1";

  const Int_t firstbin=1;
  Int_t rebin=10;
  Bool_t optSmooth=kTRUE;


  TString pCharge="Pos";
  TString pcharge="pos";
  if(cSign==kNegative){
    pCharge="Neg";
    pcharge="neg";
  }
  TString species[3]={"Pi","K","P"};
  TString speciesRoot[6]={"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"};
  Int_t idPart=iSpecies*2+cSign;
  TString partname=Form("%s%s",pCharge.Data(),species[iSpecies].Data());
  printf("%s\n",partname.Data());

  // MC histograms 
  TFile *fMC=new TFile(fnameMC,"READ");  
  TH1F *fHistDCAMC[nptbins]; 
  TH1F *fHistDCAPrim[nptbins]; 
  TH1F *fHistDCAsec[nptbins]; 
  TH1F *fHistDCAsecSt[nptbins]; 
  TDirectory *dirFileMC=(TDirectory*)fMC->Get("PWG2SpectraITSsa");
  TList *liMC = (TList*)dirFileMC->Get(lnameMC.Data());
  for(Int_t m=0;m<nptbins;m++){
    fHistDCAMC[m] = (TH1F*)liMC->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
    fHistDCAPrim[m] = (TH1F*)liMC->FindObject(Form("fHistMCPrimDCA%s%i",partname.Data(),m));
    fHistDCAsec[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecMatDCA%s%i",partname.Data(),m));
    fHistDCAsecSt[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecStDCA%s%i",partname.Data(),m));
    fHistDCAMC[m]->Rebin(rebin);
    fHistDCAPrim[m]->Rebin(rebin);
    fHistDCAsec[m]->Rebin(rebin);
    fHistDCAsecSt[m]->Rebin(rebin);
  }
  TList *liCutMC = (TList*)dirFileMC->Get(lnameCutMC.Data());
  Bool_t setDCA=kFALSE;

  Float_t DCAcut=7.;
  Double_t xyPMC[3];   //parameters of DCA cut
  xyPMC[0]=36.; //MC LHC10d1
  xyPMC[1]=43.9;
  xyPMC[2]=1.3;
  TF1* fDCAxyCutMC=0x0;
  TF1* fDCAzCutMC=0x0;
  Float_t nSigmaDCAxyMC=0.;
  Float_t nSigmaDCAzMC=0.;
  if(liCutMC){
    fDCAxyCutMC=(TF1*)liCutMC->FindObject("fDCAxyCutFunc");
    fDCAzCutMC=(TF1*)liCutMC->FindObject("fDCAzCutFunc");
    nSigmaDCAxyMC=fDCAxyCutMC->GetParameter(3);
    nSigmaDCAzMC=fDCAzCutMC->GetParameter(3);
    printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutMC->GetParameter(0),fDCAxyCutMC->GetParameter(1),fDCAxyCutMC->GetParameter(2),nSigmaDCAxyMC);
    setDCA=kTRUE;
  }
  if(!setDCA){
    printf("DCA cut values for MC not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
    fDCAxyCutMC = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
    for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutMC->SetParameter(ipar,xyPMC[ipar]);
    fDCAxyCutMC->SetParameter(3,DCAcut);
    nSigmaDCAxyMC=DCAcut;
  }

  //DATA histograms
  TFile *fDATA=new TFile(fnameDATA,"READ");
  TDirectory *dirFileDATA=(TDirectory*)fDATA->Get("PWG2SpectraITSsa");
  TList *liDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
  TH1F *fHistDCADATA[nptbins]; 
  for(Int_t m=0;m<nptbins;m++){
    fHistDCADATA[m] = (TH1F*)liDATA->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
    fHistDCADATA[m]->Rebin(rebin);
  }
  TList *liCutDATA = (TList*)dirFileDATA->Get(lnameCutDATA.Data());
  setDCA=kFALSE;
  Float_t nSigmaDCAxyDATA=0.;
  Float_t nSigmaDCAzDATA=0.;
  Double_t xyP[3];   // parameters of DCA cut
  xyP[0]=32.7;  //DATA 7 TeV pass2
  xyP[1]=44.8;
  xyP[2]=1.3;
  TF1* fDCAxyCutDATA=0x0;
  TF1* fDCAzCutDATA=0x0;
  if(liCutDATA){
    fDCAxyCutDATA=(TF1*)liCutDATA->FindObject("fDCAxyCutFunc");
    fDCAzCutDATA=(TF1*)liCutDATA->FindObject("fDCAzCutFunc");
    nSigmaDCAxyDATA=fDCAxyCutDATA->GetParameter(3);
    nSigmaDCAzDATA=fDCAzCutDATA->GetParameter(3);
    printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutDATA->GetParameter(0),fDCAxyCutDATA->GetParameter(1),fDCAxyCutDATA->GetParameter(2),nSigmaDCAxyDATA);
    setDCA=kTRUE;
  }
  if(!setDCA){
    printf("DCA cut values for DATA not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
    fDCAxyCutDATA = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
    for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutDATA->SetParameter(ipar,xyP[ipar]);
    fDCAxyCutDATA->SetParameter(3,DCAcut);
    nSigmaDCAxyDATA=DCAcut;
  }
  
  if(TMath::Abs(nSigmaDCAxyDATA-nSigmaDCAxyMC)<0.01) DCAcut=nSigmaDCAxyDATA;
  else{
    printf("ERROR: DCAxy cuts do not match between data and MC\n");
    return;
  }
  if(TMath::Abs(nSigmaDCAzDATA-nSigmaDCAzMC)>0.01){
    printf("ERROR: DCAz cuts do not match between data (%f) and MC (%f) \n",nSigmaDCAzDATA,nSigmaDCAzMC);
    return;
  }

  TCanvas *cfitDATA=new TCanvas("cfitDATA","cfitDATA");
  cfitDATA->Divide(6,4);
  TCanvas *cfitMC=new TCanvas("cfitMC","cfitMC");
  cfitMC->Divide(6,4);

  for(Int_t ibin=firstbin;ibin<nptbins;ibin++){    
    
    Double_t xyMax =fDCAxyCutDATA->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;
    Double_t xyMaxMC =fDCAxyCutMC->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;

    TObjArray *mcTemplates = 0x0;
    if(partname=="PosP"){
      mcTemplates = new TObjArray(3);        // MC histograms are put in this array
      mcTemplates->Add(fHistDCAPrim[ibin]);
      mcTemplates->Add(fHistDCAsecSt[ibin]);
      mcTemplates->Add(fHistDCAsec[ibin]);
    }else{
      mcTemplates = new TObjArray(2);        // MC histograms are put in this array
      mcTemplates->Add(fHistDCAPrim[ibin]);
      mcTemplates->Add(fHistDCAsecSt[ibin]);
    }
    //    if(fHistDCADATA[ibin]->Integral()==0) continue;

    ////////////////// Fit on DATA
    cfitDATA->cd(ibin);
    gPad->SetLogy();
    Double_t primAllDATA=PerformFit(fHistDCADATA[ibin],mcTemplates,partname,xyMax);
    
    ////////////////// Fit on MC
    cfitMC->cd(ibin);
    gPad->SetLogy();
    Double_t primAllMC=PerformFit(fHistDCAMC[ibin],mcTemplates,partname,xyMaxMC);
    
    if(hPrimAllMC!=0 && hPrimAllDATA!=0){
      hPrimAllDATA->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA);
      hPrimAllMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllMC);
      hPrimAllDATAMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA/primAllMC);
    }
  }
  
  //Adding MC truth
  TH1F* hAll=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaMean%d",pCharge.Data(),iSpecies));
  TH1F* hPrim=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaPrimMean%d",pCharge.Data(),iSpecies));
  TH1F* hPrimRecMC=(TH1F*)liMC->FindObject(Form("fHistPrimMC%sReco%d",pcharge.Data(),iSpecies));
  TH1F* hSecMatRecMC=(TH1F*)liMC->FindObject(Form("fHistSecMatMC%sReco%d",pcharge.Data(),iSpecies));
  TH1F* hSecStrRecMC=(TH1F*)liMC->FindObject(Form("fHistSecStrMC%sReco%d",pcharge.Data(),iSpecies));
  hSecStrRecMC->Add(hSecMatRecMC);
  hSecStrRecMC->Add(hPrimRecMC);
  hPrim->SetLineColor(8);
  hPrim->SetLineWidth(2);
  hPrim->SetTitle("Prim/All from MC Truth");
  hPrim->Divide(hAll);
  hPrimRecMC->SetLineColor(11);
  hPrimRecMC->SetLineWidth(2);
  hPrimRecMC->SetTitle("Prim/All from MC Truth Reco");
  hPrimRecMC->Divide(hSecStrRecMC);

  hPrimAllDATAMC->GetXaxis()->SetRangeUser(0.09,0.79);
  hPrimAllMC->GetXaxis()->SetRangeUser(0.09,0.79);
  hPrimAllDATA->GetXaxis()->SetRangeUser(0.09,0.79);
  if(optSmooth) hPrimAllDATAMC->Smooth(1,"R");
  
  hPrimAllDATA->SetTitle("Prim/all data");
  hPrimAllDATA->SetLineColor(2);
  hPrimAllDATA->SetLineWidth(2);
  hPrimAllMC->SetTitle("Prim/all mc");
  hPrimAllMC->SetLineColor(4);
  hPrimAllMC->SetLineWidth(2);
  hPrimAllDATAMC->SetTitle("Prim/all data/mc");
  hPrimAllDATAMC->SetLineColor(1);
  hPrimAllDATAMC->SetLineWidth(2);

  TCanvas *cPrimAll=new TCanvas("cPrimAll","cPrimAll");
  hPrimAllDATA->SetMinimum(0.7);
  hPrimAllDATA->SetMaximum(1.1);
  hPrimAllDATA->Draw("l");
  hPrimAllMC->Draw("lsame");
  hPrimAllDATAMC->Draw("lsame");
  hPrim->Draw("lsame");
  hPrimRecMC->Draw("lsame");
  gPad->BuildLegend();
  TLatex* tsp=new TLatex(0.4,0.75,speciesRoot[idPart].Data());
  tsp->SetTextFont(43);
  tsp->SetTextSize(28);
  tsp->SetNDC();
  tsp->Draw();
  cPrimAll->Update();

  TString fout;
  fout=Form("DCACorr%s_%s_%.0fDCA_%s_TFraction.root",period.Data(),MCname.Data(),DCAcut,partname.Data());
  TFile *out=new TFile(fout.Data(),"recreate");
  hPrimAllDATA->Write();
  hPrimAllMC->Write();
  hPrimAllDATAMC->Write();
  out->Close();
  delete out;
  
}

Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax){
  
  Double_t prim=0,secSt=0,sec=0;
  TFractionFitter* fit = new TFractionFitter(fHistDCA, mc); // initialise
  fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  if(partname=="PosP")fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
  
  //fit->SetRangeX(200,1000);                    // use only the first 15 bins in the fit
  Int_t status = fit->Fit();               // perform the fit
  cout << "fit status: " << status << endl;
  
  if (status == 0) {                       // check on fit status
    TH1F* result = (TH1F*) fit->GetPlot();
    TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
    TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
    TH1F* secMCPred=0x0;
    if(partname=="PosP"){
      secMCPred=(TH1F*)fit->GetMCPrediction(2);
      secMCPred->SetLineColor(4);
    }
    //Drawing section
    PrimMCPred->SetLineColor(2);
    secStMCPred->SetLineColor(6);
    fHistDCA->SetTitle("DCA distr - data");
    fHistDCA->DrawCopy("Ep");
    result->SetTitle("Fit result");
    result->DrawCopy("same");
    Double_t value,error;
    
    fit->GetResult(0,value,error);
    PrimMCPred->Scale(fHistDCA->GetSumOfWeights()*value/PrimMCPred->GetSumOfWeights());
    PrimMCPred->SetTitle("Primaries");
    PrimMCPred->DrawCopy("same");
    fit->GetResult(1,value,error);
    secStMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secStMCPred->GetSumOfWeights());
    secStMCPred->SetTitle("Sec from strangeness");
    secStMCPred->DrawCopy("same");
    if(partname=="PosP" && secMCPred){
      fit->GetResult(2,value,error);
      secMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secMCPred->GetSumOfWeights());
      secMCPred->SetTitle("Sec from material");
      secMCPred->DrawCopy("same");
    }
    prim=PrimMCPred->Integral(PrimMCPred->FindBin(-xyMax),PrimMCPred->FindBin(xyMax));
    secSt=secStMCPred->Integral(secStMCPred->FindBin(-xyMax),secStMCPred->FindBin(xyMax));
    if(partname=="PosP")sec=secMCPred->Integral(secMCPred->FindBin(-xyMax),secMCPred->FindBin(xyMax));
  }
  else{
    prim=1;
    secSt=1;
    if(partname=="PosP")sec=1;
  }
  Printf("Yields:  primary=%f material=%f strange=%f\n",prim,sec,secSt);
  return prim/(prim+secSt+sec);
}
Commit	Line	Data
a1c5e4a1	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TFile.h>
	3	#include <TF1.h>
	4	#include <TH1F.h>
	5	#include <TCanvas.h>
	6	#include <TLatex.h>
	7	#include <TROOT.h>
	8	#include <TStyle.h>
	9	#include <TFractionFitter.h>
	10	#endif
	11
	12	enum species {kPion, kKaon, kProton};
	13	enum chargesign {kPositive,kNegative};
	14
	15	Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax);
	16
	17	void FitDCADistributions(TString period="LHC10d",
982ff3fe	18	TString MCname="LHC10f6a",
a1c5e4a1	19	Int_t iSpecies=2,
	20	Int_t cSign=1
	21	)
	22	{
	23	gROOT->SetStyle("Plain");
	24	gStyle->SetFillColor(0);
	25	gStyle->SetOptStat(0000);
	26	//binning
	27	const Int_t nptbins = 22;
	28	Double_t ptbins[nptbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
	29	//Final correction
	30	TH1F *hPrimAllDATA=new TH1F("hPrimAllDATA","hPrimAllDATA",nptbins,ptbins);
	31	TH1F *hPrimAllMC=new TH1F("hPrimAllMC","hPrimAllMC",nptbins,ptbins);
	32	TH1F *hPrimAllDATAMC=new TH1F("hPrimAllDATAMC","hPrimAllDATAMC",nptbins,ptbins);
	33
	34
982ff3fe	35	TString fnameMC=Form("%s/AnalysisResults.root",MCname.Data());
a1c5e4a1	36	TString lnameMC="clistITSsaMult-1to-1";
a1c5e4a1	37	TString lnameCutMC="DCACutMult-1to-1";
982ff3fe	38	TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
a1c5e4a1	39	TString lnameDATA="clistITSsaMult-1to-1";
	40	TString lnameCutDATA="DCACutMult-1to-1";
	41
	42	const Int_t firstbin=1;
	43	Int_t rebin=10;
	44	Bool_t optSmooth=kTRUE;
	45
	46
	47	TString pCharge="Pos";
	48	TString pcharge="pos";
	49	if(cSign==kNegative){
	50	pCharge="Neg";
	51	pcharge="neg";
	52	}
	53	TString species[3]={"Pi","K","P"};
	54	TString speciesRoot[6]={"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"};
	55	Int_t idPart=iSpecies*2+cSign;
	56	TString partname=Form("%s%s",pCharge.Data(),species[iSpecies].Data());
	57	printf("%s\n",partname.Data());
	58
	59	// MC histograms
	60	TFile *fMC=new TFile(fnameMC,"READ");
	61	TH1F *fHistDCAMC[nptbins];
	62	TH1F *fHistDCAPrim[nptbins];
	63	TH1F *fHistDCAsec[nptbins];
	64	TH1F *fHistDCAsecSt[nptbins];
	65	TDirectory dirFileMC=(TDirectory)fMC->Get("PWG2SpectraITSsa");
	66	TList liMC = (TList)dirFileMC->Get(lnameMC.Data());
	67	for(Int_t m=0;m<nptbins;m++){
	68	fHistDCAMC[m] = (TH1F*)liMC->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
	69	fHistDCAPrim[m] = (TH1F*)liMC->FindObject(Form("fHistMCPrimDCA%s%i",partname.Data(),m));
	70	fHistDCAsec[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecMatDCA%s%i",partname.Data(),m));
	71	fHistDCAsecSt[m] = (TH1F*)liMC->FindObject(Form("fHistMCSecStDCA%s%i",partname.Data(),m));
	72	fHistDCAMC[m]->Rebin(rebin);
	73	fHistDCAPrim[m]->Rebin(rebin);
	74	fHistDCAsec[m]->Rebin(rebin);
	75	fHistDCAsecSt[m]->Rebin(rebin);
	76	}
	77	TList liCutMC = (TList)dirFileMC->Get(lnameCutMC.Data());
	78	Bool_t setDCA=kFALSE;
	79
	80	Float_t DCAcut=7.;
	81	Double_t xyPMC[3]; //parameters of DCA cut
	82	xyPMC[0]=36.; //MC LHC10d1
	83	xyPMC[1]=43.9;
	84	xyPMC[2]=1.3;
	85	TF1* fDCAxyCutMC=0x0;
	86	TF1* fDCAzCutMC=0x0;
	87	Float_t nSigmaDCAxyMC=0.;
	88	Float_t nSigmaDCAzMC=0.;
	89	if(liCutMC){
	90	fDCAxyCutMC=(TF1*)liCutMC->FindObject("fDCAxyCutFunc");
	91	fDCAzCutMC=(TF1*)liCutMC->FindObject("fDCAzCutFunc");
	92	nSigmaDCAxyMC=fDCAxyCutMC->GetParameter(3);
	93	nSigmaDCAzMC=fDCAzCutMC->GetParameter(3);
	94	printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutMC->GetParameter(0),fDCAxyCutMC->GetParameter(1),fDCAxyCutMC->GetParameter(2),nSigmaDCAxyMC);
	95	setDCA=kTRUE;
	96	}
	97	if(!setDCA){
	98	printf("DCA cut values for MC not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
	99	fDCAxyCutMC = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
	100	for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutMC->SetParameter(ipar,xyPMC[ipar]);
	101	fDCAxyCutMC->SetParameter(3,DCAcut);
	102	nSigmaDCAxyMC=DCAcut;
982ff3fe	103	}
a1c5e4a1	104
	105	//DATA histograms
	106	TFile *fDATA=new TFile(fnameDATA,"READ");
	107	TDirectory dirFileDATA=(TDirectory)fDATA->Get("PWG2SpectraITSsa");
	108	TList liDATA = (TList)dirFileDATA->Get(lnameDATA.Data());
	109	TH1F *fHistDCADATA[nptbins];
	110	for(Int_t m=0;m<nptbins;m++){
	111	fHistDCADATA[m] = (TH1F*)liDATA->FindObject(Form("fHistDCA%s%i",partname.Data(),m));
	112	fHistDCADATA[m]->Rebin(rebin);
	113	}
	114	TList liCutDATA = (TList)dirFileDATA->Get(lnameCutDATA.Data());
	115	setDCA=kFALSE;
	116	Float_t nSigmaDCAxyDATA=0.;
	117	Float_t nSigmaDCAzDATA=0.;
	118	Double_t xyP[3]; // parameters of DCA cut
	119	xyP[0]=32.7; //DATA 7 TeV pass2
	120	xyP[1]=44.8;
	121	xyP[2]=1.3;
	122	TF1* fDCAxyCutDATA=0x0;
	123	TF1* fDCAzCutDATA=0x0;
	124	if(liCutDATA){
	125	fDCAxyCutDATA=(TF1*)liCutDATA->FindObject("fDCAxyCutFunc");
	126	fDCAzCutDATA=(TF1*)liCutDATA->FindObject("fDCAzCutFunc");
	127	nSigmaDCAxyDATA=fDCAxyCutDATA->GetParameter(3);
	128	nSigmaDCAzDATA=fDCAzCutDATA->GetParameter(3);
	129	printf("DCA cut parameters from TF1: %f %f %f --- DCA cut = %f\n",fDCAxyCutDATA->GetParameter(0),fDCAxyCutDATA->GetParameter(1),fDCAxyCutDATA->GetParameter(2),nSigmaDCAxyDATA);
	130	setDCA=kTRUE;
	131	}
	132	if(!setDCA){
	133	printf("DCA cut values for DATA not found in the output file, use dafault values. DCAcut=%f with default resolution\n",DCAcut);
	134	fDCAxyCutDATA = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
	135	for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutDATA->SetParameter(ipar,xyP[ipar]);
	136	fDCAxyCutDATA->SetParameter(3,DCAcut);
	137	nSigmaDCAxyDATA=DCAcut;
	138	}
982ff3fe	139
a1c5e4a1	140	if(TMath::Abs(nSigmaDCAxyDATA-nSigmaDCAxyMC)<0.01) DCAcut=nSigmaDCAxyDATA;
	141	else{
	142	printf("ERROR: DCAxy cuts do not match between data and MC\n");
	143	return;
	144	}
982ff3fe	145	if(TMath::Abs(nSigmaDCAzDATA-nSigmaDCAzMC)>0.01){
982ff3fe	146	printf("ERROR: DCAz cuts do not match between data (%f) and MC (%f) \n",nSigmaDCAzDATA,nSigmaDCAzMC);
a1c5e4a1	147	return;
a1c5e4a1	148	}
982ff3fe	149
a1c5e4a1	150	TCanvas *cfitDATA=new TCanvas("cfitDATA","cfitDATA");
	151	cfitDATA->Divide(6,4);
	152	TCanvas *cfitMC=new TCanvas("cfitMC","cfitMC");
	153	cfitMC->Divide(6,4);
	154
	155	for(Int_t ibin=firstbin;ibin<nptbins;ibin++){
	156
	157	Double_t xyMax =fDCAxyCutDATA->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;
	158	Double_t xyMaxMC =fDCAxyCutMC->Eval(TMath::Abs((ptbins[ibin+1]+ptbins[ibin])/2))/10000.;
	159
	160	TObjArray *mcTemplates = 0x0;
	161	if(partname=="PosP"){
	162	mcTemplates = new TObjArray(3); // MC histograms are put in this array
	163	mcTemplates->Add(fHistDCAPrim[ibin]);
	164	mcTemplates->Add(fHistDCAsecSt[ibin]);
	165	mcTemplates->Add(fHistDCAsec[ibin]);
	166	}else{
	167	mcTemplates = new TObjArray(2); // MC histograms are put in this array
	168	mcTemplates->Add(fHistDCAPrim[ibin]);
	169	mcTemplates->Add(fHistDCAsecSt[ibin]);
	170	}
982ff3fe	171	// if(fHistDCADATA[ibin]->Integral()==0) continue;
a1c5e4a1	172
	173	////////////////// Fit on DATA
	174	cfitDATA->cd(ibin);
	175	gPad->SetLogy();
	176	Double_t primAllDATA=PerformFit(fHistDCADATA[ibin],mcTemplates,partname,xyMax);
	177
	178	////////////////// Fit on MC
	179	cfitMC->cd(ibin);
	180	gPad->SetLogy();
	181	Double_t primAllMC=PerformFit(fHistDCAMC[ibin],mcTemplates,partname,xyMaxMC);
	182
	183	if(hPrimAllMC!=0 && hPrimAllDATA!=0){
	184	hPrimAllDATA->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA);
	185	hPrimAllMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllMC);
	186	hPrimAllDATAMC->Fill((ptbins[ibin]+ptbins[ibin+1])/2,primAllDATA/primAllMC);
	187	}
	188	}
	189
	190	//Adding MC truth
	191	TH1F* hAll=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaMean%d",pCharge.Data(),iSpecies));
	192	TH1F* hPrim=(TH1F*)liMC->FindObject(Form("hHist%sNSigmaPrimMean%d",pCharge.Data(),iSpecies));
	193	TH1F* hPrimRecMC=(TH1F*)liMC->FindObject(Form("fHistPrimMC%sReco%d",pcharge.Data(),iSpecies));
	194	TH1F* hSecMatRecMC=(TH1F*)liMC->FindObject(Form("fHistSecMatMC%sReco%d",pcharge.Data(),iSpecies));
	195	TH1F* hSecStrRecMC=(TH1F*)liMC->FindObject(Form("fHistSecStrMC%sReco%d",pcharge.Data(),iSpecies));
	196	hSecStrRecMC->Add(hSecMatRecMC);
	197	hSecStrRecMC->Add(hPrimRecMC);
	198	hPrim->SetLineColor(8);
	199	hPrim->SetLineWidth(2);
	200	hPrim->SetTitle("Prim/All from MC Truth");
	201	hPrim->Divide(hAll);
	202	hPrimRecMC->SetLineColor(11);
	203	hPrimRecMC->SetLineWidth(2);
	204	hPrimRecMC->SetTitle("Prim/All from MC Truth Reco");
	205	hPrimRecMC->Divide(hSecStrRecMC);
	206
	207	hPrimAllDATAMC->GetXaxis()->SetRangeUser(0.09,0.79);
	208	hPrimAllMC->GetXaxis()->SetRangeUser(0.09,0.79);
	209	hPrimAllDATA->GetXaxis()->SetRangeUser(0.09,0.79);
	210	if(optSmooth) hPrimAllDATAMC->Smooth(1,"R");
	211
	212	hPrimAllDATA->SetTitle("Prim/all data");
	213	hPrimAllDATA->SetLineColor(2);
	214	hPrimAllDATA->SetLineWidth(2);
	215	hPrimAllMC->SetTitle("Prim/all mc");
	216	hPrimAllMC->SetLineColor(4);
	217	hPrimAllMC->SetLineWidth(2);
	218	hPrimAllDATAMC->SetTitle("Prim/all data/mc");
	219	hPrimAllDATAMC->SetLineColor(1);
	220	hPrimAllDATAMC->SetLineWidth(2);
	221
	222	TCanvas *cPrimAll=new TCanvas("cPrimAll","cPrimAll");
	223	hPrimAllDATA->SetMinimum(0.7);
	224	hPrimAllDATA->SetMaximum(1.1);
	225	hPrimAllDATA->Draw("l");
	226	hPrimAllMC->Draw("lsame");
	227	hPrimAllDATAMC->Draw("lsame");
	228	hPrim->Draw("lsame");
	229	hPrimRecMC->Draw("lsame");
	230	gPad->BuildLegend();
	231	TLatex* tsp=new TLatex(0.4,0.75,speciesRoot[idPart].Data());
	232	tsp->SetTextFont(43);
	233	tsp->SetTextSize(28);
	234	tsp->SetNDC();
	235	tsp->Draw();
236	cPrimAll->Update();
237
238	TString fout;
239	fout=Form("DCACorr%s_%s_%.0fDCA_%s_TFraction.root",period.Data(),MCname.Data(),DCAcut,partname.Data());
240	TFile *out=new TFile(fout.Data(),"recreate");
241	hPrimAllDATA->Write();
242	hPrimAllMC->Write();
243	hPrimAllDATAMC->Write();
244	out->Close();
245	delete out;
246
247	}
248
249	Float_t PerformFit(TH1F* fHistDCA, TObjArray *mc,TString partname,Float_t xyMax){
250
251	Double_t prim=0,secSt=0,sec=0;
252	TFractionFitter* fit = new TFractionFitter(fHistDCA, mc); // initialise
253	fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1
254	fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1
255	if(partname=="PosP")fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1
256
257	//fit->SetRangeX(200,1000); // use only the first 15 bins in the fit
258	Int_t status = fit->Fit(); // perform the fit
259	cout << "fit status: " << status << endl;
260
261	if (status == 0) { // check on fit status
262	TH1F* result = (TH1F*) fit->GetPlot();
263	TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
264	TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
265	TH1F* secMCPred=0x0;
266	if(partname=="PosP"){
267	secMCPred=(TH1F*)fit->GetMCPrediction(2);
268	secMCPred->SetLineColor(4);
269	}
270	//Drawing section
271	PrimMCPred->SetLineColor(2);
272	secStMCPred->SetLineColor(6);
273	fHistDCA->SetTitle("DCA distr - data");
274	fHistDCA->DrawCopy("Ep");
275	result->SetTitle("Fit result");
276	result->DrawCopy("same");
277	Double_t value,error;
278
279	fit->GetResult(0,value,error);
280	PrimMCPred->Scale(fHistDCA->GetSumOfWeights()*value/PrimMCPred->GetSumOfWeights());
281	PrimMCPred->SetTitle("Primaries");
282	PrimMCPred->DrawCopy("same");
283	fit->GetResult(1,value,error);
284	secStMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secStMCPred->GetSumOfWeights());
285	secStMCPred->SetTitle("Sec from strangeness");
286	secStMCPred->DrawCopy("same");
287	if(partname=="PosP" && secMCPred){
288	fit->GetResult(2,value,error);
289	secMCPred->Scale(fHistDCA->GetSumOfWeights()*value/secMCPred->GetSumOfWeights());
290	secMCPred->SetTitle("Sec from material");
291	secMCPred->DrawCopy("same");
292	}
293	prim=PrimMCPred->Integral(PrimMCPred->FindBin(-xyMax),PrimMCPred->FindBin(xyMax));
294	secSt=secStMCPred->Integral(secStMCPred->FindBin(-xyMax),secStMCPred->FindBin(xyMax));
295	if(partname=="PosP")sec=secMCPred->Integral(secMCPred->FindBin(-xyMax),secMCPred->FindBin(xyMax));
296	}
297	else{
298	prim=1;
299	secSt=1;
300	if(partname=="PosP")sec=1;
301	}
302	Printf("Yields: primary=%f material=%f strange=%f\n",prim,sec,secSt);
303	return prim/(prim+secSt+sec);
304	}