#include #include #include //ROOT headers #include #include #include #include #include #include //#include //G4VMC headers #include //VMC headers #include "PamVMCApplication.h" #include "PamVMCOptMgr.h" //This is main VMC routine using namespace std; #define DEBUG 0 void PackDump(TString path, TObjString* os); void UnPackDump(TString path, TObjString* os); int main(int argc, char **argv) { if(argc < 2) { cout<<"ERROR: PamVMCmain - no input file specified"<Load("libg4root.so"); cout<<" Input card file: "<ParseFile(argv[1]); omg->Print(); //Checking if cookie and temp_path TString working_path = omg->GetOutPath(); TString cookie = ""; if( omg->IsUseCookie() ){ cookie = gSystem->Getenv("PBS_JOBCOOKIE"); working_path = omg->GetTmpPath(); working_path += "/"; working_path += cookie; } TSystemFile *ff = new TSystemFile(); if( !ff->IsDirectory(working_path)) gSystem->mkdir(working_path); delete ff; //if I'm going to store random stuff I need temporary directory to dump g4 engine status TString tmpseeddir; Int_t save_ev; TObjString* ev_dump = NULL; TFile * dmp_file = NULL; TTree * t_dump = NULL; tmpseeddir = working_path; tmpseeddir += "/"; tmpseeddir += omg->GetOutPattern(); tmpseeddir+="_g4rndm_tmp"; if (omg->DoWriteRandom()){ cout<<"Making Geant-4 RNDM tmp directory: "<mkdir(tmpseeddir); //Changing directory gSystem->Exec("cd "+working_path); gSystem->cd(working_path); //gDirectory->cd(working_path); //Making files and trees for saved random statuses ev_dump = new TObjString(); TString file_dump_name_tmp = working_path; file_dump_name_tmp += "/"; file_dump_name_tmp += omg->GetRandomOutFilename(); dmp_file = new TFile(file_dump_name_tmp,"recreate"); t_dump = new TTree("dumps","G4 rndm for events"); t_dump->Branch("evno", &save_ev,"save_ev/I"); t_dump->Branch("dump_cont","TObjString", &ev_dump); } //INITIALIZATION OF MC TString out_tmp_fname = working_path; out_tmp_fname += "/"; out_tmp_fname += omg->GetOutPattern(); cout<SetVerboseLevel( omg->GetVerboseLevel() ); //initialize physycs with specific RunConfiguration, defined in macro appl->InitMC( omg->GetG4ConfigC() ); cout<<"MACRO READED"<ProcessGeantMacro( omg->GetG4ConfigIn() ); //set maximum allowed steps ((TGeant4*)gMC)->SetMaxNStep( omg->GetMaxSteps() ); //generate initial random seed if any TFile * read_rndm_file = NULL; TTree * read_rndm_tree = NULL; Int_t inpev; TObjString * rnd_dump = NULL; TString tmpfile; if ( omg->IsReadRandomMode() ){ //read random seed from file cout<<"Making Geant-4 RNDM tmp directory: "<mkdir(tmpseeddir); tmpfile = tmpseeddir+"/currEvent.rndm"; gSystem->Exec("cd "+working_path); gSystem->cd(working_path); read_rndm_file = new TFile( omg->GetRandomInFilename() ); read_rndm_tree = (TTree*)read_rndm_file->Get("dumps"); Int_t inpev; read_rndm_tree->SetBranchAddress("evno", &inpev); read_rndm_tree->SetBranchAddress("dump_cont", &rnd_dump); } else { Int_t seed1, seed2; if( omg->IsUserRandomMode() ){ omg->GetUserRandomSeed( seed1, seed2 ); } else { TRandom3 rinit1(0); TRandom3 rinit2(0); seed1 = rinit1.Integer(numeric_limits::max()); seed2 = rinit2.Integer(numeric_limits::max()); } cout<<"Initial RND seeds for G4 core are: " << seed1<<" , " << seed2<ProcessGeantCommand(command); } /*seting tmp directory to save status of random gen if any*/ if( omg->DoWriteRandom() ){ ((TGeant4*)gMC)->ProcessGeantCommand("/random/setDirectoryName "+tmpseeddir); } Double_t x,y,z,theta,phi,momentum; Int_t id,pdg; Double_t xl, yl, zl, xh, yh, zh, tl, th, pl, ph; Int_t nev_max = omg->GetEvMax(); //* if I have to read primary kinematics from a file do it *// TFile* rfin = NULL; TTree* primtree = NULL; TClonesArray * prims = NULL; PamVMCPrimary *p = NULL; if ( omg->GetReadPrimMode() ){ rfin = new TFile ( omg->GetReadPrimPath(), "read"); primtree = (TTree*) rfin->Get("hits"); if(!primtree) cout<<"ERROR: tree 'hits' with primaries kinematics not found"<SetBranchAddress("PRIM", &prims); nev_max = primtree->GetEntries(); } cout<<"OK1"<GetPrimaryGenerator(); gpr->SetRandom(aux_rnd); //gObjectTable->Print(); gSystem->Exec(" ps aux | grep PamVMC.exe" ); while (!stop){ if ( omg->GetReadPrimMode() ){ //read kinematics from file rfin->cd(); primtree->GetEntry(proc_ev); p = (PamVMCPrimary*)prims->At(0); x = p->fX0; y = p->fY0; z = p->fZ0; theta = p->fTHETA; phi = p->fPHI; pdg = p->fPDG; momentum = p->fP0; id = p->fID; gpr->SetParticle(pdg); gpr->SetParticleID(id); gpr->SetDirection(theta,phi); gpr->SetMomentum(momentum); gpr->SetPosition(x,y,z); } else { //generate kinematics pdg = omg->GetPrimPDG(); id = proc_ev; gpr->SetParticle(pdg); gpr->SetParticleID(id); //generate spatial distribution if( omg->IsCustomKinematics() ){ //vertex if( omg->IsBoxVertex() ){ //box omg->GetBoxVertex(xl,xh,yl,yh,zl,zh); gpr->GenPosition(xl, xh, yl, yh, zl, zh); } else { //point omg->GetFixedPointVertex(x,y,z); gpr->SetPosition(x,y,z); } //angles if( omg->IsIsotropicPlane() ){ //isoplane angles omg->GetIsotropicAngles(tl, th, pl, ph); gpr->GenDirection(tl, th, pl, ph); } else { //fixed angles omg->GetFixedAngles(theta, phi); gpr->SetDirection(theta,phi); } } else { // GOD mode gpr->Gen2PiPosAng(); } gpr->GetPosition(x,y,z); gpr->GetDirection(theta, phi); //generate MOMENTUM Momentum_units mu = omg->GetPrimMomUnits(); if( omg->GetPrimMomMode() == FIXED ){ momentum = omg->GetMomFixed(); switch (mu){ case GeV_T: gpr->SetMomentum(gpr->KinEToMomentum(momentum)); break; case GV_p: gpr->SetMomentum(momentum); case GV_R: gpr->SetMomentum( gpr->RigToMomentum(momentum)); default: break; } } if( omg->GetPrimMomMode() == FLAT ){ Double_t m1, m2; omg->GetMomLimits(m1,m2); switch (mu){ case GeV_T: gpr->GenSpe( m1, m2, kTRUE); break; case GV_p: gpr->GenSpe( m1, m2, kFALSE); case GV_R: gpr->GenSpe( gpr->RigToMomentum(m1), gpr->RigToMomentum(m2), kFALSE); default: break; } } if( omg->GetPrimMomMode() == POWERLAW ){ Double_t m1, m2; omg->GetMomLimits(m1,m2); Double_t gamma = omg->GetMomGamma(); switch (mu){ case GeV_T: gpr->GenSpe( m1, m2, gamma, kTRUE); break; case GV_p: gpr->GenSpe( m1, m2, gamma, kFALSE); case GV_R: gpr->GenSpe( gpr->RigToMomentum(m1), gpr->RigToMomentum(m2), gamma, kFALSE); default: break; } } } /*if (DEBUG){ cout<<"EVENT: "<IsReadRandomMode() ){ read_rndm_file->cd(); UnPackDump(tmpfile, rnd_dump); TString rndfile = "/random/resetEngineFrom "+tmpfile; //cout<<"reading... "<ProcessGeantCommand(rndfile); } if( omg->DoWriteRandom() ){ //cout<<"saving RND"<ProcessGeantCommand("/random/setSavingFlag 1"); } gpr->ExtrapolateTrajectory(); //Added by V.Formato on 2012/06/19 if(!gpr->ToBeSimulated()){ proc_ev++; continue;} /////RUN/////// appl->RunMC(1); ///////////// proc_ev ++; accept_ev = gpr->GetNgood(); if (appl->GetTrigger()) trig_ev++; // if(DEBUG) cout<<"TOT_EV: "<Print(); //saving dumps if any if( omg->DoWriteRandom() ){ switch( omg->GetSaveCond() ){ case EVERYTHING: dmp_file->cd(); PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); save_ev = id; t_dump->Fill(); break; case TRIG_ONLY: if(appl->GetTrigger()){ dmp_file->cd(); PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); save_ev = id; t_dump->Fill(); } break; case ACCEPT_ONLY: if(gpr->Getgood()){ dmp_file->cd(); PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); save_ev = id; t_dump->Fill(); } break; default: break; } } //condition to stop switch( omg->GetStopMode() ){ case NEV_TOT: if( proc_ev >= nev_max) stop = kTRUE; break; case NEV_ACCEPT: if( accept_ev >= nev_max) stop = kTRUE; break; case NEV_TRIGG: if( trig_ev >= nev_max) stop = kTRUE; break; default: break; } if( (omg->GetReadPrimMode()) && ( proc_ev >= nev_max) ) stop = kTRUE; //gSystem->Exec(" ps aux | grep PamVMC.exe" ); } gSystem->Exec(" ps aux | grep PamVMC.exe" ); appl->FinishRun(); if( omg->DoWriteRandom() ){ dmp_file->cd(); t_dump->Write(); dmp_file->Close(); } cout<<"done...."<Close(); delete read_rndm_file; delete appl; //delete tmp directory with dumps gSystem->Exec("rm -rf "+tmpseeddir); //moving good files to good directory if( omg->IsUseCookie() ){ TString source_dir = working_path; TString dest_path = omg->GetOutPath(); TString out_root_fname = omg->GetOutPattern()+".root"; TString out_pam_fname = omg->GetOutPattern()+".pam"; TString out_dump_fname = omg->GetRandomOutFilename(); cout<<"--->Move "<Exec("mv " + source_dir+"/"+out_root_fname + " " + dest_path +"/"); cout<<"--->Move "<Exec("mv " + source_dir+"/"+out_pam_fname + " " + dest_path +"/"); if( omg->DoWriteRandom() ){ cout<<"--->Move "<Exec("mv " + out_dump_fname + " " + dest_path +"/"); } cout<<"->Cleaning up tmpdir: "<Exec("cd "+dest_path); gSystem->Exec("rm -rf " + source_dir ); } return 0; } void PackDump(TString path, TObjString* os){ ifstream ffile; ffile.open(path.Data(), ios::binary); if(!ffile) { cout<<"ERROR! wrong filename: "<SetString(buffer); delete[] buffer; } void UnPackDump(TString path, TObjString* os){ ofstream ffile; ffile.open(path.Data(), ios::binary); if(!ffile) { cout<<"ERROR! wrong filename: "<Sizeof(); // allocate memory: const char* buffer = 0; TString tmp = os->GetString(); buffer = tmp.Data(); // read data as a block: ffile.write (buffer,length-1); cout<<"->Unpack"<