C++ TTree::GetEntriesFast方法代码示例

void MyFunctions::BuildGradeNormMap() {
	// This function build the gradeNormMap if it isn't already filled
	if (gradeNormMap.size() > 0) return;  // Already filled
	
	// This is build from the Students root file.  Maybe this should be from a course file instead?
	TFile* f = new TFile("Students.root");
	TTree* studentTree = (TTree*)f->Get("Students");
	Student* student = 0;
	studentTree->SetBranchAddress("student", &student);
	
	Long64_t nentries = studentTree->GetEntriesFast();
	for (Long64_t jentry = 0; jentry < nentries; jentry++) {
		studentTree->GetEntry(jentry);
		//Loop over all terms for this student and accumulate grade distribution info
		for (auto const& grade : student->Grades()) {
			// Only look at Fall and Spring Terms
			int term = grade.term;
			if (!regularSemester(term)) continue;
			if (!ValidGrade(grade.grade)) continue;
			auto thisKey = std::make_pair(grade.course, 0);
			auto thisKeyAll = std::make_pair("AllCourses", 0);
			MyFunctions::gradeNormMap[thisKey].AddGrade(grade.grade);  // Will automatically create entry in map if it doesn't exist
			MyFunctions::gradeNormMap[thisKeyAll].AddGrade(grade.grade);
			auto termKey = std::make_pair(grade.course, term);
			auto termKeyAll = std::make_pair("AllCourses", term);
			MyFunctions::gradeNormMap[termKey].AddGrade(grade.grade);
			MyFunctions::gradeNormMap[termKeyAll].AddGrade(grade.grade);
		}			
	}
	
	f->Close();
}

void CalibSort::Loop() {
  //   In a ROOT session, you can do:
  //      Root > .L CalibSort.C
  //      Root > CalibSort t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries
  //

  //   This is the loop skeleton where:
  //      jentry is the global entry number in the chain
  //      ientry is the entry number in the current Tree
  //   Note that the argument to GetEntry must be:
  //      jentry for TChain::GetEntry
  //      ientry for TTree::GetEntry and TBranch::GetEntry
  //
  //       To read only selected branches, Insert statements like:
  // METHOD1:
  //    fChain->SetBranchStatus("*",0);  // disable all branches
  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
  // METHOD2: replace line
  //    fChain->GetEntry(jentry);       //read all branches
  //by  b_branchname->GetEntry(ientry); //read only this branch
  if (fChain == 0) return;

  std::ofstream fileout;
  if ((flag_%10)==1) {
    fileout.open("events.txt", std::ofstream::out);
    std::cout << "Opens events.txt in output mode\n";
  } else {
    fileout.open("events.txt", std::ofstream::app);
    std::cout << "Opens events.txt in append mode\n";
  }
  fileout << "Input file: " << fname_ << " Directory: " << dirnm_ 
	  << " Prefix: " << prefix_ << "\n";
  Long64_t nbytes(0), nb(0), good(0);
  Long64_t nentries = fChain->GetEntriesFast();
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    double cut = (t_p > 20) ? 10.0 : 0.0;
    if ((flag_/10)%10 > 0) 
      std::cout << "Entry " << jentry << " p " << t_p << " Cuts " << t_qltyFlag
		<< "|" << t_selectTk << "|" << (t_hmaxNearP < cut) << "|" 
		<< (t_eMipDR < mipCut_) << std::endl;
    if (t_qltyFlag && t_selectTk && (t_hmaxNearP<cut) && (t_eMipDR<mipCut_)) {
      good++;
      fileout << good << " " << jentry << " " << t_Run  << " " << t_Event 
	      << " " << t_ieta << " " << t_p << std::endl;
    }
  }
  fileout.close();
  std::cout << "Writes " << good << " events in the file events.txt from "
	    << nentries << " entries" << std::endl;
}

void spectraMaker(std::string file)
{
  Double_t fidCut = 50.;
  
  std::string filePath = "analyzed_files/"+file;

  TFile *fout = new TFile("outputHistsTest.root","RECREATE");
  TH1D *hEreconALL = new TH1D("hEreconALL","Reconstructed Energy Spectra for All events",120,0., 1200.);
  TH1D *hErecon0 = new TH1D("hErecon0","Reconstructed Energy Spectra for All events",120,0., 1200.);
  TH1D *hErecon1 = new TH1D("hErecon1","Reconstructed Energy Spectra for All events",120,0., 1200.);
  TH1D *hErecon23 = new TH1D("hErecon23","Reconstructed Energy Spectra for All events",120,0., 1200.);

  hEreconALL->GetXaxis()->SetTitle("E_{recon} (keV)");
  hErecon0->GetXaxis()->SetTitle("E_{recon} (keV)");
  hErecon1->GetXaxis()->SetTitle("E_{recon} (keV)");
  hErecon23->GetXaxis()->SetTitle("E_{recon} (keV)");
  
  TFile *fin = new TFile(filePath.c_str(),"READ");
  TTree *tin = (TTree*)(fin->Get("SimAnalyzed"));

  //variables to be read in
  Int_t PID, side, type;
  Double_t Erecon;
  Double_t mwpcPosW[3], mwpcPosE[3];
  
  tin->SetBranchAddress("PID",&PID);
  tin->SetBranchAddress("side",&side);
  tin->SetBranchAddress("type",&type);
  tin->SetBranchAddress("Erecon",&Erecon);
  tin->GetBranch("MWPCPosAdjusted")->GetLeaf("MWPCPosAdjE")->SetAddress(mwpcPosE);
  tin->GetBranch("MWPCPosAdjusted")->GetLeaf("MWPCPosAdjW")->SetAddress(mwpcPosW);

  UInt_t nevents = tin->GetEntriesFast();

  for (UInt_t evt=0; evt<nevents; evt++) {
    tin->GetEvent(evt);

    if (PID!=1) continue;
    if ((mwpcPosW[0]*mwpcPosW[0]+mwpcPosW[1]*mwpcPosW[1]>fidCut*fidCut) || (mwpcPosE[0]*mwpcPosE[0]+mwpcPosE[1]*mwpcPosE[1]>fidCut*fidCut)) continue;

    //Fill appropriate histograms if cuts are passed.
    if (type<4) hEreconALL->Fill(Erecon);
    if (type==0) hErecon0->Fill(Erecon);
    else if (type==1) hErecon1->Fill(Erecon);
    else if (type==2 || type==3) hErecon23->Fill(Erecon);  
  }

  fin->Close();

  fout->Write();
  fout->Close();
}

