MollerAnalysis.cc

Go to the documentation of this file.

#include "MollerAnalysis.hh"
#include "MollerDetectorHit.hh"
#include "RootAnalysis.hh"

#include "G4SystemOfUnits.hh"
#include "G4ios.hh"

#include "G4RunManager.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Event.hh"
#include "G4Run.hh"
#include "G4Track.hh"
#include "G4Material.hh"
#include "G4ClassificationOfNewTrack.hh"
#include "G4TrackStatus.hh"
#include "G4Step.hh"
#include "G4Types.hh"

#include "G4HCofThisEvent.hh"
#include "G4VHitsCollection.hh"
#include "G4TrajectoryContainer.hh"
#include "G4Trajectory.hh"
#include "G4VVisManager.hh"
#include "G4SDManager.hh"
#include "G4UImanager.hh"
#include "G4UIcommand.hh"

#include "MollerAnalysisMessenger.hh"
#include "MollerDetectorConstruction.hh"

#include <time.h>
#include <iostream>

#include "TROOT.h"
#include "TNtuple.h"
#include "TApplication.h"
#include "TSystem.h"

using namespace std;

MollerAnalysis::MollerAnalysis()
{
  if (gSystem) gSystem->ProcessEvents();
  
  ev_num=0;
  partType=-1;
  splitProcess=-1;
  creatorProcess=-1;
  index=-1;

  //  for (int i=0; i<NUM_DETS; i++) {
  //    hitsCollID[i] = -1;
  //  }

  verboseLevel=0;

  kineE0=0;  kineE1=0;  kineE2=0;

  px0=0;  py0=0;  pz0=0;
  px1=0;  py1=0;  pz1=0;
  px2=0;  py2=0;  pz2=0;

  theta0=0;  theta1=0;  theta2=0;

  x0=0;  y0=0;  z0=0;

  analysisMessenger = new MollerAnalysisMessenger(this);
  //G4cout << "Root file name is "<<rootfileName << G4endl;
}
MollerAnalysis::MollerAnalysis(MollerDetectorConstruction* detector)
:myDetector(detector)
{
  if (gSystem) gSystem->ProcessEvents();

  G4int NUM_COLS = myDetector->getNumColl();

  NUM_DETS = NUM_COLS;
  G4cout<<G4endl<<G4endl<<"NUM_COLS = "<<NUM_COLS<<G4endl<<G4endl;

  ev_num=0;
  partType=-1;
  splitProcess=-1;
  creatorProcess=-1;
  index=-1;

  //  for (int i=0; i<NUM_DETS; i++) {
  //    hitsCollID[i] = -1;
  //  }

  verboseLevel=0;

  kineE0=0;  kineE1=0;  kineE2=0;

  px0=0;  py0=0;  pz0=0;
  px1=0;  py1=0;  pz1=0;
  px2=0;  py2=0;  pz2=0;

  theta0=0;  theta1=0;  theta2=0;

  x0=0;  y0=0;  z0=0;

  analysisMessenger = new MollerAnalysisMessenger(this);
  //G4cout << "Root file name is "<<rootfileName << G4endl;
}

MollerAnalysis::MollerAnalysis(const MollerAnalysis &right)
  : RootAnalysis()
{
  if (gSystem) gSystem->ProcessEvents();
  
  ev_num=right.ev_num;
  partType=right.partType;
  splitProcess=right.splitProcess;
  creatorProcess=right.creatorProcess;
  index=right.index;
  

  //  for (int i=0; i<NUM_DETS; i++) {
  //    hitsCollID[i]=right.hitsCollID[i];
  //  }
  
  verboseLevel=right.verboseLevel;

  kineE0=right.kineE0;  
  kineE1=right.kineE1;  
  kineE2=right.kineE2;
 
  px0=right.px0;  
  py0=right.py0;  
  pz0=right.pz0;
  
  px1=right.px1;  
  py1=right.py1;  
  pz1=right.pz1;

  px2=right.px2;  
  py2=right.py2;  
  pz2=right.pz2;

  theta0=right.theta0;  
  theta1=right.theta1;  
  theta2=right.theta2;

  x0=right.x0;  
  y0=right.y0;  
  z0=right.z0;

}

