PamVMC/examples/run_g4_GP.C

// $Id: run_g4.C,v 1.0 2007/06/01 
//
// Macro for running G4PAM_VMC with Geant4. 
void run_g4_GP(const TString& configMacro = "../macro/g4Config.C") {

    // Load basic libraries
    gROOT->LoadMacro("../macro/basiclibs.C");
    basiclibs();

    // Load PAMELA libs
    TString PAMLIB=gSystem->Getenv("PAM_LIB");
    gSystem->Load(PAMLIB+"/libyoda.so");
    gSystem->Load(PAMLIB+"/libDarthVader.so");
    gSystem->Load(PAMLIB+"/libPamLevel2.so");
    
    // Load Geant4 libraries
    
    gROOT->LoadMacro("../macro/PamG4Libs.C");
    PamG4Libs();
    
    // Load this example libraries
    TString G4WD=gSystem->Getenv("G4WORKDIR");
    TString PLAT=gSystem->Getenv("PLATFORM");
    gSystem->Load(G4WD+"/lib/tgt_"+PLAT+"/libPamVMC_fc");
    
    // MC application
    PamVMCApplication* appl 
      = new PamVMCApplication("PAMG4_VMC", "PAMELA GEANT4 VMC application");  
     
    
    // Initialize MC
    appl->InitMC(configMacro);
    ((TGeant4*)gMC)->ProcessGeantMacro("../macro/g4config2.in");

    //START READING INPUT FILE

    Int_t Ipa;
    Float_t X0, Y0, Z0, Theta, Phi, P0;

    TTree * h20 = 0;

    TString Filename = "protons_N0_1-2.Beam.root";

    TString HOME=gSystem->Getenv("HOME");

    TSystemDirectory *readdir=new TSystemDirectory("readdir",HOME+"/RES");
    if ((TList*)(readdir->GetListOfFiles())->FindObject(Filename) ){
      cout<< Filename<<" found..."<<endl;
      TFile* inpf = new TFile(HOME+"/RES/"+Filename);
      h20 = (TTree*)inpf->Get("h20");
      if (h20->GetBranch("Ipa")){
        h20->SetBranchAddress("Ipa",&Ipa);
        h20->SetBranchStatus("Ipa",1);
      }
      if (h20->GetBranch("X0")){
        h20->SetBranchAddress("X0",&X0);
        h20->SetBranchStatus("X0",1);
      }
      if (h20->GetBranch("Y0")){
        h20->SetBranchAddress("Y0",&Y0);
        h20->SetBranchStatus("Y0",1);
      }
      if (h20->GetBranch("Z0")){
        h20->SetBranchAddress("Z0",&Z0);
        h20->SetBranchStatus("Z0",1);
      }
      if (h20->GetBranch("P0")){
        h20->SetBranchAddress("P0",&P0);
        h20->SetBranchStatus("P0",1);
      }
      if (h20->GetBranch("Theta")){
        h20->SetBranchAddress("Theta",&Theta);
        h20->SetBranchStatus("Theta",1);
      }
      if (h20->GetBranch("Phi")){
        h20->SetBranchAddress("Phi",&Phi);
        h20->SetBranchStatus("Phi",1);
      }
    }

    Int_t nevents = (Int_t)h20->GetEntries();

     nevents = 1000; 
    cout<<"PROCESS "<<nevents<<" EVENTS"<<endl;

    TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(kProton); //correct
    Double_t mass = particlePDG->Mass();
    Double_t PX, PY, PZ;
    for (Int_t i=0; i<nevents; i++){
      h20->GetEntry(i);
      cout<<"Ipa="<<Ipa<<endl;
      appl->GetPrimaryGenerator()->SetParticle(kProton);
      appl->GetPrimaryGenerator()->SetMomentum(P0);
      cout<<"EKIN "<<appl->GetPrimaryGenerator()->GetKinEnergy()<<" GeV"<<endl;
      cout<<"applP0="<<appl->GetPrimaryGenerator()->GetMomentum()<<" GV"<<endl;
      cout<<"applRIG "<<appl->GetPrimaryGenerator()->GetRigidity()<<" GV"<<endl;
      PX=P0*sin(Theta)*cos(Phi);
      PY=P0*sin(Theta)*sin(Phi);
      PZ=-P0*cos(Theta);
      cout<<"(PX0,PY0,PZ0)="<<PX<<","<<PY<<","<<PZ<<endl;
      appl->GetPrimaryGenerator()->SetDirection(Theta,Phi);
      cout<<"(X0,Y0,Z0)"<<X0<<","<<Y0<<","<<Z0<<endl;
      appl->GetPrimaryGenerator()->SetPosition(X0,Y0,Z0);

      
            appl->RunMC(1);
     
    }

  
   // Visualization settings
    // ((TGeant4*)gMC)->ProcessGeantMacro("../macro/g4vis.in");
    // inpf.cd();
    //   delete inpf;
    // Run MC
    //  appl->RunMC(1000);
      appl->FinishRun();
      ((TGeant4*)gMC)->StartGeantUI();
      //      delete appl;

  }  
