StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StIstSlowSimMaker.cxx
1 /***************************************************************************
2  *
3  * $Id: StIstSlowSimMaker.cxx,v 1.6 2018/04/18 18:10:19 jwebb Exp $
4  *
5  * Author: Leszek Kosarzewski, March 2014
6  ****************************************************************************
7  * Description:
8  * See header file.
9  ****************************************************************************
10  * StIstSlowSimMaker.h,v 1.1
11  * Revision 1.1 2014/08/05 10:54:12 ypwang
12  * Update to integrated into IST offline chain
13  * 1) process IST MC hits, and discretize to individual IST channel
14  * 2) translate to ADC by simple GeV-2-ADC factor, raw ADC output
15  ****************************************************************************
16  * StIstSlowSimMaker.cxx,v 1.0
17  * Revision 1.0 2014/03/07 11:25:30 lkosarz
18  * Initial version
19  ****************************************************************************/
20 
21 //#include <Stiostream.h>
22 //#include <StHit.h>
23 //#include <StEventTypes.h>
24 //#include <StEvent.h>
25 #include <StMcEvent.hh>
26 #include <StMcEvent/StMcHit.hh>
27 #include <StMcIstHit.hh>
28 //#include <StMcEvent/StMcIstHit.hh>
29 #include <StMcEvent/StMcIstHitCollection.hh>
30 //#include <StMcEventTypes.hh>
31 
32 //#include <stdio.h>
33 //#include <map>
34 //#include <exception>
35 using namespace std;
36 //#include <stdexcept>
37 #include "tables/St_g2t_ist_hit_Table.h"
38 #include "tables/St_HitError_Table.h"
39 #include "TGeoManager.h"
40 //#include "TGeoMatrix.h"
41 //#include "TDataSet.h"
42 #include "TF1.h"
43 
44 #include "StIstDbMaker/StIstDb.h"
45 #include "StIstUtil/StIstCollection.h"
46 #include "StIstUtil/StIstRawHitCollection.h"
47 #include "StIstUtil/StIstRawHit.h"
48 #include "StIstUtil/StIstConsts.h"
49 #include "StIstSlowSimMaker.h"
50 #include "TRandom3.h"
51 #include "tables/St_istSimPar_Table.h"
52 
53 #include "tables/St_istMapping_Table.h"
54 #include "tables/St_istControl_Table.h"
55 
56 #include <SystemOfUnits.h>
57 
58 ClassImp(StIstSlowSimMaker)
59 
60 StIstSlowSimMaker::StIstSlowSimMaker(const char* name): StMaker(name), mIstDb(NULL), mBuildIdealGeom(kFALSE), mIstCollectionPtr(NULL), mHitEffMode(0), mMomCut(0.), mHitEff(1.0), mRndGen(nullptr)
61 {
62  mMappingGeomVec.resize( kIstNumElecIds );
63  mDefaultTimeBin = 9;
64  mCurrentTimeBinNum = 9;
65  fAdc = new TF1("fAdc","landau",0,9);
66  fAdc->SetParameters(1.14288,3.8,1.168);
67 }
68 
69 Int_t StIstSlowSimMaker::Init()
70 {
71  LOG_DEBUG<<"StIstSlowSimMaker::Init()"<<endm;
72 
73  //prepare output dataset
74  mIstCollectionPtr = new StIstCollection();
75  ToWhiteConst("istRawAdcSimu",mIstCollectionPtr);
76 
77  if(!mIstCollectionPtr ) {
78  LOG_WARN << "Error constructing istCollection" << endm;
79  return kStErr;
80  }
81 
82  return kStOk;
83 }
84 
85 //____________________________________________________________
86 Int_t StIstSlowSimMaker::InitRun(Int_t runnumber)
87 {
88  LOG_DEBUG<<"StIstSlowSimMaker::InitRun "<<runnumber<<endm;
89 
90  Int_t ierr = kStOk;
91 
92  TObjectSet *istDbDataSet = (TObjectSet *)GetDataSet("ist_db");
93  if (istDbDataSet) {
94  mIstDb = (StIstDb *)istDbDataSet->GetObject();
95  assert(mIstDb);
96  }
97  else {
98  LOG_ERROR << "InitRun : no istDb" << endm;
99  return kStErr;
100  }
101 
102  // IST control parameters
103  const istControl_st *istControlTable = mIstDb->getControl() ;
104  if (!istControlTable) {
105  LOG_ERROR << "Pointer to IST control table is null" << endm;
106  ierr = kStErr;
107  }
108  else {
109  mDefaultTimeBin = istControlTable[0].kIstDefaultTimeBin;
110  mCurrentTimeBinNum = istControlTable[0].kIstCurrentTimeBinNum;
111  }
112 
113  // IST mapping table
114  const istMapping_st *gM = mIstDb->getMapping();
115  if( !gM ) {
116  LOG_ERROR << "Pointer to IST mapping table is null" << endm;
117  ierr = kStErr;
118  }
119  else {
120  for(Int_t i=0; i<kIstNumElecIds; i++) {
121  LOG_DEBUG<<Form(" PrInt_t entry %d : geoId=%d ",i,gM[0].mapping[i])<<endm;
122  Int_t geomIdTemp = gM[0].mapping[i];
123  mMappingGeomVec[geomIdTemp-1] = i;
124  }
125  }
126 
128  mRndGen = (TRandom3 *)gRandom;
129  // MC->RC hit efficiency (momentum dependence - default)
130  istSimPar_st const* const istSimParTable = mIstDb->istSimPar();
131  mHitEffMode = istSimParTable[0].mode;
132  mMomCut = istSimParTable[0].pCut;
133  mHitEff = istSimParTable[0].effIst; // best knowledge from ZF cosmic ray study - tunable
134 
135  LOG_INFO << " IST MC hit efficiency mode used for IST slow simulator: " << mHitEffMode << endm;
136  LOG_INFO << " +++ Hit Efficiency at p > " << mMomCut << " GeV/c = " << mHitEff << endm;
137 
138  return ierr;
139 }
140 
141 //______________________________________________________________________________
142 Int_t StIstSlowSimMaker::Make()
143 {
144  // Get the input data structures from StMcEvent
145  StMcEvent* mcEvent = (StMcEvent *) GetInputDS("StMcEvent");
146  if (! mcEvent) {LOG_INFO << "No StMcEvent on input" << endl; return kStWarn;}
147 
148  TDataSetIter geant(GetInputDS("geant"));
149  if ( mBuildIdealGeom && !gGeoManager ) {
150  GetDataBase("VmcGeometry");
151  }
152 
153  StThreeVectorF mHitError(0.,0.,0.);
154 
155  // Get MC Ist hit collection. This contains all ist hits.
156  const StMcIstHitCollection* istMcHitCol = mcEvent->istHitCollection();
157 
158  mIstCollectionPtr->setNumTimeBins(mCurrentTimeBinNum);
159 
160  // Calculate pad center position from mc hits. Convert dE to ADC value. Save hits to raw hit collection.
161  if(istMcHitCol){
162  LOG_DEBUG<<"ist MC hit collection found"<<endm;
163  Int_t nIsthits=istMcHitCol->numberOfHits();
164  if(nIsthits){
165  if(istMcHitCol->layer(0)){
166  StSPtrVecMcIstHit mcHitVec = istMcHitCol->layer(0)->hits();
167  for(std::vector<StMcIstHit*>::iterator mcHitIt = mcHitVec.begin();mcHitIt!=mcHitVec.end(); ++mcHitIt){
168  LOG_DEBUG << "IST MC hit found ...... " << endm;
169  StMcIstHit* mcIstHit = (*mcHitIt);
170  if(!mcIstHit) continue;
171  float hitEff = 1.0;
172  switch (mHitEffMode) {
173  case 0: // ideal case 100% efficiency
174  break;
175  case 1:
176  if(mcIstHit->parentTrack()) {
177  float const ptot = mcIstHit->parentTrack()->momentum().mag();
178  hitEff = 1.0 - (1. - mHitEff)*ptot/mMomCut;
179  if(hitEff<mHitEff) hitEff = mHitEff;
180  }
181  break;
182  case 2:
183  hitEff = mHitEff;
184  break;
185  default:
186  break;
187  }
188 
189  if( mRndGen->Rndm()>hitEff ) continue;
190 
191  //mcIstHits stored local position
192  Double_t localIstHitPos[3]={mcIstHit->position().x(),mcIstHit->position().y(),mcIstHit->position().z()};
193  LOG_DEBUG << "StMcIstHit local position (before discretization) = " << localIstHitPos[0] << " " << localIstHitPos[1] << " " << localIstHitPos[2] << endm;
194  Int_t ladder = mcIstHit->ladder(); //coded from 1 to 24
195  Int_t sensor = mcIstHit->wafer(); //coded from 1 to 6
196  UShort_t meanColumn, meanRow;
197 
198  //get column and row value from mc hit local position
199  getMCHitRowAndColumn(mcIstHit, meanColumn, meanRow);
200 
201  //discretize hit local position (2D structure of IST sensor pads)
202  Float_t rPhiPos_mean = (meanRow-1) * kIstPadPitchRow + 0.5 * kIstPadPitchRow; //unit: cm
203  Float_t zPos_mean = (meanColumn-1) * kIstPadPitchColumn + 0.5 * kIstPadPitchColumn; //unit: cm
204  localIstHitPos[0] = kIstSensorActiveSizeRPhi/2.0 - rPhiPos_mean;
205  localIstHitPos[2] = zPos_mean - kIstSensorActiveSizeZ/2.0;
206 
207  LOG_DEBUG << "geometry position ladder = " << ladder << "\tsensor = " << sensor << "\tcolumn = " << meanColumn << "\trow = " << meanRow << endm;
208  LOG_DEBUG << "StMcIstHit local position (after discretization) = " << localIstHitPos[0] << " " << localIstHitPos[1] << " " << localIstHitPos[2] << endm;
209 
210  StMcTrack *trackMC = mcIstHit->parentTrack();
211  Long_t pdgId = trackMC->pdgId();
212  Long_t idTruth = trackMC->key();
213  LOG_DEBUG << "trackMC pdgId = " << pdgId << "\tidTruth = " << idTruth << endm;
214 
215  //
216  generateRawHits(mcIstHit);
217  }//MC hits loop over
218  }//end layer=0 cut
219  }//end MC hits number cut
220  LOG_DEBUG << "StIstSlowSimMaker::Make() -I- Loaded " << nIsthits << " IST MC raw hits. \n";
221  }
222  else{
223  LOG_DEBUG <<"No Ist MC hits found."<<endm;
224  }
225 
226  LOG_DEBUG << "Number of IST MC raw hits produced by slow simulator: " << mIstCollectionPtr->getNumRawHits() << endm;
227 
228  return kStOK;
229 }
230 
231 void StIstSlowSimMaker::Clear( Option_t *opts )
232 {
233  if(mIstCollectionPtr ) {
234  for ( UChar_t i = 0; i < kIstNumLadders; ++i ) {
235  mIstCollectionPtr->getRawHitCollection(i)->Clear(opts);
236  }
237  }
238  return StMaker::Clear();
239 }
240 
246 void StIstSlowSimMaker::getMCHitRowAndColumn(const StMcIstHit *istMChit, UShort_t &meanColumn, UShort_t &meanRow) const
247 {
248  Double_t localIstHitPos[3]={istMChit->position().x(),istMChit->position().y(),istMChit->position().z()};
249  Float_t rPhiPos = kIstSensorActiveSizeRPhi/2.0 - localIstHitPos[0];
250  Float_t zPos = localIstHitPos[2] + kIstSensorActiveSizeZ/2.0;
251  if(rPhiPos<0||zPos<0){
252  LOG_WARN<<"StIstSlowSimMaker::getMCHitRowAndColumn Wrong local position rPhiPos = "<<rPhiPos<<" zPos = "<<zPos<<endm;
253  }
254  meanColumn = (UShort_t)floor( zPos/kIstPadPitchColumn ) + 1;
255  meanRow = (UShort_t)floor( rPhiPos/kIstPadPitchRow ) + 1;
256 
257 }
258 
259 void StIstSlowSimMaker::generateRawHits(const StMcIstHit *istMChit) const
260 {
261  LOG_DEBUG << "StIstSlowSimMaker::generateRawHits - Calculating pad row and column ... " << endm;
262  Int_t ladderId = istMChit->ladder();
263  Int_t sensorId = istMChit->wafer();
264  UShort_t meanColumn, meanRow;
265  getMCHitRowAndColumn(istMChit, meanColumn, meanRow);
266 
267  Double_t dS = istMChit->dS(); //distance traveled in sensor active volume
268  StThreeVectorD midPos(istMChit->position().x(),istMChit->position().y(),istMChit->position().z());
269  Double_t mcLocPx = istMChit->localMomentum().x();
270  Double_t mcLocPy = istMChit->localMomentum().y();
271  Double_t mcLocPz = istMChit->localMomentum().z();
272 
273  StThreeVectorD mcLocP(mcLocPx,mcLocPy,mcLocPz);
274  Double_t mcLocTot = mcLocP.mag();
275  StThreeVectorD mcLocalDir = mcLocP.unit();
276 
277  StThreeVectorD inPos = midPos - mcLocalDir * dS / 2.0;
278  StThreeVectorD outPos = midPos + mcLocalDir * dS / 2.0;
279 
280  LOG_DEBUG<<"istMcHit geometry position ladder = " << ladderId << "\tsensor = " << sensorId << "\tcolumn = " << meanColumn << "\trow = " << meanRow << endm;
281  LOG_DEBUG<<"istMcHit local position x = "<<midPos.x()<< "\ty = "<<midPos.y()<<"\tz = "<<midPos.z()<<endm;
282  LOG_DEBUG<<"istMcHit local position enter x = "<<inPos.x()<< "\ty = "<<inPos.y()<<"\tz = "<<inPos.z()<<endm;
283  LOG_DEBUG<<"istMcHit local position exit x = "<<outPos.x()<< "\ty = "<<outPos.y()<<"\tz = "<<outPos.z()<<endm;
284  LOG_DEBUG<<"istMcHit local momentum x = "<<mcLocP.x()<< "\ty = "<<mcLocP.y()<<"\tz = "<<mcLocP.z()<<endm;
285  LOG_DEBUG<<"istMcHit local direction x = "<<mcLocalDir.x()<< "\ty = "<<mcLocalDir.y()<<"\tz = "<<mcLocalDir.z()<<endm;
286  LOG_DEBUG<<"istMcHit dS = "<<dS<<"\ttotMom = "<<mcLocTot<<endl;
287 
288  vector<StThreeVectorD> crossVec;
289  LOG_DEBUG<<"transforming to sensor coordinates"<<endm;
290  transformToSensor(inPos);
291  transformToSensor(outPos);
292  LOG_DEBUG<<"inPos x = "<<inPos.x()<<"\ty = "<<inPos.y()<<"\tz = "<<inPos.z()<<endm;
293  LOG_DEBUG<<"outPos x = "<<outPos.x()<<"\ty = "<<outPos.y()<<"\tz = "<<outPos.z()<<endm;
294 
295  checkPadCrossing(inPos, outPos, mcLocalDir, dS, crossVec);
296  LOG_DEBUG<<"StIstSlowSimMaker::generateRawHits crossVec.size() = "<<crossVec.size()<<endm;
297  if(crossVec.size()==1) {
298  LOG_WARN<<"StIstSlowSimMaker: McHit is outside of active area -> skip!"<<endm;
299  return;
300  }
301 
302  StThreeVectorD totalPath = crossVec[crossVec.size()-1]-crossVec[0];
303  Double_t pathLengthTotal = totalPath.mag();
304  LOG_DEBUG<<"inPos x = "<<crossVec[0].x()<<"\ty = "<<crossVec[0].y()<<"\tz = "<<crossVec[0].z()<<endm;
305 
306  //Energy-to-ADC translation factor + signal time bin dependence factor
307  //silicon dE/dx = 3.87 Mev/cm, thickness of IST sensor = 300 micros, ADC MIP = 441.0 ADC counts
308  const Double_t sidEdx = 3.87; //MeV/cm
309  const Double_t sensorThickness = 0.03; //cm
310  const Double_t mip = 441.0; //MIP adc
311  const Double_t kIstMPV = mip / (sidEdx * 1e-3 * sensorThickness);
312  Double_t kIstTimeBinFrac[mCurrentTimeBinNum];
313 
314  Int_t maxTB = mCurrentTimeBinNum/2;
315  Double_t mean = fAdc->GetParameter(1);
316  Double_t sum = fAdc->Integral(mean-0.5-maxTB,mean-0.5+(mCurrentTimeBinNum-maxTB));
317  for(Int_t iTB=0; iTB<mCurrentTimeBinNum; iTB++){
318  if(fabs(sum)>1e-6){
319  kIstTimeBinFrac[iTB] = fAdc->Integral(mean-0.5-maxTB+iTB,mean-0.5-maxTB+iTB+1)/sum; //time-bin dependence
320  }else{
321  LOG_WARN<<"StIstSlowSimMaker: integral of adc spectrum = 0"<<endm;
322  }
323  }
324 
325  StIstRawHitCollection *istRawHitCollection = mIstCollectionPtr->getRawHitCollection(ladderId-1);
326  if(istRawHitCollection) {
327  for (UInt_t i = 0; i < crossVec.size()-1; ++i) {
328  LOG_DEBUG<<"i = "<<i+1<<"\tpos x = "<<crossVec[i+1].x()<<"\ty = "<<crossVec[i+1].y()<<"\tz = "<<crossVec[i+1].z()<<endm;
329  StThreeVectorD path = crossVec[i+1]-crossVec[i];
330  StThreeVectorD meanPos = (crossVec[i+1]+crossVec[i])/2.0;
331 
332  Double_t pathLength = path.mag();
333  Double_t rPhiPos_mean, zPos_mean; //unit: cm
334  findPad(meanPos, meanColumn, meanRow, rPhiPos_mean, zPos_mean);
335 
336  //Nrow = 64 = kIstNumApvChannels/2, Ncolumn = 12 = kIstNumApvsPerArm/2
337  //One APV reads 64 rows x 12 columns
338  Int_t geoId = (ladderId-1)*kIstNumApvChannels*kIstApvsPerLadder+ (sensorId-1)*(kIstNumApvsPerArm/2)*(kIstNumApvChannels/2) + (meanColumn-1)*kIstNumApvChannels/2 + meanRow;
339  Int_t elecId = mMappingGeomVec[geoId-1];
340  LOG_DEBUG << "elecID = " << elecId << "\tgeoId = " << geoId << endm;
341 
342  Int_t idTruth = istMChit->parentTrack()->key();
343 
344  StThreeVectorD normal(0,1,0);
345  Float_t angleCorr = fabs(mcLocalDir.dot(normal));
346  LOG_DEBUG<<"incident angle w.r.t. sensor surface normal direction: "<<angleCorr<<endm;
347 
348  StIstRawHit *rawHit = istRawHitCollection->getRawHit(elecId);
349 
350  // calculate ADC by dE for every timebin
351  for (UChar_t t = 0; t < mCurrentTimeBinNum; t++) {
352  Float_t adcSum = 0.0;
353  if(rawHit->getChannelId()>=0) {
354  adcSum += rawHit->getCharge(t);
355  }
356 
357  LOG_DEBUG<<"dE = "<<1e6*istMChit->dE()<<" keV \tpathLength = "<<pathLength<<"\tpathLengthTotal = "<<pathLengthTotal<<endm;
358  Float_t charge = 0;
359  if(pathLengthTotal>1e-6) charge = istMChit->dE() * pathLength / pathLengthTotal;
360  if(kIstTimeBinFrac[maxTB]>0) charge *= kIstTimeBinFrac[t]*kIstMPV/kIstTimeBinFrac[maxTB]; //translate energy Int_to ADC with GeV-to-ADC factor
361  if ( charge > adcSum ) {
362  rawHit->setIdTruth(idTruth);
363  }
364  adcSum += charge;
365 
366  rawHit->setCharge(adcSum, t);
367  LOG_DEBUG<<"charge = "<<adcSum<<"\tat TB"<<(Int_t)t<<endm;
368  }
369  rawHit->setChannelId( elecId );
370  rawHit->setGeoId( geoId );
371  rawHit->setDefaultTimeBin( mDefaultTimeBin );
372 
373  LOG_DEBUG << "rawHitPosition ladder = " <<(Int_t)rawHit->getLadder()<< "\tsensor = " <<(Int_t)rawHit->getSensor()<<"\tcolumn = "<<(Int_t)rawHit->getColumn()<<"\trow = "<<(Int_t)rawHit->getRow()<< endm;
374  }
375  }
376  else {
377  LOG_WARN << "StIstSlowSimMaker: No rawHitCollection found for ladder " << ladderId << endm;
378  }
379 }
380 
385 void StIstSlowSimMaker::checkPadCrossing(const StThreeVectorD inPos, const StThreeVectorD outPos, StThreeVectorD mcLocalDir, Double_t dS, vector<StThreeVectorD> &cross_vec) const
386 {
387  UShort_t column;
388  UShort_t row;
389  UShort_t column_out;
390  UShort_t row_out;
391  Double_t rPhiPosMean;
392  Double_t zPosMean;
393  Double_t rPhiPosMean_out;
394  Double_t zPosMean_out;
395 
396  cross_vec.clear();
397  cross_vec.push_back(inPos);
398  StThreeVectorD current = inPos;
399 
400  findPad(outPos, column_out, row_out, rPhiPosMean_out, zPosMean_out);
401  findPad(inPos, column, row, rPhiPosMean, zPosMean);
402 
403  LOG_DEBUG<<"StIstSlowSimMaker::checkPadCrossing"<<endm;
404  LOG_DEBUG<<"\tcolumn_out = "<<column_out<<"\trow_out = "<<row_out<<endm;
405  LOG_DEBUG<<"\tcolumn_in = "<<column<<"\trow_in = "<<row<<endm;
406 
407  if(column==65535||column_out==65535) return;
408 
409  mcLocalDir.setX(-mcLocalDir.x());
410 
411  Short_t row_dist = abs(row-row_out);
412  Short_t column_dist = abs(column-column_out);
413 
414  LOG_DEBUG<<"StIstSlowSimMaker::checkPadCrossing column_dist = "<<column_dist<<"\trow_dist = "<<row_dist<<endm;
415 
416  StThreeVectorD mean(rPhiPosMean, 0.0, zPosMean);
417 
418  unsigned short row_next = row;
419  unsigned short column_next = column;
420 
421  while(row_dist>=1 || column_dist>=1)
422  {
423  StThreeVectorD distance;
424 
425  StThreeVectorD halfPad(direction(mcLocalDir.x())*kIstPadPitchRow/2.0,
426  direction(mcLocalDir.y())*scaleFromYvsX(mcLocalDir, kIstPadPitchRow/2.0),
427  direction(mcLocalDir.z())*kIstPadPitchColumn/2.0);
428 
429  LOG_DEBUG<<"current x = "<<current.x()<<"\ty = "<<current.y()<<"\tz = "<<current.z()<<endm;
430 
431  distanceToBorder(current, mcLocalDir, mean, distance);
432  LOG_DEBUG<<"current x = "<<current.x()<<"\ty = "<<current.y()<<"\tz = "<<current.z()<<endm;
433 
434  current+=distance;
435  StThreeVectorD current_next = current+distance/100.0;
436  LOG_DEBUG<<"current_next x = "<<current_next.x()<<"\ty = "<<current_next.y()<<"\tz = "<<current_next.z()<<endm;
437 
438 // findPad(current_next, column, row, rPhiPosMean, zPosMean);
439  findPad(current_next, column_next, row_next, rPhiPosMean, zPosMean);
440 
441  if(row_next==row && column_next==column) {
442  LOG_WARN << " distance calculating not yielding to the next row/column! Break! " << endm;
443  break;
444  }
445 
446  row = row_next;
447  column = column_next;
448 
449  cross_vec.push_back(current);
450  mean = StThreeVectorD(rPhiPosMean, 0.0, zPosMean);
451  row_dist = abs(row-row_out);
452  column_dist = abs(column-column_out);
453  }
454 
455  cross_vec.push_back(outPos);
456 }
457 
462 void StIstSlowSimMaker::findPad(const StThreeVectorD hitPos, UShort_t &column, UShort_t &row, Double_t &rPhiPos_mean, Double_t &zPos_mean) const
463 {
464  Double_t rPhiPos = hitPos.x();
465  Double_t zPos = hitPos.z();
466  //allow +/- 0.5 pad at edge
467  if(rPhiPos<-kIstPadPitchRow/2.||rPhiPos>kIstSensorActiveSizeRPhi+kIstPadPitchRow/2.||zPos<-kIstPadPitchColumn/2.||zPos>kIstSensorActiveSizeZ+kIstPadPitchColumn/2.){
468  LOG_ERROR<<"StIstSlowSimMaker::findPad Wrong local position rPhiPos = "<<rPhiPos<<" zPos = "<<zPos<<endm;
469  column = row = rPhiPos_mean = zPos_mean = 65535;
470  return;
471  }
472  row = (UShort_t)floor( rPhiPos/kIstPadPitchRow ) + 1;
473  if(row<1) row = 1;
474  if(row>kIstNumRowsPerSensor) row = kIstNumRowsPerSensor;
475  column = (UShort_t)floor( zPos/kIstPadPitchColumn ) + 1;
476  if(column<1) column = 1;
477  if(column>kIstNumRowsPerSensor) column = kIstNumColumnsPerSensor;
478  rPhiPos_mean = (row-1) * kIstPadPitchRow + 0.5 * kIstPadPitchRow; //unit: cm
479  zPos_mean = (column-1) * kIstPadPitchColumn + 0.5 * kIstPadPitchColumn; //unit: cm
480 }
481 
482 void StIstSlowSimMaker::transformToSensor(StThreeVectorD &hitPos) const
483 {
484  hitPos.setX(kIstSensorActiveSizeRPhi/2.0 - hitPos.x());
485  hitPos.setZ(hitPos.z() + kIstSensorActiveSizeZ/2.0);
486 }
487 
488 Double_t StIstSlowSimMaker::distanceToBorder(const StThreeVectorD hitPos, const StThreeVectorD dir, const StThreeVectorD mean, StThreeVectorD &dist) const
489 {
490  StThreeVectorD halfPad(kIstPadPitchRow/2.0, 0.0150, kIstPadPitchColumn/2.0);
491  LOG_DEBUG<<"distanceToBorder() mean x = "<<mean.x()<<"\ty = "<<mean.y()<<"\tz = "<<mean.z()<<endm;
492 
493  // yz plane at x
494  Double_t plane_x = mean.x()+direction(dir.x())*halfPad.x();
495  Double_t diff_x = plane_x-hitPos.x();
496 
497  // xy plane at z
498  Double_t plane_z = mean.z()+direction(dir.z())*halfPad.z();
499  Double_t diff_z = plane_z-hitPos.z();
500 
501  LOG_DEBUG<<"yz plane x = "<<plane_x<<"\txy plane z = "<<plane_z<<endm;
502  LOG_DEBUG<<"diff_x = "<<diff_x<<"\tdiff_z = "<<diff_z<<endm;
503  LOG_DEBUG<<"dir x = "<<dir.x()<<"\ty = "<<dir.y()<<"\tz = "<<dir.z()<<endm;
504 
505  if(fabs(dir.x())<1e-6)
506  {
507  dist.setZ(diff_z);
508  dist.setY(scaleFromYvsZ(dir, diff_z));
509  return dist.mag(); // 3D distance
510  }
511  if(fabs(dir.z())<1e-6)
512  {
513  dist.setX(diff_x);
514  dist.setY(scaleFromYvsX(dir, diff_x));
515  return dist.mag(); // 3D distance
516  }
517  if(diff_x/dir.x()<=diff_z/dir.z())
518  {
519  dist.setX(diff_x);
520  dist.setZ(diff_x*dir.z()/dir.x());
521  }
522  else
523  {
524  dist.setX(diff_z*dir.x()/dir.z());
525  dist.setZ(diff_z);
526  }
527 
528  dist.setY(scaleFromYvsX(dir, dist.x()));
529  LOG_DEBUG<<"dist x = "<<dist.x()<<"\ty = "<<dist.y()<<"\tz = "<<dist.z()<<endm;
530  return dist.mag(); // 3D distance
531 }
532 
533 
534 Double_t StIstSlowSimMaker::direction(const Double_t x) const
535 {
536  return (x>0) - (x<0);
537 }
538 
539 Double_t StIstSlowSimMaker::scaleFromYvsX(const StThreeVectorD vec, const Double_t a) const
540 {
541  if(fabs(vec.x())>1e-6) return a*vec.y()/vec.x();
542  else{
543  LOG_WARN<<"StIstSlowSimMaker::scaleFromYvsX x component is 0"<<endm;
544  return 0.;
545  }
546 }
547 
548 Double_t StIstSlowSimMaker::scaleFromYvsZ(const StThreeVectorD vec, const Double_t a) const
549 {
550  if(fabs(vec.z())>1e-6) return a*vec.y()/vec.z(); // scale a by y/z
551  else{
552  LOG_WARN<<"StIstSlowSimMaker::scaleFromYvsZ z component is 0"<<endm;
553  return 0.;
554  }
555 }
StIstRawHit::getLadder
unsigned char getLadder() const
1-24
Definition: StIstRawHit.cxx:46
StIstRawHit::getColumn
unsigned char getColumn() const
1-12
Definition: StIstRawHit.cxx:62
StMaker::Clear
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
StIstConsts::kIstNumLadders
const int kIstNumLadders
24 IST Ladders
Definition: StEnumerations.h:711
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StIstRawHit::getRow
unsigned char getRow() const
1-64
Definition: StIstRawHit.cxx:56
StIstSlowSimMaker::getMCHitRowAndColumn
void getMCHitRowAndColumn(const StMcIstHit *istMChit, UShort_t &meanColumn, UShort_t &meanRow) const
Definition: StIstSlowSimMaker.cxx:246
StIstSlowSimMaker::findPad
void findPad(const StThreeVectorD hitPos, UShort_t &column, UShort_t &row, Double_t &rPhiPos_mean, Double_t &zPos_mean) const
Definition: StIstSlowSimMaker.cxx:462
StIstRawHit::getChannelId
int getChannelId() const
0-110591
Definition: StIstRawHit.cxx:40
StIstConsts::kIstNumColumnsPerSensor
const int kIstNumColumnsPerSensor
12 columns in beam direction per each sensor
Definition: StEnumerations.h:716
kStWarn
Definition: Stypes.h:42
kStOK
Definition: Stypes.h:40
StIstSlowSimMaker::checkPadCrossing
void checkPadCrossing(const StThreeVectorD inPos, const StThreeVectorD outPos, StThreeVectorD mcLocalDir, Double_t dS, vector< StThreeVectorD > &cross_vec) const
Definition: StIstSlowSimMaker.cxx:385
TObjectSet::GetObject
virtual TObject * GetObject() const
The depricated method (left here for the sake of the backward compatibility)
Definition: TObjectSet.h:56
kStErr
Definition: Stypes.h:44
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
StIstRawHit::getSensor
unsigned char getSensor() const
1-6
Definition: StIstRawHit.cxx:51
StIstConsts::kIstNumRowsPerSensor
const int kIstNumRowsPerSensor
64 rows in r-phi direction per each sensor
Definition: StEnumerations.h:715
kStOk
Definition: Stypes.h:41
StIstSlowSimMaker::Clear
void Clear(Option_t *opts="")
User defined functions.
Definition: StIstSlowSimMaker.cxx:231
StIstSlowSimMaker::InitRun
Int_t InitRun(Int_t runNumber)
Definition: StIstSlowSimMaker.cxx:86
Generated by   doxygen 1.8.5