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

void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
{
	int n_smc  = 0;
	int n_strk = 0;
	int n_both = 0;
	double d = 0.00001;
	
	for (int i=0;i<l.GetLength()-1;++i)
	{
		for (int j=i+1;j<l.GetLength();++j)
		{
			TLorentzVector dl = l[i]->P4() - l[j]->P4();
		
			bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth());
			bool chktrk = (fabs(dl.X())<d) && (fabs(dl.Y())<d) && (fabs(dl.Z())<d) && (fabs(dl.E())<d);
			if (chkmc) n_smc++;
			if (chktrk) n_strk++;
			if (chktrk && chkmc) n_both++;
		}	
	}
	n1 = n_strk;
	n2 = n_smc;
	n3 = n_both;
}

//.........这里部分代码省略.........
	std::cout << "Got input file(s)." << std::endl;

	//Setup access to the Recon tree
	int NPIDs(0);  // This variable counts the number of particles per event
	int NVtxFGD1(0), NVtxFGD2(0);
        // Declare a TClonesArray to hold objects of type TGlobalPID
 	TClonesArray *globalPIDs = new TClonesArray("ND::TGlobalReconModule::TGlobalPID",50);
	TClonesArray *VtxFGD1 = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex",50);    
	TClonesArray *VtxFGD2 = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex",50);    
    // Associate the right branch in the TTree to the right local variable
	gRecon->SetBranchAddress("NPIDs",&NPIDs);
    gRecon->SetBranchAddress("PIDs",&globalPIDs);
	gGenVtx->SetBranchAddress("VtxFGD1", &VtxFGD1);
	gGenVtx->SetBranchAddress("NVtxFGD1", &NVtxFGD1);
	gGenVtx->SetBranchAddress("VtxFGD2", &VtxFGD2);
	gGenVtx->SetBranchAddress("NVtxFGD2", &NVtxFGD2);
	//check that truthdir and recon have the same number of entries
	if(gRecon->GetEntries() != gGenVtx->GetEntries()) 
		cout<<"not equal entries, probably wrong"<<endl;
	// Loop over the entries in the TChain.

	//========================================================
	//			Declare Graphs n stuff here
	//========================================================

	//adding tclones arrays for use with detectors
	Int_t NTPCs;
	TClonesArray *TPC;

	Int_t NFDGs;
	TClonesArray *FDG;

	Int_t NECALs;
	TClonesArray *ECAL;

	Int_t NPODs;
	TClonesArray *POD;

	Int_t NSMRDs;
	TClonesArray *SMRD;

	//adding a 2d graph general purpose, change titles each time!
	TH1D *graph1 = new TH1D("graph1","Momenta of TPC PIDs from FDGs", 100, -50.0 , 50.0);
	
	Int_t ninteract(0), nfgd(0),nfgdqes(0),fgdqesneu(0);
	//========================================================
	//	end		Declare Graphs n stuff here
	//========================================================

	// Loop over the entries in the TChain. (only 1/1000 of whole entries atm)
	for(unsigned int i = 0; i < gRecon->GetEntries()/10; ++i) {
		if((i+1)%10000 == 0) std::cout << "Processing event: " << (i+1) << std::endl;
		//display status every 10,000 th entry

	// Get an entry for the Recon tree
		gRecon->GetEntry(i);
		gGenVtx->GetEntry(i);
		ND::TGlobalReconModule::TGlobalPID *gTrack = NULL;

		//added new loop for truth vertex
		gGenVtx->GetEntry(i);
		
		for (int j=0; j<NPIDs; j++) {
			// Get a specific track from the TClonesArray
			gTrack = (ND::TGlobalReconModule::TGlobalPID*)globalPIDs->At(j);
			//get truevertex (in example, also gets trueparticle, can add in later)
			ND::TTrueVertex vtx = gTrack->TrueParticle.Vertex;
			//get position lorrentz vector
			TLorentzVector vec = vtx.Position;
			ninteract++;
			if(ABS(vec.X())<832.2 && ABS(vec.Y()-55)<832.2 && ((vec.Z()>123.45&&vec.Z()<446.95)||(vec.Z()>1481.45&&vec.Z()<1807.95))){	//is it in one of the FGDs?
				nfgd++;
				if(vtx.ReactionCode.find("Weak[NC],QES;",0)!=-1){
					unsigned long det;
					det = gTrack->Detectors;
					string detstr;
					stringstream stream;
					stream << det;
					detstr = stream.str();
					if(detstr.find("1",0)||detstr.find("2",0)||detstr.find("3",0)){
						graph1->Fill((Double_t)gTrack->FrontMomentum);	
					} 	
					nfgdqes++;
				}	
			}
			TClonesArray *TPCObjects = new TClonesArray("ND::TGlobalReconModule::TTPCObject",gTrack->NTPCs);
			ND::TGlobalReconModule::TObject *tpcTrack = NULL;
			for ( int k = 0 ; k < gTrack->NTPCs; k++) {
				tpcTrack = (ND::TGlobalReconModule::TObject*) TPCObjects->At(k);
				//now we can access  variables through tpcTrack->PullEle for example
			}
		}

	} // End loop over events

