StRoot: StPixelFastSimMaker/StPixelFastSimMaker.cxx Source File

00001 /*
00002  *
00003  * Author: A. Rose, LBL, Y. Fisyak, BNL, M. Miller, MIT
00004  *
00005  * 
00006  **********************************************************
00007  * $Log: StPixelFastSimMaker.cxx,v $
00008  * Revision 1.43  2010/09/01 20:31:47  fisyak
00009  * clean up unused varaibles
00010  *
00011  * Revision 1.42  2009/02/06 20:48:48  wleight
00012  * UPGR15 Update
00013  *
00014  * Revision 1.41  2009/01/26 14:50:46  fisyak
00015  * Clean up
00016  *
00017  * Revision 1.40  2007/12/03 20:42:43  wleight
00018  * Replaced couts with LOG_DEBUGs
00019  *
00020  * Revision 1.37  2007/11/13 19:09:51  wleight
00021  * Corrected bug causing pixel fast simulator to crash when there were no pixel and/or ist hits in the event
00022  *
00023  * Revision 1.36  2007/11/06 16:20:06  wleight
00024  * Digitized Pixel, removed all hit smearing, and implemented idTruth
00025  *
00026  * Revision 1.34  2007/10/18 16:31:44  fisyak
00027  * Add pile-up from weixie
00028  *
00029  * Revision 1.33  2007/10/18 14:25:13  didenko
00030  * updates for pile-up events
00031  * 
00032  * Revision 1.32  2007/10/16 19:53:08  fisyak
00033  * rename Hft => Pxl, remove Hpd, Igt and Fst
00034  *
00035  * Revision 1.31  2007/10/16 19:50:46  fisyak
00036  * rename Hft => Pxl, remove Hpd, Igt and Fst
00037  *
00038  * Revision 1.30  2007/09/09 17:00:32  fisyak
00039  * Fix bug 1056
00040  *
00041  * Revision 1.29  2007/09/08 00:33:05  andrewar
00042  * Modifications for pileup hit read in.
00043  *
00044  * Revision 1.28  2007/05/17 13:18:52  andrewar
00045  * Removed cout in shiftHit.
00046  *
00047  * Revision 1.27  2007/05/16 15:06:55  andrewar
00048  * Switched cout's to LOG_INFO.
00049  *
00050  * Revision 1.26  2007/04/28 17:56:36  perev
00051  * Redundant StChain.h removed
00052  *
00053  * Revision 1.25  2007/04/27 18:41:29  wleight
00054  * Removed smearing of the coordinate not controlled by the strips in the 17cm layer
00055  *
00056  * Revision 1.24  2007/04/27 14:59:10  wleight
00057  * Corrected another error in the creation of new hits
00058  *
00059  * Revision 1.23  2007/04/26 04:08:41  perev
00060  * Remove StBFChain dependency
00061  *
00062  * Revision 1.22  2007/04/25 17:44:59  wleight
00063  * Corrected error in assignment of reconstructed IST hits
00064  *
00065  * Revision 1.21  2007/04/23 18:11:30  andrewar
00066  * Removed references to Hpd (includes were obsolete)
00067  *
00068  * Revision 1.19  2007/04/23 16:32:47  wleight
00069  * Added explicit casting for double to int in calculating strip number
00070  *
00071  * Revision 1.18  2007/04/22 22:57:23  wleight
00072  * The two hits in the 17 cm layer are no longer combined into 1
00073  *
00074  * Revision 1.17  2007/04/16 19:10:52  wleight
00075  * Added IST simulation (digitization but no clustering)
00076  *
00077  * Revision 1.16  2007/04/13 19:17:15  andrewar
00078  * Removed misleading errors. Changed cout and printf to gMessMgr.
00079  *
00080  * Revision 1.15  2007/04/06 21:46:36  andrewar
00081  * Removed some debug messages.
00082  *
00083  * Revision 1.14  2007/04/06 14:55:11  andrewar
00084  * Shift of HFT hit to face of ladder.
00085  *
00086  * Revision 1.13  2007/03/28 13:33:45  mmiller
00087  * Removed cout/printf's.
00088  *
00089  * Revision 1.12  2006/12/21 18:11:59  wleight
00090  * Fixed UPGR09 compatibility so it works with all versions
00091  *
00092  * Revision 1.11  2006/12/20 16:50:21  wleight
00093  * Added fix for UPGR09 problem with layer number mismatch
00094  *
00095  * Revision 1.10  2006/12/15 02:17:20  wleight
00096  * Ist now gets hit smearing parameters from the database
00097  *
00098  * Revision 1.9  2006/12/14 23:52:51  andrewar
00099  * Added Sevil's hit error db loader.
00100  *
00101  * Revision 1.7  2006/11/29 21:42:07  andrewar
00102  * Update with Pixel resolution smearing.
00103  *
00104  * Revision 1.6  2006/11/28 22:37:42  wleight
00105  * Fixed minor smearing bug
00106  *
00107  * Revision 1.4  2006/10/13 20:15:45  fisyak
00108  * Add Hpd fast simulation (Sevil)
00109  *
00110  * Revision 1.3  2006/02/17 21:44:29  andrewar
00111  * Remover streaming of each Pixel hit.
00112  *
00113  * Revision 1.2  2006/02/08 20:57:33  fisyak
00114  * Set proper Detector Id
00115  *
00116  * Revision 1.1  2006/02/03 20:11:56  fisyak
00117  * The initial revision
00118  *
00119  *
00120  */
00121 
00122 #include "Stiostream.h"
00123 #include "StPixelFastSimMaker.h"
00124 #include "StHit.h"
00125 #include "StEventTypes.h"
00126 #include "StEvent.h"
00127 #include "StRnDHit.h"
00128 #include "StMcEvent.hh"
00129 #include "StMcHit.hh"
00130 #include "StMcIstHit.hh"
00131 
00132 #include "StMcPixelHit.hh"
00133 #include "StMcEventTypes.hh"
00134 #ifdef bug1056
00135 #include "Sti/Base/Factory.h"
00136 #include "Sti/StiPlanarShape.h"
00137 #include "Sti/StiCylindricalShape.h"
00138 #include "Sti/StiMaterial.h"
00139 #include "Sti/StiPlacement.h"
00140 #include "Sti/StiDetector.h"
00141 #include "Sti/StiToolkit.h"
00142 #include "Sti/StiDetectorBuilder.h"
00143 #endif
00144 #include <stdio.h>
00145 #include <map>
00146 #include <exception>
00147 using namespace std;
00148 #include <stdexcept>
00149 #include "tables/St_g2t_ist_hit_Table.h"
00150 #include "tables/St_g2t_pix_hit_Table.h"
00151 #include "tables/St_HitError_Table.h"
00152 #ifdef bug1056
00153 #include "Sti/StiHitErrorCalculator.h"
00154 #include "Sti/StiIsActiveFunctor.h"
00155 #include "Sti/StiNeverActiveFunctor.h"
00156 #include "Sti/StiElossCalculator.h"
00157 #include "Sti/StiVMCToolKit.h"
00158 #else
00159 #include "TGeoManager.h"
00160 #include "TGeoMatrix.h"
00161 #endif
00162 #include "StarClassLibrary/StRandom.hh"
00163 #include "tables/St_HitError_Table.h"
00164 #include <fstream>
00165 #include "TFile.h"
00166 #include "TTree.h"
00167 
00168 ClassImp(StPixelFastSimMaker)
00169 
00170   
00171   StPixelFastSimMaker::~StPixelFastSimMaker(){ /*noop*/ }
00172 
00173 int StPixelFastSimMaker::Init()
00174 {
00175   LOG_INFO<<"StPixelFastSimMaker::Init()"<<endm;
00176   int seed=time(NULL);
00177   myRandom=new StRandom();
00178   myRandom->setSeed(seed);
00179 
00180  
00181   mSmear=1;
00182   return kStOk;
00183 }
00184 
00185 //____________________________________________________________
00186 int StPixelFastSimMaker::InitRun(int RunNo)
00187 {
00188   LOG_INFO<<"StPixelFastSimMaker::InitRun"<<endm;
00189 
00190   TDataSet *set = GetDataBase("Calibrations/tracker");
00191   St_HitError *istTableSet = (St_HitError *)set->Find("ist1HitError");
00192   St_HitError *pixelTableSet = (St_HitError *)set->Find("PixelHitError");
00193   //cout<<"found St_HitErrors"<<endl;
00194   HitError_st* istHitError = istTableSet->GetTable();
00195   resXIst1 = sqrt(istHitError->coeff[0]);
00196   resZIst1 = sqrt(istHitError->coeff[3]);
00197   HitError_st* pixelHitError = pixelTableSet->GetTable();
00198   resXPix = sqrt(pixelHitError->coeff[0]);
00199   resZPix = sqrt(pixelHitError->coeff[3]);
00200 
00201   LoadPixPileUpHits(); //.. load the pile up hits for PIXEL
00202   //cout<<"done with StPixelFastSimMaker::InitRun"<<endl;
00203 
00204   return kStOk;
00205 }
00206 
00207 //___________________________
00208 void StPixelFastSimMaker::LoadPixPileUpHits()
00209 {
00210   LOG_INFO<<"+++ loading the PIXEL pileup files +++"<<endm;
00211 
00212   pileup_on = true;
00213 
00214   TFile f_pileup("pileup.root");
00215   if (f_pileup.IsZombie()) {
00216 
00217     pileup_on = false;
00218 
00219     LOG_INFO << "no PIXEL pileup file found. Will run with regular setup" << endm;
00220     return;
00221   }
00222 
00223   LOG_INFO<<"+++ Loaded pileup.root for PIXEL pileup simulation +++"<<endm;
00224 
00225   TTree* pileup_tree = (TTree*)f_pileup.Get("pileup_tree");
00226 
00227   const int maxhit = 200000;
00228   float x[maxhit], y[maxhit], z[maxhit], px[maxhit], py[maxhit], pz[maxhit], de[maxhit], ds[maxhit];
00229   long key[maxhit], vid[maxhit];
00230   int layer[maxhit], nhits;
00231 
00232   TBranch *b_x = pileup_tree->GetBranch("x");
00233   TBranch *b_y = pileup_tree->GetBranch("y");
00234   TBranch *b_z = pileup_tree->GetBranch("z");
00235   TBranch *b_px = pileup_tree->GetBranch("px");
00236   TBranch *b_py = pileup_tree->GetBranch("py");
00237   TBranch *b_pz = pileup_tree->GetBranch("pz");
00238   TBranch *b_de = pileup_tree->GetBranch("de");
00239   TBranch *b_ds = pileup_tree->GetBranch("ds");
00240   TBranch *b_key = pileup_tree->GetBranch("key");
00241   TBranch *b_vid = pileup_tree->GetBranch("vid");
00242   TBranch *b_layer = pileup_tree->GetBranch("layer");
00243   TBranch *b_nhits = pileup_tree->GetBranch("nhits");
00244   b_x->SetAddress(x);
00245   b_y->SetAddress(y);
00246   b_z->SetAddress(z);
00247   b_px->SetAddress(px);
00248   b_py->SetAddress(py);
00249   b_pz->SetAddress(pz);
00250   b_de->SetAddress(de);
00251   b_ds->SetAddress(de);
00252   b_key->SetAddress(key);
00253   b_vid->SetAddress(vid);
00254   b_layer->SetAddress(layer);
00255   b_nhits->SetAddress(&nhits);
00256 
00257   pileup_tree->GetEntry(0); //.. just one events
00258 
00259   for(int ihit = 0; ihit<nhits; ihit++) {
00260     pileup_x.push_back(x[ihit]);
00261     pileup_y.push_back(y[ihit]);
00262     pileup_z.push_back(z[ihit]);
00263 
00264     pileup_px.push_back(px[ihit]);
00265     pileup_py.push_back(py[ihit]);
00266     pileup_pz.push_back(pz[ihit]);
00267 
00268     pileup_key.push_back(key[ihit]);
00269     pileup_vid.push_back(vid[ihit]);
00270 
00271     pileup_de.push_back(de[ihit]);
00272     pileup_ds.push_back(ds[ihit]);
00273   }
00274 }
00275 //____________________________________________________________
00276 void StPixelFastSimMaker::AddPixPileUpHit(StMcPixelHitCollection* pixHitCol)
00277 {
00278   for(unsigned int i = 0; i<pileup_x.size(); i++) {
00279 
00280     StThreeVectorD pos(pileup_x[i], pileup_y[i], pileup_z[i]);
00281     StThreeVectorF mom(pileup_px[i], pileup_py[i], pileup_pz[i]);
00282 
00283     float de = pileup_de[i]; 
00284     float ds = pileup_ds[i];
00285 
00286     int key = pileup_key[i];
00287     int vid = pileup_vid[i];
00288 
00289     StMcPixelHit* pixhit = new StMcPixelHit(pos, mom, de, ds, key, vid, 0);
00290     pixHitCol->addHit(pixhit);
00291   }
00292 }
00293 //____________________________________________________________
00294 
00295 Int_t StPixelFastSimMaker::Make()
00296 {
00297   LOG_INFO<<"StPixelFastSimMaker::Make()"<<endm;
00298 
00299   // Get the input data structures from StEvent and StMcEvent
00300   StEvent* rcEvent =  (StEvent*) GetInputDS("StEvent");
00301   if (! rcEvent) {LOG_INFO << "No StEvent on input" << endl; return kStWarn;}
00302   StMcEvent* mcEvent = (StMcEvent *) GetInputDS("StMcEvent");
00303   if (! mcEvent) {LOG_INFO << "No StMcEvent on input" << endl; return kStWarn;}
00304   TDataSetIter geant(GetInputDS("geant"));
00305   if (! gGeoManager) GetDataBase("VmcGeometry");
00306   g2t_ist_hit_st* g2tIst=0;
00307   g2t_pix_hit_st* g2tPix=0;
00308   St_g2t_ist_hit *g2t_ist_hit=(St_g2t_ist_hit *)geant("g2t_ist_hit");
00309   St_g2t_pix_hit *g2t_pix_hit=(St_g2t_pix_hit *)geant("g2t_pix_hit");
00310   if(g2t_ist_hit) g2tIst=g2t_ist_hit->GetTable();
00311   if(g2t_pix_hit) g2tPix=g2t_pix_hit->GetTable();
00312 
00313   // Store hits into RnD Hit Collection until we have our own
00314   StRnDHitCollection *col = new StRnDHitCollection;
00315   if (!col )
00316     {
00317       gMessMgr->Info()<<"StPixelFastSimMaker -E- no RnDHitCollection!\n";
00318       abort();
00319     }
00320 
00321   // Don't use realistic hit errors for now. When we transit to smeared
00322   // hits, this would be a good place to store offset info
00323   StThreeVectorF mHitError(0.,0.,0.);
00324 
00325   //Get MC Pixel hit collection. This contains all pixel hits.
00326   StMcPixelHitCollection* pixHitCol = mcEvent->pixelHitCollection();
00327   unsigned int nPixLadders[2]={9,24};
00328   //int pixels[6400][640];
00329   vector<int> pixels;
00330   vector<StMcPixelHit*> pixLadderHits;
00331   multimap<int,int> pixelToKey;
00332   float smearedX,smearedZ;
00333 
00334   if (pixHitCol){
00335     int counter_layer1 = 0;
00336     int counter_layer2 = 0;
00337     int counter        = 0;
00338     float rad          = 0;
00339     if(pileup_on) AddPixPileUpHit(pixHitCol); //.. add the pileup hits into the collection
00340     
00341     Int_t nhits = pixHitCol->numberOfHits();
00342     LOG_DEBUG<<"There are "<<nhits<<" pixel hits"<<endm;
00343     if (nhits){
00344       Int_t id = col->numberOfHits();
00345       if(g2tPix){
00346         for(int jj=0;jj<g2t_pix_hit->GetTableSize();jj++)
00347           {
00348             rad  = TMath::Sqrt(g2tPix[jj].x[0]*g2tPix[jj].x[0] + g2tPix[jj].x[1]*g2tPix[jj].x[1]);  
00349             if((rad>2.2) && (rad<2.7)) {
00350               counter_layer1++;
00351               LOG_DEBUG <<"radius="<<rad<<" id="<<g2tPix[jj].id <<" x="<<g2tPix[jj].x[0]<<" y="<<g2tPix[jj].x[1]<<" z="<<g2tPix[jj].x[2]<<" truth is "<<g2tPix[jj].track_p<<" dE= "<<g2tPix[jj].de <<" counter layer 1 =" << counter_layer1 << endm;
00352             }
00353             else if((rad>7.9) && (rad<8.2)) {
00354               counter_layer2++;
00355               LOG_DEBUG <<"radius="<<rad<<" id="<<g2tPix[jj].id <<" x="<<g2tPix[jj].x[0]<<" y="<<g2tPix[jj].x[1]<<" z="<<g2tPix[jj].x[2]<<" truth is "<<g2tPix[jj].track_p<<" dE= "<<g2tPix[jj].de <<" counter layer 2 =" << counter_layer2 <<endm;
00356             }
00357           }
00358       }
00359       LOG_DEBUG << " there is " << counter_layer1 << " hits in the first layer " << endm;
00360       LOG_DEBUG << " there is " << counter_layer2 << " hits in the second layer " << endm;
00361       //for(int k=0;k<g2t_pix_hit->GetNRows();k++){
00362       for (UInt_t k=0; k<pixHitCol->numberOfLayers(); k++){
00363         if (pixHitCol->layer(k)){
00364           LOG_DEBUG<<"Layer "<<k+1<<endm;
00365           //simple simulator for perfect hits that just converts StMcPixelHit to StRnDHit
00366           //as of 11/21/08, hits are now smeared with resolution taken from hit error table
00367           UInt_t nh = pixHitCol->layer(k)->hits().size();
00368           LOG_DEBUG << " Number of hits in layer "<< k+1 <<" =" << nh << endm;
00369           for (UInt_t i = 0; i < nh; i++){
00370             counter++;
00371             int vid=g2tPix[k].volume_id;
00372             int layer=vid/1000000;
00373             int ladder=(vid%1000000)/10000;
00374             StMcHit *mcH = pixHitCol->layer(k)->hits()[i];
00375             StMcPixelHit* mcPix=dynamic_cast<StMcPixelHit*>(mcH);
00376             TString Path("");
00377             if(k==0) Path= Form("/HALL_1/CAVE_1/PXMO_1/PXLA_1/PLMI_%i/PLAC_1",mcPix->ladder());
00378             else Path=Form("/HALL_1/CAVE_1/PXMO_1/PXL1_2/PLM1_%i/PLA1_1",mcPix->ladder());
00379             //LOG_DEBUG<<"mc pixel hit location x: "<<g2tPix[k].x[0]<<"; y: "<<g2tPix[k].x[1]<<"; z: "<<g2tPix[k].x[2]<<endm;
00380             LOG_DEBUG<<"pixel hit layer/ladder is "<<layer<<"/"<<ladder<<" and volume id "<<vid<<endm;
00381             gGeoManager->RestoreMasterVolume();
00382             gGeoManager->CdTop();
00383             gGeoManager->cd(Path);
00384             //double globalPixHitPos[3]={g2tPix[k].x[0],g2tPix[k].x[1],g2tPix[k].x[2]};
00385             double globalPixHitPos[3]={mcPix->position().x(),mcPix->position().y(),mcPix->position().z()};
00386             double localPixHitPos[3]={0,0,0};
00387             gGeoManager->GetCurrentMatrix()->MasterToLocal(globalPixHitPos,localPixHitPos);
00388             smearedX=distortHit(localPixHitPos[0],resXPix,100);
00389             smearedZ=distortHit(localPixHitPos[2],resZPix,100);
00390             localPixHitPos[0]=smearedX;
00391             localPixHitPos[2]=smearedZ;
00392             double smearedGlobalPixHitPos[3]={0,0,0};
00393             gGeoManager->GetCurrentMatrix()->LocalToMaster(localPixHitPos,smearedGlobalPixHitPos);
00394             StThreeVectorF gpixpos(smearedGlobalPixHitPos);
00395             StThreeVectorD mRndHitError(0.,0.,0.);
00396             //StRnDHit* tempHit = new StRnDHit(mcPix->position(), mRndHitError, 1, 1., 0, 1, 1, id++, kPxlId);
00397             StRnDHit* tempHit = new StRnDHit(gpixpos, mRndHitError, 1, 1., 0, 1, 1, id++, kPxlId);
00398             tempHit->setVolumeId(mcPix->volumeId());
00399             tempHit->setLayer(k+1);
00400             tempHit->setLadder(mcPix->ladder());
00401             tempHit->setKey(mcPix->key());
00402             Int_t truth =0;
00403             if((k==0)&&(counter<=counter_layer1)) {
00404               truth = g2tPix[mcPix->key()-1].track_p;
00405             }
00406             else if((k==1)&&(counter<=(counter_layer2+pixHitCol->layer(0)->hits().size())))
00407               {
00408                 truth = g2tPix[mcPix->key()-1].track_p;
00409               }
00410             tempHit->setIdTruth(truth,100);
00411             LOG_DEBUG<<"key() : "<< mcPix->key()-1 << " idTruth: "<< truth <<endm;
00412             //LOG_DEBUG <<"from g2t : x= " << g2tPix[k].x[0] <<"  y= " << g2tPix[k].x[1] <<"  z= " << g2tPix[k].x[2] << endm;
00413             LOG_DEBUG<<"pixel rnd hit location x: "<<tempHit->position().x()<<"; y: "<<tempHit->position().y()<<"; z: "<<tempHit->position().z()<<endm;
00414             //LOG_DEBUG<<"pixel rnd hit location x: "<<tempHit->position().x()<<"; y: "<<tempHit->position().y()<<"; z: "<<tempHit->position().z()<<endm;
00415             col->addHit(tempHit);
00416           }
00417         }
00418       }
00419     }
00420     gMessMgr->Info() <<"StPixelFastSimMaker::Make() -I- Loaded "<<nhits<<"pixel hits. \n";
00421   }
00422   else{
00423     gMessMgr->Info() <<"No pixel hits found.\n";
00424   }
00425   
00426   
00427   const StMcIstHitCollection* istHitCol = mcEvent->istHitCollection();
00428   
00429   //new simulator for new 1-layer design
00430   int id=0;
00431   //float smearedX,smearedZ;
00432   TString Path("");
00433   if(istHitCol){
00434     LOG_INFO<<"ist hit collection found"<<endm;
00435     int nIsthits=istHitCol->numberOfHits();
00436     LOG_DEBUG<<"there are "<<nIsthits<<" ist hits"<<endm;
00437     vector<StMcIstHit*> ladderHits;
00438     if(nIsthits){
00439       if(istHitCol->layer(0)){
00440         //simple simulator for perfect hits that merely converts from StMcIstHit to StRnDHit
00441         //as of 11/21/08, the simulator now smears hits by hit resolutions taken from hit error tables
00442         for(unsigned int kk=0;kk<istHitCol->layer(0)->hits().size();kk++){
00443           StMcHit* mcH=istHitCol->layer(0)->hits()[kk];
00444           StMcIstHit* mcI=dynamic_cast<StMcIstHit*>(mcH);
00445           LOG_DEBUG<<"mc ist hit location x: "<<mcI->position().x()<<"; y: "<<mcI->position().y()<<"; z: "<<mcI->position().z()<<endm;
00446           TString Path("");
00447           Path=Form("/HALL_1/CAVE_1/IBMO_1/IBMY_1/IBAM_%i/IBLM_%i/IBSS_1",mcI->ladder(),mcI->wafer());
00448           gGeoManager->RestoreMasterVolume();
00449           gGeoManager->CdTop();
00450           gGeoManager->cd(Path);
00451           double globalIstHitPos[3]={mcI->position().x(),mcI->position().y(),mcI->position().z()};
00452           double localIstHitPos[3]={0,0,0};
00453           gGeoManager->GetCurrentMatrix()->MasterToLocal(globalIstHitPos,localIstHitPos);
00454           smearedX=distortHit(localIstHitPos[0],resXIst1,100);
00455           smearedZ=distortHit(localIstHitPos[2],resZIst1,100);
00456           localIstHitPos[0]=smearedX;
00457           localIstHitPos[2]=smearedZ;
00458           double smearedGlobalIstHitPos[3]={0,0,0};
00459           gGeoManager->GetCurrentMatrix()->LocalToMaster(localIstHitPos,smearedGlobalIstHitPos);
00460           StThreeVectorF gistpos(smearedGlobalIstHitPos);
00461           StRnDHit* tempHit = new StRnDHit(gistpos, mHitError, 1, 1., 0, 1, 1, id++, kIstId);  
00462           tempHit->setDetectorId(kIstId); 
00463           tempHit->setVolumeId(mcI->volumeId());
00464           tempHit->setKey(mcI->key());
00465           tempHit->setLayer(1);           
00466           tempHit->setLadder(mcI->ladder());           
00467           tempHit->setWafer(mcI->wafer());
00468           tempHit->setIdTruth(g2tIst[kk].track_p,100);
00469           LOG_DEBUG<<"ist hit volume id: "<<tempHit->volumeId()<<endm;
00470           LOG_DEBUG<<"ist hit ladder/wafer: "<<tempHit->ladder()<<"/"<<tempHit->wafer()<<endm;
00471           col->addHit(tempHit);
00472           LOG_DEBUG<<"ist rnd hit location x: "<<tempHit->position().x()<<"; y: "<<tempHit->position().y()<<"; z: "<<tempHit->position().z()<<" idTruth = " << g2tIst[kk].track_p << endm;
00473         }
00474       }
00475     }
00476   }
00477   else{
00478     LOG_INFO <<"No Ist hits found."<<endm;
00479   }
00480   
00481   
00482   const StMcFgtHitCollection* fgtHitCol = mcEvent->fgtHitCollection();
00483   if (fgtHitCol)
00484     {
00485       Int_t nhits = fgtHitCol->numberOfHits();
00486       if (nhits)
00487         {
00488           Int_t id = 0;
00489           //StSPtrVecHit *cont = new StSPtrVecHit();
00490           //rcEvent->addHitCollection(cont, # Name );
00491           for (UInt_t k=0; k<fgtHitCol->numberOfLayers(); k++){
00492             if (fgtHitCol->layer(k))
00493               {
00494                 UInt_t nh = fgtHitCol->layer(k)->hits().size();
00495                 for (UInt_t i = 0; i < nh; i++) {
00496                   StMcHit *mcH = fgtHitCol->layer(k)->hits()[i];
00497                   StRnDHit* tempHit = new StRnDHit(mcH->position(), mHitError, 1, 1., 0, 1, 1, id++);
00498                   tempHit->setVolumeId(mcH->volumeId());
00499                   tempHit->setKey(mcH->key());
00500 
00501                   StMcIstHit *mcI = dynamic_cast<StMcIstHit*>(mcH);
00502                   if(mcI){
00503                     tempHit->setLayer(mcI->layer());
00504                     tempHit->setLadder(mcI->ladder());
00505                     tempHit->setWafer(mcI->wafer());
00506                     //tempHit->setExtraByte0(mcI->side());
00507                   }
00508                   col->addHit(tempHit);
00509                 }
00510               }
00511           }
00512         }
00513     }
00514   
00515   rcEvent->setRnDHitCollection(col);
00516   return kStOK;
00517     }
00518 //________________________________________________________________________________
00519 Bool_t StPixelFastSimMaker::accept(StEvent* event){
00520   return event ? true : false;
00521 }
00522 //________________________________________________________________________________
00523 Bool_t StPixelFastSimMaker::accept(StMcEvent* event){
00524   return event ? true : false;
00525 }
00526 
00527 double StPixelFastSimMaker::distortHit(double x, double res, double detLength){
00528   double test;
00529   if(mSmear){
00530     test = x + myRandom->gauss(0,res);
00531     while( fabs(test) > detLength){
00532       test = x + myRandom->gauss(0,res);
00533     }
00534     //cout << " x was " <<x<< " and is now " << test<< endl;
00535     return test;
00536   }
00537   else return x;
00538 }
00539 //____________________________________________________________
00540 
00541 
00542 //____________________________________________________________
00543 void StPixelFastSimMaker::smearGaus(StThreeVectorD &mError, 
00544                                     double sigma1, double sigma2)
00545 {
00546 
00547   // smear hit in transverse plane, 
00548   // sigma's are in microns
00549   double u1=-1;
00550   double u2=-1.; 
00551   double v1=-1.;
00552   double v2=-1.;;
00553   double r = 2.;
00554   double z1 = 10.;
00555   double z2 = 10.;
00556   while(fabs(z1)>2. || fabs(z2)>2.) // sigma
00557     {
00558       r = 2.;
00559       while(r>1.)
00560         {
00561           u1 = rand()/double(RAND_MAX); 
00562           u2 = rand()/double(RAND_MAX);
00563           v1 = 2*u1 - 1.;  
00564           v2 = 2*u2 - 1.; 
00565           r = pow(v1,2) + pow(v2,2);
00566         }
00567       z1 = v1*sqrt(-2.*log(r)/r);  z2 = v2*sqrt(-2.*log(r)/r);
00568     }
00569 
00570   //set to be cumulative with other transforms 
00571   mError.setX(mError.x()+z1*sigma1/1.e04);
00572   mError.setY(mError.y()+z1*sigma1/1.e04);
00573   mError.setZ(mError.z()+z2*sigma2/1.e04);
00574 }
00575 
00576 int StPixelFastSimMaker::sector(int layer, int ladder)
00577 {
00578  
00579   if (layer ==1)
00580     {
00581       if ( ladder < 4 ) return 1;
00582       if ( ladder < 7 ) return 2;
00583       return 3;
00584     }
00585   else
00586     {
00587       if ( ladder < 9 ) return 1;
00588       if ( ladder < 17 ) return 2;
00589       return 3;
00590     }
00591 
00592   
00593 }
00594 
00595 int StPixelFastSimMaker::secLadder(int layer, int ladder)
00596 {
00597   if (layer==1)
00598     return ladder - 3*(sector(layer,ladder)-1);
00599   else
00600     return ladder - 8*(sector(layer,ladder)-1)+3;
00601 }
00602 
00603 
00604 
00605 double StPixelFastSimMaker::phiForLadder(int layer, int ladder)
00606 {
00607   int sec = sector(layer,ladder);
00608   int secLad = secLadder(layer,ladder);
00609   //double phi=0.;
00610   double secPhi=0.;
00611   double ladPhi=0.;
00612   switch (sec)
00613     {
00614     case 1:
00615       secPhi=0.;
00616       break;
00617     case 2:
00618       secPhi=120.;
00619       break;
00620     case 3:
00621       secPhi=240.;
00622       break;
00623     }
00624 
00625   switch (secLad)
00626     {
00627     case 1:
00628       ladPhi=100.;
00629       break;
00630     case 2:
00631       ladPhi=60.;
00632       break;
00633     case 3:
00634       ladPhi=20.;
00635       break;
00636     case 4:
00637       ladPhi=105.;
00638       break;
00639     case 5:
00640       ladPhi=90.;
00641       break;
00642     case 6:
00643       ladPhi=75.;
00644       break;
00645     case 7:
00646       ladPhi=60.;
00647       break;
00648     case 8:
00649       ladPhi=45.;
00650       break;
00651     case 9:
00652       ladPhi=30.;
00653       break;
00654     case 10:
00655       ladPhi=15.;
00656       break;
00657     case 11:
00658       ladPhi=0.;
00659       break;
00660     }
00661 
00662   return secPhi+ladPhi;
00663 }
00664 
00665  
00666 void StPixelFastSimMaker::shiftHit(StThreeVectorF &position,StThreeVectorF &mom, int layer, int ladder)
00667 {
00668 
00669   //printf("Entering hit shift code. %i %i\n",
00670   // sector(layer,ladder), secLadder(layer,ladder));
00671   LOG_DEBUG<<"Pixel Hit in layer "<<layer<<" and ladder "<<ladder<<" or sector "<<sector(layer,ladder)<<" and sector ladder "<<secLadder(layer,ladder)<<endm;
00672 
00673   TString Path("");
00674   Path = Form("/HALL_1/CAVE_1/PXMO_1/PSEC_%i/PLMO_%i/PLAC_1",
00675               sector(layer,ladder),secLadder(layer,ladder));
00676   gGeoManager->RestoreMasterVolume();
00677   //printf("Master volume.\n");
00678  
00679   gGeoManager->CdTop();
00680   gGeoManager->cd(Path);
00681   TGeoPhysicalNode* node= (TGeoPhysicalNode*)(gGeoManager->GetCurrentNode());
00682   if (!node ) LOG_ERROR<< "Failed to get node for sector(layer,ladder), secLadder(layer, ladder)"<<endm;
00683  
00684   double pos[3]={position.x(),position.y(),position.z()};
00685   double localpos[3]={0,0,0};
00686   gGeoManager->GetCurrentMatrix()->MasterToLocal(pos,localpos);
00687   //printf("old hit: %g %g %g\n",localpos[0],localpos[1],localpos[2]);
00688  
00689   //get ladder phi 
00690   //TGeoHMatrix  *hmat   = (TGeoHMatrix*)(node->GetMatrix());
00691   TGeoHMatrix  *hmat   =  gGeoManager->GetCurrentMatrix();
00692  
00693   if (! hmat )
00694     {
00695       LOG_ERROR<< "Can't shift hit - no hmat."<<endm;
00696     }
00697   
00698   Double_t     *rot    = hmat->GetRotationMatrix();
00699   if (! rot )
00700     {
00701       LOG_ERROR <<"Can't shift hit - no rotation matrix."<<endm;
00702     }
00703  
00704   StThreeVectorD normalVector(rot[1],rot[4],rot[7]);
00705   //double momentum = mom.magnitude();
00706   StThreeVectorF momUnit(mom);
00707   momUnit/=momUnit.magnitude();
00708  
00709   //momentum mag, pt and z stay the same; angle changes 
00710   momUnit.setPhi( momUnit.phi() - phiForLadder(layer,ladder)*3.141592654/180.);
00711   //printf("Mom dir: %g\n",momUnit.phi());
00712   
00713   // shift = x + y * tan(phi)
00714   //.006 is the half thickness of the active area. Hardcoded; not good; blah, blah.
00715   double dx = (.006)*(momUnit.y()/momUnit.x());
00716   double dz = (.006)*(momUnit.y()/momUnit.z());
00717   localpos[0] = localpos[0] - dx;
00718   localpos[2] = localpos[2] - dz;
00719   localpos[1] = localpos[1] - .006; //this isn't exactly right. The local position is just off of radius -.006, but it's close (~1-5 um).
00720    
00721   gGeoManager->GetCurrentMatrix()->LocalToMaster(localpos,pos);
00722 
00723   //printf("exit shift code.\n");
00724 
00725 }