C++ TBenchmark::Stop方法代码示例

void AuxFchTemplateFitsData_SIHIH(TString strChannel, TString strVGamma)
{
  TBenchmark time;
  time.Start("time");
  std::cout<<"CPU time = "<<time.GetCpuTime("time")<<", Real time = "<<time.GetRealTime("time")<<std::endl;  

  FullChain fch;
  TConfiguration::AnalysisParameters anPars;
  fch.SetDefaultFullChainParameters(anPars,"phoEt");

  TConfiguration conf;
  int channel;
  int vgamma;
  if (strChannel=="MUON") anPars.channel=conf.MUON;
  if (strChannel=="ELECTRON") anPars.channel=conf.ELECTRON;
  if (strVGamma=="WGamma") anPars.vgamma=conf.W_GAMMA;
  if (strVGamma=="ZGamma") anPars.vgamma=conf.Z_GAMMA;

  anPars.templFits=conf.TEMPL_SIHIH;

  AuxTemplates(anPars,0);

  time.Stop("time");
  std::cout<<"CPU time = "<<time.GetCpuTime("time")<<", Real time = "<<time.GetRealTime("time")<<std::endl;  
}

void testJetCorrectorParameters::benchmarkBinIndex(bool is3D) {
    float oldCPU = 0, newCPU = 0, oldReal = 0, newReal = 0;
    unsigned int ntests = (is3D) ? 100000 : 1000000;
    if(is3D) fX = {5.0,50.0,100.0};
    else     fX = {5.0};
    setupCorrector(is3D);

    cout << endl << "testJetCorrectorParameters::benchmarkBinIndex NTests = " << ntests << endl;
    cout << "testJetCorrectorParameters::benchmarkBinIndex benchmarking binIndex with file " << filename << " ... " << flush;
    m_benchmark->Reset();
    m_benchmark->Start("event");
    for(unsigned int i=0; i<ntests; i++) {
       L1JetPar->binIndex(fX);
    }
    m_benchmark->Stop("event");
    cout << "DONE" << endl << "testJetCorrectorParameters::benchmarkBinIndex" << endl 
         << "\tCPU time = " << m_benchmark->GetCpuTime("event")/double(ntests) << " s" << endl
         << "\tReal time = " << m_benchmark->GetRealTime("event")/double(ntests) << " s" << endl;
    oldCPU = m_benchmark->GetCpuTime("event")/double(ntests);
    oldReal = m_benchmark->GetRealTime("event")/double(ntests);

    cout << "testJetCorrectorParameters::benchmarkBinIndex benchmarking binIndexN with file " << filename << " ... " << flush;
    m_benchmark->Reset();
    m_benchmark->Start("event");
    for(unsigned int i=0; i<ntests; i++) {
       L1JetPar->binIndexN(fX);
    }
    m_benchmark->Stop("event");
    cout << "DONE" << endl << "testJetCorrectorParameters::benchmarkBinIndex" << endl
         << "\tCPU time = " << m_benchmark->GetCpuTime("event")/double(ntests) << " s" << endl
         << "\tReal time = " << m_benchmark->GetRealTime("event")/double(ntests) << " s" << endl;
    newCPU = m_benchmark->GetCpuTime("event")/double(ntests);
    newReal = m_benchmark->GetRealTime("event")/double(ntests);

    cout << "testJetCorrectorParameters::benchmarkBinIndex" << endl
         << "\tCPU speedup = " << oldCPU/newCPU << endl
         << "\tReal speedup = " << oldReal/newReal << endl;

    CPPUNIT_ASSERT(oldCPU/newCPU >= 1.0);
    //CPPUNIT_ASSERT(oldReal/newReal >= 1.0); //this might fail not due to longer L1JetPar->binIndexN(fX) execution
    if (oldReal/newReal >= 1.0)
    cout << "newReal value increased oldReal, which might be due to system load" << endl;

    destroyCorrector();
}

void AuxFchBkgEtoGamma_NoWMtCut()
{
  TBenchmark time;
  time.Start("time");
  std::cout<<"CPU time = "<<time.GetCpuTime("time")<<", Real time = "<<time.GetRealTime("time")<<std::endl;  
 
  TConfiguration conf;
  FullChain fch;
  TConfiguration::AnalysisParameters anPars;
  fch.SetDefaultFullChainParameters(anPars,"phoEt");

  anPars.channel=conf.ELECTRON;
  anPars.vgamma=conf.W_GAMMA;
  TEtoGamma etg(anPars,conf.PRELIMINARY_FOR_E_TO_GAMMA_WITH_PSV_NO_WMT_CUT,0); 
  etg.ComputePlotSave();

  time.Stop("time");
  std::cout<<"CPU time = "<<time.GetCpuTime("time")<<", Real time = "<<time.GetRealTime("time")<<std::endl;  
}

void mergeHadoopFiles(const TString& indir, const TString& outpath) {
  // Benchmark
  TBenchmark *bmark = new TBenchmark();
  bmark->Start("benchmark");
  
  //gSystem->Load("../libMiniFWLite.so");
  TChain *chain = new TChain("Events");
  chain->SetMaxTreeSize(5000000000LL); //default is 100000000000LL

  std::cout << "Merging files from dir: " << indir << std::endl
	    << "Ouputting to: " << outpath << std::endl;

  chain->Add(indir + "/*.root");
  chain->Merge(outpath, "fast");

  std::cout << "Total events merged: " << chain->GetEntries() << std::endl;
  bmark->Stop("benchmark");
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f s", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f s", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
}

