C++ TH1F::Fill方法代码示例

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

  //Chain
  TChain* chain = new TChain("MiBiCommonNTTwoPhotons/SimpleNtuple");
  // Summer 12 DY sample
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_1_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_2_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_3_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_4_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_5_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_6_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_7_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_8_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_9_*.root");
  chain->Add("root://eoscms//eos/cms/store/cmst3/user/malberti/HIGGS/VERTEX/2012/MC/DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball_Summer12-PU_S7_START52_V9-v2_AODSIM/MiBiCommonNT_1*_*.root");



  treeReader reader((TTree*)(chain));

  // histos
  TH1F *hgen = new TH1F("hgen","hgen",60,0,60);
  hgen  -> GetXaxis()->SetTitle("N_{gen} + 1");

  TH1F *hnpumc = new TH1F("hnpumc","hnpumc",60,0,60);
  hnpumc  -> GetXaxis()->SetTitle("N_{gen}");

  TH1F *htruenpumc = new TH1F("htruenpumc","htruenpumc",60,0,60);
  htruenpumc  -> GetXaxis()->SetTitle("N_{gen}");

  TH2F *hnvtxreco_vs_npugen     = new TH2F("hnvtxreco_vs_npugen","hnvtxreco_vs_npugen",60,0,60,60,0,60);
  hnvtxreco_vs_npugen -> GetXaxis()->SetTitle("N_{gen} + 1");
  hnvtxreco_vs_npugen -> GetYaxis()->SetTitle("N_{reco}");

  TH2F *hnpugen_vs_nvtxreco     = new TH2F("hnpugen_vs_nvtxreco","hnpugen_vs_nvtxreco",60,0,60,60,0,60);
  hnpugen_vs_nvtxreco -> GetXaxis()->SetTitle("N_{reco}"); 
  hnpugen_vs_nvtxreco -> GetYaxis()->SetTitle("N_{gen} + 1");
 


  int entryMax = reader.GetEntries();
  std::cout<< "Total number of entries : " << entryMax << std::endl;
  
  // estimate the response curve NPU vs NVXT
  for (int u = 0; u < entryMax; u++ )
    {
      if(u%10000 == 0) std::cout<<"reading event "<< u <<std::endl;
      reader.GetEntry(u);
      
     
      std::vector<int>* mc_PUit_NumInteractions  = reader.GetInt("mc_PUit_NumInteractions");
      int npu = mc_PUit_NumInteractions->at(0);

      std::vector<float>* mc_PUit_TrueNumInteractions  = reader.GetFloat("mc_PUit_TrueNumInteractions");
      int truenpu = mc_PUit_TrueNumInteractions->at(0);

      std::vector<float>* PV_z       = reader.GetFloat("PV_z");
      int nvtx = PV_z->size(); // number of sim PU vtx
      
      hgen->Fill(npu+1);
      hnpumc->Fill(npu);
      htruenpumc->Fill(truenpu);

      hnvtxreco_vs_npugen -> Fill(npu+1,nvtx);
      hnpugen_vs_nvtxreco -> Fill(nvtx,npu+1);
  
    }
  
  TFile ff("weights/HistosForPUReweighting_DYJetsToLL_Summer12_S9-v2.root","recreate");

  hgen->Write();
  hnpumc->Write();
  htruenpumc->Write();
  hnvtxreco_vs_npugen -> Write();
  hnpugen_vs_nvtxreco -> Write();

  ff.Close();
  return 0;
  

}

