00001
00002
00003
00004
00005
00006
00007
00008 #include "McTruthSvc.hh"
00009
00010 #include <string>
00011 #include <iostream>
00012 #include <fstream>
00013 #include <sstream>
00014
00015 #include "myglobals.hh"
00016 #include "G4Event.hh"
00017 #include "G4Track.hh"
00018 #include "G4VProcess.hh"
00019
00020 #include "MyString2Anything.hh"
00021
00022 #include "MyRoot.hh"
00023
00024 McTruthSvc* McTruthSvc::fMcTruthSvc = 0;
00025
00026 McTruthSvc::McTruthSvc()
00027 {
00028 if (fMcTruthSvc){
00029 G4Exception("McTruthSvc::McTruthSvc()","Run0031",
00030 FatalException, "McTruthSvc constructed twice.");
00031 }
00032 fMcTruthSvc = this;
00033 ReSet();
00034 }
00035
00036 McTruthSvc::~McTruthSvc()
00037 {
00038 printf("~McTruthSvc\n");
00039 }
00040
00041 McTruthSvc* McTruthSvc::GetMcTruthSvc(){
00042 if ( !fMcTruthSvc ){
00043 fMcTruthSvc = new McTruthSvc;
00044 }
00045 return fMcTruthSvc;
00046 }
00047
00048 void McTruthSvc::Initialize(){
00049 m_nTracksAll = 0;
00050 m_nTracks = 0;
00051 m_dictpid.clear();
00052 m_dicttid.clear();
00053 m_dicttime.clear();
00054 m_pid.clear();
00055 m_tid.clear();
00056 m_ptid.clear();
00057 m_ppid.clear();
00058 m_time.clear();
00059 m_px.clear();
00060 m_py.clear();
00061 m_pz.clear();
00062 m_e.clear();
00063 m_ekin.clear();
00064 m_x.clear();
00065 m_y.clear();
00066 m_z.clear();
00067 m_particleName.clear();
00068 m_charge.clear();
00069 m_process.clear();
00070 m_volume.clear();
00071 }
00072
00073 void McTruthSvc::SetBranch(){
00074 if (!m_Switch) return;
00075 MyRoot* myRoot = MyRoot::GetMyRoot();
00076 if( flag_nTracks ) myRoot->SetBranch("McTruth_nTracks", &m_nTracks);
00077 if( flag_pid ) myRoot->SetBranch("McTruth_pid", &m_pid);
00078 if( flag_tid ) myRoot->SetBranch("McTruth_tid", &m_tid);
00079 if( flag_ptid ) myRoot->SetBranch("McTruth_ptid", &m_ptid);
00080 if( flag_ppid ) myRoot->SetBranch("McTruth_ppid", &m_ppid);
00081 if( flag_time ) myRoot->SetBranch("McTruth_time", &m_time);
00082 if( flag_px ) myRoot->SetBranch("McTruth_px", &m_px);
00083 if( flag_py ) myRoot->SetBranch("McTruth_py", &m_py);
00084 if( flag_pz ) myRoot->SetBranch("McTruth_pz", &m_pz);
00085 if( flag_e ) myRoot->SetBranch("McTruth_e", &m_e);
00086 if( flag_e ) myRoot->SetBranch("McTruth_ekin", &m_ekin);
00087 if( flag_x ) myRoot->SetBranch("McTruth_x", &m_x);
00088 if( flag_y ) myRoot->SetBranch("McTruth_y", &m_y);
00089 if( flag_z ) myRoot->SetBranch("McTruth_z", &m_z);
00090 if( flag_charge ) myRoot->SetBranch("McTruth_charge", &m_charge);
00091 if( flag_particleName ) myRoot->SetBranch("McTruth_particleName", &m_particleName);
00092 if( flag_process ) myRoot->SetBranch("McTruth_process", &m_process);
00093 if( flag_volume ) myRoot->SetBranch("McTruth_volume", &m_volume);
00094 }
00095
00096 void McTruthSvc::ReadOutputCard(G4String filename){
00097 ReSet();
00098 std::ifstream fin_card(filename);
00099 if(!fin_card){
00100 std::cout<<"In McTruthSvc::ReadOutputCard, cannot open "<<filename<<"!!!"<<std::endl;
00101 G4Exception("McTruthSvc::ReadOutputCard()",
00102 "InvalidSetup", FatalException,
00103 "cannot find output card");
00104 }
00105 std::stringstream buf_card;
00106 std::string s_card;
00107 int n_output_section_symbol = 0;
00108 int n_filter_section_symbol = 0;
00109 while(getline(fin_card,s_card)){
00110 buf_card.str("");
00111 buf_card.clear();
00112 buf_card<<s_card;
00113
00114
00115 const char* c_card = s_card.c_str();
00116 int length = strlen(c_card);
00117 int offset = 0;
00118 for ( ; offset < length; offset++ ){
00119 if ( c_card[offset] != ' ' ) break;
00120 }
00121 if ( c_card[offset] == '#' || (c_card[offset] == '/' && c_card[offset+1] == '/') || length - offset == 0 ){
00122 continue;
00123 }
00124
00125 std::string name;
00126 buf_card>>name;
00127 if ( n_output_section_symbol == 0 ){
00128 if ( name == "MCTRUTH_SECTION" ){
00129 n_output_section_symbol++;
00130 }
00131 }
00132 else if ( n_output_section_symbol == 1 ){
00133 if ( name == "MCTRUTH_SECTION" ){
00134 n_output_section_symbol++;
00135 }
00136 else if( name == "nTracks" ) flag_nTracks = true;
00137 else if( name == "pid" ) flag_pid = true;
00138 else if( name == "tid2pid" ) flag_tid2pid= true;
00139 else if( name == "tid2time" ) flag_tid2time= true;
00140 else if( name == "tid" ) flag_tid = true;
00141 else if( name == "ptid" ) flag_ptid = true;
00142 else if( name == "ppid" ) flag_ppid = true;
00143 else if( name == "time" ) {flag_time = true; buf_card>>unitName_time; unit_time = MyString2Anything::get_U(unitName_time);}
00144 else if( name == "px" ) {flag_px = true; buf_card>>unitName_px; unit_px = MyString2Anything::get_U(unitName_px);}
00145 else if( name == "py" ) {flag_py = true; buf_card>>unitName_py; unit_py = MyString2Anything::get_U(unitName_py);}
00146 else if( name == "pz" ) {flag_pz = true; buf_card>>unitName_pz; unit_pz = MyString2Anything::get_U(unitName_pz);}
00147 else if( name == "e" ) {flag_e = true; buf_card>>unitName_e; unit_e = MyString2Anything::get_U(unitName_e);}
00148 else if( name == "x" ) {flag_x = true; buf_card>>unitName_x; unit_x = MyString2Anything::get_U(unitName_x);}
00149 else if( name == "y" ) {flag_y = true; buf_card>>unitName_y; unit_y = MyString2Anything::get_U(unitName_y);}
00150 else if( name == "z" ) {flag_z = true; buf_card>>unitName_z; unit_z = MyString2Anything::get_U(unitName_z);}
00151 else if( name == "particleName" ) flag_particleName = true;
00152 else if( name == "charge" ) flag_charge = true;
00153 else if( name == "process" ) flag_process = true;
00154 else if( name == "volume" ) flag_volume = true;
00155 else{
00156 std::cout<<"In McTruthSvc::ReadOutputCard, unknown name: "<<name<<" in file "<<filename<<std::endl;
00157 std::cout<<"Will ignore this line!"<<std::endl;
00158 }
00159 }
00160
00161 if ( n_filter_section_symbol == 0 ){
00162 if ( name == "MCTRUTHFILTER_SECTION" ){
00163 n_filter_section_symbol++;
00164 }
00165 }
00166 else if ( n_filter_section_symbol == 1 ){
00167 if ( name == "MCTRUTHFILTER_SECTION" ){
00168 n_filter_section_symbol++;
00169 }
00170 else if( name == "Switch" ) m_Switch = true;
00171 else if( name == "nTracks" ) buf_card>>m_maxnTracks;
00172 else if( name == "WL" ){
00173 int pid = 0;
00174 buf_card>>pid;
00175 white_list.push_back(pid);
00176 }
00177 else if( name == "BL" ){
00178 int pid = 0;
00179 buf_card>>pid;
00180 black_list.push_back(pid);
00181 }
00182 else{
00183 G4double para;
00184 std::string unit;
00185 buf_card>>para>>unit;
00186 para *= MyString2Anything::get_U(unit);
00187 if( name == "minp" ) m_minp = para;
00188 else if( name == "mine" ) m_mine = para;
00189 else if( name == "mint" ) m_mint = para;
00190 else if( name == "maxt" ) m_maxt = para;
00191 else{
00192 std::cout<<"In McTruthSvc::ReadOutputCard, unknown name: "<<name<<" in file "<<filename<<std::endl;
00193 std::cout<<"Will ignore this line!"<<std::endl;
00194 }
00195 }
00196 }
00197
00198 if ( n_output_section_symbol > 1 && n_filter_section_symbol > 1 ){
00199 break;
00200 }
00201 }
00202 buf_card.str("");
00203 buf_card.clear();
00204 if ( n_output_section_symbol <= 1 ){
00205 std::cout<<"*****************WARNING********************"<<std::endl;
00206 std::cout<<"In McTruthSvc::ReadOutputCard, failed to find enough section seperators for output in file "<<filename<<std::endl;
00207 std::cout<<"Will use default settings."<<std::endl;
00208 std::cout<<"********************************************"<<std::endl;
00209 }
00210 if ( n_filter_section_symbol<= 1 ){
00211 std::cout<<"*****************WARNING********************"<<std::endl;
00212 std::cout<<"In McTruthSvc::ReadOutputCard, failed to find enough section seperators for filter in file "<<filename<<std::endl;
00213 std::cout<<"Will use default settings."<<std::endl;
00214 std::cout<<"********************************************"<<std::endl;
00215 }
00216 fin_card.close();
00217 ShowOutCard();
00218 }
00219
00220 void McTruthSvc::ReSet(){
00221 flag_nTracks = false;
00222 flag_pid = false;
00223 flag_tid2pid = false;
00224 flag_tid2time = false;
00225 flag_tid = false;
00226 flag_ptid = false;
00227 flag_ppid = false;
00228 flag_time = false;
00229 flag_px = false;
00230 flag_py = false;
00231 flag_pz = false;
00232 flag_e = false;
00233 flag_x = false;
00234 flag_y = false;
00235 flag_z = false;
00236 flag_charge = false;
00237 flag_particleName = false;
00238 flag_process = false;
00239 flag_volume = false;
00240 m_Switch = false;
00241 m_minp = 0;
00242 m_mine = 0;
00243 m_maxnTracks = 0;
00244 m_mint = 0;
00245 m_maxt = 0;
00246 white_list.clear();
00247 black_list.clear();
00248 unitName_time="s";
00249 unitName_px ="GeV";
00250 unitName_py ="GeV";
00251 unitName_pz ="GeV";
00252 unitName_e ="GeV";
00253 unitName_x ="cm";
00254 unitName_y ="cm";
00255 unitName_z ="cm";
00256 unit_time=s;
00257 unit_px =GeV;
00258 unit_py =GeV;
00259 unit_pz =GeV;
00260 unit_e =GeV;
00261 unit_x =cm;
00262 unit_y =cm;
00263 unit_z =cm;
00264 }
00265
00266 void McTruthSvc::ShowOutCard(){
00267 std::cout<<"*************************Output settings for McTruthSvc***************************"<<std::endl;
00268 std::cout<<"output nTracks? "<<(flag_nTracks?" yes":" no")<<std::endl;
00269 std::cout<<"output pid? "<<(flag_pid?" yes":" no")<<std::endl;
00270 std::cout<<"enable tid2pid? "<<(flag_tid2pid?" yes":" no")<<std::endl;
00271 std::cout<<"enable tid2time? "<<(flag_tid2time?" yes":" no")<<std::endl;
00272 std::cout<<"output tid? "<<(flag_tid?" yes":" no")<<std::endl;
00273 std::cout<<"output ptid? "<<(flag_ptid?" yes":" no")<<std::endl;
00274 std::cout<<"output ppid? "<<(flag_ppid?" yes":" no")<<std::endl;
00275 std::cout<<"output time? "<<(flag_time?" yes":" no")<<", unit: "<<unitName_time<<std::endl;
00276 std::cout<<"output px? "<<(flag_px?" yes":" no")<<", unit: "<<unitName_px<<std::endl;
00277 std::cout<<"output py? "<<(flag_py?" yes":" no")<<", unit: "<<unitName_py<<std::endl;
00278 std::cout<<"output pz? "<<(flag_pz?" yes":" no")<<", unit: "<<unitName_pz<<std::endl;
00279 std::cout<<"output e? "<<(flag_e?" yes":" no")<<", unit: "<<unitName_e<<std::endl;
00280 std::cout<<"output x? "<<(flag_x?" yes":" no")<<", unit: "<<unitName_x<<std::endl;
00281 std::cout<<"output y? "<<(flag_y?" yes":" no")<<", unit: "<<unitName_y<<std::endl;
00282 std::cout<<"output z? "<<(flag_z?" yes":" no")<<", unit: "<<unitName_z<<std::endl;
00283 std::cout<<"output particleName? "<<(flag_particleName?" yes":" no")<<std::endl;
00284 std::cout<<"output charge? "<<(flag_charge?" yes":" no")<<std::endl;
00285 std::cout<<"output process? "<<(flag_process?" yes":" no")<<std::endl;
00286 std::cout<<"output volume? "<<(flag_volume?" yes":" no")<<std::endl;
00287 std::cout<<"Switch on? "<<(m_Switch?"yes":"no")<<std::endl;
00288 std::cout<<"minp = "<<m_minp/MeV<<"MeV"<<std::endl;
00289 std::cout<<"mine = "<<m_mine/MeV<<"MeV"<<std::endl;
00290 std::cout<<"maxnTracks = "<<m_maxnTracks<<std::endl;
00291 std::cout<<"mint = "<<m_mint/ns<<"ns"<<std::endl;
00292 std::cout<<"maxt = "<<m_maxt/ns<<"ns"<<std::endl;
00293 std::cout<<"white list: "<<std::endl;
00294 for ( int i = 0; i< white_list.size(); i++){
00295 std::cout <<" Only tracks with these following PDGCodes will be recorded:"<<std::endl;
00296 std::cout<<" "<<i<<": "<<white_list[i]<<std::endl;
00297 }
00298 if ( white_list.size() == 0 ){
00299 std::cout <<" Empty! So all tracks will be recorded!"<<std::endl;
00300 }
00301 std::cout<<"black list: "<<std::endl;
00302 for ( int i = 0; i< black_list.size(); i++){
00303 std::cout <<" Tracks with these following PDGCodes will NOT be recorded:"<<std::endl;
00304 std::cout<<" "<<i<<": "<<black_list[i]<<std::endl;
00305 }
00306 std::cout<<"******************************************************************************"<<std::endl;
00307 }
00308
00309 void McTruthSvc::SetValuePre(const G4Track* aTrack){
00310 G4int pid = aTrack->GetParticleDefinition()->GetPDGEncoding();
00311 G4double globalT=aTrack->GetGlobalTime();
00312 G4int trackID= aTrack->GetTrackID();
00313 if (flag_tid2pid||flag_tid2time){
00314 m_dicttid.push_back(trackID);
00315 if (flag_tid2pid){
00316 m_dictpid.push_back(pid);
00317 }
00318 if (flag_tid2time){
00319 m_dicttime.push_back(globalT);
00320 }
00321 }
00322 m_nTracksAll++;
00323
00324 if (!m_Switch) return;
00325
00326 if (m_maxnTracks&&m_nTracks>=m_maxnTracks) return;
00327
00328 if (m_minp&&aTrack->GetMomentum().mag()<m_minp) return;
00329
00330 if (m_mine&&aTrack->GetTotalEnergy()<m_mine) return;
00331
00332 if (m_mint&&globalT<m_mint) return;
00333 if (m_maxt&&globalT>m_maxt) return;
00334
00335 bool foundit = false;
00336 for (int i = 0; i<white_list.size(); i++){
00337 if (pid == white_list[i]) foundit=true;
00338 }
00339 if (white_list.size()==1&&white_list[0]==0){
00340 if (pid<1e7) foundit = true;
00341 }
00342 if (!foundit&&white_list.size()) return;
00343
00344 foundit = false;
00345 for (int i = 0; i<black_list.size(); i++){
00346 if (pid == black_list[i]) foundit=true;
00347 }
00348 if (foundit) return;
00349
00350 G4String processName;
00351 const G4VProcess* process = aTrack->GetCreatorProcess();
00352 if (process) {
00353 processName = process->GetProcessName();
00354 }
00355 else{
00356 processName = "NULL";
00357 }
00358 G4int ptid = aTrack->GetParentID();
00359 G4ThreeVector mom_3vec = aTrack->GetMomentum();
00360 G4ThreeVector pos_3vec = aTrack->GetVertexPosition();
00361 std::string volume = aTrack->GetLogicalVolumeAtVertex()?aTrack->GetLogicalVolumeAtVertex()->GetName():"NULL";
00362 std::string particleName = aTrack->GetParticleDefinition()->GetParticleName();
00363 int charge = aTrack->GetParticleDefinition()->GetPDGCharge();
00364 G4double energy = aTrack->GetTotalEnergy();
00365 G4double kinenergy = aTrack->GetKineticEnergy();
00366
00367 m_nTracks++;
00368 if(flag_pid) m_pid.push_back(pid);
00369 if(flag_tid) m_tid.push_back(trackID);
00370 if(flag_ptid) m_ptid.push_back(ptid);
00371 if(flag_ppid){
00372 int ptid = aTrack->GetParentID();
00373 int ppid = McTruthSvc::GetMcTruthSvc()->tid2pid(ptid);
00374 m_ppid.push_back(ppid);
00375 }
00376 if(flag_time) m_time.push_back(globalT/unit_time);
00377 if(flag_px) m_px.push_back(mom_3vec.x()/unit_px);
00378
00379 if(flag_py) m_py.push_back(mom_3vec.y()/unit_py);
00380 if(flag_pz) m_pz.push_back(mom_3vec.z()/unit_pz);
00381 if(flag_e) m_e.push_back(energy/unit_e);
00382 if(flag_e) m_ekin.push_back(kinenergy/unit_e);
00383 if(flag_x) m_x.push_back(pos_3vec.x()/unit_x);
00384 if(flag_y) m_y.push_back(pos_3vec.y()/unit_y);
00385 if(flag_z) m_z.push_back(pos_3vec.z()/unit_z);
00386 if(flag_charge) m_charge.push_back(charge);
00387 if(flag_particleName) m_particleName.push_back(particleName);
00388 if(flag_process) m_process.push_back(processName);
00389 if(flag_volume) m_volume.push_back(volume);
00390 }
00391
00392 void McTruthSvc::SetValuePost(const G4Track* aTrack){
00393
00394
00395 }
