PamVMC_update/G4main/PamVMCmain.cxx

#include <stdio.h>
#include <stdlib.h>
#include <iostream>

//ROOT headers
#include <TSystem.h>
#include <TSystemFile.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>
#include <TObjString.h>
//#include <TObjectTable.h>

//G4VMC headers
#include <TGeant4.h>

//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"<<endl;
    abort();
  }
  //Option1 montecarlo selection,  default G4
  gSystem->Load("libg4root.so");
  cout<<" Input card file: "<<argv[1]<<endl;

  PamVMCOptMgr * omg = PamVMCOptMgr::Instance();
  omg->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: "<<tmpseeddir<<endl;
    gSystem->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<<out_tmp_fname<<endl;
  PamVMCApplication* appl 
    = new PamVMCApplication("PAMG4_VMC", "PAMELA GEANT4 VMC application",out_tmp_fname,0);
  //setting verbose level
  appl->SetVerboseLevel( omg->GetVerboseLevel() );
  //initialize physycs with specific RunConfiguration, defined in macro
  appl->InitMC( omg->GetG4ConfigC() );
  cout<<"MACRO READED"<<endl;
  //additional flags and controls related to G4 core
  ((TGeant4*)gMC)->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: "<<tmpseeddir<<endl;
    gSystem->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<int>::max());
      seed2 = rinit2.Integer(numeric_limits<int>::max());
    }
    cout<<"Initial RND seeds for G4 core are: " << seed1<<" , " << seed2<<endl;
    TString command = "/random/setSeeds ";
    command += seed1;
    command += " ";
    command += seed2;
    ((TGeant4*)gMC)->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"<<endl;
    primtree->SetBranchAddress("PRIM", &prims);
    nev_max = primtree->GetEntries();
  }

  cout<<"OK1"<<endl;

  //MAIN LOOP
  Int_t proc_ev = 0;
  Int_t accept_ev = 0;
  Int_t trig_ev = 0;
  Bool_t stop = kFALSE;
  TRandom3 * aux_rnd = new TRandom3(0);
  PamVMCPrimaryGenerator* gpr = appl->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: "<<proc_ev<<endl;
      cout<<"ID: "<<id<<endl;
      cout<<"PDG:"<<pdg<<endl
          <<"X0: "<<x<<" (cm)"<<endl
          <<"Y0: "<<y<<" (cm)"<<endl
          <<"Z0: "<<z<<" (cm)"<<endl
          <<"Theta:"<<theta*180./TMath::Pi()<<" (deg)"<<endl
          <<"Phi:"<<phi*180./TMath::Pi()<<" (deg)"<<endl
          <<"Momentum: "<<momentum<<" (GV)"<<endl;
          } */  

    if ( omg->IsReadRandomMode() ){
      read_rndm_file->cd();
      UnPackDump(tmpfile, rnd_dump);
      TString rndfile = "/random/resetEngineFrom "+tmpfile;
      //cout<<"reading... "<<rndfile<<endl;
      ((TGeant4*)gMC)->ProcessGeantCommand(rndfile);
    }

    if( omg->DoWriteRandom() ){
      //cout<<"saving RND"<<endl;
      ((TGeant4*)gMC)->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: "<<proc_ev<<" TRIG_EV: "<<trig_ev<<" ACCEPT_EV:"<<accept_ev<<endl;
    
    //gObjectTable->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...."<<endl;
  
  if( read_rndm_file ) read_rndm_file->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 "<<source_dir<<"/"<<out_root_fname<<" to "<< dest_path<<endl;
    gSystem->Exec("mv " + source_dir+"/"+out_root_fname + " " + dest_path +"/");
    cout<<"--->Move "<<source_dir<<"/"<<out_pam_fname<<" to "<< dest_path<<endl;
    gSystem->Exec("mv " + source_dir+"/"+out_pam_fname + " " + dest_path +"/");
    if( omg->DoWriteRandom() ){
      cout<<"--->Move "<<out_dump_fname<<" to "<< dest_path<<endl;
      gSystem->Exec("mv " + out_dump_fname + " " + dest_path +"/");
    }
    cout<<"->Cleaning up tmpdir: "<<source_dir<<endl;
    gSystem->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: "<<path<<endl;
    return;
  } 

  // get length of file:
  ffile.seekg (0, ios::end);
  Int_t length = ffile.tellg();
  ffile.seekg (0, ios::beg);

  // allocate memory:
  char* buffer = new char [length];

  // read data as a block:
  ffile.read (buffer,length);
  ffile.close();

  os->SetString(buffer);

  delete[] buffer;
}

void UnPackDump(TString path, TObjString* os){

  ofstream ffile;
  ffile.open(path.Data(), ios::binary);
  if(!ffile) {
    cout<<"ERROR! wrong filename: "<<path<<endl;
    return;
  } 
  
  // get length of file:
  Int_t length = os->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"<<endl;
  ffile.close();
}
1	formato	1.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			}