makePythiagamma.C

Go to the documentation of this file.

// Program to use Pythia to generate H->gamm,gamma from within ROOT.
// To make an event sample (of size 100) do 
//
//    shell> root
//    root [0] .L makePyEvtsgammagamma.C
//    root [1] makePyEvts(100,800)
//
//
#define ECM 14000
// using namespace TMath;

#include <iostream>
using namespace std;

#include "TStopwatch.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TPythia6.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TRandom2.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TVector2.h"
#include "TFile.h"
#include "TROOT.h"

//double higgswidth;

// nEvents is how many events we want. 
int makePyEvts(int nEvents,float Higgsmass = 300 ) 
{ 


  // Create an instance of the Pythia event generator ... 
  TPythia6* pythia = new TPythia6; 

 
  // Initialise it to run p+p at ECM GeV in CMS 

  pythia->SetMSEL(0);        //Full user control 
  pythia->SetMSUB(102,1);    //Select gg->Higgs
  pythia->SetPMAS(25,1,Higgsmass); //Set Higgs mass
  pythia->SetMSTP(61,0);     //Turn off initial-state radiation
  //  pythia->SetMSTP(71,0);     //Turn off final-state QCD and QED radiation
  //  pythia->SetMSTP(81,0);
  //  pythia->SetMSTP(111,0);
 
 // Turn off all Higgs (h0) decay modes; then turn on correct mode
 
  int KF = 25;               // 
  int idcmin = pythia->GetMDCY(KF,2);
  int idcmax = idcmin + pythia->GetMDCY(KF,3);
  cout<<"MyOutput: Higgs idcmin="<<idcmin<<" idcmax="<<idcmax<<endl;

  for (int idc = idcmin; idc < idcmax; idc++)
    {
      pythia->SetMDME(idc,1,0);
    }   
  pythia->SetMDME(223,1,1);
  
  //Set collision
  pythia->Initialize("cms", "p", "p", ECM);

 

  float higgswidth = pythia->GetPMAS(25,2);
  float higgs_generated_mass = 0;
  cout<<"The higgs width from Pythia is "<<higgswidth<<endl;


  TStopwatch watcher;
  watcher.Reset();
  watcher.Start();


  float Higgsmass_max = 3000;
  float Higgsmaxplot = 3000;
  float Higgsminplot = 0.0;

  int Nbins = 500;
      if (Higgsmass<191)
    {
      Nbins = 500;
      Higgsmaxplot = Higgsmass + 5 * higgswidth;
      Higgsminplot = Higgsmass - 5 * higgswidth;
      }


  TH1F *Hmass = new TH1F("Hmass","Higgs Mass",Nbins,Higgsminplot,Higgsmaxplot);
  TH1F *Hpt = new TH1F("Hpt","Higgs Tranverse Momentum",500,0.,500.0);
  TH1F *Hgeneratedmass = new TH1F("Hgeneratedmass","Higgs Generated Mass",Nbins,Higgsminplot,Higgsmaxplot);

  // Now generate events

  int Num = 0;
  //loop over the events
  for (Int_t ievt = 0; ievt < nEvents; ievt++) 
    {
      // Show how far we got every 100'th event. 
      if (ievt % 100 == 0) 
        cout << "Event # " << ievt << endl;


      
      // Make one event. 
      pythia->GenerateEvent();
   
      
      //displays the particle listing for the first three events
      if(Num<1)
        {pythia->Pylist(1);
      }

      Num++; 
        
      //Get some event information!
      //Total number of Pythia tracks
      Int_t N = pythia->GetN();
      //cout << "MyOutput: Number of Pythia tracks ="<<N<<endl;
      //Loop over all tracks

      TLorentzVector gammasum;

      int N_gammas = 0;
  
      for (Int_t itk =1; itk<=N; itk++) //loop over the particles in one event
        {

          Int_t KS = pythia->GetK(itk,1); //KS code
          
          if(KS!=21) continue;

          Int_t KF = pythia->GetK(itk,2); //KF code
          
          if ( KF == 25) //higgs

            {
              TLorentzVector higgs = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4)); 
              Hgeneratedmass->Fill(higgs.M());  
            } 
   
          if( KF == 22) //photon
            {
              
              TLorentzVector lvGamma = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4));
      

              

              if(KS==21)
                {
                  N_gammas++;
                  gammasum+=lvGamma;
                  
                  
                  
                    
                }
            }

        }
        //itk //end of the loop over particles in one event
      //cout << "MyOutput: Number of muons =" << N_muons<<endl; 
      // cout<<"Number of photons that decayed from Higgs: "<<N_gammas<<endl;       


 
      
      Hmass->Fill(gammasum.M());
      Hpt->Fill(gammasum.Pt());
    }//ievt /end of loop over the events
  
 

 

  char name[100];
 
   // sprintf(name,"IVplots_Hm%5.2f_Nevt%d.eps",Higgsmass,nEvents);
   // c1->Print(name);
  
 

  // sprintf(name,"AnglePlots_Hm%5.2f_Nevt%d.eps",Higgsmass,nEvents);
  // c3->Print(name);


   TCanvas * c5 = new TCanvas("c5","**",800,800);
   c5->SetFillColor(9);
   Hmass->Draw("e");
   TF1 *f_fit = new TF1("f_fit","[0]*TMath::BreitWigner(x,[1],[2])",0,500);
   TF1 *f_theory = new TF1("f_theory","[0]*TMath::BreitWigner(x,[1],[2])",0,500);
   f_fit->SetParNames("Constant","Mass","Width");
   f_theory->SetParNames("Constant","Mass","Width");

   if (Higgsmass>180)
     {
       f_fit->SetParameters(Hmass->GetEntries(),Higgsmass,higgswidth);
       f_fit->SetParLimits(2, higgswidth - higgswidth/2, higgswidth + higgswidth/2); 
       f_fit->SetParLimits(0,0,100*Hmass->GetEntries());
       f_fit->SetParLimits(1,Higgsmass - Higgsmass/2,Higgsmass + Higgsmass/2);
     }
   else
     {
       f_fit->SetParameters(Hmass->GetEntries(),Higgsmass,higgswidth);
       f_fit->SetParLimits(2, higgswidth - higgswidth/5, higgswidth + higgswidth/5); 
       f_fit->SetParLimits(0,0,100*Hmass->GetEntries());
       f_fit->SetParLimits(1,Higgsmass - 2*higgswidth,Higgsmass + 2*higgswidth);
     }
   for(int number = 0;number<3;number++)
     { 
       Hmass->Fit("f_fit","L","",Higgsminplot,Higgsmaxplot);
       
     }
   float fit_higgswidth = f_fit->GetParameter(2);
   float fit_higgsmass = f_fit->GetParameter(1);
   
   f_theory->SetParameter(0,Hmass->GetEntries());
   f_theory->FixParameter(1,Higgsmass);
   f_theory->FixParameter(2,higgswidth);
   f_theory->SetParLimits(0,0,10000);
   
   Hmass->Fit("f_theory","L+","",Higgsminplot,Higgsmaxplot);
   
   cout<<"Fitted Higgs mass: "<<fit_higgsmass<<", Fitted Higgs width: "<<fit_higgswidth<<endl;
   cout<<"Theory Higgs mass: "<<Higgsmass<<", Theory Higgs width: "<<higgswidth<<endl;
   
   sprintf(name,"histograms_Htogammagamma_Hm%5.2f_Nevt%d.root",Higgsmass,nEvents);
   TFile f(name,"recreate");
   f.cd();
   int N=gROOT->GetList()->GetSize();
   for (int key=0; key<N ; key++)
     {
       gROOT->GetList()->At(key)->Write();
       //      gROOT->GetList()->At(key)->Delete();
     }
   f.Close();
   
   watcher.Stop();
   cout <<"Time spent for "<<nEvents<<" events :"<< watcher.RealTime()<<" [s]"<<endl;
   
   return -97534;
   */
     }//makePyEvts

---