makePyEvts.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

// nEvents is how many events we want. 
int makePyEvts(int nEvents) 
{ 

  double P4TOT[4];
  // Create an instance of the Pythia event generator ... 
 
  TPythia6* pythia = new TPythia6; 

  TH1F *h1 = new TH1F("h1","Higgs Mass",50,100.0,200.0);
  
  // 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,175); //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
  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(225,1,1);
  
  //Set collision
  pythia->Initialize("cms", "p", "p", ECM);
  #define FILENAME   "pythia.root"
  #define TREENAME   "tree"
  #define BRANCHNAME "particles"

 // Open an output file 
  TFile* file = TFile::Open(FILENAME, "RECREATE"); 
  if (!file || !file->IsOpen()) {
    Error("makeEventSample", "Couldn;t open file %s", FILENAME);
    return 1;
  }  


  // Make a tree in that file ... 
  TTree* tree = new TTree(TREENAME, "Pythia 6 tree"); 
   
  // Register a the cache of pythia on a branch (It's a 
  // TClonesArray of TMCParticle objects. )
  TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
  tree->Branch(BRANCHNAME, &particles); 
  int KF_muon =13;
  // Now generate 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();
    //
        
    for (Int_t i=0; i<4; i++){P4TOT[i]=0.0;}
    //P4TOT[0]=0.0;
    //P4TOT[1]=0.0;
    //P4TOT[2]=0.0;
    //P4TOT[3]=0.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++){
      Int_t KS = pythia->GetK(itk,1); //KS code
      Int_t KF = pythia->GetK(itk,2); //KF code
      if ((abs(KF)== KF_muon)&&(KS!=21)){
        N_muons++;
        P4TOT[0]=P4TOT[0]+pythia->GetP(itk,1); // get track 4-vector
        P4TOT[1]=P4TOT[1]+pythia->GetP(itk,2);  
        P4TOT[2]=P4TOT[2]+pythia->GetP(itk,3);
        P4TOT[3]=P4TOT[3]+pythia->GetP(itk,4);
      }//if             
    }//itk
    //cout << "MyOutput: Number of muons =" << N_muons<<endl; 
    Float_t Hmass2 = pow(P4TOT[3],2)-pow(P4TOT[0],2)-pow(P4TOT[1],2)-pow(P4TOT[2],2); 
    Float_t Hmass=0.0;
    if (Hmass2 >= 0.0) {Hmass = sqrt(Hmass2);}
    else {Hmass = -sqrt(abs(Hmass2));}
    //    cout<<"Mass(#mu^+#mu^-#mu^+#mu^-) "<<Hmass<<endl;
    h1->Fill(Hmass);
 }//ievt



   gROOT->Reset();
   c1 = new TCanvas("c1","pp -> H-> ZZ",200,10,700,500);
   c1->SetFillColor(42);
   c1->SetGrid();

  title = new TPaveLabel(0.1,0.94,0.9,0.98,
          "Pythia MC: pp -> H -> ZZ");
  title->SetFillColor(16);
  title->SetTextFont(52);
  title->Draw();

  h1->GetXaxis()->SetTitle("CM energy (GeV)");
  h1->GetXaxis()->CenterTitle();
  h1->GetYaxis()->SetTitle("Number Events/(2 GeV)");
  h1->GetYaxis()->CenterTitle();
  h1->GetYaxis()->SetRangeUser(0,1400);
  h1->GetXaxis()->SetRangeUser(140,200);
  h1->SetMarkerStyle(21);
  h1->SetMarkerColor(31);
  h1->SetMarkerSize(1);
  h1->SetLineColor(2);
  h1->SetLineWidth(4);
  h1->SetTitle("pp #rightarrow H #rightarrow ZZ");
  h1->Draw("e1p");

   c1->Update();
   c1->GetFrame()->SetFillColor(21);
   c1->GetFrame()->SetBorderSize(3);
   c1->Modified();  


   //  h1->SetMarkerStyle(21);
   //h1->SetMarkerSize(0.7);
   //h1->Draw("e1p");
   pythia->Pystat(2);
   pythia->Pystat(4);
   pythia->Pylist(11);
   pythia->Pylist(12);
   pythia->Draw();

}//makePyEvts


//____________________________________________________________________ 
//  
// EOF
//

---