dointegral.C

Go to the documentation of this file.

#include "TGraph.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TAxis.h"
#include <iostream>
#include "TFile.h"
using namespace std;

const double GeV = 1.;
const double TeV = 1000.;

TGraph * PDF;
TGraph  * ECM;
TF1 * higgs_lineshape;

double example(double *var ,double * par)
{
  double x = var[0];
  double x2 = x*x;
  double normalisation = par[0];
  double alpha = par[1];
  double beta= 1-alpha;

  return normalisation*( alpha*3/4.*(1-x2) + (beta)*3/8.*(1+x2) );
  //  return par[0]*( par[1]*3/4.*(1-var[0]*var[0]) + (1-par[1])*3/8.*(1+var[0]*var[0]) );
}

double pdf_func(double * var , double * par)
{
  return PDF->Eval(var[0]);
}

const double smax = 4*7.*TeV * 7*TeV;

double integrand_func(double * var , double * par)
{
  double s = par[0]*par[0];
  double xmin = par[1];
  double x1 = var[0];
  double x2 = s/(smax*x1);

  if (x1<xmin) return 0.;
  else return (1/x1)*pdf_func(&x1,par)*pdf_func(&x2,par);
}

double modified_higgs_lineshape_func(double  * var , double * par)
{
  double roots=var[0];
  return ECM->Eval(roots) * higgs_lineshape->Eval(roots);
}




void  dointegral()
{
  if (!PDF)
     PDF =(TGraph*) (TFile::Open("gluonpdf.root")->Get("PDF"));
  
  //  TF1 * f2 = new TF1("f2",example,-1,1,2);
  // equivalent to 
  // TF1 * f2 = new TF1("f2","[0]*( [1]*3/4.*(1-x*x) + (1-[1])*3/8.*(1+x*x))",1-,1);

  //  double epsilon = 1e-8;
  double epsilon = 1e-10;
  
  TF1 * pdf = new TF1("pdf",pdf_func,epsilon,1,0);
  
  TCanvas * c = new TCanvas("pdfcanvas","pdf distribution of gluon",400,400);
  c->Divide(1,1);
  c->cd()->SetGrid();
  c->cd()->SetLogy();
  pdf->GetXaxis()->SetTitle("x_{gluon}");
  pdf->GetYaxis()->SetTitle("xF(x,Q^{2})");
  pdf->Draw();

  TF1 * integrand = new TF1("integrand",integrand_func,epsilon,1,2);
  integrand->SetParName(0,"#sqrt{s}");
  integrand->SetParName(1,"integral lower limit");
  integrand->SetNpx(5000);
  integrand->SetParameter(0,7*TeV); //7Tev
  integrand->SetParameter(1,0);
  integrand->GetXaxis()->SetTitle("x_{gluon 1}");

  ECM = new TGraph;
  ECM->SetName("ECM");
  

  double roots=10*GeV;
  double roots_step = 10*GeV;
  double roots_max = 14*TeV;
  double 

  int nx =(int) floor( sqrt((roots_max - roots)/roots_step) )+1;

  //c = new TCanvas("integrand","integrand",300,300);
  c = new TCanvas("integrand","integrand",nx*300,nx*300);
  c->Divide(nx,nx);

  TFile * f = TFile::Open("/tmp/ECM.root","recreate");
  f->cd();
  int i=0;

  while(roots<roots_max)
    {
      c->cd(i+1)->SetGrid();
      //      c->cd(i+1)->SetLogy();

      integrand->SetParameter(0,roots);
      double integrand_lower_limit = roots**2/smax;
      integrand->SetParameter(1,integrand_lower_limit);

      cout <<"integrating at "<<roots<<endl;

      double value=integrand->Integral(epsilon,1,(double*)0,1e-7);
      cout <<"done "<<endl;
      ECM->SetPoint(i,roots,value);
      char title[100];
      sprintf(title,"#sqrt{s}_{g-g} = %5.0f [GeV]",roots);
      integrand->SetTitle(title);

      char grname[100];
      sprintf(grname,"integrand_%-5.0f",roots);
      integrand->SetRange(integrand_lower_limit,1.);
      TH1F * copy = (TH1F*) integrand->DrawCopy()->GetHistogram();
      copy->SetName(grname);
      copy->SetXTitle("x_{gluon 1}");
      copy->SetMinimum(integrand->GetMinimum());
      copy->Write();

      cout <<"done for "<<roots<<endl;
      roots+=roots_step;
      i++;
    }
   

  c = new TCanvas("ECMcanvas","gluon-gluon center of mass distribution",400,400);
  c->cd()->SetLogy();
  c->cd()->SetGrid();
  ECM->GetXaxis()->SetTitle("#sqrt{s}_{gluon-gluon}");
  //  ECM->GetYaxis()->SetTitle("P(#sqrt{s}_{gluon-gluon}");
  ECM->GetYaxis()->SetTitle("Arbitrary Unit");
  ECM->Draw("apl");

  ECM->Write();
  f->Close();
}


void multiply(  double massH = 500 *GeV, double widthH = 100 *GeV)
{
  if (!ECM)
    ECM = (TGraph * )TFile::Open("ECM.root")->Get("ECM");

  TCanvas * c = new TCanvas("ECMcanvas","gluon-gluon center of mass distribution",350,1050);
  // c->Divide(1,3);
  // c->cd(1)->SetLogy(); 
  c->cd()->SetLogy();
  // c->cd(1)->SetGrid();
  c->cd()->SetGrid();
  ECM->GetXaxis()->SetTitle("#sqrt{s}_{gluon-gluon}");
  ECM->Draw("apl");

  // c->cd(2)->SetGrid();

  // higgs_lineshape = new TF1("higgs_lineshape","TMath::BreitWigner(x,[0],[1])",massH-4*widthH,massH+4*widthH);
  // higgs_lineshape->SetParameters(massH,widthH);
  // higgs_lineshape->SetNpx(5000);
  // higgs_lineshape->GetXaxis()->SetTitle("Higgs Mass [GeV]");
  // higgs_lineshape->Draw();

  // c->cd(3)->SetGrid();
  
  // TF1 * modified_higgs_lineshape = new TF1("modified_higgs_lineshape",modified_higgs_lineshape_func,massH-4*widthH,massH+4*widthH,0);
  // modified_higgs_lineshape->SetNpx(5000);
  // modified_higgs_lineshape->GetXaxis()->SetTitle("Generated Higgs Mass [GeV]");
  // modified_higgs_lineshape->Draw();
}



   
  

---