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

// // //
double analysisClass::METlepDeltaPhi(TString Lep, int iLep){
  TLorentzVector lep;
  lep.SetPtEtaPhiM( 0, 0, 0, 0 );
  TLorentzVector met;
  met.SetPtEtaPhiM( METcorr("Pt"), 0, METcorr("Phi"), 0 );
  //std::cout<< "analysisClass::METlepDeltaPhi  iLep: "<<iLep<<std::endl;
  //
  if( Lep != "Mu" && Lep != "Tau" ){ cout<<" WRONG LEPTON TYPE SPECIFIED!! analysisClass::METlepDeltaPhi(TString Lep, int iLep) "<<endl; return 0; }
  //
  if( iLep>-1 ){
    if( Lep=="Mu" ) lep.SetPtEtaPhiM( muPtcorr(iLep), MuonEta->at(iLep), MuonPhi->at(iLep), 0 );
    if( Lep=="Tau") lep.SetPtEtaPhiM( tauPtcorr(iLep), HPSTauEta->at(iLep), HPSTauPhi->at(iLep), 0 );
  }
  //
  if( iLep==-1 ){
    if( Lep=="Mu"){
      for(unsigned int iMuR=0;  iMuR<MuonPt->size();     iMuR++){
	if(!muRisoCheck(iMuR))continue;
	if( lep.Pt()<muPtcorr(iMuR)     ) lep.SetPtEtaPhiM( muPtcorr(iMuR), MuonEta->at(iMuR), MuonPhi->at(iMuR), 0 );
      }
    }
    //------ ------
    if( Lep=="Tau"){
      for(unsigned int iTauR=0; iTauR<HPSTauPt->size(); iTauR++){
	if(!tauRisoCheck(iTauR))continue;
	if( lep.Pt()<tauPtcorr(iTauR)  ) lep.SetPtEtaPhiM( tauPtcorr(iTauR), HPSTauEta->at(iTauR), HPSTauPhi->at(iTauR), 0 );
      }
    }
  }
  //
  return fabs(met.DeltaPhi(lep));
}

// This is the pt corrected delta phi between the 2 leptons
// P and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the 
// leptons for you
double HWWKinematics::CalcDeltaPhiRFRAME(){
  // first calculate pt-corrected MR
  float mymrnew = CalcMRNEW();
  
  // Now, boost lepton system to rest in z
  // (approximate accounting for longitudinal boost)
  TVector3 BL = L1.Vect()+L2.Vect();
  BL.SetX(0.0);
  BL.SetY(0.0);
  BL = (1./(L1.P()+L2.P()))*BL;
  L1.Boost(-BL);
  L2.Boost(-BL);
  
  // Next, calculate the transverse Lorentz transformation
  // to go to Higgs approximate rest frame
  TVector3 B = L1.Vect()+L2.Vect()+MET;
  B.SetZ(0.0);
  B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
  
  L1.Boost(B);
  L2.Boost(B);
  
  //Now, re-calculate the delta phi
  // in the new reference frame:
  return L1.DeltaPhi(L2);
  
}

double SUSYLooperHists::DiTau_InvMass( TLorentzVector Met, TLorentzVector L1, TLorentzVector L2, float Al ){   
  TLorentzVector T1,T2;
  double DiTauMass; 
  TMatrixF A(2,2); 
  TVectorF C(2),X(2); 
  A(0,0)=L1.Px();
  A(0,1)=L2.Px();
  A(1,0)=L1.Py();
  A(1,1)=L2.Py();
  A=A.Invert();
  C(0)=(Met+L1+L2).Px();
  C(1)=(Met+L1+L2).Py();
  X=A*C;// double X0i=X(0), X1i=X(1);
  //---------------[ MET ReAlignement subsection ]------------------------------
  if(X(0)<0||X(1)<0){ 
    if     ( fabs(L1.DeltaPhi(Met))>Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) {/*DO NOTHING just normaly a non-Z event!*/}
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
    else if( fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))<fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
    else if( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
  }//---------------------------------------------------------------------------
  C(0)=(Met+L1+L2).Px(); C(1)=(Met+L1+L2).Py(); X=A*C;
  T1.SetPxPyPzE( L1.Px()*X(0), L1.Py()*X(0), L1.Pz()*X(0), sqrt( 3.1571 +L1.P()*L1.P()*X(0)*X(0) ) );
  T2.SetPxPyPzE( L2.Px()*X(1), L2.Py()*X(1), L2.Pz()*X(1), sqrt( 3.1571 +L2.P()*L2.P()*X(1)*X(1) ) );
  if( X(0)>0 && X(1)>0 ) DiTauMass=(T1+T2).M();  else DiTauMass=-(T1+T2).M();  return DiTauMass;
  //if((X(0)!=X0i||X(1)!=X1i))std::cout<<X(0)<<" "<<X(1)<<" <--"<<X0i<<" "<<X1i<<" RMETal.phi="<<(T1+T2-L1-L2).Phi()<<" RMETal.eta"<<(T1+T2-L1-L2).Eta()<<" MZ="<<DiTauMass<<endl; 
}

