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