Control on a sum

faca87 · March 9, 2021, 1:14pm

Hello, in this macro I calculate the emittance
emit.cpp (4.9 KB)

ts->SetBranchAddress("SecondaryParticleVert.x", &xVv[0]);
			ts->SetBranchAddress("PositionDirection.x", &xv[0]);
			ts->SetBranchAddress("SecondaryParticleAng.x", &xpv[0]);
 for(jentry=0; jentry<nentries;jentry++) {
				ts->GetEntry(jentry);
				x+=k*xv[0]-k*xVv[0];
				x2+=(k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
				xp+=xpv[0];
				xp2+=xpv[0]*xpv[0];
				xxp+=(k*xv[0]-k*xVv[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/=1e6;

I must you see, I get the values from the TTree, I save them in an array the I calculate the sums


   x+=k*xv[0]-k*xVv[0];
   x2+=(k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
   xp+=xpv[0];
   xp2+=xpv[0]*xpv[0];
   xxp+=(k*xv[0]-k*xVv[0])*xpv[0];

Currently I use all the values in the ROOT file, but now I need a check, i.e. I must sum the values only if the new value is lower than N times the previous one…
for example, lets set N=3
I start from x=0…then if using the first entry the new x value is 0.5, I use the entry, if the new value is 3.2 I don’t use that entry…
Moreover, I must count how many entries I used because currently I calulcate

   x2*=1./nentries;
   xp*=1./nentries;
   xp2*=1./nentries;
   xxp*=1./nentries;

but if I have 1000 entries and I use only 990, I must divide for 990 and not for 1000…
is there a way to add this control please?
Thank you

couet · March 9, 2021, 1:40pm

int n = 1;
double previous_value = value[0];
double sum = previous_value;
for (int i = 1; i <nentries; i++) {
   if (value[i] < previous_value/N) {
      previous_value = value[i] ;
      sum = sum + previous_value ;
      n++;
   } 
}

That’s n in the previous code.

faca87 · March 9, 2021, 2:39pm

Thank you @couet .
I maybe didn’t understand something, because I wrote this code but it doesn’t use any event

emit.cpp (5.5 KB)

couet · March 9, 2021, 2:52pm

Very likely be the test:

            if (( k*xv[i]-k*xVv[i])<previous_x/N  && ((k*xv[i]-k*xVv[i])*(k*xv[i]-k*xVv[i]))<previous_x2/N && xpv[i]<previous_xp/N && (xpv[i]*xpv[i])<previous_xp2/N && ((k*xv[i]-k*xVv[i])*xpv[i])<previous_xxp/N) {

filters all events.

faca87 · March 10, 2021, 7:57am

Thank you @couet I was thinking about this line

Maybe in my first message I didn’t explain what I need.

I start from the first value of the entries…and I assign the x, x2,xp, xp2 and xxp values.
Then I check the second entry if the new x,x2, xp,xp2 and xxp2 (ie. the sum of the first and secon entries) are lower than N times the old x,x2, xp, xp2 I use the second entry and I sum i.e.

x+=...
x2+=...
etc

otherwise I don’t use that entry. Same for the third entry… if the the sum of previous value and third entry is lower that N times the previous value I use the third entry

so maybe
if (value[i] < previous_value/N) {
isn’t the right check…or didn’t I understand the code?

I write an example

**First enrtry: **
x=0.5
xp=0.03

Second enrty
x=0.3
xp=0.04

I set the control N=3

the new x is 0.5+0.3=0.8
the new xp is 0.03+0.04=0.07
…then because 0.8<3*0.5 and 0.07< 3 * 0.03
I can use this entry

Second example
**First enrtry: **
x=0.5
xp=0.03

Second enrty
x=0.3
xp=0.4

the new x is 0.5+0.3=0.8 but the new xp is 0.03+0.4=0.43>3*0.03 then I don’t use the second entrry

couet · March 12, 2021, 8:43am

That’s what I understood from your first message but if your algorithm is a bit different simply program it . You have described it. You should be able to transcribe it into C++.

faca87 · March 14, 2021, 5:36pm

Thank you @couet I think that I did it!

faca87 · March 23, 2021, 6:54pm

Hello @couet sorry to write again regarding this problem.
I modified the macro to have this algorithm:

a. I’ve a vector xpv[0]
b. I get the first entry of xpv[0] and I set xp=xpv[0]
c. Regarding of netxt entries, If xp+xpv[0] is in the range [-cxp;cxp] I use the entry otherwise I dont’ use it. For example, I set c=1
First entry is 2 then I set xp=2
The second entry is -1, because 2-1=1 and 1 is in the range [-2;2] I use this value and my new xp is 1
The third entry is 1 because 1+1=2 and 2 isn’t in the range [-1;1], I don’t use this value and my xp is still 1.

I wrote this macro
emitcont.cpp (5.1 KB)

i.e. you see, I wrote the line
if(xp+xpv[0] <= abs(c*xp) ){

but I have a problem…you see here

for entry 3881 I’ve xp=-0.00411267 , the entry 3882 is -1.4384 and the macro get it, but it should reject this entry, because -0.00411267-1.4384=-1.44251267 and -1.44251267 is out of the range [-0.00411267;0.00411267]

Can you help to fix it please?

couet · March 24, 2021, 6:17am

Please check your macro before sending it:

root [0] 
Processing emiter.cpp...
In file included from input_line_9:1:
/Users/couet/Downloads/emiter.cpp:104:8: error: expected ';' after expression
                                n++
                                   ^
                                   ;

couet · March 24, 2021, 6:49am

Not really good either:

root emitcut.cpp
   ------------------------------------------------------------------
  | Welcome to ROOT 6.25/01                        https://root.cern |
  | (c) 1995-2021, The ROOT Team; conception: R. Brun, F. Rademakers |
  | Built for macosx64 on Mar 09 2021, 09:09:23                      |
  | From heads/master@v6-25-01-72-gb4566d45f8                        |
  | With Apple clang version 12.0.0 (clang-1200.0.32.29)             |
  | Try '.help', '.demo', '.license', '.credits', '.quit'/'.q'       |
   ------------------------------------------------------------------

root [0] 
Processing emitcut.cpp...
Error in <TUrl::TUrl>: /Users/couet/Downloads/C:/B1.root malformed, URL must contain "://"
Error in <TFile::TFile>: file /Users/couet/Downloads/C:/B1.root does not exist
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling_runtime_internal_throwIfInvalidPointer (no debug info)
[<unknown binary>] (no debug info)
[<unknown binary>] (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::IncrementalExecutor::executeWrapper(llvm::StringRef, cling::Value*) const (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::Interpreter::RunFunction(clang::FunctionDecl const*, cling::Value*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::Interpreter::EvaluateInternal(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, cling::CompilationOptions, cling::Value*, cling::Transaction**, unsigned long) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::MetaSema::actOnxCommand(llvm::StringRef, llvm::StringRef, cling::Value*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::MetaParser::isXCommand(cling::MetaSema::ActionResult&, cling::Value*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::MetaParser::isCommand(cling::MetaSema::ActionResult&, cling::Value*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] cling::MetaProcessor::process(llvm::StringRef, cling::Interpreter::CompilationResult&, cling::Value*, bool) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] HandleInterpreterException(cling::MetaProcessor*, char const*, cling::Interpreter::CompilationResult&, cling::Value*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] TCling::ProcessLine(char const*, TInterpreter::EErrorCode*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCling.so] TCling::ProcessLineSynch(char const*, TInterpreter::EErrorCode*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libCore.so] TApplication::ExecuteFile(char const*, int*, bool) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libRint.so] TRint::ProcessLineNr(char const*, char const*, int*) (no debug info)
[/Users/couet/git/roottrunk-bin/lib/libRint.so] TRint::Run(bool) (no debug info)
[/Users/couet/git/roottrunk-bin/bin/root.exe] main (no debug info)
[/usr/lib/system/libdyld.dylib] start (no debug info)
Error in <HandleInterpreterException>: Trying to dereference null pointer or trying to call routine taking non-null arguments.
Execution of your code was aborted.
In file included from input_line_9:1:
/Users/couet/Downloads/emitcut.cpp:23:4: warning: null passed to a callee that requires a non-null argument [-Wnonnull]
                fin->ls ();
                ^~~
root [1]

faca87 · March 24, 2021, 6:52am

@couet I see you miss the root file… I’m uploading it…but it’s 800MB…so it will need time! I will give you the link as soon as it will finish

Axel · March 24, 2021, 11:17am

emitcut.cpp does not have any cout << "jentry = " so I don’t think that’s the code that produces the output you are showing?

faca87 · March 24, 2021, 11:19am

Yes sorry… @Axel you are right

@couet… sorry I uploaded for the second time the bad macro!

this is the right one

emitcont.cpp (5.1 KB)

Axel · March 24, 2021, 11:54am

I cannot follow your reasoning. The code is

if (xp+xpv[0] <= abs(c*xp) ){

Here, c is 1, xp is -1.44251 and xpv[0] is -1.4384. So the if evaluates

if (-1.44251 + (-1.4384) <= abs(1.*(-1.44251)))

i.e.

if (-2.88091 <= 1.44251)

and the answer is: yes, totally.

Look - it’s better if you use printf to trace things you don’t understand in your code, or you use ACLiC (.L emitcont.cpp+) and use gdb to debug: it allows you to see all the values while things happen! It’s not really possible for us to help all our many many users get over their non-ROOT-related programming issues…

faca87 · March 26, 2021, 8:20am

Don’t worry @axel and @couet
I think I solved my self
writing

ts->GetEntry(0);
				x+=k*xv[0]-k*xVv[0];
				x2+=(k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
				xp+=xpv[0];
				xp2+=xpv[0]*xpv[0];
				xxp+=(k*xv[0]-k*xVv[0])*xpv[0];
		    	for(jentry=1; jentry<nentries;jentry++) {
				ts->GetEntry(jentry);
				if(abs(xp+xpv[0] )<= abs(c*xp)  && abs(x+(k*xv[0]-k*xVv[0])) <= abs(c*x) ){
				x+=k*xv[0]-k*xVv[0];
				x2+=(k*xv[0]-k*xVv[0])*(k*xv[0]-k*xVv[0]);
				xp+=xpv[0];
				xp2+=xpv[0]*xpv[0];
				xxp+=(k*xv[0]-k*xVv[0])*xpv[0];

system · April 9, 2021, 8:20am

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