Problem calculating the mean values

Hello, I wrote this code to calculate the emittance

Double_t emittance=0.;
			Double_t x=0;
			Double_t x2=0;	
			Double_t xp=0;
			Double_t xp2=0;
			Double_t xxp=0;
			Double_t sigmax2=0;
			Double_t sigmaxp2=0;
			Double_t sigmaxsigmaxp=0;
			int jentry;
			int nentries = ts->GetEntries();
			ts->SetBranchAddress("PositionDirection.x", &x);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xp);
		    for(jentry=0; jentry<nentries;jentry++) {
				ts->GetEntry(jentry);
				x+=x;
				x2+=x*x;
				xp+=xp;
				xp2+=xp*xp;
				xxp+=x*xp;
			}
			x*=1./nentries;
			x2*=1./nentries;
			xp*=1./nentries;
			xp2*=1./nentries;
			xxp*=1./nentries;
			sigmax2=x2-x*x;
			sigmaxp2=xp2-xp*xp;
			sigmaxsigmaxp=xxp-x*xp;
			std::cout << "N = " << nentries << "\t x = " << x << "\t x2 = " << 2 << " \t xp = " << xp <<  " \t xp2= " << xp2 << " \t xxp= " <<  xxp << std::endl;
			emittance=sigmax2*sigmaxp2-sigmaxsigmaxp;
			emittance=TMath::Sqrt(emittance);
			emittance*=1000000;

You see, that I’ve to calculate some mean values (as , etc) but these values are very different than the values calculated by Root. For example, you see here my =-3,765328*10^-7

Instead Root calculated = -0.000227

and surely the Root value is fine because using my value, I get a very big emittance (impossible)…but I can’t understand what is wrong in my code.

this is the macro emit.cpp (3.0 KB)
and this is the Root file https://we.tl/t-TMHdHaqL5D

Thank you


Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided


Hi,
The difference could be caused by the fact that in the ROOT mean value maybe the underflow and overflow;s you are having in y are not used. You should call before filling the histogram TH1::StatOverflows(true) and you should get the same value

Lorenzo

Hello @moneta thank you. Unfortunately I tried the code, and the problem was due to the code and not to the underflow and overflow!

Anyway I solvede myself
I had to use the array!

Double_t emittance=0.;
		Double_t x=0;
		Double_t x2=0;	
		Double_t xp=0;
		Double_t xp2=0;
		Double_t xxp=0;
		Double_t sigmax2=0;
		Double_t sigmaxp2=0;
		Double_t sigmaxsigmaxp=0;
		int jentry;
		int nentries = ts->GetEntries();
		double xv[10];
		double xpv[10];
		ts->SetBranchAddress("PositionDirection.x", &xv[0]);
		ts->SetBranchAddress("SecondaryParticleAng.x", &xpv[0]);
	    for(jentry=0; jentry<nentries;jentry++) {
			ts->GetEntry(jentry);
			x+=xv[0];
			x2+=xv[0]*xv[0];
			xp+=xpv[0];
			xp2+=xpv[0]*xpv[0];
			xxp+=xv[0]*xpv[0];
		}
		x*=1./nentries;
		x2*=1./nentries;
		xp*=1./nentries;
		xp2*=1./nentries;
		xxp*=1./nentries;
		sigmax2=x2-x*x;
		sigmaxp2=xp2-xp*xp;
		sigmaxsigmaxp=xxp-x*xp;
		emittance=sigmax2*sigmaxp2-sigmaxsigmaxp;
		emittance=TMath::Sqrt(emittance);
		emittance*=1000000;
		std::cout << "N = " << nentries << "\t x = " << x << "\t x2 = " << x2 << " \t xp = " << xp <<  " \t xp2= " << xp2 << " \t xxp= " <<  xxp << std::endl;

Now it works fine! Here the results and they are equal to the Root results!

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