gStyle->SetPalette(1);
gStyle->SetErrorX(0);
gStyle->SetOptStat(0);

TFile *trans = TFile::Open("SDInnerScint.root");
TTree *tree= (TTree*)trans->Get("tree");

TFile *foil = TFile::Open("SDCryoE0.root");
TTree *plane= (TTree*)foil->Get("tree");



TCanvas *ec1;
TCanvas *ec2;


double meanMU, meanmuplus, ie ;

  double  implE;

double Bfield=0;

TGraph *g1, *g2;


void esquisse(){
gStyle->SetPalette(1);
  ec1= new TCanvas();
  ec2= new TCanvas();

  ec1->SetGrid();
  ec1->Divide(2,2);

  TString poscut="positionz==14";

  ec1->cd(1);
  plane->Draw("positionz>>muPOS(1500,14,14.00003)");

  ec1->cd(2);
  plane->Draw("muon",poscut);

  ec1->cd(3);
  plane->Draw("globaltime>>hTOF(150)",poscut,"E1");

  ec1->cd(4);
  plane->Draw("kenergy+edeposit>>hKE(150)",poscut,"E1");



  ec2->SetGrid();
  ec2->Divide(2,2);
  
  TFile *foil = TFile::Open("SDCryoE0.root");
  TTree *plane= (TTree*)foil->Get("tree");
  
  ec2->cd(1);
  plane->Draw("positiony:positionx>>mxy(100,100)",poscut+"&&muon==1","zcol");
  
  ec2->cd(2);
  plane->Draw("positiony:positionx>>Mxy(100,100)",poscut+"&&muonium==1","zcol");
  
  ec2->cd(3);
  plane->Draw("globaltime>>mpg(100)",poscut+"&&muon==1","E1");
  
  ec2->cd(4);
  plane->Draw("globaltime>>mg(100)",poscut+"&&muonium==1","E1");
 
  
    
}


void moyennes(int ID=-1)
{

gStyle->SetPalette(1);
gStyle->SetErrorX(0);
gStyle->SetOptStat(11);

TCanvas *ec1;
TCanvas *ec2;


double meanMU, meanmuplus, ie ;

  double  implE;

  TString cut = "positionz==14";

  if(ID!=-1)
    {
      cut=cut+"&&runID=="+FS(ID);
    }

  plane->Draw("globaltime>>hmuplus(1000,0,500)",cut+"&&muon==1");
  meanmuplus = hmuplus->GetMean();

  plane->Draw("globaltime>>hMU(1000,0,500)",cut+"&&muonium==1");
  meanMU = hMU->GetMean();
   
  //  plane->Draw("Init_Energy>>hIE");
  plane->Draw("kenergy+edeposit>>hIE(500,0,5)",cut+"&&muon==1");

   ie=hIE->GetMean();


  ofstream outf;

  outf.open("means.dat");
  
  if(outf.is_open())
    {
      std::cout<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
       outf<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
      outf.close();
      std::cout<<"Output file means.dat created."<<std::endl;
    }
  else
    {
      cout << "Error opening file";
      return 0;
    }
  
}


