C++ TChain::GetListOfFiles方法代码示例

TChain * GetAnalysisChain(const char * incollection){
  // Builds a chain of esd files
  // incollection can be
  // - a single root file
  // - an xml collection of files on alien
  // - a ASCII containing a list of local root files
  TChain* analysisChain = 0;
  // chain
  analysisChain = new TChain("esdTree");
  if (TString(incollection).Contains(".root")){
    analysisChain->Add(incollection);
  }
  else if (TString(incollection).Contains("xml")){
    TGrid::Connect("alien://");
    TAlienCollection * coll = TAlienCollection::Open (incollection);
    while(coll->Next()){
      analysisChain->Add(TString("alien://")+coll->GetLFN());
    }
  } else {
    ifstream file_collect(incollection);
    TString line;
    while (line.ReadLine(file_collect) ) {
      analysisChain->Add(line.Data());
    }
  }
  analysisChain->GetListOfFiles()->Print();

  return analysisChain;
}

void formChain300(){
  cout << "Begin loading..." << endl;
  const TString treedirectory = "demo/T_weights";
  TChain* InputChain = new TChain(treedirectory);

  InputChain->Add("rfio:/castor/cern.ch/user/k/kreis/CUSusy/QCD/T_weight_QCD_Pt300_*.root");
  
  InputChain->Merge("/cu1/kreis/QCD/T_weights_QCD_Pt300.root");

  cout << "Loaded " << InputChain->GetListOfFiles()->GetEntries() << " files" << endl;

 return;

}

//__________________________________________________________
TChain* CreateChainCAF(Int_t nFilesMax, TFileCollection* coll, char* treeName)
{
  TIter iter(coll->GetList());
  
  TChain* target = new TChain(treeName);
  
  Int_t nFiles = 0;
  
  TFileInfo* fileInfo = 0;
  while ((fileInfo = dynamic_cast<TFileInfo*> (iter())) && (nFiles<nFilesMax || nFilesMax == 0)){
    if (fileInfo->GetFirstUrl()) {
      target->Add(fileInfo->GetFirstUrl()->GetUrl());
      nFiles++;
    }
  }
  
  Printf("Added %d files to chain", target->GetListOfFiles()->GetEntries());
  
  return target;
}

