00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00034
00041 #include <Stiostream.h>
00042 #include "StBTofSimMaker.h"
00043
00044
00045 #include <math.h>
00046 #include "TRandom.h"
00047 #include "SystemOfUnits.h"
00048 #include "phys_constants.h"
00049 #include "StThreeVectorD.hh"
00050 #include "Random.h"
00051 #include "RanluxEngine.h"
00052 #include "RandGauss.h"
00053 #include "TH1.h"
00054 #include "TH2.h"
00055 #include "TNtuple.h"
00056 #include "TFile.h"
00057
00058
00059 #include "tables/St_g2t_ctf_hit_Table.h"
00060 #include "tables/St_g2t_vpd_hit_Table.h"
00061 #include "tables/St_g2t_track_Table.h"
00062 #include "StMcTrack.hh"
00063
00064 #include "StBTofUtil/StBTofDaqMap.h"
00065 #include "StTofUtil/tofPathLength.hh"
00066 #include "StTofUtil/StTofSimParam.h"
00067 #include "StBTofUtil/StBTofGeometry.h"
00068 #include "StEventTypes.h"
00069 #include "StEvent/StBTofCollection.h"
00070
00071
00072 static RanluxEngine engine;
00073 static RandGauss ranGauss(engine);
00074
00075 ClassImp(StBTofSimMaker)
00076
00077
00078 StBTofSimMaker::StBTofSimMaker(const char *name):StMaker(name)
00079 {
00080
00081 mBookHisto=kFALSE;
00082 mSlow=kTRUE;
00083 mCellXtalk=kTRUE;
00084 mWriteStEvent=kTRUE;
00085 mDaqMap=0;
00086 Reset();
00087
00088 }
00089
00090
00091 StBTofSimMaker::~StBTofSimMaker()
00092 {
00093 delete mSimDb;
00094 delete mDaqMap;
00095
00096 }
00097
00098
00099
00100 Int_t StBTofSimMaker::Init()
00101 {
00102 Reset();
00103 mSimDb = new StTofSimParam();
00104 if (Debug()) mSimDb->print();
00105 if(mBookHisto) bookHistograms();
00106
00107 return StMaker::Init();
00108 }
00109
00110
00111 void StBTofSimMaker::Reset()
00112 {
00113 mGeantData = 0;
00114 mEvent = 0;
00115 mMcEvent = 0;
00116
00117 if (mWriteStEvent) delete mBTofCollection;
00118 delete mMcBTofHitCollection;
00119 mSimDb = 0;
00120
00121 if(mDaqMap){delete mDaqMap; mDaqMap = 0;}
00122
00123 ResetFlags();
00124 }
00125
00126
00127 Int_t StBTofSimMaker::ResetFlags()
00128 {
00130 memset(mTofHitFlag, 0, sizeof(mTofHitFlag));
00131 memset(mVpdHitFlag, 0, sizeof(mVpdHitFlag));
00132 return kStOk;
00133 }
00134
00135
00136 Int_t StBTofSimMaker::InitRun(Int_t runnumber)
00137 {
00138 LOG_INFO << "StBTofSimMaker::InitRun -- initializing BTOF DAQ map --" << endm;
00139
00141 mDaqMap = new StBTofDaqMap();
00142 mDaqMap->Init(this);
00143
00144 return kStOK;
00145 }
00146
00147
00148 Int_t StBTofSimMaker::FinishRun(Int_t runnumber)
00149 {
00150 LOG_INFO << "StBTofSimMaker::FinishRun -- cleaning up BTOF DAQ map --" << endm;
00151 if (mDaqMap){delete mDaqMap; mDaqMap = 0;}
00152 return kStOk;
00153 }
00154
00155
00156 Int_t StBTofSimMaker::Finish()
00157 {
00158 if(mBookHisto){
00159 LOG_INFO << "StBTofSimMaker::Finish writing tofsim.root ...";
00160 TFile theFile(mHistFile.c_str(),"RECREATE","tofsim");
00161 theFile.cd();
00162 writeHistograms();
00163 theFile.Close();
00164 }
00165
00166
00167 return kStOK;
00168 }
00169
00170
00171
00172 Int_t StBTofSimMaker::Make()
00173 {
00174 LOG_INFO << "StBTofSimMaker Make() starts" << endm;
00175
00176 ResetFlags();
00177
00178 mMcBTofHitCollection = new StMcBTofHitCollection();
00179
00180
00181 mGeantData = GetInputDS("geant");
00182 if(!mGeantData) {
00183 mGeantData = GetInputDS("geantBranch");
00184 }
00185 if(!mGeantData) {
00186 LOG_WARN << " No GEANT data loaded. Exit! " << endm;
00187 return kStWarn;
00188 }
00189 LOG_INFO << " Found GEANT data -- loading VPD/TOF hits... " << endm;
00190
00191
00192 St_g2t_vpd_hit* g2t_vpd_hits = 0;
00193 g2t_vpd_hits = dynamic_cast<St_g2t_vpd_hit*>(mGeantData->Find("g2t_vpd_hit"));
00194 if(!g2t_vpd_hits){
00195 LOG_WARN << " No VPD hits in GEANT" << endm; }
00196 else {
00197 Int_t nvpdhits = g2t_vpd_hits->GetNRows();
00198 LOG_DEBUG << " Found VPD hits: " << nvpdhits << endm;
00199 g2t_vpd_hit_st* vpd_hit = g2t_vpd_hits->begin();
00200 for(Int_t i=0; i<nvpdhits; i++, vpd_hit++) {
00201 VpdResponse(vpd_hit);
00202 }
00203 }
00204 LOG_INFO << " McBTofHit Size (VPD) = " << mMcBTofHitCollection->hits().size() << endm;
00205
00206
00207
00208 St_g2t_ctf_hit* g2t_tfr_hits = 0;
00209 g2t_tfr_hits = dynamic_cast<St_g2t_ctf_hit*> (mGeantData->Find("g2t_tfr_hit"));
00210 if(!g2t_tfr_hits) {
00211 LOG_WARN << " No TOF hits in GEANT" << endm; }
00212 else {
00213 Int_t nhits = g2t_tfr_hits->GetNRows();
00214 LOG_DEBUG << " Found TOF hits: " << nhits << endm;
00215 g2t_ctf_hit_st* tofHitsFromGeant = g2t_tfr_hits->begin();
00216
00217 if(mSlow) {
00218
00219 TrackVec tofResponseVec;
00220 tofResponseVec.clear();
00221 for (Int_t i=0; i<nhits; i++, tofHitsFromGeant++) {
00222
00223 CellResponse(tofHitsFromGeant,tofResponseVec);
00224 }
00225
00226 CellTimePassTh(tofResponseVec);
00227 }
00228 else if(!mSlow) {
00229 for (Int_t i=0; i<nhits; i++, tofHitsFromGeant++) FastCellResponse(tofHitsFromGeant);
00230 }
00231 else
00232 LOG_WARN << " No TOF simulator specified " << endm;
00233 }
00234 LOG_INFO << " McBTofHit Size (TOF) = " << mMcBTofHitCollection->hits().size() << endm;
00235
00236
00237 fillEvent();
00238
00239
00240 LOG_INFO << " Leaving StBTofSimMaker ... done " << endm;
00241 return kStOK;
00242 }
00243
00244
00247 Int_t StBTofSimMaker::VpdResponse(g2t_vpd_hit_st* vpd_hit)
00248 {
00249 if(!vpd_hit) {
00250 LOG_WARN << " No VPD hit!" << endm;
00251 return kStWarn;
00252 }
00253
00254 Int_t vId = vpd_hit->volume_id;
00255 Int_t itray = vId/1000 + 120;
00256 Int_t itube = vId%100;
00257
00258 Double_t tof = (vpd_hit->tof + ranGauss.shoot()*mSimDb->timeres_vpd())*1000./nanosecond;
00259 Double_t t0 = vpd_hit->tof*1000./nanosecond;
00260 Double_t de = vpd_hit->de;
00261 Double_t pathL = -9999.;
00262 Double_t q = 0.;
00263
00264
00265 StMcBTofHit *mcHit = new StMcBTofHit(itray,0,itube,de,pathL,t0,tof,q);
00266 storeMcBTofHit(mcHit);
00267 mVpdHitFlag[(itray-121)*mNVPD+(itube-1)] = 1;
00268
00269 return kStOk;
00270 }
00271
00272
00274 Int_t StBTofSimMaker::CellResponse(g2t_ctf_hit_st* tofHitsFromGeant,
00275 TrackVec& tofResponseVec)
00276 {
00285
00286
00287 if(tofHitsFromGeant->s_track<=0.0 || tofHitsFromGeant->de <=0.0) {
00288 LOG_WARN << " No deposited energy in this TOF hit!" << endm;
00289 return kStWarn;
00290 }
00291
00292 IntVec cellId = CalcCellId(tofHitsFromGeant->volume_id, tofHitsFromGeant->x[1]);
00293 Int_t icell, imodule, itray;
00294 itray = cellId[0];
00295 imodule = cellId[1];
00296 icell = cellId[2];
00297 if (itray==-1 || imodule==-1 || icell==-1) {
00298 LOG_WARN << " Not hit the sensitive MRPC volume!" << endm;
00299 return kStWarn;
00300 }
00301 if(mBookHisto) {
00302 mDeGeant->Fill(tofHitsFromGeant->de / keV);
00303 mTofGeant->Fill(tofHitsFromGeant->tof / nanosecond);
00304 }
00305
00306
00307 St_g2t_track *g2t_track = static_cast<St_g2t_track *>(mGeantData->Find("g2t_track"));
00308 if (!g2t_track) {
00309 LOG_WARN << " No G2T track table!" << endm;
00310 return kStWarn;
00311 };
00312 g2t_track_st *tof_track = g2t_track->GetTable();
00313 Int_t no_tracks= g2t_track->GetNRows();
00314
00315 Double_t beta;
00316 Int_t trackId = -1;
00317 for(Int_t j=0;j<no_tracks;j++){
00318 if(tofHitsFromGeant->track_p==tof_track[j].id){
00319 trackId = j;
00320 beta = tof_track[j].ptot/tof_track[j].e;
00321 break;
00322 }
00323 }
00324
00325 Double_t qtot=-1;
00326 Double_t tof=-1;
00327 Double_t t0 = tofHitsFromGeant->tof;
00328 Float_t wt=1.0;
00329
00330 Double_t clusterDensity = mSimDb->nclus(beta);
00331 Double_t gapLength = mSimDb->dg();
00332 Double_t alpha = mSimDb->alpha();
00333 Double_t ka = mSimDb->ka();
00334 Double_t kaa = ka/(alpha*gapLength);
00335
00336 const Int_t maxClusters=mSimDb->nmaxclus();
00337 const Int_t nTimeBins = mSimDb->ndt();
00338 Double_t driftVelocity[maxClusters],nElectrons[maxClusters],startPositionOfCluster[maxClusters],sa[maxClusters];
00339 Double_t s[maxClusters][nTimeBins];
00340
00341 Double_t chargeDepositedPerTimeBin[nTimeBins];
00342 for(Int_t j=0;j<nTimeBins;j++) {chargeDepositedPerTimeBin[j] = 0.0;}
00343
00344 Int_t nElectronClusters=-1;
00345 while(nElectronClusters<1) {nElectronClusters=gRandom->Poisson(gapLength*clusterDensity);}
00346 if(nElectronClusters>maxClusters) nElectronClusters = maxClusters;
00347
00348 for(Int_t m=0;m<nElectronClusters;m++) {
00349 driftVelocity[m] = mSimDb->vd_mean()*(0.90+0.20*gRandom->Rndm(1));
00350
00351 Int_t nElectrons_temp=-1;
00352 while(nElectrons_temp<1){nElectrons_temp = gRandom->Poisson(mSimDb->nmeane());}
00353 nElectrons[m] = Double_t(nElectrons_temp);
00354
00355 startPositionOfCluster[m]=gapLength+1;
00356 while(startPositionOfCluster[m]>gapLength) { startPositionOfCluster[m] = gRandom->Exp(1.0/clusterDensity); }
00357 }
00358
00359 Double_t ytmp=0.0;
00360 for(Int_t m=0;m<nElectronClusters;m++) {
00361 sa[m] = (exp(alpha*(gapLength-startPositionOfCluster[m]))-1)*nElectrons[m]*GammaRandom();
00362 if (sa[m]>mSimDb->nmaxe()) sa[m] = mSimDb->nmaxe();
00363 ytmp += kaa*sa[m];
00364 }
00365 qtot = ytmp*1.e+12*e_SI;
00366
00367 t0 = tofHitsFromGeant->tof*1000/nanosecond;
00368 tof=tofHitsFromGeant->tof*1000/nanosecond;
00369
00370
00371 for(Int_t j=0;j<nTimeBins;j++) {
00372 Double_t ts = t0+ mSimDb->dt()*double(j);
00373 Double_t ytmp1 = 0.;
00374 for(Int_t m=0;m<nElectronClusters;m++) {
00375 Double_t tx = (startPositionOfCluster[m])/(C_C_LIGHT*1.e-3*nanosecond/millimeter);
00376 Double_t t_drift = (gapLength-startPositionOfCluster[m])/driftVelocity[m];
00377 if( ts>=t0 + tx && ts<=t0+ tx+t_drift) {
00378 s[m][j]=(exp(alpha*driftVelocity[m]*(ts-t0-tx))-1)*nElectrons[m]*GammaRandom();
00379 if(s[m][j]>mSimDb->nmaxe()) { s[m][j] = mSimDb->nmaxe(); }
00380 } else {
00381 s[m][j]=0.0;
00382 }
00383 ytmp1 += kaa*s[m][j];
00384 }
00385 chargeDepositedPerTimeBin[j] = ytmp1*1.e+12*e_SI;
00386 }
00387
00388
00389 Int_t icellx = -1;
00390 wt = 1.0;
00391 StThreeVectorF local(tofHitsFromGeant->x[0], tofHitsFromGeant->x[1], tofHitsFromGeant->x[2]);
00392 if(mCellXtalk){ CellXtalk(icell, local.y(), wt, icellx); }
00393 TrackHit trackhit;
00394 trackhit.tray = itray;
00395 trackhit.module = imodule;
00396 trackhit.cell = icell;
00397 trackhit.trkId = trackId;
00398 trackhit.dE = tofHitsFromGeant->de * wt;
00399 trackhit.dQ = qtot * wt;
00400 for(Int_t j=0;j<nTimeBins;j++) {
00401 trackhit.dQdt[j] = chargeDepositedPerTimeBin[j] * wt;
00402 }
00403 trackhit.tof = tof;
00404 trackhit.s_track = tofHitsFromGeant->s_track;
00405 trackhit.position = local;
00406 trackhit.t0 = t0/1000.;
00407 tofResponseVec.push_back(trackhit);
00408 mTofHitFlag[itray-1][(imodule-1)*mNCell+(icell-1)] = 1;
00409
00410 if(icellx>0 && icellx<=mNCell){
00411 TrackHit trackhitx=trackhit;
00412 trackhitx.cell = icellx;
00413 trackhitx.dE = tofHitsFromGeant->de * (1.-wt);
00414 trackhitx.dQ = qtot * (1.-wt);
00415 for(Int_t j=0;j<nTimeBins;j++) {
00416 trackhitx.dQdt[j] = chargeDepositedPerTimeBin[j] * (1.-wt);
00417 }
00418 tofResponseVec.push_back(trackhitx);
00419 mTofHitFlag[itray-1][(imodule-1)*mNCell+(icellx-1)] = 1;
00420 }
00421
00422 return kStOk;
00423 }
00424
00425
00426 Int_t StBTofSimMaker::CellTimePassTh(TrackVec& tofResponseVec)
00427
00428
00429
00430
00431 {
00432 TrackVec trackSumVec;
00433 trackSumVec.clear();
00434
00435 Int_t eraseId[500000];
00436 Int_t nhits = tofResponseVec.size();
00437 Int_t nTimeBins = mSimDb->ndt();
00438 for(Int_t i=0;i<nhits;i++) eraseId[i] = 0;
00439
00440
00441 for(Int_t i=0;i<nhits;i++){
00442 if(eraseId[i]) continue;
00443 TrackHit sumhit=tofResponseVec[i];
00444
00445 for(Int_t j=i+1;j<nhits;j++) {
00446 if(eraseId[j]) continue;
00447
00448 if(tofResponseVec[j].tray != sumhit.tray ||
00449 tofResponseVec[j].module != sumhit.module ||
00450 tofResponseVec[j].cell != sumhit.cell) continue;
00451
00452 if(tofResponseVec[j].tof < sumhit.tof) {
00453 sumhit.tof = tofResponseVec[j].tof;
00454 sumhit.s_track = tofResponseVec[j].s_track;
00455 sumhit.position = tofResponseVec[j].position;
00456 if(tofResponseVec[j].trkId != sumhit.trkId) {
00457 LOG_WARN << " Two tracks match to one cell." << endm;
00458 sumhit.trkId = tofResponseVec[j].trkId;
00459 }
00460 }
00461 sumhit.dE += tofResponseVec[j].dE;
00462 sumhit.dQ += tofResponseVec[j].dQ;
00463
00464 Double_t dQdt[nTimeBins]; for(Int_t aa=0;aa<nTimeBins;aa++){dQdt[aa]=sumhit.dQdt[aa];}
00465
00466 if(sumhit.t0 == tofResponseVec[j].t0) {
00467 for(Int_t m=0;m<nTimeBins;m++) { sumhit.dQdt[m] += tofResponseVec[j].dQdt[m];}
00468 }
00469 else if(sumhit.t0 > tofResponseVec[j].t0) {
00470 Int_t nbinoffset = (int)((sumhit.t0 - tofResponseVec[j].t0) / 25.);
00471 for(Int_t m=0;m<nTimeBins;m++) {
00472 if(m<nbinoffset) {
00473 sumhit.dQdt[m] = tofResponseVec[j].dQdt[m];
00474 } else{
00475 sumhit.dQdt[m] = tofResponseVec[j].dQdt[m] + dQdt[m-nbinoffset];
00476 }
00477 }
00478 sumhit.t0 = tofResponseVec[j].t0;
00479 } else {
00480 Int_t nbinoffset = (int)((tofResponseVec[j].t0 - sumhit.t0) / 25.);
00481 for(Int_t m=0;m<nTimeBins;m++) {
00482 if(m<nbinoffset) {
00483
00484 } else{
00485 sumhit.dQdt[m] = tofResponseVec[j].dQdt[m-nbinoffset] + dQdt[m];
00486 }
00487 }
00488 }
00489
00490 eraseId[j] = 1;
00491 }
00492 trackSumVec.push_back(sumhit);
00493 }
00494
00495
00496
00498 St_g2t_track *g2t_track = static_cast<St_g2t_track *>(mGeantData->Find("g2t_track"));
00499 if (!g2t_track) {
00500 LOG_WARN << " No g2t track table !!! " << endm;
00501 return kStWarn;
00502 }
00503 g2t_track_st *tof_track = g2t_track->GetTable();
00504 for(size_t i=0;i<trackSumVec.size();i++) {
00505 Float_t tof = 0.;
00506 bool pass = kFALSE;
00507 for(Int_t m=0;m<nTimeBins;m++) {
00508 if(trackSumVec[i].dQdt[m]>(mSimDb->adc_thre()*0.001) && !pass) {
00509 tof = trackSumVec[i].tof;
00510 pass = kTRUE;
00511 break;
00512 }
00513 }
00514
00515
00516 Float_t deltaMRPC=ranGauss.shoot()*85.;
00517 tof+=deltaMRPC;
00518
00519
00520 if(pass) {
00521 StMcBTofHit *mcHit = new StMcBTofHit();
00522 StMcTrack *partnerTrk = new StMcTrack(&(tof_track[trackSumVec[i].trkId]));
00523 Int_t truthId=partnerTrk->key();
00524 mcHit->setTray(trackSumVec[i].tray);
00525 mcHit->setModule(trackSumVec[i].module);
00526 mcHit->setCell(trackSumVec[i].cell);
00527 mcHit->setdE(trackSumVec[i].dE);
00528 float pathLength=trackSumVec[i].s_track;
00529 mcHit->setPathLength(pathLength);
00530 mcHit->setTime(trackSumVec[i].tof);
00531 mcHit->setTof(tof);
00532 mcHit->setCharge(trackSumVec[i].dQ);
00533 mcHit->setPosition(trackSumVec[i].position);
00534 mcHit->setParentTrack(partnerTrk);
00535 mcHit->setParentTrackId(truthId);
00536 mMcBTofHitCollection->addHit(mcHit);
00537
00538
00539 if (mBookHisto){
00540 Float_t beta=pathLength/tof/3e-2;
00541 mBetaHist->Fill(beta);
00542 mPathLHist->Fill(pathLength);
00543 mTofHist->Fill(tof);
00544 double momentum=partnerTrk->momentum().mag();
00545 double mass=sqrt(beta*beta*momentum*momentum/(1.-beta*beta));
00546 if(beta!=1.0 && pathLength>150){ mRecMass->Fill(mass);}
00547 mTofResReco->Fill( (tof - trackSumVec[i].t0*1000.) );
00548 }
00549 }
00550 }
00551
00552 return kStOk;
00553 }
00554
00555
00556 Int_t StBTofSimMaker::fillEvent()
00557 {
00558 LOG_DEBUG << "Filling McEvent and Event"<<endm;
00559
00560
00561 if(mBookHisto) {
00562 for(Int_t i=0;i<mNTray;i++) {
00563 Int_t ncell = 0;
00564 for(Int_t j=0;j<mNTOF;j++) {
00565 if(mTofHitFlag[i][j]) {
00566 mCellGeant->Fill(j,i);
00567 ncell++;
00568 }
00569 }
00570 mNCellGeant->Fill(ncell,i);
00571 }
00572 Int_t ne = 0;
00573 Int_t nw = 0;
00574 for(Int_t i=0;i<2*mNVPD;i++) {
00575 if(mVpdHitFlag[i]) {
00576 mVpdGeant->Fill(i);
00577 if(i<mNVPD) nw++;
00578 else ne++;
00579 }
00580 }
00581 mNVpdGeant->Fill(nw,ne);
00582 }
00583
00585 mMcEvent = (StMcEvent*)GetInputDS("StMcEvent");
00586 if (!mMcEvent) {
00587 LOG_ERROR << "No StMcEvent! Bailing out ..." << endm;
00588 } else {
00589 mMcEvent->setBTofHitCollection(mMcBTofHitCollection);
00590 LOG_INFO << " ... StMcTofCollection stored in StMcEvent" << endm;
00591 }
00592
00594 if (mWriteStEvent){
00595 mBTofCollection= new StBTofCollection();
00596 mEvent = (StEvent*)GetInputDS("StEvent");
00597 if (!mEvent) {
00598 LOG_ERROR << "No StEvent! Bailing out ..." << endm;
00599 }
00600
00602 for(Int_t jj = 0; jj < (Int_t)mMcBTofHitCollection->hits().size(); jj++) {
00603 StMcBTofHit *aMcBTofHit = mMcBTofHitCollection->hits()[jj];
00604
00605 if(!aMcBTofHit) continue;
00606 Int_t trayid = aMcBTofHit->tray();
00607 Int_t moduleid = aMcBTofHit->module();
00608 Int_t cellid = aMcBTofHit->cell();
00609
00610 Int_t chn = 255;
00611 if(trayid>0&&trayid<=120) {
00612 chn = mDaqMap->Cell2TDIGChan(moduleid,cellid);
00613 } else if(trayid==121) {
00614 chn = mDaqMap->WestPMT2TDIGLeChan(cellid);
00615 } else if(trayid==122) {
00616 chn = mDaqMap->EastPMT2TDIGLeChan(cellid);
00617 }
00618
00619 Int_t index;
00620 if(trayid>0&&trayid<=120) {
00621 index = (trayid-1)*mNTOF + (moduleid-1)*mNCell + (cellid-1);
00622 } else if(trayid==121||trayid==122) {
00623 index = (trayid-1)*mNTOF + (cellid-1);
00624 }
00625
00626
00628 Float_t eff = 1.;
00629 if(trayid>0&&trayid<=120) eff = mSimDb->eff_tof(trayid, moduleid, cellid);
00630 else if(trayid==121||trayid==122) eff = mSimDb->eff_vpd(trayid, cellid);
00631 if (gRandom->Uniform(1.0) > eff){cout<<"REMOVED"<<endl; continue; }
00632
00633
00634
00635 StBTofHit aBTofHit;
00636 aBTofHit.Clear();
00637
00638 Float_t mcTof=aMcBTofHit->tof()/1000.;
00639
00640 aBTofHit.setHardwarePosition(kBTofId);
00641 aBTofHit.setTray((Int_t)aMcBTofHit->tray());
00642 aBTofHit.setModule((unsigned char)aMcBTofHit->module());
00643 aBTofHit.setCell((Int_t)aMcBTofHit->cell());
00644 aBTofHit.setLeadingEdgeTime((Double_t)mcTof);
00645 aBTofHit.setTrailingEdgeTime((Double_t)mcTof);
00646 aBTofHit.setAssociatedTrack(NULL);
00647 aBTofHit.setIdTruth(aMcBTofHit->parentTrackId(), 0);
00648 mBTofCollection->addHit(new StBTofHit(aBTofHit));
00649
00650
00651 StBTofRawHit aBTofRawHit;
00652 aBTofRawHit.Clear();
00653 aBTofRawHit.setTray((Int_t)aMcBTofHit->tray());
00654 aBTofRawHit.setChannel(6*(aMcBTofHit->module() - 1) + (Int_t)aMcBTofHit->cell());
00655 aBTofRawHit.setFlag(1);
00656 mBTofCollection->addRawHit(new StBTofRawHit(aBTofRawHit));
00657
00658 }
00659
00660
00661 StBTofHeader aHead;
00662 mBTofCollection->setHeader(new StBTofHeader(aHead));
00663
00664
00665 mEvent->setBTofCollection(mBTofCollection);
00666
00667 LOG_INFO << "... StBTofCollection Stored in StEvent! " << endm;
00668 }
00669
00671 if(Debug()) {
00672 LOG_DEBUG << " ==== Test McBTofHitCollection ==== " << endm;
00673 StSPtrVecMcBTofHit& mcBTofHits = mMcEvent->btofHitCollection()->hits();
00674 Int_t nCell[mNTray];
00675 for(Int_t i=0;i<mNTray;i++) nCell[i] = 0;
00676 Int_t nEast=0;
00677 Int_t nWest=0;
00678 for(Int_t i=0;i<(Int_t)mcBTofHits.size();i++) {
00679 LOG_DEBUG << (*mcBTofHits[i]) << endm;
00680
00681 if(mBookHisto) {
00682 Int_t itray = mcBTofHits[i]->tray();
00683 Int_t imodule = mcBTofHits[i]->module();
00684 Int_t icell = mcBTofHits[i]->cell();
00685 Float_t t0 = mcBTofHits[i]->time();
00686 Float_t tof = mcBTofHits[i]->tof();
00687 Float_t de = mcBTofHits[i]->dE();
00688
00689
00690 LOG_DEBUG << "tray# "<<itray << endm;
00691
00692
00693 if(itray>0&&itray<=120) {
00694 mCellSeen->Fill((imodule-1)*mNCell+(icell-1),itray-1);
00695 mDeSeen->Fill( de / keV );
00696 mT0Seen->Fill( t0 /1000 );
00697 mTofSeen->Fill( tof / 1000 );
00698 mTofResSeen->Fill( (tof-t0) );
00699 nCell[itray-1]++;
00700 }
00701
00702 else if(itray==121 || itray==122) {
00703 mVpdSeen->Fill((itray-121)*mNVPD + (icell-1));
00704 mVpdResSeen->Fill( (tof-t0) / nanosecond );
00705 if(itray==121) nWest++;
00706 else nEast++;
00707 }
00708 }
00709 }
00710 if(mBookHisto) {
00711 for(Int_t i=0;i<mNTray;i++) mNCellSeen->Fill(nCell[i],i);
00712 mNVpdSeen->Fill(nWest,nEast);
00713 }
00714
00715 LOG_INFO << " ==== Test TofRawDataCollection ==== " << endm;
00716 for(Int_t i=0;i<mNTray;i++) nCell[i] = 0;
00717 nEast=0;
00718 nWest=0;
00719
00720 if (mWriteStEvent){
00721 StSPtrVecBTofHit& bTofHits=mEvent->btofCollection()->tofHits();
00722 StBTofHit* bHit;
00723 for(Int_t aa=0;aa<(int)bTofHits.size();aa++){
00724 bHit=bTofHits[aa];
00725 Int_t itray=bHit->tray();
00726 Int_t imodule=bHit->module();
00727 Int_t icell=bHit->cell();
00728 if(mBookHisto) {mCellReco->Fill((imodule-1)*mNCell+(icell-1),itray-1);}
00729 }
00730
00731 if(mBookHisto) {
00732 for(Int_t i=0;i<mNTray;i++) mNCellReco->Fill(nCell[i],i);
00733 mNVpdReco->Fill(nWest,nEast);
00734 }
00735
00736 }
00737 }
00738
00739 if(Debug()) cout<<"leaving fill event"<<endl;
00740
00741 return kStOk;
00742 }
00743
00744
00745
00746
00747 IntVec StBTofSimMaker::CalcCellId(Int_t volume_id, Float_t ylocal)
00748 {
00749 IntVec cellId;
00750 Int_t ires = volume_id;
00751
00752 Int_t rileft = Int_t(ires/10/100/100);
00753 ires = ires-rileft*100*100*10;
00754 Int_t itray = Int_t(ires/10/100);
00755 ires = ires-itray*100*10;
00756 Int_t imodule = Int_t(ires/10);
00757 itray = itray + (rileft-1)*mNTray/2;
00758
00759 Int_t icell = Int_t((ylocal + mBTofPadWidth * mNCell/2) / mBTofPadWidth) + 1;
00760
00761 if(itray<=0 || itray>mNTray) itray = -1;
00762 if(imodule<=0 || imodule>mNModule) imodule = -1;
00763 if(icell<=0 || icell>mNCell) icell = -1;
00764
00765 cellId.push_back(itray);
00766 cellId.push_back(imodule);
00767 cellId.push_back(icell);
00768
00769 return cellId;
00770 }
00771
00772
00773 Double_t StBTofSimMaker::GammaRandom()
00774 {
00775 Double_t xmax,ymin,x,y,x1;
00776 xmax = 10.0;
00777 ymin = exp(-xmax);
00778
00779 back:
00780 y = ymin+(1-ymin)*gRandom->Rndm();
00781 x = -log(y);
00782 x1 = sqrt(xmax)*gRandom->Rndm();
00783 if(x1>sqrt(x)) goto back;
00784 return x/1.5;
00785
00786 }
00787
00788
00789 Int_t StBTofSimMaker::CellXtalk(Int_t icell, Float_t ylocal, Float_t& wt, Int_t& icellx)
00790 {
00791 Float_t yc = (icell-1-2.5)*mBTofPadWidth;
00792 Float_t dy = ylocal - yc;
00793
00794 wt = 1.;
00795 icellx = -1;
00796 Float_t dyCut = mSimDb->dy_xtalk();
00797 if(fabs(dy)<dyCut) return kStOk;
00798
00799 wt = 1. - (fabs(dy) - dyCut)/(mBTofPadWidth - 2.0*dyCut);
00800
00801 if(dy>0) icellx = icell + 1;
00802 else icellx = icell - 1;
00803
00804 if(icellx>mNCell) icellx = -1;
00805 if(icellx<=0) icellx = -1;
00806
00807 return kStOk;
00808 }
00809
00810
00811
00812 Int_t StBTofSimMaker::FastCellResponse(g2t_ctf_hit_st* tofHitsFromGeant)
00813 {
00814
00815 if(tofHitsFromGeant->s_track<=0.0 || tofHitsFromGeant->de <=0.0) {
00816 LOG_WARN << " No deposit energy in this tof hit! " << endm;
00817 return kStWarn;
00818 }
00819
00820 if(mBookHisto) {
00821 mDeGeant->Fill(tofHitsFromGeant->de / keV);
00822 mTofGeant->Fill(tofHitsFromGeant->tof / nanosecond);
00823 }
00824
00825 IntVec cellId = CalcCellId(tofHitsFromGeant->volume_id, tofHitsFromGeant->x[1]);
00826
00827 Int_t icell, imodule, itray;
00828 itray = cellId[0];
00829 imodule = cellId[1];
00830 icell = cellId[2];
00831 if (itray==-1 || imodule==-1 || icell==-1) {
00832 LOG_WARN << " Not hit the sensitive MRPC volume !!! " << endm;
00833 return kStWarn;
00834 }
00835
00836 StThreeVectorF local(tofHitsFromGeant->x[0], tofHitsFromGeant->x[1], tofHitsFromGeant->x[2]);
00837
00838 St_g2t_track *g2t_track = static_cast<St_g2t_track *>(mGeantData->Find("g2t_track"));
00839 if (!g2t_track) {
00840 LOG_WARN << " No g2t track table !!! " << endm;
00841 return kStWarn;
00842 }
00843 g2t_track_st *tof_track = g2t_track->GetTable();
00844 Int_t no_tracks= g2t_track->GetNRows();
00845
00846 StMcTrack *partnerTrk = 0;
00847 Int_t partnerTrkId;
00848 for(Int_t j=0;j<no_tracks;j++){
00849 if(tofHitsFromGeant->track_p==tof_track[j].id){
00850 partnerTrk = new StMcTrack(&(tof_track[j]));
00851 partnerTrkId=partnerTrk->key();
00852 }
00853 }
00854
00856 Int_t icellx = -1;
00857 Float_t wt = 1.0;
00858 if(mCellXtalk) CellXtalk(icell, local.y(), wt, icellx);
00859
00860 Double_t tof= tofHitsFromGeant->tof*1000./nanosecond + ranGauss.shoot()*mSimDb->timeres_tof()*1000./nanosecond;
00861 Double_t t0 = tofHitsFromGeant->tof*1000./nanosecond;
00862 Double_t de = tofHitsFromGeant->de * wt;
00863 Double_t pathL = tofHitsFromGeant->s_track;
00864 Double_t q = 0.;
00865
00866 StMcBTofHit *mcBTofHit = new StMcBTofHit(itray,imodule,icell,de,pathL,t0,tof,q);
00867 mcBTofHit->setPosition(local);
00868 mcBTofHit->setParentTrack(partnerTrk);
00869 storeMcBTofHit(mcBTofHit);
00870
00871 mTofHitFlag[itray-1][(imodule-1)*mNCell+(icell-1)] = 1;
00872
00873 if(icellx <= 0 || icellx > mNCell) return kStOk;
00874
00875
00876
00877 Double_t tofx = tofHitsFromGeant->tof + ranGauss.shoot()*mSimDb->timeres_tof();
00878 Double_t dex = tofHitsFromGeant->de * (1. - wt);
00879 Double_t qx = 0.*(1.-wt);
00880
00881 StMcBTofHit *mcBTofHitx = new StMcBTofHit(itray,imodule,icellx,dex,pathL,t0,tofx,qx);
00882 mcBTofHitx->setPosition(local);
00883 mcBTofHitx->setParentTrack(partnerTrk);
00884 mcBTofHitx->setParentTrackId(partnerTrkId);
00885 storeMcBTofHit(mcBTofHitx);
00886
00887 mTofHitFlag[itray-1][(imodule-1)*mNCell+(icellx-1)] = 1;
00888
00889 return kStOk;
00890 }
00891
00892
00893 Int_t StBTofSimMaker::storeMcBTofHit(StMcBTofHit* mcBTofHit)
00894 {
00895
00896
00897 Bool_t hitFound = kFALSE;
00898
00899
00900 for(size_t j=0;j<mMcBTofHitCollection->hits().size();j++) {
00901 StMcBTofHit *tempHit = mMcBTofHitCollection->hits()[j];
00902 if(!tempHit) continue;
00903 if(mcBTofHit->sameCell(*tempHit)) {
00904 hitFound = kTRUE;
00905 Float_t t1 = mcBTofHit->time();
00906 Float_t t2 = tempHit->time();
00907 Float_t tof1 = mcBTofHit->tof();
00908 Float_t dE1 = mcBTofHit->dE();
00909 Float_t dE2 = tempHit->dE();
00910 Float_t s1 = mcBTofHit->pathLength();
00911 Float_t q1 = mcBTofHit->charge();
00912 Float_t q2 = tempHit->charge();
00913 StThreeVectorF x1 = mcBTofHit->position();
00914 StThreeVectorF x2 = tempHit->position();
00915 StMcTrack *trk1 = mcBTofHit->parentTrack();
00916 if(t1>t2) {
00917
00918 } else {
00919 tempHit->setTime(t1);
00920 tempHit->setTof(tof1);
00921 tempHit->setPathLength(s1);
00922 tempHit->setPosition(x1);
00923 tempHit->setParentTrack(trk1);
00924 }
00925 tempHit->setdE(dE1+dE2);
00926 tempHit->setCharge(q1+q2);
00927 }
00928 }
00929
00930 if(!hitFound) {
00931 mMcBTofHitCollection->addHit(mcBTofHit);
00932 } else {
00933 delete mcBTofHit;
00934 }
00935 return kStOk;
00936 }
00937
00939 Int_t StBTofSimMaker::fillRaw(){
00940
00941
00942 return kStOk;
00943
00944 }
00945
00946
00948 Int_t StBTofSimMaker::electronicNoise(){
00949 return kStOk;
00950 }
00951
00952
00953 void StBTofSimMaker::setHistFileName(string s){
00954 mHistFile=s;
00955 return;
00956 }
00957
00958 Int_t StBTofSimMaker::bookHistograms()
00959 {
00960
00961 mBetaHist=new TH1F("mBetaHist","mBetaHist", 100, -2, 2);
00962 mPathLHist=new TH1F("mPathLHist","mPathLHist", 100, -2, 500);
00963 mTofHist=new TH1F("mTofHist","mTofHist", 1000, -10, 10000);
00964 mRecMass=new TH1F("mRecMass","mRecMass", 1000, -2, 4);
00965
00966 mCellGeant = new TH2F("CellGeant","CellGeant",192,0.,192.,120,1.,120.);
00967 mVpdGeant = new TH1F("VpdGeant","VpdGeant",38,0.,38.);
00968 mNCellGeant = new TH2F("NCellGeant","NCellGeant",192,0.,192.,120,1.,120.);
00969 mNVpdGeant = new TH2F("NVpdGeant","NVpdGeant",19,0.,19.,19,0.,19.);
00970 mDeGeant = new TH1F("DeGeant","DeGeant",1000,0.,10.);
00971 mTofGeant = new TH1F("TofGeant","TofGeant",1000,0.,20.);
00972
00973 mCellSeen = new TH2F("CellSeen","CellSeen",192,0.,192.,120,1.,120.);
00974 mVpdSeen = new TH1F("VpdSeen","VpdSeen",38,0.,38.);
00975 mNCellSeen = new TH2F("NCellSeen","NCellSeen",192,0.,192.,120,1.,120.);
00976 mNVpdSeen = new TH2F("NVpdSeen","NVpdSeen",19,0.,19.,19,0.,19.);
00977 mDeSeen = new TH1F("DeSeen","DeSeen",1000,0.,10.);
00978 mT0Seen = new TH1F("T0Seen","T0Seen",1000,0.,20.);
00979 mTofSeen = new TH1F("TofSeen","TofSeen",1000,0.,20.);
00980
00981 mTofResSeen = new TH1F("TofResSeen","TofResSeen",1001,-500.,500.);
00982 mVpdResSeen = new TH1F("VpdResSeen","VpdResSeen",1001,-500.,500.);
00983
00984 mCellReco = new TH2F("CellReco","CellReco",192,0.,192.,120,1.,120.);
00985 mVpdReco = new TH1F("VpdReco","VpdReco",38,0.,38.);
00986 mNCellReco = new TH2F("NCellReco","NCellReco",192,0.,192.,120,1.,120.);
00987 mNVpdReco = new TH2F("NVpdReco","NVpdReco",19,0.,19.,19,0.,19.);
00988 mADCReco = new TH1F("ADCReco","ADCReco",4096,0.,4096.);
00989 mTDCReco = new TH1F("TDCReco","TDCReco",4096,0.,4096.);
00990 mT0Reco = new TH1F("T0Reco","T0Reco",1000,0.,20.);
00991 mTofResReco = new TH1F("TofResReco","TofResReco",1000,-300.,300.);
00992 mVpdResReco = new TH1F("VpdResReco","VpdResReco",1000,-100.,100.);
00993 mTACorr = new TH2F("TACorr","TACorr",512,0.,4096.,512,0.,4096.);
00994
00995 mModHist = new TH1F("ModuleHist","ModuleHist",201,-100,100);
00996 return kStOk;
00997
00998 }
00999
01000
01001 Int_t StBTofSimMaker::writeHistograms()
01002 {
01003
01004
01005 mBetaHist->Write();
01006 mPathLHist->Write();
01007 mTofHist->Write();
01008 mRecMass->Write();
01009
01010 mCellGeant->Write();
01011 mVpdGeant->Write();
01012 mNCellGeant->Write();
01013 mNVpdGeant->Write();
01014 mDeGeant->Write();
01015 mTofGeant->Write();
01016
01017 mCellSeen->Write();
01018 mVpdSeen->Write();
01019 mNCellSeen->Write();
01020 mNVpdSeen->Write();
01021 mDeSeen->Write();
01022 mT0Seen->Write();
01023 mTofSeen->Write();
01024 mTofResSeen->Write();
01025 mVpdResSeen->Write();
01026
01027 mCellReco->Write();
01028 mVpdReco->Write();
01029 mNCellReco->Write();
01030 mNVpdReco->Write();
01031 mADCReco->Write();
01032 mTDCReco->Write();
01033 mT0Reco->Write();
01034 mTofResReco->Write();
01035 mVpdResReco->Write();
01036 mTACorr->Write();
01037
01038 mModHist->Write();
01039 return kStOk;
01040 }
01041
