Hi Everyone,
I'm having a weird problem which probably has an obvious solution, but
which i can't figure out. My reconstruction code is giving me a 'divide
by zero' error in a place where it shouldn't. i've sent the code along
below if anyone is interested. If I run my reconstruction on my coda
file, (which has > 10000 events), this is what happens:
root j.C
Processing j.C...
Reading in File: bgrid.blast
form 1
xmax -10.000000 32 10.000000 = 300.000000
field check 32 x 17 x 52 x 3 = 84864, i= 84864
File will be used as the source of data. filename = jevent.coda
Opening file (jevent.coda) for data events.
event: 1000
event: 2000
event: 3000
event: 4000
event: 5000
Error: operator '/' divided by zero FILE:j.C LINE:50
*** Interpreter error recovered ***
this same thing happens, but after a different number of events, for
different coda files. line 50 is just:
while(readFlag == 0){
(where readFlag = event.ReadNext() )
i can't figure out why this would give me a 'divde be zero' error. if
i take out the lines where actually divide by things that may be very
small (where i assign values to ntpIn[8], etc), then things *seem* to run
smoothly, but how
would this affect the readFlag. any ideas?
thanks,
jason
here's the code:
{
gROOT->Reset();
gSystem.Load("libBlast.so");
const char CODAFILE[]="jevent.coda";
const char ELECFILE[]="electronics.map";
const char GEOMFILE[]="blast.geom";
const char BFLDFILE[]="bgrid.blast";
const float PI=4.0*atan(1.0);
const float R2D=180.0/PI;
const float D2R=PI/180.0;
const float PMIN=0.1; //GeV
const float PMAX=1.0; //GeV
const float THETAMIN=20.0; //deg
const float THETAMAX=110.0; //deg
const float PHIMIN=-17.5; //deg
const float PHIMAX=17.5; //deg
int readFlag;
long int goodEvents = 0;
long int totalEvents = 0;
long int eventMax = 1000000;
TBLEvent event;
TBLRecon recon;
// Declare the ntuple:
TNtuple *ntp = new
TNtuple("n","n","p:th:phi:z0:q2:pf:thf:phif:q2f:dpf:dthf:dphif:dq2f");
// 0 1 2 3 4 5 6 7 8 9
10 11 12
float ntpIn[13];
float ebeam = 0.880;
float eprime;
// now get the data file and read it:
event.Open(CODAFILE);
event.ReadElectronicsMap(ELECFILE);
recon.Initialize(BFLDFILE,GEOMFILE);
recon.setXp();
recon.setSectorZ();
event.ReadNext();
readFlag=0;
while(readFlag == 0){
// All of the next 5 lines must be present to eliminate "tails"...
readFlag=event.ReadNext();
recon.WChits2clus(event.raw);
recon.WCclus2stubs(event.raw);
recon.WCstubs2seg2();
recon.WClink();
// since only electrons are in the event, just take events that
// have only one possible link for simplicity.
if(recon.first.nLink==1){
totalEvents++;
if(totalEvents%1000==0){
cout << "event: " << totalEvents << endl;
}
recon.first.nFast=0;
recon.WCfastFit(0);
recon.WCfastFitCorrect();
recon.first.nBetter=0;
recon.WCfit1(0,recon.first.fast[0][0],recon.first.fast[0][1],
recon.first.fast[0][2],recon.first.fast[0][3],-1);
ntpIn[0] = event.mc[0].MCp;
ntpIn[1] = event.mc[0].MCtheta;
if(event.mc[0].MCphi < -90.0){
ntpIn[2] = -1*event.mc[0].MCphi;
}
else{
ntpIn[2] = event.mc[0].MCphi;
}
ntpIn[3] = event.mc[0].MCz0;
// next calculate the reconstructed for mc, fast fit (f), and best
fit (b)
eprime = ebeam/(1 +
2*ebeam*sin(ntpIn[1]*D2R/2.0)*sin(ntpIn[1]*D2R/2.0)/0.9383);
ntpIn[4] =
4.0*ebeam*eprime*sin(ntpIn[1]*D2R/2.0)*sin(ntpIn[1]*D2R/2.0);
ntpIn[5] = recon.first.better[0][0];
ntpIn[6] = recon.first.better[0][1]*R2D;
if(recon.first.better[0][2]*R2D < -90.0){
ntpIn[7] = -1*recon.first.better[0][2]*R2D;
}
else{
ntpIn[7] = recon.first.better[0][2]*R2D;
}
eprime = ebeam/(1 +
2*ebeam*sin(ntpIn[6]*D2R/2.0)*sin(ntpIn[6]*D2R/2.0)/0.9383);
ntpIn[8] =
4.0*ebeam*eprime*sin(ntpIn[6]*D2R/2.0)*sin(ntpIn[6]*D2R/2.0);
ntpIn[9] = 100*(ntpIn[5]-ntpIn[0])/ntpIn[0];
ntpIn[10] = 100*(ntpIn[6]-ntpIn[1])/ntpIn[1];
ntpIn[11] = 100*(ntpIn[7]-ntpIn[2])/ntpIn[2];
ntpIn[12] = 100*(ntpIn[8]-ntpIn[4])/ntpIn[4];
ntp->Fill(ntpIn);
}
// cout << readFlag << " " << totalEvents << " " << ntpIn[0]
// << " " << ntpIn[1] << " " << ntpIn[2] << " " << ntpIn[3]
// << " " << ntpIn[4] << " " << ntpIn[5] << " " << ntpIn[6]
// << " " << ntpIn[7] << endl;
}
event.CloseCoda();
TFile *f = new TFile("jevent.root", "recreate");
ntp->Write();
f->Close();
}
This archive was generated by hypermail 2.1.2 : Mon Feb 24 2014 - 14:07:28 EST