void MyClass::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L MyClass.C
//      Root > MyClass 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;

    Float_t fJetEtCut = 4.2;
    TH1F * fDebugSVTX = new TH1F ("h1", "", 20, 0, 20);

    fChain->SetBranchStatus("*",0);  // disable all branches
    fChain->SetBranchStatus("Siq2da",1);  // activate branchname

    // create a histogram - for graphical visualization
    TH1F * histo_q2da = new TH1F ("h1", "Q^{2} distribution (double-angle method)", 100, 0, 1000);
    histo_q2da -> SetXTitle ("Q^{2}_{da}, GeV");

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

      // selection of the secondary vertices
       for (int vertex=0; vertex < Kt_njet_a; vertex++) {

            fDebugSVTX->Fill(0);

            // try to match one this jet (from block A) to one of the jets in block B
            // helping variables
            Bool_t      jetAinB=false;
            Int_t       correspjetB=-1;
            // loop over jets in block B
            for (int jetB=0; jetB<Kt_njet_b; jetB++) {
                // check whether Eta, Et and Phi are the same for both jets
                Float_t     deltaEta=Kt_etajet_a[vertex]-Kt_etajet_b[jetB];
                Float_t     deltaPhi=Kt_phijet_a[vertex]-Kt_phijet_b[jetB];
                Float_t     deltaEt=Kt_etjet_a[vertex]-Kt_etjet_b[jetB];
                if ( (deltaEta==0) && (deltaPhi==0) &&  (deltaEt==0) ) {
                    // found a match!
                    jetAinB=true;
                    correspjetB=jetB;
                }
            }

            // skip this block-A-jet if couldn't find a corresponding block-B-jet.
            if ( !jetAinB ) continue;
            fDebugSVTX->Fill(1);

            // JET ET CUT
            if ( Kt_etjet_b[correspjetB] < fJetEtCut) continue;
            fDebugSVTX->Fill(2);

            // JET ETA CUT
            if ( ( Kt_etajet_b[correspjetB] > 2.2 ) || ( Kt_etajet_b[correspjetB] < -1.6 ) ) continue;
            fDebugSVTX->Fill(3);

            // calculate significance of the vertex
            Float_t Significance = CalculateSignificance(vertex);

            // check whether it was successful
            if (Significance==(-999)) continue;
            fDebugSVTX->Fill(4);
            if (Significance==(-998)) continue;
            fDebugSVTX->Fill(5);

            // calculate also projected decay length
            //Float_t     ProjDecayLength = CalculateProjDecayLength(vertex);

            // calculate chi2 over NDOF (in the ntuples we have only chi2)
            Float_t     chi2ndf=Vtxsec_chi2[vertex]/Vtxsec_ndf[vertex];
            // chi2/NDOF cut
            if (chi2ndf>6.) continue;
            fDebugSVTX->Fill(6);

            // create a vertex object to be stored for later processing
    
            // ok, a good vertex was found, store it to array and raise the good-vtx-found flag
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	       ||
	       (dielectron->hltMatchBits_1 & trailingTriggerObjectBit && 
		dielectron->hltMatchBits_2 & leadingTriggerObjectBit ) ) ) continue;
	
	// The Smurf electron ID package is the same as used in HWW analysis
	// and contains cuts like VBTF WP80 for pt>20, VBTF WP70 for pt<10
	// with some customization, plus impact parameter cuts dz and dxy
	if(!passSmurf(dielectron)) continue;  

        // We have a Z candidate! HURRAY! 

	// Apply extra smearing to MC reconstructed dielectron mass
	// to better resemble the data
	/*
	  outdated lines. kept for reference
	double smear1 = escale.extraSmearingSigma(dielectron->scEta_1);
        double smear2 = escale::extraSmearingSigma(dielectron->scEta_2);
	// In systematics mode, overwrite the smear values with
	// shifted ones.
	if(systematicsMode==DYTools::RESOLUTION_STUDY){
	  smear1 = escale::extraSmearingSigmaShifted(dielectron->scEta_1,shift);
	  smear2 = escale::extraSmearingSigmaShifted(dielectron->scEta_2,shift);
	}
        double smearTotal = sqrt(smear1*smear1 + smear2*smear2);
        double massResmeared = dielectron->mass + random.Gaus(0.0,smearTotal);
	*/
	/* lines based on new features -- begin */
	double massSmear = (systematicsMode == DYTools::RESOLUTION_STUDY) ?
	  escale.generateMCSmearRandomized(dielectron->scEta_1,dielectron->scEta_2) :
	  escale.generateMCSmear(dielectron->scEta_1,dielectron->scEta_2);
        double massResmeared = dielectron->mass + massSmear;
	/* lines based on new features -- end */

	hZMassv[ifile]->Fill(massResmeared,scale * gen->weight);

	// Fill structures for response matrix and bin by bin corrections

	if(ibinFsr != -1 && ibinFsr < yieldsMcFsrOfRec.GetNoElements()){
	    yieldsMcFsrOfRec[ibinFsr] += reweight * scale * gen->weight;
	    DetCorrFactorDenominator(ibinFsr) += reweight * scale * gen->weight;
          }

	else if(ibinFsr >= yieldsMcFsrOfRec.GetNoElements())
	  cout << "ERROR: binning problem" << endl;

	int ibinRec = DYTools::findMassBin(massResmeared);
	double shape_weight = (shapeWeights && (ibinRec!=-1)) ? 
	  shapeWeights->GetBinContent(ibinRec+1) : 1.0;
	if(ibinRec != -1 && ibinRec < yieldsMcRec.GetNoElements()){
	    yieldsMcRec[ibinRec] += reweight * scale * gen->weight;
	    DetCorrFactorNumerator(ibinRec) += reweight * scale * gen->weight * shape_weight;
          }

	else if(ibinRec >= yieldsMcRec.GetNoElements())
	  cout << "ERROR: binning problem" << endl;
	
	if( ibinRec != -1 && ibinRec < yieldsMcRec.GetNoElements() 
	    && ibinFsr != -1 && ibinFsr < yieldsMcRec.GetNoElements() ){
          DetMigration(ibinFsr,ibinRec) += reweight * scale * gen->weight * shape_weight;
	}
	
        Bool_t isB1 = (fabs(dielectron->scEta_1)<kGAP_LOW);
        Bool_t isB2 = (fabs(dielectron->scEta_2)<kGAP_LOW);         
	hMassDiff->Fill(massResmeared - gen->mass);
	if( isB1 && isB2 )
	  hMassDiffBB->Fill(massResmeared - gen->mass);

void readTree_background() {
  Char_t *filename = "Background.root";
  // Retrieve the TTree
  TFile* myFile = TFile::Open(filename);
  TTree* tree = (TTree*)(myFile->Get("tree"));
  Double_t Et1, eta1, phi1, Et2, eta2, phi2;
  tree->SetBranchAddress("Et1" ,&Et1 );
  tree->SetBranchAddress("eta1",&eta1);
  tree->SetBranchAddress("phi1",&phi1);
  tree->SetBranchAddress("Et2" ,&Et2 );
  tree->SetBranchAddress("eta2",&eta2);
  tree->SetBranchAddress("phi2",&phi2);
  // Book histograms
  TH1F* hmass = new TH1F("hmass","m_{#gamma#gamma}",60,100.,160.);
  hmass->GetXaxis()->SetTitle("Invariant mass [GeV]");
  hmass->GetYaxis()->SetTitle("Events");
  // Loop over the events
  Long64_t events = tree->GetEntries();  
  for (int i=0; i<events; i++) {
    tree->GetEntry(i);
    TLorentzVector g1,g2;
    g1.SetPtEtaPhiM(Et1,eta1,phi1,0.);
    g2.SetPtEtaPhiM(Et2,eta2,phi2,0.);
    TLorentzVector gg=g1+g2;
    hmass->Fill( gg.M() );
  }
  // Test of different background options
  hmass->DrawClone("e");
  TCanvas* myCanvas = new TCanvas("myCanvas","Background fits",800,800);
  myCanvas->Divide(2,2);
  // Linear background
  TF1* myBack1 = new TF1("myBack1","[0]+[1]*x",100.,160.);
  myBack1->SetParameter(0,events);
  myBack1->SetParameter(1,-100.);
  myCanvas->cd(1);
  hmass->Fit(myBack1);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack1);
  // Quadratic background
  TF1* myBack2 = new TF1("myBack2","[0]+[1]*x+[2]*x**2",100.,160.);
  myBack2->SetParameter(0,events);
  myBack2->SetParameter(1,-100.);
  myBack2->SetParameter(1,0.);
  myCanvas->cd(2);
  hmass->Fit(myBack2);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack2);
  // Exponential background
  TF1* myBack3 = new TF1("myBack3","[0]*exp(-x/[1])",100.,160.);
  myBack3->SetParameter(0,events);
  myBack3->SetParameter(1,100.);
  myCanvas->cd(3);
  hmass->Fit(myBack3);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack3);
  // Cubic background
  TF1* myBack4 = new TF1("myBack4","[0]+[1]*x+[2]*x**2+[3]*x**3",100.,160.);
  myBack4->SetParameter(0,0.);
  myBack4->SetParameter(1,0.);
  myBack4->SetParameter(2,0.);
  myBack4->SetParameter(3,0.);
  myCanvas->cd(4);
  TFitResultPtr fit4 = hmass->Fit(myBack4,"S");
  hmass->DrawClone("e");
  EvaluatePvalue(myBack4);
}

