#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() {

gROOT->Reset();
TChain *run1 = new TChain("run1","");

run1->Add("/home/aleksandra/b10-runs/b10-136.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-137.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-138.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-139.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-140.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-141.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-142.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-143.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-154.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-155.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-156.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-157.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-159.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-160.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-161.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-163.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-164.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-165.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-166.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-173.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-174.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-175.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-176.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-177.root/run1");

run1->Add("/home/aleksandra/b10-runs/b10-178.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-179.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-180.root/run1");
run1->Add("/home/aleksandra/b10-runs/b10-181.root/run1");
/////////////////////////////////////////

//Declaration of leaves types

   Double_t           p1;
   Double_t           p2;
   Double_t           p3;
   Double_t           p4;
   Double_t           de1;
   Double_t           de2;
   Double_t           e1;
   Double_t           e2;
   Double_t           e3;
   Double_t           e4;
   Double_t           run_n;

// Set branch addresses.
   run1->SetBranchAddress("p1",&p1);
   run1->SetBranchAddress("p2",&p2);
   run1->SetBranchAddress("p3",&p3);
   run1->SetBranchAddress("p4",&p4);
   run1->SetBranchAddress("de1",&de1);
   run1->SetBranchAddress("de2",&de2);
   run1->SetBranchAddress("e1",&e1);
   run1->SetBranchAddress("e2",&e2);
   run1->SetBranchAddress("e3",&e3);
   run1->SetBranchAddress("e4",&e4);
   run1->SetBranchAddress("run_n",&run_n);

struct data_t{
Float_t p1x;
Float_t p2x;
Float_t p3x;
Float_t p4x;
Float_t de1x;
Float_t de2x;
Float_t e1x;
Float_t e2x;
Float_t e3x;
Float_t e4x;
Float_t run_nx;
Float_t theta1x;
Float_t ene1x;
Float_t r4delimy15x;
Float_t ar4delimy15x;
Float_t r4deliisox;
Float_t ar4deliisox;
Float_t r4delimy152x;
Float_t ar4delimy152x;
Float_t r4delitarx;
Float_t e4lix;
Float_t theta4x;
Float_t ene4x;
Float_t r1dealtarx;
Float_t e1alx;
Float_t th3x;
Float_t pp3xx;
Float_t ensx;
Float_t e12x;
Float_t e13x;
Float_t e23x;
Float_t q3x;
Float_t t12x;
Float_t romxx;
Float_t romyx;
Float_t ximp3x;
Float_t argox;
Float_t pal3vx;
Float_t eal3vx;
Float_t e12v3x;
Float_t e13v3x;
Float_t e23v3x;
Float_t q3v3x;
Float_t e4li7x;
};
data_t data;

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

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


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

