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

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

void fill(int const kf, TLorentzVector* b, double weight, TLorentzVector const& p1Mom, TLorentzVector const& p2Mom, TVector3 v00)
{
   int const centrality = floor(nCent * gRandom->Rndm());

   TVector3 const vertex = getVertex(centrality);
   // smear primary vertex
   // float const sigmaVertex = sigmaVertexCent[cent];
   // TVector3 const vertex(gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex));

   v00 += vertex;

   // smear momentum
   TLorentzVector const p1RMom = smearMom(0, p1Mom);
   TLorentzVector const p2RMom = smearMom(0, p2Mom);

   // smear position
   TVector3 const p1RPos = smearPosData(0, vertex.z(), centrality, p1RMom, v00);
   TVector3 const p2RPos = smearPosData(0, vertex.z(), centrality, p2RMom, v00);
   // TVector3 const kRPos = smearPos(kMom, kRMom, v00);
   // TVector3 const pRPos = smearPos(pMom, pRMom, v00);

   // reconstruct
   TLorentzVector const rMom = p1RMom + p2RMom;
   float const p1Dca = dca(p1Mom.Vect(), v00, vertex);
   float const p2Dca = dca(p2Mom.Vect(), v00, vertex);
   float const p1RDca = dca(p1RMom.Vect(), p1RPos, vertex);
   float const p2RDca = dca(p2RMom.Vect(), p2RPos, vertex);

   TVector3 v0;
   float const dca12 = dca1To2(p1RMom.Vect(), p1RPos, p2RMom.Vect(), p2RPos, v0);
   float const decayLength = (v0 - vertex).Mag();
   float const dcaD0ToPv = dca(rMom.Vect(), v0, vertex);
   float const cosTheta = (v0 - vertex).Unit().Dot(rMom.Vect().Unit());
   float const angle12 = p1RMom.Vect().Angle(p2RMom.Vect());

   TLorentzVector p1RMomRest = p1RMom;
   TVector3 beta;
   beta.SetMagThetaPhi(rMom.Beta(), rMom.Theta(), rMom.Phi());
   p1RMomRest.Boost(-beta);
   float const cosThetaStar = rMom.Vect().Unit().Dot(p1RMomRest.Vect().Unit());

   // save
   float arr[100];
   int iArr = 0;
   arr[iArr++] = centrality;
   arr[iArr++] = vertex.X();
   arr[iArr++] = vertex.Y();
   arr[iArr++] = vertex.Z();

   arr[iArr++] = kf;
   arr[iArr++] = b->M();
   arr[iArr++] = b->Perp();
   arr[iArr++] = b->PseudoRapidity();
   arr[iArr++] = b->Rapidity();
   arr[iArr++] = b->Phi();
   arr[iArr++] = v00.X();
   arr[iArr++] = v00.Y();
   arr[iArr++] = v00.Z();

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

   arr[iArr++] = dca12;
   arr[iArr++] = decayLength;
   arr[iArr++] = dcaD0ToPv;
   arr[iArr++] = cosTheta;
   arr[iArr++] = angle12;
   arr[iArr++] = cosThetaStar;

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

   arr[iArr++] = p1RMom.M();
   arr[iArr++] = p1RMom.Perp();
   arr[iArr++] = p1RMom.PseudoRapidity();
   arr[iArr++] = p1RMom.Rapidity();
   arr[iArr++] = p1RMom.Phi();
   arr[iArr++] = p1RPos.X();
   arr[iArr++] = p1RPos.Y();
   arr[iArr++] = p1RPos.Z();
   arr[iArr++] = p1RDca;
   arr[iArr++] = tpcReconstructed(0,1,centrality,p1RMom);

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

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

void recursive_inter_with_air(){
	TFile *f1 = new TFile("muon.root", "recreate");
	TCanvas *c1 = new TCanvas("test", "test", 600, 700);

	TH1F * mu_c_1 = new TH1F("mu_c_1", "mu_c_1", 100, 0, 500);
	TH1F * mu_c_2 = new TH1F("mu_c_2", "mu_c_2", 100, 0, 500);
	TH1F * mu_c_3 = new TH1F("mu_c_3", "mu_c_3", 100, 0, 500);

	fn_muon_dxde = new TF1("f1", dxde_muon, M_mu, 1000, 1);
	TNtuple * ntuple = new TNtuple("ntuple", "ntuple", "nobs:detz:theta:press");
	// TNtuple *ntuple = new TNtuple("ntuple", "ntuple", "id:flag:e:px:py:pz:x:y:z");
	r = new TRandom();
	// r ->SetSeed(12232);
	for (int ievertex_nt = 0; ievertex_nt < 500; ievertex_nt++)
	{
		n_particle = 0;

		// init muon conunter
		int mu_counter1 = 0;
		int mu_counter2 = 0;
		int mu_counter3 = 0;

		double E = 1000; // GeV
		double theta = 0;
		// double theta = acos(r->Rndm());
		double phi = r->Rndm()*pi*2;
		double x = 0, y = 0, z = 1.0e10; //infinity
		double vertex_z, vertex_x, vertex_y; //km

		

		Get_Int_Posi(pdg_proton_p, theta, phi, x, y, z, vertex_x, vertex_y, vertex_z);

		Hillas_Split(pdg_proton_p, E, theta, phi, vertex_x, vertex_y, vertex_z);

		double proton_threshold = 10;
		double pion_threshold = 10;
		int i_particle = 0;
		while(i_particle < n_particle){
			int id = Ptcl_bank[i_particle].id;
			int flag = Ptcl_bank[i_particle].flag;
			double e = Ptcl_bank[i_particle].e;
			double x = Ptcl_bank[i_particle].vertex_x;
			double y = Ptcl_bank[i_particle].vertex_y;
			double z = Ptcl_bank[i_particle].vertex_z;
			if (id == pdg_proton_p && flag == 1 && e > proton_threshold)
			{
				Get_Int_Posi(id, theta, phi, x, y, z, vertex_x, vertex_y, vertex_z);
				if (vertex_z > 0)
				{
					Hillas_Split(id, e, theta, phi, vertex_x, vertex_y, vertex_z);
					Ptcl_bank[i_particle].flag = 0;
				}
			}
			if (id == pdg_pion_p &&flag == 1){

				Get_Int_Posi(id, theta, phi, x, y, z, vertex_x, vertex_y, vertex_z)	;			
				double l_to_interaction = sqrt((vertex_x-x)*(vertex_x-x)+(vertex_y-y)*(vertex_y-y)+(vertex_z-z)*(vertex_z-z));
				double l_to_decay = get_decay_length(pdg_pion_p, tau_pion, e);
				
				if (l_to_decay>l_to_interaction){
					if (vertex_z > 0){
						Hillas_Split(id, e, theta, phi, vertex_x, vertex_y, vertex_z);
						Ptcl_bank[i_particle].flag = 0;
					}
				}else{ // pion decay
					Ptcl_bank[i_particle].flag = 0;
					double px = Ptcl_bank[i_particle].px;
					double py = Ptcl_bank[i_particle].py;
					double pz = Ptcl_bank[i_particle].pz; // GeV
					TLorentzVector pion(px, py, pz, e);
					double decay_particle_mass[2] = {M_mu, M_neu};
					TGenPhaseSpace event;
					event.SetDecay(pion, 2, decay_particle_mass);
					event.Generate();
					TLorentzVector Muon = *(event.GetDecay(0));
					TLorentzVector Neu = *(event.GetDecay(1));

					double mu_theta = Muon.Theta();
					double mu_phi = Muon.Phi();
					
					if (mu_theta>pi/2.)	{
						double mu_E = Muon.E();
						double muon_decay_time = r ->Exp(tau_muon);
						double muon_decay_length = get_decay_length(pdg_muon_n, tau_muon, mu_E);

						float det_x = 0;
						float det_y = 0;
						float det_z = 0;
						mu_counter1++;
						double l_mu_vtx_to_det = sqrt(pow(vertex_x-det_x,2)+pow(vertex_y-det_y,2)+pow(vertex_z-det_z,2));
						// cout << "MuE" << mu_E << " Mass" <<M_mu<< " dl" << muon_decay_length << " lvd" << l_mu_vtx_to_det << endl;
						if (muon_decay_length > l_mu_vtx_to_det)	
						{
							mu_counter2++;

							double mu_stop_x;
							double mu_stop_y;
							double mu_stop_z;

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

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;
} 

//.........这里部分代码省略.........
        vector<HGeantKine *> vep;
        vector<HGeantKine *> vem;
        vector<HGeantKine *> vep_eff;
        vector<HGeantKine *> vem_eff;
        vector<HGeantKine *> vep_eff_multbin;
        vector<HGeantKine *> vem_eff_multbin;
	for(Int_t j = 0; j < size; j ++){
	    kine1 = HCategoryManager::getObject(kine1,kineCat,j);
            Float_t vx,vy,vz;
            kine1->getVertex(vx,vy,vz);
	    Float_t vr = TMath::Sqrt(vx*vx+vy*vy);
            if (vz < -60 || vz > 0) continue;
            if (vr > 2) continue;
            Int_t mamaNum = kine1->getParentTrack();

            if (kine1->isInAcceptance()) {
                if (kine1->getTotalMomentum() > 100 && kine1->getTotalMomentum() < 1000) {
                    Float_t px,py,pz;
                    kine1->getMomentum(px,py,pz);
                    TH2F *smear_matr;
                    if (kine1->getID() == 2) {
                        smear_matr = smear_pos;
                    } else {
                        smear_matr = smear_ele;
                    }
                    Float_t mom_ideal = kine1->getTotalMomentum();
                    Float_t mom_reco = smear(mom_ideal,smear_matr,random);
                    Float_t reco_over_ideal = mom_reco / mom_ideal;
                    kine1->setMomentum(px*reco_over_ideal,py*reco_over_ideal,pz*reco_over_ideal);
                    TLorentzVector vec;
                    HParticleTool::getTLorentzVector(kine1,vec,kine1->getID());

                    Float_t mom = vec.Vect().Mag();
                    Float_t the = vec.Theta()*TMath::RadToDeg();
                    Float_t phi = (vec.Phi()+TMath::Pi())*TMath::RadToDeg();
                    Float_t chg = (kine1->getID() == 2) ? 1 : -1;

                    if (kine1->getID() == 3) {
                        vem.push_back(new HGeantKine(*kine1));
                    }
                    if (kine1->getID() == 2) {
                        vep.push_back(new HGeantKine(*kine1));
                    }
                    Float_t eff = 1./getEfficiencyFactor(p3DEffEle[0][6],p3DEffPos[0][6],mom_ideal,the,phi,chg,false,false); // don't debug, don't check min value
                    if (isinf(eff) || isnan(eff)) eff = 0.;
                    if (random.Uniform(1) > eff) {
                        if (kine1->getID() == 3) {
                            vem_eff.push_back(new HGeantKine(*kine1));
                        }
                        if (kine1->getID() == 2) {
                            vep_eff.push_back(new HGeantKine(*kine1));
                        }
                    }
                    Float_t eff_multbin = 1./getEfficiencyFactor(p3DEffEle[multbin][6],p3DEffPos[multbin][6],mom,the,phi,chg,false,false);
                    if (isinf(eff_multbin) || isnan(eff_multbin)) eff_multbin = 0.;
                    if (random.Uniform(1) > eff_multbin) {
                        if (kine1->getID() == 3) {
                            vem_eff_multbin.push_back(new HGeantKine(*kine1));
                        }
                        if (kine1->getID() == 2) {
                            vep_eff_multbin.push_back(new HGeantKine(*kine1));
                        }
                    }
                }
            }
        }