//.........这里部分代码省略.........
      vector<NTMuon> candMuon;


      //Load event for the selection
      sel.LoadEvent(event);

      int finalCut = sel.doFullSelection(&(datasets[d]), channelName, false, /*print*/ verbosity,
	false, false, -1., -1., applyJES,  JESParam,
	 applyEEScale,  EEScaleParam,  applyEEResol,  EEResolParam,
	 applyMEScale,  MEScaleParam,  
	 applyJER,  JERFactor,  applyMETS,  METScale);




        TLorentzVector lvTop1_gen;
        TLorentzVector lvTop2_gen;
        TLorentzVector lvTop1plusTop2_gen;
        float yt_gen = 0.;
        float ytbar_gen = 0.;
        float Delta_y_gen = 0.;

        if(datasets[d].Name() == "TTJets_Sig"){


        lvTop1_gen = (event->topAndDecays)[0].p4_t_gen;
        lvTop2_gen = (event->topAndDecays)[1].p4_t_gen;
        lvTop1plusTop2_gen = lvTop1_gen + lvTop2_gen;

        /////////////////////////
        //// GENERATOR LEVEL ////
        /////////////////////////

        histoGenTopRapidity1->Fill(lvTop1_gen.Rapidity());
        histoGenTopRapidity2->Fill(lvTop2_gen.Rapidity());
        histoGenTopMass1->Fill(lvTop1_gen.M());
        histoGenTopMass2->Fill(lvTop2_gen.M());
        histoGenTopMass1VsTopMass2->Fill(lvTop1_gen.M(),lvTop2_gen.M());

        yt_gen = lvTop1_gen.Rapidity();
        if(yt_gen<0.) yt_gen = (-1)*yt_gen;
        ytbar_gen = lvTop2_gen.Rapidity();
        if(ytbar_gen<0.) ytbar_gen = (-1)*ytbar_gen;
        Delta_y_gen = yt_gen - ytbar_gen;
	//        cout << "Delta_y_gen:" << Delta_y_gen << endl;
        histoGenDelta_y->Fill(Delta_y_gen);
        histoDelta_yEfficiencyD->Fill(Delta_y_gen);
        histoFineBinning_Delta_yEfficiencyD->Fill(Delta_y_gen);

        if(Delta_y_gen>0.){
          N_plus_gen++;
          histoGenN_plusTTRapidity->Fill(lvTop1plusTop2_gen.Rapidity());
          histoGenN_plusTTPt->Fill(lvTop1plusTop2_gen.Perp());
          histoGenN_plusTTMass->Fill(lvTop1plusTop2_gen.M());
          histoFineBinning_GenN_plusTTRapidity->Fill(lvTop1plusTop2_gen.Rapidity());
          histoFineBinning_GenN_plusTTPt->Fill(lvTop1plusTop2_gen.Perp());
          histoFineBinning_GenN_plusTTMass->Fill(lvTop1plusTop2_gen.M());
	}

        if(Delta_y_gen<0.){
          N_minus_gen++;
          histoGenN_minusTTRapidity->Fill(lvTop1plusTop2_gen.Rapidity());
          histoGenN_minusTTPt->Fill(lvTop1plusTop2_gen.Perp());
          histoGenN_minusTTMass->Fill(lvTop1plusTop2_gen.M());
          histoFineBinning_GenN_minusTTRapidity->Fill(lvTop1plusTop2_gen.Rapidity());
          histoFineBinning_GenN_minusTTPt->Fill(lvTop1plusTop2_gen.Perp());

//.........这里部分代码省略.........
    for(int i=0; i<nentries; i++) {
        nt->GetEntry(i);
        if(coll == "PPb" && run>211256) continue;//only for pPb runs at moment
        if(coll == "PbP" && run<=211256) continue;//only for pPb runs at moment
        if((TMath::Abs(vz)>15) || (!pPAcollisionEventSelectionPA) || (!pprimaryVertexFilter)|| (!pHBHENoiseFilter)) continue;
        if(!HLT_PAJet20_noJetID_v1 && !HLT_PAJet40_noJetID_v1 && !HLT_PAJet60_noJetID_v1 && !HLT_PAJet80_noJetID_v1 && !HLT_PAJet100_noJetID_v1 ) continue;

        if(i%10000==1)    cout<<"analyzing "<< i <<" th event"<<endl;
        for(int j4i = 0; j4i < nref; j4i++) {
            double jet_pt = jtpt[j4i];
            double jet_eta = jteta[j4i];
            double jet_phi = jtphi[j4i];
            double raw_pt = rawpt[j4i];
            double jetweight = 1.;
//       jetweight*=(fUE->Eval(jet_pt))*C_rel->GetBinContent(C_rel->FindBin(jet_eta));
            jetweight*=C_rel->GetBinContent(C_rel->FindBin(jet_eta));

            double chargedMax = t_chargedMax[j4i];
            double chargedSum = t_chargedSum[j4i];
            double neutralMax = t_neutralMax[j4i];
            double neutralSum = t_neutralSum[j4i];
            double photonMax = t_photonMax[j4i];
            double photonSum = t_photonSum[j4i];
            int chargedN = t_chargedN[j4i];
            int neutralN = t_neutralN[j4i];
            int photonN = t_photonN[j4i];
            double muSum = t_muSum[j4i];
            double eSum = t_eSum[j4i];
            //      if((chargedN == 0 || chargedSum == 0) && TMath::Abs(jet_eta)< 2.4) continue;  // jet id selection
            double PPTighter0 = (double)(neutralSum/jet_pt < 0.8 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 1.0 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4));
            double PPTighter1 = (double)(neutralSum/jet_pt < 0.8 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 1.0 && chargedSum/jet_pt < 0.95 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) );
            double PPTighter2 = (double)(neutralSum/jet_pt < 0.9 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 1.0 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) );
            double PPTighter3 = (double)(neutralSum/jet_pt < 0.9 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 1.0 && chargedSum/jet_pt < 0.95 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) );
            double PPTighter4 = (double)(neutralSum/jet_pt < 0.9 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 0.95 && chargedSum/jet_pt < 0.95 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) );
            double PPTighter5 = (double)(neutralSum/jet_pt < 0.8 && eSum/jet_pt < 1.0 && photonSum/jet_pt < 1.0 && (chargedSum+neutralSum+muSum+eSum)/jet_pt < 0.95 && chargedSum/jet_pt < 0.95 && ((chargedSum/jet_pt>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) );
            double PPTighter = PPTighter0*TMath::Power(2,5)+PPTighter1*TMath::Power(2,4)+PPTighter2*TMath::Power(2,3)+PPTighter3*TMath::Power(2,2)+PPTighter4*TMath::Power(2,1)+PPTighter5*TMath::Power(2,0)+0.5;
            double jetidv[nJetID]= {chargedMax,chargedSum,neutralMax,neutralSum,photonMax,photonSum,chargedMax/jet_pt,chargedSum/jet_pt,neutralMax/jet_pt,neutralSum/jet_pt,photonMax/jet_pt,photonSum/jet_pt,eSum/jet_pt,(chargedSum+neutralSum+photonSum+muSum+eSum)/jet_pt,(chargedSum+neutralSum+photonSum+muSum+eSum)/raw_pt,neutralMax/TMath::Max(chargedSum,neutralSum),(double)chargedN,(double)neutralN,(double)photonN,(double)(neutralSum/jet_pt<1.0 && eSum/jet_pt<1.0 && photonSum/jet_pt<1.0 && ((chargedSum>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) ), (double)(neutralSum/jet_pt<0.9 && eSum/jet_pt<1.0 && photonSum/jet_pt<0.9 && ((chargedSum>0 && chargedN>0 && TMath::Abs(jet_eta)<2.4) || TMath::Abs(jet_eta) >=2.4) ),PPTighter,(chargedSum+neutralSum+muSum+eSum)/jet_pt};
            if(raw_pt<22 || fabs(jet_eta)>5) continue;
            if(jet_pt>4*pt) continue;
            //if(neutralSum/jet_pt > 0.8 && chargedSum/jet_pt < 0.1) {
            if(neutralSum/jet_pt > 0.8) {
                Spikeetaphi->Fill(jet_eta,jet_phi, weight);
                Spikeetapt->Fill(jet_eta,jet_pt, weight);
            }

            jet_pt = jet_pt*jetweight;

            int dEtaBin = -1;
            if(coll=="PPb")	jet_eta=jet_eta+0.465;
            if(coll=="PbP")	jet_eta=jet_eta-0.465;


            if(TMath::Abs(jet_eta)<=1.) {
                jetpt->Fill(jet_pt,weight);
                for(int ijetid=0; ijetid<nJetID; ijetid++) {
                    jetptjetid[ijetid]->Fill(jet_pt, jetidv[ijetid],weight);     //Added
                }
            }
            for(int ieta =0; ieta<netabin; ieta++) {
                if(jet_eta>deta[ieta]&&jet_eta<=deta[ieta+1]) dEtaBin=ieta;
            }//assign the eta bin for jets
            if(dEtaBin!=-1) {
                jetptEtaBin[dEtaBin]->Fill(jet_pt,weight);
                for(int ijetid=0; ijetid<nJetID; ijetid++) {
                    if(JetIDName[ijetid].Contains("pt") || JetIDName[ijetid].Contains("Maxr") || JetIDName[ijetid].Contains("PP") || JetIDName[ijetid].Contains("N")) {
                        jetptjetidEtaBin[dEtaBin][ijetid]->Fill(jet_pt,jetidv[ijetid],weight);
                    }
                }
            }
        } //loop over jet

    }	//loop over tree

    TString dataType;
    TString out_name;
    dataType = "DATA";
    out_name=Form("%s%s%sskimJetID.root",dataType.Data(),coll.Data(),algo.Data());

    TFile *out_file = new TFile(Form("/tmp/xuq7/%s",out_name.Data()),"RECREATE");

    jetpt->Write();
    Spikeetaphi->Write();
    Spikeetapt->Write();
    for(int ijetid=0; ijetid<nJetID; ijetid++) {
        jetptjetid[ijetid]->Write();
    }
    for(int ieta=0; ieta<netabin; ieta++) {
        for(int ijetid=0; ijetid<nJetID; ijetid++) {
            if(JetIDName[ijetid].Contains("pt") || JetIDName[ijetid].Contains("Maxr") || JetIDName[ijetid].Contains("PP") || JetIDName[ijetid].Contains("N")) {
                jetptjetidEtaBin[ieta][ijetid]->Write();
            }
        }
        jetptEtaBin[ieta]->Write();
    }
    out_file->Close();
    cout<<"Output file: "<<Form("%s",out_name.Data())<<endl;

    cout<<"working done\n";

}

