C++ Jet类代码示例

double BTagWeight::getUDSGScaleFactor(const Jet& jet) const {
	double pt = getSmearedJetPtScale(jet, 0)*jet.pt();
	double eta = jet.eta();
	double SF_udsg_mean(0), SF_udsg_min(0), SF_udsg_max(0);

	if (pt < 20) {
		return 0;
	} else if (pt > 800) {
		SF_udsg_mean = getMeanUDSGScaleFactor(800., eta);
		SF_udsg_min = getMinUDSGScaleFactor(800, eta);
		SF_udsg_max = getMaxUDSGScaleFactor(800, eta);
		//use twice the uncertainty
		SF_udsg_min -= (SF_udsg_mean - SF_udsg_min);
		SF_udsg_max += (SF_udsg_max - SF_udsg_mean);
	} else {
		SF_udsg_mean = getMeanUDSGScaleFactor(pt, eta);
		SF_udsg_min = getMinUDSGScaleFactor(pt, eta);
		SF_udsg_max = getMaxUDSGScaleFactor(pt, eta);
	}
	if (LightJetSystematic_ == -1)
		return SF_udsg_min;
	else if (LightJetSystematic_ == 1)
		return SF_udsg_max;

	return SF_udsg_mean;
}

void LHCOWriter::AnalyseJets()
{
  Jet *element;
  Track *track;
  Int_t counter;

  fItJet->Reset();
  while((element = static_cast<Jet*>(fItJet->Next())))
  {
    if(element->TauTag != 0) continue;

    Reset();

    counter = 0;
    fItTrack->Reset();
    while((track = static_cast<Track*>(fItTrack->Next())))
    {
      if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
    }

    fIntParam[1] = 4;

    fDblParam[0] = element->Eta;
    fDblParam[1] = element->Phi;
    fDblParam[2] = element->PT;
    fDblParam[3] = element->Mass;
    fDblParam[4] = counter;
    fDblParam[5] = element->BTag;
    fDblParam[6] = element->EhadOverEem;

    Write();
  }
}

void timer(int v) {
	moveBullets();
	jet.calcDir();
	jet.moveJet();
	placeCamera();
	glutPostRedisplay();
	// Vector3d axis = jet.up.cross(jet.forward);
	// debugLog(jet.forward, jet.up, axis);
	glutTimerFunc(50, timer, v);
}

Jet MuXboostedBTagging::Result(std::vector<Lepton> const &muons,  Jet const &jet, Jet &core) const
{
    INFO0;
    // Keep track of the smallest x for any muon
    auto min_x = std::numeric_limits<double>::max(); // Large enough, simple significand
    // Iterate through each muon, add it to the core, and calculate the resulting boost invariant
    auto p4_neutrino_correction = Jet {};
    for (auto const &muon : muons) {
        // Add back the muon and the neutrino
        core += 2 * muon;
        // Set a hard ceiling on subjet mass, for poor reconstruction
        auto dot = muon.Spatial().Dot(core.Spatial());
        auto min = std::min(core.Mass(), MaxSubJetMass());
        if (dot == 0_eV * eV || min == 0_eV) continue;
        auto cross = muon.Spatial().Cross(core.Spatial()).Mag();
        auto x_core = static_cast<double>(core.Energy() * cross / dot / min);
        CHECK(x_core >= 0, x_core);
        if (x_core < 0) return jet;
        // Add the neutrino to original jet IFF the muon passes the boosted test
        if (x_core <= XMax()) p4_neutrino_correction += muon;
        min_x = std::min(min_x, x_core);
    }

    auto result = jet + p4_neutrino_correction;
    CHECK(result.Pt() > 0_eV, result.Pt());
    CHECK(!(result.Pt() != result.Pt()), result.Pt())
    CHECK(!(core.Pt() != core.Pt()), core.Pt())
    CHECK(!(min_x != min_x), min_x)
    if (result.Pt() <= 0_eV) return jet;
    result.SetInfo(jet.Info());
    result.Info().SetMuBTag(min_x, core.Pt() / result.Pt());
    return result;
}

double BTagWeight::getUDSGScaleFactor(const Jet& jet, std::string MCSampleTag) const {
	double pt = getSmearedJetPtScale(jet, 0)*jet.pt();
	double eta = jet.eta();
	double SF_udsg_mean(0), SF_udsg_min(0), SF_udsg_max(0);

	if (MCSampleTag == "Summer12") { // 2012
		if (pt < 20) {
			return 0;
		} else if (pt > 850 && eta >= 1.6 && eta <= 2.4) {
			SF_udsg_mean = getMeanUDSGScaleFactor(850., eta, MCSampleTag);
			SF_udsg_min = getMinUDSGScaleFactor(850, eta, MCSampleTag);
			SF_udsg_max = getMaxUDSGScaleFactor(850, eta, MCSampleTag);
			//use twice the uncertainty
			SF_udsg_min -= (SF_udsg_mean - SF_udsg_min);
			SF_udsg_max += (SF_udsg_max - SF_udsg_mean);
		} else if (pt > 1000) {
			SF_udsg_mean = getMeanUDSGScaleFactor(1000., eta, MCSampleTag);
			SF_udsg_min = getMinUDSGScaleFactor(1000, eta, MCSampleTag);
			SF_udsg_max = getMaxUDSGScaleFactor(1000, eta, MCSampleTag);
			//use twice the uncertainty
			SF_udsg_min -= (SF_udsg_mean - SF_udsg_min);
			SF_udsg_max += (SF_udsg_max - SF_udsg_mean);
		} else {
			SF_udsg_mean = getMeanUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_min = getMinUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_max = getMaxUDSGScaleFactor(pt, eta, MCSampleTag);
		}
	} else if (MCSampleTag == "Summer11Leg") { // 2011
		if (pt < 20) {
			return 0;
		} else if (pt > 670 && eta >= 0. && eta <= 2.4) {
			// Use integrated over all eta
			SF_udsg_mean = getMeanUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_min = getMinUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_max = getMaxUDSGScaleFactor(pt, eta, MCSampleTag);
			//use twice the uncertainty
			SF_udsg_min -= (SF_udsg_mean - SF_udsg_min);
			SF_udsg_max += (SF_udsg_max - SF_udsg_mean);
		} else {
			SF_udsg_mean = getMeanUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_min = getMinUDSGScaleFactor(pt, eta, MCSampleTag);
			SF_udsg_max = getMaxUDSGScaleFactor(pt, eta, MCSampleTag);
		}
	}

	if (LightJetSystematic_ == -1) {
		return SF_udsg_min;
	} else if (LightJetSystematic_ == 1) {
		return SF_udsg_max;
	}
	return SF_udsg_mean;
}

int main()
{
	// Wyrażenie sin( (x + 2) * x). Rozwijam w x0 = 3.5; liczę wartość w 4.0
	Jet<double> x ( 5);
	x.set(0, 3.5);
	x.set(1, 1.0);
	
	Jet<double> c2 ( 5);
	c2.set(0, 2.0);
	
	
	cout<<"wynik: "<<getTaylor(sin( (x+c2) * x) , x.get(0), 4.0)<<endl;
	cout<<"wynik oczekiwany: "<< sin( (4.0 + 2) * 4.0) <<endl;
	
	// Wyrażenie pow( 2.0/x + 2, 2.72). Rozwijam w x0 = 3.5; liczę wartość w 4.0
	cout<<"wynik: "<<getTaylor(pow( c2/x + c2, 2.72) , x.get(0), 4.0)<<endl;
	cout<<"wynik oczekiwany: "<<  pow( 2.0/4.0 + 2.0, 2.72) <<endl;
	
	// Wyrażenie log( x/(x+2) +2 + sin(cos(x))) . Rozwijam w x0 = 3.5; liczę wartość w 4.0
	cout<<"wynik: "<<getTaylor(log(x/(x+c2) +c2 + sin(cos(x))) , x.get(0), 4.0)<<endl;
	cout<<"wynik oczekiwany: "<<  log(4.0/(4.0+2.0) +2.0 + sin(cos(4.0))) <<endl;
	
	
	
	// Wyrażenie sin(y) całuję od 0 do 2 PI. 
	Jet<double> y ( 21);
	y.set(0, 0.0);
	y.set(1, 1.0);
	
	cout<<"wynik: "<< integrate(sin(y),y.get(0), 0.0, 2.0 * M_PI, 200, 34)<<endl;
	cout<<"oczekiwnany: 0"<<endl;
	
	return 0;
}

