StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructFlat.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructFlat.cxx,v 1.7 2012/11/16 21:23:18 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: EStructEventReader which generates events Flat in eta,phi.
10  *
11  **********************************************************************/
12 #include "StEStructFlat.h"
13 
14 #include "StEStructPool/AnalysisMaker/StEStructEventCuts.h"
15 #include "StEStructPool/AnalysisMaker/StEStructTrackCuts.h"
16 #include "StEStructPool/EventMaker/StEStructEvent.h"
17 #include "StEStructPool/EventMaker/StEStructTrack.h"
18 #include "StEStructPool/EventMaker/StEStructCentrality.h"
19 
20 StEStructFlat::StEStructFlat() {
21  mFlatEvent = 0;
22  mECuts = 0;
23  mTCuts = 0;
24  mInChain = false;
25  mEventsToDo = 0;
26  mCentBin = 0;
27  mEventCount = 0;
28  mAmDone = false;
29  mgRand2Good = false;
30 }
31 StEStructFlat::StEStructFlat( StEStructEventCuts* ecuts,
32  StEStructTrackCuts* tcuts,
33  bool inChain,
34  int centBin,
35  int eventsToDo) {
36  mFlatEvent = 0;
37  mECuts = ecuts;
38  mTCuts = tcuts;
39  mInChain = inChain;
40  mEventsToDo = eventsToDo;
41  mCentBin = centBin;
42  mEventCount = 0;
43  mAmDone = false;
44  mgRand2Good = false;
45 };
46 
47 bool StEStructFlat::hasGenerator() { return true; };
48 
49 //-------------------------------------------------------------------------
50 void StEStructFlat::setSeed(int iseed) {
51  cout << " Calling srand48(" << iseed << ")" << endl;
52  srand48(iseed);
53 }
54 
55 //-------------------------------------------------------------------------
56 StEStructEvent* StEStructFlat::next() {
57  if(mEventCount==mEventsToDo){
58  mAmDone = true;
59  return (StEStructEvent*) NULL;
60  }
61 
62  if (!mInChain ) {
63  mFlatEvent = new StEStructEvent();
64  generateEvent();
65  mEventCount++;
66  }
67 
68  float z = mFlatEvent->Vz();
69  int nTracks = countGoodTracks();
70  mFlatEvent->SetCentrality( (double) nTracks );
71  StEStructCentrality* cent=StEStructCentrality::Instance();
72  int jCent = cent->centrality( mFlatEvent->Centrality() );
73  if (((mCentBin >= 0) && (jCent != mCentBin)) ||
74  !mECuts->goodPrimaryVertexZ(z) ||
75  !mECuts->goodCentrality(mFlatEvent->Centrality())) {
76  mECuts->fillHistogram(mECuts->centralityName(),(float)nTracks,false);
77  mECuts->fillHistogram(mECuts->primaryVertexZName(),z,false);
78  return (StEStructEvent*) NULL;
79  } else {
80  mFlatEvent->FillChargeCollections();
81  mECuts->fillHistogram(mECuts->centralityName(),(float)nTracks,true);
82  mECuts->fillHistogram(mECuts->primaryVertexZName(),z,true);
83  return mFlatEvent;
84  }
85 }
86 
87 //--------------------------------------------------------------------------
88 void StEStructFlat::generateEvent() {
89  mFlatEvent->SetVx(0);
90  mFlatEvent->SetVy(0);
91  mFlatEvent->SetVz(30*gRand48());
92  mFlatEvent->SetBField(0.5);
93 
94  fillTracks(mFlatEvent);
95 }
96 
97 //--------------------------------------------------------------------------
98 void StEStructFlat::fillTracks(StEStructEvent* estructEvent) {
99 
100  mnumTracks = 0;
101  StEStructTrack* eTrack = new StEStructTrack();
102  int pid, sign;
103 
104 // int numCharge = int( -5*log(drand48()) );
105  int numCharge, numInEta = 850;
106  double v2 = 0.05, pi = 3.1415926;
107  double eta, quadEta = 0.2, etaMax = 2;
108  double phi, phiOff, sectorWidth;
109  int numSector = 12;
110 
111  sectorWidth = 360 / numSector;
112  phiOff = 360*drand48();
113 
114  float p[3], pt, pz, v[3];
115 
116  double etaMaxAmp = 1 + quadEta*etaMax*etaMax;
117  numCharge = int( numInEta*etaMax * (3+quadEta*etaMax*etaMax) / (3+quadEta) );
118  for(int i=0;i<numCharge;i++) {
119 
120  eTrack->SetInComplete();
121  if (i < numCharge/2) {
122  sign = +1;
123  pid = 211;
124  } else {
125  sign = -1;
126  pid = -211;
127  }
128  double etaAmp = 0, randAmp = 1;
129  while (randAmp > etaAmp) {
130  eta = etaMax*(2*drand48() - 1);
131  etaAmp = 1 + quadEta*eta*eta;
132  randAmp = drand48()*etaMaxAmp;
133  }
134 
135  // Put flow into phi.
136  phi = 360*drand48();
137  double h = (1+v2)*drand48();
138  while ( h > (1+v2*cos(2*phi*pi/180)) ) {
139  phi = 360*drand48();
140  h = (1+v2)*drand48();
141  }
142  phi += phiOff;
143  while (phi > 180) {
144  phi -= 360;
145  }
146 
147  pt = 0.1 - 0.5 * log(drand48());
148  pz = sqrt( pt*pt + 0.139*0.139) * (exp(eta)-exp(-eta)) / 2;
149 
150  p[0] = pt*cos(phi*pi/180);
151  p[1] = pt*sin(phi*pi/180);
152  p[2] = pz;
153  v[0] = 0;
154  v[1] = 0;
155  v[2] = estructEvent->Vz();
156 
157  // Check outer radius of track or radius of track at pad-plane,
158  // which ever is bigger. Reject track if this is less than 137.
159  if (maxRadius(eta,pt,v[2]) < 155.0) {
160  continue;
161  }
162 
163  if (!isTrackGood(v,p,eta)) {
164  mTCuts->fillHistograms(false);
165  continue;
166  }
167  mTCuts->fillHistograms(true);
168  if (pt<0.15) continue;
169 
170  mnumTracks++;
171 
172  float *gdca = globalDCA(p,v);
173  eTrack->SetBx(gdca[0]);
174  eTrack->SetBy(gdca[1]);
175  eTrack->SetBz(gdca[2]);
176  eTrack->SetBxPrimary(gdca[0]);
177  eTrack->SetByPrimary(gdca[1]);
178  eTrack->SetBzPrimary(gdca[2]);
179  eTrack->SetBxGlobal(gdca[0]);
180  eTrack->SetByGlobal(gdca[1]);
181  eTrack->SetBzGlobal(gdca[2]);
182  delete [] gdca;
183 
184  eTrack->SetPx(p[0]);
185  eTrack->SetPy(p[1]);
186  eTrack->SetPz(p[2]);
187  eTrack->SetEta(eta);
188  eTrack->SetPhi(phi*pi/180);
189  eTrack->SetDedx(0);
190  eTrack->SetChi2(1);
191  eTrack->SetTopologyMapData(0, 0xffffff80);
192  eTrack->SetTopologyMapData(1, 0x3fff);
193  eTrack->SetTopologyMapTPCNHits(45);
194  eTrack->SetNMaxPoints(45);
195  eTrack->SetNFoundPoints(45);
196  eTrack->SetNFitPoints(45);
197 
198  if (pid<0) {
199  eTrack->SetCharge(-1);
200  } else {
201  eTrack->SetCharge(1);
202  }
203  estructEvent->AddTrack(eTrack);
204  }
205 
206  estructEvent->SetCentrality( (double) mnumTracks );
207  delete eTrack;
208  return;
209 }
210 
211 
212 //-------------------------------------------------------------
213 // This method calculates maximum distance of track from beam axis.
214 // This point might be at the endplane or at twice the radius of curvature.
215 double StEStructFlat::maxRadius(double eta, double pt, double vz) {
216  double pi = 3.14159265359;
217  double lambda = pi/2 - 2*atan(exp(-eta));
218  double s;
219  if (lambda > 0) {
220  s = (+200 - vz) / sin(lambda);
221  } else {
222  s = (-200 - vz) / sin(lambda);
223  }
224  double r = pt / 0.0015;
225  double phi = s * cos(lambda) / r;
226  if (phi < pi) {
227  return r * sqrt(2 - 2*cos(phi));
228  }
229  return 2*r;
230 }
231 
232 //-------------------------------------------------------------
233 // This method checks all track cuts.
234 // No histogramming or copying data around.
235 bool StEStructFlat::isTrackGood(float *v, float *p, float eta) {
236  float pt = sqrt(p[0]*p[0]+p[1]*p[1]);
237  if (pt < 0.15) {
238  return false;
239  }
240 
241  float phi = atan2((double)p[1], (double)p[0]);
242  float *gdca = globalDCA(p,v);
243  float _r = pt/0.139;
244  float yt = log(sqrt(1+_r*_r)+_r);
245 
246  bool useTrack = true;
247  useTrack = (mTCuts->goodGlobalDCA(gdca[3]) && useTrack);
248  useTrack = (mTCuts->goodEta(eta) && useTrack);
249  useTrack = (mTCuts->goodPhi(phi) && useTrack);
250  useTrack = (mTCuts->goodPt(pt) && useTrack);
251  useTrack = (mTCuts->goodYt(yt) && useTrack);
252  delete [] gdca;
253 
254  return useTrack;
255 }
256 int StEStructFlat::countGoodTracks() {
257  return mnumTracks;
258 }
259 //--------------------------------------------------------------------------
260 double StEStructFlat::gRand48() {
261  double x1, x2, w;
262 
263  if (mgRand2Good) {
264  mgRand2Good = false;
265  return mgRand2;
266  }
267  do {
268  x1 = 2.0 * drand48() - 1.0;
269  x2 = 2.0 * drand48() - 1.0;
270  w = x1 * x1 + x2 * x2;
271  } while ( w >= 1.0 );
272 
273  w = sqrt( (-2.0 * log( w ) ) / w );
274  mgRand2 = x2 * w;
275  mgRand2Good = true;
276  return x1*w;
277 }
278 
279 
280 
281 /**********************************************************************
282  *
283  * $Log: StEStructFlat.cxx,v $
284  * Revision 1.7 2012/11/16 21:23:18 prindle
285  * EventCuts and TrackCuts were moved to EventReader. Remove that code from
286  * these readers.
287  *
288  * Revision 1.6 2006/04/06 01:03:30 prindle
289  *
290  * Rationalization of centrality binning, as described in AnalysisMaker checkin.
291  *
292  * Revision 1.5 2006/02/22 22:05:35 prindle
293  * Removed all references to multRef (?)
294  *
295  * Revision 1.4 2005/09/23 23:37:18 prindle
296  *
297  * Starting to add vertex distribution and track acceptance dependance on
298  * number of possible hits.
299  * Make Pythia interface look like Hijing interface so it now works within
300  * my Fluctuation and Correlation framework.
301  *
302  * Revision 1.3 2005/09/07 20:22:47 prindle
303  *
304  *
305  * Flat: Random changes to eta and phi distributions (which don't have to be flat).
306  *
307  * Revision 1.2 2003/11/25 22:45:14 prindle
308  * Commiting changes so I can move code to rhic
309  *
310  * Revision 1.1 2003/11/21 23:48:00 prindle
311  * Include my toy event generator in cvs
312  *
313  *
314  *********************************************************************/
Generated by   doxygen 1.8.5