source: huonglan/ANA2011/chi2calculation.C

#include <iostream>
#include <math.h>
#include <utility>

#include <TLorentzVector.h>
#include <TVector3.h>

using namespace std;

#include "TFitter.h"
#include "TMath.h"

#include "chi2class.h"

#define PRINT 0
#define mjet0 false
#define ksigma 1
#define sigmaW 2000
#define sigmadmt 4000

int bqq[] = {1,1,1,2,2,2,3,3,3,4,4,4,1,1,1,2,2,2,3,3,3,4,4,4};
int q1W[] = {2,2,3,1,1,3,1,1,2,1,1,2,2,2,3,1,1,3,1,1,2,1,1,2};
int q2W[] = {3,4,4,3,4,4,2,4,4,2,3,3,3,4,4,3,4,4,2,4,4,2,3,3};
int ben[] = {4,3,2,4,3,1,4,2,1,3,2,1,4,3,2,4,3,1,4,2,1,3,2,1};

double Tab[30];

double min(double a, double b){
  if (a<b)
    return a;
  else return b;
}

vector<double> neutrinosolution(double E1, double E2, double E3, double E4, double E5){
  vector<double> solution;

  double EXm = Tab[9] + Tab[10]*(1-E1/Tab[0]) + Tab[11]*(1-E2/Tab[1]) + Tab[12]*(1-E3/Tab[2]) + Tab[13]*(1-E4/Tab[3]) + Tab[20]*(1-E5/Tab[23]);
  double EYm = Tab[14] + Tab[15]*(1-E1/Tab[0]) + Tab[16]*(1-E2/Tab[1]) + Tab[17]*(1-E3/Tab[2]) + Tab[18]*(1-E4/Tab[3]) + Tab[21]*(1-E5/Tab[23]);

  double ETm_corr = sqrt(EXm*EXm + EYm*EYm);

  double k  = 80.4e3*80.4e3 + 2.*(Tab[20]*EXm+Tab[21]*EYm);
  double a  = 4.*Tab[19]*Tab[19];
  double b  = -4.*k*Tab[22];
  double c  = 4.*Tab[23]*Tab[23]*ETm_corr*ETm_corr-k*k;

  TVector3 pben(Tab[13],Tab[18],Tab[24]);

  double delta = b*b - 4.*a*c;
  double Eneu, costhetabenn;

  if (delta < 0)
    delta = 0;
      
  double pZmis = (-b + Tab[25]*sqrt(delta))/2./a;
  TVector3 pmis3(EXm, EYm, pZmis);

  Eneu = pmis3.Mag();
  costhetabenn = pmis3*pben/Eneu/pben.Mag();

  solution.push_back(Eneu);
  solution.push_back(costhetabenn);
  solution.push_back(EXm);
  solution.push_back(EYm);
  solution.push_back(pZmis);
  return solution;
}

