BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ExtSteppingAction Class Reference

#include <ExtSteppingAction.h>

Inheritance diagram for ExtSteppingAction:

Public Member Functions

 ExtSteppingAction ()
 ~ExtSteppingAction ()
void Reset ()
void MucReset ()
void UserSteppingAction (const G4Step *currentStep)
void SetInitialPath (double aPath)
void SetInitialTof (double aTof)
void SetBetaInMDC (double aBeta)
void SetXpErrPointer (Ext_xp_err *xpErr)
void SetMsgFlag (bool aMsgFalg)
void SetMucKalFlag (bool aMucKalFlag)
void SetMucWindow (int aMucWindow)
void SetExtTrackPointer (RecExtTrack *aExtTrack)
void CalculateEmcEndThetaPhi (int npart, int sector, int nb, int &ntheta, int &nphi)
int CalculateEmcEndPhiNb (int num)
int CalculateEmcEndCopyNb (int num)
void Set_which_tof_version (int version)
int Get_which_tof_version (void)
void InfmodMuc (Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Hep3Vector GetGapID (G4String vol)
bool TrackStop ()
void SetMucDigiColPointer (MucDigiCol *rawdigicol)

Detailed Description

Definition at line 24 of file ExtSteppingAction.h.

Constructor & Destructor Documentation

◆ ExtSteppingAction()

ExtSteppingAction::ExtSteppingAction ( )

Definition at line 33 of file ExtSteppingAction.cxx.

34 : chicc( 0.0 )
35 , initialPath( 0. )
36 , myPathIntoCrystal( 0. )
37 , myPathOutCrystal( 0. )
38 , myPathInCrystal( 0. )
39 , myBetaInMDC( 0. )
40 , extXpErr( 0 )
41 , myExtTrack( 0 )
42 , msgFlag( true )
43 , myUseMucKalFlag( true )
44 , m_trackstop( false )
45 , myMucnfit_( 0 )
46 , myMucchisq_( 0. )
47 , myMucdepth_( -99. )
48 , myMucbrLastLay_( -1 )
49 , myMucecLastLay_( -1 )
50 , myMucnhits_( -1 )
51 , myMucWindow( 6 ) {
52 myTof1R = 814.001;
53 myTof1Z = 1330.0;
54 myTof2R = 871.036;
55 myEmcR1 = 945.0;
56 myEmcR2 = 983.9;
57 myEmcZ = 1416.8;
58 myMucR = 1700.0 - 0.01;
59 myMucZ = 2050.0 - 0.01;
60 m_which_tof_version = 2;
61}

◆ ~ExtSteppingAction()

ExtSteppingAction::~ExtSteppingAction ( )

Definition at line 63 of file ExtSteppingAction.cxx.

63{}

Member Function Documentation

◆ CalculateEmcEndCopyNb()

int ExtSteppingAction::CalculateEmcEndCopyNb ( int num)

Definition at line 1059 of file ExtSteppingAction.cxx.

1059 {
1060 int copyNb;
1061 switch ( num )
1062 {
1063 case 30: copyNb = 5; break;
1064 case 31: copyNb = 6; break;
1065 case 32: copyNb = 14; break;
1066 case 33: copyNb = 15; break;
1067 case 34: copyNb = 16; break;
1068 default: copyNb = num; break;
1069 }
1070 return copyNb;
1071}

Referenced by UserSteppingAction().

◆ CalculateEmcEndPhiNb()

int ExtSteppingAction::CalculateEmcEndPhiNb ( int num)

Definition at line 1052 of file ExtSteppingAction.cxx.

1052 {
1053 if ( num == 0 || num == 1 ) { return 15 - num * 8; }
1054 else if ( num == 2 || num == 3 ) { return 30 - num * 8; }
1055 else if ( num <= 9 ) { return 17 - num; }
1056 else { return 15 - num; }
1057}

Referenced by UserSteppingAction().

◆ CalculateEmcEndThetaPhi()

void ExtSteppingAction::CalculateEmcEndThetaPhi ( int npart,
int sector,
int nb,
int & ntheta,
int & nphi )

Definition at line 983 of file ExtSteppingAction.cxx.

984 {
985 if ( ( sector >= 0 ) && ( sector < 3 ) ) sector += 16;
986 // if((sector!=7)&&(sector!=15))
987 {
988 if ( ( nb >= 0 ) && ( nb < 4 ) )
989 {
990 ntheta = 0;
991 nphi = ( sector - 3 ) * 4 + nb;
992 }
993 else if ( ( nb >= 4 ) && ( nb < 8 ) )
994 {
995 ntheta = 1;
996 nphi = ( sector - 3 ) * 4 + nb - 4;
997 }
998 else if ( ( nb >= 8 ) && ( nb < 13 ) )
999 {
1000 ntheta = 2;
1001 nphi = ( sector - 3 ) * 5 + nb - 8;
1002 }
1003 else if ( ( nb >= 13 ) && ( nb < 18 ) )
1004 {
1005 ntheta = 3;
1006 nphi = ( sector - 3 ) * 5 + nb - 13;
1007 }
1008 else if ( ( nb >= 18 ) && ( nb < 24 ) )
1009 {
1010 ntheta = 4;
1011 nphi = ( sector - 3 ) * 6 + nb - 18;
1012 }
1013 else if ( ( nb >= 24 ) && ( nb < 30 ) )
1014 {
1015 ntheta = 5;
1016 nphi = ( sector - 3 ) * 6 + nb - 24;
1017 }
1018 }
1019
1020 if ( npart == 2 )
1021 {
1022 switch ( ntheta )
1023 {
1024 case 0:
1025 if ( nphi < 32 ) nphi = 31 - nphi;
1026 else nphi = 95 - nphi;
1027 break;
1028 case 1:
1029 if ( nphi < 32 ) nphi = 31 - nphi;
1030 else nphi = 95 - nphi;
1031 break;
1032 case 2:
1033 if ( nphi < 40 ) nphi = 39 - nphi;
1034 else nphi = 119 - nphi;
1035 break;
1036 case 3:
1037 if ( nphi < 40 ) nphi = 39 - nphi;
1038 else nphi = 119 - nphi;
1039 break;
1040 case 4:
1041 if ( nphi < 48 ) nphi = 47 - nphi;
1042 else nphi = 143 - nphi;
1043 break;
1044 case 5:
1045 if ( nphi < 48 ) nphi = 47 - nphi;
1046 else nphi = 143 - nphi;
1047 default:;
1048 }
1049 }
1050}

Referenced by UserSteppingAction().

◆ Get_which_tof_version()

int ExtSteppingAction::Get_which_tof_version ( void )
inline

Definition at line 49 of file ExtSteppingAction.h.

49{ return m_which_tof_version; }

◆ GetGapID()

Hep3Vector ExtSteppingAction::GetGapID ( G4String vol)

Definition at line 1079 of file ExtSteppingAction.cxx.

1079 {
1080 int par( -1 ), se( -1 ), ga( -1 );
1081 G4String strPart = vol.substr( 5, 1 );
1082 // G4String strSeg = m_MotherVolumeName.substr(7,1);
1083
1084 G4String strSeg = vol.substr( 7, 1 );
1085 G4String strGap;
1086 if ( vol.contains( "G" ) ) strGap = vol.substr( 9, 1 );
1087 if ( vol.contains( "Ab" ) ) strGap = vol.substr( 10, 1 );
1088 std::istrstream partBuf( strPart.c_str(), strlen( strPart.c_str() ) );
1089 std::istrstream segBuf( strSeg.c_str(), strlen( strSeg.c_str() ) );
1090 std::istrstream gapBuf( strGap.c_str(), strlen( strGap.c_str() ) );
1091
1092 partBuf >> par;
1093 segBuf >> se;
1094 gapBuf >> ga;
1095 if ( vol.contains( "Ab" ) && par == 1 ) ga = ga + 1;
1096 // panelBuf >> m_panel;
1097 Hep3Vector vec;
1098 vec[0] = par;
1099 vec[1] = se;
1100 vec[2] = ga;
1101 return vec;
1102}
vec
dble_vec_t vec[12]
Definition BesEvtGen/phokhara/PHOKHARA/ranlxd.c:330

Referenced by UserSteppingAction().

◆ InfmodMuc()

void ExtSteppingAction::InfmodMuc ( Hep3Vector & pos,
Hep3Vector & mom,
HepSymMatrix & err )

Definition at line 1073 of file ExtSteppingAction.cxx.

1073 {
1074 pos = m_pos_mod;
1075 mom = m_mom_mod;
1076 err = m_err_mod;
1077}

Referenced by TrkExtAlg::execute().

◆ MucReset()

void ExtSteppingAction::MucReset ( )

Definition at line 105 of file ExtSteppingAction.cxx.

105 {
106 RemStep = 0;
107 RemDist = 99999.;
108 RemDepth = 0.;
109 RemID[0] = -1;
110 RemID[1] = -1;
111 RemID[2] = -1;
112 RemVol = "";
113}

◆ Reset()

void ExtSteppingAction::Reset ( )

Definition at line 65 of file ExtSteppingAction.cxx.

65 {
66 chicc = 0.;
67 myExtTrack = 0;
68 myPathIntoCrystal = 0.;
69 myPathOutCrystal = 0.;
70 myPathInCrystal = 0.;
71
72 myPathIntoTof1 = 0.0;
73 myPathOutTof1 = 0.0;
74 // myPathInTof1 = 0.0;
75 myPathInTof1.clear();
76
77 myPathIntoTof2 = 0.0;
78 myPathOutTof2 = 0.0;
79 // myPathInTof2 = 0.0;
80 myPathInTof2.clear();
81
82 myTofFlag = false;
83 myTof1Flag = false;
84 myTof2Flag = false;
85 myInTof1 = false;
86 myInTof2 = false;
87 myOutTof1 = false;
88 myOutTof2 = false;
89 myPhyEmcFlag = false;
90 myEmcFlag = false;
91 myEmcPathFlag = false;
92 myMucFlag = false;
93
94 m_trackstop = false;
95 myMucchisq_ = 0.;
96 myMucnfit_ = 0;
97 myMucdepth_ = -99.;
98 myMucbrLastLay_ = -1;
99 myMucecLastLay_ = -1;
100 RememberID[0] = -1;
101 RememberID[1] = -1;
102 RememberID[2] = -1;
103}

Referenced by TrkExtAlg::execute().

◆ Set_which_tof_version()

void ExtSteppingAction::Set_which_tof_version ( int version)
inline

Definition at line 48 of file ExtSteppingAction.h.

48{ m_which_tof_version = version; }

Referenced by TrkExtAlg::execute().

◆ SetBetaInMDC()

void ExtSteppingAction::SetBetaInMDC ( double aBeta)
inline

Definition at line 36 of file ExtSteppingAction.h.

36{ myBetaInMDC = aBeta; };

◆ SetExtTrackPointer()

void ExtSteppingAction::SetExtTrackPointer ( RecExtTrack * aExtTrack)
inline

Definition at line 42 of file ExtSteppingAction.h.

42{ myExtTrack = aExtTrack; };

Referenced by TrkExtAlg::execute().

◆ SetInitialPath()

void ExtSteppingAction::SetInitialPath ( double aPath)
inline

Definition at line 34 of file ExtSteppingAction.h.

34{ initialPath = aPath; };

◆ SetInitialTof()

void ExtSteppingAction::SetInitialTof ( double aTof)
inline

Definition at line 35 of file ExtSteppingAction.h.

35{ initialTof = aTof; };

◆ SetMsgFlag()

void ExtSteppingAction::SetMsgFlag ( bool aMsgFalg)
inline

Definition at line 39 of file ExtSteppingAction.h.

39{ msgFlag = aMsgFalg; };

Referenced by Ext_track::Initialization().

◆ SetMucDigiColPointer()

void ExtSteppingAction::SetMucDigiColPointer ( MucDigiCol * rawdigicol)
inline

Definition at line 54 of file ExtSteppingAction.h.

54{ m_mucdigicol = rawdigicol; }

Referenced by TrkExtAlg::execute().

◆ SetMucKalFlag()

void ExtSteppingAction::SetMucKalFlag ( bool aMucKalFlag)
inline

Definition at line 40 of file ExtSteppingAction.h.

40{ myUseMucKalFlag = aMucKalFlag; };

◆ SetMucWindow()

void ExtSteppingAction::SetMucWindow ( int aMucWindow)
inline

Definition at line 41 of file ExtSteppingAction.h.

41{ myMucWindow = aMucWindow; };

◆ SetXpErrPointer()

void ExtSteppingAction::SetXpErrPointer ( Ext_xp_err * xpErr)
inline

Definition at line 37 of file ExtSteppingAction.h.

37{ extXpErr = xpErr; };

◆ TrackStop()

bool ExtSteppingAction::TrackStop ( )
inline

Definition at line 53 of file ExtSteppingAction.h.

53{ return m_trackstop; }

Referenced by TrkExtAlg::execute().

◆ UserSteppingAction()

void ExtSteppingAction::UserSteppingAction ( const G4Step * currentStep)

Definition at line 115 of file ExtSteppingAction.cxx.

115 {
116
117 // Get track and its status
118 G4Track* currentTrack = currentStep->GetTrack();
119 if ( !currentTrack )
120 {
121 cout << "Can't get currentTrack" << endl;
122 return;
123 }
124 G4TrackStatus currentTrackStatus = currentTrack->GetTrackStatus();
125
126 int stepNumber = currentTrack->GetCurrentStepNumber();
127 if ( msgFlag ) cout << "STEP " << stepNumber << ":" << endl;
128
129 // Get current position and momentum
130 Hep3Vector currentPosition = currentTrack->GetPosition();
131 Hep3Vector currentMomentum = currentTrack->GetMomentum();
132
133 // Get current Volume
134 G4VPhysicalVolume* oldVolume = nullptr;
135 G4VPhysicalVolume* newVolume = nullptr;
136 if ( !currentTrack->GetTouchableHandle() )
137 cout << "Can't get currentTrack->GetTouchableHandle()" << endl;
138 else oldVolume = currentTrack->GetTouchableHandle()->GetVolume();
139 if ( !currentTrack->GetNextTouchableHandle() )
140 cout << "Can't get currentTrack->GetNextTouchableHandle()" << endl;
141 else newVolume = currentTrack->GetNextTouchableHandle()->GetVolume();
142 if ( !oldVolume ) cout << "Can't get oldVolume!" << endl;
143
144 //****added by LI Chunhua
145 if ( stepNumber > 50000 )
146 {
147 cout << "infinite steps in the track " << endl;
148 currentTrack->SetTrackStatus( fStopAndKill );
149 m_trackstop = true;
150 }
151
152 G4String ParticleName = currentTrack->GetDefinition()->GetParticleName();
153 double x_ = currentPosition.x();
154 double y_ = currentPosition.y();
155 double z_ = currentPosition.z();
156 bool inner =
157 ( oldVolume != newVolume ) &&
158 ( !( ( fabs( x_ ) >= myMucR ) || ( fabs( y_ ) >= myMucR ) ||
159 ( ( fabs( x_ - y_ ) / sqrt( 2. ) ) >= myMucR ) ||
160 ( ( fabs( x_ + y_ ) / sqrt( 2. ) ) >= myMucR ) || ( fabs( z_ ) >= myMucZ ) ) );
161 bool mucdec = ( fabs( x_ ) >= myMucR ) || ( fabs( y_ ) >= myMucR ) ||
162 ( ( fabs( x_ - y_ ) / sqrt( 2. ) ) >= myMucR ) ||
163 ( ( fabs( x_ + y_ ) / sqrt( 2. ) ) >= myMucR ) || ( fabs( z_ ) >= myMucZ );
164
165 // if current particle is alive.
166 if ( currentTrackStatus == fAlive )
167 {
168 // if(oldVolume!=newVolume)//if so, update Error Matrix and chicc.
169 if ( inner || mucdec ) // update Error Matrix for newvolume in MDC,TOF,EMC; update Error
170 // Matrix step by step in MUC;
171 {
172 // Get current B field
173 double currentPoint[3] = { currentPosition.x(), currentPosition.y(),
174 currentPosition.z() };
175 double currentBfield[3] = { 0.0, 0.0, 0.0 };
176 Hep3Vector currentB( 0.0, 0.0, 1.0 );
177 if ( G4TransportationManager::GetTransportationManager() )
178 {
179 G4FieldManager* fieldManager =
180 G4TransportationManager::GetTransportationManager()->GetFieldManager();
181 if ( fieldManager )
182 {
183 if ( fieldManager->GetDetectorField() )
184 {
185 fieldManager->GetDetectorField()->GetFieldValue( currentPoint, currentBfield );
186 currentB.set( currentBfield[0] / tesla, currentBfield[1] / tesla,
187 currentBfield[2] / tesla );
188 if ( msgFlag )
189 cout << "B:" << currentB.x() << "," << currentB.y() << "," << currentB.z()
190 << endl;
191 }
192 }
193 }
194
195 // Get chicc
196 G4Material* oldMaterial = oldVolume->GetLogicalVolume()->GetMaterial();
197 if ( !oldMaterial ) std::cout << "Can't get oldMaterial" << std::endl;
198 else CalculateChicc( oldMaterial );
199
200 // Update Error Matrix
201 if ( !( extXpErr->move( currentPosition, currentMomentum, currentB, 1, chicc ) ) )
202 if ( msgFlag ) cout << "can not update Error Matrix!!" << endl;
203
204 // Verbose
205 if ( msgFlag )
206 {
207 cout << "Volume name:" << newVolume->GetName() << endl;
208 cout << "Volume number:" << newVolume->GetCopyNo() << endl;
209 cout << "x:" << currentPoint[0] << "//y:" << currentPoint[1]
210 << "//z:" << currentPoint[2] << "||px:" << currentMomentum.x()
211 << "//py:" << currentMomentum.y() << "//pz:" << currentMomentum.z() << endl;
212 cout << "ptot = " << currentMomentum.mag() << endl;
213 cout << "trk len = " << currentTrack->GetTrackLength() << endl;
214 cout << "Error matrix is:" << extXpErr->get_err() << endl;
215 cout << "phi:" << atan( currentPoint[1] / currentPoint[0] ) << endl;
216 Hep3Vector nz( 0., 0., 1. );
217 cout << "Projected z error:" << extXpErr->get_plane_err( currentMomentum.unit(), nz )
218 << endl;
219 double x, y, R;
220 x = currentPoint[0];
221 y = currentPoint[1];
222 R = sqrt( x * x + y * y );
223 Hep3Vector nt( -y / R, x / R, 0. );
224 cout << "Projected phi error:"
225 << ( extXpErr->get_plane_err( currentMomentum.unit(), nt ) ) / R << endl
226 << endl;
227 }
228
229 // some often used things
230 Hep3Vector xVector( 1.0, 0, 0 );
231 Hep3Vector yVector( 0, 1.0, 0 );
232 Hep3Vector zVector( 0, 0, 1.0 );
233 Hep3Vector tzVector;
234 tzVector.setRThetaPhi( 1.0, M_PI / 2.0, currentPosition.phi() );
235 double r = currentPosition.perp();
236 double x = currentPosition.x();
237 double y = currentPosition.y();
238 double z = currentPosition.z();
239 G4String newVolumeName = newVolume->GetName();
240 G4String oldVolumeName = oldVolume->GetName();
241 G4StepPoint* postStepPoint = currentStep->GetPostStepPoint();
242 G4TouchableHistory* theTouchable =
243 (G4TouchableHistory*)( postStepPoint->GetTouchable() );
244 int newVolumeNumber = newVolume->GetCopyNo();
245 // int newVolumeNumber=theTouchable->GetReplicaNumber(2);
246
247 G4String name1;
248 if ( newVolumeName == "logical_gasLayer" )
249 {
250 name1 = theTouchable->GetVolume( 3 )->GetName();
251 newVolumeNumber = theTouchable->GetReplicaNumber( 3 );
252 }
253
254 // Get PhyTof data
255
256 // std::cout << "ExtSteppingAction NewVolumeName: " << newVolumeName <<
257 // std::endl; std::cout << "ExtSteppingAction OldVolumeName: " << oldVolumeName
258 // << std::endl; std::cout << "ExtSteppingAction name1: " <<name1<<std::endl;
259 // Comment this part to remove the airtight Tof. Most hits would be replaced by the
260 // following specific sensitive detector. The remianing few hits are abondoned in Tof
261 // reconstruction.
262 /*
263 if( (!myTofFlag) && (!myTof1Flag) && (newVolumeName.contains("Tof") ))
264 {
265 newVolumeNumber = -2;
266 double currentTrackLength = currentTrack->GetTrackLength();
267 double totalTrackLength = currentTrackLength + initialPath;
268 //double initialTof = initialPath/(myBetaInMDC*299.792458);
269 //cout<<"initialTof = "<<initialTof<<endl;
270 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
271 double localTime = currentTrack->GetLocalTime();
272 //cout<<"localTime = "<<localTime<<endl;
273 double totalTOF = localTime + initialTof;
274 //cout<<"totalTOF= "<<totalTOF<<endl;
275
276 //std::cout << "ExtSteppingAction Volumename contains one of the marker words" <<
277 std::endl;
278
279 if(myExtTrack)
280 {
281 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
282 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
283 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
284 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
285 myExtTrack->SetTof1Data(currentPosition/10.0,currentMomentum/1000.0,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
286 myTofFlag = true;
287 }
288 return;
289 }//close if (!myTofFlag)
290 */
291
292 // Get Tof layer1 Ext data
293 if ( ( !myTof1Flag ) &&
294 ( newVolumeName == "logicalScinBr1" || newVolumeName.contains( "ScinEc" ) ||
295 ( newVolumeName == "logical_gasLayer" &&
296 ( name1 == "logical_container_m0" || name1 == "logical_container_m1" ) ) ) )
297 {
298 double currentTrackLength = currentTrack->GetTrackLength();
299 double totalTrackLength = currentTrackLength + initialPath;
300 double localTime = currentTrack->GetLocalTime();
301 double totalTOF = localTime + initialTof;
302 myInTof1 = true;
303 myPathIntoTof1 = currentTrackLength;
304 if ( msgFlag ) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
305
306 newVolumeNumber = theTouchable->GetReplicaNumber( 2 );
307
308 if ( newVolumeName.contains( "ScinEc" ) )
309 {
310 newVolumeNumber = ( 95 - newVolumeNumber ) / 2;
311 newVolumeNumber = newVolumeNumber + 176;
312
313 if ( newVolumeName.contains( "East" ) ) newVolumeNumber = newVolumeNumber + 48;
314 if ( myExtTrack )
315 {
316 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
317 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
318 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
319 double tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector );
320 myExtTrack->SetTof1Data(
321 currentPosition / 10, currentMomentum / 1000, newVolumeName, newVolumeNumber,
322 totalTOF, totalTrackLength / 10, myOutputSymMatrix( extXpErr->get_err() ),
323 zError / 10, tzError / 10, xError / 10, yError / 10 );
324 myTof1Flag = true;
325 }
326 }
327
328 else if ( newVolumeName == "logical_gasLayer" )
329 {
330 G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
331 G4ThreeVector localPosition =
332 theTouchable->GetHistory()->GetTopTransform().TransformPoint(
333 currentPostPosition );
334 // cout<<"localPosition: "<<localPosition.x()<<" "<<localPosition.y()<<"
335 // "<<localPosition.z()<<endl;
336
337 // Here we assume the spread of one strip is 24+3=27mm.
338 // Distance between strips and lower large glass margin: 6, top: 8,
339 // z_strip=z_small_glass-0.5
340 double z_mm = localPosition.z() - 0.5 + ( 24 + 3 ) * 6;
341 int strip;
342 if ( z_mm <= 0 ) { strip = 0; }
343 else if ( z_mm > 0 && z_mm < 12 * 27 ) { strip = floor( z_mm / 27 ); }
344 else { strip = 11; }
345 if ( strip < 0 ) strip = 0;
346 if ( strip > 11 ) strip = 11;
347
348 newVolumeNumber = theTouchable->GetReplicaNumber( 3 );
349 newVolumeNumber = 35 - newVolumeNumber;
350 if ( currentPosition.z() < 0.0 ) { newVolumeNumber = newVolumeNumber + 36; }
351 newVolumeNumber = 12 * newVolumeNumber + strip;
352 newVolumeNumber = newVolumeNumber + 176 + 96;
353 if ( myExtTrack )
354 {
355 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
356 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
357 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
358 double tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector );
359 Hep3Vector locP =
360 Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
361 // cout<<"locP: "<<locP.x()<<" "<<locP.y()<<" "<<locP.z()
362 //<<" currentMomentum: "<<currentMomentum.x()<<" "<<currentMomentum.y()<<"
363 //"<<currentMomentum.z()
364 //<<" newVolumeName: "<<newVolumeName<<" newVolumeNumber: "<<newVolumeNumber
365 //<<" totalTOF: "<<totalTOF<<" totalTrackLength: "<<totalTrackLength<<endl;
366 myExtTrack->SetTof1Data( locP / 10, currentMomentum / 1000, newVolumeName,
367 newVolumeNumber, totalTOF, totalTrackLength / 10,
368 myOutputSymMatrix( extXpErr->get_err() ), zError / 10,
369 tzError / 10, xError / 10, yError / 10 );
370 myTof1Flag = true;
371 }
372 }
373 else
374 {
375 newVolumeNumber = ( 527 - newVolumeNumber ) / 3;
376 if ( myExtTrack )
377 {
378 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
379 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
380 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
381 double tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector );
382 myExtTrack->SetTof1Data(
383 currentPosition / 10, currentMomentum / 1000, newVolumeName, newVolumeNumber,
384 totalTOF, totalTrackLength / 10, myOutputSymMatrix( extXpErr->get_err() ),
385 zError / 10, tzError / 10, xError / 10, yError / 10 );
386 myTof1Flag = true;
387 }
388 }
389
390 return;
391
392 } // close if (!myTof1Flag)
393
394 if ( myInTof1 &&
395 ( oldVolumeName == "logicalScinBr1" || oldVolumeName.contains( "ScinEc" ) ||
396 ( oldVolumeName == "logical_gasLayer" &&
397 ( name1 == "logical_container_m0" || name1 == "logical_container_m1" ) ) ) )
398 {
399 myInTof1 = false;
400 myOutTof1 = true;
401 myPathOutTof1 = currentTrack->GetTrackLength();
402 myPathInTof1.push_back( myPathOutTof1 - myPathIntoTof1 );
403 if ( msgFlag ) cout << "myPathOutTof1 = " << myPathOutTof1 << endl;
404 }
405
406 if ( myOutTof1 &&
407 ( newVolumeName == "logicalScinBr1" || newVolumeName.contains( "ScinEc" ) ||
408 ( newVolumeName == "logical_gasLayer" &&
409 ( name1 == "logical_container_m0" || name1 == "logical_container_m1" ) ) ) )
410 {
411 myInTof1 = true;
412 myOutTof1 = false;
413 myPathIntoTof1 = currentTrack->GetTrackLength();
414 if ( msgFlag ) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
415 }
416
417 // Get Tof layer2 Ext data
418 if ( ( !myTof2Flag ) &&
419 ( newVolumeName == "logicalScinBr2" ||
420 ( newVolumeName == "logical_gasLayer" &&
421 ( name1 == "logical_container_m2" || name1 == "logical_container_m3" ) ) ) )
422 {
423 double currentTrackLength = currentTrack->GetTrackLength();
424 double totalTrackLength = currentTrackLength + initialPath;
425 // double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
426 double localTime = currentTrack->GetLocalTime();
427 double totalTOF = localTime + initialTof;
428 myInTof2 = true;
429 myPathIntoTof2 = currentTrackLength;
430 if ( msgFlag ) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
431 newVolumeNumber = theTouchable->GetReplicaNumber( 2 );
432
433 if ( newVolumeName == "logical_gasLayer" )
434 {
435 G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
436 G4ThreeVector localPosition =
437 theTouchable->GetHistory()->GetTopTransform().TransformPoint(
438 currentPostPosition );
439 // cout<<"localPosition: "<<localPosition.x()<<" "<<localPosition.y()<<"
440 // "<<localPosition.z()<<endl;
441
442 // Here we assume the spread of one strip is 24+3=27mm.
443 // Distance between strips and lower large glass margin: 6, top: 8,
444 // z_strip=z_small_glass-0.5
445 double z_mm = localPosition.z() - 0.5 + ( 24 + 3 ) * 6;
446 int strip;
447 if ( z_mm <= 0 ) { strip = 0; }
448 else if ( z_mm > 0 && z_mm < 12 * 27 ) { strip = floor( z_mm / 27 ); }
449 else { strip = 11; }
450 if ( strip < 0 ) strip = 0;
451 if ( strip > 11 ) strip = 11;
452
453 newVolumeNumber = theTouchable->GetReplicaNumber( 3 );
454 newVolumeNumber = 35 - newVolumeNumber;
455 if ( currentPosition.z() < 0.0 ) { newVolumeNumber = newVolumeNumber + 36; }
456 newVolumeNumber = 12 * newVolumeNumber + strip;
457 newVolumeNumber = newVolumeNumber + 176 + 96;
458
459 if ( myExtTrack )
460 {
461 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
462 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
463 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
464 double tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector );
465 Hep3Vector locP =
466 Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
467 // cout<<"locP: "<<locP.x()<<" "<<locP.y()<<" "<<locP.z()
468 //<<" currentMomentum: "<<currentMomentum.x()<<" "<<currentMomentum.y()<<"
469 //"<<currentMomentum.z()
470 //<<" newVolumeName: "<<newVolumeName<<" newVolumeNumber: "<<newVolumeNumber
471 //<<" totalTOF: "<<totalTOF<<" totalTrackLength: "<<totalTrackLength<<endl;
472 myExtTrack->SetTof2Data( locP / 10, currentMomentum / 1000, newVolumeName,
473 newVolumeNumber, totalTOF, totalTrackLength / 10,
474 myOutputSymMatrix( extXpErr->get_err() ), zError / 10,
475 tzError / 10, xError / 10, yError / 10 );
476 myTof2Flag = true;
477 }
478 }
479 else
480 {
481 newVolumeNumber = ( 527 - newVolumeNumber ) / 3;
482 if ( myExtTrack )
483 {
484 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
485 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
486 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
487 double tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector );
488 myExtTrack->SetTof2Data(
489 currentPosition / 10, currentMomentum / 1000, newVolumeName, newVolumeNumber,
490 totalTOF, totalTrackLength / 10, myOutputSymMatrix( extXpErr->get_err() ),
491 zError / 10, tzError / 10, xError / 10, yError / 10 );
492 myTof2Flag = true;
493 }
494 }
495 return;
496 } // close if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" )
497
498 if ( myInTof2 &&
499 ( oldVolumeName == "logicalScinBr2" ||
500 ( oldVolumeName == "logical_gasLayer" &&
501 ( name1 == "logical_container_m2" || name1 == "logical_container_m3" ) ) ) )
502 {
503 myInTof2 = false;
504 myOutTof2 = true;
505 myPathOutTof2 = currentTrack->GetTrackLength();
506 myPathInTof2.push_back( myPathOutTof2 - myPathIntoTof2 );
507 if ( msgFlag ) cout << "myPathOutTof2 = " << myPathOutTof2 << endl;
508 }
509
510 if ( myOutTof2 &&
511 ( newVolumeName == "logicalScinBr2" ||
512 ( newVolumeName == "logical_gasLayer" &&
513 ( name1 == "logical_container_m2" || name1 == "logical_container_m3" ) ) ) )
514 {
515 myInTof2 = true;
516 myOutTof2 = false;
517 myPathIntoTof2 = currentTrack->GetTrackLength();
518 if ( msgFlag ) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
519 }
520
521 // Get PhyEmc data.
522 if ( ( !myPhyEmcFlag ) && ( !myEmcFlag ) &&
523 ( newVolumeName == "EMC" || newVolumeName.contains( "BSC" ) ||
524 newVolumeName == "EndPhi" ) )
525 {
526 newVolumeNumber = -2;
527 if ( myExtTrack )
528 {
529 // myPathIntoCrystal = currentTrack->GetTrackLength();
530 Hep3Vector nPhi( -y / r, x / r, 0. );
531 double errorPhi = ( extXpErr->get_plane_err( currentMomentum.unit(), nPhi ) ) / r;
532
533 Hep3Vector aPosition = currentPosition;
534 double R = aPosition.r();
535 double aTheta = aPosition.theta();
536 aPosition.setTheta( aTheta + M_PI_2 );
537 double errorTheta;
538 errorTheta =
539 ( extXpErr->get_plane_err( currentMomentum.unit(), aPosition.unit() ) ) / R;
540 if ( msgFlag ) cout << "Theta direction: " << aPosition << endl;
541 myExtTrack->SetEmcData( currentPosition / 10, currentMomentum / 1000, newVolumeName,
542 newVolumeNumber, errorTheta, errorPhi,
543 myOutputSymMatrix( extXpErr->get_err() ) );
544 }
545 myPhyEmcFlag = true;
546 return;
547 }
548
549 if ( !myEmcPathFlag && newVolumeName.contains( "Crystal" ) )
550 {
551 myPathIntoCrystal = currentTrack->GetTrackLength();
552 // cout<<"Enter Crystal, path: "<<myPathIntoCrystal<<endl;
553 myEmcPathFlag = true;
554 }
555
556 // Get Emc Ext data.
557 if ( ( !myEmcFlag ) && newVolumeName.contains( "Crystal" ) )
558 {
559 if ( myExtTrack )
560 {
561 //------------------- record crystal number
562 int npart, ntheta, nphi;
563 if ( currentTrack->GetNextTouchableHandle()->GetVolume( 1 )->GetName().contains(
564 "BSC" ) )
565 { // barrel
566 npart = 1;
567 std::istringstream thetaBuf(
568 ( currentTrack->GetNextTouchableHandle()->GetVolume( 1 )->GetName() )
569 .substr( 16, 2 ) );
570 thetaBuf >> ntheta;
571 nphi = 308 - currentTrack->GetNextTouchableHandle()->GetCopyNumber( 2 );
572 }
573 else
574 { // endcap
575 npart = 2 - 2 * currentTrack->GetNextTouchableHandle()->GetCopyNumber( 3 );
576 int endSector = currentTrack->GetNextTouchableHandle()->GetCopyNumber( 2 );
577 int endNb, endNbGDML;
578 std::istringstream thetaBuf(
579 ( currentTrack->GetNextTouchableHandle()->GetVolume( 0 )->GetName() )
580 .substr( 20, 2 ) );
581 thetaBuf >> endNb;
582 endNbGDML = CalculateEmcEndCopyNb( endNb );
583 int endSectorGDML = CalculateEmcEndPhiNb( endSector );
584 CalculateEmcEndThetaPhi( npart, endSectorGDML, endNbGDML, ntheta, nphi );
585 }
586 ostringstream strEmc;
587 if ( ntheta < 10 )
588 { strEmc << "EmcPart" << npart << "Theta0" << ntheta << "Phi" << nphi; }
589 else
590 {
591 strEmc << "EmcPart" << npart << "Theta" << ntheta << "Phi" << nphi;
592 } //------------------------------------------
593
594 // myPathIntoCrystal = currentTrack->GetTrackLength();
595 Hep3Vector nPhi( -y / r, x / r, 0. );
596 double errorPhi = ( extXpErr->get_plane_err( currentMomentum.unit(), nPhi ) ) / r;
597
598 Hep3Vector aPosition = currentPosition;
599 double R = aPosition.r();
600 double aTheta = aPosition.theta();
601 aPosition.setTheta( aTheta + M_PI_2 );
602 double errorTheta;
603 errorTheta =
604 ( extXpErr->get_plane_err( currentMomentum.unit(), aPosition.unit() ) ) / R;
605 if ( msgFlag ) cout << "Theta direction: " << aPosition << endl;
606 myExtTrack->SetEmcData( currentPosition / 10, currentMomentum / 1000, strEmc.str(),
607 newVolumeNumber, errorTheta, errorPhi,
608 myOutputSymMatrix( extXpErr->get_err() ) );
609 }
610 myEmcFlag = true;
611 return;
612 }
613
614 // Get path in Emc
615 if ( myEmcPathFlag && oldVolumeName.contains( "Crystal" ) )
616 {
617 myPathOutCrystal = currentTrack->GetTrackLength();
618 myPathInCrystal = myPathInCrystal + myPathOutCrystal - myPathIntoCrystal;
619 // cout<<"Go out of crystal, path: "<<myPathOutCrystal<<endl;
620 myEmcPathFlag = false;
621 }
622
623 // Get Muc Ext Data.
624 if ( ( !myMucFlag ) &&
625 ( ( fabs( x ) >= myMucR ) || ( fabs( y ) >= myMucR ) ||
626 ( ( fabs( x - y ) / sqrt( 2. ) ) >= myMucR ) ||
627 ( ( fabs( x + y ) / sqrt( 2. ) ) >= myMucR ) || ( fabs( z ) >= myMucZ ) ) )
628 {
629 int newVolumeNumber = newVolume->GetCopyNo();
630 if ( myExtTrack )
631 {
632 Hep3Vector xVector( 1.0, 0, 0 );
633 Hep3Vector yVector( 0, 1.0, 0 );
634 Hep3Vector zVector( 0, 0, 1.0 );
635 double xError = extXpErr->get_plane_err( currentMomentum.unit(), xVector );
636 double yError = extXpErr->get_plane_err( currentMomentum.unit(), yVector );
637 double zError = extXpErr->get_plane_err( currentMomentum.unit(), zVector );
638 double tzError = 0.0;
639 double phi = currentPosition.phi();
640 if ( phi < 0. ) phi += M_PI;
641 int i = int( 8.0 * phi / M_PI );
642 if ( i == 0 || i == 7 || i == 8 ) { tzError = yError; }
643 if ( i == 1 || i == 2 )
644 {
645 Hep3Vector tzVector( -1.0, 1.0, 0. );
646 tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector.unit() );
647 }
648 if ( i == 3 || i == 4 ) { tzError = xError; }
649 if ( i == 5 || i == 6 )
650 {
651 Hep3Vector tzVector( 1.0, 1.0, 0. );
652 tzError = extXpErr->get_plane_err( currentMomentum.unit(), tzVector.unit() );
653 }
654 myExtTrack->SetMucData( currentPosition / 10, currentMomentum / 1000, newVolumeName,
655 newVolumeNumber, myOutputSymMatrix( extXpErr->get_err() ),
656 zError / 10, tzError / 10, xError / 10, yError / 10 );
657 }
658 myMucFlag = true;
659 if ( !( ParticleName.contains( "mu" ) && myUseMucKalFlag ) )
660 {
661 // currentTrack->SetKineticEnergy(0.0);
662 currentStep->GetTrack()->SetTrackStatus( fStopAndKill );
663 m_trackstop = true;
664 return;
665 }
666 }
667
668 //**************************************************
669 //************inset KalFilter Algorithm in MUC******
670 //**************************************************
671 HepSymMatrix XpErr = extXpErr->get_err();
672 int namesize = oldVolumeName.size();
673 bool Flag1( false ), Flag2( false ), Flag3( false ), Flag4( false ), Flag5( false );
674 Flag1 = m_mucdigicol->size() > 0;
675 Flag2 = myUseMucKalFlag;
676 Flag3 = ParticleName.contains( "mu" ) && oldVolumeName.contains( "lMuc" ) &&
677 oldVolumeName.contains( "P" ) && oldVolumeName.contains( "S" ) &&
678 ( oldVolumeName.contains( "G" ) || oldVolumeName.contains( "Ab" ) );
679 Flag4 = oldVolumeName.contains( "lMuc" ) && oldVolumeName.contains( "P" ) &&
680 oldVolumeName.contains( "S" ) &&
681 ( ( namesize == 10 && oldVolumeName.contains( "G" ) ) ||
682 ( namesize == 11 && oldVolumeName.contains( "Ab" ) ) );
683 Flag5 = !( ( RemID[0] == 1 && RemID[2] == 9 ) ||
684 ( ( RemID[0] == 0 || RemID[0] == 2 ) && RemID[2] == 8 ) );
685 // if(!Flag5) {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
686 if ( Flag1 && Flag2 && Flag3 && Flag5 )
687 {
688 // get depth in Ab
689 double depth = currentStep->GetStepLength();
690 if ( oldVolumeName.contains( "Ab" ) ) RemDepth = RemDepth + depth;
691 if ( RemStep == 0 && Flag4 )
692 {
693 Hep3Vector ID_1 = GetGapID( oldVolumeName );
694 RemID = ID_1;
695
696 bool LastLay = ( ID_1[0] == 1 && ID_1[2] < 9 ) ||
697 ( ( ID_1[0] == 0 || ID_1[0] == 2 ) && ID_1[2] < 8 );
698 if ( RememberID[2] != ID_1[2] && LastLay )
699 {
700 Hep3Vector pos_loc = MucGeoGeneral::Instance()
701 ->GetGap( ID_1[0], ID_1[1], ID_1[2] )
702 ->TransformToGap( currentPosition );
703 double dist = fabs( pos_loc.z() );
704 RemPositon = currentPosition;
705 RemMomentum = currentMomentum;
706 RemXpErr = XpErr;
707 RemDist = dist;
708 RemStep++;
709 }
710 }
711 if ( RemStep > 0 )
712 {
713 bool WillFilter = false;
714 bool LastLay_ = false;
715 double dist_b = 99999.;
716 Hep3Vector ID_2;
717 if ( Flag4 )
718 {
719 ID_2 = GetGapID( oldVolumeName );
720 LastLay_ = ( ID_2[0] == 1 && ID_2[2] < 9 ) ||
721 ( ( ID_2[0] == 0 || ID_2[0] == 2 ) && ID_2[2] < 8 );
722 if ( LastLay_ )
723 dist_b = fabs( MucGeoGeneral::Instance()
724 ->GetGap( ID_2[0], ID_2[1], ID_2[2] )
725 ->TransformToGap( currentPosition )
726 .z() );
727 if ( RemID != ID_2 ) WillFilter = true;
728 }
729
730 Hep3Vector pos_loc = MucGeoGeneral::Instance()
731 ->GetGap( RemID[0], RemID[1], RemID[2] )
732 ->TransformToGap( currentPosition );
733 double dist = fabs( pos_loc.z() );
734 // get the nearest point from the center of gap
735 if ( !WillFilter && RemDist > dist )
736 {
737 RemPositon = currentPosition;
738 RemMomentum = currentMomentum;
739 RemXpErr = XpErr;
740 RemDist = dist;
741 RemVol = oldVolumeName;
742 }
743
744 // if find the nearest point in the gap, Open Fillter
745 if ( WillFilter )
746 {
747 ExtMucKal* muckal = new ExtMucKal();
748 muckal->SetGapID( RemID );
749 muckal->SetPosMomErr( RemPositon, RemMomentum, RemXpErr );
750 muckal->SetMucDigiCol( m_mucdigicol );
751 muckal->SetMucWindow( myMucWindow );
752 bool iniOK = muckal->MucKalIniti();
753 vector<MucRecHit*> tmp = muckal->TrackHit();
754 bool filter = muckal->ExtMucFilter();
755 // bool filter = muckal->GetFilterStatus();
756 double chi2 = muckal->GetChi2();
757 bool samestrip = muckal->GetSameStrip();
758 if ( filter && iniOK )
759 {
760
761 myMucnfit_++;
762 // cout<<"myMucnfit_: "<<myMucnfit_<<endl;
763 myMucchisq_ = myMucchisq_ + chi2;
764 muckal->XPmod( m_pos_mod, m_mom_mod, m_err_mod );
765 RememberID = muckal->GetHitGap();
766 if ( RememberID[0] == 1 ) myMucbrLastLay_ = RememberID[2];
767 if ( RememberID[0] == 0 || RememberID[0] == 2 ) myMucecLastLay_ = RememberID[2];
768 if ( oldVolumeName.contains( "Ab" ) ) RemDepth = RemDepth - depth;
769 if ( myMucnfit_ == 1 ) myMucdepth_ = RemDepth;
770 else myMucdepth_ = myMucdepth_ + RemDepth;
771
772 if ( !samestrip )
773 {
774
775 extXpErr->put_err( m_err_mod );
776 extXpErr->set_pos( m_pos_mod );
777 extXpErr->set_mom( m_mom_mod );
778 currentStep->GetTrack()->SetTrackStatus( fStopAndKill );
779 }
780 else
781 {
782 // cout<<"-------hit exist, but no modify-------"<<endl;
783 RemPositon = currentPosition;
784 RemMomentum = currentMomentum;
785 RemXpErr = XpErr;
786 RemDist = dist_b;
787 RemID = ID_2;
788 if ( oldVolumeName.contains( "Ab" ) ) RemDepth = depth;
789 else RemDepth = 0.;
790 }
791 if ( msgFlag ) cout << "---------Filter is OK---------- " << endl;
792 }
793 else
794 {
795 // if(LastLay_)
796 //{
797 RemPositon = currentPosition;
798 RemMomentum = currentMomentum;
799 RemXpErr = XpErr;
800 RemDist = dist_b;
801 RemID = ID_2;
802 // }
803 }
804 delete muckal;
805 }
806 }
807 if ( myExtTrack )
808 myExtTrack->SetMucKalData( myMucchisq_, myMucnfit_, myMucdepth_, myMucbrLastLay_,
809 myMucecLastLay_, myMucnhits_ );
810 }
811
812 /*
813 if(Flag1&&Flag2&&Flag3)
814 {
815 ExtMucKal* muckal = new ExtMucKal();
816 muckal->SetVolume(oldVolume,newVolume);
817 muckal->SetMomPosErr(currentPosition,currentMomentum,XpErr);
818 muckal->SetMucDigiCol(mucdigicol);
819 muckal->ExtMucFilter();
820 muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
821 bool filter = muckal->GetFilterStatus();
822 double chi2 = muckal->GetChi2();
823 if(filter)
824 {
825 myMucnfit_++;
826 myMucchisq_ = myMucchisq_+chi2;
827 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
828 RememberVol.assign(oldVolumeName,10);
829 extXpErr->put_err(m_err_mod);
830 extXpErr->set_pos(m_pos_mod);
831 extXpErr->set_mom(m_mom_mod);
832 }
833 delete muckal;
834 }
835 */
836 // if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
837 /* //Get Muc Ext Hits Data.
838 if(newVolumeName.contains("Strip"))
839 {
840 int newVolumeNumber = newVolume->GetCopyNo();
841 if(myExtTrack)
842 {
843 ExtMucHit aExtMucHit;
844 Hep3Vector xVector(1.0,0,0);
845 Hep3Vector yVector(0,1.0,0);
846 Hep3Vector zVector(0,0,1.0);
847 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
848 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
849 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
850 double tzError;
851 double phi = currentPosition.phi();
852 if(phi<0.) phi+=M_PI;
853 int i = int(8.0*phi/M_PI);
854 if( i==0 || i==7 || i==8 )
855 {
856 tzError = yError;
857 }
858 if( i==1 || i==2 )
859 {
860 Hep3Vector tzVector(-1.0,1.0,0.);
861 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
862 }
863 if( i==3 || i==4 )
864 {
865 tzError = xError;
866 }
867 if( i==5 || i==6 )
868 {
869 Hep3Vector tzVector(1.0,1.0,0.);
870 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
871 }
872 aExtMucHit.SetExtMucHit(currentPosition,currentMomentum,newVolumeName,newVolumeNumber,extXpErr->get_err(),zError,tzError,xError,yError);
873 myExtTrack->AddExtMucHit(aExtMucHit);
874 }
875 }
876 */
877 }
878 else
879 {
880 if ( msgFlag )
881 cout << "x:" << currentPosition.x() << "//y:" << currentPosition.y()
882 << "//z:" << currentPosition.z() << "||px:" << currentMomentum.x()
883 << "//py:" << currentMomentum.y() << "//pz:" << currentMomentum.z() << endl;
884 double x = currentPosition.x();
885 double y = currentPosition.y();
886 double z = currentPosition.z();
887 if ( ( fabs( x ) >= 2 * myMucR ) || ( fabs( y ) >= 2 * myMucR ) ||
888 ( fabs( z ) >= 2 * myMucZ ) )
889 // currentTrack->SetKineticEnergy(0.0);// protection in case that a very energetic track
890 // travels without touching BESIII
891 {
892 currentStep->GetTrack()->SetTrackStatus( fStopAndKill );
893 m_trackstop = true;
894 }
895 }
896 }
897 else if ( currentTrackStatus == fStopAndKill )
898 {
899 m_trackstop = true;
900 if ( myEmcFlag ) myExtTrack->SetEmcPath( myPathInCrystal / 10. );
901 if ( myTof1Flag ) myExtTrack->setPathInTof1( myPathInTof1 );
902 if ( myTof2Flag ) myExtTrack->setPathInTof2( myPathInTof2 );
903 if ( msgFlag )
904 {
905 cout << "myPathInTof1 = ";
906 for ( int i = 0; i != myPathInTof1.size(); i++ ) cout << myPathInTof1[i] << " ";
907 cout << endl;
908 cout << "myPathInTof2 = ";
909 for ( int i = 0; i != myPathInTof2.size(); i++ ) cout << myPathInTof2[i] << " ";
910 cout << endl;
911 }
912
913 if ( msgFlag )
914 {
915 if ( newVolume != 0 )
916 {
917 std::cout << "x:" << currentPosition.x() << "//y:" << currentPosition.y()
918 << "//z:" << currentPosition.z() << "||px:" << currentMomentum.x()
919 << "//py:" << currentMomentum.y() << "//pz:" << currentMomentum.z()
920 << "//stoped" << endl;
921 cout << "Error matrix is:" << extXpErr->get_err() << endl;
922 }
923 else
924 {
925 cout << "x:" << currentPosition.x() << "//y:" << currentPosition.y()
926 << "//z:" << currentPosition.z() << "||px:" << currentMomentum.x()
927 << "//py:" << currentMomentum.y() << "//pz:" << currentMomentum.z() << "//escaped"
928 << std::endl;
929 std::cout << "Error matrix is:" << extXpErr->get_err() << std::endl;
930 }
931 }
932 }
933}
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
nPhi
const int nPhi
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TofSim/BesTofDigitizerEcV3.hh:36
M_PI
#define M_PI
Definition TConstant.h:4
ExtMucKal::ExtMucFilter
bool ExtMucFilter()
Definition ExtMucKal.cxx:22
ExtMucKal::SetGapID
void SetGapID(Hep3Vector id)
Definition ExtMucKal.h:30
ExtMucKal::GetSameStrip
bool GetSameStrip()
Definition ExtMucKal.h:44
ExtMucKal::SetPosMomErr
void SetPosMomErr(Hep3Vector pos, Hep3Vector mom, HepSymMatrix err)
Definition ExtMucKal.h:23
ExtMucKal::SetMucWindow
void SetMucWindow(int aMucWindow)
Definition ExtMucKal.h:29
ExtMucKal::XPmod
void XPmod(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition ExtMucKal.cxx:87
ExtMucKal::SetMucDigiCol
void SetMucDigiCol(MucDigiCol *amucdigi)
Definition ExtMucKal.h:28
ExtMucKal::GetChi2
double GetChi2()
Definition ExtMucKal.h:37
ExtMucKal::TrackHit
vector< MucRecHit * > TrackHit()
Definition ExtMucKal.cxx:93
ExtMucKal::MucKalIniti
bool MucKalIniti()
Definition ExtMucKal.cxx:169
ExtMucKal::GetHitGap
Hep3Vector GetHitGap()
Definition ExtMucKal.cxx:252
ExtSteppingAction::GetGapID
Hep3Vector GetGapID(G4String vol)
Definition ExtSteppingAction.cxx:1079
ExtSteppingAction::CalculateEmcEndThetaPhi
void CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
Definition ExtSteppingAction.cxx:983
ExtSteppingAction::CalculateEmcEndCopyNb
int CalculateEmcEndCopyNb(int num)
Definition ExtSteppingAction.cxx:1059
ExtSteppingAction::CalculateEmcEndPhiNb
int CalculateEmcEndPhiNb(int num)
Definition ExtSteppingAction.cxx:1052
MucGeoGap::TransformToGap
HepPoint3D TransformToGap(const HepPoint3D gPoint) const
Transform a point from global coordinate to gap coordinate.
Definition MucGeoGap.cxx:580
MucGeoGeneral::GetGap
MucGeoGap * GetGap(const int part, const int seg, const int gap) const
Get a pointer to the gap identified by (part,seg,gap).
Definition MucGeoGeneral.cxx:435
MucGeoGeneral::Instance
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
Definition MucGeoGeneral.cxx:428
rb::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:22
The documentation for this class was generated from the following files: