PamVMC/src/PamVMCApplication.cxx

// $Id: PamVMCApplication.cxx,v 1.0 2007/06/01 
//
// Geant4 PAMELA VMC 
#include <cstdlib>
#include <string>
#include <iostream>
#include <TROOT.h>
#include <TInterpreter.h>
#include <TVirtualMC.h>
#include <Riostream.h>
#include <TGeoManager.h>
#include <TVirtualGeoTrack.h>
#include <TLorentzVector.h>
#include "PamVMCApplication.h"
#include "PamVMCStack.h"
#include "PamVMCPrimaryGenerator.h"

#include "PamVMCFieldMgr.h"

#include "PamVMCSDMgr.h"

#include "PamVMCVolCrossMgr.h"

#include "PamVMCDetPamela.h"

//#include "TFluka.h"


ClassImp(PamVMCApplication)

PamVMCApplication::PamVMCApplication(const char *name, const char *title, const char* filename, Int_t seed) 
  : TVirtualMCApplication(name,title),
    fEventNo(0),
    fVerbose(2),
    fStack(0),
    fdetector(0),
    fPrimaryGenerator(0),
    fRootManager(filename,kWrite)
{
// Standard constructor

  // Create a user stack
  fStack = new PamVMCStack(1000);
    
  // Create a primary generator
  fPrimaryGenerator = new PamVMCPrimaryGenerator(fStack);

  //Create random object
  frandom = new TRandom3(seed);
  std::cout<<"!!! THE SEED of this MONTE CARLO is: "<<frandom->GetSeed()<<std::endl;

  fdetector = new PamVMCDetPamela();

  PamVMCFieldMgr::Instance()->LoadField();
  
  TString fname=filename; fname+=".pam";

  PamVMCRawMgr::Instance()-> CreateOutFile(fname.Data());
}


PamVMCApplication::~PamVMCApplication() 
{
// Destructor  
  delete frandom;
  delete fStack;
  delete fdetector;
  delete fPrimaryGenerator;
  delete gMC;
  gMC = 0;
}


void PamVMCApplication::InitMC(const char* setup)
{    
// Initialize MC.

  fVerbose.InitMC();

  gROOT->LoadMacro(setup);
  gInterpreter->ProcessLine("Config()");
  gMC->SetStack(fStack);
  gMC->SetRandom(frandom);
  
  std::cout<<"!!! BEFORE INIT GMC"<<std::endl;
  
  gMC->Init();
  cout<<"INIT MC"<<endl;
  gMC->BuildPhysics();
  cout<<"PHYSICS BUILDED"<<endl;

#if ROOT_VERSION_CODE >= 333572    
  gMC->SetMagField(PamVMCFieldMgr::Instance());
  cout<<"Magnetic field setted"<<endl;
#endif

  fdetector->InitMC();
  
  // fdetector->Print();
  
  PamVMCSDMgr::Instance()->Register();
  
  PamVMCDigMgr::Instance()->Initialize(gMC->GetRandom());

  PamVMCRawMgr::Instance()->WriteToFile();
  
  fPrimaryGenerator->SetRandom(gMC->GetRandom());
  fPrimaryGenerator->Register();
  
  std::cout<<"LIST OF ACTIVE SD'S:"<<std::endl;
  
  PamVMCSDMgr::Instance()->Print();
}

void PamVMCApplication::RunMC(Int_t nofEvents)
{    
   // MC run.

  fVerbose.RunMC(nofEvents);

  gMC->ProcessRun(nofEvents);

  fVerbose.FinishRun();
}

void PamVMCApplication::AddIons(){

  fVerbose.AddIons();
  //Arguments are:
  //name
  //Z (atomic number)
  //A (atomic mass)
  //Q (charge eplus)
  //exitation exitation energy [GeV] 0. if in ground state
  //mass [GeV], if not specified, calculated approx

  gMC->DefineIon("He3",2,3,2,0.,2.808);
  gMC->DefineIon("He4",2,4,2,0.,3.727); 

  //gMC->DefineIon("He7",2,7,2,0.);
  //gMC->DefineIon("Boron",5,10,5,0.);  
  //gMC->DefineIon("Carbon",6,12,6,0.);
}


