Skip to content
Snippets Groups Projects
Commit eeda116a authored by Eric Torrence's avatar Eric Torrence
Browse files

Merge remote-tracking branch 'upstream/master'

parents 8e04003e a580aa08
No related branches found
No related tags found
No related merge requests found
Showing
with 633 additions and 20 deletions
......@@ -25,7 +25,7 @@ build_image:
- mkdir build
- cd build
- set +e && source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh; set -e
- set +e && asetup --input=../../calypso/asetup.faser Athena,22.0.40; set -e
- set +e && asetup --input=../../calypso/asetup.faser Athena,22.0.49; set -e
- cmake ../../calypso
- make -j 3
artifacts:
......@@ -41,7 +41,7 @@ test_unittest:
- yum -y install man
- cd build
- set +e && source ${ATLAS_LOCAL_ROOT_BASE}/user/atlasLocalSetup.sh; set -e
- set +e && asetup --input=../../calypso/asetup.faser Athena,22.0.40; set -e
- set +e && asetup --input=../../calypso/asetup.faser Athena,22.0.49; set -e
- set +e && source `find . -name 'setup.sh'`; set -e
- ctest -j3
dependencies:
......
......@@ -18,12 +18,16 @@ CaloHitCollection_PERS* CaloHitCollectionCnv::createPersistent(CaloHitCollection
CaloHitCollection* CaloHitCollectionCnv::createTransient() {
MsgStream mlog(msgSvc(), "CaloHitCollectionConverter" );
CaloHitCollectionCnv_p1 converter_p1;
CaloHitCollectionCnv_p1a converter_p1a;
static const pool::Guid p1_guid("134E8045-AB99-43EF-9AD1-324C48830B64");
// static const pool::Guid p1_guid("B2573A16-4B46-4E1E-98E3-F93421680779");
static const pool::Guid p1a_guid("02C108EE-0AD9-4444-83BD-D5273FCDEF6F");
CaloHitCollection *trans_cont(0);
if( this->compareClassGuid(p1_guid)) {
if( this->compareClassGuid(p1a_guid)) {
std::unique_ptr< CaloHitCollection_p1a > col_vect( this->poolReadObject< CaloHitCollection_p1a >() );
trans_cont = converter_p1a.createTransient( col_vect.get(), mlog );
} else if( this->compareClassGuid(p1_guid)) {
std::unique_ptr< CaloHitCollection_p1 > col_vect( this->poolReadObject< CaloHitCollection_p1 >() );
trans_cont = converter_p1.createTransient( col_vect.get(), mlog );
} else {
......
......@@ -8,12 +8,14 @@
#include "FaserCaloSimEvent/CaloHitCollection.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollection_p1.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollectionCnv_p1.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollection_p1a.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollectionCnv_p1a.h"
#include "AthenaPoolCnvSvc/T_AthenaPoolCustomCnv.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
// typedef to the latest persistent version
typedef CaloHitCollection_p1 CaloHitCollection_PERS;
typedef CaloHitCollectionCnv_p1 CaloHitCollectionCnv_PERS;
typedef CaloHitCollection_p1a CaloHitCollection_PERS;
typedef CaloHitCollectionCnv_p1a CaloHitCollectionCnv_PERS;
class CaloHitCollectionCnv : public T_AthenaPoolCustomCnv<CaloHitCollection, CaloHitCollection_PERS > {
friend class CnvFactory<CaloHitCollectionCnv>;
......
......@@ -22,5 +22,5 @@ atlas_add_dictionary( FaserCaloSimEventTPCnvDict
FaserCaloSimEventTPCnv/CaloSimEventTPCnvDict.h
FaserCaloSimEventTPCnv/selection.xml
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthenaPoolCnvSvcLib GaudiKernel GeneratorObjectsTPCnv FaserCaloSimEvent TestTools StoreGateLib SGtests Identifier FaserCaloSimEventTPCnv )
LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthenaPoolCnvSvcLib GaudiKernel GeneratorObjectsTPCnv FaserCaloSimEvent TestTools StoreGateLib SGtests Identifier FaserCaloSimEventTPCnv AthenaKernel)
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef CALOHITCOLLECTIONCNV_P1A_H
#define CALOHITCOLLECTIONCNV_P1A_H
// CaloHitCollectionCnv_p1, T/P separation of Calo Hits
// author D.Costanzo <davide.costanzo@cern.ch>
// author O.Arnaez <olivier.arnaez@cern.ch>
#include "AthenaPoolCnvSvc/T_AthenaPoolTPConverter.h"
#include "FaserCaloSimEvent/CaloHitCollection.h"
#include "CaloHitCollection_p1a.h"
class CaloHitCollectionCnv_p1a : public T_AthenaPoolTPCnvBase<CaloHitCollection, CaloHitCollection_p1a>
{
public:
CaloHitCollectionCnv_p1a() {};
virtual CaloHitCollection* createTransient(const CaloHitCollection_p1a* persObj, MsgStream &log);
virtual void persToTrans(const CaloHitCollection_p1a* persCont,
CaloHitCollection* transCont,
MsgStream &log) ;
virtual void transToPers(const CaloHitCollection* transCont,
CaloHitCollection_p1a* persCont,
MsgStream &log) ;
private:
static const double m_persEneUnit;
static const double m_persLenUnit;
static const double m_persAngUnit;
static const double m_2bHalfMaximum;
static const int m_2bMaximum;
};
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef CALOHITCOLLECTION_P1A_H
#define CALOHITCOLLECTION_P1A_H
/*
Authors: Davide Costanzo Rob Duxfield
Modified for FASER by D. Casper - change internal unsigned short counters to unsigned long
p2, p3, etc reserved for possible future ATLAS changes
*/
#include <vector>
#include <string>
class CaloHitCollection_p1a
{
public:
/// Default constructor
CaloHitCollection_p1a ();
//private:
std::vector<float> m_hit1_meanTime; // 1 element per string
std::vector<float> m_hit1_x0; //
std::vector<float> m_hit1_y0; //
std::vector<float> m_hit1_z0; //
std::vector<float> m_hit1_theta; //
std::vector<float> m_hit1_phi; //
std::vector<unsigned long> m_nHits; //
std::vector<unsigned short> m_hitEne_2b; // 1 element per hit
std::vector<unsigned short> m_hitLength_2b; //
std::vector<unsigned short> m_dTheta; // 1 element per hit except for first hit in string
std::vector<unsigned short> m_dPhi; //
std::vector<float> m_hitEne_4b; // 1 element per hit with m_hitEne_2b[i] == 2**16
std::vector<float> m_hitLength_4b; // 1 element per hit with m_hitLength_2b[i] == 2**16
std::vector<unsigned long> m_barcode;
std::vector<unsigned long> m_mcEvtIndex;
std::vector<char> m_evtColl;
std::vector<unsigned long> m_nBC;
std::vector<unsigned long> m_id;
std::vector<unsigned long> m_nId;
};
// inlines
inline
CaloHitCollection_p1a::CaloHitCollection_p1a () {}
#endif
......@@ -15,6 +15,8 @@
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCnv_p1.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollectionCnv_p1.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollection_p1.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollectionCnv_p1a.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollection_p1a.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHit_p1.h"
#endif // CALOEVENTTPCNV_CALOSIMEVENTTPCNVDICT_H
......@@ -4,4 +4,6 @@
<class name="CaloHit_p1" />
<class name="std::vector<CaloHit_p1>" />
<class name="CaloHitCollection_p1" id="134E8045-AB99-43EF-9AD1-324C48830B64" />
<class name="CaloHitCollection_p1a" id="02C108EE-0AD9-4444-83BD-D5273FCDEF6F" />
</lcgdict>
......@@ -14,7 +14,10 @@
#include "CLHEP/Geometry/Point3D.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/ThreadLocalContext.h"
// Athena
#include "AthenaKernel/ExtendedEventContext.h"
#include "StoreGate/StoreGateSvc.h"
// * * * stolen from eflowRec * * * //
......@@ -60,7 +63,7 @@ const double CaloHitCollectionCnv_p1::m_2bHalfMaximum = pow(2.0, 15.0);
const int CaloHitCollectionCnv_p1::m_2bMaximum = (unsigned short)(-1);
void CaloHitCollectionCnv_p1::transToPers(const CaloHitCollection* transCont, CaloHitCollection_p1* persCont, MsgStream &/*log*/)
void CaloHitCollectionCnv_p1::transToPers(const CaloHitCollection* transCont, CaloHitCollection_p1* persCont, MsgStream &log)
{
// Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and
// persistifies the end point of each hit plus the start point of the first hit in each string.
......@@ -77,6 +80,8 @@ void CaloHitCollectionCnv_p1::transToPers(const CaloHitCollection* transCont, Ca
static const double dRcut = 1.0e-7;
static const double dTcut = 1.0;
const EventContext& ctx = Gaudi::Hive::currentContext();
const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy();
const HepMcParticleLink * lastLink=nullptr;
int lastId = -1;
double stringFirstTheta = 0.0;
......@@ -98,11 +103,23 @@ void CaloHitCollectionCnv_p1::transToPers(const CaloHitCollection* transCont, Ca
if ( !lastLink || (siHit->particleLink() != *lastLink) ) {
// store barcode once for set of consecutive hits with same barcode
// store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode
lastLink = &(siHit->particleLink());
persCont->m_barcode.push_back(lastLink->barcode());
persCont->m_mcEvtIndex.push_back(lastLink->eventIndex());
log << MSG::ALWAYS << "TP: Event collection from link: " << lastLink->getEventCollection() << " ; char value: " << lastLink->getEventCollectionAsChar() << " barcode: " << lastLink->barcode() << endmsg;
unsigned short index{0};
const HepMcParticleLink::index_type position =
HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(),
lastLink->getEventCollection(),
proxy).at(0);
if (position!=0) {
index = lastLink->eventIndex();
if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) {
log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg;
}
}
persCont->m_mcEvtIndex.push_back(index);
persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar());
if (idx > 0) {
......@@ -238,18 +255,29 @@ void CaloHitCollectionCnv_p1::transToPers(const CaloHitCollection* transCont, Ca
persCont->m_nBC.push_back(idx - endBC);
persCont->m_nId.push_back(idx - endId);
persCont->m_nHits.push_back(idx - endHit);
log << MSG::ALWAYS << "nBC: " << persCont->m_nBC << endmsg;
log << MSG::ALWAYS << "nId: " << persCont->m_nId << endmsg;
log << MSG::ALWAYS << "nHits:" << persCont->m_nHits << endmsg;
}
CaloHitCollection* CaloHitCollectionCnv_p1::createTransient(const CaloHitCollection_p1* persObj, MsgStream &log) {
for (size_t i = 0; i < persObj->m_evtColl.size(); i++)
{
if (persObj->m_evtColl[i] != 'a')
log << MSG::ALWAYS << "Corrupted in createTransient: " << persObj->m_evtColl[i] << endmsg;
}
std::unique_ptr<CaloHitCollection> trans(std::make_unique<CaloHitCollection>("DefaultCollectionName",persObj->m_nHits.size()));
persToTrans(persObj, trans.get(), log);
return(trans.release());
}
void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont, CaloHitCollection* transCont, MsgStream &/*log*/)
void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont, CaloHitCollection* transCont, MsgStream &log)
{
const EventContext& ctx = Gaudi::Hive::currentContext();
unsigned int hitCount = 0;
unsigned int angleCount = 0;
unsigned int idxBC = 0;
......@@ -260,12 +288,15 @@ void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont,
unsigned int endBC = 0;
unsigned int endId = 0;
log << MSG::ALWAYS << "nHits.size(): " << persCont->m_nHits.size() << " nBC.size(): " << persCont->m_nBC.size() << endmsg;
log << MSG::ALWAYS << " nBC: " << persCont->m_nBC << endmsg;
for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) {
if (persCont->m_nHits[i]) {
const unsigned int start = endHit;
endHit += persCont->m_nHits[i];
log << MSG::ALWAYS << "i: " << i << " persCont->m_nHits[i]: " << persCont->m_nHits[i] << " start: " << start << " endHit: " << endHit << endmsg;
const double t0 = persCont->m_hit1_meanTime[i];
const double theta0 = persCont->m_hit1_theta[i];
......@@ -273,7 +304,6 @@ void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont,
HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]);
for (unsigned int j = start; j < endHit; j++) {
if (j >= endBC + persCont->m_nBC[idxBC])
endBC += persCont->m_nBC[idxBC++];
......@@ -299,7 +329,13 @@ void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont,
HepGeom::Point3D<double> endThis( endLast + r );
HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), HepMcParticleLink::IS_INDEX );
HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX;
if (persCont->m_mcEvtIndex[idxBC] == 0) {
flag = HepMcParticleLink::IS_POSITION;
}
log << MSG::ALWAYS << "PT: Event collection from persCont: " << persCont->m_evtColl[idxBC] << " ; idxBC: " << idxBC << " evtColl.size(): " << persCont->m_evtColl.size() << " evtIndex.size(): " << persCont->m_mcEvtIndex.size() << endmsg;
HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx );
transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]);
endLast = endThis;
......@@ -309,4 +345,5 @@ void CaloHitCollectionCnv_p1::persToTrans(const CaloHitCollection_p1* persCont,
}
}
}
log << MSG::ALWAYS << "hitCount: " << hitCount << endmsg;
}
/*
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
#include "FaserCaloSimEvent/CaloHit.h"
#include "FaserCaloSimEvent/CaloHitCollection.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollection_p1a.h"
#include "FaserCaloSimEventTPCnv/CaloHits/CaloHitCollectionCnv_p1a.h"
#include "GeneratorObjects/HepMcParticleLink.h"
#include <cmath>
//CLHEP
#include "CLHEP/Geometry/Point3D.h"
// Gaudi
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/ThreadLocalContext.h"
// Athena
#include "AthenaKernel/ExtendedEventContext.h"
#include "StoreGate/StoreGateSvc.h"
// * * * stolen from eflowRec * * * //
inline double phicorr(double a)
{
if (a <= -M_PI)
{
return a+(2*M_PI*floor(-(a-M_PI)/(2*M_PI)));
}
else if (a > M_PI)
{
return a-(2*M_PI*floor((a+M_PI)/(2*M_PI)));
}
else
{
return a;
}
}
// * * * stolen from eflowRec * * * //
inline double cycle(double a, double b)
{
double del = b-a;
if (del > M_PI)
{
return a+2.0*M_PI;
}
else if (del < -M_PI)
{
return a-2.0*M_PI;
}
else
{
return a;
}
}
const double CaloHitCollectionCnv_p1a::m_persEneUnit = 1.0e-5;
const double CaloHitCollectionCnv_p1a::m_persLenUnit = 1.0e-5;
const double CaloHitCollectionCnv_p1a::m_persAngUnit = 1.0e-5;
const double CaloHitCollectionCnv_p1a::m_2bHalfMaximum = pow(2.0, 15.0);
const int CaloHitCollectionCnv_p1a::m_2bMaximum = (unsigned short)(-1);
void CaloHitCollectionCnv_p1a::transToPers(const CaloHitCollection* transCont, CaloHitCollection_p1a* persCont, MsgStream &log)
{
// Finds hits belonging to a "string" (in which the end point of one hit is the same as the start point of the next) and
// persistifies the end point of each hit plus the start point of the first hit in each string.
//
// Further compression is achieved by optimising the storage of the position vectors:- start (x,y,z) and (theta,phi) of
// first hit are stored as floats, (delta_theta,delta_phi) relative to the fisrst hit are stored as 2 byte numbers and
// used to specify the hit direction. All hit lengths are stored as 2 byte numbers.
//
// Additional savings are achieved by storing the energy loss for each hit as a 2 byte number and only storing the mean
// time of the first hit per string.
//
// See http://indico.cern.ch/getFile.py/access?contribId=11&resId=2&materialId=slides&confId=30893 for more info.
static const double dRcut = 1.0e-7;
static const double dTcut = 1.0;
const EventContext& ctx = Gaudi::Hive::currentContext();
const IProxyDict* proxy = Atlas::getExtendedEventContext(ctx).proxy();
const HepMcParticleLink * lastLink=nullptr;
int lastId = -1;
double stringFirstTheta = 0.0;
double stringFirstPhi = 0.0;
double lastT = 0.0;
double persSumE = 0.0;
double transSumE = 0.0;
unsigned int idx = 0;
unsigned int endBC = 0;
unsigned int endId = 0;
unsigned int endHit = 0;
HepGeom::Point3D<double> lastTransEnd(0.0, 0.0, 0.0);
HepGeom::Point3D<double> lastPersEnd(0.0, 0.0, 0.0);
for (CaloHitCollection::const_iterator it = transCont->begin(); it != transCont->end(); ++it) {
CaloHitCollection::const_iterator siHit = it;
if ( !lastLink || (siHit->particleLink() != *lastLink) ) {
// store barcode, eventIndex and McEventCollection once for set of consecutive hits with same barcode
lastLink = &(siHit->particleLink());
persCont->m_barcode.push_back(lastLink->barcode());
unsigned short index{0};
const HepMcParticleLink::index_type position =
HepMcParticleLink::getEventPositionInCollection(lastLink->eventIndex(),
lastLink->getEventCollection(),
proxy).at(0);
if (position!=0) {
index = lastLink->eventIndex();
if(lastLink->eventIndex()!=static_cast<HepMcParticleLink::index_type>(index)) {
log << MSG::WARNING << "Attempting to persistify an eventIndex larger than max unsigned short!" << endmsg;
}
}
persCont->m_mcEvtIndex.push_back(index);
persCont->m_evtColl.push_back(lastLink->getEventCollectionAsChar());
if (idx > 0) {
persCont->m_nBC.push_back(idx - endBC);
endBC = idx;
}
}
if ( (int)siHit->identify() != lastId ) {
// store id once for set of consecutive hits with same barcode
lastId = siHit->identify();
persCont->m_id.push_back(lastId);
if (idx > 0) {
persCont->m_nId.push_back(idx - endId);
endId = idx;
}
}
HepGeom::Point3D<double> st = siHit->localStartPosition();
HepGeom::Point3D<double> en = siHit->localEndPosition();
const double dx = st.x() - lastTransEnd.x();
const double dy = st.y() - lastTransEnd.y();
const double dz = st.z() - lastTransEnd.z();
const double t = siHit->meanTime();
const double dRLast = sqrt(dx * dx + dy * dy + dz * dz); // dR between end of previous hit and start of current one
const double dTLast = fabs(t - lastT);
CLHEP::Hep3Vector direction(0.0, 0.0, 0.0);
double theta = 0.0;
double phi = 0.0;
bool startNewString = false;
if (dRLast < dRcut && dTLast < dTcut) {
// hit is part of existing string
direction = CLHEP::Hep3Vector( en.x() - lastPersEnd.x(), en.y() - lastPersEnd.y(), en.z() - lastPersEnd.z() );
theta = direction.theta();
phi = phicorr( direction.phi() );
const int dTheta_2b = (int)( (theta - stringFirstTheta) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
const int dPhi_2b = (int)( (cycle(phi, stringFirstPhi) - stringFirstPhi) / m_persAngUnit + m_2bHalfMaximum + 0.5 );
if ( dTheta_2b < m_2bMaximum && dTheta_2b >= 0 && dPhi_2b < m_2bMaximum && dPhi_2b >= 0) {
persCont->m_dTheta.push_back(dTheta_2b);
persCont->m_dPhi.push_back(dPhi_2b);
theta = stringFirstTheta + ( (double)dTheta_2b - m_2bHalfMaximum ) * m_persAngUnit;
phi = stringFirstPhi + ( (double)dPhi_2b - m_2bHalfMaximum ) * m_persAngUnit;
phi = phicorr(phi);
}
else {
startNewString = true;
}
}
if (startNewString || dRLast >= dRcut || dTLast >= dTcut) {
// begin new hit string
direction = CLHEP::Hep3Vector( en.x() - st.x(), en.y() - st.y(), en.z() - st.z() );
theta = direction.theta();
phi = phicorr( direction.phi() );
persCont->m_hit1_meanTime.push_back(t);
persCont->m_hit1_x0.push_back(st.x());
persCont->m_hit1_y0.push_back(st.y());
persCont->m_hit1_z0.push_back(st.z());
persCont->m_hit1_theta.push_back(theta);
persCont->m_hit1_phi.push_back(phi);
lastPersEnd = HepGeom::Point3D<double>(st.x(), st.y(), st.z());
stringFirstTheta = theta;
stringFirstPhi = phi;
if (idx > 0) {
persCont->m_nHits.push_back(idx - endHit);
endHit = idx;
}
}
lastTransEnd = HepGeom::Point3D<double>(en.x(), en.y(), en.z());
transSumE += siHit->energyLoss();
const int eneLoss_2b = (int)((transSumE - persSumE) / m_persEneUnit + 0.5); // calculated to allow recovery sum over
// whole hit string to chosen precision
const int hitLength_2b = (int)(direction.mag() / m_persLenUnit + 0.5); // calculated to give the correct position to
// the chosen precision, NOT the length of the
// hit (small difference in practice).
double eneLoss = 0.0;
if (eneLoss_2b >= m_2bMaximum) {
eneLoss = siHit->energyLoss();
persCont->m_hitEne_2b.push_back(m_2bMaximum);
persCont->m_hitEne_4b.push_back(eneLoss);
}
else {
eneLoss = eneLoss_2b * m_persEneUnit;
persCont->m_hitEne_2b.push_back(eneLoss_2b);
}
double length = 0.0;
if (hitLength_2b >= m_2bMaximum) {
length = direction.mag();
persCont->m_hitLength_2b.push_back(m_2bMaximum);
persCont->m_hitLength_4b.push_back(direction.mag());
}
else {
length = hitLength_2b * m_persLenUnit;
persCont->m_hitLength_2b.push_back(hitLength_2b);
}
CLHEP::Hep3Vector persDir(length, 0.0, 0.0);
persDir.setTheta(theta);
persDir.setPhi(phi);
lastPersEnd = (CLHEP::Hep3Vector)lastPersEnd + persDir;
persSumE += eneLoss;
lastT = t;
++idx;
}
persCont->m_nBC.push_back(idx - endBC);
persCont->m_nId.push_back(idx - endId);
persCont->m_nHits.push_back(idx - endHit);
}
CaloHitCollection* CaloHitCollectionCnv_p1a::createTransient(const CaloHitCollection_p1a* persObj, MsgStream &log) {
for (size_t i = 0; i < persObj->m_evtColl.size(); i++)
{
if (persObj->m_evtColl[i] != 'a')
log << MSG::ALWAYS << "Corrupted in createTransient: " << persObj->m_evtColl[i] << endmsg;
}
std::unique_ptr<CaloHitCollection> trans(std::make_unique<CaloHitCollection>("DefaultCollectionName",persObj->m_nHits.size()));
persToTrans(persObj, trans.get(), log);
return(trans.release());
}
void CaloHitCollectionCnv_p1a::persToTrans(const CaloHitCollection_p1a* persCont, CaloHitCollection* transCont, MsgStream &/*log*/)
{
const EventContext& ctx = Gaudi::Hive::currentContext();
unsigned int hitCount = 0;
unsigned int angleCount = 0;
unsigned int idxBC = 0;
unsigned int idxId = 0;
unsigned int idxEne4b = 0;
unsigned int idxLen4b = 0;
unsigned int endHit = 0;
unsigned int endBC = 0;
unsigned int endId = 0;
for (unsigned int i = 0; i < persCont->m_nHits.size(); i++) {
if (persCont->m_nHits[i]) {
const unsigned int start = endHit;
endHit += persCont->m_nHits[i];
const double t0 = persCont->m_hit1_meanTime[i];
const double theta0 = persCont->m_hit1_theta[i];
const double phi0 = persCont->m_hit1_phi[i];
HepGeom::Point3D<double> endLast(persCont->m_hit1_x0[i], persCont->m_hit1_y0[i], persCont->m_hit1_z0[i]);
for (unsigned int j = start; j < endHit; j++) {
if (j >= endBC + persCont->m_nBC[idxBC])
endBC += persCont->m_nBC[idxBC++];
if (j >= endId + persCont->m_nId[idxId])
endId += persCont->m_nId[idxId++];
const double eneLoss_2b = persCont->m_hitEne_2b[hitCount];
const double hitLength_2b = persCont->m_hitLength_2b[hitCount];
const double eneLoss = (eneLoss_2b < m_2bMaximum) ? eneLoss_2b * m_persEneUnit : persCont->m_hitEne_4b[idxEne4b++];
const double length = (hitLength_2b < m_2bMaximum) ? hitLength_2b * m_persLenUnit : persCont->m_hitLength_4b[idxLen4b++];
const double dTheta = (j > start) ? ((double)persCont->m_dTheta[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
const double dPhi = (j > start) ? ((double)persCont->m_dPhi[angleCount] - m_2bHalfMaximum) * m_persAngUnit : 0.0;
const double meanTime = t0;
const double theta = theta0 + dTheta;
const double phi = phicorr(phi0 + dPhi);
CLHEP::Hep3Vector r(length, 0.0, 0.0);
r.setTheta(theta);
r.setPhi(phi);
HepGeom::Point3D<double> endThis( endLast + r );
HepMcParticleLink::PositionFlag flag = HepMcParticleLink::IS_INDEX;
if (persCont->m_mcEvtIndex[idxBC] == 0) {
flag = HepMcParticleLink::IS_POSITION;
}
HepMcParticleLink partLink( persCont->m_barcode[idxBC], persCont->m_mcEvtIndex[idxBC], HepMcParticleLink::ExtendedBarCode::eventCollectionFromChar(persCont->m_evtColl[idxBC]), flag, ctx );
transCont->Emplace( endLast, endThis, eneLoss, meanTime, partLink, persCont->m_id[idxId]);
endLast = endThis;
++hitCount;
if (j > start) ++angleCount;
}
}
}
}
......@@ -17,5 +17,7 @@ def EcalSensorSDCfg(ConfigFlags, name="EcalSensorSD", **kwargs):
kwargs.setdefault("LogicalVolumeNames", ["Ecal::_dd_Geometry_DownstreamRegion_Ecal_Modules_OutCell"])
kwargs.setdefault("OutputCollectionNames", [bare_collection_name])
# result.merge(acc)
return result, EcalSensorSDTool(name, **kwargs)
result = ComponentAccumulator()
result.setPrivateTools(CompFactory.EcalSensorSDTool(name, **kwargs))
return result
......@@ -89,7 +89,7 @@ asetup
source build/x8*/setup.sh
#
# Do this by hand
# asetup --input="$release_directory/calypso/asetup.faser" Athena,22.0.40
# asetup --input="$release_directory/calypso/asetup.faser" Athena,22.0.49
# source "$release_directory/build/x8*/setup.sh"
#
#
......
......@@ -21,5 +21,9 @@
<service name="sqlite_file:///cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/sqlite200/ALLP200.db" accessMode="read" />
</logicalservice>
<logicalservice name="COOLOFL_TDAQ">
<service name="sqlite_file:data/sqlite200/CABP200.db" accessMode="read" />
<service name="sqlite_file:///cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/sqlite200/CABP200.db" accessMode="read" />
</logicalservice>
</servicelist>
......@@ -50,6 +50,10 @@ def FaserByteStreamCnvSvcCfg(configFlags, **kwargs):
from WaveByteStream.WaveByteStreamConfig import WaveByteStreamCfg
result.merge(WaveByteStreamCfg(configFlags))
# Configure TrackerByteStream converter
from TrackerByteStream.TrackerByteStreamConfig import TrackerByteStreamCfg
result.merge(TrackerByteStreamCfg(configFlags))
# Load ByteStreamEventStorageInputSvc
bsInputSvc = CompFactory.FaserByteStreamInputSvc
result.addService(bsInputSvc(name = "FaserByteStreamInputSvc"))
......
......@@ -2,7 +2,7 @@
## Setup
In Athena 22.0.40, the scipy module is missing from the LCG environment for some reason.
In Athena 22.0.49, the scipy module is missing from the LCG environment for some reason.
To use the generator, the following command is required after the usual steps to setup the release:
......
################################################################################
# Package: ReadGenie
################################################################################
# Declare the package name:
atlas_subdir( GenieReader )
# Install files from the package:
atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
atlas_install_joboptions( share/*.py )
\ No newline at end of file
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
from AthenaCommon.AppMgr import ServiceMgr as svcMgr
from GeneratorModules.EvgenAlg import EvgenAlg
from AthenaPython.PyAthena import StatusCode
from AthenaCommon.SystemOfUnits import GeV, m
import ROOT
__author__ = "Dave Caser <dcasper@uci.edu>"
class GenieReader(EvgenAlg):
def __init__(self, name="GenieReader", MCEventKey="BeamTruthEvent"):
super(GenieReader,self).__init__(name=name)
self.McEventKey = MCEventKey
return
def fillEvent(self, evt):
try:
from AthenaPython.PyAthena import HepMC3 as HepMC
except ImportError:
from AthenaPython.PyAthena import HepMC as HepMC
evt.weights().push_back(1.0)
pos = HepMC.FourVector(self.evtStore["vx"]*m, self.evtStore["vy"]*m, self.evtStore["vz"]*m, 0)
gv = HepMC.GenVertex(pos)
ROOT.SetOwnership(gv, False)
evt.add_vertex(gv)
nParticles = self.evtStore["n"]
pdgc = list(self.evtStore["pdgc"])
status = list(self.evtStore["status"])
px = list(self.evtStore["px"])
py = list(self.evtStore["py"])
pz = list(self.evtStore["pz"])
E = list(self.evtStore["E"])
M = list(self.evtStore["M"])
for i in range(nParticles):
gp = HepMC.GenParticle()
mom = HepMC.FourVector(px[i]*GeV, py[i]*GeV, pz[i]*GeV, E[i]*GeV)
gp.set_momentum(mom)
gp.set_generated_mass(M[i]*GeV)
gp.set_pdg_id(pdgc[i])
genie_status = status[i]
if (genie_status == 0): # initial particle
hepmc_status = 4
elif (genie_status == 1): # stable final particle
hepmc_status = 1
elif (genie_status == 3): # decayed particle
hepmc_status = 2
else: # catch-all informational particle
hepmc_status = 3
gp.set_status(hepmc_status)
ROOT.SetOwnership(gp, False)
if (hepmc_status == 4):
gv.add_particle_in(gp)
else:
gv.add_particle_out(gp)
return StatusCode.Success
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
#
# Usage: athena.py -c'INPUT=["faser.1.gfaser.root","faser.2.gfaser.root"]; OUTPUT="gFaser.EVNT.pool.root"; EVTMAX=100' GenieReader/GenieReader_jobOptions.py >& GenieReader.log
#
# INPUT and OUTPUT are mandatory (INPUT can be a list, as shown above)
# EVTMAX can be omitted; then all input files are read to the end
#
if not 'INPUT' in dir():
print("Missing INPUT parameter")
exit()
if not isinstance(INPUT, (list,tuple)):
INPUT = [INPUT]
pass
if not 'OUTPUT' in dir():
print("Missing OUTPUT parameter")
exit()
import AthenaRootComps.ReadAthenaRoot
svcMgr.EventSelector.InputCollections = INPUT
svcMgr.EventSelector.TupleName = "gFaser"
from GenieReader.GenieReaderAlg import GenieReader
from AthenaCommon.AlgSequence import AlgSequence
job = AlgSequence()
job += GenieReader()
from AthenaPoolCnvSvc.WriteAthenaPool import AthenaPoolOutputStream
ostream = AthenaPoolOutputStream( "StreamEVGEN" , OUTPUT, noTag=True )
ostream.ItemList += [ "xAOD::EventInfo#*", "xAOD::EventAuxInfo#*",
"EventInfo#*",
"McEventCollection#*" ]
if not 'EVTMAX' in dir():
EVTMAX = -1
theApp.EvtMax = EVTMAX
......@@ -4,6 +4,7 @@
#include "NeutrinoSimEventTPCnv/NeutrinoHits/NeutrinoHitCollectionCnv_p1.h"
#include "NeutrinoSimEventTPCnv/NeutrinoHits/NeutrinoHit_p1.h"
#include "NeutrinoSimEventTPCnv/NeutrinoHits/NeutrinoHitCollectionCnv_p1a.h"
#include "NeutrinoHitCollectionCnv.h"
......@@ -18,12 +19,16 @@ NeutrinoHitCollection_PERS* NeutrinoHitCollectionCnv::createPersistent(NeutrinoH
NeutrinoHitCollection* NeutrinoHitCollectionCnv::createTransient() {
MsgStream mlog(msgSvc(), "NeutrinoHitCollectionConverter" );
NeutrinoHitCollectionCnv_p1 converter_p1;
NeutrinoHitCollectionCnv_p1a converter_p1a;
static const pool::Guid p1_guid("2CA7AF71-1494-4378-BED4-AFB5C54398AA");
// static const pool::Guid p1_guid("B2573A16-4B46-4E1E-98E3-F93421680779");
static const pool::Guid p1a_guid("0A3CD37D-64CE-405D-98CF-595D678B14B7");
NeutrinoHitCollection *trans_cont(0);
if( this->compareClassGuid(p1_guid)) {
if( this->compareClassGuid(p1a_guid)) {
std::unique_ptr< NeutrinoHitCollection_p1a > col_vect( this->poolReadObject< NeutrinoHitCollection_p1a >() );
trans_cont = converter_p1a.createTransient( col_vect.get(), mlog );
} else if( this->compareClassGuid(p1_guid)) {
std::unique_ptr< NeutrinoHitCollection_p1 > col_vect( this->poolReadObject< NeutrinoHitCollection_p1 >() );
trans_cont = converter_p1.createTransient( col_vect.get(), mlog );
} else {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment