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