//.........这里部分代码省略.........
      // Progress
      lepfilter::progress( nEventsTotal, nEventsChain );
	  
      // Analysis Code
	  float weight = ss.scale1fb()*10.0;
	  if(ss.scale1fb() > 100000.) continue;  //excludes 5to10 and 10to20 EM Enriched
	  
	  bool jetptcut = false;
	  int jetidx = 0;

	  while( (jetidx < ss.jets().size()) && !jetptcut) //check to see if at least one jet w/ pt > 40
	  	{
	  	  if( ss.jets()[jetidx].pt() > 40. )
	  		{jetptcut = true;}
	  	  jetidx++;
	  	}
	  
	  if( !(jetptcut && ss.met() < 20. && ss.mt() < 20) )
	  	{continue;}

	  // if(ss.nFOs() != 1) //if more than 1 FO in event
	  // 	{continue;}
   
	  int nbtags = 0;
	  for(int i = 0; i < ss.jets().size(); i++){
		  if(ss.jets_disc()[i] > 0.814) nbtags++;
		}
	  if (nbtags > 2.) nbtags = 2;

	  if(ss.p4().pt() > 100. || ss.p4().pt() < 10. || fabs(ss.p4().eta()) > 2.4) //What do we want here? 
	  	{continue;}

	  //------------------------------------------------------------------------------------------
	  //---------------------------------Find e = f(Pt,eta)---------------------------------------
	  //------------------------------------------------------------------------------------------

	  //Find ratio of nonprompt leps passing tight to nonprompt leps passing at least loose.  This is the fake rate 
	  // Use lep_passes_id to see if num.  Use FO to see if Fakable Object (denom)
	  //Calculate e=Nt/(Nl) where l->loose  (as opposed to loose-not-tight).
	  
	  //Using gen level info to see if prompt -> no prompt contamination in measurement region
	  //everything else is RECO (p4, id, passes_id, FO, etc.)
	  

	  if( ss.motherID() != 1 && (doBonly==0 || ss.motherID() == -1) && (doConly==0 || ss.motherID() == -2) && (doLightonly==0 || ss.motherID() == 0) )  //if el is nonprompt (GEN info)
		{
		  if( abs( ss.id() ) == 11 ) // it's an el
			{
			  if(ss.FO_NoIso())  pTrelvsIso_histo_el->Fill( ss.iso(), ss.ptrelv1() );
			}
		  if( abs( ss.id() ) == 13 ) // it's a mu
			{
 			  if(ss.FO_NoIso())  pTrelvsIso_histo_mu->Fill( ss.iso(), ss.ptrelv1() );
			}
		} 

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

	}//end event loop
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }

  pTrelvsIso_histo_el->GetXaxis()->SetTitle("Iso");
  pTrelvsIso_histo_el->GetYaxis()->SetTitle("pTrel");
  pTrelvsIso_histo_mu->GetXaxis()->SetTitle("Iso");
  pTrelvsIso_histo_mu->GetYaxis()->SetTitle("pTrel");

  gStyle->SetOptStat(0);
  gStyle->SetPaintTextFormat("1.3f");

  TCanvas *c9=new TCanvas("c9","B pTrel vs Iso (el)",800,800);
  pTrelvsIso_histo_el->Draw("colz,texte");
  TCanvas *c10=new TCanvas("c10","B pTrel vs Iso (mu)",800,800);
  pTrelvsIso_histo_mu->Draw("colz,texte");

  //---save histos-------//
  TFile *OutputFile = new TFile(outfile,"recreate");
  pTrelvsIso_histo_el->Write();
  pTrelvsIso_histo_mu->Write();

  OutputFile->Close();

  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//.........这里部分代码省略.........
      LepTree::progress(nEventsTotal, nEventsChain);

      // Analysis Code

      if (evt_isRealData() && evt_run() < 273423) continue;  // Get runs after fixing the L1 interface problem

      int nevt = evt_event();
      if (nevt != evt_num) {
        h_muonCount->Fill(nMuonCount);
        nMuonCount = 0;
        isTriggerMuon = 0;
      }

      // --- New Tag & Probe ---
      if (abs(id()) != 13) continue;
      if (p4().pt() < 10) continue;
      if (fabs(p4().eta()) > 2.4 ) continue;
      if (charge()*tag_charge() > 0) continue;
      if (tag_p4().pt() < 25) continue;
      if (tag_RelIso03EA() > 0.1) continue;
      if (dilep_mass() < 81.2 || dilep_mass() > 101.2) continue;
      if (!pid_PFMuon()) continue;

      ++nMuonCount;
      // if (nMuonCount > 2) continue;

      for (unsigned int i=0; i<triggerNames.size(); i++) {
        if (getTriggerValue(tagTrigBranches[i], event) <= 0) continue;

        fillTagMuonHists(muonHists[i],  triggerPtCuts[i]);
        if (passes_IsoCut() && passes_leptonID())
          fillProbeMuonHists(muonHists[i], trigBranches[i], event, triggerPtCuts[i]);

        fillTagMuonHists(muonHists1[i], triggerPtCuts[i]);
        fillTagMuonHists(muonHists2[i], triggerPtCuts[i]);
        fillTagMuonHists(muonHists3[i], triggerPtCuts[i]);

        fillDenMuonHists(muonHists1[i], triggerPtCuts[i]);
        if (passes_IsoCut() && passes_leptonID())
          fillNumMuonHists(muonHists1[i], triggerPtCuts[i]);

        if (passes_IsoCut()) {
          fillDenMuonHists(muonHists2[i], triggerPtCuts[i]);
          if (passes_leptonID())
            fillNumMuonHists(muonHists2[i], triggerPtCuts[i]);
        }

        if (passes_leptonID()) {
          fillDenMuonHists(muonHists3[i], triggerPtCuts[i]);
          if (passes_IsoCut())
            fillNumMuonHists(muonHists3[i], triggerPtCuts[i]);
        }
      }

      // End of Analysis Code
    }

    // Clean Up
    delete tree;
    file.Close();
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }

  // return
  bmark->Stop("benchmark");

  TFile* outfile = new TFile("hists.root", "RECREATE");

  for (unsigned int i=0; i<triggerNames.size(); i++) {
    TDirectory * dir = (TDirectory*) outfile->mkdir(triggerNames[i].c_str());
    dir->cd();
    TDirectory * dir2 = (TDirectory*) dir->mkdir("trigeff");
    dir2->cd();
    writeEfficiencyPlots(muonHists[i], triggerNames[i], outfile);
    dir2 = (TDirectory*) dir->mkdir("ID+ISO");
    dir2->cd();
    writeEfficiencyPlots(muonHists1[i], triggerNames[i], outfile);
    dir2 = (TDirectory*) dir->mkdir("ID");
    dir2->cd();
    writeEfficiencyPlots(muonHists2[i], triggerNames[i], outfile);
    dir2 = (TDirectory*) dir->mkdir("ISO");
    dir2->cd();
    writeEfficiencyPlots(muonHists3[i], triggerNames[i], outfile);
  }

  h_muonCount->Write();

  outfile->Close();

  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//.........这里部分代码省略.........
	      histos["NB0_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["NB0_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["NB0_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["NB0_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    } else {
	      histos["NBge1_MT_"+skimFilePrefix]->Fill(mu[0].mt,weight);
	      histos["NBge1_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["NBge1_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["NBge1_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["NBge1_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    } if(HT>100.){
	      histos["HTge100_MT_"+skimFilePrefix]->Fill(mu[0].mt,weight);
	      histos["HTge100_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["HTge100_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["HTge100_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["HTge100_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    } if(HT>150.){
	      histos["HTge150_MT_"+skimFilePrefix]->Fill(mu[0].mt,weight);
	      histos["HTge150_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["HTge150_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["HTge150_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["HTge150_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    } if(p.evt_pfmet>50.){
	      histos["METge50_MT_"+skimFilePrefix]->Fill(mu[0].mt,weight);
	      histos["METge50_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["METge50_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["METge50_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["METge50_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    } if(p.evt_pfmet>25.){
	      histos["METge25_MT_"+skimFilePrefix]->Fill(mu[0].mt,weight);
	      histos["METge25_MyMuPt_"+skimFilePrefix]->Fill(mu[0].p4.Pt(),weight);
	      histos["METge25_MyMuEta_"+skimFilePrefix]->Fill(mu[0].p4.Eta(),weight);
	      histos["METge25_MyMuPhi_"+skimFilePrefix]->Fill(mu[0].p4.Phi(),weight);
	      histos["METge25_MyMuCharge_"+skimFilePrefix]->Fill(mucharge,weight);	      
	    }
	  }
	}
      }
      if(!sameEvent(p,c)){
	//finally start new event.
	if(abs(p.id)==13&&p.passes_POG_tightID&&p.miniisoDB<0.2&&fabs(p.dZ)<0.1&&fabs(p.dxyPV)<0.02&&p.dilep_mass>75.&&p.dilep_mass<105&&p.p4.Pt()>20&&fabs(p.p4.Eta())<2.1) {
	  if((p.evt_isRealData&&(p.HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ || p.HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ || p.HLT_IsoMu24_eta2p1 || p.HLT_IsoTkMu24_eta2p1 || p.HLT_IsoMu27 || p.HLT_IsoTkMu27 || p.HLT_IsoMu20 || p.HLT_IsoMu20))||(!p.evt_isRealData)){
	    histos["Z_Pt_"+skimFilePrefix]->Fill(p.dilep_p4.Pt(),weight);
	  }
	}
	mu.clear();
	nleps = 1;
      }
      if(currentrun!=c.evt_run) { cout << "This is run " << c.evt_run << endl; currentrun = c.evt_run; }
      bool skip = false;
      //if(c.evt_run!=251244&&c.evt_run!=251251&&c.evt_run!=251252) skip = true;
      if(abs(c.id)!=13) skip = true;
      else if(c.p4.Pt()<20.) skip = true;
      else if(fabs(c.p4.Eta())>2.1) skip = true;
      else if(!c.passes_POG_looseID) skip = true;
      //else if(fabs(c.dZ)>0.5) skip = true;
      //else if(fabs(c.dxyPV)>0.2) skip = true;
      //else if(c.miniisoDB>0.5) skip = true;
      if(!skip) { mu.push_back(c); }
      p = c;
      resetEvent(c);
    }
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }
  
  // Example Histograms
  for(map<string,TH1F*>::iterator h=histos.begin(); h!=histos.end();++h){
    h->second->SetBinContent(h->second->GetNbinsX(), h->second->GetBinContent(h->second->GetNbinsX() )+ h->second->GetBinContent(h->second->GetNbinsX()+1) );
    h->second->SetBinError(h->second->GetNbinsX(), sqrt(pow(h->second->GetBinError(h->second->GetNbinsX() ),2)+pow(h->second->GetBinError(h->second->GetNbinsX()+1),2) ) );
  }
  for(map<string,TH1F*>::iterator h=histos.begin(); h!=histos.end();++h){
    h->second->SetBinContent(1, h->second->GetBinContent(1 )+ h->second->GetBinContent(0) );
    h->second->SetBinError(1, sqrt(pow(h->second->GetBinError(1 ),2)+pow(h->second->GetBinError(0),2) ) );
  }
  //string filename = "rootfiles/test/Histos_"+skimFilePrefix+".root";
  string filename = "rootfiles/first_20150727/Histos3_"+skimFilePrefix+".root";
  if(DCSonly) filename = "rootfiles/first_20150727/Histos3DCS_"+skimFilePrefix+".root";
  TFile *f = new TFile(filename.c_str(),"RECREATE");
  f->cd();
  for(map<string,TH1F*>::iterator h=    histos.begin(); h!=    histos.end();++h) h->second->Write();
  f->Close();
  cout << "Saved histos in " << f->GetName() << endl;  
  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//.........这里部分代码省略.........
  printf("CREATE MC KINE CUTS\n");
  TObjArray* mcList = new TObjArray(0) ;
  mcList->AddLast(mcKineCuts);
  mcList->AddLast(mcGenCuts);

  printf("CREATE ACCEPTANCE CUTS\n");
  TObjArray* accList = new TObjArray(0) ;
  accList->AddLast(mcAccCuts);

  printf("CREATE RECONSTRUCTION CUTS\n");
  TObjArray* recList = new TObjArray(0) ;
  recList->AddLast(recKineCuts);
  recList->AddLast(recQualityCuts);
  recList->AddLast(recIsPrimaryCuts);

  printf("CREATE PID CUTS\n");
  TObjArray* fPIDCutList = new TObjArray(0) ;
  fPIDCutList->AddLast(cutPID);

  //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
  printf("CREATE INTERFACE AND CUTS\n");
  AliCFManager* man = new AliCFManager() ;

  man->SetNStepEvent(1);
  man->SetEventCutsList(0,evtList);

  man->SetParticleContainer(container);
  man->SetParticleCutsList(0,mcList);
  man->SetParticleCutsList(1,accList);
  man->SetParticleCutsList(2,recList);
  man->SetParticleCutsList(3,fPIDCutList);


  //CREATE THE TASK
  printf("CREATE TASK\n");
  // create the task
  AliCFSingleTrackTask *task = new AliCFSingleTrackTask("AliSingleTrackTask");
  task->SetCFManager(man); //here is set the CF manager
  task->SetQAList(qaList);
  if (readAOD)       task->SetReadAODData() ;
  if (readTPCTracks) task->SetReadTPCTracks();

  //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
  printf("CREATE ANALYSIS MANAGER\n");
  // Make the analysis manager
  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");

  if (useGrid) mgr->SetAnalysisType(AliAnalysisManager::kGridAnalysis);
  else mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);


  AliMCEventHandler*  mcHandler = new AliMCEventHandler();
  mgr->SetMCtruthEventHandler(mcHandler);
 
  AliInputEventHandler* dataHandler ;
  
  if   (readAOD) dataHandler = new AliAODInputHandler();
  else           dataHandler = new AliESDInputHandler();
  mgr->SetInputEventHandler(dataHandler);

  // Create and connect containers for input/output

  //------ input data ------
  AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);

  // ----- output data -----
  
  //slot 0 : default output tree (by default handled by AliAnalysisTaskSE)
  AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root");

  //now comes user's output objects :
  
  // output TH1I for event counting
  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root");
  // output Correction Framework Container (for acceptance & efficiency calculations)
  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root");
  // output QA histograms 
  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("clist0", TList::Class(),AliAnalysisManager::kOutputContainer,"output.root");

  cinput0->SetData(analysisChain);

  mgr->AddTask(task);
  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
  mgr->ConnectOutput(task,0,coutput0);
  mgr->ConnectOutput(task,1,coutput1);
  mgr->ConnectOutput(task,2,coutput2);
  mgr->ConnectOutput(task,3,coutput3);
 
  printf("READY TO RUN\n");
  //RUN !!!
  if (mgr->InitAnalysis()) {
    mgr->PrintStatus();
    mgr->StartAnalysis("local",analysisChain);
  }

  benchmark.Stop("AliSingleTrackTask");
  benchmark.Show("AliSingleTrackTask");

  return kTRUE ;
}

//.........这里部分代码省略.........
				// If any negative yields were found in any decay mode, recalculate the total yield
				if( j==0 && negsFound ) {
					double newYield, newErr;
					newYield = h_bkgtype[i][0]->IntegralAndError( 0, -1, newErr );
					h_yields->SetBinContent(i+1, newYield);
					h_yields->SetBinError(i+1, newErr);
				}
				// Add zeroed histograms to total histograms
				h_bkgtype_sum[i][j]->Add( h_bkgtype[i][j] );
				h_evttype_sum[i][j]->Add( h_evttype[i][j] );
			}
		}
		h_yields_sum->Add( h_yields );

	} // End loop over files in the chain

	cout << "Cutflow yields:                        (yield)  (gen evts)" << endl;

	printf("Total number of events:             %10.2f %9i\n",   yield_total    , yGen_total    );
	if( mySample->IsData() ) {
		printf("Events passing duplicate removal:   %10.2f %9i\n", yield_unique   , yGen_unique   );
		printf("Events passing filters and trigger: %10.2f %9i\n", yield_filter   , yGen_filter   );
	}
	printf("Events with 1st vertex good:        %10.2f %9i\n",   yield_vtx      , yGen_vtx      );
	printf("Events with at least 1 good lepton: %10.2f %9i\n",   yield_1goodlep , yGen_1goodlep );
	printf("Events passing lepton selection:    %10.2f %9i\n",   yield_lepSel   , yGen_lepSel   );

	printf("\nEvents passing 2-lep requirement:   %10.2f %9i\n", yield_2lepCR   , yGen_2lepCR   );
	printf("   Events with veto lepton:         %10.2f %9i\n",   yield_2lepveto , yGen_2lepveto );
	printf("   Events with isolated track:      %10.2f %9i\n",   yield_trkVeto  , yGen_trkVeto  );
	printf("   Events with identified tau:      %10.2f %9i\n\n", yield_tauVeto  , yGen_tauVeto  );

	printf("Events with at least 2 jets:        %10.2f %9i\n",   yield_njets    , yGen_njets    );
	printf("Events with at least 1 b-tag:       %10.2f %9i\n",   yield_1bjet    , yGen_1bjet    );
	printf("Events with MET > 250 GeV:          %10.2f %9i\n",   yield_METcut   , yGen_METcut   );
	printf("Events with MT > 150 GeV:           %10.2f %9i\n",   yield_MTcut    , yGen_MTcut    );
	printf("Events with min dPhi > 0.5:         %10.2f %9i\n",   yield_dPhi     , yGen_dPhi     );
	// printf("Events with chi2 < 10:              %10.2f %9i\n",   yield_chi2     , yGen_chi2     );
	printf("Yield after preselection:           %10.2f %9i\n",   yield_chi2     , yGen_chi2     );

	if ( nEventsChain != nEventsTotal ) {
		cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
	}


	///////////////////////////////////////////////////////////////////////////////
	// Store histograms and clean them up
	TFile* plotfile = new TFile( myAnalysis->GetPlotFileName(), "READ");
	TFile* systfile = new TFile( myAnalysis->GetSystFileName(), "READ");
	TFile* sourcefile;

	// Certain histograms are cumulative across multiple samples. For those histograms, add what the
	// looper has just collected to the cumulative version stored in our output files
	for( int j=0; j<=nVariations; j++ ) {

		if( j==0 ) sourcefile = plotfile;
		else       sourcefile = systfile;

		for( int i=0; i<nSigRegs; i++ ) {

			// Build up cumulative histo of SUSY scan yields
			TH2D* hTemp2 = (TH2D*)sourcefile->Get( h_sigyields[i][j]->GetName() );
			if( hTemp2 != 0 ) h_sigyields[i][j]->Add( hTemp2 );

			// Build up cumulative histo of yields by signal/background type
			TH1D* hTemp = (TH1D*)sourcefile->Get( h_evttype_sum[i][j]->GetName() );
			if( hTemp != 0 ) h_evttype_sum[i][j]->Add( hTemp );
		}
	}
	delete plotfile;
	delete systfile;

	// Take all histograms in histdir and write them to plotfile
	plotfile = new TFile( myAnalysis->GetPlotFileName(), "UPDATE");
	plotfile->cd();
	histdir->GetList()->Write( "", TObject::kOverwrite );
	delete plotfile;

	// Take all histograms in systdir and write them to systfile
	systfile = new TFile( myAnalysis->GetSystFileName(), "UPDATE");
	systfile->cd();
	systdir->GetList()->Write( "", TObject::kOverwrite );
	delete systfile;

	// Cleanup
	zerodir->Close();
	histdir->Close();
	systdir->Close();

	// return
	bmark->Stop("benchmark");
	cout << endl;
	cout << nEventsTotal << " Events Processed" << endl;
	cout << "------------------------------" << endl;
	cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
	cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
	cout << endl;
	delete bmark;
	return 0;
}