void fillPFCandVsCount ( Jet & jet, TH2F* hist ){
  
  const int n_ids = 8;
  int pfcand_ids[n_ids + 1];
  memset ( pfcand_ids, 0, sizeof (pfcand_ids));
  
  int nPFCands = jet.getNPFCandidates();
  for (int iPFCand = 0; iPFCand < nPFCands; ++iPFCand){
    PFCand cand = jet.getPFCandidate( iPFCand );
    pfcand_ids[cand.Id()]++;
    pfcand_ids[n_ids]++;
  }
  for (int i = 0; i <= n_ids; ++i) hist -> Fill(i, pfcand_ids[i]);

}

std::vector<Jet> MuXboostedBTagging::CoreCandidates(std::vector<Lepton> const &muons, Jet const &jet) const
{
    INFO0;
    auto recluster_input = muons;
    for (auto const &consituent : jet.Constituents()) {
        if (consituent.Info().ContainsDetectorPart(DetectorPart::tower) && consituent.Pt() < min_tower_pt_ratio_ * jet.Pt()) continue; // Don't use
        recluster_input.emplace_back(consituent);
    }
    // Recluster the jet, to find core candidates
    auto cluster_sequence = boca::ClusterSequence {recluster_input, fastjet::JetDefinition(fastjet::antikt_algorithm, core_jet_radius / rad, &Settings::Recombiner())};
    // Get the core candidates (NOT sorted by pT until we remove muons)
    auto core_candidates = cluster_sequence.InclusiveJets();
    if (core_candidates.empty()) return core_candidates;
    cluster_sequence.NoLongerNeeded();

    // If a taggable muon is inside a core candidate, remove the muon p4
    for (auto const &muon : muons) for (auto &core_candidate : core_candidates) if (muon.DeltaRTo(core_candidate) < core_jet_radius) {
                core_candidate -= muon;
                break;
            }

    // Sort the core candidates by pT (highest to lowest)
    core_candidates = SortedByPt(core_candidates);
    return core_candidates;
}

Type getTaylor(const Jet<Type> &jet, const Type &x0, const Type &x, int n = -1)
{
	Type out = 0;
	Type powerAcc = 1;
	
	if(n == -1) n = jet.length();
	
	for(int i = 0; i < n; i++)
	{
		out += jet.get(i) * powerAcc;
		powerAcc *= (x - x0);
		
	}
	
	return out;
}

void display() {

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

	glColor3ub(0,0,0);

	glPushMatrix();
	glScaled(100, .05, .05);
	DrawCube(0, 0, 0, 4, true);
	glPopMatrix();

	glPushMatrix();
	glScaled(.05, 100, .05);
	DrawCube(0, 0, 0, 4, true);
	glPopMatrix();

	glPushMatrix();
	glScaled(.05, .05, 100);
	DrawCube(0, 0, 0, 4, true);
	glPopMatrix();

	paintGridAround(400);

	//DrawCube(0, 0, 0, 4, true);

	paintBullets();

	jet.paintJet();

	glutSwapBuffers();
}

bool Analysis::DphiCuts(){
	bool IsCutPass = false; 
	if (branchJet && branchMissingET){
	  MissingET *met = (MissingET*) branchMissingET->At(0);	
	  int nJets = branchJet->GetEntriesFast();	
	  for (int i=0; i < nJets; i++){
	    Jet *jet = (Jet*) branchJet->At(i);	
		if (jet->PT > 200.){
			if (jet->P4().DeltaR(met->P4()) > 1.) 
				IsCutPass=true; 
			else {IsCutPass=false; break;}
		}
	  }
	}
	IsEventPass = IsEventPass && IsCutPass;
	return IsCutPass;
}

