Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

ExtSteppingAction Class Reference

#include <ExtSteppingAction.h>

List of all members.

Public Member Functions

int CalculateEmcEndCopyNb (int num)
int CalculateEmcEndCopyNb (int num)
int CalculateEmcEndPhiNb (int num)
int CalculateEmcEndPhiNb (int num)
void CalculateEmcEndThetaPhi (int npart, int sector, int nb, int &ntheta, int &nphi)
void CalculateEmcEndThetaPhi (int npart, int sector, int nb, int &ntheta, int &nphi)
 ExtSteppingAction ()
 ExtSteppingAction ()
void Reset ()
void Reset ()
void SetBetaInMDC (double aBeta)
void SetBetaInMDC (double aBeta)
void SetExtTrackPointer (RecExtTrack *aExtTrack)
void SetExtTrackPointer (RecExtTrack *aExtTrack)
void SetInitialPath (double aPath)
void SetInitialPath (double aPath)
void SetInitialTof (double aTof)
void SetInitialTof (double aTof)
void SetMsgFlag (bool aMsgFalg)
void SetMsgFlag (bool aMsgFalg)
void SetXpErrPointer (Ext_xp_err *xpErr)
void SetXpErrPointer (Ext_xp_err *xpErr)
void UserSteppingAction (const G4Step *currentStep)
void UserSteppingAction (const G4Step *currentStep)
 ~ExtSteppingAction ()
 ~ExtSteppingAction ()

Private Member Functions

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

Private Attributes

double chicc
Ext_xp_err * extXpErr
Ext_xp_err * extXpErr
double initialPath
double initialTof
bool msgFlag
double myBetaInMDC
bool myEmcFlag
bool myEmcPathFlag
double myEmcR1
double myEmcR2
double myEmcZ
RecExtTrack * myExtTrack
RecExtTrack * myExtTrack
bool myInTof1
bool myInTof2
bool myMucFlag
double myMucR
double myMucZ
HepSymMatrix myOutputSM
bool myOutTof1
bool myOutTof2
double myPathInCrystal
double myPathIntoCrystal
vector< double > myPathInTof1
vector< double > myPathInTof1
vector< double > myPathInTof2
vector< double > myPathInTof2
double myPathIntoTof1
double myPathIntoTof2
double myPathOutCrystal
double myPathOutTof1
double myPathOutTof2
bool myPhyEmcFlag
bool myTof1Flag
double myTof1R
double myTof1Z
bool myTof2Flag
double myTof2R
bool myTofFlag

Constructor & Destructor Documentation

ExtSteppingAction::ExtSteppingAction (   ) 
  

00024                                     :chicc(0.0),initialPath(0.),myPathIntoCrystal(0.),myPathOutCrystal(0.),myPathInCrystal(0.),myBetaInMDC(0.),extXpErr(0),myExtTrack(0),msgFlag(true)
00025 {
00026   myTof1R = 814.001;
00027   myTof1Z = 1330.0;
00028   myTof2R = 871.036;
00029   myEmcR1 = 945.0; 
00030   myEmcR2 = 983.9;
00031   myEmcZ  = 1416.8;
00032   myMucR  = 1700.0-0.01;
00033   myMucZ  = 2050.0-0.01;
00034 }

ExtSteppingAction::~ExtSteppingAction (   ) 
  

00036 {}

ExtSteppingAction::ExtSteppingAction (   ) 
 

ExtSteppingAction::~ExtSteppingAction (   ) 
 

Member Function Documentation

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

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

00495 {       
00496   int n = currentMaterial->GetNumberOfElements();
00497   const double *p = currentMaterial->GetFractionVector();
00498   double density = (currentMaterial->GetDensity())/(g/cm3);
00499   double temp=0.0;
00500   for(int i=0;i<n;i++)
00501   {
00502     const G4Element* mElement = currentMaterial->GetElement(i);
00503     double z = mElement->GetZ();
00504     double a = mElement->GetN();
00505     //  std::cout<<"Fraction: "<<*p<<" z: "<<z<<" a: "<<a<<std::endl;
00506     temp += *(p++)/a*z*(z+1);
00507   }
00508   chicc = 0.39612e-3*std::sqrt(density*temp);           
00509   //   std::cout<<"chicc:"<<chicc<<std::endl;
00510 }

