Commit 8171d164 authored by Riccardo Longo's avatar Riccardo Longo
Browse files

MonteCarlo folder now available. Same as from Mike Phipps repository, except...

MonteCarlo folder now available. Same as from Mike Phipps repository, except for few adaptions for Geant4 compatibility with newer versions and CMakeLists.txt changes to have the folder as part of JCaPA
parent 96279e30
message("${BoldGreen}----------------------------")
message("Preparing MonteCarlo libraries...")
message("----------------------------${ColourReset}")
include_directories(${CMAKE_SOURCE_DIR})
#----------------------------------------------------------------------------
# Find Geant4 package, activating all available UI and Vis drivers by default
# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
# to build a batch mode only executable
#
option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
if(WITH_GEANT4_UIVIS)
find_package(Geant4 REQUIRED ui_all vis_all)
else()
find_package(Geant4 REQUIRED)
endif(WITH_GEANT4_UIVIS)
message("Geant4 use file variable is pointing at: ${Geant4_USE_FILE}")
include(${Geant4_USE_FILE})
#----------------------------------------------------------------------------
# Setup Geant4 include directories and compile definitions
# Setup include directory for this project
# Find ROOT
find_package (ROOT REQUIRED)
include_directories (${ROOT_INCLUDE_DIR})
# Add ROOTs header paths
#message(${CMAKE_MODULE_PATH} )
option(ANALYSIS_HH_ "Build with ROOT" 1)
add_definitions(-DG4UI_USE -DG4VIS_USE -DANALYSIS_HH_)
if(useROOT)
EXECUTE_PROCESS(COMMAND root-config --cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${ROOT_CXX_FLAGS}")
endif(useROOT)
include_directories(${PROJECT_SOURCE_DIR}/MonteCarlo/include ${ROOT_INCLUDE_DIR})
#----------------------------------------------------------------------------
# Locate sources and headers for this project
# NB: headers are included so they will show up in IDEs
#
file(GLOB sources ${PROJECT_SOURCE_DIR}/MonteCarlo/src/*.cc)
file(GLOB headers ${PROJECT_SOURCE_DIR}/MonteCarlo/include/*.hh)
#----------------------------------------------------------------------------
# Add the executable, and link it to the Geant4 libraries
#
add_executable(zdc zdc.cc ${sources} ${headers})
if(useROOT)
EXECUTE_PROCESS(COMMAND root-config --libs OUTPUT_VARIABLE ROOT_LD_FLAGS OUTPUT_STRIP_TRAILING_WHITESPACE)
endif(useROOT)
target_link_libraries(zdc ${Geant4_LIBRARIES} ${ROOT_LIBRARIES})
#----------------------------------------------------------------------------
# Copy all scripts to the build directory, i.e. the directory in which we
# build . This is so that we can run the executable directly because it
# relies on these scripts being in the current working directory.
#
set(EXAMPLE_SCRIPTS
#zdc.in
#zdc.out
init_vis.mac
run1.mac
run2.mac
vis.mac
)
foreach(_script ${EXAMPLE_SCRIPTS})
configure_file(
${PROJECT_SOURCE_DIR}/MonteCarlo/${_script}
${PROJECT_BINARY_DIR}/MonteCarlo/${_script}
COPYONLY
)
endforeach()
#----------------------------------------------------------------------------
# For internal Geant4 use - but has no effect if you build this
# example standalone
#
add_custom_target( DEPENDS zdc)
#----------------------------------------------------------------------------
# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
#
install(TARGETS zdc DESTINATION bin)
$Id:$
-------------------------------------------------------------------
=========================================================
atlasZDC
=========================================================
ChangeLog file
-----------------------
This file should be used to briefly summarize all major
modifications introduced in the code and keep track of
all tags.
----------------------------------------------------------
* Reverse chronological order (last date on top), please *
----------------------------------------------------------
07/31/18 M. Phipps
- Fixed global geometry bug. Fixed rod numbering issue. Cleaned up code
07/31/18 M. Phipps
- Created.
/run/verbose 2
/run/initialize
/control/execute vis.mac
/vis/open OGL 600x600-0+0
/vis/sceneHandler/create OGL
/vis/viewer/create ! ! 600x600-0+0
/vis/viewer/refresh
/vis/viewer/set/autoRefresh false
/vis/verbose errors
/vis/drawVolume
/vis/scene/create
/vis/scene/add/volume world
/vis/sceneHandler/attach
/hits/verbose 1
/vis/viewer/set/viewpointVector -1 0 0
/vis/viewer/set/lightsVector -1 0 0
/vis/viewer/set/style wireframe
/vis/viewer/set/auxiliaryEdge true
/vis/viewer/set/lineSegmentsPerCircle 100
/vis/scene/add/trajectories smooth
/tracking/storeTrajectory 2
/vis/scene/notifyHandlers
/vis/modeling/trajectories/create/drawByCharge
/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true
/vis/scene/notifyHandlers scene-0
/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/create/chargeFilter
/vis/filtering/trajectories/chargeFilter-0/add 1
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/chargeFilter-0/invert true
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/create/particleFilter
/vis/filtering/trajectories/particleFilter-0/add opticalphoton
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/particleFilter-0/add e-
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/particleFilter-0/add e+
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/particleFilter-0/add proton
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/particleFilter-0/add alpha
/vis/scene/notifyHandlers scene-0
/vis/filtering/trajectories/particleFilter-0/add GenericIon
/vis/scene/notifyHandlers scene-0
/vis/ogl/set/displayListLimit 5000000
/vis/scene/endOfEventAction accumulate
/vis/set/textColour green
/vis/set/textLayout right
/vis/scene/add/text2D 0.9 -.9 24 ! ! ZDC
/vis/scene/notifyHandlers
/vis/set/textLayout
/vis/set/textColour
/vis/scene/add/scale
/vis/scene/notifyHandlers
/vis/scene/add/axes
/vis/scene/notifyHandlers
/vis/scene/add/eventID
/vis/scene/notifyHandlers
/vis/scene/add/date
/vis/scene/notifyHandlers
/vis/set/colour red
/vis/set/lineWidth 2
/vis/scene/add/frame
/vis/scene/notifyHandlers
/vis/set/colour
/vis/set/lineWidth
/vis/geometry/set/visibility World 0 false
/vis/scene/notifyHandlers
/vis/viewer/set/hiddenMarker true
/vis/viewer/set/viewpointThetaPhi 90 90
/vis/viewer/set/autoRefresh true
/vis/viewer/refresh
/vis/verbose warnings
/run/initialize
/gun/particle proton
/gun/energy 450 GeV
# $Id: GNUmakefile 68058 2013-03-13 14:47:43Z gcosmo $
# --------------------------------------------------------------
# GNUmakefile for examples module. Gabriele Cosmo, 06/04/98.
# --------------------------------------------------------------
name := zdc
G4TARGET := $(name)
G4EXLIB := true
ifndef G4INSTALL
G4INSTALL = ../../..
endif
.PHONY: all
all: lib bin
include $(G4INSTALL)/config/binmake.gmk
visclean:
rm -f g4*.prim g4*.eps g4*.wrl
rm -f .DAWN_*
#include <iostream>
void analysis() {
char fileName[256];
sprintf(fileName,"/home/mike/ZDC/100Neutrons75GeV.root");
TFile *file = TFile::Open(fileName);
if (!file) {cerr << "file not open " << endl; return 1;}
TTree *hitTree = new TTree("Hit","Hit");
hitTree = (TTree*) file->Get("tree");
Int_t treeSize = (Int_t) hitTree->GetEntries();
/*
char h1Name[256];
sprintf(h1Name,"maxTB%d",z);
TH1F *h1 = new TH1F(h1Name,h1Name,27, -0.5, 26.5);
maxTBHisto.push_back(h1);
char h2Name[256];
sprintf(h2Name,"meanTime%d",z);
TH1F *h2 = new TH1F(h2Name,h2Name,115,10 , 18);
meanTimeHisto.push_back(h2);
*/
int ModNb = 0; int RadNb = 0; int Pid = 0; int EventNb = 0;
double EDep = 0.; double X = 0.; double Y = 0.; double Z = 0.; double Charge = 0.;
TCanvas *c1 = new TCanvas("c1","c1");
hitTree->Draw("EDep:X");
c1->SaveAs("EDepVsX.root");
c1->SaveAs("EDepVsX.png");
hitTree->Draw("EDep:Y");
c1->SaveAs("EDepVsY.root");
c1->SaveAs("EDepVsY.png");
hitTree->Draw("EDep:Z");
c1->SaveAs("EDepVsZ.root");
c1->SaveAs("EDepVsZ.png");
hitTree->Draw("X","Pid == 0");
c1->SaveAs("CherenkovVsX.root");
c1->SaveAs("CherenkovVsX.png");
hitTree->Draw("Y","Pid == 0");
c1->SaveAs("CherenkovVsY.root");
c1->SaveAs("CherenkovVsY.png");
hitTree->Draw("Z","Pid == 0");
c1->SaveAs("CherenkovVsZ.root");
c1->SaveAs("CherenkovVsZ.png");
hitTree->Draw("X","Pid == 22");
c1->SaveAs("GammaVsX.root");
c1->SaveAs("GammaVsX.png");
hitTree->Draw("X","Pid == 111");
c1->SaveAs("Pi0VsX.root");
c1->SaveAs("Pi0VsX.png");
hitTree->Draw("X","Pid == 211");
c1->SaveAs("Pi+VsX.root");
c1->SaveAs("Pi+VsX.png");
hitTree->Draw("X");
c1->SaveAs("TotalX.root");
c1->SaveAs("TotalX.png");
hitTree->Draw("Y");
c1->SaveAs("TotalY.root");
c1->SaveAs("TotalY.png");
hitTree->Draw("Z");
c1->SaveAs("TotalZ.root");
c1->SaveAs("TotalZ.png");
delete c1;
hitTree->SetBranchAddress("ModNb",&ModNb);
hitTree->SetBranchAddress("RadNb",&RadNb);
hitTree->SetBranchAddress("EDep",&EDep);
hitTree->SetBranchAddress("Pid",&Pid);
hitTree->SetBranchAddress("X",&X);
hitTree->SetBranchAddress("Y",&Y);
hitTree->SetBranchAddress("Z",&Z);
hitTree->SetBranchAddress("EventNo",&EventNb);
hitTree->SetBranchAddress("Charge",&Charge);
TH1F *photonHist = new TH1F("cherenkovYield","cherenkovYield",10,1000,100000);
TF1 *f1 = new TF1("f1","[0]+[1]*x+[2]*x*x",-45,45);
f1->SetParameter(0,0.162133);
f1->SetParameter(1,-4.27359e-05);
f1->SetParameter(2,6.55089e-05);
float prob = 0;
float aggregateprob = 0;
int prevEventNb = 0;
for (int i=0; i < treeSize; ++i) {
hitTree->GetEntry(i);
while (eventNb == prevEventNb) {
prob = f1->Eval(X);
aggregateprob += prob;
}
float totalProb = aggregateProb/nOptPhotons;
float eventYield = totalProb * nOptPhotons;
photonHist->Fill(eventYield);
std::cout << " yield " << eventYield << std::endl;
prevEventNb = eventNb;
}
delete hitTree;
}
}
#################################
###### Module Attributes ########
# currently 5 mods are maximum allowed
# mod types: 1 (EM-A), 2 (H1), 3 (standardMod), 4 (custom)
# if mod type 1-3 chosen, geometry from ATHENA simulation used
# if mod type 4 (custom) chosen, dimensions and materials must be set
# if you want just a block of absorber (eg for a testbeam) use type 4 and set core diameter to 0
# for now, casing of pixel mods (types 1 and 2) are assumed to have a 2 mm steel casing (in reality it's a little longer)
# Width of radiator gap assumed to be 0.5 mm greater than fiber diameter
# Units in mm
nModules: 4
mod1Type: 1
mod2Type: 1
mod3Type: 2
mod4Type: 2
#Custom categorizes (use same syntax for modN -- w/ 5 being the max # of mods supported)
mod1CasingThickness: 0.605
mod1NStripsPerGap: 70
# W and Pb are two currently enabled
mod1AbsorberMat: W
mod1AbsorberThickness: 3.
mod1AbsorberHeight: 180.0
# AbsorberWdith is only used if NStripsPerGap = 0 -- ie absorber only mode
mod1AbsorberWidth: 90.
# CoreDiameter should never be 0 -- even in absorber only mode
mod1CoreDiameter: 0.75
# Number of absorber layers assumed 1 fewer than # of radiators, unless 0 Radiators is entered, in which case the entire "module" is assumed to be an absorber block (like in test beam case when a Pb block is inserted)
mod1NRadiators: 5
# make sure you account for core + 2*cladding thickness + any empty buffer space
mod1RadiatorGapLength: 1.
mod1CoreIndexRefraction: 1.46
mod1Cladding: NO
# only if cladding enabled
mod1CladdingThickness: 0.
mod1CladdingIndexRefraction: 1.43
#Custom categorizes (use same syntax for modN -- w/ 5 being the max # of mods supported)
mod2CasingThickness: 0.605
mod2NStripsPerGap: 52
# W and Pb are two currently enabled
mod2AbsorberMat: W
mod2AbsorberThickness: 16.0
mod2AbsorberHeight: 180.0
# AbsorberWidth is only used if NStripsPerGap = 0 -- ie absorber only mode
mod2AbsorberWidth: 90.
mod2CoreDiameter: 1.5
# Number of absorber layers assumed 1 fewer than # of radiators
mod2NRadiators: 12
mod2RadiatorGapLength: 2
mod2CoreIndexRefraction: 1.46
mod2Cladding: YES
# only if cladding enabled
mod2CladdingThickness: 0.1
mod2CladdingIndexRefraction: 1.43
#Custom categorizes (use same syntax for modN -- w/ 5 being the max # of mods supported)
mod3CasingThickness: 0.605
mod3NStripsPerGap: 52
# W and Pb are two currently enabled
mod3AbsorberMat: W
mod3AbsorberThickness: 10.0
# AbsorberWidth is only used if NStripsPerGap = 0 -- ie absorber only mode
mod3AbsorberWidth: 90.
mod3AbsorberHeight: 180.0
mod3CoreDiameter: 1.5
# Number of absorber layers assumed 1 fewer than # of radiators
mod3NRadiators: 12
mod3RadiatorGapLength: 2
mod3CoreIndexRefraction: 1.46
mod3Cladding: NO
# only if cladding enabled
mod3CladdingThickness: 0.1
mod3CladdingIndexRefraction: 1.43
#Custom categorizes (use same syntax for modN -- w/ 5 being the max # of mods supported)
mod4CasingThickness: 0.605
mod4NStripsPerGap: 52
# W and Pb are two currently enabled
mod4AbsorberMat: W
mod4AbsorberThickness: 10.0
mod4AbsorberHeight: 180.0
# AbsorberWidth is only used if NStripsPerGap = 0 -- ie absorber only mode
mod4AbsorberWidth: 90.
mod4CoreDiameter: 1.5
# Number of absorber layers assumed 1 fewer than # of radiators
mod4NRadiators: 12
mod4RadiatorGapLength: 2
mod4CoreIndexRefraction: 1.46
mod4Cladding: NO
# only if cladding enabled
mod4CladdingThickness: 0.1
mod4CladdingIndexRefraction: 1.43
#################################
########### Physics #############
# FTFP_Bert and QGSP_Bert are two lists currently enabled
physicsList: FTFP_BERT
# enter these in units of nm
cherenkovMinWavelength: 250
cherenkovMaxWavelength: 600
# FTFP_BERT used as standard hadronic list. Other physics modules used include: G4RadioactiveDecayPhysics, G4HadronElasticPhysicsXS, G4StoppingPhysics, and G4IonPhysics
#################################
############ Output #############
# TTree with every step in radiator recorded
fullEventListing: NO
# TTree with # of Cherenkovs captured
cherenkovYield: NO
# standard: 1 channel/module. Non-standard output: every strip read out individually
standardSegmentation: YES
#################################
############ Debugging ##########
checkOverlaps: false
\ No newline at end of file
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: ActionInitialization.hh 68058 2013-03-13 14:47:43Z gcosmo $
//
/// \file ActionInitialization.hh
/// \brief Definition of the ActionInitialization class
#ifndef ActionInitialization_h
#define ActionInitialization_h 1
#include "G4VUserActionInitialization.hh"
#include "SharedData.hh"
/// Action initialization class.
class ActionInitialization : public G4VUserActionInitialization
{
public:
ActionInitialization(SharedData *sharedData);
virtual ~ActionInitialization();
virtual void BuildForMaster() const;
virtual void Build() const;
private:
SharedData *fSharedData;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// January 22, 2012
// Yakov Kulinich
// Stony Brook & BNL
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef CherenkovHit_h
#define CherenkovHit_h 1
#include "G4VHit.hh"
#include "G4THitsCollection.hh"
#include "G4Allocator.hh"
#include "G4ThreeVector.hh"
#include "G4ParticleDefinition.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class CherenkovHit : public G4VHit
{
public:
CherenkovHit();
~CherenkovHit();
CherenkovHit(const CherenkovHit&);
const CherenkovHit& operator=(const CherenkovHit&);
G4int operator==(const CherenkovHit&) const;
inline void* operator new(size_t);
inline void operator delete(void*);
void Draw();
void Print();
public:
void setTrackID (G4int track) { trackID = track; };
void setModNb (G4int mod) { modNb = mod; };
void setRadNb (G4int rad) { radNb = rad; };
void setEdep (G4double de) { edep = de; };
void setPos (G4ThreeVector xyz) { pos = xyz; };
void setMomentum (G4ThreeVector mom) { momentum = mom;};
void setParticle (G4ParticleDefinition *part) {particle = part;};
void setEnergy (G4double e) {energy = e;};
void setCharge (G4double c) {charge = c;};
void setEventNo (G4int eventNo) {fEventNo = eventNo;}
void setVelocity (G4double vel) {velocity = vel;}
void setBeta (G4double b) {beta = b;}
G4int getTrackID() { return trackID; };
G4int getModNb() { return modNb; };
G4int getRadNb() { return radNb; };
G4ParticleDefinition* getParticle() { return particle; };
G4double getEdep() { return edep; };
G4double getEnergy() { return energy; };
G4ThreeVector getPos() { return pos; };
G4ThreeVector getMomentum() { return momentum; };
G4double getCharge() {return charge;};
G4int getEventNo() {return fEventNo;}
G4double getVelocity() {return velocity;}
G4double getBeta() {return beta;}
private:
G4int trackID;
G4int modNb;
G4int radNb;
G4double edep;