BOSS6.5.5: BesMdcSD Class Reference

00031 :BesSensitiveDetector(name)
00032 {
00033     collectionName.insert("BesMdcHitsCollection");
00034     collectionName.insert("BesMdcTruthCollection");
00035   
00036     mdcGeoPointer=BesMdcGeoParameter::GetGeo();
00037     mdcCalPointer=new BesMdcCalTransfer;
00038   
00039     IMdcGeomSvc* ISvc;
00040     StatusCode sc=Gaudi::svcLocator()->service("MdcGeomSvc", ISvc);
00041     if (!sc.isSuccess()) 
00042       std::cout<<"BesMdcSD::Could not open MdcGeomSvc"<<std::endl;
00043     mdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc);
00044   
00045     IG4Svc* tmpSvc;
00046     sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
00047     if (!sc.isSuccess())
00048       G4cout <<" MdcSD::Error,could not open G4Svc"<<G4endl;
00049     m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc); 
00050   
00051     if(m_G4Svc->GetMdcDedxFlag()==1){
00052       G4cout <<" MdcSD: Use sampled dedx instead of Geant4 value"<<G4endl;
00053       dedxFuncInti();
00054     }
00055 
00057     
00058     //get DedxSimData
00059     std::string dedxTDSPath = "/Calib/DedxSim";
00060     IDataProviderSvc* pCalibDataSvc;
00061     sc = Gaudi::svcLocator()->service("CalibDataSvc", pCalibDataSvc, true);
00062     if (!sc.isSuccess())
00063         std::cout << "BesMdcSD::Could not open CalibDataSvc" << std::endl;
00064     m_calibDataSvc = dynamic_cast<CalibDataSvc*>(pCalibDataSvc);
00065     if (!sc.isSuccess())
00066     {
00067         std::cout << "Could not get CalibDataSvc"
00068                   << std::endl;
00069     }
00070     SmartDataPtr<CalibData::DedxSimData> pDedxSimData(m_calibDataSvc, dedxTDSPath);
00071     m_numDedxHists = pDedxSimData->gethistNo();
00072     m_numBg = pDedxSimData->getRangeNo();
00073     m_dedx_hists = new TH1F[m_numDedxHists];
00074     for (G4int i = 0; i < m_numBg; i++)
00075     {
00076         m_bgRange.push_back(pDedxSimData->getRange(i));
00077     }
00078     for (G4int i = 0; i < m_numDedxHists; i++)
00079     {
00080         m_dedx_hists[i] = pDedxSimData->getHist(i);
00081     }
00082 
00083     //get CalibCurSvc
00084     IDedxCurSvc* tmp_dedxCurSvc;
00085     sc = Gaudi::svcLocator()->service("DedxCurSvc", tmp_dedxCurSvc, true);
00086     if (!sc.isSuccess())
00087     {
00088         std::cout << "Could not get DedxCurSvc"
00089                   << std::endl;
00090     }
00091     m_pDedxCurSvc = dynamic_cast<DedxCurSvc*>(tmp_dedxCurSvc);
00092 
00093     if(m_G4Svc->MdcRootFlag())
00094     {
00095       m_tupleMdc = m_G4Svc->GetTupleMdc();
00096       sc = m_tupleMdc->addItem("betaGamma",m_betaGamma);
00097       sc = m_tupleMdc->addItem("fitval",m_fitval);
00098       sc = m_tupleMdc->addItem("dedx",m_dedx);
00099       sc = m_tupleMdc->addItem("de",m_de);
00100       //sc = m_tupleMdc->addItem("length",m_length);
00101       //sc = m_tupleMdc->addItem("random",m_random);
00102       sc = m_tupleMdc->addItem("charge", m_charge);
00103       sc = m_tupleMdc->addItem("costheta", m_costheta);
00104     }
00105     mdc_hit_length.clear();
00106 }

00108                    {
00109   mdc_hit_length.clear();
00110   delete []m_dedx_hists;
00111 }

