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