void test_complex(void)
{
// 
// script to check out complex numbers using TComplex, example of Breit-Wigner amplitude.
//

// #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 npts1 1001;
    Int_t npts = npts1-1;

    TComplex I(0,1);
    TComplex One(1,0);
     TComplex Amp[npts1];
     Double_t E[npts1];
     Double_t Mag[npts1];
     Double_t Phase[npts1];
     Double_t Real[npts1];
     Double_t Imaginary[npts1];

     Double_t E0 = 1.232;
     Double_t Gamma=1.18;
     Double_t Norm=1;
     Double_t Emin=-8;
     Double_t Emax=10;

    for (jj=0;jj<npts1;jj++) {
            E[jj] = Emin + (Emax-Emin)*jj/npts;
           Amp[jj] = (Gamma/2)/((E[jj]-E0) - I*Gamma/2);
           Mag[jj] = Amp[jj].Rho();
           Phase[jj] = Amp[jj].Theta();
           Real[jj] = Amp[jj].Re();
           Imaginary[jj] = Amp[jj].Im();
           // cout << " jj=" << jj << " E=" << E[jj] << " Amp[jj]=" << Mag[jj] << " Phase=" << Phase[jj] << endl;
           }

   //
   TCanvas *c1 = new TCanvas("c1","c1 test complex",200,10,700,700);
   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);  
   c1_1->SetGridx();
   c1_1->SetGridy();

   Double_t ymin=0;
   Double_t ymax=1;

   TGraph *res = new TGraph(npts,E,Mag);
   res->SetTitle("");
   res->GetXaxis()->SetRangeUser(Emin,Emax);
   res->GetYaxis()->SetRangeUser(ymin,ymax);
   res->GetXaxis()->SetTitleSize(0.05);
   res->GetYaxis()->SetTitleSize(0.05);
   res->GetYaxis()->SetTitleOffset(1.5);
   res->GetXaxis()->SetTitle("Energy (GeV)");
   res->GetYaxis()->SetTitle("Intensity");
   res->GetXaxis()->SetNdivisions(5);
   res->SetMarkerStyle(20);
   res->SetMarkerColor(2);
   res->SetMarkerSize(0.5);
   res->Draw("Ap");


   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);  
   c1_2->SetGridx();
   c1_2->SetGridy();

   Double_t ymin=0;
   Double_t ymax=4;

   TGraph *phaz = new TGraph(npts,E,Phase);
   phaz->SetTitle("");
   phaz->GetXaxis()->SetRangeUser(Emin,Emax);
   phaz->GetYaxis()->SetRangeUser(ymin,ymax);
   phaz->GetXaxis()->SetTitleSize(0.05);
   phaz->GetYaxis()->SetTitleSize(0.05);
   phaz->GetYaxis()->SetTitleOffset(1.5);
   phaz->GetXaxis()->SetTitle("Energy (GeV)");
   phaz->GetYaxis()->SetTitle("Phase (rad)");
   phaz->GetXaxis()->SetNdivisions(5);
   phaz->SetMarkerStyle(20);
   phaz->SetMarkerColor(2);
   phaz->SetMarkerSize(0.5);
   phaz->Draw("Ap");

   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);  
   c1_3->SetGridx();
   c1_3->SetGridy();

   Double_t ymin=0;
   Double_t ymax=4;

   /*TGraph *track= new TGraph(npts,Mag,Phase);
   track->SetTitle("");
   track->GetXaxis()->SetRangeUser(Emin,Emax);
   track->GetXaxis()->SetTitleSize(0.05);
   track->GetYaxis()->SetTitleSize(0.05);
   track->GetXaxis()->SetTitle("Intensity");
   track->GetYaxis()->SetTitle("Phase (rad)");
   track->GetXaxis()->SetNdivisions(5);
   track->SetMarkerStyle(20);
   track->SetMarkerColor(2);
   track->SetMarkerSize(0.5);
   track->Draw("Ap");*/

   Double_t ymin=-0.5;
   Double_t ymax=1;

   TGraph *ampR = new TGraph(npts,E,Real);
   ampR->SetTitle("");
   ampR->GetXaxis()->SetRangeUser(Emin,Emax);
   ampR->GetYaxis()->SetRangeUser(ymin,ymax);
   ampR->GetXaxis()->SetTitleSize(0.05);
   ampR->GetYaxis()->SetTitleSize(0.05);
   ampR->GetYaxis()->SetTitleOffset(1.5);
   ampR->GetXaxis()->SetTitle("Energy (GeV)");
   // ampR->GetYaxis()->SetTitle("Real Part");
   ampR->GetXaxis()->SetNdivisions(5);
   ampR->SetMarkerStyle(20);
   ampR->SetMarkerColor(2);
   ampR->SetMarkerSize(0.5);
   ampR->Draw("Ap");


   Double_t ymin=0;
   Double_t ymax=2;

   TGraph *ampIM = new TGraph(npts,E,Imaginary);
   /*ampIM->SetTitle("");
   ampIM->GetXaxis()->SetRangeUser(Emin,Emax);
   ampIM->GetYaxis()->SetRangeUser(ymin,ymax);
   ampIM->GetXaxis()->SetTitleSize(0.05);
   ampIM->GetYaxis()->SetTitleSize(0.05);
   ampIM->GetXaxis()->SetTitle("Energy (GeV)");
   ampIM->GetYaxis()->SetTitle("Imaginary Part");
   ampIM->GetXaxis()->SetNdivisions(5);*/
   ampIM->SetMarkerStyle(20);
   ampIM->SetMarkerColor(4);
   ampIM->SetMarkerSize(0.5);
   ampIM->Draw("Samep");
   // ampIM->Draw("Ap");

   TLegend *leg = new TLegend(0.15,0.8,0.45,0.9);
   leg->AddEntry(ampR,"Real Part","p");
   leg->AddEntry(ampIM,"Imaginary part","p");
   leg->Draw();

   c1->cd(4);
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);  
   c1_4->SetGridx();
   c1_4->SetGridy();

   Double_t xmin=-0.5;
   Double_t xmax=0.5;
   Double_t ymin=0;
   Double_t ymax=1;

   TGraph *RIM= new TGraph(npts,Real,Imaginary);
   RIM->SetTitle("");
   RIM->GetXaxis()->SetRangeUser(xmin,xmax);
   RIM->GetYaxis()->SetRangeUser(ymin,ymax);
   RIM->GetXaxis()->SetTitleSize(0.05);
   RIM->GetYaxis()->SetTitleSize(0.05);
   RIM->GetYaxis()->SetTitleOffset(1.5);
   RIM->GetXaxis()->SetTitle("Real Part");
   RIM->GetYaxis()->SetTitle("Imaginary Part");
   RIM->GetXaxis()->SetNdivisions(5);
   RIM->SetMarkerStyle(20);
   RIM->SetMarkerColor(4);
   RIM->SetMarkerSize(0.5);
   RIM->Draw("Ap");


}