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

void MomScaleCorrection (
  rochcor2012* rmcor, TLorentzVector& muMinus, TLorentzVector& muPlus,
  double rpx1, double rpy1, double rpz1, int charge1, 
  double rpx2, double rpy2, double rpz2, int charge2,
  float qter, int runOpt = 0, bool isData = true) {

  double mu_mass = 0.105658;
  double rp1 = rpx1*rpx1 + rpy1*rpy1 + rpz1*rpz1;
  double rp2 = rpx2*rpx2 + rpy2*rpy2 + rpz2*rpz2;
  double re1 = sqrt(mu_mass*mu_mass + rp1);
  double re2 = sqrt(mu_mass*mu_mass + rp2);

  if( charge1 > 0 ) {
    muPlus.SetPxPyPzE(rpx1, rpy1, rpz1, re1);
    muMinus.SetPxPyPzE(rpx2, rpy2, rpz2, re2);
  }
  else {
    muPlus.SetPxPyPzE(rpx2, rpy2, rpz2, re2);
    muMinus.SetPxPyPzE(rpx1, rpy1, rpz1, re1);
  }

  if( isData ) {
    rmcor->momcor_data(muPlus, 1, runOpt,qter);
    rmcor->momcor_data(muMinus, -1, runOpt,qter);
  }
  else {
    rmcor->momcor_mc(muPlus, 1, runOpt,qter);
    rmcor->momcor_mc(muMinus, -1, runOpt,qter);
  }
  //cout << "Muon = " << muPlus.Pt() << " " << muMinus.Pt() << endl;
}

 TLorentzVector GenLeptonPair::get4v(){
   TLorentzVector v;
   if(this->type == "ee") v.SetPxPyPzE((this->l1->px+this->l2->px),(this->l1->py+this->l2->py),(this->l1->pz+this->l2->pz),(sqrt(pow(this->l1->px,2)+pow(this->l1->py,2)+pow(this->l1->pz,2)+pow(MASS_EL,2))+sqrt(pow(this->l2->px,2)+pow(this->l2->py,2)+pow(this->l2->pz,2)+pow(MASS_EL,2))));
   if(this->type == "em") v.SetPxPyPzE((this->l1->px+this->l2->px),(this->l1->py+this->l2->py),(this->l1->pz+this->l2->pz),(sqrt(pow(this->l1->px,2)+pow(this->l1->py,2)+pow(this->l1->pz,2)+pow(MASS_EL,2))+sqrt(pow(this->l2->px,2)+pow(this->l2->py,2)+pow(this->l2->pz,2)+pow(MASS_MU,2))));
   if(this->type == "me") v.SetPxPyPzE((this->l1->px+this->l2->px),(this->l1->py+this->l2->py),(this->l1->pz+this->l2->pz),(sqrt(pow(this->l1->px,2)+pow(this->l1->py,2)+pow(this->l1->pz,2)+pow(MASS_EL,2))+sqrt(pow(this->l1->px,2)+pow(this->l1->py,2)+pow(this->l1->pz,2)+pow(MASS_EL,2))));
   if(this->type == "mm") v.SetPxPyPzE((this->l1->px+this->l2->px),(this->l1->py+this->l2->py),(this->l1->pz+this->l2->pz),(sqrt(pow(this->l1->px,2)+pow(this->l1->py,2)+pow(this->l1->pz,2)+pow(MASS_EL,2))+sqrt(pow(this->l2->px,2)+pow(this->l2->py,2)+pow(this->l2->pz,2)+pow(MASS_MU,2))));
   return v;
 }

