DEMPgen/src/Asymmetry.cxx

0001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0002   Asymmetry class implementation
0003   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
0004
0005 #include "Asymmetry.hxx"
0006
0007 #include <string>
0008 #include <algorithm>
0009 #include <vector>
0010 #include <array>
0011 #include <stdlib.h>
0012 #include <stdio.h>
0013
0014 #include "TF1.h"
0015 #include "TFile.h"
0016 #include "TTree.h"
0017 #include "TH1.h"
0018 #include "TH2.h"
0019 #include "TMath.h"
0020 #include "TGraph.h"
0021
0022 using namespace std;
0023 using namespace TMath;
0024
0025 extern TFile * AsymmFile;
0026 extern char* DEMPgen_Path;
0027 TTree * GK_Raw;
0028
0029 Asymmetry::Asymmetry(char * in_AsyName, char * in_Func,
0030                      vector<double> in_Qsq, bool refit)
0031 {
0032   AsyNameStr = in_AsyName;
0033   FuncForm = in_Func;
0034   if (refit) Parameterize(in_Qsq);
0035   else SetPars(in_Qsq);
0036 }
0037
0038 int Asymmetry::Parameterize(vector<double> in_Qsq)
0039 {
0040
0041   //Go to default work file if not extern not available
0042   if (AsymmFile->IsZombie()){
0043     AsymmFile = new TFile(Form("%s/data/input/Asymmetries.root", DEMPgen_Path));
0044   }
0045
0046   GK_Raw = (TTree*)AsymmFile->Get("GK_Raw");
0047
0048   nQsq = in_Qsq.size();
0049   if (nQsq == 0) {
0050     return Parameterize();
0051   }
0052
0053   else{
0054     Qsq_Vec = in_Qsq;
0055   }
0056
0057   char tfnamestr[100] = "%s_%d";
0058   char plotstr[100] = "tp:%s";
0059   char cutstr[100] = "q2==%g";
0060
0061   int i = 0;
0062   int j = 0;
0063   char intstr[2];
0064   char tempname1[100];
0065   char tempname2[100];
0066   char tempname3[100];
0067   TGraph *gtemp;
0068
0069   int n;
0070
0071   double param2;
0072
0073
0074   for (i = 0; i < nQsq; i++){
0075
0076     sprintf(tempname1, tfnamestr, AsyNameStr, i);
0077     AsyFunction.push_back(new TF1(tempname1, FuncForm));
0078     sprintf(tempname2, plotstr, AsyNameStr);
0079     sprintf(tempname3, cutstr, Qsq_Vec[i]);
0080
0081     n = GK_Raw->Draw(tempname2, tempname3, "goff");
0082
0083     // The previous version of this parameterization used
0084     // only one TF1. As a result, the parameters were
0085     // initialized by the previous fit.
0086     // The following code attempts to replicate that.
0087     if (i>0){
0088       nPars = AsyFunction.at(i)->GetNpar();
0089       for (j=0; j<nPars; j++){
0090         AsyFunction.at(i)->SetParameter(j, AsyFunction.at(i-1)
0091                                         ->GetParameter(j));
0092       }
0093     }
0094
0095     param2 = GK_Raw->GetV2()[n-1];
0096     AsyFunction.at(i)->FixParameter(2,param2);
0097
0098     gtemp = new TGraph(n, GK_Raw->GetV1(), GK_Raw->GetV2());
0099     gtemp -> Fit(tempname1, "MQ");
0100   }
0101
0102   return 0;
0103
0104
0105 }
0106
0107 int Asymmetry::Parameterize()
0108 {
0109   int n_all_Qsq = GK_Raw->Draw("q2","","goff");
0110   double * all_Qsq = GK_Raw->GetV1();
0111   Qsq_Vec = *(new vector<double>(all_Qsq, all_Qsq + n_all_Qsq));
0112   sort(Qsq_Vec.begin(), Qsq_Vec.end());
0113   vector<double>::iterator it = unique(Qsq_Vec.begin(),
0114                                        Qsq_Vec.end());
0115   Qsq_Vec.resize(distance(Qsq_Vec.begin(), it));
0116   nQsq = Qsq_Vec.size();
0117   Parameterize(Qsq_Vec);
0118   return 0;
0119 }
0120
0121 double Asymmetry::Extrap(double x0, double x1, double x2,
0122                          double y1, double y2)
0123 {
0124   if (x1==x2) {
0125     return y2;
0126   }
0127   else if (y1==y2) {
0128     return y2;
0129   }
0130
0131   else {
0132     Double_t m = (y1-y2)/(x1-x2);
0133     Double_t b = y1-(m*x1);
0134
0135     Double_t y0 = (m*x0)+b;
0136
0137     return y0;
0138
0139   }
0140
0141 }
0142
0143 double Asymmetry::GetAsyAmp(double Qsq, double tp)
0144 {
0145   double x0, x1=0, x2;
0146   double y0, y1, y2;
0147   x0 = Qsq;
0148   // Special case where only one Qsq was supplied to
0149   // Parameterize
0150   if (nQsq == 0){
0151     cerr << "Error: Qsq list not initialized" << endl;
0152     exit(EXIT_FAILURE);
0153   }
0154   if (nQsq == 1){
0155     cerr << "Warning: Only only one Qsq value in use." << endl;
0156     cerr << "No extrapolation between fits." << endl;
0157     return AsyFunction[0]->Eval(-tp);
0158   }
0159
0160   // Find nearest value in Qsq_Vec < Qsq
0161   int n_low = 0;
0162   int n_high = 0;
0163   for (n_low = 0; n_low < nQsq; n_low++){
0164     if (Qsq < Qsq_Vec[n_low+1]) break;
0165   }
0166
0167   // Special case where Qsq < any(Qsq_Vec)
0168   // (also cover Qsq_Vec[0] < Qsq < Qsq_Vec[1]
0169   if (n_low == 0)
0170     n_high = 1;
0171   // Special case where Qsq > any(Qsq_Vec)
0172   else if (n_low == nQsq){
0173     n_low = nQsq-2;
0174     n_high = nQsq-1;
0175   }
0176   //Usual case
0177   else
0178     n_high = n_low + 1;
0179
0180   x1 = Qsq_Vec[n_low];
0181   x2 = Qsq_Vec[n_high];
0182   y1 = AsyFunction[n_low]->Eval(-tp);
0183   //cout << "x1:\t" << x1 << endl;
0184   //cout << "y1:\t" << y1 << endl;
0185   y2 = AsyFunction[n_high]->Eval(-tp);
0186   //cout << "x2:\t" << x2 << endl;
0187   //cout << "y2:\t" << y2 << endl;
0188
0189   //cout<<x1<<"\t"<<x2<<"\t"<<y1<<"\t"<<y2<<endl;
0190
0191   y0 = Extrap(x0, x1, x2, y1, y2);
0192   return y0;
0193
0194 }
0195
0196 int Asymmetry::SetPars(vector<double> in_Qsq)
0197 {
0198   //Go to default work file if not extern not available
0199   if (AsymmFile->IsZombie()){
0200     AsymmFile = new TFile(Form("%s/data/input/Asymmetries.root", DEMPgen_Path));
0201   }
0202
0203   nQsq = in_Qsq.size();
0204   if (nQsq == 0) {
0205     return Parameterize();
0206   }
0207
0208   else{
0209     Qsq_Vec = in_Qsq;
0210   }
0211
0212   TTree * Pars = (TTree*)AsymmFile->Get(AsyNameStr);
0213
0214   char tfnamestr[100] = "%s_%d";
0215   char tempname1[100];
0216
0217   for (int i = 0; i < nQsq; i++){
0218     sprintf(tempname1, tfnamestr, AsyNameStr, i);
0219     AsyFunction.push_back(new TF1(tempname1, FuncForm));
0220   }
0221
0222   nPars = AsyFunction.at(0)->GetNpar();
0223
0224   double pars[nPars];
0225
0226   Pars->SetBranchAddress("pars",&pars);
0227
0228   for (int i=0; i<nQsq; i++){
0229     Pars->GetEntry(i);
0230     AsyFunction.at(i)->SetParameters(&pars[0]);
0231   }
0232
0233 }
0234
0235 void Asymmetry::PrintPars()
0236 {
0237   nPars = AsyFunction.at(0)->GetNpar();
0238   for (int i = 0; i < nQsq; i++){
0239     cout << "Qsq:\t"<< Qsq_Vec.at(i) << endl;
0240     for (int j = 0; j < nPars; j++){
0241       cout << "Par "<<j<<":\t"
0242            <<AsyFunction.at(i)->GetParameter(j)
0243            << endl;
0244     }
0245   }
0246 }