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
// 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 = 1000 ) 
{ 

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

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


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


  float Higgsmass_max = 3000;

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


  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 *Hpt = new TH1F("Hpt","Higgs Tranverse Momentum",500,0.,500.0);
  /*
  TH1F *hmuonpt_initial = new TH1F("hmuonpt_initial","Muon Pt in initial state",200,0.,200.0);
  TH1F *hmuonpt_final = new TH1F("hmuonpt_final","Muon Pt in final state",200,0.,200.0);

  TH1F *hgamma = new TH1F("hgamma","Gamma Pt",200,0.,200.0);

  TH1F * mcZ1mass = new TH1F("mcZ1mass","generated Z mass",NbinIVplot,0,Higgsmass_max);
  TH1F * mcZ2mass = new TH1F("mcZ2mass","generated Z mass",NbinIVplot,0,Higgsmass_max);
  
  TH1F * mu1p_pt = new TH1F("mu1p_pt","#mu^{+}_{1} p_{T} [GeV]",NbinIVplot,0,Higgsmass_max);
  TH1F * mu2p_pt = new TH1F("mu2p_pt","#mu^{+}_{2} p_{T} [GeV]",NbinIVplot,0,Higgsmass_max);
  TH1F * mu1m_pt = new TH1F("mu1m_pt","#mu^{-}_{1} p_{T} [GeV]",NbinIVplot,0,Higgsmass_max);
  TH1F * mu2m_pt = new TH1F("mu2m_pt","#mu^{-}_{2} p_{T} [GeV]",NbinIVplot,0,Higgsmass_max);

  TH1F *Cosine1plot = new TH1F("Cosine1graph","Theta angle of muon1 in Z1 rest frame",NbinIVplot,-1,1);
  TH1F *Cosine2plot = new TH1F("Cosine2graph","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());

  TH1F *IV1 = new TH1F("IV1","Invariant mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV1wc = new TH1F("IV1wc","wrong invariant mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max); 
  TH1F *IV1rc = new TH1F("IV1rc","right invariant mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max); 
  IV1wc->SetLineColor(2);
  IV1rc->SetLineColor(4);

  TH1F *IV2 = new TH1F("IV2","Invariant mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV2wc = new TH1F("IV2wc","wrong Invariant mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV2rc = new TH1F("IV2rc","right Invariant mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  IV2wc->SetLineColor(2);
  IV2rc->SetLineColor(4);

  TH1F *IV3 = new TH1F("IV3","Invariant mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV3wc = new TH1F("IV3wc","Invariant mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV3rc = new TH1F("IV3rc","Invariant mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max);
  IV3wc->SetLineColor(6);
  IV3rc->SetLineColor(3);

  TH1F *IV4 = new TH1F("IV4","Invariant mass #mu^{+}_{2} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV4wc = new TH1F("IV4wc","Invariant mass #mu^{+}_{2} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max);
  TH1F *IV4rc = new TH1F("IV4rc","Invariant mass #mu^{+}_{2} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max);
  
  IV4wc->SetLineColor(6);
  IV4rc->SetLineColor(3);

  IV1->SetXTitle("IV1 = Invariant mass of (#mu_{1}^{+},#mu_{1}^{-})");
  IV2->SetXTitle("IV2 = Invariant mass of (#mu_{2}^{+},#mu_{2}^{-})");
  IV3->SetXTitle("IV3 = Invariant mass of (#mu_{1}^{+},#mu_{2}^{-})");
  IV4->SetXTitle("IV4 = Invariant mass of (#mu_{2}^{+},#mu_{1}^{-})");

  TH2F *IV21 = new TH2F("IV21","mass #mu^{+}_{2} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV2w1w = new TH2F("IV2w1w","mass #mu^{+}_{2} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV2r1r = new TH2F("IV2r1r","mass #mu^{+}_{2} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  IV2w1w->SetMarkerStyle(5);
  IV2w1w->SetMarkerColor(2);//red
  //  IV2w1w->SetMarkerSize(0.5);
  IV2r1r->SetMarkerStyle(2);
  IV2r1r->SetMarkerColor(4);
  //  IV2r1r->SetMarkerSize(0.5);
  
  TH2F *IV43 = new TH2F("IV43","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV4w3w = new TH2F("IV4w3w","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV4r3r = new TH2F("IV4r3r","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{1} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  IV4w3w->SetMarkerStyle(5);
  IV4w3w->SetMarkerColor(6);//pink
  //  IV4w3w->SetMarkerSize(0.5);
  IV4r3r->SetMarkerStyle(2);
  IV4r3r->SetMarkerColor(5);//green
  //  IV4r3r->SetMarkerSize(0.5);


  TH2F *IV31 = new TH2F("IV31","mass #mu^{+}_{1} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV3r1w = new TH2F("IV3r1w","mass #mu^{+}_{1} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV3w1r = new TH2F("IV3w1r","mass #mu^{+}_{1} #mu_{2}^{-} versus mass #mu^{+}_{1} #mu_{1}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  //  IV3r1w->SetMarkerSize(0.5);
  //  IV3w1r->SetMarkerSize(0.5);


  TH2F *IV42 = new TH2F("IV42","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV4r2w = new TH2F("IV4r2w","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  TH2F *IV4w2r = new TH2F("IV4w2r","mass #mu^{+}_{2} #mu_{1}^{-} versus mass #mu^{+}_{2} #mu_{2}^{-}",NbinIVplot,0,Higgsmass_max,NbinIVplot,0,Higgsmass_max);
  //  IV4r2w->SetMarkerSize(0.5);
  //  IV4w2r->SetMarkerSize(0.5);

  IV21->SetXTitle("IV2 = Invariant mass of (#mu_{2}^{+},#mu_{2}^{-})");
  IV21->SetYTitle("IV1 = Invariant mass of (#mu_{1}^{+},#mu_{1}^{-})");

  IV43->SetXTitle("IV4 = Invariant mass of (#mu_{2}^{+},#mu_{1}^{-})");
  IV43->SetYTitle("IV3 = Invariant mass of (#mu_{1}^{+},#mu_{2}^{-})");

  IV31->SetXTitle("IV3 = Invariant mass of (#mu_{1}^{+},#mu_{2}^{-})");
  IV31->SetYTitle("IV1 = Invariant mass of (#mu_{1}^{+},#mu_{1}^{-})");

  IV42->SetXTitle("IV4 = Invariant mass of (#mu_{2}^{+},#mu_{1}^{-})");
  IV42->SetYTitle("IV2 = Invariant mass of (#mu_{2}^{+},#mu_{2}^{-})");

  */
  Int_t KF_muon =13;



  // Now generate events

  TRandom2 R;
  int counting_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();
      // 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

      
      if(counting_num<1)
        {pythia->Pylist(1);
        }
      counting_num++;
      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 ( (KF==23) && (KS==21))  //Z boson
            {
              TLorentzVector lvZ = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4));
              if (Nz==0)
                {               mcZ1mass->Fill(lvZ.M());              Nz=1;}
              else
                mcZ2mass->Fill(lvZ.M());
            }

           if((abs(KF) == 22) &&(KS!=21)) //photon
            {
              
              TLorentzVector lvGamma = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4));
              if(lvGamma.Pt() >= 10.)
                hgamma->Fill(lvGamma.Pt());
            }
          */
          
          
          if ((abs(KF)== KF_muon))
            {
              TLorentzVector lvMuon = TLorentzVector(pythia->GetP(itk,1), pythia->GetP(itk,2), pythia->GetP(itk,3), pythia->GetP(itk,4));
              
              cout<<"Hello"<<endl;

              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
      /*      TVector3 boostZ1,boostZ2;
      TLorentzVector Z1,Z2;
      Z1 = muon1plus + muon1minus;
      Z2 = muon2plus + muon2minus;
      //get the opposite of the 3-vector of the Zs
      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();

      TVector3 uZ1=-boostZ1.Unit();
      TVector3 uZ2=-boostZ2.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;
        }
          
      float IV_1=0,IV_2=0,IV_3=0,IV_4=0;

      /*      IV_1 = (muon1plus+muon1minus).M();
      IV_2 = (muon2plus+muon2minus).M();
      IV_3 = (muon1plus+muon2minus).M();
      IV_4 = (muon2plus+muon1minus).M();
     

      mu1p_pt->Fill(muon1plus.Pt());
      mu2p_pt->Fill(muon2plus.Pt());
      mu1m_pt->Fill(muon1minus.Pt());
      mu2m_pt->Fill(muon2minus.Pt());

      IV1->Fill(IV_1); 
      IV2->Fill(IV_2);      
      IV3->Fill(IV_3);      
      IV4->Fill(IV_4);
      */

      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;

      /*  if (condition)
        {
          //blue
          IV1wc->Fill(IV_1);
          IV2wc->Fill(IV_2);
          IV3rc->Fill(IV_3);
          IV4rc->Fill(IV_4);
          IV2w1w->Fill(IV_2,IV_1);
          IV4r3r->Fill(IV_4,IV_3);
          IV3r1w->Fill(IV_3,IV_1);
          IV4r2w->Fill(IV_4,IV_2);
        }
      else
        {
          //black
          IV1rc->Fill(IV_1);
          IV2rc->Fill(IV_2);
          IV3wc->Fill(IV_3);
          IV4wc->Fill(IV_4);
          IV2r1r->Fill(IV_2,IV_1);
          IV4w3w->Fill(IV_4,IV_3);
          IV3w1r->Fill(IV_3,IV_1);
          IV4w2r->Fill(IV_4,IV_2);
        }

      IV21->Fill(IV_2,IV_1);
      IV43->Fill(IV_4,IV_3);
      IV31->Fill(IV_3,IV_1);
      IV42->Fill(IV_4,IV_2);
      */

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

 
  
  
  /* TCanvas *c1 = new TCanvas("c1","**",800,800);
  c1->Divide(2,2);
  c1->SetFillColor(9); 
  c1->cd(1);
  IV1->Draw();
  IV1wc->Draw("same");
  IV1rc->Draw("same");

  (c1->cd(2))->SetLogy();
  IV2->Draw();
  IV2wc->Draw("same");
  IV2rc->Draw("same");

  c1->cd(3);
  IV3->Draw();
  IV3wc->Draw("same");
  IV3rc->Draw("same");

  c1->cd(4)->SetLogy();
  IV4->Draw();
  IV4wc->Draw("same");
  IV4rc->Draw("same");
  */
   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();
  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 *c2=new TCanvas("c2","**",800,800);
  c2->Divide(2,2);
  c2->SetFillColor(8);
  c2->cd(1);
  IV21->Draw();
  IV2r1r->Draw("same");
  IV2w1w->Draw("same");

  c2->cd(2);
  IV43->Draw();
  IV4r3r->Draw("same");
  IV4w3w->Draw("same");
  
  c2->cd(3);
  TH2F * IV3w1r_c = new TH2F(*IV3w1r);
  IV3w1r->SetMarkerStyle(2);
  IV3w1r_c->SetMarkerStyle(5);
  IV3w1r->SetMarkerColor(4); //blue
  IV3w1r_c->SetMarkerColor(6);//yellow

  TH2F * IV3r1w_c = new TH2F(*IV3r1w);
  IV3r1w->SetMarkerStyle(2);
  IV3r1w_c->SetMarkerStyle(5);
  IV3r1w->SetMarkerColor(2); //red
  IV3r1w_c->SetMarkerColor(3); //green

  IV3w1r->Draw();
  IV3r1w->Draw("same");
  IV3w1r_c->Draw("same");
  IV3r1w_c->Draw("same");
  IV31->Draw("same");

  c2->cd(4);
  TH2F * IV4w2r_c = new TH2F(*IV4w2r);
  IV4w2r->SetMarkerStyle(2);
  IV4w2r_c->SetMarkerStyle(5);
  IV4w2r->SetMarkerColor(4); //blue
  IV4w2r_c->SetMarkerColor(6); //yellow

  TH2F * IV4r2w_c = new TH2F(*IV4r2w);
  IV4r2w->SetMarkerStyle(2);
  IV4r2w_c->SetMarkerStyle(5);
  IV4r2w->SetMarkerColor(2); //red
  IV4r2w_c->SetMarkerColor(3); //green


  IV4w2r->Draw("");
  IV4r2w->Draw("same");
  IV4w2r_c->Draw("same");
  IV4r2w_c->Draw("same");
  IV42->Draw("same");
   */
  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");

  f_fit->SetParameters(Hmass->GetEntries(),Higgsmass,10);
  
  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);

  Hmass->Fit("f_theory","L+","",Higgsminplot,Higgsmaxplot);

  //  higgswidth = 3.5;

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

  sprintf(name,"histograms_HtoZZ_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-1 ; 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

---