int ExtSteppingAction::CalculateEmcEndCopyNb (  int  num  ) 
 

int ExtSteppingAction::CalculateEmcEndCopyNb (  int  num  ) 
  

00637 {
00638   int copyNb;
00639   switch(num){
00640     case 30:
00641       copyNb = 5;
00642       break;
00643     case 31:
00644       copyNb = 6;
00645       break;
00646     case 32:
00647       copyNb = 14;
00648       break;
00649     case 33:
00650       copyNb = 15;
00651       break;
00652     case 34:
00653       copyNb = 16;
00654       break;
00655     default:
00656       copyNb = num;
00657       break;
00658   }
00659   return copyNb;
00660 }

int ExtSteppingAction::CalculateEmcEndPhiNb (  int  num  ) 
 

int ExtSteppingAction::CalculateEmcEndPhiNb (  int  num  ) 
  

00624 {
00625   if(num==0||num==1) {
00626     return 15-num*8;
00627   } else if(num==2||num==3) {
00628     return 30-num*8;
00629   } else if(num<=9) {
00630     return 17-num;
00631   } else {
00632     return 15-num;
00633   }
00634 }

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

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

00541 {
00542   if((sector>=0)&&(sector<3))
00543     sector+=16;
00544   //if((sector!=7)&&(sector!=15))
00545   {
00546     if((nb>=0)&&(nb<4))
00547     {
00548       ntheta = 0;
00549       nphi = (sector-3)*4+nb;
00550     }
00551     else if((nb>=4)&&(nb<8))
00552     {
00553       ntheta = 1;
00554       nphi = (sector-3)*4+nb-4;
00555     }
00556     else if((nb>=8)&&(nb<13))
00557     {
00558       ntheta = 2;
00559       nphi = (sector-3)*5+nb-8;
00560     }
00561     else if((nb>=13)&&(nb<18))
00562     {
00563       ntheta = 3;
00564       nphi = (sector-3)*5+nb-13;
00565     }
00566     else if((nb>=18)&&(nb<24))
00567     {
00568       ntheta = 4;
00569       nphi = (sector-3)*6+nb-18;
00570     }
00571     else if((nb>=24)&&(nb<30))
00572     {
00573       ntheta = 5;
00574       nphi = (sector-3)*6+nb-24;
00575     }
00576   }
00577 
00578   if(npart==2)
00579   {
00580     switch(ntheta)
00581     {
00582       case 0:
00583         if(nphi<32)
00584           nphi = 31-nphi;
00585         else
00586           nphi = 95-nphi;
00587         break;
00588       case 1:
00589         if(nphi<32)
00590           nphi = 31-nphi;
00591         else
00592           nphi = 95-nphi;
00593         break;
00594       case 2:
00595         if(nphi<40)
00596           nphi = 39-nphi;
00597         else
00598           nphi = 119-nphi;
00599         break;
00600       case 3:
00601         if(nphi<40)
00602           nphi = 39-nphi;
00603         else
00604           nphi = 119-nphi;
00605         break;
00606       case 4:
00607         if(nphi<48)
00608           nphi = 47-nphi;
00609         else
00610           nphi = 143-nphi;
00611         break;
00612       case 5:
00613         if(nphi<48)
00614           nphi = 47-nphi;
00615         else
00616           nphi = 143-nphi;
00617       default:
00618         ;
00619     }
00620   }
00621 }

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

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

00514 {
00515   //  std::cout<<"myOutputSymMatrix:1"<<std::endl;
00516   HepSymMatrix m;
00517   //  std::cout<<"myOutputSymMatrix:2"<<std::endl;
00518   m=m1;
00519   //  std::cout<<"myOutputSymMatrix:3"<<std::endl;
00520   for(int i=0;i<6;i++)
00521   {  for(int j=0;j<=i;j++)
00522     {  if(i<3) {
00523                  m[i][j]=m[i][j]/100;     //mm*mm --> cm*cm
00524                }
00525     else if(j<3) {
00526       m[i][j]=m[i][j]/10000;   //mm*MeV --> cm*GeV
00527     }
00528     else {
00529       m[i][j]=m[i][j]/1000000; //MeV*MeV --> GeV*GeV
00530     }
00531     }   
00532   }
00533   //  std::cout<<"myOutputSymMatrix:4"<<std::endl;
00534   //  std::cout<<"m1="<<m1<<std::endl;
00535   //  std::cout<<"m="<<m<<std::endl;
00536   myOutputSM=m;
00537   return myOutputSM;
00538 }

