errorpropagation/errProp/errProp.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 //
0027 /// \file errProp/errProp.cc
0028 /// \brief Main program of the errorpropagation example
0029 //
0030 // ------------------------------------------------------------
0031 //      GEANT 4 example main
0032 // ------------------------------------------------------------
0033 //
0034 // History:
0035 // - Created:   P. Arce   May 2007
0036 //
0037
0038 #include "ExErrorDetectorConstruction.hh"
0039
0040 #include "G4ErrorCylSurfaceTarget.hh"
0041 #include "G4ErrorFreeTrajState.hh"
0042 #include "G4ErrorGeomVolumeTarget.hh"
0043 #include "G4ErrorPlaneSurfaceTarget.hh"
0044 #include "G4ErrorPropagator.hh"
0045 #include "G4ErrorPropagatorData.hh"
0046 #include "G4ErrorPropagatorManager.hh"
0047 #include "G4ErrorTrackLengthTarget.hh"
0048 #include "G4SteppingVerbose.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4UImanager.hh"
0051
0052 void Initialize();
0053 G4ErrorTarget* BuildTarget(G4int iTarget);
0054 void ProcessEvent(G4int iProp, size_t nEv);
0055 void Finalize();
0056
0057 G4ErrorTarget* theTarget;
0058 G4ErrorMode theG4ErrorMode;
0059
0060 //-------------------------------------------------------------
0061 int main()
0062 {
0063   Initialize();
0064
0065   //----- Choose propagation mode
0066   // 0: propagate until target, all steps in one go
0067   // 1: propagate until target, returning control to the user at each step
0068   G4int iProp = 0;
0069   char* prop = std::getenv("G4ERROR_PROP");
0070   if (prop) {
0071     if (G4String(prop) == G4String("UNTIL_TARGET")) {
0072       iProp = 0;
0073     }
0074     else if (G4String(prop) == G4String("STEP_BY_STEP")) {
0075       iProp = 1;
0076     }
0077     else {
0078       G4Exception("exG4eReco", "Fatal error in Argument", FatalErrorInArgument,
0079                   G4String("Variable G4ERROR_PROP = " + G4String(prop)
0080                            + "   It must be: UNTIL_TARGET or STEP_BY_STEP")
0081                     .c_str());
0082     }
0083   }
0084   else {
0085     G4Exception("exG4eReco", "Fatal error in Argument", JustWarning,
0086                 "Variable G4ERROR_PROP not defined, taking it = UNTIL_TARGET");
0087   }
0088
0089   size_t nEvents = 3;
0090   for (size_t ii = 0; ii < nEvents; ii++) {
0091     ProcessEvent(iProp, ii);
0092   }
0093
0094   Finalize();
0095 }
0096
0097 //-------------------------------------------------------------
0098 void Initialize()
0099 {
0100   G4VSteppingVerbose::SetInstance(new G4SteppingVerbose);
0101
0102   // Initialize the GEANT4e manager
0103   G4ErrorPropagatorManager* g4emgr = G4ErrorPropagatorManager::GetErrorPropagatorManager();
0104   G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
0105
0106   g4emgr->SetUserInitialization(new ExErrorDetectorConstruction);
0107
0108   G4UImanager::GetUIpointer()->ApplyCommand("/exerror/setField -10. kilogauss");
0109
0110   g4emgr->InitGeant4e();
0111
0112   G4UImanager::GetUIpointer()->ApplyCommand("/control/verbose 1");
0113   G4UImanager::GetUIpointer()->ApplyCommand("/tracking/verbose 1");
0114   G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 100 mm");
0115
0116   //----- Choose target type:
0117   // 1: PlaneSurfaceTarget
0118   // 2: CylSurfaceTarget
0119   // 3: GeomVolumeTarget
0120   // 4: TrackLengthTarget
0121   G4int iTarget = 1;
0122   char* target = std::getenv("G4ERROR_TARGET");
0123   if (target) {
0124     if (G4String(target) == G4String("PLANE_SURFACE")) {
0125       iTarget = 1;
0126     }
0127     else if (G4String(target) == G4String("CYL_SURFACE")) {
0128       iTarget = 2;
0129     }
0130     else if (G4String(target) == G4String("VOLUME")) {
0131       iTarget = 3;
0132     }
0133     else if (G4String(target) == G4String("TRKLEN")) {
0134       iTarget = 4;
0135     }
0136     else {
0137       G4Exception("exG4eReco", "Fatal error in Argument", FatalErrorInArgument,
0138                   G4String("Variable G4ERROR_TARGET = " + G4String(target)
0139                            + "   It must be:  PLANE_SURFACE, CYL_SURFACE, VOLUME, TRKLEN")
0140                     .c_str());
0141     }
0142   }
0143   else {
0144     G4Exception("exG4eReco", "Fatal error in Argument", JustWarning,
0145                 "Variable G4ERROR_TARGET not defined, taking it = PLANE_SURFACE");
0146   }
0147
0148   theTarget = BuildTarget(iTarget);
0149   g4edata->SetTarget(theTarget);
0150
0151   theG4ErrorMode = G4ErrorMode_PropBackwards;
0152   char* mode = std::getenv("G4ERROR_MODE");
0153   if (mode) {
0154     if (G4String(mode) == G4String("FORWARDS")) {
0155       theG4ErrorMode = G4ErrorMode_PropForwards;
0156     }
0157     else if (G4String(mode) == G4String("BACKWARDS")) {
0158       theG4ErrorMode = G4ErrorMode_PropBackwards;
0159     }
0160     else {
0161       G4Exception("exG4eReco", "Fatal error in Argument", FatalErrorInArgument,
0162                   G4String("Variable G4ERROR_MODE = " + G4String(mode)
0163                            + "   It must be:  FORWARDS or BACKWARDS")
0164                     .c_str());
0165     }
0166   }
0167   else {
0168     G4Exception("exG4eReco", "Fatal error in Argument", JustWarning,
0169                 "Variable G4ERROR_MODE not defined, taking it = BACKWARDS");
0170   }
0171 }
0172
0173 void ProcessEvent(G4int iProp, size_t)
0174 {
0175   // Set the starting trajectory.
0176   G4ThreeVector xv3(0, 0, 0);
0177   G4ThreeVector pv3(20.0 * GeV, 0.0, 0.0);
0178   G4ErrorTrajErr error(5, 0);
0179   G4ErrorFreeTrajState* g4ErrorTrajState = new G4ErrorFreeTrajState("mu-", xv3, pv3, error);
0180
0181   G4ErrorPropagatorManager* g4emgr = G4ErrorPropagatorManager::GetErrorPropagatorManager();
0182
0183   // int ierr = 0;
0184
0185   G4Point3D surfPos(224. * cm, 0., 0.);
0186   G4Normal3D surfNorm(1., 0., 0.);
0187   //-  G4ErrorTarget* theG4ErrorTarget
0188   //     = new G4ErrorPlaneSurfaceTarget(surfNorm, surfPos );
0189
0190   if (iProp == 0) {
0191     // Propagate until G4ErrorTarget is found all in one go
0192     // ierr =
0193     g4emgr->Propagate(g4ErrorTrajState, theTarget, theG4ErrorMode);
0194   }
0195   else if (iProp == 1) {
0196     // Propagate until G4ErrorTarget is reached step by step
0197
0198     g4emgr->InitTrackPropagation();
0199
0200     //    G4Track* aTrack
0201     //      = G4EventManager::GetEventManager()->GetTrackingManager()->GetTrack();
0202     bool moreEvt = TRUE;
0203     while (moreEvt) {
0204       // ierr =
0205       g4emgr->PropagateOneStep(g4ErrorTrajState, theG4ErrorMode);
0206
0207       //---- Check if target is reached
0208       if (g4emgr->GetPropagator()->CheckIfLastStep(g4ErrorTrajState->GetG4Track())) {
0209         g4emgr->GetPropagator()->InvokePostUserTrackingAction(g4ErrorTrajState->GetG4Track());
0210         moreEvt = 0;
0211         G4cout << "STEP_BY_STEP propagation: Last Step " << G4endl;
0212       }
0213     }
0214   }
0215
0216   G4cout << " $$$ PROPAGATION ENDED " << G4endl;
0217   // extract current state info
0218   G4Point3D posEnd = g4ErrorTrajState->GetPosition();
0219   G4Normal3D momEnd = g4ErrorTrajState->GetMomentum();
0220   G4ErrorTrajErr errorEnd = g4ErrorTrajState->GetError();
0221
0222   G4cout << " Position: " << posEnd << G4endl << " Momentum: " << momEnd << G4endl
0223          << " Error: " << errorEnd << G4endl;
0224 }
0225
0226 //-------------------------------------------------------------
0227 G4ErrorTarget* BuildTarget(G4int iTarget)
0228 {
0229   G4ErrorTarget* target = 0;
0230   if (iTarget == 1) {
0231     G4Point3D surfPos(221. * cm, 0., 0.);
0232     G4Normal3D surfNorm(1., 0., 0.);
0233     target = new G4ErrorPlaneSurfaceTarget(surfNorm, surfPos);
0234   }
0235   else if (iTarget == 2) {
0236     G4double radius = 222 * cm;
0237     target = new G4ErrorCylSurfaceTarget(radius);
0238   }
0239   else if (iTarget == 3) {
0240     target = new G4ErrorGeomVolumeTarget("MUON");
0241   }
0242   else if (iTarget == 4) {
0243     target = new G4ErrorTrackLengthTarget(223. * cm);
0244   }
0245   else {
0246     G4Exception("exG4eReco", "Fatal error in Argument", FatalErrorInArgument,
0247                 "Target type has to be between 1 and 4");
0248   }
0249   return target;
0250 }
0251
0252 //-------------------------------------------------------------
0253 void Finalize()
0254 {
0255   G4ErrorPropagatorManager::GetErrorPropagatorManager()->CloseGeometry();
0256 }