PWG/muon/upsilonCORRFW.C

   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(&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
 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 }