TwoGamma.cxx

Go to the documentation of this file.

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/TwoVector.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "TH1F.h"
#include "TMath.h"
#include <fstream>
#include <vector>
 
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
#include "VertexFit/Helix.h"
 
using CLHEP::Hep2Vector;
using CLHEP::Hep3Vector;
using CLHEP::HepLorentzVector;
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
 
#include "TwoGamma.h"
 
typedef std::vector<int>              Vint;
typedef std::vector<HepLorentzVector> Vp4;
const double                          ecms = 3.686;
const double                          velc = 299.792458;
const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
const double pai      = 3.1415926;
 
DECLARE_COMPONENT( TwoGamma )
////////////////////////////////////////////////////////////////////
TwoGamma::TwoGamma( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  //  m_tuple1 = 0;
  for ( int i = 0; i < 10; i++ ) { m_pass[i] = 0; }
  m_event = 0;
 
  Ndata1  = 0;
  Ndata2  = 0;
  m_runNo = 0;
 
  // Declare the properties
  declareProperty( "CmsEnergy", m_ecms = 3.686 );
 
  declareProperty( "max1", m_max1 = 1.2 );
  declareProperty( "max2", m_max2 = 0.8 );
  declareProperty( "costheta", m_costheta = 0.8 );
  declareProperty( "dphi1", m_dphi1 = 2.5 );
  declareProperty( "dphi2", m_dphi2 = 5 );
  declareProperty( "eff", m_eff = 0.07 );
  declareProperty( "sec", m_sec = 184.1 );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode TwoGamma::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  StatusCode status;
 
  NTuplePtr nt1( ntupleSvc(), "FILE1/DiGam" );
  if ( nt1 ) m_tuple1 = nt1;
  else
  {
    m_tuple1 =
        ntupleSvc()->book( "FILE1/DiGam", CLID_ColumnWiseTuple, "DiGam N-Tuple signal" );
    if ( m_tuple1 )
    {
      status = m_tuple1->addItem( "run", m_run );
      status = m_tuple1->addItem( "rec", m_rec );
      status = m_tuple1->addItem( "time", m_time );
 
      status = m_tuple1->addItem( "ngood", m_ngood );
      status = m_tuple1->addItem( "nchrg", m_nchrg );
 
      status = m_tuple1->addItem( "Cms_e1", m_e1 );
      status = m_tuple1->addItem( "Cms_e2", m_e2 );
      status = m_tuple1->addItem( "Cms_e", m_e );
      status = m_tuple1->addItem( "Cms_costheta1", m_costheta1 );
      status = m_tuple1->addItem( "Cms_costheta2", m_costheta2 );
      status = m_tuple1->addItem( "Cms_dltphi", m_dltphi );
      status = m_tuple1->addItem( "Cms_dltphi_1", m_dltphi_1 );
      status = m_tuple1->addItem( "Cms_dlttheta", m_dlttheta );
      status = m_tuple1->addItem( "Cms_phi1", m_phi1 );
      status = m_tuple1->addItem( "Cms_phi2", m_phi2 );
 
      status = m_tuple1->addItem( "Lab_e1", m_e1_lab );
      status = m_tuple1->addItem( "Lab_e2", m_e2_lab );
      status = m_tuple1->addItem( "Lab_e", m_e_lab );
      status = m_tuple1->addItem( "Lab_costheta1", m_costheta1_lab );
      status = m_tuple1->addItem( "Lab_costheta2", m_costheta2_lab );
      status = m_tuple1->addItem( "Lab_dltphi", m_dltphi_lab );
      status = m_tuple1->addItem( "Lab_dlttheta", m_dlttheta_lab );
      status = m_tuple1->addItem( "Lab_phi1", m_phi1_lab );
      status = m_tuple1->addItem( "Lab_phi2", m_phi2_lab );
 
      status = m_tuple1->addItem( "xBoost", m_xBoost );
      status = m_tuple1->addItem( "yBoost", m_yBoost );
      status = m_tuple1->addItem( "zBoost", m_zBoost );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple1:" << long( m_tuple1 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
  Ndata1  = 0;
  Ndata2  = 0;
  m_runNo = 0;
 
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
//*********************************************************************************************************************
StatusCode TwoGamma::execute() {
  StatusCode sc = StatusCode::SUCCESS;
  m_event++;
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  SmartDataPtr<EvtRecEvent>        evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  m_run     = eventHeader->runNumber();
  m_rec     = eventHeader->eventNumber();
  int runNo = m_run;
  int event = m_rec;
  int time  = eventHeader->time();
  if ( m_event > 1 && runNo != m_runNo ) return sc;
  m_runNo = runNo;
  m_time  = time;
 
  if ( m_rec % 1000 == 0 ) cout << "Run   " << m_run << "     Event   " << m_rec << endl;
 
  HepLorentzVector p4psip( 0.011 * m_ecms, 0.0, 0.005, m_ecms );
  double           psipBeta = ( p4psip.vect() ).mag() / ( p4psip.e() );
 
  m_pass[0] += 1;
 
  int NCharge = evtRecEvent->totalCharged();
  if ( NCharge != 0 ) return sc;
 
  m_pass[1] += 1;
 
  double Emax1 = 0.0;
  double Emax2 = 0.0;
  int    Imax[2];
 
  if ( ( ( evtRecEvent->totalTracks() - evtRecEvent->totalCharged() ) < 2 ) ||
       ( ( evtRecEvent->totalTracks() - evtRecEvent->totalCharged() ) > 15 ) )
    return sc;
  m_pass[2] += 1;
 
  HepLorentzVector p4Gam_1;
  HepLorentzVector p4Gam_2;
 
  for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
 
    HepLorentzVector p4Gam;
    double           Phi   = emcTrk->phi();
    double           Theta = emcTrk->theta();
    double           Ener  = emcTrk->energy();
 
    p4Gam.setPx( Ener * sin( Theta ) * cos( Phi ) );
    p4Gam.setPy( Ener * sin( Theta ) * sin( Phi ) );
    p4Gam.setPz( Ener * cos( Theta ) );
    p4Gam.setE( Ener );
    p4Gam.boost( -0.011, 0, -0.005 * 1.0 / 3.686 );
 
    if ( Ener > Emax2 )
    {
      Emax2   = Ener;
      Imax[1] = i;
      p4Gam_2 = p4Gam;
    }
    if ( Ener > Emax1 )
    {
      Emax2   = Emax1;
      p4Gam_2 = p4Gam_1;
      Imax[1] = Imax[0];
      Emax1   = Ener;
      p4Gam_1 = p4Gam;
      Imax[0] = i;
    }
  }
  m_e1_lab = Emax1;
  m_e2_lab = Emax2;
  m_e_lab  = Emax1 + Emax2;
  log << MSG::INFO << "Emax1 = " << Emax1 << "Emax2= " << Emax2 << endmsg;
  if ( Emax1 < m_max1 || Emax2 < m_max2 ) return sc;
  m_pass[3] += 1;
 
  // P4 in Lab
  double emcphi[2], emctht[2];
  for ( int i = 0; i < 2; i++ )
  {
    if ( i >= evtRecTrkCol->size() ) break;
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + Imax[i];
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
    emcphi[i]            = emcTrk->phi();
    emctht[i]            = emcTrk->theta();
  }
  double dltphi_lab   = ( fabs( emcphi[0] - emcphi[1] ) - pai ) * 180.0 / pai;
  double dlttheta_lab = ( fabs( emctht[0] + emctht[1] ) - pai ) * 180.0 / pai;
  m_costheta1_lab     = cos( emctht[0] );
  m_costheta2_lab     = cos( emctht[1] );
  m_phi1_lab          = emcphi[0] * 180.0 / pai;
  m_phi2_lab          = emcphi[1] * 180.0 / pai;
  m_dlttheta_lab      = dlttheta_lab;
  m_dltphi_lab        = dltphi_lab;
 
  if ( fabs( m_costheta1_lab ) > m_costheta || fabs( m_costheta2_lab ) > m_costheta )
    return sc;
  m_pass[4] += 1;
  // P4 in Lab
 
  // P4 in Cms
  double px1  = p4Gam_1.px();
  double py1  = p4Gam_1.py();
  double pz1  = p4Gam_1.pz();
  double pxy1 = sqrt( px1 * px1 + py1 * py1 );
  double e1   = p4Gam_1.e();
 
  double px2  = p4Gam_2.px();
  double py2  = p4Gam_2.py();
  double pz2  = p4Gam_2.pz();
  double pxy2 = sqrt( px2 * px2 + py2 * py2 );
  double e2   = p4Gam_2.e();
 
  double phi1 = 0;
  double phi2 = 0;
  if ( atan( py1 * 1.0 / px1 ) > 0 )
  {
    if ( px1 > 0 ) phi1 = atan( py1 * 1.0 / px1 );
    else phi1 = -( pai - atan( py1 * 1.0 / px1 ) );
  }
  else
  {
    if ( px1 > 0 ) phi1 = atan( py1 * 1.0 / px1 );
    else phi1 = pai + atan( py1 * 1.0 / px1 );
  }
 
  if ( atan( py2 * 1.0 / px2 ) > 0 )
  {
    if ( px2 > 0 ) phi2 = atan( py2 * 1.0 / px2 );
    else phi2 = -( pai - atan( py2 * 1.0 / px2 ) );
  }
  else
  {
    if ( px2 > 0 ) phi2 = atan( py2 * 1.0 / px2 );
    else phi2 = pai + atan( py2 * 1.0 / px2 );
  }
 
  double dltphi = ( fabs( phi1 - phi2 ) - pai ) * 180.0 / pai;
 
  double theta1 = 0;
  double theta2 = 0;
 
  if ( pz1 > 0 ) theta1 = atan( pxy1 * 1.0 / pz1 );
  else theta1 = ( pai + atan( pxy1 * 1.0 / pz1 ) );
 
  if ( pz2 > 0 ) theta2 = atan( pxy2 * 1.0 / pz2 );
  else theta2 = ( pai + atan( pxy2 * 1.0 / pz2 ) );
 
  double dlttheta = ( theta1 + theta2 - pai ) * 180.0 / pai;
 
  double xBoost = p4Gam_1.px() + p4Gam_2.px();
  double yBoost = p4Gam_1.py() + p4Gam_2.py();
  double zBoost = p4Gam_1.pz() + p4Gam_2.pz();
  m_xBoost      = xBoost;
  m_yBoost      = yBoost;
  m_zBoost      = zBoost;
 
  m_costheta1 = cos( theta1 );
  m_costheta2 = cos( theta2 );
  m_phi1      = phi1 * 180.0 / pai;
  m_phi2      = phi2 * 180.0 / pai;
  m_dlttheta  = dlttheta;
  m_dltphi    = dltphi;
  m_e1        = e1;
  m_e2        = e2;
  m_e         = e1 + e2;
  // P4 in Cms
 
  if ( fabs( m_dltphi ) < m_dphi1 ) Ndata1++;
  if ( fabs( m_dltphi ) > m_dphi1 && fabs( m_dltphi ) < m_dphi2 ) Ndata2++;
 
  m_tuple1->write().ignore();
  //    runNo = 0;
 
  return StatusCode::SUCCESS;
}
 
//*************************************************************************************************************
StatusCode TwoGamma::finalize() {
  char head[100];
  char foot[100] = ".txt";
  sprintf( head, "OffLineLum_%04d", m_runNo );
  strcat( head, foot );
  ofstream fout( head, ios::out );
  fout << "===============================================" << endl;
  fout << "DESIGNED For OffLine Luminosity BY 2Gam Events" << endl << endl;
  fout << "MANY THANKS to Prof. C.Z Yuan & Z.Y Wang" << endl << endl;
  fout << "                                        2009/05" << endl;
  //    fout<<"                               Charmonium Group"<<endl;
  fout << "===============================================" << endl << endl << endl;
 
  double lum = ( Ndata1 - Ndata2 ) * 1.0 / ( m_eff * m_sec );
  fout << "                  Table  List                  " << endl << endl;
  fout << "-----------------------------------------------" << endl;
  fout << "Firstly Energy Gamma " << m_max1 << " Gev" << endl;
  fout << "Secondly Energy Gamma    " << m_max2 << " Gev" << endl;
  fout << "Solid Angle Acceptance   " << m_costheta << endl;
  fout << "QED XSection     " << m_sec << " NB" << endl;
  fout << "Monte Carto Efficiency   " << m_eff * 100 << "%" << endl << endl;
  fout << "runNo                      Luminosity          " << endl << endl;
  fout << m_runNo << "                     " << lum << " nb^(-1)" << endl;
  fout << "-----------------------------------------------" << endl;
  fout.close();
 
  MsgStream log( msgSvc(), name() );
  cout << "in finalize()" << endl;
  cout << "all events   " << m_pass[0] << endl;
  cout << "NCharged==0  " << m_pass[1] << endl;
  cout << "Number of Gam    " << m_pass[2] << endl;
  cout << "N events " << m_pass[3] << endl;
  cout << "N events 0.8 " << m_pass[4] << endl;
  return StatusCode::SUCCESS;
}