int ScanChain( TChain* chain, bool fast = true, int nEvents = -1, string skimFilePrefix = "test") {

  // Benchmark
  TBenchmark *bmark = new TBenchmark();
  bmark->Start("benchmark");

  // Example Histograms
  TDirectory *rootdir = gDirectory->GetDirectory("Rint:");
  TH1F *samplehisto = new TH1F("samplehisto", "Example histogram", 200,0,200);
  samplehisto->SetDirectory(rootdir);

  // Loop over events to Analyze
  unsigned int nEventsTotal = 0;
  unsigned int nEventsChain = chain->GetEntries();
  if( nEvents >= 0 ) nEventsChain = nEvents;
  TObjArray *listOfFiles = chain->GetListOfFiles();
  TIter fileIter(listOfFiles);
  TFile *currentFile = 0;

  // File Loop
  while ( (currentFile = (TFile*)fileIter.Next()) ) {

    // Get File Content
    TFile *file = new TFile( currentFile->GetTitle() );
    TTree *tree = (TTree*)file->Get("t");
    if(fast) TTreeCache::SetLearnEntries(10);
    if(fast) tree->SetCacheSize(128*1024*1024);
    cms3.Init(tree);
    
    // Loop over Events in current file
    if( nEventsTotal >= nEventsChain ) continue;
    unsigned int nEventsTree = tree->GetEntriesFast();
    for( unsigned int event = 0; event < nEventsTree; ++event) {
    
      // Get Event Content
      if( nEventsTotal >= nEventsChain ) continue;
      if(fast) tree->LoadTree(event);
      cms3.GetEntry(event);
      ++nEventsTotal;
    
      // Progress
      StopBabies10012015::progress( nEventsTotal, nEventsChain );

      // Analysis Code

    }
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }
  
  // Example Histograms
  samplehisto->Draw();
  
  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//.........这里部分代码省略.........
	if(tmp<value["MTb"]) value["MTb"] = tmp;  else if(value["MTb"]<0) value["MTb"] = tmp;
	tmp = JetUtil::deltaR(btaggedjets[n],leplv);
	tmp = JetUtil::deltaPhi(btaggedjets[n],metlv);
	if(tmp<value["minDPhiB"]) value["minDPhiB"] = tmp;  else if(value["minDPhiB"]<0) value["minDPhiB"] = tmp;
	for(unsigned int m = n+1; m<btaggedjets.size();++m){
	  temp = btaggedjets[n]+btaggedjets[m]+leplv;
	  if(temp.M()<value["Mlbb"]) value["Mlbb"] = temp.M();  else if(value["Mlbb"]<0) value["Mlbb"] = temp.M();
	  if(temp.Pt()<value["pTlbb"]) value["pTlbb"] = temp.Pt();  else if(value["pTlbb"]<0) value["pTlbb"] = temp.Pt();
	  tmp = JetUtil::deltaR(btaggedjets[n],btaggedjets[m]);
	  if(tmp<value["dRbb"]) value["dRbb"] = tmp;  else if(value["dRbb"]<0) value["dRbb"] = tmp;
	}
      }
      float myleadjpt = -1.;
      int bj1(-1), bj2(-1),bj3(-1);
      for(unsigned int n = 0; n<jetlv.size();++n){
	float tmp = getMT(jetlv[n],metlv);
  	if(tmp<value["MTq"]) value["MTq"] = tmp;  else if(value["MTq"]<0) value["MTq"] = tmp;
  	if(tmp>value["MTqmax"]) value["MTqmax"] = tmp;  else if(value["MTqmax"]<0) value["MTqmax"] = tmp;
	if(jetlv[n].Pt()>250.){ if(tmp>value["MTq_boostmax"]) value["MTq_boostmax"] = tmp;  else if(value["MTq_boostmax"]<0) value["MTq_boostmax"] = tmp; }
	if(jetlv[n].Pt()>300.){ if(tmp<value["MTq_boost300"]) value["MTq_boost300"] = tmp;  else if(value["MTq_boost300"]<0) value["MTq_boost300"] = tmp; }
	if(jetlv[n].Pt()>myleadjpt){ value["MTq_boostLeadJet"] = tmp; myleadjpt = jetlv[n].Pt(); }// else if(value["MTq_boostLeadJet"]<0) value["MTq_boostLeadJet"] = tmp;
	for(unsigned int m = n+1; m<jetlv.size();++m){
	  tmp = getMT(jetlv[n]+jetlv[m],metlv);
	  if(tmp<value["MTqq"]) value["MTqq"] = tmp;  else if(value["MTqq"]<0) value["MTqq"] = tmp;
	}
	if(n<3){
	  tmp = JetUtil::deltaPhi(jetlv[n],metlv);
	  if(tmp<value["minDPhiJ3"]) value["minDPhiJ3"] = tmp;  else if(value["minDPhiJ3"]<0) value["minDPhiJ3"] = tmp;
	}
	float dP1l(-1), dP2l(-1), dP3l(-1);
	if(bj1>0) dP1l = JetUtil::deltaPhi(jetlv[bj1],leplv);
	if(bj2>0) dP2l = JetUtil::deltaPhi(jetlv[bj2],leplv);
	if(bj3>0) dP3l = JetUtil::deltaPhi(jetlv[bj3],leplv);
	if(JetUtil::deltaPhi(jetlv[n],leplv)>dP1l){ bj3 = bj2; bj2 = bj1; bj1 = n;}
	else if(JetUtil::deltaPhi(jetlv[n],leplv)>dP2l){ bj3 = bj2; bj2 = n;}
	else if(JetUtil::deltaPhi(jetlv[n],leplv)>dP3l){ bj3 = n;}
      }

      LorentzVector bjsumlep = jetlv[bj1] + jetlv[bj2] + jetlv[bj3];
      value["M3b"] = bjsumlep.M();

      float prune=0.;
      for(unsigned int n = 0; n<cms3.ak8pfjets_pruned_mass().size();++n) prune += cms3.ak8pfjets_pruned_mass()[n];
      
      //cout << __LINE__<<endl;
      value["sumak8prunedmass"]   = prune;
      value["MT"] = MT;
      value["MET"] = MET;
      value["HT"] = HT;
      if(HT>0) {
	value["METoverSqrtHT"] = MET/TMath::Sqrt(HT);
      }
      value["HTratio"] = HTratio;
      value["dRLepBJet"] = dRLepBJet;
      value["chi2"] = chi2;
      value["NBJets"] = NGBJets;
      value["NJets"] = NGJets;
      value["minDPhi"] = minDPhi;
      value["pTl"] = leplv.Pt();
      value["pTleadj"] = myleadjpt;
      value["pTleadb"] = btaggedjets[leadb].Pt();
      value["pTtrailb"] = btaggedjets[trailb].Pt();
      value["DeltaPhiWl"] = JetUtil::deltaPhi(leplv,metlv+leplv);

      outtree[sampleid]->Fill();


      for(unsigned int i = 0; i<histonames.size(); ++i){
	value[histonames[i] ] = -99;//reset values
      }
    }
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }

  for(unsigned int b = 0; b<3; ++b){
    if(skimFilePrefix!="TTbar"&&b>0) continue;
    outfile[b]->cd();
    outtree[b]->Write();
    outfile[b]->Close();
    cout << "Tree for making BDT saved in " << outfile[b]->GetName() << endl;
  }

  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//_____________________________________________________________________
