PamVMC_update/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 "PamVMCDetPamela.h"

#include "PamVMCRndMgr.h"
#include "PamVMCOptMgr.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),
    ftrig(0)
{
  cout<<"---> PamVMCApplication::PamVMCApplication"<<endl;
  // 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();

  cout<<"OK1"<<endl;
  TString fname=filename; fname+=".pam";
  PamVMCOptMgr::Instance();

  if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
    cout<<"OK2"<<endl;
    PamVMCRawMgr::Instance()-> CreateOutFile(fname.Data());
  }
  cout<<"<--- PamVMCApplication::PamVMCApplication"<<endl;
}


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);
  
#if ROOT_VERSION_CODE >= 333572
  gMC->SetMagField(PamVMCFieldMgr::Instance());
  cout<<"Magnetic field set"<<endl;
#endif

  std::cout<<"!!! BEFORE INIT GMC"<<std::endl;
  
  gMC->Init();
  cout<<"INIT MC"<<endl;
  gMC->BuildPhysics();
  cout<<"PHYSICS BUILDED"<<endl;

  fdetector->InitMC();
  
  // fdetector->Print();
  
  if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
    PamVMCSDMgr::Instance()->Register();
    PamVMCDigMgr::Instance()->Initialize((TRandom3*)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("H2",1,2,1,0., 1.875613);
  gMC->DefineIon("H3",1,3,1,0., 2.808921);
  gMC->DefineIon("He3",2,3,2,0.,2.808391);
  gMC->DefineIon("He4",2,4,2,0.,3.727379); 
  gMC->DefineIon("Li6",3,6,3,0.);
  gMC->DefineIon("Li7",3,7,3,0.);
  gMC->DefineIon("Be7",4,7,4,0.);
  gMC->DefineIon("Be9",4,9,4,0.);
  gMC->DefineIon("Be10",4,10,4,0.);
  gMC->DefineIon("B10",5,10,5,0.);
  gMC->DefineIon("B11",5,11,5,0.);
  gMC->DefineIon("C12",6,12,6,0.);
  gMC->DefineIon("N14",7,14,7,0.);
  gMC->DefineIon("O16",8,16,8,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");
  cout<<"--->CONSTRUCT GEOMETRY"<<endl;
  gGeoManager->CloseGeometry();
  gMC->SetRootGeometry();
  cout<<"<---GEOMETRY CONSTRUCTED"<<endl;
  
  //cout<<"--->CONSTRUCT CUTS"<<endl;
  //fdetector->ConstructCuts();
  //cout<<"<---CUTS CONSTRUCTED"<<endl;

}

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
if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
   PamVMCSDMgr::Instance()->PreEvent();
 }
 ftrig = 0;
 fEventNo++;
}

void PamVMCApplication::BeginPrimary()
{    
// User actions at beginning of a primary track
  fVerbose.BeginPrimary();
  //making a clones of random objects states
  if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) PamVMCRndMgr::Instance()->MakeRandomClones();
  Double_t x,y,z;
  gMC->TrackPosition(x,y,z);
  fPrimaryGenerator->SetZeroStep(z,x,y);
}

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


void PamVMCApplication::Stepping()
{    
  if (fStack->GetCurrentTrack()->IsPrimary()){    
    //write some code to call cross manager for primary
    Double_t x,y,z;
    gMC->TrackPosition(x,y,z);
    fPrimaryGenerator->CheckCrossing(z,x,y);
  }
  
// User actions at each step
  fVerbose.Stepping();
  PamVMCDetectorSD* temp=(PamVMCDetectorSD*) PamVMCSDMgr::Instance()->GetSD(gMC->CurrentVolName()); //old

  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
  fVerbose.FinishPrimary();
  //Here I should add stuff to check if primary was inside acceptance
  fPrimaryGenerator->CheckInsideAcceptance();
}

