1

私は C++ の専門家ではありませんが、.lhe ファイルからイベントを取得して、pp 衝突 (pp -> bbbar) で pythia による bbbar ペア質量不変分布を計算してプロットする簡単なプログラムを作成しようとしています。

ここに添付されたプログラムについて、次のように手順を実行しましたが、ヒストグラムには何も表示されません..ループに問題がある可能性があります..手順：

  • ヒストグラムを予約しました
  • イベントループを開始する
  • b 減衰を考慮して iBottom を定義 - 粒子ループを開始
  • defined b , bbar -パーティクル ループを開始します
  • ペアを見つけるために輪郭を作る
  • invM を計算します。

それについての助けはありますか？ピティアについて知っている人がいれば、そのようなプログラムの作成方法を知る..

// 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++ で実装する方法だけを尋ねるのが最善です。

手始めに、奇妙なセマンティクスにつながる不必要にネストされた多くのループに waaay があるように見えます (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 に答える