CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
ExtSteppingAction.cxx
Go to the documentation of this file.
1//
2//File: ExtSteppingAction.cc
3//Date: 2005.3.17
4//Author: L.L.Wang
5//History: 2005.3.17 created by Wang Liangliang
6// 2005.8.12 find a mistake in caculation of error direction of tof and corrected by Wang Liangliang
7
8//#include "CLHEP/config/CLHEP.h"
9
10#include "TrkExtAlg/ExtSteppingAction.h"
11
12#include "G4Track.hh"
13#include "G4TrackStatus.hh"
14#include "G4VPhysicalVolume.hh"
15#include "G4TransportationManager.hh"
16#include "G4FieldManager.hh"
17#include "G4Field.hh"
18#include "G4LogicalVolume.hh"
19#include "globals.hh"
20
21#include "TofSim/BesTofDigitizerEcV4.hh"
22#include "TofSim/BesTofDigitizerEcV4_dbs.hh"
23#include "MucRecEvent/MucRecHit.h"
24#include "TrkExtAlg/ExtMucKal.h"
25#include "TrkExtAlg/Ext_xp_err.h"
26#include <iostream>
27#include <string>
28#include "strstream"
29
30using namespace std;
31
32
33ExtSteppingAction::ExtSteppingAction():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}
45
46ExtSteppingAction::~ExtSteppingAction(){}
47
48void ExtSteppingAction::Reset()
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}
89
90void ExtSteppingAction::MucReset()
91{
92RemStep =0;
93RemDist = 99999.;
94RemDepth = 0.;
95RemID[0]=-1;
96RemID[1]=-1;
97RemID[2]=-1;
98RemVol = "";
99}
100
101
102
103void ExtSteppingAction::UserSteppingAction(const G4Step* currentStep)
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}
851
852
853
854void ExtSteppingAction::CalculateChicc(G4Material* currentMaterial)
855{
856 int n = currentMaterial->GetNumberOfElements();
857 const double *p = currentMaterial->GetFractionVector();
858 double density = (currentMaterial->GetDensity())/(g/cm3);
859 double temp=0.0;
860 for(int i=0;i<n;i++)
861 {
862 const G4Element* mElement = currentMaterial->GetElement(i);
863 double z = mElement->GetZ();
864 double a = mElement->GetN();
865 // std::cout<<"Fraction: "<<*p<<" z: "<<z<<" a: "<<a<<std::endl;
866 temp += *(p++)/a*z*(z+1);
867 }
868 chicc = 0.39612e-3*std::sqrt(density*temp);
869 // std::cout<<"chicc:"<<chicc<<std::endl;
870}
871
872
873HepSymMatrix & ExtSteppingAction::myOutputSymMatrix(const HepSymMatrix & m1)
874{
875 // std::cout<<"myOutputSymMatrix:1"<<std::endl;
876 HepSymMatrix m;
877 // std::cout<<"myOutputSymMatrix:2"<<std::endl;
878 m=m1;
879 // std::cout<<"myOutputSymMatrix:3"<<std::endl;
880 for(int i=0;i<6;i++)
881 { for(int j=0;j<=i;j++)
882 { if(i<3) {
883 m[i][j]=m[i][j]/100; //mm*mm --> cm*cm
884 }
885 else if(j<3) {
886 m[i][j]=m[i][j]/10000; //mm*MeV --> cm*GeV
887 }
888 else {
889 m[i][j]=m[i][j]/1000000; //MeV*MeV --> GeV*GeV
890 }
891 }
892 }
893 // std::cout<<"myOutputSymMatrix:4"<<std::endl;
894 // std::cout<<"m1="<<m1<<std::endl;
895 // std::cout<<"m="<<m<<std::endl;
896 myOutputSM=m;
897 return myOutputSM;
898}
899
900void ExtSteppingAction::CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
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}
982
983int ExtSteppingAction::CalculateEmcEndPhiNb(int num)
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}
995
996int ExtSteppingAction::CalculateEmcEndCopyNb(int num)
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}
1021
1022void ExtSteppingAction::InfmodMuc(Hep3Vector &pos,Hep3Vector &mom,HepSymMatrix &err)
1023{
1024 pos = m_pos_mod;
1025 mom = m_mom_mod;
1026 err = m_err_mod;
1027}
1028
1029
1030Hep3Vector ExtSteppingAction::GetGapID(G4String vol)
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}
1055
1056
nPhi
const int nPhi
Definition: BesTofDigitizerEcV3.hh:37
BesTofDigitizerEcV4_dbs.hh
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
Ext_xp_err.h
help
std::string help()
Definition: HelloServer.cpp:35
M_PI
#define M_PI
Definition: TConstant.h:4
BesTofDigitizerEcV4_dbs::Calculate_Readoutstrip_number_continuum
static G4int Calculate_Readoutstrip_number_continuum(G4double x_mm, G4double y_mm, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4_dbs.cc:461
BesTofDigitizerEcV4_dbs::Produce_unique_identifier
static G4int Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
Definition: BesTofDigitizerEcV4_dbs.cc:876
BesTofDigitizerEcV4::Produce_unique_identifier
static G4int Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
Definition: BesTofDigitizerEcV4.cc:841
BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum
static G4int Calculate_Readoutstrip_number_continuum(G4double x_mm, G4double y_mm, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4.cc:455
DstExtTrack::SetTof1Data
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.)
Definition: DstExtTrack.cxx:162
DstExtTrack::SetEmcData
void SetEmcData(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, double aThetaSigma, double aPhiSigma, HepSymMatrix aErrorMatrix)
Definition: DstExtTrack.cxx:196
DstExtTrack::SetTof2Data
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.)
Definition: DstExtTrack.cxx:179
DstExtTrack::SetMucData
void SetMucData(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, HepSymMatrix aErrorMatrix, double aZSigma=0., double aTSigma=0., double aXSigma=0., double aYSigma=0.)
Definition: DstExtTrack.cxx:209
DstExtTrack::SetEmcPath
void SetEmcPath(double path)
Definition: DstExtTrack.h:214
DstExtTrack::SetMucKalData
void SetMucKalData(double chi2, int dof, double depth, int brLastLay, int ecLastLay, int nhits)
Definition: DstExtTrack.cxx:226
ExtMucKal::ExtMucFilter
bool ExtMucFilter()
Definition: ExtMucKal.cxx:16
ExtMucKal::SetGapID
void SetGapID(Hep3Vector id)
Definition: ExtMucKal.h:28
ExtMucKal::GetSameStrip
bool GetSameStrip()
Definition: ExtMucKal.h:42
ExtMucKal::SetPosMomErr
void SetPosMomErr(Hep3Vector pos, Hep3Vector mom, HepSymMatrix err)
Definition: ExtMucKal.h:25
ExtMucKal::SetMucWindow
void SetMucWindow(int aMucWindow)
Definition: ExtMucKal.h:27
ExtMucKal::XPmod
void XPmod(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtMucKal.cxx:71
ExtMucKal::SetMucDigiCol
void SetMucDigiCol(MucDigiCol *amucdigi)
Definition: ExtMucKal.h:26
ExtMucKal::GetChi2
double GetChi2()
Definition: ExtMucKal.h:35
ExtMucKal::TrackHit
vector< MucRecHit * > TrackHit()
Definition: ExtMucKal.cxx:78
ExtMucKal::MucKalIniti
bool MucKalIniti()
Definition: ExtMucKal.cxx:150
ExtMucKal::GetHitGap
Hep3Vector GetHitGap()
Definition: ExtMucKal.cxx:228
ExtSteppingAction::Get_which_tof_version
int Get_which_tof_version(void)
Definition: ExtSteppingAction.h:51
ExtSteppingAction::GetGapID
Hep3Vector GetGapID(G4String vol)
Definition: ExtSteppingAction.cxx:1030
ExtSteppingAction::CalculateEmcEndThetaPhi
void CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
Definition: ExtSteppingAction.cxx:900
ExtSteppingAction::InfmodMuc
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtSteppingAction.cxx:1022
ExtSteppingAction::UserSteppingAction
void UserSteppingAction(const G4Step *currentStep)
Definition: ExtSteppingAction.cxx:103
ExtSteppingAction::CalculateEmcEndCopyNb
int CalculateEmcEndCopyNb(int num)
Definition: ExtSteppingAction.cxx:996
ExtSteppingAction::CalculateEmcEndPhiNb
int CalculateEmcEndPhiNb(int num)
Definition: ExtSteppingAction.cxx:983
Ext_errmx::put_err
void put_err(const double error[])
Definition: Ext_errmx.cxx:45
Ext_errmx::get_err
const HepSymMatrix & get_err() const
Definition: Ext_errmx.h:157
Ext_errmx::get_plane_err
double get_plane_err(const Hep3Vector &np, const Hep3Vector &nr) const
Definition: Ext_errmx.cxx:85
Ext_xp_err::set_mom
void set_mom(const Hep3Vector &mom)
Definition: Ext_xp_err.h:90
Ext_xp_err::move
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
Ext_xp_err::set_pos
void set_pos(const Hep3Vector &pos)
Definition: Ext_xp_err.h:85
MucGeoGap::TransformToGap
HepPoint3D TransformToGap(const HepPoint3D gPoint) const
Transform a point from global coordinate to gap coordinate.
Definition: MucGeoGap.cxx:627
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:450
MucGeoGeneral::Instance
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
Definition: MucGeoGeneral.cxx:439
RecExtTrack::setPathInTof1
void setPathInTof1(vector< double > x)
Definition: RecExtTrack.h:61
RecExtTrack::setPathInTof2
void setPathInTof2(vector< double > x)
Definition: RecExtTrack.h:62
vec
dble_vec_t vec[12]
Definition: ranlxd.c:372
num
int num[96]
Definition: ranlxd.c:373