#include <iostream>
#include <iomanip>
#include <fstream>
#include <fnmatch.h>
#include <stdio.h>
#include <time.h>
#include "TROOT.h"
#include "TStyle.h"
#include "TChain.h"
#include "TFile.h"
#include "TDirectory.h"
//#include "qvalue_data.h"
#include "qvalue_fitter.h"
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH1I.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "RooAddPdf.h"
#include "RooPolynomial.h"
#include "RooArgusBG.h"
#include "RooChebychev.h"
#include "RooVoigtian.h"
#include "RooGaussian.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "TTree.h"
#include "TGaxis.h"
#include "TCanvas.h"

#include "TMarker.h"
#include "TPaveText.h"
#include "THStack.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TAxis.h"

#include "TLatex.h"
#include "TColor.h"
#include "RooTrace.h"

#include "analysis.h"
#include "histos.h"


//  //  //  //  //  //  //  //  //  //  //  //  //  //  //
//
// Define functions
//
void initialize();
void defineHistos();
void defineBranches(TTree*);
void defineBranches2(TTree*);
void finaltreefill();
void defineTarget(int);
void writeRoot(TFile*, TCanvas*, int, int, int, float, int, int);
TLorentzVector lorentzBoost(TLorentzVector,TLorentzVector);


//  //  //  //  //  //  //  //  //  //  //  //  //  //  //
//
//  Define variables
//
#define numberOfSituations 5
#define numberOfPeriods 7
#define targetType 4
#define numberOfBeamPol 3
#define numberOfThetaBins 11
#define numberOfPhiBins 21
#define numberOfNewPhiBins 1
#define numberOfWBins 22
#define numberOfEnergyBins 29

#define SQR(a) ((a)*(a))

using namespace std; 
using namespace RooFit ;
float m_p = 938.27203;   // mass of the proton
int numberofhistostostore = 50;


