#define E03104_cxx
#include "run.h"
#include "bin.h"
#include "const.h"
#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <TDirectory.h>
#include <TEventList.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TCut.h>
#include <TROOT.h>
#include "TObjArray.h"
#include <TLorentzVector.h>
#include <TBranch.h>
#include <TLeaf.h>
using std::string;
using namespace std;
Run::Run(string filename) {
rootfile = filename;
RunNumber = 0;
run_parameter_ok = false;
beam_parameter_ok = false;
calc_ok = true;
// check for existence of file
ifstream inp;
inp.open(rootfile.c_str(), ifstream::in);
inp.close();
if (inp.fail()) {
cout << "ERROR: can't open rootfile: `" << rootfile << "'" << endl;
exit (1);
}
int iroot = filename.find(".root");
if (iroot > 4)
RunNumber = atoi((filename.substr(iroot-4,4)).c_str());
if (RunNumber == 0) {
cout << "Error: Can't find run number for file: " << filename << endl;
exit(1);
}
}
Run::~Run() {
}
ostream& operator<< (ostream &os, const Run* &obj) {
os << "Root Files:" << endl
<< "----------" << endl << endl;
if (obj == NULL) {return os;};
Run* list_of_runs = (Run*) obj;
Run* this_run = list_of_runs;
int i = 1;
while (this_run != NULL) {
os << setw(4) << i
<< setw(6) << this_run->kinematics
<< " " << this_run->rootfile << endl;
this_run = this_run->get_next_run();
i++;
};
this_run = list_of_runs;
os << endl;
os << "Run E0 h th_p p0_p th_e p0_e C_x effi chi chi chi " << endl;
os << "-----------------------------------------------------------------------------------" << endl;
while (this_run != NULL) {
os <<
setw(4) << this_run->RunNumber <<
fixed <<
setw(7) << setprecision(3) << this_run->e0 <<
setw(7) << setprecision(3) << this_run->h <<
setw(6) << setprecision(1) << this_run->thLeft / DEGTORAD <<
setw(8) << setprecision(4) << this_run->pLeft <<
setw(6) << setprecision(1) << this_run->thRight / DEGTORAD<<
setw(8) << setprecision(4) << this_run->pRight <<
setw(6) << setprecision(1) << this_run->xfpp <<
setw(7) << setprecision(3) <<
(this_run->nEventsBeforeCuts > 0 ?
(double) this_run->nEventsAfterCuts / (double) this_run->nEventsBeforeCuts : 0) <<
setw(8) << setprecision(2) << this_run->stat_chi.get_min() / DEGTORAD <<
setw(8) << setprecision(2) << this_run->stat_chi.get_max() / DEGTORAD <<
setw(8) << setprecision(2) << this_run->stat_chi.get_mean() / DEGTORAD <<
setw(8) << setprecision(3) << this_run->stat_Q2.get_min() <<
setw(8) << setprecision(3) << this_run->stat_Q2.get_max() <<
setw(8) << setprecision(3) << this_run->stat_Q2.get_mean() <<
endl;
this_run = this_run->get_next_run();
};
os << endl
<< "Counter:" << endl
<< "-------" << endl << endl;
Long_t total = 0;
Long_t after = 0;
Long_t heli_plus = 0;
Long_t heli_minus = 0;
this_run = list_of_runs;
while (this_run != NULL) {
os <<
setw(4) << this_run->RunNumber <<
setw(10) << this_run->nEventsBeforeCuts <<
setw(9) << this_run->nEventsAfterCuts <<
setw(10) << this_run->helicity_plus <<
setw(10) << this_run->helicity_minus <<
setw(10) << fixed << setprecision(5) <<
(this_run->nEventsAfterCuts > 0 ?
(double) (this_run->helicity_plus-this_run->helicity_minus) /
(double) (this_run->helicity_plus+this_run->helicity_minus) : 0) <<
endl;
total += this_run->nEventsBeforeCuts;
after += this_run->nEventsAfterCuts;
heli_plus += this_run->helicity_plus;
heli_minus += this_run->helicity_minus;
this_run = this_run->get_next_run();
};
os <<
"sum:" <<
setw(10) << total <<
setw(9) << after <<
setw(10) << heli_plus <<
setw(10) << heli_minus <<
setw(10) << fixed << setprecision(5) <<
(after > 0 ?
(double) (heli_plus-heli_minus) /
(double) (heli_plus+heli_minus) : 0) <<
endl;
return os;
}
// ----------------------------------------------------------------------------
void Run::open(TCut GlobalCut) {
// Phil: temp stuff
/* TFile tempf("/u/home/pcarter/palmetto/fpp_gep_new_2818.root", "READ");
if (tempf.IsZombie()) {
cout << "Error opening tempf" << endl;
exit (-1);
}
TTree* temptree = (TTree *)tempf.Get("h2");
float tmpvar1, tmpvar2;
temptree->SetBranchAddress("Helicite",&tmpvar1);
temptree->SetBranchAddress("Thtra",&tmpvar2);
temptree->GetEvent(10);
cout << "Helicite: " << tmpvar1 << "\nThtra: " << tmpvar2 << "\n";
tempf.Close();
exit(0);
*/
// load cosy file
if (spin_transport_model == COSY)
spin_matrix.load_cosy_file(spin_transport_file);
// open root file
RunFile = new TFile(rootfile.c_str(), "READ");
if (RunFile->IsZombie()) {
cout << "Error opening file" << endl;
exit (-1);
}
// open a tree named T
// if it doesn't exist, open a tree named h2
// if no tree by either name exists, we're sunk!
this_tree = (TTree *)RunFile->Get("T");
if (!this_tree) this_tree = (TTree *)RunFile->Get("h2");
if (!this_tree) {
cout << "Error: Can't open the Root file tree. Tried names T and h2.\n";
exit(1);
}
nEventsBeforeCuts = (Long_t) this_tree->GetEntries();
//this_tree->Draw(">>eList", GlobalCut, 0, 7000);
this_tree->Draw(">>eList", GlobalCut);
evtList = (TEventList*)gDirectory->Get("eList");
//evtList->Print("all");
nEventsAfterCuts = evtList->GetN();
cout << "Reading file: " << rootfile << endl
<< " with " << nEventsBeforeCuts << " events before cuts"
<< " and " << nEventsAfterCuts << " events after cuts" << endl;
this_tree->SetEventList(evtList);
event_data = new E03104(this_tree);
event_count = 0;
helicity_plus = 0;
helicity_minus = 0;
}
// ----------------------------------------------------------------------------
void Run::close() {
RunFile->Close(); // closes file and delete the tree object
delete event_data;
}
// ----------------------------------------------------------------------------
bool Run::is_next_event() {
if (event_count < nEventsAfterCuts) {
GetEvent(evtList->GetEntry(event_count));
// cout << "#: " << evtList->GetEntry(event_count) << " " << p_p << endl;
// cout << "#: " << event_count << " " << evtList->GetEntry(event_count) << " " << p_p << endl;
event_count++;
// if (evtList->GetEntry(event_count) > 6800) {exit(1);}
if (event_count % 50000 == 0) cout << event_count << endl;
return true;
}
return false;
}
// ----------------------------------------------------------------------------
void Run::calc_kinematics(ReactionType reaction) {
// note: this function doesn't handle the GAMMAD reaction type
double ef = 0;
double omega = 0;
double gamma = sqrt(p_p*p_p + m_p*m_p)/m_p;
precession_angle = (BendAngle * DEGTORAD + Th_tgt_p - Th_tra_p) * kappa * gamma;
// cout << "chi = " << (precession_angle / DEGTORAD) << endl;
prec = 0;
prmag = 0;
theta_pq = 0;
phi_x = 0;
kin_Ebeam = 0;
kin_Ee = 0;
kin_theta_e = 0;
kin_tau = 0;
kin_epsilon = 0;
px_prime = 0.;
pz_prime = 0.;
if (reaction == EEP || reaction == EP_COINC) {
// k_tgt : incident electron beam direction
// q_tgt : proton momentum (= three momentum-transfer vector) direction
Double_t Psi_tgt_e = atan(tan(Th_tgt_e) * cos(Ph_tgt_e));
Double_t cos_phi_e = cos(Ph_tgt_e + thRight);
Double_t sin_phi_e = sin(Ph_tgt_e + thRight);
Double_t cos_psi_e = cos(-Psi_tgt_e);
Double_t sin_psi_e = sin(-Psi_tgt_e);
Double_t Psi_tgt_p = atan(tan(Th_tgt_p) * cos(Ph_tgt_p));
Double_t cos_phi_p = cos(Ph_tgt_p + thLeft);
Double_t sin_phi_p = sin(Ph_tgt_p + thLeft);
Double_t cos_psi_p = cos(-Psi_tgt_p);
Double_t sin_psi_p = sin(-Psi_tgt_p);
Double_t m_target = 0;
if (reaction == EEP) m_target = m_4he;
else m_target = m_p;
TVector3 v3_ki(0, 0, e0);
TVector3 v3_kf(p_e*sin_psi_e, p_e*cos_psi_e * sin_phi_e, p_e*cos_psi_e * cos_phi_e);
TVector3 v3_pf(p_p*sin_psi_p, p_p*cos_psi_p * sin_phi_p, p_p*cos_psi_p * cos_phi_p);
TVector3 v3_q = v3_ki - v3_kf;
TVector3 y = (v3_ki.Cross(v3_kf)).Unit();
TVector3 z = v3_q.Unit();
TVector3 x = y.Cross(z);
// Get recoil momentum
TVector3 v3_pr = v3_q - v3_pf;
prmag = v3_pr.Mag();
if (v3_q.Dot(v3_pr) > 0) prec = -prmag;
else prec = +prmag;
// Get angles between proton and q
theta_pq = v3_q.Angle(v3_pf) / DEGTORAD;
Double_t argnum = v3_pf.Dot(y);
Double_t argden = v3_pf.Dot(x);
phi_x = 0;
if (argden != 0.) phi_x = atan(argnum/argden);
if (argden < 0) phi_x += PI;
else if (argnum < 0 && argden > 0) phi_x += 2 * PI;
else if (argnum < 0 && argden == 0) phi_x = 1.5 * PI;
else if (argnum > 0 && argden == 0) phi_x = 0.5 * PI;
else if (argnum == 0 && argden == 0) phi_x = 0.0;
phi_x /= DEGTORAD;
TLorentzVector target(0, 0, 0, m_target);
TLorentzVector proton(v3_pf, sqrt(m_p * m_p + p_p * p_p));
TLorentzVector ki(v3_ki, e0);
TLorentzVector kf(v3_kf, p_e);
Q2 = -(ki - kf).Mag2(); // -1 * magnitude squared of (ki - kf)
ef = kf.E(); // == p_e
omega = e0 - ef;
} else if (reaction == EP_SINGLE) {
Double_t Psi_tgt_p = atan(tan(Th_tgt_p) * cos(Ph_tgt_p));
Double_t cos_phi_p = cos(Ph_tgt_p + thLeft);
Double_t sin_phi_p = sin(Ph_tgt_p + thLeft);
Double_t cos_psi_p = cos(-Psi_tgt_p);
Double_t sin_psi_p = sin(-Psi_tgt_p);
TLorentzVector target(0, 0, 0, m_p);
TLorentzVector proton(p_p * sin_psi_p,
p_p * cos_psi_p * sin_phi_p,
p_p * cos_psi_p * cos_phi_p,
sqrt(m_p * m_p + p_p * p_p));
Q2 = -(proton - target).Mag2();
omega = Q2 / (2 * m_p);
ef = e0 - omega;
}
if (reaction == EP_COINC || reaction == EP_SINGLE) {
if (Q2 / ef / e0 > 0.) {
kin_Ebeam = e0;
kin_Ee = ef;
kin_theta_e = 2.* asin(sqrt(Q2 / 4. / e0 / ef));
kin_tau = Q2 / 4 / (m_p*m_p);
kin_epsilon = 1./(1 + 2.*(1+kin_tau)*pow(kin_theta_e/2,2));
// Parameterization of the proton form factor ratio from
// O. Gayou, et al., Phys. Rev. C 64, 038202 (2001), Eq. (2)
// Double_t ratio = (1. - 0.14 * (Q2 - 0.3))/mu_p;
// Parameterization of the proton form factor ratio from
// O. Gayou, et al., Phys. Rev. Lett. 88, 092301 (2002), Eq. (4)
Double_t ratio = (1. - 0.13 * (Q2 - 0.04))/mu_p;
px_prime = -2*sqrt(kin_tau*(1+kin_tau))*ratio/(ratio*ratio + kin_tau/kin_epsilon)
* tan(kin_theta_e/2);
pz_prime = (e0+ef)/m_p * sqrt(kin_tau*(1+kin_tau)) /
(ratio*ratio+kin_tau/kin_epsilon) * pow(tan(kin_theta_e/2),2);
// cout << " ========================================== " << endl;
// cout << " Ratio = " << ratio*2.79 << endl;
// cout << " Px_prime = " << px_prime << endl;
// cout << " Pz_prime = " << pz_prime << endl;
// cout << " ========================================== " << endl;
} else calc_ok = false;
}
stat_Q2.add(Q2);
stat_chi.add(precession_angle);
}
// ----------------------------------------------------------------------------
void Run::calc_spin_transport(ReactionType reaction) {
// note: this function doesn't handle the GAMMAD reaction type
spin_matrix.set_to_unit_matrix();
switch (reaction) {
case EP_SINGLE:
spin_matrix.lab_to_tra_EP(thLeft, Th_tgt_p, Ph_tgt_p);
break;
case EP_COINC:
spin_matrix.lab_to_tra_EP(thLeft, Th_tgt_p, Ph_tgt_p);
// spin_matrix.lab_to_tra_EEP(e0, p_e, thRight, thLeft, Th_tgt_e, Ph_tgt_e);
break;
case EEP:
spin_matrix.lab_to_tra_EEP(e0, p_e, thRight, thLeft, Th_tgt_e, Ph_tgt_e);
break;
default:
cout << "Error: unknown reaction: " << reaction << endl;
exit (1);
}
// cout << "real spin_cosy(" << pLeft << ", " << y_tgt_p << ", " << Th_tgt_p << ", " << Ph_tgt_p << ", " << p_p << ");\n";
switch (spin_transport_model) {
case DIPOLE:
spin_matrix.spin_dipole(precession_angle);
break;
case PENTCHEV:
spin_matrix.spin_pentchev(p_p, y_tgt_p, Th_tgt_p, Ph_tgt_p,
Th_tra_p, Ph_tra_p, 0);
break;
case PENTCHEV1:
spin_matrix.spin_pentchev(p_p, y_tgt_p, Th_tgt_p, Ph_tgt_p,
Th_tra_p, Ph_tra_p, 1);
break;
case PENTCHEV2:
spin_matrix.spin_pentchev(p_p, y_tgt_p, Th_tgt_p, Ph_tgt_p,
Th_tra_p, Ph_tra_p, 2);
break;
case COSY:
spin_matrix.spin_cosy(pLeft, y_tgt_p, Th_tgt_p, Ph_tgt_p, p_p);
break;
default:
cout << "ERROR: unknown spin transport model" << endl;
exit(1);
}
spin_matrix.tra_to_fpp(Th_tra_p, Ph_tra_p);
}
// ----------------------------------------------------------------------------
void Run::calc_weights() {
calc_weights(xfpp, h);
}
// ----------------------------------------------------------------------------
Run* Run::get_next_run() {
return next_run;
}
// ----------------------------------------------------------------------------
void Run::set_next_run(Run* run) {
next_run = run;
}
// ----------------------------------------------------------------------------
Double_t* Run::get_address(string name) {
if (name.compare("delta_p") == 0) return &Delta_p;
if (name.compare("y_tgt_p") == 0) return &y_tgt_p;
if (name.compare("Th_tgt_p") == 0) return &Th_tgt_p;
if (name.compare("Ph_tgt_p") == 0) return &Ph_tgt_p;
if (name.compare("Th_tra_p") == 0) return &Th_tra_p;
if (name.compare("Ph_tra_p") == 0) return &Ph_tra_p;
if (name.compare("Q2") == 0) return &Q2;
if (name.compare("prmag") == 0) return &prmag;
if (name.compare("prec") == 0) return ≺
if (name.compare("theta_pq") == 0) return &theta_pq;
if (name.compare("phi_x") == 0) return &phi_x;
if (name.compare("chi") == 0) return &precession_angle;
if (name.compare("ph_fpp") == 0) return &Ph_fpp;
if (name.compare("th_fpp") == 0) return &Th_fpp;
cout << "Error: get_address " << name << endl;
exit(1);
}
// ----------------------------------------------------------------------------
void Run::GetEvent(int i) {
// read the ith Event of the present run
// this puts data in the variables that were designated using SetBranchAddress
this_tree->GetEvent(i);
// change of variables and units
if (event_data->newstyle) {
Ph_tgt_e = atan(event_data->Ph_tgte/1000.0);
Th_tgt_e = atan(event_data->Th_tgte/1000.0);
Ph_tgt_p = atan(event_data->Ph_tgt/1000.0);
Th_tgt_p = atan(event_data->Th_tgt/1000.0);
Ph_tra_p = atan(event_data->Phtra/1000.0);
Th_tra_p = atan(event_data->Thtra/1000.0);
// round Helicity to the nearest integer ("unknown helicity" is something like 1e-15)
Helicity = -1.0*double(int(event_data->Helicite));
if(Helicity*event_data->Helicite != -1.0){
cout << "Helicity mismatch!!!!" << endl;
cout << "Helicity = " << Helicity << " ... Helicite = " << event_data->Helicite << endl;
}
Th_fpp = event_data->Th_fpp;
Ph_fpp = event_data->Ph_fpp;
// Ph_fpp = 90. - event_data->Ph_fpp;
// if (Ph_fpp < 0.) Ph_fpp += 360.;
p_p = event_data->Pp/1000.0;
// Phil: this next line works only for a hydrogen target (p_i = 0)
// TODO: fix that
p_e = e0 + m_p - sqrt(p_p*p_p + m_p*m_p) - event_data->E_miss/1000.0;
Delta_p = event_data->Dppc; // this variable is not used
y_tgt_p = event_data->Y_tgt/1000.0;
} else {
Ph_tgt_e = atan(event_data->R_tr_tg_ph[0]);
Th_tgt_e = atan(event_data->R_tr_tg_th[0]);
Ph_tgt_p = atan(event_data->L_tr_tg_ph[0]);
Th_tgt_p = atan(event_data->L_tr_tg_th[0]);
Ph_tra_p = atan(event_data->L_tr_ph[0]);
Th_tra_p = atan(event_data->L_tr_th[0]);
Helicity = event_data->Beam_HL_helicity;
Th_fpp = event_data->L_fpp_th_az;
Ph_fpp = 90. - event_data->L_fpp_ph_az;
if (Ph_fpp < 0.) Ph_fpp += 360.;
p_e = event_data->R_gold_p;
p_p = event_data->L_gold_p;
Delta_p = event_data->L_tr_tg_dp[0]; // this variable is not used
y_tgt_p = event_data->L_tr_tg_y[0];
}
}
// ----------------------------------------------------------------------------
void Run::calc_weights(Double_t xfpp, Double_t pol_e) {
Ac = Ac_McNaughton(p_p, Th_fpp * DEGTORAD, xfpp);
// Phil: cos1 and sin1 were formerly floats
double cos1 = cos(Ph_fpp * DEGTORAD);
double sin1 = sin(Ph_fpp * DEGTORAD);
// cout << "=========== spin ============" << endl;
// cout << spin_matrix(1,0) << " " << spin_matrix(1,1) << " " << spin_matrix(1,2) << endl;
// cout << spin_matrix(0,0) << " " << spin_matrix(0,1) << " " << spin_matrix(0,2) << endl;
// cout << "=========== spin end============" << endl;
lambda[0] = spin_matrix(1,0) * cos1 + spin_matrix(0,0) * sin1;
lambda[1] = spin_matrix(1,1) * cos1 + spin_matrix(0,1) * sin1;
lambda[2] = spin_matrix(1,2) * cos1 + spin_matrix(0,2) * sin1;
lambda[3] = pol_e * Helicity * lambda[0];
lambda[4] = pol_e * Helicity * lambda[1];
lambda[5] = pol_e * Helicity * lambda[2];
if (Helicity > 0) helicity_plus++;
else if (Helicity < 0) helicity_minus++;
}
// ----------------------------------------------------------------------------
// Calculate energy loss of the proton in carbon
double Run::CalcEnergyLoss(double spec_p /* In GeV */, double DThickness /* in inch */) {
if (DThickness==0) return spec_p;
double rho=1.7; // carbon density
double d=DThickness; // door thickness
double Z=6.; //For carbon
double A=12.011; // For carbon
double pp = spec_p; // Everything below in units of GeV.
double tp = sqrt(pp*pp + m_p*m_p) - m_p;
double Tp1, Tp2;
// add calculation of Energy loss
// First from Ed's MC he determined A function for
// the energy loss in the material before the carbon
double eloss = (59.052-0.1103*pp*1000+6.223e+1*pp*pp)/1000.;
double Tp0 = tp - eloss;
pp = sqrt( (tp+m_p)*(tp+m_p) - m_p*m_p);
// energy loss in carbon
for ( Int_t j=1;j<=2;j++) {
if(j == 2) tp=Tp1;
double bet2=pow((pp/(m_p+tp)),2);
double gam2=1./(1.-bet2);
double gam=sqrt(gam2);
double dedx=0.0003071*(Z/(A*bet2));
dedx=dedx*(log(2.*511000*gam2*bet2/78.)-bet2-log(gam));
if (j == 1) Tp1=Tp0-(d*rho*dedx/2.);
else if (j == 2) Tp2=Tp0-(d*rho*dedx/2.);
// cout << " de/dx = " << dedx << endl;
}
tp = (Tp1+Tp2)/2;
pp = sqrt( (tp+m_p)*(tp+m_p) - m_p*m_p);
// cout << " pcent = " << pp << " Tp = " << tp << endl ;
return pp; // Return the energy of the proton in the center of the carbon
}
// ----------------------------------------------------------------------------
double Run::Ac_McNaughton(double fPP, double fThFPP, double DoorThickness) {
// fPP (GeV/c)
// fThFPP (rad)
// DoorThickness (inch)
/*
Use the N-D/McNaughton parametrization
*/
// Proton mass in GeV2
// const double M = 0.938271998;
// High energy coefficients
const double a0 = 0.974;
const double a1 = -2.298;
const double a2 = 43.21;
const double a3 = 19.44;
const double b0 = -26.49;
const double b1 = -30.29;
const double b2 = 1048.7;
const double b3 = 2826.;
const double c0 = 750.;
const double c1 = -833.;
const double c2 = -861.;
const double c3 = 7205.;
const double d0 = 0.1402;
const double d1 = 0.0414;
const double d2 = -0.6052;
const double d3 = 0.42;
const double p17 = -1.2831;
// Low energy coefficients
const double a0l = 4.576;
const double a1l = -0.78;
const double a2l = 14.6;
const double a3l = 98.7;
const double b0l = -11.9;
const double b1l = -131.7;
const double b2l = 1682.;
const double b3l = -5574.;
const double c0l = 673.;
const double c1l = 1866.;
const double c2l = 7290.;
const double c3l = -38800.;
const double a4l = -532.;
const double b4l = 6946.;
const double c4l = -55800.;
const double p17l = -0.686;
// high energy coefficients mkj NIM A241 (1985) 435
const double ahipar[4]={1.6575,1.3855,9.7700,-149.27};
const double bhipar[4]={-16.346,152.53,139.16,-3231.1};
const double chipar[4]={1052.2,-3210.8,-2293.3,60327.0};
const double dhipar[4]={0.13887,-0.19266,-0.45643,8.1528};
// Low energy coefficients mkj NIM A241 (1985) 435
const double apar[5]={5.3346,-5.5361,2.8353,61.915,-145.54};
const double bpar[5]={-12.774,-68.339,1333.5,-3713.5,3738.3};
const double cpar[5]={1095.3,949.5,-28012.0,96833.0,-118830.0};
// Low energy coeff April -Gibroni nim 215 (1983) 147
//const double alpha[5]={3.6991,0.26957,-0.0012157,0.17072};
//const double beta[5]={-8.7225,-3.7161,12.869,-2.6088,1.6024};
//const double gamma[5]={351.97,271.44,-113.71,-10.407,20.331};
//const double capC=75.383;
//const double C[2]={0.12,0.18472};
double pp = CalcEnergyLoss(fPP, DoorThickness);
double t = fThFPP;
double r = pp*sin(t);
double tp = sqrt(pp*pp + m_p*m_p) - m_p;
double fAnalyzingPower = 1;
if (tp >= 0.450) {
double p = pp + p17;
double psq = p*p;
double psqp = psq*p;
// double psqsq = psq*psq;
double a = a0 + a1*p + a2*psq + a3*psqp;
double b = b0 + b1*p + b2*psq + b3*psqp;
double c = c0 + c1*p + c2*psq + c3*psqp;
double d = d0 + d1*p + d2*psq + d3*psqp;
fAnalyzingPower = d*pp*sin(5*t) + a*r / (1 + b*r*r + c*r*r*r*r);
p = pp -1.2;
a = 0; b = 0; c = 0; d = 0;
for (Int_t j=0;j<4;j++) {
a = a + ahipar[j]*TMath::Power(p,j);
b = b + bhipar[j]*TMath::Power(p,j);
c = c + chipar[j]*TMath::Power(p,j);
d = d + dhipar[j]*TMath::Power(p,j);
}
fAnalyzingPower=d*pp*sin(5*t)+a*r/(1.+b*TMath::Power(r,2)+c*TMath::Power(r,4));
} else {
double p = pp + p17l;
double psq = p*p;
double psqp = psq*p;
double psqsq = psq*psq;
double a = a0l + a1l*p + a2l*psq + a3l*psqp + a4l*psqsq;
double b = b0l + b1l*p + b2l*psq + b3l*psqp + b4l*psqsq;
double c = c0l + c1l*p + c2l*psq + c3l*psqp + c4l*psqsq;
fAnalyzingPower = a*r / (1 + b*r*r + c*r*r*r*r);
p = pp -0.7;
a = 0; b = 0; c = 0;
for (Int_t j=0;j<5;j++) {
a = a + apar[j]*TMath::Power(p,j);
b = b + bpar[j]*TMath::Power(p,j);
c = c + cpar[j]*TMath::Power(p,j);
}
fAnalyzingPower=a*r/(1.+b*TMath::Power(r,2)+c*TMath::Power(r,4));
}
// fAnalyzingPower=1;
if (fAnalyzingPower > -10) {
if (fAnalyzingPower < 10) {
return fAnalyzingPower;
}
}
// todo: this should never happen!
cout << "ERROR: Ac = " << fAnalyzingPower << endl;
fAnalyzingPower = 0;
return fAnalyzingPower;
}