Avant toute chose, veuillez installer ROOT ou configurer VirtualBox pour vous placer dans un environnement où vous pouvez l'utiliser. Pour installer VirtualBox, vous pouvez télécharger les deux fichiers suivants (attention, le fichier d'image disque virtuelle .vhd fait 7 Go).
Sauvegardez les deux fichiers dans le même répertoire. Si vous avez installé VirtualBox (en téléchargeant le programme à partir du site), il suffit de lancer le fichier .vbox pour démarrer votre machine virtuelle.
Dans tous les cas, testez que votre installation de ROOT fonctionne. Pour ce faire, veuillez suivre les instructions contenues dans ce document.
#include <iostream>
#include <TH1.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
using namespace std;
void test(){
cout << "Hello world!" << endl;
TH1F *histoTest = new TH1F("histoTest","histo;x;y",10,0,10);
histoTest->Fill(2.53);
histoTest->Fill(5.24);
histoTest->Fill(2.02);
histoTest->Fill(8.93);
TF1 *funcTest = new TF1("funcTest","[0] + [1]*x",0.,10.);
funcTest->SetParameters(1,1);
histoTest->Fit("funcTest");
histoTest->Draw();
}
#include <iostream>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
using namespace std;
// Nous allons générer des nombres aléatoires grâce à
// l algorithme de congruence linéaire.
// La structure du programme sera la suivante:
// - une fonction principale: genCL
// - déclaration de toutes les variables et histogrammes (voir http://root.cern.ch/root/html/TH1F.html)
// - une boucle pour générer les nombres aléatoires
// - on dessine les histos et on les sauve dans un fichier .root
// On définit la fonction genCL qui prend un argument: npoints
// Par défaut, npoints = 10000.
// Lorsqu on appelle la fonction, on peut lui donner l argument
// que l on veut: genCL(1000000) par exemple.
void genCL(int npoints = 10000)
{
// Pour afficher du texte à l ecran, on utilise: cout << "Mon texte" << endl;
cout << "Generateur congruentiel lineaire" << endl;
// Déclarer toutes les variables nécessaires ici:
int nBinX = 100, nBinY = 100;
double xlow = 0., xup = 1., ylow = 0., yup = 1.;
// Déclarer les histogrammes ici.
// Pour déclarer un histogramme, la syntaxe est la suivante:
TH1F *monHisto = new TH1F("monHisto", "Titre complet de mon histo;Titre axe x;Titre axe y", nBinX, xlow, xup);
// Pour un histogramme a 2 dim, la syntaxe devient:
TH2F *monHisto2D = new TH2F("monHisto2D", "Titre complet de mon histo;Titre axe x;Titre axe y", nBinX, xlow, xup, nBinY, ylow, yup);
// On introduit la boucle qui va générer les npoints nombres aléatoires
for (int i = 0; i < npoints; i++) {
// Ecrivez ici la formule de congruence linéaire
// Remplissez les histogrammes. La syntaxe est la suivante:
monHisto->Fill(0.5);
// Pour un histo 2D, cela devient:
monHisto2D->Fill(0.5, 0.5);
}// fin de la boucle
// Le reste du code est destiné à afficher les histogrammes.
// On réserve un canvas que l on divise en deux.
TCanvas *can = new TCanvas();
can->Divide(1, 2);
// Dans la 1ère partie, dessinez la distribution des nombres aléatoires
can->cd(1);
monHisto->SetMinimum(0.0);
monHisto->Draw();
// Dans la 2ème partie, dessinez l histogramme de corrélation
can->cd(2);
monHisto2D->Draw();
}
#include <iostream>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
using namespace std;
// Nous allons générer des nombres aléatoires grâce à
// l algorithme de congruence linéaire.
// La structure du programme sera la suivante:
// - une fonction principale: genCL
// - déclaration de toutes les variables et histogrammes (voir http://root.cern.ch/root/html/TH1F.html)
// - une boucle pour générer les nombres aléatoires
// - on dessine les histos et on les sauve dans un fichier .root
// On définit la fonction genCL qui prend un argument: npoints
// Par défaut, npoints = 10000.
// Lorsqu on appelle la fonction, on peut lui donner l argument
// que l on veut: genCL(1000000) par exemple.
void genCL(int npoints = 10000)
{
// Pour afficher du texte à l ecran, on utilise: cout << "Mon texte" << endl;
cout << "Generateur congruentiel lineaire" << endl;
// Déclarer toutes les variables nécessaires ici:
int nBinX = 100, nBinY = 100;
double xlow = 0., xup = 1., ylow = 0., yup = 1.;
double xCoor = 0., yCoor = 0.;
int a = 205, c = 29573, m = 139968, last = 4711; // La variable last représente I_0, la dernière valeur calculée.
// Déclarer les histogrammes ici.
// Pour déclarer un histogramme, la syntaxe est la suivante:
// TH1F *monHisto = new TH1F("monHisto", "Titre complet de mon histo;Titre axe x;Titre axe y", nBin, xmin, xmax);
// Pour un histogramme a 2 dim, la syntaxe devient:
// TH2F *monHisto2D = new TH2F("monHisto2D", "Titre complet de mon histo;Titre axe x;Titre axe y", nBinX, xmin, xmax, nBinY, ymin, ymax);
TH1F *h_Dist = new TH1F("h_Dist", "Distribution selon un generateur congruentiel lineaire;R_{j};N", nBinX, xlow, xup);
TH2F *h_Corr = new TH2F("h_Corr", "Correlation n, n+1;R_{j};R_{j+1}", nBinX, xlow, xup, nBinY, ylow, yup);
// On introduit la boucle qui va générer les npoints nombres aléatoires
for (int i = 0; i < npoints; i++) {
// Ecrivez ici la formule de congruence linéaire
last = (a * last + c) % m;
xCoor = last * 1. / m;
last = (a * last + c) % m;
yCoor = last * 1. / m;
// Remplissez les histogrammes. La syntaxe est la suivante:
h_Dist->Fill(xCoor);
// Pour un histo 2D, cela devient:
h_Corr->Fill(xCoor, yCoor);
}// fin de la boucle
// Le reste du code est destiné à afficher les histogrammes.
// On réserve un canvas que l on divise en deux.
TCanvas *can = new TCanvas();
can->Divide(1, 2);
// Dans la 1ère partie, dessinez la distribution des nombres aléatoires
can->cd(1);
h_Dist->SetMinimum(0.0);
h_Dist->Draw("e");
// Dans la 2ème partie, dessinez l histogramme de corrélation
can->cd(2);
h_Corr->Draw();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// Nous allons generer des nombres aleatoires grace au
// generateur de la classe TRandom3.
// La structure du programme reste la meme:
// - une fonction principale: genTRandom3
// - declaration de toutes les variables et histogrammes (voir http://root.cern.ch/root/html/TH1F.html)
// - une boucle pour generer les nombres aleatoires
// - on dessine les histos
// On definit la fonction genTRandom3 qui prend un argument: npoints
// Par defaut, npoints = 10000.
// Lorsqu on appelle la fonction, on peut lui donner l argument
// que l on veut: genTRandom3(1000000) par exemple.
void genTRandom3(int npoints = 10000)
{
// Declarez toutes les variables necessaires ici:
// Declarez les histogrammes, pour la distribution sur [0,1],
// pour la correlation (2D), et pour la distribution [-1000,1000]
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
// NB: on genere un nombre aleatoire uniforme sur [0,1]
// grace a la methode Rndm() de la classe TRandom3
//generateur.SetSeed(time(NULL));
// On introduit la boucle qui va generer les npoints nombres aleatoires
for (int i=0; i < npoints; i++) {
// Remplissez les histogrammes pour l exercice a)
// Generer entre -1000 et 1000 et remplissez l histogramme pour l exercice b)
}// fin de la boucle
// Reservez un 1er canvas pour la distribution sur [0, 1] et la correlation.
cout << "Generateur TRandom3 " << endl;
TCanvas *can1 = new TCanvas();
can1->Divide(1, 2);
// Dessinez le 1er histo
can1->cd(1);
// Dessinez l histo de correlation
can1->cd(2);
// Reservez un 2eme canvas pour la distribution sur (-1000, 1000).
TCanvas *can2 = new TCanvas("can2", "can2", 100, 0, 700, 500);
// Dessinez la distribution sur (-1000, 1000)
can2->cd();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// Nous allons generer des nombres aleatoires grace au
// generateur de la classe TRandom3.
// La structure du programme reste la meme:
// - une fonction principale: genTRandom3
// - declaration de toutes les variables et histogrammes (voir http://root.cern.ch/root/html/TH1F.html)
// - une boucle pour generer les nombres aleatoires
// - on dessine les histos
// On definit la fonction genTRandom3 qui prend un argument: npoints
// Par defaut, npoints = 10000.
// Lorsqu on appelle la fonction, on peut lui donner l argument
// que l on veut: genTRandom3(1000000) par exemple.
void genTRandom3(int npoints = 10000)
{
// Declarez toutes les variables necessaires ici:
int nBinX = 100, nBinY = 100;
double xlow = 0., xup = 1., ylow = 0., yup = 1., xCoor = 0., yCoor = 0.;
// Declarez les histogrammes, pour la distribution sur [0,1],
// pour la correlation (2D)
TH1F *h_Dist = new TH1F("h_Dist", "Distribution selon TRandom3;R_{j};N", nBinX, xlow, xup);
TH2F *h_Corr = new TH2F("h_Corr", "Correlation entre 2 nombres TRandom3;R_{j};R_{j+1}",nBinX, xlow, xup, nBinY, ylow, yup);
// Il faut maintenant declarer le generateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
// NB: on genere un nombre aleatoire uniforme sur [0,1]
// grace a la methode Rndm() de la classe TRandom3
//generateur.SetSeed(time(NULL));
generateur.SetSeed(123456);
// On introduit la boucle qui va generer les npoints nombres aleatoires
for (int i=0; i < npoints; i++) {
// Remplissez les histogrammes pour l exercice a)
xCoor = generateur.Rndm();
yCoor = generateur.Rndm();
h_Dist->Fill(xCoor);
h_Corr->Fill(xCoor, yCoor);
}// fin de la boucle
// Reservez un 1er canvas pour la distribution sur [0, 1] et la correlation.
cout << "Generateur TRandom3 " << endl;
TCanvas *can1 = new TCanvas();
can1->Divide(1, 2);
// Dessinez le 1er histo
can1->cd(1);
h_Dist->SetMinimum(0.0);
h_Dist->Draw("e");
// Dessinez l histo de correlation
can1->cd(2);
h_Corr->Draw();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// Toujours grace a la classe TRandom3, nous allons generer
// des nombres aleatoires suivant une loi Normale de moyenne mu
// et de variance sigma2.
// La structure du programme est la suivante:
// - une fonction principale: loiNormale (avec deux arguments: mu et sigma2)
// - declaration de toutes les variables et histogrammes
// - une boucle pour generer les nombres aleatoires
// - on dessine les histos
// On definit la fonction loiNormale qui prend deux arguments: mu ( = 0 par defaut) et sigma2 ( = 1 par defaut)
void loiNormale(double mu = 0, double sigma2 = 1)
{
// Declarez toutes les variables necessaires ici:
// nombre de points, xlow, xup...
int npoints = 5000;
// Declarez les histogrammes
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// La boucle commence ici, pour chaque passage, deux nombres
// suivant la loi Normale seront generes
for (int i = 0; i < npoints; i++) {
// Generez deux nombres aleatoires uniformes sur [0,1]
// Construisez les deux nombres gaussiens a partir des deux precedents
// Modifiez le code pour une gaussienne de moyenne mu et de variance sigma2
// Remplissez l'histogramme
}// fin de la boucle
cout << "Loi Normale - N(" << mu << ", " << sigma2 << ")" << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez l'histogramme
// Pour la partie b), fittez l'histogramme avec une gaussienne
// grace a la methode Fit() de la classe TH1.
// ( https://root.cern.ch/root/htmldoc/guides/users-guide/FittingHistograms.html )
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// Toujours grace a la classe TRandom3, nous allons generer
// des nombres aleatoires suivant une loi Normale de moyenne mu
// et de variance sigma2.
// La structure du programme est la suivante:
// - une fonction principale: loiNormale (avec deux arguments: mu et sigma2)
// - declaration de toutes les variables et histogrammes
// - une boucle pour generer les nombres aleatoires
// - on dessine les histos
// On definit la fonction loiNormale qui prend deux arguments: mu ( = 0 par defaut) et sigma2 ( = 1 par defaut)
void loiNormale(double mu = 0, double sigma2 = 1)
{
// Declarez toutes les variables necessaires ici:
// nombre de points, xlow, xup...
int npoints = 5000;
int nbins = 100;
double pi = 3.141592, xlow = -8., xup = 8., U1 = 0., U2 = 0., T1 = 0., T2 = 0.;
// Declarez les histogrammes
TH1F *h_Dist = new TH1F("h_Dist", "Distribution Normale;R_{J};N", nbins, xlow, xup);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// La boucle commence ici, pour chaque passage, deux nombres
// suivant la loi Normale seront generes
for (int i = 0; i < npoints; i++) {
// Generez deux nombres aleatoires uniformes sur [0,1]
U1 = generateur.Rndm();
U2 = generateur.Rndm();
// Construisez les deux nombres gaussiens a partir des deux precedents
T1 = sqrt(-2 * log(U1)) * cos(2 * pi * U2);
T2 = sqrt(-2 * log(U1)) * sin(2 * pi * U2);
// Modifiez le code pour une gaussienne de moyenne mu et de variance sigma2
T1 = mu + T1 * sqrt(sigma2);
T2 = mu + T2 * sqrt(sigma2);
// Remplissez l'histogramme
h_Dist->Fill(T1);
h_Dist->Fill(T2);
}// fin de la boucle
cout << "Loi Normale - N(" << mu << ", " << sigma2 << ")" << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez l'histogramme
h_Dist->SetMinimum(0.0);
h_Dist->DrawCopy();
// Pour la partie b), fittez l'histogramme avec une gaussienne
// grace a la methode Fit() de la classe TH1.
// ( https://root.cern.ch/root/htmldoc/guides/users-guide/FittingHistograms.html )
h_Dist->Fit("gaus");
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
// result = .......
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[2];
arg[0] = x;
par[0] = 1.;
par[1] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
// Elle prend en argument tau, qui est par défaut mis à 1.
void invTransform(double tau = 1.)
{
cout << "Methode de la transformée inverse" << endl;
// Nombre de points
int npoints = 100000;
double u = 0., x = 0., xlow = 0., xup = 10.;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exp(-x);x;N", 200, xlow, xup);
// Il faut maintenant declarer le generateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle for.
for (int i = 0; i < npoints; i++) {
// Generez ici les variables aleatoires necessaires
}// fin de la boucle for
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Draw();
// On declare un objet function TF1 basé sur notre
// fonction myFunc. La fonction TF1 est définie sur (0, 10)
// et possède 2 parametres: normalisation et temps de vie
TF1 *f = new TF1("f", myFunc, 0, 10, 2);
// On donne la valeur 1 a la normalisation et tau au temps de vie
f->SetParameters(1., tau);
// On repondère notre histogramme pour que
// l'intégrale soit égale à 1
// (plus facile d'interpréter les résultats du fit)
h_Dist->Scale(1. / h_Dist->Integral("width"));
// ("width" indique qu'on multiplie par la largeur du bin.)
h_Dist->Fit("f");
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
result = par[0] * (1. / par[1]) * exp(-x[0] / par[1]);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[2];
arg[0] = x;
par[0] = 1.;
par[1] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
// Elle prend en argument tau, qui est par défaut mis à 1.
void invTransform(double tau = 1.)
{
cout << "Methode de la transformée inverse" << endl;
// Nombre de points
int npoints = 100000;
double u = 0., x = 0., xlow = 0., xup = 10.;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exp(-x);x;N", 200, xlow, xup);
// Il faut maintenant declarer le generateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle for.
for (int i = 0; i < npoints; i++) {
// Generez ici les variables aleatoires necessaires
u = generateur.Rndm();
x = - tau * log(1-u);
h_Dist->Fill(x);
}// fin de la boucle for
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Draw();
// On declare un objet function TF1 basé sur notre
// fonction myFunc. La fonction TF1 est définie sur (0, 10)
// et possède 2 parametres: normalisation et temps de vie
TF1 *f = new TF1("f", myFunc, 0, 10, 2);
// On donne la valeur 1 a la normalisation et tau au temps de vie
f->SetParameters(1., tau);
// On repondère notre histogramme pour que
// l'intégrale soit égale à 1
// (plus facile d'interpréter les résultats du fit)
h_Dist->Scale(1. / h_Dist->Integral("width"));
// ("width" indique qu'on multiplie par la largeur du bin.)
h_Dist->Fit("f");
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
// result = .......
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[2];
arg[0] = x;
par[0] = 0.375;
par[1] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 100000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., xlow = -1., xup = 1., fmax = 0.75;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution 1 + x^{2} - methode de rejection", 200, xlow, xup);
// Il faut maintenant declarer le generateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
}// fin du while
cout << "Efficacite: " << "......" << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->Draw();
TF1 *f = new TF1("f", myFunc, -1, 1, 2);
f->SetParameters(0.375, 1.);
h_Dist->Scale(1. / h_Dist->Integral("width"));
h_Dist->Fit("f");
h_Dist->SetMinimum(0.);
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
result = par[0] * (1. + par[1] * x[0] * x[0]);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[2];
arg[0] = x;
par[0] = 0.375;
par[1] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 100000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., xlow = -1., xup = 1., fmax = 0.75;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution 1 + x^{2} - methode de rejection", 200, xlow, xup);
// Il faut maintenant declarer le generateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
x = xlow + generateur.Rndm() * (xup - xlow);
y = generateur.Rndm() * fmax;
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
if (y < myFunc(x)) {
h_Dist->Fill(x);
nAccepted++;
}
nTotal++;
}// fin du while
cout << "Efficacite : " << 1.*nAccepted/nTotal << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->Draw();
TF1 *f = new TF1("f", myFunc, -1, 1, 2);
f->SetParameters(0.375, 1.);
h_Dist->Scale(1. / h_Dist->Integral("width"));
h_Dist->Fit("f");
h_Dist->SetMinimum(0.);
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
// result = .......
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x, double y)
{
double arg[2], par[4];
arg[0] = x;
arg[1] = y;
par[0] = 1./12.;
par[1] = 1.;
par[2] = 3.;
par[3] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., z = 0., xlow = -1., xup = 1., ylow = 0., yup = 3., fmax = 0.5;
TH2F *h_Dist = new TH2F("h_Dist", "Distribution (3-y)(1+x^{2}) - meth rej;x;y", 200, xlow, xup, 30, ylow, yup);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
// Testez si le point (x, y, z) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
}// fin du while
cout << "Efficacite : " << "......" << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
TF2 *f = new TF2("f", myFunc, -1., 1., 0., 3., 4);
f->SetParameters(800, 1., 3., 1.);
h_Dist->Fit("f");
h_Dist->Draw("lego");
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = -1.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
result = par[0] * (par[2] - par[3] * x[1]) * (1. + par[1] * x[0] * x[0]);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x, double y)
{
double arg[2], par[4];
arg[0] = x;
arg[1] = y;
par[0] = 1./12.;
par[1] = 1.;
par[2] = 3.;
par[3] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., z = 0., xlow = -1., xup = 1., ylow = 0., yup = 3., fmax = 0.5;
TH2F *h_Dist = new TH2F("h_Dist", "Distribution (3-y)(1+x^{2}) - meth rej;x;y", 200, xlow, xup, 30, ylow, yup);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
x = xlow + generateur.Rndm() * (xup - xlow);
y = ylow + generateur.Rndm() * (yup - ylow);
z = generateur.Rndm() * fmax;
// Testez si le point (x, y, z) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
if (z < myFunc(x, y)) {
h_Dist->Fill(x, y);
nAccepted++;
}
nTotal++;
}// fin du while
cout << "Efficacite : " << nAccepted * 1. / nTotal << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
TF2 *f = new TF2("f", myFunc, -1., 1., 0., 3., 4);
f->SetParameters(800, 1., 3., 1.);
h_Dist->Fit("f");
h_Dist->Draw("lego");
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = 0.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
// result = .......
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[1];
arg[0] = x;
par[0] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., xlow = 0, xup = 4, fmax = 7./3.;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exercice 3 - meth rej", 200, xlow, xup);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
}// fin du while
cout << "Efficacite : " << "......" << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Scale(1./h_Dist->Integral("width"));
TF1 *f = new TF1("f", myFunc, 0., 4., 1);
f->SetParameter(0, 1.);
f->SetNpx(10000);
h_Dist->Fit("f");
h_Dist->Draw();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = 0.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
if (0 <= x[0] && x[0] <= 1) result = par[0] * (1./2.) * pow(x[0], 2);
else if (x[0] <= 2) result = par[0] * (1./2.) * pow(x[0] - 2, 2);
else if (x[0] <= 3) result = par[0] * (7./3.) * pow(x[0] - 2, 6);
else if (x[0] <= 4) result = par[0] * (7./3.) * pow(x[0] - 4, 6);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[1];
arg[0] = x;
par[0] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void rejMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
double x = 0., y = 0., xlow = 0, xup = 4, fmax = 7./3.;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exercice 3 - meth rej", 200, xlow, xup);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generateur;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generateur.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Generez ici les variables aleatoires necessaires
x = xlow + generateur.Rndm() * (xup - xlow);
y = generateur.Rndm() * fmax;
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
if (y < myFunc(x)) {
h_Dist->Fill(x);
nAccepted++;
}
nTotal++;
}// fin du while
cout << "Efficacite : " << nAccepted * 1. / nTotal << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Scale(1./h_Dist->Integral("width"));
TF1 *f = new TF1("f", myFunc, 0., 4., 1);
f->SetParameter(0, 1.);
f->SetNpx(10000);
h_Dist->Fit("f");
h_Dist->Draw();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = 0.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
if (0 <= x[0] && x[0] <= 1) result = par[0] * (1./2.) * pow(x[0], 2);
else if (x[0] <= 2) result = par[0] * (1./2.) * pow(x[0] - 2, 2);
else if (x[0] <= 3) result = par[0] * (7./3.) * pow(x[0] - 2, 6);
else if (x[0] <= 4) result = par[0] * (7./3.) * pow(x[0] - 4, 6);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[1];
arg[0] = x;
par[0] = 1.;
return myFunc(arg, par);
}
// Voici la fonction principale.
void twoBoxesMethode()
{
cout << "Methode de rejection avec deux boites" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exercice 3 - meth rej", 200, 0, 4);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generator;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generator.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
// Il faut d'abord tirer un nombre pour selectionner
// la boite (petite ou grande boite)
double u0 = generator.Rndm();
double x, xlow, xhigh, fmax, y;
// En fonction de u0, on choisit la petite ou la grande boite
// et on assigne les valeurs correspondantes à xlow, xhigh et fmax
if (u0 < 3./17.) { // alors on est dans le petit rectangle sur (0, 2)
xlow = 0;
xhigh = 2;
fmax = 1./2.;
}
else { // alors on est dans le grand rectangle sur (2, 4)
xlow = 2;
xhigh = 4;
fmax = 7./3.;
}
// on peut maintenant parcourir uniformément
// la boite selectionnée ce tour-ci
x = xlow + (xhigh - xlow) * generator.Rndm();
// nombre aleatoire pour y entre 0 et fmax
y = fmax * generator.Rndm();
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
if (y < myFunc(x)) {
h_Dist->Fill(x);
nAccepted++;
}
nTotal++;
}// fin du while
cout << "Efficacite: " << nAccepted * 1. / nTotal << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Scale(1./h_Dist->Integral("width"));
TF1 *f = new TF1("f", myFunc, 0., 4., 1);
f->SetParameter(0, 1.);
f->SetNpx(10000);
h_Dist->Fit("f");
h_Dist->Draw();
}
#include <iostream>
#include <time.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TRandom3.h>
using namespace std;
// La fonction myFunc sera utilisee pour fitter la distribution
// que nous allons generer. Il faut donc y inserer la forme de
// la fonction a fitter
// ( http://root.cern.ch/root/html/TF1.html#F3 )
double myFunc(double *x, double *par)
{
double result = 0.;
// Ecrivez ici la forme de la pdf avec les parametres a fitter:
// result = par[0] * (1 + par[1] * x[0]) par exemple
if (0 <= x[0] && x[0] <= 1) result = par[0] * (1./2.) * pow(x[0], 2);
else if (x[0] <= 2) result = par[0] * (1./2.) * pow(x[0] - 2, 2);
else if (x[0] <= 3) result = par[0] * (7./3.) * pow(x[0] - 2, 6);
else if (x[0] <= 4) result = par[0] * (7./3.) * pow(x[0] - 4, 6);
return result;
}
// Les lignes suivantes permettent d'utiliser la fonction
// myFunc avec des parametres deja initialises.
// VOUS NE DEVEZ RIEN Y CHANGER
double myFunc(double x)
{
double arg[1], par[1];
arg[0] = x;
par[0] = 1.;
return myFunc(arg, par);
}
// fonction qui retourne les valeurs
// triangulaires pour un triangle sur (0, 2)
double h1(double x)
{
double result;
if (x < 0) result = 0;
else if (x < 1) result = x;
else if (x < 2) result = 2 - x;
else result = 0;
return result;
}
// fonction qui retourne les valeurs
// triangulaires pour un triangle sur (2, 4)
double h2(double x)
{
double result;
if (x < 2) result = 0;
else if (x < 3) result = x - 2;
else if (x < 4) result = 4 - x;
else result = 0;
return result;
}
// Voici la fonction principale.
void twoTrianglesMethode()
{
cout << "Methode de rejection" << endl;
// Nombre de points et nombre de points acceptes
int npoints = 10000000, nAccepted = 0, nTotal = 0;
TH1F *h_Dist = new TH1F("h_Dist", "Distribution exercice 3 - meth rej", 200, 0, 4);
// Il faut maintenant declarer le genrateur TRandom3.
TRandom3 generator;
// Decommentez l instruction suivante pour avoir differents
// resultats a chaque fois que vous lancez le programme.
//generator.SetSeed(time(NULL));
// Debut de la boucle. Attention, c'est une boucle while.
// On effectue les instructions JUSQU'A CE QUE le nombre
// de points acceptes = npoints
while (nAccepted < npoints) {
double u0 = generator.Rndm();
double x, hx, xlow, xhigh, fmax, y;
if (u0 < 3./17.) {
xlow = 0;
xhigh = 2;
fmax = 1./2.;
// on genere un chiffre triangulairement distribué sur
// [0, 2]. La somme de deux nombres uniformes a cette propriété.
x = generator.Rndm() + generator.Rndm();
// on calcul la hauteur d'une distribution triangulaire
// calculé au x tiré
hx = h1(x);
}
else {
xlow = 2;
xhigh = 4;
fmax = 7./3.;
// on genere un chiffre triangulairement distribué sur
// [2, 4]. La somme de deux nombres uniformes, décalée de 2, a cette propriété.
x = 2 + generator.Rndm() + generator.Rndm();
// on calcule la hauteur d'une distribution triangulaire
// calculée au x tiré
hx = h2(x);
}
// nombre aleatoire pour y entre 0 et fmax
y = fmax * hx * generator.Rndm();
// Testez si le point (x, y) genere est accepte ou non.
// Si il est accepte, remplissez l histogramme et incrementez
// le compteur
if (y < myFunc(x)) {
h_Dist->Fill(x);
nAccepted++;
}
nTotal++;
}// fin du while
cout << "Efficacite: " << nAccepted * 1. / nTotal << endl;
TCanvas *can = new TCanvas();
can->cd();
// Dessinez et fittez la distribution obtenue
h_Dist->SetMinimum(0.);
h_Dist->Scale(1./h_Dist->Integral("width"));
TF1 *f = new TF1("f", myFunc, 0., 4., 1);
f->SetParameter(0, 1);
f->SetNpx(10000);
h_Dist->Fit("f");
h_Dist->Draw();
}
#include <iostream>
#include <time.h>
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TRandom3.h"
using namespace std;
// Dans cet exercice nous allons simuler la reaction e+ e- -> q qbar
// Dans ce programme nous allons construire plusieurs fonctions:
// - myFunc qui servira a fitter
// - dsigma qui calculera la valeur de la section efficace
// différentielle à partir des valeurs de theta générées
// - rejMethode la fonction principale pour simuler et remplir
// des histogrammes
// On commence par definir quelques constantes.
#define eBeam 5.
#define pi 3.141592
#define alphaEM 0.007297352
// On definit ici la fonction qui sera utilisee pour le fit
// (voir exercices 2 a et b de la seance 2)
double myFunc(double *x, double *par)
{
double result = 0.;
//result = ... ;
return result;
}
// On definit ici la formule pour la section efficace différentielle,
// obtenue après toutes les simplifications
double dsigma(double theta, double Qq)
{
// Codez ici la formule de la section efficace.
double s = 0; // à changer
int Nc = 0; // à changer
double dSigma = 0;
//dSigma = ... ;
return dSigma;
}
// On definit ici la fonction de generation des evenements.
void rejMethode()
{
// Introduisez la charge du quark ici (on suppose les masses nulles).
// Vous pouvez aussi paramétriser la fonction rejMethode pour qu'elle
// prenne en argument la saveur du quark.
double Qq = 2./3;
// On definit le nombre d evenements souhaites et les histos
int npoints = 100000;
TH1F *h_theta = new TH1F("h_theta", "d#sigma /d#Omega;#theta;d#sigma/d#Omega", 100, 0., pi);
//TH1F ..... ;
// Declarez les variables a generer et autres variables necessaires
double costheta;
// Initialisez le generateur
TRandom3 generateur;
// Décommentez l instruction suivante pour changer la seed
// du générateur de facon aleatoire.
//generateur.SetSeed(time(NULL));
// Calculez dsigmaMax
double dsigmaMax = dsigma(0., Qq);
cout << "dsigmaMax = " << dsigmaMax << endl;
// Compteurs
int nAccepted = 0;
int nTotal = 0;
// On commence la boucle ici.
while (nAccepted < npoints) {
// Generez les variables aleatoires uniformes necessaires.
// Testez si l evenement genere est acceptable
// et remplissez ce qu il faut
if ( true ) { //à changer
h_theta->Fill(1.); //à changer
nAccepted++;
}
nTotal++;
}// fin de la boucle while
// Imprimez l efficacite de la methode de rejection
cout << "Efficacité : " << nAccepted << "/" << nTotal << " = " << nAccepted * 1. / nTotal << endl;
// Calculez et imprimez la section efficace intégrée
double sigmaTot = 0;
// sigmaTot = .....;
// Declarez une fonction et fittez la distribution
// en cos(theta)
TF1 *f1 = new TF1("f1", myFunc, -1., 1., 2);
f1->SetParameters(100, 1.);
// Affichez les histogrammes dans un canvas
TCanvas *can = new TCanvas();
// Lorsque vous aurez d'autres histogrammes, divisez le canevas pour en dessiner d'autres
//can->Divide(2,2);
can->cd(1);
h_theta->SetMinimum(0.);
h_theta->DrawCopy();
//h_theta->Fit("f1"); // Ne fonctionnera pas (pourquoi ?)
// Sauvez le canevas en format pdf
can->SaveAs("test.pdf");
}
#include <iostream>
#include <time.h>
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TRandom3.h"
using namespace std;
// Dans cet exercice nous allons simuler la reaction e+ e- -> q qbar
// Dans ce programme nous allons construire plusieurs fonctions:
// - myFunc qui servira a fitter
// - dsigma qui calculera la valeur de la section efficace
// différentielle à partir des valeurs de theta générées
// - rejMethode la fonction principale pour simuler et remplir
// des histogrammes
// On commence par definir quelques constantes.
#define eBeam 5.
#define pi 3.141592
#define alphaEM 0.007297352
// On definit ici la fonction qui sera utilisee pour le fit
// (voir exercices 2 a et b de la seance 2)
double myFunc(double *x, double *par)
{
double result = 0.;
double s = 4. * eBeam * eBeam;
result = par[0] * (1. + par[1] * x[0] * x[0]);
return result;
}
// On definit ici la fonction obtenue apres toutes les simplifications.
double dsigma(double theta, double Qq)
{
double s = 4. * eBeam * eBeam;
int Nc = 3;
double dSigma = alphaEM * alphaEM * Nc * Qq * Qq / (4 * s);
dSigma *= (1 + cos(theta) * cos(theta) );
return dSigma;
}
// On definit ici la fonction de generation des evenements.
void rejMethode()
{
// Introduisez la charge du quark ici (on suppose les masses nulles).
// Vous pouvez aussi paramétriser la fonction rejMethode pour qu'elle
// prenne en argument la saveur du quark.
double Qq = 2./3;
// On definit le nombre d evenements souhaites et les histos
int npoints = 100000;
TH1F *h_theta = new TH1F("h_theta", "d#sigma /d#Omega;#theta;N", 100, 0., pi);
TH1F *h_cosTheta = new TH1F("h_cosTheta", "d#sigma /d#Omega;cos#theta;N", 100, -1., 1.);
//TH1F *h_theta0 = new TH1F("h_theta0", "d#sigma /d#Omega;#theta0;N", 100, 0., pi);
//TH1F *h_cosTheta0 = new TH1F("h_cosTheta0", "d#sigma /d#Omega;cos#theta0;N", 100, -1., 1.);
TH1F *h_phi = new TH1F("h_phi", "d#sigma /d#Omega;#phi;N", 100, 0., 2*pi);
TH1F *h_pT = new TH1F("h_pT", "d#sigma /d#Omega;p_{T}[GeV];N", 100, 0., 1.2*eBeam);
// Declarez les variables a generer et autres variables necessaires
double costheta, phi, z1;
double theta, pT;
// Initialisez le generateur
TRandom3 generateur;
// Décommentez l instruction suivante pour changer la seed
// du générateur de facon aleatoire.
//generateur.SetSeed(time(NULL));
// Calculez dsigmaMax
double dsigmaMax = dsigma(0., Qq);
cout << "dsigmaMax = " << dsigmaMax << endl;
// Compteurs
int nAccepted = 0;
int nTotal = 0;
// On commence la boucle ici.
while (nAccepted < npoints) {
// Generez les variables aleatoires uniformes necessaires.
costheta = -1. + generateur.Rndm() * 2.;
theta = acos(costheta);
//h_theta0->Fill(theta);
//h_cosTheta0->Fill(costheta);
pT = eBeam*sin(theta); // Quantité de mouvement transverse du quark
phi = generateur.Rndm() * 2 * pi;
z1 = generateur.Rndm() * dsigmaMax;
// Testez si l evenement genere est acceptable
// et remplissez ce qu il faut
if ( z1 < dsigma(theta, Qq) ) {
h_theta->Fill(theta);
h_cosTheta->Fill(costheta);
h_phi->Fill(phi);
h_pT->Fill(pT);
nAccepted++;
}
nTotal++;
}// fin de la boucle while
// Imprimez l efficacite de la methode de rejection
cout << "Efficacité : " << nAccepted << "/" << nTotal << " = " << nAccepted * 1. / nTotal << endl;
// Calculez et imprimez la section efficace intégrée
double sigmaTot = 4./3. * pi * alphaEM * alphaEM * 3. * Qq * Qq / 4 / eBeam / eBeam;
cout << "Section efficace intégrée : " << sigmaTot << " GeV^{-2} ou " << sigmaTot*3.8707e8 << " pb. Comme on a généré " << nAccepted << " événements, cela donne une luminosité générée de " << nAccepted / sigmaTot / 3.8707e8 << " pb^{-1}. C'est " << nAccepted / sigmaTot / 3.8707e8 / 162.1 << " fois la luminosité totale enregistrée par l'expérience DELPHI entre 1989 et 1995." << endl;
// Declarez une fonction et fittez la distribution
// en cos(theta)
TF1 *f1 = new TF1("f1", myFunc, -1., 1., 2);
f1->SetParameters(100, 1.);
// Affichez les histogrammes dans un canvas
TCanvas *can = new TCanvas();
can->Divide(2,2);
can->cd(1);
h_theta->SetMinimum(0.);
h_theta->DrawCopy();
can->cd(2);
h_cosTheta->SetMinimum(0.);
h_cosTheta->DrawCopy();
h_cosTheta->Fit("f1");
can->cd(3);
//h_theta0->SetMinimum(0.);
//h_theta0->DrawCopy();
h_phi->SetMinimum(0.);
h_phi->DrawCopy();
can->cd(4);
//h_cosTheta0->SetMinimum(0.);
//h_cosTheta0->DrawCopy();
h_pT->SetMinimum(0.);
h_pT->DrawCopy();
// Sauvez le canevas en format pdf
can->SaveAs("eeqqbar_0mass.pdf");
//can->SaveAs("eeqqbar_0mass_test21mars2023.pdf");
}
#include <iostream>
#include <time.h>
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TRandom3.h"
using namespace std;
// Dans cet exercice nous allons simuler la reaction e+ e- -> q qbar
// Dans ce programme nous allons construire plusieurs fonctions:
// - myFunc qui servira a fitter
// - dsigma qui calculera la valeur de la section efficace
// a partir des elements de matrices et des quadri-vecteurs
// que nous allons generer
// - rejMethode la fonction principale pour simuler et remplir
// des trees et histogrammes
// On commence par definir quelques constantes.
#define eBeam 5.
#define pi 3.141592
#define alphaEM 0.007297352
// On definit ici la fonction qui sera utilisee pour le fit
// (voir exercices 2 a et b de la seance 2)
double myFunc(double *x, double *par)
{
double result = 0.;
result = ...;
return result;
}
// On definit ici la formule de la section efficace, obtenue apres
// toutes les simplifications, simple fonction de qMass, Qq et theta.
// On ne l'utilisera que pour calculer sigma max.
double dsigma( ... )
{
// Codez ici la formule de la section efficace en fonction de qMass, Qq, theta, ...
double s = ...;
int Nc = ...;
double dSigma = ...;
return dSigma;
}
// On definit ici la fonction qui sera utilisee pour calculer
// la section efficace en fonction des quadri-vecteurs qu elle
// recoit en argument.
double dsigma(TLorentzVector ele, TLorentzVector pos, TLorentzVector q, TLorentzVector qbar, ... )
{
// Codez ici la formule de la section efficace en fonction de quadri-vecteurs, ...
double s = ...;
int Nc = ...;
double M2 = ...;
double dSigma = ...;
return dSigma;
}
// On definit ici la fonction de generation des evenements.
void rejMethode(TString quark = "")
{
// Les quelques lignes suivantes permettent de tenir compte de
// la masse des quarks.
// quarks table
// u ~ 0.003 GeV d ~ 0.006 GeV
// s ~ 0.105 GeV c ~ 1.25 GeV
// b ~ 4.25 GeV t ~ 173 GeV
double qMass = 0.0;
double Qq = 2./3;
if (quark == "u"){
qMass = 0.003;
Qq = 2./3;
}
else if (quark == "d"){
qMass = 0.006;
Qq = -1./3;
}
else if (quark == "s"){
qMass = 0.105;
Qq = -1./3;
}
else if (quark == "c"){
qMass = 1.25;
Qq = 2./3;
}
else if (quark == "b"){
qMass = 4.25;
Qq = -1./3;
}
else if (quark == "t"){
qMass = 173.;
Qq = 2./3;
}
// On definit le nombre d evenements souhaites.
// Le fichier .root, le TTree et un histo
int npoints = 100000;
TString fileName = "exseance3_" + quark + "_quark.root";
TFile f(fileName, "RECREATE");
TTree *tree = new TTree("tree", "MC root Tree");
TH1F *h_dsigma = new TH1F("h_dsigma", "d#sigma /d#Omega", 100, -1., 1.);
// On definit ici les variables qui seront sauvees dans le TTree
// Suivez l exemple pour ajouter d autres variables
float q_Px, q_Py, q_Pz, q_E, ele_Px, ele_Py, ele_Pz, ele_E;
float qbar_Px, qbar_Py, qbar_Pz, qbar_E, pos_Px, pos_Py, pos_Pz, pos_E;
float q_Phi, q_CosTheta;
tree->Branch("ele_Px", &ele_Px, "ele_Px/F");
tree->Branch("ele_Py", &ele_Py, "ele_Py/F");
tree->Branch("ele_Pz", &ele_Pz, "ele_Pz/F");
tree->Branch("ele_E", &ele_E, "ele_E/F");
tree->Branch("pos_Px", &pos_Px, "pos_Px/F");
tree->Branch("pos_Py", &pos_Py, "pos_Py/F");
tree->Branch("pos_Pz", &pos_Pz, "pos_Pz/F");
tree->Branch("pos_E", &pos_E, "pos_E/F");
tree->Branch("q_Px", &q_Px, "q_Px/F");
tree->Branch("q_Py", &q_Py, "q_Py/F");
tree->Branch("q_Pz", &q_Pz, "q_Pz/F");
tree->Branch("q_E", &q_E, "q_E/F");
tree->Branch("qbar_Px", &qbar_Px, "qbar_Px/F");
tree->Branch("qbar_Py", &qbar_Py, "qbar_Py/F");
tree->Branch("qbar_Pz", &qbar_Pz, "qbar_Pz/F");
tree->Branch("qbar_E", &qbar_E, "qbar_E/F");
tree->Branch("q_Phi", &q_Phi, "q_Phi/F");
tree->Branch("q_CosTheta", &q_CosTheta, "q_CosTheta/F");
// On definit ici les quadri-vecteurs des particules de la reaction.
// Suivez l exemple pour les autres particules
TLorentzVector ele, pos, q, qbar;
ele.SetPxPyPzE(0., 0., eBeam, eBeam);
pos.SetPxPyPzE(0., 0., -eBeam, eBeam);
q.SetPxPyPzE(0., 0., 0., 0.);
qbar.SetPxPyPzE(0., 0., 0., 0.);
// Declarez les variables a generer et autres variables necessaires
// Initialisez le generateur
TRandom3 generateur;
// Décommentez l instruction suivante pour changer la seed
// du gnérateur de facon aleatoire.
//generateur.SetSeed(time(NULL));
// Calculez dsigmaMax
double dsigmaMax = dsigma(0., qMass, Qq);
cout << "dsigmaMax = " << dsigmaMax << endl;
// Compteurs
int nAccepted = 0;
int nTotal = 0;
// On commence la boucle ici.
while (nAccepted < npoints) {
// Generez les variables aleatoires uniformes necessaires.
// Construisez les quadri-vecteurs des quarks
// Testez si l evenement genere est acceptable
// et remplissez ce qu il faut
if ( ........ ){
h_dsigma->Fill( ...... );
tree->Fill();
nAccepted++;
}
nTotal++;
}// fin de la boucle while
// Imprimez l efficacite de la methode de rejection
cout << nAccepted << "/" << nTotal << " = " << nAccepted * 1. / nTotal << endl;
// Declarez une fonction et fittez la distribution
// en cos(theta)
TF1 *f1 = new TF1("f1", myFunc, -1., 1., 2);
f1->SetParameters(100, 1.);
h_dsigma->Fit("f1");
// Affichez l histogramme dans un canvas
TCanvas *can = new TCanvas();
can->cd();
h_dsigma->SetMinimum(0.);
h_dsigma->DrawCopy();
// Sauvez le TTree ainsi que l histogramme
tree->Write();
h_dsigma->Write();
f.Close();
}
#include <iostream>
#include <time.h>
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TRandom3.h"
using namespace std;
// Dans cet exercice nous allons simuler la reaction e+ e- -> q qbar
// Dans ce programme nous allons construire plusieurs fonctions:
// - myFunc qui servira a fitter
// - dsigma qui calculera la valeur de la section efficace
// a partir des elements de matrices et des quadri-vecteurs
// que nous allons generer
// - rejMethode la fonction principale pour simuler et remplir
// des trees et histogrammes
// On commence par definir quelques constantes.
#define eBeam 5.
#define pi 3.141592
#define alphaEM 0.007297352
// On definit ici la fonction qui sera utilisee pour le fit
// (voir exercices 2 a et b de la seance 2)
double myFunc(double *x, double *par)
{
double result = 0.;
double s = 4. * eBeam * eBeam;
result = par[0] * ((1. + 4 * par[1] / s) + (1. - 4 * par[1] / s) * x[0] * x[0]);
return result;
}
// On definit ici la formule de la section efficace, obtenue apres
// toutes les simplifications, simple fonction de qMass, Qq et theta.
// On ne l'utilisera que pour calculer sigma max.
double dsigma(double theta, double qMass, double Qq)
{
double s = 4. * eBeam * eBeam;
int Nc = 3;
double dSigma = alphaEM * alphaEM * Nc * Qq * Qq / (4 * s);
dSigma *= sqrt(1 - (4 * qMass * qMass / s));
dSigma *= ((1 + (4 * qMass * qMass / s)) + (1 - (4 * qMass * qMass / s)) * cos(theta) * cos(theta));
return dSigma;
}
// On definit ici la fonction qui sera utilisee pour calculer
// la section efficace en fonction des quadri-vecteurs qu elle
// recoit en argument.
double dsigma(TLorentzVector ele, TLorentzVector pos, TLorentzVector q, TLorentzVector qbar, double qMass, double Qq)
{
// Codez ici la formule de la section efficace en fonction de quadri-vecteurs, ...
double s = 4.*eBeam*eBeam;
int Nc = 3;
double M2 = 8 * 4 * pi * 4 * pi * alphaEM * alphaEM * Nc * Qq * Qq / (s * s);
double bracket = (ele * q) * (pos * qbar) + (ele * qbar) * (pos * q) + qMass * qMass * (ele * pos);
M2 *= bracket;
double dSigma = 2 * sqrt( s / 4. - qMass * qMass) / sqrt(s) / (8 * pi * 8 * pi * s);
dSigma *= M2;
return dSigma;
}
// On definit ici la fonction de generation des evenements.
void rejMethode(TString quark = "u")
{
// Les quelques lignes suivantes permettent de tenir compte de
// la masse des quarks.
// quarks table
// u ~ 0.003 GeV d ~ 0.006 GeV
// s ~ 0.105 GeV c ~ 1.25 GeV
// b ~ 4.25 GeV t ~ 173 GeV
double qMass = 0.0;
double Qq = 2./3;
if (quark == "u"){
qMass = 0.003;
Qq = 2./3;
}
else if (quark == "d"){
qMass = 0.006;
Qq = -1./3;
}
else if (quark == "s"){
qMass = 0.105;
Qq = -1./3;
}
else if (quark == "c"){
qMass = 1.25;
Qq = 2./3;
}
else if (quark == "b"){
qMass = 4.25;
Qq = -1./3;
}
else if (quark == "t"){
qMass = 173.;
Qq = 2./3;
}
// On definit le nombre d evenements souhaites.
// Le fichier .root, le TTree et un histo
int npoints = 100000;
TString fileName = "exseance3_" + quark + "_quark.root";
TFile f(fileName, "RECREATE");
TTree *tree = new TTree("tree", "MC root Tree");
TH1F *h_theta = new TH1F("h_theta", "d#sigma /d#Omega;#theta;d#sigma/d#Omega", 100, 0., 3.141592);
TH1F *h_cosTheta = new TH1F("h_cosTheta", "d#sigma /d#Omega;cos#theta;d#sigma/d#Omega", 100, -1., 1.);
TH1F *h_phi = new TH1F("h_phi", "d#sigma /d#Omega;#phi;d#sigma/d#Omega", 100, 0., 2*3.141592);
TH1F *h_pT = new TH1F("h_pT", "d#sigma /d#Omega;p_{T}[GeV];d#sigma/d#Omega", 100, 0., 1.*eBeam);
// On definit ici les variables qui seront sauvees dans le TTree
// Suivez l exemple pour ajouter d autres variables
float q_Px, q_Py, q_Pz, q_E, ele_Px, ele_Py, ele_Pz, ele_E;
float qbar_Px, qbar_Py, qbar_Pz, qbar_E, pos_Px, pos_Py, pos_Pz, pos_E;
float q_Phi, q_CosTheta;
tree->Branch("ele_Px", &ele_Px, "ele_Px/F");
tree->Branch("ele_Py", &ele_Py, "ele_Py/F");
tree->Branch("ele_Pz", &ele_Pz, "ele_Pz/F");
tree->Branch("ele_E", &ele_E, "ele_E/F");
tree->Branch("pos_Px", &pos_Px, "pos_Px/F");
tree->Branch("pos_Py", &pos_Py, "pos_Py/F");
tree->Branch("pos_Pz", &pos_Pz, "pos_Pz/F");
tree->Branch("pos_E", &pos_E, "pos_E/F");
tree->Branch("q_Px", &q_Px, "q_Px/F");
tree->Branch("q_Py", &q_Py, "q_Py/F");
tree->Branch("q_Pz", &q_Pz, "q_Pz/F");
tree->Branch("q_E", &q_E, "q_E/F");
tree->Branch("qbar_Px", &qbar_Px, "qbar_Px/F");
tree->Branch("qbar_Py", &qbar_Py, "qbar_Py/F");
tree->Branch("qbar_Pz", &qbar_Pz, "qbar_Pz/F");
tree->Branch("qbar_E", &qbar_E, "qbar_E/F");
tree->Branch("q_Phi", &q_Phi, "q_Phi/F");
tree->Branch("q_CosTheta", &q_CosTheta, "q_CosTheta/F");
// On definit ici les quadri-vecteurs des particules de la reaction.
// Suivez l exemple pour les autres particules
TLorentzVector ele, pos, q, qbar;
ele.SetPxPyPzE(0., 0., eBeam, eBeam);
pos.SetPxPyPzE(0., 0., -eBeam, eBeam);
q.SetPxPyPzE(0., 0., 0., 0.);
qbar.SetPxPyPzE(0., 0., 0., 0.);
// Declarez les variables a generer et autres variables necessaires
double costheta, phi, z1;
double theta, k2, kx, ky, kz;
k2 = eBeam * eBeam - qMass * qMass;
// Initialisez le generateur
TRandom3 generateur;
// Décommentez l instruction suivante pour changer la seed
// du gnérateur de facon aleatoire.
//generateur.SetSeed(time(NULL));
// Calculez dsigmaMax
double dsigmaMax = dsigma(0., qMass, Qq);
cout << "dsigmaMax = " << dsigmaMax << endl;
// Compteurs
int nAccepted = 0;
int nTotal = 0;
// On commence la boucle ici.
while (nAccepted < npoints) {
// Generez les variables aleatoires uniformes necessaires.
// Construisez les quadri-vecteurs des quarks
costheta = -1. + generateur.Rndm() * 2.;
theta = acos(costheta);
//theta = 3.141592*generateur.Rndm();
//costheta = cos(theta);
phi = generateur.Rndm() * 2 * pi;
z1 = generateur.Rndm() * dsigmaMax;
kx = sqrt(k2) * sin(theta) * cos(phi);
ky = sqrt(k2) * sin(theta) * sin(phi);
kz = sqrt(k2) * cos(theta);
q.SetPxPyPzE(kx, ky, kz, eBeam);
qbar.SetPxPyPzE(-kx, -ky, -kz, eBeam);
ele_Px = ele.Px(); ele_Py = ele.Py(); ele_Pz = ele.Pz(); ele_E = ele.E();
pos_Px = pos.Px(); pos_Py = pos.Py(); pos_Pz = pos.Pz(); pos_E = pos.E();
q_Px = q.Px(); q_Py = q.Py(); q_Pz = q.Pz(); q_E = q.E();
qbar_Px = qbar.Px(); qbar_Py = qbar.Py(); qbar_Pz = qbar.Pz(); qbar_E = qbar.E();
q_Phi = phi;
q_CosTheta = costheta;
// Testez si l evenement genere est acceptable
// et remplissez ce qu il faut
if ( z1 < dsigma(theta, qMass, Qq) ) {
//if ( z1 < dsigma(ele, pos, q, qbar, qMass, Qq)){
h_theta->Fill(theta);
h_cosTheta->Fill(cos(theta));
h_phi->Fill(phi);
h_pT->Fill(sqrt(pow(q_Px,2)+pow(q_Py,2))); //Quantité de mouvement transverse du quark
tree->Fill();
nAccepted++;
}
nTotal++;
}// fin de la boucle while
// Imprimez l efficacite de la methode de rejection
cout << "Efficacité : " << nAccepted << "/" << nTotal << " = " << nAccepted * 1. / nTotal << endl;
//Section efficace intégrée
double sigmaTot = 4./3. * pi * alphaEM * alphaEM * 3. * Qq * Qq / 4 / eBeam / eBeam * (1. + 2. * qMass * qMass / 4 / eBeam / eBeam);
cout << "Section efficace intégrée : " << sigmaTot << " GeV^{-2} ou " << sigmaTot*3.8707e8 << " pb. Comme on a généré " << nAccepted << " événements, cela donne une luminosité générée de " << nAccepted / sigmaTot / 3.8707e8 << " pb^{-1}. C'est " << nAccepted / sigmaTot / 3.8707e8 / 162.1 << " fois la luminosité totale enregistrée par l'expérience DELPHI entre 1989 et 1995." << endl;
// Declarez une fonction et fittez la distribution
// en cos(theta)
TF1 *f1 = new TF1("f1", myFunc, -1., 1., 2);
f1->SetParameters(100, 1.);
// Affichez les histogrammes dans un canvas
TCanvas *can = new TCanvas();
can->Divide(2,2);
can->cd(1);
h_theta->SetMinimum(0.);
h_theta->DrawCopy();
can->cd(2);
h_cosTheta->SetMinimum(0.);
h_cosTheta->Fit("f1");
h_cosTheta->DrawCopy();
can->cd(3);
h_phi->SetMinimum(0.);
h_phi->DrawCopy();
can->cd(4);
h_pT->SetMinimum(0.);
h_pT->DrawCopy();
// Sauvez le TTree ainsi que les histogrammes
tree->Write();
h_theta->Write();
h_cosTheta->Write();
h_phi->Write();
h_pT->Write();
f.Close();
}
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <iostream>
#include <TRandom3.h>
#include <time.h>
#include <cmath>
// On definit des constantes valables
// dans tout le code.
#define pi 3.141592
#define Ebeam 10.
#define alphaEM 0.0072973525
#define alphaS 0.1184
#define Nc 3
using namespace std;
// Codez ici une formule permettant de faire le passage
// du referentiel alpha, beta, gamma au referentiel du
// laboratoire.
// Cette fonction doit prendre les 3 angles en argument
// plus une reference au quadrivecteur a modifier.
void myRotate(double alpha, double beta, double gamma, TLorentzVector &particle)
{
// Pour acceder au trivecteur du quadrivecteur
// on utilise la methode Vect()
TVector3 myVect = particle.Vect();
// La methode Rotate(alpha, V) effectue
// une rotation d un angle alpha autour de
// l axe defini par le vecteur V
particle.SetVect(myVect);
}
// Libre a vous de tester la fonction
// de rotation que vous venez de coder
void testMyRotate()
{
TLorentzVector ele(1.,0.,0.,0.); //(px, py, pz, E) dans le ref. alpha, beta, gamma
myRotate(pi/4, pi/4, 0, ele);
cout << ele.Px() << " " << ele.Py() << " " << ele.Pz() << endl;
TLorentzVector pos(0.,1.,0.,0.);
myRotate(pi/4, 0, 0, pos);
cout << pos.Px() << " " << pos.Py() << " " << pos.Pz() << endl;
TLorentzVector posbis(1.,1.,1.,0.);
myRotate(pi/10, pi/6, pi/4, posbis);
cout << posbis.Px() << " " << posbis.Py() << " " << posbis.Pz() << endl;
}
// Codez ici l expression de la section efficace
// differentielle en termes des quadrivecteurs.
double dsigma(TLorentzVector ele, TLorentzVector pos, TLorentzVector quark, TLorentzVector antiq, TLorentzVector gluon)
{
double result = 1;
//result = factor * bracket / deno;
return result;
}
// Codez ici l expression de la section efficace
// differentielle en termes des variables x1, x2, beta, gamma
double dsigma(double x1, double x2, double beta, double gamma)
{
double result = 1;
//result = factor*bracket/deno;
return result;
}
void testDsigma()
{
cout << dsigma(0.8, 0.8, pi/4, pi/4) << endl;
}
// Codez ici la fonction principale qui genere les
// evenements et qui applique la methode de rejection.
void rejMethode()
{
// Determinez le sigmaMax en faisant une premiere
// boucle et en sauvant la valeur maximale de la
// section efficace
double sigmaMax = 1.E-15;
// Declarez toutes les variables necessaires pour
// votre programme
int npoints = 500000, nAccepted = 0;
double alpha, cosBeta, beta, gamma, x1, x2, x3, cosTheta1, cosTheta2, fmax;
double x1_min = 0.00005;
// Declarez les 5 quadri vecteurs ainsi que les
// variables a sauver
TLorentzVector ele, pos, quark, antiq, gluon;
double x1_Rec, x2_Rec, x3_Rec;
// Declarez une serie d histogrammes vous permettant
// de controler vos resultats
TH1F *h_cosThetaQuark = new TH1F("h_cosThetaQuark","h_cosThetaQuark",100, -1., 1.);
TH1F *h_cosThetaAntiq = new TH1F("h_cosThetaAntiq","h_cosThetaAntiq",100, -1., 1.);
TH1F *h_alpha = new TH1F("h_alpha","h_alpha",100, 0., 6.29);
TH1F *h_gamma = new TH1F("h_gamma","h_gamma",100, 0., 6.29);
TH1F *h_cosBeta = new TH1F("h_cosBeta","h_cosBeta",100, -1., 1.);
TH1F *h_beta = new TH1F("h_beta","h_beta",100, 0., 3.142);
TH3F *h_xxx = new TH3F("h_xxx", "h_xxx", 100, 0, 1, 100, 0, 1, 100, 0, 1);
TH1F *h_x1 = new TH1F("h_x1","h_x1",100,0,1);
TH1F *h_x2 = new TH1F("h_x2","h_x2",100,0,1);
TH1F *h_x3 = new TH1F("h_x3","h_x3",100,0,1);
TH2F *h_x1x2 = new TH2F("h_x1x2","h_x1x2;x_{1};x_{2}",100,0,1,100,0,1);
TH1F *h_cosThetaEK = new TH1F("h_cosThetaEK","h_cosThetaEK;cos#theta_{EK}",100,0,1);
TH2F *h_sigmax1 = new TH2F("h_sigmax1","h_sigmax1",1000, 0, 0.000001, 100, 0, 1);
// Initialisez le generateur de nombres aleatoires
TRandom3 generateur;
//generateur.SetSeed(time(NULL));
// Ici commence la boucle sur les evenements a generer
for (int i=0; i<npoints; i++)
{
// Generez les differents nombres aleatoires necessaires
// sur les intervalles adequats
// Construisez toute la cinematique dans le referentiel alpha, beta, gamma
// Passage au referentiel du labo
myRotate(alpha, beta, gamma, quark);
myRotate(alpha, beta, gamma, antiq);
myRotate(alpha, beta, gamma, gluon);
h_alpha->Fill(alpha); h_gamma->Fill(gamma); h_beta->Fill(beta); h_cosBeta->Fill(cosBeta);
// Inserez la condition de la methode de rejection
double Sigma = 0;
if (Sigma > sigmaMax) {
sigmaMax = Sigma;
cout << " !!!!!! " << endl;
}
if ( Sigma < fmax ) continue;
nAccepted++;
// Remplissez vos Histogrammes
h_cosThetaQuark->Fill(cos(quark.Theta())); h_cosThetaAntiq->Fill(cos(antiq.Theta()));
}//fin de la boucle de generation d evenements
// Hors de la boucle, dessinez vos histogrammes dans
// des canevas. Faites des fit sur certaines distributions,
// sauvez les histogrammes.
TF1 *f = new TF1("f","[0]*(1.+[1]*x*x)", -0.9, 0.9);
f->SetParameters(100, 1);
TCanvas *can = new TCanvas(); can->Divide(2,3);
can->cd(1); h_cosThetaQuark->Fit("f"); h_cosThetaQuark->DrawCopy();
can->cd(2); h_cosThetaAntiq->Fit("f"); h_cosThetaAntiq->DrawCopy();
can->cd(3); h_alpha->DrawCopy();
can->cd(4); h_beta->DrawCopy();
can->cd(5); h_gamma->DrawCopy();
can->cd(6); h_cosBeta->DrawCopy();
TCanvas *can2 = new TCanvas(); can2->Divide(2,2);
can2->cd(1); h_xxx->DrawCopy();
can2->cd(2); h_x1->DrawCopy();
can2->cd(3); h_x2->DrawCopy();
can2->cd(4); h_sigmax1->DrawCopy();
can2->cd(5); h_x1x2->DrawCopy();
can2->cd(6); h_cosThetaEK->DrawCopy();
cout << "Sigma Max = " << sigmaMax << endl;
cout << "Eff = " << nAccepted*1./npoints << endl;
}
Remarque: le code ci-dessous est inutilement complexe pour la partie algorithme de Jade + algorithme anti-kT. En effet, dans votre cas vous pouvez y arriver avec un code beaucoup plus simple puisque vous n'avez que 3 particules. J'ai tenu à laisser la version complète pour vous montrer un exemple plus "complexe" de code C++, mais il n'était évidemment pas attendu que vous trouviez cela pendant le TP.
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <iostream>
#include <TRandom3.h>
#include <time.h>
#include <cmath>
#include <vector>
#include <algorithm>
// On definit des constantes valables
// dans tout le code.
#define pi 3.141592
#define Ebeam 10.
#define alphaEM 0.0072973525
#define alphaS 0.1184
#define Nc 3
using namespace std;
// Codez ici une formule permettant de faire le passage
// du referentiel alpha, beta, gamma au referentiel du
// laboratoire.
// Cette fonction doit prendre les 3 angles en argument
// plus une reference au quadrivecteur a modifier.
void myRotate(double alpha, double beta, double gamma, TLorentzVector &particle)
{
// Pour acceder au trivecteur du quadrivecteur
// on utilise la methode Vect()
TVector3 myVect = particle.Vect();
TVector3 X(1,0,0);
TVector3 Z(0,0,1);
// La methode Rotate(alpha, V) effectue
// une rotation d un angle alpha autour de
// l axe defini par le vecteur V
myVect.Rotate(-gamma, Z); X.Rotate(-gamma, Z);
myVect.Rotate(-beta, X); Z.Rotate(-beta, X);
myVect.Rotate(-alpha, Z); X.Rotate(-alpha, Z);
particle.SetVect(myVect);
}
// Libre a vous de tester la fonction
// de rotation que vous venez de coder
void testMyRotate()
{
TLorentzVector ele(1.,0.,0.,0.); //(px, py, pz, E) dans le ref. alpha, beta, gamma
myRotate(pi/4, pi/4, 0, ele);
cout << ele.Px() << " " << ele.Py() << " " << ele.Pz() << endl;
TLorentzVector pos(0.,1.,0.,0.);
myRotate(pi/4, 0, 0, pos);
cout << pos.Px() << " " << pos.Py() << " " << pos.Pz() << endl;
TLorentzVector posbis(1.,1.,1.,0.);
myRotate(pi/10, pi/6, pi/4, posbis);
cout << posbis.Px() << " " << posbis.Py() << " " << posbis.Pz() << endl;
}
// Codez ici l expression de la section efficace
// differentielle en termes des quadrivecteurs.
double dsigma(TLorentzVector ele, TLorentzVector pos, TLorentzVector quark, TLorentzVector antiq, TLorentzVector gluon)
{
double result;
double s = 4*Ebeam*Ebeam;
double Qq = 2./3;
double factor = 8*4./3*(4*pi)*(4*pi)*(4*pi)*alphaEM*alphaEM*alphaS*Nc*Qq*Qq*s;
double bracket = (ele*quark) * (ele*quark) + (ele*antiq)*(ele*antiq) + (pos*quark)*(pos*quark) + (pos*antiq)*(pos*antiq);
double deno = (ele*pos) * (quark*gluon) * (antiq*gluon) * 2 * s * 32 * pow(2*pi, 5);
result = factor * bracket / deno;
return result;
}
// Codez ici l expression de la section efficace
// differentielle en termes des variables x1, x2, beta, gamma
double dsigma(double x1, double x2, double beta, double gamma)
{
double result;
double s = 4*Ebeam*Ebeam;
double Qq = 2./3;
double x3 = 2 - x1 - x2;
double theta1 = acos( -(x1*x1 + x3*x3 - x2*x2) /2./x1/x3 );
double theta2 = acos( -(x2*x2 + x3*x3 - x1*x1) /2./x2/x3 );
double factor = alphaEM*alphaEM*alphaS*4./3*Nc*Qq*Qq/(4*pi*pi*s);
double bracket = x1*x1 + x2*x2;
bracket += x1*x1*sin(beta)*sin(beta)*sin(gamma+theta1)*sin(gamma+theta1);
bracket += x2*x2*sin(beta)*sin(beta)*sin(gamma-theta2)*sin(gamma-theta2);
double deno = (1-x1)*(1-x2);
result = factor*bracket/deno;
return result;
}
void testDsigma()
{
cout << "dsigma(0.8, 0.8, pi/4, pi/4) = " << dsigma(0.8, 0.8, pi/4, pi/4) << endl;
cout << dsigma(0.5, 0.8, 0.5, 1.2) << endl;
}
// Codez ici la fonction principale qui genere les
// evenements et qui applique la methode de rejection.
void rejMethode(int nMax = 500000)
{
// Determinez le sigmaMax en faisant une premiere
// boucle et en sauvant la valeur maximale de la
// section efficace
//double sigmaMax = 1.0e-15;
double sigmaMax = 2.83888e-09;
// Declarez toutes les variables necessaires pour
// votre programme
int npoints = 0, nAccepted = 0;
double alpha, cosBeta, beta, gamma, x1, x2, x3, cosTheta1, cosTheta2, fmax;
double x1_min = 0.00005, y1, y2;
// Declarez les 5 quadri vecteurs ainsi que les
// variables a sauver
TLorentzVector ele, pos, quark, antiq, gluon;
double x1_Rec, x2_Rec, x3_Rec;
// Declarez une serie d histogrammes vous permettant
// de controler vos resultats
TH1F *h_cosThetaQuark = new TH1F("h_cosThetaQuark","h_cosThetaQuark",100, -1., 1.);
TH1F *h_cosThetaAntiq = new TH1F("h_cosThetaAntiq","h_cosThetaAntiq",100, -1., 1.);
TH1F *h_alpha = new TH1F("h_alpha","h_alpha",100, 0., 6.29);
TH1F *h_gamma = new TH1F("h_gamma","h_gamma",100, 0., 6.29);
TH1F *h_cosBeta = new TH1F("h_cosBeta","h_cosBeta",100, -1., 1.);
TH1F *h_beta = new TH1F("h_beta","h_beta",100, 0., 3.142);
TH3F *h_xxx = new TH3F("h_xxx", "h_xxx;x_{1};x_{2};x_{3}", 100, 0, 1, 100, 0, 1, 100, 0, 1);
TH1F *h_x1 = new TH1F("h_x1","h_x1",100,0,1);
TH1F *h_x2 = new TH1F("h_x2","h_x2",100,0,1);
TH1F *h_x3 = new TH1F("h_x3","h_x3",100,0,1);
TH2F *h_x1x2 = new TH2F("h_x1x2","h_x1x2;x_{1};x_{2}",100,0,1,100,0,1);
TH1F *h_cosThetaEK = new TH1F("h_cosThetaEK","h_cosThetaEK;cos#theta_{EK}",100,0,1);
TH2F *h_sigmax1 = new TH2F("h_sigmax1","h_sigmax1",1000, 0, 0.000000003, 100, 0, 1); //mettre sigmaMax comme borne sup
TH1F *h_nbJetsJade = new TH1F("h_nbJetsJade","h_nbJetsJade;Nombre de jets (selon l'algorithme de Jade)",5,0,5);
TH1F *h_nbJetsAntiKT = new TH1F("h_nbJetsAntiKT","h_nbJetsAntiKT;Nombre de jets (selon l'algorithme d'anti-kT)",5,0,5);
// Initialisez le generateur de nombres aleatoires
TRandom3 generateur;
//generateur.SetSeed(time(NULL));
generateur.SetSeed(123456);
int previous_nAccepted = 0;
// Ici commence la boucle sur les evenements a generer
while (nAccepted < nMax)
{
if (nAccepted % 100000 == 0 && previous_nAccepted!=nAccepted) cout << nAccepted << " événements." << endl;
npoints++;
previous_nAccepted = nAccepted;
// Generez les differents nombres aleatoires necessaires
// sur les intervalles adequats
y1 = -log(1-x1_min) + generateur.Rndm()*(-log(1-(1.-x1_min)) + log(1. - x1_min));
y2 = -log(1-x1_min) + generateur.Rndm()*(-log(1-(1.-x1_min)) + log(1. - x1_min));
x1 = 1. - exp(-y1);
x2 = 1. - exp(-y2);
//x1 = x1_min + generateur.Rndm()*(1.- x1_min - x1_min);
//x2 = x1_min + generateur.Rndm()*(1.- x1_min - x1_min);
x3 = 2. - x1 - x2;
if (x1+x2 < 1. + x1_min) continue;
alpha = 2.*pi*generateur.Rndm();
cosBeta = -1. + 2*generateur.Rndm();
beta = acos(cosBeta);
gamma = 2.*pi*generateur.Rndm();
fmax = generateur.Rndm()*sigmaMax;
// Construisez toute la cinematique dans le referentiel alpha, beta, gamma
cosTheta1 = - ( x1*x1 + x3*x3 - x2*x2 )/( 2*x1*x3 );
cosTheta2 = - ( x2*x2 + x3*x3 - x1*x1 )/( 2*x2*x3 );
ele.SetPxPyPzE(0.,0., Ebeam, Ebeam); pos.SetPxPyPzE(0.,0., -Ebeam, Ebeam);
quark.SetPxPyPzE(x1*Ebeam*cosTheta1, x1*Ebeam*sqrt(1-cosTheta1*cosTheta1), 0., x1*Ebeam);
antiq.SetPxPyPzE(x2*Ebeam*cosTheta2, -x2*Ebeam*sqrt(1-cosTheta2*cosTheta2), 0., x2*Ebeam);
gluon.SetPxPyPzE(x3*Ebeam, 0., 0., x3*Ebeam);
// Passage au referentiel du labo
myRotate(alpha, beta, gamma, quark);
myRotate(alpha, beta, gamma, antiq);
myRotate(alpha, beta, gamma, gluon);
// Inserez la condition de la methode de rejection
//double Sigma = dsigma(ele, pos, quark, antiq, gluon)*exp(-y1)*exp(-y2);
double Sigma = dsigma(x1, x2, beta, gamma)*exp(-y1)*exp(-y2);
//double Sigma = dsigma(x1, x2, beta, gamma);
if (Sigma > sigmaMax) {
sigmaMax = Sigma;
cout << " !!!!!! " << endl;
cout << sigmaMax << " pour x1 = " << x1 << " et x2 = " << x2 << endl;
}
if ( Sigma < fmax ) continue;
nAccepted++;
h_alpha->Fill(alpha); h_gamma->Fill(gamma); h_beta->Fill(beta); h_cosBeta->Fill(cosBeta);
// Remplissez vos Histogrammes, Tree, etc ...
h_cosThetaQuark->Fill(cos(quark.Theta())); h_cosThetaAntiq->Fill(cos(antiq.Theta()));
x1_Rec = max( max(x1, x2), max(x1, x3) );
x3_Rec = min( min(x1, x2), min(x1, x3) );
x2_Rec = 2 - x1_Rec - x3_Rec; // De telle manière que x1_Rec > x2_Rec > x3_Rec.
h_xxx->Fill(x1_Rec, x2_Rec, x3_Rec);
h_x1->Fill(x1_Rec); h_x2->Fill(x2_Rec); h_x3->Fill(x3_Rec); h_x1x2->Fill(x1_Rec,x2_Rec);
h_cosThetaEK->Fill((x2_Rec-x3_Rec)/x1_Rec);
h_sigmax1->Fill(Sigma, x1_Rec);
// Algorithme de Jade
double ycut = 0.1;
vector<TLorentzVector> finalStateParticles;
vector<TLorentzVector> finalJets;
finalStateParticles.push_back(quark);
finalStateParticles.push_back(antiq);
finalStateParticles.push_back(gluon);
vector<int> skipEntries;
int initialSize = finalStateParticles.size();
for(int h = 0; h < initialSize ; h++){
finalJets.clear();
for(unsigned int i = 0; i < finalStateParticles.size() ; i++){
if(find(skipEntries.begin(),skipEntries.end(),i) != skipEntries.end()) continue;
bool isAJet = true;
for(unsigned int j = i+1 ; j < finalStateParticles.size() ; j++){
if((finalStateParticles[i]+finalStateParticles[j]).M() < sqrt(ycut) * Ebeam * 2.) {
finalJets.push_back(finalStateParticles[i]+finalStateParticles[j]); //Si le "jet" est trop proche d'un autre, on les combine en une seule particule.
isAJet = false;
skipEntries.push_back(j);
}
}
if(isAJet) finalJets.push_back(finalStateParticles[i]);
}
finalStateParticles = finalJets;
}
h_nbJetsJade->Fill(finalJets.size());
// Algorithme anti-kT (pour info)
double R = 0.4;
finalJets.clear();
finalStateParticles.clear();
finalStateParticles.push_back(quark);
finalStateParticles.push_back(antiq);
finalStateParticles.push_back(gluon);
skipEntries.clear();
while(skipEntries.size() < finalStateParticles.size()){
double best_diB = 1e200;
double best_dij = 1e200;
int best_diB_i = -1;
int best_dij_i = -1;
int best_dij_j = -1;
for(unsigned int i = 0; i < finalStateParticles.size() ; i++){
if(find(skipEntries.begin(),skipEntries.end(),i) != skipEntries.end()) continue;
double current_diB = 1./pow(finalStateParticles[i].Pt(),2);
if((current_diB < best_diB) || ((current_diB > 1e200) && (skipEntries.size() == finalStateParticles.size()-1)) ){
best_diB = current_diB;
best_diB_i = i;
}
for(unsigned int j = i+1 ; j < finalStateParticles.size() ; j++){
if(find(skipEntries.begin(),skipEntries.end(),j) != skipEntries.end()) continue;
double current_dij = std::min(1/pow(finalStateParticles[i].Pt(),2),1/pow(finalStateParticles[j].Pt(),2))* pow (finalStateParticles[i].Eta() - finalStateParticles[j].Eta(),2) * pow(finalStateParticles[i].Phi() - finalStateParticles[j].Phi(),2)/pow(R,2.);
if(current_dij < best_dij){
best_dij = current_dij;
best_dij_i = i;
best_dij_j = j;
}
}
}
if((best_diB < best_dij) || ((best_diB > 1e200) && (skipEntries.size() == finalStateParticles.size()-1)) ){
finalJets.push_back(finalStateParticles[best_diB_i]); //Si la particule est plus proche du faisceau que d'un autre jet, on la garde
skipEntries.push_back(best_diB_i);
}
else{
finalStateParticles[best_dij_i] = finalStateParticles[best_dij_i] + finalStateParticles[best_dij_j]; //Sinon, on combine avec le jet le plus proche et on recommence
skipEntries.push_back(best_dij_j);
}
}
h_nbJetsAntiKT->Fill(finalJets.size());
}//fin de la boucle de generation d evenements
// Hors de la boucle, dessinez vos histogrammes dans
// des canevas. Faites des fit sur certaines distributions,
// sauvez les histogrammes.
TF1 *f = new TF1("f","[0]*(1.+[1]*x*x)", -0.9, 0.9);
f->SetParameters(100, 1);
TCanvas *can = new TCanvas(); can->Divide(2,3);
can->cd(1); h_cosThetaQuark->Fit("f"); h_cosThetaQuark->DrawCopy();
can->cd(2); h_cosThetaAntiq->Fit("f"); h_cosThetaAntiq->DrawCopy();
can->cd(3); h_alpha->DrawCopy();
can->cd(4); h_beta->DrawCopy();
can->cd(5); h_gamma->DrawCopy();
can->cd(6); h_cosBeta->DrawCopy();
TCanvas *can2 = new TCanvas(); can2->Divide(2,3);
can2->cd(1); h_xxx->DrawCopy();
can2->cd(2); h_x1->DrawCopy();
can2->cd(3); h_x2->DrawCopy();
can2->cd(4); h_x3->DrawCopy();
can2->cd(5); h_x1x2->DrawCopy();
can2->cd(6); h_cosThetaEK->DrawCopy();
TCanvas *can3 = new TCanvas(); can3->Divide(2,1);
can3->cd(1); h_nbJetsJade->DrawCopy();
can3->cd(2); h_nbJetsAntiKT->DrawCopy();
cout << "Sigma Max = " << sigmaMax << endl;
cout << "Eff = " << nAccepted*1./npoints << endl;
can->SaveAs("Result_angles.pdf");
can2->SaveAs("Result_x.pdf");
can3->SaveAs("Result_nbJets.pdf");
can->SaveAs("Result_angles.root");
can2->SaveAs("Result_x.root");
can3->SaveAs("Result_nbJets.root");
}