CGEM BOSS 6.6.5.g
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 25 of file ExtSteppingAction.h.

Constructor & Destructor Documentation

◆ ExtSteppingAction()

ExtSteppingAction::ExtSteppingAction ( )

Definition at line 33 of file ExtSteppingAction.cxx.

33 :chicc(0.0),initialPath(0.),myPathIntoCrystal(0.),myPathOutCrystal(0.),myPathInCrystal(0.),myBetaInMDC(0.),extXpErr(0),myExtTrack(0),msgFlag(true),myUseMucKalFlag(true),m_trackstop(false),myMucnfit_(0),myMucchisq_(0.),myMucdepth_(-99.),myMucbrLastLay_(-1),myMucecLastLay_(-1),myMucnhits_(-1),myMucWindow(6)
34{
35 myTof1R = 814.001;
36 myTof1Z = 1330.0;
37 myTof2R = 871.036;
38 myEmcR1 = 945.0;
39 myEmcR2 = 983.9;
40 myEmcZ = 1416.8;
41 myMucR = 1700.0-0.01;
42 myMucZ = 2050.0-0.01;
43 m_which_tof_version=2;
44}

◆ ~ExtSteppingAction()

ExtSteppingAction::~ExtSteppingAction ( )

Definition at line 46 of file ExtSteppingAction.cxx.

46{}

Member Function Documentation

◆ CalculateEmcEndCopyNb()

int ExtSteppingAction::CalculateEmcEndCopyNb ( int  num)

Definition at line 996 of file ExtSteppingAction.cxx.

997{
998 int copyNb;
999 switch(num){
1000 case 30:
1001 copyNb = 5;
1002 break;
1003 case 31:
1004 copyNb = 6;
1005 break;
1006 case 32:
1007 copyNb = 14;
1008 break;
1009 case 33:
1010 copyNb = 15;
1011 break;
1012 case 34:
1013 copyNb = 16;
1014 break;
1015 default:
1016 copyNb = num;
1017 break;
1018 }
1019 return copyNb;
1020}
int num[96]
Definition: ranlxd.c:373

Referenced by UserSteppingAction().

◆ CalculateEmcEndPhiNb()

int ExtSteppingAction::CalculateEmcEndPhiNb ( int  num)

Definition at line 983 of file ExtSteppingAction.cxx.

984{
985 if(num==0||num==1) {
986 return 15-num*8;
987 } else if(num==2||num==3) {
988 return 30-num*8;
989 } else if(num<=9) {
990 return 17-num;
991 } else {
992 return 15-num;
993 }
994}

Referenced by UserSteppingAction().

◆ CalculateEmcEndThetaPhi()

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

Definition at line 900 of file ExtSteppingAction.cxx.

901{
902 if((sector>=0)&&(sector<3))
903 sector+=16;
904 //if((sector!=7)&&(sector!=15))
905 {
906 if((nb>=0)&&(nb<4))
907 {
908 ntheta = 0;
909 nphi = (sector-3)*4+nb;
910 }
911 else if((nb>=4)&&(nb<8))
912 {
913 ntheta = 1;
914 nphi = (sector-3)*4+nb-4;
915 }
916 else if((nb>=8)&&(nb<13))
917 {
918 ntheta = 2;
919 nphi = (sector-3)*5+nb-8;
920 }
921 else if((nb>=13)&&(nb<18))
922 {
923 ntheta = 3;
924 nphi = (sector-3)*5+nb-13;
925 }
926 else if((nb>=18)&&(nb<24))
927 {
928 ntheta = 4;
929 nphi = (sector-3)*6+nb-18;
930 }
931 else if((nb>=24)&&(nb<30))
932 {
933 ntheta = 5;
934 nphi = (sector-3)*6+nb-24;
935 }
936 }
937
938 if(npart==2)
939 {
940 switch(ntheta)
941 {
942 case 0:
943 if(nphi<32)
944 nphi = 31-nphi;
945 else
946 nphi = 95-nphi;
947 break;
948 case 1:
949 if(nphi<32)
950 nphi = 31-nphi;
951 else
952 nphi = 95-nphi;
953 break;
954 case 2:
955 if(nphi<40)
956 nphi = 39-nphi;
957 else
958 nphi = 119-nphi;
959 break;
960 case 3:
961 if(nphi<40)
962 nphi = 39-nphi;
963 else
964 nphi = 119-nphi;
965 break;
966 case 4:
967 if(nphi<48)
968 nphi = 47-nphi;
969 else
970 nphi = 143-nphi;
971 break;
972 case 5:
973 if(nphi<48)
974 nphi = 47-nphi;
975 else
976 nphi = 143-nphi;
977 default:
978 ;
979 }
980 }
981}

Referenced by UserSteppingAction().

◆ Get_which_tof_version()

int ExtSteppingAction::Get_which_tof_version ( void  )
inline

Definition at line 51 of file ExtSteppingAction.h.

51{return m_which_tof_version;}

Referenced by UserSteppingAction().

◆ GetGapID()

Hep3Vector ExtSteppingAction::GetGapID ( G4String  vol)

Definition at line 1030 of file ExtSteppingAction.cxx.

1031{
1032 int par(-1),se(-1),ga(-1);
1033 G4String strPart = vol.substr(5,1);
1034 //G4String strSeg = m_MotherVolumeName.substr(7,1);
1035
1036 G4String strSeg = vol.substr(7,1);
1037 G4String strGap;
1038 if(vol.contains("G")) strGap= vol.substr(9,1);
1039 if(vol.contains("Ab")) strGap= vol.substr(10,1);
1040 std::istrstream partBuf(strPart.c_str(), strlen(strPart.c_str()));
1041 std::istrstream segBuf(strSeg.c_str(), strlen(strSeg.c_str()));
1042 std::istrstream gapBuf(strGap.c_str(), strlen(strGap.c_str()));
1043
1044 partBuf >> par ;
1045 segBuf >> se;
1046 gapBuf >> ga;
1047 if(vol.contains("Ab")&&par==1) ga = ga+1;
1048 //panelBuf >> m_panel;
1049 Hep3Vector vec;
1050 vec[0]=par;
1051 vec[1]=se;
1052 vec[2]=ga;
1053 return vec;
1054}
dble_vec_t vec[12]
Definition: ranlxd.c:372

