C++ AnalysisEnvironmentLoader类代码示例

int main (int argc, char *argv[])
{
 
  cout << "#########################" << endl;
  cout << "Beginning of the program" << endl;
  cout << "#########################" << endl;

  //////////////////////
  //Global variables
  //////////////////////
  vector < Dataset > datasets;
  SSDiLeptonSelection sel;
  float Luminosity = 0;
  float LumiError = 0.;
  // 0: MC - 1: Data - 2 Data & MC
  int DataType = 0;
  int verbosity = -1;

  int AnaStep = 6;//which defines the cuts that the events should pass to be considered as selected

  float Nobs_mumu = 0.;
  //////////////////////

  //////////////////////
  // Initialisation
  //////////////////////
  string xmlFileName;
  cout<<"argc "<<argc<<" "<<argv[0]<<endl;
  if (argc>1 ) xmlFileName = string(argv[1]);
  else xmlFileName = string ("../../config/TTbarMETAnalysis.xml");
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets);	// now the list of datasets written in the xml file is known
  anaEL.LoadSSDiLeptonSelection (sel);	// now the parameters for the selection are given to the selection
  anaEL.LoadGeneralInfo (DataType, Luminosity, LumiError, verbosity);
  int flagOriginal=sel.GetFlagb();
  int methodOriginal=sel.GetMethodb();
  int systOriginal= sel.GetSystb();
  std::cout << " For btag : flag " << flagOriginal << ", method " << methodOriginal << ", syst " << systOriginal << std::endl;
  IPHCTree::NTEvent * event = 0;
  //Selection table
  SelectionTable selTable_allChannels (sel.GetCutList (), datasets, string ("*"));
  SelectionTable selTable_ee (sel.GetCutList (), datasets, string ("ee"));
  SelectionTable selTable_emu (sel.GetCutList (), datasets, string ("emu"));
  SelectionTable selTable_mumu (sel.GetCutList (), datasets, string ("mumu"));
  
