.rst .pdf

Contents

Soluzioni: Lezione 8

8.1 generazione di toy experiment

• esercizi01_02.cpp

/*
Esercizio 01: Si scriva un programma che, fissato NMAX numero di eventi di un campione, generi NToys toy experiment di quel campione e ne calcoli la media.
Esercizio 02: Si Aggiunge al programma precedente un oggetto di tipo TH1F che visualizzi la distribuzione delle medie al variare dei toy experiment.

c++ -o esercizi01_02 esercizi01_02.cpp ../../Lezione_04/Esercizi/statistiche.cc `root-config --glibs --cflags` 
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"

#include "../../Lezione_04/Esercizi/statistiche.h"  

float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  }
  
//----------------- MAIN -----------------  
  
int main (int argc, char ** argv)
  {
    srand(time(NULL));
    
    // numero di eventi generati per un singolo toy experiment
    int NMAX = 500 ;
    // numero di toy experiment per construire la pdf della media per campionamento
    int NToys = 10000 ;
            
    //definisco due oggetti della classe statistiche
    statistiche* s_singleToy = new statistiche(NMAX);
    statistiche* s_media = new statistiche(NToys);
    
    TH1F h_medie ("h_medie", "Distribuzione delle medie", 41, -0.5, 0.5) ;
    
    // loop su tutti i toy experiment
    for (int iToy = 0 ; iToy < NToys ; ++iToy)
      {
        int i = 0 ;
        //loop sul singolo esperimento
        while (i++ < NMAX) s_singleToy->aggiungiNumero (rand_range (-3., 3.)) ;
                       
        h_medie.Fill(s_singleToy->media());
        s_media->aggiungiNumero (s_singleToy->media ()) ;        
        s_singleToy->reset () ;
      } 

     TCanvas c1 ("c1", "c1", 100, 100, 1000, 1000);
     h_medie.GetXaxis()->SetTitle("media del campione");
     h_medie.GetYaxis()->SetTitle("numero di eventi");
     h_medie.Draw();
     
     c1.Print("esercizio2.png", "png");

    delete s_singleToy;
    delete s_media;


    return 0 ;
  }

8.2 disegno di toy experiment

• esercizi01_02.cpp

/*
Esercizio 01: Si scriva un programma che, fissato NMAX numero di eventi di un campione, generi NToys toy experiment di quel campione e ne calcoli la media.
Esercizio 02: Si Aggiunge al programma precedente un oggetto di tipo TH1F che visualizzi la distribuzione delle medie al variare dei toy experiment.

c++ -o esercizi01_02 esercizi01_02.cpp ../../Lezione_04/Esercizi/statistiche.cc `root-config --glibs --cflags` 
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"

#include "../../Lezione_04/Esercizi/statistiche.h"  

float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  }
  
//----------------- MAIN -----------------  
  
int main (int argc, char ** argv)
  {
    srand(time(NULL));
    
    // numero di eventi generati per un singolo toy experiment
    int NMAX = 500 ;
    // numero di toy experiment per construire la pdf della media per campionamento
    int NToys = 10000 ;
            
    //definisco due oggetti della classe statistiche
    statistiche* s_singleToy = new statistiche(NMAX);
    statistiche* s_media = new statistiche(NToys);
    
    TH1F h_medie ("h_medie", "Distribuzione delle medie", 41, -0.5, 0.5) ;
    
    // loop su tutti i toy experiment
    for (int iToy = 0 ; iToy < NToys ; ++iToy)
      {
        int i = 0 ;
        //loop sul singolo esperimento
        while (i++ < NMAX) s_singleToy->aggiungiNumero (rand_range (-3., 3.)) ;
                       
        h_medie.Fill(s_singleToy->media());
        s_media->aggiungiNumero (s_singleToy->media ()) ;        
        s_singleToy->reset () ;
      } 

     TCanvas c1 ("c1", "c1", 100, 100, 1000, 1000);
     h_medie.GetXaxis()->SetTitle("media del campione");
     h_medie.GetYaxis()->SetTitle("numero di eventi");
     h_medie.Draw();
     
     c1.Print("esercizio2.png", "png");

    delete s_singleToy;
    delete s_media;


    return 0 ;
  }

8.3 calcolo delle statistiche di un campione di toy experiment

• esercizio03.cpp

/*
Esercizio 03: Si utilizzi la classe statistiche sviluppata durante la Lezione 4 per confrontare la deviazione standard della media calcolata per ogni singolo toy con la deviazione standard del campione delle medie.

c++ -o esercizio03 esercizio03.cpp ../../Lezione_04/Esercizi/statistiche.cc  
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "../../Lezione_04/Esercizi/statistiche.h" 

float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  }
  
//----------------- MAIN -----------------  
  
int main (int argc, char ** argv)
  {
    srand(time(NULL));   
  
    // numero di eventi generati per un singolo toy experiment
    int NMAX = 500 ;
    // numero di toy experiment per construire la pdf della media per campionamento
    int NToys = 10000 ;
            
    statistiche* s_singleToy = new statistiche(NMAX);
    statistiche* s_media = new statistiche(NToys);
    statistiche* s_incertezzaMedia = new statistiche(NToys);
           
    // loop su tutti i toy experiment
    for (int iToy = 0 ; iToy < NToys ; ++iToy)
      {
        int i = 0 ;
        //loop sul singolo esperimento
        while (i++ < NMAX) s_singleToy->aggiungiNumero (rand_range (-3., 3.)) ;

        s_media->aggiungiNumero (s_singleToy->media ()) ;                      
        s_incertezzaMedia->aggiungiNumero(s_singleToy->dev_standard_media ());
                       
        s_singleToy->reset () ;
      } 


    //valori finali
    std::cout << "media delle deviazioni standard della media per i singoli toy experiment: " 
              << s_incertezzaMedia->media () << " \n" ;
    std::cout << "deviazione standard della distribuzione delle medie dei singoli toy experiment: " 
              << s_media->dev_standard_SQ () << " \n" ;


    delete s_singleToy;
    delete s_media;
    delete s_incertezzaMedia;


    return 0 ;
  }

8.4 studio dell’evoluzione delle statistiche con il numero di eventi

• esercizio04.cpp

/*
Esercizio 04: Si utilizzino due TGraph per confrontare l'evoluzione della deviazione standard della media calcolata per ogni singolo toy con la deviazione standard del campione delle medie al variare del numero di eventi generati in un singolo toy experiment.

c++ -o esercizio04 esercizio04.cpp ../../Lezione_04/Esercizi/statistiche.cc `root-config --glibs --cflags` 
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"
#include "TGraph.h"

#include "../../Lezione_04/Esercizi/statistiche.h"

float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  }
  
//----------------- MAIN -----------------  
  
int main (int argc, char ** argv)
  {
    srand(time(NULL)); 
  
    // numero massimo di eventi generati per un singolo toy experiment
    int NMAX = 40000 ;

    statistiche* s_singleToy = new statistiche(NMAX);
    statistiche* s_media = new statistiche(NMAX);
    statistiche* s_incertezzaMedia = new statistiche(NMAX);
           
    TGraph g_sigma ;
    TGraph g_sigmaMedia ;
          
    //loop sul numero di eventi da generare per singolo esperimento  
    for (int n = 10 ; n < NMAX ; n = n*2)
      {
        int i = 0 ;
                       
        //loop sul singolo esperimento
        while (i++ < n) s_singleToy->aggiungiNumero (rand_range (-3., 3.)) ;

        s_media->aggiungiNumero (s_singleToy->media ()) ;                
        g_sigma.SetPoint (g_sigma.GetN(), n, s_media->dev_standard ( )) ;        
                
        s_incertezzaMedia->aggiungiNumero(s_singleToy->dev_standard_media ( ));        
        g_sigmaMedia.SetPoint (g_sigmaMedia.GetN(), n, s_incertezzaMedia->media ()) ;        
                
        s_singleToy->reset () ;        
      }   
      
      

    TCanvas c1 ("c1", "c1", 100, 100, 1000, 1000) ;
    c1.SetLogx () ; 

    g_sigma.SetMarkerStyle (20) ;
    g_sigma.SetMarkerColor (kBlue + 1) ;
    g_sigma.SetLineColor (kBlue - 9) ;
    g_sigma.SetMarkerSize (2) ;
    g_sigma.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi") ;
    g_sigma.GetHistogram ()->GetYaxis ()->SetTitle ("media della deviazione standard") ;
    g_sigma.Draw ("ALP") ;

    c1.Print ("esercizio4_sigma.png", "png") ;

    TCanvas c2 ("c2", "c2", 100, 100, 1000, 1000) ;
    c2.SetLogx () ;

    g_sigmaMedia.SetMarkerStyle (20) ;
    g_sigmaMedia.SetMarkerColor (kGreen + 1) ;
    g_sigmaMedia.SetLineColor (kGreen -5) ;
    g_sigmaMedia.SetMarkerSize (2) ;
    g_sigmaMedia.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi") ;
    g_sigmaMedia.GetHistogram ()->GetYaxis ()->SetTitle ("deviazione standard dalla media") ;
    g_sigmaMedia.Draw ("ALP") ;    

    c2.Print ("esercizio4_sigmaMedia.png", "png") ;

    delete s_singleToy;
    delete s_media;
    delete s_incertezzaMedia;


    return 0 ;
  }

8.5 metodo di integrazione hit-or-miss

• esercizio05.cpp

/*
Esercizio 05: Si implementi il metodo di integrazione hit-or-miss con la funzione di esempio f(x) = sin (x).
Si scriva l'algoritmo che calcola l'integrale come una funzione esterna al programma main, facendo in modo che prenda come parametri di ingresso, oltre agli estremi lungo l'asse x e l'asse y, anche il numero di punti pseudo-casuali da generare.
Si faccia in modo che l'algoritmo ritorni un contenitore contenente due elementi: il primo elemento sia il valore dell'integrale, il secondo sia la sua incertezza.

c++ -o esercizio05 esercizio05.cpp ../../Lezione_03/Esercizi/esercizio04/myarray.cc
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fsin (double x) 
  {
    return sin (x) ; 
  }
  
mioArray HitOrMiss (int Nrand, double g(double), double xMin, double xMax,
                    double yMin, double yMax)
   {
     int nHit = 0;
     double area  = (xMax - xMin) * (yMax - yMin) ;
     
     for (int i = 0 ; i < Nrand ; ++i) 
      {
        double x = rand_range (xMin, xMax) ;
        double y = rand_range (yMin, yMax) ; 
        if (y < g (x)) ++nHit ; 
      }   
      
     double p = nHit / static_cast<double> (Nrand) ;
     double varianza = area * area / static_cast<double> (Nrand) * p * (1. - p) ;
           
     mioArray risultato (2) ;      
     risultato.fill (0, nHit * area / static_cast<double> (Nrand)) ;
     risultato.fill (1, sqrt (varianza)) ;
      
    return risultato ;       
   }           

//----------------- MAIN ----------------- 
int main (int argc, char ** argv)
{ 
  srand (time (NULL)) ;
  
  int N = 10000 ;
  double x_min = 0. ;
  double x_max = M_PI ; 
  double y_min = 0. ; 
  double y_max = 1. ;
   
  mioArray risultato = HitOrMiss(N,fsin,x_min,x_max,y_min,y_max);
  
  std::cout << "Integrale = " << risultato.get1 (0) << " +/- " << risultato.get1 (1) << std::endl ;
            
  return 0 ;
}

8.6 incetezza del metodo di integrazione hit-or-miss

• esercizio06.cpp

/*
Esercizio 06: Si inserisca il calcolo dell'integrale dell'esercizio precedente in un ciclo che, al variare del numero N di punti generati, mostri il valore dell'integrale e della sua incertezza.
Si utilizzi un TGraph per disegnare gli andamenti del valore dell'integrale e della sua incertezza, al variare di N con ragione logaritmica.

c++ -o esercizio06 esercizio06.cpp ../../Lezione_03/Esercizi/esercizio04/myarray.cc `root-config --glibs --cflags`
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"
#include "TGraph.h"

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fsin (double x) 
  {
    return sin (x) ; 
  }
  
mioArray HitOrMiss (int Nrand, double g(double), double xMin, double xMax,
              double yMin, double yMax )
   {
     int nHit = 0;
     double area  = (xMax - xMin) * (yMax - yMin) ;
     
     for (int i = 0 ; i < Nrand ; ++i) 
      {
        double x = rand_range (xMin, xMax) ;
        double y = rand_range (yMin, yMax) ; 
        if (y < g (x)) ++nHit ; 
      }   
      
     double p = nHit / static_cast<double> (Nrand) ;
     double varianza = area * area / static_cast<double> (Nrand) * p * (1. - p) ;
     
     mioArray risultato (2) ;      
     risultato.fill (0, nHit * area / static_cast<double> (Nrand)) ;
     risultato.fill (1, sqrt (varianza)) ;
      
    return risultato ; 
                 
   }           

//----------------- MAIN ----------------- 
int main (int argc, char ** argv)
{
 
  srand (time (NULL)) ;
  int N = 10000 ;
  double x_min = 0. ;
  double x_max = M_PI ; 
  double y_min = 0. ; 
  double y_max = 1. ;
  
  TGraph g_integrale ;
  TGraph g_incertezza;
  
  for (int iN = 10; iN < N; iN = iN*2)
  {
   mioArray risultato = HitOrMiss(iN,fsin,x_min,x_max,y_min,y_max);
  
  std::cout << "Integrale = " << risultato.get1 (0) << " +/- " << risultato.get1 (1) << std::endl ;
    
    g_integrale.SetPoint (g_integrale.GetN(), iN, risultato.get1 (0)) ;  
    g_incertezza.SetPoint (g_incertezza.GetN(), iN, risultato.get1 (1)) ;                              
  }
    
   TCanvas c1 ("c1", "c1", 100, 100, 1000, 1000) ;
   c1.SetLogx () ; 

   g_integrale.SetMarkerStyle (20) ;
   g_integrale.SetMarkerColor (kRed + 1) ;
   g_integrale.SetLineColor (kRed - 9) ;
   g_integrale.SetMarkerSize (2) ;
   g_integrale.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi generati") ;
   g_integrale.GetHistogram ()->GetYaxis ()->SetTitle ("integrale") ;
   g_integrale.Draw ("ALP") ;

   c1.Print ("esercizio6_integrale.png", "png") ; 
   
   TCanvas c2 ("c2", "c2", 100, 100, 1000, 1000) ;
   c2.SetLogx () ; 

   g_incertezza.SetMarkerStyle (20) ;
   g_incertezza.SetMarkerColor (kBlue + 1) ;
   g_incertezza.SetLineColor (kBlue - 9) ;
   g_incertezza.SetMarkerSize (2) ;
   g_incertezza.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi generati") ;
   g_incertezza.GetHistogram ()->GetYaxis ()->SetTitle ("incertezza") ;
   g_incertezza.Draw ("ALP") ;

   c2.Print ("esercizio6_incertezza.png", "png") ;    
                                  
  return 0 ;
}

8.7 metodo di integrazone crude Montecarlo

• esercizio07.cpp

/*
Esercizio 07:Si implementi il metodo di integrazione crude-MC con la funzione di esempio f(x) = sin (x).
Si scriva l'algoritmo che calcola l'integrale come una funzione esterna al programma main, facendo in modo che prenda come parametri di ingresso, oltre agli estremi lungo l'asse x, anche il numero di punti pseudo-casuali da generare.
Si faccia in modo che l'algoritmo ritorni un contenitore contenente due elementi: il primo elemento sia il valore dell'integrale, il secondo sia la sua incertezza.

c++ -o esercizio07 ../../Lezione_03/Esercizi/esercizio04/myarray.cc esercizio07.cpp 
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fsin (double x) 
  {
    return sin (x) ; 
  }
  
 
mioArray CrudeMC (int Nrand, double g(double), double xMin, double xMax )
   {
     double lunghezza  = (xMax - xMin);
     double somma = 0;
     double sommaQ = 0;
     double media = 0;
     double varianza = 0;
     
     for (int i = 0; i<Nrand; i++)
     {
        double x = rand_range(xMin, xMax);
        somma += g(x);
        sommaQ += g(x)*g(x);     
     }
      
     media = somma/static_cast<double>(Nrand);
     varianza = sommaQ/static_cast<double>(Nrand) - media*media ;
          
     mioArray risultato (2) ;      
     risultato.fill (0, media*lunghezza) ;
     risultato.fill (1, sqrt (varianza/static_cast<double>(Nrand)) * lunghezza) ;
      
    return risultato ; 
                                         
   }           

//----------------- MAIN -----------------
int main (int argc, char ** argv)
{ 
  srand (time (NULL)) ;
  int N = 10000 ;
  double x_min = 0. ;
  double x_max = M_PI ; 
  
  mioArray risultato = CrudeMC(N,fsin,x_min,x_max);  
  std::cout << "Integrale = " << risultato.get1 (0) << " +/- " << risultato.get1 (1) << std::endl ;
                                    
  return 0 ;
}

8.8 incetezza del metodo di integrazione crude Montecarlo

• esercizio08.cpp

/*
Esercizio 08:Si inserisca il calcolo dell'integrale dell'esercizio precedente in un ciclo che, al variare del numero N di punti generati, mostri il valore dell'integrale e della sua incertezza.
Si utilizzi un TGraph per disegnare gli andamenti del valore dell'integrale e della sua incertezza, al variare di N con ragione logaritmica.
Si sovrapponga questo TGraph a quello ottenuto dallo svolgimento dell'Esercizio 6.6.

c++ -o esercizio08 esercizio08.cpp ../../Lezione_03/Esercizi/esercizio04/myarray.cc `root-config --glibs --cflags`
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"
#include "TGraph.h"

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fsin (double x) 
  {
    return sin (x) ; 
  }
  
 
mioArray CrudeMC (int Nrand, double g(double), double xMin, double xMax )
   {
     double lunghezza  = (xMax - xMin);
     double somma = 0;
     double sommaQ = 0;
     double media = 0;
     double varianza = 0;
     
     for (int i = 0; i<Nrand; i++)
     {
        double x = rand_range(xMin, xMax);
        somma += g(x);
        sommaQ += g(x)*g(x);     
     }
      
     media = somma/static_cast<double>(Nrand);
     varianza = sommaQ/static_cast<double>(Nrand) - media*media ;
                                   
     mioArray risultato_MC (2) ;      
     risultato_MC.fill (0, media*lunghezza) ;
     risultato_MC.fill (1, sqrt (varianza/static_cast<double>(Nrand)) * lunghezza) ;
      
    return risultato_MC ;       
   }  
   
mioArray HitOrMiss (int Nrand, double g(double), double xMin, double xMax,
              double yMin, double yMax)
   {
     int nHit = 0;
     double area  = (xMax - xMin) * (yMax - yMin) ;
     
     for (int i = 0 ; i < Nrand ; ++i) 
      {
        double x = rand_range (xMin, xMax) ;
        double y = rand_range (yMin, yMax) ; 
        if (y < g (x)) ++nHit ; 
      }   
      
     double p = nHit / static_cast<double> (Nrand) ;
     double varianza = area * area / static_cast<double> (Nrand) * p * (1. - p) ;
           
     mioArray risultato_HoM (2) ;      
     risultato_HoM.fill (0, nHit * area / static_cast<double> (Nrand)) ;
     risultato_HoM.fill (1, sqrt (varianza)) ;
      
    return risultato_HoM ;       
   }           
            
//----------------- MAIN -----------------
int main (int argc, char ** argv)
{ 
  srand (time (NULL)) ;
  int N = 50000 ;
  double x_min = 0. ;
  double x_max = M_PI ; 
  double y_min = 0. ; 
  double y_max = 1. ;

  TGraph g_integrale_MC ;
  TGraph g_incertezza_MC;
  TGraph g_integrale_HoM ;
  TGraph g_incertezza_HoM;
  
 for (int iN = 10; iN < N; iN = iN*2)
  {  
    mioArray risultato_MC  = CrudeMC(iN,fsin,x_min,x_max); 
    mioArray risultato_HoM = HitOrMiss(iN,fsin,x_min,x_max,y_min,y_max);    

    g_integrale_MC.SetPoint (g_integrale_MC.GetN(), iN, risultato_MC.get1 (0)) ;  
    g_incertezza_MC.SetPoint (g_incertezza_MC.GetN(), iN, risultato_MC.get1 (1)) ; 
    
    g_integrale_HoM.SetPoint (g_integrale_HoM.GetN(), iN, risultato_HoM.get1 (0)) ;  
    g_incertezza_HoM.SetPoint (g_incertezza_HoM.GetN(), iN, risultato_HoM.get1 (1)) ;                             
  }
  
  TCanvas c1 ("c1", "c1", 100, 100, 1000, 1000) ;
   c1.SetLogx () ; 

   g_integrale_HoM.SetMarkerStyle (21) ;
   g_integrale_HoM.SetMarkerColor (kBlue + 1) ;
   g_integrale_HoM.SetLineColor (kBlue - 9) ;
   g_integrale_HoM.SetMarkerSize (2) ;
   g_integrale_HoM.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi generati") ;
   g_integrale_HoM.GetHistogram ()->GetYaxis ()->SetTitle ("integrale") ;
   g_integrale_HoM.GetHistogram ()->GetYaxis ()->SetRangeUser (1.1, 2.6) ; //definisco il range dell'asse y
   g_integrale_HoM.Draw ("ALP") ;

   g_integrale_MC.SetMarkerStyle (20) ;
   g_integrale_MC.SetMarkerColor (kRed + 1) ;
   g_integrale_MC.SetLineColor (kRed - 9) ;
   g_integrale_MC.SetMarkerSize (2) ;
   g_integrale_MC.Draw ("LPSAME") ;
   
   c1.Print ("esercizio8_integrale.png", "png") ; 
   
   TCanvas c2 ("c2", "c2", 100, 100, 1000, 1000) ;
   c2.SetLogx () ; 

   g_incertezza_HoM.SetMarkerStyle (21) ;
   g_incertezza_HoM.SetMarkerColor (kBlue + 1) ;
   g_incertezza_HoM.SetLineColor (kBlue - 9) ;
   g_incertezza_HoM.SetMarkerSize (2) ;
   g_incertezza_HoM.GetHistogram ()->GetXaxis ()->SetTitle ("numero di eventi generati") ;
   g_incertezza_HoM.GetHistogram ()->GetYaxis ()->SetTitle ("integrale") ;
   g_incertezza_HoM.Draw ("ALP") ;

   g_incertezza_MC.SetMarkerStyle (20) ;
   g_incertezza_MC.SetMarkerColor (kRed + 1) ;
   g_incertezza_MC.SetLineColor (kRed - 9) ;
   g_incertezza_MC.SetMarkerSize (2) ;
   g_incertezza_MC.Draw ("LPSAME") ;
   
   c2.Print ("esercizio8_incertezza.png", "png") ;    
                   
  return 0 ;
}

8.9 confronto fra crude Montecarlo e hit-or-miss

• esercizio09.cpp

/*
Esercizio 09:Si disegnino in due TGraph gli andamenti della precisione del calcolo dell'integrale, per i due algoritmi di hit-or-miss e crude-MC, in funzione del tempo di calcolo corrispondente alle varie scelte del numero totale N di eventi pseudo-casuali generati.

c++ -o esercizio09 esercizio09.cpp ../../Lezione_03/Esercizi/esercizio04/myarray.cc `root-config --glibs --cflags`
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "TH1F.h"
#include "TCanvas.h"
#include "TGraph.h"

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fsin (double x) 
  {
    return sin (x) ; 
  }
  
 
mioArray CrudeMC (int Nrand, double g(double), double xMin, double xMax )
   {
     double lunghezza  = (xMax - xMin);
     double somma = 0;
     double sommaQ = 0;
     double media = 0;
     double varianza = 0;
     
     for (int i = 0; i<Nrand; i++)
     {
        double x = rand_range(xMin, xMax);
        somma += g(x);
        sommaQ += g(x)*g(x);     
     }
      
     media = somma/static_cast<double>(Nrand);
     varianza = sommaQ/static_cast<double>(Nrand) - media*media ;
                                   
     mioArray risultato_MC (2) ;      
     risultato_MC.fill (0, media*lunghezza) ;
     risultato_MC.fill (1, sqrt (varianza/static_cast<double>(Nrand)) * lunghezza) ;
      
    return risultato_MC ;       
   }  
   
mioArray HitOrMiss (int Nrand, double g(double), double xMin, double xMax,
              double yMin, double yMax)
   {
     int nHit = 0;
     double area  = (xMax - xMin) * (yMax - yMin) ;
     
     for (int i = 0 ; i < Nrand ; ++i) 
      {
        double x = rand_range (xMin, xMax) ;
        double y = rand_range (yMin, yMax) ; 
        if (y < g (x)) ++nHit ; 
      }   
      
     double p = nHit / static_cast<double> (Nrand) ;
     double varianza = area * area / static_cast<double> (Nrand) * p * (1. - p) ;
           
     mioArray risultato_HoM (2) ;      
     risultato_HoM.fill (0, nHit * area / static_cast<double> (Nrand)) ;
     risultato_HoM.fill (1, sqrt (varianza)) ;
      
    return risultato_HoM ;       
   }           
            
//----------------- MAIN -----------------
int main (int argc, char ** argv)
{ 
  srand (time (NULL)) ;
  int N = 50000 ;
  double x_min = 0. ;
  double x_max = M_PI ; 
  double y_min = 0. ; 
  double y_max = 1. ;
  double tempoCalcolo = 0.;

  TGraph g_incertezza_MC;
  TGraph g_incertezza_HoM;
  
 for (int iN = 10; iN < N; iN = iN*2)
  {  
    mioArray risultato_MC  = CrudeMC(iN,fsin,x_min,x_max); 
    mioArray risultato_HoM = HitOrMiss(iN,fsin,x_min,x_max,y_min,y_max); 
    
    tempoCalcolo = clock();   
 
    g_incertezza_MC.SetPoint (g_incertezza_MC.GetN(), iN, risultato_MC.get1 (1)) ; 
    g_incertezza_HoM.SetPoint (g_incertezza_HoM.GetN(), iN, risultato_HoM.get1 (1)) ;                             
  }
  
  TCanvas *c1 = new TCanvas ("c1", "c1", 100, 100, 1000, 1000) ;
  c1->Divide(1,2);
  c1->cd(1);
  g_incertezza_HoM.SetMarkerStyle (20) ;
  g_incertezza_HoM.SetMarkerColor (kBlue + 1) ;
  g_incertezza_HoM.SetLineColor (kBlue - 9) ;
  g_incertezza_HoM.SetMarkerSize (2) ;
  g_incertezza_HoM.SetTitle("Metodo Hit or Miss");
  g_incertezza_HoM.GetHistogram ()->GetXaxis ()->SetTitle ("tempo di calcolo") ;
  g_incertezza_HoM.GetHistogram ()->GetYaxis ()->SetTitle ("incertezza") ;
  g_incertezza_HoM.Draw ("ALP") ;
   
  c1->cd(2);
  g_incertezza_MC.SetMarkerStyle (20) ;
  g_incertezza_MC.SetMarkerColor (kRed + 1) ;
  g_incertezza_MC.SetLineColor (kRed - 9) ;
  g_incertezza_MC.SetMarkerSize (2) ;
  g_incertezza_MC.SetTitle("Metodo Crude Monte Carlo");
  g_incertezza_MC.GetHistogram ()->GetXaxis ()->SetTitle ("tempo di calcolo") ;
  g_incertezza_MC.GetHistogram ()->GetYaxis ()->SetTitle ("incertezza") ;
  g_incertezza_MC.Draw ("ALP") ;
   
  c1->Print ("esercizio9.png", "png") ;   
  delete c1;
                     
  return 0 ;
}

8.10 integrale della Gaussiana

• esercizio10.cpp

/*
Esercizio 10: Si utilizzi il metodo hit-or-miss per stimare l’integrale sotteso ad una disrtibuzione di probabilita' Gaussiana con μ=0 e *σ=1 in un generico intervallo [a,b].
Si calcoli l'integrale contenuto entro gli intervalli [-kσ,kσ] al variare di k da 1 a 5.

c++ -o esercizio10 esercizio10.cpp ../../Lezione_03/Esercizi/esercizio04/myarray.cc
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

#include "../../Lezione_03/Esercizi/esercizio04/myarray.h"


float rand_range (float min, float max)
  {
    return min + (max - min) * rand () / static_cast<float> (RAND_MAX) ;
  } 


double fgaus (double x)  
  {
    return exp (-0.5 * x * x) ; 
  }
  
mioArray HitOrMiss (int Nrand, double g(double), double xMin, double xMax,
                    double yMin, double yMax)
   {
     int nHit = 0;
     double area  = (xMax - xMin) * (yMax - yMin) ;
     
     for (int i = 0 ; i < Nrand ; ++i) 
      {
        double x = rand_range (xMin, xMax) ;
        double y = rand_range (yMin, yMax) ; 
        if (y < g (x)) ++nHit ; 
      }   
      
     double p = nHit / static_cast<double> (Nrand) ;
     double varianza = area * area / static_cast<double> (Nrand) * p * (1. - p) ;
           
     mioArray risultato (2) ;      
     risultato.fill (0, nHit * area / static_cast<double> (Nrand)) ;
     risultato.fill (1, sqrt (varianza)) ;
      
    return risultato ;       
   }           

//----------------- MAIN -----------------
int main (int argc, char ** argv)
{ 
  srand (time (NULL)) ;
  int N = 10000 ;
  double x_min = 0.;
  double x_max = 0. ; 
  double y_min = 0. ; 
  double y_max = 1. ;

  int k = 1;
  
  while (k <= 5)
   {
    x_min = - k;
    x_max = k;
    mioArray risultato = HitOrMiss(N,fgaus,x_min,x_max,y_min,y_max);
    std::cout << "Intervallo [" << x_min << " , " << x_max << "]" << std::endl;
    std::cout << "Integrale = " << risultato.get1 (0) << " +/- " << risultato.get1 (1) << std::endl ; 
    k++;                             
  }
                   
  return 0 ;
}

previous

Esercizi per la lezione 8

next

Lezione 9: la verosimiglianza