const MollerAnalysis& MollerAnalysis::operator=(const MollerAnalysis &right)
{
  if (gSystem) gSystem->ProcessEvents();
  
  ev_num=right.ev_num;
  partType=right.partType;
  splitProcess=right.splitProcess;
  creatorProcess=right.creatorProcess;
  index=right.index;
  
  //  for (int i=0; i<NUM_DETS; i++) {
  //    hitsCollID[i]=right.hitsCollID[i];
  //  }
  
  verboseLevel=right.verboseLevel;

  kineE0=right.kineE0;  
  kineE1=right.kineE1;  
  kineE2=right.kineE2;
  
  px0=right.px0;  
  py0=right.py0;  
  pz0=right.pz0;
  
  px1=right.px1;  
  py1=right.py1;  
  pz1=right.pz1;

  px2=right.px2;  
  py2=right.py2;  
  pz2=right.pz2;

  theta0=right.theta0;  
  theta1=right.theta1;  
  theta2=right.theta2;

  x0=right.x0;  
  y0=right.y0;  
  z0=right.z0;

  return *this;

}

MollerAnalysis::~MollerAnalysis()
{
  delete analysisMessenger;
}
void MollerAnalysis::SetRootFileName(const G4String& nam) 
{
  rootfileName = nam;
}


void MollerAnalysis::BeginOfRunAction(const G4Run* /*aRun*/)
{
  if (gSystem) gSystem->ProcessEvents();

  G4cout << "Run Action has been entered\n" << G4endl;

  G4cout << "Root file name is "<<rootfileName << G4endl;


  // This works:
  //G4String filenamechar=("~/scratch/ROOTfiles/"+rootfileName+".root");
  G4String filenamechar=(rootfileName+".root");

  // This works:
  /*
  G4String filenamechar;
  filenamechar="~/scratch/ROOTfiles/";
  filenamechar += rootfileName;
  filenamechar += ".root";
  G4cout<<filenamechar<<G4endl;
  */

  //char filenamechar[200];
  /*
  time_t rawtime;
  tm * ptm;
  time ( &rawtime );
  ptm = gmtime ( &rawtime );
  snprintf(filenamechar,200,"~/scratch/ROOTfiles/coll_test0_%02i%02i%02i_%02i%02i%02i.root",
       ptm->tm_year-100,ptm->tm_mon+1,ptm->tm_mday,(ptm->tm_hour+19)%24,ptm->tm_min,ptm->tm_sec);
  //  snprintf(filenamechar,200,"~/scratch/ROOTfiles/moller_%02i%02i%02i_%02i%02i%02i.root",
  //     ptm->tm_year-100,ptm->tm_mon+1,ptm->tm_mday,(ptm->tm_hour+19)%24,ptm->tm_min,ptm->tm_sec);
  printf(filenamechar);
  */
  //snprintf(filenamechar,200,"~/scratch/ROOTfiles/%s.root",rootfileName);

  hfile = new TFile(filenamechar,"RECREATE","Root/G4 analysis");
  //hfile = new TFile(rootfileName,"RECREATE","Root/G4 analysis");


  string myntuple_string;
  myntuple_string="ion:x:y:z:x0:y0:z0:kineE:px:py:pz:kineE0:px0:py0:pz0:kineE1:px1:py1:pz1:kineE2:px2:py2:pz2:type:volume:theta0:theta1:theta2:ev_num:process:event:creator:hit:kineE_org:theta_org:track:diffXS:totXS:rate";

  char *string_ntup = new char[myntuple_string.length()+1];
  strncpy (string_ntup,myntuple_string.c_str(),myntuple_string.length()+1);

  ntup = new TNtuple("geant","G4 Results",(char *)string_ntup);

  G4cout << "Run Action has been completed\n" << G4endl;
}

