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
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
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
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 #include <string>
00139 #include <cmath>
00140 using namespace std;
00141
00142 #include "St_DataSetIter.h"
00143 #include "St_ObjectSet.h"
00144 #include "StMessMgr.h"
00145 #include "TH1.h"
00146 #include "TH2.h"
00147 #include "TString.h"
00148
00149 #include "TFile.h"
00150 #include "TDirectory.h"
00151 #include "TNtuple.h"
00152
00153 #include "StSequence.hh"
00154 #include "StDbUtilities/StSvtCoordinateTransform.hh"
00155 #include "StDbUtilities/StGlobalCoordinate.hh"
00156 #include "StDbUtilities/StSvtLocalCoordinate.hh"
00157 #include "StDbUtilities/StSvtWaferCoordinate.hh"
00158
00159 #include "StSvtClassLibrary/StSvtHybridData.hh"
00160 #include "StSvtClassLibrary/StSvtData.hh"
00161 #include "StSvtClassLibrary/StSvtHybridPixelsD.hh"
00162 #include "StSvtClassLibrary/StSvtHybridPixelsC.hh"
00163 #include "StSvtClassLibrary/StSvtConfig.hh"
00164 #include "StSvtClassLibrary/StSvtGeometry.hh"
00165 #include "StSvtClassLibrary/StSvtWaferGeometry.hh"
00166 #include "StSvtClassLibrary/StSvtT0.hh"
00167 #include "StSvtClassLibrary/StSvtHybridDriftVelocity.hh"
00168 #include "StSvtClassLibrary/StSvtHybridPed.hh"
00169 #include "StSvtClassLibrary/StSvtDaq.hh"
00170 #include "StSvtDaqMaker/StSvtHybridDaqPed.hh"
00171 #include "StSvtCalibMaker/StSvtPedMaker.h"
00172 #include "StSvtAngles.hh"
00173 #include "StSvtElectronCloud.hh"
00174 #include "StSvtSignal.hh"
00175 #include "StSvtSimulation.hh"
00176 #include "StSvtGeantHits.hh"
00177
00178 #include "StThreeVectorD.hh"
00179 #include "StPhysicalHelixD.hh"
00180 #include "SystemOfUnits.h"
00181
00182 #include "tables/St_g2t_svt_hit_Table.h"
00183 #include "tables/St_g2t_track_Table.h"
00184
00185
00186 #include "StSvtSimulationMaker.h"
00187 #include "StDbUtilities/St_svtCorrectionC.h"
00188 #include "StDbUtilities/St_svtHybridDriftVelocityC.h"
00189 #include "StMcEvent/StMcSvtHitCollection.hh"
00190 #include "StParticleDefinition.hh"
00191 #include "StMcEvent.hh"
00192 #include "StMcTrack.hh"
00193 #include "StMcSvtHit.hh"
00194 #define gufld gufld_
00195 extern "C" {void gufld(Float_t *, Float_t *);}
00196
00197 ClassImp(StSvtSimulationMaker)
00198
00199
00207 #define cTimeBinSize 0.04
00208 #define cDefaultDriftVelocity 1.E-5*675000
00209 #define cDiffusionConst 0.0035
00210 #define cLifeTime 1000000.0
00211
00212 #define cTrapConst 5.0e-5 // Stephens tuning to real data
00213
00215 StSvtSimulationMaker::StSvtSimulationMaker(const char *name):StMaker(name)
00216 , mLifeTime(cLifeTime)
00217 , mTrapConst(cTrapConst)
00218 , mDiffusionConst(cDiffusionConst)
00219 , mTimeBinSize(cTimeBinSize)
00220 , mAnodeSize(-1956)
00221 , mPedOffset(-1956)
00222 , mSigOption(0)
00223 , mDefaultDriftVelocity(-1956)
00224 , mBField(-1956)
00225 , mConfig(0)
00226 , mDoEmbedding(kFALSE)
00227 , mSvtAngles(0)
00228 , mElectronCloud(0)
00229 , mSvtSimulation(0)
00230 , mCoordTransform(0)
00231 , mSvtGeom(0)
00232 , mDriftSpeedColl(0)
00233 , mT0(0)
00234 , mSvtSimPixelColl(0)
00235 , mSvt8bitPixelColl(0)
00236 , mSvtGeantHitColl(0)
00237 , counter(0)
00238 , mNtFile(0)
00239 , mNTuple(0)
00240 {
00241 LOG_DEBUG << "StSvtSimulationMaker::constructor"<<endm;
00242
00243
00244 mNumOfHybrids = 0;
00245 LOG_DEBUG << "StSvtSimulationMaker::constructor...END"<<endm;
00246 }
00247
00248
00249 StSvtSimulationMaker::~StSvtSimulationMaker()
00250 {
00251 LOG_DEBUG << "StSvtSimulationMaker::destructor"<<endm;
00252
00253 SafeDelete(mSvtAngles);
00254 SafeDelete(mSvtSimulation);
00255 SafeDelete(mElectronCloud);
00256 SafeDelete(mCoordTransform);
00257
00258 LOG_DEBUG << "StSvtSimulationMaker::destructor...END"<<endm;
00259 }
00260
00261
00262
00263 Int_t StSvtSimulationMaker::setConst(double timBinSize, double anodeSize)
00264 {
00265 mTimeBinSize = timBinSize ;
00266 mAnodeSize = anodeSize;
00267 return kStOK;
00268 }
00269
00270
00271
00272
00273 void StSvtSimulationMaker::setElectronLifeTime(double tLife)
00274 {
00275 mLifeTime = tLife;
00276
00277 }
00278
00279 void StSvtSimulationMaker::setTrappingConst(double trapConst)
00280 {
00281 mTrapConst = trapConst;
00282 }
00283
00284
00285 void StSvtSimulationMaker::setDiffusionConst(double diffConst)
00286 {
00287 mDiffusionConst = diffConst;
00288 }
00289
00290
00291
00292 Int_t StSvtSimulationMaker::setOptions(int SigOption)
00293 {
00294 mSigOption = SigOption;
00295
00296 return kStOK;
00297 }
00298
00299
00300
00301 Int_t StSvtSimulationMaker::Init()
00302 {
00303
00304 LOG_DEBUG << "In StSvtSimulationMaker::Init() ..."<<endm;
00305
00306
00307 LOG_DEBUG << "In StSvtSimulationMaker::Init() -End"<<endm;
00308
00309 return StMaker::Init();
00310 }
00311
00312
00314 Int_t StSvtSimulationMaker::InitRun(int runumber)
00315 {
00316 LOG_DEBUG << "StSvtSimulationMaker::InitRun()"<<endm;
00317
00318
00319 Int_t res;
00320 getConfig();
00321 if ((res=getSvtGeometry())!=kStOk) return res;
00322 #if 0
00323 if ((res=getSvtDriftSpeeds())!=kStOk) return res;
00324 #endif
00325 if ((res=getSvtT0())!=kStOk) return res;
00326 if ((res=getPedestalOffset())!=kStOk) return res;
00327
00328 setSvtPixelData();
00329
00330 if (! mCoordTransform) mCoordTransform=new StSvtCoordinateTransform();
00331
00332 mCoordTransform->setParamPointers(mSvtGeom, mConfig,mDriftSpeedColl,mT0);
00333
00334
00335
00336
00337 mTimeBinSize = cTimeBinSize;
00338 mAnodeSize = mSvtGeom->getAnodePitch()*10;
00339 mDefaultDriftVelocity = cDefaultDriftVelocity;
00340
00341 if (! mElectronCloud) mElectronCloud = new StSvtElectronCloud();
00342 mElectronCloud->setDriftVelocity(mDefaultDriftVelocity);
00343 mElectronCloud->setElectronLifeTime(mLifeTime);
00344 mElectronCloud->setDiffusionConst(mDiffusionConst);
00345 mElectronCloud->setTrappingConst(mTrapConst);
00346 if (! mSvtSimulation) mSvtSimulation = new StSvtSimulation();
00347 mSvtSimulation->setOptions(mSigOption);
00348 mSvtSimulation->setAnodeTimeBinSizes(mTimeBinSize , mAnodeSize);
00349 mSvtSimulation->setDriftVelocity(mDefaultDriftVelocity);
00350 mSvtSimulation->setElCloud(mElectronCloud);
00351
00352 LOG_INFO <<"StSvtSimulationMaker::InitRun info:"<<endm;
00353 LOG_INFO <<" Anode size="<<mAnodeSize<<" ,time bin size="<<mTimeBinSize<<endm;
00354 LOG_INFO <<" default drift velocity="<<mDefaultDriftVelocity<<endm;
00355 LOG_INFO <<" pedestal offset(from database)="<<mPedOffset<<endm;
00356 LOG_INFO <<" T0(from database)= "<<mT0->getT0(1)<<endm;
00357
00358
00359 Float_t x[3] = {0,0,0};
00360 Float_t b[3];
00361 gufld(x,b);
00362 mBField = b[2]*tesla;
00363
00364 LOG_DEBUG << "StSvtSimulationMaker::InitRun()-END"<<endm;
00365
00366 return StMaker::InitRun(runumber);
00367 }
00368
00369
00370
00371 Int_t StSvtSimulationMaker:: FinishRun(int oldrunumber){
00372 LOG_INFO <<"StSvtSimulationMaker::FinishRun()"<<endm;
00373
00374 TDataSet *set;
00375 if ((set=GetDataSet("StSvtPixelData"))) delete set;
00376 if ((set=GetDataSet("StSvt8bitPixelData"))) delete set;
00377
00378 LOG_INFO <<"StSvtSimulationMaker::FinishRun() - END"<<endm;
00379 return StMaker::FinishRun(oldrunumber);
00380 }
00381
00382
00383
00384
00385 void StSvtSimulationMaker::resetPixelData(){
00386
00387
00388 StSvtHybridPixelsD* tmpPixels;
00389 StSvtHybridPixelsC* tmp8bitPixels;
00390 if (mSvtSimPixelColl) {
00391 for(int Barrel = 1;Barrel <= mSvtSimPixelColl->getNumberOfBarrels();Barrel++) {
00392 for (int Ladder = 1;Ladder <= mSvtSimPixelColl->getNumberOfLadders(Barrel);Ladder++) {
00393 for (int Wafer = 1;Wafer <= mSvtSimPixelColl->getNumberOfWafers(Barrel);Wafer++) {
00394 for( int Hybrid = 1;Hybrid <= mSvtSimPixelColl->getNumberOfHybrids();Hybrid++){
00395
00396 int index = mSvtSimPixelColl->getHybridIndex(Barrel,Ladder,Wafer,Hybrid);
00397 if( index < 0) continue;
00398
00399 tmpPixels = (StSvtHybridPixelsD*)mSvtSimPixelColl->at(index);
00400 tmp8bitPixels = (StSvtHybridPixelsC*)mSvt8bitPixelColl->at(index);
00401
00402 if(!tmpPixels) {
00403 tmpPixels = new StSvtHybridPixelsD(Barrel, Ladder, Wafer, Hybrid);
00404 mSvtSimPixelColl->put_at(tmpPixels,index);
00405 }
00406 if(!tmp8bitPixels) {
00407 tmp8bitPixels = new StSvtHybridPixelsC(Barrel, Ladder, Wafer, Hybrid);
00408 mSvt8bitPixelColl->put_at(tmp8bitPixels,index);
00409 }
00410
00411 tmpPixels->setPedOffset(mPedOffset);
00412 tmpPixels->reset();
00413
00414 }
00415 }
00416 }
00417 }
00418 }
00419 }
00420
00421
00422
00424 void StSvtSimulationMaker::setSvtPixelData()
00425 {
00426 if (GetDataSet("StSvtPixelData"))
00427 LOG_ERROR <<"Error: Found StSvtSimPIxels in the chain - should have been deleted"<<endm;
00428
00429 St_ObjectSet *set = new St_ObjectSet("StSvtPixelData");
00430 AddConst(set);
00431 mSvtSimPixelColl = new StSvtData(mConfig->getConfiguration());
00432 set->SetObject((TObject*)mSvtSimPixelColl);
00433
00434 set = new St_ObjectSet("StSvt8bitPixelData");
00435 AddConst(set);
00436 mSvt8bitPixelColl = new StSvtData(mConfig->getConfiguration());
00437 set->SetObject((TObject*)mSvt8bitPixelColl);
00438
00439 mNumOfHybrids = mSvtSimPixelColl->getTotalNumberOfHybrids();
00440 }
00441
00442
00443
00444 void StSvtSimulationMaker::setGeantData()
00445 {
00446 St_ObjectSet* set=(St_ObjectSet*)GetDataSet("StSvtGeantHits");
00447
00448 if (set) {
00449 LOG_DEBUG <<"Found StSvtGeantHits in the chain- replacing"<<endm;
00450 set->SetObject(0);
00451 }
00452 else{
00453 set = new St_ObjectSet("StSvtGeantHits");
00454 AddData(set);
00455 }
00456
00457 if (mSvtGeantHitColl) {
00458 LOG_ERROR <<"!!!!!!m SvtGeantHitColl already exists in SvtSimulationMaker.cxx:setEval"<<endm;
00459 } else{
00460
00461 mSvtGeantHitColl = new StSvtData(mConfig->getConfiguration());
00462 set->SetObject((TObject*)mSvtGeantHitColl);
00463 }
00464
00465
00466
00467 if(!counter){
00468 counter = new int[mNumOfHybrids];
00469 for( int ii=0; ii<mNumOfHybrids; ii++){
00470 counter[ii] = 0;
00471 }
00472 }
00473 }
00474
00475
00476 Int_t StSvtSimulationMaker::getSvtGeometry()
00477 {
00478
00479 St_DataSet* dataSet;
00480 dataSet = GetDataSet("StSvtGeometry");
00481 if (!dataSet){
00482 LOG_ERROR <<"BIG TROUBLE:No SVT geometry -impossible to run!!!!"<<endm;
00483 return kStFatal;
00484 }
00485
00486
00487 mSvtGeom = (StSvtGeometry*)dataSet->GetObject();
00488 if (!mSvtGeom){
00489 LOG_ERROR << "BIG TROUBLE:No SVT geometry -impossible to run!!!!"<<endm;
00490 return kStFatal;
00491 }
00492
00493
00494 return kStOk;
00495 }
00496
00497
00498 Int_t StSvtSimulationMaker::getPedestalOffset()
00499 {
00500
00501 St_DataSet* dataSet;
00502 dataSet = GetDataSet("StSvtDaq");
00503 if (!dataSet){
00504 LOG_ERROR <<"BIG TROUBLE:No DAQ parameters for SVT!!!!"<<endm;
00505 return kStFatal;
00506 }
00507
00508
00509 StSvtDaq *daq = (StSvtDaq*)dataSet->GetObject();
00510 if (daq==NULL){
00511 LOG_ERROR << "BIG TROUBLE:No DAQ parameters are empty!!!!"<<endm;
00512 return kStFatal;
00513 }
00514
00515
00516 mPedOffset=daq->getPedOffset();
00517
00518 return kStOk;
00519 }
00520 #if 0
00521
00522 Int_t StSvtSimulationMaker::getSvtDriftSpeeds()
00523 {
00524 mDriftSpeedColl =NULL;
00525 St_DataSet* dataSet;
00526 dataSet = GetDataSet("StSvtDriftVelocity");
00527 if (!dataSet){
00528 cout<<"Warning: no SVT drift velocity data available - using default drift speed:"<<mDefaultDriftVelocity<<endl;
00529 return kStErr;
00530 }
00531
00532 mDriftSpeedColl = (StSvtHybridCollection*)dataSet->GetObject();
00533 if (! mDriftSpeedColl) cout<<"Warning: SVT drift velocity data empty - using default drift speed:"<<mDefaultDriftVelocity<<endl;
00534
00535 return kStOk;
00536 }
00537 #endif
00538
00539
00540 Int_t StSvtSimulationMaker::getSvtT0()
00541 {
00542 mT0=NULL;
00543 St_DataSet* dataSet;
00544 dataSet = GetDataSet("StSvtT0");
00545 if (!dataSet){
00546 LOG_WARN <<"no SVT T0 data available -using defalt T0 = 0"<<endm;
00547 return kStErr;
00548 }
00549
00550 mT0 = (StSvtT0*)dataSet->GetObject();
00551 if (! mT0) LOG_WARN <<" SVT T0 data empty - using default T0 = 0"<<endm;
00552
00553 return kStOk;
00554 }
00555
00556
00557 Int_t StSvtSimulationMaker::getConfig()
00558 {
00559 mConfig=NULL;
00560 St_DataSet *dataSet = NULL;
00561 dataSet = GetDataSet("StSvtConfig");
00562
00563 if (!dataSet)
00564 {
00565 LOG_WARN <<" No SvtConfig data set" << endm;
00566 dataSet = new St_ObjectSet("StSvtConfig");
00567 AddConst(dataSet);
00568 mConfig=NULL;
00569 }
00570
00571 mConfig=((StSvtConfig*)(dataSet->GetObject()));
00572
00573 if (!mConfig) {
00574 LOG_WARN <<"SvtConfig data set is empty- seting default full configuration" << endm;
00575 mConfig=new StSvtConfig();
00576 mConfig->setConfiguration("FULL");
00577 dataSet->SetObject(mConfig);
00578 }
00579
00580 return kStOk;
00581 }
00582
00583
00584 Int_t StSvtSimulationMaker::ideal2RealTranslation( StThreeVector<double> *pos, StThreeVector<double> *mtm, double charge, int *wafId){
00585 StGlobalCoordinate globalCor(0,0,0);
00586 StThreeVector<double> x;
00587
00588 int index = mSvtGeom->getWaferIndex(*wafId);
00589 if (index < 0) {
00590 LOG_DEBUG << "StSvtSimulationMaker::ideal2RealTranslation: illegal wafId = " << *wafId << endl;
00591 return kStSkip;
00592 }
00593
00594
00595
00596 StSvtWaferGeometry* waferGeom = (StSvtWaferGeometry*)mSvtGeom->at(index);
00597 if (! waferGeom) return kStSkip;
00598 StThreeVectorD wafCent(waferGeom->x(0),waferGeom->x(1),waferGeom->x(2));
00599 StThreeVectorD norm(waferGeom->n(0),waferGeom->n(1),waferGeom->n(2));
00600
00601
00602
00603 StPhysicalHelixD tHelix( *mtm, *pos, mBField, charge);
00604 double s = tHelix.pathLength(wafCent, norm);
00605 x = tHelix.at(s);
00606 pos->setX(x.x());
00607 pos->setY(x.y());
00608 pos->setZ(x.z());
00609 globalCor.setPosition(*pos);
00610
00611
00612 if( mCoordTransform->IsOnWaferZ(*pos,*wafId) && mCoordTransform->IsOnWaferR(*pos,*wafId)){
00613
00614 x = tHelix.momentumAt(s,mBField);
00615 mtm->setX(x.x());
00616 mtm->setY(x.y());
00617 mtm->setZ(x.z());
00618 return kStOk;
00619 }
00620
00621
00622
00623
00624
00625 int ladder = *wafId%100;
00626 int layer = *wafId/1000;
00627 int barrel = (layer-1)/2 +1;
00628 int iladder;
00629 for( iladder = ladder-1; iladder <= ladder+1; iladder++){
00630 if( iladder==0) continue;
00631 if( iladder > mConfig->getNumberOfLadders(barrel)) break;
00632 for( int iwaf = 1; iwaf <= mConfig->getNumberOfWafers(barrel); iwaf++){
00633
00634 index = mSvtGeom->getWaferIndex(barrel, iladder, iwaf);
00635
00636
00637 waferGeom = (StSvtWaferGeometry*)mSvtGeom->at(index);
00638 if (! waferGeom) continue;
00639 wafCent.setX(waferGeom->x(0));
00640 wafCent.setY(waferGeom->x(1));
00641 wafCent.setZ(waferGeom->x(2));
00642 norm.setX(waferGeom->n(0));
00643 norm.setY(waferGeom->n(1));
00644 norm.setZ(waferGeom->n(2));
00645
00646
00647 s = tHelix.pathLength(wafCent, norm);
00648 x = tHelix.at(s);
00649 pos->setX(x.x());
00650 pos->setY(x.y());
00651 pos->setZ(x.z());
00652 globalCor.setPosition(*pos);
00653 *wafId = 1000*mConfig->getLayer(index) + 100*iwaf + iladder;
00654 if( mCoordTransform->IsOnWaferZ(*pos,*wafId) && mCoordTransform->IsOnWaferR(*pos,*wafId)){
00655
00656 x = tHelix.momentumAt(s,mBField);
00657 mtm->setX(x.x());
00658 mtm->setY(x.y());
00659 mtm->setZ(x.z());
00660 return kStOk;
00661 }
00662 }
00663 }
00664
00665 if( ladder ==1){
00666 iladder = mConfig->getNumberOfLadders(barrel);
00667 for( int iwaf = 1; iwaf <= mConfig->getNumberOfWafers(barrel); iwaf++){
00668
00669 index = mSvtGeom->getWaferIndex(barrel, iladder, iwaf);
00670
00671
00672 waferGeom = (StSvtWaferGeometry*)mSvtGeom->at(index);
00673 if (! waferGeom) continue;
00674 wafCent.setX(waferGeom->x(0));
00675 wafCent.setY(waferGeom->x(1));
00676 wafCent.setZ(waferGeom->x(2));
00677 norm.setX(waferGeom->n(0));
00678 norm.setY(waferGeom->n(1));
00679 norm.setZ(waferGeom->n(2));
00680
00681
00682 s = tHelix.pathLength(wafCent, norm);
00683 x = tHelix.at(s);
00684 pos->setX(x.x());
00685 pos->setY(x.y());
00686 pos->setZ(x.z());
00687 globalCor.setPosition(*pos);
00688 *wafId = 1000*mConfig->getLayer(index) + 100*iwaf + iladder;
00689 if( mCoordTransform->IsOnWaferZ(*pos,*wafId) &&
00690 mCoordTransform->IsOnWaferR(*pos,*wafId)){
00691
00692 x = tHelix.momentumAt(s,mBField);
00693 mtm->setX(x.x());
00694 mtm->setY(x.y());
00695 mtm->setZ(x.z());
00696 return kStOk;
00697 }
00698 }
00699 }
00700
00701 return kStSkip;
00702 }
00703
00704 Int_t StSvtSimulationMaker::Make()
00705 {
00706 LOG_DEBUG << "In StSvtSimulationMaker::Make()" << endm;
00707
00708 int volId ,barrel, layer, ladder, wafer, hybrid;
00709
00710 StThreeVector<double> VecG(0,0,0);
00711 StThreeVector<double> VecL(0,0,0);
00712 StThreeVector<double> mtm(0,0,0);
00713
00714
00715 StSvtHybridPixelsD *svtSimDataPixels;
00716
00717
00718 resetPixelData();
00719 setGeantData();
00720
00721 StSvtWaferCoordinate waferCoord (0,0,0,0,0,0);
00722 StSvtLocalCoordinate localCoord (0,0,0), localCoordb(0,0,0);
00723 StGlobalCoordinate globalCor(0,0,0);
00724
00725 St_svtHybridDriftVelocityC *driftVel = St_svtHybridDriftVelocityC::instance();
00726 assert(driftVel);
00727
00728 int tmpBadCount=0;
00729
00730 StMcEvent* mcEvent = 0;
00731 mcEvent = (StMcEvent*) GetDataSet("StMcEvent");
00732 StMcSvtHitCollection *mcCol = 0;
00733
00734 if(mcEvent)
00735 {
00736 mcCol = mcEvent->svtHitCollection();
00737 LOG_INFO <<"mNumOfGeantHits = "<<mcCol->numberOfHits()<<endm;
00738 for (unsigned int iBarrel=0; iBarrel< mcCol->numberOfBarrels(); iBarrel++) {
00739 for (unsigned int iLadder=0; iLadder<mcCol->barrel(iBarrel)->numberOfLadders(); iLadder++) {
00740 for (unsigned int iWafer = 0; iWafer < mcCol->barrel(iBarrel)->ladder(iLadder)->numberOfWafers(); iWafer++) {
00741 for (StMcSvtHitIterator iter = mcCol->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().begin();
00742 iter != mcCol->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().end();
00743 iter++) {
00744 StMcSvtHit *mchit = dynamic_cast<StMcSvtHit *> (*iter);
00745
00746 VecG.setX( mchit->position().x());
00747 VecG.setY( mchit->position().y());
00748 VecG.setZ( mchit->position().z());
00749 mtm.setX(mchit->localMomentum().x());
00750 mtm.setY(mchit->localMomentum().y());
00751 mtm.setZ(mchit->localMomentum().z());
00752 volId = 1000*mchit->layer()+100*mchit->wafer()+mchit->ladder();
00753
00754
00755
00756
00757
00758
00759 Int_t iok = ideal2RealTranslation(&VecG, &mtm, (double)mchit->parentTrack()->particleDefinition()->charge(), &volId);
00760 if (iok != kStOK) continue;
00761 mchit->setPosition(VecG);
00762 mchit->setLocalMomentum(mtm);
00763
00764 globalCor.setPosition(VecG);
00765 mCoordTransform->operator()(globalCor,localCoord,volId);
00766
00767
00768 if (! driftVel->p(localCoord.barrel(),localCoord.ladder(),localCoord.wafer(),localCoord.hybrid())) continue;
00769
00770
00771
00772 mCoordTransform->operator()(localCoord,waferCoord);
00773 mCoordTransform->operator()(waferCoord,localCoordb);
00774 double tnew, told;
00775 int mtry=0;
00776 double scale = 1;
00777 double diffx = localCoordb.position().x()- localCoord.position().x();
00778 while( fabs(diffx) > 0.002){
00779
00780 told = waferCoord.timebucket();
00781 tnew =waferCoord.timebucket()+ scale*waferCoord.timebucket()*(localCoordb.position().x()- localCoord.position().x())/localCoord.position().x();
00782 if( told-tnew > 4) tnew = told-0.5;
00783 else if( told-tnew < -4) tnew = told+0.5;
00784 waferCoord.setTimeBucket(tnew);
00785 mCoordTransform->operator()(waferCoord,localCoordb);
00786 if( fabs(diffx) < fabs(localCoordb.position().x()- localCoord.position().x())){
00787 waferCoord.setTimeBucket(told);
00788 scale *= 0.5;
00789 }
00790 else{
00791 diffx = localCoordb.position().x()- localCoord.position().x();
00792 }
00793
00794 mtry++;
00795 if( mtry > 10000) diffx = 0;
00796
00797 }
00798
00799
00800 VecL.setX(localCoord.position().x());
00801 VecL.setY(localCoord.position().y());
00802 VecL.setZ(localCoord.position().z());
00803
00804
00805 double anode,time;
00806 barrel = waferCoord.barrel();
00807 layer = waferCoord.layer(); ladder = waferCoord.ladder();
00808 wafer = waferCoord.wafer(); hybrid = waferCoord.hybrid();
00809 time = waferCoord.timebucket();
00810 double driftTime=time - mT0->getT0();
00811 anode = waferCoord.anode();
00812
00813 if(driftTime < 0.0 || time > 128.0 || anode < 0.0 || anode > 240.0)
00814 { tmpBadCount++; continue;}
00815 #if 0
00816
00817
00818
00819
00820 static const int barrels[]={3,1,1,2,2};
00821 barrel = 3;
00822 if (layer<=4) barrel = barrels[layer];
00823
00824 if ( !strncmp(mConfig->getConfiguration(), "Y1L", strlen("Y1L")) ) {
00825 if ((wafer == 1) || (wafer == 2) || (wafer == 3))
00826 ladder = 2;
00827 }
00828
00829
00830 #endif
00831 if( 1000*layer+100*wafer+ladder != volId){
00832 LOG_INFO << "trouble - skipping hit" << volId <<"\t"<< mchit->position().x() << "\t"
00833 << mchit->position().y() <<"\t and our calc"<<"\t" << layer << " "
00834 << wafer << "\t" << ladder <<endm;
00835 continue;
00836 }
00837
00838 int index = mSvtSimPixelColl->getHybridIndex(barrel,ladder,wafer,hybrid);
00839 if( index < 0) continue;
00840 svtSimDataPixels = (StSvtHybridPixelsD*)mSvtSimPixelColl->at(index);
00841 if (! mSvtAngles) mSvtAngles = new StSvtAngles();
00842 mSvtAngles->calcAngles(mSvtGeom,mtm.x(),mtm.y(),mtm.z(),layer,ladder,wafer);
00843 double theta = mSvtAngles->getTheta();
00844 double phi = mSvtAngles->getPhi();
00845
00846 double vd=-1;
00847 #if 0
00848 if (mDriftSpeedColl && mDriftSpeedColl->at(index)){
00849 vd = ((StSvtHybridDriftVelocity*)mDriftSpeedColl->at(index))->getV3(1);
00850 if (vd<=0) vd=mDefaultDriftVelocity;
00851 else vd=vd*1e-5;
00852 }
00853
00854 #else
00855 vd = driftVel->DriftVelocity(barrel,ladder,wafer,hybrid)*1e-5;
00856 #endif
00857 if (vd <= 0) continue;
00858
00859 double PasaShift=0.22/vd*25.;
00860
00861 mSvtSimulation->setDriftVelocity(vd);
00862 double energy = mchit->dE()*1e9;
00863 mSvtSimulation->doCloud(driftTime,energy,theta,phi,mchit->parentTrack()->key());
00864 mSvtSimulation->fillBuffer(anode,time-PasaShift,svtSimDataPixels);
00865
00866 if (Debug()) {
00867 FillGeantHit(barrel,ladder,wafer,hybrid,&waferCoord,&VecG,&VecL,mSvtSimulation->getPeak(),mchit->parentTrack()->key());
00868 }
00869
00870
00871 }
00872 }
00873 }
00874 }
00875 }
00876 else
00877 LOG_INFO << "There's a problem no MC event" << endm;
00878
00879 int nSimDataPixels = mSvtSimPixelColl->size();
00880 for (int index=0;index<nSimDataPixels;index++) {
00881 svtSimDataPixels = (StSvtHybridPixelsD*)mSvtSimPixelColl->at(index);
00882 if (!svtSimDataPixels) continue;
00883 svtSimDataPixels->updateTruth();
00884 }
00885 LOG_DEBUG <<"bad hits:"<<tmpBadCount<<endm;
00886 LOG_DEBUG << "In StSvtSimulationMaker::Make()...END" << endm;
00887 return kStOK;
00888 }
00889
00890
00891
00892 void StSvtSimulationMaker::FillGeantHit(int barrel, int ladder, int wafer, int hybrid,
00893 StSvtWaferCoordinate* waferCoord,StThreeVector<double>* VecG,
00894 StThreeVector<double>* VecL, double peak,int idtrk)
00895 {
00896 StSvtGeantHits* geantHit;
00897
00898 int index = mSvtGeantHitColl->getHybridIndex(barrel,ladder,wafer,hybrid);
00899 if (index < 0) return;
00900
00901 geantHit = (StSvtGeantHits*)mSvtGeantHitColl->at(index);
00902 if(!geantHit) {
00903 geantHit = new StSvtGeantHits(barrel,ladder,wafer,hybrid);
00904 mSvtGeantHitColl->put_at(geantHit,index);
00905 }
00906
00907 geantHit->setGeantHit(counter[index],waferCoord);
00908 geantHit->setLocalCoord(counter[index],VecL);
00909 geantHit->setGlobalCoord(counter[index],VecG);
00910 geantHit->setPeak(counter[index],peak);
00911 geantHit->setTrackId(counter[index],idtrk);
00912 ++counter[index];
00913 geantHit->setNumOfHits(counter[index]);
00914 }
00915
00916
00917
00918 void StSvtSimulationMaker::Clear(const char*)
00919 {
00920
00921 LOG_DEBUG << "In StSvtSimulationMaker::Clear" << endm;
00922
00923
00924
00925 mSvtGeantHitColl = NULL;
00926
00927 if (counter) delete counter; counter=NULL;
00928
00929 StMaker::Clear();
00930
00931 LOG_DEBUG << "In StSvtSimulationMaker::Clear..END" << endm;
00932 }
00933
00934
00935
00936
00937 Int_t StSvtSimulationMaker::Finish()
00938 {
00939 LOG_DEBUG << "In StSvtSimulationMaker::Finish()"<< endm;
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950 LOG_DEBUG << "In StSvtSimulationMaker::Finish() ...END"<< endm;
00951 return kStOK;
00952 }
00953
00954
00955
00956
00957
