Newer
Older
/*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* *
* Original Author : Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : January 2013 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with Beam: *
* User can enter various parameter, such as emittance or use ASCII or root *
* TH1F distribution *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
* *
*****************************************************************************/
// STL
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <sstream>
#include <cmath>
#include <vector>
// NPL header
#include "NPBeam.h"
#include "NPFunction.h"
#include "NPOptionManager.h"
// ROOT Header
#include "TDirectory.h"
// Use CLHEP System of unit and Physical Constant
#include "NPGlobalSystemOfUnits.h"
#include "NPPhysicalConstants.h"
using namespace NPL;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
fEnergy = 0;
fSigmaEnergy = -1 ;
fMeanX = 0 ;
fMeanY = 0 ;
fSigmaX = -1;
fSigmaY = -0;
fMeanThetaX = 0 ;
fMeanPhiY = 0 ;
fSigmaThetaX = -1 ;
fSigmaPhiY = -1 ;
fTargetSize = 0 ;
fEffectiveTargetSize = 0 ;
fTargetThickness = 0 ;
fEffectiveTargetThickness = 0 ;
fTargetAngle = 0 ;
fTargetZ = 0 ;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
// do that to avoid warning from multiple Hist with same name...
int offset = 0;
while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
++offset;
fEnergyHist = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Beam::Beam(string isotope){
SetUp(isotope);
fEnergy = 0;
fSigmaEnergy = -1 ;
fMeanX = 0 ;
fMeanY = 0 ;
fSigmaX = -1;
fSigmaY = -0;
fMeanThetaX = 0 ;
fMeanPhiY = 0 ;
fSigmaThetaX = -1 ;
fSigmaPhiY = -1 ;
fTargetSize = 0 ;
fEffectiveTargetSize = 0 ;
fTargetThickness = 0 ;
fEffectiveTargetThickness = 0 ;
fTargetAngle = 0 ;
fTargetZ = 0 ;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
// case of user given distribution
// do that to avoid warning from multiple Hist with same name...
int offset = 0;
while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
++offset;
fEnergyHist = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
delete fEnergyHist ;
delete fXThetaXHist ;
delete fYPhiYHist ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::ReadConfigurationFile(string Path){
////////General Reading needs////////
string LineBuffer;
string DataBuffer;
//////////////////////////////////////////////////////////////////////////////////////////
ifstream ReactionFile;
bool ReadingStatus = false ;
bool check_BeamName = false ;
bool check_Energy = false ;
bool check_SigmaEnergy = false ;
bool check_MeanX = false ;
bool check_MeanY = false ;
bool check_SigmaX = false ;
bool check_SigmaY = false ;
bool check_MeanThetaX = false ;
bool check_MeanPhiY = false ;
bool check_SigmaThetaX = false ;
bool check_SigmaPhiY = false ;
bool check_EnergyProfilePath = false ;
bool check_XThetaXPath = false ;
bool check_YPhiYPath = false ;
bool check_AllEnergy = false;
bool check_AllEmittance = false;
string GlobalPath = getenv("NPTOOL");
string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
ReactionFile.open(StandardPath.c_str());
if (ReactionFile.is_open()) {cout << "Reading Reaction File " << Path << endl ;}
// In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
else{
ReactionFile.open( Path.c_str() );
if(ReactionFile.is_open()) { if(fVerboseLevel==1) cout << "Reading Reaction File " << Path << endl;}
else {cout << "Reaction File " << Path << " not found" << endl;exit(1);}
while (!ReactionFile.eof()) {
//Pick-up next line
getline(ReactionFile, LineBuffer);
if (LineBuffer.compare(0, 4, "Beam") == 0) {
if(fVerboseLevel==1) cout << "Beam Found" << endl ;
ReadingStatus = true ;
}
while(ReadingStatus){
ReactionFile >> DataBuffer;
//Search for comment Symbol %
if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
else if (DataBuffer == "Particle=") {
check_BeamName = true ;
ReactionFile >> DataBuffer;
SetUp(DataBuffer);
if(fVerboseLevel==1) cout << "Beam Particle: " << GetName() << endl;
else if (DataBuffer == "ExcitationEnergy=") {
ReactionFile >> DataBuffer;
fExcitationEnergy = atof(DataBuffer.c_str()) * MeV;
if(fVerboseLevel==1) cout << "Excitation Energy: " << fExcitationEnergy << " MeV" << endl;
}
else if (DataBuffer == "Energy=") {
check_Energy = true ;
ReactionFile >> DataBuffer;
fEnergy = atof(DataBuffer.c_str()) * MeV;
if(fVerboseLevel==1) cout << "Beam Energy: " << fEnergy / MeV << " MeV" << endl;
}
else if (DataBuffer == "SigmaEnergy=") {
check_SigmaEnergy = true ;
ReactionFile >> DataBuffer;
fSigmaEnergy= atof(DataBuffer.c_str()) * MeV;
if(fVerboseLevel==1) cout << "Beam Energy Sigma: " << fSigmaEnergy / MeV << " MeV" << endl;
}
else if (DataBuffer=="MeanX=") {
check_MeanX = true ;
ReactionFile >> DataBuffer;
fMeanX = atof(DataBuffer.c_str()) * mm;
if(fVerboseLevel==1) cout << "Mean X: " << fMeanX / mm << " mm" << endl;
}
else if (DataBuffer=="MeanY=") {
check_MeanY = true ;
ReactionFile >> DataBuffer;
fMeanY = atof(DataBuffer.c_str()) * mm;
if(fVerboseLevel==1) cout << "Mean Y: " << fMeanY / mm << " mm" << endl;
}
else if (DataBuffer=="SigmaX=") {
check_SigmaX = true ;
ReactionFile >> DataBuffer;
fSigmaX = atof(DataBuffer.c_str()) * mm;
if(fVerboseLevel==1) cout << "Sigma X: " << fSigmaX / mm << " mm" << endl;
}
else if (DataBuffer=="SigmaY=") {
check_SigmaY = true ;
ReactionFile >> DataBuffer;
fSigmaY = atof(DataBuffer.c_str()) * mm;
if(fVerboseLevel==1) cout << "Sigma Y: " << fSigmaY / mm << " mm" << endl;
}
else if (DataBuffer == "MeanThetaX=" ) {
check_MeanThetaX = true ;
ReactionFile >> DataBuffer;
fMeanThetaX = atof(DataBuffer.c_str()) * deg;
if(fVerboseLevel==1) cout << "Mean Theta X: " << fMeanThetaX / deg << " deg" << endl;
}
else if (DataBuffer == "MeanPhiY=") {
check_MeanPhiY = true ;
ReactionFile >> DataBuffer;
fMeanPhiY = atof(DataBuffer.c_str()) * deg;
if(fVerboseLevel==1) cout << "Mean Phi Y: " << fMeanPhiY / deg << " deg" << endl;
}
else if (DataBuffer == "SigmaThetaX=" ) {
check_SigmaThetaX = true ;
ReactionFile >> DataBuffer;
fSigmaThetaX = atof(DataBuffer.c_str()) * deg;
if(fVerboseLevel==1) cout << "Sigma Theta X: " << fSigmaThetaX / deg << " deg" << endl;
}
else if (DataBuffer == "SigmaPhiY=") {
check_SigmaPhiY = true ;
ReactionFile >> DataBuffer;
fSigmaPhiY = atof(DataBuffer.c_str()) * deg;
if(fVerboseLevel==1) cout << "Sigma Phi Y: " << fSigmaPhiY / deg << " deg" << endl;
}
else if (DataBuffer == "EnergyProfilePath=") {
check_EnergyProfilePath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Energy profile file: " << FileName << endl;
SetEnergyHist( Read1DProfile(FileName, HistName ));
}
else if (DataBuffer == "XThetaXProfilePath=") {
check_XThetaXPath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading X-ThetaX profile file: " << FileName << endl;
SetXThetaXHist(Read2DProfile(FileName, HistName ) );
}
else if (DataBuffer == "YPhiYProfilePath=") {
check_YPhiYPath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Y-ThetaY profile file: " << FileName << endl;
SetYPhiYHist( Read2DProfile(FileName, HistName ));
}
///////////////////////////////////////////////////
// If no Beam Token and no comment, toggle out
else{
cout << "ERROR : Wrong Token Sequence: Getting out " << endl ;
exit(1);
}
///////////////////////////////////////////////////
if( ( check_MeanX && check_MeanY && check_SigmaX && check_SigmaY && check_SigmaThetaX && check_SigmaPhiY && check_MeanThetaX && check_MeanPhiY) || ( check_XThetaXPath && check_YPhiYPath ) ){
check_AllEmittance = true ;
}
if( ( check_Energy && check_SigmaEnergy ) || ( check_EnergyProfilePath ) ){
check_AllEnergy = true ;
}
// If all Token found toggle out
if( check_BeamName && check_AllEnergy && check_AllEmittance )
ReadingStatus = false ;
}
}
if( !check_BeamName || !check_AllEnergy || !check_AllEnergy ){
cout << "ERROR : Token Sequence Incomplete, Beam definition could not be Fonctionnal" << endl ;
exit(1);
}
if(fVerboseLevel==1) cout << "///////////////////////////////////////////////////" << endl << endl ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY ){
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
X = Y = 1000000*cm;
if(fSigmaEnergy!=-1)
E = gRandom->Gaus(fEnergy,fSigmaEnergy);
else
E = fEnergyHist->GetRandom();
if(fSigmaX!=-1){
// Shoot within the target unless target size is null (no limit)
while(sqrt(X*X+Y*Y>fEffectiveTargetSize) || fEffectiveTargetSize == 0){
NPL::RandomGaussian2D(fMeanX, fMeanThetaX, fSigmaX, fSigmaThetaX, X, ThetaX);
NPL::RandomGaussian2D(fMeanY, fMeanPhiY, fSigmaY, fSigmaPhiY, Y, PhiY);
}
}
else{
while(sqrt(X*X+Y*Y>fEffectiveTargetSize) || fEffectiveTargetSize == 0){
fXThetaXHist->GetRandom2(X,ThetaX);
fYPhiYHist->GetRandom2(Y,PhiY);
}
}
Z = fTargetZ + fEffectiveTargetThickness*(gRandom->Uniform()-0.5);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::Print() const {
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetSize(double TargetSize){
fTargetSize = TargetSize;
fEffectiveTargetSize = fTargetSize*cos(fTargetAngle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetThickness(double TargetThickness){
fTargetThickness = TargetThickness;
fEffectiveTargetThickness = fTargetThickness/cos(fTargetAngle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetAngle(double TargetAngle){
fTargetAngle = TargetAngle;
fEffectiveTargetSize = fTargetSize*cos(fTargetAngle);
fEffectiveTargetThickness = fTargetThickness/cos(fTargetAngle);
}