00001 #include <vector>
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/AlgFactory.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "GaudiKernel/SmartDataPtr.h"
00007 #include "GaudiKernel/IDataProviderSvc.h"
00008 #include "GaudiKernel/PropertyMgr.h"
00009 #include "VertexFit/IVertexDbSvc.h"
00010 #include "GaudiKernel/Bootstrap.h"
00011 #include "GaudiKernel/ISvcLocator.h"
00012
00013 #include "EventModel/EventModel.h"
00014 #include "EventModel/Event.h"
00015
00016 #include "EvtRecEvent/EvtRecEvent.h"
00017 #include "EvtRecEvent/EvtRecTrack.h"
00018 #include "DstEvent/TofHitStatus.h"
00019 #include "EventModel/EventHeader.h"
00020
00021 #include "TMath.h"
00022 #include "GaudiKernel/INTupleSvc.h"
00023 #include "GaudiKernel/NTuple.h"
00024 #include "GaudiKernel/Bootstrap.h"
00025 #include "GaudiKernel/IHistogramSvc.h"
00026 #include "CLHEP/Vector/ThreeVector.h"
00027 #include "CLHEP/Vector/LorentzVector.h"
00028 #include "CLHEP/Vector/TwoVector.h"
00029
00030 using CLHEP::Hep3Vector;
00031 using CLHEP::Hep2Vector;
00032 using CLHEP::HepLorentzVector;
00033 #include "CLHEP/Geometry/Point3D.h"
00034
00035 #include "VertexFit/KinematicFit.h"
00036 #include "VertexFit/VertexFit.h"
00037 #include "VertexFit/IVertexDbSvc.h"
00038 #include "ParticleID/ParticleID.h"
00039
00040
00041 #include "DQASelBhabha/DQASelBhabha.h"
00042
00043 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00044 typedef HepGeom::Point3D<double> HepPoint3D;
00045 #endif
00046 using CLHEP::HepLorentzVector;
00047 const double mpsi2s=3.68609;
00048 const double mpi = 0.13957;
00049 const double mk = 0.493677;
00050 const double xmass[5] = {0.000511, 0.105658, 0.139570,0.493677, 0.938272};
00051 const double velc = 299.792458;
00052 typedef std::vector<int> Vint;
00053 typedef std::vector<HepLorentzVector> Vp4;
00054
00055 static int counter[10]={0,0,0,0,0,0,0,0,0,0};
00056 static int nbhabha=0;
00058
00059 DQASelBhabha::DQASelBhabha(const std::string& name, ISvcLocator* pSvcLocator) :
00060 Algorithm(name, pSvcLocator) {
00061
00062
00063 declareProperty("writentuple",m_writentuple = false);
00064 declareProperty("ecms",m_ecms = 3.097);
00065 declareProperty("beamangle",m_beamangle = 0.022);
00066 declareProperty("Vr0cut", m_vr0cut=1.0);
00067 declareProperty("Vz0cut", m_vz0cut=8.0);
00068 declareProperty("Coscut", m_coscut=0.93);
00069
00070 declareProperty("EnergyThreshold", m_energyThreshold=0.04);
00071 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
00072 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
00073 declareProperty("GammaTrkCut", m_gammaTrkCut=20.0);
00074 declareProperty("GammaTLCut", m_gammatlCut=0);
00075 declareProperty("GammaTHCut", m_gammathCut=60);
00076
00077
00078 declareProperty ("acoll_e_cut", m_acoll_e_cut=6.);
00079 declareProperty ("acopl_e_cut", m_acopl_e_cut=6.);
00080 declareProperty ("poeb_e_cut", m_poeb_e_cut=0.5);
00081 declareProperty ("dtof_e_cut", m_dtof_e_cut=4.);
00082 declareProperty ("eoeb_e_cut", m_eoeb_e_cut=0.4);
00083 declareProperty ("etotal_e_cut", m_etotal_e_cut=0.8);
00084 declareProperty ("tpoeb_e_cut", m_tpoeb_e_cut=0.95);
00085 declareProperty ("tptotal_e_cut", m_tptotal_e_cut=0.16);
00086 declareProperty ("tetotal_e_cut", m_tetotal_e_cut=0.65);
00087
00088
00089
00090 declareProperty ("m_useEMConly", m_useEMConly= false);
00091 declareProperty ("m_usePID", m_usePID= false);
00092 declareProperty ("m_useMDC", m_useMDC= true);
00093 declareProperty ("m_useDEDX", m_useDEDX= false);
00094 declareProperty ("m_useTOF", m_useTOF= false);
00095 declareProperty ("m_useEMC", m_useEMC= true);
00096 declareProperty ("m_useMUC", m_useMUC= false);
00097
00098 }
00099
00100
00101
00102 StatusCode DQASelBhabha::initialize(){
00103 MsgStream log(msgSvc(), name());
00104
00105 log << MSG::INFO << "in initialize()" << endmsg;
00106 StatusCode status;
00107 status = service("THistSvc", m_thistsvc);
00108 if(status.isFailure() ){
00109 log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
00110 return status;
00111 }
00112
00113
00114 m_ee_mass = new TH1F( "ee_mass", "ee_mass", 80, m_ecms-0.3, m_ecms+0.5 );
00115 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_mass", m_ee_mass);
00116 m_ee_acoll = new TH1F( "ee_acoll", "ee_acoll", 60, 0, 6 );
00117 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_acoll", m_ee_acoll);
00118 m_ee_eop_ep = new TH1F( "ee_eop_ep", "ee_eop_ep", 100,0.4,1.4 );
00119 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eop_ep", m_ee_eop_ep);
00120 m_ee_eop_em = new TH1F( "ee_eop_em", "ee_eop_em", 100,0.4,1.4 );
00121 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eop_em", m_ee_eop_em);
00122 m_ee_costheta_ep = new TH1F( "ee_costheta_ep", "ee_costheta_ep", 100,-1,1 );
00123 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_costheta_ep", m_ee_costheta_ep);
00124 m_ee_costheta_em = new TH1F( "ee_costheta_em", "ee_costheta_em", 100,-1,1 );
00125 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_costheta_em", m_ee_costheta_em);
00126
00127 m_ee_phi_ep = new TH1F( "ee_phi_ep", "ee_phi_ep", 120,-3.2,3.2 );
00128 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_phi_ep", m_ee_phi_ep);
00129 m_ee_phi_em = new TH1F( "ee_phi_em", "ee_phi_em", 120,-3.2,3.2 );
00130 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_phi_em", m_ee_phi_em);
00131
00132 m_ee_nneu = new TH1I( "ee_nneu", "ee_nneu", 5,0,5 );
00133 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_nneu", m_ee_nneu);
00134
00135 m_ee_eemc_ep=new TH1F("ee_eemc_ep","ee_eemc_ep",100,1.0,2.0);
00136 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eemc_ep", m_ee_eemc_ep);
00137 m_ee_eemc_em=new TH1F("ee_eemc_em","ee_eemc_em",100,1.0,2.0);
00138 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_eemc_em", m_ee_eemc_em);
00139 m_ee_x_ep=new TH1F("ee_x_ep","ee_x_ep",100,-1.0,1.0);
00140 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_x_ep", m_ee_x_ep);
00141 m_ee_y_ep=new TH1F("ee_y_ep","ee_y_ep",100,-1.0,1.0);
00142 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_y_ep", m_ee_y_ep);
00143 m_ee_z_ep=new TH1F("ee_z_ep","ee_z_ep",100,-10.0,10.0);
00144 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_z_ep", m_ee_z_ep);
00145 m_ee_x_em=new TH1F("ee_x_em","ee_x_em",100,-1.0,1.0);
00146 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_x_em", m_ee_x_em);
00147 m_ee_y_em=new TH1F("ee_y_em","ee_y_em",100,-1.0,1.0);
00148 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_y_em", m_ee_y_em);
00149 m_ee_z_em=new TH1F("ee_z_em","ee_z_em",100,-10.0,10.0);
00150 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_z_em", m_ee_z_em);
00151
00152 m_ee_px_ep=new TH1F("ee_px_ep","ee_px_ep",200,-2.0,2.0);
00153 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_px_ep", m_ee_px_ep);
00154 m_ee_py_ep=new TH1F("ee_py_ep","ee_py_ep",200,-2.0,2.0);
00155 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_py_ep", m_ee_py_ep);
00156 m_ee_pz_ep=new TH1F("ee_pz_ep","ee_pz_ep",200,-2.0,2.0);
00157 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pz_ep", m_ee_pz_ep);
00158 m_ee_p_ep=new TH1F("ee_p_ep","ee_p_ep",100,1.0,2.0);
00159 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_p_ep", m_ee_p_ep);
00160 m_ee_px_em=new TH1F("ee_px_em","ee_px_em",100,-2.0,2.0);
00161 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_px_em", m_ee_px_em);
00162 m_ee_py_em=new TH1F("ee_py_em","ee_py_em",100,-2.0,2.0);
00163 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_py_em", m_ee_py_em);
00164 m_ee_pz_em=new TH1F("ee_pz_em","ee_pz_em",100,-2.0,2.0);
00165 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pz_em", m_ee_pz_em);
00166 m_ee_p_em=new TH1F("ee_p_em","ee_p_em",100,1.0,2.0);
00167 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_p_em", m_ee_p_em);
00168 m_ee_deltatof=new TH1F("ee_deltatof","ee_deltatof",50,0.0,10.0);
00169 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_deltatof", m_ee_deltatof);
00170
00171 m_ee_pidchidedx_ep=new TH1F("ee_pidchidedx_ep","ee_pidchidedx_ep",160,-4,4);
00172 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchidedx_ep", m_ee_pidchidedx_ep);
00173 m_ee_pidchidedx_em=new TH1F("ee_pidchidedx_em","ee_pidchidedx_em",160,-4,4);
00174 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchidedx_em", m_ee_pidchidedx_em);
00175 m_ee_pidchitof1_ep=new TH1F("ee_pidchitof1_ep","ee_pidchitof1_ep",160,-4,4);
00176 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof1_ep", m_ee_pidchitof1_ep);
00177 m_ee_pidchitof1_em=new TH1F("ee_pidchitof1_em","ee_pidchitof1_em",160,-4,4);
00178 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof1_em", m_ee_pidchitof1_em);
00179 m_ee_pidchitof2_ep=new TH1F("ee_pidchitof2_ep","ee_pidchitof2_ep",160,-4,4);
00180 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof2_ep", m_ee_pidchitof2_ep);
00181 m_ee_pidchitof2_em=new TH1F("ee_pidchitof2_em","ee_pidchitof2_em",160,-4,4);
00182 status = m_thistsvc->regHist("/DQAHist/Bhabha/ee_pidchitof2_em", m_ee_pidchitof2_em);
00183
00184
00185
00186
00187
00188
00189
00190 NTuplePtr nt1(ntupleSvc(), "DQAFILE/Bhabha");
00191 if ( nt1 ) m_tuple1 = nt1;
00192 else {
00193 m_tuple1 = ntupleSvc()->book ("DQAFILE/Bhabha", CLID_ColumnWiseTuple, "N-Tuple");
00194 if ( m_tuple1 ) {
00195 status = m_tuple1->addItem ("run", m_run);
00196 status = m_tuple1->addItem ("rec", m_rec);
00197 status = m_tuple1->addItem ("Nchrg", m_ncharg);
00198 status = m_tuple1->addItem ("Nneu", m_nneu,0,40);
00199 status = m_tuple1->addItem ("NGch", m_ngch, 0, 40);
00200 status = m_tuple1->addItem ("NGam", m_nGam);
00201
00202
00203 status = m_tuple1->addItem ("bhabhatag", m_bhabhatag);
00204
00205 status = m_tuple1->addItem ("acoll", m_acoll);
00206 status = m_tuple1->addItem ("acopl", m_acopl);
00207 status = m_tuple1->addItem ("deltatof", m_deltatof);
00208 status = m_tuple1->addItem ("eop1", m_eop1);
00209 status = m_tuple1->addItem ("eop2", m_eop2);
00210 status = m_tuple1->addItem ("eoeb1", m_eoeb1);
00211 status = m_tuple1->addItem ("eoeb2", m_eoeb2);
00212 status = m_tuple1->addItem ("poeb1", m_poeb1);
00213 status = m_tuple1->addItem ("poeb2", m_poeb2);
00214 status = m_tuple1->addItem ("etoeb1", m_etoeb1);
00215 status = m_tuple1->addItem ("etoeb2", m_etoeb2);
00216 status = m_tuple1->addItem ("mucinfo1", m_mucinfo1);
00217 status = m_tuple1->addItem ("mucinfo2", m_mucinfo2);
00218
00219
00220 status = m_tuple1->addIndexedItem ("delang",m_nneu, m_delang);
00221 status = m_tuple1->addIndexedItem ("delphi",m_nneu, m_delphi);
00222 status = m_tuple1->addIndexedItem ("delthe",m_nneu, m_delthe);
00223 status = m_tuple1->addIndexedItem ("npart",m_nneu, m_npart);
00224 status = m_tuple1->addIndexedItem ("nemchits",m_nneu, m_nemchits);
00225 status = m_tuple1->addIndexedItem ("module",m_nneu, m_module);
00226 status = m_tuple1->addIndexedItem ("x",m_nneu, m_x);
00227 status = m_tuple1->addIndexedItem ("y",m_nneu, m_y);
00228 status = m_tuple1->addIndexedItem ("z",m_nneu, m_z);
00229 status = m_tuple1->addIndexedItem ("px",m_nneu, m_px);
00230 status = m_tuple1->addIndexedItem ("py",m_nneu, m_py);
00231 status = m_tuple1->addIndexedItem ("pz",m_nneu, m_pz);
00232 status = m_tuple1->addIndexedItem ("theta",m_nneu, m_theta);
00233 status = m_tuple1->addIndexedItem ("phi",m_nneu, m_phi);
00234 status = m_tuple1->addIndexedItem ("dx",m_nneu, m_dx);
00235 status = m_tuple1->addIndexedItem ("dy",m_nneu, m_dy);
00236 status = m_tuple1->addIndexedItem ("dz",m_nneu, m_dz);
00237 status = m_tuple1->addIndexedItem ("dtheta",m_nneu, m_dtheta);
00238 status = m_tuple1->addIndexedItem ("dphi",m_nneu, m_dphi);
00239 status = m_tuple1->addIndexedItem ("energy",m_nneu, m_energy);
00240 status = m_tuple1->addIndexedItem ("dE",m_nneu, m_dE);
00241 status = m_tuple1->addIndexedItem ("eSeed",m_nneu, m_eSeed);
00242 status = m_tuple1->addIndexedItem ("nSeed",m_nneu, m_nSeed);
00243 status = m_tuple1->addIndexedItem ("e3x3",m_nneu, m_e3x3);
00244 status = m_tuple1->addIndexedItem ("e5x5",m_nneu, m_e5x5);
00245 status = m_tuple1->addIndexedItem ("secondMoment",m_nneu, m_secondMoment);
00246 status = m_tuple1->addIndexedItem ("latMoment",m_nneu, m_latMoment);
00247 status = m_tuple1->addIndexedItem ("a20Moment",m_nneu, m_a20Moment);
00248 status = m_tuple1->addIndexedItem ("a42Moment",m_nneu, m_a42Moment);
00249 status = m_tuple1->addIndexedItem ("getTime",m_nneu, m_getTime);
00250 status = m_tuple1->addIndexedItem ("getEAll",m_nneu, m_getEAll);
00251
00252
00253
00254 status = m_tuple1->addIndexedItem("charge", m_ngch, m_charge);
00255 status = m_tuple1->addIndexedItem ("vx", m_ngch, m_vx0);
00256 status = m_tuple1->addIndexedItem ("vy", m_ngch, m_vy0);
00257 status = m_tuple1->addIndexedItem ("vz", m_ngch, m_vz0);
00258
00259
00260 status = m_tuple1->addIndexedItem ("px", m_ngch, m_px) ;
00261 status = m_tuple1->addIndexedItem ("py", m_ngch, m_py) ;
00262 status = m_tuple1->addIndexedItem ("pz", m_ngch, m_pz) ;
00263 status = m_tuple1->addIndexedItem ("p", m_ngch, m_p) ;
00264
00265
00266
00267 status = m_tuple1->addIndexedItem ("kal_vx", m_ngch, m_kal_vx0);
00268 status = m_tuple1->addIndexedItem ("kal_vy", m_ngch, m_kal_vy0);
00269 status = m_tuple1->addIndexedItem ("kal_vz", m_ngch, m_kal_vz0);
00270
00271
00272 status = m_tuple1->addIndexedItem ("kal_px", m_ngch, m_kal_px) ;
00273 status = m_tuple1->addIndexedItem ("kal_py", m_ngch, m_kal_py) ;
00274 status = m_tuple1->addIndexedItem ("kal_pz", m_ngch, m_kal_pz) ;
00275 status = m_tuple1->addIndexedItem ("kal_p", m_ngch, m_kal_p) ;
00276
00277
00278 status = m_tuple1->addIndexedItem ("probPH" , m_ngch, m_probPH) ;
00279 status = m_tuple1->addIndexedItem ("normPH" , m_ngch, m_normPH) ;
00280 status = m_tuple1->addIndexedItem ("chie" , m_ngch, m_chie) ;
00281 status = m_tuple1->addIndexedItem ("chimu" , m_ngch, m_chimu) ;
00282 status = m_tuple1->addIndexedItem ("chipi" , m_ngch, m_chipi) ;
00283 status = m_tuple1->addIndexedItem ("chik" , m_ngch, m_chik) ;
00284 status = m_tuple1->addIndexedItem ("chip" , m_ngch, m_chip) ;
00285 status = m_tuple1->addIndexedItem ("ghit" , m_ngch, m_ghit) ;
00286 status = m_tuple1->addIndexedItem ("thit" , m_ngch, m_thit) ;
00287
00288 status = m_tuple1->addIndexedItem ("e_emc" , m_ngch, m_e_emc) ;
00289 status = m_tuple1->addIndexedItem ("phi_emc" , m_ngch, m_phi_emc) ;
00290 status = m_tuple1->addIndexedItem ("theta_emc" , m_ngch, m_theta_emc) ;
00291
00292 status = m_tuple1->addIndexedItem ("nhit_muc" , m_ngch, m_nhit_muc) ;
00293 status = m_tuple1->addIndexedItem ("nlay_muc" , m_ngch, m_nlay_muc) ;
00294 status = m_tuple1->addIndexedItem ("t_btof" , m_ngch, m_t_btof );
00295 status = m_tuple1->addIndexedItem ("t_etof" , m_ngch, m_t_etof );
00296 status = m_tuple1->addIndexedItem ("qual_etof" , m_ngch, m_qual_etof );
00297 status = m_tuple1->addIndexedItem ("tof_etof" , m_ngch, m_tof_etof );
00298 status = m_tuple1->addIndexedItem ("te_etof" , m_ngch, m_te_etof );
00299 status = m_tuple1->addIndexedItem ("tmu_etof" , m_ngch, m_tmu_etof );
00300 status = m_tuple1->addIndexedItem ("tpi_etof" , m_ngch, m_tpi_etof );
00301 status = m_tuple1->addIndexedItem ("tk_etof" , m_ngch, m_tk_etof );
00302 status = m_tuple1->addIndexedItem ("tp_etof" , m_ngch, m_tp_etof );
00303
00304 status = m_tuple1->addIndexedItem ("qual_btof1", m_ngch, m_qual_btof1 );
00305 status = m_tuple1->addIndexedItem ("tof_btof1" , m_ngch, m_tof_btof1 );
00306 status = m_tuple1->addIndexedItem ("te_btof1" , m_ngch, m_te_btof1 );
00307 status = m_tuple1->addIndexedItem ("tmu_btof1" , m_ngch, m_tmu_btof1 );
00308 status = m_tuple1->addIndexedItem ("tpi_btof1" , m_ngch, m_tpi_btof1 );
00309 status = m_tuple1->addIndexedItem ("tk_btof1" , m_ngch, m_tk_btof1 );
00310 status = m_tuple1->addIndexedItem ("tp_btof1" , m_ngch, m_tp_btof1 );
00311
00312 status = m_tuple1->addIndexedItem ("qual_btof2", m_ngch, m_qual_btof2 );
00313 status = m_tuple1->addIndexedItem ("tof_btof2" , m_ngch, m_tof_btof2 );
00314 status = m_tuple1->addIndexedItem ("te_btof2" , m_ngch, m_te_btof2 );
00315 status = m_tuple1->addIndexedItem ("tmu_btof2" , m_ngch, m_tmu_btof2 );
00316 status = m_tuple1->addIndexedItem ("tpi_btof2" , m_ngch, m_tpi_btof2 );
00317 status = m_tuple1->addIndexedItem ("tk_btof2" , m_ngch, m_tk_btof2 );
00318 status = m_tuple1->addIndexedItem ("tp_btof2" , m_ngch, m_tp_btof2 );
00319 status = m_tuple1->addIndexedItem ("pidcode" , m_ngch, m_pidcode);
00320 status = m_tuple1->addIndexedItem ("pidprob" , m_ngch, m_pidprob);
00321 status = m_tuple1->addIndexedItem ("pidchiDedx" , m_ngch, m_pidchiDedx);
00322 status = m_tuple1->addIndexedItem ("pidchiTof1" , m_ngch, m_pidchiTof1);
00323 status = m_tuple1->addIndexedItem ("pidchiTof2" , m_ngch, m_pidchiTof2);
00324
00325 status = m_tuple1->addItem ("px_cms_ep", m_px_cms_ep);
00326 status = m_tuple1->addItem ("py_cms_ep", m_py_cms_ep);
00327 status = m_tuple1->addItem ("pz_cms_ep", m_pz_cms_ep);
00328 status = m_tuple1->addItem ("e_cms_ep", m_e_cms_ep);
00329 status = m_tuple1->addItem ("cos_ep", m_cos_ep);
00330 status = m_tuple1->addItem ("px_cms_em", m_px_cms_em);
00331 status = m_tuple1->addItem ("py_cms_em", m_py_cms_em);
00332 status = m_tuple1->addItem ("pz_cms_em", m_pz_cms_em);
00333 status = m_tuple1->addItem ("e_cms_em", m_e_cms_em);
00334 status = m_tuple1->addItem ("cos_em", m_cos_em);
00335 status = m_tuple1->addItem ("mass_ee", m_mass_ee);
00336 status = m_tuple1->addItem ("emax", m_emax);
00337 status = m_tuple1->addItem ("esum", m_esum);
00338 status = m_tuple1->addItem ( "npip", m_npip );
00339 status = m_tuple1->addItem ( "npim", m_npim );
00340 status = m_tuple1->addItem ( "nkp", m_nkp );
00341 status = m_tuple1->addItem ( "nkm", m_nkm );
00342 status = m_tuple1->addItem ( "np", m_np );
00343 status = m_tuple1->addItem ( "npb", m_npb );
00344
00345 status = m_tuple1->addItem ( "nep", m_nep );
00346 status = m_tuple1->addItem ( "nem", m_nem );
00347 status = m_tuple1->addItem ( "nmup", m_nmup );
00348 status = m_tuple1->addItem ( "nmum", m_nmum );
00349
00350 }
00351 else {
00352 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
00353 return StatusCode::FAILURE;
00354 }
00355 }
00356
00357
00358
00359
00360
00361 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00362 return StatusCode::SUCCESS;
00363
00364
00365
00366
00367 }
00368
00369
00370 StatusCode DQASelBhabha::execute() {
00371
00372 const double beamEnergy = m_ecms/2.;
00373 const HepLorentzVector p_cms(m_ecms*sin(m_beamangle*0.5),0.0,0.0,m_ecms);
00374 const Hep3Vector u_cms = -p_cms.boostVector();
00375 MsgStream log(msgSvc(), name());
00376 log << MSG::INFO << "in execute()" << endreq;
00377
00378 setFilterPassed(false);
00379
00380 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00381 if (!eventHeader)
00382 {
00383 log << MSG::FATAL << "Could not find Event Header" << endreq;
00384 return StatusCode::SUCCESS;
00385 }
00386
00387 m_run = eventHeader->runNumber();
00388 m_rec = eventHeader->eventNumber();
00389
00390
00391
00392
00393 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00394 if (!evtRecEvent)
00395 {
00396 log << MSG::FATAL << "Could not find EvtRecEvent" << endreq;
00397 return StatusCode::SUCCESS;
00398 }
00399 log << MSG::INFO <<"ncharg, nneu, tottks = "
00400 << evtRecEvent->totalCharged() << " , "
00401 << evtRecEvent->totalNeutral() << " , "
00402 << evtRecEvent->totalTracks() <<endreq;
00403
00404 m_ncharg = evtRecEvent->totalCharged();
00405
00406 m_nneu = evtRecEvent->totalNeutral();
00407
00408
00409
00410 HepPoint3D vx(0., 0., 0.);
00411 HepSymMatrix Evx(3, 0);
00412 IVertexDbSvc* vtxsvc;
00413 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00414 if(vtxsvc->isVertexValid()){
00415 double* dbv = vtxsvc->PrimaryVertex();
00416 double* vv = vtxsvc->SigmaPrimaryVertex();
00417
00418
00419
00420 vx.setX(dbv[0]);
00421 vx.setY(dbv[1]);
00422 vx.setZ(dbv[2]);
00423 Evx[0][0]=vv[0]*vv[0];
00424 Evx[0][1]=vv[0]*vv[1];
00425 Evx[1][1]=vv[1]*vv[1];
00426 Evx[1][2]=vv[1]*vv[2];
00427 Evx[2][2]=vv[2]*vv[2];
00428 }
00429
00430 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00431 if (!evtRecTrkCol)
00432 {
00433 log << MSG::FATAL << "Could not find EvtRecTrackCol" << endreq;
00434 return StatusCode::SUCCESS;
00435 }
00436 Vint iGood;
00437 iGood.clear();
00438
00439 int nCharge = 0;
00440
00441 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00442 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00443 if(!(*itTrk)->isMdcTrackValid()) continue;
00444 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00445
00446 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00447 double pch=mdcTrk->p();
00448 double x0=mdcTrk->x();
00449 double y0=mdcTrk->y();
00450 double z0=mdcTrk->z();
00451 double phi0=mdcTrk->helix(1);
00452 double xv=vx.x();
00453 double yv=vx.y();
00454 double zv=vx.z();
00455 double Rxy=(x0-xv)*cos(phi0)+(y0-yv)*sin(phi0);
00456 double m_vx0 = x0;
00457 double m_vy0 = y0;
00458 double m_vz0 = z0;
00459 double m_vr0 = Rxy;
00460 if(fabs(z0) >= m_vz0cut) continue;
00461 if(fabs(Rxy) >= m_vr0cut) continue;
00462
00463
00464 if(fabs(m_vz0) >= m_vz0cut) continue;
00465 if(m_vr0 >= m_vr0cut) continue;
00466
00467
00468
00469
00470 iGood.push_back(i);
00471 nCharge += mdcTrk->charge();
00472
00473 }
00474
00475
00476
00477
00478
00479
00480
00481
00482 int nGood = iGood.size();
00483 m_ngch=nGood;
00484 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
00485
00486 if((nGood != 2)||(nCharge!=0)){
00487 return StatusCode::SUCCESS;
00488 }
00489
00490 counter[1]++;
00491
00492
00493
00494
00495 Vint ipip, ipim, iep,iem,imup,imum;
00496 ipip.clear();
00497 ipim.clear();
00498 iep.clear();
00499 iem.clear();
00500 imup.clear();
00501 imum.clear();
00502
00503
00504 ParticleID *pid = ParticleID::instance();
00505 for(int i = 0; i < m_ngch; i++) {
00506 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00507
00508 pid->init();
00509 pid->setMethod(pid->methodProbability());
00510 pid->setChiMinCut(4);
00511 pid->setRecTrack(*itTrk);
00512 pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2());
00513 pid->identify(pid->onlyElectron()|pid->onlyMuon()|pid->onlyPion());
00514 pid->calculate();
00515 if(!(pid->IsPidInfoValid())) continue;
00516 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00519 double prob_pi = pid->probPion();
00520 double prob_K = pid->probKaon();
00521 double prob_p = pid->probProton();
00522 double prob_e = pid->probElectron();
00523 double prob_mu = pid->probMuon();
00524
00525 HepLorentzVector ptrk;
00526 ptrk.setPx(mdcTrk->px()) ;
00527 ptrk.setPy(mdcTrk->py()) ;
00528 ptrk.setPz(mdcTrk->pz()) ;
00529 double p3 = ptrk.mag() ;
00530
00531 m_pidcode[i]=0;
00532 m_pidprob[i]=pid->prob(0);
00533 m_pidchiDedx[i]=pid->chiDedx(0);
00534 m_pidchiTof1[i]=pid->chiTof1(0);
00535 m_pidchiTof2[i]=pid->chiTof2(0);
00536 if(mdcTrk->charge() > 0) {
00537 iep.push_back(iGood[i]);
00538
00539 }
00540 if (mdcTrk->charge() < 0) {
00541 iem.push_back(iGood[i]);
00542
00543 }
00544
00545
00546
00547
00548 }
00549
00550 m_nep = iep.size() ;
00551 m_nem = iem.size() ;
00552 m_nmup = imup.size() ;
00553 m_nmum = imum.size() ;
00554
00555 counter[2]++;
00556
00557
00558
00559
00560 Vint iGam;
00561 iGam.clear();
00562 int iphoton=0;
00563 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
00564 if(i>=evtRecTrkCol->size())break;
00565 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00566 if(!(*itTrk)->isEmcShowerValid()) continue;
00567 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00568 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
00569
00570 RecEmcID showerId = emcTrk->getShowerId();
00571 unsigned int npart = EmcID::barrel_ec(showerId);
00572 int n = emcTrk->numHits();
00573 int module=emcTrk->module();
00574 double x = emcTrk->x();
00575 double y = emcTrk->y();
00576 double z = emcTrk->z();
00577 double dx = emcTrk->dx();
00578 double dy = emcTrk->dy();
00579 double dth = emcTrk->dtheta();
00580 double dph = emcTrk->dphi();
00581 double dz = emcTrk->dz();
00582 double energy = emcTrk->energy();
00583 double dE = emcTrk->dE();
00584 double eSeed = emcTrk->eSeed();
00585 double e3x3 = emcTrk->e3x3();
00586 double e5x5 = emcTrk->e5x5();
00587 double secondMoment = emcTrk->secondMoment();
00588 double latMoment = emcTrk->latMoment();
00589 double getTime = emcTrk->time();
00590 double getEAll = emcTrk->getEAll();
00591 double a20Moment = emcTrk->a20Moment();
00592 double a42Moment = emcTrk->a42Moment();
00593
00594
00595
00596
00597
00598 double nseed=0;
00599 HepPoint3D EmcPos(x,y,z);
00600 m_nemchits[iphoton]=n;
00601 m_npart[iphoton]=npart;
00602 m_module[iphoton]=module;
00603 m_theta[iphoton]=EmcPos.theta();
00604 m_phi[iphoton]=EmcPos.phi();
00605 m_x[iphoton]=x;
00606 m_y[iphoton]=y;
00607 m_z[iphoton]=z;
00608 m_dx[iphoton]=dx;
00609 m_dy[iphoton]=dy;
00610 m_dz[iphoton]=dz;
00611 m_dtheta[iphoton]=dth;
00612 m_dphi[iphoton]=dph;
00613 m_energy[iphoton]=energy;
00614 m_dE[iphoton]=dE;
00615 m_eSeed[iphoton]=eSeed;
00616 m_nSeed[iphoton]=nseed;
00617 m_e3x3[iphoton]=e3x3;
00618 m_e5x5[iphoton]=e5x5;
00619 m_secondMoment[iphoton]=secondMoment;
00620 m_latMoment[iphoton]=latMoment;
00621 m_getTime[iphoton]=getTime;
00622 m_getEAll[iphoton]=getEAll;
00623 m_a20Moment[iphoton]=a20Moment;
00624 m_a42Moment[iphoton]=a42Moment;
00625
00626
00627
00628
00629
00630
00631 double dthe = 200.;
00632 double dphi = 200.;
00633 double dang = 200.;
00634
00635
00636 for(int j = 0; j < nGood; j++) {
00637
00638
00639 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() +iGood[j];
00640 if (!(*jtTrk)->isMdcTrackValid()) continue;
00641 RecMdcTrack *jtmdcTrk = (*jtTrk)->mdcTrack();
00642 double jtcharge = jtmdcTrk->charge();
00643 if(!(*jtTrk)->isExtTrackValid()) continue;
00644 RecExtTrack *extTrk = (*jtTrk)->extTrack();
00645 if(extTrk->emcVolumeNumber() == -1) continue;
00646 Hep3Vector extpos = extTrk->emcPosition();
00647
00648 double angd = extpos.angle(emcpos);
00649 double thed = extpos.theta() - emcpos.theta();
00650 double phid = extpos.deltaPhi(emcpos);
00651 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00652 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
00653
00654 if(fabs(thed) < fabs(dthe)) dthe = thed;
00655 if(fabs(phid) < fabs(dphi)) dphi = phid;
00656 if(angd < dang) dang = angd;
00657
00658 }
00659
00660
00661
00662
00663
00664
00665
00666 dthe = dthe * 180 / (CLHEP::pi);
00667 dphi = dphi * 180 / (CLHEP::pi);
00668 dang = dang * 180 / (CLHEP::pi);
00669 double eraw = emcTrk->energy();
00670 double phi = emcTrk->phi();
00671 double the = emcTrk->theta();
00672
00673 m_delphi[iphoton]=dphi;
00674 m_delthe[iphoton]=dthe;
00675 m_delang[iphoton]=dang;
00676 if(energy < m_energyThreshold) continue;
00677 if(getTime>m_gammathCut||getTime<m_gammatlCut)continue;
00678
00679 if(dang< m_gammaTrkCut) continue;
00680 iphoton++;
00681 iGam.push_back(i);
00682 if(iphoton>=40)return StatusCode::SUCCESS;
00683 }
00684
00685 int nGam = iGam.size();
00686 m_nGam=nGam;
00687
00688
00689 counter[3]++;
00690
00691 double egam_ext=0;
00692 double ex_gam=0;
00693 double ey_gam=0;
00694 double ez_gam=0;
00695 double et_gam=0;
00696 double e_gam=0;
00697 for(int i = 0; i < m_nGam; i++) {
00698 EvtRecTrackIterator itTrk = evtRecTrkCol->begin()+ iGam[i];
00699 if(!(*itTrk)->isEmcShowerValid()) continue;
00700 RecEmcShower* emcTrk = (*itTrk)->emcShower();
00701 double eraw = emcTrk->energy();
00702 double phi = emcTrk->phi();
00703 double the = emcTrk->theta();
00704 HepLorentzVector ptrk;
00705 ex_gam+=eraw*sin(the)*cos(phi);
00706 ey_gam+=eraw*sin(the)*sin(phi);
00707 ez_gam+=eraw*cos(the);
00708 et_gam+=eraw*sin(the);
00709 e_gam+=eraw ;
00710 if(eraw>=egam_ext)
00711 {
00712 egam_ext=eraw;
00713 }
00714
00715 }
00716
00717
00718
00719
00720
00721 double px_had=0;
00722 double py_had=0;
00723 double pz_had=0;
00724 double pt_had=0;
00725 double p_had=0;
00726 double e_had=0;
00727
00728
00729
00730 int ii ;
00731 m_e_emc[0]=-0.1;
00732 m_e_emc[1]=-0.1;
00733 for(int i = 0; i < m_ngch; i++ ){
00734
00735 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
00736
00737 if(!(*itTrk)->isMdcTrackValid()) continue;
00738 if(!(*itTrk)->isMdcKalTrackValid()) continue;
00739
00740 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
00741 RecMdcKalTrack *mdcKalTrk = (*itTrk)->mdcKalTrack();
00742 ii=i;
00743
00744
00745 m_charge[ii] = mdcTrk->charge();
00746 m_vx0[ii] = mdcTrk->x();
00747 m_vy0[ii] = mdcTrk->y();
00748 m_vz0[ii] = mdcTrk->z();
00749
00750
00751 m_px[ii] = mdcTrk->px();
00752 m_py[ii] = mdcTrk->py();
00753 m_pz[ii] = mdcTrk->pz();
00754 m_p[ii] = mdcTrk->p();
00755
00756
00757 mdcKalTrk->setPidType(RecMdcKalTrack::electron);
00758
00759
00762 m_kal_vx0[ii] = mdcKalTrk->x();
00763 m_kal_vy0[ii] = mdcKalTrk->y();
00764 m_kal_vz0[ii] = mdcKalTrk->z();
00765
00766
00767 m_kal_px[ii] = mdcKalTrk->px();
00768 m_kal_py[ii] = mdcKalTrk->py();
00769 m_kal_pz[ii] = mdcKalTrk->pz();
00770 m_kal_p[ii] = mdcKalTrk->p();
00771
00772
00773 px_had+=mdcKalTrk->px();
00774 py_had+=mdcKalTrk->py();
00775 pz_had+=mdcKalTrk->pz();
00776 pt_had+=mdcKalTrk->pxy();
00777 p_had+=mdcKalTrk->p();
00778 e_had+=sqrt(mdcKalTrk->p()*mdcKalTrk->p()+mdcKalTrk->mass()*mdcKalTrk->mass());
00779
00780 double ptrk = mdcKalTrk->p() ;
00781
00782
00783 if((*itTrk)->isMdcDedxValid()) {
00784
00785 RecMdcDedx* dedxTrk = (*itTrk)->mdcDedx();
00786 m_probPH[ii]= dedxTrk->probPH();
00787 m_normPH[ii]= dedxTrk->normPH();
00788
00789 m_chie[ii] = dedxTrk->chiE();
00790 m_chimu[ii] = dedxTrk->chiMu();
00791 m_chipi[ii] = dedxTrk->chiPi();
00792 m_chik[ii] = dedxTrk->chiK();
00793 m_chip[ii] = dedxTrk->chiP();
00794 m_ghit[ii] = dedxTrk->numGoodHits();
00795 m_thit[ii] = dedxTrk->numTotalHits();
00796 }
00797
00798 if((*itTrk)->isEmcShowerValid()) {
00799
00800 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00801 m_e_emc[ii] = emcTrk->energy();
00802 m_phi_emc[ii] = emcTrk->phi();
00803 m_theta_emc[ii] = emcTrk->theta();
00804 }
00805
00806
00807 if((*itTrk)->isMucTrackValid()){
00808
00809 RecMucTrack* mucTrk = (*itTrk)->mucTrack() ;
00810 m_nhit_muc[ii] = mucTrk->numHits() ;
00811 m_nlay_muc[ii] = mucTrk->numLayers() ;
00812
00813 }
00814
00815 if((*itTrk)->isTofTrackValid()) {
00816
00817 SmartRefVector<RecTofTrack> tofTrkCol = (*itTrk)->tofTrack();
00818
00819 SmartRefVector<RecTofTrack>::iterator iter_tof = tofTrkCol.begin();
00820
00821 for(;iter_tof != tofTrkCol.end(); iter_tof++ ) {
00822 TofHitStatus *status = new TofHitStatus;
00823 status->setStatus((*iter_tof)->status());
00824
00825 if(!(status->is_barrel())){
00826 if( (status->is_cluster()) ) m_t_etof[ii] = (*iter_tof)->tof();
00827 if( !(status->is_counter()) ){if(status) delete status; continue; }
00828 if( status->layer()!=0 ) {if(status) delete status;continue;}
00829 double path=(*iter_tof)->path();
00830 double tof = (*iter_tof)->tof();
00831 double ph = (*iter_tof)->ph();
00832 double rhit = (*iter_tof)->zrhit();
00833 double qual = 0.0 + (*iter_tof)->quality();
00834 double cntr = 0.0 + (*iter_tof)->tofID();
00835 double texp[5];
00836 for(int j = 0; j < 5; j++) {
00837 double gb = ptrk/xmass[j];
00838 double beta = gb/sqrt(1+gb*gb);
00839 texp[j] = path /beta/velc;
00840 }
00841
00842 m_qual_etof[ii] = qual;
00843 m_tof_etof[ii] = tof ;
00844 }
00845 else {
00846 if( (status->is_cluster()) ) m_t_btof[ii] = (*iter_tof)->tof();
00847 if( !(status->is_counter()) ){if(status) delete status; continue;}
00848 if(status->layer()==1){
00849 double path=(*iter_tof)->path();
00850 double tof = (*iter_tof)->tof();
00851 double ph = (*iter_tof)->ph();
00852 double rhit = (*iter_tof)->zrhit();
00853 double qual = 0.0 + (*iter_tof)->quality();
00854 double cntr = 0.0 + (*iter_tof)->tofID();
00855 double texp[5];
00856 for(int j = 0; j < 5; j++) {
00857 double gb = ptrk/xmass[j];
00858 double beta = gb/sqrt(1+gb*gb);
00859 texp[j] = path /beta/velc;
00860 }
00861
00862 m_qual_btof1[ii] = qual;
00863 m_tof_btof1[ii] = tof ;
00864 }
00865
00866 if(status->layer()==2){
00867 double path=(*iter_tof)->path();
00868 double tof = (*iter_tof)->tof();
00869 double ph = (*iter_tof)->ph();
00870 double rhit = (*iter_tof)->zrhit();
00871 double qual = 0.0 + (*iter_tof)->quality();
00872 double cntr = 0.0 + (*iter_tof)->tofID();
00873 double texp[5];
00874 for(int j = 0; j < 5; j++) {
00875 double gb = ptrk/xmass[j];
00876 double beta = gb/sqrt(1+gb*gb);
00877 texp[j] = path /beta/velc;
00878 }
00879
00880 m_qual_btof2[ii] = qual;
00881 m_tof_btof2[ii] = tof ;
00882 }
00883 }
00884 if(status) delete status;
00885 }
00886 }
00887
00888 }
00889 counter[4]++;
00890
00891
00892
00893
00894 m_bhabhatag=0;
00895
00896 if(m_ngch != 2 || nCharge != 0 ) return StatusCode::SUCCESS;
00897 EvtRecTrackIterator itTrk1;
00898
00899 EvtRecTrackIterator itTrk2;
00900
00901 RecMdcKalTrack *mdcKalTrk1;
00902
00903 RecMdcKalTrack *mdcKalTrk2;
00904
00905 HepLorentzVector p41e,p42e,p4le;
00906 Hep3Vector p31e,p32e,p3le;
00907 HepLorentzVector p41m,p42m,p4lm;
00908 Hep3Vector p31m,p32m,p3lm;
00909 HepLorentzVector p41h,p42h,p4lh;
00910 Hep3Vector p31h,p32h,p3lh;
00911 WTrackParameter w1_ini;
00912 WTrackParameter w2_ini;
00913 int iip=-1;
00914 int iim=-1;
00915 for(int i = 0; i < m_ngch; i++ ){
00916 if(m_charge[i]>0)itTrk1= evtRecTrkCol->begin() + iGood[i];
00917 if(m_charge[i]<0) itTrk2= evtRecTrkCol->begin() + iGood[i];
00918 if(m_charge[i]>0) mdcKalTrk1 = (*itTrk1)->mdcKalTrack();
00919 if(m_charge[i]<0) mdcKalTrk2 = (*itTrk2)->mdcKalTrack();
00920 if(m_charge[i]>0)iip=i;
00921 if(m_charge[i]<0)iim=i;
00922
00923
00924 if(m_charge[i]>0) w1_ini=WTrackParameter (xmass[0],mdcKalTrk1->getZHelixE(),mdcKalTrk1->getZErrorE());
00925 if(m_charge[i]<0) w2_ini=WTrackParameter (xmass[0],mdcKalTrk2 ->getZHelixE(),mdcKalTrk2 ->getZErrorE());
00926 if(m_charge[i]>0) p41e =w1_ini.p();
00927 if(m_charge[i]<0) p42e =w2_ini.p();
00928 if(m_charge[i]>0) p41e.boost(u_cms);
00929 if(m_charge[i]<0) p42e.boost(u_cms);
00930 if(m_charge[i]>0) p31e = p41e.vect();
00931 if(m_charge[i]<0) p32e = p42e.vect();
00932
00933
00934
00935
00936 if(m_charge[i]>0){
00937 m_px_cms_ep=p41e.px();
00938 m_py_cms_ep=p41e.py();
00939 m_pz_cms_ep=p41e.pz();
00940 m_e_cms_ep=p41e.e();
00941 }
00942 if(m_charge[i]<0){
00943 m_px_cms_em=p42e.px();
00944 m_py_cms_em=p42e.py();
00945 m_pz_cms_em=p42e.pz();
00946 m_e_cms_em=p42e.e();
00947 }
00948
00949
00950 }
00951
00952 double e01=(iip!=-1)?m_e_emc[iip]:0;
00953 double e02=(iim!=-1)?m_e_emc[iim]:0;
00954 int ilarge=( e01 > e02 ) ?iip:iim;
00955 p4le=( e01 > e02 ) ?p41e:p42e;
00956 p4lm=( e01 > e02 ) ?p41m:p42m;
00957
00958 p3le=( e01 > e02 ) ?p31e:p32e;
00959 p3lm=( e01 > e02 ) ?p31m:p32m;
00960
00961 double acolle= 180.-p31e.angle(p32e)* 180.0 / CLHEP::pi;
00962 double acople= 180.- (p31e.perpPart()).angle(p32e.perpPart ())* 180.0 / CLHEP::pi;
00963 double poeb1e=p41e.rho()/beamEnergy;
00964 double poeb2e=p42e.rho()/beamEnergy;
00965 double poeble=p4le.rho()/beamEnergy ;
00966
00967
00968 double eemc1=m_e_emc[iip];
00969 double eemc2=m_e_emc[iim];
00970
00971
00972 double ex1=m_kal_vx0[iip];
00973 double ey1=m_kal_vy0[iip];
00974 double ez1=m_kal_vz0[iip];
00975 double epx1=m_kal_px[iip];
00976 double epy1=m_kal_py[iip];
00977 double epz1=m_kal_pz[iip];
00978 double epp1=m_kal_p[iip];
00979 double ex2=m_kal_vx0[iim];
00980 double ey2=m_kal_vy0[iim];
00981 double ez2=m_kal_vz0[iim];
00982 double epx2=m_kal_px[iim];
00983 double epy2=m_kal_py[iim];
00984 double epz2=m_kal_pz[iim];
00985 double epp2=m_kal_p[iim];
00986
00987 double pidchidedx1=m_pidchiDedx[iip];
00988 double pidchitof11= m_pidchiTof1[iip];
00989 double pidchitof21= m_pidchiTof2[iip];
00990 double pidchidedx2=m_pidchiDedx[iim];
00991 double pidchitof12= m_pidchiTof1[iim];
00992 double pidchitof22= m_pidchiTof2[iim];
00993
00994 double eoeb1=m_e_emc[iip]/beamEnergy;
00995 double eoeb2=m_e_emc[iim]/beamEnergy;
00996
00997 double eope1=0;
00998 if(p41e.rho()>0)eope1=m_e_emc[iip]/p41e.rho();
00999 double eope2=0;
01000 if(p42e.rho()>0)eope2=m_e_emc[iim]/p42e.rho();
01001
01002
01003
01004 double exoeb1= m_e_emc[iip]*sin(m_theta_emc[iip])*cos(m_phi_emc[iip])/beamEnergy;
01005 double eyoeb1= m_e_emc[iip]*sin(m_theta_emc[iip])*sin(m_phi_emc[iip])/beamEnergy;
01006 double ezoeb1=m_e_emc[iip]*cos(m_theta_emc[iip])/beamEnergy;
01007 double etoeb1=m_e_emc[iip]*sin(m_theta_emc[iip])/beamEnergy;
01008
01009 double exoeb2= m_e_emc[iim]*sin(m_theta_emc[iim])*cos(m_phi_emc[iim])/beamEnergy;
01010 double eyoeb2= m_e_emc[iim]*sin(m_theta_emc[iim])*sin(m_phi_emc[iim])/beamEnergy;
01011 double ezoeb2=m_e_emc[iim]*cos(m_theta_emc[iim])/beamEnergy;
01012 double etoeb2=m_e_emc[iim]*sin(m_theta_emc[iim])/beamEnergy;
01013
01014 double eoebl=m_e_emc[ilarge]/beamEnergy;
01015
01016 double eopl=0;
01017 if(p4le.rho()>0)eopl=m_e_emc[ilarge]/p4le.rho();
01018
01019 double exoebl= m_e_emc[ilarge]*sin(m_theta_emc[ilarge])*cos(m_phi_emc[ilarge])/beamEnergy;
01020 double eyoebl= m_e_emc[ilarge]*sin(m_theta_emc[ilarge])*sin(m_phi_emc[ilarge])/beamEnergy;
01021 double ezoebl=m_e_emc[ilarge]*cos(m_theta_emc[ilarge])/beamEnergy;
01022 double etoebl=m_e_emc[ilarge]*sin(m_theta_emc[ilarge])/beamEnergy;
01023
01024 int mucinfo1=(m_nhit_muc[iip]>=2&&m_nlay_muc[iip]>=2 ) ? 1 : 0;
01025 int mucinfo2=(m_nhit_muc[iim]>=2&&m_nlay_muc[iim]>=2) ? 1 : 0;
01026 int mucinfol=(m_nhit_muc[ilarge]>=2&&m_nlay_muc[ilarge]>=2) ? 1 : 0;
01027 int pidel=( e01 > e02 ) ? m_nep : m_nem;
01028 int pidmul=( e01 > e02 ) ? m_nmup : m_nmum;
01029 double deltatof=0;
01030
01031 if(m_t_btof[iip]*m_t_btof[iim]!=0) deltatof+=fabs(m_t_btof[iip]-m_t_btof[iim]);
01032 if(m_t_etof[iip]*m_t_etof[iim]!=0) deltatof+=fabs(m_t_etof[iip]-m_t_etof[iim]);
01033
01034
01035
01036
01037
01038
01039
01040 if (acolle<m_acoll_e_cut) m_bhabhatag+=1;
01041 if (acople<m_acopl_e_cut)m_bhabhatag+=10;
01042 if (fabs(m_ecms-mpsi2s)<0.001){
01043 if (sqrt((eoeb1-1)*(eoeb1-1)+(eoeb2-1)*(eoeb2-1))<m_tetotal_e_cut)m_bhabhatag+=100;
01044 }
01045 else{
01046 if ((fabs(poeb1e-1)<m_poeb_e_cut)&&(fabs(poeb2e-1)<m_poeb_e_cut))m_bhabhatag+=100;
01047 }
01048 if (!m_useTOF||(deltatof<m_dtof_e_cut))m_bhabhatag+=1000;
01049 if (poeb1e>m_tpoeb_e_cut||poeb2e>m_tpoeb_e_cut||(sqrt((poeb1e-1)*(poeb1e-1)+(poeb2e-1)*(poeb2e-1))<m_tptotal_e_cut))m_bhabhatag+=10000;
01050 if (!m_useMUC||(mucinfo1==0||mucinfo2==0))m_bhabhatag+=100000;
01051 if (!m_usePID||(m_nep==1&&m_nem==1))m_bhabhatag+=1000000;
01052
01053
01054
01055
01056
01057 m_acoll=acolle;
01058 m_acopl=acople;
01059 m_poeb1=poeb1e;
01060 m_poeb2=poeb2e;
01061 m_eop1=eope1;
01062 m_eop2=eope2;
01063 m_cos_ep=p41e.cosTheta ();
01064 m_cos_em=p42e.cosTheta ();
01065 m_mass_ee=(p41e+p42e).m();
01066 m_deltatof=deltatof;
01067
01068 m_eoeb1=eoeb1;
01069 m_eoeb2=eoeb2;
01070
01071 m_etoeb1=etoeb1;
01072 m_etoeb2=etoeb2;
01073 m_mucinfo1=mucinfo1;
01074 m_mucinfo2=mucinfo2;
01075
01076 if(m_bhabhatag==1111111){
01077 nbhabha++;
01078 if(m_writentuple)m_tuple1 -> write();
01080
01081
01082
01083
01084
01085 (*itTrk1)->setPartId(1);
01086 (*itTrk2)->setPartId(1);
01087
01088
01089
01090
01091
01092 (*itTrk1)->setQuality(0);
01093 (*itTrk2)->setQuality(0);
01094
01095
01096
01097
01098 setFilterPassed(true);
01100 m_ee_mass->Fill((p41e+p42e).m());
01101 m_ee_acoll->Fill(acolle);
01102 m_ee_eop_ep->Fill(eope1);
01103 m_ee_eop_em->Fill(eope2);
01104 m_ee_costheta_ep->Fill(p41e.cosTheta ());
01105 m_ee_costheta_em->Fill(p42e.cosTheta ());
01106 m_ee_phi_ep->Fill(p41e.phi ());
01107 m_ee_phi_em->Fill(p42e.phi ());
01108 m_ee_nneu->Fill(m_nGam);
01109
01110
01111
01112 m_ee_eemc_ep->Fill(eemc1);
01113 m_ee_eemc_em->Fill(eemc2);
01114
01115 m_ee_x_ep->Fill(ex1);
01116 m_ee_y_ep->Fill(ey1);
01117 m_ee_z_ep->Fill(ez1);
01118
01119 m_ee_x_em->Fill(ex2);
01120 m_ee_y_em->Fill(ey2);
01121 m_ee_z_em->Fill(ez2);
01122
01123
01124 m_ee_px_ep->Fill(epx1);
01125 m_ee_py_ep->Fill(epy1);
01126 m_ee_pz_ep->Fill(epz1);
01127 m_ee_p_ep->Fill(epp1);
01128
01129 m_ee_px_em->Fill(epx2);
01130 m_ee_py_em->Fill(epy2);
01131 m_ee_pz_em->Fill(epz2);
01132 m_ee_p_em->Fill(epp2);
01133
01134 m_ee_deltatof->Fill(deltatof);
01135
01136 m_ee_pidchidedx_ep->Fill(pidchidedx1);
01137 m_ee_pidchidedx_em->Fill(pidchidedx2);
01138 m_ee_pidchitof1_ep->Fill(pidchitof11);
01139 m_ee_pidchitof1_em->Fill(pidchitof12);
01140 m_ee_pidchitof2_ep->Fill(pidchitof21);
01141 m_ee_pidchitof2_em->Fill(pidchitof22);
01142
01143
01144
01145
01146 }
01147
01148
01149
01150
01151
01152
01153
01154
01155 return StatusCode::SUCCESS;
01156
01157
01158
01159
01160
01161
01162 }
01163
01164
01165 StatusCode DQASelBhabha::finalize() {
01166
01167 MsgStream log(msgSvc(), name());
01168 log << MSG::INFO << "in finalize()" << endmsg;
01169 return StatusCode::SUCCESS;
01170 }
01171
01172