void plot_double2(void)
{
// 
// plot the intensity distribution expected from the double slit experiment
//

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 100000;

   // define histograms

       Double_t ymin = -2;
       Double_t ymax = 2;

       //  dimensions are mm

       Double_t R0=1000.;             // reference distance
       Double_t lambda=0.0006;
       Double_t b=0.7;
       Double_t pupillary=1.2;
       Double_t nslits=2;
	
       TF1 *doubleslit= new TF1 ("doubleslit",func_doubleslit,ymin,ymax,4);
	// the parameters are 4. Unit are mm
	// R0 distance between slits and screen 
	// lambda - wavelength of light 
	// b - diameter of aperture
                    // pupillary - distance between apertures.
        doubleslit->SetParameters(R0,lambda,b,pupillary);
        doubleslit->SetParNames("R0","lambda","b","pupillary");	

       TF1 *singleslit= new TF1 ("singleslit",func_singleslit,ymin,ymax,3);
	// the parameters are 4. Unit are mm
	// R0 distance between slits and screen 
	// lambda - wavelength of light 
	// b - diameter of aperture
        singleslit->SetParameters(R0,lambda,b);
        singleslit->SetParNames("R0","lambda","b");

   //
   TCanvas *c1 = new TCanvas("c1","c1 Test random",200,10,700,600);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->SetGridx();
   // c1->SetGridy();
   // c1->SetLogy();

     /*
   c1->Divide(2,2);
   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);
   */
 
   doubleslit->SetTitle("");
   doubleslit->GetXaxis()->SetRangeUser(ymin,ymax);
   doubleslit->GetXaxis()->SetTitleSize(0.05);
   doubleslit->GetYaxis()->SetTitleSize(0.05);
   doubleslit->GetXaxis()->SetTitle("Position (mm)");
   doubleslit->GetYaxis()->SetTitle("Intensity");
   doubleslit->SetLineColor(2);
   doubleslit->Draw("c");   

    sprintf(string,"R_{0}= %.0f\n",R0);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
    sprintf(string,"#lambda= %.4f\n",lambda);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
    sprintf(string,"iris= %.2f\n",b);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
    sprintf(string,"pupillary= %.2f\n",pupillary);
   t1 = new TLatex(0.2,0.7,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();

   singleslit->SetTitle("");
   singleslit->GetXaxis()->SetRangeUser(ymin,ymax);
   singleslit->GetXaxis()->SetTitleSize(0.05);
   singleslit->GetYaxis()->SetTitleSize(0.05);
   singleslit->GetXaxis()->SetTitle("Position (mm)");
   singleslit->GetYaxis()->SetTitle("Intensity");
   singleslit->SetLineColor(1);
   singleslit->Draw("csame");  

   sprintf(filename,"plot_double2.pdf");
   c1->SaveAs(filename);
 

}
Double_t func_doubleslit (Double_t *x, Double_t *par){

	// double slit intensity profile from amplitudes
	
	Double_t R0 = par[0];
	Double_t  lambda = par[1];
	Double_t  b = par[2];
	Double_t  pupillary = par[3];
	Double_t y = x[0];

                    // Double_t pupillary =-1.2;                    
                    TComplex I(0,1);
                    TComplex One(1,0);
                    TComplex a1= 144.85*One + 0*I;
                    TComplex a2= 144.85*One -0.84*I;
                    a1 = a1 / a1.Rho();
                    a2 = a2 / a2.Rho();
                   
                    Double_t sinthe = y/sqrt(y*y + R0*R0);
                    Double_t arg2 = TMath::Pi()*pupillary*sinthe/lambda;
                    Double_t arg1 = TMath::Pi()*b*sinthe/lambda;

                    TComplex amp1=cos(arg2)*One + sin(arg2)*I;
                    TComplex amp2=cos(arg2)*One -sin(arg2)*I;
                    TComplex amp = a1*amp1 + a2*amp2;
                    Double_t interference = amp.Rho(); 

                    Double_t interference1 =1 ;
                    if (fabs(arg1) > 0) interference1 = sin(arg1)/arg1;
                    Double_t f = interference*interference*interference1*interference1;

                    // cout << " a1.rho=" << a1.Rho() << " a1.Theta()=" << a1.Theta() << " a2.rho=" << a2.Rho() << " a2.Theta()=" << a2.Theta() << endl;
                    // cout << "R0 =" << R0 << " lambda=" << lambda << " pupillary=" << pupillary << " b=" << b  <<  " interference1=" << interference1 << endl;
	
	return f;
}
Double_t func_singleslit (Double_t *x, Double_t *par){

	// single  slit intensity profile
	
	Double_t R0 = par[0];
	Double_t  lambda = par[1];
	Double_t  b = par[2];
	Double_t y = x[0];
                    
                    Double_t sinthe = y/sqrt(y*y + R0*R0);
                    Double_t arg1 = TMath::Pi()*b*sinthe/lambda;

                    Double_t interference1 =1 ;
                    if (fabs(arg1) > 0) interference1 = sin(arg1)/arg1;
                    Double_t f = 4*interference1*interference1;

                    // cout << "R0 =" << R0 << " lambda=" << lambda << " pupillary=" << pupillary << " b=" << b  << " option=" << option << endl;
	
	return f;
}