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

void smearEmEnergy( const EnergyScaleCorrection_class& egcor, TLorentzVector& p, const ZGTree& myTree, float r9 ) {

    float smearSigma =  egcor.getSmearingSigma(myTree.run, fabs(p.Eta())<1.479, r9, p.Eta(), p.Pt(), 0., 0.);
    float cor = myRandom_.Gaus( 1., smearSigma );
    p.SetPtEtaPhiM( p.Pt()*cor, p.Eta(), p.Phi(), p.M() );

}

void fillCorrectEtaBin(TH1D* hist[nEtaBins], TLorentzVector sumVect, double phiSR, TH1D* polHist[nEtaBins][nPhisrBins], double polarization, TH1D* etaHist[nEtaBins]){
	
	for (int iEta = 0; iEta < nEtaBins; iEta++)
	{
		if (sumVect.Eta() >= etaBinStart[iEta] && sumVect.Eta() < etaBinEnd[iEta])
		{
			hist[iEta]->Fill(phiSR);
			etaHist[iEta]->Fill(sumVect.Eta());

			int phiSRbin = hist[iEta]->FindBin(phiSR);			
			polHist[iEta][phiSRbin-1]->Fill(polarization);
		}
	}
}

void ZgSelectData::printEvent(int index_e1, int index_e2, int index_pho, vector<Float_t> elePt_){
	cout << "PassedZeeGamma: run " << run << " event " << event << " rho " << rho << endl ;
	cout << "electron :index,pt,eta,phi,En " << endl ;
	cout << index_e1 << " , " << elePt_[index_e1] << " , " << eleEta[index_e1] ;
	cout << " , " << elePhi[index_e1] << " , " << eleEn[index_e1] << endl ;
	cout << index_e2 << " , " << elePt_[index_e2] << " , " << eleEta[index_e2] ;
	cout << " , " << elePhi[index_e2] << " , " << eleEn[index_e2] << endl ;

	cout << "Electron Pass " << endl ;
	cout << "TrackIso " << eleIsoTrkDR03[index_e1] << " EcalIso " << eleIsoEcalDR03[index_e1] << " HcalSolidIso ";
	cout << eleIsoHcalSolidDR03[index_e1] << " passID = " << elePassID(index_e1,3,elePt_) << endl ;
	cout << "TrackIso " << eleIsoTrkDR03[index_e2] << " EcalIso " << eleIsoEcalDR03[index_e2] << " HcalSolidIso ";
	cout << eleIsoHcalSolidDR03[index_e2] << " passID = " << elePassID(index_e2,3,elePt_) <<  endl ;
	TLorentzVector ve1, ve2, vZ ;
	ve1.SetPtEtaPhiE(elePt[index_e1],eleEta[index_e1],elePhi[index_e1],eleEn[index_e1]);
	ve2.SetPtEtaPhiE(elePt[index_e2],eleEta[index_e2],elePhi[index_e2],eleEn[index_e2]);
	vZ = ve1 + ve2 ;
	
	cout << "ZMass = " << vZ.M() << endl ;

	cout << "Photon : index, pt, eta, phi, dr1, dr2 " << endl ;
	TLorentzVector vp ;
	vp.SetPtEtaPhiM(phoEt[index_pho],phoEta[index_pho],phoPhi[index_pho],0);
	cout << index_pho << " , " << vp.Pt() << " , " << vp.Eta() << " , " << vp.Phi() ;
	cout << " , " << vp.DeltaR(ve1) << " , " << vp.DeltaR(ve2) << endl ;
	cout << "=================================================================" << endl ;

}

