10 #include "TPaveText.h"
16 #include "TGraphErrors.h"
17 #include "TVirtualFitter.h"
20 #include "TSpectrum.h"
23 #include "RooGenericPdf.h"
24 #include "RooFitResult.h"
25 #include "RooRealVar.h"
28 ///////////////////////////
29 // Default Options (Global)
30 ///////////////////////////
35 // display container contents
36 Bool_t IsShowContOn = kFALSE;
38 Bool_t IsCheckBinError = kFALSE;
40 //_______________________
41 // User setting variables
43 // enumerate variables
44 enum { NEVENT, Y, PT, IMASS, TRIG, PTMU, PMU, TRIGSIDE, RABSMU, CHARGE, ETAMU, THABSMU, VZMU, DCAMU };
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
63 Bool_t UseRooFit = kFALSE;
64 if(UseRooFit) { using namespace RooFit; }
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;
78 ////////////////////////////
79 // UPSILON ANALYSIS MACRO //
80 ////////////////////////////
82 //______________________________________________________
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
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");
101 char *file_effmat = ""; // efficiency matrix
102 char *file_datmat = ""; // data matrix
103 char *file_result = ""; // correction result
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);
113 // call a function to perform fit to extract #nignal
114 if(IsExtSig) file_datmat = extSignal(dataname);
116 file_datmat = "data.container.CFMuonResUpsilon.PDC09.local.AOD.MC.root";
117 printf("Default data container: %s\n",file_datmat);
120 // call a function to perform a correction on data
121 if(IsDoCorrection) file_result = doCorrection(dataname, file_effmat, file_datmat);
123 // print out closing messages
124 printf("******************************************************\n");
125 printf(" Finishing Correction Framework for Upsilon...\n");
126 printf("******************************************************\n");
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);
137 //////////////////////////////
138 // Create Efficiency matrix //
139 //////////////////////////////
141 char* createEffMat(const char *dataname)
144 // print out starting messages
145 printf("\nNow starting efficiency matrix creation...\n");
149 // open input container
150 printf("Opening container from %s...",dataname);
151 TFile *file = new TFile(dataname);
152 AliCFContainer *cont = (AliCFContainer*) (file->Get("container"));
155 // print out container step and variables
158 Int_t nstep = cont->GetNStep();
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
175 TCanvas *csc[stepnum];
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();
189 // GridSparse from container
190 AliCFGridSparse *genSpar = (AliCFGridSparse*)cont->GetGrid(0); // GEN
191 AliCFGridSparse *recSpar = (AliCFGridSparse*)cont->GetGrid(1); // REC
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
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);
214 //____________________
215 // Loop on matrix bins
217 // cut on single muon
219 varmin[THABSMU] = 171+varEpsilon[THABSMU];
220 varmax[THABSMU] = 178-varEpsilon[THABSMU];
222 varmin[ETAMU] = -4.0+varEpsilon[ETAMU];
223 varmax[ETAMU] = -2.5-varEpsilon[ETAMU];
225 //varmin[RABSMU] = 17.6+varEpsilon[RABSMU];
226 //varmin[RABSMU] = 80.0-varEpsilon[RABSMU];
228 // canvas for fitting
229 TCanvas *cfit = new TCanvas("cfit","cfit",0,0,1600,1125);
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]);
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]);
243 TH1D* gmass = (TH1D*)genSpar->Slice(IMASS,-1,-1,varmin,varmax);
244 TH1D* rmass = (TH1D*)recSpar->Slice(IMASS,-1,-1,varmin,varmax);
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
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);
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]);
264 func->SetParameters(par);
265 rmass->Fit(func,"RL+");
266 // retreive fit parameters
267 func->GetParameters(&par[0]);
268 parErr = func->GetParErrors();
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;
277 printf("genValue = %.0f +/- %.0f\n",genValue,genValueErr);
278 printf("recValue = %.0f +/- %.0f\n",recValue,recValueErr);
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);
288 printf("saving fit plots as %s...",Form("recfitplots.%s",dataname));
289 cfit->SaveAs(Form("recfitplots.%s",dataname));
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);
304 cvars->cd(1); gvar0->Draw();
305 cvars->cd(2); gvar1->Draw();
306 cvars->cd(3); rvar0->Draw();
307 cvars->cd(4); rvar1->Draw();
309 // save efficiency container
310 char *outfile = Form("efficiency.%s",dataname);
311 printf("saving efficiency container as %s...",outfile);
312 effCont->Save(outfile);
315 //_________________________
316 // create efficiency matrix
317 printf("creating efficiency matrix...");
318 AliCFEffGrid *matEff = new AliCFEffGrid("matEff",Form("%s efficiency matrix",resname),*effCont);
320 // calcualte efficiency
321 printf("calculating efficiency...");
322 matEff->CalculateEfficiency(1,0); // REC over MC
325 //________________________
326 // display efficiency plot
327 TCanvas *ceff = new TCanvas("ceff",Form("%s efficiency",resname),0,0,1200,400);
329 // efficiency over var0
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
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)
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");
358 printf("\nEfficiency matrix created...\n");
366 char* extSignal(const char *dataname)
369 // print out starting messages
370 printf("\nNow extracting #signal by fitting...\n");
374 // open input container
375 printf("Opening container from %s...",dataname);
376 TFile *file = new TFile(dataname);
377 AliCFContainer *cont = (AliCFContainer*) (file->Get("container"));
380 // check dataname : PDC, LHC or Unknown (1, 2, 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; }
386 // print out container step and variables
388 Int_t nstep = cont->GetNStep();
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
405 TCanvas *csc[stepnum];
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();
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);
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
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
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);
454 //_____________________
455 // Loops on matrix bins
457 if(dFlag == 2) { // LHC data
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);
467 // common cut on events
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];
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);
479 TH1D *hmassi[3], *hmassm[3];
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++) {
488 varmin[TRIG] = 0+varEpsilon[TRIG];
491 varmin[TRIG] = 3+varEpsilon[TRIG];
494 varmin[TRIG] = 3.3+varEpsilon[TRIG];
501 varmin[TRIG] = 0+varEpsilon[TRIG];
504 varmin[TRIG] = 2+varEpsilon[TRIG];
507 varmin[TRIG] = 2.2+varEpsilon[TRIG];
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);
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);
530 // draw N matching mass plot and fit
536 for(Int_t i=0; i<3; i++) {
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;
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;
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]);
550 TF1 *gl = new TF1("gl",global,fitrange[0],fitrange[1],11);
551 gl->SetParameters(¶m[0]);
553 // fix yield 2S and 3S normalized from CMS result
554 gl->FixParameter(5,1.8); gl->FixParameter(8,1);
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]);}
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]);}
564 gl->SetLineColor(kBlue);
566 hmassm[i]->Fit(gl,"NR");
567 gl->GetParameters(¶m[0]);
568 parErr = gl->GetParErrors();
570 bkg->SetParameters(param[0],param[1]);
571 bkg->SetParErrors(parErr);
572 bkg->SetLineColor(kRed);
573 bkg->SetLineStyle(2);
574 bkg->SetLineWidth(2);
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);
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);
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);
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);
600 c1[i] = new TCanvas(Form("c1%d",i),Form("c1%d",i),0,0,800,800);
602 c1[i]->cd(1); //c1[i]->cd(1)->SetLogy();
603 hmassm[i]->SetMinimum(0);
604 hmassm[i]->DrawCopy();
606 //bkg_2->Draw("same");
613 chi2 = gl->GetChisquare();
615 mean = sig->GetParameter(1);
616 mean_err = sig->GetParError(1);
617 sigma = sig->GetParameter(2);
618 sigma_err = sig->GetParError(2);
621 nsig = sig->Integral(imin, imax);
622 nsig_err = TMath::Sqrt(nsig);
623 //nsig = sig->GetParameter(0);
624 //nsig_err = sig->GetParError(0);
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;
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;
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)));
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));
653 info->SetFillColor(0);
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);
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);
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);
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);
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);
702 } // dFlag ==2 (LHC data)
704 if(dFlag == 1) { // PDC09
705 // cut on single muon
707 //varmin[THABSMU] = 171+varEpsilon[THABSMU];
708 //varmax[THABSMU] = 178-varEpsilon[THABSMU];
710 //varmin[ETAMU] = -4.0+varEpsilon[ETAMU];
711 //varmax[ETAMU] = -2.5-varEpsilon[ETAMU];
713 // canvas for fitting
714 TCanvas *cfit = new TCanvas("cfit","cfit",0,0,1600,1125);
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]);
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]);
729 cfit->cd((bin0*7)+bin1+1);
731 TH1D* gmass = (TH1D*)genSpar->Slice(IMASS,-1,-1,varmin,varmax);
732 TH1D* rmass = (TH1D*)cintSpar->Slice(IMASS,-1,-1,varmin,varmax);
734 gmass->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]);
735 rmass->GetXaxis()->SetRangeUser(usrrange[0],usrrange[1]);
737 Double_t genValue, genValueErr, recValue, recValueErr;
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};
745 TF1 *func_sig[3], *func_bkg, *func_gl;
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]);
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]);
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]);
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]);
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();
801 func_bkg->Draw("same");
803 func_gl->Draw("same");
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;
812 // Define RooFit variable of the x axis value
813 RooRealVar mass("mass","Mass [GeV/c^{2}]",fitrange[0],fitrange[1]);
815 // Convert TH1D class to RooFit
816 RooDataHist hst("hst","hst",RooArgList(mass),rmass);
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);
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);
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);
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);
840 // Define RooFit variable of exponential for background
841 RooRealVar c("c","c",-20,20);
842 RooExponential exp("exp","exp",mass,c);
844 RooAddPdf glpdf("glpdf","glpdf",RooArgList(gauss0,gauss1,gauss2,exp),RooArgList(Nsig0,Nsig1,Nsig2,Nbg));
846 RooFitResult *fitres = glpdf.fitTo(hst,"mhe0");
848 // Get a frame from RooRealVar
849 RooPlot *frame = mass.frame(Title("#Upsilon mass plot"));
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;
861 // Add RooDataHist in the frame
863 glpdf.plotOn(frame,Components(gauss0),LineColor(kRed),LineStyle(kDashed));
869 printf("genValue = %.0f +/- %.0f\n",genValue,genValueErr);
870 printf("recValue = %.0f +/- %.0f\n",recValue,recValueErr);
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);
882 // set data container
883 dataCont->SetGrid(0, cintGrid);
884 dataCont->SetGrid(1, cmusGrid);
886 TCanvas *cdat = new TCanvas("cdat",Form("%s data container",resname),0,0,1200,400);
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");
892 // save data container
893 char *outfile = Form("data.%s",dataname);
894 printf("saving data container as %s...",outfile);
895 dataCont->Save(outfile);
905 char* doCorrection(char *dataname, char* file_effmat, char* file_datmat)
908 // print out starting messages
909 printf("\nNow extracting #signal by fitting...\n");
913 TFile *fEffCont = new TFile(file_effmat);
914 TFile *fDataCont = new TFile(file_datmat);
916 AliCFContainer *effCont = (AliCFContainer*)fEffCont->Get("effCont");
917 AliCFContainer *dataCont = (AliCFContainer*)fDataCont->Get("dataCont");
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);
923 // calculate efficiency
924 matEff->CalculateEfficiency(1,0);
926 //________________________
927 // display efficiency plot
928 TCanvas *ceff = new TCanvas("ceff",Form("%s efficiency",resname),0,0,1200,400);
930 // efficiency over var0
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
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)
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");
960 TCanvas *ccor = new TCanvas("ccor",Form("%s correction",resname),0,0,1200,400);
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");
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");
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");
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");
1010 //___________________
1011 // perform correction
1012 matData->ApplyEffCorrection(*matEff);
1014 //________________________
1015 // display correction plot
1016 // correction over var0
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
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
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");
1052 printf("\nCorrection matrix created...\n");
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");
1070 //////////////////////////
1071 // DEFINE FIT FUNCTIONS //
1072 //////////////////////////
1074 // single exponential
1075 Double_t expo(Double_t *x, Double_t *par) {
1076 return par[0] * TMath::Exp((par[1]*x[0]));
1079 // double exponential with exclusive region - UPSILON
1080 Double_t dexpo(Double_t *x, Double_t *par) {
1082 if (x[0] > 9.0 && x[0] < 10.0) {
1086 return expo(x, &par[0]) + expo(x, &par[2]);
1089 // double exponential
1090 Double_t dexpo2(Double_t *x, Double_t *par) {
1091 return expo(x, &par[0]) + expo(x, &par[2]);
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])));
1100 // fit two resonance
1101 Double_t normGauss2(Double_t *x, Double_t *par) {
1102 return normGauss(x, &par[0]) + normGauss(x,&par[3]);
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]);
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)
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;
1130 Double_t xMax = 2*mG;
1131 // Integral convolution of a Gaussian with an inverted Landau
1133 Double_t step = (xMax-xMin) / Nstep;
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);
1140 // Normalized result
1141 return N*sum*step/(sqrt2pi*sG*wL);
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]);
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]);
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
1159 Double_t alpha=par[0];
1162 Double_t sigma=par[3];
1165 Double_t t = (m-m0)/sigma;
1168 Double_t absAlpha = TMath::Abs((Double_t)alpha);
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);
1176 Double_t fracn=modf(n, &intn);
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;
1182 return N/(sigma*TMath::Sqrt(2*TMath::Pi()))*(a/TMath::Power(b - t, n));
1183 //return N*(a/TMath::Power(b - t, n));
1187 // Loading libraries
1189 printf("Loading libraries...");
1190 // loading libraries
1191 gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
1192 gSystem->Load("libCore");
1193 gSystem->Load("libTree");
1194 gSystem->Load("libGeom");
1195 gSystem->Load("libVMC");
1196 gSystem->Load("libPhysics");
1197 gSystem->Load("libSTEERBase");
1198 gSystem->Load("libESD");
1199 gSystem->Load("libAOD") ;
1200 gSystem->Load("libANALYSIS");
1201 gSystem->Load("libANALYSISalice");
1202 gSystem->Load("libCORRFW");
1204 gROOT->SetStyle("Plain");
1205 gStyle->SetPalette(1);
1206 gStyle->SetOptFit(1);
1207 gStyle->SetOptStat(0);