Skip to content

Commit

Permalink
Merge pull request #235 from PADME-Experiment/feature/mctruth
Browse files Browse the repository at this point in the history 
Feature/mctruth
  • Loading branch information
@eleonardi
eleonardi authored Jul 22, 2021
2 parents 6f79bd8 + 78a87aa commit 3675c50
Show file tree
Hide file tree
Showing 39 changed files with 1,173 additions and 46 deletions.
25 changes: 23 additions & 2 deletions PadmeAnalysis/PadmeAnalysis.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
#include "TSACRecoEvent.hh"
#include "THEPVetoRecoEvent.hh"
#include "TRecoVHit.hh"
#include "TMCTruthEvent.hh"

#include "HistoSvc.hh"
#include "SACAnalysis.hh"
Expand Down Expand Up @@ -234,6 +235,7 @@ int main(Int_t argc, char **argv)
TRecoVClusCollection* fPVetoRecoCl =0;
TRecoVClusCollection* fEVetoRecoCl =0;
TRecoVClusCollection* fHEPVetoRecoCl =0;
TMCTruthEvent* fMCTruthEvent =0;

TTree::SetMaxTreeSize(190000000000);

Expand Down Expand Up @@ -315,12 +317,13 @@ int main(Int_t argc, char **argv)
} else if (branchName=="HEPVeto_Clusters") {
fHEPVetoRecoCl = new TRecoVClusCollection();
fRecoChain->SetBranchAddress(branchName.Data(),&fHEPVetoRecoCl);
} else if (branchName=="MCTruth") {
fMCTruthEvent = new TMCTruthEvent();
fRecoChain->SetBranchAddress(branchName.Data(),&fMCTruthEvent);
}

}



//////////// You come here if a Chain with >=0 events has been found
Int_t jevent = 0;

Expand Down Expand Up @@ -406,6 +409,7 @@ int main(Int_t argc, char **argv)
event->PVetoRecoCl =fPVetoRecoCl ;
event->EVetoRecoCl =fEVetoRecoCl ;
event->HEPVetoRecoCl =fHEPVetoRecoCl ;
event->MCTruthEvent =fMCTruthEvent ;

//UserAn->InitHistos();
//gTimeAn->InitHistos();
Expand All @@ -432,6 +436,23 @@ int main(Int_t argc, char **argv)
std::cout<<"----------------------------------------------------Run/Event n. = "<<fRecoEvent->GetRunNumber()<<" "<<fRecoEvent->GetEventNumber()<<std::endl;
}

//// Show MCTruth information (example)
//if (fMCTruthEvent) {
// printf("MCTruthEvent - Run %d Event %d Weight %8.3f Vertices %d\n",fMCTruthEvent->GetRunNumber(),fMCTruthEvent->GetEventNumber(),fMCTruthEvent->GetEventWeight(),fMCTruthEvent->GetNVertices());
// for(Int_t ii=0;ii<fMCTruthEvent->GetNVertices();ii++) {
// TMCVertex* vtx = fMCTruthEvent->Vertex(ii);
// printf("\tVertex %d Type %s Time %8.3f ns Position (%8.3f,%8.3f,%8.3f) mm Particles in %d out %d\n",ii,vtx->GetProcess().Data(),vtx->GetTime(),vtx->GetPosition().X(),vtx->GetPosition().Y(),vtx->GetPosition().Z(),vtx->GetNParticleIn(),vtx->GetNParticleOut());
// for(Int_t j=0;j<vtx->GetNParticleIn();j++) {
// TMCParticle* p = vtx->ParticleIn(j);
// printf("\t\tParticle In %d PDGCode %d Energy %8.3f MeV Momentum (%8.3f,%8.3f,%8.3f) MeV\n",j,p->GetPDGCode(),p->GetEnergy(),p->GetMomentum().X(),p->GetMomentum().Y(),p->GetMomentum().Z());
// }
// for(Int_t j=0;j<vtx->GetNParticleOut();j++) {
// TMCParticle* p = vtx->ParticleOut(j);
// printf("\t\tParticle Out %d PDGCode %d Energy %8.3f MeV Momentum (%8.3f,%8.3f,%8.3f) MeV\n",j,p->GetPDGCode(),p->GetEnergy(),p->GetMomentum().X(),p->GetMomentum().Y(),p->GetMomentum().Z());
// }
// }
//}

TTimeStamp timevent=fRecoEvent->GetEventTime();
utc_time=timevent.GetSec();

Expand Down
4 changes: 2 additions & 2 deletions PadmeAnalysis/UserAnalysis/include/PadmeAnalysisEvent.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#include "TSACRecoEvent.hh"
#include "THEPVetoRecoEvent.hh"
#include "TRecoVHit.hh"

#include "TMCTruthEvent.hh"

