C++ TChain类代码示例

//.........这里部分代码省略.........
    listTrees = list_Spikes_Run2011B_Newkill();
  }
  /*
  if(data=="Spikes_Run2011A_HLTEG") {
    listTrees = list_Spikes_Run2011A_HLTEG();
  }
  else if(data=="Spikes_Run2011B_HLTEG") {
    listTrees = list_Spikes_Run2011B_HLTEG();
  }
  else if(data=="Spikes_Run2011B_highPU_HLTEG") {
    listTrees = list_Spikes_highPU_HLTEG();
  }
  else if(data=="Spikes_2011A_HLTEG5") {
    listTrees = list_Spikes_2011A_HLTEG5();
  }
  else if(data=="Spikes_2011B_HLTEG5") {
    listTrees = list_Spikes_2011B_PRV1_HLTEG5();
  }
  else if(data=="Spikes_2011AB_HLTEG5") {
    listTrees = list_Spikes_2011AB();
  }
  else if(data=="Spikes_Run2011A") {
    listTrees = list_Spikes_Run2011A();
  }
  else if(data=="Spikes2011A_PRV6") {
    listTrees = list_Spikes_2011A_PromptRecoV6();
  }
  else if(data=="Spikes2011A") {
    listTrees = list_Spikes_2011A();
  }
  else if(data=="Spikes2011A_PromptRecoV4") {
    listTrees = list_Spikes_2011A_PromptRecoV4();
  }
  else if(data=="Spikes2011A_1fb") {
    listTrees = list_Spikes_2011A_1fb();
  }
  else if(data=="Spikes2011A_1fb_try2") {
    listTrees = list_Spikes_2011A_1fb_try2();
  }
  */
  /*
  else if(data=="Spikes2011A_HLTEG5") {
    listTrees = list_Spikes_2011A_HLTEG5();
  }
  */
  else {
    cout << "problem of dataset tag" << endl;
    return 0;
  }
  
  if(listTrees.size()==0) {
    cout << "tree list empty" << endl;
    return 0;
  }
  
  if(iHalf<0 || nHalf<1 || iHalf>=nHalf) {
    cout << "problem with iHalf and nHalf" << endl;
    return 0;
  }

  if(nHalf>listTrees.size()) {
    cout << "requested more halfes than available trees ; set nHalf to listTrees.size()=" << listTrees.size() << endl;
    nHalf = listTrees.size();
  }

  int nTrees, nStart, nEnd;
  nTrees = listTrees.size() / nHalf;
  if(iHalf < nHalf-1) {
    nStart = iHalf*nTrees;
    nEnd = (iHalf + 1)*nTrees;
  }
  else if(iHalf==nHalf-1) {
    nStart = iHalf*nTrees;
    nEnd = listTrees.size();
  }

  if(debug) cout << "gonna add the trees to the chain" << endl;

  TChain * myChain = new TChain ("produceNtuple/eIDSimpleTree");
  for(int i=nStart ; i<nEnd ; i++)
    myChain->Add(listTrees[i]);

  if(debug) cout << "added" << endl;

  cout << "tree list size : " << listTrees.size() << endl
	 << "nStart=" << nStart << " nEnd=" << nEnd << " nTrees=" << nTrees << endl; 
  //for(int i=nStart ; i<nEnd ; i++)
  //outlog << listTrees[i] << endl;

  if(nEntries==-999) {
    cout << "no tree processing requested ; exiting" << endl;
    return 0;
  }

  // Process the tree
  if(debug) cout << "process the tree" << endl;
  TreeToHistos(myChain, data_sum, dirOut, outlog, nEntries, halftag, spikeCut, debug, graph_style,sev_user,hlt_sel);

  return 1;
}

void runMakePhotonBabies(char* prefix , bool isData = true, float kFactor = 1.){

  TChain* ch = new TChain("Events");

  //-----------------------------------------------------------------------------------

  if( strcmp( prefix , "data_53X_2012A" ) == 0 ){    
    //pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012A-13Jul2012-v1_AOD/V05-03-18_slim/merged_ntuple_999999_9_1_skim.root");
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012A-13Jul2012-v1_AOD/V05-03-18_slim/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012A-recover-06Aug2012-v1_AOD/V05-03-18_slim/merged*root");
  }

  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "data_53X_2012B" ) == 0 ){    
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012B-13Jul2012-v1_AOD/V05-03-18_slim/merged*root");
  }

  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "data_53X_2012C" ) == 0 ){    
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012C-24Aug2012-v1_AOD/V05-03-18_slim/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012C-PromptReco-v2_AOD/V05-03-18_slim/merged*root");
  }

  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "data_53X_2012D" ) == 0 ){    
    pickSkimIfExists(ch,"PhotonTriggerSkim_slim/DoubleElectron_Run2012D-PromptReco-v1_AOD/V05-03-18_slim/merged*root");
  }

  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "Photon" ) == 0 ){
    pickSkimIfExists(ch,"/hadoop/cms/store/user/cwelke/CMSSW_5_2_3_patch4_V05-02-27/Photon_Run2012A-PromptReco-v1_AOD/unmerged/store*root");
    pickSkimIfExists(ch,"/hadoop/cms/store/user/cwelke/CMSSW_5_2_3_patch4_V05-02-27/SinglePhoton_Run2012B-PromptReco-v1_AOD/unmerged/store*root");
  }

  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "DoubleElectron" ) == 0 ){
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012A-13Jul2012-v1_AOD/V05-03-13/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012A-recover-06Aug2012-v1_AOD/V05-03-13/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012B-13Jul2012-v1_AOD/V05-03-13/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012C-PromptReco-v1_AOD/V05-03-13/merged*root");
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012C-PromptReco-v2_AOD/V05-03-13/merged*root");
  }
  
  //-----------------------------------------------------------------------------------

  else if( strcmp( prefix , "DoubleElectron_2012Cv2" ) == 0 ){
    pickSkimIfExists(ch,"PhotonTriggerSkim/DoubleElectron_Run2012C-PromptReco-v2_AOD/V05-03-13/merged*root");
  }
  
  //-----------------------------------------------------------------------------------

  else{
    cout << "ERROR: cannot find sample " << prefix << endl;
    exit(0);
  }

  //-----------------------------------------------------------------------------------
    
  bool calculateTCMET = false;  //recalculate tcmet on-the-fly?
  
  cout << endl;
  cout << "Checking for corrupt files" << endl;
  cout << "Entries " << ch->GetEntries() << endl;
  ch->Draw("evt_run");
  cout << endl;
  
  makePhotonBabies* myLooper = new makePhotonBabies();
  
  cout << "Running on sample " << prefix << endl;
  myLooper->ScanChain(ch, prefix, isData, calculateTCMET, -1 ,kFactor);
  
}

//
// Fucntions
//
void skimThisBaby(TString inPath, TString inFileName, TString inTreeName, TString outPath){

  // Talk to user
  cout << "  Skimming: "<< endl;
  cout << "     " << inPath+inFileName <<endl;
  
  // Load input TChain
  TChain *ch = new TChain(inTreeName);
  TString inFile = inPath;
  inFile += inFileName;
  ch->Add(inFile);
  Long64_t nentries = ch->GetEntries();
  TTreeCache::SetLearnEntries(10);
  ch->SetCacheSize(128*1024*1024);

  // Initialize Branches
  babyAnalyzer.Init(ch->GetTree());
  babyAnalyzer.LoadAllBranches();
  
  // Setup output file name and path
  TString outFileName = outPath;
  outFileName += inFileName;
  outFileName.ReplaceAll("*", "");
  outFileName.ReplaceAll(".root", "_skimmed.root");

  cout << "  Output will be written to: " << endl;
  cout << "      "  << outFileName << endl << endl;

  // Open outputfile and Clone input TTree
  TFile *newfile = new TFile(outFileName, "recreate");
  TTree *newtree  = (TTree*)ch->GetTree()->CloneTree(0);
  if(!newtree) cout << "BAD TTREE CLONE" << endl;

  TH1D *newCounter=NULL;
  
  // Get nEntries
  unsigned int nEventsTotal = 0;
  unsigned int nEventsChain = ch->GetEntries();
  
  // Grab list of files
  TObjArray *listOfFiles = ch->GetListOfFiles();
  TIter fileIter(listOfFiles);
  TFile *currentFile = 0;

  // File Loop
  int iFile=0;
  while ( (currentFile = (TFile*)fileIter.Next()) ) {

    // Get File Content
    TFile *file = new TFile( currentFile->GetTitle() );
    TTree *tree = (TTree*)file->Get("t");
    TTreeCache::SetLearnEntries(10);
    tree->SetCacheSize(128*1024*1024);
    babyAnalyzer.Init(tree);

    if(iFile==0){
      TH1D *temp = (TH1D*)file->Get("h_counter");
      newCounter = (TH1D*)temp->Clone("h_counter");
      newCounter->SetDirectory(newfile);
    }
    //else{
    //TH1D *temp = (TH1D*)file->Get("h_counter");
    //newCounter->Add(temp);
    //}

    // Loop over Events in current file 
    if( nEventsTotal >= nEventsChain ) continue;
    unsigned int nEventsTree = tree->GetEntriesFast();
    for( unsigned int event = 0; event < nEventsTree; ++event) {
      
      // Progress
      stop_1l_babyAnalyzer::progress( nEventsTotal, nEventsChain );
    
      // Load Tree
      tree->LoadTree(event);
      babyAnalyzer.GetEntry(event);
      ++nEventsTotal;

      // Selection
      if(nvetoleps()<1) continue;
      if(nvetoleps()<2 && PassTrackVeto_v3() && PassTauVeto()) continue;
      if(ngoodjets()<2) continue;
      if(mt_met_lep()<150.0) continue;
      if(pfmet()<200.0) continue;

      // Turn on all branches on input
      babyAnalyzer.LoadAllBranches();      
      
      // Fill output tree
      newtree->Fill();
    
    } // end loop over entries
   
    iFile++;
  } // end loop over files in TChain

  // Clean up
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  
  RooDataSet* Background= new RooDataSet("BackgroundData","BackgroundData",RooArgSet(Mass,GammaCL,BDTResponse),Import(*MinimalTreeA),Cut("gamma_CL>0.1"));

  std::cout<<"Data file A loaded into Roo DataSet"<<std::endl;
  
  RooDataSet* BackgroundB= new RooDataSet("BackgroundDataB","BackgroundDataB",RooArgSet(Mass,GammaCL,BDTResponse),Import(*MinimalTreeB),Cut("gamma_CL>0.1"));

  std::cout<<"Data File B loaded into RooDataSet"<<std::endl;
  
  Background->append(*BackgroundB);

  RooPlot* massFrame =Mass.frame(Title("Data Import Check"),Bins(50));
  Background->plotOn(massFrame);

  TCanvas C;
  massFrame->Draw();
  C.SaveAs("DataImportCheck.eps");

  //================================== Make MC TChain=======================================

  TFile* TFMCA= new TFile((MCResponseA.GetFP()).data());
  TTree* MCATree=(TTree*)TFMCA->Get("DecayTree");
  if(!MCATree){
    std::cout<<"MC A Tree Null"<<std::endl;
    return 1;
  }
  TFile* TFMCB= new TFile((MCResponseB.GetFP()).data());
  TTree* MCBTree=(TTree*)TFMCB->Get("DecayTree");
  if(!MCBTree){
    std::cout<<"MC B Tree Null"<<std::endl;
    return 1;
  }

  TChain* MCChain = new TChain("DecayTree");
  MCChain->Add((MCResponseA.GetFP()).data());
  MCChain->Add((MCResponseB.GetFP()).data());

  Long64_t N= MCChain->GetEntries();

  std::cout<<"Entries In Chain = "<<N<<std::endl;

  //=============================== Loop over several cuts===============

  LoopTimer LP(0.10);
  //PunziMVACut* TestCut= new PunziMVACut(MCChain,Background,0.2,5);
  std::vector<PunziMVACut*> Cuts;
  RooDataSet* Conductor=Background;
  int nCuts=15;
  int nFineCuts=25;
  
  FillableGraph* FG= new FillableGraph("Punzi FoM Optimisation;BDT Cut;Punzi FoM");
  FillableGraph* FGThree = new FillableGraph("Punzi FoM Optimisation Sigma=3;BDT Cut; Punzi FoM");
  FillableGraph* FGFour = new FillableGraph("Punzi FoM Optimisation Sigma=4;BDT Cut; Punzi FoM");

  FillableGraph* SignalEfficiency= new FillableGraph("Signal Efficiency;BDT Cut;Signal Efficiency");

  std::vector<double> CutVals;
  for(int i=0;i<nCuts;++i){
    CutVals.push_back(0.0+(i*0.02));
  }
  for(int f=0;f<nFineCuts;++f){
    CutVals.push_back(0.30+(f*0.005));
  }
  //CutVals.push_back(0.20);                                                  
  std::sort(CutVals.begin(),CutVals.end());

void ProcessingTime(const char *inputFile)
{
  TChain *chain = new TChain("Delphes");
  chain->Add(inputFile);

  TH1F hist("time", "time", 50, 0, 0.01);
  Int_t i;

  TDirectory *currentDirectory = gDirectory;

  // Graphics style parameters to avoid grey background on figures
  gStyle->SetCanvasColor(kExRootBackgroundColor);
  gStyle->SetStatColor(kExRootBackgroundColor);
  //  gStyle->SetTitleColor(kExRootBackgroundColor);
  gStyle->SetPadColor(kExRootBackgroundColor);

  gStyle->SetPadTopMargin(0.10);
  gStyle->SetPadRightMargin(0.10);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadLeftMargin(0.15);

  gStyle->SetStatFont(kExRootFont);
  gStyle->SetStatFontSize(kExRootFontSize);

  gStyle->SetTitleFont(kExRootFont, "");
  gStyle->SetTitleFont(kExRootFont, "X");
  gStyle->SetTitleFont(kExRootFont, "Y");
  gStyle->SetTitleFont(kExRootFont, "Z");
  gStyle->SetTitleSize(kExRootFontSize, "");
  gStyle->SetTitleSize(kExRootFontSize, "X");
  gStyle->SetTitleSize(kExRootFontSize, "Y");
  gStyle->SetTitleSize(kExRootFontSize, "Z");

  gStyle->SetLabelFont(kExRootFont, "X");
  gStyle->SetLabelFont(kExRootFont, "Y");
  gStyle->SetLabelFont(kExRootFont, "Z");
  gStyle->SetLabelSize(kExRootFontSize, "X");
  gStyle->SetLabelSize(kExRootFontSize, "Y");
  gStyle->SetLabelSize(kExRootFontSize, "Z");

  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);

  gStyle->SetTextFont(kExRootFont);
  gStyle->SetTextSize(kExRootFontSize);

  gStyle->SetOptStat(111110);
  // gStyle->SetOptFit(101);

  canvas = static_cast<TCanvas*>(gROOT->FindObject("c1"));
  if(canvas)
  {
    canvas->Clear();
    canvas->UseCurrentStyle();
    canvas->SetWindowSize(800, 650);
  }
  else
  {
    canvas = new TCanvas("c1", "c1", 800, 650);
  }
  canvas->SetGrid();
  canvas->SetHighLightColor(kExRootBackgroundColor);

  currentDirectory->cd();

  for(i = 0; i < 9; ++i)
  {
    chain->Draw("Event.ProcTime >> time", TString::Format("Jet_size == %d", i+2));
    gr.SetPoint(i, i+2, hist.GetMean()*1000);
    grerr.SetPoint(i, i+2, hist.GetMean()*1000);
    grerr.SetPointError(i, 0, hist.GetRMS()*1000);
  }

  grerr.GetXaxis()->SetLimits(1.0, 11.0);
  grerr.GetXaxis()->SetTitleOffset(1.5);
  grerr.GetYaxis()->SetTitleOffset(1.75);
  grerr.GetXaxis()->SetTitle("jet multiplicity");
  grerr.GetYaxis()->SetTitle("processing time per event, ms");
  gr.SetMarkerStyle(kFullCircle);
  gr.SetMarkerColor(kBlack);
  gr.SetMarkerSize(1);
  gr.SetLineColor(kBlack);
  gr.SetLineWidth(2);
  grerr.SetFillStyle(1001);
  grerr.SetFillColor(17);
  grerr.Draw("A3");
  gr.Draw("P");

  comment.SetTextSize(kExRootFontSize);
  comment.SetTextFont(kExRootFont);
  comment.SetTextAlign(22);
  comment.SetFillColor(kExRootBackgroundColor);
  comment.SetBorderSize(0);
  comment.AddText("ttbar + jets events");
  comment.Draw();
}

bool dt_RunDrawTest(const char* from, Int_t mode = 0, Int_t verboseLevel = 0) {
  // This launch a test a TTree::Draw.
  // The mode currently available are:
  //    0: Do not load the shared library
  //    1: Load the shared library before opening the file
  //    2: Load the shared library after opening the file
  //    3: Simple TChain test with shared library
  //    4: Simple Friend test with shared library
  // The verboseLeve currently available:
  //    0: As silent as possible, only report errors and overall speed results.
  //    1: Output 0 + label for the start of each phase
  //    2: Output 1 + more details on the different phase being done
  //    3: Output 2 + stop at the first and draw a canvas showing the differences

//gDebug = 5;
   SetVerboseLevel(verboseLevel);

   if (mode == 1) {
      if (!TClassTable::GetDict("Event")) {
         gSystem->Load("Event_cxx");
     }     
      gHasLibrary = kTRUE;
   }

   TFile *hfile = 0;
   TTree *tree = 0;
   if (mode <3) {
      hfile = new TFile(from);
      tree = (TTree*)hfile->Get("T");
   }

   if (mode >= 2 && mode <= 4) {
      if (!TClassTable::GetDict("Event")) {
         gSystem->Load("Event_cxx");
      } else {
         cerr << "Since libEvent.so has already been loaded, mode 2 can not be tested!";
         cerr << endl;
      }
      gHasLibrary = kTRUE;
   }

   if (mode == 3) {
      // Test Chains.
      TChain * chain = new TChain("T");
      chain->Add(from);
      chain->Add(from);
      tree = chain;
   }

   if (mode == 4) {
      // Test friends.
      tree = new TTree("T","Base of friendship");
      tree->AddFriend("T",from);
   }

   TBranch *eb = tree->GetBranch("event");
   gBranchStyle = (int) eb->InheritsFrom(TBranchElement::Class());
   // cerr << "Branch style is " << gBranchStyle << endl;

   if (gQuietLevel<2) cout << "Generating histograms from TTree::Draw" << endl;
   TDirectory* where = GenerateDrawHist(tree,2,gQuietLevel);
 
   if (gQuietLevel<2) cout << "Comparing histograms" << endl;
   if (Compare(where)>0) {
     cout << "DrawTest: Comparison failed" << endl;
     return false;
   }
   DrawMarks();

   if (gQuietLevel<2) cout << "DrawTest: Comparison was successfull" << endl;
   if (hfile) delete hfile;
   else delete tree;
   gROOT->GetList()->Delete();

   return true;
}   

// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv) {

  // check number of inpt parameters
  if(argc < 2){
    cerr<<"Forgot to parse the cfg file --> exit "<<endl;
    return -1;
  }

  // Set Root style from global enviroment path                                                                                                                               
  string ROOTStyle;
  if(getenv ("ROOTStyle")!=NULL){
    ROOTStyle = getenv ("ROOTStyle");
    gROOT->ProcessLine((".x "+ROOTStyle+"/setTDRStyle.C").c_str());
  }

  gStyle->SetOptStat(0);
  gStyle->SetPadTopMargin(0.09);
  gStyle->SetPadLeftMargin(0.13);
  gStyle->SetErrorX(0.5);


  // parse config file parameter
  if (gConfigParser) return 1 ;
  gConfigParser = new ConfigParser();

  TString config ; 
  config.Form("%s",argv[1]);

  if(!(gConfigParser->init(config))){
      cout << ">>> parseConfigFile::Could not open configuration file " << config << endl;
      return -1;
  }

  // import base directory where samples are located and txt file with the directory name + other info
  string InputBaseDirectory  = gConfigParser -> readStringOption("Input::InputBaseDirectory");

  // import from cfg file the cross section value for this sample
  float CrossSection   = gConfigParser -> readFloatOption("Input::CrossSection");

  // total number of events
  int maxEventNumber   = gConfigParser -> readFloatOption("Input::EventsNumber");

  // treeName
  string treeName      = gConfigParser -> readStringOption("Input::TreeName");

  // take the cut list
  string InputCutList  = gConfigParser -> readStringOption("Input::InputCutList");

  // Read the cut file
  vector <cutContainer> CutList;
  if(ReadInputCutFile(InputCutList,CutList) <= 0){
    cerr<<" Empty Cut List File or not Exisisting --> Exit "<<endl; return -1;}

  // take the variable list to be plotted
  string InputVariableList  = gConfigParser -> readStringOption("Input::InputVariableList");
  vector<variableContainer> variableList;

  if(ReadInputVariableFile(InputVariableList,variableList) <= 0 ){
    cerr<<" Empty Variable List File or not Exisisting --> Exit "<<endl; return -1;}

  // take lumi and other parameters
  float lumi  =  gConfigParser -> readFloatOption("Option::Lumi"); // fb^(-1)
  lumi *= 1000. ;   // transform into pb^(-1)

  finalStateString    = gConfigParser -> readStringOption("Option::finalStateString");

  matchingCone        = gConfigParser -> readFloatOption("Option::matchingCone"); 
  minLeptonCleaningPt = gConfigParser -> readFloatOption("Option::minLeptonCleaningPt"); 
  minLeptonCutPt      = gConfigParser -> readFloatOption("Option::minLeptonCutPt");
  minJetCutPt         = gConfigParser -> readFloatOption("Option::minJetCutPt");
  usePuppiAsDefault   = gConfigParser -> readBoolOption("Option::usePuppiAsDefault");
  leptonIsoCut_mu     = gConfigParser -> readFloatOption("Option::leptonIsoCutMu");
  leptonIsoCut_el     = gConfigParser -> readFloatOption("Option::leptonIsoCutEl");
  leptonIsoCutLoose   = gConfigParser -> readFloatOption("Option::leptonIsoCutLoose");

  // output directory
  string outputPlotDirectory = gConfigParser -> readStringOption("Output::outputPlotDirectory");
  system(("mkdir -p output/"+outputPlotDirectory).c_str());
  system(("rm -r output/"+outputPlotDirectory+"/*").c_str());
  system(("mkdir -p output/"+outputPlotDirectory+"/xs").c_str());
  system(("mkdir -p output/"+outputPlotDirectory+"/norm").c_str());


  ///// Start the analysis  
  map<string,TH1F*> histoCutEff ;

  TChain* chain = new TChain (treeName.c_str()) ;  
  chain->Add ((InputBaseDirectory+"/*.root").c_str()) ;
  int totEvent = chain->GetEntries();

  readTree* reader  = new readTree((TTree*)(chain));

  cout<<"Lumi (fb-1) "<<lumi/1000<<" entries before "<<totEvent<<" cross section "<<CrossSection<<" Nevents before selections "<<lumi*CrossSection<<" weight "<<lumi*CrossSection/float(totEvent)<<endl;

  float weight = 1.0*lumi*CrossSection/float(totEvent) ;

  // make the plot container 
  vector<histoContainer> plotVector;

