]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muon/upsilonCORRFW.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adding cuts for muons
[u/mrichter/AliRoot.git] / PWG / muon / upsilonCORRFW.C
CommitLineData
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
36Bool_t IsShowContOn = kFALSE;
37// check bin error
38Bool_t IsCheckBinError = kFALSE;
39
40//_______________________
41// User setting variables
42
43// enumerate variables
44enum { NEVENT, Y, PT, IMASS, TRIG, PTMU, PMU, TRIGSIDE, RABSMU, CHARGE, ETAMU, THABSMU, VZMU, DCAMU };
45
46// set efficiency matrix bin
47const Char_t *var0name="y";
48const Char_t *var1name="pt";
49const Int_t var0bin = 5;
50const Int_t var1bin = 7;
51Double_t var0lim[var0bin+1] = {-4.0,-3.7,-3.4,-3.1,-2.8,-2.5};
52Double_t var1lim[var1bin+1] = {0,2,4,6,8,10,15,20};
53const Double_t binnorm = (Double_t)5/2; // for the last two pt bin
54
55//______________
56// Analysis cuts
57
58// Use pt trigger cut
59Bool_t IsHpt = kTRUE;
60
61//_______________
62// fit parameters
63Bool_t UseRooFit = kFALSE;
64if(UseRooFit) { using namespace RooFit; }
65
66// rebinning
67Bool_t rebinned = kTRUE;
68// 3 gauss for UPSILON family
69Char_t *resname = "Upsilon";
70Double_t fitrange[2] = {7, 12}; // fit range min, max
71Double_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)
74Double_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
76Bool_t IsFixSigma = kTRUE;
77
78////////////////////////////
79// UPSILON ANALYSIS MACRO //
80////////////////////////////
81
82//______________________________________________________
83
84Bool_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
141char* 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
366char* 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(&param[0]);
549
550 TF1 *gl = new TF1("gl",global,fitrange[0],fitrange[1],11);
551 gl->SetParameters(&param[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(&param[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(&param_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(&param_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(&param_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(&param_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(&param_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
905char* 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
1075Double_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
1080Double_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
1090Double_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
1096Double_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
1101Double_t normGauss2(Double_t *x, Double_t *par) {
1102 return normGauss(x, &par[0]) + normGauss(x,&par[3]);
1103}
1104
1105// fit three resonances
1106Double_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
1111Double_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
1145Double_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
1150Double_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
1155Double_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
1188void 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}