//plotting bits at the end :D
	graph1->Draw();
	App->Run();
	return 0;
}

bool leptonic_fitter_algebraic::fit( const TLorentzVector& B, const TH1& BITF, const TF1& Beff, 
				     const TLorentzVector& lep, 
				     double MEX, double MEY, const TF1& dnuPDF )
{
  if( _dbg > 19 ) cout<<"DBG20 Entered leptonic_fitter_algebraic::fit with B mass: "<<B.M()<<", l_m:"<<lep.M()<<", MET: "<<MEX<<" "<<MEY<<endl;
  if( B.M() <= 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a b-jet with an illegal (non-positive) mass!"); 
  if( lep.M() < 0 ) throw std::runtime_error( "leptonic_fitter_algebraic was given a lepton with an illegal (negative) mass!"); 
  _converged = _swapped = false;
  _obsB = B;
  _obsL = lep;

  _BITF = &BITF;
  _Beff = &Beff;
  _dnuPDF = dnuPDF;

  _b_m2 = B.M2();

  double lep_b_angle = lep.Angle( B.Vect() );
  double cos_lep_b = TMath::Cos( lep_b_angle );
  double sin_lep_b = TMath::Sin( lep_b_angle );
  double b_p = B.P();
  double b_e = B.E();
  _denom = b_e - cos_lep_b * b_p;
  
  _lep_p = lep.P();
  _x0 = - _W_m2 / ( 2 * _lep_p );
  _y1 = - sin_lep_b * _x0 * b_p / _denom;
  _x1_0 = _x0 * b_e / _denom  -  _y1*_y1 / _x0;
  _Z2_0 = _x0*_x0 - _W_m2 - _y1*_y1;
  if( _dbg > 219 ) cout<<"DBG220 lfa updated lepton with: "<<lv2str( lep )<<" -> x0:"<<_x0<<", y1: "<<_y1<<", x1_0: "<<_x1_0<<", Z2_0: "<<_Z2_0<<endl;

  static double bnums[3];
  bnums[0] = B.X();
  bnums[1] = B.Y();
  bnums[2] = B.Z();
  TMatrixD bXYZ( 3, 1, bnums );
  _R_T = rotation( 2, lep.Phi() ); // R_z^T
  _R_T *= rotation( 1, lep.Theta() - 0.5*TMath::Pi() ); // R_z^T R_y^T
  TMatrixD rotation_vect( _R_T, TMatrixD::kTransposeMult, bXYZ ); // R_y R_z
  double* rotation_array = rotation_vect.GetMatrixArray();
  double phi_x = - TMath::ATan2( rotation_array[2], rotation_array[1] );
  if( _dbg > 99 ) cout<<"DBG100 lfa x rotation vector is:"<<rotation_array[0]<<" "<<rotation_array[1]<<" "<<rotation_array[2]<<" -> phi_x:"<<phi_x<<endl;
  _R_T *= rotation( 0, - phi_x ); // R_z^T R_y^T R_x^T

  // set up _Nu's non-zero elements so that \vec{nu} = Nu \vec{t} for any \vec{t} (since only t's 3nd component is used, and its always 1).
  _Nu[0][2] = MEX;
  _Nu[1][2] = MEY;

  double iVarMET = TMath::Power( TMath::Max( 1., dnuPDF.GetHistogram()->GetRMS() ), -2 );
  _invFlatVar[0][0] = _invFlatVar[1][1] = iVarMET; // set up the chi^2 distance with the right order of magnitude (generalizes to rotated covariance matrix)
  if( _dbg > 209 ) cout<<"DBG210 lfa "<<dnuPDF.GetName()<<" --> iVarMET:"<<iVarMET<<endl;

  // (re)define fit parameter, so all fits start off on an equal footing
  _mini->SetPrintLevel( _minimizer_print_level );
  _mini->Clear();
  _mini->SetFunction( _functor );
  leptonic_fitter_algebraic_object = this; // set the function in the functor pointing back to this object. Doubtfull that all this redirection is needed...
  _mini->SetTolerance( _tolerance );
  bool OK = _mini->SetLimitedVariable( 0, "sB", 1.0, 0.4, 0.1, 6.0 );
  //bool OK = _mini->SetVariable( 0, "sB", 1.0, 0.4 );
  if( ! OK ) {cerr<<"minimizer (@lfa) failed to SetVariable."<<endl; return false;}

  // define 1 sigma in terms of the function
  _mini->SetErrorDef( 0.5 ); // since this is a likelihood fit

  // do the minimization
  OK = _mini->Minimize(); 
  if( _dbg > 19 && ( ! OK || _dbg > 59 ) ) cout<<"DBG INFO: initial fit @lfa returned OK: "<<OK<<", has status: "<<_mini->Status()<<endl;

  _converged = OK; // use status somehow? depends on fitter?

  // read parameters
  const double *xs = _mini->X();
  for( int ip = 0; ip < 1; ++ip ) _params[ ip ] = xs[ ip ];

  // return all intermediate results to the minimum, in particular, the discriminant
  calc_MLL( _params, true );
  TMatrixD nu_vec( _Emat, TMatrixD::kMult, _tvec );
  update_nu_and_decay_chain( nu_vec );
  if( _dbg > 203 ) cout<<"DBG204 lfa finalized _genN: "<<lv2str(_genN)<<", _W: "<<lv2str(_W)<<", & _t: "<<lv2str(_T)<<endl;

  _MLL = _mini->MinValue();
  return true;
} 

//.........这里部分代码省略.........
		CCQEScount[NCuts],
		CCREScount[NCuts],
		NCREScount[NCuts],
		DIScount[NCuts],
		otherCount[NCuts];
	memset(NCEScount, 0, NCuts*sizeof(int));
	memset(CCQEScount, 0, NCuts*sizeof(int));
	memset(CCREScount, 0, NCuts*sizeof(int));
	memset(NCREScount, 0, NCuts*sizeof(int));
	memset(DIScount, 0, NCuts*sizeof(int));
	memset(otherCount, 0, NCuts*sizeof(int));
	memset(totCount, 0, NCuts*sizeof(int));
	
cout<< "haven't segfaulted yet" << endl;
	int fgd1count(0),fgd2count(0), fgd1nces(0), fgd2nces(0);
	for(unsigned int i = 0; i < NTOT; ++i) {
		if((i+1)%1000 == 0) std::cout << "Processing event: " << (i+1) << std::endl;
		//display status every 1,000 th entry
		// Get an entry for the tree
		tree->GetEntry(i);
		int keep(1), j(0), isNCES = TrueParticle->Vertex.ReactionCode.find("Weak[NC],QES;",0)!=-1;
		if (inFGD1(FrontPosition)){
			fgd1count++;
			if (isNCES)
				fgd1nces++;
		}
		else if (inFGD2(FrontPosition)){
			fgd2count++;
			if(isNCES)
				fgd2nces++;
		}
			
		initialNCES += isNCES;

		//apply cuts here
		//cout<<TrueParticle->Vertex.ReactionCode<<endl;
		//APPLY CUTS HERE!!!
		//This program cuts from our tree, not original data.
		//Currently doesn't output anything, or cut anything
		// but cuts are defined in cuts.C
		for( ; j<NCuts; j++){
			keep = cuts[j]->apply();
			if(!keep)
				break;
			totCount[j]++;
			NCEScount[j] += isNCES;
			if (!isNCES){
				
				CCQEScount[j] += TrueParticle->Vertex.ReactionCode.find("Weak[CC],QES;",0)!=-1;
				CCREScount[j]  += TrueParticle->Vertex.ReactionCode.find("Weak[CC],RES;",0)!=-1;
				NCREScount[j]  += TrueParticle->Vertex.ReactionCode.find("Weak[NC],RES;",0)!=-1;
				DIScount[j]  += TrueParticle->Vertex.ReactionCode.find("DIS;",0)!=-1;	
			}
		}
		if(keep){	
			Double_t recon_mom = ((ND::TGlobalReconModule::TTPCObject*)TPC->At(0))->FrontMomentum;
			TLorentzVector v = TrueParticle->InitMom;
			Double_t truth_mom = sqrt(v.X()*v.X()+v.Y()*v.Y()+v.Z()*v.Z());
		}
//			if (!isNCES&&keep){
//				string type;
//				if(TrueParticle->Vertex.ReactionCode.find("Weak[CC],QES;",0)!=-1){
//					CCQEScount[j]++;
//					type = string("CCQES");
//				}
//				else if(TrueParticle->Vertex.ReactionCode.find("Weak[CC],RES;",0)!=-1){
//					CCREScount[j]++;
//					type = string("RES");
//				}
//				else if(TrueParticle->Vertex.ReactionCode.find("Weak[NC],RES;",0)!=-1){
//					NCREScount[j]++; 
//					type = string("RES");
//				}
//				else if(TrueParticle->Vertex.ReactionCode.find("DIS;",0)!=-1){	
//					DIScount[j]++;
//					type = string("DIS");
//				}
//				Outfile << type << "," << FOName->GetString().Data() << "," << EventID << endl;
//			}		
		//after cuts applied, keep will be = 1 if it is to be kept	
	} // End loop over events
	//printf("initially: eff = %6.5f, pur = %6.5f\n", (double)initialNCES/(double)NNCES, (double)initialNCES/(double)NTOT);
	Double_t mcPOT = 5.e17 * 1554;
	Double_t scale = 5.e21/mcPOT;
	cout << "fgd1: " << fgd1count*scale << "	fgd2: " << fgd2count*scale << endl;
	cout << "fgd1nces: " << fgd1nces*scale << "	fgd2nces: " << fgd2nces*scale << endl;
	cout<<"After FGD cuts: eff="<<(double)NNCES/(double)NNCES<<"  pur="<<(double)NNCES/(double)NTOT<<endl;
	cout<<"After in beam time cut: eff="<<(double)NNCESatb/(double)NNCES<<"  pur="<<(double)NNCESatb/(double)NTOTatb<<endl;
	for (int n(0); n < NCuts; n++){
		double eff = (double)NCEScount[n]/(double)NNCES, pur = (double)NCEScount[n]/(double)totCount[n];
		printf("%12s:     eff = %6.5f, pur = %6.5f, e*p = %6.5f\n",cuts[n]->name.c_str(),eff, pur, eff*pur);
	}
	cout << "Background components:" << endl;
	printf("%6.5f%% CCQES\n", (double)CCQEScount[NCuts-1]/((double)totCount[NCuts-1] - (double)NCEScount[NCuts-1]));
	printf("%6.5f%% NCRES\n", (double)NCREScount[NCuts-1]/((double)totCount[NCuts-1] - (double)NCEScount[NCuts-1]));
	printf("%6.5f%% CCRES\n", (double)CCREScount[NCuts-1]/((double)totCount[NCuts-1] - (double)NCEScount[NCuts-1]));
	printf("%6.5f%% DIS\n", (double)DIScount[NCuts-1]/((double)totCount[NCuts-1] - (double)NCEScount[NCuts-1]));
	treefile->Close();
	return 0;
}

