Home » PANDA » PandaRoot » Tracking » still trouble with Geane
| still trouble with Geane [message #8202] |
Sat, 11 April 2009 16:49  |
|
Anonymous Poster
|
|
From: *pool.einsundeins.de
|
|
Hi Geane experts,
I still hit numerical errors in Geane. There could be more, but here is a list for the moment which trigger floating point exceptions:
1) erpremc/trscsd.F line 148
2) erpremc/trsdsc.F line 76
3) erpremc/trsdsc.F line 138
We need to figure out a way how to avoid these. I want to point out that this is not something that just happens in here at TUM. It will happen to everybody who will use Geane with tracks of the same geometry we tested in the last weeks.
I solved these cases locally at TUM in the following way:
For cases 1) and 3) I check whether ABS(S(J)).GT.1E30 and then I set IERR=1 and RETURN to ertrak. In ertrak I check for IERR.NE.0 and I set x2(1)=-1.E30 and I jump to the end with GOTO 99.
For case number 2) I check whether the argument for the ASIN call in line 76 is ABS(TN(3)).GT.1 and I abort in the same fashion.
In FairGeanePro I check for x2[0]<-1.E29 and I return kFALSE from all Propagates which call Ertrak.
This is not an elegant solution. However it is one I could do without changing the interface.
This has to be fixed in our global copy of geant3. Since the ertrak routine, does not have an error flag for returning, I do not know what the best solution is.
Please discuss this with me here on the forum and it would be great if you could send fixes to our group, because we have set up a lot of code for extensive testing.
This is one of the last steps we need to go for having a stable tracking.
Lia, could you please also communicate this message to Profs. Rotondi and Fontana in case they are not on the forum?
Cheers, Christian
|
|
|
|
|
|
|
|
| Re: still trouble with Geane [message #8263 is a reply to message #8257] |
Mon, 20 April 2009 16:09  |
|
Anonymous Poster
|
|
From: 141.39.167*
|
|
Hi Lia,
thanks for your effort. I dont think the SNGL solution will help. I tried it and then the SNGL causes the FPE, if the argument is out of range.
Another way of limiting the size was also tried out by me:
If I make a check before assigning S(J) and limit it in size, e.g.:
IF(S(J).GT.1.E30)S(J)=1.E30
IF(S(J).LT.-1.E30)S(J)=-1.E30
I dont get any errors anymore at this place. But what happens is that the error just moves to another place.
I want to say this: If the value becomes GT.1.E30 the track fit will be useless anyway. So we should abort the propagation in this case.
I want to suggest that you bring the an IERR flag in ertrack to the interface, so that we can abort the fit in the fortran code and than just detect the error in FairGeanePro. This is the only clean solution in my mind.
I produced an example for you to reproduce the error. Here is what needs to be done to run it:
Copy from genfit/benchmark the files benchmark.geo and media_pnd.geo to your geometry directory. The media file contains the vacuum2 material suggested by you some months ago. Then from your work directory, run
genfit/benchmark/runTpcMC.C
which will just do one event. This will only be used for material and B-field description in the actual test.
Then run:
genfit/benchmark/runGeaneTest.C
This should crash with an FPE at event#7. I made my own signal handler in GeaneTestTask.cxx to allow core dumping. So, if you do
ulimit -c unlimited
before, you will get a core file in your workdir. You can debug this with
gdb root.exe core
and then do 'where' to see the place where it happens.
In this example you will get the crash first in erpremc/trsdsc.F in
#6 0xb6e7f230 in asin () from /lib/tls/i686/cmov/libm.so.6
#7 0xb281ca8f in trsdsc_ (pd=0xb2bac2b0, rd=0xb2babbb0, pc=0xbf82e1c4,
rc=0xb2babbb0, h=0xb2bac770, ch=@0xb2a2beb8, ierr=@0xbf82e214,
spu=@0xbf82e200, dj=0xb2bac800, dk=0xb2bac80c) at erpremc/trsdsc.F:76
#8 0xb2812a59 in ertrak_ (x1=0xb7b3228, p1=0xb7b3234, x2=0xb7b319c,
p2=0xb7b31a8, ipa=@0xbf82e384, chopt=@0xb7a75b8, _chopt=2)
at erdecks/ertrak.F:167
#9 0xb28ff5f5 in TGeant3::Ertrak (this=0xb522058, x1=0xb7b3228, p1=0xb7b3234,
x2=0xb7b319c, p2=0xb7b31a8, ipa=5, chopt=0xb7a75b8 "PE")
at TGeant3/TGeant3.cxx:5392
Here the argument of the ASIN call in line 76 is greater than 1. If you limit it by
IF(TN(3).GT.1)TN(3)=1
then you will get to the error which I discussed many times now, around line 139 in the =S(J) business.
Please, if you have any trouble at all in getting the crash in the same event as me (#7), tell me what happens. I am not 100% sure, whether I forgot anything. These things in providing an exact error source to somebody are very difficult.
I will send you my diff of my erpremc/trsdsc.F and trscdc.F files and a comment on the ertrak.F file in aanother message in a minute, on my suggestions on fixing this problem.
Cheers, Christian
|
|
|
|
| Re: still trouble with Geane [message #8264 is a reply to message #8263] |
Mon, 20 April 2009 16:21 |
|
Anonymous Poster
|
|
From: 141.39.167*
|
|
Hi again,
now for my diffs of the fortran code:
--- ../../../fairsoft_orig/transport/geant3/erpremc/trscsd.F 2008-06-19 08:35:10.000000000 +0200
+++ erpremc/trscsd.F 2009-04-11 16:04:46.000000000 +0200
@@ -145,11 +145,13 @@ C
DO 25 I=1,5
DO 20 K=I,5
J=J+1
+ IF(ABS(S(J)).GT.1E30) GOTO 901
RD(J)=S(J)
20 CONTINUE
25 CONTINUE
C
RETURN
+
C
C *** ERROR EXITS
C
--- ../../../fairsoft_orig/transport/geant3/erpremc/trsdsc.F 2008-06-19 08:35:10.000000000 +0200
+++ erpremc/trsdsc.F 2009-04-11 16:05:01.000000000 +0200
@@ -73,6 +73,7 @@ C
5 CONTINUE
C
PC(1)=PD(1)
+ IF(ABS(TN(3)).GT.1) GOTO 99
PC(2)=ASIN(TN(3))
IF (ABS (TN(1)) .LT. 1.E-30) TN(1) = 1.E-30
PC(3) = ATAN2 (TN(2),TN(1))
@@ -135,11 +136,14 @@ C
DO 25 I=1,5
DO 20 K=I,5
J=J+1
+ IF(ABS(S(J)).GT.1E30) GOTO 99
RC(J)=S(J)
20 CONTINUE
25 CONTINUE
C
RETURN
+ 99 IERR=1
+ RETURN
C
C *** ERROR EXITS
C
And finally the ertrak.F file. But please remember: I dont have the possibility to bring the IERR out of the trsxsx.F routines to the outside world, this is why I encode the crash information in the x2 (final point). But I think you can solve this with the interface. Honestly, I just lack the skills....
--- ../../../fairsoft_orig/transport/geant3/erdecks/ertrak.F 2008-06-19 08:35:21.000000000 +0200
+++ erdecks/ertrak.F 2009-04-11 16:05:18.000000000 +0200
@@ -162,9 +162,17 @@
CALL VZERO (DUM,15)
CALL TRSCSD (ERPIN(1), DUM(1), ERPIN(1), DUM(1), HI(1),
+ CHARGE, IERR, SPU, ERPLI(1,1), ERPLI(1,2))
+ IF(IERR.NE.0) THEN
+ x2(1)=-1E30;
+ GOTO 99
+ ENDIF
IF (LEONLY) GOTO 35
CALL TRSDSC (ERPIN(1), ERRIN(1), DUM(1), ERRIN(1), HI(1),
+ CHARGE, IERR, SPU, ERPLI(1,1), ERPLI(1,2))
+ IF(IERR.NE.0) THEN
+ x2(1)=-1E30;
+ GOTO 99
+ ENDIF
DO 29 I = 1, 5
DO 28 J = 1, 5
ASDSC(I,J) = A(I,J)
These diffs are made with diff -rup, so you can copy pase them to a file and I apply them with the patch comman if you wish. But they are so little that it'll be better to implement them yourself.
Cheers, Christian
|
|
|
|
| Re: still trouble with Geane [message #8265 is a reply to message #8263] |
Mon, 20 April 2009 16:45 |
|
Anonymous Poster
|
|
From: 141.39.167*
|
|
Hi again,
in my instructions on how to reproduce the error I forgot something:
you need to include the directory genfit/benchmark in your global CMakeLists.txt since it isnt built by default, because nobody needs it usually.
Cheers, Christian
|
|
|
|
| Re: still trouble with Geane [message #8297 is a reply to message #8265] |
Fri, 24 April 2009 14:22 |
Lia Lavezzi
Messages: 291
Registered: May 2007
Location: Torino
|
first-grade participant |
From: *pv.infn.it
|
|
Hi Christian,
we made some changes in geane (listed below). Can you please run your tests with these changes and tell us if they are successful?
Here is the list of the changes you should do in geant3.
For erpremc/trscsd.F and trsdsc.F we just made your ones; in addition to these we also made modifications in: TGeant3/TGeant3.cxx and .h, erdecks/ertrak.F, comad/gcomad.F and geant321/ertrio.F.
The rationale is the following:
to assure a minimal impact and to not perturb the non PANDA current users of GEANE, we added to the GEANE common ERTRIO a labelled common ERTRIO1 containing three error flags. We use only IERTR for the moment.
Therefore, the unaware user will not detect the change, apart from a standard error message from ERTRAK when the FP exception happens.
The common modification required some interventions on the TGeant3 package.
In our interface we manage the ertrio1 common and skip the hit when IERTR=1, as you suggested.
For trscsd.F:
--- erpremc.orig/trscsd.F 2009-04-23 16:55:02.000000000 +0200
+++ erpremc/trscsd.F 2009-04-23 17:07:29.000000000 +0200
@@ -145,6 +145,7 @@ C
DO 25 I=1,5
DO 20 K=I,5
J=J+1
+ IF(ABS(S(J)).GT.1E30) GOTO 901
RD(J)=S(J)
20 CONTINUE
25 CONTINUE
and for trsdsc.F:
--- erpremc.orig/trsdsc.F 2009-04-23 16:55:02.000000000 +0200
+++ erpremc/trsdsc.F 2009-04-23 17:10:06.000000000 +0200
@@ -73,6 +73,7 @@ C
5 CONTINUE
C
PC(1)=PD(1)
+ IF(ABS(TN(3)).GT.1) GOTO 99
PC(2)=ASIN(TN(3))
IF (ABS (TN(1)) .LT. 1.E-30) TN(1) = 1.E-30
PC(3) = ATAN2 (TN(2),TN(1))
@@ -135,6 +136,7 @@ C
DO 25 I=1,5
DO 20 K=I,5
J=J+1
+ IF(ABS(S(J)).GT.1E30) GOTO 99
RC(J)=S(J)
20 CONTINUE
25 CONTINUE
@@ -142,6 +144,9 @@ C
RETURN
C
C *** ERROR EXITS
-C
+C
+ 99 IERR=1
+ RETURN
+*
END
*
For erdecks/ertrak.F:
--- erdecks.orig/ertrak.F 2009-04-23 16:55:31.000000000 +0200
+++ erdecks/ertrak.F 2009-04-23 17:05:04.000000000 +0200
@@ -73,6 +73,8 @@
*
* *** Decode character option
*
+ IERTR=0
+*
CHOPTI = CHOPT
CALL UOPTC (CHOPT, 'BELMOPVX', IOPT)
*
@@ -162,9 +164,25 @@
CALL VZERO (DUM,15)
CALL TRSCSD (ERPIN(1), DUM(1), ERPIN(1), DUM(1), HI(1),
+ CHARGE, IERR, SPU, ERPLI(1,1), ERPLI(1,2))
+ IF(IERR.NE.0) THEN
+*
+* *** Tracking error - floating point exception
+*
+ IERTR=1
+ WRITE (LOUT, 780)
+ GOTO 99
+ ENDIF
IF (LEONLY) GOTO 35
CALL TRSDSC (ERPIN(1), ERRIN(1), DUM(1), ERRIN(1), HI(1),
+ CHARGE, IERR, SPU, ERPLI(1,1), ERPLI(1,2))
+ IF(IERR.NE.0) THEN
+*
+* *** Tracking error - floating point exception
+*
+ IERTR=1
+ WRITE (LOUT, 780)
+ GOTO 99
+ ENDIF
DO 29 I = 1, 5
DO 28 J = 1, 5
ASDSC(I,J) = A(I,J)
For comad/gcomad.F:
--- comad.orig/gcomad.F 2009-04-23 17:10:55.000000000 +0200
+++ comad/gcomad.F 2009-04-23 17:12:16.000000000 +0200
@@ -253,6 +253,8 @@ C
IADD=GCADDI(IQUEST)
ELSE IF(CHCOMM.EQ.'ERTRIO') THEN
IADD=GCADDD(ERDTRP)
+ ELSE IF(CHCOMM.EQ.'ERTRIO1') THEN
+ IADD=GCADDD(IERTR)
ELSE IF(CHCOMM.EQ.'EROPTS') THEN
IADD=GCADDF(ERPLI)
ELSE IF(CHCOMM.EQ.'EROPTC') THEN
For TGeant3/TGeant3.cxx:
--- TGeant3.orig/TGeant3.cxx 2009-04-23 16:55:13.000000000 +0200
+++ TGeant3/TGeant3.cxx 2009-04-23 17:03:45.000000000 +0200
@@ -1262,6 +1262,7 @@ void TGeant3::LoadAddress()
// Commons for GEANE
gcomad(PASSCHARD("ERTRIO"),(int*&) fErtrio PASSCHARL("ERTRIO"));
+ gcomad(PASSCHARD("ERTRIO1"),(int*&) fErtrio1 PASSCHARL("ERTRIO1"));
gcomad(PASSCHARD("EROPTS"),(int*&) fEropts PASSCHARL("EROPTS"));
gcomad(PASSCHARD("EROPTC"),(int*&) fEroptc PASSCHARL("EROPTC"));
gcomad(PASSCHARD("ERWORK"),(int*&) fErwork PASSCHARL("ERWORK"));
For TGeant3/TGeant3.h:
--- TGeant3.orig/TGeant3.h 2009-04-23 16:55:13.000000000 +0200
+++ TGeant3/TGeant3.h 2009-04-23 17:24:25.000000000 +0200
@@ -494,7 +494,11 @@ typedef struct {
Int_t ilpred;
Int_t iepred;
} Ertrio_t;
-
+typedef struct {
+ Int_t iertr;
+ Int_t iertr1;
+ Int_t iertr2;
+} Ertrio1_t;
//-----------EROTPS
// CHARACTER*8 CHOPTI
// LOGICAL LEEXAC, LELENG, LEONLY, LEPLAN, LEPOIN, LEVOLU
@@ -774,6 +778,7 @@ public:
// Access to GEANE commons
virtual Ertrio_t* Ertrio() const {return fErtrio;}
+ virtual Ertrio1_t* Ertrio1() const {return fErtrio1;}
virtual Eropts_t* Eropts() const {return fEropts;}
virtual Eroptc_t* Eroptc() const {return fEroptc;}
virtual Erwork_t* Erwork() const {return fErwork;}
@@ -1054,6 +1059,7 @@ public:
Float_t weight, Int_t is);
Ertrio_t *fErtrio; //! ERTRIO common structure
+ Ertrio1_t *fErtrio1; //! ERTRIO1 common structure
Eropts_t *fEropts; //! EROPTS common structure
Eroptc_t *fEroptc; //! EROPTC common structure
Erwork_t *fErwork; //! ERWORK common structure
For geant321/ertrio.F:
--- geant321.orig/ertrio.inc 2009-04-23 17:15:52.000000000 +0200
+++ geant321/ertrio.inc 2009-04-23 17:22:41.000000000 +0200
@@ -32,6 +32,8 @@
+ ERTRSP(5,5,MXPRED), ERXIN( 3), ERXOUT( 3,MXPRED),
+ ERPIN(3), ERPOUT(3,MXPRED), NEPRED,INLIST,ILPRED,
+ IEPRED(MXPRED)
+ COMMON /ERTRIO1/ IERTR, IERTR1, IERTR2
+ INTEGER IERTR, IERTR1, IERTR2
#include "geant321/eropts.inc"
Then in the pandaroot/geane interface:
in FairGeanePro.cxx:
--- FairGeanePro.orig.cxx 2009-04-24 14:00:30.000000000 +0200
+++ FairGeanePro.cxx 2009-04-24 14:02:51.000000000 +0200
@@ -38,6 +38,7 @@ FairGeanePro::FairGeanePro() : TNamed("G
fdbPDG= TDatabasePDG::Instance();
fErrorMat= new TArrayD(15);
afErtrio=gMC3->fErtrio;
+ afErtrio1=gMC3->fErtrio1;
Pos=TVector3(0, 0 , 0);
PosErr = TVector3(0,0,0);
Mom=TVector3(0,0,0);
@@ -262,7 +263,9 @@ Bool_t FairGeanePro::Propagate(Float_t *
;
gMC3->Eufill(1, ein,xlf);
gMC3->Ertrak(x1,p1,x2,p2,GeantCode, "L");
- if(x2[0]<-1.E29) return kFALSE;
+
+ if(afErtrio1->iertr != 0) return kFALSE;
+ // if(x2[0]<-1.E29) return kFALSE;
}
Bool_t FairGeanePro::Propagate(Int_t PDG) {
@@ -272,7 +275,8 @@ Bool_t FairGeanePro::Propagate(Int_t PDG
cout << " FairGeanePro::Propagate ---------------------------"<< " " << x1[0]<< " "<<
x1[1]<< " "<< x1[2] << endl;
fApp->GeanePreTrack(x1, p1, PDG);
gMC3->Ertrak(x1,p1,x2,p2,GeantCode, fPropOption.Data());
- if(x2[0]<-1.E29) return kFALSE;
+ if(afErtrio1->iertr != 0) return kFALSE;
+ // if(x2[0]<-1.E29) return kFALSE;
trklength=gMC3->TrackLength();
Double_t trasp[25];
@@ -532,7 +536,8 @@ int FairGeanePro::FindPCA(Int_t pca, Int
//check needed for low momentum tracks
gMC3->Ertrak(x1,p1,x2,p2,GeantCode, fPropOption.Data());
- if(x2[0]<-1.E29) return 1;
+ if(afErtrio1->iertr != 0) return 1;
+ // if(x2[0]<-1.E29) return 1;
gMC3->GetClose(po1,po2,po3,clen);
// check on cases when only two steps are performed!
For FairGeanePro.h:
--- FairGeanePro.orig.h 2009-04-24 14:00:30.000000000 +0200
+++ FairGeanePro.h 2009-04-23 17:15:46.000000000 +0200
@@ -84,6 +84,7 @@ public:
TVector3 Mom;
TArrayD *fErrorMat;
Ertrio_t *afErtrio;
+ Ertrio1_t *afErtrio1;
Float_t x1[3];
Float_t p1[3];
Int_t GeantCode;
If also your tests will be positive, we would ask Mohammad to send these changes to CERN (or, if
Mohammad agrees we can do that). After the CERN release we will need the new external packages realease
and when it will be ready I will put my changes concerning the geane directory in the svn repository.
Best Regards,
Alberto and Lia.
|
|
|
|
| Re: still trouble with Geane [message #8309 is a reply to message #8297] |
Mon, 27 April 2009 19:16 |
|
Anonymous Poster
|
|
From: *natpool.mwn.de
|
|
Hi Lia,
thanks for your effort. Could you please send a patch file which I can apply to new copy of geant3. You can obtain it with:
diff -u --recursive --new-file pathToOriginalGeant3 pathToYourNewGeane > patch.txt
Then it will be easy to test. Like this it is a lot of work.
Thanks, Christian
|
|
|
|
|
|
| Re: still trouble with Geane [message #8353 is a reply to message #8320] |
Thu, 30 April 2009 12:19 |
|
Anonymous Poster
|
|
From: *e18.physik.tu-muenchen.de
|
|
Hi,
I tested your patch now. For me it works just fine. If it could find its way into the external packages that would be great.
- I will remove the theta<0.4deg cut from GeaneTrackRep.
- I will put Lia's changes to FairGeanePro, but I will protext them with an #ifndef #else thing, which will be default to still work with the 'old' version of geane. Then when one patches the geane, we can just change this flag and then remove the old code, once everybody has changed their Geane.
As a more general comment: I think we will have more cases where FPEs could occur in all the other trXXYY.F. So, if we could have the external packages source code in a repository that would be so useful, for future fixes as they become necessary. But now the interface is great to the C++ code and these fixes will be quick and painless!
Thanks a lot to the Pavia group. Please thank your Porfessors in case they dont read this.
Cheers, CHristian
|
|
|
|
|
|
| Re: still trouble with Geane [message #8395 is a reply to message #8393] |
Mon, 04 May 2009 13:37 |
|
Anonymous Poster
|
|
From: *pool.einsundeins.de
|
|
Sure thing Lia. You're most welcome.
One question to Mohammad:
How are your plans of handling external packages source code in the future? I am sure we will need to make some fixes here and there to that code.
Cheers, Christian
|
|
|
|
|
|
| Re: still trouble with Geane [message #8401 is a reply to message #8400] |
Mon, 04 May 2009 21:41 |
Johan Messchendorp
Messages: 693
Registered: April 2007
Location: University of Groningen
|
first-grade participant |
From: *xs4all.nl
|
|
Dear all,
Since it might take some time before the updates will be implemented and consequently in the next release of the external packages, I propose to prepare a new tarball of the external packages which only includes patches with the fixes related to Geane (no other updates). In preparation for the Torino meeting, we can place this tarball on the Wiki site, so that users and developers have enough time to install it. Meanwhile, we (.... eh, well lets say Mohammad) can work on the next official release of the external packages with -hopefully- an update from the vmc/root developers including Christian's and Lia's fixes and many more other new goodies. Lia would you be able to prepare a patched version of the external packages?
Note that the external packages will also be on subversion at the end of the month as Mohammad already announced, but with the main reason to make sure that the pandaroot releases will be coupled to the corresponding external packages. Concerning bug fixes in any code of the external packages, I would find it of utmost important that these fixes are (eventually) made in the original code in agreement with the responsible authors in order not to deviate from their developments. I realize that this does cause some extra overhead, but on the long run it will pay off.
Kind wishes and keep up the nice work!
Johan.
[Updated on: Mon, 04 May 2009 21:46] Report message to a moderator |
|
|
|
|
|
|
|
|
|
|
|
|
|
Goto Forum:
Current Time: Tue Dec 03 02:47:48 CET 2024
Total time taken to generate the page: 0.00912 seconds
|
GSI Forum
Search
Help
Members
Login
Home