void bfcZ(const Int_t First,
	 const Int_t Last,
	 const Char_t *Chain="gstar Cy2b tfs -NoHits TpcHitFilter",
	 const Char_t *infile=0,
	 const Double_t zmin=-200, const Double_t zmax=-100,
	 const Char_t *outfile=0,
	 const Char_t *TreeFile=0){ 
  gROOT->LoadMacro("bfc.C");
  bfc(-1,Chain,infile,outfile,TreeFile);
  if (chain->GetOption("TpcHitFilter")) {
    StTpcHitFilterMaker *filtMk = (StTpcHitFilterMaker *) chain->GetMaker("tpc_hit_filter");
    if (filtMk) {
      //      St_tcl_Maker *tclMk= (St_tcl_Maker *) chain->GetMaker("tpc_hits");
      if (zmin < 0 && zmax < 0) filtMk->WestOff();
      if (zmin > 0 && zmax > 0) filtMk->EastOff();
      filtMk->SetZrange(zmin,zmax);
    }
  }
  {
    TDatime t;
    printf ("QAInfo:Run is started at Date/Time %i/%i\n",t.GetDate(),t.GetTime());
    printf ("QAInfo:Run on %s in %s\n",
	    gSystem->HostName(),
	    gSystem->WorkingDirectory());
    printf ("QAInfo: with %s\n", chain->GetCVS());
  }
  // Init the chain and all its makers
  Int_t iTotal = 0, iBad = 0;
  St_XDFFile *xdf_out = 0;
  TBenchmark evnt;
  Int_t iMake = 0, i = First;
  if (Last >= 0) {
    Int_t iInit = chain->Init();
    if (iInit >=  kStEOF) {
      chain->Fatal(iInit,"on init");
      goto END;
    }
    StEvtHddr *hd = (StEvtHddr*)chain->GetDataSet("EvtHddr");
    if (hd) hd->SetRunNumber(-2); // to be sure that InitRun calls at least once
    // skip if any
    if (First > 1) {
      if (chain->GetOption("fzin")) {
	geant = (St_geant_Maker *) chain->GetMaker("geant");
	if (geant) {
	  if (geant->IsActive()) geant->Skip(First-1);
	}
      }
      else {
	StIOMaker *inpMk      = (StIOMaker *)      chain->GetMaker("inputStream");
	if (inpMk) {printf ("Skip %i Events\n",First-1);inpMk->Skip(First-1);}
      }
    }
  }
  xdf_out = chain->GetXdfOut();
  if (chain->GetOption("Event")) evMk  = (StEventMaker   *) chain->GetMaker("StEventMaker");  
  treeMk = chain->GetMaker("OutputStream");
 EventLoop: if (i <= Last && iMake != kStEOF && iMake != kStFatal) {
   evnt.Reset();
   evnt.Start("QAInfo:");
   chain->Clear();
   iMake = chain->Make(i);
   if (iMake <kStEOF) {
     if (xdf_out){
       St_DataSet *dstSet = chain->GetInputDS("dst");
       if (dstSet) xdf_out->NextEventPut(dstSet); // xdf output
     }
     iTotal++;
     if (treeMk && iMake == kStErr) {treeMk->Make(i); iBad++;}
     //    gSystem->Exec("ps ux");
     evnt.Stop("QAInfo:");
     //  evnt.Show("QAInfo:");
     printf ("QAInfo: Done with Event [no. %d/run %d/evt. %d/Date.Time %d.%d/sta %d] Real Time = %10.2f seconds Cpu Time =  %10.2f seconds \n",
	     i,chain->GetRunNumber(),chain->GetEventNumber(),chain->GetDate(), chain->GetTime(),
	     iMake,evnt.GetRealTime("QAInfo:"),evnt.GetCpuTime("QAInfo:"));
   }
   i++; goto EventLoop;
 }
 END:
  fflush(stdout);
  printf ("QAInfo:Run completed ");
  gSystem->Exec("date");
  if (evMk) Event = (StEvent *) chain->GetInputDS("StEvent");
  {
    TDatime t;
    printf ("\nQAInfo:Run is finished at Date/Time %i/%i; Total events processed :%i and not completed: %i\n",
	    t.GetDate(),t.GetTime(),iTotal,iBad);
  }
}

