diff --git a/ForwardDetectors/ALFA/ALFA_BeamTransport/src/ALFA_BeamTransport.cxx b/ForwardDetectors/ALFA/ALFA_BeamTransport/src/ALFA_BeamTransport.cxx index 8e8b7426c35f4267a5b195e3f443b689c3420497..3664258798e6dcc9d1efeab50413ce4f812d720a 100644 --- a/ForwardDetectors/ALFA/ALFA_BeamTransport/src/ALFA_BeamTransport.cxx +++ b/ForwardDetectors/ALFA/ALFA_BeamTransport/src/ALFA_BeamTransport.cxx @@ -255,8 +255,8 @@ int ALFA_BeamTransport::AddHepMCData(HepMC::GenEvent* evt) HepMC::FourVector MomentumVectorRP1 = HepMC::FourVector(m_MomRP1.x(),m_MomRP1.y(),m_MomRP1.z(),m_EnergyRP1); - HepMC::GenVertex* VertexRP1 = new HepMC::GenVertex(PositionVectorRP1); - HepMC::GenParticle* ParticleRP1 = new HepMC::GenParticle(MomentumVectorRP1,2212,1); + HepMC::GenVertexPtr VertexRP1 = HepMC::newGenVertexPtr(PositionVectorRP1); + HepMC::GenParticlePtr ParticleRP1 = HepMC::newGenParticlePtr(MomentumVectorRP1,2212,1); //Add particle to vertex VertexRP1->add_particle_out(ParticleRP1); //add new vertex to HepMC event record @@ -269,8 +269,8 @@ int ALFA_BeamTransport::AddHepMCData(HepMC::GenEvent* evt) HepMC::FourVector MomentumVectorRP3 = HepMC::FourVector(m_MomRP3.x(),m_MomRP3.y(),m_MomRP3.z(),m_EnergyRP3); - HepMC::GenVertex* VertexRP3 = new HepMC::GenVertex(PositionVectorRP3); - HepMC::GenParticle* ParticleRP3 = new HepMC::GenParticle(MomentumVectorRP3,2212,1); + HepMC::GenVertexPtr VertexRP3 = HepMC::newGenVertexPtr(PositionVectorRP3); + HepMC::GenParticlePtr ParticleRP3 = HepMC::newGenParticlePtr(MomentumVectorRP3,2212,1); VertexRP3->add_particle_out(ParticleRP3); @@ -317,8 +317,8 @@ int ALFA_BeamTransport::DoBeamTracking(int evt_number) } int ALFA_BeamTransport::TransportSelectedParticle(HepMC::GenEvent* evt, int evt_number){ - HepMC::GenParticle* p1=0; - HepMC::GenParticle* p2=0; + HepMC::GenParticlePtr p1{nullptr}; + HepMC::GenParticlePtr p2{nullptr}; std::vector<FPTracker::Point> PosAtRP1; @@ -336,41 +336,41 @@ int ALFA_BeamTransport::TransportSelectedParticle(HepMC::GenEvent* evt, int evt_ //First we have to select the final state particles from the MC - for ( HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p ) + for ( auto p: *evt) { //Simple Eta Pt cut to remove particles from BeamTransportation which have no chance to reach RP plane - mom=std::sqrt(std::pow((*p)->momentum().px(),2)+std::pow((*p)->momentum().py(),2)+std::pow((*p)->momentum().pz(),2)); - theta=std::acos(std::abs((*p)->momentum().pz())/mom); + mom=std::sqrt(std::pow(p->momentum().px(),2)+std::pow(p->momentum().py(),2)+std::pow(p->momentum().pz(),2)); + theta=std::acos(std::abs(p->momentum().pz())/mom); eta=-std::log(std::tan(theta/2)); std::cout<<"eta "<<eta<<std::endl; std::cout<<"DeptaP/p "<<std::abs((mom)/m_FPConfig.pbeam0)<<std::endl; - if( ((*p)->status() == 1) && (!(*p)->end_vertex()) ) {//TODO What is end_vertex()??? + if( (p->status() == 1) && (!p->end_vertex()) ) {//TODO What is end_vertex()??? //Change the status code from Pythia (1) to 201 //added 120124 - (*p)->set_status(201); + p->set_status(201); - int pid = (*p)->pdg_id(); + int pid = p->pdg_id(); if(eta>m_EtaCut && 1-std::abs(mom/m_FPConfig.pbeam0) < m_XiCut){ //save a copy of the particles which passed the cut - HepMC::FourVector Position = (*p)->production_vertex()->position(); + HepMC::FourVector Position = p->production_vertex()->position(); - HepMC::FourVector Momentum = (*p)->momentum(); + HepMC::FourVector Momentum = p->momentum(); - HepMC::GenVertex* Vertex = new HepMC::GenVertex(Position); //copy of the vertex - HepMC::GenParticle* Particle = new HepMC::GenParticle(Momentum,pid,202); + HepMC::GenVertexPtr Vertex = HepMC::newGenVertexPtr(Position); //copy of the vertex + HepMC::GenParticlePtr Particle = HepMC::newGenParticlePtr(Momentum,pid,202); Vertex->add_particle_out(Particle); evt->add_vertex(Vertex); //select direction of particle - if((*p)->momentum().pz()>0. && pid == 2212){//Beam1 TODO Tracking only works for protons!!!! - p1=(*p); + if(p->momentum().pz()>0. && pid == 2212){//Beam1 TODO Tracking only works for protons!!!! + p1=p; //now we want to track the final particle if it's a protons //Positions are given in mm FPTracker needs them in meter m_Particle1=FPTracker::Particle(p1->production_vertex()->position().x()/1000., @@ -400,8 +400,8 @@ int ALFA_BeamTransport::TransportSelectedParticle(HepMC::GenEvent* evt, int evt_ } } - else if((*p)->momentum().pz()<0. && pid == 2212){//beam 2 - p2=(*p); + else if(p->momentum().pz()<0. && pid == 2212){//beam 2 + p2=p; m_Particle2=FPTracker::Particle(p2->production_vertex()->position().x()/1000., p2->production_vertex()->position().y()/1000., @@ -461,8 +461,8 @@ int ALFA_BeamTransport::TransportSelectedParticle(HepMC::GenEvent* evt, int evt_ HepMC::FourVector MomentumVectorRP1 = HepMC::FourVector(MomAtPR1.at(i).x(),MomAtPR1.at(i).y(),MomAtPR1.at(i).z(),EnergyRP1.at(i)); - HepMC::GenVertex* VertexRP1 = new HepMC::GenVertex(PositionVectorRP1); - HepMC::GenParticle* ParticleRP1 = new HepMC::GenParticle(MomentumVectorRP1,2212,1); //save the transported particle with status code 1 (added 120124) preview was 201 + HepMC::GenVertexPtr VertexRP1 = HepMC::newGenVertexPtr(PositionVectorRP1); + HepMC::GenParticlePtr ParticleRP1 = HepMC::newGenParticlePtr(MomentumVectorRP1,2212,1); //save the transported particle with status code 1 (added 120124) preview was 201 //Add particle to vertex VertexRP1->add_particle_out(ParticleRP1); //add new vertex to HepMC event record @@ -481,8 +481,8 @@ int ALFA_BeamTransport::TransportSelectedParticle(HepMC::GenEvent* evt, int evt_ HepMC::FourVector MomentumVectorRP3 = HepMC::FourVector(MomAtPR3.at(i).x(),MomAtPR3.at(i).y(),MomAtPR3.at(i).z(),EnergyRP3.at(i)); - HepMC::GenVertex* VertexRP3 = new HepMC::GenVertex(PositionVectorRP3); - HepMC::GenParticle* ParticleRP3 = new HepMC::GenParticle(MomentumVectorRP3,2212,1);//save the transported particle with status code 1 (added 120124) preview was 201 + HepMC::GenVertexPtr VertexRP3 = HepMC::newGenVertexPtr(PositionVectorRP3); + HepMC::GenParticlePtr ParticleRP3 = HepMC::newGenParticlePtr(MomentumVectorRP3,2212,1);//save the transported particle with status code 1 (added 120124) preview was 201 VertexRP3->add_particle_out(ParticleRP3); diff --git a/ForwardDetectors/ALFA/ALFA_Ntuple/src/ALFA_Ntuple.cxx b/ForwardDetectors/ALFA/ALFA_Ntuple/src/ALFA_Ntuple.cxx index f092ac8b5de910cb3cb623d9206ded3d322163aa..570c789455cd35cb68a22ce144bc33607d9632d8 100644 --- a/ForwardDetectors/ALFA/ALFA_Ntuple/src/ALFA_Ntuple.cxx +++ b/ForwardDetectors/ALFA/ALFA_Ntuple/src/ALFA_Ntuple.cxx @@ -856,10 +856,10 @@ StatusCode ALFA_Ntuple::TruthInfo() int vtx_counter[8]; memset(&vtx_counter, 0, sizeof(vtx_counter)); - HepMC::GenParticle* pi1=0; - HepMC::GenParticle* pi2=0; - HepMC::GenParticle* p1=0; - HepMC::GenParticle* p2=0; + HepMC::ConstGenParticlePtr pi1{nullptr}; + HepMC::ConstGenParticlePtr pi2{nullptr}; + HepMC::ConstGenParticlePtr p1{nullptr}; + HepMC::ConstGenParticlePtr p2{nullptr}; int pint = 0; int pcount = 0; @@ -887,37 +887,40 @@ StatusCode ALFA_Ntuple::TruthInfo() break; } - // comment to fix coverity 29013 -// if (coll_counter == 0 ) -// { -// ATH_MSG_DEBUG("no collection for the event "); -// return sc; -// } +#ifdef HEPMC3 + auto begGenVtxItr = (**mcTruBeg).vertices().begin(); + auto endGenVtxItr = (**mcTruBeg).vertices().end(); +#else HepMC::GenEvent::vertex_const_iterator begGenVtxItr = (**mcTruBeg).vertices_begin(); HepMC::GenEvent::vertex_const_iterator endGenVtxItr = (**mcTruBeg).vertices_end(); +#endif //loop over verteces belonging to one event for(;begGenVtxItr!=endGenVtxItr;++begGenVtxItr) { -// std::cout << "coll_counter-1 = " << coll_counter-1 << endl; vtx_counter[coll_counter-1]++; -// std::cout << " vtx_counter - 1 = " << vtx_counter[coll_counter-1]-1 << endl; ATH_MSG_DEBUG(" * vertex no: " << vtx_counter[coll_counter-1]); ATH_MSG_DEBUG(" * position x = " << (**begGenVtxItr).position().x() << ", y = " << (**begGenVtxItr).position().y()<< ", z =" << (**begGenVtxItr).position().z()); +#ifdef HEPMC3 + auto children=(*begGenVtxItr)->particles_out(); + auto parents=(*begGenVtxItr)->particles_in(); + children.insert( children.end(), parents.begin(), parents.end() ); + auto child = children.begin(); + auto child_end = children.end(); +#else HepMC::GenVertex::particle_iterator child; child = (*begGenVtxItr)->particles_begin(HepMC::family); HepMC::GenVertex::particle_iterator child_end; child_end = (*begGenVtxItr)->particles_end(HepMC::family); +#endif for(; child != child_end; ++child) { -// if ((*child == genEvt->beam_particles().first) || (*child == genEvt->beam_particles().second)) -// { px = (*child)->momentum().px(); py = (*child)->momentum().py(); pz = (*child)->momentum().pz(); @@ -938,7 +941,6 @@ StatusCode ALFA_Ntuple::TruthInfo() { pint++; -// std::cout << "pint = " << pint << endl; if(pz > 0) { @@ -954,20 +956,16 @@ StatusCode ALFA_Ntuple::TruthInfo() { pi2=(*child); } -// if(pint > 2){ATH_MSG_DEBUG("Strange: More than two incoming protons in this event!");} } } -// std::cout << "pip1" << endl; // outgoing protons at the interaction point if( (*child)->status() == 201 && (pcount < 2)) -// if( (*child)->status() == 212 && (pcount < 2)) { if( (*child)->pdg_id() == 2212) { pcount++; -// std::cout << "pip2" << endl; if(pz > 0) { @@ -987,7 +985,6 @@ StatusCode ALFA_Ntuple::TruthInfo() m_fvtx_beam2_f[2] = (**begGenVtxItr).position().y(); m_fvtx_beam2_f[3] = (**begGenVtxItr).position().z(); } -// if(pcount > 2){ATH_MSG_DEBUG("Strange: More than two outcoming protons in this event (elastic scaterring)!");} } } @@ -995,7 +992,6 @@ StatusCode ALFA_Ntuple::TruthInfo() { if( (*child)->pdg_id() == 2212){ -// std::cout << "pip" << endl; if(pz > 0) { @@ -1021,11 +1017,8 @@ StatusCode ALFA_Ntuple::TruthInfo() m_fp_C[2] = py; m_fp_C[3] = pz; } -// if(pint > 2){ATH_MSG_DEBUG("Strange: More than two incoming protons in this event!");} } } -// std::cout << "pip3" << endl; -// } } } @@ -1035,37 +1028,14 @@ StatusCode ALFA_Ntuple::TruthInfo() { m_bVtxKinFillFlag = true; -// HepMC::FourVector pv1 = (pi1->momentum()); -// HepMC::FourVector pv2 = (pi2->momentum()); -// HepMC::FourVector pv3 = (p1->momentum()); -// HepMC::FourVector pv4 = (p2->momentum()); -// CLHEP::HepLorentzVector hp1(pv1.px(),pv1.py(),pv1.pz(),pv1.e()); -// CLHEP::HepLorentzVector hp2(pv2.px(),pv2.py(),pv2.pz(),pv2.e()); -// CLHEP::HepLorentzVector hp3(pv3.px(),pv3.py(),pv3.pz(),pv3.e()); -// CLHEP::HepLorentzVector hp4(pv4.px(),pv4.py(),pv4.pz(),pv4.e()); -// m_ft_13 = (hp1-hp3).m2(); -// m_ft_24 = (hp2-hp4).m2(); -// ATH_MSG_DEBUG(" ******************************************************* "); -// ATH_MSG_DEBUG(" "); -// ATH_MSG_DEBUG(" t_13 = " << m_ft_13); -// ATH_MSG_DEBUG(" t_24 = " << m_ft_24); -// ATH_MSG_DEBUG(" ***************** "); -// m_fp_beam1_i[0] = pi1->momentum().e(); -// m_fp_beam1_i[1] = pi1->momentum().px(); -// m_fp_beam1_i[2] = pi1->momentum().py(); -// m_fp_beam1_i[3] = pi1->momentum().pz(); ATH_MSG_DEBUG("initial particle 1: px = " << m_fp_beam1_i[1] << ", py = " << m_fp_beam1_i[2] << ", pz = " << m_fp_beam1_i[3] << ", E = " << m_fp_beam1_i[0]); ATH_MSG_DEBUG(" ** "); -// m_fp_beam2_i[0] = pi2->momentum().e(); -// m_fp_beam2_i[1] = pi2->momentum().px(); -// m_fp_beam2_i[2] = pi2->momentum().py(); -// m_fp_beam2_i[3] = pi2->momentum().pz(); ATH_MSG_DEBUG("initial particle 2: px = " << m_fp_beam2_i[1] << ", py = " << m_fp_beam2_i[2] << ", pz = " << m_fp_beam2_i[3] << ", E = " << m_fp_beam2_i[0]); ATH_MSG_DEBUG(" ** "); @@ -1348,9 +1318,7 @@ StatusCode ALFA_Ntuple::GetLocRecData() } -// cout << "iAlgoID, iRPot, iNumTrackPot, iNumTrackPotMD << " << iAlgoID << ", " << iRPot << ", " << iNumTrackPot[iRPot] << ", " << mapNumTrackPotMD[iAlgoID][iRPot] << endl; -// cout << "iRPot, m_iDetector, fRecPosX = " << iRPot << ", " << m_iDetector << ", " << fRecPosX << endl; } @@ -1382,7 +1350,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() fOverY = (*mcColODBeg)->getOverY(); iNumY = (*mcColODBeg)->getNumY(); -// cout << "iPot, algo, yOD = " << iRPot << ", " << iAlgoID << ", " << fRecPosY << endl; if (mapNumTrackPotOD.find(iAlgoID)==mapNumTrackPotOD.end()) { @@ -1429,7 +1396,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() mapNumTrackPotOD[iAlgoID][iRPot]++; } -// cout << "iAlgoID, iRPot, iNumTrackPot, iNumTrackPotOD << " << iAlgoID << ", " << iRPot << ", " << iNumTrackPot[iRPot] << ", " << mapNumTrackPotOD[iAlgoID][iRPot] << endl; } Int_t iNumTrackPotMax = 1; @@ -1448,7 +1414,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() { for (int iPot=0; iPot<RPOTSCNT; iPot++) { -// cout << "first, second = " << it->first << ", " << it->second[iPot] << endl; if (mapNumTrackPotMax.find(it->first)==mapNumTrackPotMax.end()) mapNumTrackPotMax.insert(std::pair<int, int>(it->first, 0)); if (it->second[iPot] > mapNumTrackPotMax[it->first]) mapNumTrackPotMax[it->first] = it->second[iPot]; @@ -1458,7 +1423,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() std::map<int, int >::iterator itNum; for (itNum=mapNumTrackPotMax.begin(); itNum!=mapNumTrackPotMax.end(); itNum++) { -// cout << "first, second = " << itNum->first << ", " << itNum->second << endl; m_MapAlgoTreeMD[itNum->first].iNumTrack = itNum->second; } @@ -1469,7 +1433,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() { for (int iPot=0; iPot<RPOTSCNT; iPot++) { -// cout << "first, second = " << it->first << ", " << it->second[iPot] << endl; if (mapNumTrackPotMax.find(it->first)==mapNumTrackPotMax.end()) mapNumTrackPotMax.insert(std::pair<int, int>(it->first, 0)); if (it->second[iPot] > mapNumTrackPotMax[it->first]) mapNumTrackPotMax[it->first] = it->second[iPot]; @@ -1478,7 +1441,6 @@ StatusCode ALFA_Ntuple::GetLocRecData() for (itNum=mapNumTrackPotMax.begin(); itNum!=mapNumTrackPotMax.end(); itNum++) { -// cout << "first, second = " << itNum->first << ", " << itNum->second << endl; m_MapAlgoTreeOD[itNum->first].iNumTrack = itNum->second; } @@ -1579,7 +1541,6 @@ StatusCode ALFA_Ntuple::GetLocRecCorrData() } -// cout << "iRPot, iNumTrackPot, fRecPosX = " << iRPot << ", " << iNumTrackPot[iRPot]-1 << ", " << fRecPosX << endl; } sc = evtStore()->retrieve(pLocRecCorrODCol, m_strLocRecCorrODCollectionName); @@ -1613,7 +1574,6 @@ StatusCode ALFA_Ntuple::GetLocRecCorrData() fxBeam = iSign*22; fyBeam = (*mcColODBeg)->getYpositionBeam(); -// cout << "LocRecCorr = " << fxLHC << ", " << fyLHC << ", " << fzLHC << endl; if (mapNumTrackPot.find(iAlgoID)==mapNumTrackPot.end()) { @@ -1695,7 +1655,6 @@ StatusCode ALFA_Ntuple::GetGloRecData() iNumGloTrack++; -// cout << "NumGloTrack, Arm, x, y, x_slope, y_slope = " << iNumGloTrack << ", " << iArm << ", " << fxPos << ", " << fyPos << ", " << fxSlope << ", " << fySlope << endl; } m_iNumGloTrack = iNumGloTrack; @@ -1936,10 +1895,6 @@ StatusCode ALFA_Ntuple::COOLUpdate(IOVSVC_CALLBACK_ARGS_P(/*I*/, keys)) // to access individual elements by name, use e.g. // float var1=(*atrlist)["T04"].data<float>(); // to get the value of a float column called T04 into var1 -// std::ostringstream atrstring; -// atrlist->print(atrstring); -// cout<< "Values for folder /TDAQ/ALFA" << endl; -// cout << atrstring.str() << endl; m_BPMALFA.bpmr_r_x_pos = (*atrlist)["bpmr_r_x_pos"].data<float>(); m_BPMALFA.bpmr_r_y_pos = (*atrlist)["bpmr_r_y_pos"].data<float>(); diff --git a/ForwardDetectors/ForwardTransportFast/src/ForwardTransportFast.cxx b/ForwardDetectors/ForwardTransportFast/src/ForwardTransportFast.cxx index 1e2f28a0e63f2ef01a36ae53aa42e4995a3bd220..5ea9de9cb5d66755194d15395ad4ec975a5ac2e4 100644 --- a/ForwardDetectors/ForwardTransportFast/src/ForwardTransportFast.cxx +++ b/ForwardDetectors/ForwardTransportFast/src/ForwardTransportFast.cxx @@ -116,8 +116,8 @@ StatusCode ForwardTransportFast::execute() { for (int i=0; i<(int)fPidVector.size(); i++) { // add vertices for G4 tracking (status code = 1) - HepMC::GenVertex* gVertex = new HepMC::GenVertex (fPosVector.at(i)); - HepMC::GenParticle* gParticle = new HepMC::GenParticle(fMomVector.at(i), fPidVector.at(i), 1); + HepMC::GenVertexPtr gVertex = HepMC::newGenVertexPtr (fPosVector.at(i)); + HepMC::GenParticlePtr gParticle = HepMC::newGenParticlePtr(fMomVector.at(i), fPidVector.at(i), 1); gVertex->add_particle_out(gParticle);