[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / fitD0InvMass.C

#include <Riostream.h>
#include <TStyle.h>
#include <TFile.h>
#include <TList.h>
#include <TDirectoryFile.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TCanvas.h>

#include <AliHFMassFitter.h>

#include <fstream>
#include <cmath>

//root[0] .x fitD0InvMass.C+
//root[1] fitD0(...)

//input: nsigma is the number of sigma in which significance, signal and background are calculated

//fit of histograms in the directory PWG3_D2H_D0InvMass of AnalysisResults.root . Must specify the list name where the histo are stored. Produce a root file with the fitted histos and a text file with signal, background and significance

void fitD0(Int_t rebin=0,TString listname="coutputmassD0mycuts",Int_t nsigma=3, TString pathin="./",TString pathout="./",Int_t btype=2){

  TString file="AnalysisResults.root"; //"D0InvMass.root";

  file.Prepend(pathin);
  cout<<"Opening "<<file<<endl;
  TFile *fin=new TFile(file.Data());
  if(!fin->IsOpen()){
    cout<<"File "<<file.Data()<<" not found"<<endl;
    return;
  }

  gStyle->SetOptFit(0111);
  gStyle->SetOptStat("nemrou");

  const Int_t npt=5;
  TString dirname="PWG3_D2H_D0InvMass";
  TDirectoryFile *dir=(TDirectoryFile*)fin->GetDirectory(dirname);

  TList *lista = (TList*) dir->Get(listname.Data());
  if(!lista) {
    cout<<listname<<" doesn't exist"<<endl;
    return;
  }

  TString fileout="table.dat";
  fileout.Prepend(pathout);
  ofstream out(fileout.Data());

  AliHFMassFitter *fitter=new AliHFMassFitter();

  cout<<"mean = "<<fitter->GetMean()<<"\tsigma= "<<fitter->GetSigma()<<endl;
  Bool_t init=kTRUE;

  out<<"Bakground type = "<<btype<<";\trebin = "<<rebin<<endl;
  Int_t i=0;
  for(Int_t ipt=0;ipt<npt;ipt++){ //npt

    out<<"************************************************************\n";
    out<<"* ptbin * entriesMass * entriesSgn * entriesBkg * entriesRfl *\n";

    //ALL w cuts
    TString nameall="histMass_";
    nameall+=(ipt+1);
    TH1F *hMass=(TH1F*)lista->FindObject(nameall);
    if(!hMass){
      cout<<nameall<<" not found"<<endl;
      return;
    }

    out<<"* "<< ipt+1 <<"  *  "<<hMass->GetEntries()<<"  *  ";
    //SIGNAL from MC
    TString namesgn="histSgn_";
    namesgn+=(ipt+1);
    TH1F *hSgn=(TH1F*)lista->FindObject(namesgn);
    if(!hSgn){
      cout<<namesgn<<" not found"<<endl;
      return;
    }
    
    out<<hSgn->GetEntries()<<"  *  ";
    //BACKGROUND from MC
    TString namebkg="histBkg_";
    namebkg+=(ipt+1);
    TH1F *hBkg=(TH1F*)lista->FindObject(namebkg);
    if(!hBkg){
      cout<<namebkg<<" not found"<<endl;
      return;
    }

    out<<hBkg->GetEntries()<<"  *  ";
    //REFLECTED SIGNAL from MC
    TString namerfl="histRfl_";
    namerfl+=(ipt+1);
    TH1F *hRfl=(TH1F*)lista->FindObject(namerfl);
    if(!hRfl){
      cout<<namerfl<<" not found"<<endl;
      return;
    }

    out<<hRfl->GetEntries()<<"  *\n";
    out<<"*************************************************\n";
    if(rebin!=0){
      hSgn->Rebin(rebin);
      hBkg->Rebin(rebin);
      hRfl->Rebin(rebin);
    }
    Double_t min, max;
    Int_t    nbin;
    Double_t sgn,errsgn,errsgn2=0,bkg,errbkg,errbkg2=0,sgnf,errsgnf,sgnfMC;
    Double_t meanfrfit,sigmafrfit, meanMC=hSgn->GetMean(),sigmaMC=hSgn->GetRMS();
    Double_t width=0;
    nbin=hMass->GetNbinsX();
    min=hMass->GetBinLowEdge(1);
    max=min+nbin*hMass->GetBinWidth(nbin);

    //FIT:

    TString namentu="ntupt3bin";
    if(init) fitter->InitNtuParam((char*)namentu.Data());
    // - all
    fitter->SetHisto(hMass);
    fitter->SetRangeFit(min, max);
    //fitter->SetRangeFit(1.83,1.89);
    fitter->SetType(btype,0);//(b,s)
    if(rebin!=0) fitter->RebinMass(rebin);
    Bool_t fitOK=fitter->MassFitter(kFALSE); //kFALSE = do not draw
    if(!fitOK) {
      out<<"Fit return kFALSE, skip "<<hMass->GetName()<<endl;
      fitter->Reset(); //delete histogram set
      continue;
    }
    width=fitter->GetHistoClone()->GetBinWidth(3);
    cout<<"\nChi^2 = "<<fitter->GetChiSquare()<<"\t Reduced Chi^2 = "<<fitter->GetReducedChiSquare()<<endl;
    meanfrfit=fitter->GetMean();
    sigmafrfit=fitter->GetSigma();
    out<<"mean = "<<meanfrfit<<" (MC "<<meanMC<<")"<<"\tsigma= "<<sigmafrfit<<" (MC "<<sigmaMC<<")"<<endl;
    //cout<<"old nbin = "<<hMass->GetNbinsX()<<"\tnew nbin = "<<fitter->GetHistoClone()->GetNbinsX()<<endl;
    if(meanfrfit<0 || sigmafrfit<0 || meanfrfit<min || meanfrfit>max) {
      cout<<"Fit failed, check"<<endl;
      out<<hMass->GetName();
      out<<": \nSgn not available"<<endl;
      out<<"Bkg  not available"<<endl;
      out<<"Sgnf not available"<<endl;
      out<<"*************************************************\n";
    }
    else{
      //cout<<"Writing..."<<endl;
      fitter->WriteHisto(pathout);
      hMass= fitter->GetHistoClone();
      //TH1F *hc=fitter->GetHistoClone();
      //TF1 *fbtest=hc->GetFunction("funcbkgRecalc"); //new version of ALiHFMassFitter
     
      fitter->FillNtuParam();
     
      Double_t limsx,limdx;
      limsx=meanfrfit-nsigma*sigmafrfit;
      limdx=meanfrfit+nsigma*sigmafrfit;
      // limsx=meanMC-nsigma*sigmaMC;
      //limdx=meanMC+nsigma*sigmaMC;


      //determine limit of nsigma in bins
      Int_t binsx,bindx;
      binsx=hMass->FindBin(limsx);
      if (limsx > hMass->GetBinCenter(binsx)) binsx++;
      bindx=hMass->FindBin(limdx);
      if (limdx < hMass->GetBinCenter(bindx)) bindx--;

      //reconvert bin in x
      Double_t sxr,dxr;
      sxr=hMass->GetBinLowEdge(binsx);
      dxr=hMass->GetBinLowEdge(bindx+1);

      fitter->Signal(sxr,dxr,sgn,errsgn);
      fitter->Background(sxr,dxr,bkg,errbkg);
      fitter->Significance(sxr,dxr,sgnf,errsgnf);
      

      Float_t inttot,intsgn,intsgnerr;

      Int_t np=-99;
      switch (btype){
      case 0: //expo
	np=2;
	break;
      case 1: //linear
	np=2;
	break;
      case 2: //pol2
	np=3;
	break;
      case 3: //no bkg
	np=1;
	break;
      }

      
      TF1 *fmass=hMass->GetFunction("funcmass");
      if (fmass){
	
	inttot=fmass->GetParameter(0);
	intsgn=fmass->GetParameter(np);
	intsgnerr=fmass->GetParError(np);
	//cout<<"i = "<<i<<"inttot = "<<inttot<<"\tintsgn = "<<intsgn<<"\tintsgnErr = "<<intsgnerr<<"\twidth = "<<width<<endl;
	Double_t errbkg2rel=errbkg/bkg*TMath::Sqrt((inttot-intsgn)/width/bkg);//intsgnerr/(inttot-intsgn)*TMath::Sqrt((inttot-intsgn)/width/bkg);
	Double_t errsgn2rel=errsgn/sgn*TMath::Sqrt(intsgn/width/sgn);
	errbkg2=errbkg2rel*bkg;
	errsgn2=errsgn2rel*sgn;
      }

           
      cout<<"bin sx = "<<binsx<<"\t xsx = "<<sxr<<endl;
      cout<<"bin dx = "<<bindx<<"\t xdx = "<<dxr<<endl;

      out<<hMass->GetName();
      Double_t sgnMC,bkgMC,rflMC;
      sgnMC=hSgn->Integral(binsx,bindx);
      bkgMC=hBkg->Integral(binsx,bindx);
      rflMC=hRfl->Integral(binsx,bindx);
      sgnfMC=sgnMC/TMath::Sqrt(sgnMC+bkgMC+rflMC);
      out<<": \nSgn  "<<sgn<<"+/-"<<errsgn<<" ("<<errsgn2<<")\tCompare with: "<<sgnMC<<endl;
      out<<"Bkg  "<<bkg<<"+/-"<<errbkg<<" ("<<errbkg2<<")\tCompare with: "<<bkgMC<<" + "<<rflMC<<" = "<<bkgMC+rflMC<<endl;
      out<<"Sgnf "<<sgnf<<"+/-"<<errsgnf<<"\tCompare with: "<<sgnfMC<<endl;
      if (sgn != 0 && sgnfMC != 0) out<<"sigmaS/S = "<<errsgn/sgn<<"\tCompare with 1/signif: "<<1./sgnfMC<<endl;
      else {
	out<<"sigmaS/S = "<<errsgn<<"/"<<sgn<<"\tCompare with 1/signif: 1./"<<sgnfMC<<endl;
      }
      out<<"nsigma considered for comparison = \ndx: "<<(dxr-meanfrfit)/sigmafrfit<<"\nsx: "<<(meanfrfit-sxr)/sigmafrfit<<endl;
      out<<"Mean = "<<meanfrfit<<"\tSigma = "<<sigmafrfit<<endl;
      out<<"*************************************************\n";
      i++;
    }
    
    sgn=0; bkg=0; sgnf=0;
    errsgn=0; errbkg=0; errsgnf=0;
    if (ipt == npt-1) fitter->WriteNtuple(pathout);

    out<<endl;
    fitter->Reset(); //delete histogram set

    init=kFALSE;

  }

  out.close();
}
Commit	Line	Data
ed68e1ea	1	#include <Riostream.h>
	2	#include <TStyle.h>
	3	#include <TFile.h>
	4	#include <TList.h>
33924cbd	5	#include <TDirectoryFile.h>
ed68e1ea	6	#include <TClonesArray.h>
	7	#include <TH1F.h>
	8	#include <TCanvas.h>
	9
	10	#include <AliHFMassFitter.h>
	11
	12	#include <fstream>
	13	#include <cmath>
	14
	15	//root[0] .x fitD0InvMass.C+
33924cbd	16	//root[1] fitD0(...)
ed68e1ea	17
33924cbd	18	//input: nsigma is the number of sigma in which significance, signal and background are calculated
ed68e1ea	19
33924cbd	20	//fit of histograms in the directory PWG3_D2H_D0InvMass of AnalysisResults.root . Must specify the list name where the histo are stored. Produce a root file with the fitted histos and a text file with signal, background and significance
ed68e1ea	21
33924cbd	22	void fitD0(Int_t rebin=0,TString listname="coutputmassD0mycuts",Int_t nsigma=3, TString pathin="./",TString pathout="./",Int_t btype=2){
	23
	24	TString file="AnalysisResults.root"; //"D0InvMass.root";
ed68e1ea	25
ed68e1ea	26	file.Prepend(pathin);
	27	cout<<"Opening "<<file<<endl;
	28	TFile *fin=new TFile(file.Data());
	29	if(!fin->IsOpen()){
33924cbd	30	cout<<"File "<<file.Data()<<" not found"<<endl;
33924cbd	31	return;
ed68e1ea	32	}
ed68e1ea	33
33924cbd	34	gStyle->SetOptFit(0111);
	35	gStyle->SetOptStat("nemrou");
	36
	37	const Int_t npt=5;
	38	TString dirname="PWG3_D2H_D0InvMass";
	39	TDirectoryFile dir=(TDirectoryFile)fin->GetDirectory(dirname);
	40
	41	TList lista = (TList) dir->Get(listname.Data());
ed68e1ea	42	if(!lista) {
	43	cout<<listname<<" doesn't exist"<<endl;
	44	return;
	45	}
	46
	47	TString fileout="table.dat";
	48	fileout.Prepend(pathout);
	49	ofstream out(fileout.Data());
	50
	51	AliHFMassFitter *fitter=new AliHFMassFitter();
	52
	53	cout<<"mean = "<<fitter->GetMean()<<"\tsigma= "<<fitter->GetSigma()<<endl;
	54	Bool_t init=kTRUE;
	55
	56	out<<"Bakground type = "<<btype<<";\trebin = "<<rebin<<endl;
	57	Int_t i=0;
	58	for(Int_t ipt=0;ipt<npt;ipt++){ //npt
	59
	60	out<<"************************************************************\n";
	61	out<<"* ptbin * entriesMass * entriesSgn * entriesBkg * entriesRfl *\n";
	62
	63	//ALL w cuts
	64	TString nameall="histMass_";
	65	nameall+=(ipt+1);
	66	TH1F hMass=(TH1F)lista->FindObject(nameall);
	67	if(!hMass){
	68	cout<<nameall<<" not found"<<endl;
	69	return;
	70	}
	71
	72	out<<"* "<< ipt+1 <<" * "<<hMass->GetEntries()<<" * ";
	73	//SIGNAL from MC
	74	TString namesgn="histSgn_";
	75	namesgn+=(ipt+1);
	76	TH1F hSgn=(TH1F)lista->FindObject(namesgn);
	77	if(!hSgn){
	78	cout<<namesgn<<" not found"<<endl;
	79	return;
	80	}
	81
	82	out<<hSgn->GetEntries()<<" * ";
	83	//BACKGROUND from MC
	84	TString namebkg="histBkg_";
	85	namebkg+=(ipt+1);
	86	TH1F hBkg=(TH1F)lista->FindObject(namebkg);
	87	if(!hBkg){
	88	cout<<namebkg<<" not found"<<endl;
	89	return;
	90	}
	91
	92	out<<hBkg->GetEntries()<<" * ";
	93	//REFLECTED SIGNAL from MC
	94	TString namerfl="histRfl_";
	95	namerfl+=(ipt+1);
	96	TH1F hRfl=(TH1F)lista->FindObject(namerfl);
	97	if(!hRfl){
	98	cout<<namerfl<<" not found"<<endl;
	99	return;
	100	}
	101
	102	out<<hRfl->GetEntries()<<" *\n";
	103	out<<"*************************************************\n";
	104	if(rebin!=0){
	105	hSgn->Rebin(rebin);
106	hBkg->Rebin(rebin);
107	hRfl->Rebin(rebin);
108	}
109	Double_t min, max;
110	Int_t nbin;
111	Double_t sgn,errsgn,errsgn2=0,bkg,errbkg,errbkg2=0,sgnf,errsgnf,sgnfMC;
112	Double_t meanfrfit,sigmafrfit, meanMC=hSgn->GetMean(),sigmaMC=hSgn->GetRMS();
113	Double_t width=0;
114	nbin=hMass->GetNbinsX();
115	min=hMass->GetBinLowEdge(1);
116	max=min+nbin*hMass->GetBinWidth(nbin);
117
118	//FIT:
119
120	TString namentu="ntupt3bin";
121	if(init) fitter->InitNtuParam((char*)namentu.Data());
122	// - all
123	fitter->SetHisto(hMass);
124	fitter->SetRangeFit(min, max);
125	//fitter->SetRangeFit(1.83,1.89);
126	fitter->SetType(btype,0);//(b,s)
127	if(rebin!=0) fitter->RebinMass(rebin);
128	Bool_t fitOK=fitter->MassFitter(kFALSE); //kFALSE = do not draw
129	if(!fitOK) {
130	out<<"Fit return kFALSE, skip "<<hMass->GetName()<<endl;
131	fitter->Reset(); //delete histogram set
132	continue;
133	}
134	width=fitter->GetHistoClone()->GetBinWidth(3);
135	cout<<"\nChi^2 = "<<fitter->GetChiSquare()<<"\t Reduced Chi^2 = "<<fitter->GetReducedChiSquare()<<endl;
136	meanfrfit=fitter->GetMean();
137	sigmafrfit=fitter->GetSigma();
138	out<<"mean = "<<meanfrfit<<" (MC "<<meanMC<<")"<<"\tsigma= "<<sigmafrfit<<" (MC "<<sigmaMC<<")"<<endl;
139	//cout<<"old nbin = "<<hMass->GetNbinsX()<<"\tnew nbin = "<<fitter->GetHistoClone()->GetNbinsX()<<endl;
140	if(meanfrfit<0 \|\| sigmafrfit<0 \|\| meanfrfit<min \|\| meanfrfit>max) {
141	cout<<"Fit failed, check"<<endl;
142	out<<hMass->GetName();
143	out<<": \nSgn not available"<<endl;
144	out<<"Bkg not available"<<endl;
145	out<<"Sgnf not available"<<endl;
146	out<<"*************************************************\n";
147	}
148	else{
149	//cout<<"Writing..."<<endl;
150	fitter->WriteHisto(pathout);
151	hMass= fitter->GetHistoClone();
152	//TH1F *hc=fitter->GetHistoClone();
153	//TF1 *fbtest=hc->GetFunction("funcbkgRecalc"); //new version of ALiHFMassFitter
154
155	fitter->FillNtuParam();
156
157	Double_t limsx,limdx;
158	limsx=meanfrfit-nsigma*sigmafrfit;
159	limdx=meanfrfit+nsigma*sigmafrfit;
160	// limsx=meanMC-nsigma*sigmaMC;
161	//limdx=meanMC+nsigma*sigmaMC;
162
163
164	//determine limit of nsigma in bins
165	Int_t binsx,bindx;
166	binsx=hMass->FindBin(limsx);
167	if (limsx > hMass->GetBinCenter(binsx)) binsx++;
168	bindx=hMass->FindBin(limdx);
169	if (limdx < hMass->GetBinCenter(bindx)) bindx--;
170
171	//reconvert bin in x
172	Double_t sxr,dxr;
173	sxr=hMass->GetBinLowEdge(binsx);
174	dxr=hMass->GetBinLowEdge(bindx+1);
175
176	fitter->Signal(sxr,dxr,sgn,errsgn);
177	fitter->Background(sxr,dxr,bkg,errbkg);
178	fitter->Significance(sxr,dxr,sgnf,errsgnf);
179
180
181	Float_t inttot,intsgn,intsgnerr;
182
183	Int_t np=-99;
184	switch (btype){
185	case 0: //expo
186	np=2;
187	break;
188	case 1: //linear
189	np=2;
190	break;
191	case 2: //pol2
192	np=3;
193	break;
194	case 3: //no bkg
195	np=1;
196	break;
197	}
198
199
200	TF1 *fmass=hMass->GetFunction("funcmass");
201	if (fmass){
202
203	inttot=fmass->GetParameter(0);
204	intsgn=fmass->GetParameter(np);
205	intsgnerr=fmass->GetParError(np);
206	//cout<<"i = "<<i<<"inttot = "<<inttot<<"\tintsgn = "<<intsgn<<"\tintsgnErr = "<<intsgnerr<<"\twidth = "<<width<<endl;
207	Double_t errbkg2rel=errbkg/bkgTMath::Sqrt((inttot-intsgn)/width/bkg);//intsgnerr/(inttot-intsgn)TMath::Sqrt((inttot-intsgn)/width/bkg);
208	Double_t errsgn2rel=errsgn/sgn*TMath::Sqrt(intsgn/width/sgn);
209	errbkg2=errbkg2rel*bkg;
210	errsgn2=errsgn2rel*sgn;
211	}
212
213
214	cout<<"bin sx = "<<binsx<<"\t xsx = "<<sxr<<endl;
215	cout<<"bin dx = "<<bindx<<"\t xdx = "<<dxr<<endl;
216
217	out<<hMass->GetName();
218	Double_t sgnMC,bkgMC,rflMC;
219	sgnMC=hSgn->Integral(binsx,bindx);
220	bkgMC=hBkg->Integral(binsx,bindx);
221	rflMC=hRfl->Integral(binsx,bindx);
222	sgnfMC=sgnMC/TMath::Sqrt(sgnMC+bkgMC+rflMC);
223	out<<": \nSgn "<<sgn<<"+/-"<<errsgn<<" ("<<errsgn2<<")\tCompare with: "<<sgnMC<<endl;
224	out<<"Bkg "<<bkg<<"+/-"<<errbkg<<" ("<<errbkg2<<")\tCompare with: "<<bkgMC<<" + "<<rflMC<<" = "<<bkgMC+rflMC<<endl;
225	out<<"Sgnf "<<sgnf<<"+/-"<<errsgnf<<"\tCompare with: "<<sgnfMC<<endl;
33924cbd	226	if (sgn != 0 && sgnfMC != 0) out<<"sigmaS/S = "<<errsgn/sgn<<"\tCompare with 1/signif: "<<1./sgnfMC<<endl;
	227	else {
	228	out<<"sigmaS/S = "<<errsgn<<"/"<<sgn<<"\tCompare with 1/signif: 1./"<<sgnfMC<<endl;
	229	}
ed68e1ea	230	out<<"nsigma considered for comparison = \ndx: "<<(dxr-meanfrfit)/sigmafrfit<<"\nsx: "<<(meanfrfit-sxr)/sigmafrfit<<endl;
	231	out<<"Mean = "<<meanfrfit<<"\tSigma = "<<sigmafrfit<<endl;
	232	out<<"*************************************************\n";
	233	i++;
	234	}
	235
	236	sgn=0; bkg=0; sgnf=0;
	237	errsgn=0; errbkg=0; errsgnf=0;
	238	if (ipt == npt-1) fitter->WriteNtuple(pathout);
	239
	240	out<<endl;
	241	fitter->Reset(); //delete histogram set
	242
	243	init=kFALSE;
	244
	245	}
	246
	247	out.close();
	248	}