StRoot: StSvtSimulationMaker/StSvtEmbeddingMaker.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StSvtEmbeddingMaker.cxx,v 1.15 2009/01/22 23:19:21 fine Exp $
00004  *
00005  * Author: Selemon Bekele
00006  ***************************************************************************
00007  *
00008  * Description: Svt Embedding Maker class
00009  *
00010  ***************************************************************************
00011  *
00012  * $Log: StSvtEmbeddingMaker.cxx,v $
00013  * Revision 1.15  2009/01/22 23:19:21  fine
00014  * Prptection against fo the crash
00015  *
00016  * Revision 1.14  2009/01/22 22:45:45  fine
00017  * Add the fatal message to stop the chain
00018  *
00019  * Revision 1.13  2007/12/24 17:38:59  fisyak
00020  * spelling error StSvtRmsPedestal => StSvtRMSPedestal
00021  *
00022  * Revision 1.12  2007/07/12 20:18:18  fisyak
00023  * read Db by deman
00024  *
00025  * Revision 1.11  2007/04/28 17:57:11  perev
00026  * Redundant StChain.h removed
00027  *
00028  * Revision 1.10  2005/07/11 19:20:56  caines
00029  * Add in shift due to pasa response that is accounted for in the real data t0 calc.
00030  *
00031  * Revision 1.9  2005/02/09 14:33:35  caines
00032  * New electron expansion routine
00033  *
00034  * Revision 1.5  2004/02/24 15:53:21  caines
00035  * Read all params from database
00036  *
00037  * Revision 1.4  2004/01/22 16:30:47  caines
00038  * Getting closer to a final simulation
00039  *
00040  * Revision 1.3  2003/11/30 20:51:46  caines
00041  * New version of embedding maker and make OnlSeqAdj a stand alone maker
00042  *
00043  * Revision 1.2  2003/09/07 03:49:06  perev
00044  * gcc 3.2 + WarnOff
00045  *
00046  * Revision 1.1  2003/07/31 19:18:09  caines
00047  * Petrs improved simulation code
00048  *
00049  **************************************************************************/
00050 
00051 #include "Stiostream.h"
00052 #include <cmath>
00053 using namespace std;
00054 
00055 #include "TH2.h"
00056 #include "TFile.h"
00057 
00058 #include "St_DataSetIter.h"
00059 #include "St_ObjectSet.h"
00060 #include "StSequence.hh"
00061 #include "StSvtEmbeddingMaker.h"
00062 #include "StMessMgr.h"
00063 
00064 #include "StSvtClassLibrary/StSvtHybridCollection.hh"
00065 #include "StSvtClassLibrary/StSvtHybridData.hh"
00066 #include "StSvtClassLibrary/StSvtData.hh"
00067 #include "StSvtClassLibrary/StSvtHybridPixelsD.hh"
00068 #include "StSvtClassLibrary/StSvtHybridPixels.hh"
00069 #include "StSvtClassLibrary/StSvtConfig.hh"
00070 #include "StSvtClassLibrary/StSvtHybridPed.hh"
00071 #include "StDAQMaker/StSVTReader.h"
00072 
00073 ClassImp(StSvtEmbeddingMaker)
00079 #define cDefaultBckgRMS 1.8
00080 
00081 
00082 //____________________________________________________________________________
00083 StSvtEmbeddingMaker::StSvtEmbeddingMaker(const char *name):StMaker(name)
00084 {
00086   mDoEmbedding=kTRUE;                  
00087   mPlainSimIfNoSVT=kFALSE;            
00088   setBackGround(kTRUE,cDefaultBckgRMS);//sets to default true and sigma 1.8
00089   SetPedRmsPreferences(kTRUE, kTRUE);  //read database
00090   
00091   
00092   mSimPixelColl = NULL;
00093   mRealDataColl = NULL;
00094   mPedColl = NULL;
00095   mPedRMSColl = NULL;
00096 }
00097 
00098 //____________________________________________________________________________
00099 StSvtEmbeddingMaker::~StSvtEmbeddingMaker()
00100 {
00101 
00102 }
00103 
00104 //____________________________________________________________________________
00105 Int_t StSvtEmbeddingMaker::Init()
00106 {
00107   return StMaker::Init();
00108 }
00109 
00110 //____________________________________________________________________________
00112 Int_t StSvtEmbeddingMaker::InitRun(int runumber)
00113 {
00114   if (mUsePixelRMS){
00115     gMessMgr->Warning()<<"StSvtEmbeddingMaker - reading individual pixel RMS values from database "<<endm;
00116     mPedRMSColl= (StSvtHybridCollection*) ((TObjectSet *) GetDataSet("StSvtRMSPedestal"))->GetObject();
00117     if (mPedRMSColl) {gMessMgr->Warning()<<"StSvtEmbeddingMaker: Found RMS values for individual pixels."<<endm;}
00118     else {gMessMgr->Warning()<<"StSvtEmbeddingMaker: NO RMS values for individual pixels."<<endm;}
00119   }
00120 
00121   if (mUseHybridRMS){
00122     gMessMgr->Warning()<<"StSvtEmbeddingMaker - reading individual hybrid RMS values from database "<<endm;
00123     mPedColl= (StSvtHybridCollection*) ((TObjectSet *) GetDataSet("StSvtRMSPedestal"))->GetObject();
00124     if (mPedColl) {gMessMgr->Warning()<<"StSvtEmbeddingMaker: Found RMS values for individual hybrids."<<endm;}
00125     else {gMessMgr->Warning()<<"StSvtEmbeddingMaker: NO RMS values for individual hybrids."<<endm;}
00126   }
00127   if ((!mPedRMSColl)&&(!mPedColl))
00128     {gMessMgr->Warning()<<"Warning: no SVT RMS information available from Chain - using default backgroung:"<<mBackGSigma<<endm;}
00129 
00130   return StMaker::InitRun(runumber);
00131 }
00132 
00133 //____________________________________________________________________________
00134 Int_t StSvtEmbeddingMaker::Make()
00135 {
00136   //now decide in which mode to realy run, based on options
00137   //this could be done in InitRun, but could there be mixed event in one run with SVT and without SVT, ie.. some trigger mix
00138   mRunningEmbedding=mDoEmbedding;
00139 
00140   if (mDoEmbedding &&NoSvt() ){//check if SVT is present if not skip simulation
00141 
00142     if (mPlainSimIfNoSVT){//run plain simulation instead of embedding
00143       mRunningEmbedding=kFALSE; //run plain simulation
00144     }
00145     else 
00146       { //clear data and get out
00147         ClearOutputData();
00148         gMessMgr->Info()<<"SVT SlowSimulation: SKIPPING THIS EVENT - no SVT in real data!!"<<endm;
00149         return kStOk;
00150       }  
00151   }
00152   
00153   //write out the state of simulation
00154   char st[20];
00155   if (mRunningEmbedding) sprintf(st,"EMBEDDING");
00156   else sprintf(st,"PLAIN SIMULATION");
00157   gMessMgr->Info()<<"SVT SlowSimulation is running in the state of :"<<st<<endm;
00158   
00159   //  if (mRunningEmbedding)
00160   //{ 
00161   Int_t res;
00162   res=GetSvtData();
00163   if (res!=kStOk) return res;
00164     //}
00165   
00166   ClearMask(); //it has to be cleared here - needed for plain simulation
00167   for(int Barrel = 1;Barrel <= mSimPixelColl->getNumberOfBarrels();Barrel++) {
00168     for (int Ladder = 1;Ladder <= mSimPixelColl->getNumberOfLadders(Barrel);Ladder++) {
00169       for (int Wafer = 1;Wafer <= mSimPixelColl->getNumberOfWafers(Barrel);Wafer++) {
00170         for( int Hybrid = 1;Hybrid <= mSimPixelColl->getNumberOfHybrids();Hybrid++){
00171           mCurrentIndex = mSimPixelColl->getHybridIndex(Barrel,Ladder,Wafer,Hybrid);
00172           if( mCurrentIndex < 0) continue; 
00173           
00174           mCurrentPixelData  = (StSvtHybridPixelsD*)mSimPixelColl->at(mCurrentIndex);
00175           
00176           if(!mCurrentPixelData) { //no data from simulation Maker
00177             gMessMgr->Error()<<"Error  in StSvtEmbeddingMaker::Make(): Something is wrong, no data from simulator for hybrid index:"<<mCurrentIndex<<endm;
00178             mCurrentPixelData = new StSvtHybridPixelsD(Barrel, Ladder, Wafer, Hybrid);
00179             mSimPixelColl->put_at(mCurrentPixelData,mCurrentIndex);
00180           }
00181           
00182           if (mRunningEmbedding){  //do the embedding
00183             ClearMask(); //just in case..
00184             AddRawData();
00185           }
00186            
00187           if (mBackGrOption) CreateBackground();
00188                     
00189         }
00190       }
00191     }
00192   }
00193   
00194   return kStOK;
00195 }
00196 
00197 //____________________________________________________________________________
00198 Int_t StSvtEmbeddingMaker::GetSvtData()
00199 {
00200   //EmbeddingMaker requires some data(at least empty) from the SimulationMaker
00201   mSimPixelColl=NULL;
00202   mRealDataColl=NULL;
00203     
00204   St_DataSet* dataSet = GetDataSet("StSvtPixelData");
00205   if (dataSet==NULL){
00206     LOG_FATAL <<"BIG TROUBLE:No data from simulator to work with!!!!"<<endm;
00207     return kStErr;
00208   }
00209 
00210   mSimPixelColl= (StSvtData*)(dataSet->GetObject());
00211   if ( mSimPixelColl==NULL){
00212     gMessMgr->Error()<<"BIG TROUBLE:Data from simulator is empty!!!!"<<endm;
00213     return kStErr;
00214   }
00215   
00216   if (!mRunningEmbedding) return kStOk; //dont read real data if you don't need them
00217   dataSet = GetDataSet("StSvtRawData");
00218   if (dataSet) mRealDataColl= (StSvtData*)(dataSet->GetObject());
00219   else gMessMgr->Error()<<"No StSvtRawData in the chain"<<endm;
00220   if (!mRealDataColl)      //switching to plain simulation, because there is no raw data
00221     {
00222       gMessMgr->Error()<<"Note: StSvtEmbeddingMaker is set to do embbeding, but found no raw data- embedding into empty event!!!"<<endm;
00223       return kStErr;
00224     }
00225   
00226   return kStOk;
00227 }
00228  
00229 //____________________________________________________________________________
00231 void StSvtEmbeddingMaker::AddRawData()
00232 { 
00233   
00234   StSequence* Sequence;
00235   
00236   int numOfSeq;
00237   int* anolist;
00238   
00239   StSvtHybridData* realData = (StSvtHybridData *)mRealDataColl->at(mCurrentIndex); 
00240     
00241   double *adcArray=mCurrentPixelData->GetArray();      
00242   double offset =  mCurrentPixelData->getPedOffset(); //this could be problem if offset differs between real and simulated data!!!
00243   
00244   anolist = NULL;
00245   if (realData!=NULL)
00246     for (int iAnode= 0; iAnode<realData->getAnodeList(anolist); iAnode++){
00247       
00248       int Anode = anolist[iAnode]; //goes from 1
00249       realData->getSequences(Anode,numOfSeq,Sequence);
00250       
00251       for (int nSeq=0; nSeq< numOfSeq ; nSeq++){ 
00252         unsigned char* adc=Sequence[nSeq].firstAdc;
00253         int length = Sequence[nSeq].length;
00254         int startTimeBin=Sequence[nSeq].startTimeBin;
00255         
00256         for(int t = 0; t<length; t++){
00257           int time = startTimeBin + t;
00258           unsigned char adcVal = adc[t];
00259           //pixelIndex = mSimPixelData->getPixelIndex(Anode,time);
00260           //less clean but faster
00261           int pixelIndex=(Anode-1)*128+time;
00262           //there already is pedestal offset
00263           adcArray[pixelIndex]=adcArray[pixelIndex]+(double)adcVal-offset;
00264           
00265         //don't do noise here
00266           mMask[pixelIndex]=kFALSE;
00267         }
00268       }
00269     }
00270   
00271   //now clear rest of the mask
00272   for(int an = 0; an < 240; an++){
00273     for(int tim = 0; tim < 128; tim++){
00274       int pIndex=an*128 + tim;
00275       if (adcArray[pIndex]==offset) mMask[pIndex]=kFALSE; //don't make extra noise outside of hits
00276     }
00277   }
00278   
00279 }
00280 
00281 //____________________________________________________________________________
00282 double StSvtEmbeddingMaker::MakeGaussDev(double sigma)
00283 {
00284   
00285   double fac,rsq,v1,v2;
00286   
00287   do {
00288     v1 = 2.0*((float)rand()/(float)RAND_MAX) - 1.0;
00289     v2 = 2.0*((float)rand()/(float)RAND_MAX) - 1.0;
00290     rsq = v1*v1 + v2*v2;
00291   } while(rsq >= 1.0 );
00292   
00293   fac = sigma*::sqrt(-2.0*::log(rsq)/rsq);
00294   
00295   // gset = v1*fac;  // gset = 3.0*::sqrt(-2.0*::log(rsq))*(v1/::sqrt(rsq))
00296    return v2*fac;
00297  
00298 }
00299 
00300 //____________________________________________________________________________
00302 void StSvtEmbeddingMaker::CreateBackground()
00303 {
00304   const double rmsScale=16.;
00305   double *adcArray=mCurrentPixelData->GetArray(); // array of [128*240]
00306   
00307   //find out what background to use
00308   StSvtHybridPixels* pedRms=NULL;
00309   if (mPedRMSColl)
00310     { 
00311       pedRms = (StSvtHybridPixels*)mPedRMSColl->at(mCurrentIndex);
00312       if (pedRms == NULL) {gMessMgr->Warning()<<"Warning: Individual pixel RMS info is empty for hybrid "<<mCurrentIndex<<" =>have to use other method "<<endm;}
00313     }
00314   
00315   if(pedRms)
00316     {// I have rms for each pixel
00317       for(int an = 0; an < 240; an++)  for(int tim = 0; tim < 128; tim++){
00318         //cout<<pedRms<<"indiv rms="<<pedRms->At(pedRms->getPixelIndex(an+1,tim))/rmsScale<<endl;
00319         int pIndex=an*128 + tim;
00320         if (mMask[pIndex]) adcArray[pIndex]+=MakeGaussDev(pedRms->At(pedRms->getPixelIndex(an+1,tim))/rmsScale);// !! mAdcArray already contains PedOffset
00321       }
00322     }
00323   else {
00324     //one value for hybrid
00325     double backgsigma;
00326     StSvtHybridPed *ped=NULL;
00327     if (mPedColl){ 
00328       ped=(StSvtHybridPed *)mPedColl->at(mCurrentIndex);
00329       if (ped == NULL) {gMessMgr->Warning()<<"Warning: hybrid  RMS info is empty for hybrid "<<mCurrentIndex<<" =>using default value "<<mBackGSigma<<endm;}
00330         }
00331     if (ped) backgsigma=ped->getRMS(); else  backgsigma=mBackGSigma; //the default value
00332     if ((backgsigma<=0.)||(backgsigma>=6.)){ //check for obviously bad values 
00333       {gMessMgr->Warning()<<"Warnig for index "<<mCurrentIndex<<" pedestal RMS is:"<<backgsigma<<" => seting to default "<<mBackGSigma<<endm;}
00334       backgsigma=mBackGSigma;
00335     }
00336       if (Debug()) {gMessMgr->Debug()<<"for index "<<mCurrentIndex<<" pedestal RMS is:"<< backgsigma<<endm;}
00337     
00338     for(int an = 0; an < 240; an++){
00339       for(int tim = 0; tim < 128; tim++){
00340         int pIndex=an*128 + tim;
00341         if (mMask[pIndex]) adcArray[pIndex]+=MakeGaussDev(backgsigma);// !! mAdcArray already contains PedOffset             
00342       }
00343     }
00344   }
00345 }
00346 
00347 //____________________________________________________________________________
00348 Int_t StSvtEmbeddingMaker::Finish()
00349 { 
00350   return kStOK;
00351 }
00352 
00353 
00354 //____________________________________________________________________________
00355 void StSvtEmbeddingMaker::SetPedRmsPreferences(Bool_t usePixelRMS, Bool_t useHybridRMS)
00356 { // allows to disable reading RMS from database
00357   mUsePixelRMS=usePixelRMS;
00358   mUseHybridRMS=useHybridRMS;
00359 }
00360 
00361 
00362 //____________________________________________________________________________
00363 void  StSvtEmbeddingMaker::setDoEmbedding(Bool_t doIt){
00364   mDoEmbedding = doIt;
00365 }
00366 
00367 //____________________________________________________________________________
00368 void StSvtEmbeddingMaker::setPlainSimEvenIfNoSVT(Bool_t doIt){
00369   mPlainSimIfNoSVT= doIt;
00370 }
00371 //____________________________________________________________________________
00372 void  StSvtEmbeddingMaker::setBackGround(Bool_t backgr,double backgSigma){
00373   mBackGrOption = backgr;
00374   mBackGSigma = backgSigma;
00375 }
00376 
00377 //____________________________________________________________________________
00378 void  StSvtEmbeddingMaker::ClearMask()
00379 {
00380   for (int  i=0;i<128*240;i++)mMask[i]=kTRUE;
00381 }
00382 
00383 //_____________________________________________________________
00384 Int_t StSvtEmbeddingMaker::NoSvt()
00385 {  
00386   St_DataSet *dataSet;
00387 
00388   dataSet = GetDataSet("StDAQReader");
00389   if(!dataSet){
00390     gMessMgr->Error()<<("BIG TROUBLE: cannot find StDAQReader in the chain and you want to run embedding")<<endm;
00391     return kTRUE;
00392   }
00393   StDAQReader* daqReader = (StDAQReader*)(dataSet->GetObject());
00394   if (!daqReader){
00395     gMessMgr->Error()<<("BIG TROUBLE: StDAQReader is empty and you want to run embedding")<<endm;
00396     return kTRUE;
00397   }
00398   
00399   if (!daqReader->SVTPresent ())
00400     {
00401       gMessMgr->Info()<<("NO SVT in DAQ")<<endm;
00402       return kTRUE; //No SVT in the datastream 
00403     }
00404   
00405   gMessMgr->Info()<<("SVT found in DAQ")<<endm;
00406   return kFALSE;
00407 }
00408 
00409 //_____________________________________________________________
00410 //this resets mSvtSimPixelColl
00411 void StSvtEmbeddingMaker::ClearOutputData()
00412 {
00413   StSvtHybridPixelsD* tmpPixels;
00414  if (mSimPixelColl) {
00415   for(int Barrel = 1;Barrel <= mSimPixelColl->getNumberOfBarrels();Barrel++) {
00416      for (int Ladder = 1;Ladder <= mSimPixelColl->getNumberOfLadders(Barrel);Ladder++) {
00417        for (int Wafer = 1;Wafer <= mSimPixelColl->getNumberOfWafers(Barrel);Wafer++) {
00418          for( int Hybrid = 1;Hybrid <= mSimPixelColl->getNumberOfHybrids();Hybrid++){
00419            
00420            int index = mSimPixelColl->getHybridIndex(Barrel,Ladder,Wafer,Hybrid);
00421            if( index < 0) continue; 
00422            
00423            tmpPixels  = (StSvtHybridPixelsD*)mSimPixelColl->at(index);
00424           
00425            if(!tmpPixels) {
00426              tmpPixels = new StSvtHybridPixelsD(Barrel, Ladder, Wafer, Hybrid);
00427              mSimPixelColl->put_at(tmpPixels,index);
00428            }
00429            tmpPixels->setPedOffset(0);
00430            tmpPixels->reset();
00431          }
00432        }
00433      }
00434    }
00435   } else {
00436      LOG_ERROR << "StSvtEmbeddingMaker::ClearOutputData(): There is no SVT data to clear" << endm;
00437   }
00438 }