Segmentation Violation in simulation-particles-code

ROOT
1

Hi,
first of all best wishes for your holidays.
My program has to simulate the generation of different types of particles. The main program is called Mainmodule and it needs of ParticleType.c+, ResonanceType.c+, Particle.c+ .

When I execute my macro Mainmodule() it gives me Segmentation Violation Error.

I would to know what could be the problem, where I can find it. I think it would be helpful to have the codes of ParticleType, ResonanceType and Particle.
I specify that .L ParticleType, .L ResonanceType, .L Particle and .L Mainmodule works with no error!

===========================================================
There was a crash (kSigSegmentationViolation).
This is the entire stack trace of all threads:
===========================================================
#0  0xb775ece5 in __kernel_vsyscall ()
#1  0xb70d3123 in __waitpid_nocancel () at ../sysdeps/unix/syscall-template.S:84
#2  0xb705d9d6 in do_system (line=0x90b1e38 "/home/nirvana/Immagini/ROOT/etc/gdb-backtrace.sh 6762 1>&2") at ../sysdeps/posix/system.c:148
#3  0xb76180e2 in TUnixSystem::Exec (shellcmd=<optimized out>, this=<optimized out>) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:2118
#4  TUnixSystem::StackTrace (this=0x8417698) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:2412
#5  0xb761aa2b in TUnixSystem::DispatchSignals (this=0x8417698, sig=kSigSegmentationViolation) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:3643
#6  0xb761ab88 in SigHandler (sig=kSigSegmentationViolation) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:408
#7  0xb76141b3 in sighandler (sig=11) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:3620
#8  0xb760cfc6 in textinput::TerminalConfigUnix::HandleSignal (this=0xb772c280 <textinput::TerminalConfigUnix::Get()::s>, signum=11) at /home/nirvana/Scaricati/root-6.10.02/core/textinput/src/textinput/TerminalConfigUnix.cpp:100
#9  0xb760cff1 in (anonymous namespace)::TerminalConfigUnix__handleSignal (signum=11) at /home/nirvana/Scaricati/root-6.10.02/core/textinput/src/textinput/TerminalConfigUnix.cpp:36
#10 <signal handler called>
#11 0xb737f5a4 in ParticleType::GetName() const () from /home/nirvana/Immagini/ROOT/macros/ParticleType_c.so
#12 0xb737518e in Particle::FindParticle(char const*) () from /home/nirvana/Immagini/ROOT/macros/Particle_c.so
#13 0xb7375213 in Particle::AddParticleType(char const*, double, int, double) () from /home/nirvana/Immagini/ROOT/macros/Particle_c.so
#14 0xb736eb71 in Mainmodule() () from /home/nirvana/Immagini/ROOT/macros/Mainmodule_C.so
#15 0xa515c083 in ?? ()
#16 0xb42b3d44 in cling::Interpreter::RunFunction(clang::FunctionDecl const*, cling::Value*) () from /home/nirvana/Immagini/ROOT/lib/libCling.so
#17 0xb42b4ee6 in cling::Interpreter::EvaluateInternal(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::CompilationOptions, cling::Value*, cling::Transaction**, unsigned int) () from /home/nirvana/Immagini/ROOT/lib/libCling.so
#18 0xb42b511b in cling::Interpreter::process(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::Value*, cling::Transaction**, bool) () from /home/nirvana/Immagini/ROOT/lib/libCling.so
#19 0xb433ded3 in cling::MetaProcessor::process(char const*, cling::Interpreter::CompilationResult&, cling::Value*, bool) () from /home/nirvana/Immagini/ROOT/lib/libCling.so
#20 0xb42605a0 in HandleInterpreterException (metaProcessor=<optimized out>, input_line=<optimized out>, compRes=
0xbfe7d3b0: (unknown: 144375800), result=0xbfe7d3b4) at /home/nirvana/Scaricati/root-6.10.02/core/metacling/src/TCling.cxx:1887
#21 0xb42733fb in TCling::ProcessLine (this=0x8451a08, line=0x9013740 "#line 1 \"ROOT_prompt_6\"\nMainmodule()", error=0xbfe7d6b8) at /home/nirvana/Scaricati/root-6.10.02/core/metacling/src/TCling.cxx:2050
#22 0xb7526a68 in TApplication::ProcessLine (this=0x8448018, line=<optimized out>, sync=false, err=0xbfe7d6b8) at /home/nirvana/Scaricati/root-6.10.02/core/base/src/TApplication.cxx:1001
#23 0xb773eb1c in TRint::ProcessLineNr (this=0x8448018, filestem=0xb77519a3 "ROOT_prompt_", line=0x89b28d0 "Mainmodule()", error=0xbfe7d6b8) at /home/nirvana/Scaricati/root-6.10.02/core/rint/src/TRint.cxx:741
#24 0xb773ee6d in TRint::HandleTermInput (this=0x8448018) at /home/nirvana/Scaricati/root-6.10.02/core/rint/src/TRint.cxx:602
#25 0xb773e31e in TTermInputHandler::Notify (this=0x89af960) at /home/nirvana/Scaricati/root-6.10.02/core/rint/src/TRint.cxx:131
#26 0xb7740e90 in TTermInputHandler::ReadNotify (this=0x89af960) at /home/nirvana/Scaricati/root-6.10.02/core/rint/src/TRint.cxx:123
#27 0xb7619db5 in TUnixSystem::CheckDescriptors (this=0x8417698) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:1321
#28 0xb761b42c in TUnixSystem::DispatchOneEvent (this=0x8417698, pendingOnly=false) at /home/nirvana/Scaricati/root-6.10.02/core/unix/src/TUnixSystem.cxx:1076
#29 0xb7500998 in TSystem::InnerLoop (this=0x8417698) at /home/nirvana/Scaricati/root-6.10.02/core/base/src/TSystem.cxx:410
#30 0xb74fedad in TSystem::Run (this=0x8417698) at /home/nirvana/Scaricati/root-6.10.02/core/base/src/TSystem.cxx:360
#31 0xb752317e in TApplication::Run (this=0x8448018, retrn=false) at /home/nirvana/Scaricati/root-6.10.02/core/base/src/TApplication.cxx:1153
#32 0xb7740742 in TRint::Run (this=0x8448018, retrn=false) at /home/nirvana/Scaricati/root-6.10.02/core/rint/src/TRint.cxx:455
#33 0x08048939 in main (argc=1, argv=0xbfe7fb14) at /home/nirvana/Scaricati/root-6.10.02/main/src/rmain.cxx:30

These are the other codes:

ParticleType.h

#ifndef PARTICLETYPE_H
#define PARTICLETYPE_H


class ParticleType {

public:
ParticleType (const char* fName_, double fMass_, int fCharge_);
const char* GetName () const;
double GetMass ()const ;
int GetCharge ()const;
double GetWidth  () const;
void Print();

protected:

private:
const char* fName;
const double fMass;
const int fCharge;

};

ParticleType.c+

#include <iostream>
#include "ParticleType.h"
using namespace std; 

ParticleType::ParticleType(const char* fName_, double fMass_, int fCharge_) : fName (fName_), fMass (fMass_), fCharge (fCharge_) {} 

const char* ParticleType::GetName() const {return fName;}
double ParticleType::GetMass() const {return fMass;}
int ParticleType::GetCharge() const {return fCharge;}
double ParticleType::GetWidth () const {return 0 ;}
void ParticleType::Print() {cout<<fName<<'\n'<<fMass<<'\n'<<fCharge<<'\n';}

ResonanceType.h

#include "ParticleType.h"
#ifndef RESONANCETYPE_H 
#define RESONANCETYPE_H


class ResonanceType : public ParticleType {

public:
ResonanceType (const char* fName_, double fMass_, int fCharge_, double fWidth_);
double GetWidth ();
void Print ();

protected:

private:
const double fWidth;


};

#endif

ResonanceType.c+

#include <iostream>
#include "ParticleType.h"
#include "ResonanceType.h"

using namespace std; 

ResonanceType::ResonanceType (const char* fName_, double fMass_, int fCharge_, double fWidth_) : ParticleType(fName_, fMass_, fCharge_), fWidth (fWidth_) {}

