00001
00002 #include "BesMdcSD.hh"
00003 #include "G4HCofThisEvent.hh"
00004 #include "G4Step.hh"
00005 #include "G4ThreeVector.hh"
00006 #include "G4SDManager.hh"
00007 #include "G4UnitsTable.hh"
00008 #include "G4ios.hh"
00009 #include "G4RunManager.hh"
00010 #include "ReadBoostRoot.hh"
00011 #include "G4Svc/IG4Svc.h"
00012 #include "G4Svc/G4Svc.h"
00013 #include "CalibDataSvc/ICalibRootSvc.h"
00014 #include "CalibData/Dedx/DedxCalibData.h"
00015 #include "CalibData/Dedx/DedxSimData.h"
00016 #include "GaudiKernel/DataSvc.h"
00017 #include "TFile.h"
00018 #include "TH1F.h"
00019 #include "TH2D.h"
00020
00021 #include "GaudiKernel/Bootstrap.h"
00022 #include "GaudiKernel/IService.h"
00023 #include "GaudiKernel/Service.h"
00024 #include "GaudiKernel/SmartDataPtr.h"
00025
00026 #include <iostream>
00027
00028
00029 BesMdcSD::BesMdcSD(G4String name)
00030 :BesSensitiveDetector(name)
00031 {
00032 collectionName.insert("BesMdcHitsCollection");
00033 collectionName.insert("BesMdcTruthCollection");
00034
00035 mdcGeoPointer=BesMdcGeoParameter::GetGeo();
00036 mdcCalPointer=new BesMdcCalTransfer;
00037
00038 IMdcGeomSvc* ISvc;
00039 StatusCode sc=Gaudi::svcLocator()->service("MdcGeomSvc", ISvc);
00040 if (!sc.isSuccess())
00041 std::cout<<"BesMdcSD::Could not open MdcGeomSvc"<<std::endl;
00042 mdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc);
00043
00044 IG4Svc* tmpSvc;
00045 sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
00046 if (!sc.isSuccess())
00047 G4cout <<" MdcSD::Error,could not open G4Svc"<<G4endl;
00048 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00049
00050 if(m_G4Svc->GetMdcDedxFlag()==1){
00051 G4cout <<" MdcSD: Use sampled dedx instead of Geant4 value"<<G4endl;
00052 dedxFuncInti();
00053 }
00054
00056
00057
00058 std::string dedxTDSPath = "/Calib/DedxSim";
00059 IDataProviderSvc* pCalibDataSvc;
00060 sc = Gaudi::svcLocator()->service("CalibDataSvc", pCalibDataSvc, true);
00061 if (!sc.isSuccess())
00062 std::cout << "BesMdcSD::Could not open CalibDataSvc" << std::endl;
00063 m_calibDataSvc = dynamic_cast<CalibDataSvc*>(pCalibDataSvc);
00064 if (!sc.isSuccess())
00065 {
00066 std::cout << "Could not get CalibDataSvc"
00067 << std::endl;
00068 }
00069 SmartDataPtr<CalibData::DedxSimData> pDedxSimData(m_calibDataSvc, dedxTDSPath);
00070 m_version = pDedxSimData->getVersion();
00071 m_numDedxHists = pDedxSimData->gethistNo();
00072 m_numBg = pDedxSimData->getRangeNo();
00073 if (m_version == 0) m_numTheta = 10;
00074 else m_numTheta = pDedxSimData->getThetaNo();
00075 m_dedx_hists = new TH1F[m_numDedxHists];
00076 for (G4int i = 0; i < m_numBg; i++)
00077 {
00078 m_bgRange.push_back(pDedxSimData->getRange(i));
00079 }
00080 for (G4int i = 0; i < m_numDedxHists; i++)
00081 {
00082 m_dedx_hists[i] = pDedxSimData->getHist(i);
00083 }
00084
00085
00086 IDedxCurSvc* tmp_dedxCurSvc;
00087 sc = Gaudi::svcLocator()->service("DedxCurSvc", tmp_dedxCurSvc, true);
00088 if (!sc.isSuccess())
00089 {
00090 std::cout << "Could not get DedxCurSvc"
00091 << std::endl;
00092 }
00093 m_pDedxCurSvc = dynamic_cast<DedxCurSvc*>(tmp_dedxCurSvc);
00094
00095 if(m_G4Svc->MdcRootFlag())
00096 {
00097 m_tupleMdc = m_G4Svc->GetTupleMdc();
00098 sc = m_tupleMdc->addItem("betaGamma",m_betaGamma);
00099 sc = m_tupleMdc->addItem("fitval",m_fitval);
00100 sc = m_tupleMdc->addItem("dedx",m_dedx);
00101 sc = m_tupleMdc->addItem("de",m_de);
00102
00103
00104 sc = m_tupleMdc->addItem("charge", m_charge);
00105 sc = m_tupleMdc->addItem("costheta", m_costheta);
00106 }
00107 }
00108
00109 BesMdcSD::~BesMdcSD(){
00110 delete []m_dedx_hists;
00111 }
00112
00113 void BesMdcSD::Initialize(G4HCofThisEvent* HCE)
00114 {
00115 hitsCollection = new BesMdcHitsCollection
00116 (SensitiveDetectorName,collectionName[0]);
00117 static G4int HCID = -1;
00118 if(HCID<0)
00119 { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); }
00120 HCE->AddHitsCollection( HCID, hitsCollection );
00121 G4int i,j;
00122 for(i=0; i<43;i++){
00123 for(j=0;j<288;j++){
00124 hitPointer[i][j]=-1;
00125 truthPointer[i][j]=-1;
00126 }
00127 }
00128 }
00129
00130
00131 void BesMdcSD::BeginOfTruthEvent(const G4Event* evt)
00132 {
00133 truthCollection = new BesMdcHitsCollection
00134 (SensitiveDetectorName,collectionName[1]);
00135
00136 }
00137
00138 void BesMdcSD::EndOfTruthEvent(const G4Event* evt)
00139 { static G4int HLID=-1;
00140 if(HLID<0)
00141 {
00142 HLID = G4SDManager::GetSDMpointer()->
00143 GetCollectionID(collectionName[1]);
00144 }
00145 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00146 HCE->AddHitsCollection(HLID,truthCollection);
00147 }
00148
00149 G4bool BesMdcSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
00150 {
00151 G4Track* gTrack = aStep->GetTrack() ;
00152
00153 G4double globalT=gTrack->GetGlobalTime();
00154 if(isnan(globalT)){
00155 G4cout<<"MdcSD:error, globalT is nan "<<G4endl;
00156 return false;
00157 }
00158 if(globalT > 2000)return false;
00159
00160
00161 G4int charge = gTrack->GetDefinition()->GetPDGCharge();
00162 if (charge == 0) return false;
00163
00164
00165 G4double stepLength=aStep->GetStepLength();
00166 if(stepLength==0){
00167
00168 return false;
00169 }
00170
00171 G4double edep = aStep->GetTotalEnergyDeposit() / stepLength;
00172 if(edep==0.) return false;
00173
00174
00175 G4StepPoint* prePoint = aStep->GetPreStepPoint() ;
00176 G4StepPoint* postPoint = aStep->GetPostStepPoint() ;
00177
00178
00179 G4ThreeVector pointIn = prePoint->GetPosition();
00180 G4ThreeVector pointOut = postPoint->GetPosition();
00181
00182
00183 const G4VTouchable *touchable = prePoint->GetTouchable();
00184 G4VPhysicalVolume *volume = touchable->GetVolume(0);
00185
00186 G4double driftD = 0;
00187 G4double driftT = 0;
00188 G4double edepTemp = 0;
00189 G4double lengthTemp = 0;
00190 G4int cellId=0;
00191 G4int layerId = touchable->GetVolume(1)->GetCopyNo();
00192 if(volume->IsReplicated()){
00193 cellId = touchable->GetReplicaNumber();
00194 }else{
00195 cellId=touchable->GetVolume(0)->GetCopyNo();
00196 }
00197 if(layerId==18&&(cellId==27||cellId==42))return false;
00198
00199 if(ReadBoostRoot::GetMdc() == 2) {
00200
00201
00202 if(layerId >= 36) layerId = 36 + (layerId-36)/2;
00203 }
00204
00205 G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
00206 if(((fabs(pointIn.z())-halfLayerLength)>-7.)
00207 &&((fabs(pointOut.z())-halfLayerLength)>-7.))return false;
00208
00209 G4int trackID= gTrack->GetTrackID();
00210 G4int truthID, g4TrackID;
00211 GetCurrentTrackIndex(truthID, g4TrackID);
00212
00213 G4double theta=gTrack->GetMomentumDirection().theta();
00214
00215 G4ThreeVector hitPosition=0;
00216 G4double transferT=0;
00217 driftD = Distance(layerId,cellId,pointIn,pointOut,hitPosition,transferT);
00218
00219 G4double posPhi, wirePhi;
00220 posPhi=hitPosition.phi();
00221 if(posPhi<0)posPhi += 2*pi;
00222 wirePhi=mdcGeoPointer->SignalWire(layerId,cellId).Phi(hitPosition.z());
00223
00224 G4int posFlag;
00225 if(posPhi<=wirePhi){
00226 posFlag = 0;
00227 }else{
00228 posFlag = 1;
00229 }
00230
00231 if(posPhi < 1. && wirePhi > 5.)posFlag = 1;
00232 if(posPhi > 5. && wirePhi < 1.)posFlag = 0;
00233
00234 G4ThreeVector hitLine=pointOut-pointIn;
00235 G4double enterAngle=hitLine.phi()-wirePhi;
00236 while(enterAngle<-pi/2.)enterAngle+=pi;
00237 while(enterAngle>pi/2.)enterAngle-=pi;
00238
00239 if(m_G4Svc->GetMdcDedxFlag()==1){
00240 G4double betaGamma=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
00241 if(betaGamma<0.01)return false;
00242
00243
00244 G4double eCount=dedxSample(betaGamma, charge, theta);
00245 edep=eCount;
00246 }
00247
00248 BesMdcHit* newHit = new BesMdcHit();
00249 newHit->SetTrackID(truthID);
00250
00251 newHit->SetLayerNo(layerId);
00252 newHit->SetCellNo(cellId);
00253 newHit->SetEdep(edep);
00254 newHit->SetPos(hitPosition);
00255 newHit->SetDriftD(driftD);
00256 newHit->SetTheta(theta);
00257 newHit->SetPosFlag(posFlag);
00258 newHit->SetEnterAngle(enterAngle);
00259
00260
00261 mdcCalPointer->SetHitPointer(newHit);
00262
00263 driftT=mdcCalPointer->D2T(driftD);
00264 globalT+=transferT;
00265 driftT+=globalT;
00266
00267 newHit->SetDriftT (driftT);
00268 newHit->SetGlobalT(globalT);
00269
00270 if (hitPointer[layerId][cellId] == -1) {
00271 hitsCollection->insert(newHit);
00272 G4int NbHits = hitsCollection->entries();
00273 hitPointer[layerId][cellId]=NbHits-1;
00274 }
00275 else
00276 {
00277 G4int pointer=hitPointer[layerId][cellId];
00278 if (g4TrackID == trackID) {
00279 G4double preDriftT=(*hitsCollection)[pointer]->GetDriftT();
00280 }
00281
00282 G4double preDriftT = (*hitsCollection)[pointer]->GetDriftT();
00283 if (driftT < preDriftT) {
00284 (*hitsCollection)[pointer]->SetTrackID(truthID);
00285
00286 (*hitsCollection)[pointer]->SetDriftD(driftD);
00287 (*hitsCollection)[pointer]->SetDriftT(driftT);
00288 (*hitsCollection)[pointer]->SetPos(hitPosition);
00289 (*hitsCollection)[pointer]->SetGlobalT(globalT);
00290 (*hitsCollection)[pointer]->SetTheta(theta);
00291 (*hitsCollection)[pointer]->SetPosFlag(posFlag);
00292 (*hitsCollection)[pointer]->SetEnterAngle(enterAngle);
00293 }
00294
00295 delete newHit;
00296 }
00297
00298
00299 if(truthCollection){
00300 if(g4TrackID==trackID){
00301 G4int pointer=truthPointer[layerId][cellId];
00302 if(pointer==-1){
00303 BesMdcHit* truthHit = new BesMdcHit();
00304 truthHit->SetTrackID (truthID);
00305 truthHit->SetLayerNo(layerId);
00306 truthHit->SetCellNo(cellId);
00307 truthHit->SetEdep (edep);
00308 truthHit->SetPos (hitPosition);
00309 truthHit->SetDriftD (driftD);
00310 truthHit->SetDriftT (driftT);
00311 truthHit->SetGlobalT(globalT);
00312 truthHit->SetTheta(theta);
00313 truthHit->SetPosFlag(posFlag);
00314 truthHit->SetEnterAngle(enterAngle);
00315
00316 truthCollection->insert(truthHit);
00317 G4int NbHits = truthCollection->entries();
00318 truthPointer[layerId][cellId]=NbHits-1;
00319 }
00320 else {
00321 if(truthID == (*truthCollection)[pointer]->GetTrackID()){
00322 G4double preDriftT=(*truthCollection)[pointer]->GetDriftT();
00323 if(driftT<preDriftT){
00324 (*truthCollection)[pointer]->SetDriftD(driftD);
00325 (*truthCollection)[pointer]->SetDriftT(driftT);
00326 (*truthCollection)[pointer]->SetPos(hitPosition);
00327 (*truthCollection)[pointer]->SetGlobalT(globalT);
00328 (*truthCollection)[pointer]->SetTheta(theta);
00329 (*truthCollection)[pointer]->SetPosFlag(posFlag);
00330 (*truthCollection)[pointer]->SetEnterAngle(enterAngle);
00331 }
00332 edepTemp=(*truthCollection)[pointer]->GetEdep();
00333 (*truthCollection)[pointer]->SetEdep(edepTemp+edep);
00334 } else {
00335 BesMdcHit* truthHit = new BesMdcHit();
00336 truthHit->SetTrackID (truthID);
00337 truthHit->SetLayerNo(layerId);
00338 truthHit->SetCellNo(cellId);
00339 truthHit->SetEdep(edep);
00340 truthHit->SetPos(hitPosition);
00341 truthHit->SetDriftD (driftD);
00342 truthHit->SetDriftT (driftT);
00343 truthHit->SetGlobalT(globalT);
00344 truthHit->SetTheta(theta);
00345 truthHit->SetPosFlag(posFlag);
00346 truthHit->SetEnterAngle(enterAngle);
00347
00348 truthCollection->insert(truthHit);
00349 G4int NbHits = truthCollection->entries();
00350 truthPointer[layerId][cellId]=NbHits-1;
00351 }
00352 }
00353 }
00354 }
00355
00356
00357
00358
00359 return true;
00360 }
00361
00362 void BesMdcSD::EndOfEvent(G4HCofThisEvent*)
00363 {
00364 if (verboseLevel>0) {
00365 hitsCollection->PrintAllHits();
00366
00367
00368
00369
00370
00371
00372 }
00373 }
00374
00375 G4double BesMdcSD::Distance(G4int layerId, G4int cellId, G4ThreeVector pointIn, G4ThreeVector pointOut,G4ThreeVector& hitPosition,G4double& transferT)
00376 {
00377
00378
00379
00380
00381
00382
00383
00384
00385 G4double t1,distance,dInOut,dHitIn,dHitOut;
00386
00387 G4ThreeVector east=mdcGeomSvc->Wire(layerId,cellId)->Backward();
00388 G4ThreeVector west=mdcGeomSvc->Wire(layerId,cellId)->Forward();
00389 G4ThreeVector wireLine=east-west;
00390 G4ThreeVector hitLine=pointOut-pointIn;
00391
00392 G4ThreeVector hitXwire=hitLine.cross(wireLine);
00393 G4ThreeVector wire2hit=east-pointOut;
00394
00395
00396 if(hitXwire.mag()==0){
00397 distance=wireLine.cross(wire2hit).mag()/wireLine.mag();
00398 hitPosition=pointIn;
00399 }else{
00400 t1=hitXwire.dot(wire2hit.cross(wireLine))/hitXwire.mag2();
00401 hitPosition=pointOut+t1*hitLine;
00402
00403 dInOut=(pointOut-pointIn).mag();
00404 dHitIn=(hitPosition-pointIn).mag();
00405 dHitOut=(hitPosition-pointOut).mag();
00406 if(dHitIn<=dInOut && dHitOut<=dInOut){
00407 distance=fabs(wire2hit.dot(hitXwire)/hitXwire.mag());
00408 }else if(dHitOut>dHitIn){
00409 distance=wireLine.cross(pointIn-east).mag()/wireLine.mag();
00410 hitPosition=pointIn;
00411 }else{
00412 distance=wireLine.cross(pointOut-east).mag()/wireLine.mag();
00413 hitPosition=pointOut;
00414 }
00415 }
00416
00417
00418 G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
00419 G4double halfWireLength=wireLine.mag()/2.;
00420 G4double transferZ=0;
00421 if(layerId%2==0){
00422 transferZ=halfLayerLength+hitPosition.z();
00423 }else{
00424 transferZ=halfLayerLength-hitPosition.z();
00425 }
00426 if(layerId<8){
00427 transferT=transferZ*halfWireLength/halfLayerLength/220;
00428 }else{
00429 transferT=transferZ*halfWireLength/halfLayerLength/240;
00430 }
00431
00432 return distance;
00433
00434 }
00435
00436 void BesMdcSD::dedxFuncInti(void)
00437 {
00438 dEdE_mylanfunc = new TF1("dEdE_mylanfunc",
00439 "[3]*TMath::Exp([2]*((x[0]-[0])/[1]+TMath::Exp(-1*((x[0]-[0])/[1]))))",
00440 0,
00441 7500);
00442
00443 dEdE_mylanfunc->SetParNames("MPV","Sigma","constant1","constant2");
00444 }
00445
00446 G4double BesMdcSD::dedxSample(G4double betagamma, G4double charge, G4double theta)
00447 {
00448 G4double x = betagamma;
00449 G4double fitval = GetValDedxCurve(x, charge);
00450 if(fitval <= 0)return 0;
00451
00452 G4double random1, random2, dedx1, dedx2, de;
00453 G4double standard1, standard2, beta_temp1, beta_temp2;
00454 G4double dedx = -1;
00455
00456 G4int range_idx, bg_idx, angle_idx, charge_idx, hist_idx;
00457 range_idx = GetBetagammaIndex(betagamma);
00458 angle_idx = GetAngleIndex(theta);
00459 charge_idx = GetChargeIndex(charge);
00460
00461 if (range_idx == -1)
00462 {
00463 while (dedx <= 0)
00464 {
00465 bg_idx = 0;
00466 hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
00467 random1 = m_dedx_hists[hist_idx].GetRandom();
00468 beta_temp1 = (m_bgRange[0] + m_bgRange[1]) / 2;
00469 standard1 = GetValDedxCurve(beta_temp1, charge);
00470 dedx = random1 + fitval - standard1;
00471 }
00472 }
00473 else if (range_idx == m_numBg - 1)
00474 {
00475 while (dedx <= 0)
00476 {
00477 bg_idx = (G4int)(range_idx / 2);
00478 hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
00479 random1 = m_dedx_hists[hist_idx].GetRandom();
00480 beta_temp1 = (m_bgRange[m_numBg - 2] +
00481 m_bgRange[m_numBg - 1]) / 2;
00482 standard1 = GetValDedxCurve(beta_temp1, charge);
00483 dedx = random1 + fitval - standard1;
00484 }
00485 }
00486 else
00487 {
00488
00489 if (range_idx % 2 == 0)
00490 {
00491 while (dedx <= 0)
00492 {
00493 bg_idx = (G4int)(range_idx / 2);
00494 hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
00495 random1 = m_dedx_hists[hist_idx].GetRandom();
00496 beta_temp1 = (m_bgRange[range_idx] +
00497 m_bgRange[range_idx + 1]) / 2;
00498 standard1 = GetValDedxCurve(beta_temp1, charge);
00499 dedx1 = random1 + fitval - standard1;
00500 dedx = dedx1;
00501 }
00502 }
00503
00504
00505 else
00506 {
00507 while (dedx <= 0)
00508 {
00509
00510 beta_temp1 = (m_bgRange[range_idx - 1] +
00511 m_bgRange[range_idx]) / 2;
00512 standard1 = GetValDedxCurve(beta_temp1, charge);
00513
00514
00515 beta_temp2 = (m_bgRange[range_idx + 1] +
00516 m_bgRange[range_idx + 2]) / 2;
00517 standard2 = GetValDedxCurve(beta_temp2, charge);
00518
00519
00520 bg_idx = (G4int)(range_idx / 2);
00521 hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
00522 random1 = m_dedx_hists[hist_idx].GetRandom();
00523
00524
00525 bg_idx++;
00526 hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
00527 random2 = m_dedx_hists[hist_idx].GetRandom();
00528
00529
00530 dedx1 = random1 + fitval - standard1;
00531 dedx2 = random2 + fitval - standard2;
00532 dedx = (dedx2 * (x - m_bgRange[range_idx]) +
00533 dedx1 * (m_bgRange[range_idx + 1] - x)) /
00534 (m_bgRange[range_idx + 1] - m_bgRange[range_idx]);
00535 }
00536 }
00537 }
00538
00539
00540 de= dedx;
00541
00542
00543 if(m_G4Svc->MdcRootFlag())
00544 {
00545 m_betaGamma= x;
00546 m_fitval= fitval;
00547
00548
00549
00550
00551 m_dedx= dedx;
00552 m_de= de;
00553
00554 m_charge = charge;
00555 m_costheta = cos(theta);
00556 m_tupleMdc->write();
00557 }
00558 return de;
00559 }
00560
00561
00562
00563
00564
00565
00566 G4double BesMdcSD::GetValDedxCurve(G4double x, G4double charge)
00567 {
00568 G4int dedxflag = -1;
00569 G4int size = -1;
00570 G4double A = 0.;
00571 G4double B = 0.;
00572 G4double C = 0.;
00573 std::vector<G4double> par;
00574 G4double val;
00575 G4int index = -1;
00576
00577 par.clear();
00578 dedxflag = m_pDedxCurSvc->getCurve(0);
00579 size = m_pDedxCurSvc->getCurveSize();
00580 for (G4int i = 1; i < size; i++) {
00581 par.push_back(m_pDedxCurSvc->getCurve(i));
00582 }
00583
00584 if (x < 4.5)
00585 A = 1;
00586 else if(x < 10)
00587 B = 1;
00588 else
00589 C = 1;
00590
00591 G4double partA = par[0] * pow(sqrt(x * x + 1), par[2]) / pow(x, par[2]) *
00592 (par[1] - par[5] * log(pow(1 / x, par[3]))) -
00593 par[4] + exp(par[6] + par[7] * x);
00594 G4double partB = par[8] * pow(x, 3) + par[9] * pow(x, 2) + par[10] * x + par[11];
00595 G4double partC = - par[12] * log(par[15] + pow(1 / x, par[13])) + par[14];
00596
00597 val = 550 * (A * partA + B * partB + C * partC);
00598 return val;
00599
00600 }
00601
00602
00603
00604
00605
00606
00607 G4int BesMdcSD::GetBetagammaIndex(G4double bg)
00608 {
00609 if (bg < m_bgRange[0]) return -1;
00610 for (G4int i = 0; i < m_numBg - 1; i++)
00611 {
00612 if (bg > m_bgRange[i] && bg < m_bgRange[i + 1])
00613 {
00614 return i;
00615 }
00616 }
00617 if (bg > m_bgRange[m_numBg - 1])
00618 return (m_numBg - 1);
00619 return -1;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628 G4int BesMdcSD::GetAngleIndex(G4double theta)
00629 {
00630 if (m_version == 0) {
00631 if (fabs(cos(theta)) >= 0.93) return 9;
00632 return (G4int)(fabs(cos(theta)) * 10 / 0.93);
00633 }
00634 else {
00635 G4double cos_min = -0.93;
00636 G4double cos_max = 0.93;
00637 G4int nbin = m_numTheta;
00638 G4double cos_step = (cos_max - cos_min) / nbin;
00639
00640 if (cos(theta) < cos_min) return 0;
00641 else if (cos_min < cos(theta) && cos(theta) < cos_max) {
00642 return (G4int)((cos(theta) - cos_min) / cos_step);
00643 }
00644 else return nbin - 1;
00645 }
00646 }
00647
00648
00649
00650
00651
00652
00653 G4int BesMdcSD::GetChargeIndex(G4int charge)
00654 {
00655 if (charge > 0) return 0;
00656 if (charge == 0) return -1;
00657 if (charge < 0) return 1;
00658 return -1;
00659 }
00660