void ExtSteppingAction::Reset (   ) 
 

void ExtSteppingAction::Reset (   ) 
  

00039 {
00040   chicc        = 0.;
00041   myExtTrack   = 0;
00042   myPathIntoCrystal = 0.;
00043   myPathOutCrystal  = 0.;
00044   myPathInCrystal   = 0.;
00045 
00046   myPathIntoTof1 = 0.0;
00047   myPathOutTof1  = 0.0;
00048   //myPathInTof1   = 0.0;
00049   myPathInTof1.clear();
00050 
00051   myPathIntoTof2 = 0.0;
00052   myPathOutTof2  = 0.0;
00053   //myPathInTof2   = 0.0;
00054   myPathInTof2.clear();
00055   
00056   myTofFlag    = false;
00057   myTof1Flag   = false;
00058   myTof2Flag   = false;
00059   myInTof1     = false;
00060   myInTof2     = false;
00061   myOutTof1    = false;
00062   myOutTof2    = false;
00063   myPhyEmcFlag = false;
00064   myEmcFlag    = false;
00065   myEmcPathFlag= false;
00066   myMucFlag    = false;
00067 }

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

00035 {myBetaInMDC = aBeta;};

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

00035 {myBetaInMDC = aBeta;};

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

00040 {myExtTrack = aExtTrack;};

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

00040 {myExtTrack = aExtTrack;};

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

00033 {initialPath = aPath;};

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

00033 {initialPath = aPath;};

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

00034 {initialTof = aTof;};

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

00034 {initialTof = aTof;};

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

00038 {msgFlag = aMsgFalg;};

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

00038 {msgFlag = aMsgFalg;};

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

00036 {extXpErr = xpErr;};

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

00036 {extXpErr = xpErr;};

void ExtSteppingAction::UserSteppingAction (  const G4Step *  currentStep  ) 
 

void ExtSteppingAction::UserSteppingAction (  const G4Step *  currentStep  ) 
  