double chisq(double Ebqq, double Eq1W, double Eq2W, double Eben, double Ee)
{
  vector<double> neutrino = neutrinosolution(Ebqq,Eq1W,Eq2W,Eben,Ee);
  double En = neutrino.at(0);
  double cosbenn = neutrino.at(1);  

  TVector3 v3e_old(Tab[20],Tab[21],Tab[22]);
  TVector3 v3n(neutrino.at(2),neutrino.at(3),neutrino.at(4));

  double tbqq = (Tab[0] - Ebqq)/(ksigma*0.5)/sqrt(Ebqq*1000);
  double tq1W = (Tab[1] - Eq1W)/(ksigma*0.5)/sqrt(Eq1W*1000);
  double tq2W = (Tab[2] - Eq2W)/(ksigma*0.5)/sqrt(Eq2W*1000);
  double tben = (Tab[3] - Eben)/(ksigma*0.5)/sqrt(Eben*1000);
  double te   = (Tab[23] - Ee)/(ksigma*0.1)/sqrt(Ee*1000);
  
  double mass0 = 0, mass1 = 0, mass2 = 0, mass3 = 0, tmW = 0, tmen = 0, tmt = 0;
  double tmt1 = 0, tmt2 = 0, mintmt = 0;
  
  if (mjet0){
    double costhetane = v3e_old*v3n/v3e_old.Mag()/v3n.Mag();
    mass0 = sqrt(2*En*Ee*(1-costhetane));
    mass1 = sqrt(2*Eq1W*Eq2W*(1-Tab[4]));
    mass2 = sqrt(2*(Eq1W*Eq2W*(1-Tab[4]) + Eq1W*Ebqq*(1-Tab[5]) + Eq2W*Ebqq*(1-Tab[6])));
    mass3 = sqrt(80400*80400 + 2*(Eben*Tab[7]*(1-Tab[8]) + En*Eben*(1-cosbenn)));
  
    tmW = (mass1 - 80400)/sigmaW;
    tmen= (mass0 - 80400)/sigmaW; // Is 0 if Delta > 0, != 0 otherwise ==> good penalty ?
    tmt = (mass2 - mass3)/sigmadmt;
  }
  //double tmt1 = (mass2 - 172e3)/4000.;
  //double tmt2 = (mass3 - 172e3)/4000.;

  else{
    //Recalculate the masses with 4-vectors
    TVector3 v3bqq_old(Tab[10],Tab[15],Tab[28]);
    
    TVector3 v3q1W_old(Tab[11],Tab[16],Tab[27]);
    TVector3 v3q2W_old(Tab[12],Tab[17],Tab[26]);
    TVector3 v3ben_old(Tab[13],Tab[18],Tab[24]);

    double kbqq = Ebqq/Tab[0];
    double kq1W = Eq1W/Tab[1];
    double kq2W = Eq2W/Tab[2];
    double kben = Eben/Tab[3];
    double ke   = Ee/Tab[23];
  
    TVector3 v3bqq_new = kbqq*v3bqq_old;
    TVector3 v3q1W_new = kq1W*v3q1W_old;
    TVector3 v3q2W_new = kq2W*v3q2W_old;
    TVector3 v3ben_new = kben*v3ben_old;
    TVector3 v3e_new   = ke*v3e_old;

    TLorentzVector vbqqnew(v3bqq_new,Ebqq);
    TLorentzVector vq1Wnew(v3q1W_new,Eq1W);
    TLorentzVector vq2Wnew(v3q2W_new,Eq2W);
    TLorentzVector vbennew(v3ben_new,Eben);
    TLorentzVector venew(v3e_new,Ee);
    TLorentzVector vn(v3n,neutrino.at(0));

    TLorentzVector vWlep = venew + vn;
    mass0 = vWlep.M();

    TLorentzVector vWhad = vq1Wnew + vq2Wnew;
    mass1 = vWhad.M();
 
    TLorentzVector vth = vWhad + vbqqnew;
    mass2 = vth.M();
  
    TLorentzVector vtl = vWlep + vbennew;
    mass3 = vtl.M();
  
    tmW  = (mass1 - 80400)/sigmaW;
    tmen = (mass0 - 80400)/sigmaW; // Is 0 if Delta > 0, != 0 otherwise ==> good penalty ?
    tmt  = (mass2 - mass3)/sigmadmt;  
    tmt1 = (mass2 - 172e3)/sigmaW;
    tmt2 = (mass3 - 172e3)/sigmaW;

    mintmt = min(fabs(tmt1),fabs(tmt2));
  }
  
  if (Tab[29] == -1) //freefit = false
    return tbqq*tbqq + tq1W*tq1W + tq2W*tq2W + tben*tben + te*te + tmW*tmW + tmen*tmen + mintmt*mintmt;//this is used for the case of fit1top
  else
    return tbqq*tbqq + tq1W*tq1W + tq2W*tq2W + tben*tben + te*te + tmW*tmW + tmen*tmen + tmt*tmt;//this is used for the case of freefit

  //return tbqq*tbqq + tq1W*tq1W + tq2W*tq2W + tben*tben + te*te + tmW*tmW + tmen*tmen + tmt1*tmt1 + tmt2*tmt2;//this is used for the case of fit2top
}