run->Branch("p1x",&data.p1x,"p1x/F");
run->Branch("p2x",&data.p2x,"p2x/F");
run->Branch("p3x",&data.p3x,"p3x/F");
run->Branch("p4x",&data.p4x,"p4x/F");
run->Branch("de1x",&data.de1x,"de1x/F");
run->Branch("de2x",&data.de2x,"de2x/F");
run->Branch("e1x",&data.e1x,"e1x/F");
run->Branch("e2x",&data.e2x,"e2x/F");
run->Branch("e3x",&data.e3x,"e3x/F");
run->Branch("e4x",&data.e4x,"e4x/F");
run->Branch("run_nx",&data.run_nx,"run_nx/F");
run->Branch("theta1x",&data.theta1x,"theta1x/F");
run->Branch("ene1x",&data.ene1x,"ene1x/F");
run->Branch("r4delimy15x",&data.r4delimy15x,"r4delimy15x/F");
run->Branch("ar4delimy15x",&data.ar4delimy15x,"ar4delimy15x/F");
run->Branch("r4deliisox",&data.r4deliisox,"r4deliisox/F");
run->Branch("ar4deliisox",&data.ar4deliisox,"ar4deliisox/F");
run->Branch("r4delimy152x",&data.r4delimy152x,"r4delimy152x/F");
run->Branch("ar4delimy152x",&data.ar4delimy152x,"ar4delimy152x/F");
run->Branch("r4delitarx",&data.r4delitarx,"r4delitarx/F");
run->Branch("e4lix",&data.e4lix,"e4lix/F");
run->Branch("theta4x",&data.theta4x,"theta4x/F");
run->Branch("ene4x",&data.ene4x,"ene4x/F");
run->Branch("r1dealtarx",&data.r1dealtarx,"r1dealtarx/F");
run->Branch("e1alx",&data.e1alx,"e1alx/F");
run->Branch("th3x",&data.th3x,"th3x/F");   
run->Branch("pp3xx",&data.pp3xx,"pp3xx/F");
run->Branch("ensx",&data.ensx,"ensx/F");   
run->Branch("e12x",&data.e12x,"e12x/F");
run->Branch("e13x",&data.e13x,"e13x/F");  
run->Branch("e23x",&data.e23x,"e23x/F");
run->Branch("q3x",&data.q3x,"q3x/F");
run->Branch("t12x",&data.t12x,"t12x/F");
run->Branch("romxx",&data.romxx,"romxx/F");
run->Branch("romyx",&data.romyx,"romyx/F");
run->Branch("ximp3x",&data.ximp3x,"ximp3x/F");
run->Branch("argox",&data.argox,"argox/F");
run->Branch("pal3vx",&data.pal3vx,"pal3vx/F");
run->Branch("eal3vx",&data.eal3vx,"eal3vx/F");
run->Branch("e12v3x",&data.e12v3x,"e12v3x/F");
run->Branch("e13v3x",&data.e13v3x,"e13v3x/F");
run->Branch("e23v3x",&data.e23v3x,"e23v3x/F");
run->Branch("q3v3x",&data.q3v3x,"q3v3x/F");
run->Branch("e4li7x",&data.e4li7x,"e4li7x/F");


////////////////////////
//// DA CUTKIN!!!
////////////////////////

Double_t u=931.49432;             //atomic mass unit (rif. 12C)
Double_t xmproi= 10.01294*u;    //massa proiettile (10B)
Double_t xmb= 2.014102*u;       //massa bersaglio (d)
Double_t xmal= 4.002602*u;      //massa (alpha)
Double_t xmli=7.016928*u;      //massa  (7Be)
Double_t xmsp= 1.008664916*u;      //massa spettatore (neu)
Double_t enproi= 27.900;        //en.fin. di 10B@ 27,980 MeV dopo 1/2 target di CD2 da 56 ug/cm2
Double_t xmtr=1.007825*u;     // massa trasferito (p)

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

////////////////////////////////////////////
////////////////////////////////////////////
////	CALIBRAZIONE 		   /////////////
////////////////////////////////////////////
////////////////////////////////////////////

/////////////////////////////////////////
// Calibrazione psd1 
/////////////////////////////////////////
// SU PSD1 VANNO alfa
//////////////////////////////
//	PARAMETRI
//////////////////////////////

Double_t a0=216.650;
Double_t a1=-6.85526;
Double_t b0=-1.30761;
Double_t b1=0.0700185;
//Double_t th0=28.625;

Double_t m0=-0.766027;
Double_t m1=-0.00234838;
Double_t m2=4.58343E-05;
Double_t n0=0.00993842;
Double_t n1=3.48261E-06;

    Double_t altar0=1.49487e-02;
    Double_t altar1=-3.18779e-01;
    Double_t altar2=-1.23679e+02;
    Double_t altar3=7.03804e+00;
    Double_t altar4=1.34098e-01;
    Double_t altar5=8.62871e-01;
    Double_t altar6=-7.14954e-01;


//////////////////////////////
//		FUNZIONI
//////////////////////////////
TF2 *th11= new TF2("th11","(x-[0]-[1]*y)/([2]+[3]*y)",100.,4000., 100., 4000.);
// x=p3, y=e3
TF2 *ene11= new TF2("ene11","[0]+([1]*x)+([2]*x*x)+([3]+([4]*x))*y",20., 40., 0., 4000.);
// x=th3, y=e3
//AGGIUNGERE DELTAE IN ALLUMINA

    TF2 *r1dealtar = new TF2("r1dealtar","[0]*pow(x,[1])-([2]*exp(-[3]*pow(x,[4])-[5]*pow(x,[6])))/cos(y*0.0174532952)",0.,35,20,40);
