]>
Commit | Line | Data |
---|---|---|
014044ba | 1 | #include "TH1.h" |
2 | #include "TH3.h" | |
3 | #include "TF1.h" | |
4 | #include "TF3.h" | |
5 | #include "TLegend.h" | |
6 | #include "TCanvas.h" | |
7 | #include "TTree.h" | |
8 | #include "TFile.h" | |
9 | #include "TLine.h" | |
10 | #include "TPaveText.h" | |
11 | #include "TStyle.h" | |
12 | #include "TROOT.h" | |
13 | #include "TSystem.h" | |
14 | #include "TNtuple.h" | |
15 | #include "TArrow.h" | |
16 | #include "TGraphErrors.h" | |
17 | #include "TVirtualFitter.h" | |
18 | #include "TMath.h" | |
19 | #include "TMinuit.h" | |
20 | #include "TSpectrum.h" | |
21 | ||
22 | // RooFit | |
23 | #include "RooGenericPdf.h" | |
24 | #include "RooFitResult.h" | |
25 | #include "RooRealVar.h" | |
26 | ||
27 | ||
28 | /////////////////////////// | |
29 | // Default Options (Global) | |
30 | /////////////////////////// | |
31 | ||
32 | //________________ | |
33 | // container check | |
34 | ||
35 | // display container contents | |
36 | Bool_t IsShowContOn = kFALSE; | |
37 | // check bin error | |
38 | Bool_t IsCheckBinError = kFALSE; | |
39 | ||
40 | //_______________________ | |
41 | // User setting variables | |
42 | ||
43 | // enumerate variables | |
44 | enum { NEVENT, Y, PT, IMASS, TRIG, PTMU, PMU, TRIGSIDE, RABSMU, CHARGE, ETAMU, THABSMU, VZMU, DCAMU }; | |
45 | ||
46 | // set efficiency matrix bin | |
47 | const Char_t *var0name="y"; | |
48 | const Char_t *var1name="pt"; | |
49 | const Int_t var0bin = 5; | |
50 | const Int_t var1bin = 7; | |
51 | Double_t var0lim[var0bin+1] = {-4.0,-3.7,-3.4,-3.1,-2.8,-2.5}; | |
52 | Double_t var1lim[var1bin+1] = {0,2,4,6,8,10,15,20}; | |
53 | const Double_t binnorm = (Double_t)5/2; // for the last two pt bin | |
54 | ||
55 | //______________ | |
56 | // Analysis cuts | |
57 | ||
58 | // Use pt trigger cut | |
59 | Bool_t IsHpt = kTRUE; | |
60 | ||
61 | //_______________ | |
62 | // fit parameters | |
63 | Bool_t UseRooFit = kFALSE; | |
64 | if(UseRooFit) { using namespace RooFit; } | |
65 | ||
66 | // rebinning | |
67 | Bool_t rebinned = kTRUE; | |
68 | // 3 gauss for UPSILON family | |
69 | Char_t *resname = "Upsilon"; | |
70 | Double_t fitrange[2] = {7, 12}; // fit range min, max | |
71 | Double_t usrrange[2] = {7, 12}; // set range user | |
72 | // setting fit paramter : single exponential(par[2]) + 3 gaussian(par[9]) | |
73 | // param[11] = { A, B, N_Y(1S), mean_Y(1S), sigma_Y(1S), N_Y(2S), mean_Y(2S), sigma_Y(2S), N_Y(3S), mean_Y(3S), sigma_Y(3S) | |
74 | Double_t param[11] = {3000,-0.3,6.6,9.46,0.230,1.8,10.023,0.230,1,10.355,0.230}; | |
75 | // fix parameter switch | |
76 | Bool_t IsFixSigma = kTRUE; | |
77 | ||
78 | //////////////////////////// | |
79 | // UPSILON ANALYSIS MACRO // | |
80 | //////////////////////////// | |
81 | ||
82 | //______________________________________________________ | |
83 | ||
84 | Bool_t upsilonCORRFW(const char *dataname, // input file name | |
85 | const Bool_t IsCreateEffMat=kFALSE, // switch for efficiency matrix creation | |
86 | const Bool_t IsExtSig=kFALSE, // switch for fitting | |
87 | const Bool_t IsDoCorrection=kFALSE // switch for correction | |
88 | ) | |
89 | { | |
90 | // print out startup messages | |
91 | printf("******************************************************\n"); | |
92 | printf(" Starting Correction Framework for Upsilon...\n"); | |
93 | printf("******************************************************\n"); | |
94 | printf("\nChecking options...\n"); | |
95 | printf("Create new efficiency matrix...%s\n",IsCreateEffMat ? "yes":"no"); | |
96 | printf("Extract Signal...%s\n",IsExtSig ? "yes":"no"); | |
97 | printf("Do correction...%s\n",IsDoCorrection ? "yes":"no"); | |
98 | printf("Checking done...\n"); | |
99 | ||
100 | // output files | |
101 | char *file_effmat = ""; // efficiency matrix | |
102 | char *file_datmat = ""; // data matrix | |
103 | char *file_result = ""; // correction result | |
104 | ||
105 | // call a fuction to construct efficiency matrix | |
106 | if(IsCreateEffMat) file_effmat = createEffMat(dataname); | |
107 | else { // use default efficiency matrix | |
108 | //file_effmat = "effCont_mat_residual_aod_50k_pdccut.root"; | |
109 | file_effmat = "efficiency.container.CFMuonResUpsilon.local.AOD.MC.root"; | |
110 | printf("Default efficiency matrix: %s\n",file_effmat); | |
111 | } | |
112 | ||
113 | // call a function to perform fit to extract #nignal | |
114 | if(IsExtSig) file_datmat = extSignal(dataname); | |
115 | else { | |
116 | file_datmat = "data.container.CFMuonResUpsilon.PDC09.local.AOD.MC.root"; | |
117 | printf("Default data container: %s\n",file_datmat); | |
118 | } | |
119 | ||
120 | // call a function to perform a correction on data | |
121 | if(IsDoCorrection) file_result = doCorrection(dataname, file_effmat, file_datmat); | |
122 | ||
123 | // print out closing messages | |
124 | printf("******************************************************\n"); | |
125 | printf(" Finishing Correction Framework for Upsilon...\n"); | |
126 | printf("******************************************************\n"); | |
127 | ||
128 | // display results | |
129 | printf("\n here are results!!\n"); | |
130 | printf("efficiency container is %s\n",file_effmat); | |
131 | printf("data container is %s\n",file_datmat); | |
132 | printf("correction result file is %s\n",file_result); | |
133 | ||
134 | return kTRUE; | |
135 | } | |
136 | ||
137 | ////////////////////////////// | |
138 | // Create Efficiency matrix // | |
139 | ////////////////////////////// | |
140 | ||
141 | char* createEffMat(const char *dataname) | |
142 | { | |
143 | ||
144 | // print out starting messages | |
145 | printf("\nNow starting efficiency matrix creation...\n"); | |
146 | ||
147 | LoadLib(); | |
148 | ||
149 | // open input container | |
150 | printf("Opening container from %s...",dataname); | |
151 | TFile *file = new TFile(dataname); | |
152 | AliCFContainer *cont = (AliCFContainer*) (file->Get("container")); | |
153 | printf("done\n"); | |
154 | ||
155 | // print out container step and variables | |
156 | //cont->Print(); | |
157 | // get N step | |
158 | Int_t nstep = cont->GetNStep(); | |
159 | // get N variables | |
160 | Int_t nvar = cont->GetNVar(); | |
161 | printf("Nstep: %d\n", nstep); | |
162 | printf("Nvar: %d\n", nvar); | |
163 | // casting constant for use of array | |
164 | const Int_t stepnum = nstep; | |
165 | const Int_t varnum = nvar; | |
166 | // set variable: bin, min and max | |
167 | Double_t varbin[varnum]={ 5, 15, 100, 300, 40, 100, 100, 4, 300, 6, 100, 200, 100, 100 }; | |
168 | Double_t varmin[varnum]={ 0, -4, 0, 0, 0, 0, 0, 0, 15, 0, -4.5, 170, -50, 0 }; | |
169 | Double_t varmax[varnum]={ 5,-2.5, 100, 300, 4, 100, 100, 4, 90, 6, -2.0, 180, 50, 100 }; | |
170 | // set variable epsilon | |
171 | Double_t varEpsilon[varnum]; | |
172 | for (Int_t ie=0; ie<varnum; ie++) varEpsilon[ie] = (varmax[ie]-varmin[ie])/(100*varbin[ie]); | |
173 | // Display container contents | |
174 | if(IsShowContOn) { | |
175 | TCanvas *csc[stepnum]; | |
176 | ||
177 | for(Int_t i=0; i<2; i++) { | |
178 | csc[i] = new TCanvas(Form("csc%d",i),Form("container(%d)",i),0,0,1000,600); | |
179 | csc[i]->Divide(nvar/2,2); | |
180 | for(Int_t j=0; j<nvar; j++) { | |
181 | csc[i]->cd(j+1); cont->ShowProjection(j,i)->Draw(); | |
182 | } | |
183 | } | |
184 | } | |
185 | ||
186 | //___________ | |
187 | // GridSparse | |
188 | ||
189 | // GridSparse from container | |
190 | AliCFGridSparse *genSpar = (AliCFGridSparse*)cont->GetGrid(0); // GEN | |
191 | AliCFGridSparse *recSpar = (AliCFGridSparse*)cont->GetGrid(1); // REC | |
192 | ||
193 | // variables for efficiency matrix | |
194 | const Int_t nStep=2; // number of steps: MC and REC | |
195 | const Int_t nVar=2; // number of variables: var0 and var1 | |
196 | const Int_t var0dim=var0bin; // number of bin in var0 dimension of efficiency matrix | |
197 | const Int_t var1dim=var1bin; // number of bin in var1 dimension of efficiency matrix | |
198 | const Int_t binMat[nVar]={var0dim,var1dim}; // summary array | |
199 | ||
200 | // create new container for efficiency | |
201 | AliCFContainer *effCont = new AliCFContainer("effCont","Upsilon efficiency matrix",nStep,nVar,binMat); | |
202 | effCont->SetBinLimits(0,var0lim); | |
203 | effCont->SetBinLimits(1,var1lim); | |
204 | // create new gridsparse for efficiency | |
205 | // gridsparse of generated events | |
206 | AliCFGridSparse *genGrid = new AliCFGridSparse("genGrid","Generated Upsilon matrix",nVar,binMat); | |
207 | genGrid->SetBinLimits(0,var0lim); | |
208 | genGrid->SetBinLimits(1,var1lim); | |
209 | // gridsparse of reconstructed events | |
210 | AliCFGridSparse *recGrid = new AliCFGridSparse("recGrid","Reconstructed Upsilon matrix",nVar,binMat); | |
211 | recGrid->SetBinLimits(0,var0lim); | |
212 | recGrid->SetBinLimits(1,var1lim); | |
213 | ||
214 | //____________________ | |
215 | // Loop on matrix bins | |
216 | ||
217 | // cut on single muon | |
218 | // theta cut | |
219 | varmin[THABSMU] = 171+varEpsilon[THABSMU]; | |
220 | varmax[THABSMU] = 178-varEpsilon[THABSMU]; | |
221 | // eta cut | |
222 | varmin[ETAMU] = -4.0+varEpsilon[ETAMU]; | |
223 | varmax[ETAMU] = -2.5-varEpsilon[ETAMU]; | |
224 | // rabs cut | |
225 | //varmin[RABSMU] = 17.6+varEpsilon[RABSMU]; | |
226 | //varmin[RABSMU] = 80.0-varEpsilon[RABSMU]; | |
227 | ||
228 | // canvas for fitting | |
229 | TCanvas *cfit = new TCanvas("cfit","cfit",0,0,1600,1125); | |
230 | cfit->Divide(7,5); | |
231 | ||
232 | printf("Loop on matrix bins..."); | |
233 | for(Int_t bin0=0; bin0<var0dim; bin0++) { // loop on rapidity (var0dim) | |
234 | varmin[Y] = var0lim[bin0]+varEpsilon[Y]; | |
235 | varmax[Y] = var0lim[bin0+1]-varEpsilon[Y]; | |
236 | printf("%s range: %.2f < %s < %.2f\n",var0name,varmin[Y],var0name,varmax[Y]); | |
237 | ||
238 | for(Int_t bin1=0; bin1<var1dim; bin1++) { // loop on pt (var1dim) | |
239 | varmin[PT] = var1lim[bin1]+varEpsilon[PT]; | |
240 | varmax[PT] = var1lim[bin1+1]-varEpsilon[PT]; | |
241 | printf("%s range: %.2f < %s < %.2f\n",var1name,varmin[PT],var1name,varmax[PT]); | |
242 | ||
243 | TH1D* gmass = (TH1D*)genSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
244 | TH1D* rmass = (TH1D*)recSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
245 | ||
246 | cfit->cd((bin0*7)+bin1+1); | |
247 | if(rebinned) rmass->Rebin(2); | |
248 | // fit reconstructed resonances to get the #_signal | |
249 | // declare array to retreive fit results | |
250 | Double_t par[4]; | |
251 | Double_t *parErr; | |
252 | // landau (x) gauss : #_signal,mean,sigma,width | |
253 | TF1 *func = new TF1("func",gaussXlandau,8,10,4); | |
254 | func->SetParameters(30,9.46,0.09,0.004); | |
255 | // simple gauss | |
256 | //TF1 *func = new TF1("func",normGauss,8,10,3); | |
257 | //func->SetParameters(30,9.46,0.09); | |
258 | func->SetLineColor(kBlue); | |
259 | // fit with likelihood method | |
260 | rmass->Fit(func,"0L"); | |
261 | // retreive fit parameters | |
262 | func->GetParameters(&par[0]); | |
263 | // fit second times | |
264 | func->SetParameters(par); | |
265 | rmass->Fit(func,"RL+"); | |
266 | // retreive fit parameters | |
267 | func->GetParameters(&par[0]); | |
268 | parErr = func->GetParErrors(); | |
269 | rmass->Draw("E"); | |
270 | cfit->Update(); | |
271 | ||
272 | Double_t genValue = bin1<5 ? gmass->Integral():gmass->Integral()/binnorm; | |
273 | Double_t genValueErr = TMath::Sqrt(genValue); | |
274 | Double_t recValue = bin1<5 ? par[0]/0.1:(par[0]/0.1)/binnorm; | |
275 | Double_t recValueErr = bin1<5 ? parErr[0]/0.1:(parErr[0]/0.1)/binnorm; | |
276 | ||
277 | printf("genValue = %.0f +/- %.0f\n",genValue,genValueErr); | |
278 | printf("recValue = %.0f +/- %.0f\n",recValue,recValueErr); | |
279 | ||
280 | Int_t binCell[nVar]={bin0+1,bin1+1}; | |
281 | genGrid->SetElement(binCell,genValue); | |
282 | genGrid->SetElementError(binCell,genValueErr); | |
283 | recGrid->SetElement(binCell,recValue); | |
284 | recGrid->SetElementError(binCell,recValueErr); | |
285 | } | |
286 | } | |
287 | printf("done\n"); | |
288 | printf("saving fit plots as %s...",Form("recfitplots.%s",dataname)); | |
289 | cfit->SaveAs(Form("recfitplots.%s",dataname)); | |
290 | printf("done\n"); | |
291 | ||
292 | ||
293 | printf("fill efficiency container..."); | |
294 | //__________________________ | |
295 | // fill efficiency container | |
296 | effCont->SetGrid(0,genGrid); | |
297 | effCont->SetGrid(1,recGrid); | |
298 | TH1D *gvar0 = effCont->ShowProjection(0,0); | |
299 | TH1D *gvar1 = effCont->ShowProjection(1,0); | |
300 | TH1D *rvar0 = effCont->ShowProjection(0,1); | |
301 | TH1D *rvar1 = effCont->ShowProjection(1,1); | |
302 | TCanvas *cvars = new TCanvas("cvars","variables",0,0,800,800); | |
303 | cvars->Divide(2,2); | |
304 | cvars->cd(1); gvar0->Draw(); | |
305 | cvars->cd(2); gvar1->Draw(); | |
306 | cvars->cd(3); rvar0->Draw(); | |
307 | cvars->cd(4); rvar1->Draw(); | |
308 | printf("done\n"); | |
309 | // save efficiency container | |
310 | char *outfile = Form("efficiency.%s",dataname); | |
311 | printf("saving efficiency container as %s...",outfile); | |
312 | effCont->Save(outfile); | |
313 | printf("done\n"); | |
314 | ||
315 | //_________________________ | |
316 | // create efficiency matrix | |
317 | printf("creating efficiency matrix..."); | |
318 | AliCFEffGrid *matEff = new AliCFEffGrid("matEff",Form("%s efficiency matrix",resname),*effCont); | |
319 | printf("done\n"); | |
320 | // calcualte efficiency | |
321 | printf("calculating efficiency..."); | |
322 | matEff->CalculateEfficiency(1,0); // REC over MC | |
323 | printf("done\n"); | |
324 | ||
325 | //________________________ | |
326 | // display efficiency plot | |
327 | TCanvas *ceff = new TCanvas("ceff",Form("%s efficiency",resname),0,0,1200,400); | |
328 | ceff->Divide(3,1); | |
329 | // efficiency over var0 | |
330 | ceff->cd(1); | |
331 | TH1D *effvar0 = matEff->Project(0); | |
332 | effvar0->SetName(Form("eff_%s",var0name)); | |
333 | effvar0->SetMinimum(0); effvar0->SetMaximum(1); | |
334 | effvar0->SetTitle(Form("%s efficiency(%s)",resname,var0name)); | |
335 | effvar0->GetXaxis()->SetTitle(var0name); | |
336 | effvar0->GetYaxis()->SetTitle("Efficiency"); | |
337 | effvar0->Draw("error"); | |
338 | // efficiency over var1 | |
339 | ceff->cd(2); | |
340 | TH1D *effvar1 = matEff->Project(1); | |
341 | effvar1->SetName(Form("eff_%s",var1name)); | |
342 | effvar1->SetMinimum(0); effvar1->SetMaximum(1); | |
343 | effvar1->SetTitle(Form("%s efficiency(%s)",resname,var1name)); | |
344 | effvar1->GetXaxis()->SetTitle(var1name); | |
345 | effvar1->GetYaxis()->SetTitle("Efficiency"); | |
346 | effvar1->Draw("error"); | |
347 | // 2-dimensional efficiency plot (var0, var1) | |
348 | ceff->cd(3); | |
349 | TH1D *eff2D = matEff->Project(0,1); | |
350 | eff2D->SetName("eff_2D"); | |
351 | eff2D->SetMinimum(0); eff2D->SetMaximum(1); | |
352 | eff2D->SetTitle(Form("%s efficiency(%s,%s)",resname,var0name,var1name)); | |
353 | eff2D->GetXaxis()->SetTitle(var0name); | |
354 | eff2D->GetYaxis()->SetTitle(var1name); | |
355 | eff2D->GetZaxis()->SetTitle("Efficiency"); | |
356 | eff2D->Draw("colz"); | |
357 | ||
358 | printf("\nEfficiency matrix created...\n"); | |
359 | return outfile; | |
360 | } | |
361 | ||
362 | //////////////////// | |
363 | // Extract Signal // | |
364 | //////////////////// | |
365 | ||
366 | char* extSignal(const char *dataname) | |
367 | { | |
368 | ||
369 | // print out starting messages | |
370 | printf("\nNow extracting #signal by fitting...\n"); | |
371 | ||
372 | LoadLib(); | |
373 | ||
374 | // open input container | |
375 | printf("Opening container from %s...",dataname); | |
376 | TFile *file = new TFile(dataname); | |
377 | AliCFContainer *cont = (AliCFContainer*) (file->Get("container")); | |
378 | printf("done\n"); | |
379 | ||
380 | // check dataname : PDC, LHC or Unknown (1, 2, 0) | |
381 | Int_t dFlag = 0; | |
382 | if(strstr(dataname,"PDC09")) { printf("this is PDC09 data...\n"); dFlag = 1; } | |
383 | else if(strstr(dataname,"LHC")) { printf("this is LHC data...\n"); dFlag = 2; } | |
384 | else { printf("unknown data...\n"); dFlag = 0; } | |
385 | ||
386 | // print out container step and variables | |
387 | // get N step | |
388 | Int_t nstep = cont->GetNStep(); | |
389 | // get N variables | |
390 | Int_t nvar = cont->GetNVar(); | |
391 | printf("Nstep: %d\n", nstep); | |
392 | printf("Nvar: %d\n", nvar); | |
393 | // casting constant for use of array | |
394 | const Int_t stepnum = nstep; | |
395 | const Int_t varnum = nvar; | |
396 | // set variable: bin, min and max | |
397 | Double_t varbin[varnum]={ 5, 15, 100, 300, 40, 100, 100, 4, 300, 6, 100, 200, 100, 100 }; | |
398 | Double_t varmin[varnum]={ 0, -4, 0, 0, 0, 0, 0, 0, 15, 0, -4.5, 170, -50, 0 }; | |
399 | Double_t varmax[varnum]={ 5,-2.5, 100, 300, 4, 100, 100, 4, 90, 6, -2.0, 180, 50, 100 }; | |
400 | // set variable epsilon | |
401 | Double_t varEpsilon[varnum]; | |
402 | for (Int_t ie=0; ie<varnum; ie++) varEpsilon[ie] = (varmax[ie]-varmin[ie])/(100*varbin[ie]); | |
403 | // Display container contents | |
404 | if(IsShowContOn) { | |
405 | TCanvas *csc[stepnum]; | |
406 | ||
407 | for(Int_t i=0; i<nstep; i++) { | |
408 | csc[i] = new TCanvas(Form("csc%d",i),Form("container(%d)",i),0,0,1000,600); | |
409 | csc[i]->Divide(nvar/2,2); | |
410 | for(Int_t j=0; j<nvar; j++) { | |
411 | csc[i]->cd(j+1); cont->ShowProjection(j,i)->Draw(); | |
412 | } | |
413 | } | |
414 | ||
415 | /* | |
416 | TH1D *hCheckMass = cont->ShowProjection(IMASS,1); | |
417 | Int_t nbins = hCheckMass->GetXaxis()->GetNbins(); | |
418 | printf("nbins = %d\n", nbins); | |
419 | for(i=1;i<=nbins;i++) printf("bin %d: %f - %f\n",i, hCheckMass->GetBinError(i), TMath::Sqrt(hCheckMass->GetBinContent(i))); | |
420 | TCanvas *ctemp = new TCanvas("ctemp","ctemp",0,0,500,500); | |
421 | hCheckMass->Draw(); | |
422 | */ | |
423 | } | |
424 | ||
425 | //___________ | |
426 | // GridSparse | |
427 | ||
428 | // GridSparse from container | |
429 | AliCFGridSparse *genSpar = (AliCFGridSparse*)cont->GetGrid(0); // GEN for PDC09 | |
430 | AliCFGridSparse *cintSpar = (AliCFGridSparse*)cont->GetGrid(1); // REC for PDC09 || CINT && !CMUS | |
431 | AliCFGridSparse *cmusSpar = (AliCFGridSparse*)cont->GetGrid(4); // CMUS only | |
432 | ||
433 | // variables for data matrix | |
434 | const Int_t nStep=2; // number of steps: REC(CINT && !CMUS) and CMUS only | |
435 | const Int_t nVar=2; // number of variables: var0 and var1 | |
436 | const Int_t var0dim=var0bin; // number of bin in var0 dimension of data matrix | |
437 | const Int_t var1dim=var1bin; // number of bin in var1 dimension of data matrix | |
438 | const Int_t binMat[nVar]={var0dim,var1dim}; // summary array | |
439 | ||
440 | // create new container for data | |
441 | AliCFContainer *dataCont = new AliCFContainer("dataCont","Upsilon data matrix",nStep,nVar,binMat); | |
442 | dataCont->SetBinLimits(0,var0lim); | |
443 | dataCont->SetBinLimits(1,var1lim); | |
444 | // create new gridsparse for data | |
445 | // gridsparse of REC (or CINT && !CMUS) events | |
446 | AliCFGridSparse *cintGrid = new AliCFGridSparse("cintGrid","Upsilon data matrix(CINT)",nVar,binMat); | |
447 | cintGrid->SetBinLimits(0,var0lim); | |
448 | cintGrid->SetBinLimits(1,var1lim); | |
449 | // gridsparse of CMUS only events | |
450 | AliCFGridSparse *cmusGrid = new AliCFGridSparse("cmusGrid","Upsilon matrix(CMUS)",nVar,binMat); | |
451 | cmusGrid->SetBinLimits(0,var0lim); | |
452 | cmusGrid->SetBinLimits(1,var1lim); | |
453 | ||
454 | //_____________________ | |
455 | // Loops on matrix bins | |
456 | ||
457 | if(dFlag == 2) { // LHC data | |
458 | // event statistics | |
459 | for(Int_t i=0; i<5; i++) { | |
460 | varmin[NEVENT] = i+varEpsilon[NEVENT]; | |
461 | varmax[NEVENT] = i+1-varEpsilon[NEVENT]; | |
462 | Int_t nEvtInt = ((TH1D*)cintSpar->Slice(IMASS,-1,-1,varmin,varmax))->Integral(); | |
463 | Int_t nEvtMus = ((TH1D*)cmusSpar->Slice(IMASS,-1,-1,varmin,varmax))->Integral(); | |
464 | printf("#Event (%d) = %d(CINT1 = %d, CMUS1 = %d)\n",i,nEvtInt+nEvtMus,nEvtInt,nEvtMus); | |
465 | } | |
466 | ||
467 | // common cut on events | |
468 | // Rabs cut | |
469 | varmin[RABSMU] = 17.6+varEpsilon[RABSMU]; varmax[RABSMU] = 80.0-varEpsilon[RABSMU]; | |
470 | // acceptance cut on dimuon | |
471 | varmin[Y] = -4+varEpsilon[Y]; varmax[Y] = -2.5-varEpsilon[Y]; | |
472 | ||
473 | // canvas for display | |
474 | //TCanvas *cmassi = new TCanvas("cmassi","mass distribution(CINT)",0,0,900,600); | |
475 | //cmassi->Divide(3,2); | |
476 | //TCanvas *cmassm = new TCanvas("cmassm","mass distribution(CMUS)",0,0,900,600); | |
477 | //cmassm->Divide(3,2); | |
478 | ||
479 | TH1D *hmassi[3], *hmassm[3]; | |
480 | ||
481 | // physics selection (0: off, 1: on) | |
482 | for(Int_t i=0;i<2;i++) { | |
483 | // track-trigger matching | |
484 | for(Int_t j=0;j<3;j++) { | |
485 | if(IsHpt) { | |
486 | switch (j) { | |
487 | case 0: | |
488 | varmin[TRIG] = 0+varEpsilon[TRIG]; | |
489 | break; | |
490 | case 1: | |
491 | varmin[TRIG] = 3+varEpsilon[TRIG]; | |
492 | break; | |
493 | case 2: | |
494 | varmin[TRIG] = 3.3+varEpsilon[TRIG]; | |
495 | break; | |
496 | } | |
497 | } | |
498 | else { | |
499 | switch (j) { | |
500 | case 0: | |
501 | varmin[TRIG] = 0+varEpsilon[TRIG]; | |
502 | break; | |
503 | case 1: | |
504 | varmin[TRIG] = 2+varEpsilon[TRIG]; | |
505 | break; | |
506 | case 2: | |
507 | varmin[TRIG] = 2.2+varEpsilon[TRIG]; | |
508 | break; | |
509 | } | |
510 | } | |
511 | // (CINT && !CMUS) | |
512 | varmin[NEVENT] = i+1+varEpsilon[NEVENT]; varmax[NEVENT] = i+2-varEpsilon[NEVENT]; | |
513 | hmassi[j] = (TH1D*)cintSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
514 | hmassi[j]->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]); | |
515 | //printf("#Event (CINT1B-%dmatch,hpt) = %d\n",j,hmassi[j]->Integral()); | |
516 | //cmassi->cd(i*3+j+1); | |
517 | //hmassi[j]->Draw(); | |
518 | ||
519 | // (CMUS only) | |
520 | varmin[NEVENT] = i+3+varEpsilon[NEVENT]; varmax[NEVENT] = i+4-varEpsilon[NEVENT]; | |
521 | hmassm[j] = (TH1D*)cmusSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
522 | if(rebinned) hmassm[j]->Rebin(4); | |
523 | hmassm[j]->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]); | |
524 | //printf("#Event (CMUS1B-%dmatch,hpt) = %d\n",j,hmassm[j]->Integral()); | |
525 | //cmassm->cd(i*3+j+1); | |
526 | //hmassm[j]->Draw(); | |
527 | } | |
528 | } | |
529 | ||
530 | // draw N matching mass plot and fit | |
531 | TCanvas *c1[3]; | |
532 | TPaveText *info; | |
533 | ||
534 | Char_t *ptext[6]; | |
535 | ||
536 | for(Int_t i=0; i<3; i++) { | |
537 | ||
538 | Double_t par[15] = {0, }; | |
539 | Double_t *parErr = 0;; | |
540 | Double_t chi2=0, ndf=0, mean=0, mean_err=0, sigma=0, sigma_err=0, imin=0, imax=0, norm=0, norm_err=0, nsig=0, nsig_err=0,nbkg=0; | |
541 | ||
542 | //cout << "Bin cont. = " << hmassm[i]->GetBinContent(20) << " +- " << hmassm[i]->GetBinError(20) << endl; | |
543 | //cout << "sqrt(Bin cont.) = " << TMath::Sqrt(hmassm[i]->GetBinContent(20)) << endl; | |
544 | ||
545 | TF1 *bkg = new TF1("bkg",expo,fitrange[0],fitrange[1],2); | |
546 | bkg->SetParameters(param[0],param[1]); | |
547 | hmassm[i]->Fit(bkg,"NR"); | |
548 | bkg->GetParameters(¶m[0]); | |
549 | ||
550 | TF1 *gl = new TF1("gl",global,fitrange[0],fitrange[1],11); | |
551 | gl->SetParameters(¶m[0]); | |
552 | ||
553 | // fix yield 2S and 3S normalized from CMS result | |
554 | gl->FixParameter(5,1.8); gl->FixParameter(8,1); | |
555 | // fix mean from PDG | |
556 | gl->FixParameter(6,param[6]); gl->FixParameter(9,param[9]); | |
557 | // fix sigma from simulation with current alignment(v2) | |
558 | if(IsFixSigma) {gl->FixParameter(4,param[4]); gl->FixParameter(7,param[7]); gl->FixParameter(10,param[10]);} | |
559 | else { | |
560 | if(i!=2) {gl->SetParLimits(4,param[4]-0.2,param[4]+0.2);gl->FixParameter(7,param[7]); gl->FixParameter(10,param[10]);} | |
561 | else {gl->SetParLimits(4,param[4]-0.2,param[4]+0.2);gl->FixParameter(7,param[7]); gl->FixParameter(10,param[10]);} | |
562 | } | |
563 | ||
564 | gl->SetLineColor(kBlue); | |
565 | gl->SetLineWidth(2); | |
566 | hmassm[i]->Fit(gl,"NR"); | |
567 | gl->GetParameters(¶m[0]); | |
568 | parErr = gl->GetParErrors(); | |
569 | ||
570 | bkg->SetParameters(param[0],param[1]); | |
571 | bkg->SetParErrors(parErr); | |
572 | bkg->SetLineColor(kRed); | |
573 | bkg->SetLineStyle(2); | |
574 | bkg->SetLineWidth(2); | |
575 | ||
576 | TF1 *sig = new TF1("sig",normGauss3,fitrange[0],fitrange[1],9); | |
577 | sig->SetParameters(param[2],param[3],param[4],param[5],param[6],param[7],param[8],param[9],param[10]); | |
578 | sig->SetParErrors(&parErr[2]); | |
579 | sig->SetLineColor(kRed); | |
580 | sig->SetLineWidth(2); | |
581 | ||
582 | TF1 *sig1 = new TF1("sig1",normGauss,fitrange[0],fitrange[1],3); | |
583 | sig1->SetParameters(param[2],param[3],param[4]); | |
584 | sig1->SetParErrors(&parErr[2]); | |
585 | sig1->SetLineColor(kMagenta); | |
586 | sig1->SetLineWidth(2); | |
587 | ||
588 | TF1 *sig2 = new TF1("sig2",normGauss,fitrange[0],fitrange[1],3); | |
589 | sig2->SetParameters(param[5],param[6],param[7]); | |
590 | sig2->SetParErrors(&parErr[5]); | |
591 | sig2->SetLineColor(kViolet); | |
592 | sig2->SetLineWidth(2); | |
593 | ||
594 | TF1 *sig3 = new TF1("sig3",normGauss,fitrange[0],fitrange[1],3); | |
595 | sig3->SetParameters(param[8],param[9],param[10]); | |
596 | sig3->SetParErrors(&parErr[8]); | |
597 | sig3->SetLineColor(kPink); | |
598 | sig3->SetLineWidth(2); | |
599 | ||
600 | c1[i] = new TCanvas(Form("c1%d",i),Form("c1%d",i),0,0,800,800); | |
601 | ||
602 | c1[i]->cd(1); //c1[i]->cd(1)->SetLogy(); | |
603 | hmassm[i]->SetMinimum(0); | |
604 | hmassm[i]->DrawCopy(); | |
605 | sig->Draw("same"); | |
606 | //bkg_2->Draw("same"); | |
607 | bkg->Draw("same"); | |
608 | gl->Draw("same"); | |
609 | sig1->Draw("same"); | |
610 | sig2->Draw("same"); | |
611 | sig3->Draw("same"); | |
612 | ||
613 | chi2 = gl->GetChisquare(); | |
614 | ndf = gl->GetNDF(); | |
615 | mean = sig->GetParameter(1); | |
616 | mean_err = sig->GetParError(1); | |
617 | sigma = sig->GetParameter(2); | |
618 | sigma_err = sig->GetParError(2); | |
619 | imin = mean-2*sigma; | |
620 | imax = mean+2*sigma; | |
621 | nsig = sig->Integral(imin, imax); | |
622 | nsig_err = TMath::Sqrt(nsig); | |
623 | //nsig = sig->GetParameter(0); | |
624 | //nsig_err = sig->GetParError(0); | |
625 | if(rebinned) { | |
626 | norm = sig->GetParameter(0)/0.2; | |
627 | norm_err = sig->GetParError(0)/0.2; | |
628 | nsig = sig1->Integral(imin,imax)/0.2; | |
629 | nsig_err = TMath::Sqrt(nsig); | |
630 | nbkg = bkg->Integral(imin, imax)/0.2; | |
631 | } | |
632 | else { | |
633 | norm = sig->GetParameter(0)/0.1; | |
634 | norm_err = sig->GetParError(0)/0.1; | |
635 | nsig = sig1->Integral(imin,imax)/0.1; | |
636 | nsig_err = TMath::Sqrt(nsig); | |
637 | nbkg = bkg->Integral(imin, imax)/0.1; | |
638 | } | |
639 | ||
640 | TPaveText *info = new TPaveText(0.58,0.48,0.96,0.70,"brNDC"); | |
641 | info->InsertText(Form("N_{%s} = %.0f #pm %.0f",resname, norm, norm_err)); | |
642 | info->InsertText(Form("Mass = %.3f #pm %.3f GeV/c^{2}", mean, mean_err)); | |
643 | if(IsFixSigma) info->InsertText(Form("#sigma_{%s} = %.0f (fixed) MeV/c^{2}",resname,sigma*1000)); | |
644 | else info->InsertText(Form("#sigma_{%s} = %.0f #pm %.0f MeV/c^{2}",resname,sigma*1000, sigma_err*1000)); | |
645 | info->InsertText(Form("S/B (2#sigma)= %.3f", nsig/nbkg)); | |
646 | info->InsertText(Form("Significance (2#sigma) = %.3f" , nsig/TMath::Sqrt(nsig+nbkg))); | |
647 | ||
648 | printf("%d #geq trigger matching(%s)\n",i,IsHpt ? "Hpt" : "Lpt"); | |
649 | printf("#Chi^{2}/ndf = %.3f\n", chi2/ndf); | |
650 | printf("N_Upsilon = %.0f +- %.0f\n",norm, norm_err); | |
651 | printf("S = %f, B = %f, S/B = %f, sqrt(S+B) = %f, S/sqrt(S+B) = %f\n", nsig, nbkg, nsig/nbkg, TMath::Sqrt(nsig+nbkg), nsig/TMath::Sqrt(nsig+nbkg)); | |
652 | ||
653 | info->SetFillColor(0); | |
654 | info->Draw("same"); | |
655 | ||
656 | TLegend *leg = new TLegend(0.5928571,0.7442857,0.96,0.9542857,NULL,"brNDC"); | |
657 | leg->SetBorderSize(0); | |
658 | leg->SetTextFont(62); | |
659 | leg->SetLineColor(1); | |
660 | leg->SetLineStyle(1); | |
661 | leg->SetLineWidth(1); | |
662 | leg->SetFillColor(0); | |
663 | leg->SetFillStyle(1001); | |
664 | ||
665 | TLegendEntry *entry=leg->AddEntry(hmassm[i],"Data","lpf"); | |
666 | entry->SetFillStyle(1001); | |
667 | entry->SetLineColor(1); | |
668 | entry->SetLineStyle(1); | |
669 | entry->SetLineWidth(1); | |
670 | entry->SetMarkerColor(1); | |
671 | entry->SetMarkerStyle(8); | |
672 | entry->SetMarkerSize(1); | |
673 | ||
674 | entry=leg->AddEntry(sig,"Signal only","lpf"); | |
675 | entry->SetFillColor(19); | |
676 | entry->SetLineColor(kRed); | |
677 | entry->SetLineStyle(1); | |
678 | entry->SetLineWidth(2); | |
679 | entry->SetMarkerColor(1); | |
680 | entry->SetMarkerStyle(1); | |
681 | entry->SetMarkerSize(1); | |
682 | ||
683 | entry=leg->AddEntry(bkg,"Background only","lpf"); | |
684 | entry->SetFillColor(19); | |
685 | entry->SetLineColor(kRed); | |
686 | entry->SetLineStyle(2); | |
687 | entry->SetLineWidth(2); | |
688 | entry->SetMarkerColor(1); | |
689 | entry->SetMarkerStyle(1); | |
690 | entry->SetMarkerSize(1); | |
691 | ||
692 | entry=leg->AddEntry(gl,"Signal + Background","lpf"); | |
693 | entry->SetFillColor(19); | |
694 | entry->SetLineColor(kBlue); | |
695 | entry->SetLineStyle(1); | |
696 | entry->SetLineWidth(2); | |
697 | entry->SetMarkerColor(1); | |
698 | entry->SetMarkerStyle(1); | |
699 | entry->SetMarkerSize(1); | |
700 | leg->Draw(); | |
701 | } | |
702 | } // dFlag ==2 (LHC data) | |
703 | ||
704 | if(dFlag == 1) { // PDC09 | |
705 | // cut on single muon | |
706 | // theta cut | |
707 | //varmin[THABSMU] = 171+varEpsilon[THABSMU]; | |
708 | //varmax[THABSMU] = 178-varEpsilon[THABSMU]; | |
709 | // eta cut | |
710 | //varmin[ETAMU] = -4.0+varEpsilon[ETAMU]; | |
711 | //varmax[ETAMU] = -2.5-varEpsilon[ETAMU]; | |
712 | ||
713 | // canvas for fitting | |
714 | TCanvas *cfit = new TCanvas("cfit","cfit",0,0,1600,1125); | |
715 | cfit->Divide(7,5); | |
716 | ||
717 | ||
718 | printf("Loop on matrix bins..."); | |
719 | for(Int_t bin0=0; bin0<var0dim; bin0++) { // loop on rapidity (var0dim) | |
720 | varmin[Y] = var0lim[bin0]+varEpsilon[Y]; | |
721 | varmax[Y] = var0lim[bin0+1]-varEpsilon[Y]; | |
722 | printf("%s range: %.2f < %s < %.2f\n",var0name,varmin[Y],var0name,varmax[Y]); | |
723 | ||
724 | for(Int_t bin1=0; bin1<var1dim; bin1++) { // loop on pt (var1dim) | |
725 | varmin[PT] = var1lim[bin1]+varEpsilon[PT]; | |
726 | varmax[PT] = var1lim[bin1+1]-varEpsilon[PT]; | |
727 | printf("%s range: %.2f < %s < %.2f\n",var1name,varmin[PT],var1name,varmax[PT]); | |
728 | ||
729 | cfit->cd((bin0*7)+bin1+1); | |
730 | ||
731 | TH1D* gmass = (TH1D*)genSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
732 | TH1D* rmass = (TH1D*)cintSpar->Slice(IMASS,-1,-1,varmin,varmax); | |
733 | rmass->Rebin(2); | |
734 | gmass->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]); | |
735 | rmass->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]); | |
736 | ||
737 | Double_t genValue, genValueErr, recValue, recValueErr; | |
738 | if(!UseRooFit) { | |
739 | // fit reconstructed resonances to get the #_signal | |
740 | // declare array to retreive fit results | |
741 | Double_t param_pdc[12] = {1,1,1,9.46,0.1,0.04,1,10.023,0.1,1,10.38,0.1}; | |
742 | Double_t *parErr; | |
743 | ||
744 | // fit functions | |
745 | TF1 *func_sig[3], *func_bkg, *func_gl; | |
746 | ||
747 | //___________ | |
748 | // background | |
749 | func_bkg = new TF1("func_bkg",expo,fitrange[0],fitrange[1],2); | |
750 | func_bkg->SetParameters(param_pdc[0],param_pdc[1]); | |
751 | func_bkg->SetLineColor(kRed); | |
752 | func_bkg->SetLineStyle(2); | |
753 | func_bkg->SetLineWidth(2); | |
754 | rmass->Fit(func_bkg,"NR"); | |
755 | func_bkg->GetParameters(¶m_pdc[0]); | |
756 | ||
757 | //_______ | |
758 | // signal | |
759 | // Y(2S)+Y(3S) | |
760 | func_sig[1] = new TF1("func_sig2",normGauss2,9.80,10.60,6); | |
761 | func_sig[1]->SetParameters(param_pdc[6],param_pdc[7],param_pdc[8],param_pdc[9],param_pdc[10],param_pdc[11]); | |
762 | func_sig[1]->SetLineColor(kMagenta); | |
763 | func_sig[1]->SetLineWidth(2); | |
764 | rmass->Fit(func_sig[1],"NR"); | |
765 | func_sig[1]->GetParameters(¶m_pdc[6]); | |
766 | ||
767 | /* | |
768 | func_sig[0] = new TF1("func_sig1",gaussXlandau,9.20,9.70,4); | |
769 | func_sig[0]->SetParameters(param_pdc[2],param_pdc[3],param_pdc[4],param_pdc[5]); | |
770 | func_sig[0]->SetLineColor(kPink); | |
771 | func_sig[0]->SetLineWidth(2); | |
772 | rmass->Fit(func_sig[0],"NR"); | |
773 | func_sig[0]->GetParameters(¶m_pdc[2]); | |
774 | ||
775 | func_sig[2] = new TF1("func_sig3",normGauss,10.32,10.57,3); | |
776 | func_sig[2]->SetParameters(param_pdc[9],param_pdc[10],param_pdc[11]); | |
777 | func_sig[2]->SetLineColor(kViolet); | |
778 | func_sig[2]->SetLineWidth(2); | |
779 | rmass->Fit(func_sig[2],"0"); | |
780 | func_sig[2]->SetParameters(param_pdc[9],param_pdc[10],param_pdc[11]); | |
781 | rmass->Fit(func_sig[2],"R"); | |
782 | func_sig[2]->GetParameters(¶m_pdc[9]); | |
783 | */ | |
784 | ||
785 | // Y(1S)+Y(2S)+Y(3S)+Background | |
786 | func_gl = new TF1("func_gl",global_pdc,fitrange[0],fitrange[1],12); | |
787 | func_gl->SetParameters(param_pdc); | |
788 | func_gl->SetParLimits(3,9.30,9.60); | |
789 | //func_gl->SetParLimits(4,0.05,0.15); | |
790 | func_gl->SetParLimits(7,9.95,10.08); | |
791 | //func_gl->SetParLimits(8,0.05,0.15); | |
792 | func_gl->SetParLimits(10,10.32,10.57); | |
793 | //func_gl->SetParLimits(11,0.05,0.15); | |
794 | func_gl->SetLineColor(kBlue); | |
795 | func_gl->SetLineWidth(2); | |
796 | rmass->Fit(func_gl,"NR"); | |
797 | func_gl->GetParameters(¶m_pdc[0]); | |
798 | parErr = func_gl->GetParErrors(); | |
799 | ||
800 | rmass->Draw("E"); | |
801 | func_bkg->Draw("same"); | |
802 | ||
803 | func_gl->Draw("same"); | |
804 | cfit->Update(); | |
805 | ||
806 | genValue = bin1<5 ? gmass->Integral():gmass->Integral()/binnorm; | |
807 | genValueErr = TMath::Sqrt(genValue); | |
808 | recValue = bin1<5 ? param_pdc[2]/0.1:(param_pdc[2]/0.1)/binnorm; | |
809 | recValueErr = bin1<5 ? parErr[2]/0.1:(parErr[2]/0.1)/binnorm; | |
810 | } | |
811 | else { // RooFit | |
812 | // Define RooFit variable of the x axis value | |
813 | RooRealVar mass("mass","Mass [GeV/c^{2}]",fitrange[0],fitrange[1]); | |
814 | ||
815 | // Convert TH1D class to RooFit | |
816 | RooDataHist hst("hst","hst",RooArgList(mass),rmass); | |
817 | ||
818 | // Define RooFit variable of signal | |
819 | RooRealVar Nsig0("Nsig0","Nsig0",0,1000); | |
820 | RooRealVar Nsig1("Nsig1","Nsig1",0,1000); | |
821 | RooRealVar Nsig2("Nsig2","Nsig2",0,1000); | |
822 | RooRealVar Nbg("Nbg","Nbg",-10,5000); | |
823 | ||
824 | // Define RooFit variable of gaussian for Y(1S) | |
825 | RooRealVar mG0("mG0","mG0",9.30,9.60); | |
826 | RooRealVar sG0("sG0","sG0",0.05,0.2); | |
827 | //RooGenericPdf genG("genG","genG","Nsig/(sG*TMath::Sqrt(2*TMath::Pi())) * TMath::Exp(-((mass-mG)*(mass-mG))/(2*(sG*sG)))",RooArgSet(mass,Nsig,sG,mG)); | |
828 | RooGaussian gauss0("gauss0","gauss0",mass,mG0,sG0); | |
829 | ||
830 | // Define RooFit variable of gaussian for Y(2S) | |
831 | RooRealVar mG1("mG1","mG1",9.80,10.10); | |
832 | RooRealVar sG1("sG1","sG1",0.05,0.2); | |
833 | RooGaussian gauss1("gauss1","gauss1",mass,mG1,sG1); | |
834 | ||
835 | // Define RooFit variable of gaussian for Y(3S) | |
836 | RooRealVar mG2("mG2","mG2",10.30,10.60); | |
837 | RooRealVar sG2("sG2","sG2",0.05,0.2); | |
838 | RooGaussian gauss2("gauss2","gauss2",mass,mG2,sG2); | |
839 | ||
840 | // Define RooFit variable of exponential for background | |
841 | RooRealVar c("c","c",-20,20); | |
842 | RooExponential exp("exp","exp",mass,c); | |
843 | ||
844 | RooAddPdf glpdf("glpdf","glpdf",RooArgList(gauss0,gauss1,gauss2,exp),RooArgList(Nsig0,Nsig1,Nsig2,Nbg)); | |
845 | ||
846 | RooFitResult *fitres = glpdf.fitTo(hst,"mhe0"); | |
847 | ||
848 | // Get a frame from RooRealVar | |
849 | RooPlot *frame = mass.frame(Title("#Upsilon mass plot")); | |
850 | ||
851 | // Get Chi2, mean, sigma and Nsig | |
852 | printf("chi2 = %6.4f\n",frame->chiSquare()); | |
853 | printf("mean (1S) = %6.2f +/- %4.2f [GeV]\n",mG0.getVal(),mG0.getError()); | |
854 | printf("sigma (1S) = %6.2f +/- %4.2f [GeV]\n",sG0.getVal(),sG0.getError()); | |
855 | printf("Nsig (1S) = %6.0f +/- %4.0f\n",Nsig0.getVal(),Nsig0.getError()); | |
856 | genValue = bin1<5 ? gmass->Integral():gmass->Integral()/binnorm; | |
857 | genValueErr = TMath::Sqrt(genValue); | |
858 | recValue = bin1<5 ? Nsig0.getVal():Nsig0.getVal()/binnorm; | |
859 | recValueErr = bin1<5 ? Nsig0.getError():Nsig0.getError()/binnorm; | |
860 | ||
861 | // Add RooDataHist in the frame | |
862 | hst.plotOn(frame); | |
863 | glpdf.plotOn(frame,Components(gauss0),LineColor(kRed),LineStyle(kDashed)); | |
864 | glpdf.plotOn(frame); | |
865 | ||
866 | frame->Draw(); | |
867 | } | |
868 | ||
869 | printf("genValue = %.0f +/- %.0f\n",genValue,genValueErr); | |
870 | printf("recValue = %.0f +/- %.0f\n",recValue,recValueErr); | |
871 | ||
872 | // fill grid | |
873 | Int_t binCell[nVar]={bin0+1,bin1+1}; | |
874 | cintGrid->SetElement(binCell,genValue); | |
875 | cintGrid->SetElementError(binCell,genValueErr); | |
876 | cmusGrid->SetElement(binCell,recValue); | |
877 | cmusGrid->SetElementError(binCell,recValueErr); | |
878 | } | |
879 | } | |
880 | } | |
881 | ||
882 | // set data container | |
883 | dataCont->SetGrid(0, cintGrid); | |
884 | dataCont->SetGrid(1, cmusGrid); | |
885 | ||
886 | TCanvas *cdat = new TCanvas("cdat",Form("%s data container",resname),0,0,1200,400); | |
887 | cdat->Divide(3,1); | |
888 | cdat->cd(1); dataCont->Project(1,0)->Draw(); | |
889 | cdat->cd(2); dataCont->Project(1,1)->Draw(); | |
890 | cdat->cd(3); dataCont->Project(1,0,1)->Draw("colz"); | |
891 | ||
892 | // save data container | |
893 | char *outfile = Form("data.%s",dataname); | |
894 | printf("saving data container as %s...",outfile); | |
895 | dataCont->Save(outfile); | |
896 | printf("done\n"); | |
897 | ||
898 | return outfile; | |
899 | } // extSignal | |
900 | ||
901 | /////////////////// | |
902 | // Do Correction // | |
903 | /////////////////// | |
904 | ||
905 | char* doCorrection(char *dataname, char* file_effmat, char* file_datmat) | |
906 | { | |
907 | ||
908 | // print out starting messages | |
909 | printf("\nNow extracting #signal by fitting...\n"); | |
910 | ||
911 | LoadLib(); | |
912 | ||
913 | TFile *fEffCont = new TFile(file_effmat); | |
914 | TFile *fDataCont = new TFile(file_datmat); | |
915 | ||
916 | AliCFContainer *effCont = (AliCFContainer*)fEffCont->Get("effCont"); | |
917 | AliCFContainer *dataCont = (AliCFContainer*)fDataCont->Get("dataCont"); | |
918 | ||
919 | AliCFEffGrid *matEff = new AliCFEffGrid("matEff",Form("%s efficiency matrix",resname),*effCont); | |
920 | AliCFDataGrid *matGen = new AliCFDataGrid("matGen",Form("%s gen matrix",resname), *dataCont,0); | |
921 | AliCFDataGrid *matData = new AliCFDataGrid("matData",Form("%s data matrix",resname),*dataCont,1); | |
922 | ||
923 | // calculate efficiency | |
924 | matEff->CalculateEfficiency(1,0); | |
925 | ||
926 | //________________________ | |
927 | // display efficiency plot | |
928 | TCanvas *ceff = new TCanvas("ceff",Form("%s efficiency",resname),0,0,1200,400); | |
929 | ceff->Divide(3,1); | |
930 | // efficiency over var0 | |
931 | ceff->cd(1); | |
932 | TH1D *effvar0 = matEff->Project(0); | |
933 | effvar0->SetName(Form("eff_%s",var0name)); | |
934 | effvar0->SetMinimum(0); effvar0->SetMaximum(1); | |
935 | effvar0->SetTitle(Form("%s efficiency(%s)",resname,var0name)); | |
936 | effvar0->GetXaxis()->SetTitle(var0name); | |
937 | effvar0->GetYaxis()->SetTitle("Efficiency"); | |
938 | effvar0->Draw("error"); | |
939 | // efficiency over var1 | |
940 | ceff->cd(2); | |
941 | TH1D *effvar1 = matEff->Project(1); | |
942 | effvar1->SetName(Form("eff_%s",var1name)); | |
943 | effvar1->SetMinimum(0); effvar1->SetMaximum(1); | |
944 | effvar1->SetTitle(Form("%s efficiency(%s)",resname,var1name)); | |
945 | effvar1->GetXaxis()->SetTitle(var1name); | |
946 | effvar1->GetYaxis()->SetTitle("Efficiency"); | |
947 | effvar1->Draw("error"); | |
948 | // 2-dimensional efficiency plot (var0, var1) | |
949 | ceff->cd(3); | |
950 | TH1D *eff2D = matEff->Project(0,1); | |
951 | eff2D->SetName("eff_2D"); | |
952 | eff2D->SetMinimum(0); eff2D->SetMaximum(1); | |
953 | eff2D->SetTitle(Form("%s efficiency(%s,%s)",resname,var0name,var1name)); | |
954 | eff2D->GetXaxis()->SetTitle(var0name); | |
955 | eff2D->GetYaxis()->SetTitle(var1name); | |
956 | eff2D->GetZaxis()->SetTitle("Efficiency"); | |
957 | eff2D->Draw("colz"); | |
958 | //_________________ | |
959 | // display gen plot | |
960 | TCanvas *ccor = new TCanvas("ccor",Form("%s correction",resname),0,0,1200,400); | |
961 | ccor->Divide(3,1); | |
962 | // cora over var0 | |
963 | ccor->cd(1); | |
964 | TH1D *genvar0 = matGen->Project(0); | |
965 | genvar0->SetName(Form("cor_%s",var0name)); | |
966 | genvar0->SetMinimum(0); | |
967 | genvar0->SetLineColor(kRed); | |
968 | genvar0->SetMarkerStyle(22); | |
969 | genvar0->SetMarkerColor(kRed); | |
970 | genvar0->SetTitle(Form("%s correction(%s)",resname,var0name)); | |
971 | genvar0->GetXaxis()->SetTitle(var0name); | |
972 | genvar0->GetYaxis()->SetTitle("dN/dy"); | |
973 | genvar0->Draw(); | |
974 | // cora over var1 | |
975 | ccor->cd(2); | |
976 | TH1D *genvar1 = matGen->Project(1); | |
977 | genvar1->SetName(Form("cor_%s",var1name)); | |
978 | genvar1->SetMinimum(0); | |
979 | genvar1->SetLineColor(kRed); | |
980 | genvar1->SetMarkerStyle(22); | |
981 | genvar1->SetMarkerColor(kRed); | |
982 | genvar1->SetTitle(Form("%s correction(%s)",resname,var1name)); | |
983 | genvar1->GetXaxis()->SetTitle(var1name); | |
984 | genvar1->GetYaxis()->SetTitle("dN/dpt"); | |
985 | genvar1->Draw(); | |
986 | ||
987 | //__________________ | |
988 | // display data plot | |
989 | // data over var0 | |
990 | ccor->cd(1); | |
991 | TH1D *datvar0 = matData->Project(0); | |
992 | datvar0->SetName(Form("dat_%s",var0name)); | |
993 | datvar0->SetMinimum(0); | |
994 | datvar0->SetMarkerStyle(22); | |
995 | datvar0->SetTitle(Form("%s before correction(%s)",resname,var0name)); | |
996 | datvar0->GetXaxis()->SetTitle(var0name); | |
997 | datvar0->GetYaxis()->SetTitle("dN/dy"); | |
998 | datvar0->Draw("same"); | |
999 | // data over var1 | |
1000 | ccor->cd(2); | |
1001 | TH1D *datvar1 = matData->Project(1); | |
1002 | datvar1->SetName(Form("dat_%s",var1name)); | |
1003 | datvar1->SetMinimum(0); | |
1004 | datvar1->SetMarkerStyle(22); | |
1005 | datvar1->SetTitle(Form("%s before correction(%s)",resname,var1name)); | |
1006 | datvar1->GetXaxis()->SetTitle(var1name); | |
1007 | datvar1->GetYaxis()->SetTitle("dN/dpt"); | |
1008 | datvar1->Draw("same"); | |
1009 | ||
1010 | //___________________ | |
1011 | // perform correction | |
1012 | matData->ApplyEffCorrection(*matEff); | |
1013 | ||
1014 | //________________________ | |
1015 | // display correction plot | |
1016 | // correction over var0 | |
1017 | ccor->cd(1); | |
1018 | TH1D *corvar0 = matData->Project(0); | |
1019 | corvar0->SetName(Form("cor_%s",var0name)); | |
1020 | corvar0->SetMinimum(0); | |
1021 | corvar0->SetLineColor(kBlue); | |
1022 | corvar0->SetMarkerStyle(22); | |
1023 | corvar0->SetMarkerColor(kBlue); | |
1024 | corvar0->SetTitle(Form("%s correction(%s)",resname,var0name)); | |
1025 | corvar0->GetXaxis()->SetTitle(var0name); | |
1026 | corvar0->GetYaxis()->SetTitle("dN/dy"); | |
1027 | corvar0->Draw("same"); | |
1028 | // correction over var1 | |
1029 | ccor->cd(2); | |
1030 | TH1D *corvar1 = matData->Project(1); | |
1031 | corvar1->SetName(Form("cor_%s",var1name)); | |
1032 | corvar1->SetMinimum(0); | |
1033 | corvar1->SetLineColor(kBlue); | |
1034 | corvar1->SetMarkerStyle(22); | |
1035 | corvar1->SetMarkerColor(kBlue); | |
1036 | corvar1->SetTitle(Form("%s correction(%s)",resname,var1name)); | |
1037 | corvar1->GetXaxis()->SetTitle(var1name); | |
1038 | corvar1->GetYaxis()->SetTitle("dN/dpt"); | |
1039 | corvar1->Draw("same"); | |
1040 | // 2D correction over var0 and var1 | |
1041 | ccor->cd(3); | |
1042 | TH2D *cor2D = matData->Project(0,1); | |
1043 | cor2D->SetName("cor_2D"); | |
1044 | cor2D->SetMinimum(0); | |
1045 | cor2D->SetTitle(Form("%s correction(%s,%s)",resname,var0name,var1name)); | |
1046 | cor2D->GetXaxis()->SetTitle(var0name); | |
1047 | cor2D->GetYaxis()->SetTitle(var1name); | |
1048 | cor2D->GetZaxis()->SetTitle("dN/dpty"); | |
1049 | cor2D->Draw("colz"); | |
1050 | ||
1051 | ||
1052 | printf("\nCorrection matrix created...\n"); | |
1053 | ||
1054 | // save correction result | |
1055 | char *outfile = Form("correction-result.%s",dataname); | |
1056 | printf("saving correction result as %s...",outfile); | |
1057 | TFile file(outfile,"RECREATE"); | |
1058 | file.cd(); | |
1059 | corvar0->Write(); | |
1060 | corvar1->Write(); | |
1061 | cor2D->Write(); | |
1062 | file.Write(); | |
1063 | file.Close(); | |
1064 | ||
1065 | printf("done\n"); | |
1066 | ||
1067 | return outfile; | |
1068 | } | |
1069 | ||
1070 | ////////////////////////// | |
1071 | // DEFINE FIT FUNCTIONS // | |
1072 | ////////////////////////// | |
1073 | ||
1074 | // single exponential | |
1075 | Double_t expo(Double_t *x, Double_t *par) { | |
1076 | return par[0] * TMath::Exp((par[1]*x[0])); | |
1077 | } | |
1078 | ||
1079 | // double exponential with exclusive region - UPSILON | |
1080 | Double_t dexpo(Double_t *x, Double_t *par) { | |
1081 | // exclusive region | |
1082 | if (x[0] > 9.0 && x[0] < 10.0) { | |
1083 | TF1::RejectPoint(); | |
1084 | return 0; | |
1085 | } | |
1086 | return expo(x, &par[0]) + expo(x, &par[2]); | |
1087 | } | |
1088 | ||
1089 | // double exponential | |
1090 | Double_t dexpo2(Double_t *x, Double_t *par) { | |
1091 | return expo(x, &par[0]) + expo(x, &par[2]); | |
1092 | } | |
1093 | ||
1094 | // signal : gauss | |
1095 | // fit each resonance | |
1096 | Double_t normGauss(Double_t *x, Double_t *par) { | |
1097 | return par[0] / (par[2]*TMath::Sqrt(2*TMath::Pi())) * TMath::Exp(-((x[0]-par[1])*(x[0]-par[1]))/(2*(par[2]*par[2]))); | |
1098 | } | |
1099 | ||
1100 | // fit two resonance | |
1101 | Double_t normGauss2(Double_t *x, Double_t *par) { | |
1102 | return normGauss(x, &par[0]) + normGauss(x,&par[3]); | |
1103 | } | |
1104 | ||
1105 | // fit three resonances | |
1106 | Double_t normGauss3(Double_t *x, Double_t *par) { | |
1107 | return normGauss(x, &par[0]) + normGauss(x, &par[3]) + normGauss(x, &par[6]); | |
1108 | } | |
1109 | ||
1110 | // convolution of gauss and landau | |
1111 | Double_t gaussXlandau(Double_t *var, Double_t *par) { | |
1112 | // Gaussian convoluted with an inverted Landau function: | |
1113 | // TMath::Gaus(x,mean,sigma) | |
1114 | // TMath::Landau(x,mean,width) | |
1115 | ||
1116 | Double_t N = par[0]; // normalization factor = number of entries | |
1117 | Double_t mG = par[1]; // mean of the gaussian | |
1118 | Double_t sG = par[2]; // sigma of the gaussian | |
1119 | Double_t wL = par[3]; // width of the Landau | |
1120 | Double_t x = var[0]; // variable | |
1121 | // Numerical constant | |
1122 | Double_t sqrt2pi = 2.506628275; // Sqrt{2*Pi} | |
1123 | Double_t mL0 = -0.22278298; // maximum Landau location for mean=0 | |
1124 | // Range of integral convolution | |
1125 | Double_t Nsigma = 5; | |
1126 | Double_t Sigma = sG + wL; | |
1127 | //Double_t xMin = mG - Nsigma*Sigma; | |
1128 | //Double_t xMax = mG + Nsigma*Sigma; | |
1129 | Double_t xMin = 0; | |
1130 | Double_t xMax = 2*mG; | |
1131 | // Integral convolution of a Gaussian with an inverted Landau | |
1132 | Int_t Nstep = 2000; | |
1133 | Double_t step = (xMax-xMin) / Nstep; | |
1134 | Double_t sum = 0; | |
1135 | Double_t xx; | |
1136 | for (Double_t i=1.0; i<Nstep; i++) { | |
1137 | xx = xMin + (i-0.5)*step; | |
1138 | sum += TMath::Gaus(xx,mG,sG) * TMath::Landau(-(x-xx),0,wL); | |
1139 | } | |
1140 | // Normalized result | |
1141 | return N*sum*step/(sqrt2pi*sG*wL); | |
1142 | } | |
1143 | ||
1144 | // global fit : UPSILON Family w/ sigle exponential | |
1145 | Double_t global(Double_t *x, Double_t *par) { | |
1146 | return expo(x,&par[0]) + normGauss3(x, &par[2]); | |
1147 | } | |
1148 | ||
1149 | // global fit for PDC09 : UPSILON Family w/ sigle exponential | |
1150 | Double_t global_pdc(Double_t *x, Double_t *par) { | |
1151 | return expo(x,&par[0]) + gaussXlandau(x, &par[2]) + normGauss2(x, &par[6]); | |
1152 | } | |
1153 | ||
1154 | // Crystal ball fit | |
1155 | Double_t CrystalBall(Double_t *x, Double_t *par) { | |
1156 | //Crystal ball function for signal | |
1157 | //parameters are 0:alpha,1:n,2:mean,3:sigma,4:number of expected events | |
1158 | Double_t m=x[0]; | |
1159 | Double_t alpha=par[0]; | |
1160 | Double_t n=par[1]; | |
1161 | Double_t m0=par[2]; | |
1162 | Double_t sigma=par[3]; | |
1163 | Double_t N=par[4]; | |
1164 | ||
1165 | Double_t t = (m-m0)/sigma; | |
1166 | if(alpha<0) t = -t; | |
1167 | ||
1168 | Double_t absAlpha = TMath::Abs((Double_t)alpha); | |
1169 | ||
1170 | if(t >= -absAlpha) { | |
1171 | return N/(sigma*TMath::Sqrt(2*TMath::Pi()))*exp(-0.5*t*t); | |
1172 | //return N*exp(-0.5*t*t); | |
1173 | } | |
1174 | else { | |
1175 | Double_t intn=0; | |
1176 | Double_t fracn=modf(n, &intn); | |
1177 | if(n>10) return 0; | |
1178 | if((n/absAlpha<0) && TMath::Abs(fracn)>0) return 0; // TMath::Power(x,y): negative x and non-integer y not supported!! | |
1179 | Double_t a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha); | |
1180 | Double_t b = n/absAlpha - absAlpha; | |
1181 | ||
1182 | return N/(sigma*TMath::Sqrt(2*TMath::Pi()))*(a/TMath::Power(b - t, n)); | |
1183 | //return N*(a/TMath::Power(b - t, n)); | |
1184 | } | |
1185 | } | |
1186 | ||
1187 | // Loading libraries | |
1188 | void LoadLib() { | |
1189 | printf("Loading libraries..."); | |
1190 | // loading libraries | |
1191 | gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include"); | |
1192 | gSystem->Load("libCore.so"); | |
1193 | gSystem->Load("libTree.so"); | |
1194 | gSystem->Load("libGeom.so"); | |
1195 | gSystem->Load("libVMC.so"); | |
1196 | gSystem->Load("libPhysics.so"); | |
1197 | gSystem->Load("libSTEERBase"); | |
1198 | gSystem->Load("libESD"); | |
1199 | gSystem->Load("libAOD") ; | |
1200 | gSystem->Load("libANALYSIS.so"); | |
1201 | gSystem->Load("libANALYSISalice.so"); | |
1202 | gSystem->Load("libCORRFW.so"); | |
1203 | ||
1204 | gROOT->SetStyle("Plain"); | |
1205 | gStyle->SetPalette(1); | |
1206 | gStyle->SetOptFit(1); | |
1207 | gStyle->SetOptStat(0); | |
1208 | printf("done\n"); | |
1209 | } |