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

void PhaseSpace() {

   if (!gROOT->GetClass("TGenPhaseSpace")) gSystem->Load("libPhysics");

   TLorentzVector target(0.0, 0.0, 0.0, 0.938);
   TLorentzVector beam(0.0, 0.0, .65, .65);
   TLorentzVector W = beam + target;

   //(Momentum, Energy units are Gev/C, GeV)
   Double_t masses[3] = { 0.938, 0.139, 0.139} ;

   TGenPhaseSpace event;
   event.SetDecay(W, 3, masses);

   TH2F *h2 = new TH2F("h2","h2", 50,1.1,1.8, 50,1.1,1.8);

   for (Int_t n=0;n<100000;n++) {
      Double_t weight = event.Generate();

      TLorentzVector *pProton = event.GetDecay(0);

      TLorentzVector *pPip    = event.GetDecay(1);
      TLorentzVector *pPim    = event.GetDecay(2);

      TLorentzVector pPPip = *pProton + *pPip;
      TLorentzVector pPPim = *pProton + *pPim;

      h2->Fill(pPPip.M2() ,pPPim.M2() ,weight);
   }
   h2->Draw();
}

//.........这里部分代码省略.........
      sortedVertex = selVtxSumPt2;

    if( sortedVertex.size() < 2 ) sortedVertex.push_back(-1);
    if( sortedVertex.size() < 3 ) sortedVertex.push_back(-1);

    if( lepCat >= 0 ) {
      /// lepton tag
      sortedVertex[0] = selVtx;  
      sortedVertex[1] = -1; 
      sortedVertex[2] = -1; 
      //      vertProb = 1;
    } 

    _minitree->mtree_rho   = rho2012;
    _minitree->mtree_rho25 = rho25;    
    _minitree->mtree_zVtx  = vtxbs[selVtx][2];
    _minitree->mtree_nVtx  = nVtxBS;
    _minitree->mtree_ivtx1 = selVtx;
    _minitree->mtree_ivtx2 = sortedVertex[1];
    _minitree->mtree_ivtx3 = sortedVertex[2];
    _minitree->mtree_vtxProb = vertProb;
    _minitree->mtree_vtxMva  = vertMVA;

    _minitree->mtree_mass  = hcand.M();
    _minitree->mtree_pt    = hcand.Pt();
    _minitree->mtree_piT   = hcand.Pt()/hcand.M();
    _minitree->mtree_y     = hcand.Rapidity();

    /// spin variables
    TLorentzVector gtmp1 = glead;  gtmp1.Boost( -hcand.BoostVector() );
    TLorentzVector gtmp2 = gtrail; gtmp2.Boost( -hcand.BoostVector() );
    _minitree->mtree_cThetaLead_heli  = cos( gtmp1.Angle(hcand.BoostVector()) );
    _minitree->mtree_cThetaTrail_heli = cos( gtmp2.Angle(hcand.BoostVector()) );
    _minitree->mtree_cThetaStar_CS    = 2*(glead.E()*gtrail.Pz() - gtrail.E()*glead.Pz())/(hcand.M()*sqrt(hcand.M2()+hcand.Pt()*hcand.Pt()));
    
    /// fill photon id variables in main tree
    _minitree->mtree_minR9      = +999;
    _minitree->mtree_minPhoIdEB = +999;
    _minitree->mtree_minPhoIdEE = +999;
    _minitree->mtree_maxSCEta   =   -1;
    _minitree->mtree_minSCEta   = +999;
    for( int i = 0 ; i < 2; i++ ) {
      int ipho = -1;
      if( i == 0 ) ipho = ilead;
      if( i == 1 ) ipho = itrail;
      fillPhotonVariablesToMiniTree( ipho, selVtx, i );      
      if( _minitree->mtree_r9[i] < _minitree->mtree_minR9 ) _minitree->mtree_minR9 = _minitree->mtree_r9[i];
      if( fabs( _minitree->mtree_sceta[i] ) <  1.5 &&  _minitree->mtree_mvaid[i] <  _minitree->mtree_minPhoIdEB ) _minitree->mtree_minPhoIdEB = _minitree->mtree_mvaid[i];
      if( fabs( _minitree->mtree_sceta[i] ) >= 1.5 &&  _minitree->mtree_mvaid[i] <  _minitree->mtree_minPhoIdEE ) _minitree->mtree_minPhoIdEE = _minitree->mtree_mvaid[i];
      if( fabs( _minitree->mtree_sceta[i] ) > _minitree->mtree_maxSCEta ) _minitree->mtree_maxSCEta =  fabs(_minitree->mtree_sceta[i]);
      if( fabs( _minitree->mtree_sceta[i] ) < _minitree->mtree_minSCEta ) _minitree->mtree_minSCEta =  fabs(_minitree->mtree_sceta[i]);
    }
    
    //------------ compute diphoton mva (add var to minitree inside function) ----------------//
    massResoCalc.setP4CalPosVtxResoSmear( glead,gtrail, 
					  TVector3(phoCaloPos[ilead ][0], phoCaloPos[ilead ][1],phoCaloPos[ilead ][2]),
					  TVector3(phoCaloPos[itrail][0], phoCaloPos[itrail][1],phoCaloPos[itrail][2]),
					  TVector3(vtxbs[selVtx][0], vtxbs[selVtx][1], vtxbs[selVtx][2]),
					  _minitree->mtree_relResOverE, _minitree->mtree_relSmearing );

    
    _minitree->mtree_massResoTot = massResoCalc.relativeMassResolutionFab_total( vertProb );    
    _minitree->mtree_massResoEng = massResoCalc.relativeMassResolutionFab_energy( );    
    _minitree->mtree_massResoAng = massResoCalc.relativeMassResolutionFab_angular();    
   
    float diphotonmva = DiPhoID_MVA( glead, gtrail, hcand, massResoCalc, vertProb, 

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

//.........这里部分代码省略.........
  displayTLorentzVector(&polvec);

  }

  if ( (ampltiude == 3) || (ampltiude == 4) || (ampltiude == 34) )
  {
  	 printf("Spin polarization vector defined in rest frame of the muon\n");
  	 printf("s^{mu} = (0,phat34)\n");

   printf("\n");
  	 printf("where the phat34 is pointing towards:\n");
  	 printf("phat34 theta: %12.6f\n", phat34_theta);
  	 printf("phat34 phi:   %12.6f\n", phat34_phi);

   printf("\n");
  	 printf("phat34 components:\n");
  	 printf("x: %12.6f\n", phat34[0]);
  	 printf("y: %12.6f\n", phat34[1]);
  	 printf("z: %12.6f\n", phat34[2]);

  }

  printf("\n");
  printf("--------------------------\n");
  printf("--- Consistency checks ---\n");
  printf("--------------------------\n");

  printf("\n");
  printf("- Lorentz scalar product (ps) = 0 (orthogonality check)\n");
  printf("(ps): %12.6e     (should give really small value)\n", (*P)*polvec);

  printf("\n");
  printf("- k0^{2} = 0 (massless auxuliary vector)\n");
  printf("(k0)^2: %12.6e     (should give really small value)\n", k0.M2());


  printf("\n");
  printf("------------------\n");
  printf("--- Ampltidude ---\n");
  printf("------------------\n");

  printf("\n");
  printf("Amplitude formula chosen: --->  %d  <--- \n", ampltiude);
  printf("\n");

  printf("Amplitude formula list:\n");
  printf("(No):  ---------------------- Formula ---------------------------- (polarization)\n");
  printf("- Default case:                              \n");
  printf(" (0): 128 * (p k') * (q k)                                         (unpolarized)\n");
  printf("(-1):  64 * (p k') * (q k) + 64 * M * (s k') * (q k)               (polarized)\n");
  printf("(+1):  64 * (p k') * (q k) - 64 * M * (s k') * (q k)               (polarized)\n");

  printf("- Custom:                              \n");
  printf(" (+3): 64 * gamma * ( 1 - beta ) *     (M * (k')^{0}) * (q k) -    (polarized)\n");
  printf("      -64 * gamma * ( 1 - beta ) * M *    (s k')      * (q k)                 \n");
  printf(" (+4): 64 * gamma * ( 1 + beta ) *     (M * (k')^{0}) * (q k) +    (polarized)\n");
  printf("      +64 * gamma * ( 1 + beta ) * M *    (s k')      * (q k)                 \n");
  printf("(+34): sum of (+3) and (+4)                                        (unpolarized)\n");

  printf("\n");
  printf("Notations:\n");
  printf("(muon)- ---> (electron)- (nu_mu) (nu_electronbar)\n");
  printf("  p     --->      q        k            k' \n");
  printf("  p     --->      q        k1           k2 \n");
  printf("  P     --->     p1        p2           p3 \n");
  printf("\n");