//altar0=[0],altar1=[1],altar2=[2],altar3=[3],x=ene11, y=th11
TF2 *e11al = new TF2("e11al","x+y",0.,35.,0.,35);
//x=r1dealtar,      y=ene11


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

th11->FixParameter(0,a0);
th11->FixParameter(1,b0);
th11->FixParameter(2,a1);
th11->FixParameter(3,b1);

ene11->FixParameter(0,m0);
ene11->FixParameter(1,m1);
ene11->FixParameter(2,m2);
ene11->FixParameter(3,n0);
ene11->FixParameter(4,n1);

r1dealtar->FixParameter(0,altar0);
r1dealtar->FixParameter(1,altar1);
r1dealtar->FixParameter(2,altar2);
r1dealtar->FixParameter(3,altar3);
r1dealtar->FixParameter(4,altar4);
r1dealtar->FixParameter(5,altar5);
r1dealtar->FixParameter(6,altar6);

/////////////////////////////////////////
// Calibrazione psd4
/////////////////////////////////////////
// SU PSD4 VANNO liti
//////////////////////////////
//	PARAMETRI
//////////////////////////////

Double_t aa0=203.788;
Double_t aaa1=-13.1443;
Double_t bb0=-1.46033;
Double_t bb1=0.144176;
//Double_t tth0=14.45;

Double_t mm0=-0.725489;
Double_t mm1=-0.0213808;
Double_t mm2=0.000898637;
Double_t nn0=0.00752321;
Double_t nn1=0.000202458;
Double_t nn2=-6.60542e-06;

Double_t limy15_0= 0.250121;
Double_t limy15_2=-25.9538;
Double_t limy15_3=1.22409;
Double_t limy15_4=0.35809;
Double_t limy15_5=2.38052;
Double_t limy15_6=-0.218362;
Double_t limy15_1=-0.138612;
// sono = ... sono la stessa funzione in 2 step!
Double_t limy152_0= 0.250121;
Double_t limy152_2=-25.9538;
Double_t limy152_3=1.22409;
Double_t limy152_4=0.35809;
Double_t limy152_5=2.38052;
Double_t limy152_6=-0.218362;
Double_t limy152_1=-0.138612;
// 110 mb di C4H10
Double_t liiso0= 1.07514;
Double_t liiso1=-0.37057;
Double_t liiso2=-0.237324;
Double_t liiso3=-7.50018;
Double_t liiso4= 0.152912;
Double_t liiso5= 3.64552e-05;
Double_t liiso6=-2.48099;
Double_t liiso7= 4.36946;
Double_t liiso8= 0.260424;

Double_t litar0=0.0106633;
Double_t litar2=-165.425;
Double_t litar3=5.78178;
Double_t litar4=0.101071;
Double_t litar5=1.33022;
Double_t litar6=-0.545036;
Double_t litar1=-0.692923;
    
Double_t li7_0=3.20693;
Double_t li7_1=-0.297216;
Double_t li7_2= 0.00692405;


//		FUNZIONI

TF2 *th44= new TF2("th44","-0.15+(x-[0]-[1]*y)/([2]+0.99*[3]*y)",0.,4000., 0., 4000.);
// x=p4, y=e4
TF2 *ene44= new TF2("ene44","[0]+([1]*x)+([2]*x*x)+([3]+([4]*x)+[5]*x*x)*y",0., 35., 0., 4000.); 
// x=th4, y=e4
//AGGIUNGERE DELTAE IN ALLUMINA

TF1 *r4delimy15 = new TF1("r4delimy15","[0]*pow(x,[1])-([2]*exp(-[3]*pow(x,[4])-[5]*pow(x,[6])))",0.,35);
//limy15_0=[0],limy15_1=[1],limy15_2=[2],limy15_3=[3],x=ene3
TF2 *ar4delimy15 = new TF2("ar4delimy15","x+y",0.,35.,0,35);
//x=ene3              y=r3delimy15