//___________________[ TAUS RECONSTRUCTION AND INVMASS CALCULATION ]_________________________________________
double SUSYLooperHistsSoftBase::DiTau_InvMass( TLorentzVector Met, TLorentzVector L1, TLorentzVector L2, float Al ) {

    TLorentzVector T1,T2;    double DiTauMass;    TMatrixF A(2,2);    TVectorF C(2),X(2);
    A(0,0)=L1.Px();
    A(0,1)=L2.Px();
    A(1,0)=L1.Py();
    A(1,1)=L2.Py();
    A=A.Invert();
    C(0)=(Met+L1+L2).Px();
    C(1)=(Met+L1+L2).Py();
    X=A*C;
    //double X0i=X(0), X1i=X(1);
    if ( (X(0)<0.||X(1)<0.) && Al>0. ) {//---[MET Alignement subsection]--------------------------------------------------------------------------------------------
        if      ( fabs(L1.DeltaPhi(Met))>Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) {}//{DO NOTHING just normaly a non-Z event!}
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
        else if ( fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>Al                                                  ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))<fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L1.Phi(),0);
        else if ( fabs(L1.DeltaPhi(Met))<Al && fabs(L2.DeltaPhi(Met))<Al && fabs(L1.DeltaPhi(Met))>fabs(L2.DeltaPhi(Met)) ) Met.SetPtEtaPhiM(Met.Pt(),0,L2.Phi(),0);
    }//-------------------------------------------------------------------------------------------------------------------------------------------------------------
    C(0)=(Met+L1+L2).Px();
    C(1)=(Met+L1+L2).Py();
    X=A*C;
    T1.SetPxPyPzE( L1.Px()*X(0), L1.Py()*X(0), L1.Pz()*X(0), sqrt( 3.1571 +L1.P()*L1.P()*X(0)*X(0) ) );
    T2.SetPxPyPzE( L2.Px()*X(1), L2.Py()*X(1), L2.Pz()*X(1), sqrt( 3.1571 +L2.P()*L2.P()*X(1)*X(1) ) );
    if (  X(0)>0.  &&  X(1)>0.   )   DiTauMass=(T1+T2).M();
    //if (  (X(0)>0. && 0.<X(1)&&X(1)<1.)  ||  (X(1)>0. && 0.<X(0)&&X(0)<1.)  )  DiTauMass=(T1+T2).M(); // B
    //if (  (X(0)>1.&& X(1)<0.) || (X(1)>1.&& X(0)<0.)    )  DiTauMass=(T1+T2).M(); // C
    //if (  X(0)<0.  &&  X(1)<0.   )  DiTauMass=(T1+T2).M(); //  D    
    else DiTauMass=-(T1+T2).M();
    return DiTauMass;
}//-----------------------------------------------------------------------------------------------------------

// This is the pt corrected delta phi between the 2 mega-jets
// P and Q are the 4-vectors for the 2 hemispheres 
// M is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the 
// leptons for you
double HWWKinematics::CalcDeltaPhiNEW(TLorentzVector P, TLorentzVector Q, TVector3 M){
    // first calculate pt-corrected MR
 float mymrnew = CalcMRNEW(L1,L2,MET);

    //Next, calculate the transverse Lorentz transformation
 TVector3 B = P.Vect()+Q.Vect()+MET;
 B.SetZ(0.0);
 B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;

 P.Boost(B);
 Q.Boost(B);

    //Now, re-calculate the delta phi
    // in the new reference frame:

 return P.DeltaPhi(Q);

}

// // //
double analysisClass::METlepMT(TString Lep, int iLep){
  TLorentzVector MET;
  TLorentzVector lep;
  MET.SetPtEtaPhiM( 0 , 0 , 0 , 0 );
  lep.SetPtEtaPhiM( 0 , 0 , 0 , 0 );
  double maxPt=0;
  //std::cout<< "analysisClass::METlepMT  iLep: "<<iLep<<std::endl;
  //
  if( Lep != "Mu" && Lep != "Tau" ){ cout<<" WRONG LEPTON TYPE SPECIFIED!! analysisClass::METlepMT(TString Lep, int iLep) "<<endl; return 0; }
  //
  if( iLep>-1 ){
    if( Lep == "Mu"  ) lep.SetPtEtaPhiM( muPtcorr(iLep), MuonEta->at(iLep), MuonPhi->at(iLep), 0 );
    if( Lep == "Tau" ) lep.SetPtEtaPhiM( tauPtcorr(iLep), HPSTauEta->at(iLep), HPSTauPhi->at(iLep), 0 );
  }
  if( iLep==-1){
    if( Lep == "Mu"  ){
      maxPt=0;
      for(unsigned int iMuR=0; iMuR<MuonPt->size(); iMuR++){
	if(!muRisoCheck(iMuR))continue;
	if( muPtcorr(iMuR)>maxPt ){ maxPt=muPtcorr(iMuR);  lep.SetPtEtaPhiM( muPtcorr(iMuR), MuonEta->at(iMuR), MuonPhi->at(iMuR), 0 ); }
      }
    }
    //------- -------
    if( Lep == "Tau" ){ 
      maxPt=0;
      for(unsigned int iTauR=0; iTauR<HPSTauPt->size(); iTauR++){
	if(!tauRisoCheck(iTauR))continue;
	if( tauPtcorr(iTauR)>maxPt ){ maxPt=tauPtcorr(iTauR);  lep.SetPtEtaPhiM( tauPtcorr(iTauR), HPSTauEta->at(iTauR), HPSTauPhi->at(iTauR), 0 ); }
      }
    }
  }
  //
  MET.SetPtEtaPhiM( METcorr("Pt"), 0 , METcorr("Phi") , 0 );
  double M_T=TMath::Sqrt(2*fabs(MET.Pt())*fabs(lep.Pt())*(1-TMath::Cos(MET.DeltaPhi(lep))));
  //
  return M_T;
}