//.........这里部分代码省略.........
					  h_met_ll_tar_njets0->Fill(zmet.pfmet()*TMath::Cos( zmet.pfmetphi() ), zmet.pfmet()*TMath::Sin( zmet.pfmetphi() ), weight_mc);
					}
				  else if(zmet.njets()>=2)
					{
					  h_met_mumu_tar_njets2->Fill(zmet.pfmet(),weight_mc);
					  h_met_ll_tar_njets2->Fill(zmet.pfmet()*TMath::Cos( zmet.pfmetphi() ), zmet.pfmet()*TMath::Sin( zmet.pfmetphi() ), weight_mc);
					}
				}

			  if (zmet.lep3().Pt() <= 10.  &&  zmet.njets()>=2  &&  zmet.nbcsvm()==0 && (zmet.mjj() > 70. && zmet.mjj() < 110.) )
   				{h_mll_mumu_tar->Fill(zmet.dilmass(),weight_mc);}
			}
 
		  else if (zmet.leptype()==2)
			{
			  // h_mll_mumu_inc->Fill(zmet.dilmass(),weight_mc);
			
			  if(zmet.dilmass() > 81. && zmet.dilmass() < 101. /* && zmet.lep3().Pt() <= 10.*/)
				{
				  h_met_emu_inc->Fill(zmet.pfmet(),weight_mc);
				  fillHist(h_sumet_emu_inc, zmet.pfsumet(),weight_mc);
				  h_met_emu2_inc->Fill(zmet.pfmet()*TMath::Cos( zmet.pfmetphi() ), zmet.pfmet()*TMath::Sin( zmet.pfmetphi() ), weight_mc);

				  if(zmet.njets()==0)
					{
					  h_met_emu_tar_njets0->Fill(zmet.pfmet(),weight_mc);
					  h_met_emu2_tar_njets0->Fill(zmet.pfmet()*TMath::Cos( zmet.pfmetphi() ), zmet.pfmet()*TMath::Sin( zmet.pfmetphi() ), weight_mc);
					}
				  else if(zmet.njets()>=2)
					{
					  h_met_emu_tar_njets2->Fill(zmet.pfmet(),weight_mc);
					  h_met_emu2_tar_njets2->Fill(zmet.pfmet()*TMath::Cos( zmet.pfmetphi() ), zmet.pfmet()*TMath::Sin( zmet.pfmetphi() ), weight_mc);
					}
				}

			  // if (zmet.lep3().Pt() <= 10.  &&  zmet.njets()>=2  &&  zmet.nbcsvm()==0 && (zmet.mjj() > 70. && zmet.mjj() < 110.) )
			  //	{h_mll_mumu_tar->Fill(zmet.dilmass(),weight_mc);}
			}
	}

	
  

 //---------------------------------------------------------------------------------------------------------
    }    //end events loop.
 
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }


  //Write to file
  TFile *OutputFile = new TFile(Form("%s.root", skimFilePrefix.c_str()),"recreate");
  OutputFile->cd();
  h_mll_ee_inc->Write();
  h_mll_mumu_inc->Write();
  h_mll_ee_tar->Write();
  h_mll_mumu_tar->Write();
  h_met_ee_inc->Write();
  h_met_mumu_inc->Write();
  h_met_ll_inc->Write();
  h_met_emu_inc->Write();
  h_met_emu2_inc->Write();
  h_met_ee_tar_njets0->Write();
  h_met_mumu_tar_njets0->Write();
  h_met_ll_tar_njets0->Write();
  h_met_emu_tar_njets0->Write();
  h_met_emu2_tar_njets0->Write();
  h_met_ee_tar_njets2->Write();
  h_met_mumu_tar_njets2->Write();
  h_met_ll_tar_njets2->Write();
  h_met_emu_tar_njets2->Write();
  h_met_emu2_tar_njets2->Write();
  h_nvtx_scaled->Write();
  h_nvtx_unscaled->Write();
  h_phi->Write();
  h_sumet_ee_inc->Write();
  h_sumet_mumu_inc->Write();
  h_sumet_emu_inc->Write();

  OutputFile->Close();


  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;

  return 0;
}

