+ if (radius<2*kHeight*tan(theta+kMaxOmega)*3/4)
+ {
+
+ if(phi==0)
+ {
+ //printf("Radius: %f, Max Radius: %f\n",radius,2*kHeight*tan(theta+kMaxOmega)*3/4);
+ //printf("Loaded digit %d with coordinates x:%f, y%f\n",dig,x,y);
+ //printf("Using digit %d, for theta %f\n",dig,theta);
+ }
+
+ counter1++;
+
+ l=kHeight/cos(theta);
+
+ //x=x*kCorr;
+ //y=y*kCorr;
+ /*if(SnellAngle(theta+omega)<999)
+ {
+ //printf("(Angle)/(Snell angle):%f\n",(theta+omega)/SnellAngle(theta+omega));
+ x=x*(theta+omega)/SnellAngle(theta+omega);
+ y=y*(theta+omega)/SnellAngle(theta+omega);
+ }
+ else
+ {
+ x=0;
+ y=0;
+ }*/
+
+ //main calculation
+
+ DigitsXY->Fill(x,y,(float) 1);
+
+ theta1=SnellAngle(theta)*1.5;
+
+ aux1=-y*sin(phi)+x*cos(phi);
+ aux2=y*cos(phi)+x*sin(phi);
+ aux3=( TMath::Power(aux1,2)+TMath::Power(cos(theta1)*aux2 ,2))/TMath::Power(sin(theta1)*aux2+l,2);
+ omega=atan(sqrt(aux3));
+
+ //omega is distorted, theta1 is distorted
+
+ if(InvSnellAngle(theta+omega)<999)
+ {
+ omega1=InvSnellAngle(omega+theta1) - theta;
+ //theta1=InvSnellAngle(omega+theta) - omega1;
+ //omega1=kCorr*omega;
+
+ kCorr=InvSnellAngle(omega+theta)/(omega+theta);
+ theta1=kCorr*theta/1.4;
+ //if(phi==0)
+ //printf("Omega:%f Theta:%f Omega1:%f Theta1:%f ISA(o+t):%f ISA(t):%f\n",omega*180/kPi,theta*180/kPi,omega1*180/kPi,theta1*180/kPi,InvSnellAngle(omega+theta)*180/kPi,InvSnellAngle(theta)*180/kPi);
+ }
+ else
+ {
+ omega1=0;
+ theta1=0;
+ }
+
+ //printf("Omega:%f\n",omega);
+
+
+ //if(SnellAngle(theta+omega)<999)
+ //printf("(Angle)/(Snell angle):%f\n",(theta+omega)/SnellAngle(theta+omega));
+ if(theta==0 && phi==0)
+ {
+ //printf("Omega: %f Corrected Omega: %f\n",omega, omega/kCorr);
+ //omega=omega/kCorr;
+ }
+
+ //cout<<"\ni="<<i<<" theta="<<Int_t(2*theta*dimension/kPi)<<" phi="<<Int_t(2*phi*dimension/kPi)<<" omega="<<Int_t(2*omega*dimension/kPi)<<endl<<endl;
+ //{Int_t lixo;cin>>lixo;}
+ if(omega1<kMaxOmega && omega1>kMinOmega)
+ {
+ //printf("Omega found:%f\n",omega);
+ omega1=omega1-kMinOmega;
+
+ //printf("Omega: %f Theta: %3.1f Phi:%3.1f\n",omega, theta*180/kPi, phi*180/kPi);
+
+ bintheta=theta1*kDimensionTheta/kMaxTheta;
+ binphi=phi*kDimensionPhi/kMaxPhi;
+ binomega=omega1*kDimensionOmega/(kMaxOmega-kMinOmega);
+
+ if(Int_t(bintheta+0.5)==Int_t(bintheta))
+ inttheta=Int_t(bintheta);
+ else
+ inttheta=Int_t(bintheta+0.5);
+
+ if(Int_t(binomega+0.5)==Int_t(binomega))
+ intomega=Int_t(binomega);
+ else
+ intomega=Int_t(binomega+0.5);
+
+ if(Int_t(binphi+0.5)==Int_t(binphi))
+ intphi=Int_t(binphi);
+ else
+ intphi=Int_t(binphi+0.5);
+
+ //printf("Point added at %d %d %d\n",inttheta,intphi,intomega);
+ point[inttheta][intphi][intomega]+=1;
+ //printf("Omega stored:%d\n",intomega);
+ Points->Fill(inttheta,intphi,intomega,(float) 1);
+ ThetaPhi->Fill(inttheta,intphi,(float) 1);
+ OmegaTheta->Fill(inttheta,intomega,(float) 1);
+ OmegaPhi->Fill(intphi,intomega,(float) 1);
+ //printf("Filling at %d %d %d\n",Int_t(theta*kDimensionTheta/kMaxTheta),Int_t(phi*kDimensionPhi/kMaxPhi),Int_t(omega*kDimensionOmega/kMaxOmega));
+ }
+ //if(omega<kMaxOmega)point[Int_t(theta)][Int_t(phi)][Int_t(omega)]+=1;
+ }