//.........这里部分代码省略.........
			int min = TMath::LocMin(4, &muSum[0]);

			if( (min == 0 || min == 1) ){
				if(mZa > mZb) { Z1.SetVectM( Za, mZa); Z2.SetVectM(Zb,mZb); 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());}	//to set the highest mass the z
				else { Z2.SetVectM( Za, mZa); Z1.SetVectM(Zb,mZb); 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());}
				
			}
			else if( (min == 2 || min == 3) ){
				if(mZc > mZd) { Z1.SetVectM( Zc, mZc); Z2.SetVectM(Zd,mZd); lep_min1.SetPtEtaPhiE(lep1.Pt(),lep1.Eta(), lep1.Phi(),lep1.Mag()); lep_plus1.SetPtEtaPhiE(lep4.Pt(),lep4.Eta(), lep4.Phi(),lep4.Mag()); lep_min2.SetPtEtaPhiE(lep3.Pt(),lep3.Eta(), lep3.Phi(),lep3.Mag()); lep_plus2.SetPtEtaPhiE(lep2.Pt(),lep2.Eta(), lep2.Phi(),lep2.Mag());}	//to set the highest mass the z
				else { Z2.SetVectM( Zc, mZc); Z1.SetVectM(Zd,mZd); lep_min2.SetPtEtaPhiE(lep1.Pt(),lep1.Eta(), lep1.Phi(),lep1.Mag()); lep_plus2.SetPtEtaPhiE(lep4.Pt(),lep4.Eta(), lep4.Phi(),lep4.Mag()); lep_min1.SetPtEtaPhiE(lep3.Pt(),lep3.Eta(), lep3.Phi(),lep3.Mag()); lep_plus1.SetPtEtaPhiE(lep2.Pt(),lep2.Eta(), lep2.Phi(),lep2.Mag());}


			}
		xxxx++;
		fillFlag = true;
		}


		if ( fillFlag == true && goodEventCounter < 25001) {	//if it fullfills the specs then fill and find angles

			rec_H = Z1 + Z2;
			double Hmass = rec_H.M();
			tHrec = Hmass;
// 			cout<<tHrec<<endl;

			double Z1mass = Z1.M();
			tZ1m = Z1mass;
			double Z2mass = Z2.M();
			tZ2m = Z2mass;
			double ptlepp1 = lep_plus1.Pt();
			double ptlepm1 = lep_min1.Pt();
			double ptlepp2 = lep_plus2.Pt();
			double ptlepm2 = lep_min2.Pt();
			double dR1 = sqrt(pow(fabs(lep_min1.Eta() - lep_plus1.Eta()),2)+pow(fabs(lep_min1.DeltaPhi(lep_plus1)),2));
			double dR2 = sqrt(pow(fabs(lep_min2.Eta() - lep_plus2.Eta()),2)+pow(fabs(lep_min2.DeltaPhi(lep_plus2)),2));

// 			if ( /*Hmass<120 || Hmass>130 || */Z1mass < 49 || Z1mass>107 || Z2mass < 12 || Z2mass> 115 ){continue;}	//constrains		
	
			//filling the simple histogram values
			h_Z1_m -> Fill(Z1.M());
			h_Z1_E -> Fill(Z1.E());
			h_Z1_Pt -> Fill(Z1.Pt());
			h_Z1_eta -> Fill(Z1.Eta());
			h_Z1_phi -> Fill(Z1.Phi());
	
			h_Z2_m -> Fill(Z2.M());
			h_Z2_E -> Fill(Z2.E());
			h_Z2_Pt -> Fill(Z2.Pt());
			h_Z2_eta -> Fill(Z2.Eta());
			h_Z2_phi -> Fill(Z2.Phi());
	
	
			h_rec_H_m	-> Fill(Hmass);
			h_rec_H_E	-> Fill(rec_H.E());
			h_rec_H_Pt	-> Fill(rec_H.Pt());
			h_rec_H_eta	-> Fill(rec_H.Eta());
			h_rec_H_phi	-> Fill(rec_H.Phi());	

			h_lep_plus1_E	-> Fill(lep_plus1.E());
			h_lep_plus1_Pt	-> Fill(ptlepp1);
			h_lep_plus1_eta	-> Fill(lep_plus1.Eta());
			h_lep_plus1_phi	-> Fill(lep_plus1.Phi());

			h_lep_min1_E	-> Fill(lep_min1.E());
			h_lep_min1_Pt	-> Fill(ptlepm1);
			h_lep_min1_eta	-> Fill(lep_min1.Eta());

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

        if(((TLorentzVector*)(*eventJets.p4)[nj])->Pt() > 20)  idJet20.push_back(nj);
        if(((TLorentzVector*)(*eventJets.p4)[nj])->Pt() > 30)  idJet30.push_back(nj);
      }

      if     (nsel == 0){ // Z->ll
        if(idLep.size() == 2) nPassCuts[2]++;
	if(idLep.size() == 2 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10 && 
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[1]])->Pt() > 10) nPassCuts[3]++;
        if(idLep.size() == 2 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10 && 
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[1]])->Pt() > 10 &&
          (int)(*eventLeptons.pdgId)[idLep[0]]*(int)(*eventLeptons.pdgId)[idLep[1]] < 0 &&
	  TMath::Abs((int)(*eventLeptons.pdgId)[idLep[0]]) == TMath::Abs((int)(*eventLeptons.pdgId)[idLep[1]]) &&
	  TMath::Abs((int)(*eventLeptons.pdgId)[idLep[0]]) == typeSel &&
	 (infilecatv[ifile] == 0 || isGenLep.size() == 2)) {
          nPassCuts[4]++;
          dilep = ( ( *(TLorentzVector*)(eventLeptons.p4->At(idLep[0])) ) + ( *(TLorentzVector*)(eventLeptons.p4->At(idLep[1])) ) ); 
          if(TMath::Abs(dilep.M()-91.1876)<15.0) passSel = kTRUE;	
	}
      }
      else if(nsel == 1){ // fake
        if(idLep.size() == 1) nPassCuts[2]++;
	if(idLep.size() == 1 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10) nPassCuts[3]++;
        if(idLep.size() == 1 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10 &&
	  ((typeSel == 11 && idJet30.size() >= 1) ||(typeSel == 13 && idJet20.size() >= 1)) &&
	  TMath::Abs((int)(*eventLeptons.pdgId)[idLep[0]]) == typeSel &&
	 (infilecatv[ifile] == 0 || isGenLep.size() == 1)) {
          dilep = ( *(TLorentzVector*)(eventLeptons.p4->At(idLep[0])) );
          nPassCuts[4]++;
          deltaPhiDileptonMet = TMath::Abs(dilep.DeltaPhi(*((TLorentzVector*)(*eventMet.p4)[0])));
          mtW = TMath::Sqrt(2.0*dilep.Pt()*((TLorentzVector*)(*eventMet.p4)[0])->Pt()*(1.0 - cos(deltaPhiDileptonMet)));
	  if(mtW < 20) passSel = kTRUE;
	}
      }
      else if(nsel == 2){ // W->ln
        if(idLep.size() == 1 && idTight[0] == 1) nPassCuts[2]++;
	if(idLep.size() == 1 && idTight[0] == 1 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10) nPassCuts[3]++;
        if(idLep.size() == 1 && idTight[0] == 1 &&
     	  ((TLorentzVector*)(*eventLeptons.p4)[idLep[0]])->Pt() > 10 &&
	  TMath::Abs((int)(*eventLeptons.pdgId)[idLep[0]]) == typeSel &&
	 (infilecatv[ifile] == 0 || isGenLep.size() == 1)) {
          dilep = ( *(TLorentzVector*)(eventLeptons.p4->At(idLep[0])) );
          nPassCuts[4]++;
          deltaPhiDileptonMet = TMath::Abs(dilep.DeltaPhi(*((TLorentzVector*)(*eventMet.p4)[0])));
          mtW = TMath::Sqrt(2.0*dilep.Pt()*((TLorentzVector*)(*eventMet.p4)[0])->Pt()*(1.0 - cos(deltaPhiDileptonMet)));
	  if(mtW > 50) passSel = kTRUE;
	}
      }

      if(passSel == kTRUE) nPassCuts[5]++;
      if(passSel == kFALSE) continue;

      if(mtW < 0){
        deltaPhiDileptonMet = TMath::Abs(dilep.DeltaPhi(*((TLorentzVector*)(*eventMet.p4)[0])));
        mtW = TMath::Sqrt(2.0*dilep.Pt()*((TLorentzVector*)(*eventMet.p4)[0])->Pt()*(1.0 - cos(deltaPhiDileptonMet)));
      }

      double mcWeight = eventMonteCarlo.mcWeight;
      if(infilecatv[ifile] == 0) mcWeight = 1.0;
      double theLumi = lumi; if(infilecatv[ifile] == 0) theLumi = 1.0;
      double puWeight = nPUScaleFactor(fhDPU, (double)eventVertex.npv); if(infilecatv[ifile] == 0) puWeight = 1.0;