Jet* JetTruthEval::get_truth_jet(PHG4Particle* particle) {

  if (_strict) {assert(particle);}
  else if (!particle) {++_errors; return NULL;}
  
  if (_do_cache) {
    std::map<PHG4Particle*,Jet*>::iterator iter =
      _cache_get_truth_jet.find(particle);
    if (iter != _cache_get_truth_jet.end()) {
      return iter->second;
    }
  }

  Jet* truth_jet = NULL;
  
  // loop over all jets and look for this particle...
  for (JetMap::Iter iter = _truthjets->begin();
       iter != _truthjets->end();
       ++iter) {
    Jet* candidate = iter->second;

    // loop over all consituents and look for this particle    
    for (Jet::ConstIter jter = candidate->begin_comp();
	 jter != candidate->end_comp();
	 ++jter) {
      unsigned int index = jter->second;
      
      PHG4Particle* constituent = _truthinfo->GetParticle( index );
      if (_strict) {assert(constituent);}
      else if (!constituent) {++_errors; continue;}

      if (get_svtx_eval_stack()->get_truth_eval()->are_same_particle(constituent,particle)) {
	truth_jet = candidate;
	break;
      }
    }
    
    if (truth_jet) break;
  }  
  
  if (_do_cache) _cache_get_truth_jet.insert(make_pair(particle,truth_jet));
  
  return truth_jet;
}

double JetInfo::CoreEnergyFraction(Jet const &jet) const
{
    DEBUG0;
    auto energy = 0_eV;
    auto core_energy = 0_eV;
    for (auto const &constituent : Constituents()) if (constituent.DetectorPart() == DetectorPart::photon) {
            energy += constituent.Momentum().Et();
            if (jet.DeltaRTo(constituent.Momentum()) < 0.2_rad) core_energy += constituent.Momentum().Et();
        }
    return energy == 0_eV ? 0. : static_cast<double>(core_energy / energy);
}

Angle JetInfo::TrackRadius(Jet const &jet) const
{
    DEBUG0;
    auto energy = 0_eV;
    auto weight = 0_eV * rad;
    for (auto const &constituent : Constituents()) if (constituent.DetectorPart() == DetectorPart::track) {
            energy += constituent.Momentum().Et();
            weight += constituent.Momentum().Et() * jet.DeltaRTo(constituent.Momentum());
        }
    return energy == 0_eV ? 0_rad : Angle(weight / energy);
}

Jet smearedJet(const Jet & origJet, unsigned mode) {
	if (origJet.getVarF("jn_genpt") <= 0)
		return origJet;

	//smearing factors are described in https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
	double eta = fabs(origJet.lorentzVector().Eta());
	double pt = origJet.lorentzVector().Pt();
	double ptSF = 1.0;
	double ptSF_err = 0.06;
	if (eta < 0.5) {
		ptSF = 1.052;
		ptSF_err = sqrt(pow(0.012, 2) + pow(0.5 * (0.062 + 0.061), 2));
	} else if (eta >= 0.5 && eta < 1.1) {
		ptSF = 1.057;
		ptSF_err = sqrt(pow(0.012, 2) + pow(0.5 * (0.056 + 0.055), 2));
	} else if (eta >= 1.1 && eta < 1.7) {
		ptSF = 1.096;
		ptSF_err = sqrt(pow(0.017, 2) + pow(0.5 * (0.063 + 0.062), 2));
	} else if (eta >= 1.7 && eta < 2.3) {
		ptSF = 1.134;
		ptSF_err = sqrt(pow(0.035, 2) + pow(0.5 * (0.087 + 0.085),2));
	} else if (eta >= 2.3 && eta < 5.0) {
		ptSF = 1.288;
		ptSF_err = sqrt(pow(0.127, 2) + pow(0.5 * (0.155 + 0.153), 2));
	}

	if (mode == 1)
		ptSF += ptSF_err;
	else if (mode == 2)
		ptSF -= ptSF_err;
	
	ptSF = max(0., (origJet.getVarF("jn_genpt") + ptSF * (pt - origJet.getVarF("jn_genpt")))) / pt;  //deterministic version
	if (ptSF <= 0)
		return origJet;

	double px = origJet.getVarF("jn_px") * ptSF;
	double py = origJet.getVarF("jn_py") * ptSF;
	double pz = origJet.getVarF("jn_pz");
	double mass = origJet.lorentzVector().M();
	double en = sqrt(mass * mass + px * px + py * py + pz * pz);

	Jet smearedJet(origJet);
	smearedJet.addFloatVar( "jn_px", px );
	smearedJet.addFloatVar( "jn_py", py );
	smearedJet.addFloatVar( "jn_pz", pz );
	smearedJet.addFloatVar( "jn_en", en );
	return smearedJet;
}