Referenced by UserSteppingAction().

◆ InfmodMuc()

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

Definition at line 1022 of file ExtSteppingAction.cxx.

1023{
1024 pos = m_pos_mod;
1025 mom = m_mom_mod;
1026 err = m_err_mod;
1027}

Referenced by TrkExtAlg::execute().

◆ MucReset()

void ExtSteppingAction::MucReset ( )

Definition at line 90 of file ExtSteppingAction.cxx.

91{
92RemStep =0;
93RemDist = 99999.;
94RemDepth = 0.;
95RemID[0]=-1;
96RemID[1]=-1;
97RemID[2]=-1;
98RemVol = "";
99}

Referenced by Ext_track::Set().

◆ Reset()

void ExtSteppingAction::Reset ( )

Definition at line 48 of file ExtSteppingAction.cxx.

49{
50 chicc = 0.;
51 myExtTrack = 0;
52 myPathIntoCrystal = 0.;
53 myPathOutCrystal = 0.;
54 myPathInCrystal = 0.;
55
56 myPathIntoTof1 = 0.0;
57 myPathOutTof1 = 0.0;
58 //myPathInTof1 = 0.0;
59 myPathInTof1.clear();
60
61 myPathIntoTof2 = 0.0;
62 myPathOutTof2 = 0.0;
63 //myPathInTof2 = 0.0;
64 myPathInTof2.clear();
65
66 myTofFlag = false;
67 myTof1Flag = false;
68 myTof2Flag = false;
69 myInTof1 = false;
70 myInTof2 = false;
71 myOutTof1 = false;
72 myOutTof2 = false;
73 myPhyEmcFlag = false;
74 myEmcFlag = false;
75 myEmcPathFlag= false;
76 myMucFlag = false;
77
78 m_trackstop =false;
79 myMucchisq_ =0.;
80 myMucnfit_ =0;
81 myMucdepth_=-99.;
82 myMucbrLastLay_=-1;
83 myMucecLastLay_=-1;
84 RememberID[0]=-1;
85 RememberID[1]=-1;
86 RememberID[2]=-1;
87
88}

Referenced by TrkExtAlg::execute().

◆ Set_which_tof_version()

void ExtSteppingAction::Set_which_tof_version ( int  version)
inline

Definition at line 50 of file ExtSteppingAction.h.

50{m_which_tof_version=version;}

Referenced by TrkExtAlg::execute().

◆ SetBetaInMDC()

void ExtSteppingAction::SetBetaInMDC ( double  aBeta)
inline

Definition at line 38 of file ExtSteppingAction.h.

38{myBetaInMDC = aBeta;};

Referenced by Ext_track::Set().

◆ SetExtTrackPointer()

void ExtSteppingAction::SetExtTrackPointer ( RecExtTrack aExtTrack)
inline

Definition at line 44 of file ExtSteppingAction.h.

44{myExtTrack = aExtTrack;};

Referenced by TrkExtAlg::execute().

◆ SetInitialPath()

void ExtSteppingAction::SetInitialPath ( double  aPath)
inline

Definition at line 36 of file ExtSteppingAction.h.

36{initialPath = aPath;};

Referenced by Ext_track::Set().

◆ SetInitialTof()

void ExtSteppingAction::SetInitialTof ( double  aTof)
inline

Definition at line 37 of file ExtSteppingAction.h.

37{initialTof = aTof;};

Referenced by Ext_track::Set().

◆ SetMsgFlag()

void ExtSteppingAction::SetMsgFlag ( bool  aMsgFalg)
inline

Definition at line 41 of file ExtSteppingAction.h.

41{msgFlag = aMsgFalg;};

Referenced by Ext_track::Initialization().

◆ SetMucDigiColPointer()

void ExtSteppingAction::SetMucDigiColPointer ( MucDigiCol rawdigicol)
inline

Definition at line 56 of file ExtSteppingAction.h.

56{ m_mucdigicol = rawdigicol;}

Referenced by TrkExtAlg::execute().

◆ SetMucKalFlag()

void ExtSteppingAction::SetMucKalFlag ( bool  aMucKalFlag)
inline

Definition at line 42 of file ExtSteppingAction.h.

42{myUseMucKalFlag=aMucKalFlag;};

Referenced by Ext_track::Initialization().

◆ SetMucWindow()

void ExtSteppingAction::SetMucWindow ( int  aMucWindow)
inline

Definition at line 43 of file ExtSteppingAction.h.

43{myMucWindow=aMucWindow;};

Referenced by Ext_track::Initialization().

◆ SetXpErrPointer()

void ExtSteppingAction::SetXpErrPointer ( Ext_xp_err xpErr)
inline

Definition at line 39 of file ExtSteppingAction.h.

39{extXpErr = xpErr;};

Referenced by Ext_track::Set().

◆ TrackStop()

bool ExtSteppingAction::TrackStop ( )
inline

Definition at line 55 of file ExtSteppingAction.h.

55{return m_trackstop;}

Referenced by TrkExtAlg::execute().

◆ UserSteppingAction()

void ExtSteppingAction::UserSteppingAction ( const G4Step *  currentStep)

Definition at line 103 of file ExtSteppingAction.cxx.

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

The documentation for this class was generated from the following files: