{
gStyle->SetOptFit(1111);
// ifstream timef("timing.txt");
// Double_t timing[132];
// for(Int_t i=0;i<132;i++){
//   timef>>timing[i];
// }
//
//   ifstream coeff("coef.txt");
//   Double_t coef[132];
//   Double_t val,v;
//   for(Int_t i=0;i<132;i++){
//     coeff>>val;coeff>>val;
//     coeff>>coef[i];
//     cout<<i<<" "<<coef[i]<<endl;
//   }
//   ifstream coeff2("coefnew.txt");
//   Double_t coefwf[132];
//   for(Int_t i=0;i<132;i++){
//     coeff2>>coefwf[i];
//     cout<<i<<" "<<coefwf[i]<<endl;
//   }
//gSystem->Load("~/franck/dvcs/soft/libDVCS.so");

struct event_dvcs {
  Float_t vz;
  Float_t xc;
  Float_t yc;
  Float_t nblock;
  Float_t xtheo;
  Float_t ytheo;
  Float_t e;
  Float_t ecalo;
  Float_t etheo;
  Float_t thetagtheo;
  Float_t thetag;
  Float_t phigtheo;
  Float_t phig;
  Float_t phi;
  Float_t pmom;
  Float_t pyef;
  Float_t hrsthetatheo;
};

event_dvcs evt;

TNtuple *ntu = new TNtuple("ntu","ntuple","vz:xc:yc:nblock:xtheo:ytheo:e:ecalo:etheo:thetagtheo:thetag:phigtheo:phig:phi:pmom:pyef:hrsthetatheo");

TDVCSEvent *event=new TDVCSEvent();
TH1F *h=new TH1F("h","h",100,-1,1);
// TH1F *htheta=new TH1F("htheta","htheta",100,-0.2,0.2);
// TH1F *hphi=new TH1F("hphi","hphi",100,-0.2,0.2);
TH1F *htheta=new TH1F("htheta","htheta",100,-2,2);
TH1F *hphi=new TH1F("hphi","hphi",100,-2,2);
TH1F *hthetab=new TH1F("hthetab","htheta",100,-20,20);
TH1F *hphib=new TH1F("hphib","hphi",100,-20,20);
TH1F *hthetaa=new TH1F("hthetaa","htheta",100,-20,20);
TH1F *hphia=new TH1F("hphia","hphi",100,-20,20);
TH2F *hXvsV=new TH2F("hXvsVz","hXvsV",100,-20,20,100,-20,20);
TH2F *hYvsV=new TH2F("hYvsVz","hYvsV",100,-20,20,100,-20,20);
TH2F *hXcvsV=new TH2F("hXcvsVz","hXcvsV",100,-20,20,100,-20,20);
TH2F *hYcvsV=new TH2F("hYcvsVz","hYcvsV",100,-20,20,100,-20,20);

TH2F *specphp=new TH2F("specphp","Phi vs P",100,0,6,100,-5,5);

TH2F *hXcvsXtheo=new TH2F("hXcvsXtheo","hXcvsYtheo",100,-20,20,100,-20,20);
TH2F *hYcvsYtheo=new TH2F("hYcvsYtheo","hYcvsYtheo",100,-20,20,100,-20,20);
TH2F *hYcvsXc=new TH2F("hYcvsXc","hYcvsXc",100,-20,20,100,-20,20);
TH1F *hResY=new TH1F("hResY","hResY",100,-20,20);
TH2F *hClusSizevsYc=new TH2F("hClusSizevsYc","hClusSizevsYc",100,-20,20,20,0,20);
TH2F *hEvsYc=new TH2F("hEvsYc","hEvsYc",200,-20,20,100,3.,5.);
TH2F *hEvsXc=new TH2F("hEvsXc","hEvsXc",200,-20,20,100,3.,5.);

// TH1F *hDthetaAn[15];
// for(Int_t ii=0;ii<15;ii++) {
//   TString hna="hDthetaAn";
//   hna+=ii;
//   char *toto=hna;
//   hDthetaAn[ii] = new TH1F(toto,toto,100,-0.02,0.02);
// }

TDVCSChain *tree=new TDVCSChain("T");
//tree->Add("/local/home/data/halla/elastic/dvcs_4998_0_elos.root");   //71855 Set 2
tree->Add("/local/home/data/halla/elastic/dvcs_4996_0_elos_vdcoff.root");   //71855 Set 2
//tree->Add("/local/home/data/halla/elastic/dvcs_3229_0.root");//Set 1
//tree->Add("/work/dvcs/carlos/elas/dvcs_3232.root");//60008 Set 2
///tree->Add("/work/dvcs/carlos/elas/dvcs_3234.root");// 39507 Set 3
//tree->Add("/work/dvcs/carlos/elas/dvcs_3235.root");//37227  Set 3
//tree->Add("/work/dvcs/carlos/elas/dvcs_3236.root");//Set 3

TFile *fcut=new TFile("cuts.root");  // machine adaq cd ~/dvcs/onlana
TCutG *cutt1=(TCutG*)fcut->Get("cut1");
TCutG *cutt2=(TCutG*)fcut->Get("cut2");
TCutG *cutt3=(TCutG*)fcut->Get("cut3");

TCaloEvent *caloev;
THRS *hrs=0;
THRS *hrseloss=0;
Double_t theta,phi,ene,ener,v;
Double_t p[100],ph[100],rth[100],rph[100],ry[100];
tree->SetBranchAddress("event_spectro.",&hrs);
tree->SetBranchAddress("event_spectroEloss.",&hrseloss);
tree->SetBranchAddress("event_calo.",&caloev);
tree->SetBranchAddress("L.tr.tg_ph",&ph);
tree->SetBranchAddress("L.tr.r_ph",&rph);
tree->SetBranchAddress("L.tr.r_th",&rth);
tree->SetBranchAddress("L.tr.r_y",&ry);
tree->SetBranchAddress("L.tr.p",&p);
tree->SetBranchAddress("rpl.z",&v);
Int_t nev=tree->GetEntries();
Int_t time,ccp=0;
nev=40000;

TARSWave *wt=new TARSWave();
// wt->SetFirstWindow(5,25,"data");           // for 1.1m data !!!
// wt->SetSecondWindow(5,25,-35,65,"data"); 
wt->SetFirstWindow(-15,15,"data");            // for 5.5m data !!!
wt->SetSecondWindow(-15,15,-50,50,"data"); 
wt->SetChi20(0.1);wt->SetChi21(180);wt->SetChi22(1e10);

//Double_t ebeam=5.767;
Double_t ebeam=5.7572;
TLorentzVector pini(0.,0.,0.,0.938272);
TLorentzVector eini(0.,0.,ebeam,ebeam);

TFile *g=new TFile("ntu.root","RECREATE");

for(Int_t i=0;i<nev;i++){


  if (!i) {
    cout << "Calo Distance from DB : " << event->GetGeometry()->GetCaloDist() << endl;
    cout << "Calo Angle from DB : " << event->GetGeometry()->GetCaloTheta()*TMath::RadToDeg() << endl;
  }

  Double_t calod=event->GetGeometry()->GetCaloDist();

  //  i+=10;
  if(i%100==0) cout<<i<<"/"<<nev<<endl;
  tree->GetEntry(i);
  ene=0,ener=0;
  Int_t nb=caloev->GetNbBlocks();

  specphp->Fill(p[0],ph[0]*TMath::RadToDeg(),1);

  for(Int_t j=0;j<nb;j++){
    TCaloBlock *block=caloev->GetBlock(j);
    TARSWave *wave=caloev->GetWave(j);
    wave->SubstractPedestalDB(block->GetBlockNumber(),"CALO"); 
    wave->SetNbChannel(block->GetBlockNumber(),"CALO");
    time=wave->GetSignal()->GetMinimumBin();
    Double_t amp=0.;
    wave->SetNbChannel(0,"CALO");
    wave->Analyze("data");
    block->Erase("");
    if(wave->GetNbPulses()>0) {
      amp=wave->GetAmplitude1();
    }
    block->SetBlockEnergy(amp*caloev->GetCalibration()->GetCalib(block->GetBlockNumber()));
    //   block->SetBlockEnergy(-(wave->GetSignal()->GetIntegral(40,75))*coef[block->GetBlockNumber()]/2976.);

  }

  TLorentzVector scatp=*hrseloss->GetParticle();
  //TLorentzVector scatp=*hrs->GetParticle();
  scatp.SetE(TMath::Sqrt(scatp.P()*scatp.P()+0.938272*0.938272));
  Double_t theo=(eini+pini-scatp).E();
  TLorentzVector theovec=(eini+pini-scatp);
  event->SetCaloEvent(caloev);
  //  v=v+0.10*rth[0]+0.15*rph[0]-0.1*ry[0];  // Bertin optimization (ELOG #12 on Software elog)
  event->SetVertex(0.,0.14,v*100.); // Vertex correction from Multifoil target analysis
                                       // among which: +1.05 from center foil position and
                                       //              -(-0.09) from center foil optics survey
  //  event->SetVertex(0.,0.,v*100.);
  //    event->SetVertex(0.,0.,0.);
  // cout<<v*100.<<endl;
  caloev->DoClustering();
  //  h->Fill(caloev->GetEnergy()-theo);
  if(caloev->GetNbClusters()==1){

    caloev->GetCluster(0)->Analyze(1,4.3); // Adjustment of Wo

    Double_t xclus=caloev->GetCluster(0)->GetX();
    Double_t yclus=caloev->GetCluster(0)->GetY();

    TLorentzVector photon=event->GetPhoton(0,7.,0,1.,1.0); // Adjustment of A and B
    Double_t thetatheo=atan2(theovec.Py(),sqrt(pow(theovec.Px(),2.)+pow(theovec.Pz(),2.)));
    Double_t phitheo=atan2(theovec.Px(),theovec.Pz());
    Double_t theta=atan2(photon.Py(),sqrt(pow(photon.Px(),2.)+pow(photon.Pz(),2.)));
    Double_t phi=atan2(photon.Px(),photon.Pz());

    //    if ( (yclus>-15.)&&(yclus<15.)&&(xclus>-10.)&&(xclus<11.)&&(cutt1->IsInside(ph[0],p[0])>0) && (caloev->GetCluster(0)->GetClusSize()==9)) {
    if ( (cutt2->IsInside(ph[0],p[0])>0)&&((-9.<v*100.)&&(v*100.<9.))&&(xclus>-5)&&(xclus<5)&&(yclus>-5)&&(yclus<5) ) {
      //    if ( (yclus>-15.)&&(yclus<15.)&&(xclus>-10.)&&(xclus<11.) ) {
      h->Fill(caloev->GetCluster(0)->GetEnergy()-theo);
      hthetab->Fill(calod*tan(thetatheo));
      hphib->Fill(calod*tan(phitheo+0.3577401));
      hthetaa->Fill(calod*tan(theta));
      hphia->Fill(calod*tan(phi+0.3577401));
      htheta->Fill((theta-thetatheo)*TMath::RadToDeg());
      hphi->Fill((phi-phitheo)*TMath::RadToDeg());

      hXvsV->Fill(v*100,calod*tan(phi+0.3577401));
      hYvsV->Fill(v*100,calod*tan(theta));
      hXcvsV->Fill(v*100,xclus);
      hYcvsV->Fill(v*100,yclus);
      hYcvsXc->Fill(xclus,yclus);

      hYcvsYtheo->Fill(calod*tan(thetatheo),yclus);
      hXcvsXtheo->Fill(calod*tan(phitheo+0.3577401),xclus);

      hResY->Fill(calod*tan(thetatheo)-yclus);
      hClusSizevsYc->Fill(yclus,float(caloev->GetCluster(0)->GetClusSize()));
      hEvsYc->Fill(yclus,caloev->GetEnergy());
      hEvsXc->Fill(xclus,caloev->GetEnergy());

      evt.vz=float(v*100.);
      evt.xc=float(xclus);
      evt.yc=float(yclus);
      evt.nblock=float(caloev->GetCluster(0)->GetClusSize());
      evt.xtheo=float(calod*tan(phitheo+0.3577401));
      evt.ytheo=float(calod*tan(thetatheo));
      evt.e=float(caloev->GetCluster(0)->GetEnergy());
      evt.ecalo=float(caloev->GetEnergy());
      evt.etheo=float(theo);
      evt.thetag=float(atan2(photon.Py(),sqrt(pow(photon.Px(),2.)+pow(photon.Pz(),2.))));
      evt.thetagtheo=float(atan2(theovec.Py(),sqrt(pow(theovec.Px(),2.)+pow(theovec.Pz(),2.))));
      evt.phig=float(atan2(photon.Px(),photon.Pz()));
      evt.phigtheo=float(atan2(theovec.Px(),theovec.Pz()));
      evt.phi=float(ph[0]);
      //      evt.pmom=float(scatp.P());
      evt.pmom=float(p[0]);
      evt.pyef=float(scatp.Py());

      evt.hrsthetatheo=TMath::ACos((ebeam+0.938272)*(TMath::Sqrt(evt.pmom*evt.pmom+0.938272*0.938272)-0.938272)/(ebeam*evt.pmom));

      ntu->Fill(&evt.vz);

     }

    // cout<<theta<<" "<<phi<<" "<<thetatheo<<" "<<phitheo<<endl;

  }
  caloev->Reset();
}

g->cd();
ntu->Write();

TCanvas toto;
toto.Divide(2,2);
toto.cd(1);
ntu->Draw("(phigtheo-phig)","abs(phigtheo-phig)<0.02");
toto.cd(2);
ntu->Draw("(thetagtheo-thetag)","abs(thetagtheo-thetag)<0.02");
toto.cd(3);
ntu->Draw("(etheo-e)","abs(etheo-e)<1");
toto.cd(4);
ntu->Draw("(phigtheo-phig):vz","abs(phigtheo-phig)<0.02");
g->Close();

}