void MollerAnalysis::BeginOfEventAction(const G4Event */*anEvent*/)
{
  ev_num++;
  if (gSystem) gSystem->ProcessEvents();

  G4SDManager * SDman = G4SDManager::GetSDMpointer();

  /*  if(hitsCollID[0]<0 || hitsCollID[1]<0 || 
     hitsCollID[2]<0 || hitsCollID[3]<0 || 
     hitsCollID[4]<0 || hitsCollID[5]<0 || 
     hitsCollID[6]<0 || hitsCollID[7]<0 || 
     hitsCollID[8]<0 || hitsCollID[9]<0 || 
     hitsCollID[10]<0 || hitsCollID[11]<0 || 
     hitsCollID[12]<0 || hitsCollID[13]<0 || 
     hitsCollID[14]<0 || hitsCollID[15]<0 || hitsCollID[16]<0
     )
  */ 

  if((int) hitsCollID.size()==0)   
    {
      G4cout<<G4endl<<G4endl<<"Setting hitsCollID"<<G4endl<<G4endl;
     
      for (int i=0; i<NUM_DETS; i++) {
  G4String colNam = "detector";
  colNam+=G4UIcommand::ConvertToString(i+1);
  colNam+="/hitsColl";
  
  G4cout << colNam << G4endl;
  
  hitsCollID.push_back (SDman->GetCollectionID(colNam));
      }

      G4cout<<G4endl<<"Finished setting hitsCollID"<<G4endl<<G4endl;
      
    }


  gRootAnalysis->SetProcess(1);
  /* These are you lines Dustin! Only these three lines. */
  //gRootAnalysis->SetMomentum0();
  //gRootAnalysis->SetMomentum1();
  //gRootAnalysis->SetMomentum2();
  /*                 ****************                    */
}


void MollerAnalysis::UserSteppingAction(const G4Step *aStep)
{
  /** Get data that will go into ROOT tree **/

  G4Track* fTrack = aStep->GetTrack();


  if (fTrack->GetTrackStatus()!=fAlive) { 
 //   printf("track is not 'alive.'\n");
    return; 
  }

  G4int tracking = 1;

  G4StepPoint* thePrePoint = aStep->GetPreStepPoint();
  G4StepPoint* thePostPoint = aStep->GetPostStepPoint();
  G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
  G4Material* material = fTrack->GetMaterial();

  const G4VProcess* originProcess = fTrack->GetCreatorProcess();
  const G4VProcess* stepProcess = thePostPoint->GetProcessDefinedStep();

  G4String process;
  if (originProcess!=0) process = originProcess->GetProcessName();
  G4String process1;
  if (stepProcess!=0) process1 = stepProcess->GetProcessName();


  G4String volume = fTrack->GetVolume()->GetName();
  G4String particle = fTrack->GetDefinition()->GetParticleName();


  if (process=="eIoni") creatorProcess=0;
  else if (process=="conv") creatorProcess=1;
  else if (process=="Decay") creatorProcess=2;
  else if (process=="annihil") creatorProcess=3;
  else if (process=="eBrem") creatorProcess=4;
  else if (process=="phot") creatorProcess=5;
  else if (process=="compt") creatorProcess=6;
  else if (process=="Transportation") creatorProcess=7;
  else if (process=="msc") creatorProcess=8;
  else creatorProcess=9;

  //G4bool ion = process=="eIoni";
  
  if ((particle == "e-") 
      && (fTrack->GetVertexKineticEnergy()==11*GeV) 
      && ((thePrePoint->GetKineticEnergy()-
     thePostPoint->GetKineticEnergy())>1.00*GeV) 
      && (process1=="eIoni"))
    { 
      if (kineE0==0 || kineE1==0) {
  kineE0 = thePrePoint->GetKineticEnergy();
  kineE1 = thePostPoint->GetKineticEnergy();
  px0 = thePrePoint->GetMomentumDirection().x();
  py0 = thePrePoint->GetMomentumDirection().y();
  pz0 = thePrePoint->GetMomentumDirection().z();
  px1 = thePostPoint->GetMomentumDirection().x();
  py1 = thePostPoint->GetMomentumDirection().y();
  pz1 = thePostPoint->GetMomentumDirection().z();
  theta0 = sqrt(px0*px0 + py0*py0);
  theta1 = sqrt(px1*px1 + py1*py1);
  x0 = thePostPoint->GetPosition().x();
  y0 = thePostPoint->GetPosition().y();
  z0 = thePostPoint->GetPosition().z();
      }
      
      if (process1=="eIoni") splitProcess=0;
      else if (process1=="conv") splitProcess=1;
      else if (process1=="Decay") splitProcess=2;
      else if (process1=="annihil") splitProcess=3;
      else if (process1=="eBrem") splitProcess=4;
      else if (process1=="photo") splitProcess=5;
      else if (process1=="compton") splitProcess=6;
      else if (process1=="Transportation") splitProcess=7;
      else splitProcess=8;
    }
  
  if (tracking) {
    
    G4ThreeVector worldPos = fTrack->GetPosition();
    G4ThreeVector worldMomentum = fTrack->GetMomentumDirection();
    G4ThreeVector vertexMomentum = fTrack->GetVertexMomentumDirection();
    
    G4double scat_ang = sqrt(vertexMomentum.x()*vertexMomentum.x() + 
           vertexMomentum.y()*vertexMomentum.y());
    
    /*    if (scat_ang<0.002 || scat_ang>0.020) { 
      return; 
      }*/
    
    G4double kineE = fTrack->GetKineticEnergy();
    G4double kineE0 = fTrack->GetVertexKineticEnergy();
    
    G4bool ion = (process=="eIoni");
    
    if (particle=="e-") partType=0;
    else if (particle=="e+") partType=1;
    else if (particle=="proton") partType=2;
    else if (particle=="antiproton") partType=3;
    else if (particle=="gamma") partType=4;
    else partType=5;
    /*   
    if (partType!=4) {
      
      fntup[0] = (Float_t)ion;
      fntup[1] = (Float_t)worldPos.x();
      fntup[2] = (Float_t)worldPos.y();
      fntup[3] = (Float_t)worldPos.z();
      fntup[4] = 0;
      fntup[5] = 0;
      fntup[6] = 0;
      fntup[7] = (Float_t)kineE;
      fntup[8] = (Float_t)worldMomentum.x();
      fntup[9] = (Float_t)worldMomentum.y();
      fntup[10] = (Float_t)worldMomentum.z();
      
      fntup[11] = (Float_t)kineE0;
      fntup[12] = (Float_t)px0;
      fntup[13] = (Float_t)py0;
      fntup[14] = (Float_t)pz0;
      fntup[15] = (Float_t)kineE1;
      fntup[16] = (Float_t)px1;
      fntup[17] = (Float_t)py1;
      fntup[18] = (Float_t)pz1;
      
      fntup[19] = (Float_t)fTrack->GetVertexKineticEnergy();
      fntup[20] = (Float_t)fTrack->GetVertexMomentumDirection().x();
      fntup[21] = (Float_t)fTrack->GetVertexMomentumDirection().y();
      fntup[22] = (Float_t)fTrack->GetVertexMomentumDirection().z();
      
      fntup[23] = (Float_t)partType;
      fntup[24] = (Float_t)atof(volume);
      fntup[25] = (Float_t)theta0;
      fntup[26] = (Float_t)theta1;
      fntup[27] = (Float_t)theta2;
      fntup[28] = (Float_t)ev_num;
      fntup[29] = (Float_t)ion;
      fntup[30] = (Float_t)ev_num;
      
      
      ntup->Fill(fntup);
    }
    */
  }
    
  /* KILL particle if it hits a detector */
  /* made of tungsten! */
  /*
  //  if(thePrePV->GetName()=="1"||thePrePV->GetName()=="2"||
     thePrePV->GetName()=="3"||thePrePV->GetName()=="4"||
     thePrePV->GetName()=="5"||thePrePV->GetName()=="6"|| 
     thePrePV->GetName()=="7"||thePrePV->GetName()=="8"|| 
     thePrePV->GetName()=="9"||thePrePV->GetName()=="10"|| 
     thePrePV->GetName()=="11"||thePrePV->GetName()=="12"|| 
     thePrePV->GetName()=="13"||thePrePV->GetName()=="14"|| 
     thePrePV->GetName()=="15"||thePrePV->GetName()=="16"||thePrePV->GetName()=="17"
     ){
  */
  //    G4cout<<"Material = "<<material->GetName();
  if((material->GetName()=="Tungsten")||(material->GetName()=="Pb")||(material->GetName()=="Copper")){ 
    if (kineE2 == 0 && process=="eIoni") {
      kineE2 = fTrack->GetVertexKineticEnergy();
      px2 = fTrack->GetVertexMomentumDirection().x();
      py2 = fTrack->GetVertexMomentumDirection().y();
      pz2 = fTrack->GetVertexMomentumDirection().z();
      theta2 = sqrt(px2*px2 + py2*py2);
      x0 = fTrack->GetVertexPosition().x();
      y0 = fTrack->GetVertexPosition().y();
      z0 = fTrack->GetVertexPosition().z();
      
    }
    //    G4cout<<" Stopping the track..."<<G4endl;
    //    fTrack->SetTrackStatus(fStopAndKill);
    fTrack->SetTrackStatus(fKillTrackAndSecondaries);
    return;
  }
  //  G4cout<<endl;
}


void MollerAnalysis::EndOfEventAction(const G4Event *anEvent)
{
  
  if (gSystem) gSystem->ProcessEvents();

  G4HCofThisEvent * HCE = anEvent->GetHCofThisEvent();
  
  event = anEvent->GetEventID();

  std::vector<MollerDetectorHitsCollection*> THC;

 
  if(HCE) {
    for (int i=0; i<NUM_DETS; i++) {
      THC.push_back ((MollerDetectorHitsCollection*)(HCE->GetHC(hitsCollID[i])));
    }
  }
  G4int checkTrack1=0;
  G4int checkTrack2=0;
  //  G4cout<<"Checking for hits..."<<G4endl;
  for (int i=0; i<NUM_DETS; i++) {
    // somehow use: aHit->GetTrackID()==1 and aHit->GetTrackID()==2 with an indicator of first entry for that event/track?
    if(THC[i]){
      int n_hit = THC[i]->entries();
      /* // Just do first hit; rate ~65GHz - missing second electron?
  if (n_hit>0) {
  MollerDetectorHit* aHit = (*THC[i])[0];
  //Decide here what kind of data we want to keep
  if (aHit->GetType()==0) AddData(aHit);
  if (aHit->GetType()==4) AddData(aHit);
  }
  }
  }*/
      //  G4cout<<"Done checking for hits..."<<G4endl;
      checkTrack1=0; //need to reset for each detector...
      checkTrack2=0;
      if (n_hit>0) {
  
  for(int i1=0;i1<n_hit;i1++) {
    MollerDetectorHit* aHit = (*THC[i])[i1];
    //aHit->Print();
    //aHit->Draw();
    //    G4cout<<"Ev "<<event<<", det: "<<i<<", hit: "<<i1<<", track: "<<aHit->GetTrackID()<<", check: "<<checkTrack1<<", "<<checkTrack2<<", "<<G4endl;
    //Try to get only the first hit, but for both primary electrons...
    //Really need to define an ntuple variable for the n_hit somehow...
    //and/or do this checking in the analysis scripts...

    if ((aHit->GetTrackID()==1)&&(checkTrack1==0)) {
      if (aHit->GetType()==0) AddData(aHit);
      if (aHit->GetType()==4) AddData(aHit);
      checkTrack1++;
    } else if ((aHit->GetTrackID()==2)&&(checkTrack2==0)){
      if (aHit->GetType()==0) AddData(aHit);
      if (aHit->GetType()==4) AddData(aHit);      
      checkTrack2++;
    }

  }

      }
    }
  }
      
  gRootAnalysis->SetMomentum0(0,0,0,0);
  gRootAnalysis->SetMomentum1(0,0,0,0);
  gRootAnalysis->SetMomentum2(0,0,0,0);
  gRootAnalysis->SetProcess(0);
  
}



void MollerAnalysis::EndOfRunAction(const G4Run */*aRun*/)
{
  if (gSystem) gSystem->ProcessEvents();
  
  hfile->Write();
  hfile->Close();
  
}