void PamVMCApplication::ConstructGeometry()
{    
  // Construct geometry using detector contruction class

  fVerbose.ConstructGeometry();
  new TGeoManager("PAM_GEOM","PAMELA GEOMETRY");
  fdetector->ConstructGeometry();
  //TGeoManager::Import("/home/pamela/PamCAL210508/gp2tgeo.root");
  gGeoManager->CloseGeometry();
  gMC->SetRootGeometry();

  cout<<"!!CONSTRUCT GEOMETRY!!"<<endl;
  fdetector->ConstructCuts();

}

void PamVMCApplication::InitGeometry()
{    
  // Initialize geometry
  
  cout<<"!!SETUP CUTS"<<endl;
  //fdetector->InitGeometry();

  fVerbose.InitGeometry(); 
}


void PamVMCApplication::GeneratePrimaries()
{    
  // Fill the user stack (derived from TVirtualMCStack) with primary particles.
  
  fPrimaryGenerator->GeneratePrimary();
}

void PamVMCApplication::BeginEvent()
{    
// User actions at beginning of event
  PamVMCSDMgr::Instance()->PreEvent();
  fEventNo++;
}

void PamVMCApplication::BeginPrimary()
{    
// User actions at beginning of a primary track
  fPrimaryGenerator->SetGood(kTRUE);
  fVerbose.BeginPrimary();
}

void PamVMCApplication::PreTrack()
{    
// User actions at beginning of each track
 
}


#define CAVITY_X 8.07
#define CAVITY_Y 6.57
#define ZCAVITY_MIN 27.399
#define ZCAVITY_MAX 71.059


Bool_t PamVMCApplication::IsInsideCavity(){
  Double_t x,y,z;
  gMC->TrackPosition(x,y,z);
  if( (z<ZCAVITY_MIN) || (z>ZCAVITY_MAX)) return kTRUE;
  if( (TMath::Abs(x) > CAVITY_X) || (TMath::Abs(y) > CAVITY_Y)) return kFALSE;
  return kTRUE;
}


void PamVMCApplication::Stepping()
{    
  if (fStack->GetCurrentTrack()->IsPrimary()){
    if(!IsInsideCavity()) fPrimaryGenerator->SetGood(kFALSE);
  }
  
// User actions at each step
  fVerbose.Stepping();
  PamVMCDetectorSD* temp=(PamVMCDetectorSD*) PamVMCSDMgr::Instance()->GetSD(gMC->CurrentVolName());
  if(temp){
    fdstatus=INSIDE;
    if(gMC->IsTrackEntering()){ fdstatus=ENTERING; }
    if(gMC->IsTrackExiting() || gMC->IsTrackOut() || gMC->IsTrackStop() ) { fdstatus=EXITING;}
    temp->FillHit(fdstatus,gMC,fStack->GetCurrentTrack()->IsPrimary());    
  }
  
}

void PamVMCApplication::PostTrack()
{    
// User actions after finishing of each track

}

void PamVMCApplication::FinishPrimary()
{    
// User actions after finishing of a primary track
  //fPrimaryGenerator->SetGood(PamVMCVolCrossMgr::Instance()->IsTrackGood());
  //PamVMCVolCrossMgr::Instance()->Reset();
  fVerbose.FinishPrimary();
}

void PamVMCApplication::FinishEvent()
{    
// User actions after finishing of an event

    
    PamVMCSDMgr::Instance()->Digitize();
    //EventNo needs for trigger, particle PDG needs for TOF..
    PamVMCDigMgr::Instance()->Digitize(fEventNo,fPrimaryGenerator->GetParticle());

    PamVMCSDMgr::Instance()->Compress();

    PamVMCRawMgr::Instance()->WriteEvent();

    fRootManager.Fill();

    fPrimaryGenerator->ClearPrimCol();

    PamVMCSDMgr::Instance()->Clear();

    fStack->Reset();
   
} 

