ExtSteppingAction Class Reference

#include <ExtSteppingAction.h>

List of all members.

Public Member Functions

 ExtSteppingAction ()
 ~ExtSteppingAction ()
void Reset ()
void MucReset ()
void UserSteppingAction (const G4Step *currentStep)
void SetInitialPath (double aPath)
void SetInitialTof (double aTof)
void SetBetaInMDC (double aBeta)
void SetXpErrPointer (Ext_xp_err *xpErr)
void SetMsgFlag (bool aMsgFalg)
void SetMucKalFlag (bool aMucKalFlag)
void SetMucWindow (int aMucWindow)
void SetExtTrackPointer (RecExtTrack *aExtTrack)
void CalculateEmcEndThetaPhi (int npart, int sector, int nb, int &ntheta, int &nphi)
int CalculateEmcEndPhiNb (int num)
int CalculateEmcEndCopyNb (int num)
void Set_which_tof_version (int version)
int Get_which_tof_version (void)
void InfmodMuc (Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Hep3Vector GetGapID (G4String vol)
bool TrackStop ()
void SetMucDigiColPointer (MucDigiCol *rawdigicol)

Private Member Functions

void CalculateChicc (G4Material *currentMaterial)
HepSymMatrix & myOutputSymMatrix (const HepSymMatrix &)

Private Attributes

double chicc
double initialPath
double initialTof
double myBetaInMDC
double myPathIntoCrystal
double myPathOutCrystal
double myPathInCrystal
double myPathIntoTof1
double myPathOutTof1
vector< double > myPathInTof1
double myPathIntoTof2
double myPathOutTof2
vector< double > myPathInTof2
int myMucWindow
Ext_xp_errextXpErr
HepSymMatrix myOutputSM
MucDigiColm_mucdigicol
bool myUseMucKalFlag
Hep3Vector RememberID
bool m_trackstop
int myMucnfit_
double myMucchisq_
double myMucdepth_
int myMucbrLastLay_
int myMucecLastLay_
int myMucnhits_
HepSymMatrix m_err_mod
Hep3Vector m_pos_mod
Hep3Vector m_mom_mod
Hep3Vector RemPositon
Hep3Vector RemMomentum
HepSymMatrix RemXpErr
int RemStep
double RemDist
double RemDepth
Hep3Vector RemID
G4String RemVol
RecExtTrackmyExtTrack
bool msgFlag
bool myTofFlag
bool myTof1Flag
bool myTof2Flag
bool myInTof1
bool myOutTof1
bool myInTof2
bool myOutTof2
bool myPhyEmcFlag
bool myEmcFlag
bool myEmcPathFlag
bool myMucFlag
double myTof1R
double myTof1Z
double myTof2R
double myEmcR1
double myEmcR2
double myEmcZ
double myMucR
double myMucZ
int m_which_tof_version


Detailed Description

Definition at line 25 of file ExtSteppingAction.h.


Constructor & Destructor Documentation

ExtSteppingAction::ExtSteppingAction (  ) 

Definition at line 33 of file ExtSteppingAction.cxx.

References m_which_tof_version, myEmcR1, myEmcR2, myEmcZ, myMucR, myMucZ, myTof1R, myTof1Z, and myTof2R.

00033                                     :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)
00034 {
00035     myTof1R = 814.001;
00036     myTof1Z = 1330.0;
00037     myTof2R = 871.036;
00038     myEmcR1 = 945.0; 
00039     myEmcR2 = 983.9;
00040     myEmcZ  = 1416.8;
00041     myMucR  = 1700.0-0.01;
00042     myMucZ  = 2050.0-0.01;
00043     m_which_tof_version=2;
00044 }

ExtSteppingAction::~ExtSteppingAction (  ) 

Definition at line 46 of file ExtSteppingAction.cxx.

00046 {}


Member Function Documentation

void ExtSteppingAction::CalculateChicc ( G4Material *  currentMaterial  )  [private]

Definition at line 850 of file ExtSteppingAction.cxx.

References chicc, genRecEmupikp::i, and subSeperate::temp.

Referenced by UserSteppingAction().

00851 {       
00852     int n = currentMaterial->GetNumberOfElements();
00853     const double *p = currentMaterial->GetFractionVector();
00854     double density = (currentMaterial->GetDensity())/(g/cm3);
00855     double temp=0.0;
00856     for(int i=0;i<n;i++)
00857     {
00858         const G4Element* mElement = currentMaterial->GetElement(i);
00859         double z = mElement->GetZ();
00860         double a = mElement->GetN();
00861         //      std::cout<<"Fraction: "<<*p<<" z: "<<z<<" a: "<<a<<std::endl;
00862         temp += *(p++)/a*z*(z+1);
00863     }
00864     chicc = 0.39612e-3*std::sqrt(density*temp);         
00865     //   std::cout<<"chicc:"<<chicc<<std::endl;
00866 }

int ExtSteppingAction::CalculateEmcEndCopyNb ( int  num  ) 

Definition at line 992 of file ExtSteppingAction.cxx.

Referenced by UserSteppingAction().

00993 {
00994     int copyNb;
00995     switch(num){
00996         case 30:
00997             copyNb = 5;
00998             break;
00999         case 31:
01000             copyNb = 6;
01001             break;
01002         case 32:
01003             copyNb = 14;
01004             break;
01005         case 33:
01006             copyNb = 15;
01007             break;
01008         case 34:
01009             copyNb = 16;
01010             break;
01011         default:
01012             copyNb = num;
01013             break;
01014     }
01015     return copyNb;
01016 }

int ExtSteppingAction::CalculateEmcEndPhiNb ( int  num  ) 

Definition at line 979 of file ExtSteppingAction.cxx.

Referenced by UserSteppingAction().

00980 {
00981     if(num==0||num==1) {
00982         return 15-num*8;
00983     } else if(num==2||num==3) {
00984         return 30-num*8;
00985     } else if(num<=9) {
00986         return 17-num;
00987     } else {
00988         return 15-num;
00989     }
00990 }

void ExtSteppingAction::CalculateEmcEndThetaPhi ( int  npart,
int  sector,
int  nb,
int &  ntheta,
int &  nphi 
)

Definition at line 896 of file ExtSteppingAction.cxx.

Referenced by UserSteppingAction().

00897 {
00898     if((sector>=0)&&(sector<3))
00899       sector+=16;
00900     //if((sector!=7)&&(sector!=15))
00901     {
00902         if((nb>=0)&&(nb<4))
00903         {
00904             ntheta = 0;
00905             nphi = (sector-3)*4+nb;
00906         }
00907         else if((nb>=4)&&(nb<8))
00908         {
00909             ntheta = 1;
00910             nphi = (sector-3)*4+nb-4;
00911         }
00912         else if((nb>=8)&&(nb<13))
00913         {
00914             ntheta = 2;
00915             nphi = (sector-3)*5+nb-8;
00916         }
00917         else if((nb>=13)&&(nb<18))
00918         {
00919             ntheta = 3;
00920             nphi = (sector-3)*5+nb-13;
00921         }
00922         else if((nb>=18)&&(nb<24))
00923         {
00924             ntheta = 4;
00925             nphi = (sector-3)*6+nb-18;
00926         }
00927         else if((nb>=24)&&(nb<30))
00928         {
00929             ntheta = 5;
00930             nphi = (sector-3)*6+nb-24;
00931         }
00932     }
00933 
00934     if(npart==2)
00935     {
00936         switch(ntheta)
00937         {
00938             case 0:
00939                 if(nphi<32)
00940                   nphi = 31-nphi;
00941                 else
00942                   nphi = 95-nphi;
00943                 break;
00944             case 1:
00945                 if(nphi<32)
00946                   nphi = 31-nphi;
00947                 else
00948                   nphi = 95-nphi;
00949                 break;
00950             case 2:
00951                 if(nphi<40)
00952                   nphi = 39-nphi;
00953                 else
00954                   nphi = 119-nphi;
00955                 break;
00956             case 3:
00957                 if(nphi<40)
00958                   nphi = 39-nphi;
00959                 else
00960                   nphi = 119-nphi;
00961                 break;
00962             case 4:
00963                 if(nphi<48)
00964                   nphi = 47-nphi;
00965                 else
00966                   nphi = 143-nphi;
00967                 break;
00968             case 5:
00969                 if(nphi<48)
00970                   nphi = 47-nphi;
00971                 else
00972                   nphi = 143-nphi;
00973             default:
00974                 ;
00975         }
00976     }
00977 }