void PileUpMaker() {
  // Access data file
  TFile *file = TFile::Open("ZmumuGammaNtuple_Full2012_MuCorr.root");
  TTree *Datatree = (TTree*)file->Get("ZmumuGammaEvent");

  // Set up histogram for data pile up
  TH1F *pileupraw = new TH1F("pileup","",40,0,80);
 
  // Fill pile up distribution
  Int_t nentries = Datatree->GetEntriesFast();
  Int_t nbytes = 0; 
  for (Int_t i=0; i<nentries;i++) {
    nbytes += Datatree->GetEvent(i);
    UInt_t NVertices;
    Datatree->SetBranchAddress("NVertices",&NVertices);
    pileupraw->Fill(float(NVertices)/0.7);
  }

  // Normalize the distribution
  TH1F* PU = (TH1F*)pileupraw->Clone("pileup");
  PU->Scale(1./PU->Integral());

  // Plot both the normalized and unnormalized pile up
  TCanvas *cv = 0;
  cv= new TCanvas("cv","cv",800, 600);
  pileupraw->Draw("");
  cv->SaveAs("DataPU.pdf");

  cv= new TCanvas("cv","cv",800,600);
  PU->Draw("");
  cv->SaveAs("DataPU_normalized.pdf");

  // Save the normalized distribution to a root file for weighting MC
  TFile *file = TFile::Open("2012_PhosphorData_PileUp.root","UPDATE");
  file->cd();
  file->WriteTObject(PU,PU->GetName(),"WriteDelete");
  file->Close();
  delete file;
}

void beamSpot() {

	gSystem->Load("libCintex.so");
	//gSystem->Load("libFWCoreFWLite.so");
	gSystem->Load("libRecoParticleFlowPFAnalyses");
	gSystem->Load("libDataFormatsParticleFlowReco");
	Cintex::Enable();


	TFile* f = TFile::Open("PFlowTB_Tree_9GeV_2k.root");
	TTree* tree = (TTree*) gDirectory->FindObjectAny("Extraction");

	//gPad->Divide(2);
	TH2F beamSpotHcal("beamSpotHcal", "beamSpotHcal;eta;phi", 100, 0, 3, 100, -1, 1);
	TH2F beamSpotEcal("beamSpotEcal", "beamSpotEcal;eta;phi", 100, 0, 3, 100, -1, 1);

	using namespace pftools;
	using namespace std;
	Calibratable* calib = new Calibratable();
	tree->SetBranchAddress("Calibratable", &calib);

	for (unsigned j(0); j < tree->GetEntriesFast(); ++j) {
			tree->GetEntry(j);
			for(unsigned k(0); k < calib->tb_numEcal_; ++k) {
				beamSpotEcal.Fill(calib->tb_ecal_[k].eta_,calib->tb_ecal_[k].phi_, calib->tb_ecal_[k].energy_);
			}
			for(unsigned k(0); k < calib->tb_numHcal_; ++k) {
				beamSpotHcal.Fill(calib->tb_hcal_[k].eta_,calib->tb_hcal_[k].phi_, calib->tb_hcal_[k].energy_);
			}
	}
	//gPad->cd(1);
	beamSpotEcal.Draw("colz");
	//gPad->cd(2);
	//beamSpotHcal.Draw("colz");

}