void rochcor2012::musclefit_data( TLorentzVector& mu, TLorentzVector& mubar){

  float dpar1 = 0.0;
  float dpar2 = 0.0;
  float epar1 = 0.0;
  float epar2 = 0.0;
  
  if(fabs(mu.PseudoRapidity())<=0.9){
    dpar1 = d0par;
    epar1 = e0par;
  }else if(mu.PseudoRapidity()>0.9){
    dpar1 = d1par;
    epar1 = e1par;
  }else if(mu.PseudoRapidity()<-0.9){
    dpar1 = d2par;
    epar1 = e2par;
  }

  if(fabs(mubar.PseudoRapidity())<=0.9){
    dpar2 = d0par;
    epar2 = e0par;
  }else if(mubar.PseudoRapidity()>0.9){
    dpar2 = d1par;
    epar2 = e1par;
  }else if(mubar.PseudoRapidity()<-0.9){
    dpar2 = d2par;
    epar2 = e2par;
  }

  float corr1 = 1.0 + bpar*mu.Pt() + (-1.0)*cpar*mu.Pt()*TMath::Sign(float(1.0),float(mu.PseudoRapidity()))*TMath::Power(mu.PseudoRapidity(),2)
    + (-1.0)*dpar1*mu.Pt()*sin(mu.Phi() + epar1);
  float corr2 = 1.0 + bpar*mubar.Pt() + (1.0)*cpar*mubar.Pt()*TMath::Sign(float(1.0),float(mubar.PseudoRapidity()))*TMath::Power(mubar.PseudoRapidity(),2)
    + (1.0)*dpar2*mubar.Pt()*sin(mubar.Phi() + epar2);
  
  float px1 = mu.Px();
  float py1 = mu.Py();
  float pz1 = mu.Pz();
  float e1 = mu.E();
  
  float px2 = mubar.Px();
  float py2 = mubar.Py();
  float pz2 = mubar.Pz();
  float e2 = mubar.E();

  px1 *= corr1;
  py1 *= corr1;
  pz1 *= corr1;
  e1 *= corr1;
  
  px2 *= corr2;
  py2 *= corr2;
  pz2 *= corr2;
  e2 *= corr2;
  
  mu.SetPxPyPzE(px1,py1,pz1,e1);
  mubar.SetPxPyPzE(px2,py2,pz2,e2);

}

double  METMHTAsys(MissingET* met,vector<Jet> jetvec,vector<Muon> muonvec,vector<Electron> electronvec,vector<Photon> photonvec){
  double Met=-99;
  double METAsys=-99;
  TVector2 PUCorMet, RawMet;
   TLorentzVector allvecsum;
  allvecsum.SetPxPyPzE(0, 0, 0, 0);
  PUCorMet.Set(0., 0.);
  RawMet.Set(0.0, 0.0);
  
  for(int i=0; i<(int)jetvec.size(); i++) {allvecsum += jetvec.at(i).P4();}
  for(int j=0; j<(int)muonvec.size(); j++) {allvecsum += muonvec.at(j).P4();}
  for(int k=0; k<(int)electronvec.size(); k++) {allvecsum += electronvec.at(k).P4();}
  for(int l=0; l<(int)photonvec.size(); l++) {allvecsum += photonvec.at(l).P4();}

  PUCorMet.Set(-allvecsum.Px(),-allvecsum.Py());
  Met= PUCorMet.Mod();
  RawMet.SetMagPhi(met->MET, met->Phi);
  
  METAsys=fabs(Met-(RawMet.Mod()))/(Met+(RawMet.Mod()));//this is funny. RawMet.Mod() must return met->MET. We didn't need to build RawMet to obtain its magnitude :):0
  //cout << "......................RawMet.Mod(): " << RawMet.Mod() << endl;
  //cout << "...................... Met: " << Met << endl; 
  //cout << "...................... METAsys: " << METAsys << endl;
  return METAsys;

}

