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

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();
//.........这里部分代码省略.........

	void Boost(Double_t CM[4],Double_t Sign,Double_t ZHad[4],
			Double_t XZMomentum[4],Int_t NumPart,TObjArray &InputCandList, 
			TObjArray &OutputCandList,Int_t Scheme)
	{
	// Boost method
	// Corresponds to old DKTFRAM subroutine 
	
	// Rotation matrix
	Double_t RMatrix[4][4],HelpMatrix[4][4];
	Double_t OldVector[4], NewVector[4];
	
	// Declare 4-vector for input and output of boost
	Double_t Cand4Vec[4], FinalCand4Vec[4];
	
	// Loop. vars
	Int_t iCol, jLin, iParticle;
	
	// Init matrix RMatrix with id.
	for (iCol = 0;iCol<4;iCol++)
	for (jLin = 0;jLin < 4;jLin++)
	{
		if (iCol != jLin)
		RMatrix[iCol][jLin] = 0;
		else RMatrix[iCol][jLin] = 1;
	}
	
	// Do a Lorentz boost to/from cm frame of XXXX
	Dijets::LorentzBoost(RMatrix,CM,0);
	
	//for (iCol = 0;iCol<4;iCol++)
	//for (jLin = 0;jLin < 4;jLin++)
	// {
	//cout << "rmatrix2 = " << iCol << ' ' << jLin << ' ' << RMatrix[iCol][jLin] << endl;
	// }
	
	// Find rotation to put boosted Z on the (sign) z axis
	if (Sign!=0.)
	{
	for (jLin = 0;jLin<4;jLin++)
		{
		OldVector[jLin] = 0.;
		for (iCol = 0;iCol < 4;iCol++)
		OldVector[jLin] += RMatrix[jLin][iCol]*ZHad[iCol]; 
		}
	
	//cout << "old1 = " 
	//  << OldVector[0] << ' ' 
	//  << OldVector[1] << ' ' 
	//  << OldVector[2] << ' ' 
	//  << OldVector[3] << ' ' 
	//  << endl;
	
	// Check if not 0 and fill in new vector
	if (OldVector[0]!=0. || OldVector[1]!=0. || OldVector[2]!=0.)
		{
		NewVector[0] = 0.;
		NewVector[1] = 0.;
		NewVector[2] = Sign;
		NewVector[3] = fabs(Sign);
		//cout << "new1 = " 
		//   << NewVector[0] << ' ' 
		//   << NewVector[1] << ' ' 
		//   << NewVector[2] << ' ' 
		//   << NewVector[3] << endl;
		Dijets::RRot(RMatrix,OldVector,NewVector);
	
		// Find rotation to put boosted Z and rotated XZMomentum into XZ-plane
		// fix here
		for (jLin = 0;jLin<4;jLin++)
		{
		OldVector[jLin] = 0.;
		for (iCol = 0;iCol < 4;iCol++)
			OldVector[jLin] += RMatrix[jLin][iCol]*XZMomentum[iCol]; 
		}
	
		/*
		for (jLin = 0;jLin<4;jLin++)
		{
		cout << " rmat = " 
		<< RMatrix[jLin][0] << ' '
		<< RMatrix[jLin][1] << ' '
		<< RMatrix[jLin][2] << ' '
		<< RMatrix[jLin][3] << endl;
		}
		cout << "xzmomentum = " 
		<<  XZMomentum[0] << ' '
		<<  XZMomentum[1] << ' '
		<<  XZMomentum[2] << ' '
		<<  XZMomentum[3] << endl;
		cout << "old2 = " 
		<< OldVector[0] << ' ' 
		<< OldVector[1] << ' ' 
		<< OldVector[2] << ' ' 
		<< OldVector[3] << ' ' 
		<< endl;
		*/
	
		// check if not negative
		if (OldVector[0]!=0. || OldVector[1]!=0.)
		{
//.........这里部分代码省略.........

//!PG main function
int 
selector (TChain * tree, histos & plots)
{
  TClonesArray * tagJets = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("tagJets", &tagJets) ;
  TClonesArray * otherJets = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("otherJets", &otherJets) ;
  TClonesArray * electrons = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("electrons", &electrons) ;
  TClonesArray * muons = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("muons", &muons) ;
  TClonesArray * MET = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("MET", &MET) ;
  TClonesArray * tracks = new TClonesArray ("TLorentzVector") ;
  tree->SetBranchAddress ("tracks", &tracks) ;
  
  int EleId[100];
  float IsolEleSumPt[100];
  int nEle;
  tree->SetBranchAddress ("nEle", &nEle) ;
  tree->SetBranchAddress ("EleId",EleId ) ;
  tree->SetBranchAddress ("IsolEleSumPt",IsolEleSumPt ) ;

  float IsolMuSumPt[100];
  int nMu ;
  tree->SetBranchAddress ("nMu", &nMu) ;
  tree->SetBranchAddress ("IsolMuSumPt",IsolMuSumPt ) ;


  int nentries = (int) tree->GetEntries () ;
  //std::cout << "nentries " << nentries << std::endl;
  //PG loop over the events
  //  int nentries = 100 ;
  for (int evt = 0 ; evt < nentries ; ++evt)
    {
      tree->GetEntry (evt) ;
      int cutId = 0 ;

      plots.increase (cutId++) ; //PG 0

      if (tagJets->GetEntries () != 2) continue ; plots.increase (cutId++) ; //PG 1 FIXME ctrl numbering
      //if (electrons->GetEntries () < 1 ||
      //  muons->GetEntries () < 1) continue ; plots.increase (cutId++) ; //PG 2
            
      std::vector<lepton> leptons ;
            
      //PG pour electrons into leptons collection
      //PG ---------------------------------------

      //PG loop over electrons
      for (int iele = 0; iele < nEle ; ++iele)
        {
          TLorentzVector * theEle = (TLorentzVector*) (electrons->At (iele)) ;
          lepton dummy (theEle, 0, iele) ;
          leptons.push_back (dummy) ;
        } //PG loop over electrons

      //PG loop over muons
      for (int imu = 0 ; imu < nMu ; ++imu)
        {
          TLorentzVector * theMu = (TLorentzVector*) (muons->At (imu)) ;
          lepton dummy (theMu, 1, imu) ;
          leptons.push_back (dummy) ;
        } //PG loop over muons

//PG this check is not necessary
//PG      if (leptons.size () < 2) continue ;

      //PG 2 LEPTONS
      //PG ---------

/* 
applied after the leptons choice: 
in this case it is possible to differentiate the selections depending on the 
position of each lepton in the pt-sorting.
the algorithm searches the first two most energetic candidates which satisfy 
the ID selections required for the first and second lepton respectively.
If the so-identified first two leptons have the same flavour, the event
is rejected, since it is expected to fall into the 2e or 2mu sub-channel.
Is it certified that we do not have overlap?
*/

      sort (leptons.rbegin (), leptons.rend (), lessThan ()) ;

      lepton primoLEP ;
      lepton secondoLEP ;

      //PG find the first lepton
      int ilep = 0 ;
      for ( ; ilep < leptons.size () ; ++ilep)
        {
          if (leptons.at (ilep).m_flav == 0) //PG electron
            {
               //PG iso check
              bool eleIso = (IsolEleSumPt[leptons.at (ilep).m_index] /  
                             leptons.at (ilep).m_kine->Pt () ) < 0.2 ; // 0.2 per il momento
              if (g_ISO1[0] == 1 && eleIso != 1) continue;
              
              //PG eleID check
//.........这里部分代码省略.........