makePyEvts_angular.C

Go to the documentation of this file.

// Program to use Pythia to generate H->ZZ; Z->mu+mu- from within ROOT.
// To make an event sample (of size 100) do 
//
//    shell> root
//    root [0] .L makePyEvts.C
//    root [1] makePyEvts(100)
//
//
#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
 

  // Turn off all Z decay modes; then turn on correct mode
  int KF = 23;               // 
  int idcmin = pythia->GetMDCY(KF,2);
  int idcmax = idcmin + pythia->GetMDCY(KF,3);
  cout<<"MyOutput: Z idcmin="<<idcmin<<" idcmax="<<idcmax<<endl;
  for (int idc = idcmin; idc < idcmax; idc++)
    {
      pythia->SetMDME(idc,1,0);
    }   
  pythia->SetMDME(184,1,1);
 
  // Turn off all Higgs (h0) decay modes; then turn on correct mode
  KF = 25;               // 
  idcmin = pythia->GetMDCY(KF,2);
  idcmax = idcmin + pythia->GetMDCY(KF,3);
  cout<<"MyOutput: Higgs idcmin="<<idcmin<<" idcmax="<<idcmax<<endl;
  for (Int_t idc = idcmin; idc < idcmax; idc++)
    {
      pythia->SetMDME(idc,1,0);
    }   
  pythia->SetMDME(225,1,1);
  
  //Set collision
  pythia->Initialize("cms", "p", "p", ECM);

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


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


  float Higgsmass_max = 3000;

  // if (Higgsmass > Higgsmass_max)
  //  {cout<<"I am not happy !"<<endl; return;}
    
  int NbinIVplot = 100;


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



  TH1F *Hmass = new TH1F("Hmass","Higgs Mass",Nbins,Higgsminplot,Higgsmaxplot);

  TH1F *Z1plot = new TH1F("Z1mass","Z1mass",Nbins,0,150);
  TH1F *Z2plot = new TH1F("Z2mass","Z2mass",Nbins,0,150);

  TH1F *Cosine1plot = new TH1F("Cosine1graph","Cos(#theta) angle of muon1 in Z1 rest frame",NbinIVplot,-1,1);
  TH1F *Cosine2plot = new TH1F("Cosine2graph","Cos(#theta) angle of muon2 in Z2 rest frame",NbinIVplot,-1,1);
  TH2F *Correlationplot = new TH2F("Correlationplot","Theta angle muon 1 versus theta angle muon 2",NbinIVplot,-1,1,NbinIVplot,-1,1);

  TH1F *Phiplot1 = new TH1F("Phiplot1","Phi angle of muon 1 in Z1 rest frame",NbinIVplot,0,2*TMath::Pi());
  TH1F *Phiplot2 = new TH1F("Phiplot2","Phi angle of muon 2 in Z2 rest frame",NbinIVplot,0,2*TMath::Pi());

 
  Int_t KF_muon =13;

  // Now generate events

  TRandom2 R;

  //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();
      // Now we're ready to fill the tree, and the event is over. 
      //      tree->Fill();
      //
      
      TLorentzVector muonsum;
      TLorentzVector muon1plus,  muon2plus, muon1minus , muon2minus;

      int Nmuonplus=0;
      int Nmuonminus=0;
      
      int Nz=0;
      
      //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


      
      Int_t N_muons =0;     // initialize number muons
      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 ((abs(KF)== KF_muon))
            {
              TLorentzVector lvMuon = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4));
              
             

              if(KS==21)
                {
                  N_muons++;
                  // hmuonpt_initial->Fill(lvMuon.Pt());
                  muonsum+=lvMuon;
                  
                  
                  if (KF<0) //anti-muon
                    {
                      if (Nmuonplus==1)
                        {
                          muon2plus = lvMuon;
                        }
                      if (Nmuonplus==0)
                        {
                          muon1plus = lvMuon;
                          Nmuonplus=1;
                        }
                    }
                  else //muon
                    {
                      if (Nmuonminus==1)
                        {
                          muon2minus = lvMuon;
                        }
                      if (Nmuonminus==0)
                        {
                          muon1minus = lvMuon;
                          Nmuonminus=1;
                        }
                    }
                }
              else
                {
                  // hmuonpt_final->Fill(lvMuon.Pt());
                }
            }
          
          
        }//itk //end of the loop over particles in one event
      //cout << "MyOutput: Number of muons =" << N_muons<<endl; 
      


      
      //Want to take the four vector of the muon and anti-muon and add them to get the four vectors of the Z
     
      TLorentzVector Z1,Z2;
      Z1 = muon1plus + muon1minus;
      Z2 = muon2plus + muon2minus;
      Z1plot->Fill(Z1.M());
      Z2plot->Fill(Z2.M());
     

      /*
      //get the opposite of the 3-vector of the Zs
      TVector3 boostZ1,boostZ2;
      boostZ1 = - Z1.BoostVector();
      boostZ2 = - Z2.BoostVector();

      //      boostZ1.Print();
      //      boostZ2.Print();

      //boost muon, anti-muon and Z1 to the Z1 rest frame
      TLorentzVector Z1_star,muon1_star,antimuon1_star;
      Z1_star = Z1;
      muon1_star = muon1minus;
      antimuon1_star = muon1plus;

      Z1_star.Boost(boostZ1);
      muon1_star.Boost(boostZ1);
      antimuon1_star.Boost(boostZ1);
      // Z1_star.Vect().Print();

      //boost muon, anti-muon and Z2 to the Z2 rest frame
      TLorentzVector Z2_star,muon2_star,antimuon2_star;
      Z2_star = Z2;
      muon2_star = muon2minus;
      antimuon2_star = muon2plus;

      Z2_star.Boost(boostZ2);
      muon2_star.Boost(boostZ2);
      antimuon2_star.Boost(boostZ2);
      // Z2_star.Vect().Print();
      */

      //boost Higgs, Z1, and Z2 into Higgs rest frame
      TVector3 boostH,boostZ1_hrest,boostZ2_hrest;
      TLorentzVector H,H_hrest,Z1_hrest,Z2_hrest,muon1_hrest,muon2_hrest,muon1_star,muon2_star,Z1_star,Z2_star;
      H = Z1+Z2;
      boostH = -H.BoostVector();
      H_hrest = H;
      Z1_hrest = Z1;
      Z2_hrest = Z2;
      muon1_hrest = muon1minus;
      muon2_hrest = muon2minus;
      H_hrest.Boost(boostH);
      Z1_hrest.Boost(boostH);
      Z2_hrest.Boost(boostH);
      muon1_hrest.Boost(boostH);
      muon2_hrest.Boost(boostH);
      
      boostZ1_hrest = -Z1_hrest.BoostVector();
      boostZ2_hrest = -Z2_hrest.BoostVector();
      
      Z1_star = Z1_hrest;
      Z2_star = Z2_hrest;
      
      muon1_star = muon1_hrest;
      muon2_star = muon2_hrest;
      
      muon1_star.Boost(boostZ1_hrest);
      muon2_star.Boost(boostZ2_hrest);

      TVector3 uZ1 = Z1_hrest.Vect().Unit();
      TVector3 uZ2 = Z2_hrest.Vect().Unit();

      float cosine1, cosine2;
      TVector3 unit_muon1_star = muon1_star.Vect().Unit();
      cosine1 = uZ1.Dot(unit_muon1_star);
      //      cout <<cosine1<<endl;
      //      cout <<" radian or degree "<<uZ1.Angle(unit_muon1_star)<<endl;
      //      cout <<cos( uZ1.Angle(unit_muon1_star) )<<endl;

      TVector3 unit_muon2_star = muon2_star.Vect().Unit();
      cosine2 = uZ2.Dot(unit_muon2_star);
      
      Cosine1plot->Fill(cosine1);
      Cosine2plot->Fill(cosine2);
      Correlationplot->Fill(cosine1,cosine2);

      TVector3 x1(1.0,0.0,0.0),y1(0.0,1.0,0.0);
      x1.RotateUz(uZ1);
      y1.RotateUz(uZ1);
      //(x1,y1,uZ1) is the rest frame of the Z, with z axis along the Z boson.
      
      TVector2 muon1_star2d(x1.Dot(unit_muon1_star),y1.Dot(unit_muon1_star));
      Phiplot1->Fill(muon1_star2d.Phi());
      
      TVector3 x2(1.0,0.0,0.0),y2(0.0,1.0,0.0);
      x2.RotateUz(uZ2);
      y2.RotateUz(uZ2);
      //(x2,y2,uZ2) is the rest frame of the Z, with z axis along the Z boson.
      
      TVector2 muon2_star2d(x2.Dot(unit_muon2_star),y2.Dot(unit_muon2_star));
      Phiplot2->Fill(muon2_star2d.Phi());
      

      bool didImessmuonplus=false;
      bool didImessmuonminus=false;

      //randomizing the anti-muons order
      if(R.Rndm()>0.5)
        {
          TLorentzVector temp;
          temp=muon1plus;
          muon1plus=muon2plus;
          muon2plus=temp;
          didImessmuonplus=true;
        }
      //randomizing the muons order
      if(R.Rndm()>0.5)
        {
          TLorentzVector temp;
          temp=muon1minus;
          muon1minus=muon2minus;
          muon2minus=temp;
          didImessmuonminus=true;
        }
         
      int MM=didImessmuonminus;
      int MP=didImessmuonplus;

      
      bool condition=(MM || MP) && (!MP || !MM);
      bool auxcond=((MP+MM)&1);
      if (auxcond!=condition)
        cout <<condition<<" "<<((MP+MM)&1)<<endl;

    

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

 
  
  
 
   char name[100];
  // sprintf(name,"IVplots_Hm%5.2f_Nevt%d.eps",Higgsmass,nEvents);
  // c1->Print(name);
  
  TCanvas * c3 = new TCanvas("c3","**",800,800);
  c3->Divide(2,2);
  c3->SetFillColor(9);
  c3->cd(1);
  Cosine1plot->Draw();
 
  Cosine1plot->SetTitle("Cos(#theta) of muon1 in Z1 rest frame");