void mtime(TString particle="mu+",TString grids="ON", TString guards="ON"){
  
 gDirectory->DeleteAll();

 int ID;

  std::cout<<"Run? (-1 for all): " <<std::endl;
  std::cin>>ID;
  std::cout<<"Magnetic Field [Gauss]?" <<std::endl;
  std::cin>>Bfield;
  //  std::cout<<"Implantation Energy?" <<std::endl;
  //  std::cin>>implE;

  int b_f=1;
  if (Bfield==0)b_f=0;

#define GMUON 0.01355
#define TAU_MUON 2.197 
#define ASYM 0.28
  Int_t i,j;
  double r;
  char  func_str[256];
  char  text_str[256];
 
 GetValues();

  
  sprintf(func_str, "exp(-x/%f)", TAU_MUON);
  TF1 *f1 = new TF1("f1", func_str, 0., 20.);
  //sprintf(func_str, "exp(-x/%f)", TAU_MUON);
  TF1 *f2 = new TF1("f2", func_str, 0., 20.);
  TF1 *fitdecay = new TF1("fitdecay", TFdecayfit, meanmuplus, 10000., 2);
  TF1 *fitasym1  = new TF1("fitasym1", TFasymuonfit, meanmuplus, 10000., 3);
  TF1 *fitasym2  = new TF1("fitasym2", TFasymuonfit, meanmuplus, 10000., 3);

  /*----------- initialise fit hforw and hback ----------------*/
  fitdecay->SetParNames("N(0)","MuonLifeTime");
  fitdecay->SetParameters(0, 450.);
  fitdecay->SetParameters(1, 2.2);

  // convert frequence from Hz 1/s to 1/ns (10^-9 Hz)
  double fz = 13.*b_f; 
  double ampl = 0.2; 
  
  if(particle=="Mu")
    {
      fz = -1.39e3*b_f;
      ampl = 0.16;
    }
  fitasym1->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
  fitasym1->SetParameters(ampl,fz,0*b_f);
  fitasym1->SetLineColor(46);


  fitasym2->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
  fitasym2->SetParameters(ampl,fz,0.5*b_f);
  fitasym2->SetLineColor(50);


  //===================================================================
  //===================================================================

 /*---------------------------------------------------------*/
 TH1F *hRight = new TH1F("hRight", "hRight",100, 0., 2400.);
 TH1F *hLeft = new TH1F("hLeft", "hLeft",100, 0., 2400.);
 TH1F *LRdiff = new TH1F("LRdiff", "LRdiff",100, 0., 2400.);
 TH1F *LRsum  = new TH1F("LRsum",  "LRsum",100, 0., 2400.);
 TH1F *LRasym = new TH1F("LRasym", "LRasym",100, 0., 2400.);
 /*---------------------------------------------------------*/
 TH1F *hTop = new TH1F("hTop", "hTop",100, 0., 2400.);
 TH1F *hBottom = new TH1F("hBottom", "hBottom",100, 0., 2400.);
 TH1F *BTdiff = new TH1F("BTdiff", "BTdiff",100, 0., 2400.);
 TH1F *BTsum  = new TH1F("BTsum",  "BTsum",100, 0., 2400.);
 TH1F *BTasym = new TH1F("BTasym", "BTasym",100, 0., 2400.);
 /*---------------------------------------------------------*/
 if(particle=="Mu")
   { 
 hRight->SetBins( 100, 0.,  meanMU+200.);
 hLeft->SetBins(  100, 0.,  meanMU+200.);
 LRdiff->SetBins( 100, 0.,  meanMU+200.);
 LRsum->SetBins(  100, 0.,  meanMU+200.);
 LRasym->SetBins( 100, 0.,  meanMU+200.);
 /*---------------------------------------------------------*/
 hTop->SetBins(100, 0., 500.);
 hBottom->SetBins(100, 0., 500.);
 BTdiff->SetBins(100, 0., 500.);
 BTsum->SetBins(100, 0., 500.);
 BTasym->SetBins(100, 0., 500.);
   }
 /*---------------------------------------------------------*/
   


 hLeft->Sumw2();
 hRight->Sumw2();
 hTop->Sumw2();
 hBottom->Sumw2();
 LRdiff->Sumw2();
 LRsum->Sumw2();
 LRasym->Sumw2();
 BTdiff->Sumw2();
 BTsum->Sumw2();
 BTasym->Sumw2();

  //===================================================================
  //===================================================================
  
  //  moyennes();
  TString cut2;
  cut2 = "";
  if(ID==-1)
    {  
      cut2 = "";
    }
  else
    {
      cut2 = "positron.runID=="+FS(ID)+"&&";
    }

  tree->Draw("positron.globaltime>>hLeft",cut2+"positron.Left==1","");
  tree->Draw("positron.globaltime>>hRight",cut2+"positron.Right==1","");
  tree->Draw("positron.globaltime>>hBottom",cut2+"positron.Bottom==1","");
  tree->Draw("positron.globaltime>>hTop",cut2+"positron.Top==1","");
   
 
  LRsum->Add(hLeft,hRight,1,1);
  LRdiff->Add(hLeft,hRight,1,-1);
  LRasym->Divide(LRdiff,LRsum);

  hLeft->SetLineColor(1);
  hTop->SetLineColor(2);
  hRight->SetLineColor(3);
  hBottom->SetLineColor(4);


 

  //=============================================================================
  //                              LEFT-RIGHT
  //=============================================================================




  c1=new TCanvas("c1","LR Counts");
  c1->SetGrid();
  c1->SetFillColor(0);
  c1->SetFrameFillColor(0);
  c1->SetFrameBorderMode(0);
  
  hLeft->SetTitle("Counts in Left and Right Scintillators");
  hLeft->GetYaxis()->SetTitle("Counts (arbitrary scale)");
  hLeft->GetXaxis()->SetTitle("Time [ns]");
  hLeft->SetMarkerStyle(3);
  hLeft->SetMarkerColor(4);
  hLeft->SetLineColor(4);

  hLeft->Draw("E1lp");



  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
  pt->SetName("title");
  pt->SetTextSize(0.032);
   pt->SetBorderSize(0);
  pt->SetFillColor(0);
  text = pt->AddText("");
  pt->Draw();

  
  hRight->SetMarkerStyle(26);
  hRight->SetMarkerColor(2);
  hRight->SetLineColor(2);
  hRight->Draw("E1lpsame");

  TOF_plot();
  
  leg = new TLegend(0.2,0.75,0.43,0.88);
  leg->SetFillColor(0);
  leg->SetBorderSize(1);
  
  leg->AddEntry(hLeft,"Left ","p");
  leg->AddEntry(hRight,"Right ","p");
    leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
  leg->Draw();
  



  //++++++++++++++++++++++++++++++++++
  //++++++++++++++++++++++++++++++++++

  c3=new TCanvas("c3","LR Asymmetry");
  c3->SetGrid();
  c3->SetFillColor(0);
  c3->SetFrameFillColor(0);
  c3->SetFrameBorderMode(0);
  TString title = "Left/Right Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";


  LRasym->SetTitle(title);
  LRasym->GetXaxis()->SetTitle("Time [ns]");
  LRasym->GetYaxis()->SetTitle("Asymmetry");
  LRasym->GetYaxis()->SetRangeUser(-1,1);
  LRasym->SetMarkerStyle(2);
  LRasym->SetMarkerColor(38);
  LRasym->SetLineColor(38);

  LRasym->Fit("fitasym1");
  //  gStyle->SetOptStat(0);
  gStyle->SetOptFit(111);  
  LRasym->Draw("e1");

 
  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
  pt->SetName("title");
  pt->SetBorderSize(0);
  pt->SetTextSize(0.032);
   pt->SetFillColor(0);
  text = pt->AddText("");
  pt->Draw();



  
  TOF_plot();

  

  leg = new TLegend(0.2,0.75,0.43,0.88);
  leg->SetFillColor(0);
  leg->SetBorderSize(1);

  leg->AddEntry(LRasym,"L-R Asym","p");
  leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
  leg->Draw();


  //=============================================================================
  //                             Bottom-Top
  //=============================================================================
gStyle->SetPalette(1);
gStyle->SetOptStat(11);
gStyle->SetErrorX(0);



  BTsum->Add(hBottom,hTop,1,1);
  BTdiff->Add(hBottom,hTop,1,-1);
  BTasym->Divide(BTdiff,BTsum);
  BTasym->SetMarkerStyle(3);


  c4=new TCanvas("c4","BT Counts");
  c4->SetGrid();
  c4->SetFillColor(0);
   c4->SetFillColor(0);
  c4->SetFrameFillColor(0);
  c4->SetFrameBorderMode(0);
  hTop->SetTitle("Counts in Top and Bottom Scintillators");

  hTop->SetMarkerStyle(26);
  hTop->SetMarkerColor(2);
  hTop->SetLineColor(2);
 hTop->GetYaxis()->SetTitle("Counts (arbitrary scale)");
  hTop->GetXaxis()->SetTitle("Time [ns]");
  hTop->Draw("lp");


  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
  pt->SetName("title");
  pt->SetBorderSize(0);
  pt->SetTextSize(0.032);
   pt->SetFillColor(0);
  text = pt->AddText("");
  pt->Draw();



  hBottom->SetMarkerStyle(3);
  hBottom->SetMarkerColor(4);
  hBottom->SetLineColor(4);
  hBottom->Draw("lpsame");

  TOF_plot();

  leg = new TLegend(0.2,0.75,0.43,0.88);
  leg->SetFillColor(0);
  leg->SetBorderSize(1);
  leg->AddEntry(hBottom,"Bottom ","p");
  leg->AddEntry(hTop,"Top ","p");
  leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
  leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
  leg->Draw();
  


  //++++++++++++++++++++++++++++++++++
  //++++++++++++++++++++++++++++++++++


  c6=new TCanvas("c6","BT Asymmetry");
  c6->SetGrid();
  c6->SetFillColor(0);
  c6->SetFrameFillColor(0);
  c6->SetFrameBorderMode(0);



  title = "Bottom/Top Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";
  BTasym->SetTitle(title);
  BTasym->GetXaxis()->SetTitle("Time [ns]");
  BTasym->GetYaxis()->SetTitle("Asymmetry");
  BTasym->GetYaxis()->SetRangeUser(-1,1);
  BTasym->SetMarkerColor(31);
  BTasym->SetLineColor(31);


  //  BTasym->GetXaxis()->SetRangeUser(0,5000);
  BTasym->Fit("fitasym2");
  gStyle->SetOptStat(11);
  gStyle->SetOptFit(111);  

  BTasym->Draw("p");
  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
  pt->SetName("title");
  pt->SetTextSize(0.032);
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  text = pt->AddText("");
  pt->Draw();

  TOF_plot();

  leg = new TLegend(0.2,0.75,0.43,0.88);
  leg->SetFillColor(0);
  leg->SetBorderSize(1);

  leg->AddEntry(BTasym,"B-T Asym","p");
    leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
 leg->Draw();
  

  //  c6->SetEditable(0);



}