void fcn(int& nDim, double* gout, double& result, double par[], int flg) {
  result = chisq(par[0], par[1], par[2], par[3], par[4]);
}

void configfunction(vector<TLorentzVector> par, double METetx, double METety, int config, bool freefit){

  //jet1
  TLorentzVector vbqq  = par.at(bqq[config] - 1);  
  //jet2
  TLorentzVector vq1W  = par.at(q1W[config] - 1);
  //jet3
  TLorentzVector vq2W  = par.at(q2W[config] - 1);
  //jet4
  TLorentzVector vben  = par.at(ben[config] - 1);

  //Change a bit to win time here! But in fact not :-D
  //We check if the inv mass of the sum of vq1W and vq2W 
  //can give the mW in a safe region
  //TLorentzVector vW = vq1W + vq2W;
  //double mW = vW.M();
  //if (fabs(mW-80.4e3)>3e4) return false;

  //jet1
  TVector3 vbqq3 = vbqq.Vect();
  double ptXbqq   = vbqq.Vect().X();
  double ptYbqq   = vbqq.Vect().Y();
  double ptZbqq   = vbqq.Vect().Z();
  double Ebqq     = vbqq.T();

  //jet2
  TVector3 vq1W3 = vq1W.Vect();
  double ptXq1W   = vq1W.Vect().X();
  double ptYq1W   = vq1W.Vect().Y();
  double ptZq1W   = vq1W.Vect().Z();
  double Eq1W     = vq1W.T();

  //jet3
  TVector3 vq2W3 = vq2W.Vect();
  double ptXq2W   = vq2W.Vect().X();
  double ptYq2W   = vq2W.Vect().Y();
  double ptZq2W   = vq2W.Vect().Z();
  double Eq2W     = vq2W.T();

  //jet4
  TVector3 vben3 = vben.Vect();
  double ptXben   = vben.Vect().X();
  double ptYben   = vben.Vect().Y();
  double ptZben   = vben.Vect().Z();
  double Eben     = vben.T();

  //electron
  TLorentzVector veW   = par.at(4);
  TVector3 veW3  = veW.Vect();
  
  double costhetaq1q2 = vq1W3*vq2W3/vq1W3.Mag()/vq2W3.Mag();
  double costhetaq1b  = vq1W3*vbqq3/vq1W3.Mag()/vbqq3.Mag();
  double costhetaq2b  = vq2W3*vbqq3/vq2W3.Mag()/vbqq3.Mag();
  double costhetabene = vben3*veW3/vben3.Mag()/veW3.Mag();
  
  double pTe = veW.Pt();
  double pXe = veW3.X();
  double pYe = veW3.Y();
  double pZe = veW3.Z();
  double Ee  = veW.E();

  Tab[0] = Ebqq;
  Tab[1] = Eq1W;
  Tab[2] = Eq2W;
  Tab[3] = Eben;
  Tab[4] = costhetaq1q2; 
  Tab[5] = costhetaq1b;
  Tab[6] = costhetaq2b;
  Tab[7] = veW.T();
  Tab[8] = costhetabene;

  Tab[9]  = METetx;
  Tab[10] = ptXbqq;
  Tab[11] = ptXq1W;
  Tab[12] = ptXq2W;
  Tab[13] = ptXben;

  Tab[14] = METety;
  Tab[15] = ptYbqq;
  Tab[16] = ptYq1W;
  Tab[17] = ptYq2W;
  Tab[18] = ptYben;

  Tab[19] = pTe;
  Tab[20] = pXe;
  Tab[21] = pYe;
  Tab[22] = pZe;
  Tab[23] = Ee;

  Tab[24] = ptZben;
  Tab[26] = ptZq2W;
  Tab[27] = ptZq1W;
  Tab[28] = ptZbqq;
  
  if (config >= 0 && config < 12)
    Tab[25] = 1;
  else 
    Tab[25] = -1;

  if (PRINT) {
    cout << "Config " << config << " initialised" << endl;
    for (int i = 0; i < 26; i++)
      cout << "Tab[" << i << "] : " << Tab[i] << endl;
  }
  
  if (freefit)
    Tab[29] = 1;
  else 
    Tab[29] = -1;

  //return true;
  
}