TFile* SUSYLooperHists::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L SUSYLooperHists.C
//      Root > SUSYLooperHists 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 NULL;
   Long64_t nentries = fChain->GetEntriesFast();
  
   TFile* myFile = new TFile(outputFileName,"recreate");
   // checks the Mtt mass ++ simple plot example ++
   TH1D* LepTwoZmass= new TH1D("LepTwoZmass",";approximate Z(#tau#tau) mass [GeV];",100,0,2000);
   LepTwoZmass->Sumw2();

   // baseline plots +++++++++++++++++++++
   // allow to fill cutflow to sychronize to others
   TH1D* CutFlow = new TH1D("CutFlow",";CutFlow [unweighted];",20,-0.5,19.5);
   CutFlow->Sumw2();

   TH1D* Sig = new TH1D("Sig",";3DSig;",25,-0.,50);
   Sig->Sumw2();
   TH1D* Dxy = new TH1D("Dxy",";Dxy;",50,-0.,0.2);
   Dxy->Sumw2();
   TH1D* Dz = new TH1D("Dz",";Dz;",50,-0.,0.2);
   Dz->Sumw2();
   TH2D* DxyDz = new TH2D("DxyDz",";3DSig;",100,-0.,0.2,100,-0.,0.2);
   DxyDz->Sumw2();


   TH1D* IP3D = new TH1D("IP3D",";IP3D;",500,-0.,.5);
   IP3D->Sumw2();
   TH1D* Sig0 = new TH1D("Sig0",";3DSig;",50,-0.,50);
   Sig0->Sumw2();
   TH1D* Pt = new TH1D("Pt",";Pt;",200,-0.,25);
   Pt->Sumw2();
   TH1D* Iso = new TH1D("Is",";Is;",200,-0.,5);
   Iso->Sumw2();
   TH1D* LooseBPt = new TH1D("LooseBPt",";Pt;",200,-0.,500);
   LooseBPt->Sumw2();



 // a scan should be filled without weight*puWeights before any cut to get the efficiency. The reason is that for a scan each points have different x-section, which are not accessable during the loop
   TH2D* scan = new TH2D("scan","scan",8,112.5,312.5,40,112.5,312.5);
   // fill after some cuts the scans to get efficiency, i.e. ->Divide(scan) after the loop
   TH2D* scanA = new TH2D("scanA","scan",8,112.5,312.5,40,112.5,312.5);
   TH2D* scanB = new TH2D("scanB","scan",8,112.5,312.5,40,112.5,312.5);
   TH2D* scanC = new TH2D("scanC","scan",8,112.5,312.5,40,112.5,312.5);

   // plot a mass
   TH1D* GenMll = new TH1D("GenMll",";GenMll;",200,-0,200);
   GenMll->Sumw2();

   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;

      if(jentry%100000==0) cout << "Event: "<<jentry <<endl;
      scan->Fill(GenSusyMStop,GenSusyMNeutralino);
      // select lepton and put it into TLorentzvector 


      TLorentzVector lep;
      TLorentzVector met;
      vector <TLorentzVector> lepVec;
      for(int i =0;i<nLepGood;i++)
	{
	  TLorentzVector aLep;
	  aLep.SetPtEtaPhiM(LepGood_pt[i],LepGood_eta[i],LepGood_phi[i],0.1);
	  lepVec.push_back(aLep);
	}
      if(nLepGood>0) lep.SetPtEtaPhiM(LepGood_pt[0],LepGood_eta[0],LepGood_phi[0],0.1);
      else lep.SetPtEtaPhiM(0,0,0,0.1);
      met.SetPtEtaPhiM(met_pt,0.,met_phi,0);
      float MT =  sqrt(lep.Pt()*2*met.Pt()*(1-cos(lep.DeltaPhi(met))));
      TLorentzVector jetSum;
//.........这里部分代码省略.........

