Bose-Einstein correlations (BEC) are widely used to gain an insight into the spatiotemporal characteristics of boson emitters. It was used for the first time in the 1950s by R. Hanbury-Brown and R. Q. Twiss[Hanbury-Brown R, Twiss R Q 1954 Phil. Mag. 45 663] in astronomy to measure the dimension of distant astronomical objects emitting photons, and hence is also known as Hanbury-Brown-Twiss effect (HBT). In nuclear and particle physics field, BEC also has important applications in the investigation of the space-time properties of subatomic reaction region, especially in elementary-particle collisions and relativistic heavy-ion collisions with large multiplicity at high energies. Its potential application in exclusive reactions with low multiplicity in the non-perturbative QCD energy region may offer complementary information like duration and size of nucleon resonances, which are generally excited by hadronic or electromagnetic probes and usually decay into the ground states accompanied by emission of identical mesons. However, the event mixing technique, which is highly adopted for BEC observations in inclusive reactions at high energies with large multiplicity cannot be directly applied to the BEC measurement in exclusive reactions with very limited multiplicity at low energies. The event mixing method produces un-correlated samples from original sample through making mixed events by randomly selecting the momenta of two bosons from different original events. It works well for the high multiplicity case because the degree of freedom of final state particles is large compared with that of the low multiplicity case. In exclusive reactions with a very limited number of identical bosons in the final state, this method is however strongly interfered by non-BEC factors such as global conservation laws and decays of resonances. Appropriate constraints are required to control the event mixing process in order to eliminate the influence of those non-BEC factors. In this study, we are trying to develop an event mixing method for BEC measurement in reactions having only three final state particles and only two identical bosons among them. For this end, five constraint modes for the event mixing are proposed and investigated via Monte Carlo simulation. Each mode employs one or a combination of the following cut conditions:1) missing mass cut (MM) that requires the missing mass of the mixed event to be equal to that of the original event; 2) polar angle consistency cut (PAC) that requires that the swapping particles should come from the same polar angle bin; 3) azimuthal angle consistency cut (AAC); 4) momentum consistency cut (MC); 5) energy upper limit cut (EU) that requires that any boson energy should not exceed a given upper limit. The double neutral pion photoproduction on the proton around 1 GeV is taken for example to demonstrate the effects of these constraints on the event mixing. In the simulation, one event sample free of BEC effects and four samples in the presence of BEC effects are generated for testing the ability for these constraints to extract BEC parameters. It is found the constraint mode using the MM and PAC cuts, and the mode employing the MM and AAC cuts, and the mode adopting the MM and the EU cuts can be used to observe BEC effects and extract BEC parameters. Among them, optimum results can be achieved by the combination of the MM and EU cuts.