double t1pfmetPhiCorr(double t1pfmet, double t1pfmetPhi, double t1pfmetSumEt, bool isData){
Double_t t1pfmetCorrX;
Double_t t1pfmetCorrY;
Double_t t1pfmetCorrE;
 if(isData){
   //t1pfmetCorrX = t1pfmet*cos(t1pfmetPhi) - (-3.56  -0.0228096*t1pfmetSumEt);                                                    
   //t1pfmetCorrY = t1pfmet*sin(t1pfmetPhi) - (3.23  -0.0191691*t1pfmetSumEt);                                                     
   t1pfmetCorrX = t1pfmet*cos(t1pfmetPhi) - (-3.54059 + -0.0243867*t1pfmetSumEt);
   t1pfmetCorrY = t1pfmet*sin(t1pfmetPhi) - (3.25566 + -0.0195691*t1pfmetSumEt);
 }else{
   //t1pfmetCorrX = t1pfmet*cos(t1pfmetPhi) - (-1.34 - 0.0073*t1pfmetSumEt);                                                       
   //t1pfmetCorrY = t1pfmet*sin(t1pfmetPhi) - (0.263249  -0.0045*t1pfmetSumEt);                                                    
   t1pfmetCorrX = t1pfmet*cos(t1pfmetPhi) - (-1.33018 + -0.00738871*t1pfmetSumEt);
   t1pfmetCorrY = t1pfmet*sin(t1pfmetPhi) - (0.275181 + -0.00550282*t1pfmetSumEt);
 }


//std::cout << "px = t1pfmet*cos(t1pfmetPhi) - (" << fMETCorr[0] << " + " << fMETCorr[1] << "*t1pfmetSumEt)" << std::endl;
//std::cout << "py = t1pfmet*sin(t1pfmetPhi) - (" << fMETCorr[2] << " + " << fMETCorr[3] << "*t1pfmetSumEt)" << std::endl;
t1pfmetCorrE = sqrt(t1pfmetCorrX*t1pfmetCorrX + t1pfmetCorrY*t1pfmetCorrY);
TLorentzVector correctedMet;
correctedMet.SetPxPyPzE(t1pfmetCorrX,t1pfmetCorrY,0,t1pfmetCorrE);
Double_t t1pfmetPhiCorrected = correctedMet.Phi(); 
Double_t t1pfmetCorrected = correctedMet.Pt();

return t1pfmetPhiCorrected;
}

void printJet(Jet *jet)
{

  GenParticle *particle;
  Muon *muon;

  Track *track;
  Tower *tower;

  TObject *object;
  TLorentzVector momentum;
      momentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
      //TRefArray constituentarray(jet->Constituents);
      TRefArray particlearray(jet->Particles);
      cout<<"Looping over jet constituents. Jet pt: "<<jet->PT<<", eta: "<<jet->Eta<<", phi: "<<jet->Phi<<endl;      

      // Loop over all jet's constituents
      for(Int_t j = 0; j < jet->Constituents.GetEntriesFast(); ++j)
      {
        object = jet->Constituents.At(j);
        // Check if the constituent is accessible
        if(object == 0) continue;

        if(object->IsA() == GenParticle::Class())
        {
          particle = (GenParticle*) object;
          cout << "    GenPart pt: " << particle->PT << ", eta: " << particle->Eta << ", phi: " << particle->Phi << endl;
          momentum += particle->P4();
        }
        else if(object->IsA() == Track::Class())
        {
          track = (Track*) object;
          cout << "    Track pt: " << track->PT << ", eta: " << track->Eta << ", phi: " << track->Phi << endl;
          momentum += track->P4();
        }
        else if(object->IsA() == Tower::Class())
        {
          tower = (Tower*) object;
          cout << "    Tower pt: " << tower->ET << ", eta: " << tower->Eta << ", phi: " << tower->Phi << endl;
          momentum += tower->P4();
        }
        else if(object->IsA() == Muon::Class())
        {
          muon = (Muon*) object;
          cout << "    Muon pt: " << muon->PT << ", eta: " << muon->Eta << ", phi: " << muon->Phi << endl;
          momentum += muon->P4();
        }
      }
      cout << " constituent sum pt:  " << momentum.Pt() <<" eta "<< momentum.Eta()  <<"  phi " << momentum.Phi() << std::endl;


      for (Int_t j =0; j<jet->Particles.GetEntries();  j++){
     		GenParticle *p_tmp = (GenParticle*) particlearray.At(j);
		printGenParticle(p_tmp);
	}
}