//--------------------------------------------------------------------------------------------------
const TPhoton* recoverFsr(const Lepton              &lep,
                          const vector<Lepton>      &lepvec,
                          const pair<Lepton,Lepton> &zcand,
                          const TClonesArray	    *photonArr)
{
    vector<TLorentzVector> photonMoms;
    vector<const TPhoton*> photonObjs;
    TLorentzVector phovec;

    for(int i=0; i<photonArr->GetEntriesFast(); i++) {
        const TPhoton *photon = (TPhoton*)photonArr->At(i);

        phovec.SetPtEtaPhiM(0,0,0,0);
        if((photon->typeBits & kPFMuonPhoton) || (photon->typeBits & kPFPhoton)) {
            // case of muon FSR where photon object is not reconstructed, but photon energy
            // is included in the ECAL energy of the muon object
            phovec.SetPtEtaPhiM(photon->pfPt, photon->pfEta, photon->pfPhi, 0);

        } else {
            continue;
        }

        if(phovec.Pt()        <= 2)   continue;
        if(fabs(phovec.Eta()) >= 2.4) continue;

        double dR = phovec.DeltaR(lep.p4);

        // veto if close to an electron supercluster
        bool flagEleSC = false;
        for(unsigned int j=0; j<lepvec.size(); j++) {
            if(abs(lepvec[j].pdgId)!=ELECTRON_PDGID) continue;

            double dPhi  = fabs(phovec.DeltaPhi(lepvec[j].p4));
            double dEta  = fabs(phovec.Eta() - lepvec[j].p4.Eta());
            if((dPhi<2. && dEta<0.05) || sqrt(dPhi*dPhi + dEta*dEta)<0.15) {
                flagEleSC = true;
                break;
            }
        }
        if(flagEleSC) continue;

        // check input lepton is the closest lepton to this photon
        bool found_closer_lepton = false;
        for(unsigned int j=0; j<lepvec.size(); j++) {
            if(lep.baconObj == lepvec[j].baconObj) continue;

            double tmp_dR = phovec.DeltaR(lepvec[j].p4);
            if(tmp_dR < dR) {
                found_closer_lepton = true;
                break;
            }
        }
        if(found_closer_lepton) continue;

        // Z mass OK?
        double oldMass = (zcand.first.p4 + zcand.second.p4).M();
        double newMass = (zcand.first.p4 + zcand.second.p4 + phovec).M();
        if( newMass <= 4.   ||
                newMass >= 100. ||
                fabs(newMass - Z_MASS) >= fabs(oldMass - Z_MASS) )
            continue;

        // "keep all photons close to one of the 4L leptons..."
        bool use = false;
        if(dR < 0.07) {
            use = true;
        } else if(dR<0.5 && phovec.Pt()>4 && photon->isoForFsr03<1.0*phovec.Pt()) {  // "need tighter cuts for other photons..."
            use = true;
        }
        if(use) {
            photonMoms.push_back(phovec);
            photonObjs.push_back(photon);
        }
    }

    // choose the best one
    double maxPt  =  0, minDR  = 999;
    int    iMaxPt = -1, iMinDR = -1;
    for(unsigned int i=0; i<photonMoms.size(); i++) {
        if(photonMoms[i].Pt() > maxPt) {
            maxPt  = photonMoms[i].Pt();
            iMaxPt = i;
        }
        double tmp_dR = photonMoms[i].DeltaR(lep.p4);
        if(tmp_dR < minDR) {
            minDR  = tmp_dR;
            iMinDR = i;
        }
    }

    if(maxPt > 4) {
        return photonObjs[iMaxPt];
    } else if(minDR < 999) {
        return photonObjs[iMinDR];
    } else {
        return 0;
    }
}

