AliPhysics  d497afb (d497afb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TotalCharmCrossSection.C
Go to the documentation of this file.
1 //Author: Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de
2 //Routine to extrapolate D meson cross sections using FONLL calculations
3 #include <TFile.h>
4 #include <TTree.h>
5 #include <TMath.h>
6 #include <TH1F.h>
7 #include <iostream>
8 #include <TH1.h>
9 #include <TH1D.h>
10 #include <TKey.h>
11 #include <TH2F.h>
12 #include <TGraph.h>
13 #include <TPave.h>
14 #include <TPaveText.h>
15 #include <TPad.h>
16 #include <TCanvas.h>
17 #include <TLegend.h>
18 #include <TMultiGraph.h>
19 #include <TROOT.h>
20 #include <TAxis.h>
21 #include <TGraphAsymmErrors.h>
22 #include <TStyle.h>
23 #include <fstream>
24 
25 void TotalCharmCrossSection() {
26 
27  // FONLL predictions from _pdf.rooT files
28  TFile *D0_05 = TFile::Open("D07TeV_FONLL_y0.5_pdf.root");
29  TFile *Dst_05 = TFile::Open("Dstar7TeV_FONLL_y0.5_pdf.root");
30  TFile *Dplus_05 = TFile::Open("Dplus7TeV_FONLL_y0.5_pdf.root");
31 
32  TFile *D0_ally = TFile::Open("D07TeV_FONLL_y12_pdf.root");
33  TFile *Dst_ally = TFile::Open("Dstar7TeV_FONLL_y12_pdf.root");
34  TFile *Dplus_ally = TFile::Open("Dplus7TeV_FONLL_y12_pdf.root");
35 
36  //Spectra for error propagation
37  TFile *specD0 = TFile::Open("HFPtSpectrum_D0Kpi_combinedFD_rebinnedth_150311_newsigma_7TeV.root");
38  TFile *specDst = TFile::Open("HFPtSpectrum_DstarD0pi_combineFD_rebinnedth_newnorm_120311.root");
39  TFile *specDplus = TFile::Open("HFPtSpectrum_DplusKpipi_combinedFD_150311_newsigma_7TeV.root");
40 
41  TH1F *hD0stat = (TH1F*)specD0->Get("histoSigmaCorr");
42  TH1F *hDststat = (TH1F*)specDst->Get("histoSigmaCorr");
43  TH1F *hDplusstat = (TH1F*)specDplus->Get("histoSigmaCorr");
44 
45 
46  //Initialising the histo arrays
47  TH1F **histosD0_ally, **histosD0_05, **histosDst_ally, **histosDst_05, **histosDplus_ally, **histosDplus_05;
48 
49  histosD0_ally = new TH1F*[11];
50  histosD0_05 = new TH1F*[11];
51  histosDst_ally = new TH1F*[11];
52  histosDst_05 = new TH1F*[11];
53  histosDplus_ally = new TH1F*[11];
54  histosDplus_05 = new TH1F*[11];
55 
56 
57 //Getting all the necessary histos for extrapolation ratios
58 
59  histosD0_ally[0] = (TH1F*)D0_ally->Get("FONLL_central_prediction");
60  histosD0_ally[1] = (TH1F*)D0_ally->Get("FONLL_fr55_prediction");
61  histosD0_ally[2] = (TH1F*)D0_ally->Get("FONLL_fr22_prediction");
62  histosD0_ally[3] = (TH1F*)D0_ally->Get("FONLL_fr21_prediction");
63  histosD0_ally[4] = (TH1F*)D0_ally->Get("FONLL_fr12_prediction");
64  histosD0_ally[5] = (TH1F*)D0_ally->Get("FONLL_fr15_prediction");
65  histosD0_ally[6] = (TH1F*)D0_ally->Get("FONLL_fr51_prediction");
66  histosD0_ally[7] = (TH1F*)D0_ally->Get("FONLL_minmass_prediction");
67  histosD0_ally[8] = (TH1F*)D0_ally->Get("FONLL_maxmass_prediction");
68  histosD0_ally[9] = (TH1F*)D0_ally->Get("FONLL_minpdf_prediction");
69  histosD0_ally[10] = (TH1F*)D0_ally->Get("FONLL_maxpdf_prediction");
70 
71  histosD0_05[0] = (TH1F*)D0_05->Get("FONLL_central_prediction");
72  histosD0_05[1] = (TH1F*)D0_05->Get("FONLL_fr55_prediction");
73  histosD0_05[2] = (TH1F*)D0_05->Get("FONLL_fr22_prediction");
74  histosD0_05[3] = (TH1F*)D0_05->Get("FONLL_fr21_prediction");
75  histosD0_05[4] = (TH1F*)D0_05->Get("FONLL_fr12_prediction");
76  histosD0_05[5] = (TH1F*)D0_05->Get("FONLL_fr15_prediction");
77  histosD0_05[6] = (TH1F*)D0_05->Get("FONLL_fr51_prediction");
78  histosD0_05[7] = (TH1F*)D0_05->Get("FONLL_minmass_prediction");
79  histosD0_05[8] = (TH1F*)D0_05->Get("FONLL_maxmass_prediction");
80  histosD0_05[9] = (TH1F*)D0_05->Get("FONLL_minpdf_prediction");
81  histosD0_05[10] = (TH1F*)D0_05->Get("FONLL_maxpdf_prediction");
82 
83  histosDst_ally[0] = (TH1F*)Dst_ally->Get("FONLL_central_prediction");
84  histosDst_ally[1] = (TH1F*)Dst_ally->Get("FONLL_fr55_prediction");
85  histosDst_ally[2] = (TH1F*)Dst_ally->Get("FONLL_fr22_prediction");
86  histosDst_ally[3] = (TH1F*)Dst_ally->Get("FONLL_fr21_prediction");
87  histosDst_ally[4] = (TH1F*)Dst_ally->Get("FONLL_fr12_prediction");
88  histosDst_ally[5] = (TH1F*)Dst_ally->Get("FONLL_fr15_prediction");
89  histosDst_ally[6] = (TH1F*)Dst_ally->Get("FONLL_fr51_prediction");
90  histosDst_ally[7] = (TH1F*)Dst_ally->Get("FONLL_minmass_prediction");
91  histosDst_ally[8] = (TH1F*)Dst_ally->Get("FONLL_maxmass_prediction");
92  histosDst_ally[9] = (TH1F*)Dst_ally->Get("FONLL_minpdf_prediction");
93  histosDst_ally[10] = (TH1F*)Dst_ally->Get("FONLL_maxpdf_prediction");
94 
95 
96  histosDst_05[0] = (TH1F*)Dst_05->Get("FONLL_central_prediction");
97  histosDst_05[1] = (TH1F*)Dst_05->Get("FONLL_fr55_prediction");
98  histosDst_05[2] = (TH1F*)Dst_05->Get("FONLL_fr22_prediction");
99  histosDst_05[3] = (TH1F*)Dst_05->Get("FONLL_fr21_prediction");
100  histosDst_05[4] = (TH1F*)Dst_05->Get("FONLL_fr12_prediction");
101  histosDst_05[5] = (TH1F*)Dst_05->Get("FONLL_fr15_prediction");
102  histosDst_05[6] = (TH1F*)Dst_05->Get("FONLL_fr51_prediction");
103  histosDst_05[7] = (TH1F*)Dst_05->Get("FONLL_minmass_prediction");
104  histosDst_05[8] = (TH1F*)Dst_05->Get("FONLL_maxmass_prediction");
105  histosDst_05[9] = (TH1F*)Dst_05->Get("FONLL_minpdf_prediction");
106  histosDst_05[10] = (TH1F*)Dst_05->Get("FONLL_maxpdf_prediction");
107 
108  histosDplus_ally[0] = (TH1F*)Dplus_ally->Get("FONLL_central_prediction");
109  histosDplus_ally[1] = (TH1F*)Dplus_ally->Get("FONLL_fr55_prediction");
110  histosDplus_ally[2] = (TH1F*)Dplus_ally->Get("FONLL_fr22_prediction");
111  histosDplus_ally[3] = (TH1F*)Dplus_ally->Get("FONLL_fr21_prediction");
112  histosDplus_ally[4] = (TH1F*)Dplus_ally->Get("FONLL_fr12_prediction");
113  histosDplus_ally[5] = (TH1F*)Dplus_ally->Get("FONLL_fr15_prediction");
114  histosDplus_ally[6] = (TH1F*)Dplus_ally->Get("FONLL_fr51_prediction");
115  histosDplus_ally[7] = (TH1F*)Dplus_ally->Get("FONLL_minmass_prediction");
116  histosDplus_ally[8] = (TH1F*)Dplus_ally->Get("FONLL_maxmass_prediction");
117  histosDplus_ally[9] = (TH1F*)Dplus_ally->Get("FONLL_minpdf_prediction");
118  histosDplus_ally[10] = (TH1F*)Dplus_ally->Get("FONLL_maxpdf_prediction");
119 
120  histosDplus_05[0] = (TH1F*)Dplus_05->Get("FONLL_central_prediction");
121  histosDplus_05[1] = (TH1F*)Dplus_05->Get("FONLL_fr55_prediction");
122  histosDplus_05[2] = (TH1F*)Dplus_05->Get("FONLL_fr22_prediction");
123  histosDplus_05[3] = (TH1F*)Dplus_05->Get("FONLL_fr21_prediction");
124  histosDplus_05[4] = (TH1F*)Dplus_05->Get("FONLL_fr12_prediction");
125  histosDplus_05[5] = (TH1F*)Dplus_05->Get("FONLL_fr15_prediction");
126  histosDplus_05[6] = (TH1F*)Dplus_05->Get("FONLL_fr51_prediction");
127  histosDplus_05[7] = (TH1F*)Dplus_05->Get("FONLL_minmass_prediction");
128  histosDplus_05[8] = (TH1F*)Dplus_05->Get("FONLL_maxmass_prediction");
129  histosDplus_05[9] = (TH1F*)Dplus_05->Get("FONLL_minpdf_prediction");
130  histosDplus_05[10] = (TH1F*)Dplus_05->Get("FONLL_maxpdf_prediction");
131 
132 //Output file for extrap. ratios
133 
134  ofstream outD0("outputD0.txt"), outDst("outputDst.txt"), outDplus("outputDplus.txt");
135 
136 //Initialise doubles for max, min, central for each meson
137 Double_t ratscalmaxD0 = 0, ratscalminD0 = 9999, ratscalmaxDst = 0, ratscalminDst = 9999, ratscalmaxDplus = 0, ratscalminDplus = 9999,ratmassmaxD0 = 0, ratmassminD0 = 9999, ratmassmaxDst = 0, ratmassminDst = 9999, ratmassmaxDplus = 0, ratmassminDplus = 9999, ratpdfmaxD0 = 0, ratpdfminD0 = 9999, ratpdfmaxDst = 0, ratpdfminDst = 9999, ratpdfmaxDplus = 0, ratpdfminDplus = 9999, ratD0cent, ratDstcent, ratDpluscent,
138 ratD0[11], ratDst[11], ratDplus[11];
139 
140 
141  //Central value is with central parameters
142  ratD0cent = histosD0_ally[0]->Integral("width")/histosD0_05[0]->Integral(21,120,"width");
143  ratDstcent = histosDst_ally[0]->Integral("width")/histosDst_05[0]->Integral(21,120,"width");
144  ratDpluscent = histosDplus_ally[0]->Integral("width")/histosDplus_05[0]->Integral(21,120,"width");
145 //Headers
146 outD0 << "Name\tpT\ty\toverall\n";
147 outDst << "Name\tpT\ty\toverall\n";
148 outDplus << "Name\tpT\ty\toverall\n";
149 
150 
151 //Determining maximum and minimum ratios, and writing all ratios to file for later analysis
152  for (Int_t i = 0; i<=10; i++) {
153 
154  ratD0[i] = histosD0_ally[i]->Integral("width")/histosD0_05[i]->Integral(21,120,"width");
155 
156  //Scale variation for histos between i=1 and i=6
157  if (ratD0[i] > ratscalmaxD0 && 0 <= i <= 6) {ratscalmaxD0 = ratD0[i];}
158  if (ratD0[i] < ratscalminD0 && 0 <= i <= 6) {ratscalminD0 = ratD0[i];}
159 
160  //Mass variation for 0, 7, 8
161  if (ratD0[i] > ratmassmaxD0 && (i == 0 || i == 7 || i == 8)) {ratmassmaxD0 = ratD0[i];}
162  if (ratD0[i] < ratmassminD0 && (i == 0 || i == 7 || i == 8)) {ratmassminD0 = ratD0[i];}
163 
164  //PDF uncertainty is 0, 9 and 10
165  if (ratD0[i] > ratpdfmaxD0 && (i == 0 || i == 9 || i == 10)) {ratpdfmaxD0 = ratD0[i];}
166  if (ratD0[i] < ratpdfminD0 && (i == 0 || i == 9 || i == 10)) {ratpdfminD0 = ratD0[i];}
167 
168  outD0 << histosD0_05[i]->GetName() << '\t'
169  << histosD0_05[i]->Integral("width")/histosD0_05[i]->Integral(21,120,"width") << '\t'
170  << histosD0_ally[i]->Integral("width")/histosD0_05[i]->Integral("width") << '\t'
171  << ratD0[i] << endl;
172 
173  ratDst[i] = histosDst_ally[i]->Integral("width")/histosDst_05[i]->Integral(21,120,"width");
174  if (ratDst[i] > ratscalmaxDst && 0 <= i <= 6) {ratscalmaxDst = ratDst[i];}
175  if (ratDst[i] < ratscalminDst && 0 <= i <= 6) {ratscalminDst = ratDst[i];}
176 
177  if (ratDst[i] > ratmassmaxDst && (i == 0 || i == 7 || i == 8)) {ratmassmaxDst = ratDst[i];}
178  if (ratDst[i] < ratmassminDst && (i == 0 || i == 7 || i == 8)) {ratmassminDst = ratDst[i];}
179 
180  if (ratDst[i] > ratpdfmaxDst && (i == 0 || i == 9 || i == 10)) {ratpdfmaxDst = ratDst[i];}
181  if (ratDst[i] < ratpdfminDst && (i == 0 || i == 9 || i == 10)) {ratpdfminDst = ratDst[i];}
182 
183  outDst << histosDst_05[i]->GetName() << '\t'
184  << histosDst_05[i]->Integral("width")/histosDst_05[i]->Integral(21,120,"width") << '\t'
185  << histosDst_ally[i]->Integral("width")/histosDst_05[i]->Integral("width") << '\t'
186  << ratDst[i] << endl;
187 
188 
189  ratDplus[i] = histosDplus_ally[i]->Integral("width")/histosDplus_05[i]->Integral(21,120,"width");
190  if (ratDplus[i] > ratscalmaxDplus && 0 <= i <= 6) {ratscalmaxDplus = ratDplus[i];}
191  if (ratDplus[i] < ratscalminDplus && 0 <= i <= 6) {ratscalminDplus = ratDplus[i];}
192 
193  if (ratDplus[i] > ratmassmaxDplus && (i == 0 || i == 7 || i == 8)) {ratmassmaxDplus = ratDplus[i];}
194  if (ratDplus[i] < ratmassminDplus && (i == 0 || i == 7 || i == 8)) {ratmassminDplus = ratDplus[i];}
195 
196  if (ratDplus[i] > ratpdfmaxDplus && (i == 0 || i == 9 || i == 10)) {ratpdfmaxDplus = ratDplus[i];}
197  if (ratDplus[i] < ratpdfminDplus && (i == 0 || i == 9 || i == 10)) {ratpdfminDplus = ratDplus[i];}
198 
199  outDplus << histosDplus_05[i]->GetName() << '\t'
200  << histosDplus_05[i]->Integral("width")/histosDplus_05[i]->Integral(21,120,"width") << '\t'
201  << histosDplus_ally[i]->Integral("width")/histosDplus_05[i]->Integral("width") << '\t'
202  << ratDplus[i] << endl;
203 
204 
205 
206  }
207  //Maximum and minimum error on centrals:
208  Double_t ratscalD0UpErr, ratscalD0LowErr, ratscalDstUpErr, ratscalDstLowErr, ratscalDplusUpErr, ratscalDplusLowErr, ratmassD0UpErr, ratmassD0LowErr, ratmassDstUpErr, ratmassDstLowErr, ratmassDplusUpErr, ratmassDplusLowErr, ratpdfD0UpErr, ratpdfD0LowErr, ratpdfDstUpErr, ratpdfDstLowErr, ratpdfDplusUpErr, ratpdfDplusLowErr, ratD0UpErr, ratD0LowErr, ratDstUpErr, ratDstLowErr, ratDplusUpErr, ratDplusLowErr;
209 
210  ratscalD0UpErr = ratscalmaxD0 - ratD0cent;
211  ratscalD0LowErr = ratD0cent - ratscalminD0;
212  ratscalDstUpErr = ratscalmaxDst - ratDstcent;
213  ratscalDstLowErr = ratDstcent - ratscalminDst;
214  ratscalDplusUpErr = ratscalmaxDplus - ratDpluscent;
215  ratscalDplusLowErr = ratDpluscent - ratscalminDplus;
216 
217  ratmassD0UpErr = ratmassmaxD0 - ratD0cent;
218  ratmassD0LowErr = ratD0cent - ratmassminD0;
219  ratmassDstUpErr = ratmassmaxDst - ratDstcent;
220  ratmassDstLowErr = ratDstcent - ratmassminDst;
221  ratmassDplusUpErr = ratmassmaxDplus - ratDpluscent;
222  ratmassDplusLowErr = ratDpluscent - ratmassminDplus;
223 
224  ratpdfD0UpErr = ratpdfmaxD0 - ratD0cent;
225  ratpdfD0LowErr = ratD0cent - ratpdfminD0;
226  ratpdfDstUpErr = ratpdfmaxDst - ratDstcent;
227  ratpdfDstLowErr = ratDstcent - ratpdfminDst;
228  ratpdfDplusUpErr = ratpdfmaxDplus - ratDpluscent;
229  ratpdfDplusLowErr = ratDpluscent - ratpdfminDplus;
230 
231  //Overall uncertainties are mass uncertainty and scale uncertainty added in quadrature
232  ratD0UpErr = sqrt(ratscalD0UpErr*ratscalD0UpErr+ratmassD0UpErr*ratmassD0UpErr+ratpdfD0UpErr*ratpdfD0UpErr);
233  ratD0LowErr = sqrt(ratscalD0LowErr*ratscalD0LowErr+ratmassD0LowErr*ratmassD0LowErr+ratpdfD0LowErr*ratpdfD0LowErr);
234  ratDstUpErr = sqrt(ratscalDstUpErr*ratscalDstUpErr+ratmassDstUpErr*ratmassDstUpErr+ratpdfDstUpErr*ratpdfDstUpErr);
235  ratDstLowErr = sqrt(ratscalDstLowErr*ratscalDstLowErr+ratmassDstLowErr*ratmassDstLowErr+ratpdfDstLowErr*ratpdfDstLowErr);
236  ratDplusUpErr = sqrt(ratscalDplusUpErr*ratscalDplusUpErr+ratmassDplusUpErr*ratmassDplusUpErr+ratpdfDplusUpErr*ratpdfDplusUpErr);
237  ratDplusLowErr = sqrt(ratscalDplusLowErr*ratscalDplusLowErr+ratmassDplusLowErr*ratmassDplusLowErr+ratpdfDplusLowErr*ratpdfDplusLowErr);
238 
239 
240  //Output ratios to file for future reference:
241 
242  ofstream outratios("ratios.txt");
243  outratios << "D0\ncentral\t" << ratD0cent << "\nscales:\t+" << ratscalD0UpErr << "\t-" << ratscalD0LowErr << "\nmass:\t+" << ratmassD0UpErr << "\t-" << ratmassD0LowErr << "\nPDF:\t+" << ratpdfD0UpErr << "\t-" << ratpdfD0LowErr << "\nOverall:\t+" << ratD0UpErr << "\t-" << ratD0LowErr
244 << "\n\nD*\ncentral\t" << ratDstcent << "\nscales:\t+" << ratscalDstUpErr << "\t-" << ratscalDstLowErr << "\nmass:\t+" << ratmassDstUpErr << "\t-" << ratmassDstLowErr << "\nPDF:\t+" << ratpdfDstUpErr << "\t-" << ratpdfDstLowErr << "\nOverall:\t+" << ratDstUpErr << "\t-" << ratDstLowErr
245 << "\n\nD+\ncentral\t" << ratDpluscent << "\nscales:\t+" << ratscalDplusUpErr << "\t-" << ratscalDplusLowErr << "\nmass:\t+" << ratmassDplusUpErr << "\t-" << ratmassDplusLowErr << "\nPDF:\t+" << ratpdfDplusUpErr << "\t-" << ratpdfDplusLowErr << "\nOverall:\t+" << ratDplusUpErr << "\t-" << ratDplusLowErr << '\n';
246  //Visible spectra:
247  Double_t D0stati[6], Dststati[6], Dplusstati[6];
248 
249  Double_t D0stat =0, Dststat=0, Dplusstat=0, D0statsq=0, Dststatsq=0, Dplusstatsq=0, D0systhi=0, D0systlo=0, Dstsysthi=0, Dstsystlo=0, Dplussysthi=0, Dplussystlo=0;
250  for (Int_t i = 4; i<=9; i++) {
251  D0stati[i-4] = hD0stat->GetBinError(i)/1e6/3.89e-2;
252  Dststati[i-4] = hDststat->GetBinError(i-1)/1e6/3.89e-2/67.7e-2;
253  Dplusstati[i-4] = hDplusstat->GetBinError(i-3)/1e6/0.0922;
254 
255  //Statistical added in quadrature (bin widths included)
256  D0statsq += pow(D0stati[i-4]*hD0stat->GetBinWidth(i),2);
257  Dststatsq += pow(Dststati[i-4]*hDststat->GetBinWidth(i-1),2);
258  Dplusstatsq += pow(Dplusstati[i-4]*hDplusstat->GetBinWidth(i-3),2);
259 
260  }
261 
262 
263  //Rooting to get actual values from delta-sigma-squared
264  D0stat = sqrt(D0statsq); Dststat = sqrt(Dststatsq); Dplusstat = sqrt(Dplusstatsq);
265 
266  //Corrected systematics are not in a ROOT file, input here manually
267  D0systhi = 967844/1e6/3.89e-2;
268  D0systlo = 1849380/1e6/3.89e-2;
269  Dstsysthi = 232479/1e6/3.89e-2/67.7e-2;
270  Dstsystlo = 379356/1e6/3.89e-2/67.7e-2;
271  Dplussysthi = 1710680/1e6/0.0922;
272  Dplussystlo = 3402420/1e6/0.0922;
273 
274 
275  //Visible D cross sections from data file, "width" switch incorporates bin widths
276  Double_t D0vis = hD0stat->Integral("width")/1e6/3.89e-2,
277  Dstvis = hDststat->Integral("width")/1e6/67.7e-2/3.89e-2,
278  Dplusvis = hDplusstat->Integral("width")/1e6/0.0922;
279 
280  //Branching ratio errors as percentage:
281  Double_t ErrBrD0Perc = 0.0002785/0.0216673, ErrBrDstPerc = 9.29153e-5/0.0062678, ErrBrDplusPerc = 0.0004746/0.0208372;
282 
283  //Luminosity errors as percentages
284 
285  Double_t ErrLumD0Perc=0.07, ErrLumDstPerc=0.07, ErrLumDplusPerc=0.07;
286 
287  //Initialising doubles for final results:
288  Double_t crossD0, crossD0stat, crossD0systhi, crossD0systlo, crossD0lum, crossD0br, crossD0extrhi, crossD0extrlo, crossDst, crossDststat, crossDstsysthi, crossDstsystlo, crossDstlum, crossDstbr, crossDstextrhi, crossDstextrlo, crossDplus, crossDplusstat, crossDplussysthi, crossDplussystlo, crossDpluslum, crossDplusbr, crossDplusextrhi, crossDplusextrlo;
289 
290  //Also initialising for final ccbar cross sec:
291  Double_t ccbarD0, ccbarD0stat, ccbarD0systhi, ccbarD0systlo, ccbarD0lum, ccbarD0extrhi, ccbarD0extrlo, ccbarDst, ccbarDststat, ccbarDstsysthi, ccbarDstsystlo, ccbarDstlum, ccbarDstextrhi, ccbarDstextrlo, ccbarDplus, ccbarDplusstat, ccbarDplussysthi, ccbarDplussystlo, ccbarDpluslum, ccbarDplusextrhi, ccbarDplusextrlo;
292 
293  //Cross section is visible multiplied by central ratio
294  crossD0 = D0vis * ratD0cent;
295  crossDst = Dstvis * ratDstcent;
296  crossDplus = Dplusvis * ratDpluscent;
297 
298  //Stat on each is statitstical of visible multiplied by central ratio
299  crossD0stat = D0stat * ratD0cent;
300  crossDststat = Dststat * ratDstcent;
301  crossDplusstat = Dplusstat * ratDpluscent;
302 
303  //Systematic the same, but upper and lower separately
304  crossD0systhi = D0systhi * ratD0cent;
305  crossD0systlo = D0systlo * ratD0cent;
306  crossDstsysthi = Dstsysthi * ratDstcent;
307  crossDstsystlo = Dstsystlo * ratDstcent;
308  crossDplussysthi = Dplussysthi * ratDpluscent;
309  crossDplussystlo = Dplussystlo * ratDpluscent;
310 
311  //Luminosity errors are percentage multiplied by central result
312  crossD0lum = ErrLumD0Perc * crossD0;
313  crossDstlum = ErrLumDstPerc * crossDst;
314  crossDpluslum = ErrLumDplusPerc * crossDplus;
315 
316  //Similarly for branching ratios
317  crossD0br = ErrBrD0Perc * crossD0;
318  crossDstbr = ErrBrDstPerc * crossDst;
319  crossDplusbr = ErrBrDplusPerc * crossDplus;
320 
321  //Extrapolation hi/lo errors are overall hi/lo ratio errors multiplied by visible
322  crossD0extrhi = ratD0UpErr * D0vis;
323  crossD0extrlo = ratD0LowErr * D0vis;
324  crossDstextrhi = ratDstUpErr * Dstvis;
325  crossDstextrlo = ratDstLowErr * Dstvis;
326  crossDplusextrhi = ratDplusUpErr * Dplusvis;
327  crossDplusextrlo = ratDplusLowErr * Dplusvis;
328 
329  //ccbar cross section: divide result by fragmentation ratio. This is calculated the same way for all uncertainties
330  //Fragmentations:
331  Double_t D0frag = 0.557, Dstfrag = 0.238, Dplusfrag = 0.226;
332 
333  //D0
334  ccbarD0 = crossD0/D0frag;
335  ccbarD0stat = crossD0stat/D0frag;
336  ccbarD0systhi = crossD0systhi/D0frag;
337  ccbarD0systlo = crossD0systlo/D0frag;
338  ccbarD0lum = crossD0lum/D0frag;
339  ccbarD0extrhi = crossD0extrhi/D0frag;
340  ccbarD0extrlo = crossD0extrlo/D0frag;
341 
342  //D*
343  ccbarDst = crossDst/Dstfrag;
344  ccbarDststat = crossDststat/Dstfrag;
345  ccbarDstsysthi = crossDstsysthi/Dstfrag;
346  ccbarDstsystlo = crossDstsystlo/Dstfrag;
347  ccbarDstlum = crossDstlum/Dstfrag;
348  ccbarDstextrhi = crossDstextrhi/Dstfrag;
349  ccbarDstextrlo = crossDstextrlo/Dstfrag;
350 
351  //D+
352  ccbarDplus = crossDplus/Dplusfrag;
353  ccbarDplusstat = crossDplusstat/Dplusfrag;
354  ccbarDplussysthi = crossDplussysthi/Dplusfrag;
355  ccbarDplussystlo = crossDplussystlo/Dplusfrag;
356  ccbarDpluslum = crossDpluslum/Dplusfrag;
357  ccbarDplusextrhi = crossDplusextrhi/Dplusfrag;
358  ccbarDplusextrlo = crossDplusextrlo/Dplusfrag;
359 
360 
361  //To combine: Add D0 + D+, divide by total frag. ratio of the two. Take weighted average of this value with D* using 1/(stat)^2 as weight.
362 
363 
364  //Initialising variables:
365  Double_t totalfrag, D0Dplusccbar, D0Dplusccbarstat, D0Dplusccbarsysthi, D0Dplusccbarsystlo, D0Dplusccbarlum, D0Dplusccbarextrhi, D0Dplusccbarextrlo;
366 
367  totalfrag = D0frag + Dplusfrag;
368 
369  //Propagating this method through:
370  D0Dplusccbar = (crossD0 + crossDplus)/totalfrag;
371  //Errors are quadratic sum of D0 and D+ errors, divided by total frag. ratio
372  D0Dplusccbarstat = sqrt(crossD0stat * crossD0stat + crossDplusstat * crossDplusstat)/totalfrag;
373  D0Dplusccbarsysthi = sqrt(crossD0systhi * crossD0systhi + crossDplussysthi * crossDplussysthi)/totalfrag;
374  D0Dplusccbarsystlo = sqrt(crossD0systlo * crossD0systlo + crossDplussystlo * crossDplussystlo)/totalfrag;
375  D0Dplusccbarlum = sqrt(crossD0lum * crossD0lum + crossDpluslum * crossDpluslum)/totalfrag;
376  D0Dplusccbarextrhi = sqrt(crossD0extrhi * crossD0extrhi + crossDplusextrhi * crossDplusextrhi)/totalfrag;
377  D0Dplusccbarextrlo = sqrt(crossD0extrlo * crossD0extrlo + crossDplusextrlo * crossDplusextrlo)/totalfrag;
378 
379 cout << "Just to test, this is D0Dplusccbar: " << D0Dplusccbar << "\n\n\n";
380 
381 
382  //Weighted average:
383  Double_t weightD0Dplus, weightDst, totalweight, overallccbar, overallccbarstat, overallccbarsysthi, overallccbarsystlo, overallccbarlum, overallccbarextrhi, overallccbarextrlo;
384 
385 
386  //Weights are 1/sigma(stat)^2 (absolute value of error, not relative errors)
387  weightD0Dplus = 1/(D0Dplusccbarstat*D0Dplusccbarstat);
388  weightDst = 1/(ccbarDststat*ccbarDststat);
389 
390  totalweight = weightD0Dplus + weightDst;
391 
392  //Statistical error is square root of 1/(total weight):
393  overallccbarstat = 1/sqrt(totalweight);
394 
395  //Weighted averages of everything else:
396  overallccbar = (D0Dplusccbar * weightD0Dplus + ccbarDst * weightDst)/totalweight;
397  overallccbarsysthi = (D0Dplusccbarsysthi * weightD0Dplus + ccbarDstsysthi * weightDst)/totalweight;
398  overallccbarsystlo = (D0Dplusccbarsystlo * weightD0Dplus + ccbarDstsystlo * weightDst)/totalweight;
399  overallccbarlum = (D0Dplusccbarlum * weightD0Dplus + ccbarDstlum * weightDst)/totalweight;
400  overallccbarextrhi = (D0Dplusccbarextrhi * weightD0Dplus + ccbarDstextrhi * weightDst)/totalweight;
401  overallccbarextrlo = (D0Dplusccbarextrlo * weightD0Dplus + ccbarDstextrlo * weightDst)/totalweight;
402 
403 
404  //Outputs all cross sections to screen
405  cout << "D0:\nVisible =\t" << D0vis << " ± " << D0stat << " (stat.) + " << D0systhi << " - " << D0systlo << " (syst.) microbarn\nExtrapolated =\t" << crossD0 << " ± " << crossD0stat << " (stat.) +" << crossD0systhi <<" -" << crossD0systlo << " (syst.) ± " << crossD0lum << " (lum) ± " << crossD0br << " (br.) +" << crossD0extrhi << " -" << crossD0extrlo << " (extr.) microbarn\nTotal sigma =\t" << ccbarD0 << " ± " << ccbarD0stat << " (stat.) +" << ccbarD0systhi <<" -" << ccbarD0systlo << " (syst.) ± " << ccbarD0lum << " (lum) +" << ccbarD0extrhi << " -" << ccbarD0extrlo << " (extr.) microbarn\n"
406  << "Dst:\nVisible =\t" << Dstvis << " ± " << Dststat << " (stat.) + " << Dstsysthi << " - " << Dstsystlo << " (syst.) microbarn\nExtrapolated =\t" << crossDst << " ± " << crossDststat << " (stat.) +" << crossDstsysthi <<" -" << crossDstsystlo << " (syst.) ± " << crossDstlum << " (lum) ± " << crossDstbr << " (br.) +" << crossDstextrhi << " -" << crossDstextrlo << " (extr.) microbarn\nTotal sigma =\t" << ccbarDst << " ± " << ccbarDststat << " (stat.) +" << ccbarDstsysthi <<" -" << ccbarDstsystlo << " (syst.) ± " << ccbarDstlum << " (lum) +" << ccbarDstextrhi << " -" << ccbarDstextrlo << " (extr.) microbarn\n"
407  << "Dplus:\nVisible =\t" << Dplusvis << " ± " << Dplusstat << " (stat.) + " << Dplussysthi << " - " << Dplussystlo << " (syst.) microbarn\nExtrapolated =\t" << crossDplus << " ± " << crossDplusstat << " (stat.) +" << crossDplussysthi <<" -" << crossDplussystlo << " (syst.) ± " << crossDpluslum << " (lum) ± " << crossDplusbr << " (br.) +" << crossDplusextrhi << " -" << crossDplusextrlo << " (extr.) microbarn\nTotal sigma =\t" << ccbarDplus << " ± " << ccbarDplusstat << " (stat.) +" << ccbarDplussysthi <<" -" << ccbarDplussystlo << " (syst.) ± " << ccbarDpluslum << " (lum) +" << ccbarDplusextrhi << " -" << ccbarDplusextrlo << " (extr.) microbarn\n\nWeighted average =\t" << overallccbar << " ± " << overallccbarstat << " (stat.) + " << overallccbarsysthi << " - " << overallccbarsystlo << " (syst.) ± " << overallccbarlum << " (lum.) + " << overallccbarextrhi << " - " << overallccbarextrlo << " (extr.) µb\n";
408 
409 
410 
411 
412 }
Double_t
double Double_t
Definition: External.C:58
Int_t
int Int_t
Definition: External.C:63
TotalCharmCrossSection
void TotalCharmCrossSection()
Definition: TotalCharmCrossSection.C:25