void rochcor2012::momcor_mc( TLorentzVector& mu, float charge, float sysdev, int runopt,bool sync=false){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  //float mptsys = sran.Gaus(0.0,sysdev);
  
  float dm = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/mmavg[mu_phibin][mu_etabin];
  float da = mcor_ma[mu_phibin][mu_etabin] + mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin];
  
  float cor = 1.0/(1.0 + dm + charge*da*ptmu);

  //for the momentum tuning - eta,phi,Q correction
  px *= cor;
  py *= cor;
  pz *= cor;
  e *= cor;
  
  float gscler = 0.0;
  float deltaer = 0.0;
  float sfer = 0.0;
  
  gscler = TMath::Sqrt( TMath::Power(mgscl_stat,2) + TMath::Power(mgscl_syst,2) );
  deltaer = TMath::Sqrt( TMath::Power(delta_stat,2) + TMath::Power(delta_syst,2) );
  sfer = TMath::Sqrt( TMath::Power(sf_stat,2) + TMath::Power(sf_syst,2) );

  float tune;
  if (sync) tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*(1.0 + (sf + sysdev*sfer)));
  else tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*eran.Gaus(1.0,(sf + sysdev*sfer)));
  
  px *= (tune);
  py *= (tune);
  pz *= (tune);
  e *= (tune);
  
  float gscl = (genm_smr/mrecm);
  
  px *= (gscl + gscler_mc_dev*gscler);
  py *= (gscl + gscler_mc_dev*gscler);
  pz *= (gscl + gscler_mc_dev*gscler);
  e *= (gscl + gscler_mc_dev*gscler);

  mu.SetPxPyPzE(px,py,pz,e);
  
}

hades_track gen_track(float b_x, float b_y, float b_z, float v_px, float v_py, float v_pz, float v_e)
{
	TVector3 base(b_x, b_y, b_z);
	TLorentzVector vec;
	vec.SetPxPyPzE(v_px, v_py, v_pz, v_e);

	hades_track ht(base, vec);

	return ht;
}

