GSI Helmholtzzentrum für Schwerionenforschung
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
Quote: Lia, could you please also communicate this message to Profs. Rotondi and Fontana in case they are not on the forum? Hi Christian, ok, I can forward your message and I will discuss this with Prof. Rotondi and will let you know Ciao, Lia.
Hi Christian, I talked to Prof. Rotondi about these problems. First of all, concerning the errors: 1) erpremc/trscsd.F line 148 and 3) erpremc/trsdsc.F line 138 he suggests to try changing the line: RD(J)=S(J) with: RD(J)=SNGL(S(J)) (and in the same way the RC(J)=S(J)) in order to convert the double precision variable to single precision one. Could you try this and tell us if the error is still there, please? In addition to this, can you please send me the macro/macros to be able to reproduce these errors? ...in this way we can make tests ourselves. Prof. Rotondi' s idea is also to modify GEANE in order to be able to handle the tracking error... do I remember correctly that you already modified a bit the FORTRAN code? If so, can you please summarize your changes? Thank you very much, Lia.
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
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
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
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.
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
Hi Christian, sorry for the late reply. I attach here the file patch.txt, I hope this is ok to make your tests... if something' s still missing please tell me! Ciao, Lia.
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
Thank you to you Christian! Nice to hear that the patch works: at this point, if there are no objections, we would send the changes to CERN in order to add them to the new geane version. Ciao, Lia.
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
Hi, It is planed for end of May. Mohammad
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.
If I remember well, there was even a sort of bug in VMC while propagating some track ids of secondaries, disappearing tracks, or something similar. Maybe Mohammad could comment on it, if I remember well the bug was quite "important" and I think it should be solved in the next ext packages release (or pre-release).
Hi, I am preparing the needed files to apply the patch to geant3: I just need a few time to test if the patching procedure works. Since, however, it seems that a new error handling might be needed (see GEANE on 64bit Ubuntu 9.04 [message #8403] in the forum) I would check if this requires some further modifications before providing the patch. Regards, Lia.
Hi, yes it sounds like a good idea to use the new error flag mechanism for these cases also. Christian
Hi, Ok, perfect. It might also be a good idea to make a QA macro which would do some rigourous tests for all the machines involved in the nightly builds. Best wishes, Johan.
Hi all, I sent the patch for GEANE to Rene Brun and he has just inserted it in geant3 SVN trunk Ciao, Lia.