segmentation fault and garbage collector

From: Mattias Davidsson (mattias@quark.lu.se)
Date: Tue Mar 07 2000 - 13:51:57 MET


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