void number_of_events_table_xi1690(TString inFile=""){

	//*** In file

	TFile * input = new TFile(inFile, "READ");

	TTree * ntpMC = (TTree*) input->Get("ntpMC");
	TTree * ntpPiMinus = (TTree*) input->Get("ntpPiMinus");
	TTree * ntpPiPlus = (TTree*) input->Get("ntpPiPlus");
	TTree * ntpKaonMinus = (TTree*) input->Get("ntpKaonMinus");
	TTree * ntpProton = (TTree*) input->Get("ntpProton");
	TTree * ntpAntiProton = (TTree*) input->Get("ntpAntiProton");
	TTree * ntpLambda0 = (TTree*) input->Get("ntpLambda0");
	TTree * ntpAntiLambda0 = (TTree*) input->Get("ntpAntiLambda0");
	TTree * ntpXiPlus = (TTree*) input->Get("ntpXiPlus");
	TTree * ntpXiMinus1690 = (TTree*) input->Get("ntpXiMinus1690");
	TTree * ntpXiSys = (TTree*) input->Get("ntpXiSys");


	double nevents_mc = ntpMC->GetEntriesFast();
	TString cuts = " McTruthMatch && VtxFit_HowGood==1 && MassFit_prob>0.01";
	TString VtxCut = " McTruthMatch && VtxFit_HowGood==1 & HitTag==1";
	TString cut4c = "McTruthMatch && 4CFit_prob>0.01";


	cout << "particle|   #evts (uncut)|    #evts (ratio in %)|   MC ratio in %|   dp/p in %" << endl;


	//**** PiMinus
	TH1D * h_piminus_tht_uncut = new TH1D("h_piminus_tht_uncut", "h_piminus_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus_tht_uncut", "piminus_tht", "McTruthMatch && Mother==3122");
	double piminus_uncut =  h_piminus_tht_uncut->GetEntries();

	TH1D * h_piminus_tht = new TH1D("h_piminus_tht", "h_piminus_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus_tht", "piminus_tht", "McTruthMatch && piminus_HitTag && Mother==3122");
	int piminus =  h_piminus_tht->GetEntries();

	TH1D * h_piminus_dp = new TH1D("h_piminus_dp", "h_piminus_dp", 250, -0.1,0.1);
	ntpPiMinus->Project("h_piminus_dp", "(piminus_p-piminus_MC_p)/piminus_MC_p", "McTruthMatch && piminus_HitTag && Mother==3122");

//	Double_t param[6] = jenny::GetFitParameterDoubleFit(h_piminus_dp, false, 0.02,0.1, true);
//	jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piminus_dp, "","", false, false, false, 0.02,0.1, true);

	double ratio_piminus_cut = piminus/piminus_uncut;
	double ratio_piminus_mc = piminus/nevents_mc;

	cout << "PiMinus|   " <<  piminus_uncut << "|   " <<  piminus << "(" << ratio_piminus_cut*100 << ")|   " << ratio_piminus_mc*100 << endl; //param[2]*100 << endl;


	//**** PiPlus (AntiLambda0)
	TH1D * h_piplus_tht_uncut = new TH1D("h_piplus_tht_uncut", "h_piplus_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus_tht_uncut", "piplus_tht", "McTruthMatch && Mother==-3122");
	double piplus_uncut =  h_piplus_tht_uncut->GetEntries();

	TH1D * h_piplus_tht = new TH1D("h_piplus_tht", "h_piplus_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus_tht", "piplus_tht", "McTruthMatch && piplus_HitTag && Mother==-3122");
	int piplus =  h_piplus_tht->GetEntries();

	double ratio_piplus_cut = piplus/piplus_uncut;
	double ratio_piplus_mc = piplus/nevents_mc;

	TH1D * h_piplus_dp = new TH1D("h_piplus_dp", "h_piplus_dp", 250, -0.1,0.1);
	ntpPiPlus->Project("h_piplus_dp", "(piplus_p-piplus_MC_p)/piplus_MC_p", "McTruthMatch && piplus_HitTag && Mother==-3122");

//	Double_t parampip[6] = jenny::GetFitParameterDoubleFit(h_piplus_dp, false, 0.02,0.1, true);
//	jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piplus_dp, "","", false, false, false, 0.02,0.1, true);

	cout << "PiPlus(AL0)|   " <<  piplus_uncut << "|   " <<  piplus << "(" << ratio_piplus_cut*100 << ")|   " << ratio_piplus_mc*100 << endl; //pip[2]*100 << endl;



	//**** piplus (Xi+)
	TH1D * h_piplus2_tht_uncut = new TH1D("h_piplus2_tht_uncut", "h_piplus2_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus2_tht_uncut", "piplus_tht", "McTruthMatch && Mother==-3312");
	double piplus2_uncut =  h_piplus2_tht_uncut->GetEntries();

	TH1D * h_piplus2_tht = new TH1D("h_piplus2_tht", "h_piplus2_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus2_tht", "piplus_tht", "McTruthMatch && piplus_HitTag && Mother==-3312");
	int piplus2 =  h_piplus2_tht->GetEntries();

	TH1D * h_piplus2_dp = new TH1D("h_piplus2_dp", "h_piplus_dp", 250, -0.1,0.1);
	ntpPiPlus->Project("h_piplus2_dp", "(piplus_p-piplus_MC_p)/piplus_MC_p", "McTruthMatch && piplus_HitTag && Mother==-3312");

//	Double_t parampip2[6] = jenny::GetFitParameterDoubleFit(h_piplus2_dp, false, 0.02,0.1, true);
//	jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piplus2_dp, "","", false, false, false, 0.02,0.1, true);

	double ratio_piplus2_cut = piplus2/piplus2_uncut;
	double ratio_piplus2_mc = piplus2/nevents_mc;

	cout << "PiPlus(Xi+)|" <<  piplus2_uncut << "|" <<  piplus2 << "(" << ratio_piplus2_cut*100 << ")|" << ratio_piplus2_mc*100 << endl; //pip2[2]*100 << endl;


	//**** kaonMinus
	TH1D * h_kaonMinus_tht_uncut = new TH1D("h_kaonMinus_tht_uncut", "h_kaonMinus_tht", 100, 0,10);
	ntpKaonMinus->Project("h_kaonMinus_tht_uncut", "kaonminus_tht", "McTruthMatch && Mother==13314");
	double kaonMinus_uncut =  h_kaonMinus_tht_uncut->GetEntries();

	TH1D * h_kaonMinus_tht = new TH1D("h_kaonMinus_tht", "h_kaonMinus_tht", 100, 0,10);
	ntpKaonMinus->Project("h_kaonMinus_tht", "kaonminus_tht", "McTruthMatch && kaonminus_HitTag && Mother==13314");
	int kaonMinus =  h_kaonMinus_tht->GetEntries();
//.........这里部分代码省略.........

void readNtuplesNoFilter_forAN_isMC(TString inputfilename="results.root", int flavor=Sig::Prompt,  std::string effmeasured="IterL3_NOHP_NOL1"){

  bool doingL1 = effmeasured.find("L1OverOffline") != std::string::npos; 
  
  TFile* outfile = TFile::Open(Form("%s_efficiency.root", effmeasured.c_str()),"RECREATE");
  std::cout << "output file: " << outfile -> GetName() << std::endl;

  //Create histograms  
  TH1F* dimuon_mass             = new TH1F("h_dimuon_mass"          ,"dimuon_mass"      , 1500,  0,  150 );
  TH1F* tagiso                  = new TH1F("h_tagiso"               ,"tagiso"           ,  100,  0,  1   );
  TH1F* tagMuonPt               = new TH1F("h_tagMuonPt"            ,"tagMuonPt"        ,  150,  0,  150 );
  TH1F* nvtx_event              = new TH1F("h_nvtx_event"           ,"nvtx_event"       ,   60,  0,   60 );
  
  // L1 over GEN:
  TEfficiency* L1muonPt         = new TEfficiency("L1muonPt"        ,"L1muonPt"         ,   16,  pt_bins ); 
  TEfficiency* L1muonEta        = new TEfficiency("L1muonEta"       ,"L1muonEta"        ,   15, eta_bins );
  TEfficiency* L1muonPhi        = new TEfficiency("L1muonPhi"       ,"L1muonPhi"        ,   20, -3.2, 3.2);
  TEfficiency* L1muonEff        = new TEfficiency("L1muonEff"       ,"L1muonEff"        ,    1,   0., 1.0);
  								    
  // L2 over L1: 						    
  TEfficiency* L2muonPt         = new TEfficiency("L2muonPt"        ,"L2muonPt"         ,   16,  pt_bins ); 
  TEfficiency* L2muonEta        = new TEfficiency("L2muonEta"       ,"L2muonEta"        ,   15, eta_bins );
  TEfficiency* L2muonPhi        = new TEfficiency("L2muonPhi"       ,"L2muonPhi"        ,   20, -3.2, 3.2);
  TEfficiency* L2muonEff        = new TEfficiency("L2muonEff"       ,"L2muonEff"        ,    1,   0., 1.0);

  // L3 over L1:
  TEfficiency* muonPt           = new TEfficiency("muonPt"          ,"muonPt"           ,   16,  pt_bins ); 
  TEfficiency* muonPtTurnOn     = new TEfficiency("muonPtTurnOn"    ,"muonPtTurnOn"     ,   16,  pt_bins ); 
  TEfficiency* muonEta          = new TEfficiency("muonEta"         ,"muonEta"          ,   15, eta_bins );
  TEfficiency* muonPhi          = new TEfficiency("muonPhi"         ,"muonPhi"          ,   20, -3.2, 3.2);
  TEfficiency* muonEff          = new TEfficiency("muonEff"         ,"muonEff"          ,    1,   0., 1.0);
  TEfficiency* muonDeltaR       = new TEfficiency("muonDeltaR"      ,"muonDeltaR"       ,   30,   0., 3.0);
  TEfficiency* muonDeltaPhi     = new TEfficiency("muonDeltaPhi"    ,"muonDeltaPhi"     ,   30,   0., 3.2); 

  TEfficiency* failingMuonPt    = new TEfficiency("failingMuonPt"   ,"failingMuonPt"    ,   16,  pt_bins ); 
  TEfficiency* failingMuonEta   = new TEfficiency("failingMuonEta"  ,"failingMuonEta"   ,   15, eta_bins );
  TEfficiency* failingMuonPhi   = new TEfficiency("failingMuonPhi"  ,"failingMuonPhi"   ,   20, -3.2, 3.2);
  TEfficiency* failingMuonEff   = new TEfficiency("failingMuonEff"  ,"failingMuonEff"   ,   1 ,   0., 1.0);

  TEfficiency* fakeMuonPt    = new TEfficiency("fakeMuonPt"   ,"fakeMuonPt"    ,   16,  pt_bins ); 
  TEfficiency* fakeMuonEta   = new TEfficiency("fakeMuonEta"  ,"fakeMuonEta"   ,   15, eta_bins );
  TEfficiency* fakeMuonPhi   = new TEfficiency("fakeMuonPhi"  ,"fakeMuonPhi"   ,   20, -3.2, 3.2);
  TEfficiency* fakeMuonEff   = new TEfficiency("fakeMuonEff"  ,"fakeMuonEff"   ,   1 ,   0., 1.0);

  TH1F* PassingProbePt          = new TH1F("h_PassingProbePt"       ,"PassingMuonPt"    ,  16,  pt_bins );
  TH1F* PassingProbeEta         = new TH1F("h_PassingProbeEta"      ,"PassingMuonEta"   ,  15, eta_bins );
  TH1F* PassingProbePhi         = new TH1F("h_PassingProbePhi"      ,"PassingMuonPhi"   ,  20, -3.2, 3.2);
  TH1F* PassingProbeMll         = new TH1F("h_PassingProbeMll"      ,"PassingMuonMll"   ,  40, 60., 120.); 
  TH1F* PassingProbeDeltaR      = new TH1F("h_PassingProbeDeltaR"   ,"PassingMuonDeltaR",  30,  0.,  3.0); 

  TH1F* FailingProbePt          = new TH1F("h_FailingProbePt"       ,"FailingMuonPt"    ,  16,  pt_bins );
  TH1F* FailingProbeEta         = new TH1F("h_FailingProbeEta"      ,"FailingMuonEta"   ,  15, eta_bins );
  TH1F* FailingProbePhi         = new TH1F("h_FailingProbePhi"      ,"FailingMuonPhi"   ,  20, -3.2, 3.2);
  TH1F* FailingProbeMll         = new TH1F("h_FailingProbeMll"      ,"FailingMuonMll"   ,  40,  60., 120.);
  TH1F* FailingProbeDeltaR      = new TH1F("h_FailingProbeDeltaR"   ,"FailingMuonDeltaR",  30,  0.,  3.0); 

  // Di-muon efficiencies 
  TEfficiency* diMuonPt         = new TEfficiency("diMuonPt"      ,"diMuonPt"       ,   16,  pt_bins  , 16,  pt_bins ); 
  TEfficiency* diMuonEta        = new TEfficiency("diMuonEta"     ,"diMuonEta"      ,   15, eta_bins  , 15, eta_bins );
  TEfficiency* diMuonPhi        = new TEfficiency("diMuonPhi"     ,"diMuonPhi"      ,   20, -3.2, 3.2 , 20, -3.2, 3.2);
  TEfficiency* diMuonEff        = new TEfficiency("diMuonEff"     ,"diMuonEff"      ,    1,   0., 1.0);
  TEfficiency* diMuonDeltaR     = new TEfficiency("diMuonDeltaR"  ,"diMuonDeltaR"   ,   30,   0., 3.0);
  TEfficiency* diMuonLeadPt     = new TEfficiency("diMuonLeadPt"  ,"diMuonLeadPt"   ,   16,  pt_bins ); 
  TEfficiency* diMuonLeadEta    = new TEfficiency("diMuonLeadEta" ,"diMuonLeadEta"  ,   15,  eta_bins );
  TEfficiency* diMuonLeadPhi    = new TEfficiency("diMuonLeadPhi" ,"diMuonLeadPhi"  ,   20, -3.2, 3.2);
  TEfficiency* diMuonTrailPt    = new TEfficiency("diMuonTrailPt" ,"diMuonTrailPt"  ,   16,  pt_bins ); 
  TEfficiency* diMuonTrailEta   = new TEfficiency("diMuonTrailEta","diMuonTrailEta" ,   15, eta_bins );
  TEfficiency* diMuonTrailPhi   = new TEfficiency("diMuonTrailPhi","diMuonTrailPhi" ,   20, -3.2, 3.2);
  
  TEfficiency* nvtx             = new TEfficiency("nvtx"             ,"nvtx"             ,   60,    0,  60);
  TEfficiency* nvtx_barrel      = new TEfficiency("nvtx_barrel"      ,"nvtx_barrel"      ,   60,    0,  60);
  TEfficiency* nvtx_endcap      = new TEfficiency("nvtx_endcap"      ,"nvtx_endcap"      ,   60,    0,  60);   
 
  TFile* inputfile = TFile::Open(inputfilename, "READ");
  std::cout << "input file: " << inputfile -> GetName() << std::endl;

  TTree *tree = (TTree*) inputfile -> Get("muonNtuples/muonTree"); 
  
  if (!tree) {
    std::cout << " *** tree not found *** " << std::endl;
    return;
  }
    
  MuonEvent* ev      = new MuonEvent(); 
  TBranch*  evBranch = tree->GetBranch("event"); 
  evBranch -> SetAddress(&ev);

  int nentries = tree->GetEntriesFast();
  std::cout << "Number of entries = " << nentries << std::endl;

  bool flagfile = false;
  offlinePtCut = getLeadingPtCut(flavor);
  float ptcut1 = getLeadingPtCut(flavor);
  float ptcut2 = getTrailingPtCut(flavor);
      
  for (Int_t eventNo=0; eventNo < nentries; eventNo++)     {
    Int_t IgetEvent   = tree   -> GetEvent(eventNo);
    if (!debug) printProgBar((int)(eventNo*100./nentries));
    if (debug && eventNo==100) break;

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

int main(int argc, char* argv[]){
  using namespace muon_pog;


  if (argc < 3) 
    {
      std::cout << "Usage : "
		<< argv[0] << " PATH_TO_INPUT_FILE PAT_TO_CONFIG_FILE(s)\n";
      exit(100);
    }

  // Input root file
  TString fileName = argv[1];

  std::cout << "[" << argv[0] << "] Processing file " << fileName.Data() << std::endl;
  
  std::vector<Plotter> plotters;
  for (int iConfig = 2; iConfig < argc; ++iConfig)
    {
        std::cout << "[" << argv[0] << "] Using config file " << argv[iConfig] << std::endl;
	plotters.push_back(std::string(argv[iConfig]));
    }
  
  // Set it to kTRUE if you do not run interactively
  gROOT->SetBatch(kTRUE); 

  // Initialize Root application
  TRint* app = new TRint("CMS Root Application", &argc, argv);

  //setTDRStyle(); what to do here?
   
  // Initialize pointers to summary and full event structure
 
  muon_pog::Event* ev = new muon_pog::Event();
  TTree* tree;
  TBranch* evBranch;

  // Open file, get tree, set branches

  TFile* inputFile = TFile::Open(fileName,"READONLY");
  tree = (TTree*)inputFile->Get("MUONPOGTREE");
  if (!tree) inputFile->GetObject("MuonPogTree/MUONPOGTREE",tree);

  evBranch = tree->GetBranch("event");
  evBranch->SetAddress(&ev);

  system("mkdir -p results");
  
  TFile* outputFile = TFile::Open("results/results.root","RECREATE"); // CB find a better name for output file  

  for (auto & plotter : plotters)
    plotter.book(outputFile);
      
  // Watch number of entries
  int nEntries = tree->GetEntriesFast();
  std::cout << "[" << argv[0] << "] Number of entries = " << nEntries << std::endl;

  int nFilteredEvents = 0;
  
  for (Long64_t iEvent=0; iEvent<nEntries; ++iEvent) 
    {
      if (tree->LoadTree(iEvent)<0) break;

      evBranch->GetEntry(iEvent);

      for (auto & plotter : plotters)
	plotter.fill(ev->muons, ev->hlt);

    }

  outputFile->Write();
  
  if (!gROOT->IsBatch()) app->Run();

  return 0;
}

//.........这里部分代码省略.........
	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();
		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
			CMS3::progress( nEventsTotal, nEventsChain );

			////////////////////////////////////////////////////////////////////////////////////////////////////////
			// Analysis Code
			// ---------------------------------------------------------------------------------------------------//


			///////////////////////////////////////////////////////////////
			// Special filters to more finely categorize background events
			if(      sampleName == "tt2l"  && gen_nfromtleps_() != 2 ) continue;  //Require 2 leps from top in "tt2l" events
			else if( sampleName == "tt1l"  && gen_nfromtleps_() != 1 ) continue;  //Require 1 lep from top in "tt1l" events

			// Stitch W+NJets samples together by removing the MET<200 events from the non-nupT samples
			if( sampleName.Contains("wjets") && filename.Contains("JetsToLNu_madgraph") && nupt()>=200. ) continue;

			//FastSim anomalous event filter
			if( isFastsim && !context::filt_fastsimjets() ) continue;

			if( !mySample->PassSelections() ) continue;


			/////////////////////////////////
			// Set event weight

void ptRatioTnP()
{
    TCanvas* cnv[3][2];

    gStyle->SetLegendBorderSize(0);
    gStyle->SetLegendFont(42);
    gStyle->SetOptStat("");
    gStyle->SetTitleSize(0.7, "X");
    gStyle->SetTitleOffset(0.7, "X");

    map<string, pair<string, int> > samples;
    map<string, TH1F*> histos[3];
    THStack* bkgStack[3];
    bkgStack[0] = new THStack("bkg_ptr", "p_T Ratio: 1 fb^{-1};p_T-ratio");
    bkgStack[1] = new THStack("bkg_m", "Invariant Mass: 1 fb^{-1};M_{e#gamma}");
    bkgStack[2] = new THStack("bkg_Xpt", "Tag and Probe pair boost: 1 fb^{-1};p_{T}");
    TH1F* voids[3];
    voids[0] = new TH1F("void_ptr", "p_T Ratio: 1 fb^{-1};p_T-ratio", 300, 0, 10);
    voids[1] = new TH1F("void_m", "Invariant Mass: 1 fb^{-1};M_{e#gamma}", 300, 0, 2000);
    voids[2] = new TH1F("void_Xpt", "Tag and Probe pair boost: 1 fb^{-1};p_{T}", 300, 0, 200);

    TLegend lg(0.6,0.5,0.95,0.8);
    lg.SetFillStyle(0);
    lg.SetTextSize(0.045);
    
    samples["TTjets_ALL"] = make_pair("tt+jets", kGreen+1);
    samples["ZZ4l_ALL"] = make_pair("ZZ", kCyan+1);
    samples["WZjets_ALL"] = make_pair("WZ", kBlue);
    samples["WW2L2nu_ALL"] = make_pair("WW", kBlue);
    samples["WJetsToLNu_HT-600To800_ALL"] = make_pair("W+jets", kRed+1);
    samples["WJetsToLNu_HT-800To1200_ALL"] = make_pair("W+jets", kRed+1);
    samples["WJetsToLNu_HT-1200To2500_ALL"] = make_pair("W+jets", kRed+1);
    samples["WJetsToLNu_HT-2500ToInf_ALL"] = make_pair("W+jets", kRed+1);
    samples["DYLL_all"] = make_pair("DY", kAzure-7);
    samples["singleEle_ALL"] = make_pair("Data", kBlack);

    TH1F* h_sumW;
    float lumiForW=1000;

    float          pu_weight=0;
    float          perEveW=0;
    float          totXsec=0;
    vector<float>* e_pt_ = new vector<float>;
    vector<float>* e_eta_ = new vector<float>;
    vector<float>* e_phi_ = new vector<float>;
    vector<float>* g_pt_ = new vector<float>;
    vector<float>* g_eta_ = new vector<float>;
    vector<float>* g_phi_ = new vector<float>;
    vector<bool>* e_match_ = new vector<bool>;
    vector<int>* g_presel_ = new vector<int>;
    vector<int>* e_idx_ = new vector<int>;
    vector<int>* g_idx_ = new vector<int>;
    vector<float>* invMass_ = new vector<float>;
    vector<float>* ptRatio_ = new vector<float>;

    for(auto& sample : samples)
    {
        cout << "/afs/cern.ch/work/c/crovelli/public/spring15_v3v8/"+sample.first+".root" << endl;
        TFile* file = TFile::Open(("/afs/cern.ch/work/c/crovelli/public/spring15_v3v8/"+sample.first+".root").c_str(), "READ");
        TTree* tree = (TTree*)file->Get("tnpAna/TaPtree");
        h_sumW = (TH1F*)file->Get("tnpAna/h_sumW");
        if(!file || !tree || !h_sumW)
            continue;

        if(!histos[0][sample.second.first])
        {
            histos[0][sample.second.first] = new TH1F((sample.second.first+"_ptr").c_str(), "", 300, 0, 10);
            histos[1][sample.second.first] = new TH1F((sample.second.first+"_m").c_str(), "", 300, 0, 2000);
            histos[2][sample.second.first] = new TH1F((sample.second.first+"_Xpt").c_str(), "", 300, 0, 200);
        }

        double sampleSumWeight = (double)h_sumW->Integral();       
        tree->SetBranchAddress("pu_weight", &pu_weight);
        tree->SetBranchAddress("perEveW", &perEveW);
        tree->SetBranchAddress("totXsec", &totXsec);
        tree->SetBranchAddress("electron_pt", &e_pt_);
        tree->SetBranchAddress("electron_eta", &e_eta_);
        tree->SetBranchAddress("electron_phi", &e_phi_);
        tree->SetBranchAddress("gamma_pt", &g_pt_);
        tree->SetBranchAddress("gamma_eta", &g_eta_);
        tree->SetBranchAddress("gamma_phi", &g_phi_);
        tree->SetBranchAddress("electron_matchHLT", &e_match_);
        tree->SetBranchAddress("gamma_presel", &g_presel_);
        tree->SetBranchAddress("eleIndex", &e_idx_);
        tree->SetBranchAddress("gammaIndex", &g_idx_);
        tree->SetBranchAddress("invMass", &invMass_);
        tree->SetBranchAddress("ptRatio", &ptRatio_);

        for(unsigned int iEntry=0; iEntry<tree->GetEntriesFast(); ++iEntry)
        {
            tree->GetEntry(iEntry);
            float ptr_tmp=-1;
            for(unsigned int iPair=0; iPair<invMass_->size(); ++iPair)
            {
                TLorentzVector ele, gam;
                ele.SetPtEtaPhiM(e_pt_->at(e_idx_->at(iPair)), e_eta_->at(e_idx_->at(iPair)),
                                 e_phi_->at(e_idx_->at(iPair)), 0);
                gam.SetPtEtaPhiM(g_pt_->at(g_idx_->at(iPair)), g_eta_->at(g_idx_->at(iPair)),
                                 g_phi_->at(g_idx_->at(iPair)), 0);
                float X_pt = (ele+gam).Pt();
//.........这里部分代码省略.........

void number_of_events_table_xi1690_cc(TString inFile=""){

	//*** In file

	TFile * input = new TFile(inFile, "READ");

	TTree * ntpMC = (TTree*) input->Get("ntpMC");
	TTree * ntpPiMinus = (TTree*) input->Get("ntpPiMinus");
	TTree * ntpPiPlus = (TTree*) input->Get("ntpPiPlus");
	TTree * ntpkaonplus = (TTree*) input->Get("ntpKaonPlus");
	TTree * ntpProton = (TTree*) input->Get("ntpProton");
	TTree * ntpAntiProton = (TTree*) input->Get("ntpAntiProton");
	TTree * ntpLambda0 = (TTree*) input->Get("ntpLambda0");
	TTree * ntpAntiLambda0 = (TTree*) input->Get("ntpAntiLambda0");
	TTree * ntpXiPlus1690 = (TTree*) input->Get("ntpXiPlus1690");
	TTree * ntpXiMinus = (TTree*) input->Get("ntpXiMinus");
	TTree * ntpXiSys = (TTree*) input->Get("ntpXiSys");


	double nevents_mc = ntpMC->GetEntriesFast();
	TString cuts = " McTruthMatch && VtxFit_HowGood==1 && MassFit_prob>0.01";
	TString VtxCut = " McTruthMatch && VtxFit_HowGood==1 & HitTag==1";
	TString cut4c = "McTruthMatch && 4CFit_prob>0.01";


	cout << "particle|   #evts (uncut)|    #evts (ratio in %)|   MC ratio in %|   dp/p in %" << endl;


	//**** PiMinus(Lambda0)
	TH1D * h_piminus_tht_uncut = new TH1D("h_piminus_tht_uncut", "h_piminus_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus_tht_uncut", "piminus_tht", "McTruthMatch && Mother==3122");
	double piminus_uncut =  h_piminus_tht_uncut->GetEntries();

	TH1D * h_piminus_tht = new TH1D("h_piminus_tht", "h_piminus_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus_tht", "piminus_tht", "McTruthMatch && piminus_HitTag && Mother==3122");
	int piminus =  h_piminus_tht->GetEntries();

	TH1D * h_piminus_dp = new TH1D("h_piminus_dp", "h_piminus_dp", 250, -0.1,0.1);
	ntpPiMinus->Project("h_piminus_dp", "(piminus_p-piminus_MC_p)/piminus_MC_p", "McTruthMatch && piminus_HitTag && Mother==3122");

	Double_t param[6] = jenny::GetFitParameterDoubleFit(h_piminus_dp, false, 0.02,0.1, true);
	//jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piminus_dp, "","", false, false, false, 0.02,0.1, true);


	double ratio_piminus_cut = piminus/piminus_uncut;
	double ratio_piminus_mc = piminus/nevents_mc;

	cout << "PiMinus(L0)|   " <<  piminus_uncut << "|   " <<  piminus << "(" << ratio_piminus_cut*100 << ")|   " << ratio_piminus_mc*100 << endl;//<< "|   " << param[2]*100 << endl;



	//**** piminus (Xi)
	TH1D * h_piminus2_tht_uncut2 = new TH1D("h_piminus2_tht_uncut2", "h_piminus2_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus2_tht_uncut2", "piminus_tht", "McTruthMatch && Mother==3312");
	double piminus_uncut2 =  h_piminus2_tht_uncut2->GetEntries();

	TH1D * h_piminus2_tht = new TH1D("h_piminus2_tht", "h_piminus2_tht", 100, 0,10);
	ntpPiMinus->Project("h_piminus2_tht", "piminus_tht", "McTruthMatch && piminus_HitTag && Mother==3312");
	int piminus2 =  h_piminus2_tht->GetEntries();

	TH1D * h_piminus2_dp = new TH1D("h_piminus2_dp", "h_piminus2_dp", 250, -0.1,0.1);
	ntpPiMinus->Project("h_piminus2_dp", "(piminus_p-piminus_MC_p)/piminus_MC_p", "McTruthMatch && piminus_HitTag && Mother==3312");

	Double_t parampip2[6] = jenny::GetFitParameterDoubleFit(h_piminus2_dp, false, 0.02,0.1, true);
	//jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piminus2_dp, "","", false, false, false, 0.02,0.1, true);

	double ratio_piminus_cut2 = piminus2/piminus_uncut2;
	double ratio_piminus_mc2 = piminus2/nevents_mc;

	cout << "PiMinus(Xi)|   " <<  piminus_uncut2 << "|   " <<  piminus2 << "(" << ratio_piminus_cut2*100 << ")|   " << ratio_piminus_mc2*100<< endl;// << "|   " << parampip2[2]*100 << endl;


	//**** PiPlus (AntiLambda0)
	TH1D * h_piplus_tht_uncut = new TH1D("h_piplus_tht_uncut", "h_piplus_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus_tht_uncut", "piplus_tht", "McTruthMatch && Mother==-3122");
	double piplus_uncut =  h_piplus_tht_uncut->GetEntries();

	TH1D * h_piplus_tht = new TH1D("h_piplus_tht", "h_piplus_tht", 100, 0,10);
	ntpPiPlus->Project("h_piplus_tht", "piplus_tht", "McTruthMatch && piplus_HitTag && Mother==-3122");
	int piplus =  h_piplus_tht->GetEntries();

	TH1D * h_piplus_dp = new TH1D("h_piplus_dp", "h_piplus_dp", 250, -0.1,0.1);
	ntpPiPlus->Project("h_piplus_dp", "(piplus_p-piplus_MC_p)/piplus_MC_p", "McTruthMatch && piplus_HitTag && Mother==-3122");

	Double_t parampip[6] = jenny::GetFitParameterDoubleFit(h_piplus_dp, false, 0.02,0.1, true);
	//jenny::CreateDrawAndSaveHistogramDoulbeFit(h_piplus_dp, "","", false, false, false, 0.02,0.1, true);

	double ratio_piplus_cut = piplus/piplus_uncut;
	double ratio_piplus_mc = piplus/nevents_mc;

	cout << "PiPlus(AL0)|   " <<  piplus_uncut << "|   " <<  piplus << "(" << ratio_piplus_cut*100 << ")|   " << ratio_piplus_mc*100 << endl;//<< "|   " << parampip[2]*100 << endl;




	//**** kaonplus
	TH1D * h_kaonplus_tht_uncut = new TH1D("h_kaonplus_tht_uncut", "h_kaonplus_tht", 100, 0,10);
	ntpkaonplus->Project("h_kaonplus_tht_uncut", "kaonplus_tht", "McTruthMatch && Mother==-13314");
	double kaonplus_uncut =  h_kaonplus_tht_uncut->GetEntries();

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

void simQdep()
{
    for(int x0=1; x0<12; x0++)
    {
	TCanvas* c1 = new TCanvas();
	c1->cd();
	char file_name[60];
	sprintf(file_name, "G4_sim/Pb_FascioDisperso_004/molteplice%d.root", x0);
	TFile* dataFile = TFile::Open("ntuples/reco_Scan4.root","r");
	TTree* nt = (TTree*)dataFile->Get("reco_tree");
	TH1F* hCharge = new TH1F("hCharge", "hCharge", 300, -25000, 2000);
	hCharge->SetLineColor(kRed);

	float thX[8]={113,95,127,118,94,134,133,160};
	float thY[8]={110,101,119,136,94,98,87,97};
	char cut_hodoX[300];
	char cut_hodoY[300];
	TCut cut_sig("-charge_MiB2 > -13.28*amp_max_MiB2 - 350");
        TCut cut_trig1("charge_MiB1 < th_SB_MiB1 && charge_Roma2 < th_SB_Roma2 && sci_front_adc > 500 && sci_front_adc < 1500");
	sprintf(cut_hodoX, "(fibreX[0] > %f || fibreX[1] > %f || fibreX[2] > %f || fibreX[3] > %f || fibreX[4] > %f || fibreX[5] > %f || fibreX[6\
] > %f || fibreX[7] > %f)", thX[0], thX[1], thX[2], thX[3], thX[4], thX[5], thX[6], thX[7]);
	sprintf(cut_hodoY, "(fibreY[0] > %f || fibreY[1] > %f || fibreY[2] > %f || fibreY[3] > %f || fibreY[4] > %f || fibreY[5] > %f || fibreY[6\
] > %f || fibreY[7] > %f)", thY[0], thY[1], thY[2], thY[3], thY[4], thY[5], thY[6], thY[7]);

	nt->Draw("charge_MiB2>>hCharge", cut_trig1 && cut_sig && cut_hodoX && cut_hodoY && "run_id == 261"); 

//	cout << hCharge->GetMean() << endl;

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

	TFile* simFile = TFile::Open(file_name, "r");
	TTree* inTree = (TTree*)simFile->Get("tree");

	TH1I* hMulti = new TH1I("hMulti","hMulti", 10, 0, 10);
	TH1F* hCorr = new TH1F("hCorr","hCorr", 300, -25000, 2000);
//	TH1F* hCorr = new TH1F("hCorr","hCorr", 10, 0, 10);
	inTree->Draw("number_0>>hMulti");	

	std::vector<float>* energy=0;
	inTree->SetBranchAddress("Energies_0", &energy);
	
	float eff=0;
	for(int iEntry=0; iEntry<inTree->GetEntriesFast(); iEntry++)
	{
	    //int count=0;
	    float corrQ=0;
	    inTree->GetEntry(iEntry);
	    //---simulazione molteplicità
	    // for(int iPart=0; iPart<energy->size(); iPart++)
	    // {
	    // 	corrQ += hCharge->GetRandom();
	    // }
	    // hCorr->Fill(corrQ);
	    //---simulazion eff
	    eff += 1-pow(0.55, energy->size());
	}

/*
	for(int iBin=1; iBin<11; iBin++)
	{
	    for(int j=0; j<hMulti->GetBinContent(iBin); j++)
	    {
		int tmp=0;
		for(int i=1; i<iBin; i++)
		{
		    if(gRandom->Uniform(0,1)<0.45)
			tmp++;
		}
		hCorr->Fill(tmp);
	    }
	}
*/
	// hCorr->Scale(hCharge->GetEntries()/hCorr->GetEntries());
	// hCorr->Draw();
	// hCharge->Draw("sames");
	// cout << x0*0.71 << "  " << hCorr->GetMean()/hCharge->GetMean() << endl;
	cout << x0 << "   " << eff/inTree->GetEntriesFast() << endl;
    }
}

Double_t CalibTree::Loop(int loop) {
  char name[500];
  bool debug=false;
  if (fChain == 0) return 0;
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
  std::map<unsigned int, std::pair<double,double> >SumW;
  std::map<unsigned int, unsigned int >nTrks;
  unsigned int mask(0xFF80), ntrkMax(0);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
  //  for (Long64_t jentry=0; jentry<1000;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    if(debug) std::cout << "***Entry (Track) Number : " << ientry 
			<< " p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal 
			<< "/" << t_eMipDR << "/" << (*t_DetIds).size() 
			<< std::endl;
    if (goodTrack()) {
      if (loop == 0) hprof_ndets->Fill(t_ieta, (*t_DetIds).size());
      double Etot=0.0;
      for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) { 
	double hitEn=0.0;
        unsigned int detid = (*t_DetIds)[idet] & mask;
	if (Cprev.find(detid) != Cprev.end()) 
	  hitEn = Cprev[detid] * (*t_HitEnergies)[idet];
	else 
	  hitEn = (*t_HitEnergies)[idet];
	Etot += hitEn;
      }
      for (unsigned int idet=0; idet<(*t_DetIds).size(); idet++) {
	unsigned int detid = (*t_DetIds)[idet] & mask;
	double hitEn=0.0;
	if (debug) std::cout << "idet " << idet << " detid/hitenergy : " 
			     << std::hex << (*t_DetIds)[idet] << ":" << detid
			     << "/" << (*t_HitEnergies)[idet] << std::endl;
	if (Cprev.find(detid) != Cprev.end()) 
	  hitEn = Cprev[detid] * (*t_HitEnergies)[idet];
	else 
	  hitEn = (*t_HitEnergies)[idet];
	double Wi = hitEn/Etot;
	double Fac = (Wi* t_p) / Etot;
	if( SumW.find(detid) != SumW.end() ) {
	  Wi  += SumW[detid].first;
	  Fac += SumW[detid].second;
	  SumW[detid] = std::pair<double,double>(Wi,Fac); 
	  nTrks[detid]++;
	} else {
	  SumW.insert( std::pair<unsigned int, std::pair<double,double> >(detid,std::pair<double,double>(Wi,Fac)));
	  nTrks.insert(std::pair<unsigned int,unsigned int>(detid, 1));
	}
	if (nTrks[detid] > ntrkMax) ntrkMax = nTrks[detid];
      }
    }
  }
  
  std::map<unsigned int, std::pair<double,double> >::iterator SumWItr = SumW.begin();
  unsigned int kount(0), mkount(0);
  double       sumfactor(0);
  double       dets[150], cfacs[150], wfacs[150], nTrk[150];
  unsigned int ntrkCut = ntrkMax/10;
  for (; SumWItr != SumW.end(); SumWItr++) {
    if (debug) 
      std::cout<< "Detid/SumWi/SumFac : " << SumWItr->first << " / "
	       << (SumWItr->second).first << " / " << (SumWItr->second).second
	       << std::endl;
    unsigned int detid = SumWItr->first;
    double factor = (SumWItr->second).second / (SumWItr->second).first;

    if(nTrks[detid]>ntrkCut) {
      if (factor > 1) sumfactor += (1-1/factor);
      else            sumfactor += (1-factor);
      mkount++;
    }
    if (Cprev.find(detid) != Cprev.end()) {
      Cprev[detid] *= factor;
      cfacs[kount] = Cprev[detid];
    } else {
      Cprev.insert( std::pair<unsigned int, double>(detid, factor) );
      cfacs[kount] = factor;
    }
    int ieta = (detid>>7) & 0x3f;
    int zside= (detid&0x2000) ? 1 : -1;
    int depth= (detid>>14)&0x1F;
    wfacs[kount]= factor;
    dets[kount] = zside*(ieta+0.1*(depth-1));
    nTrk[kount] = nTrks[detid];
    kount++;
  }
  TGraph *g_fac, *g_fac2, *g_nTrk;
  g_fac = new TGraph(kount, dets, cfacs); 
  sprintf(name, "Cfacs_detid_it%d", loop);
  fout->WriteTObject(g_fac, name);

  g_fac2 = new TGraph(kount, dets, wfacs); 
  sprintf(name, "Wfacs_detid_it%d", loop);
  fout->WriteTObject(g_fac2, name);

  g_nTrk = new TGraph(kount, dets, nTrk); 
  if (loop==0) fout->WriteTObject(g_nTrk, "nTrk_detid");
//.........这里部分代码省略.........

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;
}

//.........这里部分代码省略.........
      
      //   Int_t nBins = h_rank->GetNbinsX();
      
      //   h_rank = new TH1F("h_rank","EM Rank",nBins,h_rank->GetXaxis()->GetXmin(),h_rank->GetXaxis()->GetXmax());
      //   h_rankIntegrated = new TH1F("h_rankIntegrated","Rates by trigger threshold",nBins,h_rank->GetXaxis()->GetXmin(),h_rank->GetXaxis()->GetXmax());
      std::cout << "making histograms" << std::endl;
      
      h_rank[file] = new TH1F(makeName("h_rank",file),"EM Rank",nBins,xMin,xMax);
      h_rankIntegrated[file] = new TH1F(makeName("h_rankIntegrated",file),"Rates by trigger threshold",nBins,xMin,xMax);
      h_rankBarrel[file] = new TH1F(makeName("h_rankBarrel",file),"EM Rank (|#eta| < 2.1)",nBins,xMin,xMax);
      h_rankBarrelIntegrated[file] = new TH1F(makeName("h_rankBarrelIntegrated",file),"Rates by trigger threshold (|#eta| < 2.1)",nBins,xMin,xMax);
      
      rankCanvas->cd();

      // tree stuff
      emRank = 0;
      emIso = 0;
      emIeta = 0;
      
      tree->SetBranchAddress("emRank", &emRank, &b_emRank);
      tree->SetBranchAddress("emIso", &emIso, &b_emIso);
      tree->SetBranchAddress("emIeta", &emIeta, &b_emIeta);
      b_emRank->SetAutoDelete(kTRUE);
      b_emIso->SetAutoDelete(kTRUE);
      b_emIeta->SetAutoDelete(kTRUE);

      tree->SetBranchStatus("*",0);
      tree->SetBranchStatus("emRank",1);
      tree->SetBranchStatus("emIso",1);
      tree->SetBranchStatus("emIeta",1);

      // loop over tree
      std::cout << "Looping over tree" << std::endl;
      Long64_t nentries = tree->GetEntriesFast();
      for (Long64_t jentry = 0; jentry < nentries; jentry++) 
	{
	  //std::cout << "entry " << jentry << std::endl;
	  //Long64_t ientry = tree->LoadTree(jentry);
	  //tree->LoadTree(jentry); // Set current entry.
	  tree->GetEntry(jentry); // Read all branches of entry and return total number of bytes read.
	  
	  Int_t elecSize = emRank->size();
	  for (Int_t elec = 0; elec < elecSize; elec++)
	    {
	      if (file == 0) // isolated and non-isolated
		{
		  h_rank[file]->Fill(emRank->at(elec));
		  if (emIeta->at(elec) >= barrelIEtaBoundaryLower &&
		      emIeta->at(elec) <= barrelIEtaBoundaryUpper)
		    {
		      h_rankBarrel[file]->Fill(emRank->at(elec));
		    }
		  break;
		}
	      else // isolated only
		{
		  if (emIso->at(elec) == 1)
		    {
		      h_rank[file]->Fill(emRank->at(elec));
		      if (emIeta->at(elec) >= barrelIEtaBoundaryLower &&
			emIeta->at(elec) <= barrelIEtaBoundaryUpper)
			{
			  h_rankBarrel[file]->Fill(emRank->at(elec));
			}
		      break;
		    }