//.........这里部分代码省略.........
   TH1F* hTrgRecoVsRecoEleEB;
   TH1F* hTrgRecoVsRecoEleEE;
   TH1F* hTrgRecoVsRecoMu;

   // output file
   stringstream ssOutfile;
   ssOutfile << outfileName << ".root";
   TFile *output = new TFile(ssOutfile.str().c_str(), "recreate");

   ///////////////////////////////////////////////////////////////////////////
   //LOOP OVER FILES 
   ///////////////////////////////////////////////////////////////////////////
   for (unsigned int p = 0; p < files.size(); ++p) {
      TFile* input = new TFile(files[p]);
      TTree *thetree = (TTree*)input->Get("gsfcheckerjob/tree");
      Init(thetree);
      Long64_t nentries = fChain->GetEntriesFast();
      cout << nentries << " events" << endl;

      unsigned int evCounter = 0;
      /////////////////////////////////////////////////////////////////////////////////////////////
      //LOOP OVER EVENTS
      /////////////////////////////////////////////////////////////////////////////////////////////
      //nentries = 10000;
      for (Long64_t jentry = 0; jentry < nentries; ++jentry) {
         Long64_t ientry = LoadTree(jentry);
         if (ientry < 0) break;
         fChain->GetEntry(jentry);
         // if (Cut(ientry) < 0) continue;
         //if (jentry % 50000 == 0) cout << "Processing event " << jentry << endl;
         thetree->GetEntry(jentry);

         // fill the gen histograms
         hGenEvts->Fill(genelemom_mass[0]);

         // fill the acc histograms
         if (hardGenEle_pt[0] > elePtCut) {
            if (fabs(hardGenEle_eta[0]) < 1.442) {
               hGenEvtsEleInAcc->Fill(genelemom_mass[0]);
               hGenEvtsEleInAccEB->Fill(genelemom_mass[0]);
            }
            else if (fabs(hardGenEle_eta[0]) > 1.56 && fabs(hardGenEle_eta[0]) < 2.5) {
               hGenEvtsEleInAcc->Fill(genelemom_mass[0]);
               hGenEvtsEleInAccEE->Fill(genelemom_mass[0]);
            }
         }
         if (hardGenMu_pt[0] > muPtCut && fabs(hardGenMu_eta[0]) < 2.4) {
            hGenEvtsMuInAcc->Fill(genelemom_mass[0]);
            if (fabs(hardGenEle_eta[0]) < 1.442 || (fabs(hardGenEle_eta[0]) > 1.56 && fabs(hardGenEle_eta[0]) < 2.5)) {
               if (hardGenEle_pt[0] > elePtCut) hGenEvtsInAcc->Fill(genelemom_mass[0]);
            }
         }

         // trigger?
         if (HLT_Mu22_Photon22_CaloIdL) hTrgEvts->Fill(genelemom_mass[0]);

         // at least one gsf electron and one muon above the threshold
         if (gsf_size < 1 || muon_size < 1) continue;

         vector<int> GSF_passHEEP;
         vector<int> MU_passGOOD;
         /////////////////////////////////////////////////////////////////////////////////////////////
         //loop over electrons
         for (int j = 0; j < gsf_size; ++j) {
            //cleaning: fake electrons from muons
            bool fakeEle = false;

//.........这里部分代码省略.........
	cout << "zero yield, skipping!" << endl;
	continue; 
      }

      //--------------------------
      //uncertainty from k-factor
      //--------------------------

      float kerr_up   = fabs( ( yield_kup - yield_k ) / yield_k );
      float kerr_dn   = fabs( ( yield_kdn - yield_k ) / yield_k );
      float kerr      = TMath::Max( kerr_up , kerr_dn );

      //--------------------------     
      //uncertainty from JES
      //--------------------------

      float jerr_up   = fabs( ( yield_jup - yield_k ) / yield_k );
      float jerr_dn   = fabs( ( yield_jdn - yield_k ) / yield_k );
      float jerr      = TMath::Max( jerr_up , jerr_dn );
     
      //-----------------------------------------------------------
      //add up NLO uncertainties (including k-factor uncertainty)
      //-----------------------------------------------------------

      float err2_NLO = 0.;
      err2_NLO += kerr * kerr;                                              //k-factor
      err2_NLO += jerr * jerr;                                              //JES
      err2_NLO += lumierr * lumierr;                                        //lumi
      err2_NLO += leperr * leperr;                                          //lep efficiency
      err2_NLO += pdferr * pdferr;                                          //PDF uncertainty
      accerr_NLO[m0bin-1][m12bin-1] = err2_NLO > 0 ? sqrt( err2_NLO ) : 0.; //total

      if( yield_k > 5 && yield_k < 8 )
	htotuncertainty->Fill( accerr_NLO[m0bin-1][m12bin-1] );

      ul_NLO[m0bin-1][m12bin-1]       = getUpperLimit          ( accerr_NLO[m0bin-1][m12bin-1] );
      ul_NLO_exp[m0bin-1][m12bin-1]   = getExpectedUpperLimit  ( accerr_NLO[m0bin-1][m12bin-1] );
      ul_NLO_expp1[m0bin-1][m12bin-1] = getExpectedP1UpperLimit( accerr_NLO[m0bin-1][m12bin-1] );
      ul_NLO_expm1[m0bin-1][m12bin-1] = getExpectedM1UpperLimit( accerr_NLO[m0bin-1][m12bin-1] );

      //--------------------------
      //printout errors and UL
      //--------------------------
     
      //cout << "yield            " << yield   << endl;
      cout << "yield * K        " << yield_k << endl;
      cout << "yield * Kup      " << yield_kup << endl;
      cout << "yield * Kdn      " << yield_kdn << endl;
      cout << "yield * JESup    " << yield_jup << endl;
      cout << "yield * JESdn    " << yield_jdn << endl;
      cout << "K error          " << kerr << endl;
      cout << "JES error        " << jerr << endl;
      cout << "total error NLO  " << accerr_NLO[m0bin-1][m12bin-1] << endl;
      cout << "NLO UL           " << ul_NLO[m0bin-1][m12bin-1] << endl;
      cout << "NLO UL sig cont  " << ul_NLO_SC[m0bin-1][m12bin-1] << endl;
      cout << "NLO UL no bkg    " << ul_NLO_nobkg[m0bin-1][m12bin-1] << endl;
      cout << "NLO UL exp       " << ul_NLO_exp[m0bin-1][m12bin-1] << endl;
     
      hUL_NLO->SetBinContent( m0bin , m12bin , ul_NLO[m0bin-1][m12bin-1] );
      hUL_NLO_exp->SetBinContent( m0bin , m12bin , ul_NLO_exp[m0bin-1][m12bin-1] );
      hUL_NLO_expp1->SetBinContent( m0bin , m12bin , ul_NLO_expp1[m0bin-1][m12bin-1] );
      hUL_NLO_expm1->SetBinContent( m0bin , m12bin , ul_NLO_expm1[m0bin-1][m12bin-1] );
    }
  }

  cout << __LINE__ << endl;