TF1 *r4deliiso = new TF1("r4deliiso","[0]*pow(x,[1])-([2]*exp(-[3]*pow(x,[4])-[5]*pow(x,[6])-[7]*pow(x,[8])))",0.001,26);
//liiso0=[0],liiso1=[1],liiso2=[2],liiso3=[3], x=ar3delimy15
TF2 *ar4deliiso = new TF2("ar4deliiso","x+y",0.,35.,0,35);
//x=ar3delimy15       y=r3deliiso


TF1 *r4delimy152 = new TF1("r4my152","[0]*pow(x,[1])-([2]*exp(-[3]*pow(x,[4])-[5]*pow(x,[6])))",0.,35);
//limy152_0=[0],limy152_1=[1],limy152_2=[2],limy152_3=[3],x=ar3deliiso
TF2 *ar4delimy152 = new TF2("ar4delimy152","x+y",0.,35.,0,35);
//x=ar3deliiso        y=r3delimy152


TF2 *r4delitar = new TF2("r4delitar","[0]*pow(x,[1])-([2]*exp(-[3]*pow(x,[4])-[5]*pow(x,[6])))/cos(y*0.0174532952)",0.,35,10,20);
//litar0=[0],litar1=[1],litar2=[2],litar3=[3],x=ar3delimy152, y=th33
TF2 *e44li = new TF2("e44li","x+y",0.,35.,0.,35);
//x=ar3delimy152,      y=r3delitar
TF1 *e4litio = new TF1("e4litio","x+([0]+[1]*x+[2]*x*x)+0.1",0.,35.);
//x=e44li

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

th44->FixParameter(0,aa0);
th44->FixParameter(1,bb0);
th44->FixParameter(2,aaa1);
th44->FixParameter(3,bb1);

ene44->FixParameter(0,mm0);
ene44->FixParameter(1,mm1);
ene44->FixParameter(2,mm2);
ene44->FixParameter(3,nn0);
ene44->FixParameter(4,nn1);
ene44->FixParameter(5,nn2);

r4delimy15->FixParameter(0,limy15_0);
r4delimy15->FixParameter(1,limy15_1);
r4delimy15->FixParameter(2,limy15_2);
r4delimy15->FixParameter(3,limy15_3);
r4delimy15->FixParameter(4,limy15_4);
r4delimy15->FixParameter(5,limy15_5);
r4delimy15->FixParameter(6,limy15_6);

r4deliiso->FixParameter(0,liiso0);
r4deliiso->FixParameter(1,liiso1);
r4deliiso->FixParameter(2,liiso2);
r4deliiso->FixParameter(3,liiso3);
r4deliiso->FixParameter(4,liiso4);
r4deliiso->FixParameter(5,liiso5);
r4deliiso->FixParameter(6,liiso6);
r4deliiso->FixParameter(7,liiso7);
r4deliiso->FixParameter(8,liiso8);

r4delimy152->FixParameter(0,limy152_0);
r4delimy152->FixParameter(1,limy152_1);
r4delimy152->FixParameter(2,limy152_2);
r4delimy152->FixParameter(3,limy152_3);
r4delimy152->FixParameter(4,limy152_4);
r4delimy152->FixParameter(5,limy152_5);
r4delimy152->FixParameter(6,limy152_6);

r4delitar->FixParameter(0,litar0);
r4delitar->FixParameter(1,litar1);
r4delitar->FixParameter(2,litar2);
r4delitar->FixParameter(3,litar3);
r4delitar->FixParameter(4,litar4);
r4delitar->FixParameter(5,litar5);
r4delitar->FixParameter(6,litar6);
    
e4litio->FixParameter(0,li7_0);
e4litio->FixParameter(1,li7_1);
e4litio->FixParameter(2,li7_2);

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