//Choose 4 jets for chi2
vector<vector<int> > choose4jets(vector<TLorentzVector> jets) {
  unsigned int nJets = jets.size();
  vector<vector<int> > ind;

  if (nJets==4){
    vector<int> ind1;
    for (unsigned int i1 = 0; i1 < nJets; i1++){      
      ind1.push_back(i1);
    }
    ind.push_back(ind1);
  }
  else {
    for (unsigned int i1 = 0; i1 < nJets; i1++){
      for (unsigned int i2=i1+1; i2 < nJets; i2++){
        for (unsigned int i3=i2+1; i3 < nJets; i3++){
          for (unsigned int i4=i3+1; i4 < nJets; i4++){
            vector<int> ind1;
            ind1.push_back(i1);
            ind1.push_back(i2);
            ind1.push_back(i3);
            ind1.push_back(i4);
            ind.push_back(ind1);
          }//end of loop i4
        }//end of loop i3
      }//end of loop i2
    }//end of loop i1
  }//end of else loop
  
  return ind;
  
}

chi2class chi2calculation(int iconfig)
{
  //Initialize TFitter
  TFitter* minimizer = new TFitter(5);
  double p1 = -1;
    
  minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);
  minimizer->ExecuteCommand("SET NOWarnings",0,0);

  minimizer->SetFCN(fcn);
  // Define the parameters
  // arg1 parameter number
  // arg2 parameter name
  // arg3 first guess at parameter value
  // arg4 estimated distance to minimum
  // arg5, arg6 ignore for now
  //     minimizer->SetParameter(0,"Ebqq",Tab[0],0.5*sqrt(Tab[0]*1000),0,0);
  //     minimizer->SetParameter(1,"Eq1W",Tab[1],0.5*sqrt(Tab[1]*1000),0,0);
  //     minimizer->SetParameter(2,"Eq2W",Tab[2],0.5*sqrt(Tab[2]*1000),0,0);
  //     minimizer->SetParameter(3,"Eben",Tab[3],0.5*sqrt(Tab[3]*1000),0,0);
  //     minimizer->SetParameter(4,"Ee",Tab[23],0.1*sqrt(Tab[23]*1000),0,0);
  minimizer->SetParameter(0,"Ebqq",Tab[0],0.5*sqrt(Tab[0]*1000),max(0.,Tab[0]-3*sqrt(Tab[0]*1000)),Tab[0]+3*sqrt(Tab[0]*1000));
  minimizer->SetParameter(1,"Eq1W",Tab[1],0.5*sqrt(Tab[1]*1000),max(0.,Tab[1]-3*sqrt(Tab[1]*1000)),Tab[1]+3*sqrt(Tab[1]*1000));
  minimizer->SetParameter(2,"Eq2W",Tab[2],0.5*sqrt(Tab[2]*1000),max(0.,Tab[2]-3*sqrt(Tab[2]*1000)),Tab[2]+3*sqrt(Tab[2]*1000));
  minimizer->SetParameter(3,"Eben",Tab[3],0.5*sqrt(Tab[3]*1000),max(0.,Tab[3]-3*sqrt(Tab[3]*1000)),Tab[3]+3*sqrt(Tab[3]*1000));
  minimizer->SetParameter(4,"Ee",Tab[23],0.1*sqrt(Tab[23]*1000),max(0.,Tab[23]-0.3*sqrt(Tab[23]*1000)),Tab[23]+0.3*sqrt(Tab[23]*1000));

  //Just fit 4 energies => fix the 5th one
  minimizer->FixParameter(4);
      
  // Run the simplex minimizer to get close to the minimum
  minimizer->ExecuteCommand("SIMPLEX",0,0);
  // Run the migrad minimizer (an extended Powell's method) to improve the
  // fit.
  minimizer->ExecuteCommand("MIGRAD",0,0);
  // Get the best fit values
  double BestE1 = minimizer->GetParameter(0);
  double BestE2 = minimizer->GetParameter(1);
  double BestE3 = minimizer->GetParameter(2);
  double BestE4 = minimizer->GetParameter(3);
  double BestE5 = minimizer->GetParameter(4);
  double Minimum = chisq(BestE1,BestE2,BestE3,BestE4,BestE5);

  chi2class Result(iconfig,Minimum,BestE1,BestE2,BestE3,BestE4,BestE5);

  delete minimizer;
  minimizer = 0;

  return Result;
}