void DrawCalibrationPlotsEB( Char_t* infile1 = "/data1/rgerosa/L3_Weight/PromptSkim_Single_Double_Electron_recoFlag/EB/WZAnalysis_PromptSkim_W-DoubleElectron_FT_R_42_V21B_Z_noEP.root",
			    Char_t* infile2 = "/data1/rgerosa/L3_Weight/PromptSkim_Single_Double_Electron_recoFlag/EB/Even_WZAnalysis_PromptSkim_W-DoubleElectron_FT_R_42_V21B_Z_noEP.root",
			    Char_t* infile3 = "/data1/rgerosa/L3_Weight/PromptSkim_Single_Double_Electron_recoFlag/EB/Odd_WZAnalysis_PromptSkim_W-DoubleElectron_FT_R_42_V21B_Z_noEP.root",
			    int evalStat = 0,
			    Char_t* fileType = "png", 
			    Char_t* dirName = ".")
{

  bool  printPlots = false;

  // by TT
  int nbins = 500;

  // by xtal
  //int nbins = 500;

  // Set style options
  gROOT->Reset();
  gROOT->SetStyle("Plain");

  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetOptTitle(0); 
  gStyle->SetOptStat(1110); 
  gStyle->SetOptFit(0); 
  gStyle->SetFitFormat("6.3g"); 
  gStyle->SetPalette(1); 
 
  gStyle->SetTextFont(42);
  gStyle->SetTextSize(0.05);
  gStyle->SetTitleFont(42,"xyz");
  gStyle->SetTitleSize(0.05);
  gStyle->SetLabelFont(42,"xyz");
  gStyle->SetLabelSize(0.05);
  gStyle->SetTitleXOffset(0.8);
  gStyle->SetTitleYOffset(1.1);
  gROOT->ForceStyle();

  if ( !infile1 ) {
    cout << " No input file specified !" << endl;
    return;
  }


  if ( evalStat && (!infile2 || !infile3 )){
    cout << " No input files to evaluate statistical precision specified !" << endl;
    return;
  }

  cout << "Making calibration plots for: " << infile1 << endl;
  
  TFile *f = new TFile(infile1);
  TH2F *h_scale_EB = (TH2F*)f->Get("h_scale_EB");
  TH2F *hcmap = (TH2F*) h_scale_EB->Clone("hcmap");
  hcmap -> Reset("ICEMS");

  // Mean over phi 

  for (int iEta = 1 ; iEta < h_scale_EB->GetNbinsY()+1 ; iEta ++)
  {
   float SumIC = 0;
   int numIC = 0;
   
   for(int iPhi = 1 ; iPhi < h_scale_EB->GetNbinsX()+1 ; iPhi++)
   {
    if( h_scale_EB->GetBinContent(iPhi,iEta) !=0)
    {
      SumIC = SumIC + h_scale_EB->GetBinContent(iPhi,iEta);
      numIC ++ ;
    }
   }
   
   for (int iPhi = 1; iPhi< h_scale_EB->GetNbinsX()+1  ; iPhi++)
   {
    if(numIC!=0 && SumIC!=0)
    hcmap->SetBinContent(iPhi,iEta,h_scale_EB->GetBinContent(iPhi,iEta)/(SumIC/numIC));
   }
  }
  

  
  //-----------------------------------------------------------------
  //--- Build the precision vs ieta plot starting from the TH2F of IC
  //-----------------------------------------------------------------
  TH1F *hoccall = new TH1F("hoccall", "hoccall", 1000,0.,1000.);
  for (int jbin = 1; jbin < h_occupancy-> GetNbinsY()+1; jbin++){
    for (int ibin = 1; ibin < h_occupancy-> GetNbinsX()+1; ibin++){
      float ic = h_occupancy->GetBinContent(ibin,jbin);
      	hoccall->Fill(ic);
    }
  }

  TH1F *hspreadall = new TH1F("hspreadall", "hspreadall", 800,0.,2.);
  for (int jbin = 1; jbin < hcmap-> GetNbinsY()+1; jbin++){
    for (int ibin = 1; ibin < hcmap-> GetNbinsX()+1; ibin++){
      float ic = hcmap->GetBinContent(ibin,jbin);
      if (ic>0 && ic<2 && ic !=1)    {
	hspreadall->Fill(ic);
      }
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            
            for(int ir=0;ir<25;ir++){
                //! Response vs DeltaR
                hRspVsDeltaR[nj][icen][ir] = new TH1F(Form("hRspVsDeltaR%d_%d_%d",nj,icen,ir),Form(" <recopt/refpt> vs. #DeltaR (%d) algorithm %s cent %d",ir,cjets[nj],icen),rbins,rbinl,rbinh);
            }
        }//! icen
        //for eta clsoure in different jet pt bins
        for(int ipt=0;ipt<bins;ipt++){
            hratiocorrrefpt_etaptbin[nj][ipt]= new TH2F(Form("hratiocorrrefpt_etaptbin%d_%d",nj,ipt),Form("Gen matched jet Reco jet / Gen jet p_{T} (corr.) distribution jet vs eta in pt bin %d %s",ipt,cjets[nj]),
                                                  100,-5.0,5.0,rbins,rbinl,rbinh);
            
        }
        //for different HFplusEta4 bins
        for(int ihf=0;ihf<nhfbin;ihf++){
            hratiocorrrefpt_genhfb[nj][ihf]= new TH2F(Form("hratiocorrrefpt_genhfbin%d_%d",nj,ihf),Form("Gen matched jet Reco jet / Gen jet p_{T} (corr.) distribution jet hiHF bin %d %s",ihf,cjets[nj]),
                                                      500,0,1000,rbins,rbinl,rbinh);
        }
        
    }//! nj
  std::cout<<"Initialized the histograms " <<std::endl;
  ///////////////////////////////////////////////////////////////////////////////////////// 


  //! Centrality reweighting function

  //! vertex z reweighting




  Long64_t nentries = c->GetEntries();
  std::cout<<Form("# of entries in TTree for %s : ",ksp)<<nentries<<std::endl;
  std::cout<<std::endl;
  hEvt->Fill(2,nentries);

  //! weight  for the merging of the samples for different pT hat bins
  Float_t wxs = xsection/(nentries/100000.);
  
  Int_t iEvent=0; 
  for (Long64_t ievt=0; ievt<nentries;ievt++) {//! event loop
    //for (Long64_t ievt=0; ievt<5000;ievt++) {//! event loop
    //! load the hiForest event
    c->GetEntry(ievt);

    int hiBin       = c->evt.hiBin;
    float vz        = c->evt.vz;
    float hiHF      = c->evt.hiHF;
    int ntracks     = c->evt.hiNtracks;
    float HFplusEta = c->evt.hiHFplusEta4;
    float HFminusEta = c->evt.hiHFminusEta4 ; 
    float HFplusEta4 = HFplusEta+HFminusEta;
    //! testing
    //if(hiBin>4 && strcmp(ksp,"pbpb")==0)continue;

    //! apply vertex cut
    if(fabs(vz)>kVzcut)continue;


    //! Centrality bin
    if(hiBin<0 || hiBin>100)continue;
    int multb=GetMultBin(ntracks);
    int hiHFb=GetHFplusEta4Bin(HFplusEta4);  
      
    double wcen=1;
    double wvz=1;
    //wxs=1;

void DoubleHiggsAnalysis::Analyze(){

  Long64_t numberOfEntries = treeReader_->GetEntries();

  // Get pointers to branches used in this analysis
  TClonesArray *branchJet = treeReader_->UseBranch("Jet");
  TClonesArray *branchPhoton = treeReader_->UseBranch("Photon");
  TClonesArray *branchElectron = treeReader_->UseBranch("Electron");
  TClonesArray *branchMuon = treeReader_->UseBranch("Muon");

  TH1F *histPhoMass = new TH1F("mass", "M_{inv}(#gamma#gamma)", 100, 100.0, 180.0);

  //setXsec
  getProcessXsec();
  setGenEvents(numberOfEntries);
  setEventWeight();
  eventWeight_t=eventWeight_;

  cout<<"starting loop on "<<numberOfEntries<<" entries"<<endl;
  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
    {
      if(entry%300 == 0)cout<<"### Processing entry: "<<entry<<endl; 
    // Load selected branches with data from specified event
      treeReader_->ReadEntry(entry);
      bool passedJetSelection=false;
      bool passedPhotonSelection=false;
      bool passedLeptonSelection=false;      

      bool looseBtagWP=false;

      Photon *pho1, *pho2;
      // Filling photon variables
      pho1 = (Photon *) branchPhoton->At(0);
      pho2 = (Photon *) branchPhoton->At(1);
      
      passedPhotonSelection=PhotonSelection(branchPhoton);
      if (passedPhotonSelection){
	counters_photonSel_++;
	passedJetSelection=JetSelection(branchJet, looseBtagWP,pho1,pho2);
	passedLeptonSelection=LeptonSelection(branchElectron,branchMuon,pho1,pho2);
      }
      
      //filling histograms and tree
      if(passedJetSelection && passedPhotonSelection){
	counters_jetSel_++;

	ptPhot1_t=pho1->PT;
	ptPhot2_t=pho2->PT;
	etaPhot1_t=pho1->Eta;
	etaPhot2_t=pho2->Eta;
	if(etaPhot2_t == 0)cout<<"--------"<<endl;
	phiPhot1_t=pho1->Phi;
	phiPhot2_t=pho2->Phi;
	histPhoMass->Fill(((pho1->P4()) + (pho2->P4())).M());

	//Filling jet variables
	jet1 = (Jet*) branchJet->At(bjet_indexes.first);
	jet2 = (Jet*) branchJet->At(bjet_indexes.second);

	ptJet_t[0]=jet1->PT;
	ptJet_t[1]=jet2->PT;
	etaJet_t[0]=jet1->Eta;
	etaJet_t[1]=jet2->Eta;
	phiJet_t[0]=jet1->Phi;
	phiJet_t[1]=jet2->Phi;

	deltaRGammaGamma_t=(pho1->P4().DeltaR(pho2->P4()));
	deltaRBB_t=(jet1->P4().DeltaR(jet2->P4()));
	minDeltaRGammaB_t=(pho1->P4().DeltaR(jet1->P4()));
	if(pho1->P4().DeltaR(jet2->P4()) < minDeltaRGammaB_t )minDeltaRGammaB_t=pho1->P4().DeltaR(jet2->P4());
	if(pho2->P4().DeltaR(jet1->P4()) < minDeltaRGammaB_t )minDeltaRGammaB_t=pho2->P4().DeltaR(jet1->P4());
	if(pho2->P4().DeltaR(jet2->P4()) < minDeltaRGammaB_t )minDeltaRGammaB_t=pho2->P4().DeltaR(jet2->P4());
	ptGG_t=(pho1->P4()+pho2->P4()).Pt();
	ptBB_t=(jet1->P4()+jet2->P4()).Pt();
	ptBBGG_t=ptBB_t+ptGG_t;
	mgg_t=((pho1->P4()) + (pho2->P4())).M();
	mbb_t=(jet1->P4()+jet2->P4()).M();
	mggbb_t=((pho1->P4()+pho2->P4())+(jet1->P4()+jet2->P4())).M();

	bool passedAdditionalCuts=false;
	float deltaRGGcut=2.;
	float deltaRGBcut=1.5;
	float deltaRBBcut=2.;
	int njetsCut=4;
	passedAdditionalCuts=additionalCuts(deltaRGGcut,deltaRGBcut,deltaRBBcut,njetsCut);
	if(passedAdditionalCuts){
	  counters_additionalCuts_++;
	  tree_passedEvents->Fill();
	  if((mgg_t>120. && mgg_t<130.) && (mbb_t >105. && mbb_t<145.))counters_massWindow_++;
	}
      }
    }
  
  histPhoMass->Write();
  tree_passedEvents->Write();
  outFile_->Write();
  outFile_->Close();

}

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

   //Make input branches of variables to use for the training
   simu->Branch("dEta", &deta, "dEta/D");
   simu->Branch("dPhi",  &dphi,  "dPhi/D");
   simu->Branch("Angle",    &angle,    "Angle/D");
   simu->Branch("pBalance",  &pbalance,  "pBalance/D");
   simu->Branch("EtaH", &etah, "EtaH/D");
   simu->Branch("Sphericity",  &sphericity,  "Sphericity/D");
   simu->Branch("type",   &typeO,   "type/I");
  
   //loop over signal and select events to use for training 
   Int_t i;
   for (i = 0; i < signal->GetEntries(); i++) {
     signal->GetEntry(i);
     //only consider events with 3 jets
      if (Njets==3){
          typeO=0; 
          //select only events where the two leading pT jets come from the 
	  //Higgs decay (typeH[0]==1), and these two jets are separated by 
	  //sqrt(dEta[0]*dEta[0]+dPhi[0]*dPhi[0]) > 1, and the invariant 
          //di-jet mass is greater than 50 GeV
          if (typeH[0]==1 && sqrt(dEta[0]*dEta[0]+dPhi[0]*dPhi[0]) > 1. && MH[0]>50.){          
	    deta = dEta[0];
	    dphi =2*acos(0)- dPhi[0]; 
            angle = Angle[0];
            pbalance = pBalance[0];
            etah = EtaH[0];
            sphericity = Sphericity;
            typeO=1; 
	   }
	  //fill selected signal events (type0==1)
          if (typeO==1){
	    //fill the histograms like this:
	      detas->Fill(deta); 	      
	      dphis->Fill(dphi);
	      angles->Fill(angle);
	      pbalances->Fill(pbalance);
	      etahs->Fill(etah);
	      sphericitys->Fill(sphericity);


	      //fill the training tree   
              simu->Fill();
              nsignal++; 
	  }
      }
   }

   cout << nsignal <<endl; 
   //now loop over background events
   for (i = 0; i < background->GetEntries(); i++) {
      background->GetEntry(i);
      if (Njets==3){
          typeO=1;
	  //similar for background events
          if (sqrt(dEta[0]*dEta[0]+dPhi[0]*dPhi[0]) > 1. && MH[0]>50.){
	    deta = dEta[0];
	    dphi =2*acos(0)- dPhi[0]; 
            angle = Angle[0];
            pbalance = pBalance[0];
            etah = EtaH[0];
            sphericity = Sphericity;
            typeO=0; 
          }
	  //fill the selected background events until you have the same 
          //number as signal

void pmt_testing(char* input_file) {

  Double_t width = 1000;
  Double_t height = 1000;
  TCanvas * canvas = new TCanvas("canvas", "PMT Testing", 0, 0, width, height);
  canvas->SetWindowSize(width + (width - canvas->GetWw()), 
                        height + (height - canvas->GetWh()));

  //Parameters
  Int_t nbins = 40;
  Double_t min = 550000; //hist min/max values
  Double_t max = 600000; //Must be better way to get these!

  //Initialized variables
  Int_t num_data_points = 0;
  Int_t pulse;
  Int_t channel;
  Int_t total_channels;
  const Int_t pedestal = 2020; //how to determine this?

  TFile *file = new TFile("adc_data.root", "recreate"); 

  //Style settings
  gROOT->SetStyle("Plain");
  gStyle->SetPalette(53);
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit(1111);
  gStyle->SetStatBorderSize(0);
  gStyle->SetOptTitle(1); 

  
  //Prepare Data
  //Create nTuple
  TTree *t1 = new TTree("t1", "ADC Data");
  t1->Branch("pulse", &pulse, "pulse/I");
  t1->Branch("channel", &channel, "channel/I");

 
  ifstream in;
  in.open(input_file);
  //Fill TTree
  while ( in >> channel) {
    if (num_data_points % 2520 == 0)
      pulse++;  
    t1->Fill();
    num_data_points++;
  }
  in.close();

  TH1F * spectrum = new TH1F("spectrum", "ADC Spectrum; Accumulated Channels; Counts", nbins, min, max); 
  //spectrum->SetBit(TH1::kCanRebin);
  spectrum->Sumw2(); 
  for(Int_t i = 0; i < t1->GetEntries(); i++) {
    t1->GetEntry(i);
    if (i % 2520 == 0) {
      cout << "total_channels = " << total_channels << endl;
      spectrum->Fill(total_channels);
      total_channels = 0;
    }
    total_channels += TMath::Abs(channel - pedestal);
  }

  spectrum->SetLineWidth(2);
  spectrum->SetMarkerColor(kBlack);
  spectrum->SetLineColor(kBlack);

  // Fit Gaussian to curve and Draw
  spectrum->Fit("gaus");
  spectrum->GetFunction("gaus")->SetLineColor(kRed);
  spectrum->Draw("Same");
  t1->Write();
  spectrum->Write();
}

int
Proto2ShowerCalib::process_event(PHCompositeNode *topNode)
{

  if (verbosity > 2)
    cout << "Proto2ShowerCalib::process_event() entered" << endl;

  // init eval objects
  _eval_run.reset();
  _eval_3x3_raw.reset();
  _eval_5x5_raw.reset();
  _eval_3x3_prod.reset();
  _eval_5x5_prod.reset();
  _eval_3x3_temp.reset();
  _eval_5x5_temp.reset();
  _eval_3x3_recalib.reset();
  _eval_5x5_recalib.reset();

  Fun4AllHistoManager *hm = get_HistoManager();
  assert(hm);

  if (not _is_sim)
    {
      PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
          "RUN_INFO");

      assert(info);

      PHG4Parameters run_info_copy("RunInfo");
      run_info_copy.FillFrom(info);

      _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");

      TH1F * hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
      assert(hBeam_Mom);

      hBeam_Mom->Fill(_eval_run.beam_mom);
    }

  EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
      "EventHeader");
  if (not _is_sim)
    {
      assert(eventheader);

      _eval_run.run = eventheader->get_RunNumber();
      if (verbosity > 4)
        cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;

      _eval_run.event = eventheader->get_EvtSequence();
    }

  if (_is_sim)
    {

      PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
          PHG4TruthInfoContainer>(topNode, "G4TruthInfo");

      assert(truthInfoList);

      _eval_run.run = -1;

      const PHG4Particle * p = truthInfoList->GetPrimaryParticleRange().first->second;
      assert(p);

      const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
      assert(v);

      _eval_run.beam_mom = sqrt(
          p->get_px() * p->get_px() + p->get_py() * p->get_py()
              + p->get_pz() * p->get_pz());
      _eval_run.truth_y = v->get_y();
      _eval_run.truth_z = v->get_z();

    }

  // normalization
  TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
  assert(hNormalization);

  hNormalization->Fill("ALL", 1);

  RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
      topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
  assert(TOWER_RAW_CEMC);
  RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
      topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
  assert(TOWER_CALIB_CEMC);

  // other nodes
  RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
      RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
  if (not _is_sim)
    {
      assert(TOWER_CALIB_TRIGGER_VETO);
    }

  RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
      RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
  if (not _is_sim)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    
             
    // loop on towers
    for (unsigned int tower = 0 ; tower < treeVars.nbOfTowers ; tower++) {

      int tp = getEt(treeVars.rawTPData[tower]) ;
      int emul[5] = {getEt(treeVars.rawTPEmul1[tower]),  
		     getEt(treeVars.rawTPEmul2[tower]),
		     getEt(treeVars.rawTPEmul3[tower]),
		     getEt(treeVars.rawTPEmul4[tower]),
		     getEt(treeVars.rawTPEmul5[tower])} ;
      int maxOfTPEmul = 0 ;
      int indexOfTPEmulMax = -1 ;
      for (int i=0 ; i<5 ; i++) if (emul[i]>maxOfTPEmul) {
	maxOfTPEmul = emul[i] ; 
	indexOfTPEmulMax = i ;
      }
      int ieta = treeVars.ieta[tower] ;
      int iphi = treeVars.iphi[tower] ;
      int nbXtals = treeVars.nbOfXtals[tower] ;
      int ttf = getTtf(treeVars.rawTPData[tower]) ;


      if (verbose>9 && (tp>0 || maxOfTPEmul>0)) {
	std::cout<<"(phi,eta, Nbxtals)="<<std::dec<<iphi<<" "<<ieta<<" "<<nbXtals<<std::endl ;
	std::cout<<"Data Et, TTF: "<<tp<<" "<<ttf<<std::endl ;
	std::cout<<"Emulator: " ;
	for (int i=0 ; i<5 ; i++) std::cout<<emul[i]<<" " ;
	std::cout<<std::endl ;
      }


      // Fill TP spctrum
      TP->Fill(tp) ;
      TPEmul->Fill(emul[ref]) ;
      TPEmulMax->Fill(maxOfTPEmul) ;
      TPspectrumMap3D->Fill(iphi, ieta, tp) ;


      // Fill TP occupancy
      if (tp>occupancyCut) occupancyTP->Fill(iphi, ieta) ;
      if (emul[ref]>occupancyCut) occupancyTPEmul->Fill(iphi, ieta) ;


      // Fill TP-Emulator matching
      // comparison is meaningful when:
      if (tp>0 && nbXtals == 25) {
	bool match(false) ;
	for (int i=0 ; i<5 ; i++) {
	  if (tp == emul[i]) {
	    TPMatchEmul->Fill(i+1) ;
	    TPMatchEmul3D->Fill(iphi, ieta, i+1) ;
	    match = true ;
	  }
	}
	if (!match) {
	  TPMatchEmul->Fill(-1) ;
	  TPMatchEmul3D->Fill(iphi, ieta, -1) ;
	  if (verbose>5) {
	    std::cout<<"MISMATCH"<<std::endl ;
	    std::cout<<"(phi,eta, Nbxtals)="<<std::dec<<iphi<<" "<<ieta<<" "<<nbXtals<<std::endl ;
	    std::cout<<"Data Et, TTF: "<<tp<<" "<<ttf<<std::endl ;
	    std::cout<<"Emulator: " ;
	    for (int i=0 ; i<5 ; i++) std::cout<<emul[i]<<" " ;
	    std::cout<<std::endl ;
	  }