00132 {  
00133   truthCollection = new BesMdcHitsCollection
00134                           (SensitiveDetectorName,collectionName[1]); 
00135   //  G4cout<<" begin event "<<evt->GetEventID()<<G4endl;
00136 }

00444 {
00445     dEdE_mylanfunc = new TF1("dEdE_mylanfunc",
00446                              "[3]*TMath::Exp([2]*((x[0]-[0])/[1]+TMath::Exp(-1*((x[0]-[0])/[1]))))",
00447                              0,
00448                              7500);
00449     //dEdE_mylanfunc->SetParameters(2009.35,559.776,-1.0932,6327.38);
00450     dEdE_mylanfunc->SetParNames("MPV","Sigma","constant1","constant2");
00451 }

00454 {
00455     G4double x = betagamma;
00456     G4double fitval = GetValDedxCurve(x, charge);
00457     if(fitval <= 0)return 0;
00458    
00459     G4double random1, random2, dedx1, dedx2, de;
00460     G4double standard1, standard2, beta_temp1, beta_temp2;
00461     G4double dedx = -1;
00462 
00463     G4int range_idx, bg_idx, angle_idx, charge_idx, hist_idx;
00464     range_idx = GetBetagammaIndex(betagamma);
00465     angle_idx = GetAngleIndex(theta);
00466     charge_idx = GetChargeIndex(charge);
00467     
00468     if (range_idx == -1)
00469     {
00470         while (dedx <= 0)
00471         {
00472             bg_idx = 0;
00473             hist_idx = bg_idx * 2 * 10 + angle_idx * 2 + charge_idx;
00474             random1 = m_dedx_hists[hist_idx].GetRandom();
00475             beta_temp1 = (m_bgRange[0] + m_bgRange[1]) / 2;
00476             standard1 = GetValDedxCurve(beta_temp1, charge); 
00477             dedx = random1 + fitval - standard1;
00478         }
00479     }
00480     else if (range_idx == m_numBg - 1)
00481     {
00482         while (dedx <= 0)
00483         {
00484             bg_idx = (G4int)(range_idx / 2);
00485             hist_idx = bg_idx * 2 * 10 + angle_idx * 2 + charge_idx;
00486             random1 = m_dedx_hists[hist_idx].GetRandom();
00487             beta_temp1 = (m_bgRange[m_numBg - 2] +
00488                           m_bgRange[m_numBg - 1]) / 2;
00489             standard1 = GetValDedxCurve(beta_temp1, charge);
00490             dedx = random1 + fitval - standard1;
00491         }     
00492     }
00493     else
00494     {
00495         //case 1: Given betagamma fall in one histograph range
00496         if (range_idx % 2 == 0)
00497         {
00498             while (dedx <= 0)
00499             {
00500                 bg_idx = (G4int)(range_idx / 2);
00501                 hist_idx = bg_idx * 2 * 10 + angle_idx * 2 + charge_idx;
00502                 random1 = m_dedx_hists[hist_idx].GetRandom();
00503                 beta_temp1 = (m_bgRange[range_idx] + 
00504                               m_bgRange[range_idx + 1]) / 2;
00505                 standard1 = GetValDedxCurve(beta_temp1, charge); 
00506                 dedx1 = random1 + fitval - standard1;
00507                 dedx = dedx1;
00508             }
00509         }
00510         //case 2: Given betagamma fall in interval between
00511         //        two histographs.
00512         else
00513         {
00514             while (dedx <= 0)
00515             {
00516                 //standard1
00517                 beta_temp1 = (m_bgRange[range_idx - 1] + 
00518                               m_bgRange[range_idx]) / 2;
00519                 standard1 = GetValDedxCurve(beta_temp1, charge);
00520 
00521                 //stardard2
00522                 beta_temp2 = (m_bgRange[range_idx + 1] + 
00523                               m_bgRange[range_idx + 2]) / 2;
00524                 standard2 = GetValDedxCurve(beta_temp2, charge);
00525                
00526                 //random1 
00527                 bg_idx = (G4int)(range_idx / 2);
00528                 hist_idx = bg_idx * 2 * 10 + angle_idx * 2 + charge_idx;
00529                 random1 = m_dedx_hists[hist_idx].GetRandom();
00530 
00531                 //random2
00532                 bg_idx++;
00533                 hist_idx = bg_idx * 2 * 10 + angle_idx * 2 + charge_idx;
00534                 random2 = m_dedx_hists[hist_idx].GetRandom();
00535 
00536                 //combine dedx1 and dedx2
00537                 dedx1 = random1 + fitval - standard1;
00538                 dedx2 = random2 + fitval - standard2;
00539                 dedx = (dedx2 * (x - m_bgRange[range_idx]) + 
00540                         dedx1 * (m_bgRange[range_idx + 1] - x)) / 
00541                        (m_bgRange[range_idx + 1] - m_bgRange[range_idx]);
00542             }
00543         }
00544     }
00545    
00546    
00547     de= dedx;// * y/10./1.5;// stepLength unit is mm in Geant4, cm in Rec.& Cal. software 
00548                          // dedx counts *1.5 in dedx Cal.
00549    
00550     if(m_G4Svc->MdcRootFlag())
00551     {
00552        m_betaGamma= x;
00553        m_fitval= fitval;
00554        //m_trackId = trackId; 
00555        //m_layer = layerId; 
00556        //m_wire  = cellId; 
00557        //m_random= random;
00558        m_dedx= dedx;
00559        m_de= de;
00560        //m_length=y;
00561        m_charge = charge;
00562        m_costheta = cos(theta);
00563        m_tupleMdc->write(); 
00564     }
00565     return de;
00566 }

