Hi!
I'm not sure if my last mail reached the roottalk list. Anyway here it
comes again with some further information about the problem.
I'm using a TObjArray to host a list of 4-vectors (particles) that I run
a jet algorithm on.
The main feature of the main program, and the jet algorithm can be found
in attachement.
The vector list is passed inside the jet-finder, where there is a copy
of the vector list, that is being deleted using the member function
Delete();
After a few events, something called garbage collector sweeps the floor,
and my program crashes with a segmentation fault. Can someone help me
with this?
cheers!
/mattias
--
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Mattias Davidsson, http://www.quark.lu.se/~mattias
+46 (0)46 222 76 83
Department of elementary particle physics
Box 118 SE-221 00
Lund, Sweden
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#include <stdlib.h>
#include <iostream>
#include <cmath>
#include <fstream>
#include <list>
#include <string>
#include <vector>
#include "Pythia.h"
#include "LorentzVector.h"
#include "JetFindDict.h"
//
// This program produces DIS events from pythia, adds the stable particles
// to a TObjArray (veclist). The aim is then to feed these 4-vectors
// to the jetfinder package. See bottom for troubles.
//
TROOT demo("demo","my first Root program");
int main() {
const char* frame="CMS";
int iframe=3;
const char* beam="gamma/e+";
int ibeam=8;
const char* target="p+";
int itarget=2;
double E=297.2E0;
double y_cut=0.1E0;
// Uncomment below for e+ e-
//const char* frame="CMS";
//int iframe=3;
//const char* beam="e+";
//int ibeam=2;
//const char* target="e-";
//int itarget=2;
//double E=91.2E0;
//double y_cut=0.E0;
// Choose jet finding scheme. Jade with y_cut
pydat1_.mstu[46-1]=4; // MSTU(46)=4
pydat1_.paru[45-1]=y_cut; // PARU(45)=y_cut
// Disable preclustering and don't ignore neutrinos
pydat1_.paru[43-1]=0.E0; // PARU(43)=0
pydat1_.mstu[41-1]=1; // MSTU(41)=1
// Setup ROOT jet finder
JadeJetFinder *jetf = new JadeJetFinder(0.0005);
//DurhamJetFinder *jetf = new DurhamJetFinder(0.004); // Durham algorithm
// IktJetFinder *jetf = new IktJetFinder(0.01,'I'); // Inclusive kt algorithm
// Select DIS process in PYTHIA.
pysubs_.msel=1; // MSEL=0
pysubs_.msub[1-1]=1; // MSUB(1)=1
//Initialize Pythia
pyinit_(frame,beam,target,&E,iframe,ibeam,itarget);
TObjArray* veclist = new TObjArray();
// Main event loop
int events=20;
for (long iev=0;iev<events;iev++){
// Declare 4vector list a'la ROOT
if ( veclist ) veclist->Delete();
// if (veclist) veclist->Clear();
//Declare the 4 vector a'la ROOT
TLorentzVector* vec4 = new TLorentzVector();
if(vec4) delete vec4;
// Declare a 3 vector a'la ROOT
TVector3* vec3 = new TVector3();
if(vec3) delete vec3;
// why do I get a segmentation fault when declaring these here???
// TObjArray* veclist = new TObjArray();
// if(veclist) veclist->Delete();
// TVector3* vec3 = new TVector3();
// TLorentzVector* vec4 = new TLorentzVector();
// if(vec3) delete vec3;
// if(vec4) delete vec4;
Float_t Pvec3[3];
pyevnt_();
long nlist=2;
if (iev < 3) pylist_(&nlist);
int N=pyjets_.n;
for(int i=1; i<=N; i++) {
if (pyjets_.k[1-1][i-1]==1 && pyjets_.k[2-1][i-1]!=22) { // Stable particle (not photon)
if (pyjets_.k[2-1][i-1]!=-11 && pyjets_.k[2-1][i-1]!=2212) { // not the electron or proton
Pvec3[0] = pyjets_.p[1-1][i-1];
Pvec3[0] = pyjets_.p[2-1][i-1];
Pvec3[0] = pyjets_.p[3-1][i-1];
Float_t E_part = pyjets_.p[4-1][i-1];
vec3 = new TVector3(Pvec3);
vec4 = new TLorentzVector(*vec3,E_part);
veclist->Add(vec4);
delete vec3;
}
}
}
// Perform Jet finding
jetf->setEvent(veclist);
jetf->setYCut(0.0003);
jetf->doFindJets();
std::cout << "No'f jets: " << jetf->njets() << std::endl;
//int Njet=0;
//pyclus_(&Njet);
//std::cout << "Pythia: " << Njet << std::endl;
// print jet quantities
//delete veclist;
// if (veclist) {veclist->Delete();
//delete veclist;}
// event->Clear();
// event->Clear();
std::cout << "even #: " << iev << std::endl;
}
delete veclist;
}
// Jet finder base class routines.
// Author: Rob Shanks (derived from Java routines
// written by Gary Bower)
// Version 0.1 Mar 01/99
// Version 0.2 Jul 02/99 : M.Iwasaki Make some Modification on
// masscut and ymass calculation (in JadeEJetfonder.cxx
// or DurhamJetFinder.cxx) on doFindJets().
// Adding Jade cluster (yclus).
// Version 0.3 Jul 07/99 : M.Iwasaki Change Mod to Mag function in TVector3
// so as to use root 2.22.x
// Version 0.4 Jul 21/99 : M.Iwasaki Make temporary 4-vectors copied
// from input particle 4-vectors in doFindJets()
// to protect the initial 4-vectors.
// Version 0.4 Sep 23/99 : M.Iwasaki Fix setEvent memory leak.
//
#include <iostream>
#include "JetFinder.h"
ClassImp(JetFinder)
const Int_t JetFinder::UNASSOC = -999;
JetFinder::JetFinder(double ycut, char name)
: m_resultsValid(0), m_evis(0),m_dycut(ycut), m_name(name)
{}
JetFinder::JetFinder(double ycut)
: m_resultsValid(0), m_name('X'), m_evis(0),m_dycut(ycut)
{}
//______________________________________________________
JetFinder::~JetFinder() {
m_jet.Delete();
m_ParticlesInJet.Delete();
m_4vec->Delete();
part->Delete();
delete m_4vec;
delete part;
}
//______________________________________________________
Int_t JetFinder::njets()
{
if (!m_resultsValid) doFindJets();
return m_injets;
}
//______________________________________________________
TLorentzVector* JetFinder::jet(Int_t index)
{
if (!m_resultsValid) doFindJets();
return (TLorentzVector*)m_jet[index];
}
//_______________________________________________________
Int_t JetFinder::nParticlesPerJet(Int_t index)
{
if (!m_resultsValid) doFindJets();
return m_inparts_per_jet[index];
}
//_______________________________________________________
Int_t JetFinder::fewestTracks()
{
if (!m_resultsValid) doFindJets();
return m_ifewest_tracks;
}
// MD set the name of the algorithm
void JetFinder::setName(char name)
{
m_name = name;
}
// MD get the name of the algorithm
char JetFinder::getName()
{
return m_name;
}
//_______________________________________________________
/* Obtain the current ycut
* @return The current value of ycut */
double JetFinder::getYCut()
{
return m_dycut;
}
//_______________________________________________________
/* Set the YCut value.
* If the new value for ycut is not the same as the old ycut the
* jet finding will be rerun.
* @param ycut the new value to be set for ycut */
void JetFinder::setYCut(double ycut)
{
if (m_dycut != ycut) m_resultsValid = 0;
m_dycut = ycut;
}
//_______________________________________________________
/* Call this to input a new event to the jet finder.
* Only elements of the enumeration which are accepted
* by the predicate will be used for jet finding.
* @param e: An Enumeration of either TLorentzVector(px,py,pz,E)
* or TVector3(px,py,pz) */
void JetFinder::setEvent(TObjArray* e)
{
m_resultsValid = 0;
if(m_4vec) delete m_4vec;
m_4vec = new TObjArray();
m_jet.Delete();
m_ParticlesInJet.Delete();
m_evis = 0;
Int_t ne = e->GetEntries();
for(Int_t i=0;i<ne;i++)
{
TObject* o = e->At(i);
TString* nam = new TString(o->IsA()->GetName());
if (nam->Contains("TLorentzVector"))
{
TLorentzVector* in = (TLorentzVector*) o;
m_evis += in->T();
m_4vec->Add(in);
}
else if (nam->Contains("TVector3"))
{
TVector3* in = (TVector3*) o;
m_evis += in->Mag();
m_4vec->Add(new TLorentzVector(*in,in->Mag()));
}
else printf("JetFinder::setEvent input is not a TVector3 or a TLorentzVector\n");
delete nam;
}
}
//_______________________________________________________
void JetFinder::doFindJets()
{
// MD is the masscut variable put in the right place?
// it has to be found in the LOOP => put it as member function.
minmass = 10000000;
masscut = m_dycut * m_evis * m_evis ;
m_resultsValid = 1;
m_injets = 0;
np = m_4vec->GetEntries();
if (np<2) return;
m_ipart_jet_assoc.Set(np);
for (Int_t m=0; m<np; m++) m_ipart_jet_assoc[m] = UNASSOC;
if(part) {
part->Delete();
delete part;
}
part = new TObjArray();
if(ymass) {
delete ymass;
}
ymass = new TMatrix(np,np);
double Part_temp[3], Part_temp_4;
for (Int_t Ipart=0; Ipart <np; Ipart++)
{
Part_temp[0] = ((TLorentzVector *)(m_4vec->At(Ipart)))->X();
Part_temp[1] = ((TLorentzVector *)(m_4vec->At(Ipart)))->Y();
Part_temp[2] = ((TLorentzVector *)(m_4vec->At(Ipart)))->Z();
Part_temp_4 = ((TLorentzVector *)(m_4vec->At(Ipart)))->T();
TVector3* vec3 = new TVector3(Part_temp);
TLorentzVector* vec4 = new TLorentzVector(*vec3,Part_temp_4);
part->Add(vec4);
delete vec3;
}
std::cout << "cow 2: np " << np << std::endl;
doMatrix(*part, *ymass);
//std::cout << "cow 3: " << std::endl;
Loop(*ymass, *part);
std::cout << "cow 3: " << std::endl;
Finish(*part);
}
//_______________________________________________________
void JetFinder::doMatrix(TObjArray & particle,TMatrix & distance)
{
// create invariant mass pair array.
// This can later be global (in JetFider.cxx)
for (Int_t i = 0; i < np - 1; i++ )
{
for (Int_t j = i + 1 ; j < np ; j++ )
{
distance(i,j) =
calcinvmass((TLorentzVector*)particle.At(i),
(TLorentzVector*)particle.At(j));
}
}
}
//_______________________________________________________
void JetFinder::Loop(TMatrix &distance, TObjArray &particle)
{
for (;;)
{
im = -1;
jm = -1;
std::cout << "x:" << std::endl;
if (MinMass(distance) < masscut)
{
Merge(particle, distance);
}
std::cout <<"MinMass: " <<MinMass(distance) << std::endl;
std::cout << "masscut: "<< masscut << std::endl;
if (MinMass(distance) >= masscut) break;
}
}
//_______________________________________________________
double JetFinder::MinMass(const TMatrix &distance)
{
// find least invariant mass pair.
//
minmass = 10000000;
for(Int_t i = 0 ; i < np - 1 ; i++ )
{
for(Int_t j = i + 1 ; j < np ; j++ )
{
if (m_ipart_jet_assoc[i] != JetFinder::UNASSOC) continue;
if (m_ipart_jet_assoc[j] != JetFinder::UNASSOC) continue;
if (distance(i,j) > minmass) continue;
minmass = distance(i,j);
im = i;
jm = j;
}
}
return minmass;
}
//____________________________________________________________
void JetFinder::Merge(TObjArray &particle, TMatrix &distance)
{
//
// combine particles im and jm.
// this is simple 4-vector addidtion
*(TLorentzVector*)particle.At(im) += *(TLorentzVector*)particle.At(jm);
for(Int_t i = 0; i < np; i++ )
{
if( m_ipart_jet_assoc[i] == jm )
{
m_ipart_jet_assoc[i] = im;
}
}
// Recalculate invariant masses for newly combined particle
//
m_ipart_jet_assoc[jm] = im;
for( Int_t i = 0; i < np ; i++ )
{
if( i == im) continue;
if (m_ipart_jet_assoc[i] != UNASSOC ) continue;
Int_t imin = TMath::Min(i,im);
Int_t imax = TMath::Max(i,im);
distance(imin,imax) =
calcinvmass((TLorentzVector*)particle.At(imin),
(TLorentzVector*)particle.At(imax));
}
}
//____________________________________________________________
void JetFinder::Finish(const TObjArray &particle)
{
// finish up by filling jet array.
for( Int_t i = 0 ; i < np ; i++ )
{
if (m_ipart_jet_assoc[i] == UNASSOC) m_injets++;
}
m_inparts_per_jet.Set(m_injets);
Int_t nj = 0;
Int_t ntrk;
m_ifewest_tracks = 100000; // Starting min value
for( Int_t i = 0 ; i < np ; i++ )
{
if (m_ipart_jet_assoc[i] != UNASSOC) continue;
m_jet.Add((TLorentzVector*)particle.At(i));
ntrk = 1;
for (Int_t j = 0 ; j < np ; j++)
{
if(m_ipart_jet_assoc[j] == i)
{
m_ipart_jet_assoc[j] = nj;
ntrk++;
}
}
m_ipart_jet_assoc[i] = nj;
m_inparts_per_jet[nj] = ntrk;
if( ntrk < m_ifewest_tracks) m_ifewest_tracks = ntrk;
nj++;
}
}
//_______________________________________________________
// MD this is used to add a pseudo particle to account for
// missing pz (p-remnant)
//
void JetFinder::PseudoParticle()
{
Int_t np = m_4vec->GetEntries();
TLorentzVector* psp = new TLorentzVector;
double misspz = 0.;
for(Int_t i = 1 ; i < np ; i++ )
{
misspz += ((TLorentzVector*)m_4vec->At(i))->Pz();
}
psp->SetZ(misspz);
m_4vec->Add(psp);
}
This archive was generated by hypermail 2b29 : Tue Jan 02 2001 - 11:50:21 MET