//.........这里部分代码省略.........
					{
						Pt_genneutrino = GenParticlePt -> at (ip);
						Phi_genneutrino = GenParticlePhi-> at (ip);
						Eta_genneutrino = GenParticleEta-> at (ip);
					}
				}

			}

			Pt_genMET = GenMETTrue->at(0);
			Phi_genMET = GenMETPhiTrue->at(0);
			genWTM = 0.0;
			genWTM =  TMass(Pt_genmuon1,Pt_genMET, fabs(Phi_genmuon1 - Phi_genMET) );

			// Set the recoil variables

			U1_Z_gen = 990.0;
			U2_Z_gen = 990.0;
			U1_W_gen = 990.0;
			U2_W_gen = 990.0;

			TLorentzVector v_GenMuon1, v_GenMuon2, v_GenMet;

			v_GenMuon1.SetPtEtaPhiM(Pt_genmuon1,Eta_genmuon1,Phi_genmuon1,0);
			v_GenMuon2.SetPtEtaPhiM(Pt_genmuon2,Eta_genmuon2,Phi_genmuon2,0);
			v_GenNu.SetPtEtaPhiM( Pt_genneutrino ,Eta_genneutrino,Phi_genneutrino ,0);
			v_GenMet.SetPtEtaPhiM ( Pt_genMET, 0, Phi_genMET,0 );

			Pt_Z_gen = (v_GenMuon1 + v_GenMuon2).Pt();
			Phi_Z_gen = (v_GenMuon1 + v_GenMuon2).Phi();

			TLorentzVector UZ_gen = -(v_GenMet + v_GenMuon1 + v_GenMuon2);
			TLorentzVector BZ_gen = (v_GenMuon1+v_GenMuon2 );
			U1_Z_gen = (UZ_gen.Pt()) * (cos(UZ_gen.DeltaPhi(BZ_gen))) ;
			U2_Z_gen = (UZ_gen.Pt()) * (sin(UZ_gen.DeltaPhi(BZ_gen))) ;

			TLorentzVector pfMETtest;
			pfMETtest.SetPtEtaPhiM(PFMET->at(0),0.0,PFMETPhi->at(0),0.0);
			muon1test.SetPtEtaPhiM(MuonPt->at(v_idx_muon_final[0]),MuonEta->at(v_idx_muon_final[0]),MuonPhi->at(v_idx_muon_final[0]),0.0);

			Pt_W_gen = (muon1test + v_GenNu).Pt();
			Phi_W_gen = (muon1test + v_GenNu).Phi();

			UW_gen = -(pfMETtest + muon1test );
			BW_gen = (muon1test +v_GenNu);
			U1_W_gen = (UW_gen.Pt()) * (cos(UW_gen.DeltaPhi(BW_gen))) ;
			U2_W_gen = (UW_gen.Pt()) * (sin(UW_gen.DeltaPhi(BW_gen))) ;

			if (IsW && DoRecoilCorrections)
			{

				float U1Phi = - BW_gen.Phi();
				float U2Phi = BW_gen.Phi() + piover2;

				if ((BW_gen.DeltaPhi(UW_gen)) < 0)   U2Phi = BW_gen.Phi() - piover2;

				float U1Prime = F_U1Prime(Pt_W_gen);
				float U2Prime = F_U2Prime(Pt_W_gen);

				TLorentzVector V_UPrime, V_U1Prime, V_U2Prime, V_MetPrime;

				V_U1Prime.SetPtEtaPhiM(U1Prime,0,U1Phi,0);
				V_U2Prime.SetPtEtaPhiM(U2Prime,0,U2Phi,0);
				V_UPrime = V_U1Prime + V_U2Prime;
				V_MetPrime = -(v_GenMuon1+ V_UPrime);

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.001;
  const double dCutEtaMax  = 1.6;
  const double dJetEtaMax  = 1.;
  const double dJetAreaRef = TMath::Pi() * dJetR * dJetR;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::E_scheme, fastjet::Best);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
  fastjet::JetDefinition subjDef(fastjet::kt_algorithm, dSjeR, fastjet::E_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;

  const Double_t dCut = TMath::TwoPi() / 3.;
  const double dMass  = TDatabasePDG::Instance()->GetParticle(211)->Mass();
//=============================================================================

  enum { kWgt, kLje, kLjr, kLtk, kLtr, kJet, kAje, kMje, k1sz, k1sA, k1sm, k1sr, k2sz, k2sA, k2sm, k2sr, kDsm, kDsr, kVar };

  TFile *file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  TNtuple *nt = new TNtuple("nt", "", "fWgt:fLje:fLjr:fLtk:fLtr:fJet:fAje:fMje:f1sj:f1sA:f1sm:f1sr:f2sj:f2sA:f2sm:f2sr:fDsm:fDsr");
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);

    double dLtk = -1.;
    TLorentzVector vPar, vLtk;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      double dEta = (*p)->momentum().eta(); if (TMath::Abs(dEta)>dCutEtaMax) continue;

      double dTrk = (*p)->momentum().perp();
      double dPhi = (*p)->momentum().phi();
      vPar.SetPtEtaPhiM(dTrk, dEta, dPhi, dMass);
      fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.E()));
      if (dTrk>dLtk) { dLtk = dTrk; vLtk.SetPtEtaPhiM(dTrk, dEta, dPhi, dMass); }
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
    std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