/*std::set<double> chi2calculation (vector<TLorentzVector> vecparticles, double etx, double ety, int iconfig)
{
  //Initialize TFitter
  TFitter* minimizer = new TFitter(5);
  double p1 = -1;
  
  minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);
  minimizer->ExecuteCommand("SET NOWarnings",0,0);
  
  //if (!configfunction(vecpars2,MET_RefFinal_etx,MET_RefFinal_ety,i)) continue;
  
  configfunction(vecparticles,etx,ety,iconfig);  
      
  minimizer->SetFCN(fcn);
  // Define the parameters
  // arg1 parameter number
  // arg2 parameter name
  // arg3 first guess at parameter value
  // arg4 estimated distance to minimum
  // arg5, arg6 ignore for now
  //     minimizer->SetParameter(0,"Ebqq",Tab[0],0.5*sqrt(Tab[0]*1000),0,0);
  //     minimizer->SetParameter(1,"Eq1W",Tab[1],0.5*sqrt(Tab[1]*1000),0,0);
  //     minimizer->SetParameter(2,"Eq2W",Tab[2],0.5*sqrt(Tab[2]*1000),0,0);
  //     minimizer->SetParameter(3,"Eben",Tab[3],0.5*sqrt(Tab[3]*1000),0,0);
  //     minimizer->SetParameter(4,"Ee",Tab[23],0.1*sqrt(Tab[23]*1000),0,0);
  minimizer->SetParameter(0,"Ebqq",Tab[0],0.5*sqrt(Tab[0]*1000),max(0.,Tab[0]-3*sqrt(Tab[0]*1000)),Tab[0]+3*sqrt(Tab[0]*1000));
  minimizer->SetParameter(1,"Eq1W",Tab[1],0.5*sqrt(Tab[1]*1000),max(0.,Tab[1]-3*sqrt(Tab[1]*1000)),Tab[1]+3*sqrt(Tab[1]*1000));
  minimizer->SetParameter(2,"Eq2W",Tab[2],0.5*sqrt(Tab[2]*1000),max(0.,Tab[2]-3*sqrt(Tab[2]*1000)),Tab[2]+3*sqrt(Tab[2]*1000));
  minimizer->SetParameter(3,"Eben",Tab[3],0.5*sqrt(Tab[3]*1000),max(0.,Tab[3]-3*sqrt(Tab[3]*1000)),Tab[3]+3*sqrt(Tab[3]*1000));
  minimizer->SetParameter(4,"Ee",Tab[23],0.1*sqrt(Tab[23]*1000),max(0.,Tab[23]-0.3*sqrt(Tab[23]*1000)),Tab[23]+0.3*sqrt(Tab[23]*1000));
  
  //minimizer->FixParameter(4);
  
  // Run the simplex minimizer to get close to the minimum
  minimizer->ExecuteCommand("SIMPLEX",0,0);
  // Run the migrad minimizer (an extended Powell's method) to improve the
  // fit.
  minimizer->ExecuteCommand("MIGRAD",0,0);
  // Get the best fit values
  double bestE1 = minimizer->GetParameter(0);
  double bestE2 = minimizer->GetParameter(1);
  double bestE3 = minimizer->GetParameter(2);
  double bestE4 = minimizer->GetParameter(3);
  double bestE5 = minimizer->GetParameter(4);
  double minimum = chisq(bestE1,bestE2,bestE3,bestE4,bestE5);
  vector<double> neutrino = neutrinosolution(bestE1,bestE2,bestE3,bestE4,bestE5);
  
  double Eneu = neutrino.at(0);
  double costhetabenn = neutrino.at(1);
  double METx = neutrino.at(2);
  double METy = neutrino.at(3);
  double METz = neutrino.at(4);
  TVector3 vneu3fit(METx,METy,METz);
  TLorentzVector vneufit(vneu3fit,Eneu);
  TVector3 v3eold(Tab[20],Tab[21],Tab[22]);
  double ke   = bestE5/Tab[23];
  TVector3 v3e_fit   = ke*v3eold;
  TLorentzVector vefit(v3e_fit,bestE5);
  
  double costhetane = v3e_fit*vneu3fit/v3e_fit.Mag()/vneu3fit.Mag();
   
  double mWlepfit = 0, mWhadfit = 0, m2 = 0, m3 = 0;
  TLorentzVector vth, vtl;
  TLorentzVector vWhad, vWlep;
  
  if (mjet0) {
    mWlepfit = sqrt(2*Eneu*bestE5*(1-costhetane));
    mWhadfit = sqrt(2*bestE2*bestE3*(1-Tab[4]));
    m2 = sqrt(2*(bestE2*bestE3*(1-Tab[4]) + bestE2*bestE1*(1-Tab[5]) + bestE3*bestE1*(1-Tab[6])));
    m3 = sqrt(80400*80400 + 2*(bestE4*bestE5*(1-Tab[8]) + Eneu*bestE4*(1-costhetabenn)));
  }
  else {//m_jet is not 0
        //Recalculate the masses with 4-vectors
    TVector3 v3bqqold(Tab[10],Tab[15],Tab[28]);
    TVector3 v3q1Wold(Tab[11],Tab[16],Tab[27]);
    TVector3 v3q2Wold(Tab[12],Tab[17],Tab[26]);
    TVector3 v3benold(Tab[13],Tab[18],Tab[24]);
    
    double kbqq = bestE1/Tab[0];
    double kq1W = bestE2/Tab[1];
    double kq2W = bestE3/Tab[2];
    double kben = bestE4/Tab[3];
    
    TVector3 v3bqq_fit = kbqq*v3bqqold;
    TVector3 v3q1W_fit = kq1W*v3q1Wold;
    TVector3 v3q2W_fit = kq2W*v3q2Wold;
    TVector3 v3ben_fit = kben*v3benold;

    TLorentzVector vbqqfit(v3bqq_fit,bestE1);
    TLorentzVector vq1Wfit(v3q1W_fit,bestE2);
    TLorentzVector vq2Wfit(v3q2W_fit,bestE3);
    TLorentzVector vbenfit(v3ben_fit,bestE4);
        
    vWlep = vefit + vneufit;
    vWhad = vq1Wfit + vq2Wfit;
    vth = vWhad + vbqqfit;
    vtl = vWlep + vbenfit;
    
    //mWlepfit = vWlep.M();
    //mWhadfit = vWhad.M();
    //m2 = vth.M();//masse top hadronique
    //m3 = vtl.M();//masse top leptonique
  }

  delete minimizer;
  minimizer = 0;
}
*/