dasu/UltraFastSim/Analysis.cc

#include "Analysis.h"

#include <iostream>
#include <vector>
#include <string>
using namespace std;

#include "Pythia.h"
using namespace Pythia8;

#include "UltraFastSim.h"

#include "TROOT.h"
#include "TFile.h"
#include "TH1F.h"
#include "TTree.h"

TROOT root("Mike's ROOT Session","PYTHIA Histograms");

Analysis::Analysis(char *r, Pythia *p, UltraFastSim *u, bool v) : runType(r), pythia(p), ufs(u), verbose_(v), iEvent(0) {
  string name(runType);
  name += ".root";
  outFile = TFile::Open(name.c_str(), "recreate");
  bjet_pt = new TH1D("bjet_pt","bJet P_{T} [GeV]",20000,-.5,9999.5);
  tree = new TTree("T","All");
  tree->Branch("_px",&_px,"_px/D");
  tree->Branch("_py",&_py,"_py/D");
  tree->Branch("_pz",&_pz,"_pz/D");
  tree->Branch("_e",&_e,"_e/D");
  tree->Branch("_m",&_m,"_m/D");
  tree->Branch("_id",&_id,"_id/I");
  tree->Branch("_np",&_np,"_np/I");
  tree->Branch("_mother1",&_mother1,"_mother1/I");
  tree->Branch("_mother2",&_mother2,"_mother2/I");
  higgs_inv_mass = new TH1D("higgs_invm","Higgs Invariant Mass [GeV]", 20000, -0.5, 9999.5);
 
  cuttree = new TTree("T_btagged","after b-tagging");
  cuttree->Branch("_px",&_pxc,"_pxc/D");
  cuttree->Branch("_py",&_pyc,"_pyc/D");
  cuttree->Branch("_pz",&_pzc,"_pzc/D");
  cuttree->Branch("_e",&_ec,"_ec/D");
  cuttree->Branch("_m",&_mc,"_mc/D");
  cuttree->Branch("_id",&_idc,"_idc/I");
  cuttree->Branch("_np",&_npc,"_npc/I");
  cuttree->Branch("_mother1",&_mother1c,"_mother2c/I");
  cuttree->Branch("_mother2",&_mother1c,"_mother2c/I");
  higgs_inv_mass_cut = new TH1D("higgs_inv_mcut","Higgs Invariant Mass After Cut [GeV]", 10000, -0.5, 999.5);
}

Analysis::~Analysis() {
  // Save output 
  bjet_pt->Write();
  higgs_inv_mass->Write();
  higgs_inv_mass_cut->Write();
  tree->Write();
  cuttree->Write();
  outFile->Close();
}

