StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructPythia.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructPythia.cxx,v 1.15 2012/11/16 21:23:18 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: EStructEventReader wrapper for (T)Pythia event generator
10  *
11  **********************************************************************/
12 #include "StEStructPythia.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/StEStructTrack.h"
19 #include "StEStructPool/EventMaker/StEStructCentrality.h"
20 
21 StEStructPythia::StEStructPythia() {
22  mPythia = 0;
23  mECuts = 0;
24  mTCuts = 0;
25  mInChain = false;
26  museAllTracks = false;
27  mEventsToDo = 100;
28  mEventCount = 0;
29  mAmDone = false;
30 };
31 
32 StEStructPythia::StEStructPythia(TPythia6* pythia,
33  StEStructEventCuts* ecuts,
34  StEStructTrackCuts* tcuts,
35  bool useAllTracks,
36  int eventsToDo) {
37  mPythia = pythia;
38  mECuts = ecuts;
39  mTCuts = tcuts;
40  mInChain = false;
41  museAllTracks = useAllTracks;
42  mEventsToDo = eventsToDo;
43  mEventCount = 0;
44  mAmDone = false;
45 };
46 
47 bool StEStructPythia::hasGenerator() { return (mPythia) ? true : false ; };
48 
49 //-------------------------------------------------------------------------
50 StEStructEvent* StEStructPythia::next() {
51 
52  if (!mPythia || mEventCount==mEventsToDo) {
53  mAmDone = true;
54  return (StEStructEvent*)NULL;
55  }
56  StEStructEvent* retVal = NULL;
57 
58  if (!mInChain ) {
59  mPythia->GenerateEvent();
60  mEventCount++;
61  }
62 
63  retVal = new StEStructEvent();
64 
65  retVal->SetVx(0);
66  retVal->SetVy(0);
67  retVal->SetVz(0);
68  retVal->SetBField(0.5);
69 
70  mstarTrigger=false;
71  int nTracks = countGoodTracks();
72 
73  float centMeasure;
74  if (museAllTracks) {
75  centMeasure = getNPartonic();
76  } else {
77  centMeasure = nTracks;
78  }
79  retVal->SetCentrality(centMeasure);
80  if(!mECuts->goodCentrality(centMeasure) || !mstarTrigger){
81  delete retVal;
82  retVal=NULL;
83  mECuts->fillHistogram(mECuts->centralityName(),(float)nTracks,false);
84  } else {
85  fillTracks(retVal);
86  retVal->FillChargeCollections();
87  mECuts->fillHistogram(mECuts->centralityName(),(float)nTracks,true);
88  }
89 
90  return retVal;
91 }
92 
93 //--------------------------------------------------------------------------
94 double StEStructPythia::getNPartonic(){
95  if(!mPythia) return 0.;
96  return (double)mPythia->GetN();
97 }
98 
99 //--------------------------------------------------------------------------
100 void StEStructPythia::fillTracks(StEStructEvent* estructEvent){
101 
102  mnumTracks=0;
103  Pyjets_t* pstr= mPythia->GetPyjets();
104  int numParticles=mPythia->GetN();
105 
106  StEStructTrack* eTrack= new StEStructTrack();
107  for (int i=2;i<numParticles;i++) { // 0 & 1 for incoming protons
108 
109  if (pstr->K[0][i]==21) continue;
110  int pid = pstr->K[1][i];
111  // require no daughters & if pi,k,proton, or electron; else skip...
112  if (!(0==pstr->K[3][i])) continue;
113 
114  if (!isTrackGood(i)) {
115  mTCuts->fillHistograms(false);
116  continue;
117  }
118  mTCuts->fillHistograms(true);
119  eTrack->SetInComplete();
120 
121  float p[3];
122  float v[3];
123  for(int k=0;k<3;k++){
124  p[k]=pstr->P[k][i];
125  v[k]=pstr->V[k][i];
126  }
127 
128  float pt = sqrt(p[0]*p[0] + p[1]*p[1]);
129  float theta = acos(p[2]/ TMath::Sqrt(pt*pt + p[2]*p[2]) );
130  float eta = -1.0*log(tan(theta/2.0));
131  float *gdca = globalDCA(p,v);
132  float phi=atan2((double)p[1], (double)p[0]);
133 
134  eTrack->SetBx(gdca[0]);
135  eTrack->SetBy(gdca[1]);
136  eTrack->SetBz(gdca[2]);
137  eTrack->SetBxPrimary(gdca[0]);
138  eTrack->SetByPrimary(gdca[1]);
139  eTrack->SetBzPrimary(gdca[2]);
140  eTrack->SetBxGlobal(gdca[0]);
141  eTrack->SetByGlobal(gdca[1]);
142  eTrack->SetBzGlobal(gdca[2]);
143  delete [] gdca;
144 
145  eTrack->SetPIDe_dEdx(10);
146  eTrack->SetPIDpi_dEdx(10);
147  eTrack->SetPIDk_dEdx(10);
148  eTrack->SetPIDp_dEdx(10);
149  eTrack->SetPIDd_dEdx(10);
150  eTrack->SetPIDe_ToF(10);
151  eTrack->SetPIDpi_ToF(10);
152  eTrack->SetPIDk_ToF(10);
153  eTrack->SetPIDp_ToF(10);
154  eTrack->SetPIDd_ToF(10);
155  if ((pid == 7) || (pid == 8)) {
156  eTrack->SetPIDe_dEdx(0);
157  eTrack->SetPIDe_ToF(0);
158  } else if ((pid == -211) || (pid == 211)) {
159  eTrack->SetPIDpi_dEdx(0);
160  eTrack->SetPIDpi_ToF(0);
161  } else if ((pid == -321) || (pid == 321)) {
162  eTrack->SetPIDk_dEdx(0);
163  eTrack->SetPIDk_ToF(0);
164  } else if ((pid == -2212) || (pid == 2212)) {
165  eTrack->SetPIDk_dEdx(0);
166  eTrack->SetPIDk_ToF(0);
167  } else if (pid ==95) {
168  // No anti-deuteron. I think Hijing does not make deuterons and
169  // this pid code is only intended for use with GEANT?
170  eTrack->SetPIDd_dEdx(0);
171  eTrack->SetPIDd_ToF(0);
172  }
173 
174  eTrack->SetPx(p[0]);
175  eTrack->SetPy(p[1]);
176  eTrack->SetPz(p[2]);
177  eTrack->SetEta(eta);
178  eTrack->SetPhi(phi);
179  eTrack->SetDedx(0);
180  eTrack->SetChi2(1);
181  eTrack->SetTopologyMapTPCNHits(45);
182  eTrack->SetNMaxPoints(45);
183  eTrack->SetNFoundPoints(45);
184  eTrack->SetNFitPoints(45);
185 
186  if(pid<0){
187  eTrack->SetCharge(-1);
188  } else {
189  eTrack->SetCharge(1);
190  }
191 
192 
193  // I don't know what this next block of code does. Jeff must have put this in.
194  // djp: Sept. 13, 2005
195 
196  // rjp: Feb 2006, Ok, I'll tell you...
197  // ip array is the parent history
198  // now add fragmentation history (up to 4 lines back) in the tpcmap area
199  int ip[4] = {0,0,0,0};
200  ip[0] = pstr->K[2][i];
201  for (int k=1;k<4;k++) {
202  if (ip[k-1]<3) {
203  ip[k-1]=0;
204  break;
205  }
206  ip[k] = pstr->K[2][ip[k-1]-1];
207  }
208  unsigned int map[2];
209  map[0]=(unsigned int)((ip[1]<<16)+ip[0]);
210  map[1]=(unsigned int)((ip[3]<<16)+ip[2]);
211  eTrack->SetTopologyMapData(0,map[0]);
212  eTrack->SetTopologyMapData(1,map[1]);
213 
214  // I included this line because I have in in the Hijing interface
215  // but I don't know what it does right off. djp: Sept. 13, 2005
216  eTrack->SetTopologyMapTPCNHits(45);
217 
218  estructEvent->AddTrack(eTrack);
219  }; // particle loop
220 
221  delete eTrack;
222  return;
223 
224 }
225 
226 
227 //-------------------------------------------------------------
228 // This method checks all track cuts.
229 // No histogramming or copying data around.
230 bool StEStructPythia::isTrackGood(int i) {
231  Pyjets_t* pstr= mPythia->GetPyjets();
232  int pid = pstr->K[1][i];
233  if (!measureable(pid)) { // checks if pi,k,p or e
234  return false;
235  }
236  float p[3];
237  float v[3];
238  for (int k=0;k<3;k++) {
239  p[k] = pstr->P[k][i];
240  v[k] = pstr->V[k][i];
241  }
242 
243  float pt = sqrt(p[0]*p[0]+p[1]*p[1]);
244  if (pt < 0.15) {
245  return false;
246  }
247 
248  float theta = acos(p[2]/ sqrt(pt * pt + p[2] * p[2]) );
249  float eta = -1.0*log(tan(theta/2.0));
250  float phi = atan2((double)p[1], (double)p[0]);
251  float* gdca = globalDCA(p,v);
252  float _r = pt/0.139;
253  float yt = log(sqrt(1+_r*_r)+_r);
254 
255  bool useTrack = true;
256  useTrack = (mTCuts->goodGlobalDCA(gdca[3]) && useTrack);
257  useTrack = (mTCuts->goodEta(eta) && useTrack);
258  useTrack = (mTCuts->goodPhi(phi) && useTrack);
259  useTrack = (mTCuts->goodPt(pt) && useTrack);
260  useTrack = (mTCuts->goodYt(yt) && useTrack);
261  delete [] gdca;
262 
263  if(eta<-3.5 && eta>-5.0) mstarTrigger=true;
264  if(eta>3.5 && eta<5.0) mstarTrigger=true;
265 
266  return useTrack;
267 }
268 
269 //-------------------------------------------------------------------------
270 //-------------------------------------------------------------
271 // This method counts all good track.
272 // No histogramming or copying data around.
273 int StEStructPythia::countGoodTracks() {
274  mnumTracks = 0;
275  int numParticles = mPythia->GetN();
276  for (int i=2;i<numParticles;i++) { // 0 & 1 for incoming protons
277  if (isTrackGood(i)) {
278  mnumTracks++;
279  }
280  }
281  return mnumTracks;
282 }
283 
284 
285 
286 /**********************************************************************
287  *
288  * $Log: StEStructPythia.cxx,v $
289  * Revision 1.15 2012/11/16 21:23:18 prindle
290  * EventCuts and TrackCuts were moved to EventReader. Remove that code from
291  * these readers.
292  *
293  * Revision 1.14 2010/09/02 21:24:45 prindle
294  * Pythia: Fill in ToF pid information
295  *
296  * Revision 1.13 2010/03/02 21:46:24 prindle
297  * Option to use getNPartonic as a centrality measure
298  *
299  * Revision 1.12 2006/10/02 22:22:50 prindle
300  * Changed the values stored in PID for different particle types so I can
301  * use same cuts for Pythia (and Hijing) as I do for data.
302  *
303  * Revision 1.11 2006/04/11 17:51:39 prindle
304  * Remove inChain from constructor arguments (no longer used in macro)
305  *
306  * Revision 1.10 2006/04/06 01:03:32 prindle
307  *
308  * Rationalization of centrality binning, as described in AnalysisMaker checkin.
309  *
310  * Revision 1.9 2006/04/04 22:11:27 porter
311  * StEStructPythia now uses StEtructCentrality for selection
312  *
313  * Revision 1.8 2006/02/22 22:05:37 prindle
314  * Removed all references to multRef (?)
315  *
316  * Revision 1.7 2005/09/23 23:37:23 prindle
317  *
318  * Starting to add vertex distribution and track acceptance dependance on
319  * number of possible hits.
320  * Make Pythia interface look like Hijing interface so it now works within
321  * my Fluctuation and Correlation framework.
322  *
323  * Revision 1.6 2004/09/24 01:43:12 prindle
324  * Add call to define centrality by multiplicity.
325  *
326  * Revision 1.5 2004/07/01 00:35:29 porter
327  * modified 'startrigger'... still need to make this a switch
328  *
329  * Revision 1.4 2004/06/25 03:13:01 porter
330  * added simple trigger selection implemented like BBC-AND plus CTB
331  *
332  * Revision 1.3 2004/06/09 22:37:51 prindle
333  * Moved some variable declarations inside comment to get rid of
334  * unused variable warnings.
335  * ved
336  *
337  * Revision 1.2 2004/03/24 21:26:52 porter
338  * fixed error calculating eta as rapidity
339  *
340  * Revision 1.1 2003/11/21 06:24:56 porter
341  * Pythia event generater as an StEStructEventReader
342  *
343  *
344  *********************************************************************/
Pyjets_t
Definition: P2ATest.C:19
Generated by   doxygen 1.8.5