TString flt2string(const float& number)
{
  ostringstream oss;
  oss << number;
  return oss.str();
}

TString FS(const float& number)
{
  ostringstream oss;
  oss << number;
  return oss.str();
}

saveall(  TString save_ans="")
{

 
      c1->SaveAs("LRcounts"+save_ans+".eps");
      c3->SaveAs("LRasym"+save_ans+".eps");
      c4->SaveAs("BTcounts"+save_ans+".eps");
      c6->SaveAs("BTasym"+save_ans+".eps");
      /* c1->SaveAs("LRcounts"+FS(implE)+save_ans+".gif");
      c3->SaveAs("LRasym"+FS(implE)+save_ans+".gif");
      c4->SaveAs("BTcounts"+FS(implE)+save_ans+".gif");
      c6->SaveAs("BTasym"+FS(implE)+save_ans+".gif");*/


}

void TOF_plot()
{
    Double_t a[2],b[2];
  
  a[0]=meanmuplus;a[1]=meanmuplus;
  b[0]=-1;b[1]=2500;
  g1= new TGraph(2,a,b);
  a[0]=meanMU;a[1]=meanMU;
  b[0]=-1;b[1]=2500;
  g2= new TGraph(2,a,b);
  g1->SetLineColor(1);
  g1->SetLineWidth(2);
  g1->Draw("");
  g2->SetLineColor(2);
  g2->SetLineWidth(2);
  g2->Draw("");
}


