00001
00002
00003
00004
00005
00006
00007
00008 #include "ProcessCountingSvc.hh"
00009
00010 #include <string>
00011 #include <iostream>
00012 #include <fstream>
00013 #include <sstream>
00014
00015 #include "myglobals.hh"
00016 #include "G4Track.hh"
00017 #include "G4Step.hh"
00018 #include "G4StepPoint.hh"
00019 #include "G4VProcess.hh"
00020
00021 #include "MyString2Anything.hh"
00022
00023 #include "MyRoot.hh"
00024
00025 ProcessCountingSvc* ProcessCountingSvc::fProcessCountingSvc = 0;
00026
00027 ProcessCountingSvc::ProcessCountingSvc()
00028 {
00029 if (fProcessCountingSvc){
00030 G4Exception("ProcessCountingSvc::ProcessCountingSvc()","Run0031",
00031 FatalException, "ProcessCountingSvc constructed twice.");
00032 }
00033 fProcessCountingSvc = this;
00034 ReSet();
00035 }
00036
00037 ProcessCountingSvc::~ProcessCountingSvc()
00038 {
00039 printf("~ProcessCountingSvc\n");
00040 }
00041
00042 ProcessCountingSvc* ProcessCountingSvc::GetProcessCountingSvc(){
00043 if ( !fProcessCountingSvc ){
00044 fProcessCountingSvc = new ProcessCountingSvc;
00045 }
00046 return fProcessCountingSvc;
00047 }
00048
00049 void ProcessCountingSvc::Initialize(){
00050 m_TrackIDs.clear();
00051 m_nSteps = 0;
00052 m_pid.clear();
00053 m_tid.clear();
00054 m_nSec.clear();
00055 m_time.clear();
00056 m_stepL.clear();
00057 m_prePx.clear();
00058 m_prePy.clear();
00059 m_prePz.clear();
00060 m_postPx.clear();
00061 m_postPy.clear();
00062 m_postPz.clear();
00063 m_dTheta.clear();
00064 m_dE.clear();
00065 m_edepT.clear();
00066 m_edepN.clear();
00067 m_edepI.clear();
00068 m_e.clear();
00069 m_preX.clear();
00070 m_preY.clear();
00071 m_preZ.clear();
00072 m_postX.clear();
00073 m_postY.clear();
00074 m_postZ.clear();
00075 m_particleName.clear();
00076 m_charge.clear();
00077 m_process.clear();
00078 m_ASDI_msc.clear();
00079 m_ASDI_hPairProd.clear();
00080 m_ASDI_ProtonInelastic.clear();
00081 m_ASDI_Transportation.clear();
00082 m_ASDI_hBrems.clear();
00083 m_ASDI_eBrem.clear();
00084 m_ASDI_eIoni.clear();
00085 m_PSDI_eIoni.clear();
00086 m_PSDI_eBrem.clear();
00087 m_PSDI_msc.clear();
00088 m_PSDI_hPairProd.clear();
00089 m_PSDI_ProtonInelastic.clear();
00090 m_PSDI_hIoni.clear();
00091 m_PSDI_hadElastic.clear();
00092 m_PSDI_ionIoni.clear();
00093 m_PSDI_Decay.clear();
00094 m_PSDI_Transportation.clear();
00095 m_PSDI_hBrems.clear();
00096 m_volume.clear();
00097 VolName = "";
00098 }
00099
00100 void ProcessCountingSvc::SetBranch(){
00101 if (!m_Switch) return;
00102 MyRoot* myRoot = MyRoot::GetMyRoot();
00103 if( flag_nSteps ) myRoot->SetBranch("ProcessCounting_nSteps", &m_nSteps);
00104 if( flag_pid ) myRoot->SetBranch("ProcessCounting_pid", &m_pid);
00105 if( flag_tid ) myRoot->SetBranch("ProcessCounting_tid", &m_tid);
00106 if( flag_nSec ) myRoot->SetBranch("ProcessCounting_nSec", &m_nSec);
00107 if( flag_time ) myRoot->SetBranch("ProcessCounting_time", &m_time);
00108 if( flag_stepL ) myRoot->SetBranch("ProcessCounting_stepL", &m_stepL);
00109 if( flag_prePx ) myRoot->SetBranch("ProcessCounting_prePx", &m_prePx);
00110 if( flag_prePy ) myRoot->SetBranch("ProcessCounting_prePy", &m_prePy);
00111 if( flag_prePz ) myRoot->SetBranch("ProcessCounting_prePz", &m_prePz);
00112 if( flag_postPx ) myRoot->SetBranch("ProcessCounting_postPx", &m_postPx);
00113 if( flag_postPy ) myRoot->SetBranch("ProcessCounting_postPy", &m_postPy);
00114 if( flag_postPz ) myRoot->SetBranch("ProcessCounting_postPz", &m_postPz);
00115 if( flag_dTheta ) myRoot->SetBranch("ProcessCounting_dTheta", &m_dTheta);
00116 if( flag_dE ) myRoot->SetBranch("ProcessCounting_dE", &m_dE);
00117 if( flag_edepT ) myRoot->SetBranch("ProcessCounting_edepT", &m_edepT);
00118 if( flag_edepN ) myRoot->SetBranch("ProcessCounting_edepN", &m_edepN);
00119 if( flag_edepI ) myRoot->SetBranch("ProcessCounting_edepI", &m_edepI);
00120 if( flag_e ) myRoot->SetBranch("ProcessCounting_e", &m_e);
00121 if( flag_preX ) myRoot->SetBranch("ProcessCounting_preX", &m_preX);
00122 if( flag_preY ) myRoot->SetBranch("ProcessCounting_preY", &m_preY);
00123 if( flag_preZ ) myRoot->SetBranch("ProcessCounting_preZ", &m_preZ);
00124 if( flag_postX ) myRoot->SetBranch("ProcessCounting_postX", &m_postX);
00125 if( flag_postY ) myRoot->SetBranch("ProcessCounting_postY", &m_postY);
00126 if( flag_postZ ) myRoot->SetBranch("ProcessCounting_postZ", &m_postZ);
00127 if( flag_charge ) myRoot->SetBranch("ProcessCounting_charge", &m_charge);
00128 if( flag_particleName ) myRoot->SetBranch("ProcessCounting_particleName", &m_particleName);
00129 if( flag_process ) myRoot->SetBranch("ProcessCounting_process", &m_process);
00130 if( flag_ASDI_msc ) myRoot->SetBranch("ProcessCounting_ASDI_msc", &m_ASDI_msc);
00131 if( flag_ASDI_hPairProd ) myRoot->SetBranch("ProcessCounting_ASDI_hPairProd", &m_ASDI_hPairProd);
00132 if( flag_ASDI_ProtonInelastic ) myRoot->SetBranch("ProcessCounting_ASDI_ProtonInelastic", &m_ASDI_ProtonInelastic);
00133 if( flag_ASDI_hIoni ) myRoot->SetBranch("ProcessCounting_ASDI_hIoni", &m_ASDI_hIoni);
00134 if( flag_ASDI_hadElastic ) myRoot->SetBranch("ProcessCounting_ASDI_hadElastic", &m_ASDI_hadElastic);
00135 if( flag_ASDI_ionIoni ) myRoot->SetBranch("ProcessCounting_ASDI_ionIoni", &m_ASDI_ionIoni);
00136 if( flag_ASDI_Decay ) myRoot->SetBranch("ProcessCounting_ASDI_Decay", &m_ASDI_Decay);
00137 if( flag_ASDI_Transportation ) myRoot->SetBranch("ProcessCounting_ASDI_Transportation", &m_ASDI_Transportation);
00138 if( flag_ASDI_hBrems ) myRoot->SetBranch("ProcessCounting_ASDI_hBrems", &m_ASDI_hBrems);
00139 if( flag_ASDI_eBrem ) myRoot->SetBranch("ProcessCounting_ASDI_eBrem", &m_ASDI_eBrem);
00140 if( flag_ASDI_eIoni ) myRoot->SetBranch("ProcessCounting_ASDI_eIoni", &m_ASDI_eIoni);
00141 if( flag_PSDI_eBrem ) myRoot->SetBranch("ProcessCounting_PSDI_eBrem", &m_PSDI_eBrem);
00142 if( flag_PSDI_eIoni ) myRoot->SetBranch("ProcessCounting_PSDI_eIoni", &m_PSDI_eIoni);
00143 if( flag_PSDI_msc ) myRoot->SetBranch("ProcessCounting_PSDI_msc", &m_PSDI_msc);
00144 if( flag_PSDI_hPairProd ) myRoot->SetBranch("ProcessCounting_PSDI_hPairProd", &m_PSDI_hPairProd);
00145 if( flag_PSDI_ProtonInelastic ) myRoot->SetBranch("ProcessCounting_PSDI_ProtonInelastic", &m_PSDI_ProtonInelastic);
00146 if( flag_PSDI_hIoni ) myRoot->SetBranch("ProcessCounting_PSDI_hIoni", &m_PSDI_hIoni);
00147 if( flag_PSDI_hadElastic ) myRoot->SetBranch("ProcessCounting_PSDI_hadElastic", &m_PSDI_hadElastic);
00148 if( flag_PSDI_ionIoni ) myRoot->SetBranch("ProcessCounting_PSDI_ionIoni", &m_PSDI_ionIoni);
00149 if( flag_PSDI_Decay ) myRoot->SetBranch("ProcessCounting_PSDI_Decay", &m_PSDI_Decay);
00150 if( flag_PSDI_Transportation ) myRoot->SetBranch("ProcessCounting_PSDI_Transportation", &m_PSDI_Transportation);
00151 if( flag_PSDI_hBrems ) myRoot->SetBranch("ProcessCounting_PSDI_hBrems", &m_PSDI_hBrems);
00152 if( flag_volume ) myRoot->SetBranch("ProcessCounting_volume", &m_volume);
00153 }
00154
00155 void ProcessCountingSvc::ReadOutputCard(G4String filename){
00156 ReSet();
00157 std::ifstream fin_card(filename);
00158 if(!fin_card){
00159 std::cout<<"In ProcessCountingSvc::ReadOutputCard, cannot open "<<filename<<"!!!"<<std::endl;
00160 G4Exception("ProcessCountingSvc::ReadOutputCard()",
00161 "InvalidSetup", FatalException,
00162 "cannot find output card");
00163 }
00164 std::stringstream buf_card;
00165 std::string s_card;
00166 int n_output_section_symbol = 0;
00167 int n_filter_section_symbol = 0;
00168 while(getline(fin_card,s_card)){
00169 buf_card.str("");
00170 buf_card.clear();
00171 buf_card<<s_card;
00172
00173
00174 const char* c_card = s_card.c_str();
00175 int length = strlen(c_card);
00176 int offset = 0;
00177 for ( ; offset < length; offset++ ){
00178 if ( c_card[offset] != ' ' ) break;
00179 }
00180 if ( c_card[offset] == '#' || (c_card[offset] == '/' && c_card[offset+1] == '/') || length - offset == 0 ){
00181 continue;
00182 }
00183
00184 std::string name;
00185 buf_card>>name;
00186 if ( n_output_section_symbol == 0 ){
00187 if ( name == "PROCESSCOUNTING_SECTION" ){
00188 n_output_section_symbol++;
00189 }
00190 }
00191 else if ( n_output_section_symbol == 1 ){
00192 if ( name == "PROCESSCOUNTING_SECTION" ){
00193 n_output_section_symbol++;
00194 }
00195 else if( name == "nSteps" ) flag_nSteps = true;
00196 else if( name == "pid" ) flag_pid = true;
00197 else if( name == "tid" ) flag_tid = true;
00198 else if( name == "nSec" ) flag_nSec = true;
00199 else if( name == "time" ) {flag_time = true; buf_card>>unitName_time; unit_time = MyString2Anything::get_U(unitName_time);}
00200 else if( name == "time" ) {flag_time = true; buf_card>>unitName_time; unit_time = MyString2Anything::get_U(unitName_time);}
00201 else if( name == "stepL" ) {flag_stepL = true; buf_card>>unitName_stepL; unit_stepL = MyString2Anything::get_U(unitName_stepL);}
00202 else if( name == "prePx" ) {flag_prePx = true; buf_card>>unitName_prePx; unit_prePx = MyString2Anything::get_U(unitName_prePx);}
00203 else if( name == "prePy" ) {flag_prePy = true; buf_card>>unitName_prePy; unit_prePy = MyString2Anything::get_U(unitName_prePy);}
00204 else if( name == "prePz" ) {flag_prePz = true; buf_card>>unitName_prePz; unit_prePz = MyString2Anything::get_U(unitName_prePz);}
00205 else if( name == "postPx" ) {flag_postPx = true; buf_card>>unitName_postPx; unit_postPx = MyString2Anything::get_U(unitName_postPx);}
00206 else if( name == "postPy" ) {flag_postPy = true; buf_card>>unitName_postPy; unit_postPy = MyString2Anything::get_U(unitName_postPy);}
00207 else if( name == "postPz" ) {flag_postPz = true; buf_card>>unitName_postPz; unit_postPz = MyString2Anything::get_U(unitName_postPz);}
00208 else if( name == "dTheta" ) {flag_dTheta = true; buf_card>>unitName_dTheta; unit_dTheta = MyString2Anything::get_U(unitName_dTheta);}
00209 else if( name == "dE" ) {flag_dE = true; buf_card>>unitName_dE; unit_dE = MyString2Anything::get_U(unitName_dE);}
00210 else if( name == "edepT" ) {flag_edepT = true; buf_card>>unitName_edepT; unit_edepT = MyString2Anything::get_U(unitName_edepT);}
00211 else if( name == "edepN" ) {flag_edepN = true; buf_card>>unitName_edepN; unit_edepN = MyString2Anything::get_U(unitName_edepN);}
00212 else if( name == "edepI" ) {flag_edepI = true; buf_card>>unitName_edepI; unit_edepI = MyString2Anything::get_U(unitName_edepI);}
00213 else if( name == "e" ) {flag_e = true; buf_card>>unitName_e; unit_e = MyString2Anything::get_U(unitName_e);}
00214 else if( name == "preX" ) {flag_preX = true; buf_card>>unitName_preX; unit_preX = MyString2Anything::get_U(unitName_preX);}
00215 else if( name == "preY" ) {flag_preY = true; buf_card>>unitName_preY; unit_preY = MyString2Anything::get_U(unitName_preY);}
00216 else if( name == "preZ" ) {flag_preZ = true; buf_card>>unitName_preZ; unit_preZ = MyString2Anything::get_U(unitName_preZ);}
00217 else if( name == "postX" ) {flag_postX = true; buf_card>>unitName_postX; unit_postX = MyString2Anything::get_U(unitName_postX);}
00218 else if( name == "postY" ) {flag_postY = true; buf_card>>unitName_postY; unit_postY = MyString2Anything::get_U(unitName_postY);}
00219 else if( name == "postZ" ) {flag_postZ = true; buf_card>>unitName_postZ; unit_postZ = MyString2Anything::get_U(unitName_postZ);}
00220 else if( name == "particleName" ) flag_particleName = true;
00221 else if( name == "charge" ) flag_charge = true;
00222 else if( name == "process" ) flag_process = true;
00223 else if( name == "ASDI_msc" ) flag_ASDI_msc = true;
00224 else if( name == "ASDI_hPairProd" ) flag_ASDI_hPairProd = true;
00225 else if( name == "ASDI_ProtonInelastic" ) flag_ASDI_ProtonInelastic = true;
00226 else if( name == "ASDI_hIoni" ) flag_ASDI_hIoni = true;
00227 else if( name == "ASDI_hadElastic" ) flag_ASDI_hadElastic = true;
00228 else if( name == "ASDI_ionIoni" ) flag_ASDI_ionIoni = true;
00229 else if( name == "ASDI_Decay" ) flag_ASDI_Decay = true;
00230 else if( name == "ASDI_Transportation" ) flag_ASDI_Transportation = true;
00231 else if( name == "ASDI_hBrems" ) flag_ASDI_hBrems = true;
00232 else if( name == "ASDI_eBrem" ) flag_ASDI_eBrem = true;
00233 else if( name == "ASDI_eIoni" ) flag_ASDI_eIoni = true;
00234 else if( name == "PSDI_eBrem" ) flag_PSDI_eBrem = true;
00235 else if( name == "PSDI_eIoni" ) flag_PSDI_eIoni = true;
00236 else if( name == "PSDI_msc" ) flag_PSDI_msc = true;
00237 else if( name == "PSDI_hPairProd" ) flag_PSDI_hPairProd = true;
00238 else if( name == "PSDI_ProtonInelastic" ) flag_PSDI_ProtonInelastic = true;
00239 else if( name == "PSDI_hIoni" ) flag_PSDI_hIoni = true;
00240 else if( name == "PSDI_hadElastic" ) flag_PSDI_hadElastic = true;
00241 else if( name == "PSDI_ionIoni" ) flag_PSDI_ionIoni = true;
00242 else if( name == "PSDI_Decay" ) flag_PSDI_Decay = true;
00243 else if( name == "PSDI_Transportation" ) flag_PSDI_Transportation = true;
00244 else if( name == "PSDI_hBrems" ) flag_PSDI_hBrems = true;
00245 else if( name == "volume" ) flag_volume = true;
00246 else{
00247 std::cout<<"In ProcessCountingSvc::ReadOutputCard, unknown name: "<<name<<" in file "<<filename<<std::endl;
00248 std::cout<<"Will ignore this line!"<<std::endl;
00249 }
00250 }
00251
00252 if ( n_filter_section_symbol == 0 ){
00253 if ( name == "PROCESSCOUNTINGFILTER_SECTION" ){
00254 n_filter_section_symbol++;
00255 }
00256 }
00257 else if ( n_filter_section_symbol == 1 ){
00258 if ( name == "PROCESSCOUNTINGFILTER_SECTION" ){
00259 n_filter_section_symbol++;
00260 }
00261 else if( name == "Switch" ) m_Switch = true;
00262 else if( name == "nTracks" ) buf_card>>m_maxnTracks;
00263 else if( name == "nSteps" ) buf_card>>m_maxnSteps;
00264 else if( name == "Volume" ){
00265 G4String VolName;
00266 while (buf_card>>VolName){
00267 if(VolName.substr(0,1) == "#" || VolName.substr(0,2) == "//") break;
00268 m_Volumes.push_back(VolName);
00269 }
00270 }
00271 else{
00272 G4double para;
00273 std::string unit;
00274 buf_card>>para>>unit;
00275 para *= MyString2Anything::get_U(unit);
00276 if( name == "minp" ) m_minp = para;
00277 else if( name == "mine" ) m_mine = para;
00278 else if( name == "mint" ) m_mint = para;
00279 else if( name == "maxt" ) m_maxt = para;
00280 else{
00281 std::cout<<"In ProcessCountingSvc::ReadOutputCard, unknown name: "<<name<<" in file "<<filename<<std::endl;
00282 std::cout<<"Will ignore this line!"<<std::endl;
00283 }
00284 }
00285 }
00286
00287 if ( n_output_section_symbol > 1 && n_filter_section_symbol > 1 ){
00288 break;
00289 }
00290 }
00291 buf_card.str("");
00292 buf_card.clear();
00293 if ( n_output_section_symbol <= 1 ){
00294 std::cout<<"*****************WARNING********************"<<std::endl;
00295 std::cout<<"In ProcessCountingSvc::ReadOutputCard, failed to find enough section seperators for output in file "<<filename<<std::endl;
00296 std::cout<<"Will use default settings."<<std::endl;
00297 std::cout<<"********************************************"<<std::endl;
00298 }
00299 if ( n_filter_section_symbol<= 1 ){
00300 std::cout<<"*****************WARNING********************"<<std::endl;
00301 std::cout<<"In ProcessCountingSvc::ReadOutputCard, failed to find enough section seperators for filter in file "<<filename<<std::endl;
00302 std::cout<<"Will use default settings."<<std::endl;
00303 std::cout<<"********************************************"<<std::endl;
00304 }
00305 fin_card.close();
00306 ShowOutCard();
00307 }
00308
00309 void ProcessCountingSvc::ReSet(){
00310 flag_nSteps = false;
00311 flag_pid = false;
00312 flag_tid = false;
00313 flag_nSec = false;
00314 flag_time = false;
00315 flag_stepL = false;
00316 flag_prePx = false;
00317 flag_prePy = false;
00318 flag_prePz = false;
00319 flag_postPx = false;
00320 flag_postPy = false;
00321 flag_postPz = false;
00322 flag_dTheta = false;
00323 flag_dE = false;
00324 flag_edepT = false;
00325 flag_edepN = false;
00326 flag_edepI = false;
00327 flag_e = false;
00328 flag_preX = false;
00329 flag_preY = false;
00330 flag_preZ = false;
00331 flag_postX = false;
00332 flag_postY = false;
00333 flag_postZ = false;
00334 flag_charge = false;
00335 flag_particleName = false;
00336 flag_process = false;
00337 flag_ASDI_msc = false;
00338 flag_ASDI_hPairProd = false;
00339 flag_ASDI_ProtonInelastic = false;
00340 flag_ASDI_hIoni = false;
00341 flag_ASDI_hadElastic = false;
00342 flag_ASDI_ionIoni = false;
00343 flag_ASDI_Decay = false;
00344 flag_ASDI_Transportation = false;
00345 flag_ASDI_hBrems = false;
00346 flag_ASDI_eBrem = false;
00347 flag_ASDI_eIoni = false;
00348 flag_PSDI_eIoni = false;
00349 flag_PSDI_eBrem = false;
00350 flag_PSDI_msc = false;
00351 flag_PSDI_hPairProd = false;
00352 flag_PSDI_ProtonInelastic = false;
00353 flag_PSDI_hIoni = false;
00354 flag_PSDI_hadElastic = false;
00355 flag_PSDI_ionIoni = false;
00356 flag_PSDI_Decay = false;
00357 flag_PSDI_Transportation = false;
00358 flag_PSDI_hBrems = false;
00359 flag_volume = false;
00360 m_Switch = false;
00361 m_maxnTracks = 0;
00362 m_minp = 0;
00363 m_mine = 0;
00364 m_maxnSteps = 0;
00365 m_mint = 0;
00366 m_maxt = 0;
00367 m_Volumes.clear();
00368 unitName_time ="s";
00369 unitName_stepL ="cm";
00370 unitName_prePx ="GeV";
00371 unitName_prePy ="GeV";
00372 unitName_prePz ="GeV";
00373 unitName_postPx ="GeV";
00374 unitName_postPy ="GeV";
00375 unitName_postPz ="GeV";
00376 unitName_dTheta ="rad";
00377 unitName_dE ="GeV";
00378 unitName_edepT ="GeV";
00379 unitName_edepN ="GeV";
00380 unitName_edepI ="GeV";
00381 unitName_e ="GeV";
00382 unitName_preX ="cm";
00383 unitName_preY ="cm";
00384 unitName_preZ ="cm";
00385 unitName_postX ="cm";
00386 unitName_postY ="cm";
00387 unitName_postZ ="cm";
00388 unit_time =s;
00389 unit_stepL =cm;
00390 unit_prePx =GeV;
00391 unit_prePy =GeV;
00392 unit_prePz =GeV;
00393 unit_postPx =GeV;
00394 unit_postPy =GeV;
00395 unit_postPz =GeV;
00396 unit_dTheta =rad;
00397 unit_dE =GeV;
00398 unit_edepT =GeV;
00399 unit_edepN =GeV;
00400 unit_edepI =GeV;
00401 unit_e =GeV;
00402 unit_preX =cm;
00403 unit_preY =cm;
00404 unit_preZ =cm;
00405 unit_postX =cm;
00406 unit_postY =cm;
00407 unit_postZ =cm;
00408 }
00409
00410 void ProcessCountingSvc::ShowOutCard(){
00411 std::cout<<"*************************Output settings for ProcessCountingSvc***************************"<<std::endl;
00412 std::cout<<"output nSteps? "<<(flag_nSteps?" yes":" no")<<std::endl;
00413 std::cout<<"output pid? "<<(flag_pid?" yes":" no")<<std::endl;
00414 std::cout<<"output tid? "<<(flag_tid?" yes":" no")<<std::endl;
00415 std::cout<<"output nSec? "<<(flag_nSec?" yes":" no")<<std::endl;
00416 std::cout<<"output time? "<<(flag_time?" yes":" no")<<", unit: "<<unitName_time<<std::endl;
00417 std::cout<<"output stepL? "<<(flag_stepL?" yes":" no")<<", unit: "<<unitName_stepL<<std::endl;
00418 std::cout<<"output prePx? "<<(flag_prePx?" yes":" no")<<", unit: "<<unitName_prePx<<std::endl;
00419 std::cout<<"output prePy? "<<(flag_prePy?" yes":" no")<<", unit: "<<unitName_prePy<<std::endl;
00420 std::cout<<"output prePz? "<<(flag_prePz?" yes":" no")<<", unit: "<<unitName_prePz<<std::endl;
00421 std::cout<<"output postPx? "<<(flag_postPx?" yes":" no")<<", unit: "<<unitName_postPx<<std::endl;
00422 std::cout<<"output postPy? "<<(flag_postPy?" yes":" no")<<", unit: "<<unitName_postPy<<std::endl;
00423 std::cout<<"output postPz? "<<(flag_postPz?" yes":" no")<<", unit: "<<unitName_postPz<<std::endl;
00424 std::cout<<"output dTheta? "<<(flag_dTheta?" yes":" no")<<", unit: "<<unitName_dTheta<<std::endl;
00425 std::cout<<"output dE? "<<(flag_dE?" yes":" no")<<", unit: "<<unitName_dE<<std::endl;
00426 std::cout<<"output edepT? "<<(flag_edepT?" yes":" no")<<", unit: "<<unitName_edepT<<std::endl;
00427 std::cout<<"output edepN? "<<(flag_edepN?" yes":" no")<<", unit: "<<unitName_edepN<<std::endl;
00428 std::cout<<"output edepI? "<<(flag_edepI?" yes":" no")<<", unit: "<<unitName_edepI<<std::endl;
00429 std::cout<<"output e? "<<(flag_e?" yes":" no")<<", unit: "<<unitName_e<<std::endl;
00430 std::cout<<"output preX? "<<(flag_preX?" yes":" no")<<", unit: "<<unitName_preX<<std::endl;
00431 std::cout<<"output preY? "<<(flag_preY?" yes":" no")<<", unit: "<<unitName_preY<<std::endl;
00432 std::cout<<"output preZ? "<<(flag_preZ?" yes":" no")<<", unit: "<<unitName_preZ<<std::endl;
00433 std::cout<<"output postX? "<<(flag_postX?" yes":" no")<<", unit: "<<unitName_postX<<std::endl;
00434 std::cout<<"output postY? "<<(flag_postY?" yes":" no")<<", unit: "<<unitName_postY<<std::endl;
00435 std::cout<<"output postZ? "<<(flag_postZ?" yes":" no")<<", unit: "<<unitName_postZ<<std::endl;
00436 std::cout<<"output particleName? "<<(flag_particleName?" yes":" no")<<std::endl;
00437 std::cout<<"output charge? "<<(flag_charge?" yes":" no")<<std::endl;
00438 std::cout<<"output process? "<<(flag_process?" yes":" no")<<std::endl;
00439 std::cout<<"output ASDI_msc? "<<(flag_ASDI_msc?" yes":" no")<<std::endl;
00440 std::cout<<"output ASDI_hPairProd?"<<(flag_ASDI_hPairProd?" yes":" no")<<std::endl;
00441 std::cout<<"output ASDI_ProtonInelastic?"<<(flag_ASDI_ProtonInelastic?" yes":" no")<<std::endl;
00442 std::cout<<"output ASDI_hIoni? "<<(flag_ASDI_hIoni?" yes":" no")<<std::endl;
00443 std::cout<<"output ASDI_hadElastic? "<<(flag_ASDI_hadElastic?" yes":" no")<<std::endl;
00444 std::cout<<"output ASDI_ionIoni? "<<(flag_ASDI_ionIoni?" yes":" no")<<std::endl;
00445 std::cout<<"output ASDI_Decay? "<<(flag_ASDI_Decay?" yes":" no")<<std::endl;
00446 std::cout<<"output ASDI_Transportation? "<<(flag_ASDI_Transportation?" yes":" no")<<std::endl;
00447 std::cout<<"output ASDI_hBrems? "<<(flag_ASDI_hBrems?" yes":" no")<<std::endl;
00448 std::cout<<"output ASDI_eBrem? "<<(flag_ASDI_eBrem?" yes":" no")<<std::endl;
00449 std::cout<<"output ASDI_eIoni? "<<(flag_ASDI_eIoni?" yes":" no")<<std::endl;
00450 std::cout<<"output PSDI_eIoni? "<<(flag_PSDI_eIoni?" yes":" no")<<std::endl;
00451 std::cout<<"output PSDI_eBrem? "<<(flag_PSDI_eBrem?" yes":" no")<<std::endl;
00452 std::cout<<"output PSDI_msc? "<<(flag_PSDI_msc?" yes":" no")<<std::endl;
00453 std::cout<<"output PSDI_hPairProd?"<<(flag_PSDI_hPairProd?" yes":" no")<<std::endl;
00454 std::cout<<"output PSDI_ProtonInelastic?"<<(flag_PSDI_ProtonInelastic?" yes":" no")<<std::endl;
00455 std::cout<<"output PSDI_hIoni? "<<(flag_PSDI_hIoni?" yes":" no")<<std::endl;
00456 std::cout<<"output PSDI_hadElastic? "<<(flag_PSDI_hadElastic?" yes":" no")<<std::endl;
00457 std::cout<<"output PSDI_ionIoni? "<<(flag_PSDI_ionIoni?" yes":" no")<<std::endl;
00458 std::cout<<"output PSDI_Decay? "<<(flag_PSDI_Decay?" yes":" no")<<std::endl;
00459 std::cout<<"output PSDI_Transportation? "<<(flag_PSDI_Transportation?" yes":" no")<<std::endl;
00460 std::cout<<"output PSDI_hBrems? "<<(flag_PSDI_hBrems?" yes":" no")<<std::endl;
00461 std::cout<<"output volume? "<<(flag_volume?" yes":" no")<<std::endl;
00462 std::cout<<"Switch on? "<<(m_Switch?"yes":"no")<<std::endl;
00463 std::cout<<"minp = "<<m_minp/MeV<<"MeV"<<std::endl;
00464 std::cout<<"mine = "<<m_mine/MeV<<"MeV"<<std::endl;
00465 std::cout<<"maxnSteps = "<<m_maxnSteps<<std::endl;
00466 std::cout<<"maxnTracks = "<<m_maxnTracks<<std::endl;
00467 std::cout<<"mint = "<<m_mint/ns<<"ns"<<std::endl;
00468 std::cout<<"maxt = "<<m_maxt/ns<<"ns"<<std::endl;
00469 std::cout<<m_Volumes.size()<<" Volumes Selected:";
00470 for ( int i = 0; i< m_Volumes.size(); i++ ){
00471 std::cout<<" "<<m_Volumes[i];
00472 }
00473 std::cout<<std::endl;
00474 std::cout<<"******************************************************************************"<<std::endl;
00475 }
00476
00477 void ProcessCountingSvc::InitialStep(const G4Step* aStep){
00478 std::cout<<"In ProcessCountingSvc::InitialStep"<<std::endl;
00479
00480 const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();
00481 VolName = touchable->GetVolume(0)->GetName();
00482
00483 CheckTrack(aStep);
00484 CheckFilter(aStep);
00485 if (!PASSEDFILTER) return;
00486
00487
00488 m_ASDI_msc.push_back(0);
00489 m_ASDI_hPairProd.push_back(0);
00490 m_ASDI_ProtonInelastic.push_back(0);
00491 m_ASDI_hIoni.push_back(0);
00492 m_ASDI_hadElastic.push_back(0);
00493 m_ASDI_ionIoni.push_back(0);
00494 m_ASDI_Decay.push_back(0);
00495 m_ASDI_Transportation.push_back(0);
00496 m_ASDI_hBrems.push_back(0);
00497 m_ASDI_eBrem.push_back(0);
00498 m_ASDI_eIoni.push_back(0);
00499 m_PSDI_eBrem.push_back(0);
00500 m_PSDI_eIoni.push_back(0);
00501 m_PSDI_msc.push_back(0);
00502 m_PSDI_hPairProd.push_back(0);
00503 m_PSDI_ProtonInelastic.push_back(0);
00504 m_PSDI_hIoni.push_back(0);
00505 m_PSDI_hadElastic.push_back(0);
00506 m_PSDI_ionIoni.push_back(0);
00507 m_PSDI_Decay.push_back(0);
00508 m_PSDI_Transportation.push_back(0);
00509 m_PSDI_hBrems.push_back(0);
00510 }
00511
00512 void ProcessCountingSvc::AddASDI(G4String name){
00513 if (!PASSEDFILTER) return;
00514 int i = m_nSteps;
00515 if ( name == "msc" || name == "CoulombScat" ) m_ASDI_msc[i]++;
00516 else if ( name == "eIoni" ) m_ASDI_eIoni[i]++;
00517 else if ( name == "eBrem" ) m_ASDI_eBrem[i]++;
00518 else if ( name == "hIoni" ) m_ASDI_hIoni[i]++;
00519 else if ( name == "hadElastic" ) m_ASDI_hadElastic[i]++;
00520 else if ( name == "ionIoni" ) m_ASDI_ionIoni[i]++;
00521 else if ( name == "Decay" ) m_ASDI_Decay[i]++;
00522 else if ( name == "Transportation" ) m_ASDI_Transportation[i]++;
00523 else if ( name == "hBrems" ) m_ASDI_hBrems[i]++;
00524 else if ( name == "hPairProd" ) m_ASDI_hPairProd[i]++;
00525 else if ( name == "ProtonInelastic" ) m_ASDI_ProtonInelastic[i]++;
00526
00527 }
00528
00529 void ProcessCountingSvc::AddPSDI(G4String name){
00530 if (!PASSEDFILTER) return;
00531 int i = m_nSteps;
00532 if ( name == "msc" || name == "CoulombScat" ) m_PSDI_msc[i]++;
00533 else if ( name == "eIoni" ) m_PSDI_eIoni[i]++;
00534 else if ( name == "eBrem" ) m_PSDI_eBrem[i]++;
00535 else if ( name == "hIoni" ) m_PSDI_hIoni[i]++;
00536 else if ( name == "hadElastic" ) m_PSDI_hadElastic[i]++;
00537 else if ( name == "ionIoni" ) m_PSDI_ionIoni[i]++;
00538 else if ( name == "Decay" ) m_PSDI_Decay[i]++;
00539 else if ( name == "Transportation" ) m_PSDI_Transportation[i]++;
00540 else if ( name == "hBrems" ) m_PSDI_hBrems[i]++;
00541 else if ( name == "hPairProd" ) m_PSDI_hPairProd[i]++;
00542 else if ( name == "ProtonInelastic" ) m_PSDI_ProtonInelastic[i]++;
00543
00544 }
00545
00546 void ProcessCountingSvc::SetValue(const G4Step* aStep){
00547 const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();
00548 VolName = touchable->GetVolume(0)->GetName();
00549 CheckFilter(aStep);
00550 if (!PASSEDFILTER) return;
00551
00552
00553
00554
00555 G4Track* aTrack = aStep->GetTrack() ;
00556 G4int trackID= aTrack->GetTrackID();
00557 G4int pid = aTrack->GetParticleDefinition()->GetPDGEncoding();
00558 G4int charge = aTrack->GetParticleDefinition()->GetPDGCharge();
00559 std::string particleName = aTrack->GetParticleDefinition()->GetParticleName();
00560
00561
00562 G4StepPoint* prePoint = aStep->GetPreStepPoint() ;
00563 G4ThreeVector pointIn_pos = prePoint->GetPosition();
00564 G4ThreeVector pointIn_mom = prePoint->GetMomentum();
00565 G4double pointIn_e = prePoint->GetTotalEnergy();
00566 G4double pointIn_pa = pointIn_mom.mag();
00567 G4double pointIn_time = prePoint->GetGlobalTime();
00568 G4StepPoint* postPoint = aStep->GetPostStepPoint() ;
00569 G4double pointOut_e = postPoint->GetTotalEnergy();
00570 G4ThreeVector pointOut_pos = postPoint->GetPosition();
00571 G4ThreeVector pointOut_mom = postPoint->GetMomentum();
00572 G4String processName;
00573 const G4VProcess* process = postPoint->GetProcessDefinedStep();
00574 if (process) {
00575 processName = process->GetProcessName();
00576 }
00577 else{
00578 processName = "NULL";
00579 }
00580
00581
00582 G4double stepL = aStep->GetStepLength();
00583 G4int nSec;
00584 const std::vector<const G4Track*>* pVTracks = aStep->GetSecondaryInCurrentStep();
00585 if (pVTracks){
00586 nSec = pVTracks->size();
00587 }
00588 else{
00589 nSec = -1;
00590 }
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 G4double deltaTheta = pointOut_mom.theta(pointIn_mom);
00615 G4double deltaE = pointIn_e - pointOut_e;
00616 G4double edep_total = aStep->GetTotalEnergyDeposit();
00617 G4double edep_nonIoni = aStep->GetNonIonizingEnergyDeposit();
00618 G4double edep_Ioni = edep_total - edep_nonIoni;
00619
00620 m_nSteps++;
00621 if(flag_pid) m_pid.push_back(pid);
00622 if(flag_tid) m_tid.push_back(trackID);
00623 if(flag_nSec) m_nSec.push_back(nSec);
00624 if(flag_time) m_time.push_back(pointIn_time/unit_time);
00625 if(flag_stepL) m_stepL.push_back(stepL/unit_stepL);
00626 if(flag_prePx) m_prePx.push_back(pointIn_mom.x()/unit_prePx);
00627 if(flag_prePy) m_prePy.push_back(pointIn_mom.y()/unit_prePy);
00628 if(flag_prePz) m_prePz.push_back(pointIn_mom.z()/unit_prePz);
00629 if(flag_postPx) m_postPx.push_back(pointOut_mom.x()/unit_postPx);
00630 if(flag_postPy) m_postPy.push_back(pointOut_mom.y()/unit_postPy);
00631 if(flag_postPz) m_postPz.push_back(pointOut_mom.z()/unit_postPz);
00632 if(flag_dTheta) m_dTheta.push_back(deltaTheta/unit_dTheta);
00633 if(flag_dE) m_dE.push_back(deltaE/unit_dE);
00634 if(flag_edepT) m_edepT.push_back(edep_total/unit_edepT);
00635 if(flag_edepN) m_edepN.push_back(edep_Ioni/unit_edepN);
00636 if(flag_edepI) m_edepI.push_back(edep_nonIoni/unit_edepI);
00637 if(flag_e) m_e.push_back(pointIn_e/unit_e);
00638 if(flag_preX) m_preX.push_back(pointIn_pos.x()/unit_preX);
00639 if(flag_preY) m_preY.push_back(pointIn_pos.y()/unit_preY);
00640 if(flag_preZ) m_preZ.push_back(pointIn_pos.z()/unit_preZ);
00641 if(flag_postX) m_postX.push_back(pointOut_pos.x()/unit_postX);
00642 if(flag_postY) m_postY.push_back(pointOut_pos.y()/unit_postY);
00643 if(flag_postZ) m_postZ.push_back(pointOut_pos.z()/unit_postZ);
00644 if(flag_charge) m_charge.push_back(charge);
00645 if(flag_particleName) m_particleName.push_back(particleName);
00646 if(flag_process) m_process.push_back(processName);
00647 if(flag_volume) m_volume.push_back(VolName);
00648 }
00649
00650 void ProcessCountingSvc::CheckFilter(const G4Step* aStep){
00651 PASSEDFILTER = true;
00652
00653 if (!m_Switch) {PASSEDFILTER = false; return;}
00654
00655
00656 G4StepPoint* prePoint = aStep->GetPreStepPoint() ;
00657 G4ThreeVector pointIn_mom = prePoint->GetMomentum();
00658 G4double pointIn_e = prePoint->GetTotalEnergy();
00659 G4double pointIn_pa = pointIn_mom.mag();
00660 G4double pointIn_time = prePoint->GetGlobalTime();
00661
00662
00663
00664 if (m_maxnTracks&&m_TrackIDs.size()>m_maxnTracks) {PASSEDFILTER = false; return;}
00665
00666 if (m_maxnSteps&&m_nSteps>=m_maxnSteps) {PASSEDFILTER = false; return;}
00667
00668 if (m_minp&&pointIn_pa<m_minp) {PASSEDFILTER = false; return;}
00669
00670 if (m_mine&&pointIn_e<m_mine) {PASSEDFILTER = false; return;}
00671
00672 if (m_mint&&pointIn_time<m_mint) {PASSEDFILTER = false; return;}
00673 if (m_maxt&&pointIn_time>m_maxt) {PASSEDFILTER = false; return;}
00674 if (m_Volumes.size()>0){
00675 PASSEDFILTER = false;
00676 for ( int i = 0; i< m_Volumes.size(); i++ ){
00677 if ( VolName == m_Volumes[i] ){
00678 PASSEDFILTER = true;
00679 break;
00680 }
00681 }
00682 }
00683 }
00684
00685 void ProcessCountingSvc::CheckTrack(const G4Step* aStep){
00686 G4Track* aTrack = aStep->GetTrack() ;
00687 G4int trackID= aTrack->GetTrackID();
00688
00689 int i = 0;
00690 for ( ; i < m_TrackIDs.size(); ){
00691
00692 if ( trackID == m_TrackIDs[i] ){
00693
00694 break;
00695 }
00696 i++;
00697 }
00698
00699 if ( i == m_TrackIDs.size() ){
00700
00701 m_TrackIDs.push_back(trackID);
00702 }
00703 }
