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 #include "Stiostream.h"
00051 #include <assert.h>
00052 #include <math.h>
00053 #include "TROOT.h"
00054 #include <TRandom.h>
00055 #include <TBrowser.h>
00056 #include <StMessMgr.h>
00057 #include "StPmdSimulatorMaker.h"
00059 #include "StPhmdCollection.h"
00060 #include "StPhmdClusterCollection.h"
00061 #include "StPhmdDetector.h"
00062 #include "StPhmdHit.h"
00063 #include "StPhmdModule.h"
00064
00065 #include "StPmdUtil/StPmdCollection.h"
00066 #include "StPmdUtil/StPmdDetector.h"
00067 #include "StPmdUtil/StPmdHit.h"
00068 #include "StPmdUtil/StPmdModule.h"
00069 #include "StPmdUtil/StPmdMapUtil.h"
00070 #include "tables/St_g2t_pmd_hit_Table.h"
00071 #include "tables/St_g2t_track_Table.h"
00072 #include <TTableSorter.h>
00073
00074 #include "StEvent/StEvent.h"
00075 #include "StEvent/StEventTypes.h"
00076
00077
00078 ClassImp(StPmdSimulatorMaker)
00079
00080 TDataSet *geaIn;
00081
00082 St_g2t_pmd_hit *g2t_pmd_hit;
00083 StPmdGeom *mpmdgeom;
00084 StPmdMapUtil *mPmdMapUtil;
00085
00086
00087
00088 StPmdSimulatorMaker::StPmdSimulatorMaker(const char *name):StMaker(name)
00089 {
00091 adcconstants();
00092 mpmdgeom=new StPmdGeom();
00093 mPmdMapUtil= new StPmdMapUtil();
00094 mPmdCollection=NULL;
00095 gMessMgr->SetLimit("StPmdSimulator",100);
00096 }
00097
00098 StPmdSimulatorMaker::~StPmdSimulatorMaker()
00099 {
00100 delete mpmdgeom;
00101 delete mPmdMapUtil;
00102 }
00103
00104 Int_t StPmdSimulatorMaker::Init()
00105 {
00106 bookHistograms();
00107 return StMaker::Init();
00108 }
00109 void StPmdSimulatorMaker::bookHistograms()
00110 {
00111 m_pmdEdep2D = new TH2F("PMDEdep2D" ,"PMD Edep (2D)",150,-150.,150.,150,-150.,150.);
00112 m_cpvEdep2D = new TH2F("CPVEdep2D" ,"CPV Edep (2D)",150,-150.,150.,150,-150.,150.);
00113
00114 mEdepPmd = new TH1F("PmdEdep" ,"Edep for PMD",1000,0.,500.);
00115 mEdepPmd_part = new TH1F("PmdEdep_part" ,"Edep Part for PMD",1000,0.,1000.);
00116 mPmdAdc = new TH1F("pmdadc" ,"ADC for PMD",1000, 0.,1000.);
00117 mHitPmd = new TH1F("PMDHit" ,"Hits for PMD",9000, -0.5,8999.5);
00118
00119 mEdepCpv = new TH1F("CpvEdep" ,"Edep for CPV",400,0.,200.);
00120 mEdepCpv_part = new TH1F("CpvEdep_part" ,"Edep Part for CPV",200,0.,200.);
00121 mCpvAdc = new TH1F("cpvadc" ,"ADC for CPV",1000, 0.,1000.);
00122 mHitCpv = new TH1F("CpvHit" ,"Hits for CPV",5000, -0.5,4999.5);
00123
00124 m_pmdsuper = new TH1F("PmdSuper" ,"Super for PMD",12,0.5,12.5);
00125 m_pmdrow = new TH2F("PMD_SupervsRow" ,"Pmd super vs row",12,0.5,12.5,100,0.5,100.5);
00126 m_pmdcol = new TH2F("PMD_SupervsCol" ,"Pmd super vs col",12,0.5,12.5,100,0.5,100.5);
00127
00128 m_cpvsuper = new TH1F("CpvSuper" ,"Super for CPV",12,0.5,12.5);
00129 m_cpvrow = new TH2F("CPV_SupervsRow" ,"Cpv super vs row",12,0.5,12.5,100,0.5,100.5);
00130 m_cpvcol = new TH2F("CPV_SupervsCol" ,"Cpv super vs col",12,0.5,12.5,100,0.5,100.5);
00131
00132 }
00133
00134
00135 Int_t StPmdSimulatorMaker::Make()
00136 {
00138 geaIn = GetDataSet("geant");
00139 if (geaIn == 0) {
00140 geaIn = GetDataSet("event/geant/Event");
00141 if (geaIn == 0) {
00142 gMessMgr->Error()<<"Geant Data didn't find in event/geant/Event and geant directories"<<endm;
00143 return kStWarn;
00144 }
00145 }
00146
00147 Int_t retpmd;
00148 retpmd = GetPmd();
00149 return retpmd;
00150 }
00151
00152 Int_t StPmdSimulatorMaker::GetPmd()
00153 {
00154 g2t_pmd_hit = (St_g2t_pmd_hit *) geaIn->Find("g2t_pmd_hit");
00155
00156
00157 if ( g2t_pmd_hit && g2t_pmd_hit->GetTableSize()>0)
00158 {
00159
00160 Int_t phit=makePmdHits();
00161 if(phit!=kStOK){cout<<"problem in Hit formation"<<endl;return kStWarn;}
00162 }
00163 else gMessMgr->Warning()<<
00164 "StPmdSimulatorMaker::makePmd=>table g2t_PMD_hit isn't found in dataset or size is 0 "<<
00165 geaIn->GetName()<<endm;
00166 return kStOK;
00167 }
00168
00170 Int_t StPmdSimulatorMaker::makePmdHits()
00171 {
00175 mPmdCollection = new StPmdCollection("PmdCollection");
00176 m_DataSet->Add(mPmdCollection);
00177 StPmdDetector* det0 = mPmdCollection->detector(1);
00178 StPmdDetector* det1 = mPmdCollection->detector(0);
00179
00180 g2t_pmd_hit_st *hit = g2t_pmd_hit->GetTable();
00181 Int_t nhits = g2t_pmd_hit->GetNRows();
00182
00183
00184
00185 if(nhits<=0)
00186 {
00187 return kStWarn;
00188 }
00189 for(Int_t ihit=0; ihit<nhits; ihit++,hit++) {
00190 Int_t sector=0,super=0,subdet=0,row=0,col=0;
00191 Int_t gsuper;
00192
00194
00195 Int_t decode=Decode_VolId(hit->volume_id,sector,super,subdet,row,col);
00196
00198 mpmdgeom->NModule(Int_t(sector),Int_t(super),gsuper);
00199
00201
00202 mpmdgeom->Sim2Detmap(gsuper,row,col);
00203
00205
00206 if(decode==kStOK)
00207 {
00208 StPmdHit *phit = new StPmdHit();
00209 phit->setGsuper(Int_t(gsuper));
00210 phit->setSubDetector(Int_t(subdet));
00211 phit->setRow(Int_t(row));
00212 phit->setColumn(Int_t(col));
00213 phit->setEdep((hit->de)*1000000.);
00214
00215
00216 Float_t xPMD,yPMD;
00217 if(subdet==1){
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237 if(det0->module_hit(Int_t(gsuper))==0)
00238 {
00239 det0->addHit(phit);
00240 }
00241 else
00242 {
00243
00244 StPmdHit *hitmatched=Exist(phit,det0,Int_t(gsuper));
00245 if(hitmatched)
00246 {
00247 Float_t newEdep=hitmatched->Edep()+(hit->de)*1000000.;
00248 hitmatched->setEdep(newEdep);
00249 delete phit;
00250 }
00251 else
00252 {
00253 det0->addHit(phit);
00254 }
00255 }
00256
00257 m_pmdEdep2D->Fill(xPMD,yPMD,hit->de);
00259 m_pmdsuper->Fill(Float_t(gsuper));
00260 m_pmdrow->Fill(Float_t(gsuper),Float_t(row));
00261 m_pmdcol->Fill(Float_t(gsuper),Float_t(col));
00262 }
00263
00264 if(subdet==2){
00265 if(det1->module_hit(Int_t(gsuper))==0)
00266 {
00267 det1->addHit(phit);}
00268 else
00269 {
00270
00271 StPmdHit *hitmatched=Exist(phit,det1,Int_t(gsuper));
00272 if(hitmatched)
00273 {
00274 Float_t newEdep=hitmatched->Edep()+(hit->de)*1000000.;
00275 hitmatched->setEdep(newEdep);
00276 delete phit;
00277 }
00278 else
00279 {
00280 det1->addHit(phit);
00281 }
00282 }
00284 Float_t xCPV,yCPV,etaCPV,phiCPV;
00285 mpmdgeom->DetCell_xy(gsuper,Int_t(row),Int_t(col),xCPV,yCPV,etaCPV,phiCPV);
00286 m_cpvEdep2D->Fill(xCPV,yCPV,hit->de);
00288 m_cpvsuper->Fill(Float_t(gsuper));
00289 m_cpvrow->Fill(Float_t(gsuper),Float_t(row));
00290 m_cpvcol->Fill(Float_t(gsuper),Float_t(col));
00291 }
00292 }
00293 }
00294
00295
00296 if(det0)
00297 {
00298 for(Int_t ii=1;ii<13;ii++)
00299 {
00300 FinalEdep(det0,ii);
00301 }
00302 }
00303 if(det1)
00304 {
00305 for(Int_t ii=1;ii<13;ii++)
00306 {
00307 FinalEdep(det1,ii);
00308 }
00309 }
00310
00311 for(Int_t ii=1;ii<13;ii++)
00312 {
00313 FillHistograms(det0,det1,ii);
00314 }
00316 fillStEvent(det0,det1);
00317
00318 return kStOK;
00319 }
00320
00321 void StPmdSimulatorMaker::Browse(TBrowser* b)
00322 {
00323
00324 TDataSet::Browse(b);
00325 }
00326
00328
00329 Int_t StPmdSimulatorMaker::Decode_VolId(Int_t& vol, Int_t& sector, Int_t& super,Int_t& subdet, Int_t& row, Int_t& col)
00330 {
00331
00332 Int_t temp0=vol%1000000;
00333 Int_t temp1=temp0%100000;
00334 Int_t temp2=temp1%10000;
00335 Int_t temp3=temp2%100;
00336
00337 col=temp3;
00338 row=temp2/100;
00339 subdet=temp1/10000;
00340 super=temp0/100000;
00341 sector=vol/1000000;
00342
00343 return kStOK;
00344
00345 }
00346
00347 StPmdHit* StPmdSimulatorMaker::Exist(StPmdHit* phit,StPmdDetector* pdet,Int_t id)
00348 {
00349 Int_t xpad,ypad,super;
00350 Int_t nmh=pdet->module_hit(id);
00351 StPmdModule * pmod=pdet->module(id);
00352 TIter next(pmod->Hits());
00353 StPmdHit *spmcl;
00354
00355 for(Int_t im=0; im<nmh; im++)
00356 {
00357 spmcl = (StPmdHit*)next();
00358 if(spmcl)
00359 {
00360 ypad=spmcl->Row();
00361 xpad=spmcl->Column();
00362 super = spmcl->Gsuper();
00363 if(phit->Row()==ypad && phit->Column()==xpad && phit->Gsuper() == super) return spmcl;
00364 }
00365 }
00366 return NULL;
00367 }
00368
00369
00370 void StPmdSimulatorMaker::FillHistograms(StPmdDetector* pmd_det,StPmdDetector* cpv_det,Int_t id)
00371 {
00372 StPmdModule* pmd_mod=pmd_det->module(id);
00373 Int_t nmh1=pmd_det->module_hit(id);
00374 Float_t edep_part = 0.;
00375 if(nmh1>0){
00376 TIter next(pmd_mod->Hits());
00377 StPmdHit *spmcl1;
00378 for(Int_t im=0; im<nmh1; im++)
00379 {
00380 spmcl1 = (StPmdHit*)next();
00381 if(spmcl1)
00382 {
00383 Int_t gsuper = spmcl1->Gsuper();
00384 Int_t row=spmcl1->Row();
00385 Int_t col=spmcl1->Column();
00386 Float_t edep=spmcl1->Edep();
00387 Int_t adc=spmcl1->Adc();
00388
00389 if(edep>0.)mEdepPmd->Fill(edep);
00390 if(edep>0.) edep_part = edep_part+edep;
00391 mPmdAdc->Fill(Float_t(adc));
00392 Float_t xPMD,yPMD,etaPMD,phiPMD;
00393 mpmdgeom->DetCell_xy(gsuper,Int_t(row),Int_t(col),xPMD,yPMD,etaPMD,phiPMD);
00394 m_pmdEdep2D->Fill(xPMD,yPMD,edep);
00395 m_pmdsuper->Fill(Float_t(gsuper));
00396 m_pmdrow->Fill(Float_t(gsuper),Float_t(row));
00397 m_pmdcol->Fill(Float_t(gsuper),Float_t(col));
00398 }
00399 }
00400 }
00401 if(nmh1>0)mHitPmd->Fill(nmh1);
00402
00403 if(edep_part>0) mEdepPmd_part->Fill(edep_part);
00404
00405
00406 StPmdModule* cpv_mod=cpv_det->module(id);
00407 edep_part = 0.;
00408 Int_t nmh2=cpv_det->module_hit(id);
00409 if(nmh2>0){
00410 TIter next(cpv_mod->Hits());
00411 StPmdHit *spmcl2;
00412
00413 for(Int_t im=0; im<nmh2; im++)
00414 {
00415
00416 spmcl2 = (StPmdHit*)next();
00417 if(spmcl2){
00418 Int_t gsuper = spmcl2->Gsuper();
00419 Int_t row=spmcl2->Row();
00420 Int_t col=spmcl2->Column();
00421 Float_t edep=spmcl2->Edep();
00422
00423 Float_t xCPV,yCPV,etaCPV,phiCPV;
00424 mpmdgeom->DetCell_xy(gsuper,Int_t(row),Int_t(col),xCPV,yCPV,etaCPV,phiCPV);
00425
00426 mEdepCpv->Fill(edep);
00427 if(edep>0.4) edep_part = edep_part+edep;
00428
00429 mpmdgeom->DetCell_xy(gsuper,Int_t(row),Int_t(col),xCPV,yCPV,etaCPV,phiCPV);
00430
00431 m_cpvEdep2D->Fill(xCPV,yCPV,edep);
00432 m_cpvsuper->Fill(Float_t(gsuper));
00433 m_cpvrow->Fill(Float_t(gsuper),Float_t(row));
00434 m_cpvcol->Fill(Float_t(gsuper),Float_t(col));
00435 }
00436 }
00437 }
00438 if(nmh2>0)mHitCpv->Fill(nmh2);
00439
00440 if(edep_part>0) mEdepCpv_part->Fill(edep_part);
00441 }
00442
00443
00444
00445 Int_t StPmdSimulatorMaker::fillStEvent(StPmdDetector* pmd_det, StPmdDetector* cpv_det)
00446 {
00447 StEvent *currevent = (StEvent*)GetInputDS("StEvent");
00448 mevtPmdCollection = new StPhmdCollection();
00449 currevent->setPhmdCollection(mevtPmdCollection);
00450
00451 StPhmdDetector* evtdet0 = mevtPmdCollection->detector(StDetectorId(kPhmdId));
00452 StPhmdDetector* evtdet1 = mevtPmdCollection->detector(StDetectorId(kPhmdCpvId));
00453
00454
00455 for(Int_t id=1;id<13;id++)
00456 {
00457
00458
00459 Int_t subdet=1;
00460
00461
00462 StPmdModule * pmd_mod=pmd_det->module(id);
00463 Int_t nmh1=pmd_det->module_hit(id);
00464 if(nmh1>0)
00465 {
00466 TIter next(pmd_mod->Hits());
00467 StPmdHit *spmcl1;
00468 for(Int_t im=0; im<nmh1; im++)
00469 {
00470 spmcl1 = (StPmdHit*)next();
00471 if(spmcl1)
00472 {
00473 Int_t gsuper=spmcl1->Gsuper();
00474 Int_t col=spmcl1->Column();
00475 Int_t row=spmcl1->Row();
00476 Float_t edep=spmcl1->Edep();
00477 Int_t adc=spmcl1->Adc();
00479 StPhmdHit *phit = new StPhmdHit();
00480 phit->setSuperModule(Int_t(gsuper-1));
00481 phit->setSubDetector(Int_t(subdet));
00482 phit->setRow(Int_t(row));
00483 phit->setColumn(Int_t(col));
00484 phit->setEnergy(edep);
00485 phit->setAdc(adc);
00486 evtdet0->addHit(phit);
00487 }
00488 }
00489 }
00490
00491
00492 StPmdModule * cpv_mod=cpv_det->module(id);
00493 Int_t nmh2=cpv_det->module_hit(id);
00494 if(nmh2>0)
00495 {
00496 subdet=2;
00497 TIter next(cpv_mod->Hits());
00498 StPmdHit *spmcl2;
00499 for(Int_t im=0; im<nmh2; im++)
00500 {
00501 spmcl2 = (StPmdHit*)next();
00502 if(spmcl2)
00503 {
00504 Int_t gsuper=spmcl2->Gsuper();
00505 Int_t col=spmcl2->Column();
00506 Int_t row=spmcl2->Row();
00507 Float_t edep=spmcl2->Edep();
00508 Int_t adc=spmcl2->Adc();
00510 StPhmdHit *phit = new StPhmdHit();
00511 phit->setSuperModule(Int_t(gsuper-1));
00512 phit->setSubDetector(Int_t(subdet));
00513 phit->setRow(Int_t(row));
00514 phit->setColumn(Int_t(col));
00515 phit->setEnergy(edep);
00516 phit->setAdc(adc);
00517 evtdet1->addHit(phit);
00518 }
00519 }
00520 }
00521 }
00522
00523 return kStOK;
00524 }
00525
00526 void StPmdSimulatorMaker::FinalEdep(StPmdDetector* pdet,Int_t id)
00527 {
00528
00529 StPmdModule * mod=pdet->module(id);
00530 Int_t nmh=pdet->module_hit(id);
00531 if(nmh>0)
00532 {
00533 TIter next(mod->Hits());
00534 StPmdHit *spmcl;
00535 for(Int_t im=0; im<nmh; im++)
00536 {
00537 Float_t rawadc=0.;
00538
00539 Int_t ADC=0;
00540 spmcl = (StPmdHit*)next();
00541 if(spmcl)
00542 {
00543 Float_t rawedep=spmcl->Edep();
00544 Float_t keVedep=rawedep;
00545 keV_ADC(keVedep,rawadc);
00546 if(mResFlag) ADC_Readout(rawadc,ADC);
00547 else ADC = (int)rawadc;
00548
00549
00550
00551
00552 spmcl->setAdc(ADC);
00553 }
00554 }
00555 }
00556 }
00557
00558 Float_t StPmdSimulatorMaker::keV_ADC(Float_t edep, Float_t& adc)
00559 {
00560 adc=mlcon0 + mlcon1*edep + mlcon2*pow(edep,2);
00561 return kStOK;
00562 }
00563
00564 Float_t StPmdSimulatorMaker::ADC_Readout(Float_t adc,Int_t& ADC)
00565 {
00566 Float_t reso_percent=0., reso=0.;
00567 reso_percent=mpcon0 + mpcon1*adc + mpcon2*pow(adc,2);
00568 reso=(reso_percent*100.)/adc;
00569
00570 Float_t adcprime=gRandom->Gaus(adc,reso);
00571 if(adcprime<0)adcprime=0;
00572 ADC=Int_t(adcprime);
00573 return kStOK;
00574 }
00575
00576
00577
00578
00579
00580
00581
