CMCServer.C
//-----------------------------------------------------------------------------
// $Header:
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include "CMCServer.h"
#include "C4Momentum.h"
#include "CSidc1.h"
#include "CSidc2.h"
#include "CPadChamber.h"
#include "CRandom.h"
CMCServer::CMCServer(const char* setupFile)
{
//
// Create and read setup
//
setup = new CMCServerSetup;
if (setupFile == 0) {
CString setupName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_SETUPFILE_MCSERVER);
setup->read(setupName.data());
}
else
setup->read(setupFile);
mcvertex = 0;
mcADDevertex = 0;
indexOfLastTrack = 0;
printon = 0;
}
CMCServer::~CMCServer()
{
geantTracks.clearAndDestroy();
for (int i=0; i<MaxDetectors; i++) {
digiHits[i].clearAndDestroy();
}
delete mcvertex;
delete mcADDevertex;
delete setup;
}
CMCServer::CMCServer(const CMCServer &)
{
}
CMCServer & CMCServer::operator = (const CMCServer &)
{
return *this;
}
void CMCServer::listSetup(ostream& ost) const
{
ost << "\nSetup listing of MonteCarlo server" << endl;
CString line('=', 46);
ost << line << endl;
setup->list(ost);
}
void CMCServer::resetTracksAndHits()
{
geantTracks.clearAndDestroy(); // reset 'particle' list
for (int id=0; id<MaxDetectors; id++) {
digiHits[id].clearAndDestroy(); // ... and lists of digiHits
}
indexOfLastTrack = 0;
}
CBoolean CMCServer::update(CEventServer& es)
{
return unpack(es);
}
CBoolean CMCServer::loadEventWithCloseVertex(CEventServer& es, CVertex *vtx)
{
if (printon > 0) {
cout << "searching for vertex with z=" << vtx->getZ() << endl;
}
int tries = 0;
CBoolean found = False;
while ( !found ) {
tries++;
if ( !es.loadNextEvent() ) return False;
resetTracksAndHits();
if ( !unpack(es) ) return False;
if (printon > 2) {
cout << "\tnext event has vertex at z=" << getMCVertex()->getZ() << endl;
}
found = abs(getMCVertex()->getZ() - vtx->getZ() + setup->getCloseVertexOffset() ) < setup->getCloseVertexDz();
}
if (printon > 0) {
cout << "found close vertex with dz="
<< getMCVertex()->getZ() - vtx->getZ() + setup->getCloseVertexOffset()
<< " after " << tries << " tries." << endl;
}
return True;
}
void CMCServer::listGeantTracks(ostream& os)
{
os << "========================================== start list of GEANT Tracks" << endl;
for (int i=0; i<geantTracks.length(); i++) {
geantTracks(i)->print(os);
}
os << "========================================== end list of GEANT Tracks" << endl;
}
CBoolean CMCServer::unpack(CEventServer& es)
{
const int geantTrackLabel = setup->getGeantTrackLabel();
const int rich1HitLabel = setup->getRich1HitLabel();
const int rich2HitLabel = setup->getRich2HitLabel();
const int sidc1HitLabel = setup->getSidc1HitLabel();
const int sidc2HitLabel = setup->getSidc2HitLabel();
const int padChamberHitLabel = setup->getPadChamberHitLabel();
const CLabel * mcpart = es.getLabel(geantTrackLabel);
if (mcpart == 0) {
cerr << "CMCServer::unpack():\n";
cerr << " ERROR \n";
cerr << " label " << geantTrackLabel << " not found" << endl;
return False;
}
int *data = (int *) mcpart->getData();
if (data == 0) {
cerr << "CMCServer::unpack():\n"
<< "\tERROR\n"
<< "\tno data found for label: " << geantTrackLabel << endl;
return False;
}
// structure of this label:
// first data word is number of entries
// then follow 14 words per pparticle:
// int trackNr for primary particles: GEANT track number
// for secondaries: GEANT track number of
// parent of particle, if from stack:
// ITRAK+10000*ISTAK
// int parent for primary particle: HIJET parent track #
// special: 1: beam delta rays in target
// 2: additional MC electron
// for secondaries: GEANT id of parent
// CString volume[4] primaries: ' '
// secondaries: GEANT volume name at creation
// CString process[4] primaries: 'prim', except for
// 'pDAL' pi-DALITZ electron
// 'eDAL' eta-DALITZ
// 'eepr' ee-pair from resonances
// 'ADDe' additional MC electron
// secondaries: GEANT creation process name
// e.g. 'PAIR','HADR',...
// int idGeant GEANT particle id
// int nHitsUV1 number of hits in UV1 from this particle
// int nHitsUV2 dito UV2
// float vect[7] VECT(1:7) at creation time: X,Y,Z,Px/P,Py/P,Pz/P,P
CMCGeantTrack *track;
C3Momentum orig, p;
float p0;
int i, trackIndex, ptype;
int index = 0;
data++; // skip number of entries in list
for (i=0; i<mcpart->getSize()-1; i+=14) {
track = new CMCGeantTrack();
track->setIndex ( indexOfLastTrack + index++ ); // index for this track (for operator==)
track->setTrackNr ( *data++ );
track->setParent ( *data++ );
track->setVolume ( CString( (char *) data++, 4) );
track->setProcess ( CString( (char *) data++, 4) );
ptype = *data++;
track->getParticle().setGeantId( (CMCParticle::CMCGeantId) ptype );
track->setNHitsUV1 ( *data++ );
track->setNHitsUV2 ( *data++ );
orig.setX(*((float*)data++));
orig.setY(*((float*)data++));
orig.setZ(*((float*)data++));
p.setX(*((float*)data++));
p.setY(*((float*)data++));
p.setZ(*((float*)data++));
p0 = *((float*)data++);
//
// update the direction cosines to give real components of a momentum
//
p.setX(p.getX()*p0);
p.setY(p.getY()*p0);
p.setZ(p.getZ()*p0);
track->setOrigin(orig);
track->getParticle().set4Momentum(p);
if(ptype == 2 || ptype == 3) track->getParticle().get4Momentum().setM(ElectronMass);
if(ptype == 8 || ptype == 9) track->getParticle().get4Momentum().setM(PionMass);
geantTracks.insert(track);
}
//
// in order to keep the construction/destruction under the control of ONE class,
// the hit-lists are created here, instead in the digitize() members of the detector
// classes. Only the digitized data is filled there.
//
// general structure of these labels:
// first data word is number of entries
// then follow 4 words per hit:
// int pointer to particle list (index in label geantTrackLabel=80048)
// -1 indicates added random hit
// int type of hit: 0=radiatorCher, 1=windowCher, 2=dEdx !! not for MCS1/MCS2/MCPD
// float x-coordinate of hit (pad units)
// float y-coordinate of hit (pad units)
// float energy loss dEdx !! only MCS1/MCS2/MCPD
// float direction theta at detector !! only MCS1/MCS2/MCPD
// float direction phi at detector !! only MCS1/MCS2/MCPD
CDetector::CDetectorId detId;
CMCDigiHit *digiHit;
int hitIndex; // start with index == 1 (hitlist comes from fortran)
int dataSize;
int dataLen;
// ------------------------------------ RICH-1 ------------------------------------
if ( (data = giveDataPointer(es, rich1HitLabel)) == 0 ) return False;
detId = CDetector::Rich1;
hitIndex = 0; // start with index == 1 (hitlist comes from fortran)
dataSize = 4; // number of entries per hit
dataLen = data[0]*dataSize; // data[0] is number of entries in list
for (i=1; i<dataLen; i+=dataSize) {
trackIndex = indexOfLastTrack + data[i]-1;
digiHit = new CMCDigiHit(); // create and fill MC digi hit
digiHit->setIndex ( hitIndex++ );
digiHit->setGeantTrack ( geantTracks(trackIndex) );
digiHit->setType ( data[i+1] );
digiHit->getGeantHit().setX( ((float *) data)[i+2] );
digiHit->getGeantHit().setY( ((float *) data)[i+3] );
digiHits[detId].insert(digiHit);
geantTracks[trackIndex]->addMCDigiHit(detId, digiHit);
}
// ------------------------------------ RICH-2 ------------------------------------
if ( (data = giveDataPointer(es, rich2HitLabel)) == 0 ) return False;
detId = CDetector::Rich2;
hitIndex = 0; // start with index == 1 (hitlist comes from fortran)
dataSize = 4; // number of entries per hit
dataLen = data[0]*dataSize; // data[0] is number of entries in list
for (i=1; i<dataLen; i+=dataSize) {
trackIndex = indexOfLastTrack + data[i]-1;
digiHit = new CMCDigiHit(); // create and fill MC digi hit
digiHit->setIndex ( hitIndex++ );
digiHit->setGeantTrack ( geantTracks(trackIndex) );
digiHit->setType ( data[i+1] );
digiHit->getGeantHit().setX( ((float *) data)[i+2] );
digiHit->getGeantHit().setY( ((float *) data)[i+3] );
digiHits[detId].insert(digiHit); // insert into hit pointer list
geantTracks[trackIndex]->addMCDigiHit(detId, digiHit); // add also to particle's hit list
}
// ------------------------------------ first Si detector ------------------------------------
if ( (data = giveDataPointer(es, sidc1HitLabel)) == 0 ) return False; // first Si detector
detId = CDetector::SiDC1;
hitIndex = 0; // start with index == 1 (hitlist comes from fortran)
dataSize = 6; // number of entries per hit
dataLen = data[0]*dataSize; // data[0] is number of entries in list
for (i=1; i<dataLen; i+=dataSize) {
trackIndex = indexOfLastTrack + data[i]-1;
digiHit = new CMCDigiHit(); // create and fill MC digi hit
digiHit->setIndex ( hitIndex++ );
digiHit->setGeantTrack ( geantTracks(trackIndex) );
digiHit->getGeantHit().setX ( ((float *) data)[i+1] );
digiHit->getGeantHit().setY ( ((float *) data)[i+2] );
digiHit->getGeantHit().setAmp( ((float *) data)[i+3] );
digiHit->setTheta ( ((float *) data)[i+4] );
digiHit->setPhi ( ((float *) data)[i+5] );
digiHits[detId].insert(digiHit); // insert into hit pointer list
geantTracks[trackIndex]->addMCDigiHit(detId, digiHit); // add also to particle's hit list
}
// ------------------------------------ second Si detector ------------------------------------
if ( (data = giveDataPointer(es, sidc2HitLabel)) == 0 ) return False;
detId = CDetector::SiDC2;
hitIndex = 0; // start with index == 1 (hitlist comes from fortran)
dataSize = 6; // number of entries per hit
dataLen = data[0]*dataSize; // data[0] is number of entries in list
for (i=1; i<dataLen; i+=dataSize) {
trackIndex = indexOfLastTrack + data[i]-1;
digiHit = new CMCDigiHit(); // create and fill MC digi hit
digiHit->setIndex ( hitIndex++ );
digiHit->setGeantTrack ( geantTracks(trackIndex) );
digiHit->getGeantHit().setX ( ((float *) data)[i+1] );
digiHit->getGeantHit().setY ( ((float *) data)[i+2] );
digiHit->getGeantHit().setAmp( ((float *) data)[i+3] );
digiHit->setTheta ( ((float *) data)[i+4] );
digiHit->setPhi ( ((float *) data)[i+5] );
digiHits[detId].insert(digiHit); // insert into hit pointer list
geantTracks[trackIndex]->addMCDigiHit(detId, digiHit); // add also to particle's hit list
}
// ------------------------------------ PadChamber ------------------------------------
if ( (data = giveDataPointer(es, padChamberHitLabel)) == 0 ) return False;
detId = CDetector::PadChamber;
hitIndex = 0; // start with index == 1 (hitlist comes from fortran)
dataSize = 6; // number of entries per hit
dataLen = data[0]*dataSize; // data[0] is number of entries in list
for (i=1; i<dataLen; i+=dataSize) {
trackIndex = indexOfLastTrack + data[i]-1;
digiHit = new CMCDigiHit(); // create and fill MC digi hit
digiHit->setIndex ( hitIndex++ );
digiHit->setGeantTrack ( geantTracks(trackIndex) );
digiHit->getGeantHit().setX ( ((float *) data)[i+1] );
digiHit->getGeantHit().setY ( ((float *) data)[i+2] );
digiHit->getGeantHit().setAmp( ((float *) data)[i+3] );
digiHit->setTheta ( ((float *) data)[i+4] );
digiHit->setPhi ( ((float *) data)[i+5] );
digiHits[detId].insert(digiHit); // insert into hit pointer list
geantTracks[trackIndex]->addMCDigiHit(detId, digiHit); // add also to particle's hit list
}
indexOfLastTrack = geantTracks.length();
return True;
}
int * CMCServer::giveDataPointer(CEventServer & es, int label)
{
const CLabel * mcpart = es.getLabel(label);
if (mcpart == 0) {
cerr << "CMCServer::unpack():\n";
cerr << " ERROR \n";
cerr << " label " << label << " not found" << endl;
return 0;
}
int *data = (int *) mcpart->getData();
if (data == 0) {
cerr << "CMCServer::unpack():\n"
<< " ERROR\n"
<< " no data found for label: " << label << endl;
return 0;
}
return data;
}
CVertex* CMCServer::getMCVertex()
{
if (mcvertex == 0)
mcvertex = new CVertex();
int i=0;
while ( i<geantTracks.length() &&
geantTracks(i)->getVolume() != "prim" ) // find first primary track
i++;
if (i < geantTracks.length() ) {
mcvertex->setX( geantTracks(i)->getOrigin().getX() );
mcvertex->setY( geantTracks(i)->getOrigin().getY() );
mcvertex->setZ( geantTracks(i)->getOrigin().getZ() + getSetup()->getTargetZInNewCoords() );
mcvertex->setChi2(-1000.); // 'very good'
mcvertex->setNTracks( geantTracks.length() );
} else {
mcvertex->setX( 9999. );
mcvertex->setY( 9999. );
mcvertex->setZ( 9999. );
mcvertex->setChi2(9999.); // 'very bad'
mcvertex->setNTracks( 0 );
}
return mcvertex;
}
CVertex* CMCServer::getMCVertexFromUserLabel(CEventServer &es)
{
if (mcvertex == 0)
mcvertex = new CVertex();
const CLabel * mcuser = (CLabel*) es.getLabel(setup->getUserLabel());
if (mcuser == 0) {
cerr << "CMCServer::getMCVertex():\n";
cerr << " ERROR \n";
cerr << " label " << setup->getUserLabel() << " not found" << endl;
return False;
}
float *data = (float *) mcuser->getData();
if (data == 0) {
cerr << "CMCServer::getMCVertex():\n"
<< "\tERROR\n"
<< "\tno data found for label: " << setup->getUserLabel() << endl;
return False;
}
// data[0] is dummy
mcvertex->setX(data[1]);
mcvertex->setY(data[2]);
mcvertex->setZ(data[3] + getSetup()->getTargetZInNewCoords() );
mcvertex->setChi2(-1000.); // 'very good'
mcvertex->setNTracks(10000);
return mcvertex;
}
CVertex* CMCServer::getMCADDeVertex()
{
if (mcADDevertex == 0)
mcADDevertex = new CVertex();
int i=0;
while ( i < geantTracks.length() &&
geantTracks(i)->getProcess() != "ADDe" ) // find first ADDe
i++;
if (i < geantTracks.length() ) {
mcADDevertex->setX( geantTracks(i)->getOrigin().getX() );
mcADDevertex->setY( geantTracks(i)->getOrigin().getY() );
mcADDevertex->setZ( geantTracks(i)->getOrigin().getZ() + getSetup()->getTargetZInNewCoords() );
mcADDevertex->setChi2(-1000.); // 'very good'
mcADDevertex->setNTracks( geantTracks.length() );
} else {
mcADDevertex->setX( 9999. );
mcADDevertex->setY( 9999. );
mcADDevertex->setZ( 9999. );
mcADDevertex->setChi2(9999.); // 'very bad'
mcADDevertex->setNTracks( 0 );
}
return mcADDevertex;
}
CBoolean CMCServer::alignGeantHitsToVertex(CVertex& mcVertex, CVertex& dataVertex, CSidc1& sidc1, CSidc2& sidc2, CPadChamber& padc)
{
CMCDigiHit *currentMcHit;
CRandom random;
double zmc, zevent, x, y, r, newx, newy, newr, dxVertex, dyVertex, dzVertex;
int i;
zmc = sidc1.getSetup()->getZPosition() - mcVertex.getZ();
zevent = sidc1.getSetup()->getZPosition() - dataVertex.getZ();
dxVertex = random.gauss(0.,0.0035); // to simulate Vertex resolution
dyVertex = random.gauss(0.,0.0035);
dzVertex = random.gauss(0.,0.01);
for (i=0; i<digiHits[0].length(); i++) { // SIDC1
currentMcHit = (digiHits[0])(i);
x = currentMcHit->getGeantHit().getX() - mcVertex.getX();
y = currentMcHit->getGeantHit().getY() - mcVertex.getY();
r = sqrt(x*x+y*y);
newr = (r/zmc)*(zevent + dzVertex); // r/z stays constant => theta constant
newx = (x/r)*newr + dataVertex.getX() + dxVertex; // r/x and r/y stays constant => phi constant
newy = (y/r)*newr + dataVertex.getY() + dyVertex;
currentMcHit->getGeantHit().setX(newx);
currentMcHit->getGeantHit().setY(newy);
}
zmc = sidc2.getSetup()->getZPosition() - mcVertex.getZ();
zevent = sidc2.getSetup()->getZPosition() - dataVertex.getZ();
for (i=0; i<digiHits[1].length(); i++) { // SIDC2
currentMcHit = (digiHits[1])(i);
x = currentMcHit->getGeantHit().getX() - mcVertex.getX();
y = currentMcHit->getGeantHit().getY() - mcVertex.getY();
r = sqrt(x*x+y*y);
newr = (r/zmc)*(zevent + dzVertex); // r/z stays constant => theta constant
newx = (x/r)*newr + dataVertex.getX() + dxVertex; // r/x and r/y stays constant => phi constant
newy = (y/r)*newr + dataVertex.getY() + dyVertex;
currentMcHit->getGeantHit().setX(newx);
currentMcHit->getGeantHit().setY(newy);
}
zmc = padc.getSetup()->getZPosition() + padc.getSetup()->getMCZOffset() - mcVertex.getZ();
zevent = padc.getSetup()->getZPosition() + padc.getSetup()->getMCZOffset() - dataVertex.getZ();
for (i=0; i<digiHits[4].length(); i++) { // PADC
currentMcHit = (digiHits[4])(i);
x = currentMcHit->getGeantHit().getX() - mcVertex.getX();
y = currentMcHit->getGeantHit().getY() - mcVertex.getY();
r = sqrt(x*x+y*y);
newr = (r/zmc)*zevent; // r/z stays constant => theta constant
newx = (x/r)*newr + dataVertex.getX(); // r/x and r/y stays constant => phi constant
newy = (y/r)*newr + dataVertex.getY();
currentMcHit->getGeantHit().setX(newx);
currentMcHit->getGeantHit().setY(newy);
}
return True;
}
float CMCServer::getWeight(CEventServer &es) const
{
const CLabel * mcuser = (CLabel*) es.getLabel(setup->getUserLabel());
if (mcuser == 0) {
cerr << "CMCServer::unpack():\n";
cerr << " ERROR \n";
cerr << " label " << setup->getUserLabel() << " not found" << endl;
return False;
}
int *data = (int *) mcuser->getData();
if (data == 0) {
cerr << "CMCServer::unpack():\n"
<< "\tERROR\n"
<< "\tno data found for label: " << setup->getUserLabel() << endl;
return False;
}
data++; // dummy
data++; // decay-type
float weight = *((float*)data);
return weight;
}
int CMCServer::getDecayMode(CEventServer &es) const
{
const CLabel * mcuser = (CLabel*) es.getLabel(setup->getUserLabel());
if (mcuser == 0) {
cerr << "CMCServer::unpack():\n";
cerr << " ERROR \n";
cerr << " label " << setup->getUserLabel() << " not found" << endl;
return False;
}
int *data = (int *) mcuser->getData();
if (data == 0) {
cerr << "CMCServer::unpack():\n"
<< "\tERROR\n"
<< "\tno data found for label: " << setup->getUserLabel() << endl;
return False;
}
data++; // dummy
int val = (int) *((float*)data);
return val;
}