int ExtSteppingAction::Get_which_tof_version ( void   )  [inline]

Definition at line 51 of file ExtSteppingAction.h.

References m_which_tof_version.

00051 {return m_which_tof_version;}

Hep3Vector ExtSteppingAction::GetGapID ( G4String  vol  ) 

Definition at line 1026 of file ExtSteppingAction.cxx.

References vec.

Referenced by UserSteppingAction().

01027 {
01028     int par(-1),se(-1),ga(-1);
01029     G4String strPart       = vol.substr(5,1);
01030     //G4String strSeg        = m_MotherVolumeName.substr(7,1);
01031 
01032     G4String strSeg        = vol.substr(7,1);
01033     G4String strGap;
01034     if(vol.contains("G")) strGap= vol.substr(9,1);
01035     if(vol.contains("Ab")) strGap= vol.substr(10,1);
01036     std::istrstream partBuf(strPart.c_str(),             strlen(strPart.c_str()));
01037     std::istrstream segBuf(strSeg.c_str(),               strlen(strSeg.c_str()));
01038     std::istrstream gapBuf(strGap.c_str(),               strlen(strGap.c_str()));
01039 
01040     partBuf       >> par ;
01041     segBuf        >> se;
01042     gapBuf        >> ga;
01043     if(vol.contains("Ab")&&par==1) ga = ga+1;
01044     //panelBuf      >> m_panel;
01045     Hep3Vector vec;
01046     vec[0]=par;
01047     vec[1]=se;
01048     vec[2]=ga;
01049     return vec;
01050 }

void ExtSteppingAction::InfmodMuc ( Hep3Vector &  pos,
Hep3Vector &  mom,
HepSymMatrix &  err 
)

Definition at line 1018 of file ExtSteppingAction.cxx.

References m_err_mod, m_mom_mod, and m_pos_mod.

Referenced by TrkExtAlg::execute().

01019 {
01020     pos = m_pos_mod;
01021     mom = m_mom_mod;
01022     err = m_err_mod;
01023 }

void ExtSteppingAction::MucReset (  ) 

Definition at line 90 of file ExtSteppingAction.cxx.

References RemDepth, RemDist, RemID, RemStep, and RemVol.

Referenced by Ext_track::Set().

00091 {
00092     RemStep =0;
00093     RemDist = 99999.;
00094     RemDepth = 0.;
00095     RemID[0]=-1;
00096     RemID[1]=-1;
00097     RemID[2]=-1;
00098     RemVol = "";
00099 }

HepSymMatrix & ExtSteppingAction::myOutputSymMatrix ( const HepSymMatrix &   )  [private]

Definition at line 869 of file ExtSteppingAction.cxx.

References genRecEmupikp::i, ganga-rec::j, and myOutputSM.

Referenced by UserSteppingAction().

00870 {
00871     //  std::cout<<"myOutputSymMatrix:1"<<std::endl;
00872     HepSymMatrix m;
00873     //  std::cout<<"myOutputSymMatrix:2"<<std::endl;
00874     m=m1;
00875     //  std::cout<<"myOutputSymMatrix:3"<<std::endl;
00876     for(int i=0;i<6;i++)
00877     {  for(int j=0;j<=i;j++)
00878         {  if(i<3) {
00879                        m[i][j]=m[i][j]/100;     //mm*mm --> cm*cm
00880                    }
00881         else if(j<3) {
00882             m[i][j]=m[i][j]/10000;   //mm*MeV --> cm*GeV
00883         }
00884         else {
00885             m[i][j]=m[i][j]/1000000; //MeV*MeV --> GeV*GeV
00886         }
00887         }       
00888     }
00889     //  std::cout<<"myOutputSymMatrix:4"<<std::endl;
00890     //  std::cout<<"m1="<<m1<<std::endl;
00891     //  std::cout<<"m="<<m<<std::endl;
00892     myOutputSM=m;
00893     return myOutputSM;
00894 }

void ExtSteppingAction::Reset (  ) 

Definition at line 48 of file ExtSteppingAction.cxx.

References chicc, m_trackstop, myEmcFlag, myEmcPathFlag, myExtTrack, myInTof1, myInTof2, myMucbrLastLay_, myMucchisq_, myMucdepth_, myMucecLastLay_, myMucFlag, myMucnfit_, myOutTof1, myOutTof2, myPathInCrystal, myPathIntoCrystal, myPathInTof1, myPathInTof2, myPathIntoTof1, myPathIntoTof2, myPathOutCrystal, myPathOutTof1, myPathOutTof2, myPhyEmcFlag, myTof1Flag, myTof2Flag, myTofFlag, and RememberID.

Referenced by TrkExtAlg::execute().

00049 {
00050     chicc        = 0.;
00051     myExtTrack   = 0;
00052     myPathIntoCrystal = 0.;
00053     myPathOutCrystal  = 0.;
00054     myPathInCrystal   = 0.;
00055 
00056     myPathIntoTof1 = 0.0;
00057     myPathOutTof1  = 0.0;
00058     //myPathInTof1   = 0.0;
00059     myPathInTof1.clear();
00060 
00061     myPathIntoTof2 = 0.0;
00062     myPathOutTof2  = 0.0;
00063     //myPathInTof2   = 0.0;
00064     myPathInTof2.clear();
00065 
00066     myTofFlag    = false;
00067     myTof1Flag   = false;
00068     myTof2Flag   = false;
00069     myInTof1     = false;
00070     myInTof2     = false;
00071     myOutTof1    = false;
00072     myOutTof2    = false;
00073     myPhyEmcFlag = false;
00074     myEmcFlag    = false;
00075     myEmcPathFlag= false;
00076     myMucFlag    = false;
00077 
00078     m_trackstop =false;
00079     myMucchisq_ =0.;
00080     myMucnfit_ =0;
00081     myMucdepth_=-99.;
00082     myMucbrLastLay_=-1;
00083     myMucecLastLay_=-1;
00084     RememberID[0]=-1;
00085     RememberID[1]=-1;
00086     RememberID[2]=-1;
00087 
00088 }

void ExtSteppingAction::Set_which_tof_version ( int  version  )  [inline]

Definition at line 50 of file ExtSteppingAction.h.

References m_which_tof_version.

Referenced by TrkExtAlg::execute().

00050 {m_which_tof_version=version;}

void ExtSteppingAction::SetBetaInMDC ( double  aBeta  )  [inline]

Definition at line 38 of file ExtSteppingAction.h.

References myBetaInMDC.

Referenced by Ext_track::Set().

00038 {myBetaInMDC = aBeta;};

void ExtSteppingAction::SetExtTrackPointer ( RecExtTrack aExtTrack  )  [inline]

Definition at line 44 of file ExtSteppingAction.h.

References myExtTrack.

Referenced by TrkExtAlg::execute().

00044 {myExtTrack = aExtTrack;}; 

void ExtSteppingAction::SetInitialPath ( double  aPath  )  [inline]

Definition at line 36 of file ExtSteppingAction.h.

References initialPath.

Referenced by Ext_track::Set().

00036 {initialPath = aPath;};

void ExtSteppingAction::SetInitialTof ( double  aTof  )  [inline]

Definition at line 37 of file ExtSteppingAction.h.

References initialTof.

Referenced by Ext_track::Set().

00037 {initialTof = aTof;};

void ExtSteppingAction::SetMsgFlag ( bool  aMsgFalg  )  [inline]

Definition at line 41 of file ExtSteppingAction.h.

References msgFlag.

Referenced by Ext_track::Initialization().

00041 {msgFlag = aMsgFalg;};

void ExtSteppingAction::SetMucDigiColPointer ( MucDigiCol rawdigicol  )  [inline]

Definition at line 56 of file ExtSteppingAction.h.

References m_mucdigicol.

Referenced by TrkExtAlg::execute().

00056 { m_mucdigicol = rawdigicol;}

void ExtSteppingAction::SetMucKalFlag ( bool  aMucKalFlag  )  [inline]

Definition at line 42 of file ExtSteppingAction.h.

References myUseMucKalFlag.

Referenced by Ext_track::Initialization().

00042 {myUseMucKalFlag=aMucKalFlag;};

void ExtSteppingAction::SetMucWindow ( int  aMucWindow  )  [inline]

Definition at line 43 of file ExtSteppingAction.h.

References myMucWindow.

Referenced by Ext_track::Initialization().

00043 {myMucWindow=aMucWindow;};

void ExtSteppingAction::SetXpErrPointer ( Ext_xp_err xpErr  )  [inline]

Definition at line 39 of file ExtSteppingAction.h.

References extXpErr.

Referenced by Ext_track::Set().

00039 {extXpErr = xpErr;};

bool ExtSteppingAction::TrackStop (  )  [inline]

Definition at line 55 of file ExtSteppingAction.h.

References m_trackstop.

Referenced by TrkExtAlg::execute().

00055 {return m_trackstop;}

void ExtSteppingAction::UserSteppingAction ( const G4Step *  currentStep  ) 

Definition at line 103 of file ExtSteppingAction.cxx.

References CalculateChicc(), CalculateEmcEndCopyNb(), CalculateEmcEndPhiNb(), CalculateEmcEndThetaPhi(), chicc, check_raw_filter::dist, ExtMucKal::ExtMucFilter(), extXpErr, Ext_errmx::get_err(), Ext_errmx::get_plane_err(), ExtMucKal::GetChi2(), MucGeoGeneral::GetGap(), GetGapID(), ExtMucKal::GetHitGap(), ExtMucKal::GetSameStrip(), genRecEmupikp::i, initialPath, initialTof, MucGeoGeneral::Instance(), m_err_mod, m_mom_mod, m_mucdigicol, M_PI, m_pos_mod, m_trackstop, Ext_xp_err::move(), msgFlag, ExtMucKal::MucKalIniti(), myEmcFlag, myEmcPathFlag, myExtTrack, myInTof1, myInTof2, myMucbrLastLay_, myMucchisq_, myMucdepth_, myMucecLastLay_, myMucFlag, myMucnfit_, myMucnhits_, myMucR, myMucWindow, myMucZ, myOutputSymMatrix(), myOutTof1, myOutTof2, myPathInCrystal, myPathIntoCrystal, myPathInTof1, myPathInTof2, myPathIntoTof1, myPathIntoTof2, myPathOutCrystal, myPathOutTof1, myPathOutTof2, myPhyEmcFlag, myTof1Flag, myTof2Flag, myUseMucKalFlag, nPhi, Ext_errmx::put_err(), rb::R(), RemDepth, RemDist, RememberID, RemID, RemMomentum, RemPositon, RemStep, RemVol, RemXpErr, Ext_xp_err::set_mom(), Ext_xp_err::set_pos(), DstExtTrack::SetEmcData(), DstExtTrack::SetEmcPath(), ExtMucKal::SetGapID(), DstExtTrack::SetMucData(), ExtMucKal::SetMucDigiCol(), DstExtTrack::SetMucKalData(), ExtMucKal::SetMucWindow(), RecExtTrack::setPathInTof1(), RecExtTrack::setPathInTof2(), ExtMucKal::SetPosMomErr(), DstExtTrack::SetTof1Data(), DstExtTrack::SetTof2Data(), ExtMucKal::TrackHit(), MucGeoGap::TransformToGap(), x, and ExtMucKal::XPmod().

00104 {
00105 
00106     //Get track and its status
00107     G4Track* currentTrack = currentStep->GetTrack();
00108     if(!currentTrack) 
00109     {
00110         cout<<"Can't get currentTrack"<<endl;
00111         return;
00112     }
00113     G4TrackStatus currentTrackStatus = currentTrack->GetTrackStatus();
00114 
00115     int stepNumber = currentTrack->GetCurrentStepNumber();
00116     if(msgFlag)  cout<<"STEP "<<stepNumber<<":"<<endl;
00117 
00118     //Get current position and momentum
00119     Hep3Vector currentPosition = currentTrack->GetPosition();
00120     Hep3Vector currentMomentum = currentTrack->GetMomentum();
00121 
00122 
00123     //Get current Volume
00124     G4VPhysicalVolume* oldVolume;
00125     G4VPhysicalVolume* newVolume;
00126     if(!currentTrack->GetTouchableHandle()) cout<<"Can't get currentTrack->GetTouchableHandle()"<<endl;
00127     else  oldVolume= currentTrack->GetTouchableHandle()->GetVolume();
00128     if(!currentTrack->GetNextTouchableHandle()) cout<<"Can't get currentTrack->GetNextTouchableHandle()"<<endl;
00129     else  newVolume= currentTrack->GetNextTouchableHandle()->GetVolume();
00130     if(!oldVolume) cout<<"Can't get oldVolume!"<<endl;
00131 
00132     //****added by LI Chunhua
00133     if(stepNumber>50000) {
00134         cout<<"infinite steps in the track "<<endl;
00135         currentTrack->SetTrackStatus(fStopAndKill); m_trackstop =true;
00136     }
00137 
00138     G4String ParticleName = currentTrack->GetDefinition()->GetParticleName();
00139     double x_ = currentPosition.x();
00140     double y_ = currentPosition.y();
00141     double z_ = currentPosition.z();
00142     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)) );
00143     bool mucdec = (fabs(x_)>=myMucR) || (fabs(y_)>=myMucR) || ((fabs(x_-y_)/sqrt(2.))>=myMucR) || ((fabs(x_+y_)/sqrt(2.))>=myMucR) || (fabs(z_)>=myMucZ);
00144 
00145     //if current particle is alive.
00146     if(currentTrackStatus == fAlive)
00147     {
00148         //if(oldVolume!=newVolume)//if so, update Error Matrix and chicc.
00149         if(inner||mucdec)//update Error Matrix for newvolume in MDC,TOF,EMC; update Error Matrix step by step in MUC;
00150         {
00151             //Get current B field
00152             double currentPoint[3]={currentPosition.x(),currentPosition.y(),currentPosition.z()};
00153             double currentBfield[3]={0.0,0.0,0.0};
00154             Hep3Vector currentB(0.0,0.0,1.0);
00155             if(G4TransportationManager::GetTransportationManager())
00156             {
00157                 G4FieldManager* fieldManager=G4TransportationManager::GetTransportationManager()->GetFieldManager();
00158                 if(fieldManager)       
00159                 {         
00160                     if(fieldManager->GetDetectorField())        
00161                     {
00162                         fieldManager->GetDetectorField()->GetFieldValue(currentPoint,currentBfield);      
00163                         currentB.set(currentBfield[0]/tesla,currentBfield[1]/tesla,currentBfield[2]/tesla);
00164                         if(msgFlag) cout<<"B:"<<currentB.x()<<","<<currentB.y()<<","<<currentB.z()<<endl;
00165                     }
00166                 }
00167             }
00168 
00169             //Get chicc
00170             G4Material* oldMaterial = oldVolume->GetLogicalVolume()->GetMaterial();
00171             if(!oldMaterial) std::cout<<"Can't get oldMaterial"<<std::endl;
00172             else CalculateChicc(oldMaterial);
00173 
00174             //Update Error Matrix
00175             if(!(extXpErr->move(currentPosition,currentMomentum,currentB,1,chicc)))
00176               if(msgFlag) cout<<"can not update Error Matrix!!"<<endl;
00177 
00178             //Verbose
00179             if(msgFlag)
00180             {
00181                 cout<<"Volume name:"<<newVolume->GetName()<<endl;
00182                 cout<<"Volume number:"<<newVolume->GetCopyNo()<<endl;
00183                 cout<<"x:"<<currentPoint[0]<<"//y:"<<currentPoint[1]<<"//z:"<<currentPoint[2]
00184                     <<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"
00185                     <<currentMomentum.z()<<endl;
00186                 cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
00187                 cout<<"phi:"<<atan(currentPoint[1]/currentPoint[0])<<endl;
00188                 Hep3Vector nz(0.,0.,1.);
00189                 cout<<"Projected z error:"<<extXpErr->get_plane_err(currentMomentum.unit(),nz)
00190                     <<endl;
00191                 double x,y,R;
00192                 x = currentPoint[0];
00193                 y = currentPoint[1];
00194                 R = sqrt(x*x+y*y);
00195                 Hep3Vector nt(-y/R,x/R,0.);
00196                 cout<<"Projected phi error:"<<(extXpErr->get_plane_err(currentMomentum.unit(),nt))/R
00197                     <<endl<<endl;
00198             }
00199 
00200             //some often used things
00201             Hep3Vector xVector(1.0,0,0);
00202             Hep3Vector yVector(0,1.0,0);
00203             Hep3Vector zVector(0,0,1.0);
00204             Hep3Vector tzVector;
00205             tzVector.setRThetaPhi(1.0,M_PI/2.0,currentPosition.phi());
00206             double r = currentPosition.perp();
00207             double x = currentPosition.x();
00208             double y = currentPosition.y();
00209             double z = currentPosition.z();
00210             G4String newVolumeName = newVolume->GetName();
00211             G4String oldVolumeName = oldVolume->GetName();
00212             G4StepPoint* postStepPoint = currentStep->GetPostStepPoint();
00213             G4TouchableHistory* theTouchable = (G4TouchableHistory*)(postStepPoint->GetTouchable());
00214             int newVolumeNumber=newVolume->GetCopyNo();
00215             //          int newVolumeNumber=theTouchable->GetReplicaNumber(2);
00216 
00217             G4String name1;
00218             if(newVolumeName=="logical_gasLayer")
00219             {
00220                 name1=theTouchable->GetVolume(3)->GetName();
00221                 newVolumeNumber=theTouchable->GetReplicaNumber(3);
00222             }
00223 
00224 
00225 
00226 
00227 
00228 
00229             //Get PhyTof data
00230 
00231             //std::cout << "ExtSteppingAction        NewVolumeName:  "  <<  newVolumeName << std::endl;
00232             //std::cout << "ExtSteppingAction        OldVolumeName:  "  <<  oldVolumeName << std::endl;
00233             //std::cout << "ExtSteppingAction     name1:  " <<name1<<std::endl;
00234             //Comment this part to remove the airtight Tof. Most hits would be replaced by the following specific sensitive 
00235             //detector. The remianing few hits are abondoned in Tof reconstruction.
00236             /*
00237                if( (!myTofFlag) && (!myTof1Flag) && (newVolumeName.contains("Tof") ))
00238                {
00239                newVolumeNumber = -2;
00240                double currentTrackLength = currentTrack->GetTrackLength();
00241                double totalTrackLength = currentTrackLength + initialPath;
00242             //double initialTof = initialPath/(myBetaInMDC*299.792458);
00243             //cout<<"initialTof = "<<initialTof<<endl;
00244             //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
00245             double localTime = currentTrack->GetLocalTime();
00246             //cout<<"localTime = "<<localTime<<endl;
00247             double totalTOF = localTime + initialTof;
00248             //cout<<"totalTOF= "<<totalTOF<<endl;
00249 
00250             //std::cout << "ExtSteppingAction        Volumename contains one of the marker words" << std::endl;
00251 
00252             if(myExtTrack)
00253             {
00254             double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00255             double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00256             double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00257             double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00258             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);
00259             myTofFlag = true;
00260             }
00261             return;
00262             }//close if (!myTofFlag)
00263             */
00264 
00265 
00266 
00267             //Get Tof layer1 Ext data
00268             if((!myTof1Flag) && (newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") || 
00269                             (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))))
00270             {
00271                 double currentTrackLength = currentTrack->GetTrackLength();
00272                 double totalTrackLength = currentTrackLength+initialPath;
00273                 double localTime = currentTrack->GetLocalTime();
00274                 double totalTOF = localTime + initialTof;
00275                 myInTof1 = true;
00276                 myPathIntoTof1 = currentTrackLength;
00277                 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
00278 
00279                 newVolumeNumber=theTouchable->GetReplicaNumber(2);
00280 
00281                 if(newVolumeName.contains("ScinEc")) {
00282                     newVolumeNumber=(95-newVolumeNumber)/2;
00283                     newVolumeNumber=newVolumeNumber+176;
00284 
00285                     if(newVolumeName.contains("East")) newVolumeNumber=newVolumeNumber+48;
00286                     if(myExtTrack) {
00287                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00288                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00289                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00290                         double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00291                         myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00292                         myTof1Flag = true;
00293                     }
00294                 }
00295 
00296                 else if(newVolumeName=="logical_gasLayer")
00297                 {
00298                     G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
00299                     G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(currentPostPosition);
00300                     //cout<<"localPosition: "<<localPosition.x()<<"  "<<localPosition.y()<<"  "<<localPosition.z()<<endl;
00301 
00302                     //Here we assume the spread of one strip is 24+3=27mm. 
00303                     //Distance between strips and lower large glass margin: 6, top: 8, z_strip=z_small_glass-0.5
00304                     double z_mm = localPosition.z()-0.5+(24+3)*6;
00305                     int strip;
00306                     if(z_mm<=0)
00307                     {
00308                         strip=0;
00309                     }
00310                     else if(z_mm>0 && z_mm<12*27)
00311                     {
00312                         strip=floor(z_mm/27);
00313                     }
00314                     else
00315                     {
00316                         strip=11;
00317                     }
00318                     if(strip<0) strip=0;
00319                     if(strip>11) strip=11;
00320 
00321                     newVolumeNumber=theTouchable->GetReplicaNumber(3);
00322                     newVolumeNumber = 35-newVolumeNumber;
00323                     if( currentPosition.z() < 0.0  ) { newVolumeNumber = newVolumeNumber + 36; }
00324                     newVolumeNumber = 12*newVolumeNumber+strip;
00325                     newVolumeNumber = newVolumeNumber + 176 + 96;
00326                     if(myExtTrack) {
00327                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00328                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00329                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00330                         double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00331                         Hep3Vector locP = Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
00332                         //cout<<"locP: "<<locP.x()<<"  "<<locP.y()<<"  "<<locP.z()
00333                         //<<"  currentMomentum: "<<currentMomentum.x()<<"  "<<currentMomentum.y()<<"  "<<currentMomentum.z()
00334                         //<<"  newVolumeName: "<<newVolumeName<<"  newVolumeNumber: "<<newVolumeNumber
00335                         //<<"  totalTOF: "<<totalTOF<<"  totalTrackLength: "<<totalTrackLength<<endl;
00336                         myExtTrack->SetTof1Data(locP/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00337                         myTof1Flag = true;
00338                     }
00339                 }
00340                 else
00341                 { 
00342                     newVolumeNumber=(527-newVolumeNumber)/3;
00343                     if(myExtTrack) {
00344                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00345                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00346                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00347                         double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00348                         myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00349                         myTof1Flag = true;
00350                     }
00351                 }
00352 
00353                 return;
00354 
00355             }//close if (!myTof1Flag)
00356 
00357             if( myInTof1 && ( oldVolumeName=="logicalScinBr1" || oldVolumeName.contains("ScinEc") ||
00358                             (oldVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))) ) {
00359                 myInTof1 = false;
00360                 myOutTof1 = true;
00361                 myPathOutTof1 = currentTrack->GetTrackLength();
00362                 myPathInTof1.push_back(myPathOutTof1 - myPathIntoTof1);
00363                 if(msgFlag) cout << "myPathOutTof1 = " << myPathOutTof1 << endl;
00364             }
00365 
00366             if( myOutTof1 && ( newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") || 
00367                             (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m0" || name1=="logical_container_m1"))) ) {
00368                 myInTof1 = true;
00369                 myOutTof1 = false;
00370                 myPathIntoTof1 = currentTrack->GetTrackLength();
00371                 if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
00372             }
00373 
00374 
00375             //Get Tof layer2 Ext data
00376             if( (!myTof2Flag) && (newVolumeName=="logicalScinBr2" || 
00377                             (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3")))) 
00378             {
00379                 double currentTrackLength = currentTrack->GetTrackLength();
00380                 double totalTrackLength = currentTrackLength+initialPath;
00381                 //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
00382                 double localTime = currentTrack->GetLocalTime();
00383                 double totalTOF = localTime + initialTof;
00384                 myInTof2 = true;
00385                 myPathIntoTof2 = currentTrackLength;
00386                 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
00387                 newVolumeNumber=theTouchable->GetReplicaNumber(2);
00388 
00389                 if(newVolumeName=="logical_gasLayer")
00390                 {
00391                     G4ThreeVector currentPostPosition = currentStep->GetPostStepPoint()->GetPosition();
00392                     G4ThreeVector localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(currentPostPosition);
00393                     //cout<<"localPosition: "<<localPosition.x()<<"  "<<localPosition.y()<<"  "<<localPosition.z()<<endl;
00394 
00395                     //Here we assume the spread of one strip is 24+3=27mm. 
00396                     //Distance between strips and lower large glass margin: 6, top: 8, z_strip=z_small_glass-0.5
00397                     double z_mm = localPosition.z()-0.5+(24+3)*6;
00398                     int strip;
00399                     if(z_mm<=0)
00400                     {
00401                         strip=0;
00402                     }
00403                     else if(z_mm>0 && z_mm<12*27)
00404                     {
00405                         strip=floor(z_mm/27);
00406                     }
00407                     else
00408                     {
00409                         strip=11;
00410                     }
00411                     if(strip<0) strip=0;
00412                     if(strip>11) strip=11;
00413 
00414                     newVolumeNumber=theTouchable->GetReplicaNumber(3);
00415                     newVolumeNumber = 35-newVolumeNumber;
00416                     if( currentPosition.z() < 0.0  ) { newVolumeNumber = newVolumeNumber + 36; }
00417                     newVolumeNumber = 12*newVolumeNumber+strip;
00418                     newVolumeNumber = newVolumeNumber + 176 + 96;
00419 
00420                     if(myExtTrack) 
00421                     {
00422                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00423                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00424                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00425                         double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00426                         Hep3Vector locP = Hep3Vector( localPosition.x(), localPosition.y(), localPosition.z() );
00427                         //cout<<"locP: "<<locP.x()<<"  "<<locP.y()<<"  "<<locP.z()
00428                         //<<"  currentMomentum: "<<currentMomentum.x()<<"  "<<currentMomentum.y()<<"  "<<currentMomentum.z()
00429                         //<<"  newVolumeName: "<<newVolumeName<<"  newVolumeNumber: "<<newVolumeNumber
00430                         //<<"  totalTOF: "<<totalTOF<<"  totalTrackLength: "<<totalTrackLength<<endl;
00431                         myExtTrack->SetTof2Data(locP/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00432                         myTof2Flag = true;
00433                     }
00434                 }
00435                 else
00436                 {
00437                     newVolumeNumber=(527-newVolumeNumber)/3;
00438                     if(myExtTrack)
00439                     {
00440                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00441                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00442                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00443                         double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00444                         myExtTrack->SetTof2Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00445                         myTof2Flag = true;
00446                     }
00447                 }
00448                 return;
00449             }//close if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" )
00450 
00451             if( myInTof2 && (oldVolumeName=="logicalScinBr2" || 
00452                             (oldVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3"))) ) {
00453                 myInTof2 = false;
00454                 myOutTof2 = true;
00455                 myPathOutTof2 = currentTrack->GetTrackLength();
00456                 myPathInTof2.push_back(myPathOutTof2 - myPathIntoTof2);
00457                 if(msgFlag) cout << "myPathOutTof2 = " << myPathOutTof2 << endl;
00458             }
00459 
00460             if( myOutTof2 && (newVolumeName=="logicalScinBr2" || 
00461                             (newVolumeName=="logical_gasLayer" && (name1=="logical_container_m2" || name1=="logical_container_m3"))) ) {
00462                 myInTof2 = true;
00463                 myOutTof2 = false;
00464                 myPathIntoTof2 = currentTrack->GetTrackLength();
00465                 if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
00466             }
00467 
00468 
00469             //Get PhyEmc data.
00470             if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
00471             {
00472                 newVolumeNumber = -2;
00473                 if(myExtTrack)
00474                 {
00475                     //                          myPathIntoCrystal = currentTrack->GetTrackLength();
00476                     Hep3Vector nPhi(-y/r,x/r,0.);
00477                     double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
00478 
00479                     Hep3Vector aPosition = currentPosition;
00480                     double R = aPosition.r();
00481                     double aTheta = aPosition.theta();
00482                     aPosition.setTheta(aTheta + M_PI_2);
00483                     double errorTheta;
00484                     errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
00485                     if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
00486                     myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
00487                 }
00488                 myPhyEmcFlag = true;
00489                 return;
00490             }
00491 
00492 
00493             if(!myEmcPathFlag &&  newVolumeName.contains("Crystal")  )
00494             {
00495                 myPathIntoCrystal = currentTrack->GetTrackLength();
00496                 //                      cout<<"Enter Crystal, path: "<<myPathIntoCrystal<<endl;
00497                 myEmcPathFlag = true;
00498             }
00499 
00500             //Get Emc Ext data.
00501             if( (!myEmcFlag) && newVolumeName.contains("Crystal")  )
00502             {
00503                 if(myExtTrack)
00504                 {
00505                     //------------------- record crystal number
00506                     int npart,ntheta,nphi;
00507                     if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) { //barrel
00508                         npart=1;
00509                         std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName()).substr(16,2));
00510                         thetaBuf >> ntheta ;
00511                         nphi = 308-currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
00512                     } else {  //endcap
00513                         npart=2-2*currentTrack->GetNextTouchableHandle()->GetCopyNumber(3);
00514                         int endSector=currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
00515                         int endNb,endNbGDML;
00516                         std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(0)->GetName()).substr(20,2));
00517                         thetaBuf  >> endNb ;
00518                         endNbGDML=CalculateEmcEndCopyNb(endNb);
00519                         int endSectorGDML=CalculateEmcEndPhiNb(endSector);
00520                         CalculateEmcEndThetaPhi(npart,endSectorGDML,endNbGDML,ntheta,nphi);
00521                     }
00522                     ostringstream strEmc;
00523                     if(ntheta<10) {
00524                         strEmc<<"EmcPart"<<npart<<"Theta0"<<ntheta<<"Phi"<<nphi;
00525                     } else {
00526                         strEmc<<"EmcPart"<<npart<<"Theta"<<ntheta<<"Phi"<<nphi;
00527                     } //------------------------------------------
00528 
00529                     //                          myPathIntoCrystal = currentTrack->GetTrackLength();
00530                     Hep3Vector nPhi(-y/r,x/r,0.);
00531                     double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
00532 
00533                     Hep3Vector aPosition = currentPosition;
00534                     double R = aPosition.r();
00535                     double aTheta = aPosition.theta();
00536                     aPosition.setTheta(aTheta + M_PI_2);
00537                     double errorTheta;
00538                     errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
00539                     if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
00540                     myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,strEmc.str(),
00541                                 //myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,
00542                         newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
00543                                 }
00544                                 myEmcFlag = true;
00545                                 return;
00546                                 }
00547 
00548                                 //Get path in Emc
00549                                 if(myEmcPathFlag &&  oldVolumeName.contains("Crystal") )
00550                                 {
00551                                 myPathOutCrystal = currentTrack->GetTrackLength();
00552                                 myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
00553                                 //                      cout<<"Go out of crystal, path: "<<myPathOutCrystal<<endl;
00554                                 myEmcPathFlag=false;
00555                                 }
00556 
00557 
00558                                 //Get Muc Ext Data.
00559                                 if( (!myMucFlag) && ( (fabs(x)>=myMucR) || (fabs(y)>=myMucR) || ((fabs(x-y)/sqrt(2.))>=myMucR) || ((fabs(x+y)/sqrt(2.))>=myMucR) || (fabs(z)>=myMucZ) ) )
00560                                 {
00561                                     int newVolumeNumber = newVolume->GetCopyNo();
00562                                     if(myExtTrack)
00563                                     {
00564                                         Hep3Vector xVector(1.0,0,0); 
00565                                         Hep3Vector yVector(0,1.0,0);
00566                                         Hep3Vector zVector(0,0,1.0);
00567                                         double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00568                                         double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00569                                         double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00570                                         double tzError;
00571                                         double phi = currentPosition.phi();
00572                                         if(phi<0.) phi+=M_PI;
00573                                         int i = int(8.0*phi/M_PI);
00574                                         if( i==0 || i==7 || i==8 )
00575                                         {
00576                                             tzError = yError;
00577                                         }
00578                                         if( i==1 || i==2 )
00579                                         {
00580                                             Hep3Vector tzVector(-1.0,1.0,0.);
00581                                             tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00582                                         }
00583                                         if( i==3 || i==4 )
00584                                         {
00585                                             tzError = xError;
00586                                         } 
00587                                         if( i==5 || i==6 )
00588                                         {
00589                                             Hep3Vector tzVector(1.0,1.0,0.);
00590                                             tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00591                                         }
00592                                         myExtTrack->SetMucData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00593                                     }
00594                                     myMucFlag = true;
00595                                     if(!(ParticleName.contains("mu")&&myUseMucKalFlag))
00596                                     {
00597                                         //currentTrack->SetKineticEnergy(0.0);
00598                                         currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
00599                                         m_trackstop=true;
00600                                         return;
00601                                     }
00602                                 }
00603 
00604                     //**************************************************
00605                     //************inset KalFilter Algorithm in MUC******
00606                     //**************************************************
00607                     HepSymMatrix XpErr = extXpErr->get_err();
00608                     int namesize = oldVolumeName.size();
00609                     bool Flag1(false),Flag2(false),Flag3(false),Flag4(false),Flag5(false);
00610                     Flag1 = m_mucdigicol->size()>0;
00611                     Flag2 = myUseMucKalFlag;
00612                     Flag3 = ParticleName.contains("mu")&&oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&(oldVolumeName.contains("G")||oldVolumeName.contains("Ab"));
00613                     Flag4 = oldVolumeName.contains("lMuc")&&oldVolumeName.contains("P")&&oldVolumeName.contains("S")&&((namesize==10&&oldVolumeName.contains("G"))||(namesize==11&&oldVolumeName.contains("Ab")));
00614                     Flag5 = !((RemID[0]==1&&RemID[2]==9)||((RemID[0]==0||RemID[0]==2)&&RemID[2]==8));
00615                     // if(!Flag5) {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
00616                     if(Flag1&&Flag2&&Flag3&&Flag5)
00617                     {
00618                         //get depth in Ab
00619                         double depth = currentStep-> GetStepLength();
00620                         if(oldVolumeName.contains("Ab"))
00621                           RemDepth = RemDepth+ depth;
00622                         if(RemStep==0&&Flag4)
00623                         {
00624                             Hep3Vector ID_1 = GetGapID(oldVolumeName);
00625                             RemID = ID_1;
00626 
00627                             bool LastLay = (ID_1[0]==1&&ID_1[2]<9)||((ID_1[0]==0||ID_1[0]==2)&&ID_1[2]<8);
00628                             if(RememberID[2]!=ID_1[2]&&LastLay)
00629                             {
00630                                 Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(ID_1[0],ID_1[1],ID_1[2])->TransformToGap(currentPosition);
00631                                 double dist = fabs(pos_loc.z());
00632                                 RemPositon = currentPosition;
00633                                 RemMomentum = currentMomentum;
00634                                 RemXpErr = XpErr;
00635                                 RemDist = dist;
00636                                 RemStep++;
00637                             }
00638                         }
00639                         if(RemStep>0)
00640                         {
00641                             bool WillFilter = false;
00642                             bool LastLay_ = false;
00643                             double dist_b = 99999.;
00644                             Hep3Vector ID_2;
00645                             if(Flag4)
00646                             {
00647                                 ID_2 = GetGapID(oldVolumeName);
00648                                 LastLay_ =(ID_2[0]==1&&ID_2[2]<9)||((ID_2[0]==0||ID_2[0]==2)&&ID_2[2]<8);
00649                                 if(LastLay_)dist_b=fabs(MucGeoGeneral::Instance()->GetGap(ID_2[0],ID_2[1],ID_2[2])->TransformToGap(currentPosition).z());
00650                                 if(RemID!=ID_2)
00651                                   WillFilter=true;
00652 
00653                             }
00654 
00655                             Hep3Vector pos_loc = MucGeoGeneral::Instance()->GetGap(RemID[0],RemID[1],RemID[2])->TransformToGap(currentPosition);
00656                             double dist = fabs(pos_loc.z());
00657                             //get the nearest point from the center of gap
00658                             if(!WillFilter&&RemDist>dist)
00659                             {
00660                                 RemPositon = currentPosition;
00661                                 RemMomentum = currentMomentum;
00662                                 RemXpErr = XpErr;
00663                                 RemDist = dist;
00664                                 RemVol = oldVolumeName;
00665                             }
00666 
00667                             //if find the nearest point in the gap, Open Fillter 
00668                             if(WillFilter)
00669                             {          
00670                                 ExtMucKal* muckal = new ExtMucKal();
00671                                 muckal->SetGapID(RemID);
00672                                 muckal->SetPosMomErr(RemPositon,RemMomentum,RemXpErr);
00673                                 muckal->SetMucDigiCol(m_mucdigicol);
00674                                 muckal->SetMucWindow(myMucWindow);
00675                                 bool iniOK = muckal->MucKalIniti();
00676                                 vector<MucRecHit*> tmp = muckal->TrackHit();
00677                                 bool filter = muckal->ExtMucFilter();
00678                                 //      bool filter = muckal->GetFilterStatus();
00679                                 double chi2 = muckal->GetChi2();
00680                                 bool samestrip = muckal->GetSameStrip();
00681                                 if(filter&&iniOK)
00682                                 {        
00683 
00684                                     myMucnfit_++;
00685                                     //cout<<"myMucnfit_: "<<myMucnfit_<<endl;
00686                                     myMucchisq_ = myMucchisq_+chi2;
00687                                     muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
00688                                     RememberID = muckal->GetHitGap();
00689                                     if(RememberID[0]==1)myMucbrLastLay_=RememberID[2];
00690                                     if(RememberID[0]==0||RememberID[0]==2)myMucecLastLay_=RememberID[2];
00691                                     if(oldVolumeName.contains("Ab"))
00692                                       RemDepth = RemDepth-depth;
00693                                     if(myMucnfit_==1)
00694                                       myMucdepth_ = RemDepth;
00695                                     else
00696                                       myMucdepth_=myMucdepth_+RemDepth;
00697 
00698                                     if(!samestrip)
00699                                     {
00700 
00701                                         extXpErr->put_err(m_err_mod);
00702                                         extXpErr->set_pos(m_pos_mod);
00703                                         extXpErr->set_mom(m_mom_mod);
00704                                         currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
00705 
00706                                     }
00707                                     else
00708                                     {
00709                                         //cout<<"-------hit exist, but no modify-------"<<endl;
00710                                         RemPositon = currentPosition;
00711                                         RemMomentum = currentMomentum;
00712                                         RemXpErr = XpErr;
00713                                         RemDist = dist_b;
00714                                         RemID =ID_2;
00715                                         if(oldVolumeName.contains("Ab"))
00716                                           RemDepth = depth;
00717                                         else
00718                                           RemDepth = 0.;
00719                                     }
00720                                     if(msgFlag)cout<<"---------Filter is OK---------- "<<endl;
00721                                 }
00722                                 else
00723                                 {
00724                                     // if(LastLay_)
00725                                     //{
00726                                     RemPositon = currentPosition;
00727                                     RemMomentum = currentMomentum;
00728                                     RemXpErr = XpErr;
00729                                     RemDist = dist_b;
00730                                     RemID =ID_2;
00731                                     // }
00732                                 }
00733                                 delete muckal;
00734                             }
00735                         }
00736                         if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
00737                     }
00738 
00739                     /*
00740                        if(Flag1&&Flag2&&Flag3)
00741                        {
00742                        ExtMucKal* muckal = new ExtMucKal();
00743                        muckal->SetVolume(oldVolume,newVolume);
00744                        muckal->SetMomPosErr(currentPosition,currentMomentum,XpErr);
00745                        muckal->SetMucDigiCol(mucdigicol);
00746                        muckal->ExtMucFilter();
00747                        muckal->XPmod(m_pos_mod,m_mom_mod,m_err_mod);
00748                        bool filter = muckal->GetFilterStatus();
00749                        double chi2 = muckal->GetChi2();
00750                        if(filter)
00751                        {
00752                        myMucnfit_++;
00753                        myMucchisq_ = myMucchisq_+chi2;
00754                        currentStep->GetTrack()->SetTrackStatus(fStopAndKill);
00755                        RememberVol.assign(oldVolumeName,10);
00756                        extXpErr->put_err(m_err_mod);
00757                        extXpErr->set_pos(m_pos_mod);
00758                        extXpErr->set_mom(m_mom_mod);
00759                        }
00760                        delete muckal;
00761                        }
00762                        */
00763                     //   if(myExtTrack)myExtTrack->SetMucKalData(myMucchisq_,myMucnfit_,myMucdepth_,myMucbrLastLay_,myMucecLastLay_,myMucnhits_);
00764                     /*          //Get Muc Ext Hits Data.
00765                             if(newVolumeName.contains("Strip"))
00766                             {
00767                             int newVolumeNumber = newVolume->GetCopyNo();
00768                             if(myExtTrack)
00769                             {
00770                             ExtMucHit aExtMucHit;
00771                             Hep3Vector xVector(1.0,0,0); 
00772                             Hep3Vector yVector(0,1.0,0);
00773                             Hep3Vector zVector(0,0,1.0);
00774                             double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00775                             double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00776                             double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00777                             double tzError;
00778                             double phi = currentPosition.phi();
00779                             if(phi<0.) phi+=M_PI;
00780                             int i = int(8.0*phi/M_PI);
00781                             if( i==0 || i==7 || i==8 )
00782                             {
00783                             tzError = yError;
00784                             }
00785                             if( i==1 || i==2 )
00786                             {
00787                             Hep3Vector tzVector(-1.0,1.0,0.);
00788                             tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00789                             }
00790                             if( i==3 || i==4 )
00791                             {
00792                             tzError = xError;
00793                             } 
00794                             if( i==5 || i==6 )
00795                             {
00796                             Hep3Vector tzVector(1.0,1.0,0.);
00797                             tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00798                             }
00799                             aExtMucHit.SetExtMucHit(currentPosition,currentMomentum,newVolumeName,newVolumeNumber,extXpErr->get_err(),zError,tzError,xError,yError);
00800                             myExtTrack->AddExtMucHit(aExtMucHit);                               
00801                             }
00802                             }
00803                             */
00804         }
00805         else {
00806             if(msgFlag) cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<endl;
00807             double x = currentPosition.x();
00808             double y = currentPosition.y();
00809             double z = currentPosition.z();
00810             if( (fabs(x)>=2*myMucR) || (fabs(y)>=2*myMucR) || (fabs(z)>=2*myMucZ) )
00811               //currentTrack->SetKineticEnergy(0.0);// protection in case that a very energetic track travels without touching BESIII                                                                                        
00812             {currentStep->GetTrack()->SetTrackStatus(fStopAndKill);m_trackstop=true;}
00813         }
00814 
00815     }
00816     else if(currentTrackStatus == fStopAndKill)
00817     {
00818         m_trackstop=true;
00819         if(myEmcFlag) myExtTrack->SetEmcPath(myPathInCrystal/10.);
00820         if(myTof1Flag) myExtTrack->setPathInTof1(myPathInTof1);
00821         if(myTof2Flag) myExtTrack->setPathInTof2(myPathInTof2);
00822         if(msgFlag) {
00823             cout << "myPathInTof1 = " ;
00824             for(int i=0; i!=myPathInTof1.size(); i++)
00825               cout << myPathInTof1[i] << "  ";
00826             cout << endl;
00827             cout << "myPathInTof2 = " ;
00828             for(int i=0; i!=myPathInTof2.size(); i++)
00829               cout << myPathInTof2[i] << "  ";
00830             cout << endl;
00831         }
00832 
00833         if(msgFlag)
00834         {
00835             if(newVolume!=0)
00836             {
00837                 std::cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//stoped"<<endl;
00838                 cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
00839             }
00840             else {
00841                 cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//escaped"<<std::endl;
00842                 std::cout<<"Error matrix is:"<<extXpErr->get_err()<<std::endl;
00843             }
00844         }
00845     }
00846 }


Member Data Documentation

double ExtSteppingAction::chicc [private]

Definition at line 61 of file ExtSteppingAction.h.

Referenced by CalculateChicc(), Reset(), and UserSteppingAction().

Ext_xp_err* ExtSteppingAction::extXpErr [private]

Definition at line 80 of file ExtSteppingAction.h.

Referenced by SetXpErrPointer(), and UserSteppingAction().

double ExtSteppingAction::initialPath [private]

Definition at line 62 of file ExtSteppingAction.h.

Referenced by SetInitialPath(), and UserSteppingAction().

double ExtSteppingAction::initialTof [private]

Definition at line 63 of file ExtSteppingAction.h.

Referenced by SetInitialTof(), and UserSteppingAction().

HepSymMatrix ExtSteppingAction::m_err_mod [private]

Definition at line 96 of file ExtSteppingAction.h.

Referenced by InfmodMuc(), and UserSteppingAction().

Hep3Vector ExtSteppingAction::m_mom_mod [private]

Definition at line 98 of file ExtSteppingAction.h.

Referenced by InfmodMuc(), and UserSteppingAction().

MucDigiCol* ExtSteppingAction::m_mucdigicol [private]

Definition at line 86 of file ExtSteppingAction.h.

Referenced by SetMucDigiColPointer(), and UserSteppingAction().

Hep3Vector ExtSteppingAction::m_pos_mod [private]

Definition at line 97 of file ExtSteppingAction.h.

Referenced by InfmodMuc(), and UserSteppingAction().

bool ExtSteppingAction::m_trackstop [private]

Definition at line 89 of file ExtSteppingAction.h.

Referenced by Reset(), TrackStop(), and UserSteppingAction().

int ExtSteppingAction::m_which_tof_version [private]

Definition at line 136 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction(), Get_which_tof_version(), and Set_which_tof_version().

bool ExtSteppingAction::msgFlag [private]

Definition at line 112 of file ExtSteppingAction.h.

Referenced by SetMsgFlag(), and UserSteppingAction().

double ExtSteppingAction::myBetaInMDC [private]

Definition at line 64 of file ExtSteppingAction.h.

Referenced by SetBetaInMDC().

bool ExtSteppingAction::myEmcFlag [private]

Definition at line 121 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myEmcPathFlag [private]

Definition at line 122 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myEmcR1 [private]

Definition at line 129 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

double ExtSteppingAction::myEmcR2 [private]

Definition at line 130 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

double ExtSteppingAction::myEmcZ [private]

Definition at line 131 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

RecExtTrack* ExtSteppingAction::myExtTrack [private]

Definition at line 110 of file ExtSteppingAction.h.

Referenced by Reset(), SetExtTrackPointer(), and UserSteppingAction().

bool ExtSteppingAction::myInTof1 [private]

Definition at line 116 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myInTof2 [private]

Definition at line 118 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

int ExtSteppingAction::myMucbrLastLay_ [private]

Definition at line 93 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myMucchisq_ [private]

Definition at line 91 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myMucdepth_ [private]

Definition at line 92 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

int ExtSteppingAction::myMucecLastLay_ [private]

Definition at line 94 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myMucFlag [private]

Definition at line 123 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

int ExtSteppingAction::myMucnfit_ [private]

Definition at line 90 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

int ExtSteppingAction::myMucnhits_ [private]

Definition at line 95 of file ExtSteppingAction.h.

Referenced by UserSteppingAction().

double ExtSteppingAction::myMucR [private]

Definition at line 133 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction(), and UserSteppingAction().

int ExtSteppingAction::myMucWindow [private]

Definition at line 77 of file ExtSteppingAction.h.

Referenced by SetMucWindow(), and UserSteppingAction().

double ExtSteppingAction::myMucZ [private]

Definition at line 134 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction(), and UserSteppingAction().

HepSymMatrix ExtSteppingAction::myOutputSM [private]

Definition at line 83 of file ExtSteppingAction.h.

Referenced by myOutputSymMatrix().

bool ExtSteppingAction::myOutTof1 [private]

Definition at line 117 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myOutTof2 [private]

Definition at line 119 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathInCrystal [private]

Definition at line 68 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathIntoCrystal [private]

Definition at line 66 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

vector<double> ExtSteppingAction::myPathInTof1 [private]

Definition at line 72 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

vector<double> ExtSteppingAction::myPathInTof2 [private]

Definition at line 76 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathIntoTof1 [private]

Definition at line 70 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathIntoTof2 [private]

Definition at line 74 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathOutCrystal [private]

Definition at line 67 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathOutTof1 [private]

Definition at line 71 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myPathOutTof2 [private]

Definition at line 75 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myPhyEmcFlag [private]

Definition at line 120 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

bool ExtSteppingAction::myTof1Flag [private]

Definition at line 114 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myTof1R [private]

Definition at line 125 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

double ExtSteppingAction::myTof1Z [private]

Definition at line 126 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

bool ExtSteppingAction::myTof2Flag [private]

Definition at line 115 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

double ExtSteppingAction::myTof2R [private]

Definition at line 127 of file ExtSteppingAction.h.

Referenced by ExtSteppingAction().

bool ExtSteppingAction::myTofFlag [private]

Definition at line 113 of file ExtSteppingAction.h.

Referenced by Reset().

bool ExtSteppingAction::myUseMucKalFlag [private]

Definition at line 87 of file ExtSteppingAction.h.

Referenced by SetMucKalFlag(), and UserSteppingAction().

double ExtSteppingAction::RemDepth [private]

Definition at line 105 of file ExtSteppingAction.h.

Referenced by MucReset(), and UserSteppingAction().

double ExtSteppingAction::RemDist [private]

Definition at line 104 of file ExtSteppingAction.h.

Referenced by MucReset(), and UserSteppingAction().

Hep3Vector ExtSteppingAction::RememberID [private]

Definition at line 88 of file ExtSteppingAction.h.

Referenced by Reset(), and UserSteppingAction().

Hep3Vector ExtSteppingAction::RemID [private]

Definition at line 106 of file ExtSteppingAction.h.

Referenced by MucReset(), and UserSteppingAction().

Hep3Vector ExtSteppingAction::RemMomentum [private]

Definition at line 101 of file ExtSteppingAction.h.

Referenced by UserSteppingAction().

Hep3Vector ExtSteppingAction::RemPositon [private]

Definition at line 100 of file ExtSteppingAction.h.

Referenced by UserSteppingAction().

int ExtSteppingAction::RemStep [private]

Definition at line 103 of file ExtSteppingAction.h.

Referenced by MucReset(), and UserSteppingAction().

G4String ExtSteppingAction::RemVol [private]

Definition at line 107 of file ExtSteppingAction.h.

Referenced by MucReset(), and UserSteppingAction().

HepSymMatrix ExtSteppingAction::RemXpErr [private]

Definition at line 102 of file ExtSteppingAction.h.

Referenced by UserSteppingAction().


Generated on Tue Nov 29 23:19:34 2016 for BOSS_7.0.2 by  doxygen 1.4.7