//.........这里部分代码省略.........
      }
      if(version==13){
	if(mlb){
	  if(met>350) mybin = 6; else if(met>300) mybin = 5; else mybin = 4;
	}
	else if(met>425) mybin = 3; else if(met>325) mybin = 2; else mybin = 1;
      }
      if(version==14){
	if(m2w){
	  if(met>500) mybin = 8; else if(met>400) mybin = 7; else if(met>350) mybin = 6; else if(met>300) mybin = 5; else mybin = 4;
	}
	else if(met>350) mybin = 3; else if(met>300) mybin = 2; else mybin = 1;
      }
      if(version==15){
	if(met>500) mybin = 5; else if(met>400) mybin = 4; else if(met>350) mybin = 3; else if(met>300) mybin = 2; else mybin = 1;
      }
      //cout << __LINE__ << " mybin " << mybin << " " << samplename << endl;
      //cout << "SRyield_" + samplename << endl;
      histos["SRyield_" + samplename]->Fill(mybin,weight);

    }
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }

  //do all postrunning modifications here - maybe need some additional code if running over additional types of babies
  //i.e.if we have additional babies, run extra, hadd files and do further modifications (note you might need to store additional histos)
  for(unsigned int b = 0; b<2; ++b){
    string samplename = skimFilePrefix;
    TString helper = TString(samplename);
    if(skimFilePrefix!="TTbar"&&b>0) continue;
    if(skimFilePrefix=="TTbar"&&b==0) samplename = "TTbar1l";
    if(skimFilePrefix=="TTbar"&&b==1) samplename = "TTbar2l";
    if(skimFilePrefix=="TTbar"&&b==2) samplename = "TTbarH";
    string mapname;
    for(int i = 1; i<=histos["SRyield_"+samplename]->GetNbinsX();++i){
      if(histos["SRyield_"+samplename]->Integral()>0) {
	if(histos["SRyield_"+samplename]->GetBinContent(i)>0){
	  histos["StatUnc_"+samplename]->SetBinContent(i, histos["SRyield_"+samplename]->GetBinError(i)/histos["SRyield_"+samplename]->GetBinContent(i));
	}
	else histos["StatUnc_"+samplename]->SetBinContent(i, -999);
	histos["SystUnc_"+samplename]->SetBinContent(i, 0.3);
      }
    }
  }  

  //save all relevant histograms into some rootfiles
  for(unsigned int b = 0; b<2; ++b){
    string samplename = skimFilePrefix;
    TString helper = TString(samplename);
    if(skimFilePrefix!="TTbar"&&b>0) continue;
    if(skimFilePrefix=="TTbar"&&b==0) samplename = "TTbar1l";
    if(skimFilePrefix=="TTbar"&&b==1) samplename = "TTbar2l";
    if(skimFilePrefix=="TTbar"&&b==2) continue;
    string filename = "inputfiles/"+samplename+".root";
    if(version== 1) filename = "inputfiles/v1/"+samplename+".root";
    if(version== 2) filename = "inputfiles/v2/"+samplename+".root";
    if(version== 3) filename = "inputfiles/v3/"+samplename+".root";
    if(version== 4) filename = "inputfiles/v4/"+samplename+".root";
    if(version== 5) filename = "inputfiles/v5/"+samplename+".root";
    if(version== 6) filename = "inputfiles/v6/"+samplename+".root";
    if(version== 7) filename = "inputfiles/v7/"+samplename+".root";
    if(version== 8) filename = "inputfiles/v8/"+samplename+".root";
    if(version== 9) filename = "inputfiles/v9/"+samplename+".root";
    if(version==10) filename = "inputfiles/v10/"+samplename+".root";
    if(version==11) filename = "inputfiles/v11/"+samplename+".root";
    if(version==12) filename = "inputfiles/v12/"+samplename+".root";
    if(version==13) filename = "inputfiles/v13/"+samplename+".root";
    if(version==14) filename = "inputfiles/v14/"+samplename+".root";
    if(version==15) filename = "inputfiles/v15/"+samplename+".root";
    if(version==16) filename = "inputfiles/v16/"+samplename+".root";
    if(version==17) filename = "inputfiles/v17/"+samplename+".root";
    TFile *f = new TFile(filename.c_str(),"RECREATE");
    f->cd();
    string mapname;
    for(unsigned int i = 0; i<histonames.size();++i){
      mapname = histonames[i] + "_" + samplename;
      if(storeh) histos[mapname]->Write();
    }
    f->Close();
    cout << "Saved histos in " << f->GetName() << endl;
  }

  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}

