PWG3/vertexingHF/macros/TotalCharmCrossSection.C

   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 }