--- PamVMC/src/PamVMCApplication.cxx	2007/06/28 07:16:57	1.1
+++ PamVMC/src/PamVMCApplication.cxx	2010/09/15 07:01:55	1.6
@@ -1,75 +1,76 @@
 // $Id: PamVMCApplication.cxx,v 1.0 2007/06/01 
 //
-// Geant4 ExampleN06 adapted to Virtual Monte Carlo 
+// 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 "PamVMCDetectorConstruction.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) 
+PamVMCApplication::PamVMCApplication(const char *name, const char *title, const char* filename, Int_t seed) 
   : TVirtualMCApplication(name,title),
     fEventNo(0),
-    fVerbose(0),
+    fVerbose(2),
     fStack(0),
-    fDetConstruction(0),
-    fPrimaryGenerator(0)
+    fdetector(0),
+    fPrimaryGenerator(0),
+    fRootManager(filename,kWrite)
 {
 // Standard constructor
 
   // Create a user stack
-    fStack = new PamVMCStack(10000);
-  
-  // Create detector construction
-  fDetConstruction = new PamVMCDetectorConstruction();
-  
+  fStack = new PamVMCStack(1000);
+    
   // Create a primary generator
   fPrimaryGenerator = new PamVMCPrimaryGenerator(fStack);
-#ifdef PAMFIELD
-  // Load the PAMELA magnetic field
 
-  std::string pamcal=getenv("PAM_CALIB"); pamcal+="/trk-param/field_param-0/";
-  std::cout << "PAMELA env: " << pamcal << std::endl;
-  pamfield.LoadField(pamcal.c_str());
-#endif
-}
+  //Create random object
+  frandom = new TRandom3(seed);
+  std::cout<<"!!! THE SEED of this MONTE CARLO is: "<<frandom->GetSeed()<<std::endl;
 
+  fdetector = new PamVMCDetPamela();
 
-PamVMCApplication::PamVMCApplication()
-  : TVirtualMCApplication(),
-    fEventNo(0),
-    fVerbose(0),
-    fStack(0),
-    fDetConstruction(0),
-    fPrimaryGenerator(0)
-{    
-// Default constructor
+  PamVMCFieldMgr::Instance()->LoadField();
+  
+  TString fname=filename; fname+=".pam";
+
+  PamVMCRawMgr::Instance()-> CreateOutFile(fname.Data());
 }
 
+
 PamVMCApplication::~PamVMCApplication() 
 {
 // Destructor  
-  
+  delete frandom;
   delete fStack;
-  delete fDetConstruction;
+  delete fdetector;
   delete fPrimaryGenerator;
   delete gMC;
   gMC = 0;
 }
 
-//
-// public methods
-//
 
 void PamVMCApplication::InitMC(const char* setup)
 {    
@@ -77,117 +78,215 @@
 
   fVerbose.InitMC();
 
-    gROOT->LoadMacro(setup);
-    gInterpreter->ProcessLine("Config()");
-    gMC->SetStack(fStack);
-    gMC->Init();
-    gMC->BuildPhysics(); 
+  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.
+   // MC run.
 
   fVerbose.RunMC(nofEvents);
 
   gMC->ProcessRun(nofEvents);
 
-  ////  fVerbose.FinishRun();
+  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
+  // Construct geometry using detector contruction class
 
   fVerbose.ConstructGeometry();
-  fDetConstruction->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
+  // Initialize geometry
   
-  fVerbose.InitGeometry();
+  cout<<"!!SETUP CUTS"<<endl;
+  //fdetector->InitGeometry();
+
+  fVerbose.InitGeometry(); 
 }
 
 
 void PamVMCApplication::GeneratePrimaries()
 {    
-// Fill the user stack (derived from TVirtualMCStack) with primary particles.
+  // Fill the user stack (derived from TVirtualMCStack) with primary particles.
   
-  ///// fVerbose.GeneratePrimaries();
-
-  fPrimaryGenerator->GeneratePrimaries();
+  fPrimaryGenerator->GeneratePrimary();
 }
 
 void PamVMCApplication::BeginEvent()
 {    
 // User actions at beginning of event
-
-  ///// fVerbose.BeginEvent();
-
+  PamVMCSDMgr::Instance()->PreEvent();
   fEventNo++;
- }
+}
 
 void PamVMCApplication::BeginPrimary()
 {    
 // User actions at beginning of a primary track
-
-////  fVerbose.BeginPrimary();
+  fPrimaryGenerator->SetGood(kTRUE);
+  fVerbose.BeginPrimary();
 }
 
 void PamVMCApplication::PreTrack()
 {    
 // User actions at beginning of each track
+ 
+}
 
-  ////// fVerbose.PreTrack();
-  
+
+#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();
+  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
 
-
-  ///// fVerbose.PostTrack();
 }
 
 void PamVMCApplication::FinishPrimary()
 {    
 // User actions after finishing of a primary track
-
-  ///// fVerbose.FinishPrimary();
+  //fPrimaryGenerator->SetGood(PamVMCVolCrossMgr::Instance()->IsTrackGood());
+  //PamVMCVolCrossMgr::Instance()->Reset();
+  fVerbose.FinishPrimary();
 }
 
 void PamVMCApplication::FinishEvent()
 {    
 // User actions after finishing of an event
 
+    
 
-  /////fVerbose.FinishEvent();
-  
+    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();
 
-  fStack->Reset();
+    PamVMCSDMgr::Instance()->Clear();
+
+    fStack->Reset();
+   
 } 
 
-void PamVMCApplication::Field(const Double_t* x, Double_t* b) const
+void PamVMCApplication::FinishRun()
 {
-// Uniform magnetic field
+  fVerbose.FinishRun();
+  
+  fRootManager.WriteAll();
 
-#ifdef PAMFIELD
-     b[0] = pamfield.GetBX((float *)x);
-     b[1] = pamfield.GetBY((float *)x);
-     b[2] = pamfield.GetBZ((float *)x);
-#else
-   for (Int_t i=0; i<3; i++) b[i] = 0.0;
-#endif
+  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