//.........这里部分代码省略.........
	Float_t diMuonF_FittedVtx_vx;
	Float_t diMuonF_FittedVtx_vy;
	Float_t diMuonF_FittedVtx_vz;

	//============= Counters ===========================//

	vector<double> FakesPerSample;

	Int_t ev_all = 0;
	Int_t ev_isVtxOK = 0;
	Int_t ev_is2MuJets = 0;
	Int_t ev_is2DiMuons = 0;
	Int_t ev_is2DiMuonsFittedVtxOK = 0;
	Int_t ev_isPixelHitOK = 0;
	Int_t ev_is2DiMuonsDzOK_FittedVtx = 0;
	Int_t ev_is2DiMuonsMassOK_FittedVtx = 0;
	Int_t ev_is2DiMuonsIsoTkOK_FittedVtx = 0;
	Int_t ev_isDiMuonHLTFired = 0;

	Int_t c1genm = 0;
	Int_t c2genm = 0;
	Int_t c3genm = 0;
	Int_t c4genm = 0;
	Int_t ev_4gmlxylzcut = 0;

	Int_t c1recm = 0;
	Int_t c2recm = 0;
	Int_t c3recm = 0;
	Int_t c4recm = 0;

	Int_t uncuttableFakeCounter = 0;
	//}}}

	TObjArray *fileElements=chain->GetListOfFiles();
	TIter next(fileElements);
	TChainElement *chEl=0;

	while ((chEl=(TChainElement*)next())) {
		if (verbose) std::cout << "running on file " << chEl->GetTitle() << std::endl;
		TFile* myfile = new TFile(dirname + chEl->GetTitle());
		if (!myfile) {
			if (verbose) std::cout << "File " << chEl->GetTitle() << " does not exist" << std::endl;
			continue;
		}

		if (verbose) std::cout << "Loading directory cutFlowAnalyzer" << std::endl;
		myfile->cd("cutFlowAnalyzer");

		TTree *t = (TTree*)myfile->Get("cutFlowAnalyzer/Events");
		if (!t) {
			if (verbose) std::cout << "Tree cutFlowAnalyzer/Events does not exist" << std::endl;
			continue;
		}

		if (verbose) cout<<"  Events  "<<t->GetEntries()<<endl;

		//Pull variables from nTuple
		//{{{
		// Event info
		t->SetBranchAddress("event", &event);
		t->SetBranchAddress("run",   &run);
		t->SetBranchAddress("lumi",  &lumi);

		// GEN Level Selectors
		t->SetBranchAddress("is4GenMu",    &is4GenMu);
		t->SetBranchAddress("is1GenMu17",  &is1GenMu17);

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

	int yGen_total = 0;
	int yGen_unique = 0;
	int yGen_filter = 0;
	int yGen_vtx = 0;
	int yGen_1goodlep = 0;
	int yGen_lepSel = 0;
	int yGen_2lepveto = 0;
	int yGen_trkVeto = 0;
	int yGen_tauVeto = 0;
	int yGen_2lepCR = 0;
	int yGen_njets = 0;
	int yGen_1bjet = 0;
	int yGen_METcut = 0;
	int yGen_MTcut = 0;
	int yGen_dPhi = 0;
	int yGen_chi2 = 0;

	////////////////////////////////////////////////////////////////////
	// Set up data-specific filters

	if( mySample->IsData() ) {
		set_goodrun_file_json( "reference-files/Cert_271036-284044_13TeV_23Sep2016ReReco_Collisions16_JSON.txt" );
		duplicate_removal::clear_list();
	}


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

	// 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( currentFile->GetTitle() );
		TString filename = file.GetName();
		TTree *tree = (TTree*)file.Get("t");
		if(fast) TTreeCache::SetLearnEntries(10);
		if(fast) tree->SetCacheSize(128*1024*1024);
		cms3.Init(tree);

		// Load event weight histograms
		TH2F* hNEvts = (TH2F*)file.Get("histNEvts");
		TH3D* hCounterSMS = (TH3D*)file.Get("h_counterSMS");
		TH1D* hCounter = (TH1D*)file.Get("h_counter");
		myHelper.Setup( isFastsim, hCounter, hNEvts, hCounterSMS );

		// Reset zeroing histograms
		for( int i=0; i<nSigRegs; i++ ) {
			for( int j=0; j<=nVariations; j++ ) {
				h_bkgtype[i][j]->Reset();
				h_evttype[i][j]->Reset();
			}
		}
		h_yields->Reset();

		// Loop over Events in current file
		if( nEventsTotal >= nEventsChain ) continue;
		unsigned int nEventsTree = tree->GetEntriesFast();

// float scan(unsigned int njetsCut=0, int btagCut=9999, TString tag="", float manualScale=-1.0){
int scan(unsigned int njetsLow=0, unsigned int njetsHigh=9999, int btagCut=9999, int metLow=0, int metHigh=9999, int htLow=0, int htHigh=9999, TString tag=""){
    //njetsLow, njetsHigh, btagCut, metLow, metHigh, htLow, htHigh, tag, manualScale
    //

    // TH1F* h1D_dummy_data = new TH1F("dummy", "dummyhisto", 10, 0, 10); 

    TH1F* h1D_njets_data = new TH1F("njets", "", 15, 0, 15); 
    TH1F* h1D_ht_data = new TH1F("ht", "", 20, 0, 600); 
    TH1F* h1D_met_data = new TH1F("met", "", 20, 0, 300); 
    TH1F* h1D_mt_data = new TH1F("mt", "", 20, 0, 400); 
    TH1F* h1D_mtW_data = new TH1F("mtW", "", 40, 0, 200); 
    TH1F* h1D_zmass_data = new TH1F("zmass", "", 42, 70, 112); 
    TH1F* h1D_lepeta_data = new TH1F("lepeta", "", 50, -3.0, 3.0); 
    TH1F* h1D_leppt_data = new TH1F("leppt", "", 50, 0, 150); 
    TH1F* h1D_Wleppt_data = new TH1F("Wleppt", "", 30, 0, 150); 

    TH1F* h1D_nbtags_data = new TH1F("nbtags", "", 15, 0, 15); 
    TH1F* h1D_btagval_data = new TH1F("btagval", "", 20, 0, 1); 
    TH1F* h1D_ptZ_data = new TH1F("ptZ", "", 25, 0, 400);
    TH1F* h1D_st_data = new TH1F("st", "", 25, 100, 800);
    TH1F* h1D_minRLeadingJet_data = new TH1F("minRLeadingJet", "",10,0,5);
    TH1F* h1D_ptj1_data = new TH1F("ptj1", "", 25, 0, 300);
    TH1F* h1D_ptj2_data = new TH1F("ptj2", "", 25, 0, 300);
    TH1F* h1D_massZj1_data = new TH1F("massZj1", "", 25, 100, 600);
    TH1F* h1D_massZj2_data = new TH1F("massZj2", "", 25, 100, 600);
    TH1F* h1D_ptjj_data = new TH1F("ptjj", "", 25, 0, 400);
    TH1F* h1D_massjj_data = new TH1F("massjj", "", 25, 0, 600);

    TH1F* h1D_mtWeemu_data = new TH1F("mtWeem", "", 40, 0, 200); 
    TH1F* h1D_mtWmumue_data = new TH1F("mtWmme", "", 40, 0, 200); 
    TH1F* h1D_mtWeee_data = new TH1F("mtWeee", "", 40, 0, 200); 
    TH1F* h1D_mtWmumumu_data = new TH1F("mtWmmm", "", 40, 0, 200); 
    TH1F* h1D_lepetae_data = new TH1F("lepetae", "", 50, -3.0, 3.0); 
    TH1F* h1D_lepetamu_data = new TH1F("lepetamu", "", 50, -3.0, 3.0); 

    float luminosity = 19.407;
    float ptCut = 20; // threshold for "high-pT" in AN
    float zmassCut = 15; // line 212 of AN
    float btagDiscriminant = 0.679; // line 230 of AN

    map<int, vector<int> > runLumi;
    cout << ">>> tag is " << tag << endl;
    cout << ">>> excluding events with #btaggedjets>= " << btagCut << endl;
    cout << ">>> considering events with njets in [" << njetsLow << "," << njetsHigh << "]" << endl;
    cout << ">>> considering events with met in [" << metLow << "," << metHigh << "]" << endl;
    cout << ">>> considering events with ht in [" << htLow << "," << htHigh << "]" << endl;
    // cout << ">>> manualScale is " << manualScale << endl;
    // DATADATA
    {
        initCounter();

        clear_seen();

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

        ch->Add("/home/users/namin/sandbox/condorTest/output/baby_2012ABCD.root");

        int nEventsTotal = ch->GetEntries();
        int nEventsSoFar = 0;
        int nGoodEvents = 0;
        float scale = 1.0;

        TFile *currentFile = 0;
        TObjArray *listOfFiles = ch->GetListOfFiles();
        TIter fileIter(listOfFiles);

        // File Loop
        while ( (currentFile = (TFile*)fileIter.Next()) ) { 
            TFile *file = new TFile( currentFile->GetTitle() );
            TTree *tree = (TTree*)file->Get("tree");
            cms2.Init(tree);
            // Set Good Run List
            if(evt_isRealData()) set_goodrun_file("final_19p49fb_cms2.txt");

            TString filename(currentFile->GetTitle());

            // Loop over Events in current file
            unsigned int nEventsTree = tree->GetEntriesFast();
            for( unsigned int event = 0; event < nEventsTree; ++event) {

                // if(event > 30000) break;

                // Get Event Content
                cms2.GetEntry(event);
                nEventsSoFar++;

                // Progress
                CMS2::progress( nEventsSoFar, nEventsTotal );

                // Select Good Runs
                if( evt_isRealData() && !goodrun( evt_run(), evt_lumiBlock() ) ) continue;

                if(evt_isRealData()){
                    DorkyEventIdentifier id = { evt_run(), evt_event(), evt_lumiBlock() };
                    if ( is_duplicate(id) ){
                        continue;
                    }
                }

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

void ROC_Occupancy() {
  gStyle->SetPalette(1);

  TH2D* l1 = new TH2D("l1","Layer 1;Modules along Z;Ladders",72, -4.5, 4.5,  42, -10.5,  10.5);
  TH2D* l2 = new TH2D("l2","Layer 2;Modules along Z;Ladders",72, -4.5, 4.5,  66, -16.5,  16.5);
  TH2D* l3 = new TH2D("l3","Layer 3;Modules along Z;Ladders",72, -4.5, 4.5,  90, -22.5,  22.5);
  TH2D* fpixI = new TH2D("fpixI","FPix Inner Shells (+x);;Blades", 72, -4.5, 4.5, 144,   0.5,  12.5);
  TH2D* fpixO = new TH2D("fpixO","FPix Outer Shells (-X);;Blades", 72, -4.5, 4.5, 144, -12.5,  -0.5);
  fpixI->GetXaxis()->SetBinLabel(1, "Disk-2 Pnl2");
  fpixI->GetXaxis()->SetBinLabel(9, "Disk-2 Pnl1");
  fpixI->GetXaxis()->SetBinLabel(19, "Disk-1 Pnl2");
  fpixI->GetXaxis()->SetBinLabel(27, "Disk-1 Pnl1");
  fpixI->GetXaxis()->SetBinLabel(41, "Disk+1 Pnl1");
  fpixI->GetXaxis()->SetBinLabel(49, "Disk+1 Pnl2");
  fpixI->GetXaxis()->SetBinLabel(59, "Disk+2 Pnl1");
  fpixI->GetXaxis()->SetBinLabel(67, "Disk+2 Pnl2");
  fpixI->GetXaxis()->LabelsOption("d");
  fpixO->GetXaxis()->SetBinLabel(1, "Disk-2 Pnl2");
  fpixO->GetXaxis()->SetBinLabel(9, "Disk-2 Pnl1");
  fpixO->GetXaxis()->SetBinLabel(19, "Disk-1 Pnl2");
  fpixO->GetXaxis()->SetBinLabel(27, "Disk-1 Pnl1");
  fpixO->GetXaxis()->SetBinLabel(41, "Disk+1 Pnl1");
  fpixO->GetXaxis()->SetBinLabel(49, "Disk+1 Pnl2");
  fpixO->GetXaxis()->SetBinLabel(59, "Disk+2 Pnl1");
  fpixO->GetXaxis()->SetBinLabel(67, "Disk+2 Pnl2");
  fpixO->GetXaxis()->LabelsOption("d");

  TreeReader tr;
  TChain *filechain = new TChain("filechain");
  filechain->Add("Ntuple.root");

  
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);
  TCanvas c("c","c",1200,800);
  c.Divide(3,2);

  TObjArray* files=filechain->GetListOfFiles();
  for (int nf=0; nf<files->GetEntries(); ++nf) {
    TFile* file = TFile::Open(files->At(nf)->GetTitle());
    tr.readtrees(*file);
    
    for (Long64_t i=0; i<tr.nclu(); i++) {
      tr.clu_read(i);
      Cluster clu = tr.clu();
      EventData e = tr.clu_evt();
      
      int clu_sdpx = ((clu.mod_on.disk>0) ? 1 : -1) * (2 * (abs(clu.mod_on.disk) - 1) + clu.mod_on.panel);
      // Roc BinX number for ROC map plots (clust branch)
      int clu_roc_binx = NOVAL_I;
      if (clu.mod_on.det==0) {
	for (int j=1;j<=8;j++) if (clu.y>=((8-j)*52.0)&&clu.y<((9-j)*52.0))
          clu_roc_binx = (clu.mod_on.module+4)*8 + j;
      } else if (clu.mod_on.det==1) {
        // Roc is left (0) or right (1) on the ROC map plot (+Z side)
        int binselx = (clu.mod_on.panel==1&&(clu.mod_on.module==1||clu.mod_on.module==4)) ? (clu.mod_on.module==1)
          : ((clu.mod_on.panel==1&&clu.x<80.0)||(clu.mod_on.panel==2&&clu.x>=80.0));
        // Gives the Roc location inside a panel (0 to 5 on +Z side)
        int nperpan = 2 * clu.mod_on.module + clu.mod_on.panel - 1 + binselx;
        clu_roc_binx = ((clu.mod_on.disk>0) ? nperpan : 9 - nperpan) + (clu_sdpx + 4) * 8 
          - 2 * ((abs(clu.mod_on.disk)==1) ? clu.mod_on.disk : 0);
      }
      // Roc BinY number for ROC map plots (c branch)
      int clu_roc_biny = NOVAL_I; 
      if (clu.mod_on.det==0) {
        // Roc is in bottom (0) or top bin (1) inside a ladder on th ROC map plot
        int binsely = ((clu.mod_on.half==1&&((clu.mod_on.ladder<0&&clu.mod_on.ladder%2==0)||clu.mod_on.ladder%2==1))
        	       ||(clu.mod_on.half==0&&((clu.mod_on.ladder<0 &&((clu.mod_on.ladder%2==-1&&clu.x<80.0)
        								   ||(clu.mod_on.ladder%2==0&&clu.x>=80.0)))
        					 ||(clu.mod_on.ladder>0 &&((clu.mod_on.ladder%2==0&&clu.x<80.0)
        								     ||(clu.mod_on.ladder%2==1&&clu.x>=80.0))))));
        clu_roc_biny = (clu.mod_on.layer * 6 + clu.mod_on.ladder + 4) * 2 + 1 + binsely;
      } else if (clu.mod_on.det==1) {
        // Gives the number of ROCs along ly
        int nrocly = clu.mod_on.module + clu.mod_on.panel;
        for (int j=0; j<nrocly; j++) {
          // ROC number = nrocly - 1 - j for + LX and nrocly + j for -LX.
          int k = (clu.mod_on.disk<0) ? j : nrocly - 1 - j;
          if (clu.y>=(k*52.0)&&clu.y<((k+1)*52.0))
            clu_roc_biny = 6 - nrocly + 2 * j + ((clu.mod_on.blade>0) ? clu.mod_on.blade-1 : clu.mod_on.blade + 12)*12 + 1;
        }
      }
      
      if (clu.mod_on.layer==1) l1->SetBinContent(clu_roc_binx,clu_roc_biny, l1->GetBinContent(clu_roc_binx,clu_roc_biny)+1);
      else if (clu.mod_on.layer==2) l2->SetBinContent(clu_roc_binx,clu_roc_biny, l2->GetBinContent(clu_roc_binx,clu_roc_biny)+1);
      else if (clu.mod_on.layer==3) l3->SetBinContent(clu_roc_binx,clu_roc_biny, l3->GetBinContent(clu_roc_binx,clu_roc_biny)+1);
      else if (clu.mod_on.blade>0) {
	fpixI->SetBinContent(clu_roc_binx,clu_roc_biny, fpixI->GetBinContent(clu_roc_binx,clu_roc_biny)+1);
	fpixI->SetBinContent(clu_roc_binx,clu_roc_biny+1, fpixI->GetBinContent(clu_roc_binx,clu_roc_biny+1)+1);
      } else if (clu.mod_on.blade<0) {
	fpixO->SetBinContent(clu_roc_binx,clu_roc_biny, fpixO->GetBinContent(clu_roc_binx,clu_roc_biny)+1);
	fpixO->SetBinContent(clu_roc_binx,clu_roc_biny+1, fpixO->GetBinContent(clu_roc_binx,clu_roc_biny+1)+1);
      }
	
      if (i%500000==0) {
	c.cd(1);
	l1->Draw("COLZ");
	c.cd(2);
	l2->Draw("COLZ");
	c.cd(3);
//.........这里部分代码省略.........

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 background_rates_GMN( const char *setupfilename, const char *outfilename ){

  ifstream setupfile(setupfilename);

  TString filename;

  TChain *C = new TChain("T");
  
  while( setupfile >> filename ){
    C->Add(filename);
  }
  
  G4SBSRunData *rd;
  
  // C->Add(rootfilename);

  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 << chEl->GetTitle() << endl;
  }

  cout << "Total number of generated events = " << ngen << endl;

  gmn_tree *T = new gmn_tree(C);

  TFile *fout = new TFile( outfilename, "RECREATE" );

  TH2D *hnphe_vs_edep_HCAL = new TH2D("hnphe_vs_edep_HCAL","",500,0.0,0.5,501,-0.5,500.5);
  TH2D *hrate_vs_nphe_HCAL = new TH2D("hrate_vs_nphe_HCAL","",288,-0.5,287.5,500,0.5,500.5);
  TH2D *hrate_vs_edep_HCALscint = new TH2D("hrate_vs_edep_HCALscint","",288,-0.5,287.5,500,0.0,1.0); //GeV

  TH2D *hnphe_vs_edep_CDET = new TH2D("hnphe_vs_edep_CDET","",500,0.0,0.25,501,-0.5,500.5);
  TH2D *hrate_vs_nphe_CDET = new TH2D("hrate_vs_nphe_CDET","",2352,-0.5,2351.5,500,0.5,500.5);
  TH2D *hrate_vs_edep_CDETscint = new TH2D("hrate_vs_edep_CDETscint","",2352,-0.5,2351.5,500,0.0,1.0); //GeV

  //TH2D *hnphe_vs_edep_BBhodo = new TH2D("hnphe_vs_edep_BBhodo","",500,0.0,0.25,501,-0.5,500.5);
  //TH2D *hrate_vs_nphe_BBhodo = new TH2D("hrate_vs_nphe_BBhodo","",2352,-0.5,2351.5,500,0.5,500.5);
  TH2D *hrate_vs_edep_BBHodoScint = new TH2D("hrate_vs_edep_BBHodoScint","",90,-0.5,89.5,500,0.0,0.15); //GeV

  TH2D *hnphe_vs_edep_BBPS = new TH2D("hnphe_vs_edep_BBPS","",500,0.0,0.25,501,-0.5,500.5);
  TH2D *hrate_vs_nphe_BBPS = new TH2D("hrate_vs_nphe_BBPS","",54,-0.5,53.5,501,-0.5,500.5);
  TH2D *hrate_vs_edep_BBPSTF1 = new TH2D("hrate_vs_edep_BBPSTF1","",54,-0.5,53.5,500,0.0,0.25); //GeV

  TH2D *hnphe_vs_edep_BBSH = new TH2D("hnphe_vs_edep_BBSH","",500,0.0,0.5,501,-0.5,500.5);
  TH2D *hrate_vs_nphe_BBSH = new TH2D("hrate_vs_nphe_BBSH","",189,-0.5,188.5,501,-0.5,500.5);
  TH2D *hrate_vs_edep_BBSHTF1 = new TH2D("hrate_vs_edep_BBSHTF1","",189,-0.5,188.5,500,0.0,0.5); //GeV
  
  double pmtnum[288];
  double sumedep_HCAL[288];
  // double sum2edep_HCAL[288];
  double sumnphe_HCAL[288];
  double hitrate_HCAL[288];
  
  for( int ipmt=0; ipmt<288; ipmt++ ){
    sumedep_HCAL[ipmt] = 0.0;
    pmtnum[ipmt] = ipmt;
    sumnphe_HCAL[ipmt] = 0.0;
    hitrate_HCAL[ipmt] = 0.0;
  }

  double pmt_CDET[2352];
  double sumedep_CDET[2352];
  double sumnphe_CDET[2352];
  double hitrate_CDET[2352];
  
  for( int ipmt=0; ipmt<2352; ipmt++ ){
    pmt_CDET[ipmt] = ipmt;
    sumedep_CDET[ipmt] = 0.0;
    sumnphe_CDET[ipmt] = 0.0;
    hitrate_CDET[ipmt] = 0.0;
  }

  double pmtBBHodo[90];
  double sumedep_BBHodo[90];
  double hitrate_BBHodo[90];
  
  for( int ipmt=0; ipmt<90; ipmt++ ){
    pmtBBHodo[ipmt] = ipmt;
    sumedep_BBHodo[ipmt] = 0.0;
    hitrate_BBHodo[ipmt] = 0.0;
  }

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

// void scan(float jetPtCut=40.0, TString tag = "", bool requireLeptonic = false, float mtLow = -1.0, float mtHigh = -1.0) {
void scan(float jetPtCut=40.0, TString tag = "", bool requireLeptonic = false) {

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

    // ch->Add("/hadoop/cms/store/user/namin/mStop500to800_mLSP1/ntuple_78_2_RJ3.root");

    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_309_4_16a.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_253_1_f2L.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_343_3_l8P.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_329_3_4dB.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_336_1_yfv.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_358_1_IPX.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_346_1_ptQ.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_281_1_Qqp.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_362_1_NNy.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_361_1_Yhs.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_359_1_YrF.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop825to900_mLSP1/ntuple_360_1_sfh.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop150to475_mLSP1/*.root");
    ch->Add("/hadoop/cms/store/user/namin/mStop500to800_mLSP1/*.root");

    int nEventsTotal = ch->GetEntries();
    int nEventsSoFar = 0;
    int nGoodEvents = 0;

    TFile *currentFile = 0;
    TObjArray *listOfFiles = ch->GetListOfFiles();
    TIter fileIter(listOfFiles);

    vector<TH1F*> h1D_mindR_vec;
    TH1F* h1D_mindR = new TH1F("", "min #DeltaR between reco jet and quark;;Entries", 100, 0.0, 2.0); h1D_mindR_vec.push_back(h1D_mindR);

    vector<TH1F*> h1D_stop_mass_vec;
    TH1F* h1D_stop_mass = new TH1F("", "Stop mass (events have 2 matched jets);;Entries", 100, 100, 900); h1D_stop_mass_vec.push_back(h1D_stop_mass); 

    vector<TH1F*> h1D_stop_mass_lt2match_vec;
    TH1F* h1D_stop_mass_lt2match = new TH1F("", "Stop mass (events have <2 matched jets);;Entries", 100, 100, 900); h1D_stop_mass_lt2match_vec.push_back(h1D_stop_mass_lt2match); 

    vector<TH1F*> h1D_stop_mass_all_vec;
    TH1F* h1D_stop_mass_all = new TH1F("", "Stop mass (events passing W daughter reqs);;Entries", 100, 100, 900); h1D_stop_mass_all_vec.push_back(h1D_stop_mass_all); 

    vector<TH1F*> h1D_lonelyjet_vec;
    TH1F* h1D_lonelyjet = new TH1F("", "Just 1 jet;;Entries", 100, 100, 900); h1D_lonelyjet_vec.push_back(h1D_lonelyjet); 

    vector<TH1F*> h1D_nonlonelyjet_vec;
    TH1F* h1D_nonlonelyjet = new TH1F("", "2 jets;;Entries", 100, 100, 900); h1D_nonlonelyjet_vec.push_back(h1D_nonlonelyjet); 

    TH2F* h2D_dRqq_Wpt = new TH2F("", "#DeltaR_{qq} (no matching);W_{p_{T}};#DeltaR_{qq}", 50,0.0,900, 50,0.0,5.0);

    vector<TH1F*> h1D_njets_gen_vec;
    vector<TH1F*> h1D_njets_reco_vec;
    vector<TH1F*> h1D_dRjj_masses_vec;
    vector<TH1F*> h1D_dRqq_nomatch_vec;
    vector<TH1F*> h1D_dRcone_nomatch_vec;
    vector<TH1F*> h1D_leadjmass_masses_vec;
    vector<TH1F*> h1D_mt_masses_vec;
    for(int iStop = 0; iStop <= NSTOPMASSES; iStop++) {
        TH1F* h1D_njets_gen = new TH1F("", ";;Entries", 15, 0,15); h1D_njets_gen_vec.push_back(h1D_njets_gen);
        TH1F* h1D_njets_reco = new TH1F("", ";;Entries", 15, 0,15); h1D_njets_reco_vec.push_back(h1D_njets_reco);
        TH1F* h1D_dRjj_masses = new TH1F("", ";;Entries", 100, 0.0, 5.0); h1D_dRjj_masses_vec.push_back(h1D_dRjj_masses);
        TH1F* h1D_dRqq_nomatch = new TH1F("", ";;Entries", 100, 0.0, 5.0); h1D_dRqq_nomatch_vec.push_back(h1D_dRqq_nomatch);
        TH1F* h1D_dRcone_nomatch = new TH1F("", ";;Entries", 100, 0.0, 5.0); h1D_dRcone_nomatch_vec.push_back(h1D_dRcone_nomatch);
        TH1F* h1D_leadjmass_masses = new TH1F("", ";;Entries", 30, 0.0, 100.0); h1D_leadjmass_masses_vec.push_back(h1D_leadjmass_masses);
        TH1F* h1D_mt_masses = new TH1F("", ";;Entries", 45, 0.0, 900.0); h1D_mt_masses_vec.push_back(h1D_mt_masses);
    }

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

        // Get File Content
        TFile *file = new TFile( currentFile->GetTitle() );
        TTree *tree = (TTree*)file->Get("Events");
        cms2.Init(tree);

        TString filename(currentFile->GetTitle());

        // Loop over Events in current file
        unsigned int nEventsTree = tree->GetEntriesFast();
        for( unsigned int event = 0; event < nEventsTree; ++event) {


            // if(event > 300) break;

            addToCounter("total");

            // Get Event Content
            cms2.GetEntry(event);
            nEventsSoFar++;

            // each element of this vector is a pair of genps_ indices for the 2 quarks from a W
            vector<pair<int, int> > WtoqqIdxs; 
            // each element of this vector is a pair of goodJet indices for the 2 matched jets for the quarks
            vector<pair<int, int> > qqMatchedJets; 
            // indices of b quarks from same tops as Ws; this is in same order as WtoqqIdxs
            vector<int> bIdxs; 
            // indices of Ws
            vector<int> WIdxs;

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

void Skim(TString fname="ZnnH125", TString outputname = "")
{
    gROOT->LoadMacro("HelperFunctions.h" );  // make functions visible to TTreeFormula
    gROOT->SetBatch(1);

    TChain * chain  = new TChain("tree");
    TString dijet   = "";
    TString outdir  = "/afs/cern.ch/work/d/degrutto/public/MiniAOD/ZnnHbb_Phys14_PU20bx25/skimV11_v2/";
    TString prefix  = "skim_";
    TString suffix  = "tree*.root";
    


    TCut selection = baseline.c_str();
    // Different baseline for dimuons
    // MET filters
    selection += metfilter.c_str();

    // JSON & trigger
    if (fname.Contains("Data")) {
        TCut trigger = mettrigger.c_str();
        selection += trigger;
        //selection.Print();
      }
      /*
    } else if (process == "Data_METBTag_R" || process == "Data_METBTag_P") {
        TCut trigger = metbtagtrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleMu_R" || process == "Data_SingleMu_P") {
        TCut trigger = singlemutrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleEl_R" || process == "Data_SingleEl_P") {
        TCut trigger = singleeltrigger.c_str();
        selection += trigger;
        //selection.Print();
    }
      */


    chain->Add(fname);

    // Sum Count, CountWithPU, CountWithPUP, CountWithPUM
    TObjArray * files = chain->GetListOfFiles();
    TIter next(files);
    TChainElement * chainElem = 0;
    TFile * f2 = 0;
    TH1D * h1 = new TH1D("Count", ";Counts", 16, 0, 16);
    TH1F * htemp = 0;

    while ((chainElem = (TChainElement*) next())) {
      //#ifndef XROOTD
      //      f2 = TFile::Open("dcache:" + TString(chainElem->GetTitle()));
      //#else
      std::cout << "chainElem->GetTitle() " << chainElem->GetTitle() << std::endl; 
        f2 = TFile::Open( TString(chainElem->GetTitle()));
	//#endif
        htemp = (TH1F*) f2->Get("Count");
        h1->SetBinContent(1, h1->GetBinContent(1)+htemp->GetBinContent(1));
	/*
        htemp = (TH1F*) f2->Get("CountWithPU");
        h1->SetBinContent(2, h1->GetBinContent(2)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUP");
        h1->SetBinContent(3, h1->GetBinContent(3)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUM");
        h1->SetBinContent(4, h1->GetBinContent(4)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithMCProd");
        h1->SetBinContent(5, h1->GetBinContent(5)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUMCProd");
        h1->SetBinContent(6, h1->GetBinContent(6)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("countWithSignalQCDcorrections");
        h1->SetBinContent(7, h1->GetBinContent(7)+htemp->GetBinContent(1));
        */
        std::clog << fname << ": skimmed from " << chainElem->GetTitle() << std::endl;
    }
    
    // LHE Count
    
    TH1D * h2 = new TH1D("LHECount", ";LHE Counts", 16, 0, 16);
    TString process_lhe = fname;
    if (process_lhe.BeginsWith("WJets"))
        process_lhe = "WJets";
    else if (process_lhe.BeginsWith("ZJets"))
        process_lhe = "ZJets";
    else 
        process_lhe = "";
    const std::vector<std::string>& lhecuts = GetLHECuts(process_lhe.Data());
    for (unsigned int i=0; i < lhecuts.size(); i++) {
        TCut cut2 = lhecuts.at(i).c_str();
        h2->SetBinContent(i+1, chain->GetEntries(cut2));
    }
    
    

    // Make output directory if it doesn't exist
    if (gSystem->AccessPathName(outdir))
        gSystem->mkdir(outdir);
    TString outname = outdir + prefix + Form("%s.root", outputname.Data());
    std::cout << "outname is " << outname << std::endl;
//.........这里部分代码省略.........

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

      } // end loop over mass pts (1 pt only if not signal scan)
    } // end loop over genClassifications for histogram arrays
  } // end loop over analyzers/regions for histogram arrays


	// Set up cutflow histograms
	TH1D* h_cutflow[nAnalyzers];
	for( int iAna=0; iAna<nAnalyzers; iAna++ ) {

		analyzer* thisAnalyzer = analyzers.at(iAna);
		std::string histName  = "h_cutflow_";
		std::string histTitle = "Cutflow histogram ";
		histName  += thisAnalyzer->GetLabel();
		histTitle += thisAnalyzer->GetLabel();

		h_cutflow[iAna] = selectionInfo::get_cutflowHistoTemplate( thisAnalyzer->GetSelections(), histName, histTitle );
	}


  //////////////////////
  //                  //
  // Loop Over Events //
  //                  //
  //////////////////////
    
  // Event Counters
  cout << "    Loading files to loop over" << endl << endl;
  unsigned int nEventsTotal = 0;
  unsigned int nEventsChain = chain->GetEntries();
  if( nEvents >= 0 ) nEventsChain = nEvents;
  
  // Grab list of files
  TObjArray *listOfFiles = chain->GetListOfFiles();
  TIter fileIter(listOfFiles);
  TFile *currentFile = 0;
  while ( (currentFile = (TFile*)fileIter.Next()) ) {
    
    
    //////////////////////
    //                  //
    // Get File Content //
    //                  //
    //////////////////////
    
    // Open File and Get Tree
    TFile *file = new TFile( currentFile->GetTitle(), "read" );
    TTree *tree = (TTree*)file->Get("t");
    if(readFast) TTreeCache::SetLearnEntries(10);
    if(readFast) tree->SetCacheSize(128*1024*1024);
    babyAnalyzer.Init(tree);
    
    // Get weight histogram from baby
    wgtInfo.getWeightHistogramFromBaby( file );
    
    // 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 //
      //                   //
      ///////////////////////

//
// 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
//.........这里部分代码省略.........

void flux(const char *rootfilename, const char *outputfilename){

  gStyle->SetMarkerStyle(20);
  
  TH1::SetDefaultSumw2();
  
  TChain *C = new TChain("T");
  C->Add(rootfilename);

  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 << chEl->GetTitle() << endl;
  }

  cout << "number of generated events = " << ngen << endl;

  C16_tree_with_flux *T = new C16_tree_with_flux(C);
  
  TFile *fout = new TFile(outputfilename,"RECREATE");

  long nevent=0;

  double cutgamma = 0.05*1.e-3; //energy cutoff in GeV for photons
  double cutelectron = 5.e-3; //range cutoff in GeV for electrons:
  
  TH1D *hphoton_rate_vs_theta = new TH1D("hphoton_rate_vs_theta","",80,0.0,40.0);
  TH1D *helectron_rate_vs_theta = new TH1D("helectron_rate_vs_theta","",80,0.0,40.0);
  TH1D *hphoton_rate_vs_theta_Edep = new TH1D("hphoton_rate_vs_theta_Edep","",80,0.0,40.0);
  TH1D *helectron_rate_vs_theta_Edep = new TH1D("helectron_rate_vs_theta_Edep","",80,0.0,40.0);
  TH1D *hother_rate_vs_theta = new TH1D("hother_rate_vs_theta","",80,0.0,40.0);
  TH1D *hother_rate_vs_theta_Edep = new TH1D("hother_rate_vs_theta_Edep","",80,0.0,40.0);

  TH1D *htotal_rate_vs_theta = new TH1D("htotal_rate_vs_theta","",80,0.0,40.0);
  TH1D *htotal_rate_vs_theta_Edep = new TH1D("htotal_rate_vs_theta_Edep","",80,0.0,40.0);

  double PI = TMath::Pi();

  double binwidth_theta = 40.0/80.0*PI/180.0;

  double Ibeam = 20e-6;
  double e = 1.602e-19;
  
  while( T->GetEntry( nevent++ ) ){
    TFile *f = ( (TChain*) (T->fChain) )->GetFile();

    TString fname = f->GetName();
    
    if( bad_file_list.find( fname ) == bad_file_list.end() ){ 
    
      if( nevent%1000 == 0 ) {
	cout << nevent << endl;
	cout << "Current file = " << f->GetName() << endl;
      }
      
      for( int part=0; part<T->FLUX_npart_CAL; part++ ){
	double E = (*(T->FLUX_E))[part];
	double px = (*(T->FLUX_px))[part];
	double py = (*(T->FLUX_py))[part];
	double pz = (*(T->FLUX_pz))[part];
	double p = (*(T->FLUX_p))[part];
	double theta = acos( pz/p );
	double phi = atan2( py, px );

	int bintheta = hphoton_rate_vs_theta->FindBin( theta * 180.0/PI );

	double dcostheta = cos( PI/180.0*hphoton_rate_vs_theta->GetBinLowEdge(bintheta) ) -
	  cos( PI/180.0*hphoton_rate_vs_theta->GetBinLowEdge(bintheta+1) );
	double dOmega = 2.0*PI *dcostheta;

	double weight = Ibeam/e/double(ngen)/dOmega;
	
	int pid = (*(T->FLUX_pid))[part];

	htotal_rate_vs_theta->Fill( theta*180.0/PI, weight );
	htotal_rate_vs_theta_Edep->Fill( theta*180.0/PI, weight*E );
	
	if( pid == 22 && E >= cutgamma ){//photon
	  hphoton_rate_vs_theta->Fill( theta*180.0/PI, weight );
	  hphoton_rate_vs_theta_Edep->Fill( theta*180.0/PI, weight*E );
	} else if( fabs(pid) == 11 && E >= cutelectron ){ //e+/e-
	  helectron_rate_vs_theta->Fill( theta*180.0/PI, weight );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      ch -> Add( Form("%s/WToLNu_FOv3_testFinal_baby.root",babyPath) );
    }
    else if( strcmp( samples.at(i) , "Higgs130" ) == 0 ){
      ch -> Add( Form("%s/HToWWTo2L2NuM130_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWToLNuTauNuM130_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWTo2Tau2NuM130_PU_testFinal_baby.root",babyPath) );
    }
    else if( strcmp( samples.at(i) , "Higgs160" ) == 0 ){
      ch -> Add( Form("%s/HToWWTo2L2NuM160_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWToLNuTauNuM160_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWTo2Tau2NuM160_PU_testFinal_baby.root",babyPath) );
    }
    else if( strcmp( samples.at(i) , "Higgs200" ) == 0 ){
      ch -> Add( Form("%s/HToWWTo2L2NuM200_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWToLNuTauNuM200_PU_testFinal_baby.root",babyPath) );
      ch -> Add( Form("%s/HToWWTo2Tau2NuM200_PU_testFinal_baby.root",babyPath) );
    }
    else{
      ch -> Add( Form("%s/%s_PU_testFinal_baby.root",babyPath,samples.at(i)) );
    }

    // --- Event loop

    // Prepare the event tree
    // - here the variable names have to corresponds to your tree
    // - you can use the same variables as above which is slightly faster,
    //   but of course you can use different ones and copy the values inside the event loop
    //
  
    TTree *theTree     = (TTree*) ch;

    std::cout << "--- Using input files: -------------------" <<  std::endl;

    TObjArray *listOfFiles = ch->GetListOfFiles();
    TIter fileIter(listOfFiles);
    TChainElement* currentFile = 0;
    
    while((currentFile = (TChainElement*)fileIter.Next())) {
      std::cout << currentFile->GetTitle() << std::endl;
    }

    Float_t lephard_pt_;
    Float_t lepsoft_pt_;
    Float_t lepsoft_fr_;
    Float_t dil_dphi_;
    Float_t dil_mass_;
    Float_t event_type_;
    Float_t met_projpt_;
    Int_t   jets_num_;
    Int_t   extralep_num_;
    Int_t   lowptbtags_num_;
    Int_t   softmu_num_;
    Float_t event_scale1fb_;
    Float_t met_pt_;
    Int_t   lepsoft_passTighterId_;
    Float_t mt_lephardmet_;
    Float_t mt_lepsoftmet_;
    Float_t mthiggs_;
    Float_t dphi_lephardmet_;
    Float_t dphi_lepsoftmet_;
    Float_t lepsoft_fbrem_;
    Float_t lepsoft_eOverPIn_;
    Float_t lepsoft_q_;
    Float_t lepsoft_dPhiIn_;

    theTree->SetBranchAddress( "lephard_pt_"             ,   &lephard_pt_              );