]>
Commit | Line | Data |
---|---|---|
8c1c76e9 | 1 | #include <TF1.h> |
2 | #include <TFormula.h> | |
3 | void MakeD2eSpectra(){ | |
4 | ||
5 | // load libs | |
6 | gSystem->Load("libANALYSIS.so"); | |
7 | gSystem->Load("libANALYSISalice.so"); | |
8 | gSystem->Load("libCORRFW.so"); | |
9 | gSystem->Load("libPWG0base.so"); | |
10 | gSystem->Load("/alidata10/alice_u/minjung/032.heavy/hfe/bg4conv/util/hfe/libPWG3hfeDevel.so"); | |
11 | ||
12 | setGeneralStyle(); | |
13 | ||
14 | // read input files | |
15 | //char filename0[]="/alidata10/alice_u/minjung/train/GSI/V006.MC_pp/2011-03-08_2238.5742/LHC10f6a/HFEtask.root"; // pwg3 talk(linear binning) | |
16 | //char filename1[]="../TableOldnorm/HFPtSpectrum_D0Kpi_rebinnedth_combinedFD_021210_DecayBRSubtracted.root"; // pwg3 talk(use of 71.4mb) | |
17 | //char filename2[]="../TableOldnorm/HFPtSpectrum_DplusKpipi_combinedFD_021210_DecayBRSubtracted.root"; // pwg3 talk(use of 71.4mb) | |
18 | //char filename0[]="/alidata10/alice_u/minjung/train/GSI/V006.MC_pp/2011-03-15_2121.5954/LHC10f6a/HFEtask.root"; // preliminary plot (log binning) | |
19 | char filename0[]="./HFEtaskLHC10f6awExtendedDpt.root"; // QM plot | |
20 | char filename1[]="../TableOldnorm/HFPtSpectrum_D0Kpi_combinedFD_rebinnedth_150311_newsigma_woDecayBR.root"; // preliminary plot(use of 62.1mb) | |
21 | char filename2[]="../TableOldnorm/HFPtSpectrum_DplusKpipi_combinedFD_150311_newsigma_woDecayBR.root"; // preliminary plot(use of 62.1mb) | |
22 | ||
23 | TFile *_file[3]; | |
24 | _file[0] = TFile::Open(filename0); | |
25 | _file[1] = TFile::Open(filename1); | |
26 | _file[2] = TFile::Open(filename2); | |
27 | ||
28 | TList *tl = (TList *)_file[0]->Get("TRD_HFE_QA2")->FindObject("MCqa"); | |
29 | ||
30 | //count # of events | |
31 | TList *tl_result = (TList *)_file[0]->Get("TRD_HFE_Results2"); | |
32 | AliCFContainer *containermc = (AliCFContainer *) tl_result->FindObject("eventContainer"); | |
33 | AliCFDataGrid *mcGrid1 = new AliCFDataGrid("mcgrid1","",*containermc,AliHFEcuts::kEventStepGenerated); | |
34 | AliCFDataGrid *mcGrid2 = new AliCFDataGrid("mcgrid2","",*containermc,AliHFEcuts::kEventStepRecNoCut); | |
35 | printf("mc # of entries kEventStepGenerated= %lf kEventStepRecNoCut= %lf\n",mcGrid1->GetEntries(),mcGrid2->GetEntries()); | |
36 | ||
37 | Int_t nEvt = (Int_t) mcGrid1->GetEntries(); | |
38 | Double_t nevtnorm = 1./float(nEvt); //normalize with total number of event | |
39 | ||
40 | ||
41 | // cross section for the old D meson data | |
42 | //Double_t Xsectrion4e = 71400000000.; | |
43 | //Double_t Xsectrion4D = 71400.; // consider the data point is represented as umb | |
44 | Double_t Xsectrion4e = 62100000000.; | |
45 | Double_t Xsectrion4D = 62100.; // consider the data point is represented as umb | |
46 | ||
47 | ||
48 | TString kDType[7]; | |
49 | kDType[0]="421"; | |
50 | kDType[1]="411"; | |
51 | kDType[2]="431"; | |
52 | kDType[3]="4122"; | |
53 | kDType[4]="4132"; | |
54 | kDType[5]="4232"; | |
55 | kDType[6]="4332"; | |
56 | ||
57 | TString kDType2[4]; | |
58 | kDType2[0]=""; | |
59 | kDType2[1]="D0"; | |
60 | kDType2[2]="Dp"; | |
61 | kDType2[3]="Drest"; | |
62 | ||
63 | TString hname; | |
64 | ||
65 | // 2D histo having D, e pt correlation | |
66 | TH2D *hDeCorr[4]; // 4 groups. "", "D0", "Dp", "Drest" | |
67 | for(int i=0; i<4; i++){ | |
68 | hname = "mcqa_barrel_PtCorr"+kDType2[i]+"_c"; | |
69 | hDeCorr[i]=(TH2D*)tl->FindObject(hname); | |
70 | } | |
71 | ||
72 | // D mesons' pt, y 2D histograms from PYTHIA | |
73 | TH2D *hDPYTHIA[7]; // 7 particles. 421:D0, 411:D+, 431:Ds, 4122, 4132, 4232, 4332 | |
74 | for(int i=0; i<7; i++){ | |
75 | hname = "Dmeson"+kDType[i]; | |
76 | hDPYTHIA[i]=(TH2D*)tl->FindObject(hname); | |
77 | } | |
78 | ||
79 | // PYTHIA pt distribution in |y|<0.5 of 7 D mesons | |
80 | TH1D *hDPYTHIApt[7]; | |
81 | for(int i=0; i<7; i++){ | |
82 | hname= "hD"+kDType[i]+"PYTHIApt"; | |
83 | hDPYTHIApt[i]=(TH1D*)hDPYTHIA[i]->ProjectionX(hname,5,5); // |rapidity|<0.5 | |
84 | ||
85 | CorrectFromTheWidth(hDPYTHIApt[i]); // bin width correction | |
86 | hDPYTHIApt[i]->Scale(nevtnorm); // # of event normalization | |
87 | hDPYTHIApt[i]->Scale(0.5); // (particle + antiparticle)/2 | |
88 | hDPYTHIApt[i]->Scale(Xsectrion4D); // convert to cross section | |
89 | hDPYTHIApt[i]->SetMarkerStyle(20+i); | |
90 | hDPYTHIApt[i]->GetYaxis()->SetTitle("d#sigma/dp_{t} |_{|y|<0.5} [#mu b/GeV/c]"); | |
91 | } | |
92 | ||
93 | // D measured (be careful with binnings) | |
94 | TH1D *hD0MeasuredStat= (TH1D *)_file[1]->Get("histoSigmaCorr"); // data & statistic errors, start from bin number 4 | |
95 | TGraphAsymmErrors *gD0measuredSys= (TGraphAsymmErrors *)_file[1]->Get("gSigmaCorr"); // data & systematic errors, start from bin number 4 | |
96 | ||
97 | TH1D *hDpMeasuredStat = (TH1D *)_file[2]->Get("histoSigmaCorr"); // data & statistic errors, start from bin number 1 | |
98 | TGraphAsymmErrors *gDpMeasuredSys = (TGraphAsymmErrors *)_file[2]->Get("gSigmaCorr"); // data & statistic errors, start from bin number 1 | |
99 | ||
100 | // drawing ------------------------------------------------------- | |
101 | cPtDmeson = new TCanvas("cPtDmeson","pT of D mesons",0,0,800,500); | |
102 | cPtDmeson->Divide(2,1); | |
103 | ||
104 | cPtDmeson->cd(1); | |
105 | hDPYTHIApt[0]->Draw(""); //0: D0 | |
106 | hD0MeasuredStat->Draw("same"); | |
107 | TLegend *legD0 = new TLegend(0.60,0.75,0.75,0.85); | |
108 | legD0->AddEntry(hDPYTHIApt[0],"D^{0} PYTHIA","p"); | |
109 | legD0->AddEntry(hD0MeasuredStat,"D^{0} Measured","p"); | |
110 | setLegendStyle(*legD0,0); | |
111 | legD0->Draw("same"); | |
112 | ||
113 | cPtDmeson->cd(2); | |
114 | hDPYTHIApt[1]->Draw(""); //1: Dp | |
115 | hDpMeasuredStat->Draw("same"); | |
116 | TLegend *legDp = new TLegend(0.60,0.75,0.75,0.85); | |
117 | legDp->AddEntry(hDPYTHIApt[1],"D^{+} PYTHIA","p"); | |
118 | legDp->AddEntry(hDpMeasuredStat,"D^{+} Measured","p"); | |
119 | setLegendStyle(*legDp,0); | |
120 | legDp->Draw("same"); | |
121 | //---------------------------------------------------------------- | |
122 | ||
123 | ||
124 | Double_t ptval[6]; | |
125 | Double_t XsecD0[6], XsecD0ErrStatFrac[6], XsecD0ErrSysFracLow[6], XsecD0ErrSysFracHigh[6]; | |
126 | Double_t XsecDp[6], XsecDpErrStatFrac[6], XsecDpErrSysFracLow[6], XsecDpErrSysFracHigh[6]; | |
127 | Double_t iM2PD0[6], iM2PDp[6]; | |
128 | ||
129 | Double_t xbins[12]={0.,0.5,1.,2.,3.,4.,5.,6.,8.,12.,16.,20.}; | |
130 | TH1D *hM2PD0= new TH1D("hM2PD0","D0: Measured/Pythia",11,xbins); | |
131 | TH1D *hM2PDp= new TH1D("hM2PDp","Dp: Measured/Pythia",11,xbins); | |
132 | ||
133 | for(int i=0; i<6; i++){ // 0 -> first measured point, there are 6 points for pp for the moment | |
134 | // D0 | |
135 | iM2PD0[i]=0.; | |
136 | ptval[i]=hD0MeasuredStat->GetBinCenter(i+4); // +4: starting from the first measured point | |
137 | XsecD0[i]=hD0MeasuredStat->GetBinContent(i+4); | |
138 | if(hDPYTHIApt[0]->GetBinContent(i+4)>0) iM2PD0[i]=hD0MeasuredStat->GetBinContent(i+4)/hDPYTHIApt[0]->GetBinContent(i+4); | |
139 | if(hD0MeasuredStat->GetBinContent(i+4)>0){ | |
140 | XsecD0ErrStatFrac[i] = hD0MeasuredStat->GetBinError(i+4)/hD0MeasuredStat->GetBinContent(i+4); | |
141 | XsecD0ErrSysFracLow[i] = gD0measuredSys->GetErrorYlow(i+4)/hD0MeasuredStat->GetBinContent(i+4); | |
142 | XsecD0ErrSysFracHigh[i] = gD0measuredSys->GetErrorYhigh(i+4)/hD0MeasuredStat->GetBinContent(i+4); | |
143 | } | |
144 | ||
145 | // D+ | |
146 | iM2PDp[i]=0.; | |
147 | XsecDp[i]=hDpMeasuredStat->GetBinContent(i+1); // be careful with the fact that D+ measured starting from bin 1 | |
148 | if(hDPYTHIApt[1]->GetBinContent(i+4)>0) iM2PDp[i]=hDpMeasuredStat->GetBinContent(i+1)/hDPYTHIApt[1]->GetBinContent(i+4); //D+ | |
149 | if(hDpMeasuredStat->GetBinContent(i+1)>0) { | |
150 | XsecDpErrStatFrac[i] = hDpMeasuredStat->GetBinError(i+1)/hDpMeasuredStat->GetBinContent(i+1); | |
151 | XsecDpErrSysFracLow[i] = gDpMeasuredSys->GetErrorYlow(i+1)/hDpMeasuredStat->GetBinContent(i+1); | |
152 | XsecDpErrSysFracHigh[i] = gDpMeasuredSys->GetErrorYhigh(i+1)/hDpMeasuredStat->GetBinContent(i+1); | |
153 | } | |
154 | ||
155 | hM2PD0->Fill(ptval[i],iM2PD0[i]); | |
156 | hM2PDp->Fill(ptval[i],iM2PDp[i]); | |
157 | } | |
158 | // drawing ------------------------------------------------------- | |
159 | cM2PDmeson = new TCanvas("cM2PDmeson","Measured/Pythia",0,0,800,500); | |
160 | cM2PDmeson->Divide(2,1); | |
161 | cM2PDmeson->cd(1); | |
162 | hM2PD0->SetMarkerStyle(20); | |
163 | hM2PD0->SetXTitle("p_{t} [GeV/c]"); | |
164 | hM2PD0->SetYTitle("ratio"); | |
165 | hM2PD0->Draw("p"); | |
166 | TLegend *legM2PD0 = new TLegend(0.45,0.75,0.75,0.85); | |
167 | legM2PD0->AddEntry(hM2PD0,"D0: Measured/PYTHIA","p"); | |
168 | setLegendStyle(*legM2PD0,0); | |
169 | legM2PD0->Draw("same"); | |
170 | ||
171 | cM2PDmeson->cd(2); | |
172 | hM2PDp->SetMarkerStyle(20); | |
173 | hM2PDp->SetXTitle("p_{t} [GeV/c]"); | |
174 | hM2PDp->SetYTitle("ratio"); | |
175 | hM2PDp->Draw("p"); | |
176 | TLegend *legM2PDp = new TLegend(0.45,0.75,0.75,0.85); | |
177 | legM2PDp->AddEntry(hM2PDp,"D+: Measured/PYTHIA","p"); | |
178 | setLegendStyle(*legM2PDp,0); | |
179 | legM2PDp->Draw("same"); | |
180 | //---------------------------------------------------------------- | |
181 | ||
182 | ||
183 | // electrons pt spectra from given D meson pt bins (before any scaling) | |
184 | TString kDPtbin[14]; | |
185 | kDPtbin[0]="Dpt0"; | |
186 | kDPtbin[1]="Dpt05"; | |
187 | kDPtbin[2]="Dpt1"; | |
188 | kDPtbin[3]="Dpt2"; | |
189 | kDPtbin[4]="Dpt3"; | |
190 | kDPtbin[5]="Dpt4"; | |
191 | kDPtbin[6]="Dpt5"; | |
192 | kDPtbin[7]="Dpt6"; | |
193 | kDPtbin[8]="Dpt8"; | |
194 | kDPtbin[9]="Dpt12"; | |
195 | kDPtbin[10]="Dpt16"; | |
196 | kDPtbin[11]="Dpt20"; | |
197 | kDPtbin[12]="Dpt30"; | |
198 | kDPtbin[13]="Dpt40"; | |
199 | ||
200 | Int_t kbinl[14]={1, 6,11,21,31,41,51,61, 81,121,161,201,301,401}; | |
201 | Int_t kbinh[14]={5,10,20,30,40,50,60,80,120,160,200,300,400,500}; | |
202 | ||
203 | TH1D *hD2e[4][14]; | |
204 | TH1D *hD2eNorm[4][14]; | |
205 | TH1D *hD2eNormWm[2][14]; | |
206 | Int_t nRebin = 1; | |
207 | for(int iDtype=0; iDtype<2; iDtype++){ // 4 groups. 0:"", 1:"D0", 2:"Dp", 3:"Drest" | |
208 | for(int iDptbin=0; iDptbin<14; iDptbin++){ | |
209 | hname = "h"+kDType2[iDtype+1]+kDPtbin[iDptbin]; | |
210 | hD2e[iDtype][iDptbin] = (TH1D*)hDeCorr[iDtype+1]->ProjectionY(hname,kbinl[iDptbin],kbinh[iDptbin]); | |
211 | hD2e[iDtype][iDptbin]->Rebin(nRebin); | |
212 | CorrectFromTheWidth(hD2e[iDtype][iDptbin]); | |
213 | hname = "hNorm"+kDType2[iDtype]+kDPtbin[iDptbin]; | |
214 | hD2eNorm[iDtype][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); | |
215 | hD2eNormWm[iDtype][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); | |
216 | if(iDtype==0) hD2eNorm[2][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); // for Ds, use D0 electron spectra | |
217 | if(iDtype==0) hD2eNorm[3][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); // for Lc, use D0 electron spectra | |
218 | } | |
219 | } | |
220 | ||
221 | hD0pt0to2=(TH1D*)hD2e[0][0]->Clone("hD0pt0to2"); | |
222 | hD0pt2to12=(TH1D*)hD2e[0][3]->Clone("hD0pt2to12"); | |
223 | hD0pt12to50=(TH1D*)hD2e[0][9]->Clone("hD0pt12to50"); | |
224 | hD0pt0to2->Add(hD2e[0][1]); | |
225 | hD0pt0to2->Add(hD2e[0][2]); | |
226 | hD0pt2to12->Add(hD2e[0][4]); | |
227 | hD0pt2to12->Add(hD2e[0][5]); | |
228 | hD0pt2to12->Add(hD2e[0][6]); | |
229 | hD0pt2to12->Add(hD2e[0][7]); | |
230 | hD0pt2to12->Add(hD2e[0][8]); | |
231 | hD0pt12to50->Add(hD2e[0][10]); | |
232 | hD0pt12to50->Add(hD2e[0][11]); | |
233 | hD0pt12to50->Add(hD2e[0][12]); | |
234 | hD0pt12to50->Add(hD2e[0][13]); | |
235 | hD0ptall=(TH1D*)hD0pt0to2->Clone("hD0ptall"); | |
236 | hD0ptall->Add(hD0pt2to12); | |
237 | hD0ptall->Add(hD0pt12to50); | |
238 | ||
239 | hD0lowptcontrib=(TH1D*)hD0pt0to2->Clone("hD0lowptcontrib"); | |
240 | hD0lowptcontrib->Divide(hD0pt2to12); | |
241 | ||
242 | hD0pt0to2->Divide(hD0ptall); | |
243 | hD0pt2to12->Divide(hD0ptall); | |
244 | hD0pt12to50->Divide(hD0ptall); | |
245 | ||
246 | new TCanvas; | |
247 | hD0pt0to2->SetMarkerStyle(24); | |
248 | hD0pt2to12->SetMarkerStyle(24); | |
249 | hD0pt12to50->SetMarkerStyle(24); | |
250 | hD0pt0to2->SetMarkerColor(1); | |
251 | hD0pt0to2->SetLineColor(1); | |
252 | hD0pt2to12->SetMarkerColor(4); | |
253 | hD0pt2to12->SetLineColor(4); | |
254 | hD0pt12to50->SetMarkerColor(2); | |
255 | hD0pt12to50->SetLineColor(2); | |
256 | hD0pt0to2->Draw(); | |
257 | hD0pt2to12->Draw("same"); | |
258 | hD0pt12to50->Draw("same"); | |
259 | gPad->SetLogy(); | |
260 | gPad->SetGrid(); | |
261 | ||
262 | // drawing ------------------------------------------------------- | |
263 | cPtD2e = new TCanvas("cPtD2e","pT of e from certain D pt bin",0,0,1000,700); | |
264 | cPtD2e->Divide(2,1); | |
265 | Int_t colorcode=1; | |
266 | TLegend *legD2e[4]; | |
267 | for(int iDtype=0; iDtype<2; iDtype++){ | |
268 | legD2e[iDtype] = new TLegend(0.75,0.45,0.88,0.88); | |
269 | legD2e[iDtype]->AddEntry(hD2e[0][0],kDType2[iDtype+1],""); | |
270 | for(int iDptbin=0; iDptbin<14; iDptbin++){ | |
271 | cPtD2e->cd(iDtype+1); | |
272 | colorcode=iDptbin+1; | |
273 | if(colorcode==10) colorcode=38; | |
274 | hD2e[iDtype][iDptbin]->SetLineColor(colorcode); | |
275 | hD2e[iDtype][iDptbin]->SetMarkerColor(colorcode); | |
276 | hD2e[iDtype][iDptbin]->SetMarkerStyle(20); | |
277 | if(iDptbin==0) hD2e[iDtype][iDptbin]->Draw(); | |
278 | else hD2e[iDtype][iDptbin]->Draw("samep"); | |
279 | gPad->SetLogy(); | |
280 | gPad->SetLogx(); | |
281 | legD2e[iDtype]->AddEntry(hD2e[iDtype][iDptbin],kDPtbin[iDptbin],"p"); | |
282 | setLegendStyle(*legD2e[iDtype],0); | |
283 | } | |
284 | legD2e[iDtype]->Draw("same"); | |
285 | } | |
286 | //---------------------------------------------------------------- | |
287 | ||
288 | // electrons pt spectra from given D meson pt bins "with scaling based on measured X section" | |
289 | Double_t accfactorD[4][6] = { {1.73018, 1.72687, 1.76262, 1.94291, 1.82095, 1.69734,}, // acceptance factor of D0->e | |
290 | {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054,}, // acceptance factor of Dp->e | |
291 | {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054,}, // acceptance factor of Ds->e(copy of Dp) | |
292 | {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054} }; // acceptance factor of Lc->e(copy of Dp) | |
293 | Double_t ptbwdth[6] = {1., 1., 1., 1., 2., 4.}; // pt bin width for binwidth correction | |
294 | ||
295 | // normalize PYTHIA electron spectra with the measured D meson cross section * branching ratio | |
296 | Double_t brachratioD[4]; | |
297 | brachratioD[0] = 0.0653; // D0->e branching ratio | |
298 | brachratioD[1] = 0.16; // Dp->e branching ratio | |
299 | brachratioD[2] = 0.08; // Ds->e branching ratio | |
300 | brachratioD[3] = 0.045; // Lc->e branching ratio | |
301 | //Double_t DiffMthdScaleD0= 1./1.07; // considering difference between weighting method and X section method for D0 : new | |
302 | //Double_t DiffMthdScaleDp= 1./1.19; // considering difference between weighting method and X section method for Dp : new | |
303 | Double_t DiffMthdScaleD0= 1./1.12; // considering difference between weighting method and X section method for D0 | |
304 | Double_t DiffMthdScaleDp= 1./1.25; // considering difference between weighting method and X section method for Dp | |
305 | ||
306 | ||
307 | Double_t norm, xsec; | |
308 | Double_t xsecD[4][6]; | |
309 | for(int iptb=0; iptb<6; iptb++){ | |
310 | xsecD[0][iptb] = XsecD0[iptb]; //D0 | |
311 | xsecD[1][iptb] = XsecDp[iptb]; //Dp | |
312 | xsecD[2][iptb] = (XsecD0[iptb]+XsecDp[iptb])*2.37/(8.49+2.65+5.07); //Ds: use of ZEUS ratio | |
313 | xsecD[3][iptb] = (XsecD0[iptb]+XsecDp[iptb])*3.59/(8.49+2.65+5.07); //Lc: use of ZEUS ratio | |
314 | } | |
315 | ||
316 | for(int iDtype=0; iDtype<4; iDtype++){ | |
317 | for(int iptb=0; iptb<6; iptb++){ | |
318 | xsec = xsecD[iDtype][iptb]*1000000.; // change scale from ub to pb | |
319 | norm = hD2eNorm[iDtype][iptb+3]->Integral("width")/(xsec*brachratioD[iDtype]*ptbwdth[iptb]); //D0 | |
320 | norm = accfactorD[iDtype][iptb]/norm; | |
321 | hD2eNorm[iDtype][iptb+3]->Scale(norm); | |
322 | } | |
323 | } | |
324 | ||
325 | // not measured pt data points for D0 | |
326 | hD2eNorm[0][9]->Scale(nevtnorm); | |
327 | hD2eNorm[0][9]->Scale(0.44); // measured/pythia for this pt bin | |
328 | hD2eNorm[0][9]->Scale(0.5); | |
329 | hD2eNorm[0][9]->Scale(DiffMthdScaleD0); | |
330 | hD2eNorm[0][9]->Scale(Xsectrion4e); | |
331 | hD2eNorm[0][10]->Scale(nevtnorm); | |
332 | hD2eNorm[0][10]->Scale(0.44); // measured/pythia for this pt bin | |
333 | hD2eNorm[0][10]->Scale(0.5); | |
334 | hD2eNorm[0][10]->Scale(DiffMthdScaleD0); | |
335 | hD2eNorm[0][10]->Scale(Xsectrion4e); | |
336 | hD2eNorm[0][11]->Scale(nevtnorm); | |
337 | hD2eNorm[0][11]->Scale(0.44); // measured/pythia for this pt bin | |
338 | hD2eNorm[0][11]->Scale(0.5); | |
339 | hD2eNorm[0][11]->Scale(DiffMthdScaleD0); | |
340 | hD2eNorm[0][11]->Scale(Xsectrion4e); | |
341 | hD2eNorm[0][12]->Scale(nevtnorm); | |
342 | hD2eNorm[0][12]->Scale(0.44); // measured/pythia for this pt bin | |
343 | hD2eNorm[0][12]->Scale(0.5); | |
344 | hD2eNorm[0][12]->Scale(DiffMthdScaleD0); | |
345 | hD2eNorm[0][12]->Scale(Xsectrion4e); | |
346 | hD2eNorm[0][13]->Scale(nevtnorm); | |
347 | hD2eNorm[0][13]->Scale(0.44); // measured/pythia for this pt bin | |
348 | hD2eNorm[0][13]->Scale(0.5); | |
349 | hD2eNorm[0][13]->Scale(DiffMthdScaleD0); | |
350 | hD2eNorm[0][13]->Scale(Xsectrion4e); | |
351 | ||
352 | ||
353 | // not measured pt data points for D+ | |
354 | hD2eNorm[1][9]->Scale(nevtnorm); | |
355 | hD2eNorm[1][9]->Scale(0.41); // measured/pythia for this pt bin | |
356 | hD2eNorm[1][9]->Scale(0.5); | |
357 | hD2eNorm[1][9]->Scale(DiffMthdScaleDp); | |
358 | hD2eNorm[1][9]->Scale(Xsectrion4e); | |
359 | hD2eNorm[1][10]->Scale(nevtnorm); | |
360 | hD2eNorm[1][10]->Scale(0.41); // measured/pythia for this pt bin | |
361 | hD2eNorm[1][10]->Scale(0.5); | |
362 | hD2eNorm[1][10]->Scale(DiffMthdScaleDp); | |
363 | hD2eNorm[1][10]->Scale(Xsectrion4e); | |
364 | hD2eNorm[1][11]->Scale(nevtnorm); | |
365 | hD2eNorm[1][11]->Scale(0.41); // measured/pythia for this pt bin | |
366 | hD2eNorm[1][11]->Scale(0.5); | |
367 | hD2eNorm[1][11]->Scale(DiffMthdScaleDp); | |
368 | hD2eNorm[1][11]->Scale(Xsectrion4e); | |
369 | hD2eNorm[1][12]->Scale(nevtnorm); | |
370 | hD2eNorm[1][12]->Scale(0.41); // measured/pythia for this pt bin | |
371 | hD2eNorm[1][12]->Scale(0.5); | |
372 | hD2eNorm[1][12]->Scale(DiffMthdScaleDp); | |
373 | hD2eNorm[1][12]->Scale(Xsectrion4e); | |
374 | hD2eNorm[1][13]->Scale(nevtnorm); | |
375 | hD2eNorm[1][13]->Scale(0.41); // measured/pythia for this pt bin | |
376 | hD2eNorm[1][13]->Scale(0.5); | |
377 | hD2eNorm[1][13]->Scale(DiffMthdScaleDp); | |
378 | hD2eNorm[1][13]->Scale(Xsectrion4e); | |
379 | ||
380 | ||
381 | // estimation for Ds and Lc based on D0, D+ at high pt | |
382 | Double_t fragRatio=2.37/(8.49+2.65+5.07); | |
383 | Double_t brachRatio[2]; | |
384 | Double_t brachRatio[0]=2*brachratioD[2]/(brachratioD[0]+brachratioD[1]); | |
385 | Double_t brachRatio[1]=2*brachratioD[3]/(brachratioD[0]+brachratioD[1]); | |
386 | for(int iptb=9; iptb<14; iptb++){ | |
387 | hD2eNorm[2][iptb]->Add(hD2eNorm[0][iptb]); | |
388 | hD2eNorm[2][iptb]->Add(hD2eNorm[1][iptb]); | |
389 | hD2eNorm[2][iptb]->Scale(fragRatio); | |
390 | hD2eNorm[2][iptb]->Scale(brachRatio[0]); | |
391 | hD2eNorm[3][iptb]->Add(hD2eNorm[0][iptb]); | |
392 | hD2eNorm[3][iptb]->Add(hD2eNorm[1][iptb]); | |
393 | hD2eNorm[3][iptb]->Scale(fragRatio); | |
394 | hD2eNorm[3][iptb]->Scale(brachRatio[1]); | |
395 | } | |
396 | ||
397 | ||
398 | int kptmaxb = 11; // D0, D+ | |
399 | TH1D *hD2eNormSum[4]; | |
400 | for(int iDtype=0; iDtype<4; iDtype++){ | |
401 | hname="hD2eNormSum" + iDtype; | |
402 | hD2eNormSum[iDtype]=(TH1D*)hD2eNorm[iDtype][3]->Clone(hname); | |
403 | // if(iDtype==2 || iDtype==3) kptmaxb=6; //Ds, Lc //mmjj | |
404 | for(int iptb=1; iptb<kptmaxb; iptb++){ | |
405 | hD2eNormSum[iDtype]->Add(hD2eNorm[iDtype][iptb+3]); | |
406 | } | |
407 | } | |
408 | ||
409 | for(int i=1; i<hD2eNormSum[0]->FindBin(2.0); i++){ // consider low pt contribution | |
410 | for(int iDtype=0; iDtype<4; iDtype++){ | |
411 | Double_t icontrib =hD0lowptcontrib->GetBinContent(i)*hD2eNormSum[iDtype]->GetBinContent(i); | |
412 | Double_t ipt =hD2eNormSum[iDtype]->GetBinCenter(i); | |
413 | hD2eNormSum[iDtype]->Fill(ipt, icontrib); | |
414 | } | |
415 | } | |
416 | ||
417 | ||
418 | TH1D *hD2eNormSumStat[4]; // copy of hD2eNormSum[4] to get mc statistical errors | |
419 | for(int iDtype=0; iDtype<4; iDtype++){ | |
420 | hname="hD2eNormSumStat" + iDtype; | |
421 | hD2eNormSumStat[iDtype]=(TH1D*)hD2eNormSum[iDtype]->Clone(hname); | |
422 | } | |
423 | ||
424 | // sum | |
425 | hD2eNormTotalStat=(TH1D*)hD2eNormSumStat[0]->Clone("hD2eNormTotalStat"); | |
426 | hD2eNormTotalStat->Add(hD2eNormSumStat[1]); | |
427 | hD2eNormTotalStat->Add(hD2eNormSumStat[2]); | |
428 | hD2eNormTotalStat->Add(hD2eNormSumStat[3]); | |
429 | ||
430 | // drawing, here the error bars are only for MC statistics --------------------------- | |
431 | cPtD2eSum = new TCanvas("cPtD2eSum","pT of e from D",0,0,1000,600); | |
432 | cPtD2eSum->Divide(2,1); | |
433 | cPtD2eSum->cd(1); | |
434 | hD2eNormSumStat[0]->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
435 | hD2eNormSumStat[0]->SetMarkerStyle(20); | |
436 | hD2eNormSumStat[0]->Draw(); | |
437 | hD2eNormSumStat[1]->SetMarkerColor(2); | |
438 | hD2eNormSumStat[1]->SetMarkerStyle(20); | |
439 | hD2eNormSumStat[1]->Draw("same"); | |
440 | hD2eNormSumStat[2]->SetMarkerColor(4); | |
441 | hD2eNormSumStat[2]->SetMarkerStyle(20); | |
442 | hD2eNormSumStat[2]->Draw("same"); | |
443 | hD2eNormSumStat[3]->SetMarkerColor(3); | |
444 | hD2eNormSumStat[3]->SetMarkerStyle(20); | |
445 | hD2eNormSumStat[3]->Draw("same"); | |
446 | gPad->SetLogy(); | |
447 | TLegend *legsum = new TLegend(0.25,0.70,0.75,0.85); | |
448 | legsum->AddEntry(hD2eNormSumStat[0],"D^{0}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
449 | legsum->AddEntry(hD2eNormSumStat[1],"D^{+}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
450 | legsum->AddEntry(hD2eNormSumStat[2],"D_{s}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
451 | legsum->AddEntry(hD2eNormSumStat[3],"#Lambda_{c}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
452 | setLegendStyle(*legsum,0); legsum->Draw("same"); | |
453 | ||
454 | cPtD2eSum->cd(2); | |
455 | hD2eNormTotalStat->SetMarkerStyle(20); | |
456 | hD2eNormTotalStat->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
457 | hD2eNormTotalStat->Draw(); | |
458 | TLegend *legtotal = new TLegend(0.30,0.75,0.75,0.85); | |
459 | legtotal->AddEntry(hD2eNormTotalStat,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
460 | setLegendStyle(*legtotal,0); | |
461 | legtotal->Draw("same"); | |
462 | gPad->SetLogy(); | |
463 | //---------------------------------------------------------------- | |
464 | ||
465 | ||
466 | // calculating propagated error of D meson's statistic and systematic errors | |
467 | TGraphAsymmErrors *gD2eNormSum[4]; // store propagated error from D systematic errors | |
468 | for(iDtype=0; iDtype<4; iDtype++){ | |
469 | gD2eNormSum[iDtype]= new TGraphAsymmErrors(hD2eNormSum[iDtype]); | |
470 | } | |
471 | ||
472 | Double_t statErr[11]; //kptmaxb | |
473 | Double_t sysErrLow[11]; //kptmaxb | |
474 | Double_t sysErrHigh[11]; //kptmaxb | |
475 | Double_t statErrSum[4]; // D0, D+, Ds. Lc | |
476 | Double_t sysErrLowSum[4]; // D0, D+, Ds. Lc | |
477 | Double_t sysErrHighSum[4]; // D0, D+, Ds. Lc | |
478 | Int_t ib4err = 0; | |
479 | for(int i=0; i<hD2eNormSum[0]->GetNbinsX(); i++){ | |
480 | for(int iDtype=0; iDtype<4; iDtype++){ | |
481 | statErrSum[iDtype]=0.; | |
482 | sysErrLowSum[iDtype]=0.; | |
483 | sysErrHighSum[iDtype]=0.; | |
484 | } | |
485 | ||
486 | for(int iptb=0; iptb<11; iptb++){ | |
487 | // calculate propagated error from D statistical errors | |
488 | if(iptb<6) ib4err = iptb; | |
489 | else ib4err = 5; // last two bins from pythia, so take the error from the last measured bin | |
490 | // statistical error of D propagated to systematic error of elec | |
491 | statErr[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrStatFrac[ib4err]; //D0 | |
492 | if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5; //1.5 put conservative error | |
493 | statErrSum[0] += statErr[iptb]*statErr[iptb]; | |
494 | statErr[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]; //Dp | |
495 | if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5; | |
496 | statErrSum[1] += statErr[iptb]*statErr[iptb]; | |
497 | statErr[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]*1.5; //Ds | |
498 | if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5; | |
499 | statErrSum[2] += statErr[iptb]*statErr[iptb]; | |
500 | statErr[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]*1.5; //Lc | |
501 | if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5; | |
502 | statErrSum[3] += statErr[iptb]*statErr[iptb]; | |
503 | ||
504 | // systematic error of D propagated to systematic error of elec | |
505 | sysErrLow[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrSysFracLow[ib4err]; //D0 | |
506 | if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5; | |
507 | sysErrLowSum[0] += sysErrLow[iptb]*sysErrLow[iptb]; | |
508 | sysErrLow[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]; //Dp | |
509 | if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5; | |
510 | sysErrLowSum[1] += sysErrLow[iptb]*sysErrLow[iptb]; | |
511 | sysErrLow[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]*1.5; //Ds | |
512 | if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5; | |
513 | sysErrLowSum[2] += sysErrLow[iptb]*sysErrLow[iptb]; | |
514 | sysErrLow[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]*1.5; //Lc | |
515 | if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5; | |
516 | sysErrLowSum[3] += sysErrLow[iptb]*sysErrLow[iptb]; | |
517 | ||
518 | sysErrHigh[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrSysFracHigh[ib4err]; //D0 | |
519 | if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5; | |
520 | sysErrHighSum[0] += sysErrHigh[iptb]*sysErrHigh[iptb]; | |
521 | sysErrHigh[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]; //Dp | |
522 | if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5; | |
523 | sysErrHighSum[1] += sysErrHigh[iptb]*sysErrHigh[iptb]; | |
524 | sysErrHigh[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]*1.5; //Ds | |
525 | if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5; | |
526 | sysErrHighSum[2] += sysErrHigh[iptb]*sysErrHigh[iptb]; | |
527 | sysErrHigh[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]*1.5; //Lc | |
528 | if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5; | |
529 | sysErrHighSum[3] += sysErrHigh[iptb]*sysErrHigh[iptb]; | |
530 | } | |
531 | ||
532 | for(iDtype=0; iDtype<4; iDtype++){ | |
533 | hD2eNormSum[iDtype]->SetBinError(i+1,TMath::Sqrt(statErrSum[iDtype])); | |
534 | gD2eNormSum[iDtype]->SetPointError(i,0,0,TMath::Sqrt(sysErrLowSum[iDtype]),TMath::Sqrt(sysErrHighSum[iDtype])); | |
535 | } | |
536 | } | |
537 | ||
538 | // sum ----------------------------------------------------------- | |
539 | hD2eNormTotal=(TH1D*)hD2eNormSum[0]->Clone("hD2eNormTotal"); | |
540 | hD2eNormTotal->Add(hD2eNormSum[1]); | |
541 | hD2eNormTotal->Add(hD2eNormSum[2]); | |
542 | hD2eNormTotal->Add(hD2eNormSum[3]); | |
543 | ||
544 | // sum up errors ------------------------------------------------- | |
545 | TGraphAsymmErrors *gD2eNormTotal = new TGraphAsymmErrors(hD2eNormTotal); | |
546 | TGraphAsymmErrors *gD2eNormTotalwMCstat = new TGraphAsymmErrors(hD2eNormTotal); | |
547 | gD2eNormTotal->SetName("gD2eSys"); | |
548 | gD2eNormTotal->SetTitle("electron spectra based on D measurement with systematic errors"); | |
549 | gD2eNormTotalwMCstat->SetName("gD2eSysStat"); | |
550 | gD2eNormTotalwMCstat->SetTitle("electron spectra based on D measurement with systematic and statistic errors"); | |
551 | ||
552 | Double_t errLow[4], errHigh[4], errDStat[4], errMCStat[4]; | |
553 | Double_t errLowSum, errHighSum, errDStatSum, errMCStatSum; | |
554 | for(int i=0; i<hD2eNormTotal->GetNbinsX(); i++){ | |
555 | errLowSum = 0.; | |
556 | errHighSum = 0.; | |
557 | errDStatSum = 0.; | |
558 | errMCStatSum = 0.; | |
559 | ||
560 | for(int iDtype=0; iDtype<2; iDtype++){ | |
561 | errLow[iDtype] = gD2eNormSum[iDtype]->GetErrorYlow(i); // error from D systematic error | |
562 | errLowSum += errLow[iDtype]; | |
563 | ||
564 | errHigh[iDtype]= gD2eNormSum[iDtype]->GetErrorYhigh(i); // error from D systematic error | |
565 | errHighSum += errHigh[iDtype]; | |
566 | ||
567 | errDStat[iDtype] = hD2eNormSum[iDtype]->GetBinError(i+1); //error from D statistical error | |
568 | errDStatSum += errDStat[iDtype]; | |
569 | ||
570 | errMCStat[iDtype] = hD2eNormSumStat[iDtype]->GetBinError(i+1); //error from MC statistical error | |
571 | errMCStatSum += errMCStat[iDtype]; | |
572 | } | |
573 | if(i<hD2eNormTotal->FindBin(2.0)){ | |
574 | errDStatSum = (hD0lowptcontrib->GetBinContent(i+1)+1.)*errDStatSum; //put conservative error based on the yield fraction of low pt contribution | |
575 | } | |
576 | ||
577 | gD2eNormTotal->SetPointError(i,0,0,TMath::Sqrt(errLowSum*errLowSum + errDStatSum*errDStatSum),TMath::Sqrt(errHighSum*errHighSum + errDStatSum*errDStatSum)); | |
578 | gD2eNormTotalwMCstat->SetPointError(i,0,0,TMath::Sqrt(errLowSum*errLowSum + errDStatSum*errDStatSum + errMCStatSum*errMCStatSum),TMath::Sqrt(errHighSum*errHighSum + errDStatSum*errDStatSum + errMCStatSum*errMCStatSum)); | |
579 | ||
580 | } | |
581 | ||
582 | // drawing ---------------------------------------------------------- | |
583 | cPtD2eSumwErr = new TCanvas("cPtD2eSumwErr","pT of e from D",0,0,1000,600); | |
584 | cPtD2eSumwErr->cd(1); | |
585 | gD2eNormTotalwMCstat->SetLineColor(2); | |
586 | gD2eNormTotalwMCstat->SetMarkerColor(2); | |
587 | gD2eNormTotalwMCstat->SetMarkerStyle(20); | |
588 | gD2eNormTotalwMCstat->Draw("ACP"); | |
589 | gD2eNormTotalwMCstat->GetXaxis()->SetTitle("p_{t} [GeV/c]"); | |
590 | gD2eNormTotalwMCstat->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
591 | gD2eNormTotal->SetMarkerStyle(20); | |
592 | gD2eNormTotal->SetLineStyle(2); | |
593 | gD2eNormTotal->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
594 | gD2eNormTotal->Draw("Psame"); | |
595 | TLegend *legtotalerr = new TLegend(0.50,0.75,0.75,0.85); | |
596 | legtotalerr->AddEntry(gD2eNormTotal,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","lp"); | |
597 | legtotalerr->AddEntry(gD2eNormTotalwMCstat,"(include MC stat error)","l"); | |
598 | setLegendStyle(*legtotalerr,0); | |
599 | legtotalerr->Draw("same"); | |
600 | gPad->SetLogy(); | |
601 | // ------------------------------------------------------------------ | |
602 | ||
603 | // weighting method | |
604 | Double_t iM2PD[6]; | |
605 | for(int iDtype=0; iDtype<2; iDtype++){ | |
606 | for(int iptb=0; iptb<6; iptb++){ | |
607 | if(iDtype==0) iM2PD[iptb]=iM2PD0[iptb]; | |
608 | if(iDtype==1) iM2PD[iptb]=iM2PDp[iptb]; | |
609 | hD2eNormWm[iDtype][iptb+3]->Scale(nevtnorm); | |
610 | hD2eNormWm[iDtype][iptb+3]->Scale(iM2PD[iptb]); // measured/pythia for this pt bin | |
611 | hD2eNormWm[iDtype][iptb+3]->Scale(0.5); | |
612 | hD2eNormWm[iDtype][iptb+3]->Scale(Xsectrion4e); | |
613 | } | |
614 | } | |
615 | ||
616 | // sum up to 12 GeV/c for D0 and D+ to compare 2 different methods | |
617 | TH1D *hD2eNormSumM1[2]; | |
618 | TH1D *hD2eNormSumM2[2]; | |
619 | for(int iDtype=0; iDtype<2; iDtype++){ | |
620 | hname="hD2eNormSumM1" + iDtype; | |
621 | hD2eNormSumM1[iDtype]=(TH1D*)hD2eNorm[iDtype][3]->Clone(hname); | |
622 | hname="hD2eNormSumM2" + iDtype; | |
623 | hD2eNormSumM2[iDtype]=(TH1D*)hD2eNormWm[iDtype][3]->Clone(hname); | |
624 | for(int iptb=1; iptb<6; iptb++){ | |
625 | hD2eNormSumM1[iDtype]->Add(hD2eNorm[iDtype][iptb+3]); | |
626 | hD2eNormSumM2[iDtype]->Add(hD2eNormWm[iDtype][iptb+3]); | |
627 | } | |
628 | } | |
629 | // drawing ------------------------------------------------------- | |
630 | c2Methods = new TCanvas("c2Methods","e spectra from 2 different methods",0,0,1000,700); | |
631 | c2Methods->Divide(3,1); | |
632 | TH1D* hRatio2mtd[2]; | |
633 | TLegend *legM[3]; | |
634 | for(int i=0; i<2; i++){ | |
635 | c2Methods->cd(i+1); | |
636 | legM[i]= new TLegend(0.25,0.70,0.75,0.85); | |
637 | hD2eNormSumM1[i]->SetMarkerStyle(20); | |
638 | hD2eNormSumM1[i]->Draw(); | |
639 | hD2eNormSumM1[i]->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
640 | hD2eNormSumM2[i]->SetMarkerColor(2); | |
641 | hD2eNormSumM2[i]->SetMarkerStyle(20); | |
642 | hD2eNormSumM2[i]->Draw("same"); | |
643 | hname="hRatio2mtdD" + i; | |
644 | hRatio2mtd[i]=(TH1D*)hD2eNormSumM2[i]->Clone(hname); | |
645 | hRatio2mtd[i]->Divide(hD2eNormSumM1[i]); | |
646 | gPad->SetLogy(); | |
647 | } | |
648 | legM[2]= new TLegend(0.25,0.70,0.75,0.85); | |
649 | legM[0]->AddEntry(hD2eNormSumM1[0],"2.0<mother p_{t}<12.0 GeV/c",""); | |
650 | legM[0]->AddEntry(hD2eNormSumM1[0],"D^{0}#rightarrow e(use of Xsection)","p"); | |
651 | legM[0]->AddEntry(hD2eNormSumM2[0],"D^{0}#rightarrow e(use of weighting)","p"); | |
652 | legM[1]->AddEntry(hD2eNormSumM1[1],"2.0<mother p_{t}<12.0 GeV/c",""); | |
653 | legM[1]->AddEntry(hD2eNormSumM1[1],"D^{+}#rightarrow e(use of Xsection)","p"); | |
654 | legM[1]->AddEntry(hD2eNormSumM2[1],"D^{+}#rightarrow e(use of weighting)","p"); | |
655 | ||
656 | legM[2]->AddEntry(hRatio2mtd[0],"2.0<mother p_{t}<12.0 GeV/c",""); | |
657 | legM[2]->AddEntry(hRatio2mtd[0],"D^{0}#rightarrow e: (use of weighting)/(use of Xsection)","p"); | |
658 | legM[2]->AddEntry(hRatio2mtd[1],"D^{+}#rightarrow e: (use of weighting)/(use of Xsection)","p"); | |
659 | c2Methods->cd(1); | |
660 | legM[0]->Draw("same"); | |
661 | c2Methods->cd(2); | |
662 | legM[1]->Draw("same"); | |
663 | ||
664 | setLegendStyle(*legM[0],0); | |
665 | setLegendStyle(*legM[1],0); | |
666 | setLegendStyle(*legM[2],0); | |
667 | ||
668 | c2Methods->cd(3); | |
669 | hRatio2mtd[0]->Draw(); | |
670 | hRatio2mtd[1]->SetLineColor(1); | |
671 | hRatio2mtd[1]->SetMarkerColor(1); | |
672 | hRatio2mtd[1]->Draw("same"); | |
673 | legM[2]->Draw("same"); | |
674 | //---------------------------------------------------------------- | |
675 | ||
676 | ||
677 | // produce output file | |
678 | TFile *f = new TFile("Testoutput.root","recreate"); | |
679 | f->cd(); | |
680 | hD2eNormTotal->Write(); | |
681 | gD2eNormTotal->Write(); | |
682 | gD2eNormTotalwMCstat->Write(); | |
683 | f->Close(); | |
684 | ||
685 | ||
686 | // compare with FONLL ---------------------------------------------- | |
687 | h4FONLL= new TH1F("h4FONLL","h4FONLL",60,0.25,30.25); | |
688 | ifstream stream1("./e_from_D-FONLL-abseta.lt.0.9.txt"); | |
689 | ||
690 | double b; | |
691 | if(!stream1) | |
692 | cout << "While opening a file an error is encountered" << endl; | |
693 | ||
694 | int i=0; | |
695 | int ipt2=0; | |
696 | int k=0; | |
697 | Double_t val[8][60]; | |
698 | Double_t exl[60]; | |
699 | Double_t exh[60]; | |
700 | Double_t eyl[60]; | |
701 | Double_t eyh[60]; | |
702 | while(!stream1.eof()) | |
703 | { | |
704 | stream1 >> b; | |
705 | ||
706 | i=k%8; | |
707 | ipt2=int(double(k)/8.); | |
708 | val[i][ipt2] = b; | |
709 | k++; | |
710 | if(ipt2==59 && i==7) break; | |
711 | } | |
712 | ||
713 | for(int i=0; i<60; i++){ | |
714 | h4FONLL->Fill(val[0][i],val[1][i]); | |
715 | } | |
716 | ||
717 | for(int m=0; m<60; m++){ | |
718 | exl[m]=0.0; | |
719 | exh[m]=0.0; | |
720 | eyl[m]=val[1][m]-val[2][m]; | |
721 | eyh[m]=val[3][m]-val[1][m]; | |
722 | } | |
723 | ||
724 | TGraphAsymmErrors* gr=new TGraphAsymmErrors(60,val[0],val[1],exl,exh,eyl,eyh); | |
725 | gr->SetMarkerStyle(20); | |
726 | gr->SetMarkerColor(1); | |
727 | gr->SetMarkerSize(0.5); | |
728 | gr->SetLineColor(3); | |
729 | gr->SetFillColor(5); | |
730 | gr->GetXaxis()->SetTitle("p_{t} [GeV/c]"); | |
731 | gr->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9 [pb/GeV/c]"); | |
732 | ||
733 | //gr->Draw("ZACPe3"); | |
734 | //gr->Draw("ZCe3"); | |
735 | //gr->Draw("Psame"); | |
736 | cD2eFONLL = new TCanvas("cD2eFONLL","pT of e from D and FONLL",0,0,1000,600); | |
737 | cD2eFONLL->cd(); | |
738 | gr->Draw("ACP"); | |
739 | gD2eNormTotalwMCstat->Draw("psame"); | |
740 | TLegend *legwFON = new TLegend(0.50,0.75,0.75,0.85); | |
741 | legwFON->AddEntry(gr,"D#rightarrow e FONLL","p"); | |
742 | legwFON->AddEntry(gD2eNormTotalwMCstat,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p"); | |
743 | setLegendStyle(*legwFON,0); | |
744 | legwFON->Draw("same"); | |
745 | gPad->SetLogy(); | |
746 | ||
747 | } | |
748 | ||
749 | //-------------------------- | |
750 | void setDataStyle(TH1D &h, Int_t lc, Int_t lw, Int_t ls){ | |
751 | ||
752 | h.SetLineColor(lc); | |
753 | h.SetLineWidth(lw); | |
754 | h.SetLineStyle(ls); | |
755 | h.SetMarkerColor(lc); | |
756 | ||
757 | } | |
758 | ||
759 | //-------------------------- | |
760 | void setLegendStyle(TLegend &legend, Int_t bs){ | |
761 | ||
762 | legend.SetBorderSize(bs); | |
763 | legend.SetFillColor(0); | |
764 | legend.SetTextFont(132); | |
765 | legend.SetTextSize(0.04); | |
766 | legend.SetMargin(0.2); | |
767 | ||
768 | } | |
769 | ||
770 | //-------------------------- | |
771 | void setPadStyle(Int_t lvl, TPad *pad){ | |
772 | ||
773 | pad->SetLogy(); | |
774 | if(lvl>0) gPad->SetGridy(); | |
775 | if(lvl>1) gPad->SetGridx(); | |
776 | ||
777 | } | |
778 | ||
779 | //-------------------------- | |
780 | void setGeneralStyle(){ | |
781 | ||
782 | gStyle->SetPalette(1); | |
783 | gStyle->SetFrameBorderMode(0); | |
784 | gStyle->SetFrameFillColor(0); | |
785 | gStyle->SetPadBorderSize(0); | |
786 | gStyle->SetPadBorderMode(0); | |
787 | gStyle->SetCanvasColor(0); | |
788 | gStyle->SetCanvasBorderSize(10); | |
789 | gStyle->SetOptStat(0); | |
790 | gStyle->SetOptTitle(0); | |
791 | gStyle->SetTitleFillColor(10); | |
792 | gStyle->SetTitleBorderSize(0); | |
793 | gStyle->SetTitleSize(0.04,"X"); | |
794 | gStyle->SetTitleSize(0.04,"Y"); | |
795 | gStyle->SetTitleFont(132,"X"); | |
796 | gStyle->SetTitleFont(132,"Y"); | |
797 | gStyle->SetTitleXOffset(0.9); | |
798 | gStyle->SetTitleYOffset(0.9); | |
799 | gStyle->SetLabelFont(132,"X"); | |
800 | gStyle->SetLabelFont(132,"Y"); | |
801 | gStyle->SetLabelSize(0.04,"X"); | |
802 | gStyle->SetLabelSize(0.04,"Y"); | |
803 | gStyle->SetTitleSize(0.05,"X"); | |
804 | gStyle->SetTitleSize(0.05,"Y"); | |
805 | gStyle->SetLineWidth(2); | |
806 | ||
807 | } | |
808 | ||
809 | void CorrectFromTheWidth(TH1D *h1) const { | |
810 | // | |
811 | // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1} | |
812 | // | |
813 | ||
814 | TAxis *axis = h1->GetXaxis(); | |
815 | Int_t nbinX = h1->GetNbinsX(); | |
816 | ||
817 | for(Int_t i = 1; i <= nbinX; i++) { | |
818 | ||
819 | Double_t width = axis->GetBinWidth(i); | |
820 | Double_t content = h1->GetBinContent(i); | |
821 | Double_t error = h1->GetBinError(i); | |
822 | h1->SetBinContent(i,content/width); | |
823 | h1->SetBinError(i,error/width); | |
824 | Double_t pt = h1->GetBinCenter(i); | |
825 | } | |
826 | ||
827 | } |