StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructRQMD.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructRQMD.cxx,v 1.4 2012/11/16 21:23:19 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: EStructEventReader which generates events Flat in eta,phi.
10  *
11  **********************************************************************/
12 #include "StEStructRQMD.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 
19 StEStructRQMD::StEStructRQMD(): meventCount(0), meventsToDo(0), mAmDone(false) {
20  mNFile = 0;
21  mIFile = 0;
22  mMaxFiles = 200;
23  mFileDir = "/aztera/estruct/aya/simulations/RQMD/AuAu200/with_rescattering/out/";
24  mFile = NULL;
25 };
26 
27 StEStructRQMD::StEStructRQMD(int nevents, StEStructEventCuts* ecuts, StEStructTrackCuts* tcuts): meventCount(0), mAmDone(false){
28 
29  mNFile = 0;
30  mIFile = 0;
31  mMaxFiles = 200;
32  mFileDir = "/aztera/estruct/aya/simulations/RQMD/AuAu200/with_rescattering/out/";
33  mFile = NULL;
34 
35  meventsToDo=nevents;
36  mECuts=ecuts;
37  mTCuts=tcuts;
38 };
39 
40 bool StEStructRQMD::hasGenerator() { return true; };
41 
42 
43 //-------------------------------------------------------------------------
44 StEStructEvent* StEStructRQMD::next() {
45 
46  if(meventCount==meventsToDo){
47  mAmDone=true;
48  return (StEStructEvent*)NULL;
49  }
50  return generateEvent();
51 }
52 
53 //--------------------------------------------------------------------------
54 StEStructEvent* StEStructRQMD::generateEvent() {
55 
56  StEStructEvent* retVal=NULL;
57 
58  retVal = new StEStructEvent();
59 
60  fillTracks(retVal);
61  retVal->SetCentrality( (float) mnumTracks );
62  if (!mECuts->goodCentrality(retVal->Centrality())) {
63  delete retVal;
64  retVal=NULL;
65  } else {
66  retVal->FillChargeCollections();
67  }
68 
69  return retVal;
70 }
71 
72 //--------------------------------------------------------------------------
73 void StEStructRQMD::fillTracks(StEStructEvent* estructEvent) {
74  char buf[1024];
75  int skip[] = { 2, 244, 368, 440, 520, 698, 826, 934, 1060,
76  1061, 1108, 1143, 1292, 1458, 1495, 1585, 1710, 1766,
77  1857, 2015, 2065, 2101, 2147, 2158, 2344, 2448};
78  int nSkip = 26;
79 
80  const int kKMinus = 12;
81  const int kKPlus = 11;
82  const int kLambda = 18;
83  const int kLambdaBar = 26;
84  const int kProton = 14;
85  const int kPbar = 15;
86  const int kPhi = 51;
87  const int kPiMinus = 9;
88  const int kPiPlus = 8;
89  const int kPi0 = 7;
90  const int kXi0 = 22;
91  const int kXiMinus = 23;
92  const int kAntiXi0 = 30;
93  const int kAntiXiPlus = 31;
94  const int kOmega = 24;
95  const int kOmegaBar = 32;
96  const int kSigmaPlus = 19;
97  const int kSigma0 = 20;
98  const int kSigmaMinus = 21;
99  const int kAntiSigmaPlus = 27;
100  const int kAntiSigma0 = 28;
101  const int kAntiSigmaMinus = 29;
102 
103 
104  mnumTracks=0;
105  StEStructTrack* eTrack = new StEStructTrack();
106 
107  char tDum[50];
108 
109  // On first call file will not have been open yet.
110  // Also after reading event we check if we are at end of file. In that
111  // case we close the file.
112  if (!mFile) {
113  if (mIFile >= mMaxFiles) {
114  mAmDone = true;
115  return;
116  }
117  mNFile++;
118  mIFile++;
119  for (int i=0;i<nSkip;i++) {
120  if (skip[i] == mNFile) {
121  mNFile++;
122  mIFile++;
123  }
124  }
125  sprintf(buf, "%s%s%d", mFileDir, "file9_", mNFile);
126  cout << "opening file " << buf << endl;
127  mFile = new ifstream(buf);
128  mLine = 0;
129  if (mFile->is_open() == 0) {
130  cout << "file not found, trying next one. " << buf << endl;
131  delete mFile;
132  mFile = NULL;;
133  return;
134  }
135  // Read and discard first line of file.
136  *mFile >> tDum; mLine++;
137  if (mFile->fail()) {
138  cout << "Failed to read first line of file????, going on to next file" << endl;
139  delete mFile;
140  mFile = NULL;;
141  return;
142  }
143  }
144 
145  int mNFreezeOutPart = 0;
146 
147  // read header and extract number of tracks in event
148  *mFile >> tDum >> tDum >> tDum; mLine++;
149  int tNucleonTarget = atoi(tDum);
150  *mFile >> tDum >> tDum >> tDum; mLine++;
151  *mFile >> tDum >> tDum >> tDum; mLine++;
152  int tNucleonProj = atoi(tDum);
153  *mFile >> tDum >> tDum >> tDum; mLine++;
154  *mFile >> tDum >> tDum >> tDum; mLine++;
155  *mFile >> tDum >> tDum >> tDum; mLine++;
156  *mFile >> tDum >> tDum >> tDum; mLine++;
157  *mFile >> tDum >> tDum >> tDum; mLine++;
158  *mFile >> mNFreezeOutPart >> tDum >> tDum >> tDum; mLine++;
159  mNFreezeOutPart += (tNucleonTarget+tNucleonProj);
160  *mFile >> tDum >> tDum >> tDum >> tDum >> tDum; mLine++;
161 
162  if ((mNFreezeOutPart < 394) || (15000 < mNFreezeOutPart)) {
163  cout << "An unusual number of tracks in this event?? ";
164  cout << mNFreezeOutPart << " tracks at line " << mLine << endl;
165  }
166 
167  // read event
168  int tPid1,tPid2,tCharge,tNClnt,tlastcl;
169  float tt,tx,ty,tz,tE,tPx,tPy,tPz,tM,tDecay;
170 
171  while (mNFreezeOutPart) {
172  eTrack->SetInComplete();
173  int gPID = 0;
174  *mFile >> tPid1 >> tPid2 >> tt >> tx >> ty >> tz
175  >> tE >> tPx >> tPy >> tPz >> tM
176  >> tDecay >> tNClnt >> tlastcl ; mLine++;
177 
178  float pt = sqrt(pow(tPx,2)+pow(tPy,2));
179  float eta = getPseudoRapidity(pt, tPz);
180  float phi = atan2(tPy, tPx);
181  mNFreezeOutPart--;
182 
183  tCharge = 0;
184  if (tPid1 == 14 && tPid2 == -2) { // K-
185  gPID = kKMinus;
186  tCharge = -1;
187  }
188  if (tPid1 == 14 && tPid2 == 2) { // K+
189  gPID = kKPlus;
190  tCharge = +1;
191  }
192  if (tPid1 == 13 && tPid2 == 0) { // Lambda
193  gPID = kLambda;
194  tCharge = 0;
195  }
196  if (tPid1 == 99 && tPid2 == -57) { // Lambdabar
197  gPID = kLambdaBar;
198  tCharge = 0;
199  }
200  if (tPid1 == 2 && tPid2 == 0) { // proton
201  gPID = kProton;
202  tCharge = +1;
203  }
204  if (tPid1 == 15 && tPid2 == 1) { // SigmaPlus
205  gPID = kSigmaPlus;
206  tCharge = +1;
207  }
208  if (tPid1 == 15 && tPid2 == 0) { // Sigma0
209  gPID = kSigma0;
210  tCharge = 0;
211  }
212  if (tPid1 == 15 && tPid2 == -1) { // SigmaMinus
213  gPID = kSigmaMinus;
214  tCharge = -1;
215  }
216  if (tPid1 == 99 && tPid2 == -43) { // AntiSigmaPlus
217  gPID = kAntiSigmaPlus;
218  tCharge = -1;
219  }
220  if (tPid1 == 99 && tPid2 == -44) { // AntiSigma0
221  gPID = kAntiSigma0;
222  tCharge = 0;
223  }
224  if (tPid1 == 99 && tPid2 == -45) { // AntiSigmaMinus
225  gPID = kAntiSigmaMinus;
226  tCharge = +1;
227  }
228  if (tPid1 == 99 && tPid2 == -41) { // pbar
229  gPID = kPbar;
230  tCharge = -1;
231  }
232  if (tPid1 == 99 && tPid2 == 35) { // phi
233  gPID = kPhi;
234  tCharge = 0;
235  }
236  if (tPid1 == 7 && tPid2 == 0) { // pi-
237  gPID = kPiMinus;
238  tCharge = -1;
239  }
240  if (tPid1 == 9 && tPid2 == 0) { // pi+
241  gPID = kPiPlus;
242  tCharge = +1;
243  }
244  if (tPid1 == 8 && tPid2 == 0) { // pi0
245  gPID = kPi0;
246  tCharge = 0;
247  }
248  if (tPid1 == 99 && tPid2 == 46) { // Xi0
249  gPID = kXi0;
250  tCharge = 0;
251  }
252  if (tPid1 == 99 && tPid2 == 47) { // XiMinus
253  gPID = kXiMinus;
254  tCharge = -1;
255  }
256  if (tPid1 == 99 && tPid2 == -46) { // Anti-Xi0
257  gPID = kAntiXi0;
258  tCharge = 0;
259  }
260  if (tPid1 == 99 && tPid2 == -47) { // Anti-XiPlus
261  gPID = kAntiXiPlus;
262  tCharge = +1;
263  }
264  if (tPid1 == 99 && tPid2 == 70) { // Omega
265  gPID = kOmega;
266  tCharge = -1;
267  }
268  if (tPid1 == 99 && tPid2 == -70) { // Anti-Omega
269  gPID = kOmegaBar;
270  tCharge = +1;
271  }
272  if (0 == tCharge) {
273  continue;
274  }
275 
276  bool useTrack = true;
277  useTrack = (mTCuts->goodEta(eta) && useTrack);
278  useTrack = (mTCuts->goodPhi(phi) && useTrack);
279 
280  if (pt<0.15) continue;
281 
282  useTrack = (mTCuts->goodPt(pt) && useTrack);
283  float _r=pt/tM;
284  float yt=log(sqrt(1+_r*_r)+_r);
285  useTrack = (mTCuts->goodYt(yt) && useTrack);
286  mTCuts->fillHistograms(useTrack);
287  if (!useTrack) continue;
288  mnumTracks++;
289 
290  eTrack->SetBx(0);
291  eTrack->SetBy(0);
292  eTrack->SetBz(0);
293  eTrack->SetBxGlobal(0);
294  eTrack->SetByGlobal(0);
295  eTrack->SetBzGlobal(0);
296 
297  eTrack->SetPx(tPx);
298  eTrack->SetPy(tPy);
299  eTrack->SetPz(tPz);
300 
301  eTrack->SetEta(eta);
302  eTrack->SetPhi(phi);
303  eTrack->SetCharge(tCharge);
304 
305  estructEvent->AddTrack(eTrack);
306  }
307 
308  // Read first line of event header for next event.
309  // If we just read last event this will cause eof.
310  *mFile >> tDum; mLine++;
311  if (mFile->eof()) {
312  cout << "End of file. Closing it." << endl;
313  delete mFile;
314  mFile = NULL;
315  }
316 
317  delete eTrack;
318  return;
319 }
320 //--------------------------------------------------------------------------
321 float StEStructRQMD::getRapidity(float E, float pz) {
322  float y = 0.5*log((E + pz)/(E - pz));
323  return y;
324 }
325 
326 //--------------------------------------------------------------------------
327 float StEStructRQMD::getPseudoRapidity(float pt, float pz) {
328  float theta = fabs(atan(pt/pz));
329  float eta = -log(tan(theta/2.));
330  if (pz > 0) {
331  return eta;
332  } else {
333  return -eta;
334  }
335 }
336 
337 
338 
339 /**********************************************************************
340  *
341  * $Log: StEStructRQMD.cxx,v $
342  * Revision 1.4 2012/11/16 21:23:19 prindle
343  * EventCuts and TrackCuts were moved to EventReader. Remove that code from
344  * these readers.
345  *
346  * Revision 1.3 2006/04/06 01:03:35 prindle
347  *
348  * Rationalization of centrality binning, as described in AnalysisMaker checkin.
349  *
350  * Revision 1.2 2006/02/22 22:05:42 prindle
351  * Removed all references to multRef (?)
352  *
353  * Revision 1.1 2004/03/02 21:51:01 prindle
354  *
355  * I forgot to cvs add my EventGenerator readers.
356  *
357  * Revision 1.2 2003/11/25 22:45:14 prindle
358  * Commiting changes so I can move code to rhic
359  *
360  * Revision 1.1 2003/11/21 23:48:00 prindle
361  * Include my toy event generator in cvs
362  *
363  *
364  *********************************************************************/
Generated by   doxygen 1.8.5