00383 {
00384   //For two lines r=r1+t1.v1 & r=r2+t2.v2 
00385   //the closest approach is d=|(r2-r1).(v1 x v2)|/|v1 x v2|
00386   //the point where closest approach are
00387   //t1=(v1 x v2).[(r2-r1) x v2]/[(v1 x v2).(v1 x v2)]
00388   //t2=(v1 x v2).[(r2-r1) x v1]/[(v1 x v2).(v1 x v2)]
00389   //if v1 x v2=0 means two lines are parallel
00390   //d=|(r2-r1) x v1|/|v1| 
00391 
00392   G4double t1,distance,dInOut,dHitIn,dHitOut;
00393   //Get wirepoint @ endplate
00394   G4ThreeVector east=mdcGeomSvc->Wire(layerId,cellId)->Backward();
00395   G4ThreeVector west=mdcGeomSvc->Wire(layerId,cellId)->Forward();
00396   G4ThreeVector wireLine=east-west;
00397   G4ThreeVector hitLine=pointOut-pointIn;
00398 
00399   G4ThreeVector hitXwire=hitLine.cross(wireLine);
00400   G4ThreeVector wire2hit=east-pointOut;
00401   //Hitposition is the position on hit line where closest approach 
00402   //of two lines, but it may out the area from pointIn to pointOut
00403   if(hitXwire.mag()==0){
00404     distance=wireLine.cross(wire2hit).mag()/wireLine.mag();
00405     hitPosition=pointIn;
00406   }else{
00407     t1=hitXwire.dot(wire2hit.cross(wireLine))/hitXwire.mag2();
00408     hitPosition=pointOut+t1*hitLine;
00409     
00410     dInOut=(pointOut-pointIn).mag();
00411     dHitIn=(hitPosition-pointIn).mag();
00412     dHitOut=(hitPosition-pointOut).mag();
00413     if(dHitIn<=dInOut && dHitOut<=dInOut){ //Between point in & out
00414       distance=fabs(wire2hit.dot(hitXwire)/hitXwire.mag());
00415     }else if(dHitOut>dHitIn){ // out pointIn
00416       distance=wireLine.cross(pointIn-east).mag()/wireLine.mag();
00417       hitPosition=pointIn;
00418     }else{ // out pointOut
00419       distance=wireLine.cross(pointOut-east).mag()/wireLine.mag(); 
00420       hitPosition=pointOut;
00421     }
00422   }
00423 
00424   //Calculate signal transferT on wire
00425   G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
00426   G4double halfWireLength=wireLine.mag()/2.;
00427   G4double transferZ=0;
00428   if(layerId%2==0){
00429     transferZ=halfLayerLength+hitPosition.z(); //West readout
00430   }else{
00431     transferZ=halfLayerLength-hitPosition.z(); //East readout
00432   }
00433   if(layerId<8){
00434     transferT=transferZ*halfWireLength/halfLayerLength/220;
00435   }else{
00436     transferT=transferZ*halfWireLength/halfLayerLength/240;
00437   }
00438 
00439   return distance;
00440 
00441 }

