nanobeam/src/RunAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Please cite the following paper if you use this software
0027 // Nucl.Instrum.Meth.B260:20-27, 2007
0028
0029 // #define MATRIX_BOUND_CHECK
0030
0031 #include "RunAction.hh"
0032 #include "G4AnalysisManager.hh"
0033 #include "G4AutoLock.hh"
0034
0035 namespace
0036 {
0037   G4Mutex aMutex = G4MUTEX_INITIALIZER;
0038 }
0039
0040 using namespace std;
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0043
0044 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* pri)
0045 :fDetector(det),fPrimary(pri)
0046 {}
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0049
0050 RunAction::~RunAction()
0051 {}
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0054
0055 void RunAction::BeginOfRunAction(const G4Run* /*aRun*/)
0056 {
0057    // Vector initialization
0058
0059    fRow=0;
0060
0061    fXVector = CLHEP::HepVector(32);
0062    fYVector = CLHEP::HepVector(32);
0063    fThetaVector = CLHEP::HepVector(32);
0064    fPhiVector = CLHEP::HepVector(32);
0065
0066    // Histograms
0067
0068    // Get/create analysis manager
0069    G4cout << "##### Create analysis manager " << "  " << this << G4endl;
0070
0071    G4AnalysisManager* man = G4AnalysisManager::Instance();
0072    man->SetDefaultFileType("root");
0073
0074    G4cout << "Using " << man->GetType() << " analysis manager" << G4endl;
0075
0076    // Open an output file
0077    man->OpenFile("nanobeam");
0078    man->SetFirstHistoId(1);
0079    man->SetFirstNtupleId(1);
0080
0081    // Create 1st ntuple (id = 1)
0082    man->CreateNtuple("ntuple0", "BeamProfile");
0083    man->CreateNtupleDColumn("xIn");
0084    man->CreateNtupleDColumn("yIn");
0085    man->CreateNtupleDColumn("zIn");
0086    man->FinishNtuple();
0087    G4cout << "Ntuple-1 created" << G4endl;
0088
0089    // Create 2nd htuple (id = 2)
0090    man->CreateNtuple("ntuple1","Grid");
0091    man->CreateNtupleDColumn("xIn");
0092    man->CreateNtupleDColumn("yIn");
0093    man->CreateNtupleDColumn("e");
0094    man->FinishNtuple();
0095    G4cout << "Ntuple-2 created" << G4endl;
0096
0097    // Create 3rd ntuple (id = 3)
0098    man->CreateNtuple("ntuple2","Coef");
0099    man->CreateNtupleDColumn("xIn");
0100    man->CreateNtupleDColumn("yIn");
0101    man->CreateNtupleDColumn("thetaIn");
0102    man->CreateNtupleDColumn("phiIn");
0103    man->FinishNtuple();
0104    G4cout << "Ntuple-3 created" << G4endl;
0105
0106    return;
0107
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0111
0112 void RunAction::EndOfRunAction(const G4Run* /*aRun*/)
0113 {
0114
0115 if (fDetector->GetCoef()==1)
0116 {
0117         // 17/12/2013 - thanks to A. Dotti
0118     // CLHEP Matrix inversion (as for CLHEP version 2.1.4.1)
0119         // is not thread safe, we need to protect this code.
0120         // Since this is not performance critical we simply lock
0121         // all this part.
0122     G4AutoLock l(&aMutex);
0123         //
0124
0125     CLHEP::HepMatrix m;
0126
0127     // VECTOR READING
0128     // VECTOR READING
0129
0130     m = CLHEP::HepMatrix(32,32);
0131     m = fPrimary->GetMatrix();
0132
0133     G4cout << G4endl;
0134     G4cout << "===> NANOBEAM LINE INTRINSIC ABERRATION COEFFICIENTS (units of micrometer and mrad) :" << G4endl;
0135     G4cout << G4endl;
0136
0137     int inv;
0138
0139     m.invert(inv);
0140     CLHEP::HepVector tmp(32,0);
0141     tmp=m*fXVector;
0142     CLHEP::HepVector b;
0143     b=tmp.sub(2,2);   G4cout << "<x|theta>=" << b << G4endl;
0144     b=tmp.sub(8,8);   G4cout << "<x|theta*delta>=" << b << G4endl;
0145     b=tmp.sub(10,10); G4cout << "<x|theta^3>=" << b << G4endl;
0146     b=tmp.sub(12,12); G4cout << "<x|theta*phi^2>=" << b << G4endl;
0147     m.invert(inv);
0148
0149     m.invert(inv);
0150     tmp = m*fThetaVector;
0151     m.invert(inv);
0152     b=tmp.sub(2,2); G4cout << "<x|x>=" << b << G4endl;
0153
0154     m.invert(inv);
0155     tmp=m*fYVector;
0156     b=tmp.sub(3,3);   G4cout << "<y|phi>=" << b << G4endl;
0157     b=tmp.sub(9,9);   G4cout << "<y|phi*delta>=" << b << G4endl;
0158     b=tmp.sub(11,11); G4cout << "<y|theta^2*phi>=" << b << G4endl;
0159     b=tmp.sub(13,13); G4cout << "<y|phi^3>=" << b << G4endl;
0160     m.invert(inv);
0161
0162     m.invert(inv);
0163     tmp = m*fPhiVector;
0164     m.invert(inv);
0165     b=tmp.sub(3,3); G4cout << "<y|y>=" << b << G4endl;
0166
0167 }
0168
0169  // Save histograms
0170
0171  G4AnalysisManager* man = G4AnalysisManager::Instance();
0172  man->Write();
0173  man->CloseFile();
0174
0175  // Complete clean-up
0176  man->Clear();
0177 }