////////// FIRST STEP, DEFINING THE METRIC
//THIS IS FOR KS CHANNEL WHICH I USE 7 VARIABLE FOR METRIC
 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);
		  }

                
		  
		  // 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;

//then store the kinevar information, into array or ntuples

///////// SECOND STEP : TWO LOOP, FOR EVERY EVENT WE LOOK INTO ANOTHER EVENT IN THE SAME ENERGY BIN TO FIND THE NEIGHBOR (300-1000) DEPEND ON STATISTIC
////// REMEMBER THAT THERE IS TWO LOOP HERE

for (int counter=0; counter<nentries; counter++) {
	    vartree->GetEntry(counter);  // Get row for one event
	    

	    if (massflag != 1) continue; // Q-factor won't be calculated for this event, to make sure that event within the range of desired mass
	    else {

              
	      

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

	      float d2array[numNN_def]; //numNN_def is the number of neighbor, depend on statistic
	      float massarray[numNN_def];
	      int eventnumarray[numNN_def];

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

	      seedeventarray[0] = kinevar1; //put the kinematic variable into our event
	      seedeventarray[1] = kinevar2;
	      seedeventarray[2] = kinevar3;
	      seedeventarray[3] = kinevar4;
	      seedeventarray[4] = kinevar5;
              seedeventarray[5] = kinevar6;
              seedeventarray[6] = kinevar7;

	      seedmass = fitvar; //in this case, this is the mass of proton pi0 or Mppi0 since I do Qvalue on sigma peak

	      // // // // // // // // // // // // // // // // // // // 
	      //
	      // 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++) { //this is the second loop, where for each event we loop to the whole event to find the nearest neighbor

		  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[7];
		      float eventmass = -999.99;
		      float d2 = 0;
		      
                      float kinewidth[7];

		      eventarray[0] = kinevar1; //put the kinematic variable into our event
		      eventarray[1] = kinevar2;
		      eventarray[2] = kinevar3;
		      eventarray[3] = kinevar4;
		      eventarray[4] = kinevar5;
                      eventarray[5] = kinevar6;
                      eventarray[6] = kinevar7;

		      //this is the maximum distance for each variable for normalization
                        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<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;
		}
	      }

////////////// QValue calculation
//THIRD STEP : DEFINE SIGNAL FUNCTION DEFINE BACKGROUND FUNCTION AND DO THE FIT
//THE FIT USING MODULE QFITTER. I WILL LOOK INTO THIS FITTER MODULE AND SEND IT LATER
// SORRY THAT IT STILL CONTAIN ALOT OF FLAG. THIS IS PART OF GIANT CODE THAT CONTAIN MANY REACTION FROM MANY DATASET



/////
//  //  //  //  //  //  //  //  //  //  //  //  //  //  //
     	      //
	      //  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;
	      }
///////////