00370 {
00371   if (verboseLevel>0) { 
00372   hitsCollection->PrintAllHits();
00373   /*
00374      G4int NbHits = hitsCollection->entries();
00375      G4cout << "\n-------->Hits Collection: in this event they are " << NbHits 
00376             << " hits in the MDC chambers: " << G4endl;
00377      for (G4int i=0;i<NbHits;i++) (*hitsCollection)[i]->Print();
00378   */
00379   } 
00380 }

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 }

00634 {
00635     if (fabs(cos(theta)) >= 0.93) return 9;
00636     return (G4int)(fabs(cos(theta)) * 10 / 0.93);
00637 }

00615 {
00616     if (bg < m_bgRange[0]) return -1;
00617     for (G4int i = 0; i < m_numBg - 1; i++)
00618     {
00619         if (bg > m_bgRange[i] && bg < m_bgRange[i + 1])
00620         {
00621             return i;
00622         }
00623     }
00624     if (bg > m_bgRange[m_numBg - 1])
00625         return (m_numBg - 1);
00626 }

00645 {
00646     if (charge > 0)  return 0;
00647     if (charge == 0) return -1;  // warning: -1 is illegal, for charged tracks are expected.
00648     if (charge < 0)  return 1;
00649 }

00031 {
00032   BesSensitiveManager* manager = BesSensitiveManager::GetSensitiveManager();
00033   if(manager)
00034   {
00035     trackIndex = manager->GetCurrentTrackIndex();
00036     std::vector<BesTruthTrack*>* trackList = manager->GetTrackList();
00037     if(trackList)
00038     {
00039       G4int size = trackList->size();
00040       if(size>0)
00041       {
00042         for(G4int i=0;i<size;i++)
00043         {
00044           if( (*trackList)[i]->GetIndex() == trackIndex )
00045           {
00046             g4TrackId = (*trackList)[i]->GetG4TrackId();
00047             break;
00048           }
00049         }
00050       } 
00051     }
00052   }
00053 }

00574 {
00575     G4int dedxflag = -1;
00576     G4int size = -1;
00577     G4double A = 0.;
00578     G4double B = 0.;
00579     G4double C = 0.;
00580     std::vector<G4double> par;
00581     G4double val;
00582     G4int index = -1;
00583 
00584     par.clear();
00585     dedxflag = m_pDedxCurSvc->getCurve(0);
00586     size = m_pDedxCurSvc->getCurveSize();
00587     for (G4int i = 1; i < size; i++) {
00588         par.push_back(m_pDedxCurSvc->getCurve(i));
00589     }
00590 
00591     if (x < 4.5)
00592         A = 1;
00593     else if(x < 10)
00594         B = 1;
00595     else
00596         C = 1;
00597     
00598     G4double partA = par[0] * pow(sqrt(x * x + 1), par[2]) / pow(x, par[2]) *
00599                   (par[1] - par[5] * log(pow(1 / x, par[3]))) -
00600                    par[4] + exp(par[6] + par[7] * x);
00601     G4double partB = par[8] * pow(x, 3) + par[9] * pow(x, 2) + par[10] * x + par[11];
00602     G4double partC = - par[12] * log(par[15] + pow(1 / x, par[13])) + par[14];
00603     
00604     val = 550 * (A * partA + B * partB + C * partC);
00605     return val;
00606     
00607 }

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 }