class PadmeAnalysisEvent{
public:
Expand All @@ -32,7 +32,7 @@ public:
TRecoVClusCollection* PVetoRecoCl ;
TRecoVClusCollection* EVetoRecoCl ;
TRecoVClusCollection* HEPVetoRecoCl ;

TMCTruthEvent* MCTruthEvent ;


};
Expand Down
16 changes: 1 addition & 15 deletions PadmeAnalysis/UserAnalysis/src/PadmeAnalysisEvent.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,20 +15,6 @@ PadmeAnalysisEvent::PadmeAnalysisEvent(){
PVetoRecoCl =0;
EVetoRecoCl =0;
HEPVetoRecoCl =0;















MCTruthEvent =0;

}
4 changes: 4 additions & 0 deletions PadmeMC/Beam/include/BeamGenerator.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ class G4Event;
class BeamMessenger;
class HistoManager;
class DetectorConstruction;
class MCTruthManager;

struct BeamPrimaryPositron
{
Expand Down Expand Up @@ -64,6 +65,7 @@ public:
private:

void GeneratePrimaryPositron();
void GenerateTargetPositron();
void CreateFinalStateUboson();
void CreateFinalStateThreeGamma();
void CreateFinalStateTwoGamma();
Expand All @@ -78,6 +80,8 @@ private:

DetectorConstruction* fDetector;

MCTruthManager* fMCTruthMgr;

BeamMessenger* fBeamMessenger;
HistoManager* fHistoManager;

Expand Down
7 changes: 7 additions & 0 deletions PadmeMC/Beam/include/BeamMessenger.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,5 +80,12 @@ private:
G4UIcmdWithADoubleAndUnit* fSetCalibRunCenterYCmd;
G4UIcmdWithADoubleAndUnit* fSetCalibRunRadiusCmd;

G4UIcmdWithADoubleAndUnit* fSetBeamTargetPosZCmd;
G4UIcmdWithADoubleAndUnit* fSetBeamTargetSigmaXCmd;
G4UIcmdWithADoubleAndUnit* fSetBeamTargetSigmaYCmd;
G4UIcmdWithADoubleAndUnit* fSetBeamTargetEmittanceXCmd;
G4UIcmdWithADoubleAndUnit* fSetBeamTargetEmittanceYCmd;
G4UIcmdWithADoubleAndUnit* fSetBeamTargetPathLengthCmd;

};
#endif
26 changes: 26 additions & 0 deletions PadmeMC/Beam/include/BeamParameters.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,24 @@ public:
G4double GetCalibRunRadius() { return fCalibRunRadius; }
void SetCalibRunRadius(G4double r) { fCalibRunRadius = r; }

G4double GetBeamTargetPosZ() { return fBeamTargetPosZ; }
void SetBeamTargetPosZ(G4double v) { fBeamTargetPosZ = v; }

G4double GetBeamTargetSigmaX() { return fBeamTargetSigmaX; }
void SetBeamTargetSigmaX(G4double v) { fBeamTargetSigmaX = v; }

G4double GetBeamTargetSigmaY() { return fBeamTargetSigmaY; }
void SetBeamTargetSigmaY(G4double v) { fBeamTargetSigmaY = v; }

G4double GetBeamTargetEmittanceX() { return fBeamTargetEmittanceX; }
void SetBeamTargetEmittanceX(G4double v) { fBeamTargetEmittanceX = v; }

G4double GetBeamTargetEmittanceY() { return fBeamTargetEmittanceY; }
void SetBeamTargetEmittanceY(G4double v) { fBeamTargetEmittanceY = v; }

G4double GetBeamTargetPathLength() { return fBeamTargetPathLength; }
void SetBeamTargetPathLength(G4double v) { fBeamTargetPathLength = v; }

private:

// Average number of positrons in each bunch
Expand Down Expand Up @@ -181,5 +199,13 @@ private:
G4double fCalibRunCenterY; // Y of center of circle
G4double fCalibRunRadius; // Radius of circle

// Beam distribution at target
G4double fBeamTargetPosZ; // Z position (1um in front of target)
G4double fBeamTargetSigmaX; // Sigma of X coordinate
G4double fBeamTargetSigmaY; // Sigma of Y coordinate
G4double fBeamTargetEmittanceX; // Emittance X component
G4double fBeamTargetEmittanceY; // Emittance Y component
G4double fBeamTargetPathLength; // Length of path from beam origin to target

};
#endif
104 changes: 96 additions & 8 deletions PadmeMC/Beam/src/BeamGenerator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include "BeamMessenger.hh"
#include "HistoManager.hh"
#include "DetectorConstruction.hh"
#include "MCTruthManager.hh"