void ResonanceType::Print () {ParticleType::Print(); cout<<fWidth<<'\n';}
double ResonanceType::GetWidth() {return fWidth;}

Particle.h

#include "ParticleType.h"
#include "ResonanceType.h"
#include <cmath>
#ifndef PARTICLE_H
#define PARTICLE_H
//far ereditare non :: gli attributi massa charge e width
class Particle {

public:
Particle();
Particle(const char *fName_, double fPx_, double fPy_, double fPz_);
double Energy() const;
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
int GetFIndex();
static void AddParticleType (const char *fName_, double fMass_, int fCharge_, double fWidth_);
void SetFIndex (int a);
static void PrintArray ();
void PrintOther ();
double GetfPx () const;
double GetfPy () const;
double GetfPz () const;
double GetNewMass () const;
double GetEnergy () const;
double InvMass (Particle &p);
void SetP(double px, double py, double pz);
int Decay2Body(Particle &dau1, Particle &dau2) const;

protected: 

private:
static ParticleType* fParticleType [10];
static const int fMaxNumParticleType = 10;
static int fNParticleType;
int fIndex;
double fPx = 0;
double fPy = 0;
double fPz = 0;
static int FindParticle(const char *fName_);
void Boost (double bx, double by, double bz);
int a;

};


#endif

Particle.c+

#include <iostream>
#include "ResonanceType.h"
#include "ParticleType.h"
#include "Particle.h"
#include <cmath>
#include <cstdlib> 

using namespace std;
int Particle::fNParticleType = 10;
ParticleType* Particle:: fParticleType[] = {0x0};