void rochcor2012::momcor_mc( TLorentzVector& mu, float charge, float sysdev, int runopt, float& qter){
  
  //sysdev == num : deviation = num

  float ptmu = mu.Pt();
  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();

  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  //float mptsys = sran.Gaus(0.0,sysdev);
  
  float Mf = (mcor_bf[mu_phibin][mu_etabin] + mptsys_mc_dm[mu_phibin][mu_etabin]*mcor_bfer[mu_phibin][mu_etabin])/(mpavg[mu_phibin][mu_etabin]+mmavg[mu_phibin][mu_etabin]);
  float Af = ((mcor_ma[mu_phibin][mu_etabin]+mptsys_mc_da[mu_phibin][mu_etabin]*mcor_maer[mu_phibin][mu_etabin]) - Mf*(mpavg[mu_phibin][mu_etabin]-mmavg[mu_phibin][mu_etabin]));     
  
  float cor = 1.0/(1.0 + 2.0*Mf + charge*Af*ptmu);
  
  //for the momentum tuning - eta,phi,Q correction
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  float gscler = mgscl_stat;
  float deltaer = delta_stat;
  float sfer = sf_stat;
    
  float gscl = (genm_smr/mrecm);
  
  px *= (gscl + gscler_mc_dev*gscler);
  py *= (gscl + gscler_mc_dev*gscler);
  pz *= (gscl + gscler_mc_dev*gscler);
  e  *= (gscl + gscler_mc_dev*gscler);
  
  float momscl = sqrt(px*px + py*py)/ptmu;
  
  float tune = 1.0/(1.0 + (delta + sysdev*deltaer)*sqrt(px*px + py*py)*eran.Gaus(1.0,(sf + sysdev*sfer)));
  
  px *= (tune); 
  py *= (tune);  
  pz *= (tune);  
  e  *= (tune);   
  
  qter *= (momscl*momscl + (1.0-tune)*(1.0-tune));
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

TLorentzVector buildP (const LHEF::HEPEUP & event, int iPart)
{
  TLorentzVector dummy ;
  dummy.SetPxPyPzE (
      event.PUP.at (iPart).at (0), // px
      event.PUP.at (iPart).at (1), // py
      event.PUP.at (iPart).at (2), // pz
      event.PUP.at (iPart).at (3) // E
    ) ;
  return dummy ;  
}

void getKinematics(TLorentzVector& b, double const mass)
{
   float const pt = gRandom->Uniform(momentumRange.first, momentumRange.second);
   float const y = gRandom->Uniform(-acceptanceRapidity, acceptanceRapidity);
   float const phi = TMath::TwoPi() * gRandom->Rndm();

   float const mT = sqrt(mass * mass + pt * pt);
   float const pz = mT * sinh(y);
   float const E = mT * cosh(y);

   b.SetPxPyPzE(pt * cos(phi), pt * sin(phi) , pz, E);
}

void rochcor2012::momcor_data( TLorentzVector& mu, float charge, int runopt, float& qter){
  
  float ptmu = mu.Pt();

  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();
  
  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);

  float Mf = 0.0;
  float Af = 0.0;

  if(runopt==0){
    Mf = (dcor_bf[mu_phibin][mu_etabin]+mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfer[mu_phibin][mu_etabin])/(dpavg[mu_phibin][mu_etabin]+dmavg[mu_phibin][mu_etabin]);
    Af = ((dcor_ma[mu_phibin][mu_etabin]+mptsys_da_da[mu_phibin][mu_etabin]*dcor_maer[mu_phibin][mu_etabin]) - Mf*(dpavg[mu_phibin][mu_etabin]-dmavg[mu_phibin][mu_etabin]));     
  }else if(runopt==1){
    Mf = (dcor_bfD[mu_phibin][mu_etabin]+mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfDer[mu_phibin][mu_etabin])/(dpavgD[mu_phibin][mu_etabin]+dmavgD[mu_phibin][mu_etabin]);
    Af = ((dcor_maD[mu_phibin][mu_etabin]+mptsys_da_da[mu_phibin][mu_etabin]*dcor_maDer[mu_phibin][mu_etabin]) - Mf*(dpavgD[mu_phibin][mu_etabin]-dmavgD[mu_phibin][mu_etabin]));     
  }
  
  float cor = 1.0/(1.0 + 2.0*Mf + charge*Af*ptmu);
  
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  //after Z pt correction
  float gscler = dgscl_stat;
  float gscl = (genm_smr/drecm);
  
  px *= (gscl + gscler_da_dev*gscler);
  py *= (gscl + gscler_da_dev*gscler);
  pz *= (gscl + gscler_da_dev*gscler);
  e  *= (gscl + gscler_da_dev*gscler);
  
  float momscl = sqrt(px*px + py*py)/ptmu;
  qter *= momscl;
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

void rochcor::momcor_data( TLorentzVector& mu, float charge, float sysdev, int runopt){
  
  float ptmu = mu.Pt();

  float muphi = mu.Phi();
  float mueta = mu.Eta(); // same with mu.Eta() in Root

  float px = mu.Px();
  float py = mu.Py();
  float pz = mu.Pz();
  float e = mu.E();
  
  int mu_phibin = phibin(muphi);
  int mu_etabin = etabin(mueta);
  
  //float mptsys1 = sran.Gaus(0.0,sysdev);
  
  float dm = (dcor_bf[mu_phibin][mu_etabin] + mptsys_da_dm[mu_phibin][mu_etabin]*dcor_bfer[mu_phibin][mu_etabin])/dmavg[mu_phibin][mu_etabin];
  float da = dcor_ma[mu_phibin][mu_etabin] + mptsys_da_da[mu_phibin][mu_etabin]*dcor_maer[mu_phibin][mu_etabin];
  
  float cor = 1.0/(1.0 + dm + charge*da*ptmu);
  
  px *= cor;
  py *= cor;
  pz *= cor;
  e  *= cor;
  
  //after Z pt correction
  float gscler = 0.0;
  
  gscler = TMath::Sqrt( TMath::Power(dgscl_stat,2) + TMath::Power(dgscl_syst,2) );
  
  float gscl = (genm_smr/drecm);
  
  px *= (gscl + gscler_da_dev*gscler);
  py *= (gscl + gscler_da_dev*gscler);
  pz *= (gscl + gscler_da_dev*gscler);
  e  *= (gscl + gscler_da_dev*gscler);
  
  mu.SetPxPyPzE(px,py,pz,e);
  
}

void Output_Delphes::AnalyseParticles(ExRootTreeBranch *branch, const HepMC::GenEvent& evt)
{
  TRootC::GenParticle *element;
  TLorentzVector momentum;
  Double_t signPz;
  ReadStats();
  for(int n=1; n<=evt.particles_size(); n++) {
    int mo1, mo2, da1, da2, status, pid;
    getStatsFromTuple(mo1,mo2,da1,da2,status,pid,n);
    element = static_cast<TRootC::GenParticle*>(branch->NewEntry());
    element->PID = pid;
    element->Status = status;
    element->M1 = mo1 - 1; // added -1 as the numbering in the tree starts from 0
    element->M2 = mo2 - 1;
    element->D1 = da1 - 1;
    element->D2 = da2 - 1;
    element->E = index_to_particle[n]->momentum().e();
    element->Px = index_to_particle[n]->momentum().px();
    element->Py = index_to_particle[n]->momentum().py();
    element->Pz = index_to_particle[n]->momentum().pz();
    element->PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
    element->Eta = element->PT < 1e-6 ? signPz*999.9 : momentum.Eta();
    element->Phi = index_to_particle[n]->momentum().phi();
    HepMC::GenVertex* vrtI = (index_to_particle[n])->production_vertex();
    HepMC::GenVertex::particles_in_const_iterator partI;
    if(vrtI) {
      element->T = vrtI->position().t();
      element->X = vrtI->position().x();
      element->Y = vrtI->position().y();
      element->Z = vrtI->position().z();
    }
    else {
      element->T = 0.;
      element->X = 0.;
      element->Y = 0.;
      element->Z = 0.;
    }  
  }
}

void AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
{
  const LHEF::HEPEUP &hepeup = reader->hepeup;

  Int_t particle;
  Double_t signPz, cosTheta;

  TLorentzVector momentum;

  TRootLHEFParticle *element;

  for(particle = 0; particle < hepeup.NUP; ++particle)
  {
    element = (TRootLHEFParticle*) branch->NewEntry();

    element->PID = hepeup.IDUP[particle];
    element->Status = hepeup.ISTUP[particle];
    element->Mother1 = hepeup.MOTHUP[particle].first;
    element->Mother2 = hepeup.MOTHUP[particle].second;
    element->ColorLine1 = hepeup.ICOLUP[particle].first;
    element->ColorLine2 = hepeup.ICOLUP[particle].second;
    element->Px = hepeup.PUP[particle][0];
    element->Py = hepeup.PUP[particle][1];
    element->Pz = hepeup.PUP[particle][2];
    element->E = hepeup.PUP[particle][3];
    element->M = hepeup.PUP[particle][4];

    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);

    cosTheta = TMath::Abs(momentum.CosTheta());
    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;

    element->PT = momentum.Perp();
    element->Eta = (cosTheta == 1.0 ? signPz*999.9 : momentum.Eta());
    element->Phi = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());
    element->Rapidity = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());

    element->LifeTime = hepeup.VTIMUP[particle];
    element->Spin = hepeup.SPINUP[particle];
  }
}