]>
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; |
31 | return; | |
ed68e1ea | 32 | } |
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/bkg*TMath::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 | } |