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<<"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;
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
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
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
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
269
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
303
304
305
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
335
336
337
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 }
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
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
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
396
397
398
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
430
431
432
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 }
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
472 if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
473 {
474 newVolumeNumber = -2;
475 if(myExtTrack)
476 {
477
478 Hep3Vector
nPhi(-
y/r,x/r,0.);
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
499 myEmcPathFlag = true;
500 }
501
502
503 if( (!myEmcFlag) && newVolumeName.contains("Crystal") )
504 {
505 if(myExtTrack)
506 {
507
508 int npart,ntheta,nphi;
509 if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) {
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 {
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 ;
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
532 Hep3Vector
nPhi(-
y/r,x/r,0.);
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
550 if(myEmcPathFlag && oldVolumeName.contains("Crystal") )
551 {
552 myPathOutCrystal = currentTrack->GetTrackLength();
553 myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
554
555 myEmcPathFlag=false;
556 }
557
558
559
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
599 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
600 m_trackstop=true;
601 return;
602 }
603 }
604
605
606
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
617 if(Flag1&&Flag2&&Flag3&&Flag5)
618 {
619
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 {
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 {
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
657 double dist = fabs(pos_loc.z());
658
659 if(!WillFilter&&RemDist>dist)
660 {
661 RemPositon = currentPosition;
662 RemMomentum = currentMomentum;
663 RemXpErr = XpErr;
664 RemDist = dist;
665 RemVol = oldVolumeName;
666 }
667
668
669 if(WillFilter)
670 {
677 vector<MucRecHit*> tmp = muckal->
TrackHit();
679
680 double chi2 = muckal->
GetChi2();
682 if(filter&&iniOK)
683 {
684
685 myMucnfit_++;
686
687 myMucchisq_ = myMucchisq_+chi2;
688 muckal->
XPmod(m_pos_mod,m_mom_mod,m_err_mod);
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
705 currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
706
707 }
708 else
709 {
710
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
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
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 }
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
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.);
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}
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)