Double_t ximp1,ximp2,cos1,cos2,sin1,sin2,cf1,cf2,sf1,sf2,cf12,c12,pp1x,pp1z,pp2x,pp2z; 
Double_t pp3x,pp3z,ximp3,pc3v3,pc3zv3,pc3xv3,pc3yv3;
Double_t pp,vp,vv,aa1,aa2,aa3,pc1x,pc1y,pc1z,pc2x,pc2y,pc2z,pc3x,pc3y,pc3z,pc1;
Double_t ec1,ec2,ec3,etotdata,etot,pc2,pc3,q3_ec;
Double_t tt,ttt,ctt,stt,v1,v2,vt,pc2xv3,pc2yv3,pc2zv3,pc2v3;
Double_t xnumx,xdena,xdenb,xnum,xden,argo,xm,xnumy,qvalorea0,qvalorea1;
Double_t th3,p3f,t12,e12,e23,e13,ens,q3,romx,romy,fi3,pp1y,pp2y,pp3y;
Double_t qvalore,aeq,beq,ceq,alfam,betam,delta,sdelta,vartest,th3v3,p3fv3;
Double_t pu20,pal3v,eal3v,ximp2v3,pp2xv3,pp2zv3,pp3zv3,pp3xv3,ximp3v3,ensv3;
Double_t e12v3,e23v3,e13v3,q3v3,ec2v3,ec3v3,etotdatav3,etotv3,q3_ecv3;

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

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


Double_t theta1=th11->Eval(p1,e1);
Double_t theta4=th44->Eval(p4,e4);
Double_t ene1=ene11->Eval(theta1,e1);
Double_t ene4=ene44->Eval(theta4,e4);
Double_t r4limy15=r4delimy15->Eval(ene4);
Double_t ar4limy15=ar4delimy15->Eval(ene4,r4limy15);
Double_t r4liiso=r4deliiso->Eval(ar4limy15);
Double_t ar4liiso=ar4deliiso->Eval(ar4limy15,r4liiso);
Double_t r4limy152=r4delimy152->Eval(ar4liiso);
Double_t ar4limy152=ar4delimy152->Eval(ar4liiso,r4limy152);
Double_t r4litar=r4delitar->Eval(ar4limy152,theta4);
Double_t e4li=e44li->Eval(ar4limy152,r4litar);
Double_t e4li7=e4litio->Eval(e4li);
Double_t r1altar=r1dealtar->Eval(ene1,theta1);
Double_t e1al=e11al->Eval(ene1,r1altar);

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


// KINEMATICS!

// sia ora 4 = li -> part "1"
// e 1 = alpha -> part "2" ricostruita
// dunque "3" sarà lo spettatore

	cos2=cos(theta1*rad);
	cos1=cos(theta4*rad);
	sin2=sin(theta1*rad);
	sin1=sin(theta4*rad);

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



        pp=sqrt(2.*xmproi*enproi);
        ximp1=sqrt(2.*e4li7*xmli);
        //        ximp2v3=sqrt(2.*eal3v*xmal);
        ximp2=sqrt(2.*e1al*xmal);
        // litio
        pp1x=ximp1*sin1;
        pp1z=ximp1*cos1;
        //alfa
        pp2x=-ximp2*sin2;
        pp2z=ximp2*cos2;
        //        pp2xv3=-ximp2v3*sin2;
        //        pp2zv3=ximp2*cos2;
        //protone spettatore
        pp3x=-pp1x-pp2x;
        pp3z=pp-pp1z-pp2z;
        pp1y=ximp1*sin1*sf1;
        pp2y=-ximp2*sin2*sf2;
        pp3y=-pp1y-pp2y;
        //        pp3xv3=-pp1x-pp2xv3;
        //        pp3zv3=pp-pp1z-pp2zv3;
        
        //	( pippo=p3x*p3x;
        //	turi=p3z*p3z;
        //	gianni=pippo+turi;
        //	roby=p3x/p3z;
        //	p3fz=msp*(-vproj+vsp)
        //  102.771=p0z=(sqrt(2*enproi*xmproi))*(xmsp/xmproi)
        //	p3fz=p3z-107.047;
        //  NO: target break up  )
        
   /*
        
        if (pp3x<=0)
        {
            th3=-ar(pp3z,pp3x);
            
            ximp3=-sqrt((pp3x*pp3x)+(pp3z*pp3z));
            
            //	th3=-atan(roby);
            //	p3f=-sqrt((p3fz*p3fz)+(p3x*p3x));
            p3f=-sqrt((pp3z*pp3z)+(pp3x*pp3x));
        }
        // ximp3=p3f
        else
        {
            //	th3=atan(roby);
            th3=ar(pp3z,pp3x);
            p3f=sqrt((pp3z*pp3z)+(pp3x*pp3x));
            ximp3=sqrt((pp3x*pp3x)+(pp3z*pp3z));
        }
    
        ens=(ximp3*ximp3)/(2.*xmsp);
     */

//----------------------------------------------------------------------------
//	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/(xmproi+xmb);
// vv è quella che nell' ohlsen si chiama V, velocità del centro di massa
        aa2=vv*sqrt(xmal/2.);
        aa1=vv*sqrt(xmli/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.*xmli);
        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.*xmal);
        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.*xmli);
        ec2=pc2*pc2/(2.*xmal);
        ec3=pc3*pc3/(2.*xmsp);
        
        fi3=deg*ar(pp3y,pp3x);
        th3=deg*ar(pc3x,(pc3z+xmsp*vv));
        
        if (fi3==0) {
             ximp3=sqrt((pp3x*pp3x)+(pp3z*pp3z));
        }
        else
        {
             ximp3=-sqrt((pp3x*pp3x)+(pp3z*pp3z));
        }
        
        ens=(ximp3*ximp3)/(2.*xmsp);


        etotdata=e4li7+e1al+ens;               // energia totale a disp. nel c.m. 1,2,3
        etot=ec1+ec2+ec3;
	
        q3_ec=ec1+ec2+ec3-enproi;
        q3=e4li7+e1al+ens-enproi;

    

/////////////////////////////////////////////////////////////////////
// Calcolo a tre variabili. Escludo energia 4He su PSD4      /////////
//
        qvalorea0=-1.07877;                                        //
        qvalorea1=-1.50787;                                         //
        
        aeq=0.5*((xmsp+xmal)/xmal);                             //
        beq=sqrt(2.*xmli*e4li7)*c12-sqrt(2.*xmproi*enproi)*cos2; //
        alfam=sqrt(xmli*e4li7*xmproi*enproi);                    //
        if (q3<-0.9&&q3>-1.2) {
            qvalore=qvalorea0;
        }
        else if (q3<-1.4&&q3>-1.6) {
            qvalore=qvalorea1;
        }
        else
        {
            qvalore=0;
        }
        betam=xmsp*(e4li7-enproi-qvalore);                       //
        ceq=xmproi*enproi+xmli*e4li7-2.*alfam*cos1+betam;        //
        delta=beq*beq-4.*aeq*ceq;                               //
        sdelta=sqrt(fabs(delta));                               //
        pu20=sqrt(2.*xmal*e1al);                                //
        vartest=2*aeq*pu20+beq;                                 //
        
        if(vartest>=0.)                                         //
        {                                                       //
            pal3v=(-beq+sdelta)/(2.*aeq);                       //
            eal3v=(pal3v*pal3v)/(2.*xmal);                    //
        }                                                       //
        else                                                //
        {                                                       //
            pal3v=(-beq-sdelta)/(2.*aeq);                       //
            eal3v=(pal3v*pal3v)/(2.*xmal);                      //
        }                                                       //
////////////////////////////////////////////////////////////////////////
        //3v
        
        ximp2v3=sqrt(2.*eal3v*xmal);
        pp2xv3=-ximp2v3*sin2;
        pp2zv3=ximp2v3*cos2;
        pp3xv3=-pp1x-pp2xv3;
        pp3zv3=pp-pp1z-pp2zv3;
        
        if (pp3xv3<=0)
        {
            th3v3=-ar(pp3zv3,pp3xv3);
            ximp3v3=-sqrt((pp3xv3*pp3xv3)+(pp3zv3*pp3zv3));
            p3fv3=-sqrt((pp3zv3*pp3zv3)+(pp3xv3*pp3xv3));
        }
        else
        {
            th3v3=ar(pp3zv3,pp3xv3);
            p3fv3=sqrt((pp3zv3*pp3zv3)+(pp3xv3*pp3xv3));
            ximp3v3=-sqrt((pp3xv3*pp3xv3)+(pp3zv3*pp3zv3));
        }
        ensv3=(ximp3v3*ximp3v3)/(2.*xmsp);
        
        pc2xv3=ximp2v3*sin2*cf2;
        pc2yv3=ximp2v3*sin2*sf2;
        pc2zv3=ximp2v3*cos2-aa2*sqrt(2.*xmal);
        pc2v3=sqrt(pc2xv3*pc2xv3+pc2yv3*pc2yv3+pc2zv3*pc2zv3);
        
        pc3xv3=-pc1x-pc2xv3;
        pc3yv3=-pc1y-pc2yv3;
        pc3zv3=-pc1z-pc2zv3;
        pc3v3=sqrt(pc3xv3*pc3xv3+pc3yv3*pc3yv3+pc3zv3*pc3zv3);
        
        
        
        ec2v3=pc2v3*pc2v3/(2.*xmal);
        ec3v3=pc3v3*pc3v3/(2.*xmsp);
        etotdatav3=e4li7+eal3v+ensv3;          // energia totale a disp. nel c.m. 1,2,3
        etotv3=ec1+ec2v3+ec3v3;
        
        q3_ecv3=ec1+ec2v3+ec3v3-enproi;
        q3v3=e4li7+eal3v+ensv3-enproi;

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

//ohlsen, calcoli in fogli miei: le formule sotto sono giuste

	xm=xmli+xmal+xmsp;

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

// 3v
    e12v3=etotv3-xm*ec3v3/(xmli+xmal);
    e13v3=etotv3-xm*ec2v3/(xmli+xmsp);
	e23v3=etotv3-xm*ec1/(xmal+xmsp);
////////////////////////////////////////////////////////////

// Calcolo theta_cm 4 variabili		

       Double_t pi=3.1415; 


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

/*      b-u proj
        vp=sqrt(2.*enproi/xmproi);
        v1=sqrt(2.*e3li7/xmli);
        v2=sqrt(2.*e4al/xmal);

	vtzc= ((xmsp+xmtr)*vp/xmtr)-pp3z/xmtr;
	vtx=-pp3x/xmtr;
		
// costruzione dell'argomento del theta_cm (slaus)

        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

       vp=sqrt(2.*enproi/(xmproi));
       v1=sqrt(2.*e4li7/xmli);
       v2=sqrt(2.*e1al/xmal);
       vt=fabs(ximp3)/(xmtr);
       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=enproi-e4li7-e1al;

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


//g->cd();

data.p1x=p1;
data.p2x=p2;
data.p3x=p3;
data.p4x=p4;
data.de1x=de1;
data.de2x=de2;
data.e1x=e1;
data.e2x=e2;
data.e3x=e3;
data.e4x=e4;
data.run_nx=run_n;
data.theta1x=theta1;
data.ene1x=ene1;
data.r4delimy15x=r4limy15;
data.ar4delimy15x=ar4limy15;
data.r4deliisox=r4liiso;
data.ar4deliisox=ar4liiso;
data.r4delimy152x=r4limy152;
data.ar4delimy152x=ar4limy152;
data.r4delitarx=r4litar;
data.e4lix=e4li;
data.theta4x=theta4;
data.ene4x=ene4;
data.r1dealtarx=r1altar;
data.e1alx=e1al;
data.th3x=th3;
data.pp3xx=pp3x;
data.ensx=ens;
data.e12x=e12;
data.e13x=e13;
data.e23x=e23;
data.q3x=q3;
data.t12x=t12;
data.romxx=romx;
data.romyx=romy;
data.ximp3x=ximp3;
data.argox=argo;
data.pal3vx=pal3v;
data.eal3vx=eal3v;
data.e12v3x=e12v3;
data.e13v3x=e13v3;
data.e23v3x=e23v3;
data.q3v3x=q3v3;
data.e4li7x=e4li7;

run->Fill();

}


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

exit(0);

}