#include <iostream>
#include <cstdlib>
#include <cmath>
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TROOT.h"
#include "TCutG.h"
#include "TNtuple.h"
#include "TF1.h"
#include "TF2.h"
#include <string>
#include "Riostream.h"
using namespace std;
float ar(float y, float x)
{
  float angolo=0;
  double rad = 0.0174532952;
  if(x>0) angolo = atan(y/x);
  if(x<0 && y>0) angolo = atan(y/x) + 180*rad;
  if(x<0 && y<0) angolo = atan(y/x) - 180*rad;
  if(x>0 && y==0) angolo = 0;
  if(x==0 && y>0) angolo = 90*rad;
  if(x<0 && y==0) angolo = 180*rad;
  if(x==0 && y<0) angolo = -90*rad;
  return(angolo);
}

int main() {

TChain *thm = new TChain("thm","");

thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag95.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag96.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag97.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag98.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag99.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag100.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag101.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag102.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag103.root/thm");
thm->Add("/home/aleksandra/zagreb2017_runs/root_data/zag104.root/thm");

//----------------------------------------------------------------------------
//	Declaration of leaves types
//----------------------------------------------------------------------------

   Double_t           pu1;
   Double_t           pu2;
   Double_t           pd1;
   Double_t           pd2;
   Double_t           eu1;
   Double_t           eu2;
   Double_t           ed1;
   Double_t           ed2;
   Double_t           nst;

//----------------------------------------------------------------------------
//	Set branch addresses.
//----------------------------------------------------------------------------

   thm->SetBranchAddress("pu1",&pu1);
   thm->SetBranchAddress("pu2",&pu2);
   thm->SetBranchAddress("pd1",&pd1);
   thm->SetBranchAddress("pd2",&pd2);
   thm->SetBranchAddress("eu1",&eu1);
   thm->SetBranchAddress("eu2",&eu2);
   thm->SetBranchAddress("ed1",&ed1);
   thm->SetBranchAddress("ed2",&ed2);
   thm->SetBranchAddress("nst",&nst);

struct data_t{
Float_t pu1x;
Float_t pu2x;
Float_t pd1x;
Float_t pd2x;
Float_t eu1x;
Float_t eu2x;
Float_t ed1x;
Float_t ed2x;
Float_t nstx;
Float_t thetad2x;
Float_t thetau2x;
Float_t ened2x;
Float_t eneu2x;
Float_t decarbond2x;
Float_t ecd2x;
Float_t dehalftargetd2x;
Float_t ead2x;
Float_t dehalftargetu2x;
Float_t eau2x;
Float_t ximp1x;
Float_t ximp2x;
Float_t ximp3x;
Float_t espx;
Float_t theta3x;
Float_t etotx;
Float_t q3x;
Float_t pa1v3x;
Float_t ea1v3x;
Float_t ximp1v3x;
Float_t pa2v3x;
Float_t ea2v3x;
Float_t ximp2v3x;
Float_t theta3v3x;
Float_t ximp3v3x;
Float_t espv3x;
Float_t q3v3x;
Float_t e12x;
Float_t e13x;
Float_t e23x;
Float_t e12v3x;
Float_t e13v3x;
Float_t e23v3x;
Float_t t12x;
Float_t romxx;
Float_t romyx;
};
data_t data;

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

	Long64_t nentries = thm->GetEntries();
	cout << nentries << "\n";

	new TFile("ks.root","RECREATE");
	TTree *run = new TTree("run"," ");

run->Branch("pu1x",&data.pu1x,"pu1x/F");
run->Branch("pu2x",&data.pu2x,"pu2x/F");
run->Branch("pd1x",&data.pd1x,"pd1x/F");
run->Branch("pd2x",&data.pd2x,"pd2x/F");
run->Branch("eu1x",&data.eu1x,"eu1x/F");
run->Branch("eu2x",&data.eu2x,"eu2x/F");
run->Branch("ed1x",&data.ed1x,"ed1x/F");
run->Branch("ed2x",&data.ed2x,"ed2x/F");
run->Branch("nstx",&data.nstx,"nstx/F");
run->Branch("thetad2x",&data.thetad2x,"thetad2x/F");
run->Branch("thetau2x",&data.thetau2x,"thetau2x/F");   
run->Branch("ened2x",&data.ened2x,"ened2x/F");
run->Branch("eneu2x",&data.eneu2x,"eneu2x/F");
run->Branch("decarbond2x",&data.decarbond2x,"decarbond2x/F");
run->Branch("ecd2x",&data.ecd2x,"ecd2x/F");
run->Branch("dehalftargetd2x",&data.dehalftargetd2x,"dehalftargetd2x/F");
run->Branch("ead2x",&data.ead2x,"ead2x/F");
run->Branch("dehalftargetu2x",&data.dehalftargetu2x,"dehalftargetu2x/F");
run->Branch("eau2x",&data.eau2x,"eau2x/F");
run->Branch("ximp1x",&data.ximp1x,"ximp1x/F");
run->Branch("ximp2x",&data.ximp2x,"ximp2x/F");
run->Branch("ximp3x",&data.ximp3x,"ximp3x/F");
run->Branch("espx",&data.espx,"espx/F");
run->Branch("theta3x",&data.theta3x,"theta3x/F");
run->Branch("etotx",&data.etotx,"etotx/F");
run->Branch("q3x",&data.q3x,"q3x/F");
run->Branch("pa1v3x",&data.pa1v3x,"pa1v3x/F");
run->Branch("ea1v3x",&data.ea1v3x,"ea1v3x/F");
run->Branch("ximp1v3x",&data.ximp1v3x,"ximp1v3x/F");
run->Branch("pa2v3x",&data.pa1v3x,"pa1v3x/F");
run->Branch("ea2v3x",&data.ea1v3x,"ea1v3x/F");
run->Branch("ximp2v3x",&data.ximp1v3x,"ximp1v3x/F");
run->Branch("theta3v3x",&data.theta3v3x,"theta3v3x/F");
run->Branch("ximp3v3x",&data.ximp3v3x,"ximp3v3x/F");
run->Branch("espv3x",&data.espv3x,"espv3x/F");
run->Branch("q3v3x",&data.q3v3x,"q3v3x/F");
run->Branch("e12x",&data.e12x,"e12x/F");
run->Branch("e13x",&data.e13x,"e13x/F");
run->Branch("e23x",&data.e23x,"e23x/F");
run->Branch("e12v3x",&data.e12v3x,"e12v3x/F");
run->Branch("e13v3x",&data.e13v3x,"e13v3x/F");
run->Branch("e23v3x",&data.e23v3x,"e23v3x/F");
run->Branch("t12x",&data.t12x,"t12x/F");
run->Branch("romxx",&data.romxx,"romxx/F");
run->Branch("romyx",&data.romyx,"romyx/F");

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------
//	PARAMETERS
//----------------------------------------------------------------------------

Double_t u=931.49432;	//atomic mass unit (rif. 12C)
Double_t xmpro=3.0160293*u;	//massa proiettile (3He)
Double_t xmtar=9.0121822*u;	//massa target (9Be)
Double_t xmal1=4.002602*u;	//massa (alpha1)
Double_t xmal2=4.002602*u;	//massa  (alpha2)
Double_t xmsp=4.002602*u;	//massa spettatore (alpha3)
Double_t epro=3.9344;	//en.fin. di 3He@ 4 MeV dopo 1/2 target di 9Be da 124 ug/cm2
Double_t xmtrs=5.01222*u;	// massa trasferito (5He)

Double_t deg=180./3.1415;
Double_t rad=1./deg;

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------
//	CALIBRAZIONE 		   
//----------------------------------------------------------------------------

//////////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------
// 	Calibrazione d2 
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
//	PARAMETRI
//----------------------------------------------------------------------------

Double_t mad2=0.06065;
Double_t qad2=-3.89114;
Double_t mbd2=-0.95537;
Double_t qbd2=15.24339;

Double_t aad2=-0.11406;
Double_t bad2=0;
Double_t abd2=0.00775;
Double_t bbd2=0;
Double_t acd2=0;
Double_t bcd2=0;

Double_t dec0=0.10489;
Double_t dec1=0.24323;
Double_t dec2=2.65652;
Double_t dec3=-0.17757;
Double_t dec4=0.58807;
Double_t dec5=-1.46695;

Double_t deht0=0.15178;
Double_t deht1=0.14815;
Double_t deht2=1.16162;
Double_t deht3=-0.63179;
Double_t deht4=1.63179;
Double_t deht5=-0.07006;

//----------------------------------------------------------------------------
//	FUNZIONI
//----------------------------------------------------------------------------

TF2 *thd21= new TF2("thd21","(x-[0]*y-[2])/([1]*y+[3])",0.,4000., 400., 4000.);
	//x=pd2, y=ed2
TF2 *ened21= new TF2("ened21","([0]+[1]*x+([2]+[3]*x)*y+([4]+[5]*x)*y*y)",40.,60., 400., 4000.);
	//x=thetad2, y=ed2

TF2 *decarbond21 = new TF2("decarbond21","[0]*pow(x,[1])-(exp(-[2]*pow(x,[3])-[4]*pow(x,[5])))/cos(y*0.0174532952-45*0.0174532952)",0.,40.,35.,60.);
	//dec0=[0], dec1=[1], dec2=[2], dec3=[3], dec4=[4], dec5=[5], x=ened21, y=thd21
TF2 *ecd21 = new TF2("ecd21","x+y",0.,40,0.,40);
	//x=decarbond21, y=ened21

TF2 *dehalftargetd21 = new TF2("dehalftargetd21","[0]*pow(x,[1])-(exp(-[2]*pow(x,[3])-[4]*pow(x,[5])))/cos(y*0.0174532952-45*0.0174532952)",0.,40.,35.,60.);
	//deht0=[0], deht1=[1], deht2=[2], deht3=[3], deht4=[4], deht5=[5], x=ecd21, y=thd21
TF2 *ead21 = new TF2("ead21","x+y",0.,40,0.,40);
	//x=dehalftargetd21, y=ecd21

//----------------------------------------------------------------------------
//	ASSEGNAZIONE DEI PARAMETRI ALLE FUNZIONI
//----------------------------------------------------------------------------

thd21->FixParameter(0,qad2);
thd21->FixParameter(1,mad2);
thd21->FixParameter(2,qbd2);
thd21->FixParameter(3,mbd2);

ened21->FixParameter(0,aad2);
ened21->FixParameter(1,bad2);
ened21->FixParameter(2,abd2);
ened21->FixParameter(3,bbd2);
ened21->FixParameter(4,acd2);
ened21->FixParameter(5,bcd2);

decarbond21->FixParameter(0,dec0);
decarbond21->FixParameter(1,dec1);
decarbond21->FixParameter(2,dec2);
decarbond21->FixParameter(3,dec3);
decarbond21->FixParameter(4,dec4);
decarbond21->FixParameter(5,dec5);

dehalftargetd21->FixParameter(0,deht0);
dehalftargetd21->FixParameter(1,deht1);
dehalftargetd21->FixParameter(2,deht2);
dehalftargetd21->FixParameter(3,deht3);
dehalftargetd21->FixParameter(4,deht4);
dehalftargetd21->FixParameter(5,deht5);

//////////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------
//	Calibrazione u2
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
//	PARAMETRI
//----------------------------------------------------------------------------

Double_t mau2=0.06605;
Double_t qau2=-4.25284;
Double_t mbu2=-0.8453;
Double_t qbu2=-1.18903;

Double_t aau2=-0.11406;
Double_t bau2=0;
Double_t abu2=0.00775;
Double_t bbu2=0;
Double_t acu2=0;
Double_t bcu2=0;

//----------------------------------------------------------------------------
//	FUNZIONI
//----------------------------------------------------------------------------

TF2 *thu22= new TF2("thu22","(x-[0]*y-[2])/([1]*y+[3])",0.,4000., 0., 4000.);
	//x=pu2, y=eu2
TF2 *eneu22= new TF2("eneu22","([0]+[1]*x+([2]+[3]*x)*y+([4]+[5]*x)*y*y)",0.,4000., 0., 4000.);
	//x=thetau2, y=eu2

TF2 *dehalftargetu22 = new TF2("dehalftargetu22","[0]*pow(x,[1])-(exp(-[2]*pow(x,[3])-[4]*pow(x,[5])))/cos(y*0.0174532952+45*0.0174532952)",0.,60.,40.,90.);
	//deht0=[0], deht1=[1], deht2=[2], deht3=[3], deht4=[4], deht5=[5], x=eneu22, y=thu22
TF2 *eau22 = new TF2("eau22","x+y",0.,40.,0.,60.);
	//x=dehalftargetu22, x=eneu22

//----------------------------------------------------------------------------
//	ASSEGNAZIONE DEI PARAMETRI ALLE FUNZIONI
//----------------------------------------------------------------------------

thu22->FixParameter(0,qau2);
thu22->FixParameter(1,mau2);
thu22->FixParameter(2,qbu2);
thu22->FixParameter(3,mbu2);

eneu22->FixParameter(0,aau2);
eneu22->FixParameter(1,bau2);
eneu22->FixParameter(2,abu2);
eneu22->FixParameter(3,bbu2);
eneu22->FixParameter(4,acu2);
eneu22->FixParameter(5,bcu2);

dehalftargetu22->FixParameter(0,deht0);
dehalftargetu22->FixParameter(1,deht1);
dehalftargetu22->FixParameter(2,deht2);
dehalftargetu22->FixParameter(3,deht3);
dehalftargetu22->FixParameter(4,deht4);
dehalftargetu22->FixParameter(5,deht5);

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

Double_t cos1,cos2,sin1,sin2,cf1,cf2,sf1,sf2,cf12,c12,pp,ximp1,ximp2,pp1x,pp1z,pp2x,pp2z,pp3x,pp3z,pp1y,pp2y,pp3y;
Double_t vv,aa1,aa2,aa3,pc1x,pc1y,pc1z,pc1,pc2x,pc2y,pc2z,pc2,pc3x,pc3y,pc3z,pc3,ec1,ec2,ec3,fi3,theta3,ximp3,esp,etotdata,etot,q3_ec,q3;
Double_t qvalorea,aeq1,beq1,alfam1,qvalore,betam1,ceq1,delta1,sdelta1,pu201,vartest1,pa1v3,ea1v3;
Double_t aeq2,beq2,alfam2,betam2,ceq2,delta2,sdelta2,pu202,vartest2,pa2v3,ea2v3;
Double_t ximp1v3,pp1xv3,pp1zv3,pp3xv3,pp3zv3,theta3v3,ximp3v3,p3fv3,espv3,pc1xv3,pc1yv3,pc1zv3,pc1v3,pc3xv3,pc3yv3,pc3zv3,pc3v3,ec1v3,ec3v3;
Double_t ximp2v3,pp2xv3,pp2zv3,pc2xv3,pc2yv3,pc2zv3,pc2v3,ec2v3;
Double_t etotdatav3,etotv3,q3_ecv3,q3v3;
Double_t xm,e12,e13,e23;
Double_t e12v3,e13v3,e23v3;
Double_t tt,ttt,ctt,stt,vp,v1,v2,vt,xnumx,xnumy,xnum,xdena,xdenb,xden,argo,t12;
Double_t romx,romy;

	Long64_t nbytes = 0;
   for (Long64_t i=0; i<nentries;i++) 
	{
//	thm->cd();
	nbytes += thm->GetEntry(i);

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// thm->SetBranchStatus("*",0);  // disable all branches
// TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname


Double_t thetad2=thd21->Eval(pd2,ed2);
Double_t thetau2=thu22->Eval(pu2,eu2);
Double_t ened2=ened21->Eval(thetad2,ed2);
Double_t eneu2=eneu22->Eval(thetau2,eu2);
Double_t decarbond2=decarbond21->Eval(ened2,thetad2);
Double_t ecd2=ecd21->Eval(decarbond2,ened2);
Double_t dehalftargetd2=dehalftargetd21->Eval(ecd2,thetad2);
Double_t ead2=ead21->Eval(dehalftargetd2,ecd2);
Double_t dehalftargetu2=dehalftargetu22->Eval(eneu2,thetau2);
Double_t eau2=eau22->Eval(dehalftargetu2,eneu2);


//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------
// 	KINEMATICS
//----------------------------------------------------------------------------

	//ed2 = alpha1 -> part "1"
	//eu2 = alpha2 -> part "2"
	//dunque "3" sarà lo spettatore

	cos1=cos(thetad2*rad);
	cos2=cos(thetau2*rad);
	sin1=sin(thetad2*rad);
	sin2=sin(thetau2*rad);

	cf1=1.;
	cf2=-1.;
	sf1=0.;
	sf2=0.;
	cf12=-1.;
	c12=cos1*cos2+sin1*sin2*cf12;

        pp=sqrt(2.*xmpro*epro);
        ximp1=sqrt(2.*ead2*xmal1);
        ximp2=sqrt(2.*eau2*xmal2);

        //alpha1
        pp1x=ximp1*sin1;
        pp1z=ximp1*cos1;

        //alpha2
        pp2x=-ximp2*sin2;
        pp2z=ximp2*cos2;

        //alpha spettatore
        pp3x=-pp1x-pp2x;
        pp3z=pp-pp1z-pp2z;

        pp1y=ximp1*sin1*sf1;
        pp2y=-ximp2*sin2*sf2;

        pp3y=-pp1y-pp2y;

//----------------------------------------------------------------------------
//	calcolo dell'energia totale di ciascuna particella
//	nel centro di massa di 1+2+3
//	calcolo del q-valore
// ----------------------------------------------------------------------------

	//calcolo delle componenti di p1,p2,p3 nel c.m. di 1,2,3
        vv=pp/(xmpro+xmtar);

	//vv è quella che nell' ohlsen si chiama V, velocità del centro di massa
	aa1=vv*sqrt(xmal1/2.);        
	aa2=vv*sqrt(xmal2/2.);
        aa3=vv*sqrt(xmsp/2.);

	//a1,a2,a3 sono a1,a2,a3 dell'ohlsen
        pc1x=ximp1*sin1*cf1;
        pc1y=ximp1*sin1*sf1;
        pc1z=ximp1*cos1-aa1*sqrt(2.*xmal1);
        pc1=sqrt(pc1x*pc1x+pc1y*pc1y+pc1z*pc1z);

	//facendo i conti questo coincide con quanto riporta l'ohlsen alla formula (3)
        pc2x=ximp2*sin2*cf2;
        pc2y=ximp2*sin2*sf2;
        pc2z=ximp2*cos2-aa2*sqrt(2.*xmal2);
        pc2=sqrt(pc2x*pc2x+pc2y*pc2y+pc2z*pc2z);

        pc3x=-pc1x-pc2x;
        pc3y=-pc1y-pc2y;
        pc3z=-pc1z-pc2z;
	pc3=sqrt(pc3x*pc3x+pc3y*pc3y+pc3z*pc3z);

        ec1=pc1*pc1/(2.*xmal1);
        ec2=pc2*pc2/(2.*xmal2);
        ec3=pc3*pc3/(2.*xmsp);
        
        fi3=deg*ar(pp3y,pp3x);
        theta3=deg*ar(pc3x,(pc3z+xmsp*vv));
        
        if (fi3==0) {
        ximp3=sqrt((pp3x*pp3x)+(pp3z*pp3z));
        }
        else
        {
        ximp3=-sqrt((pp3x*pp3x)+(pp3z*pp3z));
        }
        
        esp=(ximp3*ximp3)/(2.*xmsp);

        etotdata=ead2+eau2+esp;		//energia totale a disp. nel c.m. 1,2,3
        etot=ec1+ec2+ec3;
	
        q3_ec=ec1+ec2+ec3-epro;
        q3=ead2+eau2+esp-epro;

//----------------------------------------------------------------------------
// 	Calcolo a tre variabili. Escludo energia alpha1 su D2    
//----------------------------------------------------------------------------
                                       
	qvalorea=19.005;                                         
        
	aeq1=0.5*((xmsp+xmal1)/xmal1);                             
	beq1=sqrt(2.*xmal2*eau2)*c12-sqrt(2.*xmpro*epro)*cos1; 
	alfam1=sqrt(xmal2*eau2*xmpro*epro);                    
	if (q3<19.2 && q3>18.75) {
	qvalore=qvalorea;
	}
	else
	{
	qvalore=0;
	}
	betam1=xmsp*(eau2-epro-qvalore);                      
	ceq1=xmpro*epro+xmal2*eau2-2.*alfam1*cos2+betam1;        
	delta1=beq1*beq1-4.*aeq1*ceq1;                               
	sdelta1=sqrt(fabs(delta1));                              
	pu201=sqrt(2.*xmal1*ead2);                                
	vartest1=2*aeq1*pu201+beq1;                                 
        
	if(vartest1>=0.)                                         
	{                                                       
	pa1v3=(-beq1+sdelta1)/(2.*aeq1);                       
	ea1v3=(pa1v3*pa1v3)/(2.*xmal1);                   
	}                                                     
	else                                                
	{                                                       
	pa1v3=(-beq1-sdelta1)/(2.*aeq1);                       
	ea1v3=(pa1v3*pa1v3)/(2.*xmal1);                     
	}                                    

//----------------------------------------------------------------------------
// 	Calcolo a tre variabili. Escludo energia alpha2 su U2    
//----------------------------------------------------------------------------                                        
                                          
	aeq2=0.5*((xmsp+xmal2)/xmal2);                             
	beq2=sqrt(2.*xmal1*ead2)*c12-sqrt(2.*xmpro*epro)*cos2; 
	alfam2=sqrt(xmal1*ead2*xmpro*epro);                    
	betam2=xmsp*(ead2-epro-qvalore);                      
	ceq2=xmpro*epro+xmal1*ead2-2.*alfam2*cos1+betam2;        
	delta2=beq2*beq2-4.*aeq2*ceq2;                               
	sdelta2=sqrt(fabs(delta2));                              
	pu202=sqrt(2.*xmal2*eau2);                                
	vartest2=2*aeq2*pu202+beq2;                                 
        
	if(vartest2>=0.)                                         
	{                                                       
	pa2v3=(beq2+sdelta2)/(2.*aeq2);                       
	ea2v3=(pa2v3*pa2v3)/(2.*xmal2);                   
	}                                                     
	else                                                
	{                                                       
	pa2v3=(beq2-sdelta2)/(2.*aeq2);                       
	ea2v3=(pa2v3*pa2v3)/(2.*xmal2);                     
	}                           
 
//----------------------------------------------------------------------------       
//	3v
//----------------------------------------------------------------------------  
      
	ximp1v3=sqrt(2.*ea1v3*xmal1);
	pp1xv3=ximp1v3*sin1;
	pp1zv3=ximp1v3*cos1;
	ximp2v3=sqrt(2.*ea2v3*xmal2);
	pp2xv3=ximp2v3*sin2;
	pp2zv3=ximp2v3*cos2;
	pp3xv3=-pp1xv3-pp2xv3;
	pp3zv3=pp-pp1zv3-pp2zv3;
        
	if (pp3xv3<=0)
	{
	theta3v3=-ar(pp3zv3,pp3xv3);
	ximp3v3=-sqrt((pp3xv3*pp3xv3)+(pp3zv3*pp3zv3));
	p3fv3=-sqrt((pp3zv3*pp3zv3)+(pp3xv3*pp3xv3));
	}
	else
	{
	theta3v3=ar(pp3zv3,pp3xv3);
	p3fv3=sqrt((pp3zv3*pp3zv3)+(pp3xv3*pp3xv3));
	ximp3v3=sqrt((pp3xv3*pp3xv3)+(pp3zv3*pp3zv3));
	}
	espv3=(ximp3v3*ximp3v3)/(2.*xmsp);
        
	pc1xv3=ximp1v3*sin1*cf1;
	pc1yv3=ximp1v3*sin1*sf1;
	pc1zv3=ximp1v3*cos1-aa1*sqrt(2.*xmal1);
	pc1v3=sqrt(pc1xv3*pc1xv3+pc1yv3*pc1yv3+pc1zv3*pc1zv3);

	pc2xv3=ximp2v3*sin2*cf2;
	pc2yv3=ximp2v3*sin2*sf2;
	pc2zv3=ximp2v3*cos2-aa2*sqrt(2.*xmal2);
	pc2v3=sqrt(pc2xv3*pc2xv3+pc2yv3*pc2yv3+pc2zv3*pc2zv3);
        
	pc3xv3=-pc1xv3-pc2xv3;
	pc3yv3=-pc1yv3-pc2yv3;
	pc3zv3=-pc1zv3-pc2zv3;
	pc3v3=sqrt(pc3xv3*pc3xv3+pc3yv3*pc3yv3+pc3zv3*pc3zv3);
             
	ec1v3=pc1v3*pc1v3/(2.*xmal1);
	ec2v3=pc2v3*pc2v3/(2.*xmal2);
	ec3v3=pc3v3*pc3v3/(2.*xmsp);
	etotdatav3=ea1v3+ea2v3+espv3;	//energia totale a disp. nel c.m. 1,2,3
	etotv3=ec1v3+ec2v3+ec3v3;
        
	q3_ecv3=ec1v3+ec2v3+ec3v3-epro;
	q3v3=ea1v3+ea2v3+espv3-epro;

//---------------------------------------------------------------------------
//	CALCOLO DELLE ENERGIE RELATIVE e12,e13,e23
//---------------------------------------------------------------------------

	xm=xmal1+xmal2+xmsp;

	e12=etot-xm*ec3/(xmal1+xmal2);
  	e13=etot-xm*ec2/(xmal1+xmsp);
	e23=etot-xm*ec1/(xmal2+xmsp);

//	3v
    
	e12v3=etotv3-xm*ec3v3/(xmal1+xmal2);
	e13v3=etotv3-xm*ec2v3/(xmal1+xmsp);
	e23v3=etotv3-xm*ec1v3/(xmal2+xmsp);

//---------------------------------------------------------------------------
//	Calcolo theta_cm 4 variabili
//---------------------------------------------------------------------------	

	Double_t pi=3.1415; 

	tt=abs(theta3);
	if(theta3<=0.)
	{
	ttt=(pi-tt*rad);	//tt angolo della particella trasferita
	}	
	else
	{	
	ttt=-(pi-tt*rad);
	}
        ctt=cos(ttt);
        stt=sin(ttt);

	//b-u proj	//costruzione dell'argomento del theta_cm (slaus)
//	vp=sqrt(2.*epro/xmproi);
//	v1=sqrt(2.*e3li7/xmli);
//	v2=sqrt(2.*e4al/xmal);

//	vtzc=((xmsp+xmtrs)*vp/xmtrs)-pp3z/xmtrs;
//	vtx=-pp3x/xmtrs;

//	xnumz=vtzc*(v1*cos1-v2*cos2);
//	xnumx=vtx*(v1*sin1-v2*sin2);
//	xnum=xnumx+xnumz;

//	xdena=sqrt(vtzc*vtzc+vtx*vtx);
//	xdenb=sqrt((v1*cos1-v2*cos2)*(v1*cos1-v2*cos2)+(v1*sin1-v2*sin2)*(v1*sin1-v2*sin2));
        
//	xden=xdena*xdenb;
//	argo=xnum/xden;


	//b-u target	//costruzione dell'argomento del theta_cm (slaus)
	vp=sqrt(2.*epro/(xmpro));
	v1=sqrt(2.*ead2/xmal1);
	v2=sqrt(2.*eau2/xmal2);
	vt=fabs(ximp3)/(xmtrs);

	xnumx=(vp-vt*ctt)*(v1*cos1-v2*cos2);
	xnumy=(-vt*stt)*(v1*sin1+v2*sin2);
	xnum=xnumx+xnumy;

	xdena=sqrt((vp-vt*ctt)*(vp-vt*ctt)+(-vt*stt)*(-vt*stt));
	xdenb=sqrt((v1*cos1-v2*cos2)*(v1*cos1-v2*cos2)+(v1*sin1+v2*sin2)*(v1*sin1+v2*sin2));
	xden=xdena*xdenb;

	argo=xnum/xden;

	if(argo>=-1 || argo<=1)
	{
	t12=(acos(xnum/xden))*deg;	//angolo cm. della particella 2
//	t12=180.-acos(xnum/xden)*deg;	//angolo cm. della particella 1
	}
	else
	{
	t12=0.;
	}

// -----------------------------------------------------------------------
//	calcolo romx
//	calcolo romy
//      Costanzo et al.: NIM A295 (1990) 373
// -----------------------------------------------------------------------

       romx=(ximp3*ximp3)/(2.*u);
       romy=epro-ead2-eau2;

//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------------
//	g->cd();
//---------------------------------------------------------------------------

data.pu1x=pu1;
data.pu2x=pu2;
data.pd1x=pd1;
data.pd2x=pd2;
data.eu1x=eu1;
data.eu2x=eu2;
data.ed1x=ed1;
data.ed2x=ed2;
data.nstx=nst;
data.thetad2x=thetad2;
data.thetau2x=thetau2;
data.ened2x=ened2;
data.eneu2x=eneu2;
data.decarbond2x=decarbond2;
data.ecd2x=ecd2;
data.dehalftargetd2x=dehalftargetd2;
data.ead2x=ead2;
data.dehalftargetu2x=dehalftargetu2;
data.eau2x=eau2;
data.ximp1x=ximp1;
data.ximp2x=ximp2;
data.ximp3x=ximp3;
data.espx=esp;
data.theta3x=theta3;
data.etotx=etot;
data.q3x=q3;
data.pa1v3x=pa1v3;
data.ea1v3x=ea1v3;
data.ximp1v3x=ximp1v3;
data.pa2v3x=pa1v3;
data.ea2v3x=ea1v3;
data.ximp2v3x=ximp1v3;
data.theta3v3x=theta3v3;
data.ximp3v3x=ximp3v3;
data.espv3x=espv3;
data.q3v3x=q3v3;
data.e12x=e12;
data.e13x=e13;
data.e23x=e23;
data.e12v3x=e12v3;
data.e13v3x=e13v3;
data.e23v3x=e23v3;
data.t12x=t12;
data.romxx=romx;
data.romyx=romy;

run->Fill();

}

run->Print();
run->Write();

exit(0);

}