bool Analysis::run() {
  //ml added ------------------------------------------------------------// 
  double bjet_ptcut(25.); //bjet pt cut.
  int bjet_counter25(0);
  int muon_counter(0),hadron_counter(0), tau_counter(0);
  int mh_size = pythia->event.size(); //this is not a good idea as it will crash sometimes...
  double hadron_store[mh_size][4]; //obviously not a good solution... to set this..
  double muon_store[mh_size][4];
   
  int tau_daughters[mh_size][2];

  double Einv=0.;
  double Pxinv=0.;
  double Pyinv=0.;
  double Pzinv=0.;

  for(int i=0; i<mh_size; i++){
    for(int j=0; j<4; j++){
      muon_store[i][j]=0.;
      hadron_store[i][j]=0.;
    }
    for(int j=0; j<2;j++){
      tau_daughters[i][j]=0;
    }

  }
  _np=0;
  _npc=0;

  vector<TLorentzVector> sortedbJets = ufs->bJetList();
  // below is to work out what taulist and gentaulist is
  int countertau_fin(0), countertau_nf(0), countermu_fin(0), countermu_nf(0);
  for (int i = 0; i < pythia->event.size(); ++i){
    Particle& theParticle = pythia->event[i];
    //if(abs(theParticle.id())==15){cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;}
      
    // if((theParticle.id())==15 && theParticle.isFinal()){ countertau_fin++;}
    // else  if(abs(theParticle.id())==15 && !theParticle.isFinal() ){countertau_nf++;  }
    if(abs(theParticle.id())==15){ 
      if(verbose_) cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;
      tau_daughters[tau_counter][0]=theParticle.daughter1();
      tau_daughters[tau_counter][1]=theParticle.daughter2();
      tau_counter++;
    }

    if((theParticle.id())==13 && theParticle.isFinal()){countermu_fin++;}
    else if (abs(theParticle.id())==13 && !theParticle.isFinal()){countermu_nf++;}

    //if(abs(theParticle.id())==13){cout<<theParticle.mother1()<<" mu mothers "<<theParticle.mother2()<<endl;}
  }
 
  if(verbose_){  
    cout<<" tauList.size()= "<<ufs->tauList().size()
        <<" genTauList.size()= "<<ufs->genTauList().size()
        <<" muList.size()= "<<ufs->muonList().size()
        <<" countertau_fin "<< countertau_fin
        <<" countertau_nf  "<<countertau_nf
        <<" countermu_fin "<<countermu_fin
        <<" countermu_nf "<<countermu_nf
        <<endl;
  }

    
  // int indextau(0);
  if(ufs->bJetList().size()>1 && (ufs->genTauList().size()>1 && countermu_fin>0) )//basic cuts
    {   
      //the first tree before the b tagging requirements
      for (int i = 0; i < pythia->event.size(); ++i){
        Particle& theParticle = pythia->event[i];         
        //  if(theParticle.isFinal()){ //displays only final state particles
        _np++;
        _px=theParticle.px();
        _py=theParticle.py();
        _pz=theParticle.pz();
        _e=theParticle.e();
        _m=theParticle.m();
        _id=theParticle.id();
        _mother1=theParticle.mother1();
        _mother2=theParticle.mother2();
        tree->Fill();
            
        //this stores the muons and hadrons 4mom -- checks if mother of particle is tau+- daughter by looping over taus
        for(int i=0; i<tau_counter; i++){
          //checks whether the particle is a muon or not
          if( theParticle.mother1() == tau_daughters[i][0] || theParticle.mother1()==tau_daughters[i][1] ||  (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
            {           
              if(abs(theParticle.id())==13){ //if the particle is muon+-
                if(verbose_)cout<<theParticle.mother1()<<"  muon mothers  "<<theParticle.mother2()<<endl;
                muon_store[muon_counter][1]=theParticle.px();
                muon_store[muon_counter][2]=theParticle.py();
                muon_store[muon_counter][3]=theParticle.pz();
                muon_store[muon_counter][0]=theParticle.e();
                //cout<<muon_store[muon_counter][1]<<endl; //some issue with storing the data
                muon_counter++;
              }
              else if( abs(theParticle.id())>18 && (abs(theParticle.id())<21 || abs(theParticle.id())>42) ){ //if the particle is neither lepton or sm gauge boson or single quark
                // cout<<theParticle.mother1()<<" hadron mothers "<<theParticle.mother2()<<endl;
                hadron_store[hadron_counter][1]=theParticle.px();
                hadron_store[hadron_counter][2]=theParticle.py();
                hadron_store[hadron_counter][3]=theParticle.pz();
                hadron_store[hadron_counter][0]=theParticle.e();
                hadron_counter++;      
              } 
            }
        }//end loop over tau
        //  }//end loop over final particles only
      }//end loop over particles in event

      if(verbose_) cout<<"Muon counter "<<muon_counter<<" had counter "<<hadron_counter<<endl;
      //find invariant mass for all combos
      for(int i=0; i<muon_counter; i++){
        for(int j=0; j<hadron_counter; j++){
          Einv = muon_store[i][0]+hadron_store[j][0];
          Pxinv= muon_store[i][1]+hadron_store[j][1];
          Pyinv= muon_store[i][2]+hadron_store[j][2];
          Pzinv= muon_store[i][3]+hadron_store[j][3];
          //cout<<" inv mass "<<Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv<<endl; 
          higgs_inv_mass->Fill(Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv); //fills the inv mass for the higgs  
        }
      }

      //-----------------NEXT DO FOR ADDITIONAL CUTS and repeat -------------
      for(int i=0; i<mh_size; i++){
        for(int j=0; j<4; j++){
          muon_store[i][j]=0.;
          hadron_store[i][j]=0.;
        }
      }
      hadron_counter=0;
      muon_counter=0;


      //checks the number of bjets with pt over cut
      for(unsigned int bi =0; bi<ufs->bJetList().size(); bi++){
        if(abs(sortedbJets[bi].Et())>bjet_ptcut) bjet_counter25++; //this counts the bjets that pass cutt
        bjet_pt->Fill(abs(sortedbJets[bi].Et()));
      }

      //implement the cut on event of btagging
      if(bjet_counter25>1){
        if(verbose_) cout<<" event passed btagging selection!!!"<<endl;
        for(int p=0; p< pythia->event.size(); p++){
          Particle& theParticle=pythia->event[p];     
          //if(theParticle.status()>0){
          _npc++;
          _pxc=theParticle.px();
          _pyc=theParticle.py();
          _pzc=theParticle.pz();
          _ec =theParticle.e();
          _mc = theParticle.m();
          _idc=theParticle.id();
          _mother1c=theParticle.mother1();
          _mother2c=theParticle.mother2();
          cuttree->Fill();
              
          for(int i=0; i<tau_counter; i++){
            if( theParticle.mother1() == tau_daughters[i][0] || theParticle.mother1()==tau_daughters[i][1] ||  (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
              {
                    
                if(abs(theParticle.id())==13){ //if the particle is muon+-
                  muon_store[muon_counter][1]=theParticle.px();
                  muon_store[muon_counter][2]=theParticle.py();
                  muon_store[muon_counter][3]=theParticle.pz();
                  muon_store[muon_counter][0]=theParticle.e();
                  muon_counter++;
                }
                else if( (abs(theParticle.id())<11 || abs(theParticle.id())>16) && (abs(theParticle.id())<21 || abs(theParticle.id())>25)){ //if the particle is neither lepton or sm gauge boson
                  hadron_store[hadron_counter][1]=theParticle.px();
                  hadron_store[hadron_counter][2]=theParticle.py();
                  hadron_store[hadron_counter][3]=theParticle.pz();
                  hadron_store[hadron_counter][0]=theParticle.e();
                  hadron_counter++;      
                }
              } 
          }//loop over taus
          //} //loop over only final state particles
        }//loop over particles in one event
      }//bjet cut loop
    }//basic cut loop
    
  //find invariant mass for all combos after cuts
  for(int i=0; i<muon_counter; i++){
    for(int j=0; j<hadron_counter; j++){
      Einv = muon_store[i][0]+hadron_store[j][0];
      Pxinv= muon_store[i][1]+hadron_store[j][1];
      Pyinv= muon_store[i][2]+hadron_store[j][2];
      Pzinv= muon_store[i][3]+hadron_store[j][3];
      higgs_inv_mass_cut->Fill(Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv); //fills the inv mass for the higgs  
    }
  }

  if(verbose_) {
    cout << "Number of Gen Chrgd = " << ufs->chargedHadronList().size() << endl;
    cout << "Number of Gen Neutr = " << ufs->neutralHadronList().size() << endl;
    cout << "Number of Gen Phtns = " << ufs->photonList().size() << endl;
    cout << "Number of Gen Elcns = " << ufs->electronList().size() << endl;
    cout << "Number of Gen Muons = " << ufs->muonList().size() << endl;
    cout << "Number of Gen Taus =  " << ufs->genTauList().size() << endl;
    cout << "Number of taus =      " << ufs->tauList().size() << endl;
    cout << "Number of c Quarks =  " << ufs->cQuarkList().size() << endl;
    cout << "Number of b Quarks =  " << ufs->bQuarkList().size() << endl;
    cout << "Number of Jets =      " << ufs->jetList().size() << endl;
    cout << "Number of bJets =     " << ufs->bJetList().size() << endl;
    cout << endl;
  }

  iEvent++;
  return true;
}

Revision:	1.4
Committed:	Tue Mar 15 04:35:17 2011 UTC (14 years, 2 months ago) by dasu
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+0 -0 lines
State:	*FILE REMOVED*
Log Message:	Minor changes
#	User	Rev	Content
1	dasu	1.1	#include "Analysis.h"
2
3			#include <iostream>
4			#include <vector>
5			#include <string>
6	dasu	1.2	using namespace std;
7	dasu	1.1
8	dasu	1.2	#include "Pythia.h"
9			using namespace Pythia8;
10	dasu	1.1
11			#include "UltraFastSim.h"
12
13			#include "TROOT.h"
14			#include "TFile.h"
15			#include "TH1F.h"
16			#include "TTree.h"
17
18			TROOT root("Mike's ROOT Session","PYTHIA Histograms");
19
20	dasu	1.2	Analysis::Analysis(char r, Pythia p, UltraFastSim *u, bool v) : runType(r), pythia(p), ufs(u), verbose_(v), iEvent(0) {
21	dasu	1.1	string name(runType);
22			name += ".root";
23			outFile = TFile::Open(name.c_str(), "recreate");
24			bjet_pt = new TH1D("bjet_pt","bJet P_{T} [GeV]",20000,-.5,9999.5);
25			tree = new TTree("T","All");
26			tree->Branch("_px",&_px,"_px/D");
27			tree->Branch("_py",&_py,"_py/D");
28			tree->Branch("_pz",&_pz,"_pz/D");
29			tree->Branch("_e",&_e,"_e/D");
30			tree->Branch("_m",&_m,"_m/D");
31			tree->Branch("_id",&_id,"_id/I");
32			tree->Branch("_np",&_np,"_np/I");
33			tree->Branch("_mother1",&_mother1,"_mother1/I");
34			tree->Branch("_mother2",&_mother2,"_mother2/I");
35			higgs_inv_mass = new TH1D("higgs_invm","Higgs Invariant Mass [GeV]", 20000, -0.5, 9999.5);
36
37			cuttree = new TTree("T_btagged","after b-tagging");
38			cuttree->Branch("_px",&_pxc,"_pxc/D");
39			cuttree->Branch("_py",&_pyc,"_pyc/D");
40			cuttree->Branch("_pz",&_pzc,"_pzc/D");
41			cuttree->Branch("_e",&_ec,"_ec/D");
42			cuttree->Branch("_m",&_mc,"_mc/D");
43			cuttree->Branch("_id",&_idc,"_idc/I");
44			cuttree->Branch("_np",&_npc,"_npc/I");
45			cuttree->Branch("_mother1",&_mother1c,"_mother2c/I");
46			cuttree->Branch("_mother2",&_mother1c,"_mother2c/I");
47			higgs_inv_mass_cut = new TH1D("higgs_inv_mcut","Higgs Invariant Mass After Cut [GeV]", 10000, -0.5, 999.5);
48			}
49
50			Analysis::~Analysis() {
51			// Save output
52			bjet_pt->Write();
53			higgs_inv_mass->Write();
54			higgs_inv_mass_cut->Write();
55			tree->Write();
56			cuttree->Write();
57			outFile->Close();
58			}
59
60			bool Analysis::run() {
61			//ml added ------------------------------------------------------------//
62			double bjet_ptcut(25.); //bjet pt cut.
63			int bjet_counter25(0);
64			int muon_counter(0),hadron_counter(0), tau_counter(0);
65	dasu	1.2	int mh_size = pythia->event.size(); //this is not a good idea as it will crash sometimes...
66	dasu	1.1	double hadron_store[mh_size][4]; //obviously not a good solution... to set this..
67			double muon_store[mh_size][4];
68
69			int tau_daughters[mh_size][2];
70
71			double Einv=0.;
72			double Pxinv=0.;
73			double Pyinv=0.;
74			double Pzinv=0.;
75
76			for(int i=0; i<mh_size; i++){
77			for(int j=0; j<4; j++){
78			muon_store[i][j]=0.;
79			hadron_store[i][j]=0.;
80			}
81			for(int j=0; j<2;j++){
82			tau_daughters[i][j]=0;
83			}
84
85			}
86			_np=0;
87			_npc=0;
88
89	dasu	1.3	vector<TLorentzVector> sortedbJets = ufs->bJetList();
90	dasu	1.1	// below is to work out what taulist and gentaulist is
91			int countertau_fin(0), countertau_nf(0), countermu_fin(0), countermu_nf(0);
92	dasu	1.2	for (int i = 0; i < pythia->event.size(); ++i){
93			Particle& theParticle = pythia->event[i];
94	dasu	1.1	//if(abs(theParticle.id())==15){cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;}
95
96			// if((theParticle.id())==15 && theParticle.isFinal()){ countertau_fin++;}
97			// else if(abs(theParticle.id())==15 && !theParticle.isFinal() ){countertau_nf++; }
98			if(abs(theParticle.id())==15){
99			if(verbose_) cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;
100			tau_daughters[tau_counter][0]=theParticle.daughter1();
101			tau_daughters[tau_counter][1]=theParticle.daughter2();
102			tau_counter++;
103			}
104
105			if((theParticle.id())==13 && theParticle.isFinal()){countermu_fin++;}
106			else if (abs(theParticle.id())==13 && !theParticle.isFinal()){countermu_nf++;}
107
108			//if(abs(theParticle.id())==13){cout<<theParticle.mother1()<<" mu mothers "<<theParticle.mother2()<<endl;}
109			}
110
111			if(verbose_){
112			cout<<" tauList.size()= "<<ufs->tauList().size()
113			<<" genTauList.size()= "<<ufs->genTauList().size()
114			<<" muList.size()= "<<ufs->muonList().size()
115			<<" countertau_fin "<< countertau_fin
116			<<" countertau_nf "<<countertau_nf
117			<<" countermu_fin "<<countermu_fin
118			<<" countermu_nf "<<countermu_nf
119			<<endl;
120			}
121
122
123			// int indextau(0);
124			if(ufs->bJetList().size()>1 && (ufs->genTauList().size()>1 && countermu_fin>0) )//basic cuts
125			{
126			//the first tree before the b tagging requirements
127	dasu	1.2	for (int i = 0; i < pythia->event.size(); ++i){
128			Particle& theParticle = pythia->event[i];
129	dasu	1.1	// if(theParticle.isFinal()){ //displays only final state particles
130			_np++;
131			_px=theParticle.px();
132			_py=theParticle.py();
133			_pz=theParticle.pz();
134			_e=theParticle.e();
135			_m=theParticle.m();
136			_id=theParticle.id();
137			_mother1=theParticle.mother1();
138			_mother2=theParticle.mother2();
139			tree->Fill();
140
141			//this stores the muons and hadrons 4mom -- checks if mother of particle is tau+- daughter by looping over taus
142			for(int i=0; i<tau_counter; i++){
143			//checks whether the particle is a muon or not
144			if( theParticle.mother1() == tau_daughters[i][0] \|\| theParticle.mother1()==tau_daughters[i][1] \|\| (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
145			{
146			if(abs(theParticle.id())==13){ //if the particle is muon+-
147			if(verbose_)cout<<theParticle.mother1()<<" muon mothers "<<theParticle.mother2()<<endl;
148			muon_store[muon_counter][1]=theParticle.px();
149			muon_store[muon_counter][2]=theParticle.py();
150			muon_store[muon_counter][3]=theParticle.pz();
151			muon_store[muon_counter][0]=theParticle.e();
152			//cout<<muon_store[muon_counter][1]<<endl; //some issue with storing the data
153			muon_counter++;
154			}
155			else if( abs(theParticle.id())>18 && (abs(theParticle.id())<21 \|\| abs(theParticle.id())>42) ){ //if the particle is neither lepton or sm gauge boson or single quark
156			// cout<<theParticle.mother1()<<" hadron mothers "<<theParticle.mother2()<<endl;
157			hadron_store[hadron_counter][1]=theParticle.px();
158			hadron_store[hadron_counter][2]=theParticle.py();
159			hadron_store[hadron_counter][3]=theParticle.pz();
160			hadron_store[hadron_counter][0]=theParticle.e();
161			hadron_counter++;
162			}
163			}
164			}//end loop over tau
165			// }//end loop over final particles only
166			}//end loop over particles in event
167
168			if(verbose_) cout<<"Muon counter "<<muon_counter<<" had counter "<<hadron_counter<<endl;
169			//find invariant mass for all combos
170			for(int i=0; i<muon_counter; i++){
171			for(int j=0; j<hadron_counter; j++){
172			Einv = muon_store[i][0]+hadron_store[j][0];
173			Pxinv= muon_store[i][1]+hadron_store[j][1];
174			Pyinv= muon_store[i][2]+hadron_store[j][2];
175			Pzinv= muon_store[i][3]+hadron_store[j][3];
176			//cout<<" inv mass "<<EinvEinv - PxinvPxinv - PyinvPyinv - PzinvPzinv<<endl;
177			higgs_inv_mass->Fill(EinvEinv - PxinvPxinv - PyinvPyinv - PzinvPzinv); //fills the inv mass for the higgs
178			}
179			}
180
181			//-----------------NEXT DO FOR ADDITIONAL CUTS and repeat -------------
182			for(int i=0; i<mh_size; i++){
183			for(int j=0; j<4; j++){
184			muon_store[i][j]=0.;
185			hadron_store[i][j]=0.;
186			}
187			}
188			hadron_counter=0;
189			muon_counter=0;
190
191
192			//checks the number of bjets with pt over cut
193	dasu	1.2	for(unsigned int bi =0; bi<ufs->bJetList().size(); bi++){
194	dasu	1.3	if(abs(sortedbJets[bi].Et())>bjet_ptcut) bjet_counter25++; //this counts the bjets that pass cutt
195			bjet_pt->Fill(abs(sortedbJets[bi].Et()));
196	dasu	1.1	}
197
198			//implement the cut on event of btagging
199			if(bjet_counter25>1){
200			if(verbose_) cout<<" event passed btagging selection!!!"<<endl;
201	dasu	1.2	for(int p=0; p< pythia->event.size(); p++){
202			Particle& theParticle=pythia->event[p];
203	dasu	1.1	//if(theParticle.status()>0){
204			_npc++;
205			_pxc=theParticle.px();
206			_pyc=theParticle.py();
207			_pzc=theParticle.pz();
208			_ec =theParticle.e();
209			_mc = theParticle.m();
210			_idc=theParticle.id();
211			_mother1c=theParticle.mother1();
212			_mother2c=theParticle.mother2();
213			cuttree->Fill();
214
215			for(int i=0; i<tau_counter; i++){
216			if( theParticle.mother1() == tau_daughters[i][0] \|\| theParticle.mother1()==tau_daughters[i][1] \|\| (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
217			{
218
219			if(abs(theParticle.id())==13){ //if the particle is muon+-
220			muon_store[muon_counter][1]=theParticle.px();
221			muon_store[muon_counter][2]=theParticle.py();
222			muon_store[muon_counter][3]=theParticle.pz();
223			muon_store[muon_counter][0]=theParticle.e();
224			muon_counter++;
225			}
226			else if( (abs(theParticle.id())<11 \|\| abs(theParticle.id())>16) && (abs(theParticle.id())<21 \|\| abs(theParticle.id())>25)){ //if the particle is neither lepton or sm gauge boson
227			hadron_store[hadron_counter][1]=theParticle.px();
228			hadron_store[hadron_counter][2]=theParticle.py();
229			hadron_store[hadron_counter][3]=theParticle.pz();
230			hadron_store[hadron_counter][0]=theParticle.e();
231			hadron_counter++;
232			}
233			}
234			}//loop over taus
235			//} //loop over only final state particles
236			}//loop over particles in one event
237			}//bjet cut loop
238			}//basic cut loop
239
240			//find invariant mass for all combos after cuts
241			for(int i=0; i<muon_counter; i++){
242			for(int j=0; j<hadron_counter; j++){
243			Einv = muon_store[i][0]+hadron_store[j][0];
244			Pxinv= muon_store[i][1]+hadron_store[j][1];
245			Pyinv= muon_store[i][2]+hadron_store[j][2];
246			Pzinv= muon_store[i][3]+hadron_store[j][3];
247			higgs_inv_mass_cut->Fill(EinvEinv - PxinvPxinv - PyinvPyinv - PzinvPzinv); //fills the inv mass for the higgs
248			}
249			}
250
251	dasu	1.2	if(verbose_) {
252			cout << "Number of Gen Chrgd = " << ufs->chargedHadronList().size() << endl;
253			cout << "Number of Gen Neutr = " << ufs->neutralHadronList().size() << endl;
254			cout << "Number of Gen Phtns = " << ufs->photonList().size() << endl;
255			cout << "Number of Gen Elcns = " << ufs->electronList().size() << endl;
256	dasu	1.1	cout << "Number of Gen Muons = " << ufs->muonList().size() << endl;
257			cout << "Number of Gen Taus = " << ufs->genTauList().size() << endl;
258	dasu	1.2	cout << "Number of taus = " << ufs->tauList().size() << endl;
259			cout << "Number of c Quarks = " << ufs->cQuarkList().size() << endl;
260	dasu	1.1	cout << "Number of b Quarks = " << ufs->bQuarkList().size() << endl;
261			cout << "Number of Jets = " << ufs->jetList().size() << endl;
262			cout << "Number of bJets = " << ufs->bJetList().size() << endl;
263			cout << endl;
264			}
265	dasu	1.2
266			iEvent++;
267	dasu	1.1	return true;
268			}
269