1	nikolas	1.1	// $Id: run_g4.C,v 1.0 2007/06/01
2			//
3			// Macro for running G4PAM_VMC with Geant4.
4	nikolas	1.3	void run_g4_GP(const TString& configMacro = "../macro/g4Config.C") {
5	nikolas	1.1
6			// Load basic libraries
7	nikolas	1.2	gROOT->LoadMacro("../macro/basiclibs.C");
8	nikolas	1.1	basiclibs();
9
10			// Load PAMELA libs
11			TString PAMLIB=gSystem->Getenv("PAM_LIB");
12			gSystem->Load(PAMLIB+"/libyoda.so");
13			gSystem->Load(PAMLIB+"/libDarthVader.so");
14			gSystem->Load(PAMLIB+"/libPamLevel2.so");
15
16			// Load Geant4 libraries
17
18	nikolas	1.2	gROOT->LoadMacro("../macro/PamG4Libs.C");
19	nikolas	1.1	PamG4Libs();
20
21			// Load this example libraries
22			TString G4WD=gSystem->Getenv("G4WORKDIR");
23			TString PLAT=gSystem->Getenv("PLATFORM");
24			gSystem->Load(G4WD+"/lib/tgt_"+PLAT+"/libPamVMC_fc");
25
26			// MC application
27			PamVMCApplication* appl
28			= new PamVMCApplication("PAMG4_VMC", "PAMELA GEANT4 VMC application");
29
30
31			// Initialize MC
32			appl->InitMC(configMacro);
33	nikolas	1.2	((TGeant4*)gMC)->ProcessGeantMacro("../macro/g4config2.in");
34	nikolas	1.1
35			//START READING INPUT FILE
36
37			Int_t Ipa;
38			Float_t X0, Y0, Z0, Theta, Phi, P0;
39
40			TTree * h20 = 0;
41
42			TString Filename = "protons_N0_1-2.Beam.root";
43
44			TString HOME=gSystem->Getenv("HOME");
45
46			TSystemDirectory *readdir=new TSystemDirectory("readdir",HOME+"/RES");
47			if ((TList*)(readdir->GetListOfFiles())->FindObject(Filename) ){
48			cout<< Filename<<" found..."<<endl;
49			TFile* inpf = new TFile(HOME+"/RES/"+Filename);
50			h20 = (TTree*)inpf->Get("h20");
51			if (h20->GetBranch("Ipa")){
52			h20->SetBranchAddress("Ipa",&Ipa);
53			h20->SetBranchStatus("Ipa",1);
54			}
55			if (h20->GetBranch("X0")){
56			h20->SetBranchAddress("X0",&X0);
57			h20->SetBranchStatus("X0",1);
58			}
59			if (h20->GetBranch("Y0")){
60			h20->SetBranchAddress("Y0",&Y0);
61			h20->SetBranchStatus("Y0",1);
62			}
63			if (h20->GetBranch("Z0")){
64			h20->SetBranchAddress("Z0",&Z0);
65			h20->SetBranchStatus("Z0",1);
66			}
67			if (h20->GetBranch("P0")){
68			h20->SetBranchAddress("P0",&P0);
69			h20->SetBranchStatus("P0",1);
70			}
71			if (h20->GetBranch("Theta")){
72			h20->SetBranchAddress("Theta",&Theta);
73			h20->SetBranchStatus("Theta",1);
74			}
75			if (h20->GetBranch("Phi")){
76			h20->SetBranchAddress("Phi",&Phi);
77			h20->SetBranchStatus("Phi",1);
78			}
79			}
80
81			Int_t nevents = (Int_t)h20->GetEntries();
82
83			nevents = 1000;
84			cout<<"PROCESS "<<nevents<<" EVENTS"<<endl;
85
86			TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(kProton); //correct
87			Double_t mass = particlePDG->Mass();
88			Double_t PX, PY, PZ;
89			for (Int_t i=0; i<nevents; i++){
90			h20->GetEntry(i);
91			cout<<"Ipa="<<Ipa<<endl;
92			appl->GetPrimaryGenerator()->SetParticle(kProton);
93			appl->GetPrimaryGenerator()->SetMomentum(P0);
94			cout<<"EKIN "<<appl->GetPrimaryGenerator()->GetKinEnergy()<<" GeV"<<endl;
95			cout<<"applP0="<<appl->GetPrimaryGenerator()->GetMomentum()<<" GV"<<endl;
96			cout<<"applRIG "<<appl->GetPrimaryGenerator()->GetRigidity()<<" GV"<<endl;
97			PX=P0sin(Theta)cos(Phi);
98			PY=P0sin(Theta)sin(Phi);
99			PZ=-P0*cos(Theta);
100			cout<<"(PX0,PY0,PZ0)="<<PX<<","<<PY<<","<<PZ<<endl;
101			appl->GetPrimaryGenerator()->SetDirection(Theta,Phi);
102			cout<<"(X0,Y0,Z0)"<<X0<<","<<Y0<<","<<Z0<<endl;
103			appl->GetPrimaryGenerator()->SetPosition(X0,Y0,Z0);
104
105
106			appl->RunMC(1);
107
108			}
109
110
111			// Visualization settings
112	nikolas	1.2	// ((TGeant4*)gMC)->ProcessGeantMacro("../macro/g4vis.in");
113	nikolas	1.1	// inpf.cd();
114			// delete inpf;
115			// Run MC
116			// appl->RunMC(1000);
117			appl->FinishRun();
118			((TGeant4*)gMC)->StartGeantUI();
119			// delete appl;
120
121			}