00150 {
00151   G4Track* gTrack = aStep->GetTrack() ;
00152 
00153   G4double globalT=gTrack->GetGlobalTime();//Time since the event in which the track belongs is created
00154   if(isnan(globalT)){
00155     G4cout<<"MdcSD:error, globalT is nan "<<G4endl;
00156     return false;
00157   }
00158   if(globalT > 2000)return false; //MDC T window is 2 microsecond 
00159 
00160   //skip neutral tracks
00161   G4int charge = gTrack->GetDefinition()->GetPDGCharge();
00162   if (charge == 0) return false;
00163 
00164   //skip no energy deposit tracks
00165   G4double stepLength=aStep->GetStepLength();
00166   if(stepLength==0){
00167     //    G4cout<<"step length equal 0!!"<<G4endl;
00168     return false;
00169   }
00170 
00171   G4double edep = aStep->GetTotalEnergyDeposit() / stepLength;
00172   if(edep==0.) return false;
00173 
00174   // get position of the track at the beginning and at the end of step
00175   G4StepPoint* prePoint  = aStep->GetPreStepPoint() ;
00176   G4StepPoint* postPoint = aStep->GetPostStepPoint() ;
00177 
00178   //get position coordinate
00179   G4ThreeVector pointIn  = prePoint->GetPosition();
00180   G4ThreeVector pointOut = postPoint->GetPosition();
00181 
00182   // get physical volumes
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; // no sense wire 
00198 
00199   if(ReadBoostRoot::GetMdc() == 2) { //gdml
00200     //layerId  0-35 -> CopyNo 0-35 in gdml 
00201     //layerId 36-42 -> CopyNo (36,37),(38,39),...(48,49)    
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;//Out sensitive area
00208 
00209   G4int trackID= gTrack->GetTrackID(); //G4 track ID of current track.
00210   G4int truthID, g4TrackID;
00211   GetCurrentTrackIndex(truthID, g4TrackID); //ID of current primary track.
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();//from -pi to pi
00221   if(posPhi<0)posPhi += 2*pi;
00222   wirePhi=mdcGeoPointer->SignalWire(layerId,cellId).Phi(hitPosition.z());//from 0 to 2pi
00223 
00224   G4int posFlag;
00225   if(posPhi<=wirePhi){
00226     posFlag    = 0;
00227   }else{
00228     posFlag    = 1;
00229   }
00230   // if x axis is between pos and wire, phi will has a jump of one of them.
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;//too low momentum
00242      //if (betaGamma < 10.0) betaGamma = 10.0;
00243 
00244      G4double eCount=dedxSample(betaGamma, charge, theta);
00245      edep=eCount;
00246   }
00247 
00248   BesMdcHit* newHit = new BesMdcHit();
00249   newHit->SetTrackID(truthID);
00250   //newHit->SetTrkID(trackID); 
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   //Transfer hit pointer to BesMdcCalTransfer
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       mdc_hit_length.push_back(stepLength);
00273       G4int NbHits = hitsCollection->entries();
00274       hitPointer[layerId][cellId]=NbHits-1;
00275   } 
00276   else
00277   {
00278       G4int pointer=hitPointer[layerId][cellId]; 
00279       if (g4TrackID == trackID) {
00280           //get edep and length of previous hits
00281           edepTemp = (*hitsCollection)[pointer]->GetEdep();
00282           lengthTemp = mdc_hit_length[pointer];
00283           //set edep and length to current hit
00284           (*hitsCollection)[pointer]->SetEdep((edepTemp * lengthTemp + edep * stepLength) / (lengthTemp + stepLength));
00285           mdc_hit_length[pointer] += stepLength;
00286           G4double preDriftT=(*hitsCollection)[pointer]->GetDriftT();
00287       }
00288       
00289       G4double preDriftT = (*hitsCollection)[pointer]->GetDriftT();
00290       if (driftT < preDriftT) {   
00291           (*hitsCollection)[pointer]->SetTrackID(truthID);
00292           //(*hitsCollection)[pointer]->SetTrkID(trackID);
00293           (*hitsCollection)[pointer]->SetDriftD(driftD);
00294           (*hitsCollection)[pointer]->SetDriftT(driftT);
00295           (*hitsCollection)[pointer]->SetPos(hitPosition);
00296           (*hitsCollection)[pointer]->SetGlobalT(globalT);
00297           (*hitsCollection)[pointer]->SetTheta(theta);
00298           (*hitsCollection)[pointer]->SetPosFlag(posFlag);
00299           (*hitsCollection)[pointer]->SetEnterAngle(enterAngle);
00300       }    
00301    
00302       delete newHit;
00303   }
00304 
00305   //for mc truth 
00306   if(truthCollection){
00307     if(g4TrackID==trackID){ //This track is the primary track & will appear in MC truth
00308       G4int pointer=truthPointer[layerId][cellId];  
00309       if(pointer==-1){
00310         BesMdcHit* truthHit = new BesMdcHit();
00311         truthHit->SetTrackID  (truthID);
00312         truthHit->SetLayerNo(layerId);
00313         truthHit->SetCellNo(cellId);
00314         truthHit->SetEdep     (edep);
00315         truthHit->SetPos      (hitPosition);
00316         truthHit->SetDriftD (driftD);
00317         truthHit->SetDriftT (driftT);
00318         truthHit->SetGlobalT(globalT);
00319         truthHit->SetTheta(theta);
00320         truthHit->SetPosFlag(posFlag);
00321         truthHit->SetEnterAngle(enterAngle);
00322 
00323         truthCollection->insert(truthHit);
00324         G4int NbHits = truthCollection->entries();
00325         truthPointer[layerId][cellId]=NbHits-1;
00326       } 
00327       else {
00328         if(truthID == (*truthCollection)[pointer]->GetTrackID()){
00329             G4double preDriftT=(*truthCollection)[pointer]->GetDriftT();
00330             if(driftT<preDriftT){
00331                 (*truthCollection)[pointer]->SetDriftD(driftD);
00332                 (*truthCollection)[pointer]->SetDriftT(driftT);
00333                 (*truthCollection)[pointer]->SetPos(hitPosition);
00334                 (*truthCollection)[pointer]->SetGlobalT(globalT);
00335                 (*truthCollection)[pointer]->SetTheta(theta);
00336                 (*truthCollection)[pointer]->SetPosFlag(posFlag);
00337                 (*truthCollection)[pointer]->SetEnterAngle(enterAngle);
00338                 }
00339             edepTemp=(*truthCollection)[pointer]->GetEdep();
00340             (*truthCollection)[pointer]->SetEdep(edepTemp+edep);
00341             } else {
00342                 BesMdcHit* truthHit = new BesMdcHit();
00343                 truthHit->SetTrackID  (truthID);
00344                 truthHit->SetLayerNo(layerId);
00345                 truthHit->SetCellNo(cellId);
00346                 truthHit->SetEdep(edep);
00347                 truthHit->SetPos(hitPosition);
00348                 truthHit->SetDriftD (driftD);
00349                 truthHit->SetDriftT (driftT);
00350                 truthHit->SetGlobalT(globalT);
00351                 truthHit->SetTheta(theta);
00352                 truthHit->SetPosFlag(posFlag);
00353                 truthHit->SetEnterAngle(enterAngle);
00354 
00355                 truthCollection->insert(truthHit);
00356                 G4int NbHits = truthCollection->entries();
00357                 truthPointer[layerId][cellId]=NbHits-1;
00358             }
00359       }
00360     }
00361   }
00362 
00363   //newHit->Print();
00364 //  newHit->Draw();
00365 
00366   return true;
00367 }


