C++ TLorentzVector::Phi方法代码示例

void rochcor::momcor_mc( TLorentzVector& mu, float charge, float sysdev, int runopt){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  //float mptsys = sran.Gaus(0.0,sysdev);
  
  float dm = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/mmavg[mu_phibin][mu_etabin];
  float da = mcor_ma[mu_phibin][mu_etabin] + mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin];
  
  float cor = 1.0/(1.0 + dm + charge*da*ptmu);

  //for the momentum tuning - eta,phi,Q correction
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  float gscler = 0.0;
  float deltaer = 0.0;
  float sfer = 0.0;
  
  gscler = TMath::Sqrt( TMath::Power(mgscl_stat,2) + TMath::Power(mgscl_syst,2) );
  deltaer = TMath::Sqrt( TMath::Power(delta_stat,2) + TMath::Power(delta_syst,2) );
  sfer = TMath::Sqrt( TMath::Power(sf_stat,2) + TMath::Power(sf_syst,2) );
  
  float tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*eran.Gaus(1.0,(sf + sysdev*sfer)));
  
  px *= (tune); 
  py *= (tune);  
  pz *= (tune);  
  e  *= (tune);   
  
  float gscl = (genm_smr/mrecm);
  
  px *= (gscl + gscler_mc_dev*gscler);
  py *= (gscl + gscler_mc_dev*gscler);
  pz *= (gscl + gscler_mc_dev*gscler);
  e  *= (gscl + gscler_mc_dev*gscler);

  mu.SetPxPyPzE(px,py,pz,e);
  
}

double Histogrammer::minDrPhoB(int PhoInd, EventTree* tree){
	// find the closest b-jet
	TLorentzVector b;
	TLorentzVector bBar;
	int phoGen=-1;
	double mindr = 999.0;
	for( int mcI = 0; mcI < tree->nMC_; ++mcI){
		if( tree->mcIndex->at(mcI) == tree->phoGenIndex_->at(PhoInd) ) 
			phoGen=mcI;
		if( tree->mcPID->at(mcI) == 5) 
			b.SetPtEtaPhiM(tree->mcPt->at(mcI), tree->mcEta->at(mcI), tree->mcPhi->at(mcI), tree->mcMass->at(mcI));
		if( tree->mcPID->at(mcI) == -5) 
			bBar.SetPtEtaPhiM(tree->mcPt->at(mcI), tree->mcEta->at(mcI), tree->mcPhi->at(mcI), tree->mcMass->at(mcI));
	}
	if( phoGen > 0 && b.Pt() > 0.0001 && bBar.Pt() > 0.0001 ) {
		mindr = std::min(dR(tree->mcEta->at(phoGen), tree->mcPhi->at(phoGen), b.Eta(), b.Phi()),
						 dR(tree->mcEta->at(phoGen), tree->mcPhi->at(phoGen), bBar.Eta(), bBar.Phi()));
	}
	return mindr;
}

void rochcor2012jan22::momcor_mc( TLorentzVector& mu, float charge, int runopt, float& qter){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);
  
  if(mu_phibin>=0 && mu_etabin>=0){
    
    float Mf = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/(mpavg[mu_phibin][mu_etabin]+mmavg[mu_phibin][mu_etabin]);
    float Af = ((mcor_ma[mu_phibin][mu_etabin]+mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin]) - Mf*(mpavg[mu_phibin][mu_etabin]-mmavg[mu_phibin][mu_etabin]));     
    
    float cor = 1.0/(1.0 + 2.0*Mf + charge*Af*ptmu);
    
    //for the momentum tuning - eta,phi,Q correction
    px *= cor;
    py *= cor;
    pz *= cor;
    e  *= cor;
    
    float gscler = mgscl_stat;    
    float gscl = (genm_smr/mrecm);
    
    px *= (gscl + gscler_mc_dev*gscler);
    py *= (gscl + gscler_mc_dev*gscler);
    pz *= (gscl + gscler_mc_dev*gscler);
    e  *= (gscl + gscler_mc_dev*gscler);
    
    float momscl = sqrt(px*px + py*py)/ptmu;
    
    float tune = gsf[mu_etabin]*eran.Gaus(1.0,sf[mu_etabin]);
    
    px *= (tune); 
    py *= (tune);  
    pz *= (tune);  
    e  *= (tune);   
    
    qter *= sqrt(momscl*momscl + (1.0-tune)*(1.0-tune));
  }
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

