| 1 |
#include <stdio.h> |
| 2 |
#include <stdlib.h> |
| 3 |
#include <iostream> |
| 4 |
|
| 5 |
//ROOT headers |
| 6 |
#include <TSystem.h> |
| 7 |
#include <TSystemFile.h> |
| 8 |
#include <TDirectory.h> |
| 9 |
#include <TFile.h> |
| 10 |
#include <TTree.h> |
| 11 |
#include <TObjString.h> |
| 12 |
//#include <TObjectTable.h> |
| 13 |
|
| 14 |
//G4VMC headers |
| 15 |
#include <TGeant4.h> |
| 16 |
|
| 17 |
//VMC headers |
| 18 |
#include "PamVMCApplication.h" |
| 19 |
#include "PamVMCOptMgr.h" |
| 20 |
//This is main VMC routine |
| 21 |
using namespace std; |
| 22 |
|
| 23 |
#define DEBUG 0 |
| 24 |
|
| 25 |
void PackDump(TString path, TObjString* os); |
| 26 |
void UnPackDump(TString path, TObjString* os); |
| 27 |
|
| 28 |
int main(int argc, char **argv) |
| 29 |
{ |
| 30 |
|
| 31 |
if(argc < 2) { |
| 32 |
cout<<"ERROR: PamVMCmain - no input file specified"<<endl; |
| 33 |
abort(); |
| 34 |
} |
| 35 |
//Option1 montecarlo selection, default G4 |
| 36 |
gSystem->Load("libg4root.so"); |
| 37 |
cout<<" Input card file: "<<argv[1]<<endl; |
| 38 |
|
| 39 |
PamVMCOptMgr * omg = PamVMCOptMgr::Instance(); |
| 40 |
omg->ParseFile(argv[1]); |
| 41 |
omg->Print(); |
| 42 |
|
| 43 |
//Checking if cookie and temp_path |
| 44 |
TString working_path = omg->GetOutPath(); |
| 45 |
TString cookie = ""; |
| 46 |
if( omg->IsUseCookie() ){ |
| 47 |
cookie = gSystem->Getenv("PBS_JOBCOOKIE"); |
| 48 |
working_path = omg->GetTmpPath(); |
| 49 |
working_path += "/"; |
| 50 |
working_path += cookie; |
| 51 |
} |
| 52 |
TSystemFile *ff = new TSystemFile(); |
| 53 |
if( !ff->IsDirectory(working_path)) gSystem->mkdir(working_path); |
| 54 |
delete ff; |
| 55 |
|
| 56 |
//if I'm going to store random stuff I need temporary directory to dump g4 engine status |
| 57 |
TString tmpseeddir; |
| 58 |
Int_t save_ev; |
| 59 |
TObjString* ev_dump = NULL; |
| 60 |
TFile * dmp_file = NULL; |
| 61 |
TTree * t_dump = NULL; |
| 62 |
tmpseeddir = working_path; |
| 63 |
tmpseeddir += "/"; |
| 64 |
tmpseeddir += omg->GetOutPattern(); |
| 65 |
tmpseeddir+="_g4rndm_tmp"; |
| 66 |
if (omg->DoWriteRandom()){ |
| 67 |
|
| 68 |
cout<<"Making Geant-4 RNDM tmp directory: "<<tmpseeddir<<endl; |
| 69 |
gSystem->mkdir(tmpseeddir); |
| 70 |
|
| 71 |
//Changing directory |
| 72 |
|
| 73 |
gSystem->Exec("cd "+working_path); |
| 74 |
gSystem->cd(working_path); |
| 75 |
//gDirectory->cd(working_path); |
| 76 |
|
| 77 |
//Making files and trees for saved random statuses |
| 78 |
|
| 79 |
ev_dump = new TObjString(); |
| 80 |
TString file_dump_name_tmp = working_path; |
| 81 |
file_dump_name_tmp += "/"; |
| 82 |
file_dump_name_tmp += omg->GetRandomOutFilename(); |
| 83 |
dmp_file = new TFile(file_dump_name_tmp,"recreate"); |
| 84 |
t_dump = new TTree("dumps","G4 rndm for events"); |
| 85 |
t_dump->Branch("evno", &save_ev,"save_ev/I"); |
| 86 |
t_dump->Branch("dump_cont","TObjString", &ev_dump); |
| 87 |
} |
| 88 |
|
| 89 |
//INITIALIZATION OF MC |
| 90 |
TString out_tmp_fname = working_path; |
| 91 |
out_tmp_fname += "/"; |
| 92 |
out_tmp_fname += omg->GetOutPattern(); |
| 93 |
cout<<out_tmp_fname<<endl; |
| 94 |
PamVMCApplication* appl |
| 95 |
= new PamVMCApplication("PAMG4_VMC", "PAMELA GEANT4 VMC application",out_tmp_fname,0); |
| 96 |
//setting verbose level |
| 97 |
appl->SetVerboseLevel( omg->GetVerboseLevel() ); |
| 98 |
//initialize physycs with specific RunConfiguration, defined in macro |
| 99 |
appl->InitMC( omg->GetG4ConfigC() ); |
| 100 |
cout<<"MACRO READED"<<endl; |
| 101 |
//additional flags and controls related to G4 core |
| 102 |
((TGeant4*)gMC)->ProcessGeantMacro( omg->GetG4ConfigIn() ); |
| 103 |
//set maximum allowed steps |
| 104 |
((TGeant4*)gMC)->SetMaxNStep( omg->GetMaxSteps() ); |
| 105 |
|
| 106 |
//generate initial random seed if any |
| 107 |
TFile * read_rndm_file = NULL; |
| 108 |
TTree * read_rndm_tree = NULL; |
| 109 |
Int_t inpev; |
| 110 |
TObjString * rnd_dump = NULL; |
| 111 |
TString tmpfile; |
| 112 |
if ( omg->IsReadRandomMode() ){ |
| 113 |
//read random seed from file |
| 114 |
cout<<"Making Geant-4 RNDM tmp directory: "<<tmpseeddir<<endl; |
| 115 |
gSystem->mkdir(tmpseeddir); |
| 116 |
tmpfile = tmpseeddir+"/currEvent.rndm"; |
| 117 |
gSystem->Exec("cd "+working_path); |
| 118 |
gSystem->cd(working_path); |
| 119 |
read_rndm_file = new TFile( omg->GetRandomInFilename() ); |
| 120 |
read_rndm_tree = (TTree*)read_rndm_file->Get("dumps"); |
| 121 |
Int_t inpev; |
| 122 |
read_rndm_tree->SetBranchAddress("evno", &inpev); |
| 123 |
read_rndm_tree->SetBranchAddress("dump_cont", &rnd_dump); |
| 124 |
} else { |
| 125 |
Int_t seed1, seed2; |
| 126 |
if( omg->IsUserRandomMode() ){ |
| 127 |
omg->GetUserRandomSeed( seed1, seed2 ); |
| 128 |
} else { |
| 129 |
TRandom3 rinit1(0); |
| 130 |
TRandom3 rinit2(0); |
| 131 |
seed1 = rinit1.Integer(numeric_limits<int>::max()); |
| 132 |
seed2 = rinit2.Integer(numeric_limits<int>::max()); |
| 133 |
} |
| 134 |
cout<<"Initial RND seeds for G4 core are: " << seed1<<" , " << seed2<<endl; |
| 135 |
TString command = "/random/setSeeds "; |
| 136 |
command += seed1; |
| 137 |
command += " "; |
| 138 |
command += seed2; |
| 139 |
((TGeant4*)gMC)->ProcessGeantCommand(command); |
| 140 |
} |
| 141 |
|
| 142 |
/*seting tmp directory to save status of random gen if any*/ |
| 143 |
if( omg->DoWriteRandom() ){ |
| 144 |
((TGeant4*)gMC)->ProcessGeantCommand("/random/setDirectoryName "+tmpseeddir); |
| 145 |
} |
| 146 |
|
| 147 |
|
| 148 |
Double_t x,y,z,theta,phi,momentum; |
| 149 |
Int_t id,pdg; |
| 150 |
Double_t xl, yl, zl, xh, yh, zh, tl, th, pl, ph; |
| 151 |
Int_t nev_max = omg->GetEvMax(); |
| 152 |
|
| 153 |
|
| 154 |
//* if I have to read primary kinematics from a file do it *// |
| 155 |
TFile* rfin = NULL; |
| 156 |
TTree* primtree = NULL; |
| 157 |
TClonesArray * prims = NULL; |
| 158 |
PamVMCPrimary *p = NULL; |
| 159 |
if ( omg->GetReadPrimMode() ){ |
| 160 |
rfin = new TFile ( omg->GetReadPrimPath(), "read"); |
| 161 |
primtree = (TTree*) rfin->Get("hits"); |
| 162 |
if(!primtree) cout<<"ERROR: tree 'hits' with primaries kinematics not found"<<endl; |
| 163 |
primtree->SetBranchAddress("PRIM", &prims); |
| 164 |
nev_max = primtree->GetEntries(); |
| 165 |
} |
| 166 |
|
| 167 |
cout<<"OK1"<<endl; |
| 168 |
|
| 169 |
//MAIN LOOP |
| 170 |
Int_t proc_ev = 0; |
| 171 |
Int_t accept_ev = 0; |
| 172 |
Int_t trig_ev = 0; |
| 173 |
Bool_t stop = kFALSE; |
| 174 |
TRandom3 * aux_rnd = new TRandom3(0); |
| 175 |
PamVMCPrimaryGenerator* gpr = appl->GetPrimaryGenerator(); |
| 176 |
gpr->SetRandom(aux_rnd); |
| 177 |
|
| 178 |
//gObjectTable->Print(); |
| 179 |
gSystem->Exec(" ps aux | grep PamVMC.exe" ); |
| 180 |
while (!stop){ |
| 181 |
|
| 182 |
|
| 183 |
if ( omg->GetReadPrimMode() ){ |
| 184 |
//read kinematics from file |
| 185 |
rfin->cd(); |
| 186 |
primtree->GetEntry(proc_ev); |
| 187 |
p = (PamVMCPrimary*)prims->At(0); |
| 188 |
x = p->fX0; |
| 189 |
y = p->fY0; |
| 190 |
z = p->fZ0; |
| 191 |
theta = p->fTHETA; |
| 192 |
phi = p->fPHI; |
| 193 |
pdg = p->fPDG; |
| 194 |
momentum = p->fP0; |
| 195 |
id = p->fID; |
| 196 |
|
| 197 |
|
| 198 |
gpr->SetParticle(pdg); |
| 199 |
gpr->SetParticleID(id); |
| 200 |
gpr->SetDirection(theta,phi); |
| 201 |
gpr->SetMomentum(momentum); |
| 202 |
gpr->SetPosition(x,y,z); |
| 203 |
|
| 204 |
} else { |
| 205 |
//generate kinematics |
| 206 |
|
| 207 |
pdg = omg->GetPrimPDG(); |
| 208 |
id = proc_ev; |
| 209 |
gpr->SetParticle(pdg); |
| 210 |
gpr->SetParticleID(id); |
| 211 |
|
| 212 |
//generate spatial distribution |
| 213 |
if( omg->IsCustomKinematics() ){ |
| 214 |
//vertex |
| 215 |
if( omg->IsBoxVertex() ){ |
| 216 |
//box |
| 217 |
omg->GetBoxVertex(xl,xh,yl,yh,zl,zh); |
| 218 |
gpr->GenPosition(xl, xh, yl, yh, zl, zh); |
| 219 |
} else { |
| 220 |
//point |
| 221 |
omg->GetFixedPointVertex(x,y,z); |
| 222 |
gpr->SetPosition(x,y,z); |
| 223 |
} |
| 224 |
//angles |
| 225 |
if( omg->IsIsotropicPlane() ){ |
| 226 |
//isoplane angles |
| 227 |
omg->GetIsotropicAngles(tl, th, pl, ph); |
| 228 |
gpr->GenDirection(tl, th, pl, ph); |
| 229 |
} else { |
| 230 |
//fixed angles |
| 231 |
omg->GetFixedAngles(theta, phi); |
| 232 |
gpr->SetDirection(theta,phi); |
| 233 |
} |
| 234 |
} else { |
| 235 |
// GOD mode |
| 236 |
gpr->Gen2PiPosAng(); |
| 237 |
} |
| 238 |
gpr->GetPosition(x,y,z); |
| 239 |
gpr->GetDirection(theta, phi); |
| 240 |
|
| 241 |
|
| 242 |
//generate MOMENTUM |
| 243 |
Momentum_units mu = omg->GetPrimMomUnits(); |
| 244 |
if( omg->GetPrimMomMode() == FIXED ){ |
| 245 |
momentum = omg->GetMomFixed(); |
| 246 |
|
| 247 |
switch (mu){ |
| 248 |
case GeV_T: |
| 249 |
gpr->SetMomentum(gpr->KinEToMomentum(momentum)); |
| 250 |
break; |
| 251 |
case GV_p: |
| 252 |
gpr->SetMomentum(momentum); |
| 253 |
case GV_R: |
| 254 |
gpr->SetMomentum( gpr->RigToMomentum(momentum)); |
| 255 |
default: |
| 256 |
break; |
| 257 |
} |
| 258 |
|
| 259 |
} |
| 260 |
if( omg->GetPrimMomMode() == FLAT ){ |
| 261 |
Double_t m1, m2; |
| 262 |
omg->GetMomLimits(m1,m2); |
| 263 |
switch (mu){ |
| 264 |
case GeV_T: |
| 265 |
gpr->GenSpe( m1, m2, kTRUE); |
| 266 |
break; |
| 267 |
case GV_p: |
| 268 |
gpr->GenSpe( m1, m2, kFALSE); |
| 269 |
case GV_R: |
| 270 |
gpr->GenSpe( gpr->RigToMomentum(m1), gpr->RigToMomentum(m2), kFALSE); |
| 271 |
default: |
| 272 |
break; |
| 273 |
} |
| 274 |
|
| 275 |
} |
| 276 |
if( omg->GetPrimMomMode() == POWERLAW ){ |
| 277 |
Double_t m1, m2; |
| 278 |
omg->GetMomLimits(m1,m2); |
| 279 |
Double_t gamma = omg->GetMomGamma(); |
| 280 |
switch (mu){ |
| 281 |
case GeV_T: |
| 282 |
gpr->GenSpe( m1, m2, gamma, kTRUE); |
| 283 |
break; |
| 284 |
case GV_p: |
| 285 |
gpr->GenSpe( m1, m2, gamma, kFALSE); |
| 286 |
case GV_R: |
| 287 |
gpr->GenSpe( gpr->RigToMomentum(m1), gpr->RigToMomentum(m2), gamma, kFALSE); |
| 288 |
default: |
| 289 |
break; |
| 290 |
} |
| 291 |
} |
| 292 |
|
| 293 |
} |
| 294 |
|
| 295 |
/*if (DEBUG){ |
| 296 |
cout<<"EVENT: "<<proc_ev<<endl; |
| 297 |
cout<<"ID: "<<id<<endl; |
| 298 |
cout<<"PDG:"<<pdg<<endl |
| 299 |
<<"X0: "<<x<<" (cm)"<<endl |
| 300 |
<<"Y0: "<<y<<" (cm)"<<endl |
| 301 |
<<"Z0: "<<z<<" (cm)"<<endl |
| 302 |
<<"Theta:"<<theta*180./TMath::Pi()<<" (deg)"<<endl |
| 303 |
<<"Phi:"<<phi*180./TMath::Pi()<<" (deg)"<<endl |
| 304 |
<<"Momentum: "<<momentum<<" (GV)"<<endl; |
| 305 |
} */ |
| 306 |
|
| 307 |
if ( omg->IsReadRandomMode() ){ |
| 308 |
read_rndm_file->cd(); |
| 309 |
UnPackDump(tmpfile, rnd_dump); |
| 310 |
TString rndfile = "/random/resetEngineFrom "+tmpfile; |
| 311 |
//cout<<"reading... "<<rndfile<<endl; |
| 312 |
((TGeant4*)gMC)->ProcessGeantCommand(rndfile); |
| 313 |
} |
| 314 |
|
| 315 |
if( omg->DoWriteRandom() ){ |
| 316 |
//cout<<"saving RND"<<endl; |
| 317 |
((TGeant4*)gMC)->ProcessGeantCommand("/random/setSavingFlag 1"); |
| 318 |
} |
| 319 |
|
| 320 |
gpr->ExtrapolateTrajectory(); //Added by V.Formato on 2012/06/19 |
| 321 |
if(!gpr->ToBeSimulated()){ proc_ev++; continue;} |
| 322 |
|
| 323 |
/////RUN/////// |
| 324 |
appl->RunMC(1); |
| 325 |
///////////// |
| 326 |
proc_ev ++; |
| 327 |
accept_ev = gpr->GetNgood(); |
| 328 |
if (appl->GetTrigger()) trig_ev++; |
| 329 |
|
| 330 |
// if(DEBUG) |
| 331 |
cout<<"TOT_EV: "<<proc_ev<<" TRIG_EV: "<<trig_ev<<" ACCEPT_EV:"<<accept_ev<<endl; |
| 332 |
|
| 333 |
//gObjectTable->Print(); |
| 334 |
|
| 335 |
//saving dumps if any |
| 336 |
if( omg->DoWriteRandom() ){ |
| 337 |
switch( omg->GetSaveCond() ){ |
| 338 |
case EVERYTHING: |
| 339 |
dmp_file->cd(); |
| 340 |
PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); |
| 341 |
save_ev = id; |
| 342 |
t_dump->Fill(); |
| 343 |
break; |
| 344 |
case TRIG_ONLY: |
| 345 |
if(appl->GetTrigger()){ |
| 346 |
dmp_file->cd(); |
| 347 |
PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); |
| 348 |
save_ev = id; |
| 349 |
t_dump->Fill(); |
| 350 |
} |
| 351 |
break; |
| 352 |
case ACCEPT_ONLY: |
| 353 |
if(gpr->Getgood()){ |
| 354 |
dmp_file->cd(); |
| 355 |
PackDump(tmpseeddir+"/currentEvent.rndm", ev_dump); |
| 356 |
save_ev = id; |
| 357 |
t_dump->Fill(); |
| 358 |
} |
| 359 |
break; |
| 360 |
default: |
| 361 |
break; |
| 362 |
} |
| 363 |
} |
| 364 |
|
| 365 |
//condition to stop |
| 366 |
switch( omg->GetStopMode() ){ |
| 367 |
case NEV_TOT: |
| 368 |
if( proc_ev >= nev_max) stop = kTRUE; |
| 369 |
break; |
| 370 |
case NEV_ACCEPT: |
| 371 |
if( accept_ev >= nev_max) stop = kTRUE; |
| 372 |
break; |
| 373 |
case NEV_TRIGG: |
| 374 |
if( trig_ev >= nev_max) stop = kTRUE; |
| 375 |
break; |
| 376 |
default: |
| 377 |
break; |
| 378 |
} |
| 379 |
|
| 380 |
if( (omg->GetReadPrimMode()) && ( proc_ev >= nev_max) ) stop = kTRUE; |
| 381 |
|
| 382 |
//gSystem->Exec(" ps aux | grep PamVMC.exe" ); |
| 383 |
} |
| 384 |
gSystem->Exec(" ps aux | grep PamVMC.exe" ); |
| 385 |
|
| 386 |
|
| 387 |
appl->FinishRun(); |
| 388 |
if( omg->DoWriteRandom() ){ |
| 389 |
dmp_file->cd(); |
| 390 |
t_dump->Write(); |
| 391 |
dmp_file->Close(); |
| 392 |
} |
| 393 |
|
| 394 |
cout<<"done...."<<endl; |
| 395 |
|
| 396 |
if( read_rndm_file ) read_rndm_file->Close(); |
| 397 |
delete read_rndm_file; |
| 398 |
delete appl; |
| 399 |
|
| 400 |
//delete tmp directory with dumps |
| 401 |
gSystem->Exec("rm -rf "+tmpseeddir); |
| 402 |
|
| 403 |
//moving good files to good directory |
| 404 |
if( omg->IsUseCookie() ){ |
| 405 |
TString source_dir = working_path; |
| 406 |
TString dest_path = omg->GetOutPath(); |
| 407 |
TString out_root_fname = omg->GetOutPattern()+".root"; |
| 408 |
TString out_pam_fname = omg->GetOutPattern()+".pam"; |
| 409 |
TString out_dump_fname = omg->GetRandomOutFilename(); |
| 410 |
|
| 411 |
cout<<"--->Move "<<source_dir<<"/"<<out_root_fname<<" to "<< dest_path<<endl; |
| 412 |
gSystem->Exec("mv " + source_dir+"/"+out_root_fname + " " + dest_path +"/"); |
| 413 |
cout<<"--->Move "<<source_dir<<"/"<<out_pam_fname<<" to "<< dest_path<<endl; |
| 414 |
gSystem->Exec("mv " + source_dir+"/"+out_pam_fname + " " + dest_path +"/"); |
| 415 |
if( omg->DoWriteRandom() ){ |
| 416 |
cout<<"--->Move "<<out_dump_fname<<" to "<< dest_path<<endl; |
| 417 |
gSystem->Exec("mv " + out_dump_fname + " " + dest_path +"/"); |
| 418 |
} |
| 419 |
cout<<"->Cleaning up tmpdir: "<<source_dir<<endl; |
| 420 |
gSystem->Exec("cd "+dest_path); |
| 421 |
gSystem->Exec("rm -rf " + source_dir ); |
| 422 |
} |
| 423 |
return 0; |
| 424 |
} |
| 425 |
|
| 426 |
|
| 427 |
void PackDump(TString path, TObjString* os){ |
| 428 |
|
| 429 |
ifstream ffile; |
| 430 |
ffile.open(path.Data(), ios::binary); |
| 431 |
if(!ffile) { |
| 432 |
cout<<"ERROR! wrong filename: "<<path<<endl; |
| 433 |
return; |
| 434 |
} |
| 435 |
|
| 436 |
// get length of file: |
| 437 |
ffile.seekg (0, ios::end); |
| 438 |
Int_t length = ffile.tellg(); |
| 439 |
ffile.seekg (0, ios::beg); |
| 440 |
|
| 441 |
// allocate memory: |
| 442 |
char* buffer = new char [length]; |
| 443 |
|
| 444 |
// read data as a block: |
| 445 |
ffile.read (buffer,length); |
| 446 |
ffile.close(); |
| 447 |
|
| 448 |
os->SetString(buffer); |
| 449 |
|
| 450 |
delete[] buffer; |
| 451 |
} |
| 452 |
|
| 453 |
void UnPackDump(TString path, TObjString* os){ |
| 454 |
|
| 455 |
ofstream ffile; |
| 456 |
ffile.open(path.Data(), ios::binary); |
| 457 |
if(!ffile) { |
| 458 |
cout<<"ERROR! wrong filename: "<<path<<endl; |
| 459 |
return; |
| 460 |
} |
| 461 |
|
| 462 |
// get length of file: |
| 463 |
Int_t length = os->Sizeof(); |
| 464 |
|
| 465 |
// allocate memory: |
| 466 |
const char* buffer = 0; |
| 467 |
|
| 468 |
TString tmp = os->GetString(); |
| 469 |
buffer = tmp.Data(); |
| 470 |
|
| 471 |
// read data as a block: |
| 472 |
ffile.write (buffer,length-1); |
| 473 |
|
| 474 |
cout<<"->Unpack"<<endl; |
| 475 |
ffile.close(); |
| 476 |
} |
Parent Directory
| 
Revision Log