Particle::Particle(){}
Particle::Particle (const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 

double Particle::GetfPx () const  
{
 return fPx;
}

double Particle::GetfPy () const
{
 return fPy;
}

double Particle::GetfPz () const
{
 return fPz;
}

double Particle::GetNewMass ()  const
{
 return fParticleType[fIndex]->GetMass();
}


void Particle::SetFIndex (int a) 
{
 int FIndex=a;
}

int Particle::GetFIndex () 
{
 return fIndex;
}

//void Particle::SetFIndex (char*fName) {  //!

//;}

void Particle::SetP(double px, double py, double pz) 
{
 double fPx = px; 
 double fPy = py; 
 double fPz = pz;
}


int Particle::FindParticle (const char *fName_)
{
 for (int i=0; i<fNParticleType; i++) 
 { 
   if (fName_ == fParticleType[i]->GetName())
   {
     return i;
   }
   else 
   {
     cout<<"Non c'è corrispondenza"<<'\n';
   } 
 }
return -1;
}


void Particle::AddParticleType (const char *fName_, double fMass_, int fCharge_, double fWidth_)
{
 for (int i=0; i<fNParticleType; i++) 
 {
   int b = FindParticle(fName_);
   if (b == -1) 
   {
     if (fWidth_ == 0) 
     {
       fParticleType[i]=new ParticleType (fName_, fMass_, fCharge_);
     }
     else 
     {
       fParticleType[i]=new ResonanceType (fName_, fMass_, fCharge_, fWidth_);
     }
   }
 }
}


void Particle::PrintArray() 
{
 for(int i=0; i<fNParticleType; i++)
 {
   fParticleType[i]->Print();
 }
}  

void Particle::PrintOther () 
{
cout<<"Name"<<fParticleType[fIndex]->GetName()<<'\n'<<"fPx"<<fPx<<'\n'<<"fPy"<<fPy<<'\n'<<"fPz"<<fPz;
}


double Particle::Energy() const
{ 
 double Pulse=pow(GetfPx(), 2)+pow(GetfPy(), 2)+pow(GetfPz(), 2);
 double Masstwo= GetNewMass()*GetNewMass();
 double Energy=sqrt(Masstwo+Pulse);
  
  return Energy;
} 


double Particle::InvMass(Particle &p)
{ 
  double InvMass=(pow(pow(Energy()+p.Energy(),2)-pow(GetfPx()+p.GetfPx(),2)+pow(GetfPy()+p.GetfPy(),2)+pow(GetfPz()+p.GetfPz(),2), 1/2));
  return InvMass;
} 


int Particle::Decay2Body(Particle &dau1, Particle &dau2) const {
  if(GetNewMass() == 0.0){
    printf("Decayment cannot be preformed if mass is zero\n");
    return 1;
  }
  
  double massMot = GetNewMass();
  double massDau1 = dau1.GetNewMass();
  double massDau2 = dau2.GetNewMass();

  if(fIndex > -1){ // add width effect

    // gaussian random numbers

    float x1, x2, w, y1, y2;
    
    double invnum = 1./RAND_MAX;
    do {
      x1 = 2.0 * rand()*invnum - 1.0;
      x2 = 2.0 * rand()*invnum - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
    
    w = sqrt( (-2.0 * log( w ) ) / w );
    y1 = x1 * w;
    y2 = x2 * w;

    massMot += fParticleType[fIndex]->GetWidth() * y1;

  }

  if(massMot < massDau1 + massDau2){
    printf("Decayment cannot be preformed because mass is too low in this channel\n");
    return 2;
  }
  
  double pout = sqrt((massMot*massMot - (massDau1+massDau2)*(massDau1+massDau2))*(massMot*massMot - (massDau1-massDau2)*(massDau1-massDau2)))/massMot*0.5;

  double norm = 2*M_PI/RAND_MAX;

  double phi = rand()*norm;
  double theta = rand()*norm*0.5 - M_PI/2.;
  dau1.SetP(pout*sin(theta)*cos(phi),pout*sin(theta)*sin(phi),pout*cos(theta));
  dau2.SetP(-pout*sin(theta)*cos(phi),-pout*sin(theta)*sin(phi),-pout*cos(theta));

  double energy = sqrt(fPx*fPx + fPy*fPy + fPz*fPz + massMot*massMot);

  double bx = fPx/energy;
  double by = fPy/energy;
  double bz = fPz/energy;

  dau1.Boost(bx,by,bz);
  dau2.Boost(bx,by,bz);

  return 0;
}
void Particle::Boost(double bx, double by, double bz)
{

  double energy = GetEnergy();

  //Boost this Lorentz vector
  double b2 = bx*bx + by*by + bz*bz;
  double gamma = 1.0 / sqrt(1.0 - b2);
  double bp = bx*fPx + by*fPy + bz*fPz;
  double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

  fPx += gamma2*bp*bx + gamma*bx*energy;
  fPy += gamma2*bp*by + gamma*by*energy;
  fPz += gamma2*bp*bz + gamma*bz*energy;
}

Mainmodule.c+

#include <iostream>
#include "ParticleType.h"
#include "ResonanceType.h"
#include "Particle.h"
#include "cmath"
#include "TMath.h"
#include "TRandom.h"
#include "TH1D.h"
#include "TFile.h"
#include "TCanvas.h"

using namespace std; 

int Mainmodule() 

{ 

double Px, Py, Pz;
double Energia, Impulso_Trasverso, Impulso;
double MassaInv;

Particle particella [110]; //100 particelle per volta, sistema e riempie, poi lo fa per 10^5 volte

Particle::AddParticleType ("pioni", 0.13957, 1, 0);
Particle::AddParticleType ("Kaoni", 0.49367, 1, 0);
Particle::AddParticleType ("protoni", 0.93827, 1, 0);
Particle::AddParticleType ("K*", 0.89166, 0, 0.050);
Particle::AddParticleType ("pioni", 0.13957, -1, 0);
Particle::AddParticleType ("Kaoni", 0.49367, -1, 0);
Particle::AddParticleType ("protoni", 0.93827, -1, 0);
Double_t x = 0;
Double_t y = 0;
Double_t phi = 0;
Double_t theta = 0;

TFile *file=new TFile("Particelle.root", "RECREATE");

TH1D *HTypes=new TH1D ("HTypes", "Types", 7, 0, 110);
TH1D *HAzimuth=new TH1D ("HAzimuth", "Azimuth_Angle_Distribution", 100, 0, 2*TMath::Pi());
TH1D *HPolar=new TH1D ("HPolar", "Polar_Angle_Distribution", 100, 0, TMath::Pi());
TH1D *HPulse=new TH1D ("HPulse", "Pulse_Distribution", 100, 0, 1000);
TH1D *HTrasverso_Pulse=new TH1D ("HTrasverso_Pulse", "Trasverso_Pulse_Distribution", 100, 0, 1000);
TH1D *HEnergy=new TH1D ("HEnergy", "Energy_Distribution", 100, 0, 1000);

TH1D *HInvEverybody=new TH1D ("HInvEverybody", "InvMass_Everybody", 500, 0, 1000);
TH1D *HInvOpposite_Sign=new TH1D ("HInvOpposite_Sign", "InvMass_Opposite_Sign", 500, 0, 1000);
TH1D *HInvSame_Sign=new TH1D ("HInvSame_Sign", "InvMass_Same_Sign", 500, 0, 1000);
TH1D *HInvPione_concorde=new TH1D ("HInvPione_concorde", "InvMass_Pioneplus/Kaoneminus", 500, 0, 1000);
TH1D *HInvPione_discorde=new TH1D ("HInvPione_discorde", "InvMass_Pioneminus/Kaoneplus", 500, 0, 1000);

TH1D *HInvK=new TH1D ("HInvK", "InvMass_K*", 500, 0, 1000);

 for (Int_t i=0; i<10; i++)
  {
   int n=0;
   for (Int_t j=0; j<110; j++)
   { 
	cout<<"asdasd"<<endl;
     x=gRandom->Rndm();
     y=gRandom->Rndm();
     phi=gRandom->Uniform(0, 2*TMath::Pi());
     theta=gRandom->Uniform(0, TMath::Pi());
     Impulso=gRandom->Exp(1.0);

     if (x<0.01) 
     {
       if (y<0.50) 
       {
          particella[j].SetFIndex(1);
       }
     }
     
     if (x<0.10) 
     {
          if (y<0.50) {particella[j].SetFIndex(2);}
       else {particella[j].SetFIndex(3);}
     }

     if (x<0.20) 
         {if (y<0.50) {particella[j].SetFIndex(4);}
        else {particella[j].SetFIndex(5);}}

     if (x<1.00)  
        {if (y<0.50) {particella[j].SetFIndex(6);}
        else {particella[j].SetFIndex(7);} }


HTypes->Fill(particella[j].GetFIndex()); //Rimepio Types

     Px=Impulso*sin(theta)*cos(phi);
     Py=Impulso*sin(theta)*sin(phi);
     Pz=Impulso*cos(theta);

HAzimuth->Fill(phi);
 //Rimepio Azimuth Angle

HPolar->Fill(theta);
 //Rimpio Polar Angle

HPulse->Fill(Impulso);
 // Rimepio Pulse

particella[j].SetP(Px,Py,Pz); //Impulso per ogni particella

Impulso_Trasverso = pow((pow(Px, 2)+pow(Py, 2)),1/2);
HTrasverso_Pulse->Fill(Impulso_Trasverso);
//Riempio Trasverso Pulse

Energia = particella[j].GetEnergy();
HEnergy->Fill(Energia);
 //Rimepio Energy 

double InvMass;

     if (x<0.01)
     {
       particella[n].Decay2Body(particella[n+1], particella[n+2]);
       MassaInv=particella[n+1].InvMass(particella[n+2]);
       HInvK->Fill(MassaInv);
       //Rimepio K*
       n+=2;
       if (x<0.005) 
       {
         particella[100+n].SetFIndex(6);
         particella[100+n+1].SetFIndex(5);

       
       }
        else
       {
         particella[100+n].SetFIndex(7);
         particella[100+n+1].SetFIndex(4);
       }
        particella[j].Decay2Body(particella[100+n], particella[100+n+1]);
        MassaInv=particella[100+n].InvMass(particella[100+n+1]);
        HInvK->Fill(MassaInv);
     }
   }

   for(Int_t k=0; k<110; k++)
   {
     for (Int_t h=0; h<110; h++)
     {
      MassaInv=particella[k].InvMass(particella[h]);

      HInvEverybody->Fill(MassaInv);
//Rimepio Everybody


      if( (particella[k].GetFIndex()==4 || particella[k].GetFIndex()==2 || particella[k].GetFIndex()==6) &&
      (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==2 || particella[h].GetFIndex()==6)) 
      {
       HInvSame_Sign->Fill(MassaInv);
      } 

       if((particella[k].GetFIndex()==7 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==3) &&
       (particella[h].GetFIndex()==7 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==3))        
        {
          HInvSame_Sign->Fill(MassaInv);

        }
     
       //concorde

       if((particella[k].GetFIndex()==4 || particella[k].GetFIndex()==2 || particella[k].GetFIndex()==6) &&
       (particella[h].GetFIndex()==7 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==3)) 
       {
        HInvOpposite_Sign->Fill(MassaInv);
       } 

       if((particella[k].GetFIndex()==7 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==3) && 
       (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==2 || particella[h].GetFIndex()==6))       
        {
          HInvOpposite_Sign->Fill(MassaInv);
        }
       
       //discorde

       if((particella[k].GetFIndex()==6) && (particella[h].GetFIndex()==4) || (particella[k].GetFIndex()==7) && (particella[h].GetFIndex()==5) ||    (particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==6) || (particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==7))
       {
        HInvPione_concorde->Fill(MassaInv);

       }
       //pionekaoneconcorde

       if((particella[k].GetFIndex()==6) && (particella[h].GetFIndex()==5) || (particella[k].GetFIndex()==7) && (particella[h].GetFIndex()==4) ||    (particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==6) || (particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==7))
       {
        HInvPione_discorde->Fill(MassaInv);
       }
       //pionekaonediscorde

}
}
}
 file->Write();
 file->Close();

return 0;
}

Thanks very much.

2

Rename your “*.c” files into “*.cxx”.
Then try (fix all reported bugs and warnings):

.L ParticleType.cxx++g
.L ResonanceType.cxx++g
.L Particle.cxx++g
.L Mainmodule.cxx++g
Mainmodule()
3

Hi,
Thanks for the answer.
I’ve resolved the error but I’ve a new problem. I have modified something in the codes. So I report them above. Ignore the strategic cout, please.
The problem is this:

when I execute Mainmodule() the generation of 100 elements (1st for cycle) is done, but at line 187 of Mainmodule() it stops and gives
Error in HandleInterpreterException>: Trying to dereference null pointer or trying to call routine taking non-null arguments.
Execution of your code was aborted.

Mainmodule.cxx

#include <iostream>
#include "ParticleType.h"
#include "ResonanceType.h"
#include "Particle.h"
#include "cmath"
#include "TMath.h"
#include "TRandom.h"
#include "TH1D.h"
#include "TFile.h"
#include "TCanvas.h"

using namespace std; 

int Mainmodule() 

{ 

double Px, Py, Pz;
double Energia, Impulso_Trasverso, Impulso;
double MassaInv;

Particle particella [110]; //100 particelle per volta, sistema e riempie, poi lo fa per 10^5 volte
cout<<"L'errore non è prima degli add";

Particle::AddParticleType ("pioni", 0.13957, 1, 0);
Particle::AddParticleType ("Kaoni", 0.49367, 1, 0);
Particle::AddParticleType ("protoni", 0.93827, 1, 0);
Particle::AddParticleType ("K*", 0.89166, 0, 0.050);
Particle::AddParticleType ("pioni", 0.13957, -1, 0);
Particle::AddParticleType ("Kaoni", 0.49367, -1, 0);
Particle::AddParticleType ("protoni", 0.93827, -1, 0);

Double_t x = 0;
Double_t y = 0;
Double_t phi = 0;
Double_t theta = 0;
cout<<"L'errore non è negli add";

TFile *file=new TFile("Particelle.root", "RECREATE");

TH1D *HTypes=new TH1D ("HTypes", "Types", 7, 0, 110);
TH1D *HAzimuth=new TH1D ("HAzimuth", "Azimuth_Angle_Distribution", 100, 0, 2*TMath::Pi());
TH1D *HPolar=new TH1D ("HPolar", "Polar_Angle_Distribution", 100, 0, TMath::Pi());
TH1D *HPulse=new TH1D ("HPulse", "Pulse_Distribution", 100, 0, 1000);
TH1D *HTrasverso_Pulse=new TH1D ("HTrasverso_Pulse", "Trasverso_Pulse_Distribution", 100, 0, 1000);
TH1D *HEnergy=new TH1D ("HEnergy", "Energy_Distribution", 100, 0, 1000);
cout<<"L'errore non è nella prima parte degli histo";
TH1D *HInvEverybody=new TH1D ("HInvEverybody", "InvMass_Everybody", 500, 0, 1000);
TH1D *HInvOpposite_Sign=new TH1D ("HInvOpposite_Sign", "InvMass_Opposite_Sign", 500, 0, 1000);
TH1D *HInvSame_Sign=new TH1D ("HInvSame_Sign", "InvMass_Same_Sign", 500, 0, 1000);
TH1D *HInvPione_concorde=new TH1D ("HInvPione_concorde", "InvMass_Pioneplus/Kaoneminus", 500, 0, 1000);
TH1D *HInvPione_discorde=new TH1D ("HInvPione_discorde", "InvMass_Pioneminus/Kaoneplus", 500, 0, 1000);
TH1D *HInvK=new TH1D ("HInvK", "InvMass_K*", 500, 0, 1000);
cout<<"L'errore non è negli istrogrammi";


 for (Int_t i=0; i<10e5; i++)
  {
   int n=0;
   for (Int_t j=0; j<100; j++)
   { 
	cout<<"asdasd"<<endl;
     x=gRandom->Rndm(); cout<<"Ho generato x";
     y=gRandom->Rndm(); cout<<"Ho generato y";
     phi=gRandom->Uniform(0, 2*TMath::Pi()); cout<<"Ho generato phi";
     theta=gRandom->Uniform(0, TMath::Pi()); cout<<"Ho genrato theta";
     Impulso=gRandom->Exp(1.0); cout<<"Ho generato Impulso";

     if (x<0.01) 
     {
          particella[j].SetFIndex(1);
          cout<<"Sono K*";
     }
     else 
     {
      
      if (x<0.10) 
      {
          if (y<0.50) 
         {
           particella[j].SetFIndex(2); cout<<"sono protone+";
         }
          else 
         {
           particella[j].SetFIndex(3); cout<<"sono protone -";
         }
      }
       else 
      {

        if (x<0.20) 
        {
          if (y<0.50) 
         {
           particella[j].SetFIndex(4); cout<<"sono Kaone +";
         }
          else 
         {
           particella[j].SetFIndex(5); cout<<"sono Kaone -";
         }
        }
        else 
        {

          if (x<1.00)  
          {
            if (y<0.50) 
            {
             particella[j].SetFIndex(6); cout<<"sono pione +";
            }
            else 
            {
             particella[j].SetFIndex(7); cout<<"sono pione -";
             } 
          }
         }
        }
       }
      
   

HTypes->Fill(particella[j].GetFIndex()); //Rimepio Types

     Px=Impulso*sin(theta)*cos(phi);
     Py=Impulso*sin(theta)*sin(phi);
     Pz=Impulso*cos(theta);

HAzimuth->Fill(phi);
 //Rimepio Azimuth Angle

HPolar->Fill(theta);
 //Rimpio Polar Angle

HPulse->Fill(Impulso);
 // Rimepio Pulse
cout<<"Ho riempito qualche istogramma";

particella[j].SetP(Px,Py,Pz); //Impulso per ogni particella

cout<<"Roba impulso";

Impulso_Trasverso = pow((pow(Px, 2)+pow(Py, 2)),1/2);
HTrasverso_Pulse->Fill(Impulso_Trasverso);
cout<<"Rimepito impulso";
/*Riempio Trasverso Pulse

Energia = particella[j].Energy(); //SI BLOCCA QUI 
cout<<"Assegno l'energia";
HEnergy->Fill(Energia);
cout<<"rimepita l'energia";
 Rimepio Energy 
*/cout<<"ho finito di riempire istogrammi";

double InvMass;

    if (x<0.01)
     {
       if (x<0.5) 
       { cout<<"sono k* e devo decadere";
         particella[100+n].SetFIndex(6); cout<<"ho settato il primo indice";
         particella[100+n+1].SetFIndex(5); cout<<"ho  settato il seocndo indice";
         particella[j].Decay2Body(particella[100+n], particella[100+n+1]); cout<<"sono decaduta";
      
       }
        else
       { cout<<"sono k* e devo decadere";
         particella[100+n].SetFIndex(7); cout<<"ho settato il primo indice";
         particella[100+n+1].SetFIndex(4); cout<<"ho  settato il seocndo indice";
        particella[j].Decay2Body(particella[100+n], particella[100+n+1]); cout<<"sono decaduta";
       }

        MassaInv=particella[100+n].InvMass(particella[100+n+1]);
        HInvK->Fill(MassaInv);
       n+=2;
       cout<<"sono";
     }
   
} 


   for(Int_t k=0; k<100+n; k++)
   { 
     cout<<"Prima delle masse invarianti ci sono";
     for (Int_t h=0; h<100+n; h++)
     {
      cout<<"ti sto dando una massa invariante";
      MassaInv=particella[k].InvMass(particella[h]); //SI BLOCCA QUI 

      HInvEverybody->Fill(MassaInv); cout<<"riempio 0";
//Rimepio Everybody


      if(( particella[k].GetFIndex()==4 || particella[k].GetFIndex()==2 || particella[k].GetFIndex()==6) && (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==2 || particella[h].GetFIndex()==6)) 
      {
       HInvSame_Sign->Fill(MassaInv); cout<<"riempio 1";
      } 

      if((particella[k].GetFIndex()==7 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==3) &&
      (particella[h].GetFIndex()==7 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==3))     
      {
       HInvSame_Sign->Fill(MassaInv); cout<<"riempio 2";

      }
     
       //concorde

       if((particella[k].GetFIndex()==4 || particella[k].GetFIndex()==2 || particella[k].GetFIndex()==6) &&
       (particella[h].GetFIndex()==7 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==3) )
       {
        HInvOpposite_Sign->Fill(MassaInv); cout<<"riempio 3";
       } 

       if((particella[k].GetFIndex()==7 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==3) && 
       (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==2 || particella[h].GetFIndex()==6))       
        {
          HInvOpposite_Sign->Fill(MassaInv); cout<<"riempio 4";
        }
       
       //discorde

       if(((particella[k].GetFIndex()==6) && (particella[h].GetFIndex()==4)) || ((particella[k].GetFIndex()==7) && (particella[h].GetFIndex()==5)) ||  ((particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==6)) || ((particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==7)))
       {
        HInvPione_concorde->Fill(MassaInv); cout<<"riempio 5";

       }
       //pionekaoneconcorde

       if(((particella[k].GetFIndex()==6) && (particella[h].GetFIndex()==5)) || ((particella[k].GetFIndex()==7) && (particella[h].GetFIndex()==4)) ||  ((particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==6)) || ((particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==7)))
       {
        HInvPione_discorde->Fill(MassaInv); cout<<"riempio  6";
       }
       //pionekaonediscorde

}
cout<<"Sono alla "<< i <<"esima iterazione";
}
}
 file->Write();
 file->Close();

return 0;
cout<<"Fine generazione";
}

Particle.cxx (where methods are defined)

#include <iostream>
#include "ResonanceType.h"
#include "ParticleType.h"
#include "Particle.h"
#include <cmath>
#include <cstdlib> 

using namespace std;


int Particle::fNParticleType =0;
ParticleType*Particle:: fParticleType[] = {0x0};
const int Particle::fMaxNumParticleType;

Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
Particle::Particle():fPx(0), fPy(0), fIndex(-1){}

double Particle::GetfPx() const  
{
 return fPx;
}

double Particle::GetfPy() const
{
 return fPy;
}

double Particle::GetfPz() const
{
 return fPz;
}

double Particle::GetNewMass()  const
{
 return fParticleType[fIndex]->GetMass(); //MI DICE CHE E' UN PUNTATORE NULLO
}


void Particle::SetFIndex(int a) 
{
 int fIndex=a;
}

int Particle::GetFIndex() 
{
 return fIndex;
}

//void Particle::SetFIndex (char*fName) {  //!

//;}

void Particle::SetP(double px, double py, double pz) 
{
 double fPx = px; 
 double fPy = py; 
 double fPz = pz;
}



int Particle::FindParticle(const char *fName_)
{ cout<<"Sto eseguendo FindParticle";
 for (int i=0; i<fNParticleType; i++) 
 { 
   if (fName_ == fParticleType[i]->GetName())
   { cout<<"Ecco l'indice corrispondente";
     return i;
   }
   else 
   {
     cout<<"Non c'è corrispondenza"<<'\n';
   } 
 }
return -1;
}


void Particle::AddParticleType(const char *fName_, double fMass_, int fCharge_, double fWidth_)
{ cout<<"Sto eseguendo Addparticle";
 for (int i=0; i<fNParticleType; i++) 
 {
   int b = FindParticle(fName_);
   if (b == -1) 
   {
     if (fWidth_ == 0) 
     {
       fParticleType[i]=new ParticleType (fName_, fMass_, fCharge_); cout<<"Ecco una nuova ParticleType";
     }
     else 
     {
       fParticleType[i]=new ResonanceType (fName_, fMass_, fCharge_, fWidth_); cout<<"Ecco una nuova ResonanceType";
     }
   }
 }
}



void Particle::PrintArray() 
{
 for(int i=0; i<fNParticleType; i++)
 {
   fParticleType[i]->Print();
 }
}  

void Particle::PrintOther() 
{
cout<<"Name"<<fParticleType[fIndex]->GetName()<<'\n'<<"fPx"<<fPx<<'\n'<<"fPy"<<fPy<<'\n'<<"fPz"<<fPz;
}


double Particle::TotalEnergy() const
{ 
 double Pulse=pow(GetfPx(), 2)+pow(GetfPy(), 2)+pow(GetfPz(), 2);
 double Masstwo= GetNewMass()*GetNewMass();
 double Energy=sqrt(Masstwo+Pulse);
  
  return Energy;
} 


double Particle::InvMass(Particle &p)
{ 
  double InvMass=pow(pow(TotalEnergy()+p.TotalEnergy(),2)-(pow(GetfPx()+p.GetfPx(),2)+pow(GetfPy()+p.GetfPy(),2)+pow(GetfPz()+p.GetfPz(),2)), 1/2);
  return InvMass;
} 


int Particle::Decay2Body(Particle &dau1, Particle &dau2) const 
 {
  if(GetNewMass() == 0.0){
    cout<<"Decayment cannot be preformed if mass is zero\n";
    return 1;
  }
  
  double massMot = GetNewMass();
  double massDau1 = dau1.GetNewMass();
  double massDau2 = dau2.GetNewMass();

  if(fIndex > -1){ // add width effect

    // gaussian random numbers

    float x1, x2, w, y1, y2;
    
    double invnum = 1./RAND_MAX;
    do {
      x1 = 2.0 * rand()*invnum - 1.0;
      x2 = 2.0 * rand()*invnum - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
    
    w = sqrt( (-2.0 * log( w ) ) / w );
    y1 = x1 * w;
    y2 = x2 * w;

    massMot += fParticleType[fIndex]->GetWidth() * y1;

  }

  if(massMot < massDau1 + massDau2){
    cout<<"Decayment cannot be preformed because mass is too low in this channel\n";
    return 2;
  }
  
  double pout = sqrt((massMot*massMot - (massDau1+massDau2)*(massDau1+massDau2))*(massMot*massMot - (massDau1-massDau2)*(massDau1-massDau2)))/massMot*0.5;

  double norm = 2*M_PI/RAND_MAX;

  double phi = rand()*norm;
  double theta = rand()*norm*0.5 - M_PI/2.;
  dau1.SetP(pout*sin(theta)*cos(phi),pout*sin(theta)*sin(phi),pout*cos(theta));
  dau2.SetP(-pout*sin(theta)*cos(phi),-pout*sin(theta)*sin(phi),-pout*cos(theta));

  double energy = sqrt(fPx*fPx + fPy*fPy + fPz*fPz + massMot*massMot);

  double bx = fPx/energy;
  double by = fPy/energy;
  double bz = fPz/energy;

  dau1.Boost(bx,by,bz);
  dau2.Boost(bx,by,bz);

  return 0;
}



void Particle::Boost(double bx, double by, double bz)
{

  double energy = TotalEnergy();

  //Boost this Lorentz vector
  double b2 = bx*bx + by*by + bz*bz;
  double gamma = 1.0 / sqrt(1.0 - b2);
  double bp = bx*fPx + by*fPy + bz*fPz;
  double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

  fPx += gamma2*bp*bx + gamma*bx*energy;
  fPy += gamma2*bp*by + gamma*by*energy;
  fPz += gamma2*bp*bz + gamma*bz*energy;
}

Thanks very much

4

I guess you modified “Particle.h”, too.

5

Hi,
yes. I post here all the codes because now I have Segmentation Violation Error after many thousand of generation.
Here’s the stack

===========================================================
There was a crash (kSigSegmentationViolation).
This is the entire stack trace of all threads:
===========================================================
#0  0xb7736d05 in __kernel_vsyscall ()
#1  0xb70a3213 in __waitpid_nocancel () at ../sysdeps/unix/syscall-template.S:84
#2  0xb702bc56 in do_system (line=0x9b86890 "/home/nirvana/Scrivania/root-6.08.06/etc/gdb-backtrace.sh 2885 1>&2") at ../sysdeps/posix/system.c:148
#3  0xb75f4e38 in TUnixSystem::Exec (shellcmd=<optimized out>, this=0x8eb4698) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:2118
#4  TUnixSystem::StackTrace (this=<optimized out>) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:2405
#5  0xb75f7b3b in TUnixSystem::DispatchSignals (this=0x8eb4698, sig=kSigSegmentationViolation) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:3625
#6  0xb75f7c98 in SigHandler (sig=kSigSegmentationViolation) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:408
#7  0xb75ec831 in (anonymous namespace)::TerminalConfigUnix__handleSignal (signum=11) at /home/nirvana/Scaricati/root-6.08.06/core/textinput/src/textinput/TerminalConfigUnix.cpp:36
#8  <signal handler called>
#9  0xacf3bbe8 in ?? ()
#10 0xacf31840 in ?? ()
#11 0xacf3502e in ?? ()
#12 0xb454312c in cling::Interpreter::RunFunction(clang::FunctionDecl const*, cling::Value*) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#13 0xb4544cc6 in cling::Interpreter::EvaluateInternal(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::CompilationOptions, cling::Value*, cling::Transaction**, unsigned int) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#14 0xb454501a in cling::Interpreter::process(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::Value*, cling::Transaction**) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#15 0xb45e27ba in cling::MetaProcessor::process(char const*, cling::Interpreter::CompilationResult&, cling::Value*) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#16 0xb449d741 in HandleInterpreterException (metaProcessor=<optimized out>, input_line=<optimized out>, compRes=
0xbf8e6b9c: (unknown: 150229336), result=0xbf8e6ba4) at /home/nirvana/Scaricati/root-6.08.06/core/meta/src/TCling.cxx:1883
#17 0xb44b19b0 in TCling::ProcessLine (this=0x8eee8e8, line=0x960edb0 "#line 1 \"ROOT_prompt_7\"\nMainmodule()", error=0xbf8e6ea8) at /home/nirvana/Scaricati/root-6.08.06/core/meta/src/TCling.cxx:2049
#18 0xb74ca69e in TApplication::ProcessLine (this=0x8ee5260, line=<optimized out>, sync=false, err=<optimized out>) at /home/nirvana/Scaricati/root-6.08.06/core/base/src/TApplication.cxx:1005
#19 0xb77177e3 in TRint::ProcessLineNr (this=0x8ee5260, filestem=0xb7728e96 "ROOT_prompt_", line=0x95bf9c0 "Mainmodule()", error=0xbf8e6ea8) at /home/nirvana/Scaricati/root-6.08.06/core/rint/src/TRint.cxx:749
#20 0xb7717b90 in TRint::HandleTermInput (this=<optimized out>) at /home/nirvana/Scaricati/root-6.08.06/core/rint/src/TRint.cxx:610
#21 0xb75f6e6d in TUnixSystem::CheckDescriptors (this=0x8eb4698) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:1321
#22 0xb75f851c in TUnixSystem::DispatchOneEvent (this=<optimized out>, pendingOnly=<optimized out>) at /home/nirvana/Scaricati/root-6.08.06/core/unix/src/TUnixSystem.cxx:1076
#23 0xb752c118 in TSystem::InnerLoop (this=0x8eb4698) at /home/nirvana/Scaricati/root-6.08.06/core/base/src/TSystem.cxx:408
#24 0xb752a490 in TSystem::Run (this=<optimized out>) at /home/nirvana/Scaricati/root-6.08.06/core/base/src/TSystem.cxx:358
#25 0xb74c72ce in TApplication::Run (this=0x8ee5260, retrn=false) at /home/nirvana/Scaricati/root-6.08.06/core/base/src/TApplication.cxx:1157
#26 0xb771954e in TRint::Run (this=<optimized out>, retrn=false) at /home/nirvana/Scaricati/root-6.08.06/core/rint/src/TRint.cxx:463
#27 0x08048cdf in main (argc=<optimized out>, argv=0xbf8e92c4) at /home/nirvana/Scaricati/root-6.08.06/main/src/rmain.cxx:30
===========================================================


The lines below might hint at the cause of the crash.
You may get help by asking at the ROOT forum http://root.cern.ch/forum.
Only if you are really convinced it is a bug in ROOT then please submit a
report at http://root.cern.ch/bugs. Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#9  0xacf3bbe8 in ?? ()
#10 0xacf31840 in ?? ()
#11 0xacf3502e in ?? ()
#12 0xb454312c in cling::Interpreter::RunFunction(clang::FunctionDecl const*, cling::Value*) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#13 0xb4544cc6 in cling::Interpreter::EvaluateInternal(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::CompilationOptions, cling::Value*, cling::Transaction**, unsigned int) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#14 0xb454501a in cling::Interpreter::process(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, cling::Value*, cling::Transaction**) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#15 0xb45e27ba in cling::MetaProcessor::process(char const*, cling::Interpreter::CompilationResult&, cling::Value*) () from /home/nirvana/Scrivania/root-6.08.06/lib/libCling.so
#16 0xb449d741 in HandleInterpreterException (metaProcessor=<optimized out>, input_line=<optimized out>, compRes=
0xbf8e6b9c: (unknown: 150229336), result=0xbf8e6ba4) at /home/nirvana/Scaricati/root-6.08.06/core/meta/src/TCling.cxx:1883
===========================================================

ParticleType.h and .cxx and ResonanceType.h and .cxx are the same.

Particle.h

#include <iostream>
#include "ParticleType.h"
#include "ResonanceType.h"
#include <cmath>

#ifndef PARTICLE_H
#define PARTICLE_H


class Particle {

public:
Particle();
Particle(const char *fName_, double fPx_, double fPy_, double fPz_);
double TotalEnergy() const;
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
int GetFIndex();
static void AddParticleType (const char *fName_, double fMass_, int fCharge_, double fWidth_);
void SetFIndex (int a);
static void PrintArray ();
void PrintOther ();
double GetfPx () const;
double GetfPy () const;
double GetfPz () const;
double GetNewMass () const;
double InvMass (Particle &p);
void SetP(double px, double py, double pz);
int Decay2Body(Particle &dau1, Particle &dau2) const;

protected: 

private:
static const int fMaxNumParticleType = 10;
static ParticleType*fParticleType [fMaxNumParticleType];
static int fNParticleType;
int fIndex;
double fPx = 0;
double fPy = 0;
double fPz = 0;
static int FindParticle(const char *fName_);
void Boost (double bx, double by, double bz);
int a;





};


#endif

Particle.cxx

#include <iostream>
#include "ResonanceType.h"
#include "ParticleType.h"
#include "Particle.h"
#include <cmath>
#include <cstdlib> 

using namespace std;


int Particle::fNParticleType =0;
ParticleType*Particle:: fParticleType[] = {0x0};
const int Particle::fMaxNumParticleType;

Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
Particle::Particle():fPx(0), fPy(0), fIndex(-1){}

double Particle::GetfPx() const  
{
 return fPx;
}

double Particle::GetfPy() const
{
 return fPy;
}

double Particle::GetfPz() const
{
 return fPz;
}


void Particle::SetFIndex(int a) 
{
 fIndex=a;
}

double Particle::GetNewMass()  const
{
 return fParticleType[fIndex]->GetMass(); //MI DICE CHE E' UN PUNTATORE NULLO
}

int Particle::GetFIndex() 
{
 return fIndex;
}

//void Particle::SetFIndex (char*fName) {  //!

//;}

void Particle::SetP(double px, double py, double pz) 
{
  fPx = px; 
  fPy = py; 
  fPz = pz;
}



int Particle::FindParticle(const char *fName_)
{ cout<<"Sto eseguendo FindParticle";
 for (int i=0; i<fNParticleType; i++) 
 { 
   if (fName_ == fParticleType[i]->GetName())
   { cout<<"Ecco l'indice corrispondente";
     return i;
   }
   else 
   {
     cout<<"Non c'è corrispondenza"<<'\n';
   } 
 }
return -1;
}


void Particle::AddParticleType(const char *fName_, double fMass_, int fCharge_, double fWidth_)
{ cout<<"Sto eseguendo Addparticle";

   int b = FindParticle(fName_);
   if (b == -1) 
   {
     if (fWidth_ == 0) 
     {
       fParticleType[fNParticleType]=new ParticleType (fName_, fMass_, fCharge_); cout<<"Ecco una nuova ParticleType";
     }
     else 
     {
       fParticleType[fNParticleType]=new ResonanceType (fName_, fMass_, fCharge_, fWidth_); cout<<"Ecco una nuova ResonanceType";
     }
     fNParticleType+=1;
   }
 }




void Particle::PrintArray() 
{
 for(int i=0; i<fNParticleType; i++)
 {
   fParticleType[i]->Print();
 }
}  

void Particle::PrintOther() 
{
cout<<"Name"<<fParticleType[fIndex]->GetName()<<'\n'<<"fPx"<<fPx<<'\n'<<"fPy"<<fPy<<'\n'<<"fPz"<<fPz;
}


double Particle::TotalEnergy() const
{ 
 double Pulse=pow(GetfPx(), 2)+pow(GetfPy(), 2)+pow(GetfPz(), 2);
 double Masstwo= GetNewMass()*GetNewMass();
 double Energy=sqrt(Masstwo+Pulse);
  
  return Energy;
} 


double Particle::InvMass(Particle &p)
{ 
  double InvMass=pow(pow(TotalEnergy()+p.TotalEnergy(),2)-(pow(GetfPx()+p.GetfPx(),2)+pow(GetfPy()+p.GetfPy(),2)+pow(GetfPz()+p.GetfPz(),2)), 1/2);
  return InvMass;
} 


int Particle::Decay2Body(Particle &dau1, Particle &dau2) const 
 {
  if(GetNewMass() == 0.0){
    cout<<"Decayment cannot be preformed if mass is zero\n";
    return 1;
  }
  
  double massMot = GetNewMass();
  double massDau1 = dau1.GetNewMass();
  double massDau2 = dau2.GetNewMass();
  cout<<massMot<<'\n'<<massDau1<<'\n'<<massDau2;
  if(fIndex > -1){ // add width effect

    // gaussian random numbers

    float x1, x2, w, y1, y2;
    
    double invnum = 1./RAND_MAX;
    do {
      x1 = 2.0 * rand()*invnum - 1.0;
      x2 = 2.0 * rand()*invnum - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
    
    w = sqrt( (-2.0 * log( w ) ) / w );
    y1 = x1 * w;
    y2 = x2 * w;

    massMot += fParticleType[fIndex]->GetWidth() * y1;

  }

  if(massMot < massDau1 + massDau2){
    cout<<"Decayment cannot be preformed because mass is too low in this channel\n";
    return 2;
  }
  
  double pout = sqrt((massMot*massMot - (massDau1+massDau2)*(massDau1+massDau2))*(massMot*massMot - (massDau1-massDau2)*(massDau1-massDau2)))/massMot*0.5;

  double norm = 2*M_PI/RAND_MAX;

  double phi = rand()*norm;
  double theta = rand()*norm*0.5 - M_PI/2.;
  dau1.SetP(pout*sin(theta)*cos(phi),pout*sin(theta)*sin(phi),pout*cos(theta));
  dau2.SetP(-pout*sin(theta)*cos(phi),-pout*sin(theta)*sin(phi),-pout*cos(theta));

  double energy = sqrt(fPx*fPx + fPy*fPy + fPz*fPz + massMot*massMot);

  double bx = fPx/energy;
  double by = fPy/energy;
  double bz = fPz/energy;

  dau1.Boost(bx,by,bz);
  dau2.Boost(bx,by,bz);

  return 0;
}



void Particle::Boost(double bx, double by, double bz)
{
  double energy = TotalEnergy();

  //Boost this Lorentz vector
  double b2 = bx*bx + by*by + bz*bz;
  double gamma = 1.0 / sqrt(1.0 - b2);
  double bp = bx*fPx + by*fPy + bz*fPz;
  double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

  fPx += gamma2*bp*bx + gamma*bx*energy;
  fPy += gamma2*bp*by + gamma*by*energy;
  fPz += gamma2*bp*bz + gamma*bz*energy;
}

Mainmodule.cxx

#include <iostream>
#include "ParticleType.h"
#include "ResonanceType.h"
#include "Particle.h"
#include "cmath"
#include "TMath.h"
#include "TRandom.h"
#include "TH1D.h"
#include "TFile.h"
#include "TCanvas.h"

using namespace std; 

int Mainmodule() 

{ 

double Px, Py, Pz;
double Energia, Impulso_Trasverso, Impulso;
double MassaInv;

Particle particella [110]; //100 particelle per volta, sistema e riempie, poi lo fa per 10^5 volte
cout<<"L'errore non è prima degli add";


particella[0].AddParticleType("pioni+", 0.13957, 1, 0);
particella[1].AddParticleType("Kaoni+", 0.49367, 1, 0);
particella[2].AddParticleType("protoni+", 0.93827, 1, 0);
particella[3].AddParticleType("K*", 0.89166, 0, 0.050);
particella[4].AddParticleType("pioni-", 0.13957, -1, 0);
particella[5].AddParticleType("Kaoni-", 0.49367, -1, 0);
particella[6].AddParticleType("protoni-", 0.93827, -1, 0);



Double_t x = 0;
Double_t y = 0;
Double_t phi = 0;
Double_t theta = 0;
cout<<"L'errore non è negli add";

TFile *file=new TFile("Particelle.root", "RECREATE");

TH1D *HTypes=new TH1D ("HTypes", "Types", 7, 0, 110);
TH1D *HAzimuth=new TH1D ("HAzimuth", "Azimuth_Angle_Distribution", 100, 0, 2*TMath::Pi());
TH1D *HPolar=new TH1D ("HPolar", "Polar_Angle_Distribution", 100, 0, TMath::Pi());
TH1D *HPulse=new TH1D ("HPulse", "Pulse_Distribution", 100, 0, 1000);
TH1D *HTrasverso_Pulse=new TH1D ("HTrasverso_Pulse", "Trasverso_Pulse_Distribution", 100, 0, 1000);
TH1D *HEnergy=new TH1D ("HEnergy", "Energy_Distribution", 100, 0, 1000);
cout<<"L'errore non è nella prima parte degli histo";
TH1D *HInvEverybody=new TH1D ("HInvEverybody", "InvMass_Everybody", 500, 0, 1000);
TH1D *HInvOpposite_Sign=new TH1D ("HInvOpposite_Sign", "InvMass_Opposite_Sign", 500, 0, 1000);
TH1D *HInvSame_Sign=new TH1D ("HInvSame_Sign", "InvMass_Same_Sign", 500, 0, 5);
TH1D *HInvPione_concorde=new TH1D ("HInvPione_concorde", "InvMass_Pioneplus/Kaoneminus", 500, 0, 1000);
TH1D *HInvPione_discorde=new TH1D ("HInvPione_discorde", "InvMass_Pioneminus/Kaoneplus", 500, 0, 1000);
TH1D *HInvK=new TH1D ("HInvK", "InvMass_K*", 500, 0, 1000);
cout<<"L'errore non è negli istrogrammi";


 for (Int_t i=0; i<10e5; i++)
  {
   int n=0;
   for (Int_t j=0; j<100; j++)
   { 
	cout<<"asdasd"<<endl;
     x=gRandom->Rndm(); cout<<"Ho generato x";
     y=gRandom->Rndm(); cout<<"Ho generato y";
     phi=gRandom->Uniform(0, 2*TMath::Pi()); cout<<"Ho generato phi";
     theta=gRandom->Uniform(0, TMath::Pi()); cout<<"Ho genrato theta";
     Impulso=gRandom->Exp(1.0); cout<<"Ho generato Impulso";

     if (x<0.01) 
     {
          particella[j].SetFIndex(3);
          cout<<"Sono K*";
          
     }
     else 
     {
      
      if (x<0.10) 
      {
          if (y<0.50) 
         {
           particella[j].SetFIndex(2); cout<<"sono protone+";
           
         }
          else 
         {
           particella[j].SetFIndex(6); cout<<"sono protone -";
           
         }
      }
       else 
      {

        if (x<0.20) 
        {
          if (y<0.50) 
         {
           particella[j].SetFIndex(1); cout<<"sono Kaone +";
         }
          else 
         {
           particella[j].SetFIndex(5); cout<<"sono Kaone -";
         }
        }
        else 
        {

          if (x<=1.00)  
          {
            if (y<0.50) 
            {
             particella[j].SetFIndex(0); cout<<"sono pione +";
            }
            else 
            {
             particella[j].SetFIndex(4); cout<<"sono pione -";
             } 
          }
         }
        }
       }
      
   

HTypes->Fill(particella[j].GetFIndex()); //Rimepio Types

     Px=Impulso*sin(theta)*cos(phi);
     Py=Impulso*sin(theta)*sin(phi);
     Pz=Impulso*cos(theta);

HAzimuth->Fill(phi);
 //Rimepio Azimuth Angle

HPolar->Fill(theta);
 //Rimpio Polar Angle

HPulse->Fill(Impulso);
 // Rimepio Pulse
cout<<"Ho riempito qualche istogramma";

particella[j].SetP(Px,Py,Pz);
 //Impulso per ogni particella

cout<<"Roba impulso"<<'\n';

cout<<Px<<'\n'<<Py<<'\n'<<Pz<<'\n';
cout<<particella[j].GetfPx()<<'\n'<<particella[j].GetfPy()<<'\n'<<particella[j].GetfPz()<<'\n';

Impulso_Trasverso = pow((pow(Px, 2)+pow(Py, 2)),1/2);
HTrasverso_Pulse->Fill(Impulso_Trasverso);
cout<<"Rimepito impulso";

/*Riempio Trasverso Pulse


Energia = particella[j].Energy(); //SI BLOCCA QUI 
cout<<"Assegno l'energia";
HEnergy->Fill(Energia);
cout<<"rimepita l'energia";
 Rimepio Energy 
*/cout<<"ho finito di riempire istogrammi";



    if (x<0.01)
     {
       if (x<0.5) 
       { cout<<"sono k* e devo decadere";
         particella[100+n].SetFIndex(0); cout<<"ho settato il primo indice";
         particella[100+n+1].SetFIndex(5); cout<<"ho  settato il seocndo indice";
         particella[j].Decay2Body(particella[100+n], particella[100+n+1]); cout<<"sono decaduta";
      
       }
        else
       { cout<<"sono k* e devo decadere";
         particella[100+n].SetFIndex(4); cout<<"ho settato il primo indice";
         particella[100+n+1].SetFIndex(1); cout<<"ho  settato il seocndo indice";
        particella[j].Decay2Body(particella[100+n], particella[100+n+1]); cout<<"sono decaduta";
       }

        MassaInv=particella[100+n].InvMass(particella[100+n+1]);
        HInvK->Fill(MassaInv);
       n+=2;
       particella[100].PrintOther();
       
     } 
   
} //IL FOR DELLE 100 SI CHIUDE QUI 


   for(Int_t k=0; k<100+n; k++)
   { 
     cout<<"Prima delle masse invarianti ci sono";
     for (Int_t h=0; h<100+n; h++)
     {
      cout<<"ti sto dando una massa invariante";
     cout<<pow(particella[k].TotalEnergy()+particella[h].TotalEnergy(),2)-(pow(particella[k].GetfPx()+particella[h].GetfPx(),2)+pow(particella[k].GetfPy()+particella[h].GetfPy(),2)+pow(particella[k].GetfPz()+particella[h].GetfPz(),2));
      MassaInv=particella[k].InvMass(particella[h]); //SI BLOCCA QUI 

      HInvEverybody->Fill(MassaInv); cout<<"riempio 0";
//Rimepio Everybody


      if(( particella[k].GetFIndex()==0 || particella[k].GetFIndex()==1 || particella[k].GetFIndex()==2) && (particella[h].GetFIndex()==0 || particella[h].GetFIndex()==2 || particella[h].GetFIndex()==1)) 
      {
       HInvSame_Sign->Fill(MassaInv); cout<<"riempio 1";
      } 

      if((particella[k].GetFIndex()==4 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==6) &&
      (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==6))     
      {
       HInvSame_Sign->Fill(MassaInv); cout<<"riempio 2";

      }
     
       //concorde

       if((particella[k].GetFIndex()==4 || particella[k].GetFIndex()==5 || particella[k].GetFIndex()==6) &&
       (particella[h].GetFIndex()==0 || particella[h].GetFIndex()==1 || particella[h].GetFIndex()==2) )
       {
        HInvOpposite_Sign->Fill(MassaInv); cout<<"riempio 3";
       } 

       if((particella[k].GetFIndex()==0 || particella[k].GetFIndex()==1 || particella[k].GetFIndex()==2) && 
       (particella[h].GetFIndex()==4 || particella[h].GetFIndex()==5 || particella[h].GetFIndex()==6))       
        {
          HInvOpposite_Sign->Fill(MassaInv); cout<<"riempio 4";
        }
       
       //discorde

       if(((particella[k].GetFIndex()==0) && (particella[h].GetFIndex()==1)) || ((particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==5)) ||  ((particella[k].GetFIndex()==1) && (particella[h].GetFIndex()==0)) || ((particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==4)))
       {
        HInvPione_concorde->Fill(MassaInv); cout<<"riempio 5";

       }
       //pionekaoneconcorde

       if(((particella[k].GetFIndex()==0) && (particella[h].GetFIndex()==5)) || ((particella[k].GetFIndex()==4) && (particella[h].GetFIndex()==1)) ||  ((particella[k].GetFIndex()==5) && (particella[h].GetFIndex()==0)) || ((particella[k].GetFIndex()==1) && (particella[h].GetFIndex()==4)))
       {
        HInvPione_discorde->Fill(MassaInv); cout<<"riempio  6";
       }
       //pionekaonediscorde

}
cout<<"Sono alla "<< i <<"esima iterazione";
}
}
 file->Write(); cout<<"Scrivo il file ROOT";
 file->Close();


cout<<"Fine generazione";
return 0;
}

Thanks very much.

6 
root [2] .L Particle.cxx++g
Info in <TUnixSystem::ACLiC>: creating shared library /..././Particle_cxx.so
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:31: warning: field 'fPx' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                              ^
././Particle.cxx:15:11: note: during field initialization in this constructor
Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
          ^
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:42: warning: field 'fPy' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                                         ^
././Particle.cxx:15:11: note: during field initialization in this constructor
Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
          ^
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:53: warning: field 'fPz' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                                                    ^
././Particle.cxx:15:11: note: during field initialization in this constructor
Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
          ^
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:31: warning: field 'fPx' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                              ^
././Particle.cxx:16:11: note: during field initialization in this constructor
Particle::Particle():fPx(0), fPy(0), fIndex(-1){}
          ^
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:42: warning: field 'fPy' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                                         ^
././Particle.cxx:16:11: note: during field initialization in this constructor
Particle::Particle():fPx(0), fPy(0), fIndex(-1){}
          ^
In file included from input_line_12:9:
In file included from ././Particle.cxx:4:
/.../Particle.h:16:53: warning: field 'fPz' is uninitialized when used here [-Wuninitialized]
double Pulse = pow(1/2, pow(2,fPx)+pow(2,fPy)+pow(2,fPz));
                                                    ^
././Particle.cxx:16:11: note: during field initialization in this constructor
Particle::Particle():fPx(0), fPy(0), fIndex(-1){}
          ^

[...]$ `root-config --cxx --cflags` -O3 -Wall -Wextra -c Particle.cxx
Particle.cxx:15:1: warning: unused parameter ‘fName_’ [-Wunused-parameter]
 Particle::Particle(const char*fName_, double fPx_, double fPy_, double fPz_) : fPx(fPx_), fPy(fPy_), fPz(fPz_) {} 
 ^
In file included from Particle.cxx:4:0:
Particle.h: In constructor ‘Particle::Particle()’:
Particle.h:38:14: warning: ‘Particle::fPy’ will be initialized after [-Wreorder]
 double fPy = 0;
              ^
Particle.h:36:5: warning:   ‘int Particle::fIndex’ [-Wreorder]
 int fIndex;
     ^
Particle.cxx:16:1: warning:   when initialized here [-Wreorder]
 Particle::Particle():fPx(0), fPy(0), fIndex(-1){}
 ^
Particle.cxx: In member function ‘int Particle::Decay2Body(Particle&, Particle&) const’:
Particle.cxx:146:26: warning: variable ‘y2’ set but not used [-Wunused-but-set-variable]
     float x1, x2, w, y1, y2;
                          ^
7

Thanks very much.
But I have no idea how to solve the -Wuninitialized. What you advice, please?

8 
{
  // name this ROOT macro file "RunMe.cxx" and then debug everything using:
  // valgrind --tool=memcheck --leak-check=full --suppressions=`root-config --etcdir`/valgrind-root.supp `root-config --bindir`/root.exe -b -n -q -l RunMe.cxx
  // valgrind --tool=exp-sgcheck --suppressions=`root-config --etcdir`/valgrind-root.supp `root-config --bindir`/root.exe -b -n -q -l RunMe.cxx
  // note: study valgrind's output messages from running the "Mainmodule"
  gROOT->LoadMacro("ParticleType.cxx++g");
  gROOT->LoadMacro("ResonanceType.cxx++g");
  gROOT->LoadMacro("Particle.cxx++g");
  gROOT->LoadMacro("Mainmodule.cxx++g");
  std::cout << " ... starting Mainmodule() ..." << std::endl;
  gROOT->ProcessLine("Mainmodule();");
  std::cout << " ... finished Mainmodule() ..." << std::endl;
}

Particle.h (984 Bytes)
Particle.cxx (4.3 KB)
Particelle.root (8.6 KB)

I think, the main problem, as reported by the AddressSanitizer, is in the “Mainmodule.cxx”, where you define Particle particella [110]; but then you access particella[100+n+1] without any check that n<=8. The simplest “brutal fix” could be to define Particle particella [(100 + (100 * 2))];.

BTW. In the “Mainmodule.cxx”, better use delete file;, instead of file->Close();.

1 Like
9

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.