104{
105
106
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
119 Hep3Vector currentPosition = currentTrack->GetPosition();
120 Hep3Vector currentMomentum = currentTrack->GetMomentum();
121
122
123
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
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
146 if(currentTrackStatus == fAlive)
147 {
148
149 if(inner||mucdec)
150 {
151
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
170 G4Material* oldMaterial = oldVolume->GetLogicalVolume()->GetMaterial();
171 if(!oldMaterial) std::cout<<"Can't get oldMaterial"<<std::endl;
172 else CalculateChicc(oldMaterial);
173
174
175 if(!(extXpErr->
move(currentPosition,currentMomentum,currentB,1,chicc)))
176 if(msgFlag) cout<<"can not update Error Matrix!!"<<endl;
177
178
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;
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
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
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
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
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
301
302
303
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
333
334
335
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 }
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
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
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
394
395
396
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
428
429
430
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 }
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
470 if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
471 {
472 newVolumeNumber = -2;
473 if(myExtTrack)
474 {
475
476 Hep3Vector
nPhi(-
y/r,x/r,0.);
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
497 myEmcPathFlag = true;
498 }
499
500
501 if( (!myEmcFlag) && newVolumeName.contains("Crystal") )
502 {
503 if(myExtTrack)
504 {
505
506 int npart,ntheta,nphi;
507 if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) {
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 {
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 ;
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
530 Hep3Vector
nPhi(-
y/r,x/r,0.);
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
542 newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->
get_err()));
543 }
544 myEmcFlag = true;
545 return;
546 }
547
548
549 if(myEmcPathFlag && oldVolumeName.contains("Crystal") )
550 {
551 myPathOutCrystal = currentTrack->GetTrackLength();
552 myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
553
554 myEmcPathFlag=false;
555 }
556
557
558
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
598 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
599 m_trackstop=true;
600 return;
601 }
602 }
603
604
605
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
616 if(Flag1&&Flag2&&Flag3&&Flag5)
617 {
618
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 {
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 {
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
656 double dist = fabs(pos_loc.z());
657
658 if(!WillFilter&&RemDist>dist)
659 {
660 RemPositon = currentPosition;
661 RemMomentum = currentMomentum;
662 RemXpErr = XpErr;
663 RemDist = dist;
664 RemVol = oldVolumeName;
665 }
666
667
668 if(WillFilter)
669 {
676 vector<MucRecHit*> tmp = muckal->
TrackHit();
678
679 double chi2 = muckal->
GetChi2();
681 if(filter&&iniOK)
682 {
683
684 myMucnfit_++;
685
686 myMucchisq_ = myMucchisq_+chi2;
687 muckal->
XPmod(m_pos_mod,m_mom_mod,m_err_mod);
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
704 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
705
706 }
707 else
708 {
709
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
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
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
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
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.);
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}
void SetTof1Data(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, double aTof, double aPath, HepSymMatrix aErrorMatrix, double aZSigma=0., double aTSigma=0., double aXSigma=0., double aYSigma=0.)
void SetEmcData(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, double aThetaSigma, double aPhiSigma, HepSymMatrix aErrorMatrix)
void SetTof2Data(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, double aTof, double aPath, HepSymMatrix aErrorMatrix, double aZSigma=0., double aTSigma=0., double aXSigma=0., double aYSigma=0.)
void SetMucData(Hep3Vector aPosition, Hep3Vector aMomentum, string aVolumeName, int aVolumeNumber, HepSymMatrix aErrorMatrix, double aZSigma=0., double aTSigma=0., double aXSigma=0., double aYSigma=0.)
void SetEmcPath(double path)
void SetMucKalData(double chi2, int dof, double depth, int brLastLay, int ecLastLay, int nhits)
void SetGapID(Hep3Vector id)
void SetPosMomErr(Hep3Vector pos, Hep3Vector mom, HepSymMatrix err)
void SetMucWindow(int aMucWindow)
void XPmod(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
void SetMucDigiCol(MucDigiCol *amucdigi)
vector< MucRecHit * > TrackHit()
Hep3Vector GetGapID(G4String vol)
void CalculateEmcEndThetaPhi(int npart, int sector, int nb, int &ntheta, int &nphi)
int CalculateEmcEndCopyNb(int num)
int CalculateEmcEndPhiNb(int num)
void put_err(const double error[])
const HepSymMatrix & get_err() const
double get_plane_err(const Hep3Vector &np, const Hep3Vector &nr) const
void set_mom(const Hep3Vector &mom)
bool move(const Hep3Vector &xv1, const Hep3Vector &pv1, const Hep3Vector &B, const int ms_on, const double chi_cc)
void set_pos(const Hep3Vector &pos)
HepPoint3D TransformToGap(const HepPoint3D gPoint) const
Transform a point from global coordinate to gap coordinate.
MucGeoGap * GetGap(const int part, const int seg, const int gap) const
Get a pointer to the gap identified by (part,seg,gap).
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
void setPathInTof1(vector< double > x)
void setPathInTof2(vector< double > x)
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)