StRoot: StTofSimMaker/StTofSimMaker.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StTofSimMaker.cxx,v 1.11 2007/04/17 23:02:36 dongx Exp $
00004  *
00005  * Author: Frank Geurts
00006  ***************************************************************************
00007  *
00008  * Description: StTofSimMaker class for TOFp Simulations
00009  *
00010  ***************************************************************************
00011  *
00012  * $Log: StTofSimMaker.cxx,v $
00013  * Revision 1.11  2007/04/17 23:02:36  dongx
00014  * replaced with standard STAR Loggers
00015  *
00016  * Revision 1.10  2006/12/08 18:55:26  dongx
00017  * Update to avoid zero tdc value in denominator - modified by Jing Liu
00018  *
00019  * Revision 1.9  2005/04/13 16:03:29  dongx
00020  * corresponding changes because of the update of StTofData structure
00021  *
00022  * Revision 1.8  2004/04/01 21:33:46  jeromel
00023  * More than one place where m_Mode should be used
00024  *
00025  * Revision 1.6  2003/09/17 19:49:10  geurts
00026  * zeroed pointers in constructor
00027  *
00028  * Revision 1.5  2003/09/02 17:59:10  perev
00029  * gcc 3.2 updates + WarnOff
00030  *
00031  * Revision 1.4  2003/08/08 00:22:11  geurts
00032  * changed location of header files for the local collections
00033  *
00034  * Revision 1.3  2003/07/25 04:34:44  geurts
00035  * - upper adc and tdc limits
00036  * - geometry initialization moved to InitRun()
00037  *
00038  * Revision 1.2  2002/12/12 01:43:46  geurts
00039  * Introduced InitRun() and FinishRun() members.
00040  * TofData in TofCollection is filled with adc and tdc data.
00041  * Extra checks for StEvent object to prevent null pointers.
00042  * Primitive ADC response function, disabled slatResponseExp().
00043  *
00044  * Revision 1.1  2001/09/28 19:11:11  llope
00045  * first version
00046  *
00047  **************************************************************************/
00049 
00056 #include <Stiostream.h>
00057 #include "StTofSimMaker.h"
00058 
00059 // SCL
00060 #include "StThreeVectorD.hh"
00061 #include "Random.h"
00062 #include "RanluxEngine.h"
00063 #include "RandGauss.h"
00064 #include "TH1.h"
00065 #include "TFile.h"
00066 
00067 #include "StTofUtil/StTofCalibration.h"
00068 #include "StTofUtil/StTofSimParam.h"
00069 #include "StTofUtil/StTofGeometry.h"
00070 #include "StEventTypes.h"
00071 
00072 // g2t tables and collections
00073 #include "tables/St_g2t_ctf_hit_Table.h"
00074 #include "tables/St_g2t_vpd_hit_Table.h"
00075 //#include "tables/St_g2t_track_Table.h"
00076 //#include "tables/St_g2t_tpc_hit_Table.h"
00077 #include "StTofUtil/StTofDataCollection.h"
00078 #include "StTofUtil/StTofSlatCollection.h"
00079 #include "StTofMCSlat.h"
00080 
00081 typedef vector<StTofMCSlat> tofMCSlatVector;
00082 typedef tofMCSlatVector::iterator tofMCSlatVecIter;
00083 
00084 
00085 static RanluxEngine engine;
00086 
00087 ClassImp(StTofSimMaker)
00088 
00089 
00090 StTofSimMaker::StTofSimMaker(const char *name):StMaker(name){
00091   mGeomDb  = 0;
00092   mCalibDb = 0;
00093   mSimDb   = 0;
00094 }
00095 
00097 StTofSimMaker::~StTofSimMaker(){/* nope */}
00098 
00099 
00101 Int_t StTofSimMaker::Init(){
00102   //mGeomDb = new StTofGeometry();
00103   //mGeomDb->init(this);
00104   mCalibDb = new StTofCalibration();
00105   mCalibDb->init();
00106   mSimDb = new StTofSimParam();
00107   mSimDb->init();
00108   
00109   if (m_Mode){
00110     // book histograms
00111     mdE = new TH1F("energy deposition","dE",100,0.,0.011);
00112     mdS = new TH1F("distance","ds",100,0.,10);
00113     mNumberOfPhotoelectrons = new TH1F("number of photoelectrons","nphe",1000,0,5000);
00114     mT = new TH1F("delay corrected tof","time",100,0.,12e-7);
00115     mTime = new TH1F("only hit-pos resolution added","tt",100,0.,12e-7);
00116     mTime1 = new TH1F("fully corrected tof","tt1",100,0.,120e-7);
00117     mPMlength = new TH1F("distance in slat","length",100,0,22);
00118     mAdc = new TH1F("adc","adc",1025,-0.5,1024.5);
00119     mTdc = new TH1F("tdc","tdc",2049,-0.5,2048.5);
00120   }
00121 
00122   return StMaker::Init();
00123 }
00124 
00125 
00126 
00128 Int_t StTofSimMaker::InitRun(int runnumber){
00129   LOG_INFO << "StTofSimMaker::InitRun  -- initializing TofGeometry --" << endm;
00130   mGeomDb = new StTofGeometry();
00131   mGeomDb->init(this);
00132   return kStOK;
00133 }
00134 
00135 
00136 
00138 Int_t StTofSimMaker::FinishRun(int runnumber){
00139   LOG_INFO << "StTofSimMaker::FinishRun -- cleaning up TofGeometry --" << endm;
00140   if (mGeomDb) delete mGeomDb;
00141   mGeomDb=0;
00142   return 0;
00143 }
00144 
00145 
00146 
00148 Int_t StTofSimMaker::Make(){
00149   LOG_INFO << "StTofSimMaker  Make() starts" << endm;
00150 
00151   StTofSlatCollection *mSlatCollection =  new StTofSlatCollection;
00152   tofMCSlatVector tofMC;
00153 
00154   // Check for GEANT data
00155   St_DataSet *geantData = GetInputDS("geant");
00156   if (geantData) {
00157     St_DataSetIter geantIter(geantData);
00158 
00159     // TOFp hits
00160     St_g2t_ctf_hit *g2t_tof_hit = (St_g2t_ctf_hit*) geantIter("g2t_tof_hit");
00161     if (g2t_tof_hit){
00162       g2t_ctf_hit_st* tof_hit = g2t_tof_hit->GetTable();
00163       int numberOfTofHits = g2t_tof_hit->GetNRows();
00164       LOG_INFO << "TOF #hits: " << numberOfTofHits << endm;
00165       
00166       tofMCSlatVector MCSlatVec;
00167       MCSlatVec.clear();
00168       for (int i=0;i<numberOfTofHits;i++,tof_hit++){
00169         MCSlatVec.push_back(detectorResponse(tof_hit));
00170       }
00171 
00172       // build ... from MCSlatVec which may have entries with same slatId
00173       tofMCSlatVector slatTempVec = MCSlatVec;
00174       tofMCSlatVector slatErasedVec = slatTempVec;
00175       tofMCSlatVecIter slatTempIter, slatErasedIter;
00176         
00177       while (slatTempVec.size()!=0){
00178         unsigned short fastTdc;
00179         int nFired=0, accumNPhe=0;
00180         float accumDe=0., accumDs=0., fastTof=0.;
00181         slatTempIter=slatTempVec.begin();
00182         slatErasedIter=slatErasedVec.begin();
00183         fastTof = slatTempIter->mcInfo().mTof;
00184         fastTdc = slatTempIter->tdc();
00185             
00186         while(slatErasedIter!=slatErasedVec.end()){
00187           if(slatTempIter->slatIndex() == slatErasedIter->slatIndex()){
00188             nFired++;
00189             accumDe += slatErasedIter->mcInfo().mDe;
00190             accumDs += slatErasedIter->mcInfo().mDs;
00191             accumNPhe += slatErasedIter->mcInfo().mNPhe;
00192             fastTof = min(fastTof, slatErasedIter->mcInfo().mTof);
00193             fastTdc = min(fastTdc, slatErasedIter->tdc());
00194               
00195             slatErasedVec.erase(slatErasedIter);
00196             slatErasedIter--;
00197           }
00198           slatErasedIter++;
00199         }
00200         StTofMCSlat MCSlat = *slatTempIter;
00201         MCSlat.setNHits(nFired);
00202         MCSlat.setNPhe(accumNPhe);
00203         MCSlat.setDe(accumDe);
00204         MCSlat.setDs(accumDs);
00205         MCSlat.setTof(fastTof);
00206         MCSlat.setTdc(fastTdc);
00207         tofMC.push_back(MCSlat);
00208         
00209         slatTempVec = slatErasedVec;
00210       }
00211       LOG_INFO << "StTofSimMaker::make()  vector size from " << MCSlatVec.size()
00212            << " to " << tofMC.size() << endm;
00213         
00215         
00216       fillRaw();
00217       if (mSimDb->elec_noise() < 0)  electronicNoise();
00218       fillEvent();
00219       for (unsigned int i=0;i<tofMC.size(); i++){
00220         StTofMCSlat *MCSlatPtr = new StTofMCSlat();
00221         *MCSlatPtr = tofMC[i];
00222         //LOG_INFO << *MCSlatPtr ;
00223         mSlatCollection->push_back(MCSlatPtr);
00224       }
00225     }
00226     else
00227       LOG_INFO << "StTofSimMaker Make()  no TOF hits found" << endm;
00228 
00229     // pVPD section hits
00230     St_g2t_vpd_hit *g2t_vpd_hit = (St_g2t_vpd_hit*) geantIter("g2t_vpd_hit");
00231     if (g2t_vpd_hit){
00232       //   g2t_vpd_hit_st* vpd_hit = g2t_vpd_hit->GetTable();
00233       int numberOfVpdHits = g2t_vpd_hit->GetNRows();
00234       LOG_INFO << "VPD #hits: " << numberOfVpdHits << endm;
00235     }
00236     else
00237       LOG_INFO << "StTofSimMaker Make()  no VPD hits found" << endm;    
00238   }
00239 
00240 
00241   // send off to StEvent
00242   StTofCollection *mTheTofCollection = new StTofCollection();
00243   for (size_t j=0;j<mSlatCollection->size();j++){
00244     mTheTofCollection->addSlat(dynamic_cast<StTofMCSlat*>(mSlatCollection->getSlat(j))); 
00245     // mTheTofCollection->addSlat(mSlatCollection->getSlat(j)); 
00246 }
00247 
00248   // create tofData collection
00249   StTofDataCollection *mDataCollection = new StTofDataCollection;
00250   for (int i=0;i<48;i++){    
00251     bool slatFound = false;
00252     int j=0;
00253     for (j=0;j<(int)mSlatCollection->size();j++){
00254       StTofSlat *tempSlat = mSlatCollection->getSlat(j);
00255       unsigned short indexSlat = tempSlat->slatIndex();
00256       if (indexSlat == mGeomDb->daqToSlatId(i)){
00257         slatFound = true;
00258         //      StTofData *rawTofData = new  StTofData(indexSlat,tempSlat->adc(),tempSlat->tdc(),0,0);
00259         // update for year 5 new format
00260         StTofData *rawTofData = new  StTofData(indexSlat,tempSlat->adc(),tempSlat->tdc(),0,0,0,0);
00261         if (Debug()) LOG_INFO << indexSlat << ":  A" << tempSlat->adc() << "  T" << tempSlat->tdc() << endm;
00262         mDataCollection->push_back(rawTofData);
00263       }
00264     }
00265     if (!slatFound){
00266       //      StTofData *rawTofData = new  StTofData(mGeomDb->daqToSlatId(i),0,0,0,0);
00267       // update for year 5 new format
00268       StTofData *rawTofData = new  StTofData(mGeomDb->daqToSlatId(i),0,0,0,0,0,0);
00269       mDataCollection->push_back(rawTofData);
00270     }
00271   }
00272   for(size_t jj = 0; jj < mDataCollection->size(); jj++)
00273     mTheTofCollection->addData(mDataCollection->getData(jj));
00274 
00275   mEvent = (StEvent*) GetInputDS("StEvent");
00276   if (mEvent) mEvent->setTofCollection(mTheTofCollection);
00277   else{
00278     LOG_INFO << "StTofSimMaker: Where is StEvent !?! Unable to store data" << endm;
00279     return kStWarn;
00280   }
00281 
00282 
00283   // verify existence of tofCollection in StEvent (mEvent) 
00284   LOG_INFO << "StTofSimMaker: verifying TOF StEvent data ..." << endm;
00285   StTofCollection *mmTheTofCollection = mEvent->tofCollection();
00286   if(mmTheTofCollection) {
00287     LOG_INFO << " + StEvent tofCollection Exists" << endm;
00288     if(mmTheTofCollection->slatsPresent())
00289       LOG_INFO << " + StEvent TofSlatCollection Exists" << endm;
00290     else
00291       LOG_INFO << " - StEvent TofSlatCollection DOES NOT Exist" << endm;
00292     if(mmTheTofCollection->hitsPresent())
00293       LOG_INFO << " + StEvent TofHitCollection Exists" << endm;
00294     else
00295       LOG_INFO << " - StEvent TofHitCollection DOES NOT Exist" << endm;
00296   }
00297   else {
00298     LOG_INFO << " - StEvent tofCollection DOES NOT Exist" << endm;
00299     LOG_INFO << " - StEvent TofSlatCollection DOES NOT Exist" << endm;
00300     LOG_INFO << " - StEvent TofHitCollection DOES NOT Exist" << endm;
00301   }
00302 
00303   LOG_INFO << "StTofSimMaker  Make() finished" << endm;
00304   return kStOK;
00305 }
00306 
00307 
00309 StTofMCSlat StTofSimMaker::detectorResponse(g2t_ctf_hit_st* tof_hit)
00310 {
00311   if(Debug()){
00312     // dump the g2t structure ...
00313     LOG_INFO << " " <<setw( 3) << tof_hit->id      << " " <<setw( 4) << tof_hit->next_tr_hit_p
00314          << " " <<setw( 4) << tof_hit->track_p << " " <<setw( 8) << tof_hit->volume_id
00315          << " " <<setw(13) << tof_hit->de      << " " <<setw(11) << tof_hit->ds
00316          << " " <<setw(12) << tof_hit->p[0]    << " " <<setw(12) << tof_hit->p[1]
00317          << " " <<setw(12) << tof_hit->p[2]    << " " <<setw( 7) << tof_hit->s_track
00318          << " " <<setw(13) << tof_hit->tof     << " " <<setw(10) << tof_hit->x[0]
00319          << " " <<setw(10) << tof_hit->x[1]    << " " <<setw(10) << tof_hit->x[2]
00320          << endm;
00321   }
00322     // skip the consistency checks for now,
00323 
00324     // determine eta and phi indices from hitpoint and slatId
00325     StThreeVectorD hitPoint  = StThreeVectorD(tof_hit->x[0],
00326                                               tof_hit->x[1],
00327                                               tof_hit->x[2]);
00328     //fg    StTofCross*  mTofCross = new StTofCross;
00329     int slatId= (int) mGeomDb->tofSlatCrossId(hitPoint);
00330     int volId = (int) mGeomDb->tofSlatCrossId(tof_hit->volume_id);
00331 
00332     if (slatId != volId){
00333       LOG_INFO << "StTofSimMaker::Make  Warning: volume_id ("<< volId 
00334            << ") and hit ("<<slatId<<") inconsistent. Switching to volumeid."<< endm;
00335       slatId=volId;
00336     }
00337 
00338     // retrieve dbase parameters
00339     float zmin = mGeomDb->tofSlat(slatId).z_min;
00340     float zmax = mGeomDb->tofSlat(slatId).z_max;
00341     float cosang  = mGeomDb->tofSlat(slatId).cosang;
00342     
00343     float length = (zmax-tof_hit->x[2])/cosang ;
00344     float max_distance = (zmax-zmin)/cosang ;
00345         
00346     if (length>max_distance || length<0){
00347       LOG_INFO <<  "StTofSimMaker:  length="<<length<<" max="<<max_distance
00348            << " zmin="<<zmin<<" zmax="<<zmax<<" coasng="<<cosang<<endm;
00349       mGeomDb->printSlat(slatId);
00350     }
00351 
00352         
00353     // do the slat response modelling similar to cts
00354     long numberOfPhotoelectrons;
00355     if (mSimDb->slat_para()){
00356       LOG_INFO << "StTofSimMaker  Slat Response Table not implemented yet. "
00357       << " Switching to exponential model instead" <<endm;
00358     }
00359     //numberOfPhotoelectrons = long (tof_hit->de * slatResponseExp(length));
00360     numberOfPhotoelectrons = long (tof_hit->de * mSimDb->GeV_2_n_photons()
00361                                * mSimDb->cath_eff() * mSimDb->cath_surf() * mSimDb->surf_loss());
00362 
00363     // prepare some random generator stuff
00364     RandGauss random(engine);
00365 
00366     // calculate TOFs with all kinds of resolutions
00367     float time= tof_hit->tof + length*mSimDb->delay();
00368     float resl=  mSimDb->time_res() * ::sqrt(length);
00369     if (resl<50e-12) resl=50e-12;
00370     float tt =  tof_hit->tof + (float) random.shoot()* resl;
00371     float tt1 =  time +  (float) random.shoot()* mSimDb->start_res();
00372 
00373     // fill the histograms
00374     if(m_Mode){
00375       mPMlength->Fill(length);
00376       mdE->Fill(tof_hit->de);
00377       mdS->Fill(tof_hit->ds);
00378       mNumberOfPhotoelectrons->Fill(numberOfPhotoelectrons);
00379       mT->Fill(time);
00380       mTime->Fill(tt);
00381       mTime1->Fill(tt1);
00382     }
00383 
00384     // fill or update the mcInfo structure of StTofMCslat
00385     StTofMCSlat slat;
00386     slat.setSlatIndex(slatId);
00387     StTofMCInfo slatData;
00388     slatData.mNHits = 1; // slats will be reorderd and #hits/slat recounted
00389     slatData.mNPhe  = numberOfPhotoelectrons;
00390     slatData.mDe    = tof_hit->de;
00391     slatData.mDs    = tof_hit->ds;
00392     slatData.mTof =  tof_hit->tof;// warning: this is *NOT* correct ... prop-delay is not included !!!
00393     slatData.mTime = time;
00394     slatData.mMTime = tt;
00395     slatData.mMTimeL = tt1;
00396     slatData.mPmLength = length;
00397     slatData.mSLength = tof_hit->s_track;
00398     StThreeVectorD hitMomentum  = StThreeVectorD(tof_hit->p[0],
00399                                                  tof_hit->p[1],
00400                                                  tof_hit->p[2]);
00401     slatData.mPTot = hitMomentum.mag();
00402     slatData.mGId = tof_hit->id;
00403     slatData.mTrkId = tof_hit->track_p;
00404     
00405     slat.setMCInfo(slatData);
00406 
00407     // do a rough cts-style calibration
00408     float tdcOffset = mCalibDb->slat(slatId).offset_tdc
00409                       + random.shoot() * mCalibDb->slat(slatId).ods_tdc;
00410     float cctdc = mCalibDb->slat(slatId).cc_tdc;
00411     if(cctdc==0) cctdc=0.0001;
00412     unsigned short tdc = (unsigned short)((float)tt/cctdc+ tdcOffset);
00413     float adcOffset = mCalibDb->slat(slatId).offset_adc
00414                       + random.shoot() * mCalibDb->slat(slatId).ods_adc;
00415     unsigned short adc =  (unsigned short)((float)numberOfPhotoelectrons
00416                                      * mSimDb->nphe_to_adc() + adcOffset);
00417     if (tdc>2048) tdc=2048;
00418     if (adc>1024) adc=1024;
00419     slat.setAdc(adc);
00420     slat.setTdc(tdc);
00421 
00422     if(m_Mode){
00423       mAdc->Fill(adc);
00424       mTdc->Fill(tdc);
00425     }
00426 
00427     if (Debug()){
00428       LOG_INFO << "StTofmcInfo slatId " << slatId << "  " << slatData;
00429       LOG_INFO << "    a:" << adc << " t:" << tdc << " dE:"<< slatData.mDe << " dS:"<<  slatData.mDs 
00430            << " mTof:" <<  slatData.mTof << " mTime:"<<  slatData.mTime << " mMTime:"<<slatData.mMTime
00431            << " mMTimeL:"<< slatData.mMTimeL << " mSLength:" << slatData.mSLength
00432            << " mPTot:" << slatData.mPTot  << endm;
00433     }
00434 
00435     // this part considers X-talk between slats (based on parameter below) ...
00436     //    LOG_INFO << "PHYSNOISE PARAMETER: " << mSimDb->phys_noise() << endm;
00437   
00438     return slat;
00439 }
00440 
00441 
00442 
00444 float StTofSimMaker::slatResponseExp(float& dz)
00445 {
00446   // Exponential model for slat response
00447   return mSimDb->GeV_2_n_photons() * mSimDb->cath_eff()
00448          * mSimDb->cath_surf() * mSimDb->surf_loss()
00449          * exp(-dz/mSimDb->attlen());
00450 }
00451 
00452 
00454 Int_t StTofSimMaker::Finish(){
00455 
00456   if(m_Mode){
00457     LOG_INFO << "StTofSimMaker::Finish  writing tofsim.root ...";
00458     TFile theFile("tofsim.root","RECREATE","tofsim");
00459     theFile.cd();
00460     mdE->Write();
00461     mdS->Write();
00462     mNumberOfPhotoelectrons->Write();
00463     mT->Write();
00464     mTime->Write();
00465     mTime1->Write();
00466     mPMlength->Write();
00467     mAdc->Write();
00468     mTdc->Write();
00469     LOG_INFO << "done"<<endm;
00470   }
00471   return kStOK;
00472 }
00473 
00474 
00476 void StTofSimMaker::fillRaw(){
00477   //fill the adc and tdc entries.
00478 }
00479 
00480 
00482 void StTofSimMaker::electronicNoise(){
00483 }
00484 
00485 
00487  void StTofSimMaker::fillEvent(){
00488 }