00070 {
00071 
00072   //Get track and its status
00073   G4Track* currentTrack = currentStep->GetTrack();
00074   if(!currentTrack) 
00075   {
00076     cout<<"Can't get currentTrack"<<endl;
00077     return;
00078   }
00079   G4TrackStatus currentTrackStatus = currentTrack->GetTrackStatus();
00080 
00081   if(msgFlag) 
00082   {
00083     int stepNumber = currentTrack->GetCurrentStepNumber();
00084     cout<<"STEP "<<stepNumber<<":"<<endl;
00085   }
00086 
00087   //Get current position and momentum
00088   Hep3Vector currentPosition = currentTrack->GetPosition();
00089   Hep3Vector currentMomentum = currentTrack->GetMomentum();
00090 
00091 
00092   //Get current Volume
00093   G4VPhysicalVolume* oldVolume;
00094   G4VPhysicalVolume* newVolume;
00095   if(!currentTrack->GetTouchableHandle()) cout<<"Can't get currentTrack->GetTouchableHandle()"<<endl;
00096   else  oldVolume= currentTrack->GetTouchableHandle()->GetVolume();
00097   if(!currentTrack->GetNextTouchableHandle()) cout<<"Can't get currentTrack->GetNextTouchableHandle()"<<endl;
00098   else  newVolume= currentTrack->GetNextTouchableHandle()->GetVolume();
00099   if(!oldVolume) cout<<"Can't get oldVolume!"<<endl;
00100 
00101 
00102   //if current particle is alive.
00103   if(currentTrackStatus == fAlive)
00104   {
00105     if(oldVolume!=newVolume)//if so, update Error Matrix and chicc.
00106     {
00107       //Get current B field
00108       double currentPoint[3]={currentPosition.x(),currentPosition.y(),currentPosition.z()};
00109       double currentBfield[3]={0.0,0.0,0.0};
00110       Hep3Vector currentB(0.0,0.0,1.0);
00111       if(G4TransportationManager::GetTransportationManager())
00112       {
00113         G4FieldManager* fieldManager=G4TransportationManager::GetTransportationManager()->GetFieldManager();
00114         if(fieldManager)       
00115         {         
00116           if(fieldManager->GetDetectorField())  
00117           {
00118             fieldManager->GetDetectorField()->GetFieldValue(currentPoint,currentBfield);      
00119             currentB.set(currentBfield[0]/tesla,currentBfield[1]/tesla,currentBfield[2]/tesla);
00120             if(msgFlag) cout<<"B:"<<currentB.x()<<","<<currentB.y()<<","<<currentB.z()<<endl;
00121           }
00122         }
00123       }
00124 
00125       //Get chicc
00126       G4Material* oldMaterial = oldVolume->GetLogicalVolume()->GetMaterial();
00127       if(!oldMaterial) std::cout<<"Can't get oldMaterial"<<std::endl;
00128       else CalculateChicc(oldMaterial);
00129 
00130       //Update Error Matrix
00131       if(!(extXpErr->move(currentPosition,currentMomentum,currentB,1,chicc)))
00132         if(msgFlag) cout<<"can not update Error Matrix!!"<<endl;
00133 
00134       //Verbose
00135       if(msgFlag)
00136       {
00137         cout<<"Volume name:"<<newVolume->GetName()<<endl;
00138         cout<<"Volume number:"<<newVolume->GetCopyNo()<<endl;
00139         cout<<"x:"<<currentPoint[0]<<"//y:"<<currentPoint[1]<<"//z:"<<currentPoint[2]
00140           <<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"
00141           <<currentMomentum.z()<<endl;
00142         cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
00143         cout<<"phi:"<<atan(currentPoint[1]/currentPoint[0])<<endl;
00144         Hep3Vector nz(0.,0.,1.);
00145         cout<<"Projected z error:"<<extXpErr->get_plane_err(currentMomentum.unit(),nz)
00146           <<endl;
00147         double x,y,R;
00148         x = currentPoint[0];
00149         y = currentPoint[1];
00150         R = sqrt(x*x+y*y);
00151         Hep3Vector nt(-y/R,x/R,0.);
00152         cout<<"Projected phi error:"<<(extXpErr->get_plane_err(currentMomentum.unit(),nt))/R
00153           <<endl<<endl;
00154       }
00155 
00156       //some often used things
00157       Hep3Vector xVector(1.0,0,0);
00158       Hep3Vector yVector(0,1.0,0);
00159       Hep3Vector zVector(0,0,1.0);
00160       Hep3Vector tzVector;
00161       tzVector.setRThetaPhi(1.0,M_PI/2.0,currentPosition.phi());
00162       double r = currentPosition.perp();
00163       double x = currentPosition.x();
00164       double y = currentPosition.y();
00165       double z = currentPosition.z();
00166       G4String newVolumeName = newVolume->GetName();
00167       G4String oldVolumeName = oldVolume->GetName();
00168       G4StepPoint* postStepPoint = currentStep->GetPostStepPoint();
00169       G4TouchableHistory* theTouchable = (G4TouchableHistory*)(postStepPoint->GetTouchable());
00170       int newVolumeNumber=newVolume->GetCopyNo();
00171       //                int newVolumeNumber=theTouchable->GetReplicaNumber(2);
00172 
00173       //Get PhyTof data
00174       if( (!myTofFlag) && (!myTof1Flag) && newVolumeName.contains("Tof") ) 
00175       {
00176         newVolumeNumber = -2;
00177         double currentTrackLength = currentTrack->GetTrackLength();
00178         double totalTrackLength = currentTrackLength + initialPath;
00179         //double initialTof = initialPath/(myBetaInMDC*299.792458);
00180         //cout<<"initialTof = "<<initialTof<<endl;
00181         //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
00182         double localTime = currentTrack->GetLocalTime();
00183         //cout<<"localTime = "<<localTime<<endl;
00184         double totalTOF = localTime + initialTof;
00185         //cout<<"totalTOF= "<<totalTOF<<endl;
00186 
00187         if(myExtTrack)
00188         {
00189           double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00190           double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00191           double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00192           double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00193           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);
00194           myTofFlag = true;
00195         }
00196         return;
00197       }
00198 
00199       //Get Tof layer1 Ext data
00200       if( (!myTof1Flag) && (newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ) ) 
00201       {
00202         double currentTrackLength = currentTrack->GetTrackLength();
00203         double totalTrackLength = currentTrackLength+initialPath;
00204         //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
00205         double localTime = currentTrack->GetLocalTime();
00206         double totalTOF = localTime + initialTof;
00207         myInTof1 = true;
00208         myPathIntoTof1 = currentTrackLength;
00209         if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
00210         newVolumeNumber=theTouchable->GetReplicaNumber(2);
00211         if(newVolumeName.contains("ScinEc")) {
00212           newVolumeNumber=(95-newVolumeNumber)/2;
00213           newVolumeNumber=newVolumeNumber+176;
00214           if(newVolumeName.contains("East")) newVolumeNumber=newVolumeNumber+48;
00215         }
00216         else newVolumeNumber=(527-newVolumeNumber)/3;
00217 
00218         if(myExtTrack)
00219         {
00220           double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00221           double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00222           double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00223           double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00224           myExtTrack->SetTof1Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00225           myTof1Flag = true;
00226         }
00227         return;
00228       }
00229 
00230       if( myInTof1 && ( oldVolumeName=="logicalScinBr1" || oldVolumeName.contains("ScinEc") ) ) {
00231               myInTof1 = false;
00232               myOutTof1 = true;
00233               myPathOutTof1 = currentTrack->GetTrackLength();
00234               myPathInTof1.push_back(myPathOutTof1 - myPathIntoTof1);
00235               if(msgFlag) cout << "myPathOutTof1 = " << myPathOutTof1 << endl;
00236       }
00237 
00238       if( myOutTof1 && ( newVolumeName=="logicalScinBr1" || newVolumeName.contains("ScinEc") ) ) {
00239               myInTof1 = true;
00240               myOutTof1 = false;
00241               myPathIntoTof1 = currentTrack->GetTrackLength();
00242               if(msgFlag) cout << "myPathIntoTof1 = " << myPathIntoTof1 << endl;
00243       }
00244 
00245       //Get Tof layer2 Ext data
00246       if( (!myTof2Flag) && newVolumeName=="logicalScinBr2" ) 
00247       {
00248         double currentTrackLength = currentTrack->GetTrackLength();
00249         double totalTrackLength = currentTrackLength+initialPath;
00250         //double totalTOF = totalTrackLength/(myBetaInMDC*299.792458);
00251         double localTime = currentTrack->GetLocalTime();
00252         double totalTOF = localTime + initialTof;
00253         newVolumeNumber=(527-theTouchable->GetReplicaNumber(2))/3;
00254 
00255         myInTof2 = true;
00256         myPathIntoTof2 = currentTrackLength;
00257         if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
00258 
00259         if(myExtTrack)
00260         {
00261           double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00262           double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00263           double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00264           double tzError= extXpErr->get_plane_err(currentMomentum.unit(),tzVector);
00265           myExtTrack->SetTof2Data(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,totalTOF,totalTrackLength/10,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00266           myTof2Flag = true;
00267         }
00268         return;
00269       }
00270 
00271       if( myInTof2 && oldVolumeName=="logicalScinBr2" ) {
00272               myInTof2 = false;
00273               myOutTof2 = true;
00274               myPathOutTof2 = currentTrack->GetTrackLength();
00275               myPathInTof2.push_back(myPathOutTof2 - myPathIntoTof2);
00276               if(msgFlag) cout << "myPathOutTof2 = " << myPathOutTof2 << endl;
00277       }
00278 
00279       if( myOutTof2 && newVolumeName=="logicalScinBr2" ) {
00280               myInTof2 = true;
00281               myOutTof2 = false;
00282               myPathIntoTof2 = currentTrack->GetTrackLength();
00283               if(msgFlag) cout << "myPathIntoTof2 = " << myPathIntoTof2 << endl;
00284       }
00285 
00286       //Get PhyEmc data.
00287       if( (!myPhyEmcFlag) && (!myEmcFlag) && (newVolumeName=="EMC" || newVolumeName.contains("BSC") || newVolumeName=="EndPhi"))
00288       {
00289         newVolumeNumber = -2;
00290         if(myExtTrack)
00291         {
00292           //                            myPathIntoCrystal = currentTrack->GetTrackLength();
00293           Hep3Vector nPhi(-y/r,x/r,0.);
00294           double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
00295 
00296           Hep3Vector aPosition = currentPosition;
00297           double R = aPosition.r();
00298           double aTheta = aPosition.theta();
00299           aPosition.setTheta(aTheta + M_PI_2);
00300           double errorTheta;
00301           errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
00302           if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
00303           myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
00304         }
00305         myPhyEmcFlag = true;
00306         return;
00307       }
00308 
00309 
00310       if(!myEmcPathFlag &&  newVolumeName.contains("Crystal")  )
00311       {
00312         myPathIntoCrystal = currentTrack->GetTrackLength();
00313         //                      cout<<"Enter Crystal, path: "<<myPathIntoCrystal<<endl;
00314         myEmcPathFlag = true;
00315       }
00316 
00317       //Get Emc Ext data.
00318       if( (!myEmcFlag) && newVolumeName.contains("Crystal")  )
00319       {
00320         if(myExtTrack)
00321         {
00322           //------------------- record crystal number
00323           int npart,ntheta,nphi;
00324           if(currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName().contains("BSC")) { //barrel
00325             npart=1;
00326             std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(1)->GetName()).substr(16,2));
00327             thetaBuf >> ntheta ;
00328             nphi = 308-currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
00329           } else {  //endcap
00330             npart=2-2*currentTrack->GetNextTouchableHandle()->GetCopyNumber(3);
00331             int endSector=currentTrack->GetNextTouchableHandle()->GetCopyNumber(2);
00332             int endNb,endNbGDML;
00333             std::istringstream thetaBuf((currentTrack->GetNextTouchableHandle()->GetVolume(0)->GetName()).substr(20,2));
00334             thetaBuf  >> endNb ;
00335             endNbGDML=CalculateEmcEndCopyNb(endNb);
00336             int endSectorGDML=CalculateEmcEndPhiNb(endSector);
00337             CalculateEmcEndThetaPhi(npart,endSectorGDML,endNbGDML,ntheta,nphi);
00338           }
00339           ostringstream strEmc;
00340           if(ntheta<10) {
00341             strEmc<<"EmcPart"<<npart<<"Theta0"<<ntheta<<"Phi"<<nphi;
00342           } else {
00343             strEmc<<"EmcPart"<<npart<<"Theta"<<ntheta<<"Phi"<<nphi;
00344           } //------------------------------------------
00345 
00346           //                            myPathIntoCrystal = currentTrack->GetTrackLength();
00347           Hep3Vector nPhi(-y/r,x/r,0.);
00348           double errorPhi = (extXpErr->get_plane_err(currentMomentum.unit(),nPhi))/r;
00349 
00350           Hep3Vector aPosition = currentPosition;
00351           double R = aPosition.r();
00352           double aTheta = aPosition.theta();
00353           aPosition.setTheta(aTheta + M_PI_2);
00354           double errorTheta;
00355           errorTheta =(extXpErr->get_plane_err(currentMomentum.unit(),aPosition.unit()))/R;
00356           if(msgFlag) cout<<"Theta direction: "<<aPosition<<endl;
00357           myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,strEmc.str(),
00358               //myExtTrack->SetEmcData(currentPosition/10,currentMomentum/1000,newVolumeName,
00359             newVolumeNumber,errorTheta,errorPhi,myOutputSymMatrix(extXpErr->get_err()));
00360               }
00361               myEmcFlag = true;
00362               return;
00363               }
00364 
00365               //Get path in Emc
00366               if(myEmcPathFlag &&  oldVolumeName.contains("Crystal") )
00367               {
00368               myPathOutCrystal = currentTrack->GetTrackLength();
00369               myPathInCrystal = myPathInCrystal+myPathOutCrystal-myPathIntoCrystal;
00370               //                        cout<<"Go out of crystal, path: "<<myPathOutCrystal<<endl;
00371               myEmcPathFlag=false;
00372               }
00373 
00374 
00375               //Get Muc Ext Data.
00376               if( (!myMucFlag) && ( (fabs(x)>=myMucR) || (fabs(y)>=myMucR) || ((fabs(x-y)/sqrt(2.))>=myMucR) || ((fabs(x+y)/sqrt(2.))>=myMucR) || (fabs(z)>=myMucZ) ) )
00377               {
00378               int newVolumeNumber = newVolume->GetCopyNo();
00379               if(myExtTrack)
00380               {
00381                 Hep3Vector xVector(1.0,0,0); 
00382                 Hep3Vector yVector(0,1.0,0);
00383                 Hep3Vector zVector(0,0,1.0);
00384                 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00385                 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00386                 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00387                 double tzError;
00388                 double phi = currentPosition.phi();
00389                 if(phi<0.) phi+=M_PI;
00390                 int i = int(8.0*phi/M_PI);
00391                 if( i==0 || i==7 || i==8 )
00392                 {
00393                   tzError = yError;
00394                 }
00395                 if( i==1 || i==2 )
00396                 {
00397                   Hep3Vector tzVector(-1.0,1.0,0.);
00398                   tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00399                 }
00400                 if( i==3 || i==4 )
00401                 {
00402                   tzError = xError;
00403                 } 
00404                 if( i==5 || i==6 )
00405                 {
00406                   Hep3Vector tzVector(1.0,1.0,0.);
00407                   tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00408                 }
00409                 myExtTrack->SetMucData(currentPosition/10,currentMomentum/1000,newVolumeName,newVolumeNumber,myOutputSymMatrix(extXpErr->get_err()),zError/10,tzError/10,xError/10,yError/10);
00410               }
00411               myMucFlag = true;
00412               //currentTrack->SetKineticEnergy(0.0);
00413               return;
00414               }
00415 
00416                         //Get Muc Ext Hits Data.
00417                // if(newVolumeName.contains("Strip"))
00418                 if(newVolumeName.contains("Muc"))
00419                 {
00420                 int newVolumeNumber = newVolume->GetCopyNo();
00421                 if(myExtTrack)
00422                 {
00423                 ExtMucHit aExtMucHit;
00424                 Hep3Vector xVector(1.0,0,0); 
00425                 Hep3Vector yVector(0,1.0,0);
00426                 Hep3Vector zVector(0,0,1.0);
00427                 double xError = extXpErr->get_plane_err(currentMomentum.unit(),xVector);
00428                 double yError = extXpErr->get_plane_err(currentMomentum.unit(),yVector);
00429                 double zError = extXpErr->get_plane_err(currentMomentum.unit(),zVector);
00430                 double tzError;
00431                 double phi = currentPosition.phi();
00432                 if(phi<0.) phi+=M_PI;
00433                 int i = int(8.0*phi/M_PI);
00434                 if( i==0 || i==7 || i==8 )
00435                 {
00436                 tzError = yError;
00437                 }
00438                 if( i==1 || i==2 )
00439                 {
00440                 Hep3Vector tzVector(-1.0,1.0,0.);
00441                 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00442                 }
00443                 if( i==3 || i==4 )
00444                 {
00445                 tzError = xError;
00446                 } 
00447                 if( i==5 || i==6 )
00448                 {
00449                 Hep3Vector tzVector(1.0,1.0,0.);
00450                 tzError =extXpErr->get_plane_err(currentMomentum.unit(),tzVector.unit());
00451                 }
00452                 aExtMucHit.SetExtMucHit(currentPosition,currentMomentum,newVolumeName,newVolumeNumber,extXpErr->get_err(),zError,tzError,xError,yError);
00453                 myExtTrack->AddExtMucHit(aExtMucHit);                           
00454                 }
00455                 }
00456            
00457     }
00458     else if(msgFlag) cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<endl;
00459 
00460   }
00461   else if(currentTrackStatus == fStopAndKill)
00462   {
00463     if(myEmcFlag) myExtTrack->SetEmcPath(myPathInCrystal/10.);
00464     if(myTof1Flag) myExtTrack->setPathInTof1(myPathInTof1);
00465     if(myTof2Flag) myExtTrack->setPathInTof2(myPathInTof2);
00466     if(msgFlag) {
00467             cout << "myPathInTof1 = " ;
00468             for(int i=0; i!=myPathInTof1.size(); i++)
00469                     cout << myPathInTof1[i] << "  ";
00470             cout << endl;
00471             cout << "myPathInTof2 = " ;
00472             for(int i=0; i!=myPathInTof2.size(); i++)
00473                     cout << myPathInTof2[i] << "  ";
00474             cout << endl;
00475     }
00476 
00477     if(msgFlag)
00478     {
00479       if(newVolume!=0)
00480       {
00481         std::cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//stoped"<<endl;
00482         cout<<"Error matrix is:"<<extXpErr->get_err()<<endl;
00483       }
00484       else {
00485         cout<<"x:"<<currentPosition.x()<<"//y:"<<currentPosition.y()<<"//z:"<<currentPosition.z()<<"||px:"<<currentMomentum.x()<<"//py:"<<currentMomentum.y()<<"//pz:"<<currentMomentum.z()<<"//escaped"<<std::endl;
00486         std::cout<<"Error matrix is:"<<extXpErr->get_err()<<std::endl;
00487       }
00488     }
00489   }
00490 }

Member Data Documentation

double ExtSteppingAction::chicc [private]
 

Ext_xp_err* ExtSteppingAction::extXpErr [private]
 

Ext_xp_err* ExtSteppingAction::extXpErr [private]
 

double ExtSteppingAction::initialPath [private]
 

double ExtSteppingAction::initialTof [private]
 

bool ExtSteppingAction::msgFlag [private]
 

double ExtSteppingAction::myBetaInMDC [private]
 

bool ExtSteppingAction::myEmcFlag [private]
 

bool ExtSteppingAction::myEmcPathFlag [private]
 

double ExtSteppingAction::myEmcR1 [private]
 

double ExtSteppingAction::myEmcR2 [private]
 

double ExtSteppingAction::myEmcZ [private]
 

RecExtTrack* ExtSteppingAction::myExtTrack [private]
 

RecExtTrack* ExtSteppingAction::myExtTrack [private]
 

bool ExtSteppingAction::myInTof1 [private]
 

bool ExtSteppingAction::myInTof2 [private]
 

bool ExtSteppingAction::myMucFlag [private]
 

double ExtSteppingAction::myMucR [private]
 

double ExtSteppingAction::myMucZ [private]
 

HepSymMatrix ExtSteppingAction::myOutputSM [private]
 

bool ExtSteppingAction::myOutTof1 [private]
 

bool ExtSteppingAction::myOutTof2 [private]
 

double ExtSteppingAction::myPathInCrystal [private]
 

double ExtSteppingAction::myPathIntoCrystal [private]
 

vector<double> ExtSteppingAction::myPathInTof1 [private]
 

vector<double> ExtSteppingAction::myPathInTof1 [private]
 

vector<double> ExtSteppingAction::myPathInTof2 [private]
 

vector<double> ExtSteppingAction::myPathInTof2 [private]
 

double ExtSteppingAction::myPathIntoTof1 [private]
 

double ExtSteppingAction::myPathIntoTof2 [private]
 

double ExtSteppingAction::myPathOutCrystal [private]
 

double ExtSteppingAction::myPathOutTof1 [private]
 

double ExtSteppingAction::myPathOutTof2 [private]
 

bool ExtSteppingAction::myPhyEmcFlag [private]
 

bool ExtSteppingAction::myTof1Flag [private]
 

double ExtSteppingAction::myTof1R [private]
 

double ExtSteppingAction::myTof1Z [private]
 

bool ExtSteppingAction::myTof2Flag [private]
 

double ExtSteppingAction::myTof2R [private]
 

bool ExtSteppingAction::myTofFlag [private]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:14:56 2011 for BOSS6.5.5 by  doxygen 1.3.9.1