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