StRoot: StSsdSimulationMaker/StSsdEmbeddingMaker.cxx Source File

00001  /***************************************************************************
00002  *
00003  * Author: Bouchet Jonathan
00004  ***************************************************************************
00005  *
00006  * Description: SsdEmbeddingMaker class
00007  *
00008  **************************************************************************/
00009 
00010 #include "Stiostream.h"
00011 #include <cmath>
00012 using namespace std;
00013 
00014 #include "TH2.h"
00015 #include "TFile.h"
00016 #include "TNtuple.h"
00017 
00018 #include "St_DataSetIter.h"
00019 #include "St_ObjectSet.h"
00020 #include "StSsdEmbeddingMaker.h"
00021 #include "StMessMgr.h"
00022 
00023 #include "tables/St_sls_strip_Table.h"
00024 #include "tables/St_spa_strip_Table.h"
00025 #include "StSsdDbMaker/StSsdDbMaker.h"
00026 #include "StSsdUtil/StSsdBarrel.hh"
00027 #include "StDAQMaker/StDAQReader.h"
00028 #include "StDAQMaker/StSSDReader.h"
00029 #include "StSsdSimulationMaker/St_spa_Maker.h"
00030 #include "StSsdPointMaker/StSsdPointMaker.h"
00031 
00032 ClassImp(StSsdEmbeddingMaker)
00033 
00034 //____________________________________________________________________________
00035 StSsdEmbeddingMaker::StSsdEmbeddingMaker(const char *name):StMaker(name)
00036 {
00038   mDoEmbedding     = kTRUE;                  
00039   mPlainSimIfNoSSD = kFALSE;    
00041   mSsdrealData = 0;
00042   mSsdSimuData = 0;
00043  }
00044 //____________________________________________________________________________
00045 StSsdEmbeddingMaker::~StSsdEmbeddingMaker()
00046 {
00047 }
00048 //____________________________________________________________________________
00049 Int_t StSsdEmbeddingMaker::Init()
00050 {
00051   if(IAttr(".histos")){
00052     hStripsSimu   = new TH1F("hStripsSimu","Log_{10} of the number of strip after dAQ cut in simulation",100,0,10);
00053     hStripsReal   = new TH1F("hStripsReal","Log_{10} of the number of strip after dAQ cut in real data",100,0,10);
00054     hStripsCommon = new TH1F("hStripsCommon","number of strips in common from real file and simu file",1000,0,1000);
00055   }
00056   return StMaker::Init();
00057 }
00058 //____________________________________________________________________________
00060 Int_t StSsdEmbeddingMaker::InitRun(int RunNo)
00061 {
00062   assert(StSsdBarrel::Instance());
00063   TDataSet *ssdparams = GetInputDB("Geometry/ssd");
00064   if (! ssdparams) {
00065     LOG_ERROR << "No  access to Geometry/ssd parameters" << endm;
00066     return kStFatal;
00067   }
00068   TDataSetIter    local(ssdparams);
00069   m_dimensions  = (St_ssdDimensions     *)local("ssdDimensions");
00070   ssdDimensions_st     *dimensions = m_dimensions->GetTable(); 
00071   setSsdParameters(dimensions);
00072   return kStOk;
00073 }
00074 //____________________________________________________________________________
00075 Int_t StSsdEmbeddingMaker::Make()
00076 {
00077   LOG_INFO<<"##################################"<<endm;
00078   LOG_INFO<<"####  START OF SSD EMBEDDING  ####"<<endm;
00079   //now decide in which mode to realy run, based on options
00080   //this could be done in InitRun, but could there be mixed event in one run with SSD and without SSD, ie.. some trigger mix
00081   mRunningEmbedding=mDoEmbedding;
00082   
00083   if (mDoEmbedding && NoSsd() ){//check if SSD is present if not skip simulation
00084 
00085     if (mPlainSimIfNoSSD){//run plain simulation instead of embedding
00086       mRunningEmbedding=kFALSE; //run plain simulation
00087     }
00088     else 
00089       { //clear data and get out
00090         LOG_INFO <<"SSD SlowSimulation: SKIPPING THIS EVENT - no SSD in real data!!"<<endm;
00091         return kStOk;
00092       }  
00093   }
00094   Int_t res;
00095   res = GetSsdData();
00096   if (res!=kStOk) return res;
00097   //write out the state of simulation
00098   char st[20];
00099   if (mRunningEmbedding) sprintf(st,"EMBEDDING");
00100   else sprintf(st,"PLAIN SIMULATION");
00101   LOG_INFO<<"SSD SlowSimulation is running in the state of :"<<st<<endm;
00102   if (mRunningEmbedding){  //do the embedding
00103     MODE = 1;
00104     Int_t SIZEREAL = mSsdrealData->GetNRows();
00105     Int_t SIZESIMU = mSsdSimuData->GetNRows();
00106     LOG_INFO <<"####    NUMBER OF STRIPS (SIMU) : "<<SIZESIMU <<" ####"<<endm;
00107     LOG_INFO <<"####    NUMBER OF STRIPS (REAL) : "<<SIZEREAL <<" ####"<<endm;
00108     if(Debug())CheckTables();
00109     Int_t numberofstrips_in_common = AddRawData();
00110     LOG_INFO <<"####  MERGE SIMU+REAL : "<<numberofstrips_in_common << "####" << endm;
00111     if (IAttr(".histos")){
00112       if(SIZESIMU>0){hStripsSimu->Fill(TMath::Log10(SIZESIMU));}
00113       if(SIZEREAL>0){hStripsReal->Fill(TMath::Log10(SIZEREAL));}
00114       hStripsCommon->Fill(numberofstrips_in_common);
00115     }
00116   }     
00117   return kStOk;
00118 }
00119 //____________________________________________________________________________
00120 Int_t StSsdEmbeddingMaker::GetSsdData()
00121 {
00122   //EmbeddingMaker requires some data(at least empty) from the SimulationMaker
00123   if(!mRunningEmbedding) return kStOk;//dont read real data if you don't need them
00124   mSsdSimuData = (St_spa_strip*)GetDataSet("SsdSimuData");
00125   //mSsdSimuData = (St_spa_strip*)GetDataSet("spa_strip/.data/spa_strip");
00126   if (!mSsdSimuData || mSsdSimuData->GetNRows()==0){ 
00127     LOG_ERROR << "GetSsdData : no simu input (fired strip for the SSD)"<<endm;   
00128     return kStErr;
00129   }
00130   if(mSsdSimuData) LOG_INFO << "Numbers of strips (simu) = " << mSsdSimuData->GetNRows()<< endm;
00131   
00132   mSsdrealData = (St_spa_strip*)GetDataSet("SsdRealData");
00133   if (!mSsdrealData || mSsdrealData->GetNRows()==0){ 
00134     LOG_ERROR << "GetSsdData : no real input (fired strip for the SSD)"<<endm;   
00135     return kStErr;
00136   }
00137   if(mSsdrealData) LOG_INFO << "Numbers of strips (real) = " << mSsdrealData->GetNRows()<< endm;
00138   return kStOk;
00139 }
00140 //____________________________________________________
00141 Int_t StSsdEmbeddingMaker::AddRawData(){
00142   mSsdSimuReal = new St_spa_strip("spa_strip",50000);   
00143   m_DataSet->Add(mSsdSimuReal);
00144   spa_strip_st out_strip ;
00145   Int_t idWafSimu         = 0;
00146   Int_t iWafSimu          = 0;
00147   Int_t iLadSimu          = 0;
00148   Int_t nStripSimu        = 0;
00149   Int_t iSideSimu         = 0;
00150   spa_strip_st *strip_simu = mSsdSimuData->GetTable();
00151   Int_t idWafReal         = 0;
00152   Int_t iWafReal          = 0;
00153   Int_t iLadReal          = 0;
00154   Int_t nStripReal        = 0;
00155   Int_t iSideReal         = 0;
00156   spa_strip_st *strip_real = mSsdrealData->GetTable();
00157   Int_t currRecord        = 0;
00158   Int_t currRecordSimu    = 0;
00159   Int_t currRecordReal    = 0;
00160   Int_t currRecordCommon  = 0;
00161   Int_t IdCommon[mSsdSimuData->GetNRows()];
00162   for(Int_t ll=0;ll<mSsdSimuData->GetNRows();ll++) IdCommon[ll]=0;  
00163   for (Int_t i = 0 ; i < mSsdrealData->GetNRows(); i++)
00164     {
00165       nStripReal  = (int)(strip_real[i].id_strip/100000.);
00166       idWafReal   = strip_real[i].id_strip-10000*((int)(strip_real[i].id_strip/10000.));
00167       iWafReal    = idWaferToWafer(idWafReal);
00168       iLadReal    = (int)(idWafReal - mSsdLayer*1000 - (iWafReal+1)*100 - 1);
00169       iSideReal   = (strip_real[i].id_strip - nStripReal*100000 - idWafReal)/10000;
00170       Int_t in =0 ;
00171       for (Int_t j = 0 ; j < mSsdSimuData->GetNRows(); j++)
00172         {
00173           nStripSimu  = (int)(strip_simu[j].id_strip/100000.);
00174           idWafSimu   = strip_simu[j].id_strip-10000*((int)(strip_simu[j].id_strip/10000.));
00175           iWafSimu    = idWaferToWafer(idWafSimu);
00176           iLadSimu    = (int)(idWafSimu - mSsdLayer*1000 - (iWafSimu+1)*100 - 1);
00177           iSideSimu   = (strip_simu[j].id_strip - nStripSimu*100000 - idWafSimu)/10000;
00178           
00179           if((nStripReal == nStripSimu) && (idWafReal == idWafSimu) && (iSideReal == iSideSimu)) // we found the same strip in simu and real
00180             {
00181               LOG_DEBUG<<Form("simu side=%d idWafer=%d Ladder=%d wafer=%d nstrip=%d signal=%d",iSideSimu,idWafSimu,iLadSimu,iWafSimu,nStripSimu,strip_simu[j].adc_count)<<endm;
00182               LOG_DEBUG<<Form("real side=%d idWafer=%d Ladder=%d wafer=%d nstrip=%d signal=%d",iSideReal,idWafReal,iLadReal,iWafReal,nStripReal,strip_real[i].adc_count)<<endm;
00183               out_strip.id          = currRecord + 1;
00184               out_strip.adc_count   = (int)(strip_simu[j].adc_count + strip_real[i].adc_count);
00185               out_strip.id_strip    = 10000*(10*nStripSimu + iSideSimu)+idWafSimu;
00186               for (Int_t e = 0 ; e < 5 ; e++)
00187                 {
00188                   out_strip.id_mchit[e] = strip_simu[j].id_mchit[e];
00189                   out_strip.id_mctrack[e] = strip_simu[j].id_mctrack[e];
00190                 }
00191               mSsdSimuReal->AddAt(&out_strip);
00192               IdCommon[currRecordCommon] = (int)strip_real[i].id_strip; // we keep this id 
00193               currRecord++;
00194               currRecordCommon++;
00195               in = 1;
00196               LOG_DEBUG<<Form("# of strips in common =%d id=%d",currRecordCommon,(int)strip_real[i].id_strip)<<endm;
00197               if(Debug())
00198                 {
00199                   for(Int_t ii=0;ii<5;ii++){
00200                     LOG_INFO <<"mchit["<<ii<<"]="<<strip_simu[j].id_mchit[ii]<<endm;
00201                     LOG_INFO <<"mctrack["<<ii<<"]="<<strip_simu[j].id_mctrack[ii]<<endm;
00202                     LOG_INFO << " "<< endm;
00203                   }
00204                 }
00205             }
00206           break; // break the loop
00207         }
00208       if(in==0) //no common strips, now we fill the real strip 
00209         {
00210           out_strip.id          = currRecord + 1;
00211           out_strip.adc_count   = (int)(strip_real[i].adc_count);
00212           out_strip.id_strip    = 10000*(10*nStripReal + iSideReal)+idWafReal;
00213           for (Int_t e = 0 ; e < 5 ; e++)
00214             {
00215               out_strip.id_mchit[e]   = 0;
00216               out_strip.id_mctrack[e] = 0;
00217             }
00218           mSsdSimuReal->AddAt(&out_strip);
00219           currRecordReal++;
00220           currRecord++;
00221         } 
00222     }
00223   // we continue to fill the table with simu strips not found (in common with the real)
00224   Int_t found ;
00225   for(Int_t k=0; k <mSsdSimuData->GetNRows();k++)
00226     {
00227       found = 0;
00228       for (Int_t kk=0;kk<currRecordCommon;kk++)
00229         {               
00230           if((int)strip_simu[k].id_strip == IdCommon[kk]){
00231             LOG_DEBUG<<Form("k = %d strip id =%d",k,(int)strip_simu[k].id_strip)<<endm;
00232             found = 1;
00233           }
00234         }
00235       if(found==0){
00236         nStripSimu  = (int)(strip_simu[k].id_strip/100000.);
00237         idWafSimu   = strip_simu[k].id_strip-10000*((int)(strip_simu[k].id_strip/10000.));
00238         iWafSimu    = idWaferToWafer(idWafSimu);
00239         iLadSimu    = (int)(idWafSimu - mSsdLayer*1000 - (iWafSimu+1)*100 - 1);
00240         iSideSimu   = (strip_simu[k].id_strip - nStripSimu*100000 - idWafSimu)/10000;
00241         out_strip.id          = currRecord + 1;
00242         out_strip.adc_count   = (int)(strip_simu[k].adc_count);
00243         out_strip.id_strip    = 10000*(10*nStripSimu + iSideSimu)+idWafSimu;
00244         for (Int_t e = 0 ; e < 5 ; e++)
00245           {
00246             out_strip.id_mchit[e] = strip_simu[k].id_mchit[e];
00247             out_strip.id_mctrack[e] = strip_simu[k].id_mctrack[e];
00248           }
00249         mSsdSimuReal->AddAt(&out_strip);
00250         currRecord++;
00251         currRecordSimu++;
00252       }
00253     }
00254   LOG_INFO<<Form("# in common = %d totstrips after add = %d totstrips for simu only = %d for real only = %d",currRecordCommon,currRecord,currRecordSimu,currRecordReal)<<endm;
00255   return currRecord;
00256 }
00257 //____________________________________________________________________________
00258 Int_t StSsdEmbeddingMaker::Finish()
00259 { 
00260   return kStOK;
00261 }
00262 //____________________________________________________________________________
00263 void  StSsdEmbeddingMaker::setDoEmbedding(Bool_t doIt){
00264   mDoEmbedding = doIt;
00265 }
00266 //____________________________________________________________________________
00267 void StSsdEmbeddingMaker::setPlainSimEvenIfNoSSD(Bool_t doIt){
00268   mPlainSimIfNoSSD= doIt;
00269 }
00270 //____________________________________________________________________________
00271 Int_t StSsdEmbeddingMaker::NoSsd()
00272 {  
00273   St_DataSet *dataSet;
00274 
00275   dataSet = GetDataSet("StDAQReader");
00276   if(!dataSet){
00277     LOG_ERROR<<("BIG TROUBLE: cannot find StDAQReader in the chain and you want to run embedding")<<endm;
00278     return kTRUE;
00279   }
00280   StDAQReader* daqReader = (StDAQReader*)(dataSet->GetObject());
00281   if (!daqReader){
00282     LOG_ERROR<<("BIG TROUBLE: StDAQReader is empty and you want to run embedding")<<endm;
00283     return kTRUE;
00284   }
00285     if (!daqReader->SSDPresent ())
00286     {
00287       LOG_INFO<<("NO SSD in DAQ")<<endm;
00288       return kTRUE; //No SSD in the datastream 
00289     }
00290   
00291   LOG_INFO<<("SSD found in DAQ")<<endm;
00292   return kFALSE;
00293 }
00294 //______________________________________________________________________
00295 void StSsdEmbeddingMaker::CheckTables(){
00296  spa_strip_st *strip = mSsdSimuData->GetTable();
00297   
00298   Int_t idWaf         = 0;
00299   Int_t iWaf          = 0;
00300   Int_t iLad          = 0;
00301   Int_t nStrip        = 0;
00302   Int_t iSide         = 0;
00303   Int_t idMcHit[5]    = {0,0,0,0,0};
00304   Int_t i=0 ,e = 0;
00305   LOG_DEBUG<<"check simu table :" << endm;
00306   for (i = 0 ; i < mSsdSimuData->GetNRows(); i++)
00307     {
00308       nStrip  = (int)(strip[i].id_strip/100000.);
00309       idWaf   = strip[i].id_strip-10000*((int)(strip[i].id_strip/10000.));
00310       iWaf    = idWaferToWafer(idWaf);
00311       iLad    = (int)(idWaf - mSsdLayer*1000 - (iWaf+1)*100 - 1);
00312       iSide   = (strip[i].id_strip - nStrip*100000 - idWaf)/10000;
00313       for (e = 0 ; e < 5;e++) {idMcHit[e] = strip[i].id_mchit[e];}
00314       LOG_DEBUG<<Form("side=%d idWafer=%d Ladder=%d wafer=%d nstrip=%d signal=%d",iSide,idWaf,iLad,iWaf,nStrip,strip[i].adc_count)<<endm;
00315       for(Int_t ii=0;ii<5;ii++){printf("id_mchit[%d] =%d\n",ii,strip[i].id_mchit[ii]);}
00316     }
00317   if(Debug()){
00318     LOG_DEBUG << "check real table : " << endm;
00319     strip  = mSsdrealData->GetTable();
00320     for (i = 0 ; i < mSsdrealData->GetNRows(); i++)
00321       {
00322         nStrip  = (int)(strip[i].id_strip/100000.);
00323         idWaf   = strip[i].id_strip-10000*((int)(strip[i].id_strip/10000.));
00324         iWaf    = idWaferToWafer(idWaf);
00325         iLad    = (int)(idWaf - mSsdLayer*1000 - (iWaf+1)*100 - 1);
00326         iSide   = (strip[i].id_strip - nStrip*100000 - idWaf)/10000;
00327         for (e = 0 ; e < 5;e++) {idMcHit[e] = strip[i].id_mchit[e];
00328         LOG_DEBUG<<Form("side=%d Ladder=%d wafer=%d nstrip=%d signal=%d",iSide,iLad,iWaf,nStrip,strip[i].adc_count)<<endm;
00329         }
00330       }
00331   }
00332 }
00333 //__________________________________________________________________________
00334 void StSsdEmbeddingMaker::setSsdParameters(ssdDimensions_st *geom_par){
00335   mDimensions          = geom_par;
00336   mSsdLayer            = 7; 
00337   mDetectorLargeEdge   = 2.*geom_par[0].waferHalfActLength;
00338   mDetectorSmallEdge   = 2.*geom_par[0].waferHalfActWidth;
00339   mNLadder             = 20;
00340   mNWaferPerLadder     = geom_par[0].wafersPerLadder;
00341   mNStripPerSide       = geom_par[0].stripPerSide;
00342   mStripPitch          = geom_par[0].stripPitch;
00343   mTheta               = geom_par[0].stereoAngle;
00344 }
00345 //___________________________________________________________________________
00346