//.........这里部分代码省略.........

void produce_bgtemplate(){//main  

//  TFile *efile = TFile::Open("plot_elefakepho_test.root","read");
//  TTree* etree = (TTree*)efile->Get("FakeRateTree");

 TChain *etree = new TChain("BGTree");
 etree->Add("../data/plot_bgtemplate_singleMu.root");
// etree->Add("../data/plot_elefakepho_DY_74X.root"); 


  int tracks=0; 
  int nVertex=0;
  float mass_denmg=0,mass_nummg=0,et_denmg=0,et_nummg=0, rap_denmg=0, rap_nummg=0;
   
  etree->SetBranchAddress("tracks",&tracks);
  etree->SetBranchAddress("nVertex",&nVertex);
  etree->SetBranchAddress("mass_denmg", &mass_denmg);
  etree->SetBranchAddress("mass_nummg", &mass_nummg);
  etree->SetBranchAddress("et_denmg", &et_denmg);
  etree->SetBranchAddress("et_nummg", &et_nummg);
  etree->SetBranchAddress("rap_denmg", &rap_denmg);
  etree->SetBranchAddress("rap_nummg", &rap_nummg);

  float PtBins[]={25,30,35,40,50,60,70};
  //float PtBins[]={25,30,35,40,45,50,60,70,90,110};
  float NtrkBins[] = {10,20,30,40,50,60,70,90,120};
  float EtaBins[] = {0.0,0.1,0.2,0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3};
  float VertexBins[]={0,4,8,10,11,12,13,14,15,16,17,18,20};
  unsigned nPtBins = sizeof(PtBins)/sizeof(float);
  unsigned nNtrkBins = sizeof(NtrkBins)/sizeof(float);
  unsigned nEtaBins = sizeof(EtaBins)/sizeof(float);
  unsigned nVertexBins = sizeof(VertexBins)/sizeof(float);

  TFile *outputfile = TFile::Open("histo_bgtemplate_76X.root","RECREATE");
  outputfile->cd();

  TH1F* den_mg_pt[nPtBins]; 
  TH1F* den_mg_trk[nNtrkBins];
  TH1F* den_mg_eta[nEtaBins];
  TH1F* den_mg_vtx[nVertexBins];
  TH1F* num_mg_pt[nPtBins]; 
  TH1F* num_mg_trk[nNtrkBins];
  TH1F* num_mg_eta[nEtaBins];
  TH1F* num_mg_vtx[nVertexBins];
 
  std::ostringstream binname; 
  for(unsigned iPt(0); iPt < nPtBins; iPt++){
    binname.str("");
    if(iPt != nPtBins-1)binname << "den_mg_pt-" << PtBins[iPt] << "-" << PtBins[iPt+1];
    else binname << "den_mg_pt-" << PtBins[iPt] << "-inf";
    den_mg_pt[iPt] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
    binname.str("");
    if(iPt != nPtBins-1)binname << "num_mg_pt-" << PtBins[iPt] << "-" << PtBins[iPt+1];
    else binname << "num_mg_pt-" << PtBins[iPt] << "-inf";
    num_mg_pt[iPt] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
  }
    

  for(unsigned iNtrk(0); iNtrk < nNtrkBins; iNtrk++){
    binname.str("");
    if(iNtrk != nNtrkBins-1)binname << "den_mg_trk-" << NtrkBins[iNtrk] << "-" << NtrkBins[iNtrk+1];
    else binname << "den_mg_trk-" << NtrkBins[iNtrk] << "-inf";
    den_mg_trk[iNtrk] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
    binname.str("");
    if(iNtrk != nNtrkBins-1)binname << "num_mg_trk-" << NtrkBins[iNtrk] << "-" << NtrkBins[iNtrk+1];
    else binname << "num_mg_trk-" << NtrkBins[iNtrk] << "-inf";
    num_mg_trk[iNtrk] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
  }


  for(unsigned iEta(0); iEta < nEtaBins; iEta++){
    binname.str("");
    if(iEta != nEtaBins-1)binname << "den_mg_eta-" << EtaBins[iEta] << "-" << EtaBins[iEta+1];
    else binname << "den_mg_eta-" << EtaBins[iEta] << "-inf";
    den_mg_eta[iEta] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
    binname.str("");
    if(iEta != nEtaBins-1)binname << "num_mg_eta-" << EtaBins[iEta] << "-" << EtaBins[iEta+1];
    else binname << "num_mg_eta-" << EtaBins[iEta] << "-inf";
    num_mg_eta[iEta] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
  }

  for(unsigned iVertexBins(0); iVertexBins < nVertexBins; iVertexBins++){
    binname.str("");
    if(iVertexBins != nVertexBins-1)binname << "den_mg_vtx-" << VertexBins[iVertexBins] << "-" << VertexBins[iVertexBins+1];
    else binname << "den_mg_vtx-" << VertexBins[iVertexBins] << "-inf";
    den_mg_vtx[iVertexBins] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
    binname.str("");
    if(iVertexBins != nVertexBins-1)binname << "num_mg_vtx-" << VertexBins[iVertexBins] << "-" << VertexBins[iVertexBins+1];
    else binname << "num_mg_vtx-" << VertexBins[iVertexBins] << "-inf";
    num_mg_vtx[iVertexBins] = new TH1F(binname.str().c_str(), binname.str().c_str(), 60, 60,120);
  }

  const unsigned nEvts = etree->GetEntries(); 
  for(unsigned iEvt(0); iEvt < nEvts; iEvt++){
    etree->GetEntry(iEvt);

    if(mass_denmg > 40.0){
      den_mg_pt[findIndex(PtBins,et_denmg,nPtBins)]->Fill(mass_denmg);
      den_mg_trk[findIndex(NtrkBins,tracks,nNtrkBins)]->Fill(mass_denmg);
      den_mg_eta[findIndex(EtaBins,fabs(rap_denmg),nEtaBins)]->Fill(mass_denmg); 
//.........这里部分代码省略.........

void srcPositionStudy(Int_t binWidth, TString source, TString geometry) {

  int numFiles = 200; 

  TString simLocation;
  TChain *chain = new TChain("anaTree");

  if (geometry==TString("2011-2012")) simLocation = TString(getenv("SIM_2011_2012"));
  else if (geometry==TString("2012-2013")) simLocation = TString(getenv("SIM_2012_2013"));
  else if (geometry==TString("2012-2013_ISOBUTANE")) simLocation = TString(getenv("SIM_2012_2013_ISOBUTANE"));
  else { std::cout << "BAD GEOMETRY\n"; exit(0); }

  for (int i=0; i<numFiles; i++) {
    chain->AddFile(TString::Format("%s/%s/analyzed_%i.root",simLocation.Data(),source.Data(),i));
  }
  

  Double_t maxEn = 1500.;
  Double_t minEn = 0.;
  Double_t fidMax = 55.;
  Int_t nHists = (int)(fidMax/binWidth);
  
  
  
  std::vector <TH1D*> histsE(nHists, 0);
  std::vector <TH1D*> histsW(nHists, 0);
  std::vector <Double_t> rmin(nHists,0);
  std::vector <Double_t> rmax(nHists,0);
  std::vector <Double_t> rmid(nHists,0);

  //final means and errors
  std::vector < Double_t > EastMeans(nHists,0.);
  std::vector < Double_t > WestMeans(nHists,0.);
  
  std::vector < Double_t > EastMeanErrors(nHists,0.);
  std::vector < Double_t > WestMeanErrors(nHists,0.);

  for (Int_t i=0; i<nHists; i++) {

    rmin[i] = i*binWidth;
    rmax[i] = (i+1)*binWidth;
    rmid[i] = (double)rmin[i] + (double)binWidth/2.;

    histsE[i] = new TH1D(TString::Format("his%iE",i), 
			 TString::Format("%s %i mm radius bins East", source.Data(),binWidth),
			550, 0., 1100.);
    histsW[i] = new TH1D(TString::Format("his%iW",i), 
			 TString::Format("%s %i mm radius bins West", source.Data(),binWidth),
			550, 0., 1100.);

    histsE[i]->SetLineColor(i+1);
    histsW[i]->SetLineColor(i+1);
  }

  Double_t primPos[4];

  // Set the addresses of the information read in from the simulation file
  chain->SetBranchAddress("MWPCEnergy",&mwpcE);
  chain->SetBranchAddress("time",&Time);
  chain->SetBranchAddress("Edep",&edep);
  chain->SetBranchAddress("EdepQ",&edepQ);
  chain->SetBranchAddress("MWPCPos",&mwpc_pos);
  chain->SetBranchAddress("ScintPos",&scint_pos);
  chain->SetBranchAddress("primKE",&primKE);
  chain->SetBranchAddress("primTheta",&primTheta);
  chain->SetBranchAddress("primPos",&primPos);


  //Get total number of events in TChain
  UInt_t nevents = chain->GetEntries();
  cout << "events = " << nevents << endl;

  for (Int_t i=0; i<nevents; i++) {
    
    chain->GetEvent(i);

    Int_t nBin = primPos[3]*1000./binWidth;

    if (edepQ.EdepQE>0. && primKE<maxEn && primKE>minEn && mwpcE.MWPCEnergyE>0.1 && primTheta>TMath::Pi()/2.) histsE[nBin]->Fill(edepQ.EdepQE);
    if (edepQ.EdepQW>0. && primKE<maxEn && primKE>minEn && mwpcE.MWPCEnergyW>0.1 && primTheta<TMath::Pi()/2.) histsW[nBin]->Fill(edepQ.EdepQW);

    if (i%100000==0) std::cout << "*";
  }
  
  std::cout << std::endl;

  TCanvas *c1 = new TCanvas("c1","c1",1600,1200);
  c1->Divide(2,2);
  //TCanvas *c2 = new TCanvas("c2");
  
  //histsE[1]->Draw("SAME");
  
  Double_t refMeanE = 0., refMeanW = 0.; //This will hold the mean of the center pixel

  std::cout << nHists << endl;  
  for (Int_t i=0; i<nHists; i++) {
    c1->cd(1);
    histsE[i]->Draw("SAME");
    
    EastMeans[i] = histsE[i]->GetMean();
//.........这里部分代码省略.........

void gep_trigger_analysis_elastic( const char *rootfilename, const char *logicfilename_ecal, const char *logicfilename_hcal, const char *thresholdfilename_ecal, const char *thresholdfilename_hcal, const char *outputfilename, double thetacaldeg=29.0, int pheflag=0, const char *assocfilename="ECAL_HCAL_correlations_nophe.txt", int Q2cut=0 ){

  double nominal_threshold_HCAL = 0.5;
  double nominal_threshold_ECAL = 0.9;
  
  double thetacal = thetacaldeg*PI/180.0;
  
  TFile *fout = new TFile(outputfilename,"RECREATE");
  TChain *C = new TChain("T");
  C->Add(rootfilename);

  gep_tree_elastic *T = new gep_tree_elastic( C );

  G4SBSRunData *rd;

  long ngen = 0;
  int nfiles = 0;
  
  TObjArray *FileList = C->GetListOfFiles();
  TIter next(FileList);

  TChainElement *chEl = 0;

  set<TString> bad_file_list;
  
  while( (chEl=(TChainElement*)next() )){
    TFile newfile(chEl->GetTitle());
    newfile.GetObject("run_data",rd);
    if( rd ){
      ngen += rd->fNtries;
      nfiles++;
    } else {
      bad_file_list.insert( chEl->GetTitle());
    }
  }

  cout << "number of generated events = " << ngen << endl;
  
  set<int> list_of_nodes_ecal;
  map<int, set<int> > cells_logic_sums_ecal; //mapping between node numbers and cell numbers
  map<int, double> logic_mean_ecal; //mean peak positions by node number
  map<int, double> logic_sigma_ecal; //peak width by node number
  map<int, double> threshold_ecal; //threshold by node number
  map<std::pair<int,int>, int > cell_rowcol_ecal; //cell numbers mapped by unique row and column pairs
  map<int,set<int> > nodes_cells_ecal; //mapping of nodes by cell number:
  map<int,int> rows_cells_ecal;
  map<int,int> cols_cells_ecal;
  map<int,double> xcells_ecal;
  map<int,double> ycells_ecal;
  
  //keep track of min and max x by row number:
  double ycellmin,ycellmax;
  map<int,double> ycell_rows;
  map<int,double> cellsize_rows;
  map<int,double> xcellmin_rows;
  map<int,double> xcellmax_rows;
  
  int minrow=1000,maxrow=-1;
  
  set<int> rows_ecal;
  map<int,set<int> > columns_rows_ecal;
  
  map<int,double> elastic_peak_new_ecal;
  map<int,double> sigma_new_ecal;
  map<int,double> threshold_new_ecal;
  
  ifstream logicfile_ecal(logicfilename_ecal);
  //ifstream thresholdfile(thresholdfilename);

  TString currentline;
  
  int current_node = 1;

  bool first_cell = true;
  
  while( currentline.ReadLine( logicfile_ecal ) ){
    if( !currentline.BeginsWith( "#" ) ){
      
      TObjArray *tokens = currentline.Tokenize(" ");
      int ntokens = tokens->GetEntries();
      if( ntokens >= 11 ){
	cout << currentline.Data() << ", ntokens = " << ntokens << endl;
	
	TString snode = ( (TObjString*) (*tokens)[0] )->GetString();
	int nodenumber = snode.Atoi();
	
	TString scell = ( (TObjString*) (*tokens)[1] )->GetString();
	int cellnumber = scell.Atoi();
	
	TString speakpos = ( (TObjString*) (*tokens)[8] )->GetString();
	double mean = speakpos.Atof();
	
	TString ssigma = ( (TObjString*) (*tokens)[9] )->GetString();
	double sigma = ssigma.Atof();

	TString sthreshold = ( (TObjString*) (*tokens)[10] )->GetString();
	double threshold = sthreshold.Atof();

	TString srow = ( (TObjString*) (*tokens)[2] )->GetString();
	TString scol = ( (TObjString*) (*tokens)[3] )->GetString();
//.........这里部分代码省略.........

void anaFragShape(TString infile="dj_HCPR-GoodTrkAndPixel_CleanEvt1130.root")
{
  TChain * djcalo = new TChain("djcalo/djTree");
  djcalo->Add(infile);
  aliases_dijet(djcalo);
  TChain * djcalopix = new TChain("djcalo_pxl/djTree");
  djcalopix->Add(infile);
  aliases_dijet(djcalopix);
  TString evtSel("(cent<10 && nljet>120 && abs(nljeta)<2 && aljet>50 && abs(aljeta)<2 && jdphi>2.5 && Aj>0.24)");
  TH1::SetDefaultSumw2();

  TCanvas * c0 = new TCanvas("c0","c0",500,500);
  djcalo->Draw("Aj>>hAj(20,0,1)",evtSel);
  djcalopix->Draw("Aj",evtSel,"Esame");

  TCanvas * c2 = new TCanvas("c2","c2",500,500);
  TH1D * hJDPhi = new TH1D("hJDPhi","hJDPhi",50,0,3.1416);
  TH1D * hJDPhi2 = new TH1D("hJDPhi2","hJDPhi",50,0,3.1416);
  Float_t numDJ = djcalo->Draw("jdphi>>hJDPhi",evtSel);
  Float_t numDJ2 = djcalopix->Draw("jdphi>>hJDPhi2",evtSel);
  cout << "num djs: " << numDJ << "  djs(pix)" << numDJ2 << endl;

  TH1D * hPNDRTrk = new TH1D("hPNDRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPNDRPix = new TH1D("hPNDRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRTrk = new TH1D("hPADRTrk",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  TH1D * hPADRPix = new TH1D("hPADRPix",";R(trk,jet);1/(N_{DJ} 2#piR) dp_{T}^{Trk}/dR",10,0,TMath::PiOver2());
  djcalo->Draw("pndr>>hPNDRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalopix->Draw("pndr>>hPNDRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*pndr))",evtSel.Data()));
  djcalo->Draw("padr>>hPADRTrk",Form("(%s&&ppt>=1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  djcalopix->Draw("padr>>hPADRPix",Form("(%s&&ppt>=0.3&&ppt<1.2)*(ppt/(TMath::TwoPi()*padr))",evtSel.Data()));
  hPNDRTrk->Scale(1./(numDJ*hPNDRTrk->GetBinWidth(1)));
  hPNDRPix->Scale(1./(numDJ*hPNDRPix->GetBinWidth(1)));
  hPADRTrk->Scale(1./(numDJ*hPADRTrk->GetBinWidth(1)));
  hPADRPix->Scale(1./(numDJ*hPADRPix->GetBinWidth(1)));
  hPNDRPix->SetMinimum(50);
  hPNDRPix->SetMaximum(100);
  hPNDRTrk->SetMinimum(10);
  hPNDRTrk->SetMaximum(2000);
  hPNDRPix->SetMarkerStyle(kFullCircle);
  hPNDRTrk->SetMarkerStyle(kFullSquare);
  hPNDRPix->SetMarkerColor(kRed);
  hPNDRTrk->SetMarkerColor(kRed);
  hPADRPix->SetMarkerStyle(kOpenCircle);
  hPADRTrk->SetMarkerStyle(kOpenSquare);
  hPADRPix->SetMarkerColor(kBlue);
  hPADRTrk->SetMarkerColor(kBlue);
  TCanvas * c3 = new TCanvas("c3","c3",500,500);
  hPNDRPix->Draw("E");
  hPADRPix->Draw("Esame");
  TCanvas * c3_2 = new TCanvas("c3_2","c3_2",500,500);
  c3_2->SetLogy();
  hPNDRTrk->Draw("E");
  hPADRTrk->Draw("Esame");
}

void tnpScale_IDISO( int leptype = 1, bool printplot = false ) {

  cout << endl;
  cout << "-------------------" << endl;
  if     ( leptype == 0 ) cout << "Doing electrons" << endl;
  else if( leptype == 1 ) cout << "Doing muons"      << endl;
  else{
    cout << "ERROR! unrecognized leptype " << leptype << endl;
    exit(0);
  }
  cout << "-------------------" << endl;

  //----------------------------------------
  // Files
  //----------------------------------------

  TChain *chmc   = new TChain("leptons");
  TChain *chdata = new TChain("leptons");

  char* version = (char*) "V00-00-06";
  char* suffix = "";
  //char* suffix = "_2jets";
  //char* suffix = "_probept100";


  chmc->  Add(Form("smurf/ZJetsFull_%s/merged%s.root",version,suffix));

  if( leptype == 1 ){
    chdata->Add(Form("smurf/SingleMu2012AFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleMu2012BFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleMu2012CFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleMu2012DFull_%s/merged_json%s.root",version,suffix));
  }
  else{
    chdata->Add(Form("smurf/SingleEl2012AFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleEl2012BFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleEl2012CFull_%s/merged_json%s.root",version,suffix));
    chdata->Add(Form("smurf/SingleEl2012DFull_%s/merged_json%s.root",version,suffix));
  }

  //----------------------------------------
  // bins 
  //----------------------------------------

  // float ptbin[]  = {10., 15., 20., 30., 40., 50., 7000.};
  // float ptbin[]  = { 30. , 40. , 50. , 60. , 80.0 , 100.0 , 120.0 , 150.0 , 7000.};
  // float etabin[] = {0, 0.8, 1.5, 2.1};
  // int nptbin=8;
  // int netabin=3;

  // float ptbin[]  = { 20., 30. , 40. , 50. , 60. , 80.0 , 100.0 , 150.0 , 200.0 , 300.0 , 500.0 , 1000.0 , 10000000.0};
  // int   nptbin   = 12;

  // float etabin[] = {0,2.1};
  // int   netabin  = 1;

  float ptbin[] = { 20., 30. , 40. , 50. , 60. , 80.0 , 100.0 , 150.0 , 200.0 , 300.0, 10000.0};
  int   nptbin  = 10;

  float etabin[4];
  int   netabin = 0;

  if( leptype == 1 ){
    cout << "DOING MUON ETA BINS" << endl;
    netabin=3;
    etabin[0] = 0.0;
    etabin[1] = 0.8;
    etabin[2] = 1.5;
    etabin[3] = 2.1;
  }

  if( leptype == 0 ){
    cout << "DOING ELECTRON ETA BINS" << endl;
    netabin=2;
    etabin[0] = 0.0;
    etabin[1] = 0.8;
    etabin[2] = 1.4442;

    // netabin=1;
    // etabin[0] = 0.0;
    // etabin[1] = 1.4442;
  }


  //deno
  TH2F *hmcid_deno 	= new TH2F("hmcid_deno"   , "hmcid_deno"   , nptbin, ptbin, netabin, etabin);
  TH2F *hmciso_deno 	= new TH2F("hmciso_deno"  , "hmciso_deno"  , nptbin, ptbin, netabin, etabin);
  TH2F *hdataid_deno 	= new TH2F("hdataid_deno" , "hdataid_deno" , nptbin, ptbin, netabin, etabin);
  TH2F *hdataiso_deno	= new TH2F("hdataiso_deno", "hdataiso_deno", nptbin, ptbin, netabin, etabin);
  hmcid_deno->Sumw2();
  hmciso_deno->Sumw2();
  hdataid_deno->Sumw2();
  hdataiso_deno->Sumw2();

  //num
  TH2F *hmcid_num 	= new TH2F("hmcid_num"    , "hmcid_num"    , nptbin, ptbin, netabin, etabin);
  TH2F *hmciso_num 	= new TH2F("hmciso_num"   , "hmciso_num"   , nptbin, ptbin, netabin, etabin);
  TH2F *hdataid_num 	= new TH2F("hdataid_num"  , "hdataid_num"  , nptbin, ptbin, netabin, etabin);
  TH2F *hdataiso_num 	= new TH2F("hdataiso_num" , "hdataiso_num" , nptbin, ptbin, netabin, etabin);
  hmcid_num->Sumw2();
//.........这里部分代码省略.........

void sampleComparison(){

  gROOT->SetStyle("Plain");
  
  /*
  TChain *uno = new TChain("sT");
  uno->Add("/raid/sguazz/split/MuGun_NomGeo.root"); 
  TString n1("MuGun NomGeo");
  
  TChain *due = new TChain("sT");
  due->Add("/raid/sguazz/split/minBias_nH4_v2.root"); 
  TString n2("MinBias");
  */
  TChain *uno = new TChain("sT");
  uno->Add("/raid/sguazz/split/minBias_nH4_v2.root"); 
  TString n1("minBias MC");
  
  TChain *due = new TChain("sT");
  due->Add("/raid/sguazz/split/minbias_run2012Cpart1_nH4.root"); 
  TString n2("minBias run2012C");
    
  std::cout << " Sample uno:"  << uno->GetEntries() << " due:" << due->GetEntries() << std::endl;
  
  TH1F * chi2SplitUno = new TH1F("chi2Split",";#chi^{2}_{SPLIT}/d.o.f.;",40, 0.,20.);
  TH1F * chi2SplitDue = new TH1F("chi2SplitDue",";#chi^{2}_{SPLIT}/d.o.f.;",40, 0.,20.);
  
  TH1F * etaUno = new TH1F("eta",";#eta;",50,-2.5,2.5);
  TH1F * etaDue = new TH1F("etaDue",";#eta;",50,-2.5,2.5);
  
  TH1F * phiUno = new TH1F("phi",";#phi;",50,-M_PI,M_PI);
  TH1F * phiDue = new TH1F("phiDue",";#phi;",50,-M_PI,M_PI);
  
  TH1F * ptUno = new TH1F("pt",";P_{T} [GeV];",50,0.7,1.3);
  TH1F * ptDue = new TH1F("ptDue",";P_{T} [GeV];",50,0.7,1.3);
  
  TH1F * ptNormUno = new TH1F("ptNorm",";P_{T} [GeV];",200,0.9,1.1);
  TH1F * ptNormDue = new TH1F("ptNormDue",";P_{T} [GeV];",200,0.9,1.1);
  
  TH1F * nUno = new TH1F("n",";N_{SPLIT};",7,-0.5,6.5);
  TH1F * nDue = new TH1F("nDue",";N_{SPLIT};",7,-0.5,6.5);

  TH1F * dzUno = new TH1F("dz",";dz;",40,-20.,20.);
  TH1F * dzDue = new TH1F("dzDue",";dz;",40,-20.,20.);
  
  TH1F * routUno = new TH1F("rOut",";rOut;",30,100.,115.);
  TH1F * routDue = new TH1F("rOutDue",";rOut;",30,100.,115.);
  
  // MC
  //    TString cut("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&pLossSim>0.&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&chi2Split>0.");
  //    TString cutNoPt("chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&pLossSim>0.&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&chi2Split>0.");
  // Data
  TString        cut("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&(chi2Split/(NSplit-2.))>0.&&(chi2Split/(NSplit-2.))<4.");
  TString    cutROut("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&(chi2Split/(NSplit-2.))>0.&&(chi2Split/(NSplit-2.))<4.&&eta<0.9&&eta>-0.9");
  TString    cutNoDz("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&maxchi2<4.&&(chi2Split/(NSplit-2.))>0.&&(chi2Split/(NSplit-2.))<4.");
  TString   cutNoChi("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&(chi2Split/(NSplit-2.))>0.");
  TString    cutNoPt(                "chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.&&(chi2Split/(NSplit-2.))>0.&&(chi2Split/(NSplit-2.))<4.");
  TString cutNoSplit("pt>0.9&&pt<1.1&&chi2<2.&&rIn<8.&&((rOut>100.)||(zOut>260.)||(zOut<-260.))&&(nHitIna+nHitMis+nHitBad)==0&&dz<5.&&dz>-5.&&maxchi2<4.");
  
  int maxentries = 10000000;
  //int maxentries = 999999999;
  
  //TString assoc("&&pLossSim>0.");
  TString assoc("");
  
    uno->Draw("(chi2Split/(NSplit-2.))>>chi2Split",cutNoChi+assoc,"goff",maxentries);
    due->Draw("(chi2Split/(NSplit-2.))>>chi2SplitDue",cutNoChi+assoc,"goff",maxentries);
    draw(chi2SplitUno, chi2SplitDue, n1, n2);
    
    uno->Draw("eta>>eta",cut+assoc,"goff",maxentries);
    due->Draw("eta>>etaDue",cut+assoc,"goff",maxentries);
    draw(etaUno, etaDue, n1, n2);

    uno->Draw("pt>>pt",cutNoPt+assoc,"goff",maxentries);
    due->Draw("pt>>ptDue",cutNoPt+assoc,"goff",maxentries);
    draw(ptUno, ptDue, n1, n2);

    uno->Draw("pt>>pt",cut+assoc,"goff",maxentries);
    due->Draw("pt>>ptDue",cut+assoc,"goff",maxentries);
    draw(ptUno, ptDue, n1, n2);
  
  uno->Draw("dz>>dz",cutNoDz+assoc,"goff",maxentries);
  due->Draw("dz>>dzDue",cutNoDz+assoc,"goff",maxentries);
  draw(dzUno, dzDue, n1, n2);
  
  uno->Draw("rOut>>rOut",cutROut+assoc,"goff",maxentries);
  due->Draw("rOut>>rOutDue",cutROut+assoc,"goff",maxentries);
  draw(routUno, routDue, n1, n2);
  
    uno->Draw("phi>>phi",cut+assoc,"goff",maxentries);
    due->Draw("phi>>phiDue",cut+assoc,"goff",maxentries);
    draw(phiUno, phiDue, n1, n2);
    
    uno->Draw("NSplit>>n",cutNoSplit+assoc,"goff",maxentries);
    due->Draw("NSplit>>nDue",cutNoSplit+assoc,"goff",maxentries);
    draw(nUno, nDue, n1, n2);

}

runLooper_electron(){
  
  
  //  gStyle->SetTitleX(0.1f);
  //   gStyle->SetTitleW(0.5f);
  
  TChain *Chain = new TChain("Events");
 
  // Chain->Add("/store/disk01/yanjuntu/SingleElectron83055667ceda31da68d5aa532f72919d/preprocessing/ntuple*.root");  //tas03
  Chain->Add("/hadoop/cms/store/user/yanjuntu/conversion_sample_cms2-V01-03-01/SingleElectron83055667ceda31da68d5aa532f72919d/preprocessing/ntuple*.root");   //uaf6
  //Chain->Add("/store/disk01/yanjuntu/2_2_10/singleElectron_ntuple_v3.root");
 
  
 
  TCanvas* c1 = new TCanvas("c1","c1");
  TH1F hist_singleElectron("singleElectron", "; Eta", 100,-3,3);
  Chain->Draw("els_p4.eta() >> singleElectron", susy_baseline_2);
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_eta_susy.eps");
  float ntotal = hist_singleElectron.Integral();
 
  
  TH1F hist_singleElectron("singleElectron", "; Number of valid pixel hits", 16,0,8);
  Chain->Draw("els_valid_pixelhits >> singleElectron",susy_baseline_2 && forward);
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_valid_pixelhits_susy.eps");
  float ntotal_forward = hist_singleElectron.Integral();
  float n_removed_pixelhit = hist_singleElectron.GetBinContent(0);
  
  TH1F hist_singleElectron("singleElectron", "; DetID", 10,0,5);
  Chain->Draw("els_layer1_det  >> singleElectron", susy_baseline_2 && forward &&"(els_valid_pixelhits==2||els_valid_pixelhits==1)");
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_layer1_detid_susy.eps");

  TH1F hist_singleElectron("singleElectron", "; Layer of the first valid pixel hit (Barrel)", 10,0,5);
  Chain->Draw("els_layer1_layer  >> singleElectron", susy_baseline_2 && forward &&"(els_valid_pixelhits==2||els_valid_pixelhits==1) && els_layer1_det==1");
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_layer1_layer_barrel_susy.eps");
  float n_removed_layer = hist_singleElectron.Integral(hist_singleElectron.FindBin(2),hist_singleElectron.GetNbinsX()+1);
 

  TH1F hist_singleElectron("singleElectron", "; Charge of the first valid pixel hit (Disk)", 100,10000,200000);
  Chain->Draw("els_layer1_charge  >> singleElectron", susy_baseline_2 && forward && "(els_valid_pixelhits==2||els_valid_pixelhits==1) &&els_layer1_det==2");
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_layer1_charge_disk_susy.eps");
  float n_removed_charge = hist_singleElectron.Integral(hist_singleElectron.FindBin(40000),hist_singleElectron.GetNbinsX()+1);
 
  
  TH2F hist_Electron("Electron", "; Charge of the first valid pixel hit (Disk); Eta", 100, 10000,100000, 300,1.5,2.8);
  Chain->Draw("abs(els_p4.eta()):els_layer1_charge  >> Electron",susy_baseline_2 && forward && "(els_valid_pixelhits==2||els_valid_pixelhits==1) &&els_layer1_det==2", "box");
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_layer1_chargeVsEta_disk_susy.eps");
  
  TH1F hist_singleElectron("singleElectron", "; Number of expected inner layers", 16,0,8);
  Chain->Draw("els_n_inner_layers >> singleElectron", susy_baseline_2 );
  c1->SaveAs("~/sample_plots/conversions/3_1_x/singleElectron_ele_ninner_susy.eps");
  float n_removed_m2_0 = hist_singleElectron.Integral(hist_singleElectron.FindBin(1),hist_singleElectron.GetNbinsX()+1);
  float n_removed_m2_1 = hist_singleElectron.Integral(hist_singleElectron.FindBin(2),hist_singleElectron.GetNbinsX()+1);
 
  float ineff_forward = (n_removed_pixelhit+n_removed_layer+n_removed_charge)/ntotal_forward;
  float ineff = (n_removed_pixelhit+n_removed_layer+n_removed_charge)/ntotal;
  std::cout<< "total electrons = "<<ntotal <<std::endl;
  std::cout<< "total electrons in the forward region = "<<ntotal_forward <<std::endl;
  std::cout<< "electrons removed by cutting on the number of pixel hits  = "<< n_removed_pixelhit<<std::endl;
  std::cout<< "electrons removed by cutting on the layer  = "<< n_removed_layer<<std::endl;
  std::cout<< "electrons removed by cutting on the charge  = "<< n_removed_charge<<std::endl;
  std::cout<< "the inefficeincy of prompt electrons in forward region = "<<ineff_forward <<std::endl;
  std::cout<< "the inefficeincy of prompt electrons = "<<ineff <<std::endl;
  
  float ineff_m2_0 = n_removed_m2_0/ntotal;
  float ineff_m2_1 = n_removed_m2_1/ntotal;
  std::cout<< "electrons removed by cutting ninner layers (loose)  = "<<n_removed_m2_0 <<std::endl;
  std::cout<< "electrons removed by cutting ninner layers (tight)  = "<<n_removed_m2_1 <<std::endl;
  std::cout<< "method 2 loose cut: the inefficeincy of prompt electrons = "<<ineff_m2_0 <<std::endl;
  std::cout<< "method 2 tight cut: the inefficeincy of prompt electrons = "<<ineff_m2_1 <<std::endl;
}

int main(int argc, char *argv[])
{
  char appName[] = "root2lhco";
  stringstream message;
  FILE *outputFile = 0;
  TChain *inputChain = 0;
  LHCOWriter *writer = 0;
  ExRootTreeReader *treeReader = 0;
  Long64_t entry, allEntries;

  if(argc < 2 || argc > 3)
  {
    cerr << " Usage: " << appName << " input_file" << " [output_file]" << endl;
    cerr << " input_file - input file in ROOT format," << endl;
    cerr << " output_file - output file in LHCO format," << endl;
    cerr << " with no output_file, or when output_file is -, write to standard output." << endl;
    return 1;
  }

  signal(SIGINT, SignalHandler);

  gROOT->SetBatch();

  int appargc = 1;
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  try
  {
    cerr << "** Reading " << argv[1] << endl;
    inputChain = new TChain("Delphes");
    inputChain->Add(argv[1]);

    ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);

    if(argc == 2 || strcmp(argv[2], "-") == 0)
    {
      outputFile = stdout;
    }
    else
    {
      outputFile = fopen(argv[2], "w");

      if(outputFile == NULL)
      {
        message << "can't open " << argv[2];
        throw runtime_error(message.str());
      }
    }

    fprintf(outputFile, "   #  typ      eta      phi      pt    jmas   ntrk   btag  had/em   dum1   dum2\n");

    allEntries = treeReader->GetEntries();
    cerr << "** Input file contains " << allEntries << " events" << endl;

    if(allEntries > 0)
    {
      // Create LHC Olympics converter:
      writer = new LHCOWriter(treeReader, outputFile);

      ExRootProgressBar progressBar(allEntries - 1);
      // Loop over all events
      for(entry = 0; entry < allEntries && !interrupted; ++entry)
      {
        if(!treeReader->ReadEntry(entry))
        {
          cerr << "** ERROR: cannot read event " << entry << endl;
          break;
        }

        writer->ProcessEvent();

        progressBar.Update(entry);
      }
      progressBar.Finish();

      delete writer;
    }

    cerr << "** Exiting..." << endl;

    if(outputFile != stdout) fclose(outputFile);
    delete treeReader;
    delete inputChain;
    return 0;
  }
  catch(runtime_error &e)
  {
    if(writer) delete writer;
    if(treeReader) delete treeReader;
    if(inputChain) delete inputChain;
    cerr << "** ERROR: " << e.what() << endl;
    return 1;
  }
}