void PamVMCApplication::FinishRun()
{
  fVerbose.FinishRun();
  
  fRootManager.WriteAll();

  PamVMCDigMgr::Instance()->FinishRun(); //write RunTrailer to buffer

  PamVMCRawMgr::Instance()->FinishRun(); //write RunTrailer to file and close;

  PamVMCDigMgr::Instance()->PrintCollections();

}

#if ROOT_VERSION_CODE < 333572  
void PamVMCApplication::Field(const Double_t* x, Double_t* b) const
{
  PamVMCFieldMgr::Instance()->Field(x,b);
}
#endif
1	nikolas	1.1	// $Id: PamVMCApplication.cxx,v 1.0 2007/06/01
2			//
3	nikolas	1.2	// Geant4 PAMELA VMC
4	nikolas	1.1	#include <cstdlib>
5			#include <string>
6	nikolas	1.2	#include <iostream>
7	nikolas	1.1	#include <TROOT.h>
8			#include <TInterpreter.h>
9			#include <TVirtualMC.h>
10			#include <Riostream.h>
11			#include <TGeoManager.h>
12			#include <TVirtualGeoTrack.h>
13	nikolas	1.2	#include <TLorentzVector.h>
14	nikolas	1.1	#include "PamVMCApplication.h"
15			#include "PamVMCStack.h"
16			#include "PamVMCPrimaryGenerator.h"
17
18	nikolas	1.2	#include "PamVMCFieldMgr.h"
19
20			#include "PamVMCSDMgr.h"
21
22			#include "PamVMCVolCrossMgr.h"
23
24			#include "PamVMCDetPamela.h"
25
26			//#include "TFluka.h"
27
28
29
30	nikolas	1.1	ClassImp(PamVMCApplication)
31
32	pam-rm2	1.5	PamVMCApplication::PamVMCApplication(const char name, const char title, const char* filename, Int_t seed)
33	nikolas	1.1	: TVirtualMCApplication(name,title),
34			fEventNo(0),
35	nikolas	1.2	fVerbose(2),
36	nikolas	1.1	fStack(0),
37	nikolas	1.2	fdetector(0),
38			fPrimaryGenerator(0),
39	pam-rm2	1.5	fRootManager(filename,kWrite)
40	nikolas	1.1	{
41			// Standard constructor
42
43			// Create a user stack
44	nikolas	1.2	fStack = new PamVMCStack(1000);
45	pam-rm2	1.5
46	nikolas	1.1	// Create a primary generator
47			fPrimaryGenerator = new PamVMCPrimaryGenerator(fStack);
48
49	pam-rm2	1.5	//Create random object
50			frandom = new TRandom3(seed);
51	pizzolot	1.6	std::cout<<"!!! THE SEED of this MONTE CARLO is: "<<frandom->GetSeed()<<std::endl;
52	nikolas	1.2
53			fdetector = new PamVMCDetPamela();
54
55			PamVMCFieldMgr::Instance()->LoadField();
56	pam-rm2	1.5
57			TString fname=filename; fname+=".pam";
58
59			PamVMCRawMgr::Instance()-> CreateOutFile(fname.Data());
60	nikolas	1.1	}
61
62
63			PamVMCApplication::~PamVMCApplication()
64			{
65			// Destructor
66	pam-rm2	1.5	delete frandom;
67	nikolas	1.1	delete fStack;
68	nikolas	1.2	delete fdetector;
69	nikolas	1.1	delete fPrimaryGenerator;
70			delete gMC;
71			gMC = 0;
72			}
73
74
75			void PamVMCApplication::InitMC(const char* setup)
76			{
77			// Initialize MC.
78
79			fVerbose.InitMC();
80
81	nikolas	1.3	gROOT->LoadMacro(setup);
82			gInterpreter->ProcessLine("Config()");
83			gMC->SetStack(fStack);
84	pam-rm2	1.5	gMC->SetRandom(frandom);
85
86	nikolas	1.3	std::cout<<"!!! BEFORE INIT GMC"<<std::endl;
87
88			gMC->Init();
89			cout<<"INIT MC"<<endl;
90			gMC->BuildPhysics();
91			cout<<"PHYSICS BUILDED"<<endl;
92	pam-rm2	1.5
93			#if ROOT_VERSION_CODE >= 333572
94			gMC->SetMagField(PamVMCFieldMgr::Instance());
95			cout<<"Magnetic field setted"<<endl;
96			#endif
97
98	nikolas	1.3	fdetector->InitMC();
99
100			// fdetector->Print();
101
102			PamVMCSDMgr::Instance()->Register();
103
104	pam-rm2	1.5	PamVMCDigMgr::Instance()->Initialize(gMC->GetRandom());
105
106	nikolas	1.3	PamVMCRawMgr::Instance()->WriteToFile();
107
108	pam-rm2	1.5	fPrimaryGenerator->SetRandom(gMC->GetRandom());
109	nikolas	1.3	fPrimaryGenerator->Register();
110
111			std::cout<<"LIST OF ACTIVE SD'S:"<<std::endl;
112
113			PamVMCSDMgr::Instance()->Print();
114	nikolas	1.1	}
115
116			void PamVMCApplication::RunMC(Int_t nofEvents)
117			{
118	nikolas	1.2	// MC run.
119	nikolas	1.1
120			fVerbose.RunMC(nofEvents);
121
122			gMC->ProcessRun(nofEvents);
123
124	nikolas	1.2	fVerbose.FinishRun();
125			}
126
127			void PamVMCApplication::AddIons(){
128
129			fVerbose.AddIons();
130	pam-rm2	1.5	//Arguments are:
131			//name
132			//Z (atomic number)
133			//A (atomic mass)
134			//Q (charge eplus)
135			//exitation exitation energy [GeV] 0. if in ground state
136			//mass [GeV], if not specified, calculated approx
137
138			gMC->DefineIon("He3",2,3,2,0.,2.808);
139			gMC->DefineIon("He4",2,4,2,0.,3.727);
140
141			//gMC->DefineIon("He7",2,7,2,0.);
142			//gMC->DefineIon("Boron",5,10,5,0.);
143			//gMC->DefineIon("Carbon",6,12,6,0.);
144	nikolas	1.1	}
145
146	nikolas	1.2
147	nikolas	1.1	void PamVMCApplication::ConstructGeometry()
148			{
149	nikolas	1.2	// Construct geometry using detector contruction class
150	nikolas	1.1
151			fVerbose.ConstructGeometry();
152	nikolas	1.2	new TGeoManager("PAM_GEOM","PAMELA GEOMETRY");
153			fdetector->ConstructGeometry();
154			//TGeoManager::Import("/home/pamela/PamCAL210508/gp2tgeo.root");
155			gGeoManager->CloseGeometry();
156			gMC->SetRootGeometry();
157
158			cout<<"!!CONSTRUCT GEOMETRY!!"<<endl;
159			fdetector->ConstructCuts();
160
161	nikolas	1.1	}
162
163			void PamVMCApplication::InitGeometry()
164			{
165	nikolas	1.2	// Initialize geometry
166	nikolas	1.1
167	nikolas	1.2	cout<<"!!SETUP CUTS"<<endl;
168			//fdetector->InitGeometry();
169
170			fVerbose.InitGeometry();
171	nikolas	1.1	}
172
173
174			void PamVMCApplication::GeneratePrimaries()
175			{
176	nikolas	1.2	// Fill the user stack (derived from TVirtualMCStack) with primary particles.
177	nikolas	1.1
178	nikolas	1.2	fPrimaryGenerator->GeneratePrimary();
179	nikolas	1.1	}
180
181			void PamVMCApplication::BeginEvent()
182			{
183			// User actions at beginning of event
184	nikolas	1.2	PamVMCSDMgr::Instance()->PreEvent();
185	nikolas	1.1	fEventNo++;
186	nikolas	1.2	}
187	nikolas	1.1
188			void PamVMCApplication::BeginPrimary()
189			{
190			// User actions at beginning of a primary track
191	pam-rm2	1.5	fPrimaryGenerator->SetGood(kTRUE);
192	nikolas	1.2	fVerbose.BeginPrimary();
193	nikolas	1.1	}
194
195			void PamVMCApplication::PreTrack()
196			{
197			// User actions at beginning of each track
198	nikolas	1.2
199	nikolas	1.1	}
200
201	pam-rm2	1.5
202			#define CAVITY_X 8.07
203			#define CAVITY_Y 6.57
204			#define ZCAVITY_MIN 27.399
205			#define ZCAVITY_MAX 71.059
206
207
208			Bool_t PamVMCApplication::IsInsideCavity(){
209			Double_t x,y,z;
210			gMC->TrackPosition(x,y,z);
211			if( (z<ZCAVITY_MIN) \|\| (z>ZCAVITY_MAX)) return kTRUE;
212			if( (TMath::Abs(x) > CAVITY_X) \|\| (TMath::Abs(y) > CAVITY_Y)) return kFALSE;
213			return kTRUE;
214			}
215
216
217	nikolas	1.1	void PamVMCApplication::Stepping()
218			{
219	nikolas	1.2	if (fStack->GetCurrentTrack()->IsPrimary()){
220	pam-rm2	1.5	if(!IsInsideCavity()) fPrimaryGenerator->SetGood(kFALSE);
221	nikolas	1.2	}
222	pam-rm2	1.5
223	nikolas	1.1	// User actions at each step
224	nikolas	1.2	fVerbose.Stepping();
225			PamVMCDetectorSD* temp=(PamVMCDetectorSD*) PamVMCSDMgr::Instance()->GetSD(gMC->CurrentVolName());
226			if(temp){
227			fdstatus=INSIDE;
228			if(gMC->IsTrackEntering()){ fdstatus=ENTERING; }
229			if(gMC->IsTrackExiting() \|\| gMC->IsTrackOut() \|\| gMC->IsTrackStop() ) { fdstatus=EXITING;}
230	pam-rm2	1.5	temp->FillHit(fdstatus,gMC,fStack->GetCurrentTrack()->IsPrimary());
231	nikolas	1.2	}
232
233	nikolas	1.1	}
234
235			void PamVMCApplication::PostTrack()
236			{
237			// User actions after finishing of each track
238
239			}
240
241			void PamVMCApplication::FinishPrimary()
242			{
243			// User actions after finishing of a primary track
244	pam-rm2	1.5	//fPrimaryGenerator->SetGood(PamVMCVolCrossMgr::Instance()->IsTrackGood());
245			//PamVMCVolCrossMgr::Instance()->Reset();
246	nikolas	1.2	fVerbose.FinishPrimary();
247	nikolas	1.1	}
248
249			void PamVMCApplication::FinishEvent()
250			{
251			// User actions after finishing of an event
252
253	nikolas	1.2
254
255			PamVMCSDMgr::Instance()->Digitize();
256			//EventNo needs for trigger, particle PDG needs for TOF..
257			PamVMCDigMgr::Instance()->Digitize(fEventNo,fPrimaryGenerator->GetParticle());
258
259			PamVMCSDMgr::Instance()->Compress();
260	nikolas	1.1
261	nikolas	1.2	PamVMCRawMgr::Instance()->WriteEvent();
262
263			fRootManager.Fill();
264
265			fPrimaryGenerator->ClearPrimCol();
266
267			PamVMCSDMgr::Instance()->Clear();
268
269			fStack->Reset();
270
271			}
272
273			void PamVMCApplication::FinishRun()
274			{
275			fVerbose.FinishRun();
276	nikolas	1.1
277	nikolas	1.2	fRootManager.WriteAll();
278
279			PamVMCDigMgr::Instance()->FinishRun(); //write RunTrailer to buffer
280
281			PamVMCRawMgr::Instance()->FinishRun(); //write RunTrailer to file and close;
282	nikolas	1.1
283	nikolas	1.2	PamVMCDigMgr::Instance()->PrintCollections();
284	nikolas	1.1
285	nikolas	1.2	}
286	nikolas	1.1
287	pam-rm2	1.5	#if ROOT_VERSION_CODE < 333572
288	nikolas	1.1	void PamVMCApplication::Field(const Double_t* x, Double_t* b) const
289			{
290	nikolas	1.2	PamVMCFieldMgr::Instance()->Field(x,b);
291	nikolas	1.1	}
292	pam-rm2	1.5	#endif