Public Member Functions
virtual void	BeginOfTrack (const G4Track *)
virtual void	BeginOfTrack (const G4Track *)
void	BeginOfTruthEvent (const G4Event *)
void	BeginOfTruthEvent (const G4Event *)
	BesMdcSD (G4String)
	BesMdcSD (G4String)
void	dedxFuncInti (void)
void	dedxFuncInti (void)
G4double	Distance (G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
G4double	Distance (G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
void	EndOfEvent (G4HCofThisEvent *)
void	EndOfEvent (G4HCofThisEvent *)
virtual void	EndOfTrack (const G4Track *)
virtual void	EndOfTrack (const G4Track *)
void	EndOfTruthEvent (const G4Event *)
void	EndOfTruthEvent (const G4Event *)
void	Initialize (G4HCofThisEvent *)
void	Initialize (G4HCofThisEvent *)
G4bool	ProcessHits (G4Step , G4TouchableHistory )
G4bool	ProcessHits (G4Step , G4TouchableHistory )
	~BesMdcSD ()
	~BesMdcSD ()
Protected Member Functions
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
Private Member Functions
G4double	dedxSample (G4double betagamma, G4double length, G4double theta)
G4double	dedxSample (G4double betagamma, G4double length, G4double theta)
G4int	GetAngleIndex (G4double)
G4int	GetAngleIndex (G4double)
G4int	GetBetagammaIndex (G4double bg)
G4int	GetBetagammaIndex (G4double bg)
G4int	GetChargeIndex (G4int)
G4int	GetChargeIndex (G4int)
G4double	GetValDedxCurve (G4double bg, G4double charge)
G4double	GetValDedxCurve (G4double bg, G4double charge)
Private Attributes
TF1 *	dEdE_mylanfunc
TF1 *	dEdE_mylanfunc
G4int	hitPointer [43][288]
BesMdcHitsCollection *	hitsCollection
BesMdcHitsCollection *	hitsCollection
NTuple::Item< double >	m_betaGamma
NTuple::Item< double >	m_betaGamma
std::vector< G4double >	m_bgRange
std::vector< G4double >	m_bgRange
CalibDataSvc *	m_calibDataSvc
	dedx sim ---------------------------
CalibDataSvc *	m_calibDataSvc
	dedx sim ---------------------------
NTuple::Item< double >	m_charge
NTuple::Item< double >	m_charge
NTuple::Item< double >	m_costheta
NTuple::Item< double >	m_costheta
NTuple::Item< double >	m_de
NTuple::Item< double >	m_de
NTuple::Item< double >	m_dedx
NTuple::Item< double >	m_dedx
TH1F *	m_dedx_hists
TH1F *	m_dedx_hists
NTuple::Item< double >	m_fitval
NTuple::Item< double >	m_fitval
G4Svc *	m_G4Svc
G4Svc *	m_G4Svc
G4int	m_numBg
G4int	m_numDedxHists
IDedxCurSvc *	m_pDedxCurSvc
IDedxCurSvc *	m_pDedxCurSvc
NTuple::Item< double >	m_random
NTuple::Item< double >	m_random
NTuple::Tuple *	m_tupleMdc
NTuple::Tuple *	m_tupleMdc
std::vector< G4double >	mdc_hit_length
std::vector< G4double >	mdc_hit_length
BesMdcCalTransfer *	mdcCalPointer
BesMdcCalTransfer *	mdcCalPointer
MdcGeomSvc *	mdcGeomSvc
MdcGeomSvc *	mdcGeomSvc
BesMdcGeoParameter *	mdcGeoPointer
BesMdcGeoParameter *	mdcGeoPointer
BesMdcHitsCollection *	truthCollection
BesMdcHitsCollection *	truthCollection
G4int	truthPointer [43][288]

BesMdcSD Class Reference

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation