BOSS 7.0.5
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{
92 RemStep =0;
93 RemDist = 99999.;
94 RemDepth = 0.;
95 RemID[0]=-1;
96 RemID[1]=-1;
97 RemID[2]=-1;
98 RemVol = "";
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 G4String name1;
218 if(newVolumeName=="logical_gasLayer")
219 {
220 name1=theTouchable->GetVolume(3)->GetName();
221 newVolumeNumber=theTouchable->GetReplicaNumber(3);
222 }
223
224
225
226
227
228
229 //Get PhyTof data
230
231 //std::cout << "ExtSteppingAction NewVolumeName: " << newVolumeName << std::endl;
232 //std::cout << "ExtSteppingAction OldVolumeName: " << oldVolumeName << std::endl;
233 //std::cout << "ExtSteppingAction name1: " <<name1<<std::endl;
234 //Comment this part to remove the airtight Tof. Most hits would be replaced by the following specific sensitive
235 //detector. The remianing few hits are abondoned in Tof reconstruction.
236 /*
237 if( (!myTofFlag) && (!myTof1Flag) && (newVolumeName.contains("Tof") ))
238 {
239 newVolumeNumber = -2;
240 double currentTrackLength = currentTrack->GetTrackLength();
241 double totalTrackLength = currentTrackLength + initialPath;
242 //double initialTof = initialPath/(myBetaInMDC*299.792458);
243 //cout<<"initialTof = "<<initialTof<<endl;
244 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
245 double localTime = currentTrack->GetLocalTime();
246 //cout<<"localTime = "<<localTime<<endl;
247 double totalTOF = localTime + initialTof;
248 //cout<<"totalTOF= "<<totalTOF<<endl;
249
250 //std::cout << "ExtSteppingAction Volumename contains one of the marker words" << std::endl;
251
252 if(myExtTrack)
253 {
254 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
255 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
256 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
257 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
258 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);
259 myTofFlag = true;
260 }
261 return;
262 }//close if (!myTofFlag)
263 */
264
265
266
267 //Get Tof layer1 Ext data
268 if((!myTof1Flag) && (newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ||
269 (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))))
270 {
271 double currentTrackLength = currentTrack->GetTrackLength();
272 double totalTrackLength = currentTrackLength+initialPath;
273 double localTime = currentTrack->GetLocalTime();
274 double totalTOF = localTime + initialTof;
275 myInTof1 = true;
276 myPathIntoTof1 = currentTrackLength;
277 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
278
279 newVolumeNumber=theTouchable->GetReplicaNumber(2);
280
281 if(newVolumeName.contains("ScinEc")) {
282 newVolumeNumber=(95-newVolumeNumber)/2;
283 newVolumeNumber=newVolumeNumber+176;
284
285 if(newVolumeName.contains("East")) newVolumeNumber=newVolumeNumber+48;
286 if(myExtTrack) {
287 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
288 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
289 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
290 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
291 myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
292 myTof1Flag = true;
293 }
294 }
295
296 else if(newVolumeName=="logical_gasLayer")
297 {
298 G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
299 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(currentPostPosition);
300 //cout<<"localPosition: "<<localPosition.x()<<" "<<localPosition.y()<<" "<<localPosition.z()<<endl;
301
302 //Here we assume the spread of one strip is 24+3=27mm.
303 //Distance between strips and lower large glass margin: 6, top: 8, z_strip=z_small_glass-0.5
304 double z_mm = localPosition.z()-0.5+(24+3)*6;
305 int strip;
306 if(z_mm<=0)
307 {
308 strip=0;
309 }
310 else if(z_mm>0 && z_mm<12*27)
311 {
312 strip=floor(z_mm/27);
313 }
314 else
315 {
316 strip=11;
317 }
318 if(strip<0) strip=0;
319 if(strip>11) strip=11;
320
321 newVolumeNumber=theTouchable->GetReplicaNumber(3);
322 newVolumeNumber = 35-newVolumeNumber;
323 if( currentPosition.z() < 0.0 ) { newVolumeNumber = newVolumeNumber + 36; }
324 newVolumeNumber = 12*newVolumeNumber+strip;
325 newVolumeNumber = newVolumeNumber + 176 + 96;
326 if(myExtTrack) {
327 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
328 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
329 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
330 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
331 Hep3Vector locP = Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
332 //cout<<"locP: "<<locP.x()<<" "<<locP.y()<<" "<<locP.z()
333 //<<" currentMomentum: "<<currentMomentum.x()<<" "<<currentMomentum.y()<<" "<<currentMomentum.z()
334 //<<" newVolumeName: "<<newVolumeName<<" newVolumeNumber: "<<newVolumeNumber
335 //<<" totalTOF: "<<totalTOF<<" totalTrackLength: "<<totalTrackLength<<endl;
336 myExtTrack->SetTof1Data(locP/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
337 myTof1Flag = true;
338 }
339 }
340 else
341 {
342 newVolumeNumber=(527-newVolumeNumber)/3;
343 if(myExtTrack) {
344 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
345 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
346 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
347 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
348 myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
349 myTof1Flag = true;
350 }
351 }
352
353 return;
354
355 }//close if (!myTof1Flag)
356
357 if( myInTof1 && ( oldVolumeName=="logicalScinBr1" || oldVolumeName.contains("ScinEc") ||
358 (oldVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))) ) {
359 myInTof1 = false;
360 myOutTof1 = true;
361 myPathOutTof1 = currentTrack->GetTrackLength();
362 myPathInTof1.push_back(myPathOutTof1 - myPathIntoTof1);
363 if(msgFlag) cout << "myPathOutTof1 = " << myPathOutTof1 << endl;
364 }
365
366 if( myOutTof1 && ( newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ||
367 (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))) ) {
368 myInTof1 = true;
369 myOutTof1 = false;
370 myPathIntoTof1 = currentTrack->GetTrackLength();
371 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
372 }
373
374
375 //Get Tof layer2 Ext data
376 if( (!myTof2Flag) && (newVolumeName=="logicalScinBr2" ||
377 (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3"))))
378 {
379 double currentTrackLength = currentTrack->GetTrackLength();
380 double totalTrackLength = currentTrackLength+initialPath;
381 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
382 double localTime = currentTrack->GetLocalTime();
383 double totalTOF = localTime + initialTof;
384 myInTof2 = true;
385 myPathIntoTof2 = currentTrackLength;
386 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
387 newVolumeNumber=theTouchable->GetReplicaNumber(2);
388
389 if(newVolumeName=="logical_gasLayer")
390 {
391 G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
392 G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(currentPostPosition);
393 //cout<<"localPosition: "<<localPosition.x()<<" "<<localPosition.y()<<" "<<localPosition.z()<<endl;
394
395 //Here we assume the spread of one strip is 24+3=27mm.
396 //Distance between strips and lower large glass margin: 6, top: 8, z_strip=z_small_glass-0.5
397 double z_mm = localPosition.z()-0.5+(24+3)*6;
398 int strip;
399 if(z_mm<=0)
400 {
401 strip=0;
402 }
403 else if(z_mm>0 && z_mm<12*27)
404 {
405 strip=floor(z_mm/27);
406 }
407 else
408 {
409 strip=11;
410 }
411 if(strip<0) strip=0;
412 if(strip>11) strip=11;
413
414 newVolumeNumber=theTouchable->GetReplicaNumber(3);
415 newVolumeNumber = 35-newVolumeNumber;
416 if( currentPosition.z() < 0.0 ) { newVolumeNumber = newVolumeNumber + 36; }
417 newVolumeNumber = 12*newVolumeNumber+strip;
418 newVolumeNumber = newVolumeNumber + 176 + 96;
419
420 if(myExtTrack)
421 {
422 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
423 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
424 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
425 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
426 Hep3Vector locP = Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
427 //cout<<"locP: "<<locP.x()<<" "<<locP.y()<<" "<<locP.z()
428 //<<" currentMomentum: "<<currentMomentum.x()<<" "<<currentMomentum.y()<<" "<<currentMomentum.z()
429 //<<" newVolumeName: "<<newVolumeName<<" newVolumeNumber: "<<newVolumeNumber
430 //<<" totalTOF: "<<totalTOF<<" totalTrackLength: "<<totalTrackLength<<endl;
431 myExtTrack->SetTof2Data(locP/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
432 myTof2Flag = true;
433 }
434 }
435 else
436 {
437 newVolumeNumber=(527-newVolumeNumber)/3;
438 if(myExtTrack)
439 {
440 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
441 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
442 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
443 double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
444 myExtTrack->SetTof2Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
445 myTof2Flag = true;
446 }
447 }
448 return;
449 }//close if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" )
450
451 if( myInTof2 && (oldVolumeName=="logicalScinBr2" ||
452 (oldVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3"))) ) {
453 myInTof2 = false;
454 myOutTof2 = true;
455 myPathOutTof2 = currentTrack->GetTrackLength();
456 myPathInTof2.push_back(myPathOutTof2 - myPathIntoTof2);
457 if(msgFlag) cout << "myPathOutTof2 = " << myPathOutTof2 << endl;
458 }
459
460 if( myOutTof2 && (newVolumeName=="logicalScinBr2" ||
461 (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3"))) ) {
462 myInTof2 = true;
463 myOutTof2 = false;
464 myPathIntoTof2 = currentTrack->GetTrackLength();
465 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
466 }
467
468
469 //Get PhyEmc data.
470 if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
471 {
472 newVolumeNumber = -2;
473 if(myExtTrack)
474 {
475 // myPathIntoCrystal = currentTrack->GetTrackLength();
476 Hep3Vector nPhi(-y/r,x/r,0.);
477 double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
478
479 Hep3Vector aPosition = currentPosition;
480 double R = aPosition.r();
481 double aTheta = aPosition.theta();
482 aPosition.setTheta(aTheta + M_PI_2);
483 double errorTheta;
484 errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
485 if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
486 myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
487 }
488 myPhyEmcFlag = true;
489 return;
490 }
491
492
493 if(!myEmcPathFlag && newVolumeName.contains("Crystal") )
494 {
495 myPathIntoCrystal = currentTrack->GetTrackLength();
496 // cout<<"Enter Crystal, path: "<<myPathIntoCrystal<<endl;
497 myEmcPathFlag = true;
498 }
499
500 //Get Emc Ext data.
501 if( (!myEmcFlag) && newVolumeName.contains("Crystal") )
502 {
503 if(myExtTrack)
504 {
505 //------------------- record crystal number
506 int npart,ntheta,nphi;
507 if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) { //barrel
508 npart=1;
509 std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName()).substr(16,2));
510 thetaBuf >> ntheta ;
511 nphi = 308-currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
512 } else { //endcap
513 npart=2-2*currentTrack->GetNextTouchableHandle()->GetCopyNumber(3);
514 int endSector=currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
515 int endNb,endNbGDML;
516 std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(0)->GetName()).substr(20,2));
517 thetaBuf >> endNb ;
518 endNbGDML=CalculateEmcEndCopyNb(endNb);
519 int endSectorGDML=CalculateEmcEndPhiNb(endSector);
520 CalculateEmcEndThetaPhi(npart,endSectorGDML,endNbGDML,ntheta,nphi);
521 }
522 ostringstream strEmc;
523 if(ntheta<10) {
524 strEmc<<"EmcPart"<<npart<<"Theta0"<<ntheta<<"Phi"<<nphi;
525 } else {
526 strEmc<<"EmcPart"<<npart<<"Theta"<<ntheta<<"Phi"<<nphi;
527 } //------------------------------------------
528
529 // myPathIntoCrystal = currentTrack->GetTrackLength();
530 Hep3Vector nPhi(-y/r,x/r,0.);
531 double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
532
533 Hep3Vector aPosition = currentPosition;
534 double R = aPosition.r();
535 double aTheta = aPosition.theta();
536 aPosition.setTheta(aTheta + M_PI_2);
537 double errorTheta;
538 errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
539 if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
540 myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,strEmc.str(),
541 //myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,
542 newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
543 }
544 myEmcFlag = true;
545 return;
546 }
547
548 //Get path in Emc
549 if(myEmcPathFlag && oldVolumeName.contains("Crystal") )
550 {
551 myPathOutCrystal = currentTrack->GetTrackLength();
552 myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
553 // cout<<"Go out of crystal, path: "<<myPathOutCrystal<<endl;
554 myEmcPathFlag=false;
555 }
556
557
558 //Get Muc Ext Data.
559 if( (!myMucFlag) && ( (fabs(x)>=myMucR) || (fabs(y)>=myMucR) || ((fabs(x-y)/sqrt(2.))>=myMucR) || ((fabs(x+y)/sqrt(2.))>=myMucR) || (fabs(z)>=myMucZ) ) )
560 {
561 int newVolumeNumber = newVolume->GetCopyNo();
562 if(myExtTrack)
563 {
564 Hep3Vector xVector(1.0,0,0);
565 Hep3Vector yVector(0,1.0,0);
566 Hep3Vector zVector(0,0,1.0);
567 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
568 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
569 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
570 double tzError;
571 double phi = currentPosition.phi();
572 if(phi<0.) phi+=M_PI;
573 int i = int(8.0*phi/M_PI);
574 if( i==0 || i==7 || i==8 )
575 {
576 tzError = yError;
577 }
578 if( i==1 || i==2 )
579 {
580 Hep3Vector tzVector(-1.0,1.0,0.);
581 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
582 }
583 if( i==3 || i==4 )
584 {
585 tzError = xError;
586 }
587 if( i==5 || i==6 )
588 {
589 Hep3Vector tzVector(1.0,1.0,0.);
590 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
591 }
592 myExtTrack->SetMucData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
593 }
594 myMucFlag = true;
595 if(!(ParticleName.contains("mu")&&myUseMucKalFlag))
596 {
597 //currentTrack->SetKineticEnergy(0.0);
598 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
599 m_trackstop=true;
600 return;
601 }
602 }
603
604 //**************************************************
605 //************inset KalFilter Algorithm in MUC******
606 //**************************************************
607 HepSymMatrix XpErr = extXpErr->get_err();
608 int namesize = oldVolumeName.size();
609 bool Flag1(false),Flag2(false),Flag3(false),Flag4(false),Flag5(false);
610 Flag1 = m_mucdigicol->size()>0;
611 Flag2 = myUseMucKalFlag;
612 Flag3 = ParticleName.contains("mu")&&oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&(oldVolumeName.contains("G")||oldVolumeName.contains("Ab"));
613 Flag4 = oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&((namesize==10&&oldVolumeName.contains("G"))||(namesize==11&&oldVolumeName.contains("Ab")));
614 Flag5 = !((RemID[0]==1&&RemID[2]==9)||((RemID[0]==0||RemID[0]==2)&&RemID[2]==8));
615 // if(!Flag5) {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
616 if(Flag1&&Flag2&&Flag3&&Flag5)
617 {
618 //get depth in Ab
619 double depth = currentStep-> GetStepLength();
620 if(oldVolumeName.contains("Ab"))
621 RemDepth = RemDepth+ depth;
622 if(RemStep==0&&Flag4)
623 {
624 Hep3Vector ID_1 = GetGapID(oldVolumeName);
625 RemID = ID_1;
626
627 bool LastLay = (ID_1[0]==1&&ID_1[2]<9)||((ID_1[0]==0||ID_1[0]==2)&&ID_1[2]<8);
628 if(RememberID[2]!=ID_1[2]&&LastLay)
629 {
630 Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(ID_1[0],ID_1[1],ID_1[2])->TransformToGap(currentPosition);
631 double dist = fabs(pos_loc.z());
632 RemPositon = currentPosition;
633 RemMomentum = currentMomentum;
634 RemXpErr = XpErr;
635 RemDist = dist;
636 RemStep++;
637 }
638 }
639 if(RemStep>0)
640 {
641 bool WillFilter = false;
642 bool LastLay_ = false;
643 double dist_b = 99999.;
644 Hep3Vector ID_2;
645 if(Flag4)
646 {
647 ID_2 = GetGapID(oldVolumeName);
648 LastLay_ =(ID_2[0]==1&&ID_2[2]<9)||((ID_2[0]==0||ID_2[0]==2)&&ID_2[2]<8);
649 if(LastLay_)dist_b=fabs(MucGeoGeneral::Instance()->GetGap(ID_2[0],ID_2[1],ID_2[2])->TransformToGap(currentPosition).z());
650 if(RemID!=ID_2)
651 WillFilter=true;
652
653 }
654
655 Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(RemID[0],RemID[1],RemID[2])->TransformToGap(currentPosition);
656 double dist = fabs(pos_loc.z());
657 //get the nearest point from the center of gap
658 if(!WillFilter&&RemDist>dist)
659 {
660 RemPositon = currentPosition;
661 RemMomentum = currentMomentum;
662 RemXpErr = XpErr;
663 RemDist = dist;
664 RemVol = oldVolumeName;
665 }
666
667 //if find the nearest point in the gap, Open Fillter
668 if(WillFilter)
669 {
670 ExtMucKal* muckal = new ExtMucKal();
671 muckal->SetGapID(RemID);
672 muckal->SetPosMomErr(RemPositon,RemMomentum,RemXpErr);
673 muckal->SetMucDigiCol(m_mucdigicol);
674 muckal->SetMucWindow(myMucWindow);
675 bool iniOK = muckal->MucKalIniti();
676 vector<MucRecHit*> tmp = muckal->TrackHit();
677 bool filter = muckal->ExtMucFilter();
678 // bool filter = muckal->GetFilterStatus();
679 double chi2 = muckal->GetChi2();
680 bool samestrip = muckal->GetSameStrip();
681 if(filter&&iniOK)
682 {
683
684 myMucnfit_++;
685 //cout<<"myMucnfit_: "<<myMucnfit_<<endl;
686 myMucchisq_ = myMucchisq_+chi2;
687 muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
688 RememberID = muckal->GetHitGap();
689 if(RememberID[0]==1)myMucbrLastLay_=RememberID[2];
690 if(RememberID[0]==0||RememberID[0]==2)myMucecLastLay_=RememberID[2];
691 if(oldVolumeName.contains("Ab"))
692 RemDepth = RemDepth-depth;
693 if(myMucnfit_==1)
694 myMucdepth_ = RemDepth;
695 else
696 myMucdepth_=myMucdepth_+RemDepth;
697
698 if(!samestrip)
699 {
700
701 extXpErr->put_err(m_err_mod);
702 extXpErr->set_pos(m_pos_mod);
703 extXpErr->set_mom(m_mom_mod);
704 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
705
706 }
707 else
708 {
709 //cout<<"-------hit exist, but no modify-------"<<endl;
710 RemPositon = currentPosition;
711 RemMomentum = currentMomentum;
712 RemXpErr = XpErr;
713 RemDist = dist_b;
714 RemID =ID_2;
715 if(oldVolumeName.contains("Ab"))
716 RemDepth = depth;
717 else
718 RemDepth = 0.;
719 }
720 if(msgFlag)cout<<"---------Filter is OK---------- "<<endl;
721 }
722 else
723 {
724 // if(LastLay_)
725 //{
726 RemPositon = currentPosition;
727 RemMomentum = currentMomentum;
728 RemXpErr = XpErr;
729 RemDist = dist_b;
730 RemID =ID_2;
731 // }
732 }
733 delete muckal;
734 }
735 }
736 if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
737 }
738
739 /*
740 if(Flag1&&Flag2&&Flag3)
741 {
742 ExtMucKal* muckal = new ExtMucKal();
743 muckal->SetVolume(oldVolume,newVolume);
744 muckal->SetMomPosErr(currentPosition,currentMomentum,XpErr);
745 muckal->SetMucDigiCol(mucdigicol);
746 muckal->ExtMucFilter();
747 muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
748 bool filter = muckal->GetFilterStatus();
749 double chi2 = muckal->GetChi2();
750 if(filter)
751 {
752 myMucnfit_++;
753 myMucchisq_ = myMucchisq_+chi2;
754 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
755 RememberVol.assign(oldVolumeName,10);
756 extXpErr->put_err(m_err_mod);
757 extXpErr->set_pos(m_pos_mod);
758 extXpErr->set_mom(m_mom_mod);
759 }
760 delete muckal;
761 }
762 */
763 // if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
764 /* //Get Muc Ext Hits Data.
765 if(newVolumeName.contains("Strip"))
766 {
767 int newVolumeNumber = newVolume->GetCopyNo();
768 if(myExtTrack)
769 {
770 ExtMucHit aExtMucHit;
771 Hep3Vector xVector(1.0,0,0);
772 Hep3Vector yVector(0,1.0,0);
773 Hep3Vector zVector(0,0,1.0);
774 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
775 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
776 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
777 double tzError;
778 double phi = currentPosition.phi();
779 if(phi<0.) phi+=M_PI;
780 int i = int(8.0*phi/M_PI);
781 if( i==0 || i==7 || i==8 )
782 {
783 tzError = yError;
784 }
785 if( i==1 || i==2 )
786 {
787 Hep3Vector tzVector(-1.0,1.0,0.);
788 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
789 }
790 if( i==3 || i==4 )
791 {
792 tzError = xError;
793 }
794 if( i==5 || i==6 )
795 {
796 Hep3Vector tzVector(1.0,1.0,0.);
797 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
798 }
799 aExtMucHit.SetExtMucHit(currentPosition,currentMomentum,newVolumeName,newVolumeNumber,extXpErr->get_err(),zError,tzError,xError,yError);
800 myExtTrack->AddExtMucHit(aExtMucHit);
801 }
802 }
803 */
804 }
805 else {
806 if(msgFlag) cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<endl;
807 double x = currentPosition.x();
808 double y = currentPosition.y();
809 double z = currentPosition.z();
810 if( (fabs(x)>=2*myMucR) || (fabs(y)>=2*myMucR) || (fabs(z)>=2*myMucZ) )
811 //currentTrack->SetKineticEnergy(0.0);// protection in case that a very energetic track travels without touching BESIII
812 {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
813 }
814
815 }
816 else if(currentTrackStatus == fStopAndKill)
817 {
818 m_trackstop=true;
819 if(myEmcFlag) myExtTrack->SetEmcPath(myPathInCrystal/10.);
820 if(myTof1Flag) myExtTrack->setPathInTof1(myPathInTof1);
821 if(myTof2Flag) myExtTrack->setPathInTof2(myPathInTof2);
822 if(msgFlag) {
823 cout << "myPathInTof1 = " ;
824 for(int i=0; i!=myPathInTof1.size(); i++)
825 cout << myPathInTof1[i] << " ";
826 cout << endl;
827 cout << "myPathInTof2 = " ;
828 for(int i=0; i!=myPathInTof2.size(); i++)
829 cout << myPathInTof2[i] << " ";
830 cout << endl;
831 }
832
833 if(msgFlag)
834 {
835 if(newVolume!=0)
836 {
837 std::cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//stoped"<<endl;
838 cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
839 }
840 else {
841 cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//escaped"<<std::endl;
842 std::cout<<"Error matrix is:"<<extXpErr->get_err()<<std::endl;
843 }
844 }
845 }
846}
847
848
849
850void ExtSteppingAction::CalculateChicc(G4Material* currentMaterial)
851{
852 int n = currentMaterial->GetNumberOfElements();
853 const double *p = currentMaterial->GetFractionVector();
854 double density = (currentMaterial->GetDensity())/(g/cm3);
855 double temp=0.0;
856 for(int i=0;i<n;i++)
857 {
858 const G4Element* mElement = currentMaterial->GetElement(i);
859 double z = mElement->GetZ();
860 double a = mElement->GetN();
861 // std::cout<<"Fraction: "<<*p<<" z: "<<z<<" a: "<<a<<std::endl;
862 temp += *(p++)/a*z*(z+1);
863 }
864 chicc = 0.39612e-3*std::sqrt(density*temp);
865 // std::cout<<"chicc:"<<chicc<<std::endl;
866}
867
868
869HepSymMatrix & ExtSteppingAction::myOutputSymMatrix(const HepSymMatrix & m1)
870{
871 // std::cout<<"myOutputSymMatrix:1"<<std::endl;
872 HepSymMatrix m;
873 // std::cout<<"myOutputSymMatrix:2"<<std::endl;
874 m=m1;
875 // std::cout<<"myOutputSymMatrix:3"<<std::endl;
876 for(int i=0;i<6;i++)
877 { for(int j=0;j<=i;j++)
878 { if(i<3) {
879 m[i][j]=m[i][j]/100; //mm*mm --> cm*cm
880 }
881 else if(j<3) {
882 m[i][j]=m[i][j]/10000; //mm*MeV --> cm*GeV
883 }
884 else {
885 m[i][j]=m[i][j]/1000000; //MeV*MeV --> GeV*GeV
886 }
887 }
888 }
889 // std::cout<<"myOutputSymMatrix:4"<<std::endl;
890 // std::cout<<"m1="<<m1<<std::endl;
891 // std::cout<<"m="<<m<<std::endl;
892 myOutputSM=m;
893 return myOutputSM;
894}
895
896void ExtSteppingAction::CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
897{
898 if((sector>=0)&&(sector<3))
899 sector+=16;
900 //if((sector!=7)&&(sector!=15))
901 {
902 if((nb>=0)&&(nb<4))
903 {
904 ntheta = 0;
905 nphi = (sector-3)*4+nb;
906 }
907 else if((nb>=4)&&(nb<8))
908 {
909 ntheta = 1;
910 nphi = (sector-3)*4+nb-4;
911 }
912 else if((nb>=8)&&(nb<13))
913 {
914 ntheta = 2;
915 nphi = (sector-3)*5+nb-8;
916 }
917 else if((nb>=13)&&(nb<18))
918 {
919 ntheta = 3;
920 nphi = (sector-3)*5+nb-13;
921 }
922 else if((nb>=18)&&(nb<24))
923 {
924 ntheta = 4;
925 nphi = (sector-3)*6+nb-18;
926 }
927 else if((nb>=24)&&(nb<30))
928 {
929 ntheta = 5;
930 nphi = (sector-3)*6+nb-24;
931 }
932 }
933
934 if(npart==2)
935 {
936 switch(ntheta)
937 {
938 case 0:
939 if(nphi<32)
940 nphi = 31-nphi;
941 else
942 nphi = 95-nphi;
943 break;
944 case 1:
945 if(nphi<32)
946 nphi = 31-nphi;
947 else
948 nphi = 95-nphi;
949 break;
950 case 2:
951 if(nphi<40)
952 nphi = 39-nphi;
953 else
954 nphi = 119-nphi;
955 break;
956 case 3:
957 if(nphi<40)
958 nphi = 39-nphi;
959 else
960 nphi = 119-nphi;
961 break;
962 case 4:
963 if(nphi<48)
964 nphi = 47-nphi;
965 else
966 nphi = 143-nphi;
967 break;
968 case 5:
969 if(nphi<48)
970 nphi = 47-nphi;
971 else
972 nphi = 143-nphi;
973 default:
974 ;
975 }
976 }
977}
978
979int ExtSteppingAction::CalculateEmcEndPhiNb(int num)
980{
981 if(num==0||num==1) {
982 return 15-num*8;
983 } else if(num==2||num==3) {
984 return 30-num*8;
985 } else if(num<=9) {
986 return 17-num;
987 } else {
988 return 15-num;
989 }
990}
991
992int ExtSteppingAction::CalculateEmcEndCopyNb(int num)
993{
994 int copyNb;
995 switch(num){
996 case 30:
997 copyNb = 5;
998 break;
999 case 31:
1000 copyNb = 6;
1001 break;
1002 case 32:
1003 copyNb = 14;
1004 break;
1005 case 33:
1006 copyNb = 15;
1007 break;
1008 case 34:
1009 copyNb = 16;
1010 break;
1011 default:
1012 copyNb = num;
1013 break;
1014 }
1015 return copyNb;
1016}
1017
1018void ExtSteppingAction::InfmodMuc(Hep3Vector &pos,Hep3Vector &mom,HepSymMatrix &err)
1019{
1020 pos = m_pos_mod;
1021 mom = m_mom_mod;
1022 err = m_err_mod;
1023}
1024
1025
1026Hep3Vector ExtSteppingAction::GetGapID(G4String vol)
1027{
1028 int par(-1),se(-1),ga(-1);
1029 G4String strPart = vol.substr(5,1);
1030 //G4String strSeg = m_MotherVolumeName.substr(7,1);
1031
1032 G4String strSeg = vol.substr(7,1);
1033 G4String strGap;
1034 if(vol.contains("G")) strGap= vol.substr(9,1);
1035 if(vol.contains("Ab")) strGap= vol.substr(10,1);
1036 std::istrstream partBuf(strPart.c_str(), strlen(strPart.c_str()));
1037 std::istrstream segBuf(strSeg.c_str(), strlen(strSeg.c_str()));
1038 std::istrstream gapBuf(strGap.c_str(), strlen(strGap.c_str()));
1039
1040 partBuf >> par ;
1041 segBuf >> se;
1042 gapBuf >> ga;
1043 if(vol.contains("Ab")&&par==1) ga = ga+1;
1044 //panelBuf >> m_panel;
1045 Hep3Vector vec;
1046 vec[0]=par;
1047 vec[1]=se;
1048 vec[2]=ga;
1049 return vec;
1050}
1051
1052
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
M_PI
#define M_PI
Definition: TConstant.h:4
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::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::SetPosMomErr
void SetPosMomErr(Hep3Vector pos, Hep3Vector mom, HepSymMatrix err)
Definition: InstallArea/include/TrkExtAlg/TrkExtAlg/ExtMucKal.h:25
ExtMucKal::XPmod
void XPmod(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtMucKal.cxx:71
ExtMucKal::SetMucDigiCol
void SetMucDigiCol(MucDigiCol *amucdigi)
Definition: InstallArea/include/TrkExtAlg/TrkExtAlg/ExtMucKal.h:26
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:1026
ExtSteppingAction::CalculateEmcEndThetaPhi
void CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
Definition: ExtSteppingAction.cxx:896
ExtSteppingAction::InfmodMuc
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtSteppingAction.cxx:1018
ExtSteppingAction::UserSteppingAction
void UserSteppingAction(const G4Step *currentStep)
Definition: ExtSteppingAction.cxx:103
ExtSteppingAction::CalculateEmcEndCopyNb
int CalculateEmcEndCopyNb(int num)
Definition: ExtSteppingAction.cxx:992
ExtSteppingAction::CalculateEmcEndPhiNb
int CalculateEmcEndPhiNb(int num)
Definition: ExtSteppingAction.cxx:979
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: InstallArea/include/TrkExtAlg/TrkExtAlg/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: InstallArea/include/TrkExtAlg/TrkExtAlg/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: InstallArea/include/TrkExtAlg/TrkExtAlg/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
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8