00001
00002
00003
00004
00005 #include "GaudiKernel/Bootstrap.h"
00006
00007 #include "BesMucNoise.hh"
00008 #include "G4Trap.hh"
00009 #include "G4VSolid.hh"
00010 #include <iostream>
00011 #include <fstream>
00012 #include <sstream>
00013 #include "ReadBoostRoot.hh"
00014 #include "Randomize.hh"
00015 #include "math.h"
00016 #include "strstream"
00017 #include "G4Box.hh"
00018 #include "stdlib.h"
00019
00020
00021 using namespace std;
00022
00023
00024 const int BesMucNoise::m_kSegment[m_kPart] = {4, 8, 4};
00025 const int BesMucNoise::m_kAbsorber[m_kPart] = {9, 9, 9};
00026 const int BesMucNoise::m_kGap[m_kPart] = {8, 9, 8};
00027 const int BesMucNoise::m_kPanel[m_kPart] = {4, 3, 4};
00028
00029
00030 BesMucNoise * BesMucNoise::fPointer=0;
00031 BesMucNoise * BesMucNoise::Instance(void){
00032 if(!fPointer)fPointer = new BesMucNoise();
00033 return fPointer;
00034 }
00035
00036
00037 BesMucNoise::BesMucNoise()
00038 {
00039 if(fPointer)
00040 {G4Exception("BesMucNoise constructed twice.");}
00041 fPointer=this;
00042
00043 }
00044
00045 BesMucNoise::~BesMucNoise()
00046 {
00047 if( m_ptrIdTr != NULL ) delete m_ptrIdTr;
00048
00049 }
00050
00051 void BesMucNoise::Initialize(G4String filename,G4LogicalVolume* logicalMuc,G4String temp)
00052 {
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 }
00083
00084 void BesMucNoise::Initialize(G4String filename,G4LogicalVolume* logicalMuc)
00085 {
00086 for(G4int part=0;part<3;part++){
00087 for(G4int seg=0;seg<8;seg++){
00088 for(G4int gap=0;gap<9;gap++){
00089 m_noise[part][seg][gap]=1;
00090 }
00091 }
00092 }
00093
00094 G4cout<<"filename: "<<filename<<G4endl;
00095 std::ifstream fin(filename);
00096 if(!fin){
00097 G4cout<<"error opening muc_noise data"<<G4endl;
00098 }
00099 char buffer[200];
00100 fin.getline(buffer,200,'\n');
00101 std::istringstream stringBuf(buffer);
00102
00103
00104 G4int tot_NoDaughter = logicalMuc->GetNoDaughters();
00105
00106 for(G4int i=0; i<tot_NoDaughter;i++){
00107 G4LogicalVolume* i_LogicalGap = logicalMuc->GetDaughter(i)->GetLogicalVolume();
00108 G4String i_GapName = i_LogicalGap->GetName();
00109
00110 if(i_GapName.find("G")==8){
00111 G4LogicalVolume* i_LogicalBox = i_LogicalGap->GetDaughter(0)->GetLogicalVolume();
00112 G4LogicalVolume* i_LogicalStripPlane = i_LogicalBox->GetDaughter(0)->GetLogicalVolume();
00113
00114 G4String strPart = i_GapName.substr(5,1);
00115 G4String strSeg = i_GapName.substr(7,1);
00116 G4String strGap = i_GapName.substr(9,1);
00117
00118 std::istrstream partBuf(strPart.c_str(), strlen(strPart.c_str()));
00119 std::istrstream segBuf(strSeg.c_str(), strlen(strSeg.c_str()));
00120 std::istrstream gapBuf(strGap.c_str(), strlen(strGap.c_str()));
00121
00122 G4int part,seg,gap;
00123
00124 partBuf >> part;
00125 segBuf >> seg;
00126 gapBuf >> gap;
00127
00128 G4int tot_NoStrip = i_LogicalStripPlane->GetNoDaughters();
00129 area[part][seg][gap][0]=tot_NoStrip;
00130 G4float tot_Area=0;
00131
00132 G4LogicalVolume* i_LogicalStrip1 = i_LogicalStripPlane->GetDaughter(1)->GetLogicalVolume();
00133 G4LogicalVolume* i_LogicalStrip2 = i_LogicalStripPlane->GetDaughter(2)->GetLogicalVolume();
00134 G4Box *temp1; G4Box *temp2;
00135
00136 temp1=(G4Box *)i_LogicalStrip1->GetSolid();temp2=(G4Box *)i_LogicalStrip2->GetSolid();
00137 G4float Width1 =temp1->GetXHalfLength()*2;G4float Width2 =temp2->GetXHalfLength()*2;
00138 G4float pos1 =i_LogicalStripPlane->GetDaughter(1)->GetObjectTranslation().x();
00139 G4float pos2 =i_LogicalStripPlane->GetDaughter(2)->GetObjectTranslation().x();
00140 if ( (part == 1 && gap%2 != 0) || (part != 1 && gap%2 == 0) ) {
00141 Width1=temp1->GetYHalfLength()*2; Width2 =temp2->GetYHalfLength()*2;
00142 pos1 =i_LogicalStripPlane->GetDaughter(1)->GetObjectTranslation().y();
00143 pos2 =i_LogicalStripPlane->GetDaughter(2)->GetObjectTranslation().y();
00144 }
00145 G4float width_between_strip=pos2-pos1-Width1/2-Width2/2;
00146
00147
00148 for(G4int j=0;j<tot_NoStrip;j++){
00149 G4LogicalVolume* i_LogicalStrip = i_LogicalStripPlane->GetDaughter(j)->GetLogicalVolume();
00150 G4Box *temp;
00151 temp=(G4Box *)i_LogicalStrip->GetSolid();
00152 G4float Width =temp->GetXHalfLength()*2;
00153 G4float Length=temp->GetYHalfLength()*2;
00154 if ( (part == 1 && gap%2 != 0) || (part != 1 && gap%2 == 0) ) {
00155 Width =temp->GetYHalfLength()*2;
00156 Length=temp->GetXHalfLength()*2;
00157 }
00158
00159 if(j==0||j==(tot_NoStrip-1)) Width=Width+width_between_strip/2;
00160 else Width=Width+width_between_strip;
00161 G4float Strip_Area=fabs(Width*Length);
00162 tot_Area=tot_Area+Strip_Area;
00163 area[part][seg][gap][j+1]=tot_Area;
00164 strip_area[part][seg][gap][j] = Strip_Area;
00165
00166 }
00167
00168
00169 box_area[part][seg][gap] = tot_Area;
00170
00171 for(G4int k=1;k<tot_NoStrip+1;k++){
00172 area[part][seg][gap][k]=area[part][seg][gap][k]/tot_Area;
00173
00174 }
00175
00176 }
00177 }
00178
00179
00180
00181
00182
00183
00184
00185 m_ptrIdTr = new MucIdTransform();
00186 CheckCalibSvc();
00187 InitProb();
00188 }
00189
00190 void BesMucNoise::CheckCalibSvc()
00191 {
00192 ISvcLocator* svcLocator = Gaudi::svcLocator();
00193 StatusCode sc = svcLocator->service("MucCalibConstSvc", m_ptrCalibSvc, true);
00194
00195 if( sc != StatusCode::SUCCESS){
00196 G4cout<< "Can not use MucCalibConstSvc!" << G4endl;
00197 }
00198 }
00199
00200 void BesMucNoise::InitProb()
00201 {
00202 m_HitMean = 0.08;
00203 int NHIT = 20;
00204
00205 for( int i=0; i<NHIT; i++) {
00206 m_Prob[i][0] = m_Prob[i][1] = 0.;
00207 }
00208
00209 double sum = 0.;
00210 double EU = TMath::Power(TMath::Exp(1.0), -m_HitMean);
00211 for( int i=0; i<NHIT; i++)
00212 {
00213 m_Prob[i][0] = EU*TMath::Power(m_HitMean, i)/TMath::Factorial(i);
00214 sum += m_Prob[i][0];
00215 m_Prob[i][1] = sum;
00216
00217 }
00218 }
00219
00220 G4int BesMucNoise::AddNoise(int model, BesMucHitsCollection* aMucHitCollection,BesMucHitsCollection* aMucHitList)
00221 {
00222 G4int noiseNum = 0;
00223 if( model == 1 ) noiseNum = NoiseByCnt(aMucHitCollection, aMucHitList);
00224 else noiseNum = NoiseByNosRatio(aMucHitCollection, aMucHitList);
00225
00226 return noiseNum;
00227
00228 }
00229
00230 G4int BesMucNoise::NoiseByCnt(BesMucHitsCollection* MucHitCollection,BesMucHitsCollection* MucHitList)
00231 {
00232 G4int noiseNum = 0;
00233 m_noiseLevel = m_ptrCalibSvc->getLevel();
00234
00235 for(int i = 0; i < m_kPart; i++) {
00236 for(int j = 0; j < m_kSegment[i]; j++) {
00237 for(int k = 0; k < m_kGap[i]; k++) {
00238 if(m_noiseLevel == 2)
00239 {
00240 G4int hitNum = NoiseSampling( m_noiseLevel, i, j, k, 0 );
00241 for(G4int ii=0;ii<hitNum;ii++)
00242 {
00243 G4int strip = GetStripNo(i,j,k);
00244 BesMucHit *noiseHit=new BesMucHit(i, j, k, strip, -1, -1);
00245 BesMucHit *noiseHit2=new BesMucHit(i, j, k, strip, -1, -1);
00246 bool noiseHitExist = false;
00247 noiseHitExist = IsExist(noiseHit, MucHitList);
00248 MucHitCollection->insert(noiseHit);
00249 if(!noiseHitExist) {
00250 MucHitList->insert(noiseHit2);
00251 noiseNum += 1;
00252
00253 }
00254 else delete noiseHit2;
00255 }
00256 }
00257
00258 if(m_noiseLevel == 3)
00259 {
00260 int stripNum = m_ptrIdTr->GetStripMax(i, j, k);
00261 for(int strip = 0; strip < stripNum; strip++)
00262 {
00263 G4int hitNum = NoiseSampling( m_noiseLevel, i, j, k, strip );
00264 if(hitNum > 0){
00265 BesMucHit *noiseHit=new BesMucHit(i, j, k, strip, -1, -1);
00266 BesMucHit *noiseHit2=new BesMucHit(i, j, k, strip, -1, -1);
00267 bool noiseHitExist = false;
00268 noiseHitExist = IsExist(noiseHit, MucHitList);
00269 MucHitCollection->insert(noiseHit);
00270 if(!noiseHitExist) {
00271 MucHitList->insert(noiseHit2);
00272 noiseNum += 1;
00273
00274 }
00275 else delete noiseHit2;
00276 }
00277 }
00278 }
00279 }
00280 }
00281 }
00282
00283 return noiseNum;
00284 }
00285
00286 G4int BesMucNoise::NoiseByNosRatio(BesMucHitsCollection* MucHitCollection,BesMucHitsCollection* MucHitList)
00287 {
00288 G4int noiseNum = 0;
00289 G4float random = 0;
00290
00291
00292 random = G4UniformRand();
00293 for(int i=0; i<20; i++) {
00294 if(random<m_Prob[i][1]) {noiseNum = i; break;}
00295 }
00296
00297
00298 int prt, seg, lay, str, tmp_strip;
00299 G4float nosRatio = 0.;
00300
00301 for(int i=0; i<noiseNum; i++)
00302 {
00303 do{
00304 random = G4UniformRand();
00305 tmp_strip = TMath::Nint(random*STRIP_MAX);
00306 m_ptrIdTr->SetStripPos(tmp_strip, &prt, &seg, &lay, &str);
00307 nosRatio = m_ptrCalibSvc->getNosRatio(prt, seg, lay, str);
00308 random = G4UniformRand();
00309 if( random<nosRatio ) break;
00310 }while( 1 );
00311
00312
00313
00314 BesMucHit *noiseHit = new BesMucHit(prt, seg, lay, str, -1, -1);
00315 BesMucHit *noiseHit2 = new BesMucHit(prt, seg, lay, str, -1, -1);
00316
00317 bool noiseHitExist = false;
00318 noiseHitExist = IsExist(noiseHit, MucHitList);
00319 MucHitCollection->insert(noiseHit);
00320 if(!noiseHitExist) MucHitList->insert(noiseHit2);
00321 else delete noiseHit2;
00322 }
00323
00324 return noiseNum;
00325 }
00326
00327 bool BesMucNoise::IsExist(BesMucHit* aNoiseHit, BesMucHitsCollection* aMucHitList)
00328 {
00329 bool isExist = false;
00330 G4int n_hit = aMucHitList->entries();
00331 for(G4int iNoise=0;iNoise<n_hit;iNoise++)
00332 {
00333 if ( aNoiseHit->GetTrackIndex()%1000 == (*aMucHitList)[iNoise]->GetTrackIndex()%1000 &&
00334 aNoiseHit->GetPart() == (*aMucHitList)[iNoise]->GetPart() &&
00335 aNoiseHit->GetSeg() == (*aMucHitList)[iNoise]->GetSeg() &&
00336 aNoiseHit->GetGap() == (*aMucHitList)[iNoise]->GetGap() &&
00337 aNoiseHit->GetStrip() == (*aMucHitList)[iNoise]->GetStrip() )
00338 {
00339 isExist = true;
00340 break;
00341 }
00342 }
00343
00344 return isExist;
00345 }
00346
00347 G4int BesMucNoise::NoiseSampling(int level, int prt, int seg, int lay, int strip)
00348 {
00349 G4int hitNum = 0;
00350 G4double lambda;
00351
00352 G4double t=800;
00353 G4double e=2.71828182845904590;
00354
00355 if( level == 2 )
00356 lambda = m_ptrCalibSvc->getBoxCnt(prt,seg,lay) * box_area[prt][seg][lay] * 1E-2 * 1E-9;
00357 else if( level == 3)
00358 lambda = m_ptrCalibSvc->getStripCnt(prt,seg,lay,strip) * strip_area[prt][seg][lay][strip] * 1E-2 * 1E-9;
00359 else
00360 lambda = 0.;
00361
00362
00363 G4float random=G4UniformRand();
00364 G4double prob = 0;
00365 do{
00366 prob=prob+pow(e,-lambda*t)*pow(lambda*t,hitNum)/Factorial(hitNum);
00367 if(random<prob) break;
00368 hitNum++;
00369 }while(1);
00370
00371
00372 return hitNum;
00373 }
00374
00375 G4float BesMucNoise::Factorial(G4int i)
00376 {
00377 G4float fact=1;
00378 if(i==0||i==1)return 1;
00379 else{
00380 for(G4int ii=2;ii<=i;ii++){
00381 fact=fact*ii;}
00382 return fact;
00383 }
00384
00385 }
00386
00387 G4int BesMucNoise::GetStripNo(G4int part,G4int seg,G4int gap)
00388 { G4int stripno;
00389
00390 G4float random=(rand()%100000)/100000.0;
00391 if(part==1){
00392 G4float width=area[part][seg][gap][3]-area[part][seg][gap][2];
00393 stripno=G4int((random-area[part][seg][gap][1])/width)+2;
00394 if(stripno<1)stripno=1;
00395 if(stripno>111)stripno=111;
00396
00397 G4int step = IsNearestStrip(stripno,part,seg,gap,random);
00398 while(step!=0) {
00399
00400 stripno += step;
00401 step = IsNearestStrip(stripno,part,seg,gap,random);
00402 }
00403 stripno--;
00404 return stripno;
00405 }
00406 else{
00407 G4int max,min,mid,pass;
00408 min=1;
00409 max=area[part][seg][gap][0];
00410 mid=G4int((min+max)/2);
00411
00412 do{
00413
00414 pass=0;
00415 if(random>area[part][seg][gap][mid]){
00416 min=mid;
00417 mid=G4int((min+max)/2);
00418 }
00419 else if(random<area[part][seg][gap][mid-1]){
00420 if(random<area[part][seg][gap][1]){
00421 pass=1; mid=1;max=1;
00422 }else{
00423 max=mid-1;
00424 mid=G4int((min+max)/2);
00425 }
00426 }else{pass=1;}
00427
00428 if(min==mid)mid=max;
00429
00430 }while(pass==0&&(max>mid&&mid>min));
00431
00432 return mid-1;
00433
00434 }
00435
00436 }
00437
00438 G4int BesMucNoise::IsNearestStrip(G4int stripno,G4int part,G4int seg,G4int gap,G4float random)
00439 {
00440 if(stripno==1){
00441 return 0;
00442 }else{
00443 if(area[part][seg][gap][stripno]!=0){
00444
00445 if(random<=area[part][seg][gap][stripno]&&random>area[part][seg][gap][stripno-1])return 0;
00446 if(random<=area[part][seg][gap][stripno-1]) return -1;
00447 if(random>area[part][seg][gap][stripno]) return 1;
00448 }else{
00449 return -1;
00450 }
00451 }
00452 }