source: huonglan/ANA2011/Thrust_func_collection.C

/* ************************************************************************************* */
/* Thrust_func_collection.C                                                              */
/*                                                                                       */     
/* Author:      Consorti Valerio, valerio.consorti@physik.uni-freiburg.de                */
/* created:     June 2010                                                                */
/* Description: A small collection of functions to compute the transverse thrust.        */
/*              The thrust value is normalized in the range [0,1] where 1 is for         */
/*              perfectly spherical events and 0 is for "pencil like" events             */
/*                                                                                       */
/* ************************************************************************************* */

#include "Thrust_func_collection.h"

//A modified step function used to compute the thrust axis
double Thrust_func_collection::epsilon(double x){
  if(x>=0) return 1;
  if(x<0) return -1;
  return -99;
}

//This function return the transverse thrust value for a given thrust axis
double Thrust_func_collection::thrustService(TVector2 &n, std::vector<TVector2> &obj){
  double num=0;
  double den=0;

  for(unsigned int i=0; i<obj.size(); i++){
    num+=TMath::Abs(n*obj[i]);
    den+=obj[i].Mod();
  }

  if(den==0) return -99;
  return num/den;
}


//Main function where the thrust axis which maximize the transverse thrust value is computed
vector<double> Thrust_func_collection::Calc_Thrust(std::vector<TLorentzVector> &Vector) {
  vector<double> vec0;

  if(Vector.size()==1) return vec0;

  std::vector<TVector2> obj;
 
  for(unsigned int index=0; index<Vector.size(); index++){
    obj.push_back(Vector[index].Vect().XYvector());
  }

  if(obj.size()==0) return vec0;

  TVector2 n;
  TVector2 np1(obj[0]);
  TVector2 np1T;
  np1=np1.Unit();
  bool loop_controller=true;

  int iloop=0;
  int iloop2=0;
  double thrust_old=0;
  double thrust_new=-1;
  double thrust_new_T=-1;
  double thrust_larger=0;
  double thrust_larger_T = 0.;
  TVector2 np1_choose, np1T_choose;
  double T=-99;
  double TT=-99;
  int n_lar_T_with_same_value=0;

//Initializing Thrust axis. By default it use the leading jet direction to initialize the thrust axis
  for(unsigned int i=1; i<obj.size();i++) np1+=obj[i];
  if(np1.X()==0 && np1.Y()==0) np1.Set(obj[0]);
  np1=np1.Unit();
 

//Computing transverse Thrust value. To ensure the convergence the algorithm requires to obtain 3 times the same maximum transverse Thrust value using 3 different initializations for the thrust axis
  do{
    loop_controller=true;

    //changing thrust axis initialization after the first loop
    if(iloop2>0){
      np1.Set(1.,0.);
      np1=np1.Rotate((iloop2-1));
      np1=np1.Unit();
    }
    iloop=0;

    //loop to find the axis which maximize the transverse Thrust value
    do{
      n=np1;
      iloop++;

      if(iloop==1000){
        loop_controller=false;
        thrust_new=-1;
        break;
      }

      thrust_old=thrustService(np1,obj);

      //the axis which maximize the thrust value is the result of the convergence of this series
      for(unsigned int i=0; i<obj.size(); i++){
        np1+=epsilon(n*obj[i])*obj[i];
      }

      np1=np1.Unit();
      thrust_new=thrustService(np1,obj);
      //Add by me!!!
      np1T.Set(-np1.Y(),np1.X());
      thrust_new_T = thrustService(np1T,obj);

    }while(fabs(thrust_new-thrust_old)>10e-18);


    if(thrust_new>thrust_larger && loop_controller ){
      thrust_larger=thrust_new;
      thrust_larger_T=thrust_new_T;//Add by me
      np1_choose = np1;
      np1T_choose = np1T;
      n_lar_T_with_same_value=0;
    }

    
    if(fabs((thrust_new-thrust_larger)<10e-18) && loop_controller ) n_lar_T_with_same_value++;

    iloop2++;

    loop_controller=true;

    if(iloop2==1000){
      loop_controller=false;
      std::cout<<"Thrust Computation: convergence 2 failed after: "<<iloop2<<" loop"<<std::endl;
      std::cout<<"n_lar_T_with_same_value: "<<n_lar_T_with_same_value<<std::endl;
      std::cout<<"(1-thrust_larger)/(1-2/TMath::Pi()): "<<(1-thrust_larger)/(1-2/TMath::Pi())<<std::endl;
      std::cout<<"Vector.size(): "<<Vector.size()<<std::endl;
      std::cout<<"i     pt      px      py      pxobj   pyobj"<<std::endl;
      for(unsigned int i=0; i<Vector.size();i++){
       std::cout<<i<<"  "<<Vector[i].Pt()<<"    "<<Vector[i].Px()<<"    "<<Vector[i].Py()<<"    "<<obj[i].X()<<"        "<<obj[i].Y()<<std::endl;
      }
      return vec0;
    }
//to accept the thrust value we require:
//1- it has to be into [2/pi,1] which is the range for the transverse thrust
//2- The maximum thrust value has to be obtain 3 times using 3 different initializations for the trust axis

  }while( ( thrust_larger<2/TMath::Pi() || n_lar_T_with_same_value<3 ) && loop_controller==true );

  T=thrust_larger;
  TT=thrust_larger_T;

  if(T<0 || T>1){
    if(Vector.size()>1){
    std::cout<<"Thrust Computation: ERROR, thrust value out of range: "<<T<<std::endl;

    std::cout<<"i       pt      px      py      pxobj   pyobj"<<std::endl;
      for(unsigned int i=0; i<Vector.size();i++){
       std::cout<<i<<"  "<<Vector[i].Pt()<<"    "<<Vector[i].Px()<<"    "<<Vector[i].Py()<<"    "<<obj[i].X()<<"        "<<obj[i].Y()<<std::endl;
      }
    }
    return vec0;
  }

  vector<double> res;
  res.push_back(T);
  res.push_back(TT);
  return res;
  
}