//Cosine1graph->Rebin(8);
  TF1 *f2 = new TF1("f2","[0]*([1]*(3/4)*(1-x**2)+(3/8)*(1-[1])*(1+x**2))",-1,1);
  TF1 *f3 = new TF1("f3","[0]*( (1/(1+[1]))*(3/4)*(1-x**2) + (3/8)*([1]/(1+[1]))*(1+x**2) )",-1,1);
  //This fit, f4 is incorrect, look at ~/FromJeff/gamma_Z_macros/makePyEvts_angular_incorrect_boost_1.C for the correct fit  
TF1 *f4 = new TF1("f4","[0]*( ((1-[1])/([1]+3)) * (3/4)*(1-x**2) + (3/8)*((2*[1]+2)/([1]+3))*(1+x**2) )",-1,1);

  f2->SetParNames("Normalization","alpha");
  f3->SetParNames("Normalization","R");
  f4->SetParNames("Normalization","tullyR");

  f2->SetParameter(1,0.9);
  f3->SetParameter(1,0.1);
  f4->SetParameter(1,0.8);
  
  f2->FixParameter(0,Cosine1plot->Integral("bin width"));
  f3->FixParameter(0,Cosine1plot->Integral("bin width"));
  f4->FixParameter(0,Cosine1plot->Integral("bin width"));

  Cosine1plot->Fit("f2","L");
  Cosine1plot->Fit("f3","L+");
  Cosine1plot->Fit("f4","L+");




  c3->cd(2);
  Cosine2plot->Draw();
  c3->cd(3);
  Correlationplot->Draw();

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


  TCanvas *c4 = new TCanvas("c4","**",800,800);
  c4->Divide(2,2);
  c4->SetFillColor(9);
  c4->cd(1);
  Phiplot1->Draw();
  c4->cd(2);
  Phiplot2->Draw();

 
  TCanvas * c5 = new TCanvas("c5","**",800,800);
  c5->SetFillColor(9);
  Hmass->Draw("e");

  TCanvas * c6 = new TCanvas("c6","**",800,800);
  c6->cd();
  Z1plot->Draw();

  TCanvas * c7 = new TCanvas("c6","**",800,800);
  c7->cd();
  Z2plot->Draw();

  sprintf(name,"histograms_correct_HtoZZ_angular_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++)
    {

      TObject * obj = gROOT->GetList()->At(key);
      if (obj->InheritsFrom("TH1F") || obj->InheritsFrom("TH2F"))
        {
          cout <<"write "<<obj->GetName()<<endl;
          obj->Write();

        }
    }
  f.Close();

  
  for (int key=N-1; key>=0 ; key--)
    {
      
      TObject * obj = gROOT->GetList()->At(key);
      cout <<"delete "<<obj->GetName()<<endl;
      if (obj)  delete obj;
      cout <<"deleted "<<gROOT->GetList()->GetSize()<<" objects remaining"<<endl;
      }

  watcher.Stop();
  cout <<"Time spent for "<<nEvents<<" events :"<< watcher.RealTime()<<" [s]"<<endl;

  return -97534;

}//makePyEvts

---