StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMtdSimMaker.cxx
1 /***************************************************************************
2  *
3  * $Id: StMtdSimMaker.cxx,v 1.13 2015/11/02 18:19:46 marr Exp $
4  *
5  * Author: Frank Geurts
6  *
7  * Modified: Yi Yang (yiyang@bnl.gov)
8  * Date: 2014-05-06
9  *
10  ***************************************************************************
11  *
12  * Description: StMtdSimMaker class for Muon Telescope Simulations
13  *
14  ***************************************************************************
15  *
16  *
17  **************************************************************************/
18 #include <Stiostream.h>
19 #include "StMtdSimMaker.h"
20 
21 // SCL
22 #include <math.h>
23 #include "TRandom.h"
24 #include "SystemOfUnits.h"
25 #include "phys_constants.h"
26 #include "StThreeVectorD.hh"
27 #include "Random.h"
28 #include "RanluxEngine.h"
29 #include "RandGauss.h"
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "TNtuple.h"
33 #include "TFile.h"
34 
35 // g2t tables and collections
36 #include "tables/St_g2t_mtd_hit_Table.h"
37 #include "tables/St_g2t_track_Table.h"
38 #include "tables/St_mtdGeant2BacklegIDMap_Table.h"
39 #include "StMcTrack.hh"
40 
41 #include "StEventTypes.h"
42 
43 static RanluxEngine engine;
44 static RandGauss ranGauss(engine);
45 
46 Int_t geant2backlegIDMap[30];
47 
48 ClassImp(StMtdSimMaker)
49 
50 const float StMtdSimMaker::kMtdPadWidth = 3.8 + 0.6;
51 
52 //_____________________________________________________________________________
53 StMtdSimMaker::StMtdSimMaker(const char *name):StMaker(name),
54  mGeantData(0),
55  mEvent(0),
56  mMtdCollection(0),
57  mNMtdHits(0),
58  mMtdHitsFromGeant(0),
59 
60  mBetaHist(0),
61  mPathLHist(0),
62  mTofHist(0),
63  mRecMass(0),
64 
65  mCellGeant(0),
66  mVpdGeant(0),
67  mNCellGeant(0),
68  mNVpdGeant(0),
69  mDeGeant(0),
70  mTofGeant(0),
71 
72  mCellSeen(0),
73  mVpdSeen(0),
74  mNCellSeen(0),
75  mNVpdSeen(0),
76  mDeSeen(0),
77  mT0Seen(0),
78  mTofSeen(0),
79  mTofResSeen(0),
80  mVpdResSeen(0),
81 
82  mCellReco(0),
83  mVpdReco(0),
84  mNCellReco(0),
85  mNVpdReco(0),
86  mTDCReco(0),
87  mADCReco(0),
88  mT0Reco(0),
89  mTofResReco(0),
90  mVpdResReco(0),
91  mTACorr(0),
92  mModHist(0),
93  QABacklegChannel(0),
94 
96  mdE(0),
97  mdS(0),
98  mNumberOfPhotoelectrons(0),
99  mT(0),
100  mTime(0),
101  mTime1(0),
102  mPMlength(0),
103  mAdc(0),
104  mTdc(0),
105 
106  starHall(0) {
107 
108  memset(mModuleChannel,0,5*24*sizeof(Int_t));
109 
110  mBookHisto=kTRUE;
111  mWriteHisto=kFALSE;
112  mHistFile="mtdsim.root";
113  mWriteStEvent=kTRUE;
114  Reset();
115 }
116 
117 //_____________________________________________________________________________
118 StMtdSimMaker::~StMtdSimMaker() {
119  // do nothing
120 }
121 
122 
123 //_____________________________________________________________________________
124 Int_t StMtdSimMaker::Init() {
125  if(mBookHisto) bookHistograms();
126  return StMaker::Init();
127 }
128 
129 //_____________________________________________________________________________
130 void StMtdSimMaker::Reset() {
131  mGeantData = 0;
132  mEvent = 0;
133  if (mWriteStEvent) delete mMtdCollection;
134 }
135 //_____________________________________________________________________________
136 Int_t StMtdSimMaker::InitRun(Int_t runnumber) {
137 
138  // Load geant2backlegID Map
139  // Extract MTD maps from database
140  LOG_INFO << "Retrieving geant2backlegID table from database ..." << endm;
141  TDataSet *dataset = GetDataBase("Geometry/mtd/mtdGeant2BacklegIDMap");
142  St_mtdGeant2BacklegIDMap *mtdGeant2BacklegIDMap = static_cast<St_mtdGeant2BacklegIDMap*>(dataset->Find("mtdGeant2BacklegIDMap"));
143  if ( !mtdGeant2BacklegIDMap ){
144  LOG_ERROR << "No mtdTrayToTdigMap table found in database" << endm;
145  return kStErr;
146  }
147  mtdGeant2BacklegIDMap_st *mGeant2BLTable = static_cast<mtdGeant2BacklegIDMap_st*>(mtdGeant2BacklegIDMap->GetTable());
148  for ( Int_t i = 0; i < 30; i++ ){
149  geant2backlegIDMap[i] = 0;
150  geant2backlegIDMap[i] = (Int_t)mGeant2BLTable->geant2backlegID[i];
151  }
152 
153 
155  int channel(0);
156  for (int module=0; module<5; module++){
157  for (int cell=0; cell<24; cell ++){
158  mModuleChannel[module][cell] = ++channel;
159  }
160  }
161 
162  return kStOK;
163 }
164 
165 //_____________________________________________________________________________
166 Int_t StMtdSimMaker::FinishRun(Int_t runnumber) {
167  // do nothing
168  return kStOk;
169 }
170 
171 //_____________________________________________________________________________
173  if(mWriteHisto){
174  LOG_INFO << "StMtdSimMaker::Finish writing mtdsim.root ...";
175  TFile theFile(mHistFile.c_str(),"RECREATE","mtdsim");
176  theFile.cd();
177  writeHistograms();
178  theFile.Close();
179  }
180 
181  return kStOK;
182 }
183 
184 
185 //_____________________________________________________________________________
187 
189  mGeantData = GetInputDS("geant"); // in bfc chain
190  if(!mGeantData) { // when reading the geant.root file
191  mGeantData = GetInputDS("geantBranch");
192  }
193  if(!mGeantData) {
194  LOG_WARN << " No GEANT data loaded. Exit! " << endm;
195  return kStWarn;
196  }
197  LOG_INFO << " Found GEANT data -- loading MTD hits... " << endm;
198 
200  St_g2t_mtd_hit* g2t_mtd_hits = dynamic_cast<St_g2t_mtd_hit*>(mGeantData->Find("g2t_mtd_hit"));
201  if(!g2t_mtd_hits){
202  LOG_WARN << " No MTD hits in GEANT" << endm;
203  return kStWarn;}
204 
205  mNMtdHits = g2t_mtd_hits->GetNRows();
206  if(!mNMtdHits){
207  LOG_WARN << " No MTD hits in GEANT" << endm;
208  return kStWarn;}
209 
210  LOG_INFO << " Found MTD hits: " << mNMtdHits << endm;
211 
212  mEvent = (StEvent*)GetInputDS("StEvent");
213  mMtdHitsFromGeant = g2t_mtd_hits->GetTable();
214 
216 
217  return kStOK;
218 }
219 
221 int StMtdSimMaker::CalcCellId(Int_t volume_id, Float_t ylocal, int &ibackleg,int &imodule,int &icell)
222 {
224  Int_t backlegTemp;
225  Int_t ires = volume_id/100;
226  backlegTemp = ires/10;
227  ibackleg = geant2backlegIDMap[backlegTemp - 1];
228 
229  imodule = ires%10;
230  if ( ibackleg >= 12 && ibackleg <= 20 ) imodule = imodule + 1;
231 
233  icell = Int_t((ylocal + kMtdPadWidth * kNCell/2) / kMtdPadWidth);
234  // Get the correct cell ID
235  if ( imodule > 3 ) icell = 11 - icell;
236 
238  if(ibackleg<0 || ibackleg>kNBackleg) return -3;
239  if(imodule <0 || imodule >kNModule ) return -2;
240  if(icell <0 || icell >= kNCell ) return -1;
241 
242  return icell + 100*(imodule+100*ibackleg);
243 }
244 
245 
248 {
249  std::map<int,StMtdHit*> myMap;
250  std::map<int,StMtdHit*>::const_iterator myIter;
251 
252  for (int ihit=0;ihit<mNMtdHits;ihit++) {
253  g2t_mtd_hit_st &ghit = mMtdHitsFromGeant[ihit];
254  if(ghit.s_track<=0.0 || ghit.de <=0.0) continue;
255 
256  if(mBookHisto) {
257  mDeGeant->Fill(ghit.de / keV);
258  mTofGeant->Fill(ghit.tof / nanosecond);
259  }
260 
261 
262 
263  Int_t icell, imodule, ibackleg;
264  int cellId = CalcCellId(ghit.volume_id, ghit.x[1],ibackleg,imodule,icell);
265 
266  LOG_DEBUG << "Volume ID: "<< ghit.volume_id << " Cell ID: " << cellId << endm;
267  LOG_DEBUG << " icell: " << icell << " imodule: " << imodule << " ibackleg: " << ibackleg << " ghit.tof: " << ghit.tof << endm;
268  LOG_DEBUG << " track: " << ghit.track_p << endm;
269 
270  if (cellId<0) continue;
271  StMtdHit *&sthit = myMap[cellId];
272  if (!sthit) { sthit = new StMtdHit;}
273  else { if (sthit->tof()>ghit.tof) continue; }
274  Float_t wt = 1.0;
275 
276  //fg temporarily fix MTD resolution at 99ps
277  Double_t tof= ghit.tof/nanosecond + ranGauss.shoot()*(99.e-12)/nanosecond;
278  Double_t de = ghit.de * wt;
279 
280  // Calculate MTD leading and tailing edge time
281  // Assume the velocity = 56 ps/cm, module length = 87 cm
282  double vel = 56.e-3; // ns/cm
283  double leadingW = ghit.tof/nanosecond - ghit.x[2]*vel;
284  double leadingE = ghit.tof/nanosecond + ghit.x[2]*vel;
285  double trailingW = leadingW + 15;
286  double trailingE = leadingE + 15;
287 
288  int channel1 = mModuleChannel[imodule - 1][icell];
289  int channel2 = channel1 + 12;
290 
291  LOG_INFO << "First hit for this cell, assigning channel 1: " << channel1 << " channel 2: " << channel2 << endm;
292 
293  //Fill Histograms here AJ////
294  QABacklegChannel->Fill(ibackleg, channel1);
295  QABacklegChannel->Fill(ibackleg, channel2);
296 
297  //mCellGeant->Fill(icell);
298  mDeGeant->Fill(ghit.de);
299  mTofGeant->Fill(ghit.tof);
300 
301  //mCellSeen->Fill(cellId);
302  mDeSeen->Fill(de);
303  mTofSeen->Fill(tof);
305 
306  sthit->setBackleg(ibackleg);
307  sthit->setModule(imodule);
308  sthit->setCell(icell);
309  sthit->setLeadingEdgeTime(pair<double,double>(leadingW, leadingE));
310  sthit->setTrailingEdgeTime(pair<double,double>(trailingW, trailingE));
311  sthit->setAssociatedTrack(NULL); //done in StMtdMatchMaker
312  sthit->setIdTruth(ghit.track_p, 100);
313 
314  LOG_INFO << "sthit " << sthit->tof() << " tof:" << tof << endm;
315  }
316 
317  int nRealHits = 0;
318  if (mEvent)
319  {
320  mMtdCollection = mEvent->mtdCollection();
321  if(!mMtdCollection)
322  {
323  mMtdCollection= new StMtdCollection();
324  mEvent->setMtdCollection(mMtdCollection);
325  }
326  nRealHits = mMtdCollection->mtdHits().size();
327  }
328  else mMtdCollection= new StMtdCollection();
329 
330  int nMcHits = 0;
331  for (myIter=myMap.begin(); myIter!=myMap.end(); ++myIter) {
332  mMtdCollection->addHit((*myIter).second);
333  nMcHits++;
334  }
335  if (mEvent) {
336  LOG_INFO << "... " << nMcHits << " MC hits stored in StEvent with " << nRealHits << " real hits! " << endm;
337  }
338 
339  return kStOK;
340 }
341 
342 
344 {
345  //only done if Histogram setting is turned on
346  mBetaHist=new TH1F("mBetaHist","mBetaHist", 100, -2, 2);
347  mPathLHist=new TH1F("mPathLHist","mPathLHist", 100, -2, 500);//cm's
348  mTofHist=new TH1F("mTofHist","mTofHist", 1000, -10, 10000);
349  mRecMass=new TH1F("mRecMass","mRecMass", 1000, -2, 4);
350 
351  mCellGeant = new TH2F("CellGeant","CellGeant",192,0.,192.,120,1.,120.);
352  mNCellGeant = new TH2F("NCellGeant","NCellGeant",192,0.,192.,120,1.,120.);
353  mDeGeant = new TH1F("DeGeant","DeGeant",1000,0.,10.);
354  mTofGeant = new TH1F("TofGeant","TofGeant",1000,0.,20.);
355 
356  mCellSeen = new TH2F("CellSeen","CellSeen",192,0.,192.,120,1.,120.);
357  mNCellSeen = new TH2F("NCellSeen","NCellSeen",192,0.,192.,120,1.,120.);
358  mDeSeen = new TH1F("DeSeen","DeSeen",1000,0.,10.);
359  mT0Seen = new TH1F("T0Seen","T0Seen",1000,0.,20.);
360  mTofSeen = new TH1F("TofSeen","TofSeen",1000,0.,20.);
361 
362  mTofResSeen = new TH1F("TofResSeen","TofResSeen",1001,-500.,500.);
363 
364  mCellReco = new TH2F("CellReco","CellReco",192,0.,192.,120,1.,120.);
365  mNCellReco = new TH2F("NCellReco","NCellReco",192,0.,192.,120,1.,120.);
366  mADCReco = new TH1F("ADCReco","ADCReco",4096,0.,4096.);
367  mTDCReco = new TH1F("TDCReco","TDCReco",4096,0.,4096.);
368  mT0Reco = new TH1F("T0Reco","T0Reco",1000,0.,20.);
369  mTofResReco = new TH1F("TofResReco","TofResReco",1000,-300.,300.);//ps
370  mTACorr = new TH2F("TACorr","TACorr",512,0.,4096.,512,0.,4096.);
371 
372  mModHist = new TH1F("ModuleHist","ModuleHist",201,-100,100);
373  QABacklegChannel = new TH2I("QABacklegChannel", "QABacklegChannel", 30, 1., 30., 120, 1., 120.);
374  return kStOk;
375 
376 }
377 
378 //_____________________________________________________________________________
379 Int_t StMtdSimMaker::writeHistograms()
380 {
381  //only done if Histogram setting is turned on
382 
383  mBetaHist->Write();
384  mPathLHist->Write();
385  mTofHist->Write();
386  mRecMass->Write();
387 
388  mCellGeant->Write();
389  mNCellGeant->Write();
390  mDeGeant->Write();
391  mTofGeant->Write();
392 
393  mCellSeen->Write();
394  mNCellSeen->Write();
395  mDeSeen->Write();
396  mT0Seen->Write();
397  mTofSeen->Write();
398  mTofResSeen->Write();
399 
400  mCellReco->Write();
401  mNCellReco->Write();
402  mADCReco->Write();
403  mTDCReco->Write();
404  mT0Reco->Write();
405  mTofResReco->Write();
406  mTACorr->Write();
407 
408  mModHist->Write();
409  QABacklegChannel->Write();
410  return kStOk;
411 }
TH2I * QABacklegChannel
T-A Slewing Correlation.
TH2F * mNCellSeen
Vpd tubeId after DetectorResponse.
TH1F * mRecMass
total time of flight of partilce
string mHistFile
switch to enable Maker to write out simulated hits to StEvent
Definition: StMtdSimMaker.h:97
TH1F * mPathLHist
speed of particles hitting tof
Definition: StMtdSimMaker.h:99
Definition: tof.h:15
TH1F * mDeGeant
of vpd tubes of geant hit
TH1F * mDeSeen
of vpd tubes after DetectorResponse
TH1F * mT0Seen
deposited-energy after DetectorResponse
TH1F * mTofGeant
deposited-energy in geant hit
TH2F * mTACorr
vpd time resolution after recon
TH1F * mTDCReco
of vpd tubes after recon
Int_t CalcCellId(Int_t volume_id, Float_t ylocal, Int_t &ibackleg, Int_t &imodule, Int_t &icell)
This will calculate the cell ID as well as decode the module # and backleg # to store in an MTD Colle...
TH2F * mNCellReco
Vpd tubeId after recon.
TH2F * mCellSeen
tof in geant hit
TH2F * mCellReco
vpd time resolution after DetectorResponse
static const float kMtdPadWidth
Pad Width: 38mm padwidth + 6mm innerspacing.
Definition: StMtdSimMaker.h:90
StEvent * mEvent
geant table
Definition: StMtdSimMaker.h:75
virtual Int_t Finish()
TH1F * mTofResSeen
smeared-tof after DetectorResponse
Int_t bookHistograms()
TH1F * mT0Reco
ADC recon – empty.
Definition: Stypes.h:42
Definition: Stypes.h:40
virtual Int_t Make()
TH1F * mADCReco
TDC recon.
Int_t InitRun(Int_t)
Int_t FastCellResponse()
This will define the maps to store the StMtdHits and will call the necessary functions to obtain the ...
TH1F * mTofHist
speed of particles hitting tof
TH2F * mNCellGeant
Vpd tubeId of geant hit.
TH2F * mCellGeant
reconstructed mass of particle
Definition: Stypes.h:44
Definition: Stypes.h:41
TH1F * mModHist
T-A Slewing Correlation.
virtual TDataSet * Find(const char *path) const
Definition: TDataSet.cxx:362