void rochcor2012::momcor_data( TLorentzVector& mu, float charge, int runopt, float& qter){
  
  float ptmu = mu.Pt();

  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();
  
  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  float Mf = 0.0;
  float Af = 0.0;

  if(runopt==0){
    Mf = (dcor_bf[mu_phibin][mu_etabin]+mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfer[mu_phibin][mu_etabin])/(dpavg[mu_phibin][mu_etabin]+dmavg[mu_phibin][mu_etabin]);
    Af = ((dcor_ma[mu_phibin][mu_etabin]+mptsys_da_da[mu_phibin][mu_etabin]*dcor_maer[mu_phibin][mu_etabin]) - Mf*(dpavg[mu_phibin][mu_etabin]-dmavg[mu_phibin][mu_etabin]));     
  }else if(runopt==1){
    Mf = (dcor_bfD[mu_phibin][mu_etabin]+mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfDer[mu_phibin][mu_etabin])/(dpavgD[mu_phibin][mu_etabin]+dmavgD[mu_phibin][mu_etabin]);
    Af = ((dcor_maD[mu_phibin][mu_etabin]+mptsys_da_da[mu_phibin][mu_etabin]*dcor_maDer[mu_phibin][mu_etabin]) - Mf*(dpavgD[mu_phibin][mu_etabin]-dmavgD[mu_phibin][mu_etabin]));     
  }
  
  float cor = 1.0/(1.0 + 2.0*Mf + charge*Af*ptmu);
  
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  //after Z pt correction
  float gscler = dgscl_stat;
  float gscl = (genm_smr/drecm);
  
  px *= (gscl + gscler_da_dev*gscler);
  py *= (gscl + gscler_da_dev*gscler);
  pz *= (gscl + gscler_da_dev*gscler);
  e  *= (gscl + gscler_da_dev*gscler);
  
  float momscl = sqrt(px*px + py*py)/ptmu;
  qter *= momscl;
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

TLorentzVector doCalEnergy(double BeamEnergy,
                           TLorentzVector Particle1,
                           TLorentzVector Particle2,
                           double nucleusMass,
                           double Particle2Mass,
                           double Particle3Mass)
{
    double E_Particle1   = Particle1.E();
    double p_Particle1_x = Particle1.Px();
    double p_Particle1_y = Particle1.Py();
    double p_Particle1_z = Particle1.Pz();
    double p_Particle1   = sqrt(TMath::Power(p_Particle1_x,2.0) +
                                    TMath::Power(p_Particle1_y,2.0) +
                                    TMath::Power(p_Particle1_z,2.0));
    double phi   = Particle2.Phi();
    double theta = Particle2.Theta();
    double b     = 2.0 * ( p_Particle1_x * cos(phi) * sin(theta) +
                           p_Particle1_y * sin(phi) * sin(theta) +
                           p_Particle1_z * cos(theta) -
                           BeamEnergy * cos(theta)
                         );
    double c     = p_Particle1 * p_Particle1 + BeamEnergy * BeamEnergy - 2.0 * BeamEnergy * p_Particle1_z;
    double d     = BeamEnergy + nucleusMass - E_Particle1;
    double e     = TMath::Power(Particle3Mass,2.0) - TMath::Power(Particle2Mass,2.0) - d * d + c;
    double Delta = 16.0 * TMath::Power(d,2.0) * (TMath::Power(e,2.0) +
                                                 TMath::Power(b * Particle2Mass,2.0) -
                                                 TMath::Power(d * Particle2Mass * 2.0,2.0));
    
    TLorentzVector NewParticle(0.0,0.0,0.0,0.0);
    if(Delta>0.)
    {
       double sol2     = (2.0 * e * b + sqrt(Delta)) / (2.0 * (4.0 * TMath::Power(d,2.0) - TMath::Power(b,2.0)));
       double newpxcal = sol2 * cos(phi) * sin(theta);
       double newpycal = sol2 * sin(phi) * sin(theta);
       double newpzcal = sol2 * cos(theta);
       double energy   = sqrt(TMath::Power(sol2,2.0) + TMath::Power(Particle2Mass,2.0));
      
       TLorentzVector NewParticle2(newpxcal,newpycal,newpzcal,energy);
       NewParticle = NewParticle2;
    }
  
    return NewParticle;
}

void rochcor::momcor_data( TLorentzVector& mu, float charge, float sysdev, int runopt){
  
  float ptmu = mu.Pt();

  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();
  
  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);
  
  //float mptsys1 = sran.Gaus(0.0,sysdev);
  
  float dm = (dcor_bf[mu_phibin][mu_etabin] + mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfer[mu_phibin][mu_etabin])/dmavg[mu_phibin][mu_etabin];
  float da = dcor_ma[mu_phibin][mu_etabin] + mptsys_da_da[mu_phibin][mu_etabin]*dcor_maer[mu_phibin][mu_etabin];
  
  float cor = 1.0/(1.0 + dm + charge*da*ptmu);
  
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  //after Z pt correction
  float gscler = 0.0;
  
  gscler = TMath::Sqrt( TMath::Power(dgscl_stat,2) + TMath::Power(dgscl_syst,2) );
  
  float gscl = (genm_smr/drecm);
  
  px *= (gscl + gscler_da_dev*gscler);
  py *= (gscl + gscler_da_dev*gscler);
  pz *= (gscl + gscler_da_dev*gscler);
  e  *= (gscl + gscler_da_dev*gscler);
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

double Delta(double pt1, double eta1, double phi1,double pt2, double eta2, double phi2,double t1pfmetPhi){
  TLorentzVector fLorentzVec1; 
  TLorentzVector fLorentzVec2; 
  TLorentzVector fLorentzVecgg;
  
  // Fill TLorentzVector
  fLorentzVec1.SetPtEtaPhiM(pt1,eta1,phi1,0.);
  fLorentzVec2.SetPtEtaPhiM(pt2,eta2,phi2,0.);
  fLorentzVecgg = fLorentzVec1 + fLorentzVec2;
  double phigg= fLorentzVecgg.Phi();
  double delta = phigg-t1pfmetPhi;
  if(delta>3.14){
    delta-=6.28;
  }
  if(delta<-3.14){
    delta+=6.28;
  }
  delta =fabs(delta);
  return delta;
}

void smearEmEnergy( TLorentzVector& p ) {

    float smearEBlowEta     = 0.013898;
    float smearEBhighEta    = 0.0189895;
    float smearEElowEta     = 0.027686;
    float smearEEhighEta    = 0.031312;
    float theGaussMean      = 1.;
    float pt = p.Pt();
    float eta = p.Eta();
    float phi = p.Phi();
    float mass = p.M();
    float theSmear = 0.;
    if      (fabs(eta)<1.                   )  theSmear = smearEBlowEta;
    else if (fabs(eta)>=1.  && fabs(eta)<1.5)  theSmear = smearEBhighEta;
    else if (fabs(eta)>=1.5 && fabs(eta)<2. )  theSmear = smearEElowEta;
    else if (fabs(eta)>=2.  && fabs(eta)<2.5)  theSmear = smearEEhighEta;
    float fromGauss = myRandom_.Gaus(theGaussMean,theSmear);
    p.SetPtEtaPhiM( fromGauss*pt, eta, phi, mass ); // keep mass and direction same

}

std::vector<double> CalculateDeltaEtaAndPhi(bool possibleMultW, bool possibleMultAntiW, int w1, TLorentzVector TL_munu){

  //Don't take into account the events with multiple W's. Think of something

  double deltaphi, deltaeta;
  if(possibleMultW == false && possibleMultAntiW == false){
    deltaphi = TMath::Abs(TVector2::Phi_mpi_pi(TL_munu.Phi() - genPhi[w1]) );
     deltaeta = TMath::Abs(TL_munu.Eta() - genEta[w1] );
  } else {
    deltaphi = 99.;
    deltaeta = 99.;
  }

  std::vector<double> EtaAndPhi;

  EtaAndPhi.push_back(deltaeta);
  EtaAndPhi.push_back(deltaphi);

  return EtaAndPhi;

}

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



    int nentries = tree->GetEntries();

    for( int iEntry=0; iEntry<nentries; ++iEntry ) {

        if( iEntry % 50000 == 0 ) std::cout << "    Entry: " << iEntry << " / " << nentries << std::endl;

        myTree.GetEntry(iEntry);

        float weight = (myTree.isData) ? 1. : myTree.evt_scale1fb;
        // pu reweighting:
        if( !myTree.isData ) {
            //weight *= myTree.puWeight;
        }


        // first find leptons
        //if( myTree.ngenLep!=2 ) continue;
        //TLorentzVector genLep0, genLep1;
        //genLep0.SetPtEtaPhiM( myTree.genLep_pt[0], myTree.genLep_eta[0], myTree.genLep_phi[0], myTree.genLep_mass[0] );
        //genLep1.SetPtEtaPhiM( myTree.genLep_pt[1], myTree.genLep_eta[1], myTree.genLep_phi[1], myTree.genLep_mass[1] );

        std::vector<TLorentzVector> genLeptons;
        int genLeptType = 0;
        for( int iGen=0; iGen<myTree.ngenPart && genLeptons.size()<2; ++iGen ) {

            if( myTree.genPart_pt[iGen]<1. ) continue;
            if( myTree.genPart_status[iGen]!=1 ) continue;
            if( abs(myTree.genPart_pdgId[iGen])!=11 && abs(myTree.genPart_pdgId[iGen])!=13 ) continue;
            if( myTree.genPart_motherId[iGen]!=myTree.genPart_pdgId[iGen] ) continue;
            TLorentzVector tmpLep;
            tmpLep.SetPtEtaPhiM( myTree.genPart_pt[iGen], myTree.genPart_eta[iGen], myTree.genPart_phi[iGen], myTree.genPart_mass[iGen] );
            genLeptons.push_back(tmpLep);
            genLeptType = abs(myTree.genPart_pdgId[iGen]);

        }

        if( genLeptType!=11 && genLeptType!=13 ) continue;
        if( genLeptons.size()<2 ) continue;

        TLorentzVector genLep0, genLep1;
        genLep0 = genLeptons[0];
        genLep1 = genLeptons[1];

        float maxPt = TMath::Max( genLep0.Pt(), genLep1.Pt() );
        float minPt = TMath::Min( genLep0.Pt(), genLep1.Pt() );
        if( maxPt<25. ) continue;
        if( minPt<20. ) continue;

        if( fabs(genLep0.Eta()) > 2.4 ) continue;
        if( fabs(genLep1.Eta()) > 2.4 ) continue;

        TLorentzVector genZ = genLep0 + genLep1;
        if( genZ.M()<50. ) continue;
        //if( genZ.M()<50. || genZ.M()>130. ) continue;


        TLorentzVector genPhoton;
        bool foundGenPhoton = false;

        for( int iGen=0; iGen<myTree.ngenPart && !foundGenPhoton; ++iGen ) {

            if( myTree.genPart_pdgId[iGen]!=22 ) continue;
            if( myTree.genPart_status[iGen]!=1 ) continue;

void fill(int const kf, TLorentzVector* b, double const weight, TClonesArray& daughters)
{
   TLorentzVector kMom;
   TLorentzVector p1Mom;
   TLorentzVector p2Mom;

   int nTrk = daughters.GetEntriesFast();
   for (int iTrk = 0; iTrk < nTrk; ++iTrk)
   {
      TParticle* ptl0 = (TParticle*)daughters.At(iTrk);

      switch(abs(ptl0->GetPdgCode()))
      {
        case 321:
          ptl0->Momentum(kMom);
          break;
        case 211:
          if(!p1Mom.P()) ptl0->Momentum(p1Mom);
          else ptl0->Momentum(p2Mom);
          break;
        default:
          break;
      }
   }
   daughters.Clear();

   // smear and get total momentum
   TLorentzVector kRMom = smearMom(kMom,fKaonMomResolution);
   TLorentzVector p1RMom = smearMom(p1Mom,fPionMomResolution);
   TLorentzVector p2RMom = smearMom(p2Mom,fPionMomResolution);
   TLorentzVector rMom = kRMom + p1RMom + p2RMom;

   // save
   float arr[100];
   int iArr = 0;
   arr[iArr++] = kf;
   arr[iArr++] = weight;
   arr[iArr++] = b->M();
   arr[iArr++] = b->Perp();
   arr[iArr++] = b->PseudoRapidity();
   arr[iArr++] = b->Rapidity();
   arr[iArr++] = b->Phi();

   arr[iArr++] = rMom.M();
   arr[iArr++] = rMom.Perp();
   arr[iArr++] = rMom.PseudoRapidity();
   arr[iArr++] = rMom.Rapidity();
   arr[iArr++] = rMom.Phi();

   arr[iArr++] = kMom.M();
   arr[iArr++] = kMom.Perp();
   arr[iArr++] = kMom.PseudoRapidity();
   arr[iArr++] = kMom.Rapidity();
   arr[iArr++] = kMom.Phi();

   arr[iArr++] = kRMom.M();
   arr[iArr++] = kRMom.Perp();
   arr[iArr++] = kRMom.PseudoRapidity();
   arr[iArr++] = kRMom.Rapidity();
   arr[iArr++] = kRMom.Phi();

   arr[iArr++] = p1Mom.M();
   arr[iArr++] = p1Mom.Perp();
   arr[iArr++] = p1Mom.PseudoRapidity();
   arr[iArr++] = p1Mom.Rapidity();
   arr[iArr++] = p1Mom.Phi();

   arr[iArr++] = p1RMom.M();
   arr[iArr++] = p1RMom.Perp();
   arr[iArr++] = p1RMom.PseudoRapidity();
   arr[iArr++] = p1RMom.Rapidity();
   arr[iArr++] = p1RMom.Phi();

   arr[iArr++] = p2Mom.M();
   arr[iArr++] = p2Mom.Perp();
   arr[iArr++] = p2Mom.PseudoRapidity();
   arr[iArr++] = p2Mom.Rapidity();
   arr[iArr++] = p2Mom.Phi();

   arr[iArr++] = p2RMom.M();
   arr[iArr++] = p2RMom.Perp();
   arr[iArr++] = p2RMom.PseudoRapidity();
   arr[iArr++] = p2RMom.Rapidity();
   arr[iArr++] = p2RMom.Phi();
   
   nt->Fill(arr);
}

void skim(std::string var, TTree *fullTree, TTree *newTree, std::string cutstr, bool isWeighted , bool reweighted, bool addDiMuSel, bool addMuSel, bool addPhoSel){


  std::cout << fullTree->GetName() << " -> " <<  newTree->GetName() << " " << cutstr.c_str() << std::endl;
  fullTree->Draw(">>cutlist",cutstr.c_str());
  TEventList *keep_points = (TEventList*)gDirectory->Get("cutlist");
  int nEntries = fullTree->GetEntries();

  float x,x_orig;
  float mvametPhi, metRaw;
  float mvametPhiCor ;// corected phis for deltaPhi(met,jet) cut (// NOT to be used in calculation of MT!)
  float weight, weight_wpt, weight_in, weight_in_pu; 
  float catvar1,catvar2,catvar3,catvar4,catvar5,catvar6;   // These are empty, do as we want with them 
  float tau1,tau2;
  float dimuM, mt;	// used for the W/Z selection
  int passVBF;
  unsigned int njets;
  unsigned int nphotons;
  float phopt, phoeta;
  // Why do we use these?
  float metBasicPhi, metBasic;

 
  // Used for reweighting if needed
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *genZ  = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
  TBranch *brgenZ ;

  // W/Z selection for leptions
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *lepV   = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *lep2V  = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();

  TBranch *brlep   = fullTree->GetBranch("lep1");
  TBranch *brlep_2 = fullTree->GetBranch("lep2");
  brlep->SetAddress(&lepV);
  brlep_2->SetAddress(&lep2V);

  // Check VBF selection 
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *jetV   = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *jet2V  = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();

  // Also need these for second jet veto 
  TBranch *brjet   = fullTree->GetBranch("jet1");
  TBranch *brjet_2 = fullTree->GetBranch("jet2");
  brjet->SetAddress(&jetV);
  brjet_2->SetAddress(&jet2V);


  // if addPhoSel add the two photons
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *pho1V  = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > *pho2V  = new ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> >();
  TBranch *brpho_1 ;
  TBranch *brpho_2 ;
  if (addPhoSel){
    brpho_1 = fullTree->GetBranch("pho1");
    brpho_2 = fullTree->GetBranch("pho2");    
    brpho_1->SetAddress(&pho1V);
    brpho_2->SetAddress(&pho2V);
  }


  // mu/dimu ?
  std::string origVar = var;
  if (addMuSel)   var+="CorW";
  if (addDiMuSel) var+="CorZ";

  // Add Output branches to our very skimmed down tree.
  //TTree *newTree = new TTree(newname.c_str(),newname.c_str());

  newTree->Branch("mvamet",&x,"mvamet/Float_t");
  newTree->Branch("mvamet_orig",&x_orig,"mvamet_orig/Float_t");
  newTree->Branch("weight",&weight,"weight/Float_t");
  newTree->Branch("weight_wpt",&weight_wpt,"weight_wpt/Float_t");
  newTree->Branch("metRaw",&metRaw,"metRaw/Float_t"); // useful to keep this
  newTree->Branch("jet1mprune",&catvar1,"jet1mprume/Float_t");
  newTree->Branch("jet1mtrim",&catvar2,"jet1mtrim/Float_t");
  newTree->Branch("jet1tau2o1",&catvar3,"jet1tau12o1/Float_t");
  newTree->Branch("jpt",&catvar4,"jpt/Float_t");
  newTree->Branch("jet1eta",&catvar6,"jet1eta/Float_t");
  newTree->Branch("jet1QGtag",&catvar5,"jet1QGtag/Float_t");
  newTree->Branch("passVBF",&passVBF,"passVBF/Int_t");
  newTree->Branch("njets",&njets,"njets/Int_t");

  if (addDiMuSel) {	
    newTree->Branch("dimu_m",&dimuM,"dimu_m/Float_t");
  }
  if (addMuSel) {	
    newTree->Branch("mt",&mt,"mt/Float_t");
  }
  if (addPhoSel) {	
    newTree->Branch("phopt",&phopt,"phopt/Float_t");
    newTree->Branch("phoeta",&phoeta,"phoeta/Float_t");
  }

  fullTree->SetBranchAddress(var.c_str(),&x);
  if (addDiMuSel || addMuSel) fullTree->SetBranchAddress(origVar.c_str(),&x_orig);
  fullTree->SetBranchAddress(Form("%sPhi",var.c_str()),&mvametPhiCor);
  if (addDiMuSel || addMuSel) fullTree->SetBranchAddress(Form("%sPhi",origVar.c_str()),&mvametPhi);
  if (isWeighted) fullTree->SetBranchAddress("weight",&weight_in);
  fullTree->SetBranchAddress("metRaw",&metRaw);
  fullTree->SetBranchAddress("njets",&njets);
//.........这里部分代码省略.........

void tree1r()
{
 bool muon=0;
//JetCorrectionUncertainty *jecUnc;
// jecUnc= new JetCorrectionUncertainty("Fall12_V7_DATA_Uncertainty_AK5PFchs.txt");

  TChain myTree("analyzeBasicPat/MuonTree");
  myTree.Add("/afs/cern.ch/user/b/bbilin/eos/cms/store/group/phys_smp/ZPlusJets/8TeV/ntuples/2111/skimmed/EE_hadd_2012ABCD.root");

  TH1::AddDirectory(true);



  Double_t        PU_npT;
  Double_t        PU_npIT;
  double           EvtInfo_NumVtx;
  double          MyWeight;
   double nup;
  Int_t          event;
  double          realdata;
  int          run;
//  int          lumi;
  double HLT_Elec17_Elec8;

  std::vector<double> *Dr01LepPt=0;
  std::vector<double> *Dr01LepEta=0;
  std::vector<double> *Dr01LepPhi=0;
  std::vector<double> *Dr01LepE=0;
  std::vector<double> *Dr01LepM=0;
  std::vector<double> *Dr01LepId=0;
  std::vector<double> *Dr01LepStatus=0;

  std::vector<double>  *Bare01LepPt=0;
  std::vector<double> *Bare01LepEta=0;
  std::vector<double> *Bare01LepPhi=0;
  std::vector<double> *Bare01LepE=0;
  std::vector<double> *Bare01LepM=0;
  std::vector<double> *Bare01LepId=0;
  std::vector<double> *Bare01LepStatus=0;

  std::vector<double> *St01PhotonPt=0;
  std::vector<double> *St01PhotonEta=0;
  std::vector<double> *St01PhotonPhi=0;
  std::vector<double> *St01PhotonE=0;
  std::vector<double> *St01PhotonM=0;
  std::vector<double> *St01PhotonId=0;
  std::vector<double> *St01PhotonMomId=0;
  std::vector<double> *St01PhotonNumberMom=0;
  std::vector<double> *St01PhotonStatus=0;



  std::vector<double> *St03Pt=0;
  std::vector<double> *St03Eta=0;
  std::vector<double> *St03Phi=0;
  std::vector<double> *St03E=0;
  std::vector<double> *St03M=0;
  std::vector<double> *St03Id=0;
  std::vector<double> *St03Status=0;

  std::vector<double> *GjPt=0;
  std::vector<double> *Gjeta=0;
  std::vector<double> *Gjphi=0;
  std::vector<double> *GjE=0;
  std::vector<double> *GjPx=0;
  std::vector<double> *GjPy=0;
  std::vector<double> *GjPz=0;

  vector<double>  *patMuonEn_ =0; 
  vector<double>  *patMuonCharge_ =0;
  vector<double>  *patMuonPt_ =0;
  vector<double>  *patMuonEta_=0;
  vector<double>  *patMuonPhi_ =0;
  vector<double>  *patMuonCombId_ =0;
  vector<double>  *patMuonTrig_ =0;
  vector<double>  *patMuonDetIsoRho_ =0;
  vector<double>  *patMuonPfIsoDbeta_ =0;

  std::vector<double> *patElecTrig_ =0;
  std::vector<double> *patElecCharge_ =0;
  std::vector<double> *patElecEnergy_ =0;
  std::vector<double> *patElecEta_ =0;
  std::vector<double> *patElecScEta_ =0;
  std::vector<double> *patElecPhi_ =0;
//  std::vector<double> *patElecEcalEnergy_ =0;
  std::vector<double> *patElecPt_ =0;
  std::vector<double> *patElecPfIso_ =0;
  std::vector<double> *patElecPfIsodb_ =0;
  std::vector<double> *patElecPfIsoRho_ =0;
  std::vector<double>   *patElecMediumIDOff_ =0;
//  std::vector<double> *patElecMVATrigId_ =0;
//  std::vector<double> *patElecMVANonTrigId_ =0;

  vector<double>  *patJetPfAk05En_ =0;
  vector<double>  *patJetPfAk05Pt_ =0;
  vector<double>  *patJetPfAk05Eta_ =0;
  vector<double>  *patJetPfAk05Phi_ =0;
//  vector<double>  *patJetPfAk05JesUncert_ =0;
  vector<double>  *patJetPfAk05LooseId_ =0;
  vector<double>  *patJetPfAk05Et_ =0;
//.........这里部分代码省略.........

Bool_t monojet::Process(Long64_t entry)
{
    GetEntry(entry);

    if( entry % 100000 == 0 ) cout << "Processing event number: " << entry << endl;
    //cout << "Processing event number: " << entry << endl;

    // To make the processing fast, apply a very looose selection
    if (((TLorentzVector*)((*metP4)[0]))->Pt() < 40. or jetP4->GetEntries() < 1) return kTRUE;

    //this is the type tree
    tm->SetValue("run",runNum);
    tm->SetValue("event",eventNum);
    tm->SetValue("lumi",lumiNum);

    float dR = 0.;
    TClonesArray *tightLep;
    TClonesArray *cleanJet;
    TClonesArray *cleanTau;

    tightLep = new TClonesArray("TLorentzVector",20);
    cleanJet = new TClonesArray("TLorentzVector",20);
    cleanTau = new TClonesArray("TLorentzVector",20);

    std::vector<bool>  jetMonojetId_clean;
    jetMonojetId_clean.clear();
    std::vector<bool>  jetMonojetIdLoose_clean;
    jetMonojetIdLoose_clean.clear();
    //std::vector<float> jetPuId_clean;
    //jetPuId_clean.clear();

    std::vector<float>  tauId_clean;
    tauId_clean.clear();
    std::vector<float>  tauIso_clean;
    tauIso_clean.clear();

    int n_tightlep = 0;

    // ********* Leptons ********** //
    for(int lepton = 0; lepton < lepP4->GetEntries(); lepton++) {
        TLorentzVector* Lepton = (TLorentzVector*) lepP4->At(lepton);
        // check if this is a tight lep, and check the overlap

        //iso_1 = divide(input_tree.lepIso[0],input_tree.lepP4[0].Pt())
        //if (input_tree.lepTightId[0]==0 or iso_1 > 0.12): continue

        if (Lepton->Pt() > 20. && (*lepTightId)[lepton] > 1) {
            n_tightlep +=1;
            new ( (*tightLep)[tightLep->GetEntriesFast()]) TLorentzVector(Lepton->Px(), Lepton->Py(), Lepton->Pz(), Lepton->Energy());

            //check overlap with jets
            for(int j = 0; j < jetP4->GetEntries(); j++) {
                TLorentzVector* Jet = (TLorentzVector*) jetP4->At(j);
                dR = deltaR(Lepton,Jet);
                if (dR > dR_cut) {
                    new ( (*cleanJet)[cleanJet->GetEntriesFast()]) TLorentzVector(Jet->Px(), Jet->Py(), Jet->Pz(), Jet->Energy());
                    jetMonojetId_clean.push_back((*jetMonojetId)[j]);
                    jetMonojetIdLoose_clean.push_back((*jetMonojetIdLoose)[j]);
                    //jetPuId_clean.push_back((*jetPuId)[j]);
                }
            }
            //check overlap with taus
            for(int tau = 0; tau < tauP4->GetEntries(); tau++) {
                TLorentzVector* Tau = (TLorentzVector*) tauP4->At(tau);
                dR = deltaR(Lepton,Tau);
                if (dR > dR_cut) new ( (*cleanTau)[cleanTau->GetEntriesFast()]) TLorentzVector(Tau->Px(), Tau->Py(), Tau->Pz(), Tau->Energy());
                tauId_clean.push_back((*tauId)[tau]);
                tauIso_clean.push_back((*tauIso)[tau]);
            } // tau overlap
        } // tight lepton selection
    }//lepton loop

    tm->SetValue("n_tightlep",n_tightlep);

    TLorentzVector fakeMET;

    // Z Selection
    TLorentzVector Z;
    if(lepP4->GetEntries() == 2 && n_tightlep > 0) {
        if (((*lepPdgId)[0]+(*lepPdgId)[1])==0 ) {
            Z = *((TLorentzVector*)((*lepP4)[0])) + *((TLorentzVector*)((*lepP4)[1]));
            fakeMET = *((TLorentzVector*)((*metP4)[0])) + Z;
        }
    }

    float MT = 0.0;
    //// W Selection
    if(lepP4->GetEntries() == 1 && n_tightlep == 1) {
        fakeMET = *((TLorentzVector*)((*metP4)[0])) + *((TLorentzVector*)((*lepP4)[0])) ;
        MT = transverseMass( ((TLorentzVector*)((*lepP4)[0]))->Pt(), ((TLorentzVector*)((*lepP4)[0]))->Phi(), ((TLorentzVector*)((*metP4)[0]))->Pt(), ((TLorentzVector*)((*metP4)[0]))->Phi());
    }

    tm->SetValue("mt",MT);

    // ********* Jets ********** //
    for(int jet = 0; jet < jetP4->GetEntries(); jet++) {
        TLorentzVector* Jet = (TLorentzVector*) jetP4->At(jet);
        //cout << (*jetMonojetId)[0] <<endl;
        //cout << Jet->Pt()<<endl;
    }
//.........这里部分代码省略.........

void xAna_ele_subleading_ele_pt_optimize(TString inputFile, TCanvas *c1 , TCanvas *c2 , TCanvas *c3 ,TCanvas *c4 ,TCanvas *c5 , TCanvas *c6 ,
int mass_point ,double eff,double eff_err, TString dir_name, int dir_flag ,int signal_background_flag, int background_file_index ){


  // define histograms

  TString title2 = Form("ele pT for Zprime mass = %d",mass_point);
  TString title3 = Form("lepton pairs' pT for Zprime mass = %d",mass_point);
  TString title4 = Form("SD mass for Zprime mass = %d",mass_point);
  TString title5 = Form("Z mass for Zprime mass = %d",mass_point);
  TString title6 = Form("Zprime mass for Zprime mass = %d",mass_point);
  TString title7 = Form("|dEta| of ZH for Zprime mass = %d",mass_point);
  TString title8 = Form("|dPhi| of ZH for Zprime mass = %d",mass_point);
  TString title9 = Form("dR of ZH for Zprime mass = %d",mass_point);

  TH1D* h_ele_pT = new TH1D("h_ele_pT", title2 ,3000 , 0,3000 );
  TH1D* h_lepton_pair_pT = new TH1D("h_lepton_pair_pT", title3 ,400 , 0,4000 );
  TH1D* h_SD =new TH1D("h_SD",title4 ,100,0,200);
  TH1D* h_Z_mass = new TH1D("h_Z_mass",title5 ,250,0,500);
  TH1D* h_Zprime_mass = new TH1D("h_Zprime_mass",title6 ,1000,0,5000);
  TH1D* h_abs_dEta_ZH = new TH1D("h_abs_dEta_ZH",title7 ,80,-1,3);
  TH1D* h_abs_dPhi_ZH = new TH1D("h_abs_dPhi_ZH",title8 ,80,0,4);
  TH1D* h_abs_dR_ZH   = new TH1D("h_abs_dR_ZH"  ,title9 ,80,0,4);

  //get TTree from file ...
//  TreeReader data(inputFile.data());
  TreeReader data(inputFile.Data());

  Long64_t nTotal=0;
  Long64_t nPass[20]={0};
  ofstream fout;
  fout.open("ele_Eiko.txt");

  //Event loop
  for(Long64_t jEntry=0; jEntry<data.GetEntriesFast() ;jEntry++){

    if (jEntry % 50000 == 0)
      fprintf(stderr, "Processing event %lli of %lli\n", jEntry + 1, data.GetEntriesFast());

    data.GetEntry(jEntry);
    nTotal ++;

    // 0. check the generator-level information and make sure there is a Z->e+e-

/*
    Int_t nGenPar        = data.GetInt("nGenPar");
    Int_t* genParId      = data.GetPtrInt("genParId");
    Int_t* genParSt      = data.GetPtrInt("genParSt");
    Int_t* genMomParId   = data.GetPtrInt("genMomParId");
    Int_t* genDa1      = data.GetPtrInt("genDa1");
    Int_t* genDa2      = data.GetPtrInt("genDa2");
    Int_t* genMo1      = data.GetPtrInt("genMo1");
    Int_t* genMo2      = data.GetPtrInt("genMo2");

    bool hasLepton=false;

    for(int ig=0; ig < nGenPar; ig++){

      if(genParId[ig]!=23)continue;


      int da1=genDa1[ig];
      int da2=genDa2[ig];


      if(da1<0 || da2<0)continue;
      int da1pdg = genParId[da1];
      int da2pdg = genParId[da2];


      if(abs(da1pdg)==11)
     	hasLepton=true;

      if(hasLepton)break;

    }
*/

/*
    // 1. make sure there is a h-> bb, Yu-Hsiang change it
    bool hasHadron=false;

    for(int ig=0; ig < nGenPar; ig++){

      if(genParId[ig]!=25)continue;
      int da1=genDa1[ig];
      int da2=genDa2[ig];

      if(da1<0 || da2<0)continue;
      int da1pdg = genParId[da1];
      int da2pdg = genParId[da2];

      if(abs(da1pdg)==5)
     	hasHadron=true;

      if(hasHadron)break;

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