]>
Commit | Line | Data |
---|---|---|
6bb658a2 | 1 | //gROOT->Reset(); |
2 | #include <math.h> | |
3 | #include <time.h> | |
4 | #include <stdio.h> | |
5 | #include <stdlib.h> | |
6 | #include <Riostream.h> | |
7 | #include "RooFit.h" | |
8 | #include "RooVoigtian.h" | |
9 | #include "RooAbsReal.h" | |
10 | #include "RooRealVar.h" | |
11 | #include "RooComplex.h" | |
12 | #include "RooMath.h" | |
13 | ||
14 | #include "TVector2.h" | |
15 | #include "TFile.h" | |
16 | #include "TString.h" | |
17 | #include "TF1.h" | |
18 | #include "TH1.h" | |
19 | #include "TH2.h" | |
20 | #include "TH3.h" | |
21 | #include "TProfile.h" | |
22 | #include "TProfile2D.h" | |
23 | #include "TMath.h" | |
24 | #include "TText.h" | |
25 | #include "TRandom3.h" | |
26 | #include "TArray.h" | |
27 | #include "TLegend.h" | |
28 | #include "TStyle.h" | |
29 | #include "TMinuit.h" | |
30 | ||
31 | #include "TCanvas.h" | |
32 | #include "TLatex.h" | |
33 | #include "TImage.h" | |
34 | #include "TGaxis.h" | |
35 | #include "TGraphErrors.h" | |
36 | #include "TPaveStats.h" | |
37 | ||
38 | using namespace std; | |
39 | ||
40 | ||
41 | double myLevyPtFunc(Double_t *x, Double_t *par); | |
42 | double myMtExpFunc(Double_t *x, Double_t *par); | |
43 | double myPtExpFunc(Double_t *x, Double_t *par); | |
44 | double myBoltzFunc(Double_t *x, Double_t *par); | |
45 | double myPowerLawPtFunc(Double_t *x, Double_t *par); | |
46 | double myPowerLawFunc(Double_t *x, Double_t *par); | |
47 | double myBoltsBWFunc(Double_t *x, Double_t *par); | |
48 | double IntegrandBG(const double *x, const double *p); | |
49 | ||
50 | void Plot_XiStarResults(){ | |
51 | ||
52 | TGaxis::SetMaxDigits(3); | |
53 | gStyle->SetOptStat(0); | |
54 | gStyle->SetOptFit(01111); | |
55 | ||
56 | ||
57 | TFile *filesCutVar[13]; | |
58 | TH1D *Spectra_CutVar[13]; | |
59 | TF1 *LevyFits_CutVar[13]; | |
60 | TH1D *Yields_CutVar[13]; | |
61 | TH1D *Eff_CutVar[13]; | |
62 | for(int i=0; i<13; i++){ | |
63 | TString *name=new TString("CutVariation_"); | |
64 | *name += i; | |
65 | name->Append(".root"); | |
66 | filesCutVar[i]=new TFile(name->Data()); | |
67 | Spectra_CutVar[i] = (TH1D*)filesCutVar[i]->Get("h_Spectrum"); | |
68 | Spectra_CutVar[i]->SetDirectory(0); | |
69 | LevyFits_CutVar[i] = (TF1*)filesCutVar[i]->Get("myLevy"); | |
70 | Yields_CutVar[i] = (TH1D*)filesCutVar[i]->Get("RawYields"); | |
71 | Eff_CutVar[i] = (TH1D*)filesCutVar[i]->Get("Efficiency"); | |
72 | } | |
73 | // | |
74 | TFile *filesFitVar[4]; | |
75 | TH1D *Spectra_FitVar[4]; | |
76 | TF1 *LevyFits_FitVar[4]; | |
77 | TH1D *Yields_FitVar[4]; | |
78 | for(int i=0; i<4; i++){ | |
79 | TString *name=new TString("FitVariation_"); | |
80 | *name += i+1; | |
81 | name->Append(".root"); | |
82 | filesFitVar[i]=new TFile(name->Data()); | |
83 | Spectra_FitVar[i] = (TH1D*)filesFitVar[i]->Get("h_Spectrum"); | |
84 | Spectra_FitVar[i]->SetDirectory(0); | |
85 | LevyFits_FitVar[i] = (TF1*)filesFitVar[i]->Get("myLevy"); | |
86 | Yields_FitVar[i] = (TH1D*)filesFitVar[i]->Get("RawYields"); | |
87 | } | |
88 | ||
89 | double FVavg=0; | |
90 | for(int bin=2; bin<=9; bin++){ | |
91 | // FV 1 has largest spectrum point variation | |
92 | FVavg += fabs(Spectra_FitVar[0]->GetBinContent(bin)-Spectra_FitVar[1]->GetBinContent(bin))/Spectra_FitVar[0]->GetBinContent(bin); | |
93 | } | |
94 | cout<<"Average point-by-point spectrum variation due to fitting: "<<FVavg/8<<endl; | |
95 | ||
96 | double TrackVavg=0; | |
97 | for(int bin=2; bin<=9; bin++){ | |
98 | // CutVar 5 has largest spectrum point variation in the "tracking systematic" category | |
99 | TrackVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[5]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin); | |
100 | } | |
101 | cout<<"Average point-by-point spectrum variation due to Tracking: "<<TrackVavg/8<<endl; | |
102 | ||
103 | double TopologVavg=0; | |
104 | for(int bin=2; bin<=9; bin++){ | |
105 | // CutVar 8 has largest spectrum point variation in the "topological systematic" category | |
106 | TopologVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[8]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin); | |
107 | } | |
108 | cout<<"Average point-by-point spectrum variation due to Topological Cuts: "<<TopologVavg/8<<endl; | |
109 | ||
110 | //////////////////////////////////// | |
111 | TH1D *FinalSpectrum = (TH1D*)Spectra_CutVar[0]->Clone(); | |
112 | FinalSpectrum->SetMarkerSize(1.); | |
113 | double StatError[8]={0}; | |
114 | for(int ptbin=1; ptbin<=8; ptbin++) StatError[ptbin-1] = FinalSpectrum->GetBinError(ptbin+1);// add 1 since first bin is a dummy | |
115 | ||
116 | ||
117 | TH1D *PureSysFromCutVar[12]; | |
118 | double MeanSysCutVar[12]={0}; | |
119 | double sigmaSysCutVar[12]={0}; | |
120 | for(int i=0; i<12; i++) { | |
121 | TString *name=new TString("PureSysFromCutVar"); | |
122 | *name += i+1; | |
123 | PureSysFromCutVar[i] = new TH1D(name->Data(),"",8,0.5,8.5); | |
124 | PureSysFromCutVar[i]->GetXaxis()->SetTitle("pt bin"); | |
125 | PureSysFromCutVar[i]->GetYaxis()->SetTitle("|y_{s}-y_{v}|/y_{s}"); | |
126 | } | |
127 | ||
128 | ////////////////////////////////////// | |
129 | // Sys error assignment for pt bins to be used for mean value fits | |
130 | ||
131 | int CutVarLow=1, CutVarHigh=13; | |
132 | for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13 | |
133 | for(int ptbin=2; ptbin<=9; ptbin++){ | |
134 | double SysVar = fabs(FinalSpectrum->GetBinContent(ptbin) - Spectra_CutVar[i]->GetBinContent(ptbin))/FinalSpectrum->GetBinContent(ptbin); | |
135 | double sigma = sqrt(fabs(pow(FinalSpectrum->GetBinError(ptbin),2) - pow(Spectra_CutVar[i]->GetBinError(ptbin),2)));// stat part | |
136 | PureSysFromCutVar[i-1]->SetBinContent(ptbin-1, SysVar); | |
137 | PureSysFromCutVar[i-1]->SetBinError(ptbin-1, sigma/FinalSpectrum->GetBinContent(ptbin)); | |
138 | } | |
139 | PureSysFromCutVar[i-1]->Fit("pol0","IMEQ","",0,9); | |
140 | MeanSysCutVar[i-1] = pol0->GetParameter(0); | |
141 | sigmaSysCutVar[i-1] = pol0->GetParError(0); | |
142 | cout<<"Cut Var "<<i<<" Mean Sys = "<<MeanSysCutVar[i-1]<<" += "<<pol0->GetParError(0)<<endl; | |
143 | //cout<<"Chi2/NDF = "<<pol0->GetChisquare()/pol0->GetNDF()<<endl; | |
144 | //TPaveStats *sb2 = (TPaveStats*)(PureSysFromCutVar[i-1]->GetListOfFunctions()->FindObject("stats")); | |
145 | //sb2->SetX1NDC(.2); | |
146 | } | |
147 | ||
148 | ||
149 | double SysError[9]={0}; | |
150 | // Cut Variations | |
151 | //int CutVarN=12; | |
152 | for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13 | |
153 | for(int ptbin=1; ptbin<=9; ptbin++){ | |
154 | if(MeanSysCutVar[i-1] > 0 && MeanSysCutVar[i-1]>sigmaSysCutVar[i-1]) { | |
155 | SysError[ptbin-1] += pow(MeanSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2) - pow(sigmaSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2); | |
156 | } | |
157 | } | |
158 | } | |
159 | // Fit Variations | |
160 | for(int i=2; i<3; i++){// now only use i=2 (max deviation). was 0 to 4 | |
161 | for(int ptbin=1; ptbin<=9; ptbin++){ | |
162 | SysError[ptbin-1] += pow(FinalSpectrum->GetBinContent(ptbin) - Spectra_FitVar[i]->GetBinContent(ptbin),2); | |
163 | } | |
164 | } | |
165 | ||
166 | // Add Stat errors and set Final Error | |
167 | for(int ptbin=1; ptbin<=9; ptbin++){ | |
168 | double TotError = pow(FinalSpectrum->GetBinError(ptbin),2);// Stat | |
169 | TotError += SysError[ptbin-1];// Add on Sys | |
170 | FinalSpectrum->SetBinError(ptbin, sqrt(TotError)); | |
171 | cout<<"Spectrum ptbin:"<<ptbin-1<<" "<<FinalSpectrum->GetBinContent(ptbin)<<" +- "<<FinalSpectrum->GetBinError(ptbin)<<endl; | |
172 | } | |
173 | ||
174 | for(int ptbin=1; ptbin<=9; ptbin++) SysError[ptbin-1] = sqrt(SysError[ptbin-1]); | |
175 | ||
176 | // print efficiencies | |
177 | for(int ptbin=1; ptbin<=8; ptbin++){ | |
178 | cout<<"Efficiency ptbin:"<<ptbin<<" "<<Eff_CutVar[0]->GetBinContent(ptbin)<<" +- "<<Eff_CutVar[0]->GetBinError(ptbin)<<endl; | |
179 | } | |
180 | ||
181 | TCanvas *can = new TCanvas("can","can",13,34,700,500); | |
182 | can->Range(-1.25,-0.2625,11.25,2.3625); | |
183 | can->SetFillColor(10); | |
184 | can->SetBorderMode(0); | |
185 | can->SetBorderSize(2); | |
186 | //can->SetGridx(); | |
187 | //can->SetGridy(); | |
188 | can->SetFrameFillColor(0); | |
189 | can->SetFrameBorderMode(0); | |
190 | can->SetFrameBorderMode(0); | |
191 | ||
192 | TLegend *legend = new TLegend(.1,.7,.35,.9,NULL,"brNDC"); | |
193 | legend->SetBorderSize(1); | |
194 | legend->SetTextSize(.04);// small .03; large .036 | |
195 | //legend->SetLineColor(0); | |
196 | legend->SetFillColor(0); | |
197 | TLegend *legend2 = new TLegend(.35,.62,.9,.72,NULL,"brNDC"); | |
198 | legend2->SetBorderSize(1); | |
199 | legend2->SetTextSize(.04);// small .03; large .036 | |
200 | //legend2->SetLineColor(0); | |
201 | legend2->SetFillColor(0); | |
202 | ||
203 | ||
204 | double minFitpoint = 0.8; double maxFitpoint = 5.6; | |
205 | ||
206 | ||
207 | ||
208 | TF1 *myLevy=new TF1("myLevy",myLevyPtFunc,0,8.5,3); | |
209 | myLevy->SetParName(0,"dN/dy"); | |
210 | myLevy->SetParName(1,"C"); | |
211 | myLevy->SetParName(2,"n"); | |
212 | myLevy->SetParameter(0,.008); | |
213 | myLevy->SetParameter(1,.3); | |
214 | myLevy->SetParameter(2,15); | |
215 | myLevy->SetParLimits(0,.001,.2); | |
216 | myLevy->SetParLimits(1,.1,1); | |
217 | myLevy->SetParLimits(2,1,500); | |
218 | myLevy->SetLineColor(1); | |
219 | TF1 *myPtExp = new TF1("myPtExp",myPtExpFunc,0,8.5,2); | |
220 | myPtExp->SetParName(0,"dN/dy"); | |
221 | myPtExp->SetParName(1,"C"); | |
222 | myPtExp->SetParameter(0,.003); | |
223 | myPtExp->SetParameter(1,.3); | |
224 | myPtExp->SetParLimits(0,0.0001,0.1); | |
225 | myPtExp->SetParLimits(1,.1,1); | |
226 | myPtExp->SetLineColor(2); | |
227 | TF1 *myMtExp = new TF1("myMtExp",myMtExpFunc,0,8.5,2); | |
228 | myMtExp->SetParName(0,"dN/dy"); | |
229 | myMtExp->SetParName(1,"C"); | |
230 | myMtExp->SetParameter(0,.003); | |
231 | myMtExp->SetParameter(1,.3); | |
232 | myMtExp->SetParLimits(0,0.0001,0.1); | |
233 | myMtExp->SetParLimits(1,.1,1); | |
234 | myMtExp->SetLineColor(4); | |
235 | TF1 *myBoltz = new TF1("myBoltz",myBoltzFunc,0,8.5,2); | |
236 | myBoltz->SetParName(0,"dN/dy"); | |
237 | myBoltz->SetParName(1,"C"); | |
238 | myBoltz->SetParameter(0,.003); | |
239 | myBoltz->SetParameter(1,.3); | |
240 | myBoltz->SetParLimits(0,0.0001,0.1); | |
241 | myBoltz->SetParLimits(1,.1,1); | |
242 | myBoltz->SetLineColor(6); | |
243 | TF1 *myPowerLawPt=new TF1("myPowerLawPt",myPowerLawPtFunc,0,8.5,3); | |
244 | myPowerLawPt->SetParName(0,"dN/dy"); | |
245 | myPowerLawPt->SetParName(1,"pt0"); | |
246 | myPowerLawPt->SetParName(2,"n"); | |
247 | myPowerLawPt->SetParameter(0,.008); | |
248 | myPowerLawPt->SetParameter(1,50); | |
249 | myPowerLawPt->SetParameter(2,50); | |
250 | myPowerLawPt->SetParLimits(0,.001,.01); | |
251 | myPowerLawPt->SetParLimits(1,1,500); | |
252 | myPowerLawPt->SetParLimits(2,1,500); | |
253 | myPowerLawPt->SetLineColor(7); | |
254 | TF1 *myPowerLaw=new TF1("myPowerLaw",myPowerLawFunc,0,8.5,3); | |
255 | myPowerLaw->SetParName(0,"dN/dy"); | |
256 | myPowerLaw->SetParName(1,"pt0"); | |
257 | myPowerLaw->SetParName(2,"n"); | |
258 | myPowerLaw->SetParameter(0,.008); | |
259 | myPowerLaw->SetParameter(1,50); | |
260 | myPowerLaw->SetParameter(2,50); | |
261 | myPowerLaw->SetParLimits(0,.001,.01); | |
262 | myPowerLaw->SetParLimits(1,1,500); | |
263 | myPowerLaw->SetParLimits(2,1,500); | |
264 | myPowerLaw->SetLineColor(8); | |
265 | TF1 *myBGBlastWave=new TF1("myBGBlastWave",myBoltsBWFunc,0,8.5,5); | |
266 | myBGBlastWave->SetParName(0,"mass"); | |
267 | myBGBlastWave->SetParName(1,"beta"); | |
268 | myBGBlastWave->SetParName(2,"C"); | |
269 | myBGBlastWave->SetParName(3,"n"); | |
270 | myBGBlastWave->SetParName(4,"Norm"); | |
271 | myBGBlastWave->FixParameter(0,1.5318); | |
272 | myBGBlastWave->SetParameter(1,.3); | |
273 | myBGBlastWave->SetParameter(2,.3); | |
274 | myBGBlastWave->SetParameter(3,5); | |
275 | myBGBlastWave->SetParameter(4,.3); | |
276 | myBGBlastWave->SetParLimits(1,.001,.9); | |
277 | myBGBlastWave->SetParLimits(2,.1,200); | |
278 | myBGBlastWave->SetParLimits(3,.1,200); | |
279 | myBGBlastWave->SetLineColor(9); | |
280 | ||
281 | ||
282 | FinalSpectrum->Draw(); | |
283 | legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2","p"); | |
284 | //legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2: Current Results","p"); | |
285 | FinalSpectrum->Fit(myLevy,"IME","",minFitpoint,maxFitpoint); | |
286 | FinalSpectrum->Fit(myPtExp,"IMEN","",minFitpoint,maxFitpoint); | |
287 | FinalSpectrum->Fit(myMtExp,"IMEN","",minFitpoint,maxFitpoint); | |
288 | FinalSpectrum->Fit(myBoltz,"IMEN","",minFitpoint,maxFitpoint); | |
289 | FinalSpectrum->Fit(myPowerLaw,"IMEN","",minFitpoint,maxFitpoint); | |
290 | FinalSpectrum->Fit(myPowerLawPt,"IMEN","",minFitpoint,maxFitpoint); | |
291 | //FinalSpectrum->Fit(myBGBlastWave,"IME","",minFitpoint,maxFitpoint); | |
292 | myLevy->Draw("same"); | |
293 | /*myPtExp->Draw("same"); | |
294 | myMtExp->Draw("same"); | |
295 | myBoltz->Draw("same"); | |
296 | myPowerLaw->Draw("same"); | |
297 | myPowerLawPt->Draw("same");*/ | |
298 | //myBGBlastWave->Draw("same"); | |
299 | // | |
300 | legend->AddEntry(myLevy,"Levy","l"); | |
301 | /*legend->AddEntry(myPtExp,"Pt-exp","l"); | |
302 | legend->AddEntry(myMtExp,"Mt-exp","l"); | |
303 | legend->AddEntry(myBoltz,"Pt-Boltzmann","l"); | |
304 | legend->AddEntry(myPowerLaw,"PowerLaw","l"); | |
305 | legend->AddEntry(myPowerLawPt,"Pt-PowerLaw","l"); | |
306 | */ | |
307 | double AvgLowPtExtrap = myLevy->Integral(0,0.8) + myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8) + myMtExp->Integral(0,0.8); | |
308 | AvgLowPtExtrap /= 4.; | |
309 | cout<<"Levy Extrap % = "<<myLevy->Integral(0,0.8)/myLevy->Integral(0,10.)<<" Avg Extrap % = "<<AvgLowPtExtrap/myLevy->Integral(0,10.)<<endl; | |
310 | //double LowPtExtrapUncertaintyPercentage = 0.17;// Levy underestimation of Pythia | |
311 | //LowPtExtrapUncertaintyPercentage -= (AvgLowPtExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.); | |
312 | double High_lowptExtrap = (myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8))/2.; | |
313 | double Low_lowptExtrap = myMtExp->Integral(0,0.8); | |
314 | double High_LowPtExtrapUncertaintyPercentage = (High_lowptExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.); | |
315 | double Low_LowPtExtrapUncertaintyPercentage = (myLevy->Integral(0,0.8) - Low_lowptExtrap)/myLevy->Integral(0,10.); | |
316 | cout<<"high low-pt uncertainty = "<<High_LowPtExtrapUncertaintyPercentage<<endl; | |
317 | cout<<"low low-pt uncertainty = "<<Low_LowPtExtrapUncertaintyPercentage<<endl; | |
318 | ||
319 | ||
320 | //PureSysFromCutVar->Draw(); | |
321 | ||
322 | ||
323 | double dNdY = myLevy->Integral(0,10); | |
324 | double dNdY_exp = myMtExp->Integral(0,10); | |
325 | double dNdY_Stat = LevyFits_CutVar[0]->GetParError(0); | |
326 | double dNdY_Covered = myLevy->Integral(0.8,5.6); | |
327 | double dNdY_Covered_exp = myMtExp->Integral(0.8,5.6); | |
328 | double dNdY_DataPoints = 0, dNdY_DataPoints_e=0; | |
329 | for(int ii=0; ii<8; ii++){// to explicitly integrate dNdY without fit (Christina's suggestion) | |
330 | dNdY_DataPoints += Spectra_CutVar[0]->GetBinContent(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2); | |
331 | dNdY_DataPoints_e += pow(Spectra_CutVar[0]->GetBinError(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2),2); | |
332 | } | |
333 | dNdY_DataPoints_e = sqrt(dNdY_DataPoints_e); | |
334 | double dNdYSysPlus = 0; | |
335 | double dNdYSysMinus = 0; | |
336 | double CSys = 0; | |
337 | double nSys = 0; | |
338 | ////////////////////////////////////// | |
339 | // Sys errors of Levy parameters | |
340 | double sigma_dNdY_cv[12]={0}; | |
341 | double sigma_C_cv[12]={0}; | |
342 | double sigma_n_cv[12]={0}; | |
343 | double SysVar_dNdY_cv[12]={0}; | |
344 | double SysVar_C_cv[12]={0}; | |
345 | double SysVar_n_cv[12]={0}; | |
346 | ||
347 | for(int cv=1; cv<13; cv++){ | |
348 | if(MeanSysCutVar[cv-1] > 0 && MeanSysCutVar[cv-1] > sigmaSysCutVar[cv-1]) {// significant deviation | |
349 | sigma_dNdY_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_CutVar[cv]->GetParError(0),2))); | |
350 | sigma_C_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_CutVar[cv]->GetParError(1),2))); | |
351 | sigma_n_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_CutVar[cv]->GetParError(2),2))); | |
352 | SysVar_dNdY_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_CutVar[cv]->GetParameter(0)); | |
353 | SysVar_C_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_CutVar[cv]->GetParameter(1)); | |
354 | SysVar_n_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_CutVar[cv]->GetParameter(2)); | |
355 | ||
356 | if(SysVar_dNdY_cv[cv-1] > sigma_dNdY_cv[cv-1]){ | |
357 | dNdYSysPlus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2); | |
358 | dNdYSysMinus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2); | |
359 | } | |
360 | if(SysVar_C_cv[cv-1] > sigma_C_cv[cv-1]){ | |
361 | CSys += pow(SysVar_C_cv[cv-1],2)-pow(sigma_C_cv[cv-1],2); | |
362 | } | |
363 | if(SysVar_n_cv[cv-1] > sigma_n_cv[cv-1]){ | |
364 | nSys += pow(SysVar_n_cv[cv-1],2)-pow(sigma_n_cv[cv-1],2); | |
365 | } | |
366 | ||
367 | } | |
368 | ||
369 | } | |
370 | ||
371 | ||
372 | ||
373 | double sigma_dNdY_fv[4]={0}; | |
374 | double sigma_C_fv[4]={0}; | |
375 | double sigma_n_fv[4]={0}; | |
376 | double SysVar_dNdY_fv[4]={0}; | |
377 | double SysVar_C_fv[4]={0}; | |
378 | double SysVar_n_fv[4]={0}; | |
379 | for(int fv=0; fv<4; fv++){ | |
380 | sigma_dNdY_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_FitVar[fv]->GetParError(0),2))); | |
381 | SysVar_dNdY_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2)); | |
382 | sigma_C_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_FitVar[fv]->GetParError(1),2))); | |
383 | SysVar_C_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2)); | |
384 | sigma_n_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_FitVar[fv]->GetParError(2),2))); | |
385 | SysVar_n_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2)); | |
386 | ||
387 | if(fv==2){//max case for dN/dy | |
388 | //dNdYSysPlus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2); | |
389 | //dNdYSysMinus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2); | |
390 | dNdYSysPlus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2); | |
391 | dNdYSysMinus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2); | |
392 | } | |
393 | //if(fv==2) CSys += pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2) - pow(LevyFits_CutVar[0]->GetParError(1)-LevyFits_FitVar[fv]->GetParError(1),2); | |
394 | //if(fv==1) nSys += pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2) - pow(LevyFits_CutVar[0]->GetParError(2)-LevyFits_FitVar[fv]->GetParError(2),2); | |
395 | if(fv==2) CSys += pow(SysVar_C_fv[fv],2) - pow(sigma_C_fv[fv],2); | |
396 | if(fv==1) nSys += pow(SysVar_n_fv[fv],2) - pow(sigma_n_fv[fv],2); | |
397 | } | |
398 | ||
399 | ||
400 | double dNdYSys_forRatio = dNdYSysPlus; | |
401 | dNdYSysPlus += pow(0.062*dNdY,2);// INEL Norm uncertainty | |
402 | dNdYSysPlus += pow(0.04*dNdY,2);// Material Budget uncertainty | |
403 | dNdYSysPlus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty | |
404 | ||
405 | dNdYSysPlus += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (was LowPtExtrapUncertaintyPercentage * fabs(dNdY-dNdY_Covered)) | |
406 | dNdYSys_forRatio += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty | |
407 | dNdYSysPlus = sqrt(dNdYSysPlus); | |
408 | dNdYSys_forRatio = sqrt(dNdYSys_forRatio); | |
409 | // | |
410 | dNdYSysMinus += pow(0.03*dNdY,2);// INEL Norm uncertainty | |
411 | dNdYSysMinus += pow(0.04*dNdY,2);// Material Budget uncertainty | |
412 | dNdYSysMinus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty | |
413 | dNdYSysMinus += pow(Low_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (40% from Pythia study). | |
414 | dNdYSysMinus = sqrt(dNdYSysMinus); | |
415 | ||
416 | CSys = sqrt(CSys); | |
417 | nSys = sqrt(nSys); | |
418 | ||
419 | cout<<"dN/dy = "<<dNdY<<" , Covered range = "<<dNdY_Covered<<" , DataPoint Sum = "<<dNdY_DataPoints<<" , stat = "<<dNdY_DataPoints_e<<endl; | |
420 | cout<<"dN/dy stat = "<<dNdY_Stat<<" +sys = "<<dNdYSysPlus<<" -sys = "<<dNdYSysMinus<<endl; | |
421 | ||
422 | cout<<"C = "<<myLevy->GetParameter(1)<<" , stat = "<<LevyFits_CutVar[0]->GetParError(1)<<" , sys = "<<CSys<<endl; | |
423 | cout<<"n = "<<myLevy->GetParameter(2)<<" , stat = "<<LevyFits_CutVar[0]->GetParError(2)<<" , sys = "<<nSys<<endl; | |
424 | ||
425 | // mean pt | |
426 | double ptEdges[9]={0.8, 1.2, 1.6, 2.0, 2.4, 3.2, 4.0, 4.8, 5.6}; | |
427 | double meanpt_def=0, meanpt_Sys1=0, meanpt_Sys2=0, meanpt_Sys3=0;// Sys1 = Spectrum only, Sys2 = bin centers, Sys3 = low-pt extrapolation | |
428 | double yieldSum_def=0, yieldSum_Sys1=0, yieldSum_Sys2=0, yieldSum_Sys3=0; | |
429 | double meanpt_statError=0; | |
430 | meanpt_def += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8);// low-pt extrapolation | |
431 | yieldSum_def += myLevy->Integral(0,0.8); | |
432 | meanpt_Sys1 += meanpt_def; meanpt_Sys2 += meanpt_def; | |
433 | yieldSum_Sys1 += yieldSum_def; yieldSum_Sys2 += yieldSum_def; | |
434 | meanpt_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY) * myLevy->Moment(1,0.,0.8);// low-pt extrapolation uncertainty | |
435 | yieldSum_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY);// low-pt extrapolation uncertainty | |
436 | meanpt_statError += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8) * myLevy->GetParError(0)/myLevy->GetParameter(0); | |
437 | ||
438 | ||
439 | double StatErrorPtBins=0; | |
440 | for(int bin=1; bin<=8; bin++){ | |
441 | double bin_width = fabs(ptEdges[bin-1]-ptEdges[bin]); | |
442 | meanpt_def += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
443 | yieldSum_def += FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
444 | meanpt_Sys1 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * myLevy->Integral(ptEdges[bin-1],ptEdges[bin]); | |
445 | yieldSum_Sys1 += myLevy->Integral(ptEdges[bin-1],ptEdges[bin]); | |
446 | meanpt_Sys2 += (ptEdges[bin]+ptEdges[bin-1])/2. * FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
447 | yieldSum_Sys2 += FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
448 | meanpt_Sys3 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
449 | yieldSum_Sys3 += FinalSpectrum->GetBinContent(bin+1)*bin_width; | |
450 | // | |
451 | StatErrorPtBins += pow(myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * Spectra_CutVar[0]->GetBinError(bin+1)*bin_width,2); | |
452 | } | |
453 | ||
454 | meanpt_def += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation | |
455 | yieldSum_def += myLevy->Integral(5.6,10.); | |
456 | meanpt_Sys1 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation | |
457 | yieldSum_Sys1 += myLevy->Integral(5.6,10.); | |
458 | meanpt_Sys2 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation | |
459 | yieldSum_Sys2 += myLevy->Integral(5.6,10.); | |
460 | meanpt_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.))*myLevy->Moment(1,5.6,10.);// high-pt extrapolation | |
461 | yieldSum_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.)); | |
462 | // | |
463 | meanpt_statError += sqrt(StatErrorPtBins); | |
464 | meanpt_statError += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.) * myLevy->GetParError(0)/myLevy->GetParameter(0);// high-pt extrapolation | |
465 | meanpt_statError /= yieldSum_def; | |
466 | ||
467 | meanpt_def /= yieldSum_def; | |
468 | meanpt_Sys1 /= yieldSum_Sys1; | |
469 | meanpt_Sys2 /= yieldSum_Sys2; | |
470 | meanpt_Sys3 /= yieldSum_Sys3; | |
471 | ||
472 | double meanpt_Sys = sqrt(pow(meanpt_def-meanpt_Sys1,2)+pow(meanpt_def-meanpt_Sys2,2)+pow(meanpt_def-meanpt_Sys3,2)); | |
473 | ||
474 | cout<<"<pt> = "<<meanpt_def<<" , stat = "<<meanpt_statError<<" , sys = "<<meanpt_Sys<<endl; | |
475 | ||
476 | /*double meanpt = myLevy->Moment(1,0.,10.);// 0.,30 for full range | |
477 | double meanpt_CoveredLevy = myLevy->Moment(1,0.8,5.6); | |
478 | double mPt_Stat = 0; | |
479 | double mPt_Sys = 0; | |
480 | double mPt_StatSys = 0; | |
481 | double base_C = myLevy->GetParameter(1); | |
482 | double base_n = myLevy->GetParameter(2); | |
483 | ||
484 | ||
485 | // Sys | |
486 | myLevy->SetParameter(1, base_C + myLevy->GetParError(1)); | |
487 | mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2); | |
488 | // | |
489 | // | |
490 | myLevy->SetParameter(2, base_n + myLevy->GetParError(2)); | |
491 | mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2); | |
492 | mPt_StatSys = sqrt(mPt_StatSys); | |
493 | // | |
494 | // | |
495 | // Stat | |
496 | myLevy->SetParameter(1, base_C + LevyFits_CutVar[0]->GetParError(1)); | |
497 | mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2); | |
498 | // | |
499 | // | |
500 | myLevy->SetParameter(2, base_n + LevyFits_CutVar[0]->GetParError(2)); | |
501 | mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2); | |
502 | mPt_Stat = sqrt(mPt_Stat); | |
503 | // | |
504 | myLevy->SetParameter(1, base_C); myLevy->SetParameter(1, base_n); | |
505 | ||
506 | cout<<"Levy mean pt = "<<meanpt<<" , stat = "<<mPt_Stat<<" , sys = "<<sqrt( pow(mPt_StatSys,2)-pow(mPt_Stat,2) )<<endl; | |
507 | cout<<"Levy mean pt in covered range = "<<meanpt_CoveredLevy<<endl; | |
508 | */ | |
509 | ||
510 | double XidNdY = 0.0079;// (Xi^+ + Xi^-)/2. Published | |
511 | double XiStarToXi = dNdY/XidNdY; | |
512 | double XiStarToXi_Stat = sqrt(pow(dNdY_Stat/XidNdY,2) + pow(dNdY*0.0001/pow(XidNdY,2),2));// 0.0001 is stat error of Xi^{\pm} | |
513 | double XiStartoXi_Sys = dNdYSys_forRatio/XidNdY; | |
514 | ||
515 | ||
516 | cout<<"Xi*/Xi = "<<XiStarToXi<<" , stat = "<<XiStarToXi_Stat<<" , sys = "<<XiStartoXi_Sys<<endl; | |
517 | cout<<endl; | |
518 | cout<<endl; | |
519 | ||
520 | //Old QM2011 values for (Xi* + cc)/2 | |
521 | double Old_XiStarPoints[8]={0.001259, 0.0008586, 0.0006316, 0.0003604, 1.856e-4, 5.4601e-5, 2.0187e-5, 7.6559e-6}; | |
522 | double Old_XiStarPoints_e[8]={2.44e-4, 1.224e-4, 7.136e-5, 3.85e-5, 1.59e-5, 7.49e-6, 3.51e-6, 1.822e-6}; | |
523 | double pt_points[8]={1.0, 1.4, 1.8, 2.2, 2.8, 3.6, 4.4, 5.2}; | |
524 | double pt_points_e[8]={.2, .2, .2, .2, .4, .4, .4, .4}; | |
525 | TGraphErrors *gr_OldXiStar = new TGraphErrors(8, pt_points, Old_XiStarPoints, pt_points_e, Old_XiStarPoints_e); | |
526 | gr_OldXiStar->SetMarkerStyle(21); | |
527 | gr_OldXiStar->SetMarkerColor(2); | |
528 | gr_OldXiStar->SetLineColor(2); | |
529 | double NewOldRatio[8]={0}; | |
530 | for(int ii=0; ii<8; ii++){ | |
531 | NewOldRatio[ii] = FinalSpectrum->GetBinContent(ii+2)/Old_XiStarPoints[ii]; | |
532 | //cout<<FinalSpectrum->GetBinError(ii+2)/FinalSpectrum->GetBinContent(ii+2)<<endl; | |
533 | } | |
534 | TGraph *gr_NewOldRatio = new TGraph(8, pt_points, NewOldRatio); | |
535 | gr_NewOldRatio->SetMarkerStyle(20); | |
536 | gr_NewOldRatio->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
537 | gr_NewOldRatio->GetYaxis()->SetTitle("New/Old"); | |
538 | gr_NewOldRatio->SetTitle("New Spectrum Divided by QM11 Spectrum"); | |
539 | //gr_NewOldRatio->Draw("AP"); | |
540 | ||
541 | //gr_OldXiStar->Draw("P"); | |
542 | //legend->AddEntry(gr_OldXiStar,"(#Xi(1530) + cc)/2: QM11 Results","p"); | |
543 | // To Fit Old points | |
544 | /*TH1D *h_OldXiStar = (TH1D*)h_Xispectrum->Clone(); | |
545 | h_OldXiStar->Add(h_OldXiStar,-1); | |
546 | for(int i=0; i<8; i++){ | |
547 | h_OldXiStar->SetBinContent(i+2, Old_XiStarPoints[i]); | |
548 | h_OldXiStar->SetBinError(i+2, Old_XiStarPoints_e[i]); | |
549 | } | |
550 | h_OldXiStar->SetBinContent(1,+100);// first bin has no data | |
551 | h_OldXiStar->Fit(myLevy,"IME","",minFitpoint,maxFitpoint); | |
552 | */ | |
553 | //legend->Draw("same"); | |
554 | ||
555 | double StatPercentage[8]={}; | |
556 | double SysPercentage[8]={}; | |
557 | double TotErrPercentage[8]={}; | |
558 | for(int ii=0; ii<8; ii++){ | |
559 | StatPercentage[ii] = StatError[ii]/FinalSpectrum->GetBinContent(ii+2); | |
560 | SysPercentage[ii] = SysError[ii+1]/FinalSpectrum->GetBinContent(ii+2); | |
561 | TotErrPercentage[ii] = sqrt(pow(StatError[ii],2)+pow(SysError[ii+1],2))/FinalSpectrum->GetBinContent(ii+2); | |
562 | } | |
563 | TGraph *gr_StatPercentage = new TGraph(8, pt_points, StatPercentage); | |
564 | gr_StatPercentage->SetMarkerStyle(25); | |
565 | gr_StatPercentage->SetMarkerColor(4); | |
566 | gr_StatPercentage->SetLineColor(4); | |
567 | gr_StatPercentage->GetXaxis()->SetTitle("p_{T} (GeV/c)"); | |
568 | gr_StatPercentage->GetYaxis()->SetTitle("Uncertainty (fraction)"); | |
569 | gr_StatPercentage->SetTitle(""); | |
570 | gr_StatPercentage->SetMinimum(0); | |
571 | gr_StatPercentage->SetMaximum(0.2);// 0.35 to compare with Enrico | |
572 | //gr_StatPercentage->Draw("APC"); | |
573 | TGraph *gr_SysPercentage = new TGraph(8, pt_points, SysPercentage); | |
574 | gr_SysPercentage->SetMarkerStyle(25); | |
575 | gr_SysPercentage->SetMarkerColor(2); | |
576 | gr_SysPercentage->SetLineColor(2); | |
577 | //gr_SysPercentage->Draw("PC"); | |
578 | TGraph *gr_TotErrPercentage = new TGraph(8, pt_points, TotErrPercentage); | |
579 | gr_TotErrPercentage->SetMarkerStyle(25); | |
580 | gr_TotErrPercentage->SetMarkerColor(1); | |
581 | gr_TotErrPercentage->SetLineColor(1); | |
582 | //gr_TotErrPercentage->Draw("PC"); | |
583 | //legend->AddEntry(gr_StatPercentage,"Stat. uncertainty","p"); | |
584 | //legend->AddEntry(gr_SysPercentage,"Sys. uncertainty","p"); | |
585 | //legend->AddEntry(gr_TotErrPercentage,"#sqrt{Stat^{2}+Sys^{2}}","p"); | |
586 | legend->Draw("same"); | |
587 | ||
588 | //////////////////////////////// | |
589 | // Cut Var | |
590 | for(int var=1; var<13; var++){// 0 to 4 for FitVar, 1 to 13 for CutVar | |
591 | TString *Varname = new TString(); | |
592 | if(var==1) Varname->Append("Nclusters TPC"); | |
593 | if(var==2) Varname->Append("DCA PV proton"); | |
594 | if(var==3) Varname->Append("DCA PV 1st pion"); | |
595 | if(var==4) Varname->Append("DCA PV 2nd pion"); | |
596 | if(var==5) Varname->Append("DCA PV 3rd pion"); | |
597 | if(var==6) Varname->Append("DCA PV 4th pion"); | |
598 | if(var==7) Varname->Append("DCA proton-pion"); | |
599 | if(var==8) Varname->Append("DCA Lambda-pion"); | |
600 | if(var==9) Varname->Append("Decay length xy Lambda"); | |
601 | if(var==10) Varname->Append("Decay length xy Xi"); | |
602 | if(var==11) Varname->Append("Cos PointingAngle Lambda"); | |
603 | if(var==12) Varname->Append("Cos PointingAngle Xi"); | |
604 | ||
605 | cout<<Varname->Data()<<" & "; | |
606 | cout.precision(2); | |
607 | for(int par=0; par<3; par++) { | |
608 | cout<<100.*(LevyFits_CutVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/LevyFits_CutVar[0]->GetParameter(par)<<"("; | |
609 | if(par==0) cout<<100.*sigma_dNdY_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
610 | if(par==1) cout<<100.*sigma_C_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
611 | if(par==2) cout<<100.*sigma_n_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
612 | ||
613 | if((par+1)!=3) cout<<"& "; | |
614 | else cout<<" \\\\\ \\hline"; | |
615 | } | |
616 | cout<<endl; | |
617 | } | |
618 | cout<<endl; | |
619 | ||
620 | ////////////////////////////// | |
621 | // Fit Var | |
622 | for(int var=0; var<4; var++){// 0 to 4 for FitVar, 1 to 13 for CutVar | |
623 | TString *Varname = new TString(); | |
624 | if(var==0) Varname->Append("Pure Voigtian"); | |
625 | if(var==1) Varname->Append("Pure Voigtian w/o bkg subtraction"); | |
626 | if(var==2) Varname->Append("bkg norm to the right"); | |
627 | if(var==3) Varname->Append("1/2 bin-counting range"); | |
628 | ||
629 | cout<<Varname->Data()<<" & "; | |
630 | cout.precision(2); | |
631 | for(int par=0; par<3; par++) { | |
632 | cout<<100.*(LevyFits_FitVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/(LevyFits_CutVar[0]->GetParameter(par))<<"("; | |
633 | if(par==0) cout<<100.*sigma_dNdY_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
634 | if(par==1) cout<<100.*sigma_C_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
635 | if(par==2) cout<<100.*sigma_n_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" "; | |
636 | ||
637 | if((par+1)!=3) cout<<"& "; | |
638 | else cout<<" \\\\\ \\hline"; | |
639 | } | |
640 | cout<<endl; | |
641 | } | |
642 | cout<<endl; | |
643 | ||
644 | /////////////////////////////////////////// | |
645 | ||
646 | /* | |
647 | cout<<"Yields"<<endl; | |
648 | // print yields | |
649 | for(int var=1; var<=13; var++){ | |
650 | TString *Varname = new TString(); | |
651 | if(var==1) Varname->Append("Standard "); | |
652 | else {Varname->Append("Cut Var "); *Varname += var-1;} | |
653 | ||
654 | cout<<Varname->Data()<<" & "; | |
655 | cout.precision(4);// 4 or 6 | |
656 | for(int ptbin=1; ptbin<=8; ptbin++){ | |
657 | //cout<<Yields_CutVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[var-1]->GetBinContent(ptbin)<<" "; | |
658 | cout<<Eff_CutVar[var-1]->GetBinContent(ptbin)<<" "; | |
659 | if((ptbin)!=8) cout<<"& "; | |
660 | else cout<<" \\\\\ \\hline"; | |
661 | } | |
662 | cout<<endl; | |
663 | } | |
664 | */ | |
665 | /* for(int var=1; var<5; var++){ | |
666 | TString *Varname = new TString("Fit Var "); | |
667 | *Varname += var; | |
668 | cout<<Varname->Data()<<" & "; | |
669 | cout.precision(6); | |
670 | for(int ptbin=1; ptbin<=8; ptbin++){ | |
671 | cout<<Yields_FitVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[0]->GetBinContent(ptbin)<<" "; | |
672 | if((ptbin)!=8) cout<<"& "; | |
673 | else cout<<" \\\\\ \\hline"; | |
674 | } | |
675 | cout<<endl; | |
676 | }*/ | |
677 | cout<<endl; | |
678 | ||
679 | //TString *name=new TString("../AnaNotes/XiStar/XiStar_yield_ErrorsCutVar"); | |
680 | //*name += CutVarN; | |
681 | /*TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSys"); | |
682 | *name += CutVarLow; | |
683 | name->Append(".png"); | |
684 | c1->SaveAs(name->Data());*/ | |
685 | //TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSummary"); | |
686 | //*name += CutVarLow; | |
687 | //name->Append(".png"); | |
688 | //can->SaveAs(name->Data()); | |
689 | ||
690 | } | |
691 | ||
692 | ||
693 | //________________________________________________________________________ | |
694 | double myLevyPtFunc(Double_t *x, Double_t *par) | |
695 | { | |
696 | Double_t lMass=0; | |
697 | lMass = 1.5318; //Xi* mass | |
698 | ||
699 | ||
700 | Double_t ldNdy = par[0]; // dN/dy | |
701 | Double_t l2pi = 2*TMath::Pi(); // 2pi (cancels in return statement) | |
702 | Double_t lTemp = par[1]; // Temperature | |
703 | Double_t lPower = par[2]; // power=n | |
704 | ||
705 | Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2))); | |
706 | Double_t lInPower = 1 + (TMath::Sqrt(x[0]*x[0]+lMass*lMass)-lMass) / (lPower*lTemp); | |
707 | ||
708 | return l2pi * ldNdy * x[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower); | |
709 | } | |
710 | //________________________________________________________________________ | |
711 | double myMtExpFunc(Double_t *x, Double_t *par) | |
712 | { | |
713 | Double_t lMass=0; | |
714 | lMass = 1.5318; //Xi* mass | |
715 | ||
716 | return 2*TMath::Pi()*x[0]*par[0]*exp(-TMath::Sqrt(x[0]*x[0]+lMass*lMass)/par[1]); | |
717 | } | |
718 | //________________________________________________________________________ | |
719 | double myPtExpFunc(Double_t *x, Double_t *par) | |
720 | { | |
721 | return 2*TMath::Pi()*x[0]*par[0]*exp(-x[0]/par[1]); | |
722 | } | |
723 | //________________________________________________________________________ | |
724 | double myBoltzFunc(Double_t *x, Double_t *par) | |
725 | { | |
726 | Double_t lMass=0; | |
727 | lMass = 1.5318; //Xi* mass | |
728 | Double_t Mt = sqrt(pow(x[0],2)+pow(lMass,2)); | |
729 | return par[0]*x[0]*Mt*exp(-Mt/par[1]); | |
730 | ||
731 | } | |
732 | //________________________________________________________________________ | |
733 | double myPowerLawPtFunc(Double_t *x, Double_t *par) | |
734 | { | |
735 | return ( x[0]*par[0]*pow(1+x[0]/par[1], -par[2]) ); | |
736 | } | |
737 | //________________________________________________________________________ | |
738 | double myPowerLawFunc(Double_t *x, Double_t *par) | |
739 | { | |
740 | return ( par[0]*pow(1+x[0]/par[1], -par[2]) ); | |
741 | } | |
742 | //________________________________________________________________________ | |
743 | double myBoltsBWFunc(Double_t *x, Double_t *par) | |
744 | { | |
745 | double pT = x[0]; | |
746 | double mass = par[0]; | |
747 | double beta = par[1]; | |
748 | double temp = par[2]; | |
749 | double n = par[3]; | |
750 | double norm = par[4]; | |
751 | ||
752 | static TF1 * fIntBG = 0; | |
753 | if(!fIntBG) | |
754 | fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5); | |
755 | ||
756 | fIntBG->SetParameters(mass, pT, beta, temp,n); | |
757 | double result = fIntBG->Integral(0,1); | |
758 | return result*norm; | |
759 | } | |
760 | //________________________________________________________________________ | |
761 | double IntegrandBG(const double * x, const double* par){ | |
762 | double x0 = x[0]; | |
763 | ||
764 | double mass = par[0]; | |
765 | double pT = par[1]; | |
766 | double beta_max = par[2]; | |
767 | double temp = par[3]; | |
768 | Double_t n = par[4]; | |
769 | ||
770 | // Keep beta within reasonable limits | |
771 | Double_t beta = beta_max * TMath::Power(x0, n); | |
772 | if (beta > 0.9999999999999999) beta = 0.9999999999999999; | |
773 | ||
774 | double mT = TMath::Sqrt(mass*mass+pT*pT); | |
775 | ||
776 | double rho0 = TMath::ATanH(beta); | |
777 | double arg00 = pT*TMath::SinH(rho0)/temp; | |
778 | if (arg00 > 700.) arg00 = 700.; // avoid FPE | |
779 | double arg01 = mT*TMath::CosH(rho0)/temp; | |
780 | double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01); | |
781 | ||
782 | // printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01); | |
783 | ||
784 | return f0; | |
785 | } |