void GetValues()
{
 


 ifstream inf;
  inf.open("means.dat");
  while(!inf.is_open())
    {
      std::cout<<"No 'means.dat' file.\nAttempt to create it."<<std::endl;
      
      moyennes();
      inf.open("means.dat");
    }

  inf >> meanmuplus >> meanMU >> implE;// ie;
  inf.close();

}



Double_t MuonLifeFunc(Double_t *x, Double_t *par) {
  // two fit parameters: the normalisation
  //                     the muon life time
  return par[0]*TMath::Exp(-x[0]/par[1]);
}

Double_t TFasymuonfit(Double_t *x, Double_t *par) {
  // one fit parameter: the asymmetry
  Double_t fitval, arg;
  arg=-x[0]*Bfield*6.28*par[1]*1.e-6+3.14*par[2];
  fitval = par[0]*TMath::Cos(arg);
  return fitval;
}


Double_t TFdecayfit(Double_t *x, Double_t *par) {
 // 
 Double_t fitval;
 fitval = MuonLifeFunc(x, par);
 return fitval;
}


format()
{

  for(Int_t i==0; i<3; i++)
    {
      CASYM->cd(i)->SetFrameBorderMode(0);
      CASYM->cd(i)->SetFrameBorderSize(0);
      CASYM->cd(i)->SetGrid();
      CASYM->cd(i)->SetEditable(1);
   }


format(TCanvas cv)
{
      cv->SetFrameBorderMode(0);
      cv->SetFrameBorderSize(0);
      cv->SetGrid();
      cv->SetEditable(0);
   }

}