//.........这里部分代码省略.........
      if(fast) tree->LoadTree(event);
      cms3.GetEntry(event);
      ++nEventsTotal;
    
      // Progress
      CMS3::progress( nEventsTotal, nEventsChain );
      
      float weight = cms3.scale1fb()*2.26;
      if(event==0) cout << "weight " << weight << " nEvents " << cms3.nEvents() << " filename " << currentFile->GetTitle() << endl;
      
      //ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float> > metlv;
      //metlv.SetPxPyPzE(pfmet()*TMath::Cos(pfmet_phi()),pfmet()*TMath::Sin(pfmet_phi()),0.,pfmet());
      int NSLeps = 0;
      int NAddVetoLeps = 0;
      if(lep1_is_mu()){
	if(lep1_pt()>20&&fabs(lep1_eta())<2.4) {++NSLeps;}
      } else if (lep1_is_el()){
	if(lep1_pt()>20&&fabs(lep1_eta())<1.4442) {++NSLeps; }
      } if(lep2_is_mu()){
	if(lep2_pt()>20&&fabs(lep2_eta())<2.4) {++NSLeps;}
      } else if (lep2_is_el()){
	if(lep2_pt()>20&&fabs(lep2_eta())<1.4442) {++NSLeps; }
      }
      if(lep2_is_mu()){
	if(lep2_pt()>10&&fabs(lep2_eta())<2.4) {++NAddVetoLeps;}
      } else if (lep2_is_el()){
	if(lep2_pt()>10&&fabs(lep2_eta())<2.4) {++NAddVetoLeps; }
      }


      if(nvtxs()<0)               continue;
      if(ngoodleps()!=1)          continue;
      if(nvetoleps()!=1)          continue;
      if(!PassTrackVeto_v3())     continue;
      if(!PassTauVeto())          continue;
      if(ngoodjets()<2)           continue;
      if(ngoodbtags()<1)          continue;
      if(pfmet()<100)             continue;
      if(mt_met_lep()< 50)        continue;
      //if(pfmet()<250)             continue;
      //if(mt_met_lep()<150)        continue;
      //if(mindphi_met_j1_j2()<0.8) continue;

      string samplename = skimFilePrefix;
      if(skimFilePrefix=="TTbar"){
	if(is2lep() ) samplename = "TTbar2l";
	else if(is1lepFromTop() ) samplename = "TTbar1l";
	else continue;
      }

      if(ngoodjets()>=2 && mindphi_met_j1_j2()>0.8 && pfmet()>250) histos["MT_"+samplename]->Fill(mt_met_lep(),weight);
      if(ngoodjets()>=2 && mindphi_met_j1_j2()>0.8 && pfmet()>350) histos["MT_harder_"+samplename]->Fill(mt_met_lep(),weight);
      if(ngoodjets()>=3 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150) histos["MT2W_"+samplename]->Fill(MT2W(),weight);
      if(ngoodjets()>=4 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150) histos["MT2W_4j_"+samplename]->Fill(MT2W(),weight);
      if(ngoodjets()>=2 && mindphi_met_j1_j2()>0.8 && mt_met_lep()>150) histos["MET_"+samplename]->Fill(pfmet(),weight);
      if(ngoodjets()>=3 && mindphi_met_j1_j2()>0.8 && mt_met_lep()>150) histos["MET_3j_"+samplename]->Fill(pfmet(),weight);
      if(ngoodjets()>=4 && mindphi_met_j1_j2()>0.8 && mt_met_lep()>150) histos["MET_4j_"+samplename]->Fill(pfmet(),weight);
      if(ngoodjets()>=3 && mindphi_met_j1_j2()>0.8 && mt_met_lep()>150 && MT2W()>200) histos["MET_hMT2W_"+samplename]->Fill(pfmet(),weight);
      if(ngoodjets()>=4 && mindphi_met_j1_j2()>0.8 && mt_met_lep()>150 && MT2W()>200) histos["MET_hMT2W_4j_"+samplename]->Fill(pfmet(),weight);
      if(ngoodjets()>=2 && pfmet()>250 && mt_met_lep()>150) histos["MinDPhi_"+samplename]->Fill(mindphi_met_j1_j2(),weight);
      if(ngoodjets()>=2 && pfmet()>350 && mt_met_lep()>150) histos["MinDPhi_harder_"+samplename]->Fill(mindphi_met_j1_j2(),weight);
      if(ngoodjets()>=3 && pfmet()>250 && mt_met_lep()>150) histos["MinDPhi_3j_"+samplename]->Fill(mindphi_met_j1_j2(),weight);
      if(ngoodjets()>=4 && pfmet()>250 && mt_met_lep()>150) histos["MinDPhi_4j_"+samplename]->Fill(mindphi_met_j1_j2(),weight);
      if(ngoodjets()>=2 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150) histos["NJets_"+samplename]->Fill(ngoodjets(),weight);
      if(ngoodjets()==2 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150 && topnessMod()>6.4) histos["NJets_hard_"+samplename]->Fill(ngoodjets(),weight);
      if(ngoodjets()>=3 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150 && MT2W()>200)       histos["NJets_hard_"+samplename]->Fill(ngoodjets(),weight);
      if(ngoodjets()==2 && mindphi_met_j1_j2()>0.8 && pfmet()>350 && mt_met_lep()>150 && topnessMod()>6.4) histos["NJets_harder_"+samplename]->Fill(ngoodjets(),weight);
      if(ngoodjets()>=3 && mindphi_met_j1_j2()>0.8 && pfmet()>350 && mt_met_lep()>150 && MT2W()>200)       histos["NJets_harder_"+samplename]->Fill(ngoodjets(),weight);
      if(ngoodjets()>=2 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150) histos["tMod_"+samplename]->Fill(topnessMod(),weight);
      if(ngoodjets()==2 && mindphi_met_j1_j2()>0.8 && pfmet()>250 && mt_met_lep()>150) histos["tMod_2j_"+samplename]->Fill(topnessMod(),weight);

    }//event loop
  
    // Clean Up
    delete tree;
    file->Close();
    delete file;
  }//file loop
  if ( nEventsChain != nEventsTotal ) {
    cout << Form( "ERROR: number of events from files (%d) is not equal to total number of events (%d)", nEventsChain, nEventsTotal ) << endl;
  }
  
  string filename = "rootfiles/PASfigure.root";
  TFile *f = new TFile(filename.c_str(),"update");
  f->cd();
  for(map<string,TH1F*>::iterator h=    histos.begin(); h!=    histos.end();++h) h->second->Write();
  f->Close();
  cout << "Saved histos in " << f->GetName() << endl;
  delete f;
  // return
  bmark->Stop("benchmark");
  cout << endl;
  cout << nEventsTotal << " Events Processed" << endl;
  cout << "------------------------------" << endl;
  cout << "CPU  Time:	" << Form( "%.01f", bmark->GetCpuTime("benchmark")  ) << endl;
  cout << "Real Time:	" << Form( "%.01f", bmark->GetRealTime("benchmark") ) << endl;
  cout << endl;
  delete bmark;
  return 0;
}