int main()
{
  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
  //
  //  The start time to check the running time
  //
  clock_t start = clock();

  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //                                                                                     
  //                                                                                     
  //  Setting gStyle                                                                                                                                 
  //
  gROOT -> SetBatch();
  gROOT -> SetStyle("Plain");

  gStyle -> SetOptStat(0);
  gStyle -> SetPalette(1);
  gStyle -> SetOptTitle(0);

  //  //  //  //  //  //  //  //
  //
  // Initialize variables
  //
  initialize();

  int Per_val = 9;

  int sit_limit = 9;

  int BPol_val = 9;

  float CL_cut = 0.0;

  float ebin_width = -9.0;  // in MeV
  float lowest_energy = -9.0; //in MeV

  int NewEBin_val = 9;
  int NewEBin = -1;
  int NewEBin_car = -1;

  int NewPBin_val = 9;
  int NewPBin = 9;
  int NewPBin_car = 9; 

  int TestRun = 1;
  int NumOfEvnt = 9;

  int meson_flag = -1;
  int run_flag = -1;
  omega_flag = -1;
  phi_flag = -1;
  eta_flag = -1;
  Kaon_flag = -1;
  Sigma_flag = -1;

  carbon_flag = -1;
  argusbg_flag = -1;

  binning_flag = -1;

  int finaltree_save = 0;

  //double signalFitmean = -999.0;
  //double signalFitsigma = -999.0;

  double par1 = -9.0;
  double par2 = -9.0;
  double par3 = -9.0;
  double par4 = -9.0;
  double par5 = -9.0;
  double par6 = -9.0;
  double par7 = -9.0;
  double par8 = -9.0;
  double par9 = -9.0; 

  int Emax =  int(lowest_energy + ebin_width*NewEBin_val); // in MeV

  if (omega_flag == 1) meson_flag = 1;
  else if (Kaon_flag == 1) meson_flag = 2;
  else if (eta_flag == 1) meson_flag = 3;
  else if (phi_flag == 1) meson_flag = 4;
  else if (Sigma_flag == 1) meson_flag = 5;
  else meson_flag = 0;

  if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
    if (sit_limit==1 || sit_limit==4) {
      beginfit = 0.;
      if (Emax<901) endfit = 270.0;
      else if (Emax>901 && Emax<1301) endfit = 400;
      else endfit = 450;

      if (Emax<1301) hist_Binning = 12;
      else hist_Binning = 14;
      PrintPlot = 49;
    }

    if (sit_limit==2) {
      beginfit = 0.;
      if (Emax<901) endfit = 270.0;
      else if (Emax>901 && Emax<1301) endfit = 400;
      else endfit = 450;

      if (Emax<1301) hist_Binning = 12;
      else hist_Binning = 14;
      PrintPlot = 49;
    }
  }

  else if (omega_flag == 1 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
    if (argusbg_flag == 0) {
      beginfit = 650;
      endfit= 900;
      hist_Binning = 120;
      PrintPlot = 49;
    }
    else if (argusbg_flag == 1) {
      beginfit = 650;
      endfit = 900;
      hist_Binning = 120;
      PrintPlot = 49;
    }
  }

  else if (Kaon_flag == 1 && phi_flag == 0 && omega_flag == 0 && eta_flag == 0 && Sigma_flag == 0) {
    beginfit = 473.0;
    endfit = 523.0;
    hist_Binning = 20;
    PrintPlot = 49;
  }

  else if (Kaon_flag == 0 && phi_flag == 0 && omega_flag == 0 && eta_flag == 0 && Sigma_flag == 1) {
    beginfit = 1149.0;
    endfit = 1229.0;
    hist_Binning = 32;
    PrintPlot = 49;
  }

  else if (Kaon_flag == 0 && phi_flag == 1 && omega_flag == 0 && eta_flag == 0 && Sigma_flag == 0) {
    beginfit = 980.0;
    endfit = 1080.0;
    hist_Binning = 12;
    PrintPlot = 49;
  }

  else if (Kaon_flag == 0 && phi_flag == 0 && omega_flag == 0 && eta_flag == 1 && Sigma_flag == 0) {
    //beginfit = 522.0;
    //endfit = 577.0;
    //hist_Binning = 22;
    beginfit = 455.0;
    endfit = 650.0;
    hist_Binning = 78;
    PrintPlot = 49;
  }
  else {
    cerr << "Warning: Your reaction is not properly defined." << endl
	 << "Neither (omega_flag, phi_flag, eta_flag, Kaon_flag, Sigma_flag) is set." << endl << endl;
    exit(0);
  }

  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
  //
  //  Qvalue txt
  //
  if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) numNN_def = 300;
  else if (omega_flag == 1) numNN_def = 1000;
  else if (phi_flag == 1) numNN_def = 300;
  else if (eta_flag == 1) numNN_def = 500;
  else if (Sigma_flag == 1) numNN_def = 1000;
  else if (Kaon_flag == 1) {
    if (run_flag == 2) numNN_def = 500;
    else if (run_flag == 0) numNN_def = 200;
  }

  if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) 
    sprintf(name,"outQvalue_Sit%02d_Per%02d_NewEBin%03d_BPol%02d_NewPBin%02d.txt",sit_limit,Per_val,NewEBin_val,BPol_val,NewPBin_val);
  else if (omega_flag == 1 || phi_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1)
    sprintf(name,"outQvalue_Sit%02d_Per%02d_NewEBin%03d_BPol%02d_NewPBin%02d.txt",5,Per_val,NewEBin_val,BPol_val,NewPBin_val);

  ofstream OutQvalue(name, ios::out);
  OutQvalue.setf(ios::showpoint|ios::fixed);
  OutQvalue.precision(5);


  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
  //
  // Output Root file
  //
  if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) 
    sprintf(name,"output_Sit%02d_Per%02d_NewEBin%03d_BPol%02d_NewPBin%02d.root",sit_limit,Per_val,NewEBin_val,BPol_val,NewPBin_val);
  else if (omega_flag == 1 || phi_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1)
    sprintf(name,"output_Sit%02d_Per%02d_NewEBin%03d_BPol%02d_NewPBin%02d.root",5,Per_val,NewEBin_val,BPol_val,NewPBin_val);

  outputfile = new TFile(name,"RECREATE");

  //////////////////////////////
  // Defining the directories //                                                                                               
  //////////////////////////////
  FitDir[0] = new TDirectoryFile("Miscellaneous","Miscellaneous");
  //FitDir[1] = new TDirectoryFile("2pion_analysis"," gamma p --> p pi+ pi-");
  //FitDir[2] = new TDirectoryFile("omega_analysis"," gamma p --> p omega -> p pi+ pi- (pi0)");
  //FitDir[3] = new TDirectoryFile("eta_analysis"," gamma p --> p eta -> p pi+ pi- (pi0)");
  //FitDir[4] = new TDirectoryFile("KS_analysis"," gamma p --> K Sigma -> p pi+ pi- (pi0)");
  //FitDir[5] = new TDirectoryFile("phi_analysis"," gamma p --> p phi -> p K+ K-");

  cout << "Is this a gamma p --> p omega analysis?   Yes (1), No (0): " << omega_flag << endl;
  cout << "Is this a gamma p --> p eta analysis?     Yes (1), No (0): " << eta_flag << endl;
  cout << "Is this a gamma p --> p phi analysis?     Yes (1), No (0): " << phi_flag << endl;
  cout << "Is this a gamma p --> <K> Sigma analysis? Yes (1), No (0): " << Kaon_flag << endl;
  cout << "Is this a gamma p --> K <Sigma> analysis? Yes (1), No (0): " << Sigma_flag << endl << endl;

  cout << "Are we describing the background with an Argus function?" << endl
       << "(needed for root(s) bin near omega threshold) Y/N (1/0): " << argusbg_flag << endl << endl;

  if (argusbg_flag == 1 && (eta_flag == 1 || phi_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1)) {
    cerr << "Warning: Argus (background) function currently not properly set up in combination with " << endl
	 << "either (phi_flag, eta_flag, Kaon_flag, Sigma_flag) or the two-pion channel." << endl << endl;
    exit(0);
  }

  cout << "Run group: g9a (0), g9b (1), g12 (2), g8b (3): " << run_flag << endl;

  if (omega_flag == 1 || Kaon_flag == 1 || eta_flag == 1 || Sigma_flag == 1) {
    cout << "CL cut applied for missing-pi0 hypothesis: " << CL_cut << endl << endl;
  } 
  else if (phi_flag == 1) {
    cout << "CL cut applied for K+ K- hypothesis: " << CL_cut << endl << endl;
  }

  cout << "Fit range defined from " << beginfit << " to " << endfit << " [MeV]." << endl;
  cout << "The number of nearest events used for the analysis is: " << numNN_def << "." << endl << endl;
 
  cout << "Binning scenario: E (1), W (2): " << binning_flag << endl << endl;

  cout << "Energy binning: (1) ebin_width = " << ebin_width << " [MeV], (2) lowest energy = " << lowest_energy << " [MeV]" << endl;
  cout << "       NewEBin[ " << NewEBin_val << " ] " << endl;

  if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
    cout << "E_max in the defined energy bin = " << Emax << endl;
    cout << "Define fit range from " << beginfit << " to " << endfit << endl;
    cout << "endfitrange = " << endfit << endl;
    cout << "  Period[ " << Per_val << " ] " << endl;
    cout << "  Situation[ " << sit_limit << " ] " << endl;
    cout << "  BPol type[ " << BPol_val << " ] " << endl;
  }

  cout << "       NewPBin[ " << NewPBin_val << " ]  -->  NewPBin[1] = Phi*[0,2*M_PI]" << endl << endl;
  cout << "Are we reading in any carbon data? Y/N (1/0) " << carbon_flag << endl;

  if (TestRun) cout << endl << "This is a test run with " << NumOfEvnt << " events." << endl;
   
  for (int period_limit=Per_val; period_limit<Per_val+1; period_limit++) {
    for (int newebin_limit=NewEBin_val; newebin_limit<NewEBin_val+1; newebin_limit++) {
      for (int bpol_limit=BPol_val; bpol_limit<BPol_val+1; bpol_limit++) {
	for (int newpbin_limit=NewPBin_val; newpbin_limit<NewPBin_val+1; newpbin_limit++) {
	   
          // Define histograms
	  defineHistos();
	     
	  cout << endl << "Summary: Period[" << period_limit << "] Sit[" << sit_limit << "] BPoltype[" << bpol_limit 
	       << "] NewEBin[" << newebin_limit << "] NewPBin[" << newpbin_limit << "] " << endl << endl << endl;
	  
	  counter_But = 0;
	  count_total = 0;
	  count_top = 0;
	  counter_Car = 0;
	  
	  if (run_flag == 0 || run_flag == 1) {
	    cout << " ////////////////////////////// " << endl;
	    cout << " //                          // " << endl;
	    cout << " //   Data (butanol) input   // " << endl;
	    cout << " //                          // " << endl;
	    cout << " ////////////////////////////// " << endl;
	  }
	  else {
            cout << " ////////////////////////////////////// " << endl;
            cout << " //                                  // " << endl;
            cout << " //   Data (liquid hydrogen) input   // " << endl;
            cout << " //                                  // " << endl;
            cout << " ////////////////////////////////////// " << endl;
	  }

	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
          //
          //  Define tree to hold fit variables
          //
          TTree *vartree = new TTree("Qfactree", "Qfactree");

	  // Define branches
	  defineBranches(vartree);

	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	  //
	  // Define a tree to hold variables for determining correlation matrix
	  //
	  TTree *finaltree = new TTree("finaltree", "finaltree");
	  int NNarray_nt[numNN_def];

	  if (finaltree_save == 1) {
	    for (int i=0; i<numNN_def; i++) NNarray_nt[i]=-1;

	    if (omega_flag==0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) 
	      finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[300]/I");
	    else if (phi_flag == 1 || (Kaon_flag == 1 && run_flag == 2)) 
	      finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[500]/I");
            else if (Kaon_flag == 1 && run_flag == 0)
              finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[200]/I");
	    else if (omega_flag == 1)
	      finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[1000]/I");
	    else if (eta_flag == 1)
              finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[500]/I");
            else if (Sigma_flag == 1)
              finaltree->Branch ("NNarray_nt",NNarray_nt, "NNarray_nt[1000]/I");

	    defineBranches2(finaltree);
	  }

	  // Read the input text file
	  ifstream inputfile;
	  sprintf(name,"../Qvalue_Per%02d.txt",period_limit);
	  inputfile.open(name,ios::in);
	  cout << " Get the data from " << name << endl << endl;

	  while (!inputfile.eof()) {
	    
	    for (int i=0; i<36; i++) inputfile >> temp[i]; 

  	    //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	    //
	    //  Initialize arrays for each event
            //
	    run_number = (int)temp[0];
	    event_num = (int)temp[1];  
	    topology = (int)temp[2];

	    if (temp[3]==0 || temp[3]==1 || temp[3]==2) {
	      // It means that the 4th entry in the text file is the beam pol. type.
	      // cerr << "temp[3] = " << temp[3] << endl;
	      bpol = (int)temp[3];
	      deg_bpol = 0.00;
	    }
	    else {
	      // The 4th entry in the text file is the deg. of beam pol. per event
	      bpol = 0;
	      deg_bpol = temp[3];
	    }
	    // cerr << "bpol = " << bpol << endl << "deg bpol = " << deg_bpol << endl;

	    realphi = temp[4];
	    cosrealtheta = temp[5];
            coscmstheta_pro = temp[6];

	    pho_energy = temp[7];
	    pro_px = temp[8];                              
            pro_py = temp[9];                                                           
            pro_pz = temp[10];                                                           
            pro_E = temp[11];                                                               

            pip_px = temp[12];                                                   
            pip_py = temp[13];                                                       
            pip_pz = temp[14];                                                     
            pip_E = temp[15];                                                         

            pim_px = temp[16];                                                      
            pim_py = temp[17];                                                   
            pim_pz = temp[18];       
	    pim_E = temp[19];

	    cl_2pi = temp[20];
	    cl_3pi = temp[21];

	    z_vertex = temp[22];

	    pho_energy_kinf = temp[23];                                                       

	    pro_px_kinf = temp[24];                                                              
	    pro_py_kinf = temp[25];                                                               
	    pro_pz_kinf = temp[26];                                                               
	    pro_E_kinf = temp[27];                                                                

	    pip_px_kinf = temp[28];                                     
            pip_py_kinf = temp[29];                                                                  
            pip_pz_kinf = temp[30];                                                                 
            pip_E_kinf = temp[31];                                                                   

            pim_px_kinf = temp[32];                                                                  
            pim_py_kinf = temp[33];                                                                  
            pim_pz_kinf = temp[34];                                                                  
            pim_E_kinf = temp[35];

	    if (run_flag < 2) defineTarget(run_flag);
	    else if (run_flag == 2 || run_flag == 3) target = 1; // Events are written to the g12 or g8b text file after a vertex cut.
	    //cerr << "z_vertex = " << z_vertex << endl << " target = " << target << endl;

            if (binning_flag == 1) {
	      NewEBin = int( ((pho_energy_kinf - lowest_energy) / ebin_width)+1 ); // all energies in MeV
	    }
	    else if (binning_flag == 2) {
	      Wenergy = TMath::Sqrt(2*pho_energy_kinf*938.27 + 938.27*938.27);
              NewEBin = int( ((Wenergy - lowest_energy) / ebin_width)+1 );
	    } 
	    //cerr << "pho_energy = " << pho_energy << "ebin_width = " << ebin_width << endl << "NewEBin = " << NewEBin << endl;

	    NewPBin = 5;
	    //cerr << "realphi = " << realphi << endl;
           
	    if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
	      if (realphi >= 0.0 && realphi <= 2*M_PI) NewPBin = 1; // i.e. phi* from 0 to 360 deg

	      if (realphi < 0) 
		cerr << "WARNING: realphi is negative! Check and change the range of realphi - 0 to 2pi or -pi to pi ?" << endl;
	    }
	    else if (omega_flag == 1 || phi_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1) NewPBin = 1; // realphi in the textfile may or may not be the correct one.
	                                                                                                                  // We will calculate realphi in this code.

	    //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	    //
	    //  The data (hydrogen or butanol) input
	    //
	    if (target == 1) count_total++;  
	    if (topology == sit_limit && bpol == bpol_limit && NewEBin == newebin_limit && NewPBin == newpbin_limit && target == 1) {
	      count_top++;

	      if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
	        //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	        //
	        //  The definition for the first boost
	        //
	        MyProton.SetPxPyPzE( pro_px,pro_py, pro_pz, pro_E );
	        MyPip.SetPxPyPzE( pip_px, pip_py, pip_pz, pip_E );
	        MyPim.SetPxPyPzE( pim_px,pim_py, pim_pz, pim_E );
	        MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
	      
                kinevar1 = pho_energy;
	        kinevar2 = cosrealtheta;
	        kinevar3 = realphi;
	        kinevar4 = coscmstheta_pro; 
                //kinevar5 = ;
		//kinevar6 = ;

	        if (sit_limit==1) fitvar = MyPim.M();
	        else if (sit_limit==2) fitvar = MyPip.M();
	        else if (sit_limit==3) fitvar = MyProton.M();
	        else if (sit_limit==4) {
		  My_Initial.SetPxPyPzE(0.0,0.0,pho_energy,pho_energy+938.27);
		  MM = My_Initial - MyProton - MyPip;
		  fitvar = MM.M();
		}

		butanolhisto[0]->Fill(fitvar);
		 
                if (fitvar < 0. || fitvar > endfit) massflag = 2;
		else {
		  if (fitvar < beginfit || fitvar > endfit ) massflag = 0;
		  else massflag = 1;
		}

		if (massflag < 2) vartree->Fill(); // Fill the butanol events which pass the top, ebin, phibin and mass cuts
		if (fitvar < beginfit || fitvar > endfit) {
		  totalhisto[0] -> Fill(fitvar);
		  backgroundhisto[0] -> Fill(fitvar);
		  continue;
		}
		counter_But++;
	      }
	         
	      else if (omega_flag == 1 || eta_flag == 1) {
		if (cl_3pi > CL_cut) {
		 
		  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
		  //
		  // Boosting to define kinematic variables (helicity, cos(theta), etc.)
                  // NOTE: Here the variables sometimes have 'omega' int their names. Basically it 
		  //       refers to the particle X in gamma p -> p X and applies to eta, too.
                  MyProton.SetPxPyPzE( pro_px_kinf,pro_py_kinf, pro_pz_kinf, pro_E_kinf );
                  MyPip.SetPxPyPzE( pip_px_kinf, pip_py_kinf, pip_pz_kinf, pip_E_kinf );
                  MyPim.SetPxPyPzE( pim_px_kinf, pim_py_kinf, pim_pz_kinf, pim_E_kinf );
                  MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
                  My_Initial.SetPxPyPzE(0.0,0.0,pho_energy_kinf, pho_energy_kinf+938.27);
                  MyPi0 = My_Initial - MyProton - MyPip - MyPim;
                  MM = MyPip + MyPim + MyPi0;
		  omega_pvec = MM.Vect();
		  omega_p = omega_pvec.Mag();
                  Wenergy = TMath::Sqrt(2*pho_energy_kinf*938.27 + 938.27*938.27);
                  // cerr << "omega 4-vector: (" << MM.Px() << "," << MM.Py() << "," << MM.Pz() << "," << MM.E()<< ")" << endl;
		  // cerr << "omega momentum = " << omega_p << endl;
                          
                  // Boosting to the overall c.m. frame -

		  // Method 1 (Priya) -
		  /*
		  TLorentzVector boostVector(0.0, 0.0, -My_Initial.Pz(), My_Initial.E());
		  omega_cm1 = lorentzBoost(boostVector, MM);
		  pip_cm1 = lorentzBoost(boostVector, MyPip);
		  pim_cm1 = lorentzBoost(boostVector, MyPim);
		  pro_rest_cm1 = lorentzBoost(boostVector, MyProton_Rest);
		  double costheta_omega1 = omega_cm1.CosTheta();
		  */

		  // Method 2 (Chris) -
		  TLorentzRotation cm_r;
		  TVector3 cm_v = My_Initial.BoostVector();
		  cm_r.Boost(-cm_v);
		  TLorentzVector omega_cm2 = cm_r * MM;
		  TLorentzVector pip_cm2 = cm_r * MyPip;
		  TLorentzVector pim_cm2 = cm_r * MyPim;
		  TLorentzVector pi0_cm2 = cm_r * MyPi0;
		  double costheta_omega2 = omega_cm2.CosTheta();
		  /*
		  cerr << "Results of first boost (comparing the 2 methods) -" << endl << endl;
		  cerr << "Priya's method - " << endl;
		  cerr << "       omega_cm1(" << omega_cm1.Px() << "," << omega_cm1.Py() << "," << omega_cm1.Pz() << "," << omega_cm1.E() << ")" << endl;
		  cerr << "       pip_cm1(" << pip_cm1.Px() << "," << pip_cm1.Py() << "," << pip_cm1.Pz() << "," << pip_cm1.E() << ")" << endl;
		  cerr << "       pim_cm1(" << pim_cm1.Px() << "," << pim_cm1.Py() << "," << pim_cm1.Pz() << "," << pim_cm1.E() << ")" << endl;
		  cerr << "       costheta_omega1 = " << costheta_omega1 << endl << endl; 
			     
		  cerr << "Chris's method - " << endl;
		  cerr << "       omega_cm2(" << omega_cm2.Px() << "," << omega_cm2.Py() << "," << omega_cm2.Pz() << "," << omega_cm2.E() << ")" << endl;
                  cerr << "       pip_cm2(" << pip_cm2.Px() << "," << pip_cm2.Py() << "," << pip_cm2.Pz() << "," << pip_cm2.E() << ")" << endl;
                  cerr << "       pim_cm2(" << pim_cm2.Px() << "," << pim_cm2.Py() << "," << pim_cm2.Pz() << "," << pim_cm2.E() << ")" << endl;
                  cerr << "       costheta_omega2 = " << costheta_omega2 << endl << endl;

                  cerr << "from Volker's text file -  " << endl;
		  cerr << " costheta_recoil_pro = " << coscmstheta_pro << endl;
		  cerr << " This should be same as costheta_omega " << endl;
		  */

		  // Boosting to omega rest frame - 
                            
		  // Method 1 (Priya's method) -
		  /*
		  TVector3 zPrime1 = omega_cm1.Vect();
		  TVector3 yPrime1 = zPrime1.Cross(pro_rest_cm1.Vect()); 
		  TVector3 xPrime1 = yPrime1.Cross(zPrime1);

		  TLorentzVector boostVector2(-omega_cm1.Px(), -omega_cm1.Py(), -omega_cm1.Pz(), omega_cm1.E());        
		  TLorentzVector omega_hel1 = lorentzBoost(boostVector2, omega_cm1);
		  TLorentzVector pip_hel1 = lorentzBoost(boostVector2, pip_cm1);
		  TLorentzVector pim_hel1 = lorentzBoost(boostVector2, pim_cm1);
		  TVector3 normal = ( pip_hel1.Vect() ).Cross( pim_hel1.Vect() );

		  double costheta1 = (( ( pip_hel1.Vect() ).Cross( pim_hel1.Vect() ) ).Dot(zPrime1) )/( (zPrime1).Mag() * ( ( pip_hel1.Vect() ).Cross( pim_hel1.Vect() ) ).Mag() );
		  double cosphi1 = ( (zPrime1.Cross(normal)).Dot(yPrime1) )/( (yPrime1).Mag()*(zPrime1.Cross(normal)).Mag() );
		  double lambda1 = ( (pip_hel1.Vect()).Cross(pim_hel1.Vect()).Mag2() )/(( pip_hel1.Vect() ).Mag2() * ( pim_hel1.Vect() ).Mag2() );
		  */

		  // Method 2 (Chris' method) - 
		  // Rotating the particle's mom vectors in c.m. frame in such a way that omega mom. vector falls along z
		  TLorentzRotation o_r;
		  Float_t r1 = omega_cm2.Phi();
		  Float_t r2 = omega_cm2.Theta();
		  o_r.RotateZ(-r1);
                  o_r.RotateY(-r2);

		  TLorentzVector pip_tmp = o_r * pip_cm2;
		  TLorentzVector pim_tmp = o_r * pim_cm2;
		  TLorentzVector omega_tmp = o_r * omega_cm2;
		  TLorentzVector pi0_tmp = o_r * pi0_cm2;

		  TLorentzRotation helo_r;
		  TVector3 helo_v = omega_tmp.BoostVector();
		  helo_r.Boost(-helo_v);

		  TLorentzVector pip_hel2 = helo_r * pip_tmp;
		  TLorentzVector pim_hel2 = helo_r * pim_tmp; 
	          TLorentzVector pi0_hel2 = helo_r * pi0_tmp;

		  double costheta2 = (pip_hel2.Vect().Cross(pim_hel2.Vect())).CosTheta(); 
		  double phi2 = ( pip_hel2.Vect().Cross(pim_hel2.Vect()) ).Phi();
		  double cosphi2 = cos(phi2);
		  double lambda2_uncor = (pip_hel2.Vect().Cross(pim_hel2.Vect())).Mag2();

		  Float_t T_pip = pip_hel2.E() - pip_hel2.Mag(); // K.E. of pi+
		  Float_t T_pim = pim_hel2.E() - pim_hel2.Mag(); // K.E. of pi-
		  Float_t T_pi0 = pi0_hel2.E() - pi0_hel2.Mag(); // K.E. of pi0
	          Float_t Q = T_pip + T_pim + T_pi0; // total K.E. in omega rest frame
	     	  Float_t M = pip_hel2.Mag();
		  Float_t lambda2_max = Q*Q*((Q*Q/108.0)+(Q*M/9.0)+(M*M/3.0)); // lamdba is max when the 3 pions are
			                                                       // on the same plane and are 120 deg. apart from each other.

		  //Note: lambda max is not |pi+|^2 * |pi-|^2  (i.e., when they go perpendicular to each other, because of E-p conservation constraints with pi0)
		  Float_t lambda2 = lambda2_uncor/lambda2_max;
		  total_KinE_histo -> Fill(Q);

		  /*
		  cerr << "Results from second boost - " << endl << endl;
		  cerr << "Priya's method - " << endl;
			     
		  cerr << "    pip_hel(" << pip_hel1.Px() << "," << pip_hel1.Py() << "," << pip_hel1.Pz() << "," << pip_hel1.E() << ")" << endl;
		  cerr << "    pim_hel(" << pim_hel1.Px() << "," << pim_hel1.Py() << "," << pim_hel1.Pz() << "," << pim_hel1.E() << ")" << endl;
		  cerr << "    costheta = " << costheta1 << endl;
		  cerr << "    cosphi = " << cosphi1 << endl;
		  cerr << "    lambda = " << lambda1 << endl << endl;

		  cerr << "Chris' method - " << endl;
                  cerr << "    pip_hel(" << pip_hel2.Px() << "," << pip_hel2.Py() << "," << pip_hel2.Pz() << "," << pip_hel2.E() << ")" << endl;
                  cerr << "    pim_hel(" << pim_hel2.Px() << "," << pim_hel2.Py() << "," << pim_hel2.Pz() << "," << pim_hel2.E() << ")" << endl;
                  cerr << "    costheta = " << costheta2 << endl;
                  cerr << "    cosphi = " << cosphi2 << endl;
                  cerr << "    lambda = " << lambda2 << endl << endl;
		  */

		  fitvar = MM.M();

		  kinevar1 = costheta_omega2;
                  kinevar2 = costheta2;
                  kinevar3 = phi2;
                  kinevar4 = lambda2;
		  kinevar5 = MM.Phi();

		  if (binning_flag == 1)
		    kinevar6 = pho_energy_kinf;
		  else if (binning_flag == 2)
		    kinevar6 = Wenergy;

                  kinevar1_histo->Fill(kinevar1);
                  kinevar2_histo->Fill(kinevar2);
                  kinevar3_histo->Fill(kinevar3);
                  kinevar4_histo->Fill(kinevar4);
                  kinevar5_histo->Fill(kinevar5);
                  kinevar6_histo->Fill(kinevar6);

		  if (omega_flag == 1) {
		    butanolhisto[1]->Fill(fitvar);
		    //if (fitvar < 650.0 || fitvar > 950.0) massflag = 2;
                    if (fitvar < 650.0 || fitvar > 900.0) massflag = 2;
		    else {
		      if (fitvar < beginfit || fitvar > endfit) massflag = 0;
		      else massflag = 1;
		    }

		    if (massflag < 2) vartree->Fill(); // Fill the butanol events which pass the top, ebin, phibin and mass cuts
                
		    if (fitvar < beginfit || fitvar > endfit) {
		      totalhisto[1]->Fill(fitvar);
		      backgroundhisto[1]->Fill(fitvar);

		      continue;
		    }
		  }
		  else if (eta_flag == 1) {
		    butanolhisto[3]->Fill(fitvar);
		    //if (fitvar < 522.0 || fitvar > 577.0) massflag = 2;
                    if (fitvar < 455.0 || fitvar > 650.0) massflag = 2;
		    else {
		      if (fitvar < beginfit || fitvar > endfit) massflag = 0;
		      else massflag = 1;
		    }

		    if (massflag < 2) vartree->Fill(); // Fill the butanol events which pass the top, ebin, phibin and mass cuts    

		    if (fitvar < beginfit || fitvar > endfit) {
		      totalhisto[3]->Fill(fitvar);
		      backgroundhisto[3]->Fill(fitvar);

		      continue;
		    }
		  }

		  counter_But++;
		}
	      }

	      else if (Kaon_flag == 1 || Sigma_flag == 1) {
		if (cl_3pi > CL_cut) {

                  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
                  //
                  // Boosting to define kinematic variables (helicity, cos(theta), etc.)
                  // NOTE: Here the variables sometimes have 'omega' int their names. Basically it
                  //       refers to the particle X in gamma p -> p X
		  MyProton.SetPxPyPzE( pro_px_kinf, pro_py_kinf, pro_pz_kinf, pro_E_kinf );
		  MyPip.SetPxPyPzE( pip_px_kinf, pip_py_kinf, pip_pz_kinf, pip_E_kinf );
		  MyPim.SetPxPyPzE( pim_px_kinf, pim_py_kinf, pim_pz_kinf, pim_E_kinf );
		  MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
		  My_Initial.SetPxPyPzE(0.0,0.0,pho_energy_kinf, pho_energy_kinf+938.27);
		  MyPi0 = My_Initial - MyProton - MyPip - MyPim;
		  Wenergy = TMath::Sqrt(2*pho_energy_kinf*938.27 + 938.27*938.27);
		  MM = MyPip + MyPim;
		  Kaon_pvec = MM.Vect();
		  Kaon_p = Kaon_pvec.Mag();
		  Pboost = MyProton;
		  p0boost = MyPi0;

		  TVector3 b;
		  b = My_Initial.BoostVector();
		  p0boost.Boost(-b);
		  Pboost.Boost(-b);

		  // Boosting to the overall c.m. frame -

		  // Method 2 (Chris) -
		  TLorentzRotation cm_r;
		  TVector3 cm_v = My_Initial.BoostVector();
		  cm_r.Boost(-cm_v);

                  TLorentzVector proton_cm = cm_r * MyProton;
                  TLorentzVector proton_rest_cm = cm_r * MyProton_Rest;
		  TLorentzVector Kaon_cm = cm_r * MM;
		  TLorentzVector pip_cm = cm_r * MyPip;
		  TLorentzVector pim_cm = cm_r * MyPim;
		  TLorentzVector pi0_cm = cm_r * MyPi0;

                  TVector3 pip3_cm(pip_cm.Px(), pip_cm.Py(), pip_cm.Pz());
		  TVector3 pim3_cm(pim_cm.Px(), pim_cm.Py(), pim_cm.Pz());
                  TVector3 pi03_cm(pi0_cm.Px(), pi0_cm.Py(), pi0_cm.Pz());
                  TVector3 proton3_cm(proton_cm.Px(), proton_cm.Py(), proton_cm.Pz());
                  TVector3 proton_rest3_cm(proton_rest_cm.Px(), proton_rest_cm.Py(), proton_rest_cm.Pz());

		  double costheta_Kaon = Kaon_cm.CosTheta();
//
  TLorentzVector MyPhoton;
  MyPhoton.SetPxPyPzE(0.0,0.0,pho_energy_kinf, pho_energy_kinf);
  TLorentzVector MyPhoton_cm2 = cm_r * MyPhoton;
  TVector3 zprime = (MyPhoton_cm2.Vect().Cross((proton_cm + pi0_cm).Vect())).Unit();
  TLorentzRotation cm_s;
  TVector3 cm_vs = (proton_cm + pi0_cm).BoostVector();
 
		  cm_s.Boost(-cm_vs);
		  TLorentzVector proton_cms = cm_s * proton_cm;
  TVector3 proton_angle = (proton_cms.Vect()).Unit();

double costheta_ad = proton_angle.Dot(zprime);
//

                  

                  TVector3 zPrime = proton3_cm + pi03_cm;
		  TVector3 yPrime = proton_rest3_cm.Cross(proton3_cm + pi03_cm);
		  TVector3 xPrime = yPrime.Cross(zPrime);
		  float cosRealPhi = ( yPrime.Dot(pi03_cm.Cross(proton3_cm)) ) / ( (yPrime).Mag() * (pi03_cm.Cross(proton3_cm)).Mag() );
		  float RealPhi = 0.0;

		  boostVector.SetPxPyPzE(-zPrime.X(), -zPrime.Y(), -zPrime.Z(), (proton_cm + pi0_cm).E() );
		  TLorentzVector proton_Prime = lorentzBoost(boostVector, proton_cm);
		  TLorentzVector pi0_Prime = lorentzBoost(boostVector, pi0_cm);

		  TVector3 proton3_Prime(proton_Prime.Px(), proton_Prime.Py(), proton_Prime.Pz());
		  TVector3 pi03_Prime(pi0_Prime.Px(), pi0_Prime.Py(), pi0_Prime.Pz());

		  float cosRealTheta = cos(proton3_Prime.Angle(zPrime)); // range from -1.0 to 1.0

		  if ( acos( pi03_cm.Cross(proton3_cm).Dot(xPrime) / ( (pi03_cm.Cross(proton3_cm)).Mag() * xPrime.Mag() ) ) < (M_PI/2) ) {
		    RealPhi = -acos(cosRealPhi);
		  }
		  else {
		    RealPhi =  acos(cosRealPhi);
		  }

                  // Boosting to the pi+ pi- rest frame -

		  // Method a la Alvin -
		  /*TVector3 zPrime = pip3_cm + pim3_cm;
		  TVector3 yPrime = proton_rest3_cm.Cross(proton3_cm + pi03_cm);
		  TVector3 xPrime = yPrime.Cross(zPrime);
		  float cosRealPhi = ( yPrime.Dot(pim3_cm.Cross(pip3_cm)) ) / ( (yPrime).Mag() * (pim3_cm.Cross(pip3_cm)).Mag() );
		  float RealPhi = 0.0;

		  boostVector.SetPxPyPzE(-zPrime.X(), -zPrime.Y(), -zPrime.Z(), (pip_cm + pim_cm).E() );
		  TLorentzVector pip_Prime = lorentzBoost(boostVector, pip_cm);
		  TLorentzVector pim_Prime = lorentzBoost(boostVector, pim_cm);

		  TVector3 pip3_Prime(pip_Prime.Px(), pip_Prime.Py(), pip_Prime.Pz());
		  TVector3 pim3_Prime(pim_Prime.Px(), pim_Prime.Py(), pim_Prime.Pz());

		  float cosRealTheta = cos(pip3_Prime.Angle(zPrime)); // range from -1.0 to 1.0

		  if ( acos( pim3_cm.Cross(pip3_cm).Dot(xPrime) / ( (pim3_cm.Cross(pip3_cm)).Mag() * xPrime.Mag() ) ) < (M_PI/2) ) {
		    RealPhi = -acos(cosRealPhi);
		  }
		  else {
		    RealPhi =  acos(cosRealPhi);
		  }*/
		  
		  // Defining metric variables
                  if (Kaon_flag == 1) MM = MyPip + MyPim;
                  else if (Sigma_flag == 1) MM = MyProton + MyPi0;
		  fitvar = MM.M();

		  if (binning_flag == 1)
		    kinevar1 = pho_energy_kinf;
		  else if (binning_flag == 2)
		    kinevar1 = Wenergy;

		  kinevar2 = cosRealTheta;
                  kinevar3 = RealPhi;
                  kinevar4 = costheta_Kaon;
                  kinevar5 = MM.Phi();
                  kinevar6 = TMath::Cos(Pboost.Angle(p0boost.Vect()));
                  kinevar7 = costheta_ad;

                  kinevar1_histo->Fill(kinevar1);
		  kinevar2_histo->Fill(kinevar2);
		  kinevar3_histo->Fill(kinevar3);
		  kinevar4_histo->Fill(kinevar4);
		  kinevar5_histo->Fill(kinevar5);
                  kinevar6_histo->Fill(kinevar6);

		  if (Kaon_flag == 1) {
		    butanolhisto[2]->Fill(fitvar);
		    if (fitvar < 473.0 || fitvar > 523.0) massflag = 2;
		    else {
		      if (fitvar < beginfit || fitvar > endfit) massflag = 0;
		      else massflag = 1;
		    }

		    if (massflag < 2) vartree->Fill(); // Fill the butanol events which pass the top, ebin, phibin and mass cuts

		    if (fitvar < beginfit || fitvar > endfit) {
		      totalhisto[2]->Fill(fitvar);
		      backgroundhisto[2]->Fill(fitvar);

  		      continue;
		    }
		  }
		  else if (Sigma_flag == 1) {
                    butanolhisto[5]->Fill(fitvar);
                    if (fitvar < 1149.0 || fitvar > 1229.0) massflag = 2;
                    else {
                      if (fitvar < beginfit || fitvar > endfit) massflag = 0;
                      else massflag = 1;
                    }

                    if (massflag < 2) vartree->Fill(); // Fill the butanol events which pass the top, ebin, phibin and mass cuts                                                                          

                    if (fitvar < beginfit || fitvar > endfit) {
                      totalhisto[5]->Fill(fitvar);
                      backgroundhisto[5]->Fill(fitvar);

                      continue;
                    }
		  }

		  counter_But++;
		}
	      }
	    }

	    //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	    //
	    //  The carbon input
	    //
	    if (run_flag == 1 && carbon_flag == 1) {    
	      if (topology == sit_limit && bpol == bpol_limit && NewEBin == newebin_limit && NewPBin == newpbin_limit && target == 2) {

                // // // // // // // // // // // // // // // // // // // // // // // // // // // // // //
	        //
	        // The definition for the first boost
	        //
	        MyProton.SetPxPyPzE( pro_px,pro_py, pro_pz, pro_E );
                MyPip.SetPxPyPzE( pip_px, pip_py, pip_pz, pip_E );
                MyPim.SetPxPyPzE( pim_px,pim_py, pim_pz, pim_E );
                MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );

		labphi_pip_car = MyPip.Phi();
                labphi_pim_car = MyPim.Phi();

		if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
		  if (sit_limit==1) MM_Car = MyPim.M();
	          else if (sit_limit==2) MM_Car = MyPip.M();
	          else if (sit_limit==3) MM_Car = MyProton.M();
	          else if (sit_limit==4) {
		    My_Initial.SetPxPyPzE(0.0,0.0,pho_energy,pho_energy+938.27);
		    MM = My_Initial - MyProton - MyPip;
		    MM_Car = MM.M();
		  }
		  carbonhisto[0]->Fill(MM_Car);
		  if (MM_Car > 0.0 && MM_Car < 150.0) {
                    labphi_pip_Car_histo[0]->Fill(labphi_pip_car);
                    labphi_pim_Car_histo[0]->Fill(labphi_pim_car);
                  }

                  if (MM_Car > 350.0 && MM_Car < 500.0) {
                    labphi_pip_Car_MM350up_histo[0]->Fill(labphi_pip_car);
                    labphi_pim_Car_MM350up_histo[0]-> Fill(labphi_pim_car);
                  }
		}

		else if (omega_flag == 1 || phi_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1) {
		  My_Initial.SetPxPyPzE(0.0,0.0,pho_energy,pho_energy+938.27);
		  MyPi0 = My_Initial - MyProton - MyPip - MyPim;
		  MM = MyPip + MyPim + MyPi0;
		  MM_Car = MM.M();

		  if (omega_flag == 1) carbonhisto[1]->Fill(MM_Car);
		  else if (Kaon_flag == 1) carbonhisto[2]->Fill(MM_Car);
		  else if (eta_flag == 1) carbonhisto[3]->Fill(MM_Car);
                  else if (phi_flag == 1) carbonhisto[4]->Fill(MM_Car);
                  else if (Sigma_flag == 1) carbonhisto[5]->Fill(MM_Car);
		}

	        if (MM_Car < beginfit || MM_Car > endfit) continue;
	      
	        counter_Car++; // You will notice that counter_Car will be less than the no. of events written in the
                               // text file for given W and PhiIndex. That's bec. of the above if statement ...
	      }
	    }
	  }

	  if (run_flag == 0 && carbon_flag == 1) { // Read-in g9b carbon for g9a analysis, g9a carbon has "ice" contamination                                           

	    cerr << "Read the input textfile (g9b circ. Carbon) for g9a background comparison - " << endl;

            ifstream inputfile_car;
            sprintf(name,"/d/grid11/hurley/root2text/kmkp2phitest.txt"); // hurley - read-in
	    inputfile_car.open(name,ios::in);
            cout << " Get the data from " << name << endl << endl;

	    while (!inputfile_car.eof()) {
              for (int i=0; i<36; i++) inputfile_car >> temp_car[i];

              topology_car = (int)temp_car[2];
              pho_energy_car = temp_car[7];
              pro_px_car = temp_car[8];
              pro_py_car = temp_car[9];
              pro_pz_car = temp_car[10];
              pro_E_car = temp_car[11];

              pip_px_car = temp_car[12];
              pip_py_car = temp_car[13];
              pip_pz_car = temp_car[14];
              pip_E_car = temp_car[15];

              pim_px_car = temp_car[16];
              pim_py_car = temp_car[17];
              pim_pz_car = temp_car[18];
              pim_E_car = temp_car[19];
              NewEBin_car = int( ((pho_energy_car-lowest_energy)/ebin_width)+1 ); // all energies in MeV
	      if (realphi >= 0. && realphi <= 2*M_PI) NewPBin_car = 1;
              else if (realphi < 0) 
		cerr << "WARNING: ealphi of g9b circ. carbon is negative!! Check and change the range of realphi - 0 to 2pi or -pi to pi?" << endl;

	      if (topology_car == sit_limit && NewEBin_car == newebin_limit && NewPBin_car == newpbin_limit) {

		// The definition for the first boost
		labphi_pip_car = -10.0;
		labphi_pim_car =-10.0;

		MyProton_car.SetPxPyPzE( pro_px_car,pro_py_car, pro_pz_car, pro_E_car );
		MyPip_car.SetPxPyPzE( pip_px_car, pip_py_car, pip_pz_car, pip_E_car );
		MyPim_car.SetPxPyPzE( pim_px_car,pim_py_car, pim_pz_car, pim_E_car );
		MyProton_Rest_car.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );

		labphi_pip_car = MyPip_car.Phi();
		labphi_pim_car = MyPim_car.Phi();

		if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
		  if (sit_limit==1) MM_Car = MyPim_car.M();
		  else if (sit_limit==2) MM_Car = MyPip_car.M();
		  else if (sit_limit==3) MM_Car = MyProton_car.M();
		  else if (sit_limit==4) {
		    My_Initial_car.SetPxPyPzE(0.0,0.0,pho_energy_car,pho_energy_car+938.27);
		    MM = My_Initial_car - MyProton_car - MyPip_car;
		    MM_Car = MM.M();
		  }
		  carbonhisto[0]->Fill(MM_Car);

		  if (MM_Car > 0.0 && MM_Car < 150.0) {
		    labphi_pip_Car_histo[0]->Fill(labphi_pip_car);
		    labphi_pim_Car_histo[0]->Fill(labphi_pim_car);
		  }

		  if (MM_Car > 350.0 && MM_Car < 500.0) {
		    labphi_pip_Car_MM350up_histo[0]->Fill(labphi_pip_car);
		    labphi_pim_Car_MM350up_histo[0]-> Fill(labphi_pim_car);
		  }

		  if (MM_Car < beginfit || MM_Car > endfit) continue;
		  counter_Car++;
		}
	      }
	    }
	  }

	  cerr << " //////////////////////////////////////// " << endl;
	  cerr << " //                                    // " << endl;
	  cerr << " //   Start the Q-factor calculation   // " << endl;
	  cerr << " //                                    // " << endl;
	  cerr << " //////////////////////////////////////// " << endl;
	  if (run_flag == 0 || run_flag == 1) {
  	    cerr << endl << "Total #events in the data textfile = " << count_total << endl;
	    cerr << endl << "Total #events for ( top = [" << sit_limit << "] BPol[" << bpol_limit 
                 << "] NewEBin[" << newebin_limit << "] NewPBin[" << newpbin_limit << "] ) = " << count_top << endl;  

	    cerr << endl << "Total #events in the mass range (" << beginfit << " - " << endfit << " MeV) = " << counter_But << endl 
                 << "  --> for ( Period[" << period_limit 
                 << "] Situation[" << sit_limit << "] BPol[" << bpol_limit << "] NewEBin[" << newebin_limit << "] NewPBin[" << newpbin_limit 
                 << "] ): Butanol[" << counter_But << "] Carbon[" << counter_Car << "]" << endl << endl;
	  }
	  else {
            cerr << endl << "Total #events in the data textfile = " << count_total << endl;
            cerr << endl << "Total #events for ( top = [" << sit_limit << "] BPol[" << bpol_limit
                 << "] NewEBin[" << newebin_limit << "] NewPBin[" << newpbin_limit << "] ) = " << count_top << endl;

            cerr << endl << "Total #events in the mass range (" << beginfit << " - " << endfit << " MeV) = " << counter_But << endl
                 << "  --> for ( Period[" << period_limit
                 << "] Situation[" << sit_limit << "] BPol[" << bpol_limit << "] NewEBin[" << newebin_limit << "] NewPBin[" << newpbin_limit
                 << "] ): Data[" << counter_But << "]" << endl << endl;
	  }

	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	  //  Reading data from the tree
	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	  int event = 0;
	  Int_t nentries = (Int_t)vartree->GetEntries();

     	  cerr << "Total no. of events stored in the tree = " << nentries << endl;

	  // Beginning of the Q-factor calculation loop
	  //
	  for (int counter=0; counter<nentries; counter++) {
	    vartree->GetEntry(counter);  // Get row for one event
	    //cerr << "event mass flag = " << massflag << endl << "If massflag = 0 go to the next event" << endl;

	    if (massflag != 1) continue; // Q-factor won't be calculated for this event
	    else {

              // Initialize
	      labphi_pip = -10.0;
	      labphi_pim = -10.0;

	      event++;
		  
              // Assign values to the other branches of finaltree
              finaltreefill();

              // Check whether the finaltree was assigned values properly -
	      cerr << endl << endl;
	      //cerr << " Run-number_nt[" << run_number_nt << "] " << endl;
              //cerr << " eventnum_nt[" << eventnum_nt << "] " << endl;
              //cerr << " Sit_nt[" << topology_nt << "] BPol_nt[" << bpol_nt << "]" << endl;
              //cerr << " RealPhi_nt[" << realphi_nt << "] photon-Energy_nt[" << pho_energy_nt << "] " << endl;
              //cerr << " proton_nt[ " << pro_px_nt << ", " << pro_py_nt << ", " << pro_pz_nt << ", " << pro_E_nt << " ] " << endl;
              //cerr << " pip_nt[ " << pip_px_nt << ", " << pip_py_nt << ", " << pip_pz_nt << ", " << pip_E_nt << " ] " << endl;
              //cerr << " pim_nt[ " << pim_px_nt << ", " << pim_py_nt << ", " << pim_pz_nt << ", " << pim_E_nt << " ] " << endl;
              //cerr << " CLvalue_2pi_nt[" << cl_2pi_nt << "] cosRealTheta_nt[" << cosrealtheta_nt << "] cosCMSTheta_pro_nt[" << coscmstheta_pro_nt << "] " << endl;
              //cerr << "Z.vertex_nt[" << z_vertex_nt << "]" << endl;
              //cerr << " photon-Energy-KinF_nt[" << pho_energy_kinf_nt << "] " << endl;
              //cerr << " proton-KinF_nt[ " << pro_px_kinf_nt << ", " << pro_py_kinf_nt << ", " << pro_pz_kinf_nt << ", " << pro_E_kinf_nt << " ] " << endl;
              //cerr << " pip-KinF_nt[ " << pip_px_kinf_nt << ", " << pip_py_kinf_nt << ", " << pip_pz_kinf_nt << ", " << pip_E_kinf_nt << " ] " << endl;
              //cerr << " pim-KinF_nt[ " << pim_px_kinf_nt << ", " << pim_py_kinf_nt << ", " << pim_pz_kinf_nt << ", " << pim_E_kinf_nt << " ] " << endl;

	      MyPip.SetPxPyPzE(pip_px_nt, pip_py_nt, pip_pz_nt, pip_E_nt);
	      MyPim.SetPxPyPzE(pim_px_nt, pim_py_nt, pim_pz_nt, pim_E_nt);

	      MyProton.SetPxPyPzE( pro_px_kinf_nt,pro_py_kinf_nt, pro_pz_kinf_nt, pro_E_kinf_nt );
              MyProton_Rest.SetPxPyPzE( 0.0, 0.0, 0.0, m_p );
              My_Initial.SetPxPyPzE(0.0,0.0,pho_energy_kinf_nt, pho_energy_kinf_nt+938.27);
              MyPi0 = My_Initial - MyProton - MyPip - MyPim;
              MM = MyPip + MyPim;
	      omega_pvec = MM.Vect();
              omega_p = omega_pvec.Mag();
		
	      labphi_pip = MyPip.Phi();
	      labphi_pim = MyPim.Phi();

	      //float seedeventarray[6];
              float seedeventarray[7];
	      float seedmass = -999.99;
	      int j = 0, change = 1;

	      float d2array[numNN_def];
	      float massarray[numNN_def];
	      int eventnumarray[numNN_def];

	      float max = -99.99;
	      int index = -1;
	      int iteration = 0;

	      seedeventarray[0] = kinevar1;
	      seedeventarray[1] = kinevar2;
	      seedeventarray[2] = kinevar3;
	      seedeventarray[3] = kinevar4;
	      seedeventarray[4] = kinevar5;
              seedeventarray[5] = kinevar6;
              seedeventarray[6] = kinevar7;

	      seedmass = fitvar;

	      // // // // // // // // // // // // // // // // // // // 
	      //
	      // Find nearest neighbors for this event# = counter
	      //
	      if (nentries <= numNN_def) {
		cerr << "the total # of entries is less than the no. of nearest neighbors needed." << endl;
		break;
	      }
	      else {
		for (int i=0; i<nentries; i++) {

		  if (i == counter) continue;
		  else {
		    vartree->GetEntry(i); // Read information of ith event
		    if (massflag == 2) continue; // not considering events with negative missing mass

		    else if (massflag == -1) {
                      cerr << "ERROR: massflag = -1!! Check why it was not assigned a proper value (0 - 2)" << endl;
		      continue;
		    }
		    else {
		     // float eventarray[6];
                       float eventarray[7];
		      float eventmass = -999.99;
		      float d2 = 0;
		      //float kinewidth[6];
                      float kinewidth[7];

		      eventarray[0] = kinevar1;
		      eventarray[1] = kinevar2;
		      eventarray[2] = kinevar3;
		      eventarray[3] = kinevar4;
		      eventarray[4] = kinevar5;
                      eventarray[5] = kinevar6;
                      eventarray[6] = kinevar7;

		      if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
			kinewidth[0] = ebin_width; // ebin width in MeV
			kinewidth[1] = 2;          // cosrealtheta range
			kinewidth[2] = 2*M_PI;     // phi range
			kinewidth[3] = 2;          // coscmstheta_pro range
		      }
                      else if (omega_flag == 1 && eta_flag == 0 && phi_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
                        kinewidth[0] = 2;          // cos(theta) of (pi+ pi- pi0) in center-of-mass
                        kinewidth[1] = 2;          // cos(theta*) - helicity angle
                        kinewidth[2] = 2*M_PI;     // phi* - helicity angle
                        kinewidth[3] = 1;          // lambda - 3pi decay parameter
                        kinewidth[4] = 2*M_PI;     // lab phi of (pi+ pi- pi0)
                        kinewidth[5] = ebin_width; // EBin width in MeV
		      }
		      else if (phi_flag == 1 && omega_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {

                      }
		      else if (eta_flag == 1 && omega_flag == 0 && phi_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
			kinewidth[0] = 2;          // cos(theta) of (pi+ pi- pi0) in center-of-mass
			kinewidth[1] = 2;          // cos(theta*) - helicity angle
			kinewidth[2] = 2*M_PI;     // phi* - helicity angle
			kinewidth[3] = 1;          // lambda - 3pi decay parameter
                        kinewidth[4] = 2*M_PI;     // lab phi of (pi+ pi- pi0)
			kinewidth[5] = ebin_width; // EBin width in MeV
                      }
		      else if (Kaon_flag == 1 && omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Sigma_flag == 0) {
			kinewidth[0] = ebin_width; // EBin width in MeV
			kinewidth[1] = 2;          // cos(theta*) of pi(+)
			kinewidth[2] = 2*M_PI;     // phi*
                        kinewidth[3] = 2;          // cos(theta) of (pi+ pi-) in center-of-mass
                        kinewidth[4] = 2*M_PI;     // lab phi of (pi+ pi-)
                        kinewidth[5] = 2;          // cos(angle) between proton and pi(0)
		      }
                      else if (Kaon_flag == 0 && omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Sigma_flag == 1) {
                        kinewidth[0] = ebin_width; // EBin width in MeV
			kinewidth[1] = 2;          // cos(theta*) of pi(+)
                        kinewidth[2] = 2*M_PI;     // phi*
                        kinewidth[3] = 2;          // cos(theta) of (pi+ pi-) in center-of-mass
                        kinewidth[4] = 2*M_PI;     // lab phi of (pi+ pi-)
                        kinewidth[5] = 2;          // cos(angle) between proton and pi(0)
                        kinewidth[6] = 2;
                      }
		      eventmass = fitvar;
			  
		      // Calculate the distance between the seedevent and the ith event
		      /*for (int kv=0; kv<6; kv++) {
			d2 += TMath::Power( ( seedeventarray[kv] - eventarray[kv] ) / kinewidth[kv], 2.0);*/
                       for (int kv=0; kv<7; kv++) {
			d2 += TMath::Power( ( seedeventarray[kv] - eventarray[kv] ) / kinewidth[kv], 2.0);
		      }

		      // Create arrays of distance and mass with numNN_def # of events
		      if (j<numNN_def) {
			             
			// Fill array with distance and mass
			d2array[j] = d2;
		        massarray[j] = eventmass;
			eventnumarray[j] = i; // Don't use event_num here since two same event no.s can be from diff. runs.
			j++;
		      }

		      if (j == numNN_def) { // Find distance max in the d2array
			if (change == 1) {
			  index = 0;
			  max = d2array[0];

			  for(int k=1; k<numNN_def; k++) {
			    if (d2array[k] > max) {
			      index = k;
			      max = d2array[k]; 
			    }
			  }
			}
			iteration++;
			     
			if (d2 < max && iteration>1) {
			  //cerr << "d2 of ith event is < max. d2 array changed." << endl << endl;
			  d2array[index] = d2;
			  massarray[index] = eventmass;
			  eventnumarray[index] = i;
			  change = 1;
			}
			else {
			  //cerr << "d2 array not changed." << endl << endl;
			  change = 0;
			}
		      }
		    }
		  }
		}
	      }

	      // end of finding nearest neighbors for the seed event
	      //cerr << "the NNarray to seed event # " << eventnum_nt << " is: " << endl;
	      if (finaltree_save == 1) {
		for (int a=0; a<numNN_def; a++) {
                   NNarray_nt[a]=eventnumarray[a];
		   //cerr << "NNarray[" << a  << "] = " << NNarray_nt[a] << endl;
		}
	      }

	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
     	      //
	      //  Initialize Q-value fitter
	      //
	      Qvalue_Fitter qfitter;
	 
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Data Variable
	      //
	      RooRealVar wmass("wmass","wmass",beginfit,endfit,"MeV");

	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Gaussian function (signal) parameters and pdf
	      //
	      double Mean = 0.0;
	      double Mean_Min = 0.0;
	      double Mean_Max = 0.0;

              double Mean_DS = 0.0; // DS = double signal
              double Mean_Min_DS = 0.0;
              double Mean_Max_DS = 0.0;

	      double Sigma = 0.0;
	      double Sigma_Min = 0.0;
	      double Sigma_Max = 0.0;

              double Sigma_I = 0.0;
              double Sigma_I_Min = 0.0;
              double Sigma_I_Max = 0.0;

              double Sigma_II = 0.0;
              double Sigma_II_Min = 0.0;
              double Sigma_II_Max = 0.0;

	      double c_1 = 0.0;
	      double c_1_Min = 0;
	      double c_1_Max = 0;

	      double c_2 = 0.0;
	      double c_2_Min = 0.0;
              double c_2_Max = 0.0;

	      double c_3 = 0.0;
	      double c_3_Min = 0.0;
              double c_3_Max = 0.0;

	      if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
		if (sit_limit == 1) {
		  Mean = 140;
		  Mean_Min = 120.;
		  Mean_Max = 160.;

		  c_2 = -0.1;
		  c_2_Min = -0.5;
		  c_2_Max = 0.5;
		}
		else if (sit_limit == 2) {
		  Mean = 140;
		  Mean_Min = 120.;
	          Mean_Max = 160.;

		  c_2 = -0.2;
		  c_2_Min = -1.0;
		  c_2_Max = 1.0;
		}
		else if (sit_limit == 4) {
		  Mean = 140.;
		  Mean_Min = 120.;
		  Mean_Max = 160.;

		  c_2 = -0.1;
		  c_2_Min = -0.5;
		  c_2_Max = 0.5;
		}

		Sigma = 28.5;
		Sigma_Min = 2;
		Sigma_Max = 200;

		c_1 = 0.2;
		c_1_Min = -0.5;
		c_1_Max = 1.0;

		c_3 = 0.05;
	        c_3_Min = -0.01;
	        c_3_Max = 0.1;
	      }
		
	      else if (omega_flag == 1 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
		Mean = 782.65;
                Mean_Min = 782.65;
                Mean_Max = 782.65;

                Mean_DS = 782.65;
                //Mean_Min_DS = 775.0;
		//Mean_Max_DS = 805.0;
		Mean_Min_DS = 740.0;
		Mean_Max_DS = 865.0;

                //Sigma = 8.0;
                //Sigma_Min = 3.0;
                //Sigma_Max = 23.0;

                Sigma_I = 8.0;
                //Sigma_I_Min = 2.0;
                //Sigma_I_Max = 18.0;
                Sigma_I_Min = 0.0;
                Sigma_I_Max = 30.0;

                Sigma_II = 15.0;
                //Sigma_II_Min = 14.0;
                //Sigma_II_Max = 45.0;
                Sigma_II_Min = 0.0;
		Sigma_II_Max = 55.0;

                c_1 = 0.5;
                c_1_Min = 0.0;
                c_1_Max = 1.8;

		c_2 = 0.1;
                c_2_Min = -1.2;
                c_2_Max = 1.2;

		//c_3 = -0.005;
	        //c_3_Min = -1.0;
	        //c_3_Max = 0.8;
	      }

	      else if (omega_flag == 0 && Kaon_flag == 1 && eta_flag == 0 && phi_flag == 0 && Sigma_flag == 0) {
		Mean = 497.614;
		Mean_Min = 497.614;
		Mean_Max = 497.614;

                //Mean_DS = 497.614;
                //Mean_Min_DS = 490.0;
		//Mean_Max_DS = 510.0;

		Sigma = 4.5;
	        Sigma_Min = 2.0;
                Sigma_Max = 8.0;

                //Sigma_II = 10.0;
                //Sigma_II_Min = 6.0;
		//Sigma_II_Max = 28.0;

		c_1 = 0.1;
                c_1_Min = -1.5;
                c_1_Max = 1.5;

                c_2 = 0.1;
                c_2_Min = -1.5;
                c_2_Max = 1.5;

                //c_3 = 0.5;
                //c_3_Min = -1.5;
                //c_3_Max = 1.5;
	      }

              else if (omega_flag == 0 && Kaon_flag == 0 && eta_flag == 0 && phi_flag == 0 && Sigma_flag == 1) {
                Mean = 1189.37;
                Mean_Min = 1189.37;
                Mean_Max = 1189.37;

                Sigma = 4.5;
                Sigma_Min = 0.0;
                Sigma_Max = 9.0;

                c_1 = 0.1;
                c_1_Min = -1.5;
                c_1_Max = 1.5;

                c_2 = 0.1;
                c_2_Min = -1.5;
                c_2_Max = 1.5;

                //c_3 = 0.5;
                //c_3_Min = -1.5;
                //c_3_Max = 1.5;
              }

              else if (omega_flag == 0 && Kaon_flag == 0 && eta_flag == 1 && phi_flag == 0 && Sigma_flag == 0) {
                Mean = 548.4; // Peak appears shifted by 0.5 MeV in g12 data; M = 547.862 (PDG 2014)  
		Mean_Min = 548.4;
                Mean_Max = 548.4;

                Mean_DS = 548.4; // Peak appears shifted by 0.5 MeV in g12 data
                Mean_Min_DS = 540.0;
		Mean_Max_DS = 560.0;

                //Sigma = 5.0;
                //Sigma_Min = 3.0;
                //Sigma_Max = 10.0;

                Sigma_I = 5.0;
                Sigma_I_Min = 1.5;
                Sigma_I_Max = 9.5;

                Sigma_II = 10.0;
                Sigma_II_Min = 6.0;
                Sigma_II_Max = 28.0;

		c_1 = 0.8;
                c_1_Min = 0.0;
                c_1_Max = 1.6;

                c_2 = 0.19;
                c_2_Min = -0.6;
                c_2_Max = 1.5;

                //c_3 = -0.005;
                //c_3_Min = -0.5;
                //c_3_Max = 0.5;
              }


	      ///////////////////////
	      // Signal Definition //
	      ///////////////////////

              // Fix the (first) mass of the omega, eta, Kaon, and Sigma
	      v_mean = new RooRealVar("v_mean","v_mean", Mean, Mean_Min, Mean_Max);
              v_mean_DS = new RooRealVar("v_mean_DS","v_mean_DS", Mean_DS, Mean_Min_DS, Mean_Max_DS);
	      if (omega_flag == 1 || eta_flag == 1 || phi_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1) 
		v_mean->setConstant(kTRUE);
              else if (omega_flag == 0 && eta_flag == 0 && phi_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0)
                v_mean->setConstant(kTRUE);

              // Fit the (second) mass for eta and omega (double signal)
              if (eta_flag == 1 || omega_flag == 1) {
                qfitter.SetSignalVariable(v_mean_DS);
	      }

	      v_gamma = new RooRealVar("v_gamma","v_gamma", 8.49);
              if (omega_flag == 1)
		v_gamma->setConstant(kTRUE);
	       
	      v_sigma = new RooRealVar("v_sigma","v_sigma", Sigma, Sigma_Min, Sigma_Max);
              if (omega_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && phi_flag == 0 && Sigma_flag == 0)
	        qfitter.SetSignalVariable(v_sigma);
              else if (Kaon_flag == 1 || Sigma_flag == 1) 
		qfitter.SetSignalVariable(v_sigma);

              v_sigma_I = new RooRealVar("v_sigma_I","v_sigma_II", Sigma_I, Sigma_I_Min, Sigma_I_Max);
	      if (eta_flag == 1 || omega_flag == 1) {
		qfitter.SetSignalVariable(v_sigma_I);
	      }

              v_sigma_II = new RooRealVar("v_sigma_II","v_sigma_II", Sigma_II, Sigma_II_Min, Sigma_II_Max);
	      if (eta_flag == 1 || omega_flag == 1) {
		qfitter.SetSignalVariable(v_sigma_II);
	      }

              sig1frac = new RooRealVar("sig1frac","fraction of component 1 in signal",0.8,0.0,1.0);
	      if (eta_flag == 1 || omega_flag == 1) {
		qfitter.SetSignalVariable(sig1frac);
	      }

              // Define signal pdf's
	      signalpdf = new RooGaussian("signalpdf","Signal (Gaussian)",wmass,*v_mean,*v_sigma);

	      signalpdf_Voigtian_I = new RooVoigtian("signalpdf_Voigtian","Signal (Voigtian)",wmass,*v_mean,*v_gamma,*v_sigma_I);
              signalpdf_Voigtian_II = new RooVoigtian("signalpdf_Voigtian_II","Signal (Voigtian II)",wmass,*v_mean_DS,*v_gamma,*v_sigma_II);

	      //signalpdf_eta = new RooBifurGauss("signalpdf_eta","Signal (Bifur Gaussian)",wmass,*v_mean,*v_sigmaLeft,*v_sigmaRight);
	      signalpdf_Gaussian_I = new RooGaussian("signalpdf_I","Signal I (Gaussian)",wmass,*v_mean,*v_sigma_I);
	      signalpdf_Gaussian_II = new RooGaussian("signalpdf_II","Signal II (Gaussian)",wmass,*v_mean_DS,*v_sigma_II);

	      signalpdf_eta = new RooAddPdf("signalpdf_eta","Signal (Double Gaussian)",RooArgList(*signalpdf_Gaussian_I,*signalpdf_Gaussian_II),*sig1frac);
              signalpdf_omega = new RooAddPdf("signalpdf_omega","Signal (Double Voigtian)",RooArgList(*signalpdf_Voigtian_I,*signalpdf_Voigtian_II),*sig1frac);

	      if (eta_flag == 1) qfitter.SetSignalPdf(signalpdf_eta); // Double Gaussian
	      else if (omega_flag == 1) qfitter.SetSignalPdf(signalpdf_omega); // Double Voigtian
              else qfitter.SetSignalPdf(signalpdf); // Gaussian otherwise


	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //  
	      //
	      //  Chebychev function (background) parameters and pdf
	      //
	      c1 = new RooRealVar("c1","c1", c_1, c_1_Min, c_1_Max); 
	      if (argusbg_flag == 0) qfitter.SetBackgroundVariable(c1);
 
	      c2 = new RooRealVar("c2","c2", c_2, c_2_Min, c_2_Max);
	      if (argusbg_flag == 0 && (omega_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1)) qfitter.SetBackgroundVariable(c2);
	     
	      c3 = new RooRealVar("c3","c3", c_3, c_3_Min, c_3_Max);
              //if (omega_flag == 1) qfitter.SetBackgroundVariable(c3);

	      m0 = new RooRealVar("m0","m0", 820.0, 790.0, 950.0);
              if (argusbg_flag == 1) qfitter.SetBackgroundVariable(m0);

              c = new RooRealVar("c","c", -1.0, -10.0, 0.2);
              if (argusbg_flag == 1) qfitter.SetBackgroundVariable(c); 
              //c->setConstant(kTRUE);
	      //qfitter.SetBackgroundVariable(c);

              if (omega_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1) {
                RooArgList coefficientList(*c1, *c2);
                backgroundpdf = new RooChebychev("backgroundpdf","background pdf", wmass, coefficientList);
              }
	      else {
                RooArgList coefficientList(*c1);
                backgroundpdf = new RooChebychev("backgroundpdf","background pdf", wmass, coefficientList);
	      }
	      backgroundpdf_prime = new RooArgusBG("backgroundpdf_prime","background pdf (Argus)", wmass, *m0, *c);
	      if (argusbg_flag == 0) qfitter.SetBackgroundPdf(backgroundpdf);
	      else if (argusbg_flag == 1) qfitter.SetBackgroundPdf(backgroundpdf_prime);

	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Set Fit range
	      //
      	      qfitter.SetFitRange(beginfit, endfit);
	  
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Define TCanvas to hold visual plot of the fits
	      //     
	      histoname = "EventFit_";
	      histoname += event;

	      TCanvas *fitcanvas;
              fitcanvas = new TCanvas(histoname,histoname,hist_Binning,10,700,500);

              //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Define tree to hold fit variables for a unbinned maximum likelihood fit.
	      //
	      TTree *fitvartree = new TTree("pionmass", "pion mass");
	      fitvartree->Branch("wmass",&mass,"pionmass/D");

	      if (event%50 == 0) cerr << " QVALUEINFO: eventnumber = " << event << endl;

   	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Fill ttree and plot nearest neighbors
	      //
              //if (event>1004 && event<1007) cerr << "list of NN masses to for dummy event #" << event << endl;
	      for (int i=0; i<numNN_def; i++) {
		 //pionhisto1->Fill(massarray[i]);
	         mass = massarray[i];
	         //if (event>1004 && event<1007) cerr << mass << "   ";
		 fitvartree->Fill();
	      }
	  
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Define the data set for roofit
	      //	  
	      RooDataSet treedata("treedata","treedata",fitvartree,wmass);
 
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Tell the fitter about the data
	      //
	      qfitter.SetDataSet(&wmass, &treedata, numNN_def);
	  
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Tell the fitter where to define the qvalue at
	      //
	      //cerr << "Do event-based fit for seedevent with seedmass = " << seedmass << endl;	  
	      qfitter.SetSeatEventValue(seedmass);

	      // Trace memory
	      //RooTrace::active(true);
	      //RooTrace::verbose(true);

	      if (event == 50 || event == 99 || event%3000 == 0) plotflag=1;
	      else plotflag=0;

	      //RooTrace::active(1);
	      //qfitter.Fit(&chi2, fitcanvas, &SigFrac, &NDF, hist_Binning, &SigFlag, &signalFitmean, &signalFitsigma, &C1, &Qvalue_error_total, plotflag);
	      qfitter.Fit(&chi2, fitcanvas, &SigFrac, &NDF, hist_Binning, meson_flag, &par1, &par2, &par3, &Qvalue_error_total, plotflag);

	      par1 = v_mean->getVal();
	      if (omega_flag == 1 && eta_flag == 0 && phi_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
                par2 = v_sigma_I->getVal();
                par3 = v_sigma_II->getVal();
		par4 = v_gamma->getVal();
		par7 = v_mean_DS->getVal();
                par8 = sig1frac->getVal();
                if (argusbg_flag == 0) {
		  par5 = c1->getVal();
		  par6 = c2->getVal();
		} 
                else if (argusbg_flag == 1) {
                  par5 = m0->getVal();
                  par6 = c->getVal();
		}
	      }
              else if (phi_flag == 0 && omega_flag == 0 && Kaon_flag == 1 && eta_flag == 0 && Sigma_flag == 0) {
		par2 = v_sigma->getVal();
                par3 = c1->getVal();
                par4 = c2->getVal();
	      }
              else if (phi_flag == 0 && omega_flag == 0 && Kaon_flag == 0 && eta_flag == 0 && Sigma_flag == 1) {
                par2 = v_sigma->getVal();
                par3 = c1->getVal();
                par4 = c2->getVal();
              }
              else if (phi_flag == 0 && omega_flag == 0 && Kaon_flag == 0 && eta_flag == 1 && Sigma_flag == 0) {
                par2 = v_sigma_I->getVal();
                par3 = v_sigma_II->getVal();
                par4 = c1->getVal();
		par5 = c2->getVal();
                par6 = v_mean_DS->getVal();
                par7 = sig1frac->getVal();
              }
	      else if (omega_flag == 0 && phi_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) {
                par2 = v_sigma->getVal();
		par3 = c1->getVal();
	      }

	      // Deleting all the pointers -
	      v_mean->Delete();
              v_mean_DS->Delete();
	      v_sigma->Delete();
              v_sigma_I->Delete();
              v_sigma_II->Delete();
	      v_gamma->Delete();
              sig1frac->Delete();

	      c1->Delete();
	      c2->Delete();
	      c3->Delete();

	      m0->Delete();
              c->Delete();

	      signalpdf->Delete();
	      signalpdf_omega->Delete();
              signalpdf_Voigtian_I->Delete();
              signalpdf_Voigtian_II->Delete();
              signalpdf_Gaussian_I->Delete();
              signalpdf_Gaussian_II->Delete();
              signalpdf_eta->Delete();
	      backgroundpdf->Delete();
	      backgroundpdf_prime->Delete();

	      fitvartree->Delete();

	      //RooTrace::verbose(kTRUE);
	      //cerr << endl << endl << "This dump is right after the qfitter.Fit command and after deleting v_mean, v_sigma and pdfs" << endl << endl;
	      //RooTrace::dump();

	      if (event >= 0 && event < 20) {
		cout << endl;
	        cout << " Signal Fraction:   " << SigFrac << endl;
	        cout << " numEvnt:           " << numNN_def << endl;
	        cout << " chi2:              " << chi2 << endl;
	        cout << " NDF:               " << NDF << endl;
	        cout << " signal-Fit-Mean:   " << par1 << endl;
		//cout << " c1 (Chebychev):    " << par3 << endl;
	      }

	      qvalue[0] = qfitter.GetQvalue(0,&qvalue[1]);

	      Qval_nt = qvalue[0]; // finaltree branch
	      Qfit_error_nt = qvalue[1]; // finaltree branch
		
	      if ((event >= 0 && event < 20) || (event > 1003 && event < 1007)) {
		cerr << endl;
		cerr << "The qvalues and errors for seed event #" << eventnum_nt << " - " << endl;
		cerr << "qvalue = " << qvalue[0] << endl;
		cerr << "fit error = " << qvalue[1] << endl << endl;
	      }

	      // Filling histograms		
	      if (omega_flag == 0 && eta_flag == 0 && phi_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) FitC1histo -> Fill(par3);

	      if (!eta_flag == 1 && !omega_flag == 1) FitSigmahisto -> Fill(par2);
		   
              if (phi_flag == 0 && eta_flag == 0 && omega_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) totalhisto[0] -> Fill(seedmass);
	      else if (omega_flag == 1) totalhisto[1] -> Fill(seedmass);
              else if (Kaon_flag == 1) totalhisto[2] -> Fill(seedmass);
              else if (eta_flag == 1) totalhisto[3] -> Fill(seedmass);
              else if (phi_flag == 1) totalhisto[4] -> Fill(seedmass);
              else if (Sigma_flag == 1) totalhisto[5] -> Fill(seedmass);

              if (phi_flag == 0 && eta_flag == 0 && omega_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) FitMeanhisto[0] -> Fill(par1);
              else if (omega_flag == 1) FitMeanhisto[1] -> Fill(par1);
	      else if (Kaon_flag == 1) FitMeanhisto[2] -> Fill(par1);
	      else if (eta_flag == 1) FitMeanhisto[3] -> Fill(par1);
	      else if (phi_flag == 1) FitMeanhisto[4] -> Fill(par1);
              else if (Sigma_flag == 1) FitMeanhisto[5] -> Fill(par1);

              if (phi_flag == 0 && eta_flag == 0 && omega_flag == 0 && Kaon_flag == 0 && Sigma_flag == 0) Sigma_vs_Mean_histo[0] -> Fill(par1,par2);
              else if (omega_flag == 1) Sigma_vs_Mean_histo[1] -> Fill(par1,par2);
              else if (Kaon_flag == 1) Sigma_vs_Mean_histo[2] -> Fill(par1,par2);
              else if (eta_flag == 1) Sigma_vs_Mean_histo[3] -> Fill(par1,par2);
              else if (phi_flag == 1) Sigma_vs_Mean_histo[4] -> Fill(par1,par2);
              else if (Sigma_flag == 1) Sigma_vs_Mean_histo[5] -> Fill(par1,par2);

              if (omega_flag == 1) {
                FitSigma_I_histo -> Fill(par2);
                FitSigma_II_histo -> Fill(par3);
		FitGammahisto -> Fill(par4);
                Sigma_I_vs_Sigma_II -> Fill(par2,par3);
 
                if (argusbg_flag == 0) {
		  FitC1histo -> Fill(par5);
		  FitC2histo -> Fill(par6);
                }
                else if (argusbg_flag == 1) {
                  FitM0histo -> Fill(par5);
		  FitChisto -> Fill(par6);
		}

                FitMeanDS -> Fill(par7);
                FitSig1frac -> Fill(par8);
	      }

	      if (Kaon_flag == 1) {
		FitC1histo -> Fill(par3);
                FitC2histo -> Fill(par4);
	      }

              if (Sigma_flag == 1) {
                FitC1histo -> Fill(par3);
                FitC2histo -> Fill(par4);
              }

	      if (eta_flag == 1) {
                FitSigma_I_histo -> Fill(par2);
                FitSigma_II_histo -> Fill(par3);
		Sigma_I_vs_Sigma_II -> Fill(par2,par3);

		FitC1histo -> Fill(par4);
		FitC2histo -> Fill(par5);

                FitMeanDS -> Fill(par6);
                FitSig1frac -> Fill(par7);
	      }

	      chi2histo1 -> Fill(chi2);

	      qvaluehisto -> Fill(qvalue[0]);
	      qvalueerrorhisto -> Fill(qvalue[1]);
	      qvalerror_vs_qval_histo -> Fill(qvalue[0], qvalue[1]);

              if (qvalue[0] >= 0.0) {
		if (omega_flag == 1) {
		  backgroundhisto[1] -> Fill(seedmass, 1.0 - qvalue[0]);
		  mmhisto[1] -> Fill(seedmass, qvalue[0]);
		  Chi2_vs_Sigma_histo[1] -> Fill(par2, chi2);
                  Chi2_vs_Mean_histo[1] -> Fill(par1, chi2);

                  signal_lambdahisto -> Fill(seedeventarray[3], qvalue[0]);
                  bckgrd_lambdahisto -> Fill(seedeventarray[3], 1.0 - qvalue[0]);
                  total_lambdahisto -> Fill(seedeventarray[3]);

                  signal_costhetaomegahisto -> Fill(seedeventarray[0], qvalue[0]);
                  bckgrd_costhetaomegahisto -> Fill(seedeventarray[0], 1.0 - qvalue[0]);
                  total_costhetaomegahisto -> Fill(seedeventarray[0]);

                  omega_sig_mom_histo -> Fill(omega_p, qvalue[0]);
                  omega_bckgrd_mom_histo -> Fill(omega_p, 1.0 - qvalue[0]);
		}
                else if (Kaon_flag == 1) {
		  backgroundhisto[2] -> Fill(seedmass, 1.0 - qvalue[0]);
                  mmhisto[2] -> Fill(seedmass, qvalue[0]);
                  Chi2_vs_Sigma_histo[2] -> Fill(par2, chi2);
                  Chi2_vs_Mean_histo[2] -> Fill(par1, chi2);
		}
                else if (eta_flag == 1) {
		  backgroundhisto[3] -> Fill(seedmass, 1.0 - qvalue[0]);
                  mmhisto[3] -> Fill(seedmass, qvalue[0]);
                  Chi2_vs_Sigma_histo[3] -> Fill(par2, chi2);
                  Chi2_vs_Mean_histo[3] -> Fill(par1, chi2);

                  signal_lambdahisto -> Fill(seedeventarray[3], qvalue[0]);
                  bckgrd_lambdahisto -> Fill(seedeventarray[3], 1.0 - qvalue[0]);
                  total_lambdahisto -> Fill(seedeventarray[3]);
		}
                else if (phi_flag == 1) {
		  backgroundhisto[4] -> Fill(seedmass, 1.0 - qvalue[0]);
                  mmhisto[4] -> Fill(seedmass, qvalue[0]);
		  Chi2_vs_Sigma_histo[4] -> Fill(par2, chi2);
                  Chi2_vs_Mean_histo[4] -> Fill(par1, chi2);
		}
                else if (Sigma_flag == 1) {
                  backgroundhisto[5] -> Fill(seedmass, 1.0 - qvalue[0]);
                  mmhisto[5] -> Fill(seedmass, qvalue[0]);
                  Chi2_vs_Sigma_histo[5] -> Fill(par2, chi2);
                  Chi2_vs_Mean_histo[5] -> Fill(par1, chi2);
                }

                if (omega_flag == 0 && eta_flag == 0 && Kaon_flag == 0 && phi_flag == 0 && Sigma_flag == 0) {
		  if (seedmass > 0 && seedmass < 150.0) {
		    background_labphi_pip_histo[0] -> Fill(labphi_pip, 1.-qvalue[0]);
		    background_labphi_pim_histo[0] -> Fill(labphi_pim, 1.-qvalue[0]);
		  }

                  if (seedmass > 350.0 && seedmass < 500) {
                    background_labphi_pip_MM350up_histo[0] -> Fill(labphi_pip, 1.-qvalue[0]);
                    background_labphi_pim_MM350up_histo[0] -> Fill(labphi_pim, 1.-qvalue[0]);
		  }
		}

		chi2histo -> Fill(chi2);
		SFhisto -> Fill(SigFrac);
	      }
                                                
	      //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	      //
	      //  Write canvases 
	      //		
	      /* if (qvalue[0] >=.8) plotflag=1;
	      else plotflag=0;*/

	      if (event==503 || event==6003 || event==12003 || event==20003 || event == (counter_But-1) || event == NumOfEvnt || plotflag > 0) {
		writeRoot(outputfile, fitcanvas, event, counter_But, NumOfEvnt, seedmass, NewEBin_val, period_limit);
	      }
	      //if (TestRun==1 && event == NumOfEvnt)writeRoot(outputfile, fitcanvas, event, counter_But, NumOfEvnt, seedmass);	
              delete fitcanvas;
		 
	      if (finaltree_save == 1) {
		// Filling final tree       
		finaltree->Fill();
	      }

	      // Writing the output Qvalue textfile
              if (omega_flag == 1 || phi_flag == 1 || eta_flag == 1 || Kaon_flag == 1 || Sigma_flag == 1) topology_nt = 5;
	      if (Qval_nt > 0.00000){
		if (deg_bpol_nt > 0.00) OutQvalue << run_number_nt << "  " << eventnum_nt << "  " << topology_nt << "  " << deg_bpol_nt << "  " << endl;
		else OutQvalue << run_number_nt << "  " << eventnum_nt << "  " << topology_nt << "  " << bpol_nt << "  " << endl;
		OutQvalue << realphi_nt << "  " << cosrealtheta_nt << "  " << coscmstheta_pro_nt << "  " << endl;
		OutQvalue << pho_energy_nt << "  " << endl;
		OutQvalue << pro_px_nt << "  " << pro_py_nt << "  " << pro_pz_nt << "  " << pro_E_nt << "  " << endl;
		OutQvalue << pip_px_nt << "  " << pip_py_nt << "  " << pip_pz_nt << "  " << pip_E_nt << "  " << endl;
		OutQvalue << pim_px_nt << "  " << pim_py_nt << "  " << pim_pz_nt << "  " << pim_E_nt << "  " << endl;
		OutQvalue << cl_2pi_nt << "  " << cl_3pi_nt << endl;
		OutQvalue << z_vertex_nt << "  " << endl;
		OutQvalue << pho_energy_kinf_nt << "  " << endl;
		OutQvalue << pro_px_kinf_nt << "  " << pro_py_kinf_nt << "  " << pro_pz_kinf_nt << "  " << pro_E_kinf_nt << "  " << endl;
		OutQvalue << pip_px_kinf_nt << "  " << pip_py_kinf_nt << "  " << pip_pz_kinf_nt << "  " << pip_E_kinf_nt << "  " << endl;
		OutQvalue << pim_px_kinf_nt << "  " << pim_py_kinf_nt << "  " << pim_pz_kinf_nt << "  " << pim_E_kinf_nt << "  " << endl;
                OutQvalue << Qval_nt << "  " << Qfit_error_nt << endl;
	      }

	      if (TestRun) {
	        if (event == NumOfEvnt) break;
	      }
	    }
	  }
	    
	  //
	  //  End of Qvalue fitting loop over all events
	  //
	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
	  vartree->Delete(); // Remove vartree  

	  if (finaltree_save == 0) finaltree->Delete();
	  //else finaltree->Print();

	  //  //  //  //  //  //  //  //  //  //  //  //  //  //  // 
	  // 
	  // Finding Qvalue error 
	  //
	  //
	  /*
	  cerr<<endl<<endl<<"calculating net error"<<endl;

	  Int_t Nentries = (Int_t)finaltree->GetEntries();
	  cerr<<endl<<endl<<"no. of entries in final tree = "<<Nentries<<endl;

	  qerror=0.00;

	  for (int i=0; i<Nentries; i++)
	    { 
	      finaltree->GetEntry(i); 
	      int event_num1=eventnum_nt;
	      if(i<2)cerr<<"event_num1 = " <<event_num1<<endl;

	      float Qfit_error1= Qfit_error_nt;
	      int NNarray1[numNN_def];

	      //if(i<2)cerr<<"NN array for event #"<<event_num1<<" : "<<endl<<endl;

	      for (int j=0; j<numNN_def; j++)
		{
                   NNarray1[j] = NNarray_nt[j];
		   //if(i<2)cerr<<"NNarray["<<j<<"] = "<< NNarray1[j]<<endl;
		}
	      for (int k=0; k<Nentries; k++)
		{
		  finaltree->GetEntry(k);
		  int correlation = 0;
		  int event_num2=eventnum_nt;

		  //if(event_num2 >= event_num1)
		    //{
                    // double counting is needed ! don't avoid it !

		      if(k<3)cerr<<"event_num2 = "<<event_num2<<endl;
		      float Qfit_error2= Qfit_error_nt;
		      int NNarray2[numNN_def];
		      //if(k<3)cerr<<"NN array for event #"<<event_num2<<" : "<<endl<<endl;

		      for (int j=0; j<numNN_def; j++)
			{
                           NNarray2[j] = NNarray_nt[j];
			   //if(k<3)cerr<<"NNarray["<<j<<"] = "<< NNarray2[j]<<endl;
			}

		      // find the # of common nearest neighbors between event_num1 and event_num2
		      for (int a=0; a<numNN_def; a++)
			{
			  for (int b=0; b<numNN_def; b++)
			    {
			      if(NNarray1[a]==NNarray2[b] && NNarray1[a]!=-1)correlation++;
			    }
			}
		      //cerr<<"# common events to event_num1 (#"<<event_num1<<") and event_num2(#"<<event_num2<<") = "<<correlation<<endl;
		      //cerr<<"Qfiterror1 = "<<Qfit_error1<<endl<<"Qfit error2 = "<<Qfit_error2<<endl;

		      qerror = qerror + Qfit_error1*(correlation/numNN_def)*Qfit_error2;
		 
                  //}
		  //else continue;
		}
	    }
          qerror = TMath::Sqrt(qerror);

	  cerr<<"the net error in the determination of Qvalues for each event = "<<qerror<<endl;
	    */
     
	  cerr << endl << endl << "Total no. of events for which Q factors have been extracted = " << event << endl;

	} //for (int newpbin_limit=PBin_Min; newpbin_limit<PBin_Max+1; newpbin_limit++) {
      } //for (int bpol_limit=1; bpol_limit<3; bpol_limit++) {       
    } //for (int wbin_limit=3; wbin_limit<13; wbin_limit++) {
  }//for (int period_limit=7; period_limit<8; period_limit++) {
    
  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
  //
  //  Close output
  //
  outputfile->Write();
  outputfile->Close();
  
  cerr << endl << " Bye Bye " << endl;
  //  //  //  //  //  //  //  //  //  //  //  //  //  //  //
  //
  //  The end time to check the running time
  //     
  printf(" Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
}
  
TLorentzVector lorentzBoost(TLorentzVector boostVector,TLorentzVector vector)
{
  TLorentzVector resultVector ( vector );
  TVector3 b = boostVector.BoostVector();
  //cerr<<"this is lorentzboost script"<<endl<<"3vector b = ("<<b.X()<<","<<b.Y()<<","<<b.Z()<<")"<<endl;                                                                           
  resultVector.Boost ( boostVector.BoostVector() );
  return resultVector;
}