yangyong/zSelection/testSelectionZee.C

#include "rootheader.h"
#include "roofitheader.h"
#include "RecoAnalyzer.h"


#ifdef __MAKECINT__
#pragma link C++ class vector<vector<float> >+;
#pragma link C++ class vector<vector<unsigned short> >+;
#endif


#include "TMVAGui.C"

#if not defined(__CINT__) || defined(__MAKECINT__)
#include "TMVA/Tools.h"
#include "TMVA/Factory.h"
#include "TMVA/Reader.h"
#endif
using namespace TMVA;


#include "testSelection.h"

#include "ecalGap.cc"
#include "egammMVACorrection.cc"

int debug_ = -2; 


#include "setBranchAddress_ele.cc"

#include "physUtils.cc"
#include "usefullcode.cc"

#include "getGoodLS.cc"

#include "electronSelection.cc"

///re-weighting code 
#include "LumiReWeighting.cc" 


void testSelectionZee(char *test_datasetname, int test_evtRange){
  
  
  rgen_ = new TRandom3(0);
  cout<< rgen_->Gaus(0,0.01) <<endl; 
  
  
  dataversion = "v3";
  datasetname = test_datasetname;
  evtRange = test_evtRange; 
  
  
  egammaMVACorrection_LoadWeights();
  loadecalGapCoordinates();
    
  fChain  =new TChain("Analysis");
  TString filename ; 

  /// those files are preselected with two electron Et > 25GeV
  filename = TString( Form("/castor/cern.ch/user/y/yangyong/data/Run2011A/HiggsAnalysis/dielectronSkimmed/testAnalysisZee.%s.%s.allconv1.vtxmethod2.presel5.vtxtmva0.r%d.root",dataversion.c_str(),datasetname.c_str(),evtRange));
  cout<<filename<<endl;
  fChain->Add(filename);
  
  
  TChain *fevtInfo = new TChain("evtInfo");
  fevtInfo->Add(filename);
  fevtInfo->SetBranchAddress("totalNumberofEvents", &totalNumberofEvents);
  fevtInfo->SetBranchAddress("preselectedNumberofEvents", &preselectedNumberofEvents);
  fevtInfo->SetBranchAddress("datasettag", &datasettag);
  
  fevtInfo->GetEntry(0);
  string dname = datasettag->at(0);
  string dv = datasettag->at(1);
  int ntotal = totalNumberofEvents; 
  int nskimmed = preselectedNumberofEvents; 
  for(int n=1; n< fevtInfo->GetEntries(); n++){
    fevtInfo->GetEntry(n);
    ntotal += totalNumberofEvents; 
    nskimmed += preselectedNumberofEvents; 
  }
  
  
  setBranchAddress_ele();
  
  totalEntries = fChain->GetEntries();
  
  cout<<"totalEntries " << totalEntries <<endl; 
  
  ///get json file good lumis for each run 
  vector<string> certFile; 
  certFile.push_back("JsonFile/Cert_160404-163869_7TeV_May10ReReco_Collisions11_CMSSWConfig_v2.txtv1");
  certFile.push_back("JsonFile/Cert_160404-167913_7TeV_PromptReco_Collisions11_CMSSWConfig.txt.afterMay10");
  certFile.push_back("JsonFile/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_CMSSWConfig.txtv1");
  certFile.push_back("JsonFile/Cert_160404-173692_7TeV_PromptReco_Collisions11_CMSSWConfig.txtv1.run172620up");
  getLSrangeofEachRuns(certFile);

  
  TH1D *hh_pileupmc = new TH1D("hh_pileupmc","hh_npileup mc",51,-0.5,50.5);
  vector<double> weights_pileup;
  if( datasetname.find("SIM") != string::npos){
    /////for PU reweighting (added root files together)
    filename = TString( Form("PileUPDistrubtion/testAnalysis.%s.%s.puHist.root",dataversion.c_str(),datasetname.c_str()));
    TFile *fpumc = new TFile(filename,"read");
    TH1D *hh_npileupbx0 = (TH1D*)fpumc->Get("hh_npileupbx0");
    TH1D *hh_npileupbx1 = (TH1D*)fpumc->Get("hh_npileupbx1");
    TH1D *hh_npileupbxm1 = (TH1D*)fpumc->Get("hh_npileupbxm1");
    TH1D *hh_npileupbxav = (TH1D*)fpumc->Get("hh_npileupbxav"); //right now using averaged
    TFile *fpudata = new TFile("PileUPDistrubtion/Pileup_2011_EPS_8_jul.root","read");
    
    TH1D *hh_npileupdata = (TH1D*)fpudata->Get("pileup");
    weights_pileup = LumiReWeighting_getWeights(hh_npileupdata, hh_npileupbxav);
    fpumc->Close();
    fpudata->Close();
  }
  
  
  filename = TString (Form("selres/testSelectionZee.%s.%s.r%d.root",dataversion.c_str(), datasetname.c_str(),evtRange));
  cout<<filename<<endl; 
  TFile *fnew = new TFile(filename,"recreate");
  
  TH1I *hh_nevents = new TH1I("hh_nevents","total and skimm events",2,0,2);
  hh_nevents->SetBinContent(1,ntotal);
  hh_nevents->SetBinContent(2,nskimmed);
  
  
  vector<string> mpair_var;
  mpair_var.push_back("mpair_ebeb");
  mpair_var.push_back("mpair_ebeb_highr9");
  mpair_var.push_back("mpair_ebeb_lowr9");
  mpair_var.push_back("mpair_eeee");
  mpair_var.push_back("mpair_eeee_highr9");
  mpair_var.push_back("mpair_eeee_lowr9");
  mpair_var.push_back("mpair_eeee_highr9_loweps"); //eps_escraw < 0.05
  mpair_var.push_back("mpair_eeee_highr9_higheps");//eps_escraw > 0.05
  mpair_var.push_back("mpair_eeee_lowr9_loweps"); //eps_escraw < 0.05
  mpair_var.push_back("mpair_eeee_lowr9_higheps");//eps_escraw > 0.05
  
  
  RooRealVar *rv_mass = new RooRealVar("rv_mass","mass",100,0,1000);
  RooRealVar *rv_weight = new RooRealVar("rv_weight","weight",1.0,0,1E6);
  map<string,RooDataSet*> rhs_map_mpair;
  map<string,RooDataSet*> rhs_map_mpair_corr;
  for(int j=0; j< int( mpair_var.size()); j++){
    string mpairname = mpair_var[j];
    TString rname = TString( Form("rhs_%s",mpairname.c_str()) );
    rhs_map_mpair[mpairname] = new RooDataSet(rname,rname,RooArgList(*rv_mass,*rv_weight),rv_weight->GetName());
    mpairname = string (Form("%s_corr",mpair_var[j].c_str()));
    rname = TString( Form("rhs_%s",mpairname.c_str()) );
    rhs_map_mpair_corr[mpairname] = new RooDataSet(rname,rname,RooArgList(*rv_mass,*rv_weight),rv_weight->GetName());
  }
  map<string,TH1F*> hh_map_mpair; 
  map<string,TH1F*> hh_map_mpaircorr; 
  for(int j=0; j< int( mpair_var.size()); j++){
    string mpairname = mpair_var[j];
    filename = TString(Form("hh_%s",mpairname.c_str()));
    TH1F *hhtmp = new TH1F(filename,filename,2000,0,200);
    hh_map_mpair[mpairname] = hhtmp;
    hh_map_mpair[mpairname] ->Sumw2();
    
    mpairname = string (Form("%s_corr",mpair_var[j].c_str()));
    filename = TString(Form("hh_%s",mpairname.c_str())); 
    TH1F *hhtmp2 = new TH1F(filename,filename,2000,0,200);
    hh_map_mpaircorr[mpairname] = hhtmp2;
    hh_map_mpaircorr[mpairname] ->Sumw2();
  }
  
  
  bool goodCurLumiBlock = false; 
  int curLumiBlock = -1; 
  int curRun = -1; 
  vector<int> runNotUsed; 
  
  
  float en[2];
  float encorr[2];
  float pt[2];
  float eta[2];
  float phi[2];
  float res[20];
  float scr9[2];
  float scrps[2];
  
  float mpaircorr; 
  
  
  int nsel = 0; 
  
  
  ////loop over all events ( di-photon skimmed)
  for(entry = 0; entry < totalEntries; entry ++){
    fChain->GetEntry(entry);
    if(entry % 1000 ==0 ) cout<<" entry " <<entry<<endl; 
    
    
    ///JSON file good run and LS
    if( isRealData ){
      vector<int>::iterator it = find(goodRunList.begin(),goodRunList.end(),runNumber); 
      if( it == goodRunList.end()){
        vector<int>::iterator it2 = find(runNotUsed.begin(),runNotUsed.end(),runNumber); 
        if( it2 == runNotUsed.end()){
          runNotUsed.push_back(runNumber); 
        }
        continue; 
      }
      if( curLumiBlock != lumiBlock || curRun != runNumber){ /// a new lumiBlock  or starting of a new Run
        curLumiBlock = lumiBlock; 
        curRun =  runNumber; 
        goodCurLumiBlock = checkLumiBlockofRun();  //check this lumiBlock
      }
      if( ! goodCurLumiBlock) continue; 
      
    }
    
    for(int j=0; j< nElectron; j++){
      electronecalEnergy[j] = electronscenergy[j];  //right now ecalEnerg = scEnergy  
      float corr = getElectronCorrection(j); ///correction from regression 
      electronecalEnergyCorr[j] = (electronscrawEnergy[j] + electronscpreshowerEnergy[j]) * corr;
    }
    
    
    vector<int> indsel; 
    for(int j=0; j< nElectron; j++){
      if( electronecalDrivenSeed[j] ==0) continue; 
      float et = electronscenergy[j] * sin(2*atan(exp(-electroneta[j])));
      if( et < 30 ) continue;  ///ET > 30 
      if(! passElectronID(j,3)) continue; ///right now WP80 
      indsel.push_back(j);
    }
    if( int(indsel.size()) <2) continue; 
    
    nsel ++; 
    
    
    for(int j=0; j<2;j++){
      int k = indsel[j];
      en[j] = electronecalEnergy[k];
      encorr[j] = electronecalEnergyCorr[k];
      eta[j] = electrongsfTracketa[k];
      phi[j] = electrongsfTrackphi[k];
      scr9[j] = electrone3x3[k]/electronscrawEnergy[k];
      scrps[j] = electronscpreshowerEnergy[k] / electronscrawEnergy[k];
    }
    
    calcPairObjects(11,11,en,eta,phi,res);
    mpair = res[0];
    calcPairObjects(11,11,encorr,eta,phi,res);
    mpaircorr = res[0];
    
    
    double evtweight =1; 


    ///get pile up weghit 
    if( ! isRealData){
      int npubxav = 0; 
      for(unsigned int j =0; j < pileupBunchX->size(); j++){
        int BX = pileupBunchX->at(j);
        int npv = pileupNInteraction->at(j);
        if( BX >=-1 && BX <= 1){
          npubxav += npv; 
        }
      }
      npvbxav = npubxav*1.0/3;
      int bin_weight = hh_pileupmc->FindFixBin(npvbxav) -1;
      double puwt = LumiReWeighting_weightINT_v2(bin_weight, weights_pileup);
      
      evtweight = puwt; 
      
    }
    
    
    int ind1 = indsel[0];
    int ind2 = indsel[1];
    int bothEB =  fabs(electronsceta[ind1]) < 1.482 && fabs(electronsceta[ind2]) < 1.482;
    int bothEE =  fabs(electronsceta[ind1]) > 1.482 && fabs(electronsceta[ind2]) > 1.482;
    int bothHighR9 = scr9[0] > 0.94 && scr9[1] > 0.94; 
    int bothLowR9 = scr9[0] < 0.94 && scr9[1] < 0.94; 
    int bothLowEPS = scrps[0] < 0.05 && scrps[1] < 0.05; 
    int bothHighEPS = scrps[0] > 0.05 && scrps[1] > 0.05; 
    
    
    bool fillflag[20] = {bothEB,
                         bothEB && bothHighR9,
                         bothEB && bothLowR9, 
                         bothEE, 
                         bothEE && bothHighR9, 
                         bothEE && bothLowR9,
                         bothEE && bothHighR9 && bothLowEPS, 
                         bothEE && bothHighR9 && bothHighEPS,
                         bothEE && bothLowR9 && bothLowEPS, 
                         bothEE && bothLowR9 && bothHighEPS};
    
    for(int j=0; j< int( mpair_var.size()); j++){
      if( fillflag[j]==false) continue;
      string mpairname = mpair_var[j];
      rv_mass->setVal(mpair);
      hh_map_mpair[mpairname]->Fill(mpair,evtweight);
      rhs_map_mpair[mpairname]->add(*rv_mass,evtweight);
      rv_mass->setVal(mpaircorr);
      mpairname = TString( Form("%s_corr",mpair_var[j].c_str()));
      hh_map_mpaircorr[mpairname]->Fill(mpaircorr,evtweight);
      rhs_map_mpair_corr[mpairname]->add(*rv_mass,evtweight);
    }
    
  }

  cout <<" nsel" << nsel <<endl; 
  
  
  fnew->cd();
  //RooContainer_Save();
  RooWorkspace *w = new RooWorkspace("zeeShape","workspace") ;
  for (std::map<string,RooDataSet*>::iterator it_data = rhs_map_mpair.begin()
         ;it_data != rhs_map_mpair.end();it_data++)     {
    w->import(*it_data->second);
  }
  for (std::map<string,RooDataSet*>::iterator it_data = rhs_map_mpair_corr.begin()
         ;it_data != rhs_map_mpair_corr.end();it_data++)        {
    w->import(*it_data->second);
  }
  w->Write();
  
  fnew->Write();
  fnew->Close();
}
Revision:	1.2
Committed:	Fri Sep 30 15:27:51 2011 UTC (13 years, 7 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-00-00
Changes since 1.1:	+3 -2 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	yangyong	1.1	#include "rootheader.h"
2			#include "roofitheader.h"
3			#include "RecoAnalyzer.h"
4
5
6			#ifdef __MAKECINT__
7			#pragma link C++ class vector<vector<float> >+;
8			#pragma link C++ class vector<vector<unsigned short> >+;
9			#endif
10
11
12
13			#include "TMVAGui.C"
14
15			#if not defined(__CINT__) \|\| defined(__MAKECINT__)
16			#include "TMVA/Tools.h"
17			#include "TMVA/Factory.h"
18			#include "TMVA/Reader.h"
19			#endif
20			using namespace TMVA;
21
22
23
24			#include "testSelection.h"
25
26			#include "ecalGap.cc"
27			#include "egammMVACorrection.cc"
28
29			int debug_ = -2;
30
31
32			#include "setBranchAddress_ele.cc"
33
34			#include "physUtils.cc"
35			#include "usefullcode.cc"
36
37			#include "getGoodLS.cc"
38
39			#include "electronSelection.cc"
40
41			///re-weighting code
42			#include "LumiReWeighting.cc"
43
44
45			void testSelectionZee(char *test_datasetname, int test_evtRange){
46
47
48			rgen_ = new TRandom3(0);
49			cout<< rgen_->Gaus(0,0.01) <<endl;
50
51
52			dataversion = "v3";
53			datasetname = test_datasetname;
54			evtRange = test_evtRange;
55
56
57			egammaMVACorrection_LoadWeights();
58			loadecalGapCoordinates();
59
60			fChain =new TChain("Analysis");
61			TString filename ;
62	yangyong	1.2
63			/// those files are preselected with two electron Et > 25GeV
64			filename = TString( Form("/castor/cern.ch/user/y/yangyong/data/Run2011A/HiggsAnalysis/dielectronSkimmed/testAnalysisZee.%s.%s.allconv1.vtxmethod2.presel5.vtxtmva0.r%d.root",dataversion.c_str(),datasetname.c_str(),evtRange));
65	yangyong	1.1	cout<<filename<<endl;
66			fChain->Add(filename);
67
68
69			TChain *fevtInfo = new TChain("evtInfo");
70			fevtInfo->Add(filename);
71			fevtInfo->SetBranchAddress("totalNumberofEvents", &totalNumberofEvents);
72			fevtInfo->SetBranchAddress("preselectedNumberofEvents", &preselectedNumberofEvents);
73			fevtInfo->SetBranchAddress("datasettag", &datasettag);
74
75			fevtInfo->GetEntry(0);
76			string dname = datasettag->at(0);
77			string dv = datasettag->at(1);
78			int ntotal = totalNumberofEvents;
79			int nskimmed = preselectedNumberofEvents;
80			for(int n=1; n< fevtInfo->GetEntries(); n++){
81			fevtInfo->GetEntry(n);
82			ntotal += totalNumberofEvents;
83			nskimmed += preselectedNumberofEvents;
84			}
85
86
87			setBranchAddress_ele();
88
89			totalEntries = fChain->GetEntries();
90
91			cout<<"totalEntries " << totalEntries <<endl;
92
93			///get json file good lumis for each run
94			vector<string> certFile;
95			certFile.push_back("JsonFile/Cert_160404-163869_7TeV_May10ReReco_Collisions11_CMSSWConfig_v2.txtv1");
96			certFile.push_back("JsonFile/Cert_160404-167913_7TeV_PromptReco_Collisions11_CMSSWConfig.txt.afterMay10");
97			certFile.push_back("JsonFile/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_CMSSWConfig.txtv1");
98			certFile.push_back("JsonFile/Cert_160404-173692_7TeV_PromptReco_Collisions11_CMSSWConfig.txtv1.run172620up");
99			getLSrangeofEachRuns(certFile);
100
101
102			TH1D *hh_pileupmc = new TH1D("hh_pileupmc","hh_npileup mc",51,-0.5,50.5);
103			vector<double> weights_pileup;
104			if( datasetname.find("SIM") != string::npos){
105			/////for PU reweighting (added root files together)
106			filename = TString( Form("PileUPDistrubtion/testAnalysis.%s.%s.puHist.root",dataversion.c_str(),datasetname.c_str()));
107			TFile *fpumc = new TFile(filename,"read");
108			TH1D hh_npileupbx0 = (TH1D)fpumc->Get("hh_npileupbx0");
109			TH1D hh_npileupbx1 = (TH1D)fpumc->Get("hh_npileupbx1");
110			TH1D hh_npileupbxm1 = (TH1D)fpumc->Get("hh_npileupbxm1");
111			TH1D hh_npileupbxav = (TH1D)fpumc->Get("hh_npileupbxav"); //right now using averaged
112			TFile *fpudata = new TFile("PileUPDistrubtion/Pileup_2011_EPS_8_jul.root","read");
113
114			TH1D hh_npileupdata = (TH1D)fpudata->Get("pileup");
115			weights_pileup = LumiReWeighting_getWeights(hh_npileupdata, hh_npileupbxav);
116			fpumc->Close();
117			fpudata->Close();
118			}
119
120
121			filename = TString (Form("selres/testSelectionZee.%s.%s.r%d.root",dataversion.c_str(), datasetname.c_str(),evtRange));
122			cout<<filename<<endl;
123			TFile *fnew = new TFile(filename,"recreate");
124
125			TH1I *hh_nevents = new TH1I("hh_nevents","total and skimm events",2,0,2);
126			hh_nevents->SetBinContent(1,ntotal);
127			hh_nevents->SetBinContent(2,nskimmed);
128
129
130			vector<string> mpair_var;
131			mpair_var.push_back("mpair_ebeb");
132			mpair_var.push_back("mpair_ebeb_highr9");
133			mpair_var.push_back("mpair_ebeb_lowr9");
134			mpair_var.push_back("mpair_eeee");
135			mpair_var.push_back("mpair_eeee_highr9");
136			mpair_var.push_back("mpair_eeee_lowr9");
137			mpair_var.push_back("mpair_eeee_highr9_loweps"); //eps_escraw < 0.05
138			mpair_var.push_back("mpair_eeee_highr9_higheps");//eps_escraw > 0.05
139			mpair_var.push_back("mpair_eeee_lowr9_loweps"); //eps_escraw < 0.05
140			mpair_var.push_back("mpair_eeee_lowr9_higheps");//eps_escraw > 0.05
141
142
143
144			RooRealVar *rv_mass = new RooRealVar("rv_mass","mass",100,0,1000);
145			RooRealVar *rv_weight = new RooRealVar("rv_weight","weight",1.0,0,1E6);
146			map<string,RooDataSet*> rhs_map_mpair;
147			map<string,RooDataSet*> rhs_map_mpair_corr;
148			for(int j=0; j< int( mpair_var.size()); j++){
149			string mpairname = mpair_var[j];
150			TString rname = TString( Form("rhs_%s",mpairname.c_str()) );
151			rhs_map_mpair[mpairname] = new RooDataSet(rname,rname,RooArgList(rv_mass,rv_weight),rv_weight->GetName());
152			mpairname = string (Form("%s_corr",mpair_var[j].c_str()));
153			rname = TString( Form("rhs_%s",mpairname.c_str()) );
154			rhs_map_mpair_corr[mpairname] = new RooDataSet(rname,rname,RooArgList(rv_mass,rv_weight),rv_weight->GetName());
155			}
156			map<string,TH1F*> hh_map_mpair;
157			map<string,TH1F*> hh_map_mpaircorr;
158			for(int j=0; j< int( mpair_var.size()); j++){
159			string mpairname = mpair_var[j];
160			filename = TString(Form("hh_%s",mpairname.c_str()));
161			TH1F *hhtmp = new TH1F(filename,filename,2000,0,200);
162			hh_map_mpair[mpairname] = hhtmp;
163			hh_map_mpair[mpairname] ->Sumw2();
164
165			mpairname = string (Form("%s_corr",mpair_var[j].c_str()));
166			filename = TString(Form("hh_%s",mpairname.c_str()));
167			TH1F *hhtmp2 = new TH1F(filename,filename,2000,0,200);
168			hh_map_mpaircorr[mpairname] = hhtmp2;
169			hh_map_mpaircorr[mpairname] ->Sumw2();
170			}
171
172
173
174			bool goodCurLumiBlock = false;
175			int curLumiBlock = -1;
176			int curRun = -1;
177			vector<int> runNotUsed;
178
179
180
181			float en[2];
182			float encorr[2];
183			float pt[2];
184			float eta[2];
185			float phi[2];
186			float res[20];
187			float scr9[2];
188			float scrps[2];
189
190			float mpaircorr;
191
192
193			int nsel = 0;
194
195
196
197			////loop over all events ( di-photon skimmed)
198			for(entry = 0; entry < totalEntries; entry ++){
199			fChain->GetEntry(entry);
200			if(entry % 1000 ==0 ) cout<<" entry " <<entry<<endl;
201
202
203			///JSON file good run and LS
204			if( isRealData ){
205			vector<int>::iterator it = find(goodRunList.begin(),goodRunList.end(),runNumber);
206			if( it == goodRunList.end()){
207			vector<int>::iterator it2 = find(runNotUsed.begin(),runNotUsed.end(),runNumber);
208			if( it2 == runNotUsed.end()){
209			runNotUsed.push_back(runNumber);
210			}
211			continue;
212			}
213			if( curLumiBlock != lumiBlock \|\| curRun != runNumber){ /// a new lumiBlock or starting of a new Run
214			curLumiBlock = lumiBlock;
215			curRun = runNumber;
216			goodCurLumiBlock = checkLumiBlockofRun(); //check this lumiBlock
217			}
218			if( ! goodCurLumiBlock) continue;
219
220			}
221
222			for(int j=0; j< nElectron; j++){
223			electronecalEnergy[j] = electronscenergy[j]; //right now ecalEnerg = scEnergy
224			float corr = getElectronCorrection(j); ///correction from regression
225			electronecalEnergyCorr[j] = (electronscrawEnergy[j] + electronscpreshowerEnergy[j]) * corr;
226			}
227
228
229			vector<int> indsel;
230			for(int j=0; j< nElectron; j++){
231			if( electronecalDrivenSeed[j] ==0) continue;
232			float et = electronscenergy[j] * sin(2*atan(exp(-electroneta[j])));
233			if( et < 30 ) continue; ///ET > 30
234			if(! passElectronID(j,3)) continue; ///right now WP80
235			indsel.push_back(j);
236			}
237			if( int(indsel.size()) <2) continue;
238
239			nsel ++;
240
241
242			for(int j=0; j<2;j++){
243			int k = indsel[j];
244			en[j] = electronecalEnergy[k];
245			encorr[j] = electronecalEnergyCorr[k];
246			eta[j] = electrongsfTracketa[k];
247			phi[j] = electrongsfTrackphi[k];
248			scr9[j] = electrone3x3[k]/electronscrawEnergy[k];
249			scrps[j] = electronscpreshowerEnergy[k] / electronscrawEnergy[k];
250			}
251
252			calcPairObjects(11,11,en,eta,phi,res);
253			mpair = res[0];
254			calcPairObjects(11,11,encorr,eta,phi,res);
255			mpaircorr = res[0];
256
257
258			double evtweight =1;
259
260
261			///get pile up weghit
262			if( ! isRealData){
263			int npubxav = 0;
264			for(unsigned int j =0; j < pileupBunchX->size(); j++){
265			int BX = pileupBunchX->at(j);
266			int npv = pileupNInteraction->at(j);
267			if( BX >=-1 && BX <= 1){
268			npubxav += npv;
269			}
270			}
271			npvbxav = npubxav*1.0/3;
272			int bin_weight = hh_pileupmc->FindFixBin(npvbxav) -1;
273			double puwt = LumiReWeighting_weightINT_v2(bin_weight, weights_pileup);
274
275			evtweight = puwt;
276
277			}
278
279
280			int ind1 = indsel[0];
281			int ind2 = indsel[1];
282			int bothEB = fabs(electronsceta[ind1]) < 1.482 && fabs(electronsceta[ind2]) < 1.482;
283			int bothEE = fabs(electronsceta[ind1]) > 1.482 && fabs(electronsceta[ind2]) > 1.482;
284			int bothHighR9 = scr9[0] > 0.94 && scr9[1] > 0.94;
285			int bothLowR9 = scr9[0] < 0.94 && scr9[1] < 0.94;
286			int bothLowEPS = scrps[0] < 0.05 && scrps[1] < 0.05;
287			int bothHighEPS = scrps[0] > 0.05 && scrps[1] > 0.05;
288
289
290
291			bool fillflag[20] = {bothEB,
292			bothEB && bothHighR9,
293			bothEB && bothLowR9,
294			bothEE,
295			bothEE && bothHighR9,
296			bothEE && bothLowR9,
297			bothEE && bothHighR9 && bothLowEPS,
298			bothEE && bothHighR9 && bothHighEPS,
299			bothEE && bothLowR9 && bothLowEPS,
300			bothEE && bothLowR9 && bothHighEPS};
301
302			for(int j=0; j< int( mpair_var.size()); j++){
303			if( fillflag[j]==false) continue;
304			string mpairname = mpair_var[j];
305			rv_mass->setVal(mpair);
306			hh_map_mpair[mpairname]->Fill(mpair,evtweight);
307			rhs_map_mpair[mpairname]->add(*rv_mass,evtweight);
308			rv_mass->setVal(mpaircorr);
309			mpairname = TString( Form("%s_corr",mpair_var[j].c_str()));
310			hh_map_mpaircorr[mpairname]->Fill(mpaircorr,evtweight);
311			rhs_map_mpair_corr[mpairname]->add(*rv_mass,evtweight);
312			}
313
314			}
315
316			cout <<" nsel" << nsel <<endl;
317
318
319			fnew->cd();
320			//RooContainer_Save();
321			RooWorkspace *w = new RooWorkspace("zeeShape","workspace") ;
322			for (std::map<string,RooDataSet*>::iterator it_data = rhs_map_mpair.begin()
323			;it_data != rhs_map_mpair.end();it_data++) {
324			w->import(*it_data->second);
325			}
326			for (std::map<string,RooDataSet*>::iterator it_data = rhs_map_mpair_corr.begin()
327			;it_data != rhs_map_mpair_corr.end();it_data++) {
328			w->import(*it_data->second);
329			}
330			w->Write();
331
332			fnew->Write();
333			fnew->Close();
334			}