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