////
// Nuclear Instruments and Methods in Physics Research A 521 (2004) 261-298
/////
void fback(){

  gROOT->Reset();

#define MAX_nentries 100;

  TString pstit;

  // define style here 
  // general parameters
  gStyle->SetOptDate(0);     gStyle->SetOptTitle(0);
  gStyle->SetStatColor(10);  gStyle->SetStatH(0.3);
  gStyle->SetStatW(0.3);     
  gStyle->SetOptStat(1111); 
  
  // canvas parameters
  gStyle->SetFrameBorderMode(0);
  gStyle->SetFrameBorderSize(0);
  gStyle->SetFrameFillColor(0);
  gStyle->SetCanvasColor(0);

//
 
  // pads parameters
  //  gStyle->SetPadColor(39); 
  gStyle->SetPadColor(0); 
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadBorderSize(0);
  gStyle->SetPadBottomMargin(0.18);
  gStyle->SetPadRightMargin(0.10);
  gStyle->SetPadLeftMargin(0.14);
  gStyle->SetLabelSize(0.06,"x");
  gStyle->SetLabelSize(0.06,"y");
  gStyle->SetTitleXSize(0.08);
  gStyle->SetPaperSize(10,12);

  gStyle->SetOptLogx(0);
  gStyle->SetOptLogy(0);
  
  gStyle->SetTitleYOffset(0.8);
  gStyle->SetTitleYSize(0.10);
  gROOT->ForceStyle();  

  ////
  // Random number generators:

  TRandom *rsave = gRandom;
  TRandom *r0 = new TRandom();
  TRandom *r1 = new TRandom1();
  TRandom *r2 = new TRandom2();
  TRandom *r3 = new TRandom3();

  //

  TH1F* asymh1 = new TH1F("asymh1","asymh1", 200,-10.0,10.0);
  TH1F* asymh2 = new TH1F("asymh2","asymh2", 200,-10.0,10.0);

  //

  Double_t err1[MAX_nentries] = 0.0;

  Double_t asym1[MAX_nentries], asym2[MAX_nentries];
  Double_t nfb1[MAX_nentries], nfb2[MAX_nentries];

  Double_t Aopt = 0.10; // Constant Asymmetry, real optical asymmetry
  Double_t Actrl1 = 0.0; // Control Device Asymmetry
  Double_t Actrl2[MAX_nentries+1] = 0.0; // Control Device Asymmetry

  Double_t mean1 = 0.0;
  Double_t sigma1 = 1.0;
  Double_t mean2 = 0.0;
  Double_t sigma2 = 1.0;

  Double_t ave1[MAX_nentries], ave2[MAX_nentries];
  Double_t sum1 = 0.0;
  Double_t sum2 = 0.0;

  Int_t nentries = MAX_nentries;

  ///// No Feedback 

  r3->SetSeed(04242003); // Set the random generator sequence

  for (Int_t jentry=0; jentry<nentries;jentry++) {
    
    //    cout << " No. of Try = " << jentry << endl; 
    asym1[jentry] = r3->Gaus(mean1,sigma1) + Actrl1 + Aopt; // 
    //    cout << " asym (no FB) = " << asym1[jentry] << endl; 

    nfb1[jentry] = jentry+1.0;
    asymh1->Fill(asym1[jentry]);
    // 

    sum1 = sum1 + asym1[jentry];
    ave1[jentry] = sum1/(nfb1[jentry]);

    //    cout << " average asym (no FB) = " << ave2[jentry] << endl;

  }
 
 ///// With Feedback

  r3->SetSeed(04242003); // Set the random generator sequence

  for (Int_t jentry=0; jentry<nentries;jentry++) {
    
    //    cout << " No. of Try = " << jentry << endl; 
    asym2[jentry] = r3->Gaus(mean2,sigma2) + Actrl2[jentry] + Aopt; // 
    Actrl2[jentry+1] = Actrl2[jentry] - 1.0 * asym2[jentry]; // Unity gain
    //    cout << " asym (with FB) = " << asym2[jentry] << endl;

    nfb2[jentry] = jentry+1.0;
    asymh2->Fill(asym2[jentry]);
    // 

    sum2 = sum2 + asym2[jentry];
    ave2[jentry] = sum2/(nfb2[jentry]);

    //    cout << " average asym (with FB) = " << ave2[jentry] << endl;

  }
  
  TGraphErrors *grA1 = new TGraphErrors(nentries,nfb1,asym1,err1,err1);
  TGraphErrors *grA2 = new TGraphErrors(nentries,nfb2,asym2,err1,err1);

  TGraphErrors *grAve1 = new TGraphErrors(nentries,nfb1,ave1,err1,err1);
  TGraphErrors *grAve2 = new TGraphErrors(nentries,nfb2,ave2,err1,err1);

  /////
  
  b1 = new TCanvas("b1","Asymmtery Feedback",40,40,800,800);
  b1->Divide(1,3);

  b1->cd(1); 
  gPad->SetGrid();
  
  //
 grA1->GetXaxis()->SetLimits(0.0,nentries);
 grA1->SetMaximum(10.0);
 grA1->SetMinimum(-10.0);

 grA1->SetMarkerColor(kBlue);
 grA1->SetLineColor(kBlue);
 grA1->SetLineWidth(2);
 grA1->SetMarkerSize(1.0);
 grA1->SetMarkerStyle(24);
 
 grA1->GetXaxis()->CenterTitle();
 grA1->GetXaxis()->SetTitleOffset(1.2);
 grA1->GetXaxis()->SetTitleSize(0.06);
 TString tit = "Mini Run";
 grA1->GetXaxis()->SetTitle(tit.Data());
 
 grA1->GetYaxis()->CenterTitle();
 grA1->GetYaxis()->SetTitleOffset(1.0);
 grA1->GetYaxis()->SetTitleSize(0.06);
 TString tit = "Asymmetry";
 grA1->GetYaxis()->SetTitle(tit.Data());
 
 grA1->Draw("A p");

 //

 grA2->SetMarkerColor(kRed);
 grA2->SetLineColor(kRed);
 grA2->SetLineWidth(2);
 grA2->SetMarkerSize(1.0);
 grA2->SetMarkerStyle(24);
 
 grA2->Draw("same p");

 //

 leg = new TLegend(0.6,0.2,0.85,0.45);
 leg->SetFillColor(0);
 leg->SetShadowColor(0);
 leg->AddEntry(grA1,"No Feadback ","L");
 leg->AddEntry(grA2,"With Feadback ","L");
 leg->Draw();

 // //

 b1->cd(2); 
 gPad->SetGrid();

 asymh1->SetMarkerColor(kBlue);
 asymh1->SetLineColor(kBlue);
 asymh1->SetLineWidth(1);
 asymh1->SetMarkerSize(1.0);
 asymh1->SetMarkerStyle(24);
 
 asymh1->GetXaxis()->CenterTitle();
 asymh1->GetXaxis()->SetTitleOffset(1.2);
 asymh1->GetXaxis()->SetTitleSize(0.06);
 TString tit = "Asymmetry";
 asymh1->GetXaxis()->SetTitle(tit.Data());
 
 asymh1->GetYaxis()->CenterTitle();
 asymh1->GetYaxis()->SetTitleOffset(1.0);
 asymh1->GetYaxis()->SetTitleSize(0.06);
 TString tit = "";
 asymh1->GetYaxis()->SetTitle(tit.Data());
 
 asymh1->Draw("histo");

 // //

 b1->cd(3); 
 gPad->SetGrid();

 asymh2->SetMarkerColor(kRed);
 asymh2->SetLineColor(kRed);
 asymh2->SetLineWidth(1);
 asymh2->SetMarkerSize(1.0);
 asymh2->SetMarkerStyle(24);
 
 asymh2->GetXaxis()->CenterTitle();
 asymh2->GetXaxis()->SetTitleOffset(1.2);
 asymh2->GetXaxis()->SetTitleSize(0.06);
 TString tit = "Asymmetry";
 asymh2->GetXaxis()->SetTitle(tit.Data());
 
 asymh2->GetYaxis()->CenterTitle();
 asymh2->GetYaxis()->SetTitleOffset(1.0);
 asymh2->GetYaxis()->SetTitleSize(0.06);
 TString tit = "";
 asymh2->GetYaxis()->SetTitle(tit.Data());
 
 asymh2->Draw("histo");

 /////

 b1->cd();
 pstit = "fback.gif";
 b1->Print(pstit.Data());

 /////

  /////
  
  c1 = new TCanvas("c1","Average Feedback",40,40,800,500);
  c1->Divide(1,1);

  c1->cd(1); 
  gPad->SetGrid();
  
  //
 grAve1->GetXaxis()->SetLimits(0.0,nentries);
 grAve1->SetMaximum(1.0);
 grAve1->SetMinimum(-1.0);

 grAve1->SetMarkerColor(kBlue);
 grAve1->SetLineColor(kBlue);
 grAve1->SetLineWidth(2);
 grAve1->SetMarkerSize(1.0);
 grAve1->SetMarkerStyle(24);
 
 grAve1->GetXaxis()->CenterTitle();
 grAve1->GetXaxis()->SetTitleOffset(1.2);
 grAve1->GetXaxis()->SetTitleSize(0.06);
 TString tit = "Mini Run";
 grAve1->GetXaxis()->SetTitle(tit.Data());
 
 grAve1->GetYaxis()->CenterTitle();
 grAve1->GetYaxis()->SetTitleOffset(1.0);
 grAve1->GetYaxis()->SetTitleSize(0.06);
 TString tit = "Asymmetry";
 grAve1->GetYaxis()->SetTitle(tit.Data());
 
 grAve1->Draw("A p");

 //

 grAve2->SetMarkerColor(kRed);
 grAve2->SetLineColor(kRed);
 grAve2->SetLineWidth(2);
 grAve2->SetMarkerSize(1.0);
 grAve2->SetMarkerStyle(24);
 
 grAve2->Draw("same p");

 //

 leg = new TLegend(0.6,0.3,0.85,0.45);
 leg->SetFillColor(0);
 leg->SetShadowColor(0);
 leg->AddEntry(grAve1,"No Feadback ","L");
 leg->AddEntry(grAve2,"With Feadback ","L");
 leg->Draw();

 // //

 c1->cd();
 pstit = "fbackave.gif";
 c1->Print(pstit.Data());

}