/*
  //Book keeping of standard histos
  bool doHistoManager = true;
  SSDiLepAnaHistoManager histoManager;
  if(doHistoManager){
  	histoManager.LoadDatasets (datasets);
  	histoManager.LoadSelectionSteps (sel.GetCutList ());
  	histoManager.LoadChannels (sel.GetChannelList ());
  	histoManager.CreateHistos ();
  }
*/

  TH1F* h_mult_lep_0  =  new TH1F("h_mult_lep_0",     "# of leptons (bf cuts)",   5,0.0,5.0);
  TH1F* h_mult_jet_0  =  new TH1F("h_mult_jet_0",     "# of jets (bf cuts)",     10,0.0,10.0);
  TH1F* h_mult_bjet_0 =  new TH1F("h_mult_bjet_0",    "# of b-jets (bf cuts)",    5,0.0,5.0);
  TH1F* h_pt_lep1_0   =  new TH1F("h_pt_lep1_0",      "PT of leading lepton (bf cuts)",60,0.0,300.0);
  TH1F* h_eta_lep1_0  =  new TH1F("h_eta_lep1_0",     "Eta of leading lepton (bf cuts)",66,-3.3,3.3);
  TH1F* h_pt_jet1_0   =  new TH1F("h_pt_jet1_0",      "PT of leading jet (bf cuts)",60,0.0,300.0);
  TH1F* h_pt_jet2_0   =  new TH1F("h_pt_jet2_0",      "PT of jet 2 (bf cuts)",60,0.0,300.0);
  TH1F* h_pt_jet3_0   =  new TH1F("h_pt_jet3_0",      "PT of jet 3 (bf cuts)",60,0.0,300.0);
  TH1F* h_pt_jet4_0   =  new TH1F("h_pt_jet4_0",      "PT of jet 4 (bf cuts)",60,0.0,300.0);
  TH1F* h_met_0       =  new TH1F("h_met_0",          "MET (bf cuts)",40,0.0,400.0);
  TH1F* h_mt_0        =  new TH1F("h_mt_0",           "MT(invisible+lepton) (bf cuts)",100,0.0,1000.0);
  TH1F* h_mt_0_2      =  new TH1F("h_mt_0_2",         "MT(invisible+lepton) (bf cuts)",20,0.0,1000.0);
  TH1F* h_ht_0        =  new TH1F("h_ht_0",           "HT (bf cuts)",100,0.0,500.0);
  TH1F* h_m3_0        =  new TH1F("h_m3_0",           "M3 (bf cuts)",100,0.0,1000.0);

  TH1F* h_mult_jet_1  =  new TH1F("h_mult_jet_1",     "# of jets (af cuts)",     10,0.0,10.0);
  TH1F* h_met_1       =  new TH1F("h_met_1",          "MET (af cuts)",40,0.0,400.0);
  TH1F* h_met_1_2     =  new TH1F("h_met_1_2",        "MET (af cuts)",10,0.0,400.0);
  TH1F* h_mt_1        =  new TH1F("h_mt_1",           "MT(invisible+lepton) (af cuts)",100,0.0,1000.0);
  TH1F* h_mt_1_2      =  new TH1F("h_mt_1_2",         "MT(invisible+lepton) (af cuts)",20,0.0,1000.0);
  TH1F* h_pt_lep1_1   =  new TH1F("h_pt_lep1_1",      "PT of leading lepton (af cuts)",60,0.0,300.0);
  TH1F* h_eta_lep1_1  =  new TH1F("h_eta_lep1_1",     "Eta of leading lepton (af cuts)",66,-3.3,3.3);
  TH1F* h_pt_jet1_1   =  new TH1F("h_pt_jet1_1",      "PT of leading jet (af cuts)",60,0.0,300.0);
  TH1F* h_pt_jet2_1   =  new TH1F("h_pt_jet2_1",      "PT of jet 2 (af cuts)",60,0.0,300.0);
  TH1F* h_pt_jet3_1   =  new TH1F("h_pt_jet3_1",      "PT of jet 3 (af cuts)",60,0.0,300.0);
  TH1F* h_pt_jet4_1   =  new TH1F("h_pt_jet4_1",      "PT of jet 4 (af cuts)",60,0.0,300.0);
  TH1F* h_ht_1        =  new TH1F("h_ht_1",           "HT (af cuts)",100,0.0,500.0);
  TH1F* h_m3_1        =  new TH1F("h_m3_1",           "M3 (af cuts)",100,0.0,1000.0);
  TH1F* h_mult_bjet_1 =  new TH1F("h_mult_bjet_1",    "# of b-jets (af cuts)",    5,0.0,5.0);
  TH1F* h_pt_bjet1_2  =  new TH1F("h_pt_bjet1_2",     "PT of leading b-jet (af cuts)",60,0.0,300.0);

  TH1F* h_tophad_pt_1      =  new TH1F("h_tophad_pt_1",     "PT of top1 (af cuts)",100,0.,1000.);
  TH1F* h_toplep_pt_1      =  new TH1F("h_toplep_pt_1",     "PT of top2 (af cuts)",100,0.,1000.);
  TH1F* h_wlep_pt_1        =  new TH1F("h_wlep_pt_1",       "PT of w2 (af cuts)",100,0.,1000.);
  TH1F* h_toplep_m_1       =  new TH1F("h_toplep_m_1",      "M of top2 (af cuts)",100,0.,1000.);
  TH1F* h_dphi_l_mis_1     =  new TH1F("h_dphi_l_mis_1",    "Dphi(l,MET) (af cuts)",30,-1.*pi,pi);
  TH1F* h_dphi_l_jet4th_1  =  new TH1F("h_dphi_l_jet4th_1", "Dphi(l,jet4th) (af cuts)",30,-1.*pi,pi);
  TH1F* h_dphi_2top_1      =  new TH1F("h_dphi_2top_1",     "Dphi(top1,top2) (af cuts)",30,-1.*pi,pi);
  TH1F* h_deta_l_tophad_1  =  new TH1F("h_deta_l_tophad_1", "Deta(l,top1) (af cuts)",50,0.,5.);
  TH1F* h_deta_l_jet4th_1  =  new TH1F("h_deta_l_jet4th_1", "Deta(l,jet4th) (af cuts)",50,0.,5.);
  TH1F* h_tophad_pt_2      =  new TH1F("h_tophad_pt_2",     "PT of top1 (af cuts)",100,0.,1000.);
  TH1F* h_toplep_pt_2      =  new TH1F("h_toplep_pt_2",     "PT of top2 (af cuts)",100,0.,1000.);
  TH1F* h_wlep_pt_2        =  new TH1F("h_wlep_pt_2",       "PT of w2 (af cuts)",100,0.,1000.);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TH2F * histoDelta_yVsTTMassEfficiency = new TH2F("Delta_yVsTTMassEfficiency", "Delta_yVsTTMassEfficiency", 10, -2, 2, 5, 0, 1000);
  // Fine binning
  TH1F * histoFineBinning_Delta_yEfficiencyN = new TH1F("FineBinning_Delta_yEfficiencyN", "FineBinning_Delta_yEfficiencyN", 200, -2, 2);
  TH1F * histoFineBinning_Delta_yEfficiencyD = new TH1F("FineBinning_Delta_yEfficiencyD", "FineBinning_Delta_yEfficiencyD", 200, -2, 2);
  TH1F * histoFineBinning_Delta_yEfficiency = new TH1F("FineBinning_Delta_yEfficiency", "FineBinning_Delta_yEfficiency", 200, -2, 2);

  TH2F * histoFineBinning_Delta_yVsTTRapidityEfficiency = new TH2F("FineBinning_Delta_yVsTTRapidityEfficiency", "FineBinning_Delta_yVsTTRapidityEfficiency", 200, -2, 2, 500, -2.5, 2.5);
  TH2F * histoFineBinning_Delta_yVsTTPtEfficiency = new TH2F("FineBinning_Delta_yVsTTPtEfficiency", "FineBinning_Delta_yVsTTPtEfficiency", 200, -2, 2, 250, 0, 250);
  TH2F * histoFineBinning_Delta_yVsTTMassEfficiency = new TH2F("FineBinning_Delta_yVsTTMassEfficiency", "FineBinning_Delta_yVsTTMassEfficiency", 200, -2, 2, 100, 0, 1000);




  vector < Dataset > datasets;
  DiLeptonSelection sel;
  float Luminosity = 0;
  // 0: MC - 1: Data - 2 Data & MC
  int DataType = 0;
  int verbosity = -1;
  
  // Analysis variables
  int step;
  //4-vectors
  vector < TLorentzVector > nu1;
  vector < TLorentzVector > nu2;
  //  TLorentzVector lvTop1_aux;
  //TLorentzVector lvTop2_aux;

  //////////////////////

  //////////////////////
  // Initialisation
  //////////////////////
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets); // now the list of datasets written in the xml file is known
  anaEL.LoadDiLeptonSelection (sel); // now the parameters for the selection are given to the selection
  anaEL.LoadGeneralInfo(DataType, Luminosity, verbosity );

  float weight_MC; 
  int nEvents;
  int firstEvent;
  int initialEvents;
  anaEL.LoadAnalysisInfo("Run", "firstEvent", firstEvent);
  anaEL.LoadAnalysisInfo("Run", "nEvents", nEvents);
  anaEL.LoadAnalysisInfo("Run", "initialEvents", initialEvents);
  string channelName;
  anaEL.LoadAnalysisInfo("MassAnalysis", "ChannelName", channelName);

  anaEL.LoadWeight (sel); // now the parameters for SFBweight are initialized (for b-tag!)

  //******************************************
  //Load Scale Factors for lepton efficiencies
  //******************************************
  sel.LoadElScaleFactors();
  sel.LoadMuScaleFactors();
  sel.InitJESUnc();

  bool nonskimmed = false;

  TopTree::NTEvent * event = 0;

  TFile * TupleFile = new TFile("TupleFileSkimmed_TTbarSemileptonicPowheg.root","RECREATE");
  float mttbar_rec, mttbar_gen;
  float diffysquare_rec, diffysquare_gen;
  float diffabseta_rec, diffabseta_gen;
  float weight;

