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

void Boost_To_Stop_Rest_Frame(TLorentzVector& stop4, TLorentzVector& chargino4, TLorentzVector& b4, TLorentzVector& neutralino4, TLorentzVector& W4, TLorentzVector& up4, TLorentzVector& down4, TLorentzVector& s4)
{
    TVector3 betaV(-stop4.Px()/stop4.Energy(),-stop4.Py()/stop4.Energy(),-stop4.Pz()/stop4.Energy());
    stop4.Boost(betaV);
    chargino4.Boost(betaV);
    b4.Boost(betaV);
    neutralino4.Boost(betaV);
    W4.Boost(betaV);
    up4.Boost(betaV);
    down4.Boost(betaV);
    s4.SetE(chargino4.P()/chargino4.M());
    s4.SetVect(chargino4.Vect().Unit()*chargino4.Gamma());
}

Bool_t monojet::Process(Long64_t entry)
{
    GetEntry(entry);

    if( entry % 100000 == 0 ) cout << "Processing event number: " << entry << endl;
    //cout << "Processing event number: " << entry << endl;

    // To make the processing fast, apply a very looose selection
    if (((TLorentzVector*)((*metP4)[0]))->Pt() < 40. or jetP4->GetEntries() < 1) return kTRUE;

    //this is the type tree
    tm->SetValue("run",runNum);
    tm->SetValue("event",eventNum);
    tm->SetValue("lumi",lumiNum);

    float dR = 0.;
    TClonesArray *tightLep;
    TClonesArray *cleanJet;
    TClonesArray *cleanTau;

    tightLep = new TClonesArray("TLorentzVector",20);
    cleanJet = new TClonesArray("TLorentzVector",20);
    cleanTau = new TClonesArray("TLorentzVector",20);

    std::vector<bool>  jetMonojetId_clean;
    jetMonojetId_clean.clear();
    std::vector<bool>  jetMonojetIdLoose_clean;
    jetMonojetIdLoose_clean.clear();
    //std::vector<float> jetPuId_clean;
    //jetPuId_clean.clear();

    std::vector<float>  tauId_clean;
    tauId_clean.clear();
    std::vector<float>  tauIso_clean;
    tauIso_clean.clear();

    int n_tightlep = 0;

    // ********* Leptons ********** //
    for(int lepton = 0; lepton < lepP4->GetEntries(); lepton++) {
        TLorentzVector* Lepton = (TLorentzVector*) lepP4->At(lepton);
        // check if this is a tight lep, and check the overlap

        //iso_1 = divide(input_tree.lepIso[0],input_tree.lepP4[0].Pt())
        //if (input_tree.lepTightId[0]==0 or iso_1 > 0.12): continue

        if (Lepton->Pt() > 20. && (*lepTightId)[lepton] > 1) {
            n_tightlep +=1;
            new ( (*tightLep)[tightLep->GetEntriesFast()]) TLorentzVector(Lepton->Px(), Lepton->Py(), Lepton->Pz(), Lepton->Energy());

            //check overlap with jets
            for(int j = 0; j < jetP4->GetEntries(); j++) {
                TLorentzVector* Jet = (TLorentzVector*) jetP4->At(j);
                dR = deltaR(Lepton,Jet);
                if (dR > dR_cut) {
                    new ( (*cleanJet)[cleanJet->GetEntriesFast()]) TLorentzVector(Jet->Px(), Jet->Py(), Jet->Pz(), Jet->Energy());
                    jetMonojetId_clean.push_back((*jetMonojetId)[j]);
                    jetMonojetIdLoose_clean.push_back((*jetMonojetIdLoose)[j]);
                    //jetPuId_clean.push_back((*jetPuId)[j]);
                }
            }
            //check overlap with taus
            for(int tau = 0; tau < tauP4->GetEntries(); tau++) {
                TLorentzVector* Tau = (TLorentzVector*) tauP4->At(tau);
                dR = deltaR(Lepton,Tau);
                if (dR > dR_cut) new ( (*cleanTau)[cleanTau->GetEntriesFast()]) TLorentzVector(Tau->Px(), Tau->Py(), Tau->Pz(), Tau->Energy());
                tauId_clean.push_back((*tauId)[tau]);
                tauIso_clean.push_back((*tauIso)[tau]);
            } // tau overlap
        } // tight lepton selection
    }//lepton loop

    tm->SetValue("n_tightlep",n_tightlep);

    TLorentzVector fakeMET;

    // Z Selection
    TLorentzVector Z;
    if(lepP4->GetEntries() == 2 && n_tightlep > 0) {
        if (((*lepPdgId)[0]+(*lepPdgId)[1])==0 ) {
            Z = *((TLorentzVector*)((*lepP4)[0])) + *((TLorentzVector*)((*lepP4)[1]));
            fakeMET = *((TLorentzVector*)((*metP4)[0])) + Z;
        }
    }

    float MT = 0.0;
    //// W Selection
    if(lepP4->GetEntries() == 1 && n_tightlep == 1) {
        fakeMET = *((TLorentzVector*)((*metP4)[0])) + *((TLorentzVector*)((*lepP4)[0])) ;
        MT = transverseMass( ((TLorentzVector*)((*lepP4)[0]))->Pt(), ((TLorentzVector*)((*lepP4)[0]))->Phi(), ((TLorentzVector*)((*metP4)[0]))->Pt(), ((TLorentzVector*)((*metP4)[0]))->Phi());
    }

    tm->SetValue("mt",MT);

    // ********* Jets ********** //
    for(int jet = 0; jet < jetP4->GetEntries(); jet++) {
        TLorentzVector* Jet = (TLorentzVector*) jetP4->At(jet);
        //cout << (*jetMonojetId)[0] <<endl;
        //cout << Jet->Pt()<<endl;
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
         leptons.push_back( electrons[1] );

       } else {

         leptons.push_back( electrons[1] );
         leptons.push_back( electrons[0] );

       }

     } else if( muons.size() == 2 ) {

       leptType_ = 0;

       if( muons[0].Pt() > muons[1].Pt() ) {

         leptons.push_back( muons[0] );
         leptons.push_back( muons[1] );

       } else {

         leptons.push_back( muons[1] );
         leptons.push_back( muons[0] );

       }

     } else {

     //std::cout << "There must be an error this is not possible." << std::endl;
     //exit(9101);

     }

     
     eLept1_ = (leptons.size()>0) ? leptons[0].Energy() : 0.;
     ptLept1_ = (leptons.size()>0) ? leptons[0].Pt() : 0.;
     etaLept1_ = (leptons.size()>0) ? leptons[0].Eta() : 0.;
     phiLept1_ = (leptons.size()>0) ? leptons[0].Phi() : 0.;
     
     eLept2_ = (leptons.size()>1) ? leptons[1].Energy() : 0.;
     ptLept2_ = (leptons.size()>1) ? leptons[1].Pt() : 0.;
     etaLept2_ = (leptons.size()>1) ? leptons[1].Eta() : 0.;
     phiLept2_ = (leptons.size()>1) ? leptons[1].Phi() : 0.;


     // --------------------
     // match leptons to MC:
     // --------------------
     int correctIdMc = (leptType_==0 ) ? 13 : 11;

     for( unsigned iLept=0; iLept<leptons.size(); ++iLept ) {

       float deltaRmin = 100.;
       TLorentzVector matchedLeptonMC;

       for( unsigned iMc=0; iMc<nMc; ++iMc ) {

         if( statusMc[iMc]==1 && fabs(idMc[iMc])==correctIdMc && idMc[mothMc[mothMc[iMc]]]==23 ) {

           TLorentzVector* thisParticle = new TLorentzVector();
           thisParticle->SetPtEtaPhiE( pMc[iMc]*sin(thetaMc[iMc]), etaMc[iMc], phiMc[iMc], energyMc[iMc] );
           float thisDeltaR = leptons[iLept].DeltaR( *thisParticle );
           if( thisDeltaR < deltaRmin ) {
             deltaRmin = thisDeltaR;
             matchedLeptonMC = *thisParticle;
           }

//.........这里部分代码省略.........
    elePt.clear(); eleEta.clear(); elePhi.clear(); eleEnr.clear();elescEta.clear();
    recoElec.SetPtEtaPhiE(0.,0.,0.,0.);
    eleP4->Clear();

    primVtx = nPV;

    if(!Ele23_WPLoose) continue;

    for(int i=0;i<ptElec->size();i++){

      recoElec.SetPtEtaPhiE(ptElec->at(i),etaElec->at(i), phiElec->at(i), energyElec->at(i));
      new ((*eleP4)[i]) TLorentzVector(recoElec);
      TLorentzVector* fourmom = (TLorentzVector*) eleP4->At(i);

   //  if(fourmom->Eta() < 1.4442) *fourmom *= 0.9994;
   //   if(fourmom->Eta() > 1.566) *fourmom *= 1.0034;

      if(fourmom->Eta() < 1.4442) *fourmom *= 1.0;
      if(fourmom->Eta() > 1.566) *fourmom *= 1.0;

     //   mediumId = passHEEPId->at(i);
       mediumId = passMediumId->at(i);

       passId = mediumId;

//      if (ptElec->at(i) < 500.) { passId = mediumId; }
//     else { passId = dzMaskId; }

      if(passId) {
	if(fabs(etaSC->at(i)) < 2.5 && !(fabs(etaSC->at(i)) > 1.4442 && fabs(etaSC->at(i)) < 1.566)){

	  newelePt.push_back(fourmom->Pt());
	  neweleEta.push_back(fourmom->Eta());
	  neweleEnr.push_back(fourmom->Energy());
	  newelePhi.push_back(fourmom->Phi());
          newscEta.push_back(etaSC->at(i));
	}
      }
    }

    // Sorting
    int index[newelePt.size()];
    float pt[newelePt.size()];

    for(unsigned int el=0; el<newelePt.size(); el++)
    { 
      pt[el]=newelePt.at(el);
    }

    int size = sizeof(pt)/sizeof(pt[0]);
    TMath::Sort(size,pt,index,true);

    if(newelePt.size()==2){

      for(unsigned int j = 0; j < pt_Ele23->size(); j++){

	double dR1_comp = 1000.;
	double dR2_comp = 1000.;

	dR1 = deltaR(neweleEta.at(index[0]), newelePhi.at(index[0]), eta_Ele23->at(j), phi_Ele23->at(j));
	dR2 = deltaR(neweleEta.at(index[1]), newelePhi.at(index[1]), eta_Ele23->at(j), phi_Ele23->at(j));

	if(dR1 < 0.1){
	  if (dR1 < dR1_comp)
	  {
	    dR1_comp = dR1;

//.........这里部分代码省略.........
    hGen_Onia_phi[0][0]->Fill(onia_phi);
    hGen_Onia_phi[rapIntegratedPTIndex][0]->Fill(onia_phi);
    hGen_Onia_phi[0][rapForPTIndex]->Fill(onia_phi);
      
    hGen_Onia_pt[0]->Fill(onia_pt);
    // hGen_Onia_eta[0]->Fill(onia_eta);
    // hGen_Onia_rap[0]->Fill(onia_rap);
    //b) individual pT and rap bins:
    hGen_Onia_mass[pTIndex][rapForPTIndex]->Fill(onia_mass);
    hGen_Onia_phi[pTIndex][rapForPTIndex]->Fill(onia_phi);
    hGen_Onia_pt[rapForPTIndex]->Fill(onia_pt);
    // hGen_Onia_eta[pTIndex]->Fill(onia_eta);
    // hGen_Onia_rap[pTIndex]->Fill(onia_rap);

    hGen_Onia_rap_pT->Fill(onia_rap, onia_pt);

    //=====================
    calcPol(*muPos, *muNeg);
    //=====================

    //test:
//     calcPol(*muNeg, *muPos);
//     //H: test:
//     if(jentry%2 == 0)
//       calcPol(*muPos, *muNeg);
//     else
//       calcPol(*muNeg, *muPos);

    //===================================================
    //calculate delta, the angle between the CS and HX frame
    //Formula from EPJC paper
    Double_t deltaHXToCS = TMath::ACos(onia_mass * onia->Pz() / (onia_mT * onia_P));
    //     Double_t deltaCSToHX = -deltaHXToCS;
    Double_t sin2Delta = pow((onia_pt * onia->Energy() / (onia_P * onia_mT)),2);
    //sin2Delta does not change sign when going from HX-->CS or vice versa
    hDelta[pTIndex][rapForPTIndex]->Fill(deltaHXToCS * 180./TMath::Pi());
    hSin2Delta[pTIndex][rapForPTIndex]->Fill(sin2Delta);
    //===================================================

    Double_t deltaPhi = muPos->Phi() - muNeg->Phi();
    if(deltaPhi < -TMath::Pi()) deltaPhi += 2.*TMath::Pi();
    else if(deltaPhi > TMath::Pi()) deltaPhi = 2.*TMath::Pi() - deltaPhi;

    //debugging histos
    hPhiPos_PhiNeg[pTIndex][rapForPTIndex]->Fill(180./TMath::Pi() * muNeg->Phi(), 180./TMath::Pi() * muPos->Phi());
    hPtPos_PtNeg[pTIndex][rapForPTIndex]->Fill(muNeg->Pt(), muPos->Pt());
    hEtaPos_EtaNeg[pTIndex][rapForPTIndex]->Fill(muNeg->PseudoRapidity(), muPos->PseudoRapidity());
    // hDeltaPhi[pTIndex][rapIndex]->Fill(deltaPhi);
    hDeltaPhi[pTIndex][rapForPTIndex]->Fill(deltaPhi);

    hGen_mupl_pt[pTIndex][rapForPTIndex]->Fill(muPos->Pt());
    hGen_mupl_eta[pTIndex][rapForPTIndex]->Fill(muPos->PseudoRapidity());
    hGen_mupl_phi[pTIndex][rapForPTIndex]->Fill(muPos->Phi());
      
    hGen_mumi_pt[pTIndex][rapForPTIndex]->Fill(muNeg->Pt());
    hGen_mumi_eta[pTIndex][rapForPTIndex]->Fill(muNeg->PseudoRapidity());
    hGen_mumi_phi[pTIndex][rapForPTIndex]->Fill(muNeg->Phi());

    //fill the histos for all the different frames
    for(int iFrame = 0; iFrame < jpsi::kNbFrames; iFrame++){

      thisCosPhi[iFrame] = TMath::Cos(2.*thisPhi_rad[iFrame]);

      Double_t weight = CalcPolWeight(onia_P, thisCosTh[iFrame]);

      //1a) polariztion histos - all pT

//.........这里部分代码省略.........
      	if( fTree.electrons_Pt[iele] > 20 &&
	    fabs(fTree.electrons_Eta[iele]) < 2.5 && 
	    ( fabs(fTree.electrons_Eta[iele]) < 1.4442 || fabs(fTree.electrons_Eta[iele]) > 1.566) &&
	    fTree.electrons_vidVeto[iele] > 0.5 )
	  nLooseElectrons++;
      }
      
      if( nLooseElectrons > 0 )
      	continue;

      if( printEventIds )
	(*(EventIdFiles[cutindex])) << (fTree.Event_EventNumber) << endl ;
      cutflowtable.Fill( cutindex , weight , EventsIsPSeudoData < fabs(weight) );
      cutflowtablew1.Fill( cutindex , weight/fabs(weight) , EventsIsPSeudoData < fabs(weight) );
      cutindex ++ ;

      int j1index = -1;
      int j2index = -1;
      int j3index = -1;
      int nJets = 0 ;
      int howManyBJets = 0;
      int bjIndex = -1;
      int bj2Index = -1;
      int jprimeIndex ; 

      int nJetsPt20_47 = 0;

      int nJetsPt20_24 = 0;
      int nLbJetsPt20 = 0;
      int nTbJetsPt20 = 0;

      TLorentzVector muon;
      muon.SetPtEtaPhiE( fTree.muons_Pt[tightMuIndex] , fTree.muons_Eta[tightMuIndex] , fTree.muons_Phi[tightMuIndex] , fTree.muons_E[tightMuIndex] );
      if( muon.Energy() == 0.0 )
	cout << "ZERO???" <<endl;
      for( int jid = 0 ; jid < fTree.jetsAK4_size ; jid++ ){

	float NHF = fTree.jetsAK4_JetID_neutralHadronEnergyFraction[jid] ;
	float NEMF = fTree.jetsAK4_JetID_neutralEmEnergyFraction[jid] ;
	float NumConst = fTree.jetsAK4_JetID_numberOfDaughters[jid] ;
	float eta = fTree.jetsAK4_Eta[jid] ;
	float CHF = fTree.jetsAK4_JetID_chargedHadronEnergyFraction[jid] ;
	float CHM = fTree.jetsAK4_JetID_chargedMultiplicity[jid] ;
	float CEMF = fTree.jetsAK4_JetID_chargedEmEnergyFraction[jid] ;
	float NumNeutralParticle = fTree.jetsAK4_JetID_neutralMultiplicity[jid] ;


	bool looseid ;
	if(  abs(eta)<=3.0 ) {
	  looseid = (NHF<0.99 && NEMF<0.99 && NumConst>1) && ((abs(eta)<=2.4 && CHF>0 && CHM>0 && CEMF<0.99) || abs(eta)>2.4) ;
	}
	else
	  looseid = (NEMF<0.90 && NumNeutralParticle>10) ;

      	if( fTree.jetsAK4_CorrPt[jid] > 20 &&
	    fabs( fTree.jetsAK4_Eta[jid] ) < 4.7 &&
	    looseid
	    ){
	  TLorentzVector jet;
	  jet.SetPtEtaPhiE( fTree.jetsAK4_CorrPt[jid] , fTree.jetsAK4_Eta[jid] , fTree.jetsAK4_Phi[jid] , fTree.jetsAK4_CorrE[jid] );
	  double DR = muon.DeltaR( jet );

	  if ( DR > 0.3 ){

	    if( fTree.jetsAK4_CorrPt[jid] <= 40 ){
	      nJetsPt20_47 ++ ;

//.........这里部分代码省略.........
		mlvjj = (lep + jet1 + jet2 + temp).M();
		if (mlvjj>masswindow_low&&mlvjj<masswindow_high)hlvjj3->Fill(mlvjj);
        pz_comparison3->Fill((pz-pnu->at(0)->pz())/pnu->at(0)->pz());
        lvjj_comparison3->Fill((mlvjj-lvjjmass_true)/lvjjmass_true);



		//4 ---------nu pz = far jet pz--------------------use DetaR to choose
		double dR1 = deltaR(lep.Eta(),lep.Phi(),jet1.Eta(),jet1.Phi());
		double dR2 = deltaR(lep.Eta(),lep.Phi(),jet2.Eta(),jet2.Phi());
		if(dR1>dR2)pz=jet1.Pz();
		else pz = jet2.Pz();
		temp.SetPxPyPzE(px,py,pz,sqrt(px*px+py*py+pz*pz));
		mlvjj = (lep + jet1 + jet2 + temp).M();
		if (mlvjj>masswindow_low&&mlvjj<masswindow_high)hlvjj4->Fill(mlvjj);
        pz_comparison4->Fill((pz-pnu->at(0)->pz())/pnu->at(0)->pz());
        lvjj_comparison4->Fill((mlvjj-lvjjmass_true)/lvjjmass_true);

*/
		// 5 ---------------------slove the function ------------milano

		float alpha = lep.Px()*pnu->at(0)->px() + lep.Py()*pnu->at(0)->py();

		float delta = (alpha + 0.5*80.399*80.399)*(alpha + 0.5*80.399*80.399) - lep.Pt()*lep.Pt()*pnu->at(0)->pt()*pnu->at(0)->pt();
		if( delta < 0. ) {
			delta = 0.;
		}
			//double wlv = (lep+neu).M();
			//hwlv->Fill(wlv);
			//cout<<"event with imaginary roots :"<<i<<" mwlv is  "<<wlv<<endl;
		
		//else continue;//use only imaginary root

		float pz1 = ( lep.Pz()*(alpha + 0.5*80.399*80.399) + lep.Energy()*sqrt(delta) ) / lep.Pt() / lep.Pt();
		float pz2 = ( lep.Pz()*(alpha + 0.5*80.399*80.399) - lep.Energy()*sqrt(delta) ) / lep.Pt() / lep.Pt();

		if( delta >= 0. )//choose the larger pz
		{
			if(fabs(pz1)<fabs(pz2))pz=pz2;
			else pz=pz1;   
		}
		//else continue;//discard unreal pz
		temp.SetPxPyPzE(px,py,pz,sqrt(px*px+py*py+pz*pz));
		mlvjj = (lep + jet1 + jet2 + temp).M();
		if (mlvjj>masswindow_low&&mlvjj<masswindow_high)hlvjj5->Fill(mlvjj);
        pz_comparison5->Fill((pz-pnu->at(0)->pz())/pnu->at(0)->pz());
        lvjj_comparison5->Fill((mlvjj-lvjjmass_true)/lvjjmass_true);


		//4 choose the smaller pz
		if( delta >= 0. )
        {   
            if(fabs(pz1)>fabs(pz2))pz=pz2;
            else pz=pz1;   
        }
        temp.SetPxPyPzE(px,py,pz,sqrt(px*px+py*py+pz*pz));
        mlvjj = (lep + jet1 + jet2 + temp).M();
        if (mlvjj>masswindow_low&&mlvjj<masswindow_high)hlvjj4->Fill(mlvjj);
        pz_comparison4->Fill((pz-pnu->at(0)->pz())/pnu->at(0)->pz());
        lvjj_comparison4->Fill((mlvjj-lvjjmass_true)/lvjjmass_true);

		//3 choose the one closest to lep pz
		
		if(fabs(pz1-lep.Pz())>fabs(pz2-lep.Pz()))pz=pz2;
		else pz=pz1;
        temp.SetPxPyPzE(px,py,pz,sqrt(px*px+py*py+pz*pz));

//.........这里部分代码省略.........
     if( isMC_ ) {
       
       // first look for W->qq'
       std::vector<TLorentzVector> quarksMC;

       for( unsigned iMc=0; iMc<nMc && quarksMC.size()<2; ++iMc ) {

         // quarks have status 3
         if( statusMc[iMc] != 3 ) continue;

         TLorentzVector* thisParticle = new TLorentzVector();
         thisParticle->SetPtEtaPhiE( pMc[iMc]*sin(thetaMc[iMc]), etaMc[iMc], phiMc[iMc], energyMc[iMc] );

         if( fabs(idMc[iMc])<7 && fabs(idMc[mothMc[iMc]])==24 ) {
           wIndexqq = mothMc[iMc];
           quarksMC.push_back( *thisParticle );
         }
       }

       // (checked that always 2 quarks are found)
       if( quarksMC.size()==2 && wIndexqq!=-1 ) {
	   eQuark1_=quarksMC[0].E();
	   ptQuark1_=quarksMC[0].Pt();
	   etaQuark1_=quarksMC[0].Eta();
	   phiQuark1_=quarksMC[0].Phi();
	   eQuark2_=quarksMC[1].E();
	   ptQuark2_=quarksMC[1].Pt();
	   etaQuark2_=quarksMC[1].Eta();
	   phiQuark2_=quarksMC[1].Phi();
 
         TLorentzVector WqqMC;
         WqqMC.SetPtEtaPhiE( pMc[wIndexqq]*sin(thetaMc[wIndexqq]), etaMc[wIndexqq], phiMc[wIndexqq], energyMc[wIndexqq] );
         ptWqqMC_  = WqqMC.Pt();
         eWqqMC_   = WqqMC.Energy();
         etaWqqMC_ = WqqMC.Eta();
         phiWqqMC_ = WqqMC.Phi();
	 //float deltaRqq = quarksMC[0].DeltaR(quarksMC[1]);
	 // h1_deltaRqq->Fill(deltaRqq);
       }
       
       // now look for W->lv

       std::vector<TLorentzVector> electronMC;
       std::vector<TLorentzVector> muonMC;
       std::vector<TLorentzVector> neutrinoMC;

       for( unsigned iMc=0; iMc<nMc; ++iMc ) {
	 
	 if( statusMc[iMc] != 3 ) continue;
	 
         TLorentzVector* thisParticle = new TLorentzVector();
         thisParticle->SetPtEtaPhiE( pMc[iMc]*sin(thetaMc[iMc]), etaMc[iMc], phiMc[iMc], energyMc[iMc] );

           if( fabs(idMc[iMc])==11 && fabs(idMc[mothMc[iMc]])==24 ) {electronMC.push_back( *thisParticle ); wIndexll = mothMc[iMc]; }
	   if( fabs(idMc[iMc])==13 && fabs(idMc[mothMc[iMc]])==24 ) {muonMC.push_back( *thisParticle ); wIndexll = mothMc[iMc]; }
	   if( (fabs(idMc[iMc])==12 || fabs(idMc[iMc])==14 ) && fabs(idMc[mothMc[iMc]])==24 ) {
	     neutrinoMC.push_back( *thisParticle );  wIndexll = mothMc[iMc];
	       }//for comparison with fit
 
         delete thisParticle;
         thisParticle = 0;
       }
       
       if( electronMC.size()==1 && neutrinoMC.size()==1 ) {
	 lept1MC = electronMC[0];
	 lept2MC = neutrinoMC[0];