BOSS 6.6.4.p01
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
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}
44
45ExtSteppingAction::~ExtSteppingAction(){}
46
47void ExtSteppingAction::Reset()
48{
49 chicc = 0.;
50 myExtTrack = 0;
51 myPathIntoCrystal = 0.;
52 myPathOutCrystal = 0.;
53 myPathInCrystal = 0.;
54
55 myPathIntoTof1 = 0.0;
56 myPathOutTof1 = 0.0;
57 //myPathInTof1 = 0.0;
58 myPathInTof1.clear();
59
60 myPathIntoTof2 = 0.0;
61 myPathOutTof2 = 0.0;
62 //myPathInTof2 = 0.0;
63 myPathInTof2.clear();
64
65 myTofFlag = false;
66 myTof1Flag = false;
67 myTof2Flag = false;
68 myInTof1 = false;
69 myInTof2 = false;
70 myOutTof1 = false;
71 myOutTof2 = false;
72 myPhyEmcFlag = false;
73 myEmcFlag = false;
74 myEmcPathFlag= false;
75 myMucFlag = false;
76
77 m_trackstop =false;
78 myMucchisq_ =0.;
79 myMucnfit_ =0;
80 myMucdepth_=-99.;
81 myMucbrLastLay_=-1;
82 myMucecLastLay_=-1;
83 RememberID[0]=-1;
84 RememberID[1]=-1;
85 RememberID[2]=-1;
86
87}
88
89void ExtSteppingAction::MucReset()
90{
91RemStep =0;
92RemDist = 99999.;
93RemDepth = 0.;
94RemID[0]=-1;
95RemID[1]=-1;
96RemID[2]=-1;
97RemVol = "";
98}
99
100
101
102void ExtSteppingAction::UserSteppingAction(const G4Step* currentStep)
103{
104
105 //Get track and its status
106 G4Track* currentTrack = currentStep->GetTrack();
107 if(!currentTrack)
108 {
109 cout<<"Can't get currentTrack"<<endl;
110 return;
111 }
112 G4TrackStatus currentTrackStatus = currentTrack->GetTrackStatus();
113
114 int stepNumber = currentTrack->GetCurrentStepNumber();
115 if(msgFlag) cout<<"STEP "<<stepNumber<<":"<<endl;
116
117 //Get current position and momentum
118 Hep3Vector currentPosition = currentTrack->GetPosition();
119 Hep3Vector currentMomentum = currentTrack->GetMomentum();
120
121
122 //Get current Volume
123 G4VPhysicalVolume* oldVolume;
124 G4VPhysicalVolume* newVolume;
125 if(!currentTrack->GetTouchableHandle()) cout<<"Can't get currentTrack->GetTouchableHandle()"<<endl;
126 else oldVolume= currentTrack->GetTouchableHandle()->GetVolume();
127 if(!currentTrack->GetNextTouchableHandle()) cout<<"Can't get currentTrack->GetNextTouchableHandle()"<<endl;
128 else newVolume= currentTrack->GetNextTouchableHandle()->GetVolume();
129 if(!oldVolume) cout<<"Can't get oldVolume!"<<endl;
130
131 //****added by LI Chunhua
132 if(stepNumber>50000) {
133 cout<<"infinite steps in the track "<<endl;
134 currentTrack->SetTrackStatus(fStopAndKill); m_trackstop =true;
135 }
136
137 G4String ParticleName = currentTrack->GetDefinition()->GetParticleName();
138 double x_ = currentPosition.x();
139 double y_ = currentPosition.y();
140 double z_ = currentPosition.z();
141 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)) );
142 bool mucdec = (fabs(x_)>=myMucR) || (fabs(y_)>=myMucR) || ((fabs(x_-y_)/sqrt(2.))>=myMucR) || ((fabs(x_+y_)/sqrt(2.))>=myMucR) || (fabs(z_)>=myMucZ);
143
144 //if current particle is alive.
145 if(currentTrackStatus == fAlive)
146 {
147 //if(oldVolume!=newVolume)//if so, update Error Matrix and chicc.
148 if(inner||mucdec)//update Error Matrix for newvolume in MDC,TOF,EMC; update Error Matrix step by step in MUC;
149 {
150 //Get current B field
151 double currentPoint[3]={currentPosition.x(),currentPosition.y(),currentPosition.z()};
152 double currentBfield[3]={0.0,0.0,0.0};
153 Hep3Vector currentB(0.0,0.0,1.0);
154 if(G4TransportationManager::GetTransportationManager())
155 {
156 G4FieldManager* fieldManager=G4TransportationManager::GetTransportationManager()->GetFieldManager();
157 if(fieldManager)
158 {
159 if(fieldManager->GetDetectorField())
160 {
161 fieldManager->GetDetectorField()->GetFieldValue(currentPoint,currentBfield);
162 currentB.set(currentBfield[0]/tesla,currentBfield[1]/tesla,currentBfield[2]/tesla);
163 if(msgFlag) cout<<"B:"<<currentB.x()<<","<<currentB.y()<<","<<currentB.z()<<endl;
164 }
165 }
166 }
167
168 //Get chicc
169 G4Material* oldMaterial = oldVolume->GetLogicalVolume()->GetMaterial();
170 if(!oldMaterial) std::cout<<"Can't get oldMaterial"<<std::endl;
171 else CalculateChicc(oldMaterial);
172
173 //Update Error Matrix
174 if(!(extXpErr->move(currentPosition,currentMomentum,currentB,1,chicc)))
175 if(msgFlag) cout<<"can not update Error Matrix!!"<<endl;
176
177 //Verbose
178 if(msgFlag)
179 {
180 cout<<"Volume name:"<<newVolume->GetName()<<endl;
181 cout<<"Volume number:"<<newVolume->GetCopyNo()<<endl;
182 cout<<"x:"<<currentPoint[0]<<"//y:"<<currentPoint[1]<<"//z:"<<currentPoint[2]
183 <<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"
184 <<currentMomentum.z()<<endl;
185 cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
186 cout<<"phi:"<<atan(currentPoint[1]/currentPoint[0])<<endl;
187 Hep3Vector nz(0.,0.,1.);
188 cout<<"Projected z error:"<<extXpErr->get_plane_err(currentMomentum.unit(),nz)
189 <<endl;
190 double x,y,R;
191 x = currentPoint[0];
192 y = currentPoint[1];
193 R = sqrt(x*x+y*y);
194 Hep3Vector nt(-y/R,x/R,0.);
195 cout<<"Projected phi error:"<<(extXpErr->get_plane_err(currentMomentum.unit(),nt))/R
196 <<endl<<endl;
197 }
198
199 //some often used things
200 Hep3Vector xVector(1.0,0,0);
201 Hep3Vector yVector(0,1.0,0);
202 Hep3Vector zVector(0,0,1.0);
203 Hep3Vector tzVector;
204 tzVector.setRThetaPhi(1.0,M_PI/2.0,currentPosition.phi());
205 double r = currentPosition.perp();
206 double x = currentPosition.x();
207 double y = currentPosition.y();
208 double z = currentPosition.z();
209 G4String newVolumeName = newVolume->GetName();
210 G4String oldVolumeName = oldVolume->GetName();
211 G4StepPoint* postStepPoint = currentStep->GetPostStepPoint();
212 G4TouchableHistory* theTouchable = (G4TouchableHistory*)(postStepPoint->GetTouchable());
213 int newVolumeNumber=newVolume->GetCopyNo();
214 // int newVolumeNumber=theTouchable->GetReplicaNumber(2);
215
216 int help_mrpc_nb=-2;
217
218
219
220
221
222
223 //Get PhyTof data
224
225 // std::cout << "ExtSteppingAction NewVolumeName: " << newVolumeName << std::endl;
226 // std::cout << "ExtSteppingAction OldVolumeName: " << oldVolumeName << std::endl;
227 if( (!myTofFlag) && (!myTof1Flag) && (newVolumeName.contains("Tof") ))
228 {
229 newVolumeNumber = -2;
230 double currentTrackLength = currentTrack->GetTrackLength();
231 double totalTrackLength = currentTrackLength + initialPath;
232 //double initialTof = initialPath/(myBetaInMDC*299.792458);
233 //cout<<"initialTof = "<<initialTof<<endl;
234 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
235 double localTime = currentTrack->GetLocalTime();
236 //cout<<"localTime = "<<localTime<<endl;
237 double totalTOF = localTime + initialTof;
238 //cout<<"totalTOF= "<<totalTOF<<endl;
239
240 //std::cout << "ExtSteppingAction Volumename contains one of the marker words" << std::endl;
241
242 if(myExtTrack)
243 {
244 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
245 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
246 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
247 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
248 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);
249 myTofFlag = true;
250 }
251 return;
252 }//close if (!myTofFlag)
253
254
255
256
257
258
259
260 //Get Tof layer1 Ext data
261 if( (!myTof1Flag) && (newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ||
262 newVolumeName.contains("logical_sensitive_detector_east") || newVolumeName.contains("logical_sensitive_detector_west") ) )
263 {
264 double currentTrackLength = currentTrack->GetTrackLength();
265 double totalTrackLength = currentTrackLength+initialPath;
266 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
267 double localTime = currentTrack->GetLocalTime();
268 double totalTOF = localTime + initialTof;
269 myInTof1 = true;
270 myPathIntoTof1 = currentTrackLength;
271 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
272
273
274
275 newVolumeNumber=theTouchable->GetReplicaNumber(2);
276 help_mrpc_nb = theTouchable->GetReplicaNumber(3); // for the mrpc
277
278 if(newVolumeName.contains("ScinEc")) {
279 newVolumeNumber=(95-newVolumeNumber)/2;
280 newVolumeNumber=newVolumeNumber+176;
281
282 if(newVolumeName.contains("East")) newVolumeNumber=newVolumeNumber+48;
283
284 //std::cout<< "ExtSteppingAction Endcapart! Volumenumber" << newVolumeNumber << std::endl;
285 }
286 else if( newVolumeName.contains("_west_1")) //New MRPC - Layer 1
287 {
288
289 newVolumeNumber = (help_mrpc_nb)*(-0.5)+18.5;
290
291 int partID_help =5;
292
293 int help = BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(currentPosition.x(),currentPosition.y(), partID_help, newVolumeNumber);
294 newVolumeNumber = BesTofDigitizerEcV4::Produce_unique_identifier(newVolumeNumber, help);
295 /*
296 std::cout << "ExtSteppingAction partID | ModulNumber " << partID_help << " | " << (help_mrpc_nb)*(-0.5)+18.5 << std::endl;
297 std::cout << "ExtSteppingAction Readoutstripnumber = " << help << std::endl;
298 std::cout << "ExtSteppingAction unique identifier = " << newVolumeNumber << std::endl;
299 if(help==0){std::cout << "ExtSteppingAction x | y" <<currentPosition.x()<< " | " << currentPosition.y() << std::endl;}
300 */
301 //std::cout << "ExtSteppingAction You are using the new MRPC - Layer 1. Modulenumber " << newVolumeNumber << std::endl;
302 //std::cout << "ExtSteppingAction Current Position (x,y,z) = (" << currentPosition.x() << " | " << currentPosition.y() << " | " << currentPosition.z() << " )" << std::endl;
303 //std::cout << "ExtSteppingAction Current Position (x,y,z)*mm = (" << currentPosition.x()*mm << " | " << currentPosition.y()*mm << " | " << currentPosition.z()*mm << " )" << std::endl;
304 }
305 else if(newVolumeName.contains("_east_1") ) //New MRPC - Layer 1
306 {
307
308 newVolumeNumber = (help_mrpc_nb)*(0.5)+0.5;
309
310 int partID_help =4;
311
312 int help = BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(currentPosition.x(),currentPosition.y(), partID_help, newVolumeNumber);
313 newVolumeNumber = BesTofDigitizerEcV4::Produce_unique_identifier(newVolumeNumber, help);
314
315 /*
316 std::cout << "ExtSteppingAction partID | ModulNumber " << partID_help << " | " << (help_mrpc_nb)*(0.5)+0.5 << std::endl;
317 std::cout << "ExtSteppingAction Readoutstripnumber = " << help << std::endl;
318 std::cout << "ExtSteppingAction unique identifier = " << newVolumeNumber << std::endl;
319 if(help==0){std::cout << "ExtSteppingAction x | y" <<currentPosition.x()<< " | " << currentPosition.y() << std::endl;}
320 */
321 //std::cout << "ExtSteppingAction You are using the new MRPC - Layer 1. Modulenumber " << newVolumeNumber << std::endl;
322 //std::cout << "ExtSteppingAction Current Position (x,y,z) = (" << currentPosition.x() << " | " << currentPosition.y() << " | " << currentPosition.z() << " )" << std::endl;
323 //std::cout << "ExtSteppingAction Current Position (x,y,z)*mm = (" << currentPosition.x()*mm << " | " << currentPosition.y()*mm << " | " << currentPosition.z()*mm << " )" << std::endl;
324 }
325
326 else if( newVolumeName.contains("_west_2")) //New MROC - Layer 2
327 {
328 newVolumeNumber = (help_mrpc_nb)*(-0.5)+18;
329
330 int partID_help =6;
331
332 int help = BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(currentPosition.x(),currentPosition.y(), partID_help, newVolumeNumber);
333 newVolumeNumber = BesTofDigitizerEcV4::Produce_unique_identifier(newVolumeNumber, help);
334 /*
335 std::cout << "ExtSteppingAction partID | ModulNumber " << partID_help << " | " << (help_mrpc_nb)*(-0.5)+18 << std::endl;
336 std::cout << "ExtSteppingAction Readoutstripnumber = " << help << std::endl;
337 std::cout << "ExtSteppingAction unique identifier = " << newVolumeNumber << std::endl;
338 if(help==0){std::cout << "ExtSteppingAction x | y" <<currentPosition.x()<< " | " << currentPosition.y() << std::endl;}
339 */
340 //std::cout << "ExtSteppingAction You are using the new MRPC - Layer 2. Modulenumber " << newVolumeNumber << std::endl;
341 //std::cout << "ExtSteppingAction Current Position (x,y,z) = (" << currentPosition.x()<< " | " << currentPosition.y() << " | " << currentPosition.z() << " )" << std::endl;
342 //std::cout << "ExtSteppingAction Current Position (x,y,z)*mm = (" << currentPosition.x()*mm << " | " << currentPosition.y()*mm << " | " << currentPosition.z()*mm << " )" << std::endl;
343 }
344 else if(newVolumeName.contains("_east_2") ) //New MROC - Layer 2
345 {
346 newVolumeNumber = (help_mrpc_nb)*(0.5)+1;
347
348 int partID_help =3;
349
350
351 int help = BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(currentPosition.x(),currentPosition.y(), partID_help, newVolumeNumber);
352 newVolumeNumber = BesTofDigitizerEcV4::Produce_unique_identifier(newVolumeNumber, help);
353 /*
354 std::cout << "ExtSteppingAction partID | ModulNumber " << partID_help << " | " << (help_mrpc_nb)*(0.5)+1 << std::endl;
355 std::cout << "ExtSteppingAction Readoutstripnumber = " << help << std::endl;
356 std::cout << "ExtSteppingAction unique identifier = " << newVolumeNumber << std::endl;
357 if(help==0){std::cout << "ExtSteppingAction x | y" <<currentPosition.x()<< " | " << currentPosition.y() << std::endl;}
358 */
359 //std::cout << "ExtSteppingAction You are using the new MRPC - Layer 2. Modulenumber " << newVolumeNumber << std::endl;
360 //std::cout << "ExtSteppingAction Current Position (x,y,z) = (" << currentPosition.x()<< " | " << currentPosition.y() << " | " << currentPosition.z() << " )" << std::endl;
361 //std::cout << "ExtSteppingAction Current Position (x,y,z)*mm = (" << currentPosition.x()*mm << " | " << currentPosition.y()*mm << " | " << currentPosition.z()*mm << " )" << std::endl;
362 }
363 else{ newVolumeNumber=(527-newVolumeNumber)/3;}
364
365
366
367 if(myExtTrack)
368 {
369 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
370 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
371 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
372 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
373 myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
374 myTof1Flag = true;
375 }
376 return;
377 }//close if (!myTof1Flag)
378
379
380
381
382 if( myInTof1 && ( oldVolumeName=="logicalScinBr1" || oldVolumeName.contains("ScinEc") ||
383 newVolumeName.contains("logical_sensitive_detector_east") || newVolumeName.contains("logical_sensitive_detector_west")) ) {
384 myInTof1 = false;
385 myOutTof1 = true;
386 myPathOutTof1 = currentTrack->GetTrackLength();
387 myPathInTof1.push_back(myPathOutTof1 - myPathIntoTof1);
388 if(msgFlag) cout << "myPathOutTof1 = " << myPathOutTof1 << endl;
389 }
390
391 if( myOutTof1 && ( newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ||
392 newVolumeName.contains("logical_sensitive_detector_east") || newVolumeName.contains("logical_sensitive_detector_west")) ) {
393 myInTof1 = true;
394 myOutTof1 = false;
395 myPathIntoTof1 = currentTrack->GetTrackLength();
396 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
397 }
398
399
400
401 //Get Tof layer2 Ext data
402 if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" )
403 {
404 double currentTrackLength = currentTrack->GetTrackLength();
405 double totalTrackLength = currentTrackLength+initialPath;
406 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
407 double localTime = currentTrack->GetLocalTime();
408 double totalTOF = localTime + initialTof;
409 newVolumeNumber=(527-theTouchable->GetReplicaNumber(2))/3;
410
411 myInTof2 = true;
412 myPathIntoTof2 = currentTrackLength;
413 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
414
415 if(myExtTrack)
416 {
417 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
418 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
419 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
420 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
421 myExtTrack->SetTof2Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
422 myTof2Flag = true;
423 }
424 return;
425 }//close if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" )
426
427 if( myInTof2 && oldVolumeName=="logicalScinBr2" ) {
428 myInTof2 = false;
429 myOutTof2 = true;
430 myPathOutTof2 = currentTrack->GetTrackLength();
431 myPathInTof2.push_back(myPathOutTof2 - myPathIntoTof2);
432 if(msgFlag) cout << "myPathOutTof2 = " << myPathOutTof2 << endl;
433 }
434
435 if( myOutTof2 && newVolumeName=="logicalScinBr2" ) {
436 myInTof2 = true;
437 myOutTof2 = false;
438 myPathIntoTof2 = currentTrack->GetTrackLength();
439 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
440 }
441
442 //Get PhyEmc data.
443 if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
444 {
445 newVolumeNumber = -2;
446 if(myExtTrack)
447 {
448 // myPathIntoCrystal = currentTrack->GetTrackLength();
449 Hep3Vector nPhi(-y/r,x/r,0.);
450 double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
451
452 Hep3Vector aPosition = currentPosition;
453 double R = aPosition.r();
454 double aTheta = aPosition.theta();
455 aPosition.setTheta(aTheta + M_PI_2);
456 double errorTheta;
457 errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
458 if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
459 myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
460 }
461 myPhyEmcFlag = true;
462 return;
463 }
464
465
466 if(!myEmcPathFlag && newVolumeName.contains("Crystal") )
467 {
468 myPathIntoCrystal = currentTrack->GetTrackLength();
469 // cout<<"Enter Crystal, path: "<<myPathIntoCrystal<<endl;
470 myEmcPathFlag = true;
471 }
472
473 //Get Emc Ext data.
474 if( (!myEmcFlag) && newVolumeName.contains("Crystal") )
475 {
476 if(myExtTrack)
477 {
478 //------------------- record crystal number
479 int npart,ntheta,nphi;
480 if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) { //barrel
481 npart=1;
482 std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName()).substr(16,2));
483 thetaBuf >> ntheta ;
484 nphi = 308-currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
485 } else { //endcap
486 npart=2-2*currentTrack->GetNextTouchableHandle()->GetCopyNumber(3);
487 int endSector=currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
488 int endNb,endNbGDML;
489 std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(0)->GetName()).substr(20,2));
490 thetaBuf >> endNb ;
491 endNbGDML=CalculateEmcEndCopyNb(endNb);
492 int endSectorGDML=CalculateEmcEndPhiNb(endSector);
493 CalculateEmcEndThetaPhi(npart,endSectorGDML,endNbGDML,ntheta,nphi);
494 }
495 ostringstream strEmc;
496 if(ntheta<10) {
497 strEmc<<"EmcPart"<<npart<<"Theta0"<<ntheta<<"Phi"<<nphi;
498 } else {
499 strEmc<<"EmcPart"<<npart<<"Theta"<<ntheta<<"Phi"<<nphi;
500 } //------------------------------------------
501
502 // myPathIntoCrystal = currentTrack->GetTrackLength();
503 Hep3Vector nPhi(-y/r,x/r,0.);
504 double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
505
506 Hep3Vector aPosition = currentPosition;
507 double R = aPosition.r();
508 double aTheta = aPosition.theta();
509 aPosition.setTheta(aTheta + M_PI_2);
510 double errorTheta;
511 errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
512 if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
513 myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,strEmc.str(),
514 //myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,
515 newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
516 }
517 myEmcFlag = true;
518 return;
519 }
520
521 //Get path in Emc
522 if(myEmcPathFlag && oldVolumeName.contains("Crystal") )
523 {
524 myPathOutCrystal = currentTrack->GetTrackLength();
525 myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
526 // cout<<"Go out of crystal, path: "<<myPathOutCrystal<<endl;
527 myEmcPathFlag=false;
528 }
529
530
531 //Get Muc Ext Data.
532 if( (!myMucFlag) && ( (fabs(x)>=myMucR) || (fabs(y)>=myMucR) || ((fabs(x-y)/sqrt(2.))>=myMucR) || ((fabs(x+y)/sqrt(2.))>=myMucR) || (fabs(z)>=myMucZ) ) )
533 {
534 int newVolumeNumber = newVolume->GetCopyNo();
535 if(myExtTrack)
536 {
537 Hep3Vector xVector(1.0,0,0);
538 Hep3Vector yVector(0,1.0,0);
539 Hep3Vector zVector(0,0,1.0);
540 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
541 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
542 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
543 double tzError;
544 double phi = currentPosition.phi();
545 if(phi<0.) phi+=M_PI;
546 int i = int(8.0*phi/M_PI);
547 if( i==0 || i==7 || i==8 )
548 {
549 tzError = yError;
550 }
551 if( i==1 || i==2 )
552 {
553 Hep3Vector tzVector(-1.0,1.0,0.);
554 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
555 }
556 if( i==3 || i==4 )
557 {
558 tzError = xError;
559 }
560 if( i==5 || i==6 )
561 {
562 Hep3Vector tzVector(1.0,1.0,0.);
563 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
564 }
565 myExtTrack->SetMucData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
566 }
567 myMucFlag = true;
568 if(!(ParticleName.contains("mu")&&myUseMucKalFlag))
569 {
570 //currentTrack->SetKineticEnergy(0.0);
571 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
572 m_trackstop=true;
573 return;
574 }
575 }
576
577 //**************************************************
578 //************inset KalFilter Algorithm in MUC******
579 //**************************************************
580 HepSymMatrix XpErr = extXpErr->get_err();
581 int namesize = oldVolumeName.size();
582 bool Flag1(false),Flag2(false),Flag3(false),Flag4(false),Flag5(false);
583 Flag1 = m_mucdigicol->size()>0;
584 Flag2 = myUseMucKalFlag;
585 Flag3 = ParticleName.contains("mu")&&oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&(oldVolumeName.contains("G")||oldVolumeName.contains("Ab"));
586 Flag4 = oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&((namesize==10&&oldVolumeName.contains("G"))||(namesize==11&&oldVolumeName.contains("Ab")));
587 Flag5 = !((RemID[0]==1&&RemID[2]==9)||((RemID[0]==0||RemID[0]==2)&&RemID[2]==8));
588 // if(!Flag5) {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
589 if(Flag1&&Flag2&&Flag3&&Flag5)
590 {
591 //get depth in Ab
592 double depth = currentStep-> GetStepLength();
593 if(oldVolumeName.contains("Ab"))
594 RemDepth = RemDepth+ depth;
595 if(RemStep==0&&Flag4)
596 {
597 Hep3Vector ID_1 = GetGapID(oldVolumeName);
598 RemID = ID_1;
599
600 bool LastLay = (ID_1[0]==1&&ID_1[2]<9)||((ID_1[0]==0||ID_1[0]==2)&&ID_1[2]<8);
601 if(RememberID[2]!=ID_1[2]&&LastLay)
602 {
603 Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(ID_1[0],ID_1[1],ID_1[2])->TransformToGap(currentPosition);
604 double dist = fabs(pos_loc.z());
605 RemPositon = currentPosition;
606 RemMomentum = currentMomentum;
607 RemXpErr = XpErr;
608 RemDist = dist;
609 RemStep++;
610 }
611 }
612 if(RemStep>0)
613 {
614 bool WillFilter = false;
615 bool LastLay_ = false;
616 double dist_b = 99999.;
617 Hep3Vector ID_2;
618 if(Flag4)
619 {
620 ID_2 = GetGapID(oldVolumeName);
621 LastLay_ =(ID_2[0]==1&&ID_2[2]<9)||((ID_2[0]==0||ID_2[0]==2)&&ID_2[2]<8);
622 if(LastLay_)dist_b=fabs(MucGeoGeneral::Instance()->GetGap(ID_2[0],ID_2[1],ID_2[2])->TransformToGap(currentPosition).z());
623 if(RemID!=ID_2)
624 WillFilter=true;
625
626 }
627
628 Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(RemID[0],RemID[1],RemID[2])->TransformToGap(currentPosition);
629 double dist = fabs(pos_loc.z());
630 //get the nearest point from the center of gap
631 if(!WillFilter&&RemDist>dist)
632 {
633 RemPositon = currentPosition;
634 RemMomentum = currentMomentum;
635 RemXpErr = XpErr;
636 RemDist = dist;
637 RemVol = oldVolumeName;
638 }
639
640 //if find the nearest point in the gap, Open Fillter
641 if(WillFilter)
642 {
643 ExtMucKal* muckal = new ExtMucKal();
644 muckal->SetGapID(RemID);
645 muckal->SetPosMomErr(RemPositon,RemMomentum,RemXpErr);
646 muckal->SetMucDigiCol(m_mucdigicol);
647 muckal->SetMucWindow(myMucWindow);
648 bool iniOK = muckal->MucKalIniti();
649 vector<MucRecHit*> tmp = muckal->TrackHit();
650 bool filter = muckal->ExtMucFilter();
651 // bool filter = muckal->GetFilterStatus();
652 double chi2 = muckal->GetChi2();
653 bool samestrip = muckal->GetSameStrip();
654 if(filter&&iniOK)
655 {
656
657 myMucnfit_++;
658 //cout<<"myMucnfit_: "<<myMucnfit_<<endl;
659 myMucchisq_ = myMucchisq_+chi2;
660 muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
661 RememberID = muckal->GetHitGap();
662 if(RememberID[0]==1)myMucbrLastLay_=RememberID[2];
663 if(RememberID[0]==0||RememberID[0]==2)myMucecLastLay_=RememberID[2];
664 if(oldVolumeName.contains("Ab"))
665 RemDepth = RemDepth-depth;
666 if(myMucnfit_==1)
667 myMucdepth_ = RemDepth;
668 else
669 myMucdepth_=myMucdepth_+RemDepth;
670
671 if(!samestrip)
672 {
673
674 extXpErr->put_err(m_err_mod);
675 extXpErr->set_pos(m_pos_mod);
676 extXpErr->set_mom(m_mom_mod);
677 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
678
679 }
680 else
681 {
682 //cout<<"-------hit exist, but no modify-------"<<endl;
683 RemPositon = currentPosition;
684 RemMomentum = currentMomentum;
685 RemXpErr = XpErr;
686 RemDist = dist_b;
687 RemID =ID_2;
688 if(oldVolumeName.contains("Ab"))
689 RemDepth = depth;
690 else
691 RemDepth = 0.;
692 }
693 if(msgFlag)cout<<"---------Filter is OK---------- "<<endl;
694 }
695 else
696 {
697 // if(LastLay_)
698 //{
699 RemPositon = currentPosition;
700 RemMomentum = currentMomentum;
701 RemXpErr = XpErr;
702 RemDist = dist_b;
703 RemID =ID_2;
704 // }
705 }
706 delete muckal;
707 }
708 }
709 if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
710 }
711
712 /*
713 if(Flag1&&Flag2&&Flag3)
714 {
715 ExtMucKal* muckal = new ExtMucKal();
716 muckal->SetVolume(oldVolume,newVolume);
717 muckal->SetMomPosErr(currentPosition,currentMomentum,XpErr);
718 muckal->SetMucDigiCol(mucdigicol);
719 muckal->ExtMucFilter();
720 muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
721 bool filter = muckal->GetFilterStatus();
722 double chi2 = muckal->GetChi2();
723 if(filter)
724 {
725 myMucnfit_++;
726 myMucchisq_ = myMucchisq_+chi2;
727 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
728 RememberVol.assign(oldVolumeName,10);
729 extXpErr->put_err(m_err_mod);
730 extXpErr->set_pos(m_pos_mod);
731 extXpErr->set_mom(m_mom_mod);
732 }
733 delete muckal;
734 }
735 */
736 // if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
737 /* //Get Muc Ext Hits Data.
738 if(newVolumeName.contains("Strip"))
739 {
740 int newVolumeNumber = newVolume->GetCopyNo();
741 if(myExtTrack)
742 {
743 ExtMucHit aExtMucHit;
744 Hep3Vector xVector(1.0,0,0);
745 Hep3Vector yVector(0,1.0,0);
746 Hep3Vector zVector(0,0,1.0);
747 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
748 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
749 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
750 double tzError;
751 double phi = currentPosition.phi();
752 if(phi<0.) phi+=M_PI;
753 int i = int(8.0*phi/M_PI);
754 if( i==0 || i==7 || i==8 )
755 {
756 tzError = yError;
757 }
758 if( i==1 || i==2 )
759 {
760 Hep3Vector tzVector(-1.0,1.0,0.);
761 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
762 }
763 if( i==3 || i==4 )
764 {
765 tzError = xError;
766 }
767 if( i==5 || i==6 )
768 {
769 Hep3Vector tzVector(1.0,1.0,0.);
770 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
771 }
772 aExtMucHit.SetExtMucHit(currentPosition,currentMomentum,newVolumeName,newVolumeNumber,extXpErr->get_err(),zError,tzError,xError,yError);
773 myExtTrack->AddExtMucHit(aExtMucHit);
774 }
775 }
776 */
777 }
778 else {
779 if(msgFlag) cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<endl;
780 double x = currentPosition.x();
781 double y = currentPosition.y();
782 double z = currentPosition.z();
783 if( (fabs(x)>=2*myMucR) || (fabs(y)>=2*myMucR) || (fabs(z)>=2*myMucZ) )
784 //currentTrack->SetKineticEnergy(0.0);// protection in case that a very energetic track travels without touching BESIII
785 {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
786 }
787
788 }
789 else if(currentTrackStatus == fStopAndKill)
790 {
791 m_trackstop=true;
792 if(myEmcFlag) myExtTrack->SetEmcPath(myPathInCrystal/10.);
793 if(myTof1Flag) myExtTrack->setPathInTof1(myPathInTof1);
794 if(myTof2Flag) myExtTrack->setPathInTof2(myPathInTof2);
795 if(msgFlag) {
796 cout << "myPathInTof1 = " ;
797 for(int i=0; i!=myPathInTof1.size(); i++)
798 cout << myPathInTof1[i] << " ";
799 cout << endl;
800 cout << "myPathInTof2 = " ;
801 for(int i=0; i!=myPathInTof2.size(); i++)
802 cout << myPathInTof2[i] << " ";
803 cout << endl;
804 }
805
806 if(msgFlag)
807 {
808 if(newVolume!=0)
809 {
810 std::cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//stoped"<<endl;
811 cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
812 }
813 else {
814 cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//escaped"<<std::endl;
815 std::cout<<"Error matrix is:"<<extXpErr->get_err()<<std::endl;
816 }
817 }
818 }
819}
820
821
822
823void ExtSteppingAction::CalculateChicc(G4Material* currentMaterial)
824{
825 int n = currentMaterial->GetNumberOfElements();
826 const double *p = currentMaterial->GetFractionVector();
827 double density = (currentMaterial->GetDensity())/(g/cm3);
828 double temp=0.0;
829 for(int i=0;i<n;i++)
830 {
831 const G4Element* mElement = currentMaterial->GetElement(i);
832 double z = mElement->GetZ();
833 double a = mElement->GetN();
834 // std::cout<<"Fraction: "<<*p<<" z: "<<z<<" a: "<<a<<std::endl;
835 temp += *(p++)/a*z*(z+1);
836 }
837 chicc = 0.39612e-3*std::sqrt(density*temp);
838 // std::cout<<"chicc:"<<chicc<<std::endl;
839}
840
841
842HepSymMatrix & ExtSteppingAction::myOutputSymMatrix(const HepSymMatrix & m1)
843{
844 // std::cout<<"myOutputSymMatrix:1"<<std::endl;
845 HepSymMatrix m;
846 // std::cout<<"myOutputSymMatrix:2"<<std::endl;
847 m=m1;
848 // std::cout<<"myOutputSymMatrix:3"<<std::endl;
849 for(int i=0;i<6;i++)
850 { for(int j=0;j<=i;j++)
851 { if(i<3) {
852 m[i][j]=m[i][j]/100; //mm*mm --> cm*cm
853 }
854 else if(j<3) {
855 m[i][j]=m[i][j]/10000; //mm*MeV --> cm*GeV
856 }
857 else {
858 m[i][j]=m[i][j]/1000000; //MeV*MeV --> GeV*GeV
859 }
860 }
861 }
862 // std::cout<<"myOutputSymMatrix:4"<<std::endl;
863 // std::cout<<"m1="<<m1<<std::endl;
864 // std::cout<<"m="<<m<<std::endl;
865 myOutputSM=m;
866 return myOutputSM;
867}
868
869void ExtSteppingAction::CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
870{
871 if((sector>=0)&&(sector<3))
872 sector+=16;
873 //if((sector!=7)&&(sector!=15))
874 {
875 if((nb>=0)&&(nb<4))
876 {
877 ntheta = 0;
878 nphi = (sector-3)*4+nb;
879 }
880 else if((nb>=4)&&(nb<8))
881 {
882 ntheta = 1;
883 nphi = (sector-3)*4+nb-4;
884 }
885 else if((nb>=8)&&(nb<13))
886 {
887 ntheta = 2;
888 nphi = (sector-3)*5+nb-8;
889 }
890 else if((nb>=13)&&(nb<18))
891 {
892 ntheta = 3;
893 nphi = (sector-3)*5+nb-13;
894 }
895 else if((nb>=18)&&(nb<24))
896 {
897 ntheta = 4;
898 nphi = (sector-3)*6+nb-18;
899 }
900 else if((nb>=24)&&(nb<30))
901 {
902 ntheta = 5;
903 nphi = (sector-3)*6+nb-24;
904 }
905 }
906
907 if(npart==2)
908 {
909 switch(ntheta)
910 {
911 case 0:
912 if(nphi<32)
913 nphi = 31-nphi;
914 else
915 nphi = 95-nphi;
916 break;
917 case 1:
918 if(nphi<32)
919 nphi = 31-nphi;
920 else
921 nphi = 95-nphi;
922 break;
923 case 2:
924 if(nphi<40)
925 nphi = 39-nphi;
926 else
927 nphi = 119-nphi;
928 break;
929 case 3:
930 if(nphi<40)
931 nphi = 39-nphi;
932 else
933 nphi = 119-nphi;
934 break;
935 case 4:
936 if(nphi<48)
937 nphi = 47-nphi;
938 else
939 nphi = 143-nphi;
940 break;
941 case 5:
942 if(nphi<48)
943 nphi = 47-nphi;
944 else
945 nphi = 143-nphi;
946 default:
947 ;
948 }
949 }
950}
951
952int ExtSteppingAction::CalculateEmcEndPhiNb(int num)
953{
954 if(num==0||num==1) {
955 return 15-num*8;
956 } else if(num==2||num==3) {
957 return 30-num*8;
958 } else if(num<=9) {
959 return 17-num;
960 } else {
961 return 15-num;
962 }
963}
964
965int ExtSteppingAction::CalculateEmcEndCopyNb(int num)
966{
967 int copyNb;
968 switch(num){
969 case 30:
970 copyNb = 5;
971 break;
972 case 31:
973 copyNb = 6;
974 break;
975 case 32:
976 copyNb = 14;
977 break;
978 case 33:
979 copyNb = 15;
980 break;
981 case 34:
982 copyNb = 16;
983 break;
984 default:
985 copyNb = num;
986 break;
987 }
988 return copyNb;
989}
990
991void ExtSteppingAction::InfmodMuc(Hep3Vector &pos,Hep3Vector &mom,HepSymMatrix &err)
992{
993 pos = m_pos_mod;
994 mom = m_mom_mod;
995 err = m_err_mod;
996}
997
998
999Hep3Vector ExtSteppingAction::GetGapID(G4String vol)
1000{
1001 int par(-1),se(-1),ga(-1);
1002 G4String strPart = vol.substr(5,1);
1003 //G4String strSeg = m_MotherVolumeName.substr(7,1);
1004
1005 G4String strSeg = vol.substr(7,1);
1006 G4String strGap;
1007 if(vol.contains("G")) strGap= vol.substr(9,1);
1008 if(vol.contains("Ab")) strGap= vol.substr(10,1);
1009 std::istrstream partBuf(strPart.c_str(), strlen(strPart.c_str()));
1010 std::istrstream segBuf(strSeg.c_str(), strlen(strSeg.c_str()));
1011 std::istrstream gapBuf(strGap.c_str(), strlen(strGap.c_str()));
1012
1013 partBuf >> par ;
1014 segBuf >> se;
1015 gapBuf >> ga;
1016 if(vol.contains("Ab")&&par==1) ga = ga+1;
1017 //panelBuf >> m_panel;
1018 Hep3Vector vec;
1019 vec[0]=par;
1020 vec[1]=se;
1021 vec[2]=ga;
1022 return vec;
1023}
1024
1025
nPhi
const int nPhi
Definition: BesTofDigitizerEcV3.hh:37
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::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::GetGapID
Hep3Vector GetGapID(G4String vol)
Definition: ExtSteppingAction.cxx:999
ExtSteppingAction::CalculateEmcEndThetaPhi
void CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
Definition: ExtSteppingAction.cxx:869
ExtSteppingAction::InfmodMuc
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtSteppingAction.cxx:991
ExtSteppingAction::UserSteppingAction
void UserSteppingAction(const G4Step *currentStep)
Definition: ExtSteppingAction.cxx:102
ExtSteppingAction::CalculateEmcEndCopyNb
int CalculateEmcEndCopyNb(int num)
Definition: ExtSteppingAction.cxx:965
ExtSteppingAction::CalculateEmcEndPhiNb
int CalculateEmcEndPhiNb(int num)
Definition: ExtSteppingAction.cxx:952
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