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

TLorentzVector makePtEtaPhiE(double ptCor, double pt, double eta, double phi, double e, 
                             bool rescaleEnergy=true)
{
    TLorentzVector j;
    j.SetPtEtaPhiE(ptCor, eta, phi, (rescaleEnergy ? (e * ptCor / pt) : e));
    return j;
}

void SingleTopTree::FillJets(){
  Jets.clear();
  bJets.clear();

  TLorentzVector muon;
  int tightMuIndex;
  bool mu = TightIso06Muons.size() > 0;
  
  if( mu )
    tightMuIndex = TightIso06Muons[0] ;
  if(!mu){
    mu = Tight12IsoMuons.size() > 0 ;
    if( mu )
      tightMuIndex = Tight12IsoMuons[0] ;
  }
  muon.SetPtEtaPhiE( muons_Pt[tightMuIndex] , muons_Eta[tightMuIndex] , muons_Phi[tightMuIndex] , muons_E[tightMuIndex] );

  for( int jid = 0 ; jid < jetsAK4_size ; jid++ ){

    //if( Event_EventNumber == 672782113 ) cout << (jetsAK4_CorrPt[jid] > 40) << " " <<  (fabs( jetsAK4_Eta[jid] ) < 4.7) << " " << jetsAK4_PassesID[jid];

    if( jetsAK4_CorrPt[jid] > 40 &&
	fabs( jetsAK4_Eta[jid] ) < 4.7 &&
	jetsAK4_PassesID[jid]
	){

      double DR = 1.0 ; 
      if( mu ){
	TLorentzVector jet;
	jet.SetPtEtaPhiE( jetsAK4_CorrPt[jid] , jetsAK4_Eta[jid] , jetsAK4_Phi[jid] , jetsAK4_CorrE[jid] );
	DR = muon.DeltaR( jet );
      }

      // if( Event_EventNumber == 672782113 )
      // 	cout << DR << endl;

      if ( DR > 0.3 ){
	Jets.push_back( jid );
	if( jetsAK4_IsCSVT[jid]  && fabs( jetsAK4_Eta[jid] ) <= 2.4 ){
	  bJets.push_back(jid);
	}
      }
    }
      //if( Event_EventNumber == 672782113 )
      //cout << endl;
  }
}

//Getting FourVectors
TLorentzVector D3PDSelector::GetFourVecByType(std::string ObjectType, unsigned int index)
{
	TLorentzVector par;
	float E = GetValByType<float>(ObjectType, "_E", index);
	float eta = GetValByType<float>(ObjectType, "_eta", index);
	float phi = GetValByType<float>(ObjectType, "_phi", index);
	float pt = GetValByType<float>(ObjectType, "_pt", index);
	par.SetPtEtaPhiE(pt,eta,phi,E);
	
	return par;
}

float BDTOhio_v2::Evaluate(std::string categoryLabel, const std::vector<pat::Muon>& selectedMuons, const std::vector<pat::Electron>& selectedElectrons, const std::vector<pat::Jet>& selectedJets, const std::vector<pat::Jet>& selectedJetsLoose, const pat::MET& pfMET){

  if(selectedMuons.size()+selectedElectrons.size()!=1){
    //    cerr << "BDTOhio_v2: not a SL event" << endl;
    return -2;
  }
  // ==================================================
  // construct object vectors etc
  TLorentzVector lepton_vec;
  TLorentzVector met_vec;
  vector<TLorentzVector> jet_vecs;
  vector<double> jetCSV;
  vector< vector<double> > jets_vvdouble;
  vector<TLorentzVector> jet_loose_vecs;
  vector<TLorentzVector> tagged_jet_vecs;
  vector<double> jetCSV_loose;
  vector<double> sortedCSV;
  if(selectedMuons.size()>0) lepton_vec.SetPtEtaPhiE(selectedMuons[0].pt(),selectedMuons[0].eta(),selectedMuons[0].phi(),selectedMuons[0].energy());
  if(selectedElectrons.size()>0) lepton_vec.SetPtEtaPhiE(selectedElectrons[0].pt(),selectedElectrons[0].eta(),selectedElectrons[0].phi(),selectedElectrons[0].energy());
  met_vec.SetPtEtaPhiE(pfMET.pt(),0,pfMET.phi(),pfMET.pt());
  for(auto jet=selectedJets.begin();jet!=selectedJets.end(); jet++){
    TLorentzVector jetvec;
    jetvec.SetPtEtaPhiE(jet->pt(),jet->eta(),jet->phi(),jet->energy());
    jet_vecs.push_back(jetvec);
    if(MiniAODHelper::GetJetCSV(*jet)>btagMcut){
      tagged_jet_vecs.push_back(jetvec);
    }
    vector<double> pxpypzE;
    pxpypzE.push_back(jet->px());
    pxpypzE.push_back(jet->py());
    pxpypzE.push_back(jet->pz());
    pxpypzE.push_back(jet->energy());
    jets_vvdouble.push_back(pxpypzE);
    jetCSV.push_back(MiniAODHelper::GetJetCSV(*jet));
  }
  sortedCSV=jetCSV;
  std::sort(sortedCSV.begin(),sortedCSV.end(),std::greater<float>());
  for(auto jet=selectedJetsLoose.begin();jet!=selectedJetsLoose.end(); jet++){
    TLorentzVector jetvec;
    jetvec.SetPtEtaPhiE(jet->pt(),jet->eta(),jet->phi(),jet->energy());
    jet_loose_vecs.push_back(jetvec);
    jetCSV_loose.push_back(MiniAODHelper::GetJetCSV(*jet));
  }

  // TODO loose jet and csv defintion

  // ==================================================
  // calculate variables
  // aplanarity and sphericity
  float aplanarity,sphericity;
  bdtvar.getSp(lepton_vec,met_vec,jet_vecs,aplanarity,sphericity);

  // Fox Wolfram
  float h0,h1,h2,h3,h4;
  bdtvar.getFox(jet_vecs,h0,h1,h2,h3,h4);

  // best higgs mass 1
  double minChi,dRbb;
  TLorentzVector bjet1,bjet2;
  float bestHiggsMass = bdtvar.getBestHiggsMass(lepton_vec,met_vec,jet_vecs,jetCSV,minChi,dRbb,bjet1,bjet2, jet_loose_vecs,jetCSV_loose);
  
  // study top bb system
  TLorentzVector dummy_metv;
  double minChiStudy, chi2lepW, chi2leptop, chi2hadW, chi2hadtop, mass_lepW, mass_leptop, mass_hadW, mass_hadtop, dRbbStudy, testquant1, testquant2, testquant3, testquant4, testquant5, testquant6, testquant7; 
  TLorentzVector b1,b2;
  bdtvar.study_tops_bb_syst (pfMET.pt(), pfMET.phi(), dummy_metv, lepton_vec, jets_vvdouble, jetCSV, minChiStudy, chi2lepW, chi2leptop, chi2hadW, chi2hadtop, mass_lepW, mass_leptop, mass_hadW, mass_hadtop, dRbbStudy, testquant1, testquant2, testquant3, testquant4, testquant5, testquant6, testquant7, b1, b2);
  float dEta_fn=testquant6;
  // ptE ratio
  float pt_E_ratio = bdtvar.pt_E_ratio_jets(jets_vvdouble);
  
  // etamax
  float jet_jet_etamax = bdtvar.get_jet_jet_etamax (jets_vvdouble);
  float jet_tag_etamax = bdtvar.get_jet_tag_etamax (jets_vvdouble,jetCSV);
  float tag_tag_etamax = bdtvar.get_tag_tag_etamax (jets_vvdouble,jetCSV);

  // jet variables
  float sum_pt_jets=0;
  float dr_between_lep_and_closest_jet=99;
  float mht_px=0;
  float mht_py=0;
  TLorentzVector p4_of_everything=lepton_vec;
  p4_of_everything+=met_vec;
  for(auto jetvec = jet_vecs.begin() ; jetvec != jet_vecs.end(); ++jetvec){
    dr_between_lep_and_closest_jet=fmin(dr_between_lep_and_closest_jet,lepton_vec.DeltaR(*jetvec));
    sum_pt_jets += jetvec->Pt();
    mht_px += jetvec->Px();
    mht_py += jetvec->Py();
    p4_of_everything += *jetvec;
  }
  mht_px+=lepton_vec.Px();
  mht_py+=lepton_vec.Py();
  float mass_of_everything=p4_of_everything.M();
  float sum_pt_wo_met=sum_pt_jets+lepton_vec.Pt();
  float sum_pt_with_met=pfMET.pt()+sum_pt_wo_met;
  float MHT=sqrt( mht_px*mht_px + mht_py*mht_py );

  float Mlb=0;   // mass of lepton and closest bt-tagged jet
  float minDr_for_Mlb=999.;
  for(auto tagged_jet=tagged_jet_vecs.begin();tagged_jet!=tagged_jet_vecs.end();tagged_jet++){
    float drLep=lepton_vec.DeltaR(*tagged_jet);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
				totaLlCounter++;
			}
			if ( (typ[inst] == 4 && jmas[inst] > 10. )|| (typ[inst] == 6 && pT[inst] > 10. )){
				lCounter = 0; //dont count the event
			}

		
		}//end instance loop (particles in an event

// 		cout<< "leptons in the event are "<< lCounter<<endl;
// 		if (lCounter == 4) { 

		fillFlag = false;

		// If else if loops reconstructing the particles according to the type 4e,4mu, 2e2mu
		
		if (el.size() == 1 && mu.size() == 1 && antiel.size() == 1 && antimu.size() == 1){ //2e2m
			goodEventCounter++;


			lep1.SetPtEtaPhi( pT[el[0]], eta[el[0]]	, phi[el[0]]);			//set up of lepton four-vectors
			lep2.SetPtEtaPhi( pT[antiel[0]], eta[antiel[0]]	, phi[antiel[0]]);
			lep3.SetPtEtaPhi( pT[mu[0]], eta[mu[0]]	, phi[mu[0]]);
			lep4.SetPtEtaPhi( pT[antimu[0]], eta[antimu[0]]	, phi[antimu[0]]);

			Za = lep1 + lep2;
			Zb = lep3 + lep4;

			mZ1 = sqrt(pow(lep1.Mag()+lep2.Mag(),2)-Za.Mag2());	// reconstruct z masses 
			mZ2 = sqrt(pow(lep3.Mag()+lep4.Mag(),2)-Zb.Mag2());


			//select leading Z
			if(mZ1 > mZ2) { Z1.SetVectM( Za, mZ1); Z2.SetVectM(Zb,mZ2); lep_min1.SetPtEtaPhiE(lep1.Pt(),lep1.Eta(), lep1.Phi(),lep1.Mag()); lep_plus1.SetPtEtaPhiE(lep2.Pt(),lep2.Eta(), lep2.Phi(),lep2.Mag()); lep_min2.SetPtEtaPhiE(lep3.Pt(),lep3.Eta(), lep3.Phi(),lep3.Mag()); lep_plus2.SetPtEtaPhiE(lep4.Pt(),lep4.Eta(), lep4.Phi(),lep4.Mag());eexx++;}	//to set the highest mass the z
			else { Z2.SetVectM( Za, mZ1); Z1.SetVectM(Zb,mZ2); lep_min2.SetPtEtaPhiE(lep1.Pt(),lep1.Eta(), lep1.Phi(),lep1.Mag()); lep_plus2.SetPtEtaPhiE(lep2.Pt(),lep2.Eta(), lep2.Phi(),lep2.Mag()); lep_min1.SetPtEtaPhiE(lep3.Pt(),lep3.Eta(), lep3.Phi(),lep3.Mag()); lep_plus1.SetPtEtaPhiE(lep4.Pt(),lep4.Eta(), lep4.Phi(),lep4.Mag());xxee++;}


			

		fillFlag = true;
		}



		else if (el.size() == 2 && mu.size() == 0  && antiel.size() == 2 && antimu.size() == 0){ //4e
			goodEventCounter++;

			lep1.SetPtEtaPhi( pT[el[0]], eta[el[0]]	, phi[el[0]]);	
			lep2.SetPtEtaPhi( pT[antiel[0]], eta[antiel[0]]	, phi[antiel[0]]);
			lep3.SetPtEtaPhi( pT[el[1]], eta[el[1]]	, phi[el[1]]);	
			lep4.SetPtEtaPhi( pT[antiel[1]], eta[antiel[1]]	, phi[antiel[1]]);

			Za = lep1 + lep2;
			Zb = lep3 + lep4;
			Zc = lep1 + lep4;
			Zd = lep3 + lep2;

			double mZa = sqrt(pow(lep1.Mag()+lep2.Mag(),2)-Za.Mag2());
			double mZb = sqrt(pow(lep3.Mag()+lep4.Mag(),2)-Zb.Mag2());
			double mZc = sqrt(pow(lep1.Mag()+lep4.Mag(),2)-Zc.Mag2());
			double mZd = sqrt(pow(lep2.Mag()+lep3.Mag(),2)-Zd.Mag2());

			double s1a;
			double s1b;
			double s2a;
			double s2b;

//.........这里部分代码省略.........
                    events.push_back(event);
                    it->second.insert( std::pair<int, std::vector<int> > (LS, events) );

                } else { //LS exists, look for event

                    std::vector<int>::iterator ev;
                    for( ev=it_LS->second.begin(); ev!=it_LS->second.end(); ++ev )
                        if( *ev==event ) break;


                    if( ev==it_LS->second.end() ) {

                        it_LS->second.push_back(event);

                    } else {

                        std::cout << "DISCARDING DUPLICATE EVENT!! Run: " << run << " LS: " << LS << " event: " << event << std::endl;

                        continue;

                    }
                }
            }

        } //if is mc




        h1_nvertex_PUW->Fill( nvertex, eventWeight);


        TLorentzVector lept1, lept2;
        lept1.SetPtEtaPhiE( ptLeptZ1, etaLeptZ1, phiLeptZ1, eLeptZ1 );
        lept2.SetPtEtaPhiE( ptLeptZ2, etaLeptZ2, phiLeptZ2, eLeptZ2 );

        TLorentzVector Z = lept1 + lept2;

        ptZ = Z.Pt();
        mZ = Z.M();
        etaZ = Z.Eta();

        h1_mZ->Fill( Z.M(), eventWeight );

        if( Z.M()<70. || Z.M()>110. ) continue;

        if( nJets==0 ) continue;
        if( ptJet[0] < 20. ) continue;


        // jet id:
        TLorentzVector jet;
        jet.SetPtEtaPhiE( ptJet[0], etaJet[0], phiJet[0], eJet[0] );

        if( fabs(jet.Eta())<2.4 && nChargedJet[0]==0 ) continue;
        if( (nNeutralJet[0]+nChargedJet[0])==1 ) continue;
        if( (ePhotonsJet[0]+eHFEMJet[0])/jet.E()>0.99 ) continue;
        if( (eNeutralHadronsJet[0])/jet.E()>0.99 ) continue;


        pdgIdPartJet_t = 0;

        if( isMC ) {

            // check if matched to parton/genjet
            TLorentzVector part;

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


    //Event filters at RECO level
    fillVariableWithValue( "PassBeamScraping", !isBeamScraping ) ;
    fillVariableWithValue( "PassPrimaryVertex", isPrimaryVertex ) ;
    fillVariableWithValue( "PassHBHENoiseFilter", passHBHENoiseFilter ) ;
    fillVariableWithValue( "PassBeamHaloFilterLoose", passBeamHaloFilterLoose ) ;
    fillVariableWithValue( "PassBeamHaloFilterTight", passBeamHaloFilterTight ) ;
    fillVariableWithValue( "PassTrackingFailure", !isTrackingFailure ) ;
    fillVariableWithValue( "PassCaloBoundaryDRFilter", passCaloBoundaryDRFilter ) ;
    fillVariableWithValue( "PassEcalMaskedCellDRFilter", passEcalMaskedCellDRFilter ) ;
    
    // Evaluate cuts (but do not apply them)
    evaluateCuts();


    //Basic Event Selection
    if( passedCut("PassJSON") 
	&& passedCut("run")
	&& passedCut("PassBPTX0") 
	&& passedCut("PassBeamScraping") 
	&& passedCut("PassPrimaryVertex")
	&& passedCut("PassHBHENoiseFilter")
	&& passedCut("PassBeamHaloFilterTight") 
	&& passedCut("PassHLT")
	)
      {      
	//Loop over probes
	for(int iprobe=0; iprobe<HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTSC17HEDoubleFilterPt->size(); iprobe++)
	  {

	    TLorentzVector probe;
	    probe.SetPtEtaPhiE(HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTSC17HEDoubleFilterPt->at(iprobe),
			       HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTSC17HEDoubleFilterEta->at(iprobe),
			       HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTSC17HEDoubleFilterPhi->at(iprobe),
			       HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTSC17HEDoubleFilterEnergy->at(iprobe)
			       );

	    int isProbeBarrel = 0;
	    int isProbeEndcap = 0;
	    
	    if( fabs( probe.Eta() ) < eleEta_bar )       isProbeBarrel = 1;
	    if( fabs( probe.Eta() ) > eleEta_end_min &&
		fabs( probe.Eta() ) < eleEta_end_max )   isProbeEndcap = 1;


	    CreateAndFillUserTH1D("Pt_Probe", 200, 0, 200, probe.Pt() );

	    //Loop over tags
 	    for(int itag=0; itag<HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTEle17TrackIsolFilterPt->size(); itag++)
	      {

		TLorentzVector tag;
		tag.SetPtEtaPhiE(HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTEle17TrackIsolFilterPt->at(itag),
				 HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTEle17TrackIsolFilterEta->at(itag),
				 HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTEle17TrackIsolFilterPhi->at(itag),
				 HLTEle32CaloIdTCaloIsoTTrkIdTTrkIsoTEle17TrackIsolFilterEnergy->at(itag)
				 );	    

		CreateAndFillUserTH1D("DR_ProbeVsTag", 100, 0, 10, probe.DeltaR(tag) );

		//-----------------
		//if the tag matches in DR with the probe --> move to the next tag candidate
		if( probe.DeltaR(tag) < 0.5) 
		  continue;

//.........这里部分代码省略.........
    truthpid = truth->get_pid();
    
    truth_g4particles->Fill();

  }


  /***********************************************

  GET THE EMCAL CLUSTERS

  ************************************************/


  RawClusterContainer::ConstRange begin_end = clusters->getClusters();
  RawClusterContainer::ConstIterator clusiter;
  
  for(clusiter = begin_end.first; clusiter!=begin_end.second; ++clusiter){

    RawCluster *cluster = clusiter->second;

    clus_energy = cluster->get_energy();
    clus_eta = cluster->get_eta();
    clus_theta = 2.*TMath::ATan((TMath::Exp(-1.*clus_eta)));
    clus_pt = clus_energy*TMath::Sin(clus_theta);
    clus_phi = cluster->get_phi();
    
    if(clus_pt<0.5)
      continue;

   
    
    TLorentzVector *clus = new TLorentzVector();
    clus->SetPtEtaPhiE(clus_pt,clus_eta,clus_phi,clus_energy);
    
    float dumarray[4];
    clus->GetXYZT(dumarray);
    clus_x = dumarray[0];
    clus_y = dumarray[1];
    clus_z = dumarray[2];
    clus_t = dumarray[3];

    clus_px = clus_pt*TMath::Cos(clus_phi);
    clus_py = clus_pt*TMath::Sin(clus_phi);
    clus_pz = sqrt(clus_energy*clus_energy-clus_px*clus_px-clus_py*clus_py);

    cluster_tree->Fill();

    //only interested in high pt photons to be isolated
    if(clus_pt<mincluspt)
      continue;
    if(fabs(clus_eta)>(1.0-isoconeradius))
      continue;

    float energysum = ConeSum(cluster,clusters,trackmap,isoconeradius);
    bool conecut = energysum > 0.1*clus_energy;
    if(conecut)
      continue;
   
    isolated_clusters->Fill();
    
    GetRecoHadronsAndJets(cluster, trackmap, reco_jets,recoeval);

  }

void Ntp1Analyzer_HWWlvjj::Loop(){

   DEBUG_VERBOSE_ = false;

   if (fChain == 0) return;

   Long64_t nentries;

   if( DEBUG_ ) nentries = 100000;
   else nentries = fChain->GetEntries();

   Long64_t nbytes = 0, nb = 0;

   TRandom3 rand;
   float Cont_inclusive=0., Cont_PV=0., Cont_MU=0., Cont_ELE=0., Cont_VetoMU=0., Cont_VetoELE=0., Cont_JetsELE=0., Cont_JetsMU=0.;
   Long64_t Jentry;

requiredTriggerElectron.push_back("1-164237:HLT_Ele27_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT_v");
requiredTriggerElectron.push_back("165085-166967:HLT_Ele32_CaloIdVT_CaloIsoT_TrkIdT_TrkIsoT_v");
requiredTriggerElectron.push_back("166968-999999:HLT_Ele52_CaloIdVT_TrkIdT_v");
requiredTriggerMuon.push_back("1-163261:HLT_Mu15_v");
requiredTriggerMuon.push_back("163262-164237:HLT_Mu24_v");
requiredTriggerMuon.push_back("165085-166967:HLT_Mu30_v");
requiredTriggerMuon.push_back("163262-166967:HLT_IsoMu17_v");
requiredTriggerMuon.push_back("167039-999999:HLT_IsoMu20_eta2p1_v");
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
     Long64_t ientry = LoadTree(jentry);
     if (ientry < 0) break;
     nb = fChain->GetEntry(jentry);   nbytes += nb;
     // if (Cut(ientry) < 0) continue;

     Jentry=jentry;

     if( DEBUG_VERBOSE_ ) std::cout << "entry n." << jentry << std::endl;
     

     if( (jentry%100000) == 0 ) std::cout << "Event #" << jentry  << " of " << nentries << std::endl;
  
   //HLT_Mu11_ = this->PassedHLT("HLT_Mu11");
   //HLT_Ele17_SW_EleId_L1R_ = this->PassedHLT("HLT_Ele17_SW_EleId_L1R");
   //HLT_DoubleMu3_ = this->PassedHLT("HLT_DoubleMu3");

     run_ = runNumber;
     LS_ = lumiBlock;
     event_ = eventNumber;
     eventWeight_ = -1.; //default

     Cont_inclusive++;//Other;

     if( !isGoodEvent( jentry ) ) continue; //this takes care also of integrated luminosity and trigger

       if( nPV==0 ) continue;
     bool goodVertex = (ndofPV[0] >= 4.0 && sqrt(PVxPV[0]*PVxPV[0]+PVyPV[0]*PVyPV[0]) < 2. && fabs(PVzPV[0]) < 24. );
     if( !goodVertex ) continue;
     nvertex_ = nPV;
     rhoPF_ = rhoFastjet;

     Cont_PV++; //Other

     //trigger:
     // not yet

     if( !isMC_ ){
     reloadTriggerMask(runNumber);
     bool passedElectronTrigger=true, passedMuonTrigger=true;
     std::string SingEle("SingleElectron_6july");//WW
     std::string SingMuo("SingleMu_6july");
     int elect=dataset_.compare(SingEle), muo=dataset_.compare(SingMuo);
     if( elect==0 ){ passedElectronTrigger = hasPassedHLT(0); }
     if( muo==0 ){ passedMuonTrigger = hasPassedHLT(1); }
     if( elect!=0 && muo!=0  ){ std::cout<<"Trigger doesn't works; dataset wrong"<<std::endl; }
     if( !passedElectronTrigger && elect==0 ){ continue;}
     if( !passedMuonTrigger && muo==0 ){ continue;}
     }
         ptHat_ = (isMC_) ? genPtHat : ptHat_;

     //if( isMC_ ) 
     //  if( (ptHat_ > ptHatMax_) || (ptHat_ < ptHatMin_) ) continue;

     bool noLeptons = false;
     TLorentzVector lept1MC, lept2MC;
     int wIndexqq=-1;
     int wIndexll=-1;
   
     if( isMC_ ) {
       
       // first look for W->qq'
       std::vector<TLorentzVector> quarksMC;

       for( unsigned iMc=0; iMc<nMc && quarksMC.size()<2; ++iMc ) {

         // quarks have status 3
         if( statusMc[iMc] != 3 ) continue;

         TLorentzVector* thisParticle = new TLorentzVector();
         thisParticle->SetPtEtaPhiE( pMc[iMc]*sin(thetaMc[iMc]), etaMc[iMc], phiMc[iMc], energyMc[iMc] );

         if( fabs(idMc[iMc])<7 && fabs(idMc[mothMc[iMc]])==24 ) {
           wIndexqq = mothMc[iMc];
           quarksMC.push_back( *thisParticle );
//.........这里部分代码省略.........

Bool_t Analysis::Process(Long64_t entry)
{

  jpt  ->clear();    
  jeta ->clear(); 
  jphi ->clear(); 
  jmass->clear();
  jconst->clear();
  toppt->clear();
  jhasb->clear();
  jhasq->clear();
  
  if (entry%1==0) cout << "\r Events  " << entry << "\r" << flush;
  this->GetEntry(entry);
    
  double rho = 2*173.5;
  double min_r = 0.4;
  double max_r = 1.5;
  double ptmin = 200.;

  // fill jets
  //-----------------------------------------------
  vector<fastjet::PseudoJet> pseudojets; pseudojets.clear();
  for (int ijet=0;ijet<jets_n;++ijet){
    if(jets_isBad->at(ijet)!=0) return kFALSE;
    if(jets_pt->at(ijet)<20.) continue;
    if(fabs(jets_eta->at(ijet))>2.8) continue;
    TLorentzVector tvit = TLorentzVector();
    tvit.SetPtEtaPhiE(jets_pt->at(ijet),
		      jets_eta->at(ijet),
		      jets_phi->at(ijet),
		      jets_e->at(ijet));
    pseudojets.push_back( fastjet::PseudoJet(tvit.Px(), tvit.Py(), 
					     tvit.Pz(), tvit.E()) );
  }

  // fill truth
  //----------------------------------------------
  vector<int> tops;
  vector<int> Wq1s;
  vector<int> Wq2s;
  vector<int> bs;
  for (int imc=0;imc<mc_n;++imc)
    if (fabs(mc_pdgId->at(imc))==6){
      int windex = mc_children_index->at(imc)[0];
      int bindex = mc_children_index->at(imc)[1];
      int wq1 = mc_children_index->at(windex)[0];
      int wq2 = mc_children_index->at(windex)[1];
      if(fabs(mc_pdgId->at(wq1))>5 || fabs(mc_pdgId->at(wq2))>5) continue;
      tops.push_back( imc );
      Wq1s.push_back( wq1 );
      Wq2s.push_back( wq2 );
      bs.push_back( bindex );
      toppt->push_back( mc_pt->at(imc) );
    }//endif

  for (int it=0;it<tops.size();++it){
    fastjet::PseudoJet wq1; wq1.reset_PtYPhiM(1e-10,mc_eta->at(Wq1s[it]),mc_phi->at(Wq1s[it]));
    fastjet::PseudoJet wq2; wq2.reset_PtYPhiM(1e-10,mc_eta->at(Wq2s[it]),mc_phi->at(Wq2s[it]));
    fastjet::PseudoJet bq; bq.reset_PtYPhiM(1e-10,mc_eta->at(bs[it]),mc_phi->at(bs[it]));
    wq1.set_user_info(new MyUserInfo(it, false));
    wq1.set_user_info(new MyUserInfo(it, false));
    bq.set_user_info(new MyUserInfo(it, true));
    pseudojets.push_back( wq1 );
    pseudojets.push_back( wq2 );
    pseudojets.push_back( bq );
  }// top loop

  fastjet::contrib::VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, fastjet::contrib::VariableRPlugin::AKTLIKE);
  fastjet::JetDefinition jet_def(&lvjet_pluginAKT);
  fastjet::ClusterSequence clust_seqAKT(pseudojets, jet_def);
  vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);

  for (int ijet=0;ijet<inclusive_jetsAKT.size();++ijet){
    fastjet::PseudoJet *jet = &(inclusive_jetsAKT[ijet]);
    jpt->push_back( jet->pt() );
    jeta->push_back( jet->eta() );
    jphi->push_back( jet->phi() );
    jmass->push_back( jet->m() );
    int ntruth=0;
    vector<int> hasb;
    vector<int> hasq;
    for (int iconst=0;iconst<jet->constituents().size();++iconst){
      fastjet::PseudoJet *jetconst = &(jet->constituents()[iconst]);
      if (jetconst->has_user_info<MyUserInfo>()){
	ntruth++;
	int topid = jetconst->user_info<MyUserInfo>().top();
	bool isbq = jetconst->user_info<MyUserInfo>().isbq();
	if(isbq) hasb.push_back(topid);
	else hasq.push_back(topid);
      }//endif is truth
    }//loop constituents
    jconst->push_back( jet->constituents().size()-ntruth );
    jhasb->push_back( hasb );
    jhasq->push_back( hasq );
  }// loop reclustered jets

  outtree->Fill();

  return kTRUE;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  tree->Branch("Ele2Phi", &Ele2Phi, "Ele2Phi/D");
   tree->Branch("Ele1Enr", &Ele1Enr, "Ele1Enr/D");
    tree->Branch("Ele2Enr", &Ele2Enr, "Ele2Enr/D");
  tree->Branch("ZMass", &ZMass, "ZMass/D");
  tree->Branch("ZPT", &ZPT, "ZPT/D");
  tree->Branch("ZEta", &ZEta, "ZEta/D");
  tree->Branch("ZRapidity", &ZRapidity, "ZRapidity/D");
  tree->Branch("ZPhi", &ZPhi, "ZPhi/D");
  tree->Branch("BB", &BB, "BB/B");
  tree->Branch("BE", &BE, "BE/B");
  tree->Branch("EE", &EE, "EE/B");

  vector <double> newelePt; vector <double> neweleEta; vector <double> neweleEnr; vector <double> newelePhi;vector <double> newscEta;
  vector <double> elePt; vector <double> eleEta; vector <double> eleEnr; vector <double> elePhi; vector <double> elescEta;
  TClonesArray *eleP4 = new TClonesArray("TLorentzVector", 1500);

  int nentries = T1->GetEntries();
  //int nentries = 10000;
  cout<<"entries: "<<nentries<<endl;
  for (unsigned int jentry=0; jentry < nentries; jentry++) {
    T1->GetEntry(jentry);

    if(jentry%1000000 == 0){
      cout << "Events Processed :  " << jentry << endl;
    }

    trigMatch_lead = false; trigMatch_slead = false;
    dR = 0.; dR1 = 0.; dR2 = 0.;
    BB=false; BE=false; EE=false;
    mediumId = 0; dzMaskId = 0; passId =0;

    newelePt.clear(); neweleEta.clear(); newelePhi.clear(); neweleEnr.clear();newscEta.clear();
    elePt.clear(); eleEta.clear(); elePhi.clear(); eleEnr.clear();elescEta.clear();
    recoElec.SetPtEtaPhiE(0.,0.,0.,0.);
    eleP4->Clear();

    primVtx = nPV;

    if(!Ele23_WPLoose) continue;

    for(int i=0;i<ptElec->size();i++){

      recoElec.SetPtEtaPhiE(ptElec->at(i),etaElec->at(i), phiElec->at(i), energyElec->at(i));
      new ((*eleP4)[i]) TLorentzVector(recoElec);
      TLorentzVector* fourmom = (TLorentzVector*) eleP4->At(i);

   //  if(fourmom->Eta() < 1.4442) *fourmom *= 0.9994;
   //   if(fourmom->Eta() > 1.566) *fourmom *= 1.0034;

      if(fourmom->Eta() < 1.4442) *fourmom *= 1.0;
      if(fourmom->Eta() > 1.566) *fourmom *= 1.0;

     //   mediumId = passHEEPId->at(i);
       mediumId = passMediumId->at(i);

       passId = mediumId;

//      if (ptElec->at(i) < 500.) { passId = mediumId; }
//     else { passId = dzMaskId; }

      if(passId) {
	if(fabs(etaSC->at(i)) < 2.5 && !(fabs(etaSC->at(i)) > 1.4442 && fabs(etaSC->at(i)) < 1.566)){

	  newelePt.push_back(fourmom->Pt());
	  neweleEta.push_back(fourmom->Eta());
	  neweleEnr.push_back(fourmom->Energy());

//.........这里部分代码省略.........
      for( int iele = 0 ; iele < fTree.electrons_size;iele ++ ){
      	if( fTree.electrons_Pt[iele] > 20 &&
	    fabs(fTree.electrons_Eta[iele]) < 2.5 && 
	    ( fabs(fTree.electrons_Eta[iele]) < 1.4442 || fabs(fTree.electrons_Eta[iele]) > 1.566) &&
	    fTree.electrons_vidVeto[iele] > 0.5 )
	  nLooseElectrons++;
      }
      
      if( nLooseElectrons > 0 )
      	continue;

      if( printEventIds )
	(*(EventIdFiles[cutindex])) << (fTree.Event_EventNumber) << endl ;
      cutflowtable.Fill( cutindex , weight , EventsIsPSeudoData < fabs(weight) );
      cutflowtablew1.Fill( cutindex , weight/fabs(weight) , EventsIsPSeudoData < fabs(weight) );
      cutindex ++ ;

      int j1index = -1;
      int j2index = -1;
      int j3index = -1;
      int nJets = 0 ;
      int howManyBJets = 0;
      int bjIndex = -1;
      int bj2Index = -1;
      int jprimeIndex ; 

      int nJetsPt20_47 = 0;

      int nJetsPt20_24 = 0;
      int nLbJetsPt20 = 0;
      int nTbJetsPt20 = 0;

      TLorentzVector muon;
      muon.SetPtEtaPhiE( fTree.muons_Pt[tightMuIndex] , fTree.muons_Eta[tightMuIndex] , fTree.muons_Phi[tightMuIndex] , fTree.muons_E[tightMuIndex] );
      if( muon.Energy() == 0.0 )
	cout << "ZERO???" <<endl;
      for( int jid = 0 ; jid < fTree.jetsAK4_size ; jid++ ){

	float NHF = fTree.jetsAK4_JetID_neutralHadronEnergyFraction[jid] ;
	float NEMF = fTree.jetsAK4_JetID_neutralEmEnergyFraction[jid] ;
	float NumConst = fTree.jetsAK4_JetID_numberOfDaughters[jid] ;
	float eta = fTree.jetsAK4_Eta[jid] ;
	float CHF = fTree.jetsAK4_JetID_chargedHadronEnergyFraction[jid] ;
	float CHM = fTree.jetsAK4_JetID_chargedMultiplicity[jid] ;
	float CEMF = fTree.jetsAK4_JetID_chargedEmEnergyFraction[jid] ;
	float NumNeutralParticle = fTree.jetsAK4_JetID_neutralMultiplicity[jid] ;


	bool looseid ;
	if(  abs(eta)<=3.0 ) {
	  looseid = (NHF<0.99 && NEMF<0.99 && NumConst>1) && ((abs(eta)<=2.4 && CHF>0 && CHM>0 && CEMF<0.99) || abs(eta)>2.4) ;
	}
	else
	  looseid = (NEMF<0.90 && NumNeutralParticle>10) ;

      	if( fTree.jetsAK4_CorrPt[jid] > 20 &&
	    fabs( fTree.jetsAK4_Eta[jid] ) < 4.7 &&
	    looseid
	    ){
	  TLorentzVector jet;
	  jet.SetPtEtaPhiE( fTree.jetsAK4_CorrPt[jid] , fTree.jetsAK4_Eta[jid] , fTree.jetsAK4_Phi[jid] , fTree.jetsAK4_CorrE[jid] );
	  double DR = muon.DeltaR( jet );

	  if ( DR > 0.3 ){

	    if( fTree.jetsAK4_CorrPt[jid] <= 40 ){

//.........这里部分代码省略.........
	  if (!isTagTightEle[iEle])
	    continue;
	}
      else if (tagTightnessLevel=="Medium")
	{
	  if (!isTagMediumEle[iEle])
	    continue;
	}
      else if (tagTightnessLevel=="Loose")
	{
	  if (!isTagLooseEle[iEle])
	    continue;
	}
      else
	{
	  std::cout << "Tag Selection undefined" << std::endl;
	  exit(-1);
	}

      TLorentzVector theEle;
      //      theEle.SetPtEtaPhiE(electron_pt[iEle], electron_eta[iEle], electron_phi[iEle], electron_energy[iEle]);
      theEle.SetPtEtaPhiM(electron_pt[iEle], electron_eta[iEle], electron_phi[iEle], 0.);

      for (int iPho=0; iPho<nPhot; iPho++) {
	if (isMC && mcMatch)
	  if (!isGenMatchPhot[iPho])
	    continue;

	//apply preselection as probe selection for photons (should be loosened to measure preselection SF)
	if (!isProbePreselPhot[iPho])
	  continue;

	TLorentzVector thePho;
	thePho.SetPtEtaPhiE(ptPhot[iPho], etaPhot[iPho], phiPhot[iPho], ePhot[iPho]);
      
	TLorentzVector theTaP = theEle + thePho;
	float theMass = theTaP.M();
      
	if (fabs(theMass-91.181)>DeltaMZ)
	  continue;

	// filling the tree
	mass = theMass;
	
	probe_eta    = etascPhot[iPho];
	probe_abseta = fabs(etascPhot[iPho]);
	probe_phi    = phiPhot[iPho];
	probe_pt     = ptPhot[iPho];
	probe_r9     = r9Phot[iPho];
	
	numvtx = nvtx;
	rho=rhoAllJets;

	if(isMC) {
	  puW   = pu_weight;
	  puW30 = pu_weight30;
	  puW50 = pu_weight50;
	  puW75 = pu_weight75;
	  puW90 = pu_weight90;
	} else {
	  puW   = 1.;
	  puW30 = 1.;
	  puW50 = 1.;
	  puW75 = 1.;
	  puW90 = 1.;
	}

//.........这里部分代码省略.........
	TH1F * pz_comparison5_2 = new TH1F("pz_comparison5_2","pz_comparison",200,-1,1);
	TH1F * pz_comparison5_1 = new TH1F("pz_comparison5_1","pz_comparison",200,-1,1);
	TH1F * pz_comparison6 = new TH1F("pz_comparison6","pz_comparison",200,-1,1);
	TH1F * lvjj_comparison1 = new TH1F("lvjj_comparison1","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison2 = new TH1F("lvjj_comparison2","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison2_1 = new TH1F("lvjj_comparison2_1","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison3 = new TH1F("lvjj_comparison3","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison4 = new TH1F("lvjj_comparison4","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison5 = new TH1F("lvjj_comparison5","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison5_1 = new TH1F("lvjj_comparison5_1","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison5_2 = new TH1F("lvjj_comparison5_2","lvjj_comparison",200,-1,1);
	TH1F * lvjj_comparison6 = new TH1F("lvjj_comparison6","lvjj_comparison",200,-1,1); 

	TCanvas * c1 = new TCanvas();	

	TLorentzVector lep;
	TLorentzVector jet1;
	TLorentzVector jet2;
	TLorentzVector neu;
	TLorentzVector temp;

	double mlvjj;
	double lvjjmass_true;
	cout<<chain->GetEntries()<<endl;
	for(int i=0; i<chain->GetEntries(); i++)
	{
		//if(i>5000)break;//for test
		if(i%5000==0)cout<<"Processing the events: "<<i<<endl;
		chain->GetEntry(i);
		//higgs 2D
		//h1->Fill(phig->at(0)->eta(),phig->at(0)->pt(),normal/(float)chain->GetEntries());

		//------------------------------neutrino pz
		if(!pmu->empty())lep.SetPtEtaPhiE(pmu->at(0)->pt(), pmu->at(0)->eta(),pmu->at(0)->phi(),pmu->at(0)->energy());
		else if(!pele->empty())lep.SetPtEtaPhiE(pele->at(0)->pt(), pele->at(0)->eta(),pele->at(0)->phi(),pele->at(0)->energy());
		else {
				cout<<"skip event: "<<i<<"  no lepton here"<<endl;
				//cout<<phig->size()<<endl;
				continue;
			}
		if(pjet->size()!=2)
		{
			cout<<"skip event: "<<i<<" jet sise != 2 "<<endl;
			continue;
		}
		jet1.SetPtEtaPhiE(pjet->at(0)->pt(), pjet->at(0)->eta(),pjet->at(0)->phi(),pjet->at(0)->energy());
		jet2.SetPtEtaPhiE(pjet->at(1)->pt(), pjet->at(1)->eta(),pjet->at(1)->phi(),pjet->at(1)->energy());
		if(jet1.Pt()<jet2.Pt())
		{
			temp=jet1;
			jet1=jet2;
			jet2=temp;
		}
		neu.SetPxPyPzE(pnu->at(0)->px(), pnu->at(0)->py(),pnu->at(0)->pz(),pnu->at(0)->energy());//real neutrino
		double px = pnu->at(0)->px();
		double py = pnu->at(0)->py();
		double pz=0;


		// ------------------------------- make cuts----------------------------
		//mt
		//double mt = sqrt(2. * lep.Pt() * neu.Pt() * ( 1 - cos(deltaPhi(lep.Phi(), neu.Phi()) ) ) );
		//if(mt>20)continue;

		//lepton pt,eta
		if(lep.Pt()<50)continue;

//.........这里部分代码省略.........
  tree_->SetBranchAddress("ePart", ePart);
  Float_t ptPart[20];
  tree_->SetBranchAddress("ptPart", ptPart);
  Float_t etaPart[20];
  tree_->SetBranchAddress("etaPart", etaPart);
  Float_t phiPart[20];
  tree_->SetBranchAddress("phiPart", phiPart);
  Int_t pdgIdPart[20];
  tree_->SetBranchAddress("pdgIdPart", pdgIdPart);


  //QGLikelihoodCalculator *qglikeli = new QGLikelihoodCalculator("QG_QCD_Pt_15to3000_TuneZ2_Flat_7TeV_pythia6_Fall10.root", nBins);
  //QGLikelihoodCalculator *qglikeli = new QGLikelihoodCalculator("QG_QCD_Pt_15to3000_TuneZ2_Flat_7TeV_pythia6_Spring11-PU_S1_START311_V1G1-v1.root", nBins);
  QGLikelihoodCalculator *qglikeli = new QGLikelihoodCalculator("QG_QCD_Pt-15to3000_TuneZ2_Flat_7TeV_pythia6_Summer11-PU_S3_START42_V11-v2.root", nBins);


  int nEntries = tree_->GetEntries();
  //nEntries = 10000;

  for(int iEntry=0; iEntry<nEntries; ++iEntry) {
  //for(int iEntry=0; iEntry<500000; ++iEntry) {

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

    tree_->GetEntry(iEntry);

    if( eventWeight <= 0. ) eventWeight = 1.;


    //for( unsigned iJet=0; iJet<nJet; ++iJet ) {
    for( unsigned iJet=0; (iJet<nJet && iJet<3); ++iJet ) { //only 3 leading jets considered

      TLorentzVector thisJet;
      thisJet.SetPtEtaPhiE( ptJet[iJet], etaJet[iJet], phiJet[iJet], eJet[iJet]);

      float deltaRmin=999.;
      int partFlavor=-1;
      TLorentzVector foundPart;

      for( unsigned iPart=0; iPart<nPart; iPart++ ) {

        TLorentzVector thisPart;
        thisPart.SetPtEtaPhiE( ptPart[iPart], etaPart[iPart], phiPart[iPart], ePart[iPart]);

        float thisDeltaR = thisJet.DeltaR(thisPart);

        if( thisDeltaR < deltaRmin ) {
          deltaRmin = thisDeltaR;
          partFlavor = pdgIdPart[iPart];
          foundPart = thisPart;
        }

      } //for partons

      if( deltaRmin > 0.5 ) continue;
      //if( deltaRmin > 0.5 ) partFlavor=21; //lets try this


      int thisBin=-1;
      if( thisJet.Pt() > ptBins[nBins] ) {
        thisBin = nBins-1;
      } else {
        for( unsigned int iBin=0; iBin<nBins; ++iBin ) {
          if( thisJet.Pt()>ptBins[iBin] && thisJet.Pt()<ptBins[iBin+1] ) {
            thisBin = iBin;
            break;