StRoot: StSsdFastSimMaker/StSsdFastSimMaker.cxx Source File

00001 /*
00002  *
00003  * Author: J. Bouchet, Subatech; Y. Fisyak, BNL; J. Vanfossen, Kent State
00004  *
00005  * 
00006  **********************************************************
00007  *
00008  */
00009 
00010 #include "Stiostream.h"
00011 #include "StSsdFastSimMaker.h"
00012 #include "StHit.h"
00013 #include "StEventTypes.h"
00014 #include "StChain.h"
00015 
00016 #include "Sti/StiVMCToolKit.h"
00017 #include "StarClassLibrary/StRandom.hh"
00018 #include "tables/St_HitError_Table.h"
00019 
00020 #include "TH1.h"
00021 #include "TH2.h"
00022 #include <math.h>
00023 
00024 #include "tables/St_g2t_ssd_hit_Table.h"
00025 #include "StSsdDbMaker/StSsdDbMaker.h"
00026 #include "StSsdUtil/StSsdBarrel.hh"
00027 #include "StSsdHitCollection.h"
00028 #include "StEvent.h"
00029 #include "StMcEvent.hh"
00030 ClassImp(StSsdFastSimMaker)
00031   
00032   StSsdFastSimMaker::StSsdFastSimMaker(const char *name):StMaker(name){
00033     mHit        = new StSsdHit();
00034     WaferNumb   = 0;
00035     HitsMap     = 0;
00036     dX          = 0;
00037     dY          = 0;
00038     dZ          = 0;
00039     Local       = 0;
00040     Ratio       = 0;
00041   };
00042 
00043   StSsdFastSimMaker::~StSsdFastSimMaker(){ 
00044     delete mHit;
00045     delete WaferNumb   ; 
00046     delete HitsMap     ; 
00047     delete dX          ; 
00048     delete dY          ; 
00049     delete dZ          ; 
00050     delete Local       ; 
00051     delete Ratio       ;
00052 }
00053 
00054 int StSsdFastSimMaker::Init()
00055 {
00056   int seed=time(NULL);
00057   myRandom=new StRandom();
00058   myRandom->setSeed(seed);
00059 
00060   // Define various SSD specific geom. variables.
00061   mResXSsd = 0.0000;                            
00062   mResZSsd = 0.0000;
00063   mSmear=1;     //used to turn on or off smearing
00064   
00065   LOG_DEBUG << "Init() : Defining the histograms" << endm;
00066   if(IAttr(".histos")){
00067     WaferNumb = new TH1S("WaferNumb","Hits from geant per Wafer Id",1620,7000,8620);
00068     WaferNumb->SetXTitle("from 7000 to 7000*(iW*100)+16");
00069     WaferNumb->SetLabelSize(0.03,"X");
00070     WaferNumb->SetLabelSize(0.03,"Y");
00071     HitsMap = new TH2S("HitsMap","Map of hits in ssd wafers",20,1,21,16,1,17);
00072     HitsMap->SetXTitle("Ladder");
00073     HitsMap->SetYTitle("Wafer");
00074     HitsMap->SetLabelSize(0.03,"X");
00075     HitsMap->SetLabelSize(0.03,"Y");
00076     dX = new TH1F("dX","difference in X after smearing the hit",200,-0.02,0.02);
00077     dX->SetXTitle("x_{G} - x_{G}^{smear} [cm]");
00078     dX->SetLabelSize(0.03,"X");
00079     dX->SetLabelSize(0.03,"Y");
00080     dY = new TH1F("dY","difference in Y after smearing the hit",200,-0.02,0.02);
00081     dY->SetXTitle("y_{G} - y_{G}^{smear} [cm]");
00082     dY->SetLabelSize(0.03,"X");
00083     dY->SetLabelSize(0.03,"Y");
00084     dZ = new TH1F("dZ","difference in Z after smearing the hit",800,-2,2);
00085     dZ->SetXTitle("z_{G} - z_{G}^{smear} [cm]");
00086     dZ->SetLabelSize(0.03,"X");
00087     dZ->SetLabelSize(0.03,"Y");
00088     Local = new TH2F("Local","Local x vs local y : hits ",800,-4,4,250,-2.25,2.25);
00089     Local->SetXTitle("x_{l}[cm]");
00090     Local->SetYTitle("y_{l}[cm]");
00091     Local->SetLabelSize(0.03,"X");
00092     Local->SetLabelSize(0.03,"Y");
00093     Ratio  = new TH1F("Ratio","Ratio between initial geant hits and after removal of the dead areas",120,0,1.2);
00094     Ratio->SetLabelSize(0.03,"X");
00095     Ratio->SetLabelSize(0.03,"Y");
00096   }
00097   return kStOk;
00098 }
00099 //____________________________________________________________
00100 int StSsdFastSimMaker::InitRun(int RunNo)
00101 {
00102   // Define various SSD hit errors from database
00103   TDataSet *set = GetDataBase("Calibrations/tracker");
00104   St_HitError *tableSet = (St_HitError *)set->Find("ssdHitError");
00105   HitError_st* hitError = tableSet->GetTable();
00106   mResXSsd = sqrt(hitError->coeff[0]);
00107   mResZSsd = sqrt(hitError->coeff[3]);
00108   LOG_DEBUG << "Smearing SSD hits by " << mResXSsd << " " << mResZSsd << " cm in the SSD " << endm;
00109   
00110   //            geometry parameters
00111   TDataSet *ssdparams = GetInputDB("Geometry/ssd");
00112   if (! ssdparams) {
00113     LOG_ERROR << "No  access to Geometry/ssd parameters" << endm;
00114     return kStFatal;
00115   }
00116   TDataSetIter    local(ssdparams);
00117   m_ctrl        = (St_slsCtrl           *)local("slsCtrl");
00118   m_dimensions  = (St_ssdDimensions     *)local("ssdDimensions"); 
00119   m_positions   = (St_ssdWafersPosition *)local("ssdWafersPosition");
00120   
00121   if (!m_ctrl) {
00122     LOG_ERROR << "No  access to control parameters" << endm;
00123     return kStFatal;
00124   }   
00125   if ((!m_dimensions)||(!m_positions)) {
00126     LOG_ERROR << "No  access to geometry parameters" << endm;
00127     return kStFatal;
00128   }  
00129   St_ssdConfiguration* configTable = (St_ssdConfiguration*) local("ssdConfiguration");
00130   if (!configTable) {
00131     LOG_ERROR << "InitRun : No access to ssdConfiguration database" << endm;
00132     return kStFatal;
00133   }
00134   m_config = (ssdConfiguration_st*) configTable->GetTable() ;
00135   
00136   return kStOk;
00137 }
00138 //____________________________________________________________
00139 Int_t StSsdFastSimMaker::Finish(){return kStOk;}
00140 //____________________________________________________________
00141 Int_t StSsdFastSimMaker::Make()
00142 {
00143   // Get the input data structures from StEvent
00144   mEvent =  (StEvent*) GetInputDS("StEvent");
00145   if(mEvent)
00146     {
00147       mCol = mEvent->ssdHitCollection();
00148       if (!mCol)  
00149         { 
00150           LOG_WARN <<"StSsdFastSimulatorMaker -E- no SsdHitCollection!" << endm; 
00151           mCol = new StSsdHitCollection; 
00152           mEvent->setSsdHitCollection(mCol); 
00153           LOG_WARN <<"Make() has added a non existing StSsdHitCollection" <<endm; 
00154         }   
00155     }
00156   if (! mEvent) {
00157     LOG_WARN << "No StEvent on input, bye bye" << endm; return kStWarn;
00158     mCol =0;
00159   }
00160   //StMcEvent* mevent = (StMcEvent*) GetInputDS("StMcEvent");
00161   //if(!mevent){LOG_WARN << "No StMcEvent on input "<<endm; return kStWarn;}
00162   TDataSetIter geant(GetInputDS("geant"));
00163   if (! gGeoManager) GetDataBase("VmcGeometry");
00164   LOG_DEBUG << "Geometry Loaded" << endm;
00165   St_g2t_ssd_hit *g2t_ssd_hit = (St_g2t_ssd_hit *) geant("g2t_ssd_hit");
00166   if (! g2t_ssd_hit) {
00167     LOG_WARN << "No g2t_ssd_hit on input, bye bye" << endm; return kStWarn;
00168     return kStWarn;
00169   }
00170   g2t_ssd_hit_st *g2t         = g2t_ssd_hit->GetTable();
00171 
00172   LOG_INFO<<"####      START OF SSD FAST SIM MAKER        ####"<<endm;
00173   StSsdBarrel *mySsd =gStSsdDbMaker->GetSsd();
00174   ssdDimensions_st *dimensions = m_dimensions->GetTable();
00175   setSsdParameters(dimensions);
00176   if(Debug()>1) printSsdParameters();
00177   Int_t inContainer = 0;
00178   Int_t goodHits    = 0;
00179   if(g2t_ssd_hit){
00180     LOG_INFO<<Form("NumberOfRows = %d Size of g2t_ssd table=%d",g2t_ssd_hit->GetNRows(),g2t_ssd_hit->GetTableSize())<<endm;
00181     Int_t minWaf      = mSsdLayer*1000;
00182     Int_t currWafId   = 0;
00183     Int_t currWafNumb = 0;
00184     Int_t currLadder  = 0;
00185     Int_t i           = 0;
00186     unsigned char c   = 0;
00187     Int_t iok         = 0;
00188     Int_t in          = 0;
00189     for (i = 0; i < g2t_ssd_hit->GetNRows() ; i++)    {
00190       currWafId=g2t[i].volume_id%10000;
00191       if (currWafId > minWaf)   {
00192         currLadder  = StSsdBarrel::Instance()->idWaferToLadderNumb(currWafId);
00193         currWafNumb = StSsdBarrel::Instance()->idWaferToWafer(currWafId);
00194         StThreeVectorF error(0.,0.,0.);
00195         mHit = new StSsdHit(g2t[i].x,error,currWafId,g2t[i].de,c);
00196         LOG_DEBUG<<Form("currWafId=%d currWafNum=%d currLadder=%d",currWafId,currWafNumb,currLadder)<<endm;
00197         // now we get the position of this wafer    
00198         Double_t xg[3] = {mHit->position().x(),mHit->position().y(),mHit->position().z()};
00199         LOG_DEBUG <<Form("global position x=%f y=%f z=%f",xg[0],xg[1],xg[2])<<endm;
00200         Double_t xl[3];
00201         xl[0] =0.0;  
00202         xl[1] =0.0;  
00203         xl[2] =0.0;
00204         mySsd->mLadders[currLadder]->mWafers[currWafNumb]->MasterToLocal(xg,xl);
00205         LOG_DEBUG <<Form("local position x=%f y=%f z=%f",xl[0],xl[1],xl[2])<<endm;
00206         LOG_DEBUG <<"will smear X"<<endm;
00207         Double_t xlSmear[3];
00208         xlSmear[0] =0.0; 
00209         xlSmear[1] =0.0; 
00210         xlSmear[2] =0.0;
00211         xlSmear[0] = (distortHit(xl[0], mResXSsd));
00212         LOG_DEBUG <<"will smear Y"<<endm;
00213         xlSmear[1] = (distortHit(xl[1], mResZSsd));
00214         xlSmear[2] = xl[2];
00215         LOG_DEBUG <<Form("smeared local position x=%f y=%f z=%f",xlSmear[0],xlSmear[1],xlSmear[2])<<endm;
00216         Double_t xgSmear[3];
00217         xgSmear[0] =0.0;  
00218         xgSmear[1] =0.0;  
00219         xgSmear[2] =0.0;  
00220         mySsd->mLadders[currLadder]->mWafers[currWafNumb]->LocalToMaster(xlSmear,xgSmear);
00221         LOG_DEBUG<< "After : global position are :"<<endm;
00222         LOG_DEBUG <<Form("x=%f y=%f z=%f",xgSmear[0],xgSmear[1],xgSmear[2])<<endm;
00223         LOG_DEBUG <<Form("Smearing done...")<< endm;
00224         in  = IsOnWafer(xlSmear[0],xlSmear[1]);
00225         if(in==0)continue;
00226         LOG_DEBUG << "good hit in wafer active area " << endm;
00227         iok = RemoveTriangle(xlSmear[0],xlSmear[1]);
00228         if(iok==0)continue;
00229         goodHits++;
00230         LOG_DEBUG << "good hit after triangle rejection" << endm;
00231         Local->Fill(xlSmear[0],xlSmear[1]);
00232         StSsdHit  *mSmearedhit = new StSsdHit(xgSmear,error,currWafId,g2t[i].de,c);
00233         //fill histograms
00234         if(IAttr(".histos")){ 
00235           HitsMap->Fill(currLadder+1,currWafNumb+1,1);
00236           WaferNumb->Fill(currWafId);
00237           dX->Fill(mHit->position().x()-xgSmear[0]);
00238           dY->Fill(mHit->position().y()-xgSmear[1]);
00239           dZ->Fill(mHit->position().z()-xgSmear[2]);
00240         }
00241           // encode the hardware position
00242         // 2^3  detector ID number (8) 
00243         // 2^4  4-12 num_wafer (0-319)
00244         Int_t hw  =         
00245         8                                                                             
00246         +         16 * idWaferToWaferNumb(currWafId);                                        
00247         mSmearedhit->setHardwarePosition(hw);
00248         mSmearedhit->setLocalPosition(xlSmear[0],xlSmear[1]);
00249         mSmearedhit->setIdTruth(g2t[i].track_p,100);
00250         //Float_t p[3]      = {0,0,0};
00251         //Int_t j =0 ;
00252         //for (j = 0; j<3; j++) {p[j] = g2t[i].p[j];}
00253         //mySsd->mLadders[currLadder]->mWafers[currWafNumb]->addHit(g2t[i].id, g2t[i].id, g2t[i].track_p, g2t[i].x, g2t[i].de, p);
00254         inContainer+=mCol->addHit(mSmearedhit);
00255       }
00256     }
00257   }
00258   if(inContainer){
00259     LOG_INFO<<"####   -> "<<inContainer<<" HITS WRITTEN INTO CONTAINER   ####"<<endm;
00260     Ratio->Fill(inContainer/(float)g2t_ssd_hit->GetNRows());
00261   }
00262   else {
00263     LOG_INFO<<"######### NO SSD HITS WRITTEN INTO CONTAINER  ####"<<endm;
00264   }
00265   LOG_DEBUG << "ssd col= "<<mCol->numberOfHits()<< " inContainer=" << inContainer <<" good hits = " << goodHits <<" initial # hits="<<g2t_ssd_hit->GetNRows()<<endm;
00266   mySsd->Reset();
00267   return kStOK;
00268 }
00269 //________________________________________________________________________________
00270 Bool_t StSsdFastSimMaker::accept(StEvent* event){
00271   return event ? true : false;
00272 }
00273 //_______________________________________________________________________________
00274 Float_t StSsdFastSimMaker::distortHit(Float_t x, double res){
00275   double test;
00276   LOG_DEBUG << "In distortHit " << endm;
00277   LOG_DEBUG << Form("smear me by %f",res)<<endm;
00278   if(mSmear){
00279     test = x + myRandom->gauss(0,res);
00280     
00281     LOG_DEBUG<< " x was " <<x<< " and is now " << test << endm;
00282     return test;
00283   }
00284   else return x;
00285 }
00286 //_______________________________________________________________
00287 Int_t StSsdFastSimMaker::idWaferToWaferNumb(Int_t idWafer)
00288 {
00289   // idwafer = layer*1000+waf*100+ladder
00290   Int_t iW = (int)((idWafer - mSsdLayer*1000)/100);
00291   Int_t iL = idWafer - mSsdLayer*1000 - iW*100;
00292   return ((iL-1)*mNWaferPerLadder + iW -1);
00293 }
00294 //_______________________________________________________________
00295 Int_t StSsdFastSimMaker::idWaferToLadderNumb(Int_t idWafer)
00296 {
00297   // idwafer = layer*1000+waf*100+ladder
00298   Int_t iW = (int)((idWafer - mSsdLayer*1000)/100);
00299   Int_t iL = idWafer - mSsdLayer*1000 - iW*100;
00300   return iL-1;
00301 }
00302 //_______________________________________________________________
00303 Int_t StSsdFastSimMaker::waferNumbToIdWafer(Int_t waferNumb)
00304 {
00305   Int_t iL = 1+(int)((waferNumb)/mNLadder);
00306   Int_t iW = waferNumb-((iL-1)*mNLadder)+1;
00307   return mSsdLayer*1000 + iW*100 + iL;
00308 }
00309 //________________________________________________________________________________
00310 void StSsdFastSimMaker::setSsdParameters(ssdDimensions_st *geom_par){
00311   mDimensions          = geom_par;
00312   mSsdLayer            = 7; // all layers : 1->7
00313   mDetectorLargeEdge   = 2.*geom_par[0].waferHalfActLength;
00314   mDetectorSmallEdge   = 2.*geom_par[0].waferHalfActWidth;
00315   mNLadder             = 20;
00316   mNWaferPerLadder     = geom_par[0].wafersPerLadder;
00317   mNStripPerSide       = geom_par[0].stripPerSide;
00318   mStripPitch          = geom_par[0].stripPitch;
00319   mTheta               = geom_par[0].stereoAngle;
00320 }
00321 //________________________________________________________________________________
00322 void StSsdFastSimMaker::printSsdParameters(){
00323   LOG_INFO << "###Ladders = " <<mNLadder<<"###"<< endm;
00324   LOG_INFO << "###Wafers per Ladder = " <<mNWaferPerLadder<<"###"<< endm;
00325 }
00326 void StSsdFastSimMaker::Clear(Option_t *option)
00327 {
00328   /*noop*/
00329 }
00330 //________________________________________________________________________________
00331 Int_t StSsdFastSimMaker::RemoveTriangle(Float_t x, Float_t y)
00332 {
00333   Int_t iok = 0; 
00334   Float_t beta = (TMath::Pi()/2.) - (mTheta/2.);
00335   Float_t X = (mDetectorSmallEdge/2.)*TMath::Tan(mTheta/2.);
00336   Float_t l = (mDetectorLargeEdge/2.) -X ;
00337   LOG_DEBUG <<Form("beta=%f arctan(beta)=%f X=%f l=%f mdetectorsmallEdge=%f mDetectorLargeEdge=%f\n",beta,TMath::ATan(beta),X,l,mDetectorSmallEdge,mDetectorLargeEdge)<<endm;
00338   if(x>(mDetectorLargeEdge/2.) -X){
00339     if(y < TMath::ATan(beta)*(x-l)){
00340       LOG_DEBUG<<"upper right corner , hit rejected at x=" << x <<" y=" <<y << endm;
00341       iok =0;
00342     }
00343     else 
00344       if(y>-1*TMath::ATan(beta)*(x-l)){
00345         LOG_DEBUG<<"bottom right corner , hit rejected at x=" << x <<" y=" <<y << endm;
00346           iok = 0; //bottom right corner
00347       }
00348   }
00349   if(x<((-1*mDetectorLargeEdge/2.)+X)){
00350     if(y<-1*TMath::ATan(beta)*(x+l)){
00351       LOG_DEBUG<<"upper  left corner , hit rejected  at x=" << x <<" y=" <<y << endm;
00352       iok =0; //upper left corner
00353     }
00354     else {
00355       if(y>1*TMath::ATan(beta)*(x+l)){
00356         LOG_DEBUG<<"bottom left corner , hit rejected  at x=" << x <<" y=" <<y << endm;
00357         iok = 0; //bottom left corner
00358       }
00359     }
00360   }  
00361   else{ iok= 1;}
00362 return iok; // good hit
00363 }
00364 //_______________________________________________________________________________
00365 int StSsdFastSimMaker::IsOnWafer(Float_t x , Float_t y){
00366   //Find out for a given z coord and Hardware pos is it on the wafer
00367   if((x <(mDetectorLargeEdge/2.)) && (x > (-mDetectorLargeEdge/2.)) &&  
00368      (y <(mDetectorSmallEdge/2.)) && (y > (-mDetectorSmallEdge/2.)))
00369     return 1;   
00370   return 0;
00371 }