int main (int argc, char *argv[])
{
 
  cout << "#########################" << endl;
  cout << "Beginning of the program" << endl;
  cout << "#########################" << endl;

  //////////////////////
  //Global variables
  //////////////////////
  vector < Dataset > datasets;
  SSDiLeptonSelection sel;
  float Luminosity = 0;
  float LumiError = 0.;
  // 0: MC - 1: Data - 2 Data & MC
  int DataType = 0;
  int verbosity = -1;

  int AnaStep = 6;//which defines the cuts that the events should pass to be considered as selected

  float Nobs_mumu = 0.;
  //////////////////////

  //////////////////////
  // Initialisation
  //////////////////////
  string xmlFileName;
  cout<<"argc "<<argc<<" "<<argv[0]<<endl;
  if (argc>1 ) xmlFileName = string(argv[1]);
  else xmlFileName = string ("../../config/RPVAnalysis.xml");
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets);	// now the list of datasets written in the xml file is known
  anaEL.LoadSSDiLeptonSelection (sel);	// now the parameters for the selection are given to the selection
  anaEL.LoadGeneralInfo (DataType, Luminosity, LumiError, verbosity);
  int flagOriginal=sel.GetFlagb();
  int methodOriginal=sel.GetMethodb();
  int systOriginal= sel.GetSystb();
  std::cout << " For btag : flag " << flagOriginal << ", method " << methodOriginal << ", syst " << systOriginal << std::endl;
  IPHCTree::NTEvent * event = 0;
  //Selection table
  SelectionTable selTable_allChannels (sel.GetCutList (), datasets, string ("*"));
  SelectionTable selTable_ee (sel.GetCutList (), datasets, string ("ee"));
  SelectionTable selTable_emu (sel.GetCutList (), datasets, string ("emu"));
  SelectionTable selTable_mumu (sel.GetCutList (), datasets, string ("mumu"));
  
  //Book keeping of standard histos
  bool doHistoManager = true;
  SSDiLepAnaHistoManager histoManager;
  if(doHistoManager){
  	histoManager.LoadDatasets (datasets);
  	histoManager.LoadSelectionSteps (sel.GetCutList ());
  	histoManager.LoadChannels (sel.GetChannelList ());
  	histoManager.CreateHistos ();
  }
  //////////////////////


  cout << "The verbosity mode is " << verbosity << endl;
  cout << "The luminosity is equal to " << Luminosity << endl;
  cout << "The DataType is ";
  switch (DataType) {
  case 0:
    cout << "MC" << endl;
    break;
  case 1:
    cout << "Data" << endl;
    break;
  case 2:
    cout << "Data & MC" << endl;
    break;
  default:
    cout << " unknown" << endl;
    break;
  }
  //////////////////////



  //////////////////////
  //LOOP OVER THE DATASETS
  //////////////////////
  if (verbosity > 0) {
    cout << "#########################" << endl;
    cout << " Loop over the datasets  " << endl;
    cout << "#########################" << endl;
  }

  for (unsigned int d = 0; d < datasets.size (); d++) 
  {

    if(verbosity>2) cout<<"RPVAnalysis> Dataset: "<<datasets[d].Name()<<endl;
    datasets[d].eventTree ()->SetBranchAddress ("NTEvent", &event);

    unsigned int nEvents = static_cast<unsigned int>(datasets[d].eventTree ()->GetEntries ());
    cout << "RPVAnalysis> NEvents = " << nEvents << endl;
    cout << "RPVAnalysis> NEvents to run over = "<<datasets[d].NofEvtsToRunOver()<<endl;


    //LOOP OVER THE EVENTS
    for (unsigned int ievt = 0; ievt < datasets[d].NofEvtsToRunOver(); ievt++)
//.........这里部分代码省略.........

int main(int argc, char* argv[]){
  
  //This line could be commented if you don't want display while running, by example using screen.
  //TApplication theApp("App",&argc,argv);
  
  //---------------------------------------//
  // Global variables: could be give as argument later
  //---------------------------------------//
  
  int nwnodes = 10; //8 to 10 is the optimal
  string macroName = "ProofSelectorTTBarTagAndProbe.C+"; //"+" should be put at the end to use ACLIC complication - This macro should inherit from TSelector 
  string xmlFileName = string("/opt/sbg/data/data1/cms/jandrea/TopIPHC_2012_01_16/CMSSW_4_2_8_patch7/src/MiniTreeAnalysis/NTupleAnalysis/config/TTBarTagAndProbe.xml");  
  string outputFileName = "proof.root";
  
  //---------------------------------------//
  //	Decleration of TProof
  //---------------------------------------//
  
  //to be done before colling TProof
  system("../GeneralExamples/./clean_proof.sh ; echo 'Wait few seconds ... ' ; sleep 6");
  system("rm -r $HOME/.proof");
  
  TProof *proof = TProof::Open("");
  proof->SetParallel(nwnodes);
  //you should not have any package yet
  proof->ShowPackages();
  //proof->ClearPackages();
  //Loading package related to NTupleAnalysis
  cout<<" ## Upload package NTAna.par: ";
  proof->UploadPackage("../NTAna.par");
  cout<<" DONE [don't worry with symlink error - do rm NTAna if you change NTAna.par in the meanwhile !] "<<endl;
  proof->EnablePackage("NTAna");
  //Adding histograms for feedback: must exist in the TSelector !
  //proof->AddFeedback("fHist"); //give the "name" of the histogram and not the name of the variable TH1F* (could be the same !)
  
  //This line is required to display histograms durint the process
  TDrawFeedback fb(proof);
  
  
  //---------------------------------------//
  // Xml Loading  & Dataset registration
  //---------------------------------------//
  
  vector < Dataset > datasets;
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets); // now the list of datasets written in the xml file is known
  
  cout<<" #------------------------------------# "<<endl;
  cout<<" PROOF DATASETS SUMMARY [normaly 0]"<<endl;
  proof->ShowDataSets();
  cout<<" #------------------------------------# "<<endl;
  cout<<" # Registring dataset ... "<<endl;
  cout<<" Don't be worry with the checksum error message [at least I'm not ;-) ]"<<endl;
  cout<<" #------------------------------------# "<<endl;
  //Create datasets in proof format
  TFileCollection** fileCollec = new TFileCollection*[datasets.size()];
  for(unsigned int i=0;i<datasets.size();i++){
    fileCollec[i]  = new TFileCollection(datasets[i].Name().c_str(),"");
    for(unsigned int j=0;j<datasets[i].Filenames().size();j++){
      fileCollec[i]->Add(datasets[i].Filenames()[j].c_str());
    }
    //register dataset in proof
    proof->RegisterDataSet(datasets[i].Name().c_str(),fileCollec[i]);
    proof->VerifyDataSet(datasets[i].Name().c_str());
  }
  
  //summarize the list of datasets
  cout<<" #------------------------------------# "<<endl;
  cout<<" PROOF DATASETS SUMMARY"<<endl;
  proof->ShowDataSets();
  cout<<" #------------------------------------# "<<endl;
  
  //---------------------------------------//
  // 	Loading of the xml file
  //---------------------------------------//
  //Possibility to give float ... ex:
  //Double_t f = 3.14;
  //proof->SetParameter("IN_FLOAT",f);
  
  //---------------------------------------//
  // 	Processing of the datasets
  //---------------------------------------//
  
  string outputFileNameModif = outputFileName.substr(0,outputFileName.size()-5);
  //string MergingCommand = "hadd "+outputFileNameModif+"_merged.root "+outputFileNameModif+"_*.root  ";
  
  
  for(unsigned int i=0;i<datasets.size();i++){
    proof->AddInput(new TNamed("PROOF_DATASETNAME", datasets[i].Name()));
    //---------------------------------------//
    // 	Loading of the xml file
    //---------------------------------------//
    //Possibility to give float ... ex:
    //Double_t f = 3.14;
    //proof->SetParameter("IN_FLOAT",f);
    
    proof->AddInput(new TNamed("PROOF_XMLFILENAME", xmlFileName));
    proof->AddInput(new TNamed("PROOF_OUTPUTFILE", outputFileName));
    
    cout<<"#------------------------------------# "<<endl;
//.........这里部分代码省略.........

int main (int argc, char *argv[])
{

 
  cout << "#########################" << endl;
  cout << "Beginning of the program" << endl;
  cout << "#########################" << endl;

  //////////////////////
  //Global variables
  //////////////////////
  vector < Dataset > datasets;
  DiLeptonSelection sel;
  float Lumi = 0;
  float LumiError = 0.;
  // 0: MC - 1: Data - 2 Data & MC
  int DataType = 0;
  int verbosity = 2;
  int doPseudoExp = 0;
  int NofPseudoExp = 0;
  //////////////////////

  //////////////////////
  // Initialisation
  //////////////////////
  string xmlFileName;
  cout<<"argc "<<argc<<" "<<argv[0]<<endl;
  if (argc>1 ) xmlFileName = string(argv[1]);
  else xmlFileName = string ("../../config/PLRMeas.xml");

  AnalysisEnvironmentLoader anaEL (xmlFileName);
  if(verbosity>1) cout<<" - Loading datasets ..."<<endl;
  anaEL.LoadSamples (datasets);	// now the list of datasets written in the xml file is known
  if(verbosity>1) cout<<" - Loading general info ..."<<endl;
  anaEL.LoadGeneralInfo (DataType, Lumi, LumiError, verbosity);

  if(verbosity>1) cout<<" - Initializing PLR ..."<<endl;
  PLRMeasurement plr;
  cout<<"a"<<endl;
  plr.LoadDatasets(datasets);
  cout<<"a"<<endl;
  plr.SetLumi(Lumi,LumiError);
  cout<<"a"<<endl;
  //Load additionnal spectific info
  if(verbosity>1) cout<<" - Loading PLRInformation ..."<<endl;
  anaEL.LoadPLRInformation(plr,doPseudoExp, NofPseudoExp);
  if(verbosity>1) if(doPseudoExp) cout<<"Wil perform "<<NofPseudoExp<<" pseudo experiments "<<endl;
  if(verbosity>1) cout<<" - Loading bkg info ..."<<endl;
  anaEL.LoadBgkInformation(plr);
  //Do this after loading the other info (datasets, backgrounds)
  if(verbosity>1) cout<<" - Loading histograms ..."<<endl;
  plr.LoadHistos();
  plr.AddBkgDeterminationParameter();
  if(verbosity>1) cout<<" - Loading Signal Systematics ..."<<endl;
  anaEL.LoadSystematicsEffect(plr);

 
  if(verbosity>1) cout<<" - Run the PLR ..."<<endl;
  plr.RunLikelihood(doPseudoExp,NofPseudoExp,false,true,true);
  //plr.RunLikelihood(false,1,false,true,true);
  //plr.RunLikelihood(true,10000,false,true,true);
  //plr.RunLikelihood(true,1,false,true,true);
 
  plr.TimingReport();
  
  if(verbosity>1) cout<<" - Write output tex file ..."<<endl;
  string ofoutname = string(argv[2])+".tex";
  ofstream ofout(ofoutname.c_str());
  plr.MinimizationTable(ofout);
  plr.UncertaintiesTable(ofout);
  plr.NumberOfEventsTable(ofout);

  if(verbosity>1) cout<<" - Write output root file ..."<<endl;
  TFile* fout = new TFile(argv[2],"RECREATE");
  //TFile* fout = new TFile("PLR.root","RECREATE");
  plr.Write(fout);
  fout->Close();
  delete fout;

 cout << "#########################" << endl;
  cout << "  End of the program" << endl;
  cout << "#########################" << endl;

}

int main (int argc, char *argv[])
{
  cout << endl;
  cout << "   ,-----------------------," << endl;
  cout << "   |   Starting analysis   |" << endl;
  cout << "   `-----------------------`" << endl;
  cout << endl;

  // ##################################
  // #	Some parameters for the user  #
  // ##################################
 
  bool doHistoManager = false;		// False will not produce histos
  bool doLatexOutput = true;		// False will not produce latex table
  bool doPileUpReweigthing = false;	// False will not use pileup reweighting

  // (can be overwritten by argv[1])
  string defaultConfiguration("../../config/TTbarMETAnalysis.xml");
  // (can be overwritten by argv[2])
  string defaultLatexTableName("CutflowTable");
  // (can be overwritten by argv[3])
  string defaultRootOutputName("TTbarMETanalysis.root");

  // ############################
  // #	Initializing variables  #
  // ############################
  
  INFO1_MSG << "Initializing variables..." << endl;
 
  vector < Dataset > datasets;
  TTbarMetSelection sel;
  
  float Luminosity = 0;
  float LumiError = 0.;
  int DataType = 0;				// DataType : 0: MC - 1: Data - 2 Data & MC
  int verbosity = -1;


  reweight::LumiReWeighting *LumiWeights;
  IPHCTree::NTEvent * event = 0;


  std::vector<std::string> dileptonCol;
  dileptonCol.push_back("l+jets");
  dileptonCol.push_back("ll");
  dileptonCol.push_back("llHadrTau");
  dileptonCol.push_back("llHadrTau3+");

  std::vector<std::string> dileptonRow;
  dileptonRow.push_back("baseline");
  dileptonRow.push_back("MET");
  dileptonRow.push_back("stdVeto");
  dileptonRow.push_back("stdVetoEff");
  dileptonRow.push_back("vetoTest0");
  dileptonRow.push_back("vetoTest0Eff");
  dileptonRow.push_back("vetoTest1");
  dileptonRow.push_back("vetoTest1Eff");
  dileptonRow.push_back("vetoTest2");
  dileptonRow.push_back("vetoTest2Eff");
  dileptonRow.push_back("vetoTest3");
  dileptonRow.push_back("vetoTest3Eff");

  // Reading parameters from argv
  // 	-> configuration file
  string xmlFileName;
  if (argc > 1) xmlFileName = string(argv[1]);
  else          xmlFileName = defaultConfiguration;
  // 	-> root output name
  string rootOutputName;
  if (argc > 2) rootOutputName = string(argv[2]);
  else 			rootOutputName = defaultRootOutputName;
  // 	-> latex table name
  string latexTableName;
  if (argc > 3) latexTableName = string(argv[3]);
  else latexTableName = defaultLatexTableName;
  
  // #############################
  // # 	 Loading configuration   #
  // #############################

  cout << endl;
  INFO1_MSG << "Loading configuration..." << endl;
  cout << "        (config : " << xmlFileName << ")" << endl;
  
 
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets);	// now the list of datasets written in the xml file is known
  anaEL.LoadSelection (sel);	// now the parameters for the selection are given to the selection // no specific TTbarMET parameters
  anaEL.LoadGeneralInfo (DataType, Luminosity, LumiError, verbosity);

  // #########################################
  // # 	 Creating tables and histomanagers   #
  // #########################################

  cout << endl;
  INFO1_MSG << "Creating tables and histomanagers..." << endl;

  // Stuff
  
  TFile dileptonTreeFile("dileptonTree.root","RECREATE","Tree for Dilepton background study");
//.........这里部分代码省略.........

int main ()
{
  cout<<"#########################"<<endl;
  cout<<"Beginning of the program"<<endl;
  cout<<"#########################"<<endl;
  
  //////////////////////
  //Global variables
  //////////////////////
  vector < Dataset > datasets;
  DiLeptonSelection sel; 
  float Luminosity = 0;
  // 0: MC - 1: Data - 2 Data & MC
  int DataType = 0; 
  int verbosity = -1;

  //////////////////////

  //////////////////////
  // Initialisation
  //////////////////////
  string xmlFileName = string ("../../config/MyVar.xml");
  AnalysisEnvironmentLoader anaEL (xmlFileName);
  anaEL.LoadSamples (datasets); // now the list of datasets written in the xml file is known
  anaEL.LoadDiLeptonSelection (sel); // now the parameters for the selection are given to the selection
  anaEL.LoadGeneralInfo(DataType, Luminosity, verbosity );
  anaEL.LoadWeight (sel); // now the parameters for SFBweight are initialized (for b-tag!)
  TopTree::NTEvent * event = 0;

  //dataset weights
  double Dweight[101]; 
  for(int k0=0; k0<5; ++k0) {
    for(int k1=0; k1<101; ++k1) {
      Dweight[k1] = 0.;
    }
  } 
  vector<string> CutName;
  vector<string> VecChannelName;
  


// Here define the studied channel (ee/mumu/emu)
  string ChannelName  = "mumu";  // "mumu", "ee", "emu"


  int ITypeMC	  = -1;
  int ICut	  = -1;  
  int IChannel    = -1;  
  if (  ChannelName  == "mumu" )  IChannel = 0;
  if (  ChannelName  == "ee" )    IChannel = 1;
  if (  ChannelName  == "emu" )   IChannel = 2;

  
  
  int nbins = 200;
  float minx = 0.;
  float maxx = 350;
  
  
  cout<<"The verbosity mode is "<<verbosity <<endl;
  cout<<"The luminosity is equal to "<< Luminosity<<endl;
  cout<<"The DataType is ";
  switch(DataType){
  case 0: 
    cout<<"MC"<<endl; 
    break;
  case 1:
    cout<<"Data"<<endl; 
    break;
  case 2:
    cout<<"Data & MC"<<endl; 
    break;
  default:
    cout<<" unknown"<<endl; 
    break;
  }
  //////////////////////
  
    HistoManager MyhistoManager;
    MyhistoManager.LoadDatasets(datasets);   
     

//     std::vector<TH1F> PATElecPt;
//     std::vector<TH1F> PATElecRelIso;
//     std::vector<TH1F> PATMuonPt;
//     std::vector<TH1F> PATMuonRelIso;
    
    std::vector<TH1F> PFElecPt;
    std::vector<TH1F> PFElecNeutralHadIso;
    std::vector<TH1F> PFElecChargedHadIso;
    std::vector<TH1F> PFElecPhotonIso;
    std::vector<TH1F> PFElecTrackIso;
    std::vector<TH1F> PFElecRelIso;
    std::vector<TH1F> PFElecEcalDriven;
    std::vector<TH1F> PFElecGsfEl;
    std::vector<TH1F> PFElecD0;
    std::vector<TH1F> PFElecId90;
    std::vector<TH1F> PFElecConv;
    std::vector<TH1F> PFElecRelIsoRhoCorrected;
    std::vector<TH1F> PFElecRelIsoNeutralRhoCorrected;
//.........这里部分代码省略.........

int main ()
{
  cout<<"#########################"<<endl;
  cout<<"Beginning of the program"<<endl;
  cout<<"#########################"<<endl;
  
  //////////////////////
  // Initialisation
  //////////////////////
  float Luminosity = 0;
  float LumiError = 0;
  string PUWeightFileName;
  int DataType =0;

  string xmlFileName = string ("../../config/MatrixMethod_MC.xml");
  AnalysisEnvironmentLoader anaEl (xmlFileName);
  vector < Dataset > datasets;
  anaEl.LoadSamples(datasets);
  int verbosity = -1;
  anaEl.LoadGeneralInfo(DataType, Luminosity, LumiError, PUWeightFileName, verbosity);

  ////////////////////////////////
  // Matrix Method declaration
  ///////////////////////////////
  float looseIsoMM = 0.8; float tightIsoMMe = 0.17; float tightIsoMMmu = 0.2;
  unsigned int nbinsMM = 11; float lowEdgeMM = -0.5; float highEdgeMM = 10.5;
  MMEstimation MMestEE(datasets, looseIsoMM, tightIsoMMe, nbinsMM, lowEdgeMM, highEdgeMM, "EE");
  MMEstimation MMestMuMu(datasets, looseIsoMM, tightIsoMMe, nbinsMM, lowEdgeMM, highEdgeMM, "MuMu");

  //////////////////////

  //////////////////////////////////////////
  // Running the Matrix Method in fast way
  //////////////////////////////////////////
  //MMestEE.ReadMMFileForPullTest("MatrixMethod_OutPut_EE_MC_forNoWeight.root", 0.95, 0.2);
  MMestMuMu.ReadMMFileForPullTest("MatrixMethod_OutPut_MuMu_MC_NoWeight.root", 0.95, 0.2);

  /*vector<struct MMEpsilons> valMMEpsilons_EE;
  for(unsigned int bin_index = 0; bin_index < nbinsMM; bin_index++){
    struct MMEpsilons valMMEpsilonsTmp_EE;
    valMMEpsilonsTmp_EE.EpsilonSignal = 0.95;
    valMMEpsilonsTmp_EE.EpsilonSignalErr = 0.1;
    valMMEpsilonsTmp_EE.EpsilonFake = 0.2;
    valMMEpsilonsTmp_EE.EpsilonFakeErr = 0.1;
    valMMEpsilons_EE.push_back(valMMEpsilonsTmp_EE);
  }*/

  vector<struct MMEpsilons> valMMEpsilons_MuMu;
  for(unsigned int bin_index = 0; bin_index < nbinsMM; bin_index++){
    struct MMEpsilons valMMEpsilonsTmp_MuMu;
    valMMEpsilonsTmp_MuMu.EpsilonSignal = 0.95;
    valMMEpsilonsTmp_MuMu.EpsilonSignalErr = 0.1;
    valMMEpsilonsTmp_MuMu.EpsilonFake = 0.2;
    valMMEpsilonsTmp_MuMu.EpsilonFakeErr = 0.1;
    valMMEpsilons_MuMu.push_back(valMMEpsilonsTmp_MuMu);
  }


  float weight = datasets[0].NormFactor()*Luminosity;

  unsigned int NbIterations = 100000;
  bool doStatistical = true; bool doSystematic = false;
  
  /*MMestEE.RunTheMatrixMethodForEpsilonsTest(valMMEpsilons_EE, NbIterations, doStatistical, doSystematic, weight);
  MMestEE.PrintMMEstimated();
  MMestEE.WriteMMFileFastForPullTest("MatrixMethod_OutPut_EE_MC_forNoWeight.root","MatrixMethod_OutPut_EE_MC_Fast_EpsilonsTest.root", 0.95, 0.2);*/


  MMestMuMu.RunTheMatrixMethodForEpsilonsTest(valMMEpsilons_MuMu, NbIterations, doStatistical, doSystematic, weight);
  MMestMuMu.PrintMMEstimated();
  MMestMuMu.WriteMMFileFastForPullTest("MatrixMethod_OutPut_MuMu_MC_NoWeight.root","MatrixMethod_OutPut_MuMu_MC_Fast_EpsilonsTest.root", 0.95, 0.2);


  cout<<"#########################"<<endl;
  cout<<"    End of the program   "<<endl;
  cout<<"#########################"<<endl;

  return (0);
}

void MatrixMethod (string xmlFileName, bool mc=true, string channel="") // MuMu or EE for data, do both for MC
{
  cout<<"#########################"<<endl;
  cout<<"Beginning of the program"<<endl;
  cout<<"#########################"<<endl;
  
  //////////////////////
  // Initialisation
  //////////////////////
  float Luminosity = 0;
  float LumiError = 0;
  string PUWeightFileName;
  int DataType = 0; 

  AnalysisEnvironmentLoader anaEL (xmlFileName);
  vector < Dataset > datasets;
  anaEL.LoadSamples (datasets); // now the list of datasets written in the xml file is known
  int verbosity = -1;

  // Matrix Method
  int MMselStepCut = 3;
  float looseIsoMM = 0.8; float tightIsoMMmu = 0.2; float tightIsoMMe = 0.17;
  unsigned int nbinsMM = 11; float lowEdgeMM = -0.5; float highEdgeMM = 10.5; // Histo nb events en fonction #jets
  MMEstimation MMestEE(datasets, looseIsoMM, tightIsoMMe, nbinsMM, lowEdgeMM, highEdgeMM, "EE"); 
  MMEstimation MMestMuMu(datasets,looseIsoMM, tightIsoMMmu, nbinsMM, lowEdgeMM, highEdgeMM, "MuMu"); 

  SSDiLeptonSelection sel;
  anaEL.LoadSSDiLeptonSelection (sel); // now the parameters for the selection are given to the selection
  anaEL.LoadGeneralInfo(DataType, Luminosity, LumiError, PUWeightFileName, verbosity );
  //Load for PU:
  if(mc) sel.GeneratePUWeight(PUWeightFileName);

  IPHCTree::NTEvent * event = 0;

  /*PUWeighting  thePUReweighter;
  TFile* file1  = new TFile(PUWeightFileName.c_str(),"READ"); 
  TH1D *  hPUData = 0;
  hPUData         = (TH1D*)file1->Get("pileup");
  TH1F *  hPUMC   = new TH1F("pileup_MC", "pileup_MC", hPUData->GetXaxis()->GetNbins(), hPUData->GetXaxis()->GetXmin(), hPUData->GetXaxis()->GetXmax() );
  TFile* file2  = new TFile( "../data/CrossSection_pileup.root" ,"READ");
  hPUMC           = (TH1F*)file2->Get("pileup_TTbarSig");
  // histo in data, histo in Mc, use out-of-time pu in the reweighting
  cout << "get MC histo  " << endl;
  thePUReweighter.setPUHisto( hPUData, hPUMC);
  cout << "set MC histo in thePUReweighter " << endl;
  thePUReweighter.setUseOutOfTimePU(false); // set to true to use out-of-time PU
  */
  
  //////////////////////
  //LOOP OVER THE DATASETS
  //////////////////////
  cout<<"#########################"<<endl;
  cout<<" Loop over the datasets  "<<endl;
  cout<<"#########################"<<endl;

  for (unsigned int d = 0; d < datasets.size (); d++) {
   datasets[d].eventTree ()->SetBranchAddress ("NTEvent",&event);
   cout << "dataset = " << datasets[d].Name() << endl;
   unsigned int nEvents = (int) (datasets[d].eventTree ()->GetEntries ());
    cout << "NEvents = " << nEvents << endl;
    float weight_init;
    weight_init = datasets[d].NormFactor()*Luminosity;
    cout << "weight_init = " << weight_init << endl;
    //LOOP OVER THE EVENTS
    for (unsigned int ievt = 0; ievt < nEvents; ievt++) {
      datasets[d].eventTree ()->GetEntry (ievt);
      IPHCTree::NTTransient::InitializeAfterReading(event); // Important line to read new format files
      //Load event for the selection
      sel.LoadEvent(event);
      if(ievt%10000 == 0) cout << "number of processed events " << ievt << endl;
      
      float weight = 1;
      
      if(mc){
      	
	//Manage DY samples to avoid overlaps
	double dileptInvMass = 0;
	if( (event->mc.zAndDecays).size() > 0){
          TLorentzVector dilept = (event->mc.zAndDecays)[0].p4_Lep1_gen + (event->mc.zAndDecays)[0].p4_Lep2_gen;
          dileptInvMass = dilept.M();
	}
	if(datasets[d].Name()=="Zjets" && dileptInvMass < 50) continue;
	if(datasets[d].Name()=="DYToMuMu_M-20"	   && (dileptInvMass > 50 || dileptInvMass < 20) ) continue;
	if(datasets[d].Name()=="DYToEE_M-20"        && (dileptInvMass > 50 || dileptInvMass < 20) ) continue;
	if(datasets[d].Name()=="DYToTauTau_M-20"    && (dileptInvMass > 50 || dileptInvMass < 20) ) continue;
	if(datasets[d].Name()=="DYToMuMu_M-10To20"   &&  dileptInvMass > 20) continue;
	if(datasets[d].Name()=="DYToEE_M-10To20"    &&  dileptInvMass > 20) continue;
	if(datasets[d].Name()=="DYToTauTau_M-10To20" &&  dileptInvMass > 20) continue;


        weight = weight_init;
        //float weight = weight_init*sel.GetPUWeight();
        /*if(thePUReweighter.getUseOutOfTimePU()){
           weight = weight_init*thePUReweighter.weight(event->pileup.intime_npu, event->general.runNb);
         }else{
	   weight = weight_init*thePUReweighter.weight(event->pileup.intime_npu);
       }*/
      }
      
      if(mc || (!mc && channel=="EE")) MMestEE.CountNSel(sel, datasets[d], "ee_ss", weight, MMselStepCut, &(event->mc));
//.........这里部分代码省略.........