void PamVMCApplication::FinishEvent()
{    
// User actions after finishing of an event
  ftrig = PamVMCSDMgr::Instance()->GetTrigger();

  switch( PamVMCOptMgr::Instance()->GetSaveCond()){
  case EVERYTHING:
    if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
    fRootManager.Fill();
    break;

  case TRIG_ONLY:
    if( ftrig ){
      if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
      fRootManager.Fill();
    }
    break;

  case ACCEPT_ONLY:
    if( fPrimaryGenerator->Getgood() ){
      if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
      fRootManager.Fill();
    }
    break;

  default:
    break;


  }


  fPrimaryGenerator->ClearPrimCol();
  
  //Check for trigger:
  
  
  PamVMCSDMgr::Instance()->Clear();
  if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) PamVMCRndMgr::Instance()->ClearColl();
  
  fStack->Reset(); 
} 


void PamVMCApplication::FillEv(){

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

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

  if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
    PamVMCRawMgr::Instance()->FinishRun();  
    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	formato	1.1	// $Id: PamVMCApplication.cxx,v 1.0 2007/06/01
2			//
3			// Geant4 PAMELA VMC
4			#include <cstdlib>
5			#include <string>
6			#include <iostream>
7			#include <TROOT.h>
8			#include <TInterpreter.h>
9			#include <TVirtualMC.h>
10			#include <Riostream.h>
11			#include <TGeoManager.h>
12			#include <TVirtualGeoTrack.h>
13			#include <TLorentzVector.h>
14			#include "PamVMCApplication.h"
15			#include "PamVMCStack.h"
16			#include "PamVMCPrimaryGenerator.h"
17
18			#include "PamVMCFieldMgr.h"
19
20			#include "PamVMCSDMgr.h"
21
22			#include "PamVMCDetPamela.h"
23
24			#include "PamVMCRndMgr.h"
25			#include "PamVMCOptMgr.h"
26
27
28			ClassImp(PamVMCApplication)
29
30			PamVMCApplication::PamVMCApplication(const char name, const char title, const char* filename, Int_t seed)
31			: TVirtualMCApplication(name,title),
32			fEventNo(0),
33			fVerbose(2),
34			fStack(0),
35			fdetector(0),
36			fPrimaryGenerator(0),
37			fRootManager(filename,kWrite),
38			ftrig(0)
39			{
40			cout<<"---> PamVMCApplication::PamVMCApplication"<<endl;
41			// Create a user stack
42			fStack = new PamVMCStack(1000);
43
44			// Create a primary generator
45			fPrimaryGenerator = new PamVMCPrimaryGenerator(fStack);
46
47			//Create random object
48			frandom = new TRandom3(seed);
49			std::cout<<"!!! THE SEED of this MONTE CARLO is: "<<frandom->GetSeed()<<std::endl;
50
51			fdetector = new PamVMCDetPamela();
52
53			PamVMCFieldMgr::Instance()->LoadField();
54
55			cout<<"OK1"<<endl;
56			TString fname=filename; fname+=".pam";
57			PamVMCOptMgr::Instance();
58
59			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
60			cout<<"OK2"<<endl;
61			PamVMCRawMgr::Instance()-> CreateOutFile(fname.Data());
62			}
63			cout<<"<--- PamVMCApplication::PamVMCApplication"<<endl;
64			}
65
66
67			PamVMCApplication::~PamVMCApplication()
68			{
69			// Destructor
70			delete frandom;
71			delete fStack;
72			delete fdetector;
73			delete fPrimaryGenerator;
74			delete gMC;
75			gMC = 0;
76			}
77
78
79			void PamVMCApplication::InitMC(const char* setup)
80			{
81			// Initialize MC.
82
83			fVerbose.InitMC();
84
85			gROOT->LoadMacro(setup);
86			gInterpreter->ProcessLine("Config()");
87			gMC->SetStack(fStack);
88			gMC->SetRandom(frandom);
89
90			#if ROOT_VERSION_CODE >= 333572
91			gMC->SetMagField(PamVMCFieldMgr::Instance());
92			cout<<"Magnetic field set"<<endl;
93			#endif
94
95			std::cout<<"!!! BEFORE INIT GMC"<<std::endl;
96
97			gMC->Init();
98			cout<<"INIT MC"<<endl;
99			gMC->BuildPhysics();
100			cout<<"PHYSICS BUILDED"<<endl;
101
102			fdetector->InitMC();
103
104			// fdetector->Print();
105
106			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
107			PamVMCSDMgr::Instance()->Register();
108			PamVMCDigMgr::Instance()->Initialize((TRandom3*)gMC->GetRandom());
109			PamVMCRawMgr::Instance()->WriteToFile();
110			}
111			fPrimaryGenerator->SetRandom(gMC->GetRandom());
112
113			fPrimaryGenerator->Register();
114
115			std::cout<<"LIST OF ACTIVE SD'S:"<<std::endl;
116
117			PamVMCSDMgr::Instance()->Print();
118			}
119
120			void PamVMCApplication::RunMC(Int_t nofEvents)
121			{
122			// MC run.
123
124			fVerbose.RunMC(nofEvents);
125
126			gMC->ProcessRun(nofEvents);
127
128			fVerbose.FinishRun();
129			}
130
131			void PamVMCApplication::AddIons(){
132
133			fVerbose.AddIons();
134			//Arguments are:
135			//name
136			//Z (atomic number)
137			//A (atomic mass)
138			//Q (charge eplus)
139			//exitation exitation energy [GeV] 0. if in ground state
140			//mass [GeV], if not specified, calculated approx
141			gMC->DefineIon("H2",1,2,1,0., 1.875613);
142			gMC->DefineIon("H3",1,3,1,0., 2.808921);
143			gMC->DefineIon("He3",2,3,2,0.,2.808391);
144			gMC->DefineIon("He4",2,4,2,0.,3.727379);
145			gMC->DefineIon("Li6",3,6,3,0.);
146			gMC->DefineIon("Li7",3,7,3,0.);
147			gMC->DefineIon("Be7",4,7,4,0.);
148			gMC->DefineIon("Be9",4,9,4,0.);
149			gMC->DefineIon("Be10",4,10,4,0.);
150			gMC->DefineIon("B10",5,10,5,0.);
151			gMC->DefineIon("B11",5,11,5,0.);
152			gMC->DefineIon("C12",6,12,6,0.);
153			gMC->DefineIon("N14",7,14,7,0.);
154			gMC->DefineIon("O16",8,16,8,0.);
155			}
156
157
158			void PamVMCApplication::ConstructGeometry()
159			{
160			// Construct geometry using detector contruction class
161
162			fVerbose.ConstructGeometry();
163			new TGeoManager("PAM_GEOM","PAMELA GEOMETRY");
164			fdetector->ConstructGeometry();
165			//TGeoManager::Import("/home/pamela/PamCAL210508/gp2tgeo.root");
166			cout<<"--->CONSTRUCT GEOMETRY"<<endl;
167			gGeoManager->CloseGeometry();
168			gMC->SetRootGeometry();
169			cout<<"<---GEOMETRY CONSTRUCTED"<<endl;
170
171			//cout<<"--->CONSTRUCT CUTS"<<endl;
172			//fdetector->ConstructCuts();
173			//cout<<"<---CUTS CONSTRUCTED"<<endl;
174
175			}
176
177			void PamVMCApplication::InitGeometry()
178			{
179			// Initialize geometry
180
181			cout<<"!!SETUP CUTS"<<endl;
182			//fdetector->InitGeometry();
183
184			fVerbose.InitGeometry();
185			}
186
187
188			void PamVMCApplication::GeneratePrimaries()
189			{
190			// Fill the user stack (derived from TVirtualMCStack) with primary particles.
191
192			fPrimaryGenerator->GeneratePrimary();
193			}
194
195			void PamVMCApplication::BeginEvent()
196			{
197			// User actions at beginning of event
198			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
199			PamVMCSDMgr::Instance()->PreEvent();
200			}
201			ftrig = 0;
202			fEventNo++;
203			}
204
205			void PamVMCApplication::BeginPrimary()
206			{
207			// User actions at beginning of a primary track
208			fVerbose.BeginPrimary();
209			//making a clones of random objects states
210			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) PamVMCRndMgr::Instance()->MakeRandomClones();
211			Double_t x,y,z;
212			gMC->TrackPosition(x,y,z);
213			fPrimaryGenerator->SetZeroStep(z,x,y);
214			}
215
216			void PamVMCApplication::PreTrack()
217			{
218			// User actions at beginning of each track
219
220			}
221
222
223			#define CAVITY_X 8.07
224			#define CAVITY_Y 6.57
225			#define ZCAVITY_MIN 27.399
226			#define ZCAVITY_MAX 71.059
227
228
229
230			void PamVMCApplication::Stepping()
231			{
232			if (fStack->GetCurrentTrack()->IsPrimary()){
233			//write some code to call cross manager for primary
234			Double_t x,y,z;
235			gMC->TrackPosition(x,y,z);
236			fPrimaryGenerator->CheckCrossing(z,x,y);
237			}
238
239			// User actions at each step
240			fVerbose.Stepping();
241			PamVMCDetectorSD* temp=(PamVMCDetectorSD*) PamVMCSDMgr::Instance()->GetSD(gMC->CurrentVolName()); //old
242
243			if(temp){
244
245			fdstatus=INSIDE;
246			if(gMC->IsTrackEntering()){ fdstatus=ENTERING; }
247			if(gMC->IsTrackExiting() \|\| gMC->IsTrackOut() \|\| gMC->IsTrackStop() ) { fdstatus=EXITING;}
248			temp->FillHit(fdstatus,gMC,fStack->GetCurrentTrack()->IsPrimary());
249			}
250
251			}
252
253			void PamVMCApplication::PostTrack()
254			{
255			// User actions after finishing of each track
256
257			}
258
259			void PamVMCApplication::FinishPrimary()
260			{
261			// User actions after finishing of a primary track
262			fVerbose.FinishPrimary();
263			//Here I should add stuff to check if primary was inside acceptance
264			fPrimaryGenerator->CheckInsideAcceptance();
265			}
266
267			void PamVMCApplication::FinishEvent()
268			{
269			// User actions after finishing of an event
270			ftrig = PamVMCSDMgr::Instance()->GetTrigger();
271
272			switch( PamVMCOptMgr::Instance()->GetSaveCond()){
273			case EVERYTHING:
274			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
275			fRootManager.Fill();
276			break;
277
278			case TRIG_ONLY:
279			if( ftrig ){
280			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
281			fRootManager.Fill();
282			}
283			break;
284
285			case ACCEPT_ONLY:
286			if( fPrimaryGenerator->Getgood() ){
287			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) FillEv();
288			fRootManager.Fill();
289			}
290			break;
291
292			default:
293			break;
294
295
296			}
297
298
299
300
301
302			fPrimaryGenerator->ClearPrimCol();
303
304			//Check for trigger:
305
306
307			PamVMCSDMgr::Instance()->Clear();
308			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS) PamVMCRndMgr::Instance()->ClearColl();
309
310			fStack->Reset();
311			}
312
313
314			void PamVMCApplication::FillEv(){
315
316			PamVMCSDMgr::Instance()->Digitize();
317			//EventNo needs for trigger, particle PDG needs for TOF..
318			PamVMCDigMgr::Instance()->Digitize(fEventNo,fPrimaryGenerator->GetParticle());
319
320			PamVMCSDMgr::Instance()->Compress();
321
322			PamVMCRawMgr::Instance()->WriteEvent();
323			PamVMCRndMgr::Instance()->Compress();
324			}
325
326			void PamVMCApplication::FinishRun()
327			{
328			fVerbose.FinishRun();
329
330			fRootManager.WriteAll();
331
332			if (PamVMCOptMgr::Instance()->GetSaveMode() == ALL_DETECTORS){
333			PamVMCRawMgr::Instance()->FinishRun();
334			PamVMCDigMgr::Instance()->PrintCollections();
335			}
336
337			}
338
339			#if ROOT_VERSION_CODE < 333572
340			void PamVMCApplication::Field(const Double_t* x, Double_t* b) const
341			{
342			PamVMCFieldMgr::Instance()->Field(x,b);
343			}
344			#endif