BeamGenerator::BeamGenerator(DetectorConstruction* myDC)
:fDetector(myDC)
Expand All @@ -43,8 +44,13 @@ BeamGenerator::BeamGenerator(DetectorConstruction* myDC)

// Connect to BeamMessenger
fBeamMessenger = new BeamMessenger(this);

// Connect to Histogram manager
fHistoManager = HistoManager::GetInstance();

// Connect to MCTruth manager
fMCTruthMgr = MCTruthManager::GetInstance();

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Expand All @@ -59,6 +65,9 @@ BeamGenerator::~BeamGenerator()
void BeamGenerator::GenerateBeam(G4Event* anEvent)
{

// Reset MCTruth event
if (MCTruthManager::GetInstance()->IsEnabled()) MCTruthManager::GetInstance()->Clear();

fEvent = anEvent;

BeamParameters* bpar = BeamParameters::GetInstance();
Expand Down Expand Up @@ -99,8 +108,9 @@ void BeamGenerator::GenerateBeam(G4Event* anEvent)

for(int ib = 0; ib < nUbosonDecays; ib++) {

// Generate primary e+ which will decay to Uboson+gamma
GeneratePrimaryPositron();
// Generate primary e+ in front of Target which will decay to Uboson+gamma
//GeneratePrimaryPositron();
GenerateTargetPositron();

// Generate Uboson+gamma final state
CreateFinalStateUboson();
Expand All @@ -112,8 +122,9 @@ void BeamGenerator::GenerateBeam(G4Event* anEvent)
//*********************
for(int iggg = 0; iggg < nThreePhotonDecays; iggg++) {

// Generate primary e+ which will decay to three gammas
GeneratePrimaryPositron();
// Generate primary e+ in front of Target which will decay to three gamma
//GeneratePrimaryPositron();
GenerateTargetPositron();

// Generate gamma+gamma+gamma final state
CreateFinalStateThreeGamma();
Expand All @@ -126,14 +137,15 @@ void BeamGenerator::GenerateBeam(G4Event* anEvent)
//*********************
for(int iggg = 0; iggg < nTwoPhotonDecays; iggg++) {

// Generate primary e+ which will decay to two gammas
GeneratePrimaryPositron();
// Generate primary e+ in front of Target which will decay to two gamma
//GeneratePrimaryPositron();
GenerateTargetPositron();

// Generate gamma+gamma final state
CreateFinalStateTwoGamma();

}

//******************************************************
//General BG generator particles on the target per event
//******************************************************
Expand Down Expand Up @@ -229,7 +241,7 @@ void BeamGenerator::GeneratePrimaryPositron()

// Theta is gaussian with sigma from a phi-based combination of emittance along X and Y
G4double sigma_theta = bpar->GetBeamEmittanceX()*cos(phi)+bpar->GetBeamEmittanceX()*sin(phi);
G4double theta = G4RandGauss::shoot(0,sigma_theta);
G4double theta = G4RandGauss::shoot(0.,sigma_theta);

// Compute particle direction assuming beam is directed along Z (default direction)
G4double pX = sin(theta)*cos(phi);
Expand Down Expand Up @@ -263,12 +275,56 @@ void BeamGenerator::GeneratePrimaryPositron()

}

void BeamGenerator::GenerateTargetPositron()
{

// Generate a standard positron so that timing and energy are correct
GeneratePrimaryPositron();

BeamParameters* bpar = BeamParameters::GetInstance();

// If the beam is generated at Target front face, do nothing
if (bpar->GetBeamTargetPathLength() == 0.) return;

// Correct energy for mylar window crossing (to be implemented)

// Position positron in front of Target (need distribution)
G4double sigmaX = bpar->GetBeamTargetSigmaX();
G4double sigmaY = bpar->GetBeamTargetSigmaY();
G4double x = G4RandGauss::shoot(0.,sigmaX);
G4double y = G4RandGauss::shoot(0.,sigmaY);
G4double z = bpar->GetBeamTargetPosZ();
fPositron.pos = G4ThreeVector(x,y,z);

// Define direction (need emittance distribution)
G4double emitX = bpar->GetBeamTargetEmittanceX();
G4double emitY = bpar->GetBeamTargetEmittanceY();
G4double phi = twopi*G4UniformRand();
G4double sigmaTheta = emitX*cos(phi)+emitY*sin(phi);
G4double theta = G4RandGauss::shoot(0.,sigmaTheta);

// Compute direction vector
G4double pX = sin(theta)*cos(phi);
G4double pY = sin(theta)*sin(phi);
G4double pZ = cos(theta);
fPositron.dir = G4ThreeVector(pX,pY,pZ).unit();

// Define momentum vector
fPositron.P = sqrt(fPositron.E*fPositron.E-fPositron.m*fPositron.m);
fPositron.p = G4ThreeVector(fPositron.P*pX,fPositron.P*pY,fPositron.P*pZ);

// Correct positron time for path length from beam origin to Target
fPositron.t = fPositron.t+bpar->GetBeamTargetPathLength()/(c_light*fPositron.P/fPositron.E);

}

//void BeamGenerator::CreateFinalStateUboson(G4ParticleGun* pGun)
void BeamGenerator::CreateFinalStateUboson()
{

// Get mass of the Uboson
G4double ubosonM = BeamParameters::GetInstance()->GetUbosonMass();
//G4cout << "BeamGenerator - Uboson mass " << ubosonM/MeV << " MeV" << G4endl;

//Define the process: X->Uboson+gamma

Expand Down Expand Up @@ -402,6 +458,14 @@ void BeamGenerator::CreateFinalStateUboson()

// Add primary vertex to event
fEvent->AddPrimaryVertex(vtx);

// Store Uboson event in MCTruth
if (fMCTruthMgr->IsEnabled()) {
MCTruthVertex* tvtx = fMCTruthMgr->AddVertex("UBosonGamma",G4ThreeVector(Dx,Dy,Dz),Dt);
tvtx->AddParticleIn(G4ParticleTable::GetParticleTable()->FindParticle("e+")->GetPDGEncoding(),fPositron.E,fPositron.p);
tvtx->AddParticleOut(G4ParticleTable::GetParticleTable()->FindParticle("geantino")->GetPDGEncoding(),Up[0],G4ThreeVector(Up[1],Up[2],Up[3]));
tvtx->AddParticleOut(G4ParticleTable::GetParticleTable()->FindParticle("gamma")->GetPDGEncoding(),gp[0],G4ThreeVector(gp[1],gp[2],gp[3]));
}

}

Expand Down Expand Up @@ -455,6 +519,13 @@ void BeamGenerator::CreateFinalStateThreeGamma()
// Create primary vertex at decay point
G4PrimaryVertex* vtx = new G4PrimaryVertex(G4ThreeVector(Dx,Dy,Dz),Dt);

// Store vertex and primary positron in MCTruth
MCTruthVertex* tvtx;
if (fMCTruthMgr->IsEnabled()) {
tvtx = fMCTruthMgr->AddVertex("ThreeGamma",G4ThreeVector(Dx,Dy,Dz),Dt);
tvtx->AddParticleIn(G4ParticleTable::GetParticleTable()->FindParticle("e+")->GetPDGEncoding(),fPositron.E,fPositron.p);
}

//G4cout << "BeamGenerator 3Gamma - Positron Theta " << theta << " Phi " << phi
// << " Vertex " << G4ThreeVector(Dx,Dy,Dz) << " Time " << Dt << G4endl;

Expand Down Expand Up @@ -491,6 +562,11 @@ void BeamGenerator::CreateFinalStateThreeGamma()
gamma_p.x(),gamma_p.y(),gamma_p.z(),p[0]);
vtx->SetPrimary(gamma);

// Store gamma in MCTruth
if (fMCTruthMgr->IsEnabled()) {
tvtx->AddParticleOut(G4ParticleTable::GetParticleTable()->FindParticle("gamma")->GetPDGEncoding(),p[0],gamma_p);
}

//sum[0] += gamma_p.x();
//sum[1] += gamma_p.y();
//sum[2] += gamma_p.z();
Expand Down Expand Up @@ -559,6 +635,13 @@ void BeamGenerator::CreateFinalStateTwoGamma()
// Create primary vertex at decay point
G4PrimaryVertex* vtx = new G4PrimaryVertex(G4ThreeVector(Dx,Dy,Dz),Dt);

// Store vertex and primary positron in MCTruth
MCTruthVertex* tvtx;
if (fMCTruthMgr->IsEnabled()) {
tvtx = fMCTruthMgr->AddVertex("TwoGamma",G4ThreeVector(Dx,Dy,Dz),Dt);
tvtx->AddParticleIn(G4ParticleTable::GetParticleTable()->FindParticle("e+")->GetPDGEncoding(),fPositron.E,fPositron.p);
}

// Decode input line
std::istringstream iss(Line);

Expand Down Expand Up @@ -598,6 +681,11 @@ void BeamGenerator::CreateFinalStateTwoGamma()
gamma_p.x(),gamma_p.y(),gamma_p.z(),p[0]);
vtx->SetPrimary(gamma);

// Store gamma in MCTruth
if (fMCTruthMgr->IsEnabled()) {
tvtx->AddParticleOut(G4ParticleTable::GetParticleTable()->FindParticle("gamma")->GetPDGEncoding(),p[0],gamma_p);
}

}

// Add primary vertex to event
Expand Down
Loading

0 comments on commit 3675c50

Please sign in to comment.