1

我在 C++ 方面不是那么专家，但我尝试制作一个简单的程序来计算和绘制 pythia 在 pp 碰撞（pp -> bbbar）中的 bbbar 对质量不变分布，从 .lhe 文件中获取事件。

现在关于此处附加的程序，我按照以下步骤进行了操作，但直方图中没有显示任何内容..可能是循环中有问题..步骤：

  • 预订直方图
  • 启动事件循环
  • 定义 iBottom 以解释 b 衰减 - 启动粒子循环
  • 定义 b , bbar - 启动粒子循环
  • 制作轮廓以找到对
  • 计算invM。

那么有什么帮助吗？如果有人知道pythia，就知道如何制作这样的程序..

// The program

#include <iostream>
// Header file to access Pythia 8 program elements.
#include "Pythia.h"
// ROOT, for histogramming.
#include "TH1.h"
// ROOT, for interactive graphics.
#include "TVirtualPad.h"

{#include "TApplication.h"
// ROOT, for saving file.
#include "TFile.h"
using namespace Pythia8;
int main(int argc, char* argv[]) {

// Create the ROOT application environment.
TApplication theApp("hist", &argc, argv);

// create Pythia object and set up generation
Pythia pythia;
pythia.readString("PartonLevel:MI = off");
pythia.init("events.lhe"); 
pythia.readString("Beams:eCM = 14000.");
pythia.readString("PhaseSpace:pTHatMin = 20.");
pythia.init();

// Create file on which histograms can be saved.
TFile* outFile = new TFile("hist.root", "RECREATE");

TH1F *Mbbbar = new TH1F("MbB","MbB invariant mass", 0,0.,2000.);


// Begin event loop. Genera1te event; skip if generation aborted.
for (int iEvent = 0; iEvent<100; ++iEvent)
{
int iBottom = 0;

if(!pythia.next()) continue;
for (int i = 0; i< pythia.event.size(); ++i)
{
if (pythia.event[i].id() == 5) iBottom = i;

// Look for BQuarks among decay products.
      vector<int> b, bbar;
      for (int i = 0; i< pythia.event.size(); ++i) {
        int id = pythia.event[i].id();  
        if (id ==  5) b.push_back(i);
        if (id == -5) bbar.push_back(i);



// Check whether pair(s) present.
      int N = bbar.size();
      int n = b.size();
      if (N + n > 1) {

// Fill masses of BQuarks pair.
for (int i1 = 0; i1 < n - 1; ++i1)
       for (int i2 = 0; i2 < N - 1; ++i2)
          Mbbbar->Fill(
            (pythia.event[b[i1]].p() +pythia.event[bbar[i2]].p()).mCalc() );

}
}
}
}

pythia.stat();
Mbbbar->Draw();
std::cout << "\nDouble click on the histogram window to quit.\n";
gPad->WaitPrimitive();
// Save histogram on file and close file.
Mbbbar->Write();
delete outFile;

// Done
return 0;
}

谢谢你，萨菲纳兹

4

1 回答 1

0

您的问题似乎既涉及特定领域的知识，也涉及您在此处难以找到的编程知识。如果您尝试用基本术语解释您的意图，并且只询问如何在 C++ 中实现这些语义，那将是最好的。

对于初学者来说，看起来你有很多不必要的嵌套循环，这会导致奇怪的语义（你也应该小心声明多个int i同名的循环）。我试图猜测你想要做什么并用评论稍微解释一下

// Begin event loop. Genera1te event; skip if generation aborted.
int numEvent = 0;
int maxEvents = 100;
//as long as pytha.next() generates something and we have processed less than 100 events
while(pythia.next() && numEvent++ < maxEvents){
    vector<int> b, bbar;

    //for each particle in the event
    for (int i = 0; i < pythia.event.size(); ++i){
        //add the particle index to b if the id() is 5
        //otherwise add it to bbar
        switch(pythia.event[i].id()){
            case 5: 
                b.push_back(i);
                break;
            case -5: 
                bbar.push_back(i);
                break;
        }
    }

    int bSize  = b.size();
    int bbSize = bbar.size();
    //for every possible pair of id()=5 and id()=-5 particles found during this event
    for (int i1 = 0; i1 < bSize; ++i1)
        for (int i2 = 0; i2 < bbSize; ++i2)
            Mbbbar->Fill(
            (pythia.event[b[i1]].p() + pythia.event[bbar[i2]].p()).mCalc());
}

此外，您应该单步执行并检查您认为已执行的代码是否已执行。如果由于某种原因无法在调试模式下执行此操作，您还可以打印调试消息，例如：

std::cout << "starting event loop" << std::endl;
while(pythia.next() && numEvent++ < maxEvents){
    std::cout << "starting event cycle number " << numEvent << " with " << pythia.event.size() << " particles" << std::endl;
    //...

    std::cout << "found " << b.size() << " particles with id()=5" << std::endl;
    std::cout << "found " << bbar.size() << " particles with id()=-5" << std::endl; 
    for (int i1 = 0; i1 < bSize; ++i1)
    //...
于 2013-09-10T16:20:00.883 回答