//=============================================================================

    if (selectJets.size()>0) {
      std::vector<fastjet::PseudoJet> sortedJets = fastjet::sorted_by_pt(selectJets);
      TLorentzVector vLje; vLje.SetPtEtaPhiM(sortedJets[0].pt(), sortedJets[0].eta(), sortedJets[0].phi(), sortedJets[0].m());

      TLorentzVector vJet;
      int kJrl = -1; double dJrl = -1.;
      int kTrl = -1; double dTrl = -1.;
      for (int j=0; j<sortedJets.size(); j++) {
        double dJet = sortedJets[j].pt();
        sortedJets[j].set_user_index(-1);
        vJet.SetPtEtaPhiM(dJet, sortedJets[j].eta(), sortedJets[j].phi(), sortedJets[j].m());
        if (TMath::Abs(vJet.DeltaPhi(vLje))>dCut) { sortedJets[j].set_user_index(1); if (dJet>dJrl) { dJrl = dJet; kJrl = j; } }
        if (TMath::Abs(vJet.DeltaPhi(vLtk))>dCut) { sortedJets[j].set_user_index(2); if (dJet>dTrl) { dTrl = dJet; kTrl = j; } }
      }
//=============================================================================

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

void UEAnalysisUE::ueAnalysisRECO(float weight,std::string tkpt,float etaRegion,float ptThreshold, TClonesArray* Track, TClonesArray* TracksJet, TFile* f, std::string hltBit)
{
  f->cd( hltBit.c_str() );

  // find leading jet in visible phase space
  TLorentzVector* leadingJet;
  Float_t PTLeadingTJ = -10;
  for(int j=0;j<TracksJet->GetSize();++j){
    TLorentzVector *v = (TLorentzVector*)TracksJet->At(j);
    if(fabs(v->Eta())<etaRegion){
      leadingJet = v;
      PTLeadingTJ= v->Pt();
      break;
    }
  }

  // save <pT> vs track multiplicity
  int numTracks( 0 );
  double trackPtSum( 0. );
  double averageTrackPt( 0. );
  for(int i=0;i<Track->GetSize();++i)
    {
      TLorentzVector *v = (TLorentzVector*)Track->At(i);
      if(v->Pt()>ptThreshold) 
	{
	  ++numTracks;
	  trackPtSum += v->Pt();
	}
    }
  if ( numTracks > 0 ) 
    {
      averageTrackPt = trackPtSum/numTracks;
      //std::cout << "[RECO] N(tracks)=" << numTracks << ", <pT>(tracks)=" << averageTrackPt << std::endl;
      h2d_averageTrackPt_vs_numTracks->Fill( numTracks, averageTrackPt, weight );
    }

  // catch events where no charged jet is found in the central region
  if ( PTLeadingTJ == -10 ) return;

  Float_t  PTLeadingCJ = cc->calibrationPt(PTLeadingTJ,tkpt)*PTLeadingTJ;

  for(int i=0;i<Track->GetSize();++i){
    TLorentzVector *v = (TLorentzVector*)Track->At(i);

    if(v->Pt()>ptThreshold){

      fHistPtDistRECO->Fill(v->Pt(),weight);
      fHistEtaDistRECO->Fill(v->Eta(),weight);
      fHistPhiDistRECO->Fill(v->Phi(),weight);
      temp3RECO->Fill(fabs(v->Eta()));
      temp4RECO->Fill(fabs(v->Pt()));
    }

    if(fabs(v->Eta())<etaRegion&&v->Pt()>=ptThreshold){
      
      //if (hltBit=="HLTMinBias") std::cout << "Track: pT=" << v->Pt() << ", eta=" << v->Eta() << ", phi=" << v->Phi() << std::endl;

      // use ROOT method to calculate dphi                                                                                    
      // convert dphi from radiants to degrees                                                                                
      Float_t conv = 180/piG;
      Float_t Dphi_reco = conv * leadingJet->DeltaPhi(*v);
      
      temp1RECO->Fill(Dphi_reco);
      temp2RECO->Fill(Dphi_reco,v->Pt());
    }
  }

  
  for(int i=0;i<100;i++){
    pdN_vs_etaRECO->Fill((i*0.05)+0.025,temp3RECO->GetBinContent(i+1)/0.1,weight);
  }
  for(int i=0;i<1000;i++){
    pdN_vs_ptRECO->Fill((i*0.1)+0.05,temp4RECO->GetBinContent(i+1)/0.1,weight);
  }
  
  temp3RECO->Reset();
  temp4RECO->Reset();
  
  Float_t transN1=0;
  Float_t transN2=0;
  Float_t transP1=0;
  Float_t transP2=0;
  Float_t towardN=0;
  Float_t towardP=0;
  Float_t awayN=0;
  Float_t awayP=0;
  
  for(int i=0;i<100;i++){
    if(i<=14){

      //std::cout << "[RECO] Away (" << i << "): dN=" << temp1RECO->GetBinContent(i+1) << ", dpT=" << temp2RECO->GetBinContent(i+1) << std::endl;

      awayN += temp1RECO->GetBinContent(i+1);
      awayP += temp2RECO->GetBinContent(i+1);
    }
    if(i>14 && i<33 ){

      //std::cout << "[RECO] Trans1 (" << i << "): dN=" << temp1RECO->GetBinContent(i+1) << ", dpT=" << temp2RECO->GetBinContent(i+1) << std::endl;

      transN1 += temp1RECO->GetBinContent(i+1);
//.........这里部分代码省略.........

void UEAnalysisUE::ueAnalysisMC(float weight,std::string tkpt,float etaRegion, float ptThreshold, 
				TClonesArray* MonteCarlo, TClonesArray* ChargedJet, TFile* f, std::string hltBit)
{
  f->cd( hltBit.c_str() );
  //  std::cout << "UEAnalysisUE::ueAnalysisMC(...), HLT " << hltBit << std::endl;

  TLorentzVector* leadingJet;
  Float_t PTLeadingCJ = -10;
  for(int j=0;j<ChargedJet->GetSize();++j){
    TLorentzVector *v = (TLorentzVector*)ChargedJet->At(j);
    if(fabs(v->Eta())<etaRegion){
      leadingJet = v;
      PTLeadingCJ= v->Pt();
      break;
    }
  }

  // save <pT> vs particle multiplicity
  int numParticles( 0 );
  double particlePtSum( 0. );
  double averageParticlePt( 0. );
  for(int i=0;i<MonteCarlo->GetSize();++i)
    {
      TLorentzVector *v = (TLorentzVector*)MonteCarlo->At(i);
      if(v->Pt()>ptThreshold)
        {
          ++numParticles;
          particlePtSum += v->Pt();
        }
    }
  if ( numParticles > 0 ) 
    {
      averageParticlePt = particlePtSum/numParticles;
      //  std::cout << "[MC] N(chg. part's)=" << numParticles << ", <pT>(chg. part's)=" << averageParticlePt << std::endl;
      h2d_averageParticlePt_vs_numParticles->Fill( numParticles, averageParticlePt, weight );
    }


  //std::cout << "PTLeadingCJ " << PTLeadingCJ << std::endl;

  if ( PTLeadingCJ == -10. )
    {
      //std::cout << "return" << std::endl;
      return;
    }

  h_pTChgGenJet->Fill( PTLeadingCJ, weight );
  //std::cout << "for(int i=0;i<MonteCarlo->GetSize();i++){" << std::endl;

  //   if ( hltBit == "All" )
  //     {
  //       std::cout << "[UEAnalysisUE]" << std::endl;
  //     }
  for(int i=0;i<MonteCarlo->GetSize();i++){
    TLorentzVector *v = (TLorentzVector*)MonteCarlo->At(i);    

    if(v->Pt()>=ptThreshold){
      fHistPtDistMC->Fill(v->Pt(),weight);
      fHistEtaDistMC->Fill(v->Eta(),weight);
      fHistPhiDistMC->Fill(v->Phi(),weight);
      temp3MC->Fill(fabs(v->Eta()));
      temp4MC->Fill(fabs(v->Pt()));
    }

    if(fabs(v->Eta())<etaRegion && v->Pt()>=ptThreshold){

      //if (hltBit=="HLTMinBias") std::cout << "Particle: pT=" << v->Pt() << ", eta=" << v->Eta() << ", phi=" << v->Phi() << std::endl;

      Float_t conv = 180/piG;
      Float_t Dphi_mc = conv * leadingJet->DeltaPhi(*v);

      //       if ( hltBit == "All" )
      //        	{
      //        	  std::cout << "(" << i << ") ";
      //        	  std::cout << "pT, eta, phi, dphi ";
      //        	  std::cout << v->Pt() << ", ";
      //        	  std::cout << v->Eta() << ", ";
      //        	  std::cout << v->Phi() << ", ";
      //        	  std::cout << Dphi_mc << std::endl;
      //        	}
      
      temp1MC->Fill(Dphi_mc);
      temp2MC->Fill(Dphi_mc,v->Pt());
    }
  }
  
  //std::cout << "for(int i=0;i<100;i++){" << std::endl;

  for(int i=0;i<100;i++){
    pdN_vs_etaMC->Fill((i*0.05)+0.025,temp3MC->GetBinContent(i+1)/0.1,weight);
  }
  for(int i=0;i<1000;i++){
    pdN_vs_ptMC->Fill((i*0.1)+0.05,temp4MC->GetBinContent(i+1)/0.1,weight);
  }

  temp3MC->Reset();
  temp4MC->Reset();
    
  Float_t transN1=0;
  Float_t transN2=0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if( abs(RecoSignalType())!=1010  ) return false;//this is one mu
  }
  if( Process == "FakeMuonsV2" ){
    ////double LeadingMuPt_  = LeadingMuPt();
    ////if( LeadingMuPt_<35          ) return false;//LJ filter removed
    if( TotalJet<1                   ) return false;
    //double LeadingJetPt_ = LeadingJetPt();
    //if( LeadingJetPt_<50         ) return false;
    //if( !isZToMuMu()             ) return false;
    //if( ST()>350                 ) return false;
    if( METlepMT("Mu")>20            ) return false;
    if( METcorr("Pt")>20      ) return false;
    if( abs(RecoSignalType())!=1010  ) return false;//this is one mu
  }
  if( Process == "FakeMuonsV3" ){
    if( METlepMT("Mu")>10            ) return false;
    if( METcorr("Pt")>10      ) return false;
    if( TotalJet<1                   ) return false;//
    if( abs(RecoSignalType())!=1010  ) return false;//this is one mu
    if( !isAllMuonsHLT()             ) return false;//mu has to match the trigger
    //
    //--------- Veto on 2nd Global Muon ---------------------
    int nGlobalMuons=0;
    for(unsigned int iMuR=0;  iMuR<MuonPt->size();     iMuR++){
      if( muPtcorr(iMuR)>MuonPtCut && fabs(MuonEta->at(iMuR))<2.1 && MuonIsGlobal->at(iMuR)==1 ) nGlobalMuons++;
    }
    if( nGlobalMuons>1               ) return false;
    //
    //--------- Check for back-to-back Jet-Muon -------------
    bool isBackToBackJet=false;
    TLorentzVector Mu;
    TLorentzVector Jet;
    for(unsigned int iMuR=0;  iMuR<MuonPt->size();     iMuR++){
      if( !muRisoCheck(iMuR)     ) continue;
      Mu.SetPtEtaPhiM(muPtcorr(iMuR), MuonEta->at(iMuR), MuonPhi->at(iMuR), 0);
      for(unsigned int iJetR=0;  iJetR<PFJetPt->size();     iJetR++){
	if( !jetRTightCheck(iJetR) ) continue;
	Jet.SetPtEtaPhiM( jetPtcorr(iJetR), PFJetEta->at(iJetR), PFJetPhi->at(iJetR), 0 );
	if( (fabs(fabs(Mu.DeltaPhi(Jet))-TMath::Pi())/TMath::Pi())<0.1 )  isBackToBackJet=true;
      }
    }
    //--------- Check for back-to-back Jet-Muon ------------- 
    if( !isBackToBackJet             ) return false;
  }
  if( Process == "ControlRegion1" ){
    //low in signal contamination, to check tau fakes loose->tight ratio
    if( TotalBJet<1                    ) return false;
    if( TotalJet<2                     ) return false;
    ///ST cut reverted
    ///remove MaxMuTauInvMass CUT:  /*if( MaxMuTauInvMass()<100          ) return false;*/
    if( LeadingTauPt()<50              ) return false;
    if( RecoSignalType()!=2110         ) return false;// SS selection
    if( isZToMuMu()                    ) return false;// required to exclude events in MuTrig Calculation  
    //
    if( ST()>400                       ) return false;// Revert the ST cut
  }
  if( Process == "TTBar" ){
    double ZToMuMuCentral = 90;
    double ZToMuTauCentral = 75;
    if( MaxDiLepInvMass()<65       ) return false;
    if( TotalN<5                 ) return false;
    if( ST()<250           ) return false;
    if( RecoSignalType()>0 ) return false;
    if( TotalJet<3               ) return false;
    if( fabs(MuTauInZpeak("ZToMuMu")-ZToMuMuCentral)<10  ) return false; //exclue 80-100
    if( fabs(MuTauInZpeak("ZToMuTau")-ZToMuTauCentral)<15 ) return false; //exclude 60-90
  }
  if( Process == "WJets" ){
    double ZToMuMuCentral = 90;
    double ZToMuTauCentral = 75;
    if( MaxDiLepInvMass()<65       ) return false;
    if( TotalMu!=1          ) return false;
    if( TotalJet!=0         ) return false;
    if( ST()>250      ) return false;
    if( fabs(MuTauInZpeak("ZToMuMu")-ZToMuMuCentral)<10  ) return false; //exclue 80-100
    if( fabs(MuTauInZpeak("ZToMuTau")-ZToMuTauCentral)<15 ) return false; //exclude 60-90
  }
  if( Process == "ZJets" ){
    if( MaxDiLepInvMass()<65 ) return false;
    if( METlepMT("Mu")>50  ) return false;
    if( TotalMu < 2                ) return false;
    double ZToMuMuCentral  = 90;
    double ZToMuTauCentral = 75;
    bool   dilepZpeak_     = false;
    if( fabs(MuTauInZpeak("ZToMuMu")-ZToMuMuCentral)<10  ) dilepZpeak_=true; //accept 80-100
    if( fabs(MuTauInZpeak("ZToMuTau")-ZToMuTauCentral)<15 ) dilepZpeak_=true; //accept 60-90
    if( !dilepZpeak_               ) return false;
  }
  if( Process == "LQ3M400" ){
    if( MaxDiLepInvMass()<65       ) return false;
    if( LeadMuTauDeltaR>4          ) return false;
    if( TotalN<4                   ) return false;
    if( METcorr("Pt")<30    ) return false;
    if( ST()<400             ) return false;
    if( RecoSignalType()<0   ) return false;
  }
  //
  return true;
  //
}