00001
00002
00003
00004
00005
00006
00007
00008 #include "MyAnalysisSvc.hh"
00009 #include "MyAnalysisSvcMessenger.hh"
00010
00011 #include "myglobals.hh"
00012 #include "G4Step.hh"
00013 #include "G4Event.hh"
00014 #include "G4Run.hh"
00015
00016 #include "MyRoot.hh"
00017 #include "LogSvc.hh"
00018 #include "MyDetectorManager.hh"
00019 #include "EventHeaderSvc.hh"
00020 #include "McTruthSvc.hh"
00021 #include "ProcessCountingSvc.hh"
00022 #include "MyTriggerSvc.hh"
00023 #include "PrimaryGeneratorAction.hh"
00024 #include "Randomize.hh"
00025
00026 #include "MyProcessManager.hh"
00027 #include "DEBUG.hh"
00028
00029 MyAnalysisSvc* MyAnalysisSvc::fMyAnalysisSvc = 0;
00030
00031 MyAnalysisSvc::MyAnalysisSvc()
00032 {
00033 if (fMyAnalysisSvc){
00034 G4Exception("MyAnalysisSvc::MyAnalysisSvc()","Run0031",
00035 FatalException, "MyAnalysisSvc constructed twice.");
00036 }
00037 fMyAnalysisSvc = this;
00038
00039 pMyAnalysisSvcMessenger = new MyAnalysisSvcMessenger(this);
00040 pMyRoot = MyRoot::GetMyRoot();
00041 pMyDetectorManager = MyDetectorManager::GetMyDetectorManager();
00042 pEventHeaderSvc = EventHeaderSvc::GetEventHeaderSvc();
00043 pMcTruthSvc = McTruthSvc::GetMcTruthSvc();
00044 pProcessCountingSvc = ProcessCountingSvc::GetProcessCountingSvc();
00045 pMyTriggerSvc = MyTriggerSvc::GetMyTriggerSvc();
00046 pPrimaryGeneratorAction = PrimaryGeneratorAction::GetPrimaryGeneratorAction();
00047 pMyProcessManager = MyProcessManager::GetMyProcessManager();
00048
00049
00050 run_name = getenv("RUNNAMEROOT");
00051
00052 ofile_name = getenv("OFILENAMEROOT");
00053
00054 std::string out_card = getenv("OUTCARDROOT");
00055 set_out_card(out_card.c_str());
00056
00057 m_minT = -1;
00058 m_maxT = -1;
00059 fPrintModulo = 100;
00060
00061 }
00062
00063 MyAnalysisSvc::~MyAnalysisSvc()
00064 {
00065 delete pMyAnalysisSvcMessenger;
00066 }
00067
00068 MyAnalysisSvc* MyAnalysisSvc::GetMyAnalysisSvc(){
00069 if ( !fMyAnalysisSvc ){
00070 fMyAnalysisSvc = new MyAnalysisSvc;
00071 }
00072 return fMyAnalysisSvc;
00073 }
00074
00075 void MyAnalysisSvc::set_out_card(G4String file_name){
00076 if(file_name[0] != '/'){
00077 G4String dir_name = getenv("CONFIGUREROOT");
00078 if (dir_name[dir_name.size()-1] != '/') dir_name.append("/");
00079 file_name = dir_name + file_name;
00080 }
00081
00082 pEventHeaderSvc->ReadOutputCard(file_name);
00083 pMyDetectorManager->ReadOutputCard(file_name);
00084 pMcTruthSvc->ReadOutputCard(file_name);
00085 pProcessCountingSvc->ReadOutputCard(file_name);
00086 pMyTriggerSvc->SetMyTrigger(file_name);
00087
00088 ReadOutputCard(file_name);
00089
00090 pMyRoot->SetVerbose(fVerbose);
00091 pMyRoot->SetPrintModulo(fPrintModulo);
00092 }
00093
00094 void MyAnalysisSvc::BeginOfRunAction(){
00095
00096 std::string log_name = getenv("LOGFILEROOT");
00097 if ( std::strlen(log_name.c_str()) != 0 ){
00098 LogSvc::GetLogSvc()->SetLogFile( log_name.c_str() );
00099 }
00100 run_num = LogSvc::GetLogSvc()->AddLog( run_name.c_str() );
00101 std::cout<<"logfile finished for "<<run_name<<", "<<run_num<<std::endl;
00102
00103
00104 pMyRoot->OpenFile(ofile_name);
00105 pMyRoot->CreateTree(tree_name, fCircular);
00106
00107
00108 pMyDetectorManager->SetBranch();
00109 pEventHeaderSvc->SetBranch();
00110 pMcTruthSvc->SetBranch();
00111 pProcessCountingSvc->SetBranch();
00112
00113 t_begin = (double)clock();
00114 std::cout<<"evt_num memSize time R0 R1"<<std::endl;
00115 std::cout<<" MB sec "<<std::endl;
00116 }
00117
00118 void MyAnalysisSvc::PreUserTrackingAction(const G4Track* aTrack){
00119
00120 pMcTruthSvc->SetValuePre(aTrack);
00121 }
00122
00123 void MyAnalysisSvc::PostUserTrackingAction(const G4Track* aTrack) {
00124
00125 pMcTruthSvc->SetValuePost(aTrack);
00126 }
00127
00128 void MyAnalysisSvc::EndOfRunAction(const G4Run* aRun){
00129 int NbOfEvents = aRun->GetNumberOfEvent();
00130 t_end = (double)clock();
00131 std::cout<<"\n##############################################"<<std::endl;
00132 std::cout<<"TOTAL TIME COST IS: "<<(double)(t_end-t_begin)/CLOCKS_PER_SEC*1000<<"ms"<<std::endl;
00133 std::cout<<"TIME COST PER EVENT: "<<(double)(t_end-t_begin)/CLOCKS_PER_SEC/NbOfEvents*1000<<"ms"<<std::endl;
00134 std::cout<<"##############################################\n"<<std::endl;
00135 pMyRoot->Write();
00136 pMyRoot->Close();
00137 }
00138
00139 void MyAnalysisSvc::BeginOfEventAction(){
00140 event_start_time = (double) clock();
00141
00142
00143 pMcTruthSvc->Initialize();
00144 pProcessCountingSvc->Initialize();
00145 }
00146
00147 void MyAnalysisSvc::EndOfEventAction(const G4Event* evt){
00148 int evt_num = evt->GetEventID();
00149
00150 pMyDetectorManager->Digitize();
00151
00152 double weight = 1;
00153 void *result = pPrimaryGeneratorAction->get_extra("weight");
00154 if (result) weight = *((double*)result);
00155 pEventHeaderSvc->SetValue(evt,run_num,weight);
00156
00157
00158 if ( pMyTriggerSvc->TriggerIt(evt) ) pMyRoot->Fill();
00159
00160
00161 if (fAutoSave){
00162 if ( evt_num%fAutoSave == 0 ){
00163 pMyRoot->Save();
00164 int nBytes = pMyRoot->FlushBaskets();
00165 std::cout<<"Event "<<evt_num<<", "<<nBytes<<"Bytes data written"<<std::endl;
00166 }
00167 }
00168
00169
00170 if (evt_num%fPrintModulo == 0) {
00171 std::cout <<evt_num
00172 <<" "<<pMyProcessManager->GetMemorySize()
00173 <<" "<<(((double)clock()-t_begin)/CLOCKS_PER_SEC)
00174 <<" "<<pEventHeaderSvc->get_R0()
00175 <<" "<<pEventHeaderSvc->get_R1()
00176 <<std::endl;
00177 }
00178 }
00179
00180 void MyAnalysisSvc::SteppingAction(const G4Step* aStep){
00181
00182
00183 double current_time = (double) clock();
00184
00185 pProcessCountingSvc->SetValue(aStep);
00186 G4Track* aTrack = aStep->GetTrack() ;
00187 G4int nSteps = aTrack->GetCurrentStepNumber();
00188 bool needStopAndKill = false;
00189 if (nSteps>2e5){
00190 std::cout<<"### This track is killed for that nSteps>2e5 ###"<<std::endl;
00191 needStopAndKill = true;
00192 }
00193 G4double globalT=aTrack->GetGlobalTime();
00194 if ((m_minT>=0&&m_minT>globalT)&&(m_maxT>=0&&m_maxT<globalT)){
00195 std::cout<<"### This track is killed for that globalT < "<<m_minT<<" || globalT > "<<m_maxT<<" ###"<<std::endl;
00196 needStopAndKill = true;
00197 }
00198 if (( current_time - event_start_time > 5*CLOCKS_PER_SEC)){
00199 std::cout<<"### This track is killed for that elapsed time in this event is larger than 5 sec ###"<<std::endl;
00200 needStopAndKill = true;
00201 }
00202 if (needStopAndKill){
00203 aTrack->SetTrackStatus(fStopAndKill);
00204 std::cout<<" => ["<<aTrack->GetTrackID()<<"]: "<<aTrack->GetParticleDefinition()->GetParticleName()<<", "<<G4BestUnit(aStep->GetPreStepPoint()->GetKineticEnergy(),"Energy")<<", @ \""<<aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetName()<<"\""<<std::endl;
00205 }
00206 }
00207
00208 void MyAnalysisSvc::InitialStepAction(const G4Step* aStep){
00209 pProcessCountingSvc->InitialStep(aStep);
00210 }
00211
00212 void MyAnalysisSvc::ASDI(G4String name){
00213 pProcessCountingSvc->AddASDI(name);
00214 }
00215
00216 void MyAnalysisSvc::PSDI(G4String name){
00217 pProcessCountingSvc->AddPSDI(name);
00218 }
00219
00220 void MyAnalysisSvc::ReadOutputCard(G4String file_name){
00221 std::ifstream fin_card(file_name);
00222 if(!fin_card){
00223 std::cout<<"In MyAnalysisSvc::ReadOutputCard, cannot open "<<file_name<<"!!!"<<std::endl;
00224 G4Exception("MyAnalysisSvc::ReadOutputCard()",
00225 "InvalidSetup", FatalException,
00226 "cannot find output card");
00227 }
00228 std::stringstream buf_card;
00229 std::string s_card;
00230 bool find_tree_name = false;
00231 bool find_AutoSave = false;
00232 bool find_Circular = false;
00233 bool find_Verbose = false;
00234 bool find_PrintModulo = false;
00235 while(getline(fin_card,s_card)){
00236 buf_card.str("");
00237 buf_card.clear();
00238 buf_card<<s_card;
00239
00240
00241 const char* c_card = s_card.c_str();
00242 int length = strlen(c_card);
00243 int offset = 0;
00244 for ( ; offset < length; offset++ ){
00245 if ( c_card[offset] != ' ' ) break;
00246 }
00247 if ( c_card[offset] == '#' || (c_card[offset] == '/' && c_card[offset+1] == '/') || length - offset == 0 ){
00248 continue;
00249 }
00250 std::string name;
00251 buf_card>>name;
00252 if ( name == "tree_name" ){
00253 buf_card>>tree_name;
00254 std::cout<<"In MyAnalysisSvc::ReadOutputCard, tree_name will be set to "<<tree_name<<std::endl;
00255 find_tree_name = true;
00256 }
00257 else if ( name == "AutoSave" ){
00258 buf_card>>fAutoSave;
00259 std::cout<<"In MyAnalysisSvc::ReadOutputCard, fAutoSave will be set to "<<fAutoSave<<std::endl;
00260 find_AutoSave = true;
00261 }
00262 else if ( name == "Circular" ){
00263 buf_card>>fCircular;
00264 std::cout<<"In MyAnalysisSvc::ReadOutputCard, fCircular will be set to "<<fCircular<<std::endl;
00265 find_Circular = true;
00266 }
00267 else if ( name == "Verbose" ){
00268 buf_card>>fVerbose;
00269 std::cout<<"In MyAnalysisSvc::ReadOutputCard, fVerbose will be set to "<<fVerbose<<std::endl;
00270 find_Verbose = true;
00271 }
00272 else if ( name == "PrintModulo" ){
00273 buf_card>>fPrintModulo;
00274 std::cout<<"In MyAnalysisSvc::ReadOutputCard, fPrintModulo will be set to "<<fPrintModulo<<std::endl;
00275 find_PrintModulo = true;
00276 }
00277 if ( find_Verbose && find_PrintModulo && find_Circular && find_AutoSave && find_tree_name ) break;
00278 }
00279 if (!find_tree_name){
00280 std::cout<<"In MyAnalysisSvc::ReadOutputCard, tree_name not found in card, will be set to t as default"<<std::endl;
00281 tree_name = "t";
00282 }
00283 if (!find_AutoSave){
00284 std::cout<<"In MyAnalysisSvc::ReadOutputCard, AutoSave not found in card, will be set to 0 as default"<<std::endl;
00285 fAutoSave = 0;
00286 }
00287 if (!find_Circular){
00288 std::cout<<"In MyAnalysisSvc::ReadOutputCard, Circular not found in card, will be set to 0 as default"<<std::endl;
00289 fCircular = 0;
00290 }
00291 if (!find_Verbose){
00292 std::cout<<"In MyAnalysisSvc::ReadOutputCard, Verbose not found in card, will be set to 0 as default"<<std::endl;
00293 fVerbose = 0;
00294 }
00295 if (!find_PrintModulo){
00296 std::cout<<"In MyAnalysisSvc::ReadOutputCard, PrintModulo not found in card, will be set to 0 as default"<<std::endl;
00297 fPrintModulo = 100;
00298 }
00299 fin_card.close();
00300 buf_card.str("");
00301 buf_card.clear();
00302 }