// Get corrections from table
float qq2vvEWKcorrections::getqq2WWEWKCorr(
                       float ptl1 , float etal1 , float phil1 , float idl1 , // lepton from 1st W
                       float ptl2 , float etal2 , float phil2 , float idl2 , // lepton from 2nd W
                       float ptv1 , float etav1 , float phiv1 ,             // neutrino from 1st W
                       float ptv2 , float etav2 , float phiv2 ,             // neutrino from 2nd W
                       float x1   , float x2 ,                              // parton x-Bjorken
                       float id1  , float id2                               // parton PDG id's
//                       int   TypeFirst = 1 ,                                // ??????????????????????????
//                        float Energy = 8000.
                     ) {

 float Energy = 8000.;
 // Create lepton and neutrino vectors
 TLorentzVector l1;
 TLorentzVector l2;
 TLorentzVector v1;
 TLorentzVector v2;

 float me   = 0.001*0.510998928 ;
 float mmu  = 0.001*105.6583715 ;
 float mtau = 0.001*1776.82     ;

 float ml1(0) , ml2(0) ;

 if      ( idl1 == 11 ) ml1 = me   ;
 else if ( idl1 == 13 ) ml1 = mmu  ;
 else if ( idl1 == 15 ) ml1 = mtau ;
 if      ( idl2 == 11 ) ml2 = me   ;
 else if ( idl2 == 13 ) ml2 = mmu  ;
 else if ( idl2 == 15 ) ml2 = mtau ;

 l1.SetPtEtaPhiM(ptl1,etal1,phil1,ml1);
 l2.SetPtEtaPhiM(ptl2,etal2,phil2,ml2);
 v1.SetPtEtaPhiM(ptv1,etav1,phiv1,0.0);
 v2.SetPtEtaPhiM(ptv2,etav2,phiv2,0.0);

 TLorentzVector W1 = l1+v1; // W1
 TLorentzVector W2 = l2+v2; // W2
 TLorentzVector WW = W1+W2;

 //---- FIX
 //---- swap 1 <->2 for leptons and neutrinos, to get correctly association l-v:
 //---- use as test invariatn mass of di-leptons
 //----    needed because leptons are ordered by "pt" and not by "mother" in the input! 
 if ((W1.M() - 80.385) > 0.1 && (W2.M() - 80.385) > 0.1) {
  W1 = l1+v2; // W1
  W2 = l2+v1; // W2
  WW = W1+W2;
 }
 //---- end FIX


 float M_12   = 80.385 , M_22 = 80.385 ;

 TLorentzVector p1; p1.SetPxPyPzE(0.,0.,Energy*x1,Energy*x1);
 TLorentzVector p2; p2.SetPxPyPzE(0.,0.,-Energy*x2,Energy*x2);

 //S-HAT
 float s_hat = pow(WW.E(),2)-pow(WW.Px(),2)-pow(WW.Py(),2)-pow(WW.Pz(),2); // ScalarProd(WW) (p1+p2)^2 = (p3+p4)^2 ~ +2*p1*p2
  //T_HAT
  //float t_hat2 = TypeFirst == 1 ? p1.M()*p1.M() + M_12*M_12 - 2*( p1.E()*W1.E() - p1.Px()*W1.Px() - p1.Py()*W1.Py() - p1.Pz()*W1.Pz() ) :
  //                                p1.M()*p1.M() + M_22*M_22 - 2*( p1.E()*W2.E() - p1.Px()*W2.Px() - p1.Py()*W2.Py() - p1.Pz()*W2.Pz() )   ; //T_HAT LO
 float la1 = sqrt( pow(s_hat,2) );            //la = sqrt( pow(a,2)+pow(b,2)+pow(c,2)-2*(a*b+a*c+b*c) );
 float la2 = sqrt( pow(s_hat,2) + pow(M_12,2) + pow(M_22,2) - 2*(s_hat*M_12 + s_hat*M_22 + M_12*M_22) );
  //  Boost: boost ext. momenta in CM frame of W1,W2
 TLorentzVector W1_b = W1, W2_b = W2;
 TLorentzVector p1_b = p1, p2_b = p2;
 W1_b.Boost( -WW.BoostVector() );
 W2_b.Boost( -WW.BoostVector() );
 p1_b.Boost( -WW.BoostVector() );
 p2_b.Boost( -WW.BoostVector() );

 //  Uni-vector
 TLorentzVector ee1 = p1_b*(1./sqrt( pow(p1_b.X(),2)+pow(p1_b.Y(),2)+pow(p1_b.Z(),2) ));
 TLorentzVector ee2 = p2_b*(1./sqrt( pow(p2_b.X(),2)+pow(p2_b.Y(),2)+pow(p2_b.Z(),2) ));
 TLorentzVector z1  = W1_b*(1./sqrt( pow(W1_b.X(),2)+pow(W1_b.Y(),2)+pow(W1_b.Z(),2) ));
 TLorentzVector z2  = W2_b*(1./sqrt( pow(W2_b.X(),2)+pow(W2_b.Y(),2)+pow(W2_b.Z(),2) ));
 //  "effective" beam axis
 float abse = sqrt( pow(ee1.X()-ee2.X(),2) + pow(ee1.Y()-ee2.Y(),2) + pow(ee1.Z()-ee2.Z(),2) );
 TLorentzVector ee = (ee1-ee2) * (1. / abse);
  //  "effective" scattering angle
 float costh = ee.X()*z1.X()+ee.Y()*z1.Y()+ee.Z()*z1.Z();
  //  final T_HAT
 float t_hat= M_12 - (1./2.)*(s_hat+M_12-M_22) + (1/(2*s_hat))*la1*la2*costh; //Mz-1/2*s+1/(2*s)*la(s,0,0)*la(s,mZ,mZ)*costh  or: (p1-p3)^2 = (p2-p4)^2 ~ -2*p1*p3
  //Quark Type
 int quarkType = -1.;
 if(      fabs(id1)==2 && fabs(id2)==2 ) quarkType=0; //delta_uub
 else if( fabs(id1)==1 && fabs(id2)==1 ) quarkType=1; //delta_ddb
 else if( fabs(id1)==5 && fabs(id2)==5 ) quarkType=2; //delta_bbb
 else if( fabs(id1)==4 && fabs(id2)==4 ) quarkType=0; // cc as delta_buu
 else if( fabs(id1)==3 && fabs(id2)==3 ) quarkType=1; // ss as delta_bdd

 // Extracting corrections
 float EWK_w = 1. ;
//  std::cout << " quarkType = " << quarkType << " id1 = " << id1 << " id2 = " << id2 << std::endl;
 if( quarkType!=-1 ){
  float sqrt_s_hat = sqrt(s_hat);
  std::vector<std::vector<float> > EWK_w2_vec = findCorrection( sqrt_s_hat, t_hat );
  float EWK_w2 = 1. + EWK_w2_vec[0][quarkType];
//.........这里部分代码省略.........

	void BreitBoost(TObjArray &InputCandList, TObjArray &OutputCandList, Double_t PzInLepton, Double_t PzInHadron,TLorentzVector POutLepton,Int_t Scheme)
	{
	// Boost input particle into Breit-frame
	// PInLepton is momentum of incoming lepton in +z direction 
	// POutLepton is momentum of outgoning lepton  
	// PInHadron is momentum of incoming hadron in +z direction 
	// Scheme = 1 => boost to Breit-frame  
	//        = 2 => boost from Breit-frame
	// Loop. vars
	Int_t k4Vec;
	
	// Number of entries
	Int_t NumInputEntries;
	
	// Breit-vars
	Double_t ArbFactor, PRelLepton,CenterOfMass[4];
	Double_t ZHadrons[4], XZ[4];
	
	// boost to and from Breit-frame implemented
	if (Scheme != 1 && Scheme != 2)
	return;
	
	// Get number of particles stored in the list
	if (InputCandList.GetEntries()>0)
	NumInputEntries = InputCandList.GetEntries();  
	else {
	printf("Inputlist is empty: No boost possible\n");
	return;
	}
	/*   if(m_Debug > 1)  */
	/*     cout << "In Breitboost: NumInputEntries = "  */
	/* 	 << NumInputEntries << endl; */
	
	// Init
	TLorentzVector *Boosted4Vector = new TLorentzVector(0,0,0);
	for (k4Vec = 0;k4Vec < 4;k4Vec++)
	{
	CenterOfMass[k4Vec] = 0;
	ZHadrons[k4Vec] = 0;
	XZ[k4Vec] = 0;
	}
	
	// Find 4-momentum of Breit frame (times an arbitrary factor)
	// Breit-frame 4-vector given by "CenterOfMass"
	
	/*  if(m_Debug > 1)  */
	/*     cout << PzInHadron << ' ' << PzInLepton << ' ' */
	/* 	 << POutLepton.X() << ' ' */
	/* 	 << POutLepton.Y() << ' ' */
	/* 	 << POutLepton.Z() << ' ' */
	/* 	 << POutLepton.E() << endl; */
	
	// Calculate arbitrary factor 
	ArbFactor = fabs(PzInHadron)*(fabs(PzInLepton)-POutLepton.E())-PzInHadron*(PzInLepton-POutLepton.Z());
	/*  if(m_Debug > 1)  */
	/*     cout << "In Breitboost: ArbFactor = "  */
	/* 	 << ArbFactor << endl; */
	
	// Relative momentum squared
	PRelLepton = pow((fabs(PzInLepton)-POutLepton.E()),2)-pow((PzInLepton-POutLepton.Z()),2)-pow(POutLepton.Y(),2)-pow(POutLepton.X(),2);
	/*  if(m_Debug > 1)  */
	/*     cout << "In Breitboost: PRelLepton = "  */
	/* 	 << PRelLepton << endl; */
	
	// CM-calculation
	CenterOfMass[0] = ArbFactor * ( -POutLepton.X()); 
	CenterOfMass[1] = ArbFactor * ( -POutLepton.Y());
	CenterOfMass[2] = ArbFactor * (PzInLepton - POutLepton.Z())-PRelLepton*PzInHadron;
	CenterOfMass[3] = ArbFactor * (fabs(PzInLepton)-POutLepton.E())-PRelLepton*fabs(PzInHadron);
	/*  if(m_Debug > 1)  */
	/*     cout << "In Breitboost: CenterOfMass = "  */
	/* 	 << CenterOfMass[0] << ' '  */
	/* 	 << CenterOfMass[1] << ' '  */
	/* 	 << CenterOfMass[2] << ' '  */
	/* 	 << CenterOfMass[3] << ' '  */
	/* 	 << endl; */
	
	// Find rotation to put incoming hadron back on z-axis      
	ZHadrons[0] = 0.;
	ZHadrons[1] = 0.;
	ZHadrons[2] = PzInHadron;
	ZHadrons[3] = fabs(PzInHadron);
	/*   if(m_Debug > 1)  */
	/*     cout << "In Breitboost: ZHadrons = "  */
	/* 	 << ZHadrons[0] << ' '  */
	/* 	 << ZHadrons[1] << ' '  */
	/* 	 << ZHadrons[2] << ' '  */
	/* 	 << ZHadrons[3] << ' '  */
	/* 	 << endl; */
	
	// Init XZ
	XZ[0] = 0;
	XZ[1] = 0;
	XZ[2] = 0;
	XZ[3] = 0;
	
	// introduce a 4-array to get the POutLept into the Boost method
	Double_t helpPOut[4];
	// fix here - this is 2 or 3 in ktbrei - do the rotation to xz
	helpPOut[0] = POutLepton.Px();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		*/
	
		// check if not negative
		if (OldVector[0]!=0. || OldVector[1]!=0.)
		{
		NewVector[0] = 1.;
		NewVector[1] = 0.;
		NewVector[2] = 0.;
		NewVector[3] = 1.;
		OldVector[2] = 0.;
		/*cout << "new2 = " 
			<< NewVector[0] << ' ' 
			<< NewVector[1] << ' ' 
			<< NewVector[2] << ' ' 
			<< NewVector[3] << ' ' 
			<< OldVector[0] << ' ' 
			<< OldVector[1] << ' ' 
			<< OldVector[2] << ' ' 
			<< OldVector[3] << ' ' 
			<< endl; */
		
		//---NOTE THAT A POTENTIALLY AWKWARD SPECIAL CASE IS AVERTED, BECAUSE IF
		//   OLD AND NEW ARE EXACTLY BACK-TO-BACK, THE ROTATION AXIS IS UNDEFINED
		//   BUT IN THAT CASE DKTRROT WILL USE THE Z AXIS, AS REQUIRED
		Dijets::RRot(RMatrix, OldVector, NewVector);
	
		}
		} 
	}
	
	// Invert the transformation if necessary
	if (Scheme==2)
	{
	Dijets::MatrixInv(RMatrix,HelpMatrix);
	for (iCol = 0;iCol<4;iCol++)
		for (jLin = 0;jLin < 4;jLin++)
		RMatrix[iCol][jLin] = HelpMatrix[iCol][jLin];
	}
	
	// Clear output list
	OutputCandList.Clear();
	
	// Apply the result for all vectors
	for (iParticle=0;iParticle<NumPart;iParticle++)
	{
	// get actual 4-vector candidate
	// For Energy --> mass assumption = 0 => massless
	TLorentzVector Candidate;
	if(InputCandList.At(iParticle)->InheritsFrom(TLorentzVector::Class()))
		{
		Candidate = *((TLorentzVector*)InputCandList.At(iParticle));
		}
	else
		{
		cout << "ERROR in Boost: InputCand of unknown type !" << endl; 
		}
	
	Cand4Vec[0] = Candidate.X();
	Cand4Vec[1] = Candidate.Y();
	Cand4Vec[2] = Candidate.Z();
	Cand4Vec[3] = Candidate.T();
	/*   if(m_Debug > 1) */
	/* 	cout << "Cand = "  */
	/* 	     << Cand4Vec[0] << ' ' << Cand4Vec[1] << ' ' */
	/* 	     << Cand4Vec[2] << ' ' << Cand4Vec[3] << endl; */
	
	// Boost Cand4Vec --> FinalCand4Vec
	for (jLin = 0;jLin < 4;jLin++)
		{
		FinalCand4Vec[jLin] = 0.;
		for (iCol = 0;iCol<4;iCol++)
		FinalCand4Vec[jLin] += RMatrix[jLin][iCol]*Cand4Vec[iCol]; 
		}
	
	//cout << "Final= " 
	//   << FinalCand4Vec[0] << ' '
	//   << FinalCand4Vec[1] << ' '
	//   << FinalCand4Vec[2] << ' '
	//   << FinalCand4Vec[3] 
	//   << endl;
	
	// Create particle with boosted 4-vector 
	TLorentzVector* Cand5 = new TLorentzVector( FinalCand4Vec[0],
							FinalCand4Vec[1],
							FinalCand4Vec[2],
							FinalCand4Vec[3]);
	OutputCandList.Add(Cand5);
// 	delete Cand5;
	// Fill in outputlist
	/* if(m_Debug > 1)  */
	/* 	{ */
	/* 	  cout << "-------- Output particle list (Boost) --------" << endl; */
	/* 	  cout << ((TLorentzVector*)Cand5)->Px() << ' ' */
	/* 	       << ((TLorentzVector*)Cand5)->Py() << ' ' */
	/* 	       << ((TLorentzVector*)Cand5)->Pz() << ' ' */
	/* 	       << ((TLorentzVector*)Cand5)->E()  */
	/* 	       << endl; */
	/* 	} */
	}  
	}