void UEAnalysisMPI::mpiAnalysisRECO(float weight,float etaRegion,float ptThreshold,TClonesArray * TracksJet)
{
  std::vector<TLorentzVector*> JetRECO;
  JetRECO.clear();
  
  for(int j=0;j<TracksJet->GetSize();++j){
    TLorentzVector *v = (TLorentzVector*)TracksJet->At(j);
    if(fabs(v->Eta())<etaRegion){
      JetRECO.push_back(v);
    }
  }
  
  std::vector<AssociatedObject> assoJetRECO;
  assoJetRECO.clear();

  while(JetRECO.size()>1){
    int oldSize = JetRECO.size();
    std::vector<TLorentzVector*>::iterator itH = JetRECO.begin();
    if((*itH)->Pt()>=ptThreshold){
      for(std::vector<TLorentzVector*>::iterator it=JetRECO.begin();it!=JetRECO.end();it++){
	float azimuthDistanceJet = fabs( (*itH)->Phi() - (*it)->Phi() );
	if((*it)->Pt()/(*itH)->Pt()>=0.3){
	  if( (piG - rangePhi) <  azimuthDistanceJet && azimuthDistanceJet < (piG + rangePhi)) {
	    AssociatedObject tmpPair((*itH),(*it));
	    assoJetRECO.push_back(tmpPair);
	    JetRECO.erase(it);
	    int newSize = oldSize -1;
	    oldSize = newSize;
	    JetRECO.resize(newSize);
	    break;
	  }
	}
      }
    }
    JetRECO.erase(itH);
    int newSize = oldSize -1;
    JetRECO.resize(newSize);
  }
  
  if(assoJetRECO.size()){
    fNumbMPIRECO->Fill(assoJetRECO.size());
    std::vector<AssociatedObject>::iterator at= assoJetRECO.begin();
    
    const TLorentzVector* leadingJet((*at).first);
    const TLorentzVector* secondJet((*at).second);

    pPtRatio_vs_PtJleadRECO->Fill(leadingJet->Pt(),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_EtaJleadRECO->Fill(fabs(leadingJet->Eta()),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_PhiJleadRECO->Fill(leadingJet->Phi(),(secondJet->Pt()/leadingJet->Pt()));
    
    fdEtaLeadingPairRECO->Fill(leadingJet->Eta()-secondJet->Eta());
    float dPhiJet = fabs(leadingJet->Phi()-secondJet->Phi());
    if(dPhiJet> piG) dPhiJet = 2*piG -dPhiJet;
    dPhiJet = (180*dPhiJet)/piG;
    fdPhiLeadingPairRECO->Fill(dPhiJet);
    fptRatioLeadingPairRECO->Fill(secondJet->Pt()/leadingJet->Pt());
  }
}

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

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

double Histogrammer::minDrPhoB(int PhoInd, EventTree* tree){
	// find the closest b-jet
	TLorentzVector b;
	TLorentzVector bBar;
	int phoGen=-1;
	double mindr = 999.0;
	for( int mcI = 0; mcI < tree->nMC_; ++mcI){
		if( tree->mcIndex->at(mcI) == tree->phoGenIndex_->at(PhoInd) ) 
			phoGen=mcI;
		if( tree->mcPID->at(mcI) == 5) 
			b.SetPtEtaPhiM(tree->mcPt->at(mcI), tree->mcEta->at(mcI), tree->mcPhi->at(mcI), tree->mcMass->at(mcI));
		if( tree->mcPID->at(mcI) == -5) 
			bBar.SetPtEtaPhiM(tree->mcPt->at(mcI), tree->mcEta->at(mcI), tree->mcPhi->at(mcI), tree->mcMass->at(mcI));
	}
	if( phoGen > 0 && b.Pt() > 0.0001 && bBar.Pt() > 0.0001 ) {
		mindr = std::min(dR(tree->mcEta->at(phoGen), tree->mcPhi->at(phoGen), b.Eta(), b.Phi()),
						 dR(tree->mcEta->at(phoGen), tree->mcPhi->at(phoGen), bBar.Eta(), bBar.Phi()));
	}
	return mindr;
}

TLorentzVector selectPhoton( const ZGConfig& cfg, const ZGTree& myTree, int index, const TLorentzVector& lept0, const TLorentzVector& lept1, EnergyScaleCorrection_class egcor ) {

    TLorentzVector photon;
    if( myTree.ngamma<=index ) {
        photon.SetPtEtaPhiM( 0.01, 0., 0., 0. );
        return photon;
    }
    photon.SetPtEtaPhiM( myTree.gamma_pt[index], myTree.gamma_eta[index], myTree.gamma_phi[index], myTree.gamma_mass[index] );

    // photon energy corrections/smearing
    if( !myTree.isData ) {
        if( cfg.smearing() ) {
            smearEmEnergy( egcor, photon, myTree, myTree.gamma_r9[0] );
        }
    } else {
        if( cfg.smearing() ) {
            applyEmEnergyScale( egcor, photon, myTree, myTree.gamma_r9[0] );
        }
    }

    bool goodPhoton = true;
    if( photon.Pt()<40. ) goodPhoton=false;
    if( fabs(photon.Eta())>1.4442 && fabs(photon.Eta())<1.566 ) goodPhoton=false;
    if( fabs(photon.Eta())>2.5 ) goodPhoton=false;
    if( myTree.gamma_idCutBased[index]==0 ) goodPhoton=false;
    if( myTree.gamma_chHadIso[index]>2.5 ) goodPhoton=false;
    if( fabs(myTree.gamma_eta[index])<1.44 ) {
        if( myTree.gamma_sigmaIetaIeta[index]>0.0102 ) goodPhoton=false;
    } else {
        if( myTree.gamma_sigmaIetaIeta[index]>0.0274 ) goodPhoton=false;
    }
    float deltaR_thresh = 0.4;
    if( photon.DeltaR(lept0)<deltaR_thresh || photon.DeltaR(lept1)<deltaR_thresh ) goodPhoton=false;


    if( !goodPhoton )
        photon.SetPtEtaPhiM( 0.01, 0., 0., 0. );

    return photon;

}

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 smearEmEnergy( TLorentzVector& p ) {

    float smearEBlowEta     = 0.013898;
    float smearEBhighEta    = 0.0189895;
    float smearEElowEta     = 0.027686;
    float smearEEhighEta    = 0.031312;
    float theGaussMean      = 1.;
    float pt = p.Pt();
    float eta = p.Eta();
    float phi = p.Phi();
    float mass = p.M();
    float theSmear = 0.;
    if      (fabs(eta)<1.                   )  theSmear = smearEBlowEta;
    else if (fabs(eta)>=1.  && fabs(eta)<1.5)  theSmear = smearEBhighEta;
    else if (fabs(eta)>=1.5 && fabs(eta)<2. )  theSmear = smearEElowEta;
    else if (fabs(eta)>=2.  && fabs(eta)<2.5)  theSmear = smearEEhighEta;
    float fromGauss = myRandom_.Gaus(theGaussMean,theSmear);
    p.SetPtEtaPhiM( fromGauss*pt, eta, phi, mass ); // keep mass and direction same

}

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

std::vector<double> CalculateDeltaEtaAndPhi(bool possibleMultW, bool possibleMultAntiW, int w1, TLorentzVector TL_munu){

  //Don't take into account the events with multiple W's. Think of something

  double deltaphi, deltaeta;
  if(possibleMultW == false && possibleMultAntiW == false){
    deltaphi = TMath::Abs(TVector2::Phi_mpi_pi(TL_munu.Phi() - genPhi[w1]) );
     deltaeta = TMath::Abs(TL_munu.Eta() - genEta[w1] );
  } else {
    deltaphi = 99.;
    deltaeta = 99.;
  }

  std::vector<double> EtaAndPhi;

  EtaAndPhi.push_back(deltaeta);
  EtaAndPhi.push_back(deltaphi);

  return EtaAndPhi;

}