makePyEvts_angular_wphi_axial.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
  pythia->SetMSTP(25,2);

  // 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());
     

      //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);

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


      
      double projection_z_1, projection_z_2;
      TVector3 muon1_hrest_vector3, muon1_x1;
      TVector3 muon2_hrest_vector3, muon2_T;      

      muon1_hrest_vector3 = muon1_hrest.Vect();
      projection_z_1 = -uZ1.Dot(muon1_hrest_vector3);
      muon1_x1 = muon1_hrest_vector3 + projection_z_1*uZ1;
      double check_1 = muon1_x1.Dot(uZ1);
      //cout<<"The dot product between the x-axis and z-axis for muon 1 is:"<<check_1<<endl;
      
      muon2_hrest_vector3 = muon2_hrest.Vect();
      projection_z_2 = uZ1.Dot(muon2_hrest_vector3);
      muon2_T = muon2_hrest_vector3 - projection_z_2*uZ1;
      double check_2 = muon2_T.Dot(uZ2);
      //cout<<"The dot product between the x-axis and z-axis for muon 2 is:"<<check_2<<endl;      

      double phi, phi_0_2pi;

      TVector3 y_axis_1, muon1_x1_unit;
      muon1_x1_unit = muon1_x1.Unit();
      y_axis_1 = muon1_x1.Dot(uZ1);
      

      double dot_y_mu2;
      double sign;
      dot_y_mu2 = y_axis_1.Dot(muon2_T);
      if(dot_y_mu2 > 0)
        {sign = +1;
        }
      else
        {sign = -1;}
        
      phi = sign*(muon2_T.Angle(muon1_x1));
      if(phi<0)
        {phi_0_2pi = 2*TMath::Pi()+phi;
        }
      else
        {phi_0_2pi = phi;
        }
      Phiplot1->Fill(phi_0_2pi);
          
      

      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("c7","**",800,800);
  c7->cd();
  Z2plot->Draw();

  TCanvas * c8 = new TCanvas("c8","**",800,800);
  c8->cd();
  Phiplot1->Draw();

  TF1 *f7 = new TF1("f7","[0]*(1+[1]*[2]*cos(x)+[2]*cos(2*x))",0,2*TMath::Pi());
  
  f7->SetParNames("Normalization","a1","a2");
 
  //f7->SetParameter(1,0.9);
  
  
  
  f7->FixParameter(0,(Phiplot1->Integral("bin width"))/(2*TMath::Pi()));
  
  Phiplot1->Fit("f7","L");
  
  


  sprintf(name,"histograms_correct_HtoZZ_angular_wphi_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

---