//.........这里部分代码省略.........
      //correct for muons
      for (int i = 0 ; i < nMu ; i++)
        {
           TLorentzVector * mu_v = (TLorentzVector*) (muons->At (i)) ;
   
           if (mu_v->Pt () > 3)
             {
               met->SetPx (met->Px () - mu_v->Px ()) ;
               met->SetPy (met->Py () - mu_v->Py ()) ;
             }
         }      
      //      if (met->Pt () < g_METMin) continue ; plots.increase (cutId++) ; //PG 11
//      if (((TLorentzVector*) (MET->At (0)))->Pt () < g_METMin) continue ; plots.increase (cutId++) ; //PG 10
      
      //PG 2 TAGS
      //PG ------

      TLorentzVector * primoTAG = (TLorentzVector*) (tagJets->At (0)) ;
      TLorentzVector * secondoTAG = (TLorentzVector*) (tagJets->At (1)) ; 
      //PG get the first two in pt
      if (primoTAG->Pt () < secondoTAG->Pt ())
        {
          primoTAG = (TLorentzVector*) (tagJets->At (1)) ;
          secondoTAG = (TLorentzVector*) (tagJets->At (0)) ; 
        }

      if (met->Pt () < g_METMin) continue ; plots.increase (cutId++) ; //PG 11                                                                                                           

      TLorentzVector total = *primoTAG + *secondoTAG + sumLEP ;
      if (total.Pt () < g_PtTotMax) continue ; plots.increase (cutId++) ; //PG 11                                                                                                 
      total += *met ;
      if (total.Pt () > g_PtTotMetMax) continue ; plots.increase (cutId++) ; //PG 11                                                                                                                                       

      if (primoTAG->Pt () < g_hardTAGPtMin) continue ; plots.increase (cutId++) ; //PG 12
      if (secondoTAG->Pt () < g_softTAGPtMin) continue ; plots.increase (cutId++) ; //PG 13
      
      if (primoTAG->Eta () * secondoTAG->Eta () > g_TAGDProdEtaMax) continue ; plots.increase (cutId++) ; //PG 14
      if (fabs (primoTAG->Eta () - secondoTAG->Eta ()) < g_TAGDetaMin) continue ; plots.increase (cutId++) ; //PG 15
      TLorentzVector sumTAG = *primoTAG + *secondoTAG ;
      if (sumTAG.M () < g_TAGMinv) continue ; plots.increase (cutId++) ; //PG 16

      //PG JET VETO
      //PG --------

      //PG loop over other jets
      double etaMean = 0.5*(primoTAG->Eta () + secondoTAG->Eta ()); 
      int ojetsNum = 0 ;
      for (int ojetIt = 0 ; ojetIt < otherJets->GetEntries () ; ++ojetIt)
        {
          if ( ((TLorentzVector*) (otherJets->At (ojetIt)))->Pt () < g_ojetPtMin) continue ;
          if ( ((TLorentzVector*) (otherJets->At (ojetIt)))->Eta () < primoTAG->Eta () ||
               ((TLorentzVector*) (otherJets->At (ojetIt)))->Eta () > secondoTAG->Eta ()) continue ;
	  //if ((((TLorentzVector*) (otherJets->At (ojetIt)))->Eta () - etaMean) > g_ojetEtaFromMean) continue;
          ++ojetsNum ;        
        } //PG loop over other jets
      if (ojetsNum > g_ojetsMaxNum) continue ; plots.increase (cutId++) ; //PG 17


      /*
      double primoDR = primoELE->DeltaR (*primoTAG) ;
      if (primoDR < primoELE->DeltaR (*secondoTAG)) primoDR = fabs (primoELE->Eta () - secondoTAG->Eta ()) ;
      
      double secondoDR = secondoELE->DeltaR (*primoTAG) ;
      if (secondoDR < secondoELE->DeltaR (*secondoTAG)) secondoDR = fabs (secondoELE->Eta () - secondoTAG->Eta ()) ;

      if (secondoDR < g_eleTagDR || primoDR < g_eleTagDR) continue ; plots.increase (cutId++) ; //PG 12
      */

      //PG Ztautau vetos
      //PG -------------
      
      //PG the two electrons should not be opposite to each other
      
      TVector2 primoLEPT (primoLEP.m_kine->X (), primoLEP.m_kine->Y ()) ;
      TVector2 secondoLEPT (secondoLEP.m_kine->X (), secondoLEP.m_kine->Y ()) ;
      TVector2 METT (met->X (), met->Y ()) ;

      double Sum = METT * primoLEPT + METT * secondoLEPT / (1 + primoLEPT * secondoLEPT) ;
      double Dif = METT * primoLEPT - METT * secondoLEPT / (1 - primoLEPT * secondoLEPT) ;
      
      TVector2 METT1 = 0.5 * (Sum + Dif) * primoLEPT ;
      TVector2 METT2 = 0.5 * (Sum - Dif) * secondoLEPT ;
      
      double ptNu1 = METT1.Mod () / cos (primoLEP.m_kine->Theta ()) ;
      double ptNu2 = METT2.Mod () / cos (secondoLEP.m_kine->Theta ()) ;

    } //PG loop over the events

  plots.normalize () ;

  delete tagJets  ;  
  delete otherJets  ;
  delete electrons  ;
  delete muons  ;    
  delete MET  ;      
  delete tracks  ;   

  return 0;
  
}