StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructHijing.cxx
1 int sortCompare ( const void * elem1, const void * elem2 ) {
2  int *e1 = (int *) elem1;
3  int *e2 = (int *) elem2;
4  return (e1[1] - e2[1]);
5 }
6 int rSortCompare ( const void * elem1, const void * elem2 ) {
7  double *e1 = (double *) elem1;
8  double *e2 = (double *) elem2;
9  return (e1[2] - e2[2]);
10 }
11 /**********************************************************************
12  *
13  * $Id: StEStructHijing.cxx,v 1.10 2012/11/16 21:28:32 prindle Exp $
14  *
15  * Author: Chunhui Han
16  *
17  **********************************************************************
18  *
19  * Description: EStructEventReader wrapper for (T)Hijing event generator
20  *
21  **********************************************************************/
22 #include <stdlib.h>
23 #include "TRandom3.h"
24 
25 #include "StEStructHijing.h"
26 
27 #include "StEStructPool/AnalysisMaker/StEStructEventCuts.h"
28 #include "StEStructPool/AnalysisMaker/StEStructTrackCuts.h"
29 #include "StEStructPool/EventMaker/StEStructEvent.h"
30 #include "StEStructPool/EventMaker/StEStructTrack.h"
31 #include "StEStructPool/EventMaker/StEStructCentrality.h"
32 
33 StEStructHijing::StEStructHijing() {
34  mHijing = 0;
35  mInChain = false;
36  mEventsToDo = 100;
37  museImpactParameter = true;
38  mEventCount = 0;
39  mAmDone = false;
40  mTrackList = 0;
41  rTrackList = 0;
42 };
43 
44 StEStructHijing::StEStructHijing(THijing* hijing,
45  StEStructEventCuts* ecuts,
46  StEStructTrackCuts* tcuts,
47  bool useImpactParameter,
48  int eventsToDo) {
49  mHijing = hijing;
50  mInChain = false;
51  mEventsToDo = eventsToDo;
52  museImpactParameter = useImpactParameter;
53  mEventCount = 0;
54  mAmDone = false;
55  mTrackList = 0;
56  rTrackList = 0;
57 };
58 
59 bool StEStructHijing::hasGenerator() { return (mHijing) ? true : false; }
60 
61 //-------------------------------------------------------------------------
62 StEStructEvent* StEStructHijing::next() {
63  if(!mHijing || mEventCount== mEventsToDo) {
64  mAmDone=true;
65  return (StEStructEvent*)NULL;
66  }
67  StEStructEvent* retVal=NULL;
68 
69  if (!mInChain ) {
70  mHijing->GenerateEvent();
71  mEventCount++;
72  }
73 // int rej = filterTracks();
74 // int rej = filterTrackArea();
75 
76  retVal = new StEStructEvent();
77 
78  retVal->SetVx(0);
79  retVal->SetVy(0);
80  retVal->SetVz(0);
81  retVal->SetBField(0.5);
82 
83  int nTracks = countGoodTracks();
84  float centMeasure;
85  mImpact = mHijing->GetImpactParameter();
86  if (museImpactParameter) {
87  centMeasure = mImpact;
88  } else {
89  centMeasure = nTracks;
90  }
91  retVal->SetCentrality(centMeasure);
92  if (!mECuts->goodCentrality(centMeasure)) {
93  delete retVal;
94  retVal=NULL;
95  mECuts->fillHistogram(mECuts->centralityName(),centMeasure,false);
96  } else {
97  fillTracks(retVal);
98  retVal->FillChargeCollections();
99  mECuts->fillHistogram(mECuts->centralityName(),centMeasure,true);
100  }
101  return retVal;
102 }
103 
104 //--------------------------------------------------------------------------
105 void StEStructHijing::fillTracks(StEStructEvent* estructEvent) {
106  int numParticles=mHijing->GetNParticles();
107 
108  StEStructTrack* eTrack= new StEStructTrack();
109  for (int i=0;i<numParticles;i++) {
110  if (!isTrackGood(i)) {
111  mTCuts->fillHistograms(false);
112  continue;
113  }
114  mTCuts->fillHistograms(true);
115  eTrack->SetInComplete();
116 
117  int pid = mHijing->GetPdg(i);
118  float p[3];
119  float v[3];
120  for(int k=0;k<3;k++){
121  p[0] = mHijing->GetPx(i);
122  p[1] = mHijing->GetPy(i);
123  p[2] = mHijing->GetPz(i);
124  v[0] = mHijing->GetVx(i);
125  v[1] = mHijing->GetVy(i);
126  v[2] = mHijing->GetVz(i);
127  }
128 
129  float pt = sqrt(p[0]*p[0] + p[1]*p[1]);
130  float theta = acos(p[2]/ TMath::Sqrt(pt*pt + p[2]*p[2]) );
131  float eta = -1.0*log(tan(theta/2.0));
132  float *gdca = globalDCA(p,v);
133  float phi=atan2((double)p[1], (double)p[0]);
134 
135  eTrack->SetBx(gdca[0]);
136  eTrack->SetBy(gdca[1]);
137  eTrack->SetBz(gdca[2]);
138  eTrack->SetBxPrimary(gdca[0]);
139  eTrack->SetByPrimary(gdca[1]);
140  eTrack->SetBzPrimary(gdca[2]);
141  eTrack->SetBxGlobal(gdca[0]);
142  eTrack->SetByGlobal(gdca[1]);
143  eTrack->SetBzGlobal(gdca[2]);
144  delete [] gdca;
145 
146  eTrack->SetPIDe_dEdx(10);
147  eTrack->SetPIDpi_dEdx(10);
148  eTrack->SetPIDk_dEdx(10);
149  eTrack->SetPIDp_dEdx(10);
150  eTrack->SetPIDd_dEdx(10);
151  eTrack->SetPIDe_ToF(10);
152  eTrack->SetPIDpi_ToF(10);
153  eTrack->SetPIDk_ToF(10);
154  eTrack->SetPIDp_ToF(10);
155  eTrack->SetPIDd_ToF(10);
156  if ((pid == 7) || (pid == 8)) {
157  eTrack->SetPIDe_dEdx(0);
158  eTrack->SetPIDe_ToF(0);
159  } else if ((pid == -211) || (pid == 211)) {
160  eTrack->SetPIDpi_dEdx(0);
161  eTrack->SetPIDpi_ToF(0);
162  } else if ((pid == -321) || (pid == 321)) {
163  eTrack->SetPIDk_dEdx(0);
164  eTrack->SetPIDk_ToF(0);
165  } else if ((pid == -2212) || (pid == 2212)) {
166  eTrack->SetPIDp_dEdx(0);
167  eTrack->SetPIDp_ToF(0);
168  } else if ((pid == -2213) || (pid == 2213)) {
169  // These numbers aren't right, but I don't care about dueterons right now.
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->SetBeta(0);
180  eTrack->SetDedx(0);
181  eTrack->SetChi2(1);
182  eTrack->SetTopologyMapData(0, 0xffffff80);
183  eTrack->SetTopologyMapData(1, 0x3fff);
184  eTrack->SetTopologyMapTPCNHits(45);
185  eTrack->SetNMaxPoints(45);
186  eTrack->SetNFoundPoints(45);
187  eTrack->SetNFitPoints(45);
188 
189  if(pid<0){
190  eTrack->SetCharge(-1);
191  } else {
192  eTrack->SetCharge(1);
193  }
194 
195  estructEvent->AddTrack(eTrack);
196  }
197  delete eTrack;
198  return;
199 }
200 
201 //-------------------------------------------------------------
202 // This method checks all track cuts.
203 // No histogramming or copying data around.
204 bool StEStructHijing::isTrackGood(int it) {
205  int i = it;
206  bool useTrack = true;
207  if (mTrackList) {
208  i = mTrackList[2*it];
209  if (mTrackList[2*it+1] < 0) { // reject because of two-track efficiency
210  useTrack = false;
211  }
212  }
213  if (rTrackList) {
214  i = (int) rTrackList[4*it];
215  if (rTrackList[4*it+3] < 0) { // reject because of general track efficiency
216  useTrack = false;
217  }
218  }
219  int pid = mHijing->GetPdg(i);
220  if (!measureable(pid)) { // checks if pi,k,p or e
221  useTrack = false;
222  }
223  float p[3];
224  float v[3];
225  p[0] = mHijing->GetPx(i);
226  p[1] = mHijing->GetPy(i);
227  p[2] = mHijing->GetPz(i);
228  v[0] = mHijing->GetVx(i);
229  v[1] = mHijing->GetVy(i);
230  v[2] = mHijing->GetVz(i);
231 
232  int q;
233  pid < 0 ? q=-1 : q=+1;
234  float pt = sqrt(p[0]*p[0]+p[1]*p[1]);
235  float theta = acos(p[2]/ TMath::Sqrt(pt * pt + p[2] * p[2]) );
236  float eta = -1.0*log(tan(theta/2.0));
237  float phi = atan2((double)p[1], (double)p[0]);
238  float* gdca = globalDCA(p,v);
239  float _r = pt/0.139;
240  float yt = log(sqrt(1+_r*_r)+_r);
241  int iFrag = mHijing->GetOrigin(i);
242 
243  useTrack = (mTCuts->goodGlobalDCA(gdca[3]) && useTrack);
244  useTrack = (mTCuts->goodCharge(q) && useTrack);
245  useTrack = (mTCuts->goodPt(pt) && useTrack);
246  useTrack = (mTCuts->goodEta(eta) && useTrack);
247  useTrack = (mTCuts->goodPhi(phi) && useTrack);
248  useTrack = (mTCuts->goodYt(yt) && useTrack);
249  if (useTrack) {
250  useTrack = mTCuts->goodFragment(iFrag);
251  }
252  delete [] gdca;
253 
254  return useTrack;
255 }
256 
257 //-------------------------------------------------------------------------
258 //-------------------------------------------------------------
259 // This method bins tracks in (\eta,\phi) bins and applies track
260 // efficiencies to bins with more than one track.
261 int StEStructHijing::filterTracks() {
262  int nRejBadEta = 0;
263  int nRejNext = 0;
264  int nRejLast = 0;
265  if (mTrackList) {
266  delete [] mTrackList;
267  }
268  double pi = acos(-1);
269  int numParticles = mHijing->GetNParticles();
270 
271  // Determine bin number in (\eta,\phi) for each track.
272  mTrackList = new int[2*numParticles];
273  for (int i=0;i<numParticles;i++) {
274  double px = mHijing->GetPx(i);
275  double py = mHijing->GetPy(i);
276  double pz = mHijing->GetPz(i);
277  double phi = atan2(py,px);
278  if (phi < 0) {
279  phi += 2*pi;
280  }
281  int iphi = 150 * phi / (2*pi);
282  double theta = acos(pz/ TMath::Sqrt(px*px + py*py + pz*pz) );
283  double eta = -1.0*log(tan(theta/2.0));
284  int ieta = 75 * (1 + eta) / 2;
285  int ibin = -1;
286  if (0 <= ieta && ieta < 75) {
287  // Only use tracks in -1 <= eta <= 1
288  ibin = 150 * ieta + iphi;
289  }
290  mTrackList[2*i] = i;
291  mTrackList[2*i+1] = ibin;
292  }
293 
294  // Sort by bin number
295  qsort(mTrackList, numParticles, 2*sizeof(int), sortCompare);
296 
297  // Randomly reject tracks from bins that have more than one.
298  TRandom3 ran;
299  for (int i=0;i<numParticles-1;i++) {
300  if (mTrackList[2*i+1] < 0) {
301  mTrackList[2*i] = -1;
302  nRejBadEta++;
303  } else if (mTrackList[2*i+1] == mTrackList[2*i+3]) {
304  if (ran.Rndm() > 1.0) {
305  mTrackList[2*i+1] *= -1;
306  nRejNext++;
307  }
308  } else if (mTrackList[2*i+1] == mTrackList[2*i-1]) {
309  if (ran.Rndm() > 1.0) {
310  mTrackList[2*i+1] *= -1;
311  nRejLast++;
312  }
313  }
314  }
315  return nRejNext + nRejLast;
316 }
317 //-------------------------------------------------------------------------
318 //-------------------------------------------------------------
319 // This method applies tracking efficiency depending on track density in
320 // a broader area than filterTracks (that is essentially a two-track resolution).
321 // Sort list in eta. Step through tracks both toward +phi and -phi until we
322 // exceed DeltaPhi. For each track see if we are within DeltaEta. Count total
323 // number, apply efficiency and flag if it fails.
324 int StEStructHijing::filterTrackArea() {
325  int nRejEta = 0;
326  int nRejEff = 0;
327  if (rTrackList) {
328  delete [] rTrackList;
329  }
330  double pi = acos(-1);
331  int numParticles = mHijing->GetNParticles();
332 
333  // Determine bin number in (\eta,\phi) for each track.
334  rTrackList = new double[4*numParticles];
335  for (int i=0;i<numParticles;i++) {
336  double px = mHijing->GetPx(i);
337  double py = mHijing->GetPy(i);
338  double pz = mHijing->GetPz(i);
339  double phi = atan2(py,px);
340  double theta = acos(pz/ TMath::Sqrt(px*px + py*py + pz*pz) );
341  double eta = -1.0*log(tan(theta/2.0));
342  rTrackList[4*i] = i;
343  rTrackList[4*i+1] = phi;
344  rTrackList[4*i+2] = eta;
345  rTrackList[4*i+3] = 1;
346  }
347 
348  // Sort by bin number
349  qsort(rTrackList, numParticles, 4*sizeof(double), rSortCompare);
350 
351  // Consider each track.
352  // Count number of other tracks in DeltaPhi, DeltaEta
353  double DeltaEta = 2.0/15;
354  double DeltaPhi = 2.0*pi/15;
355  int j;
356  TRandom3 ran;
357  for (int i=0;i<numParticles;i++) {
358  int iDense = 0;
359  double eta = rTrackList[4*i+2];
360  if (fabs(eta)-1 > DeltaEta) {
361  nRejEta++;
362  continue;
363  }
364  double phi = rTrackList[4*i+1];
365  double dPhi;
366  j = i + 1;
367  while ((j < numParticles) && (rTrackList[4*j+2]-eta < DeltaEta)) {
368  // Think this works. Is it faster?
369  dPhi = fmod(fabs(phi-rTrackList[4*j+1])+pi,2*pi) - pi;;
370  if (fabs(dPhi) < DeltaPhi) {
371  iDense++;
372  }
373  j++;
374  }
375  j = i - 1;
376  while ((j >= 0) && (rTrackList[4*j+2]-eta < DeltaEta)) {
377  dPhi = fmod(fabs(phi-rTrackList[4*j+1])+pi,2*pi) - pi;;
378  if (fabs(dPhi) < DeltaPhi) {
379  iDense++;
380  }
381  j--;
382  }
383  // Local track density is now iDense / (DeltaEta*DeltaPhi)
384  double eff = 0.9 - 0.0025*iDense/(DeltaEta*DeltaPhi);
385  eff = 2.0;
386  if (ran.Rndm() > eff) {
387  rTrackList[4*i+3] = -1;
388  nRejEff++;
389  }
390  }
391  return nRejEff;
392 }
393 //--------------------------------------------------------------------------
394 //-------------------------------------------------------------------------
395 //-------------------------------------------------------------
396 // This method counts all good track.
397 // No histogramming or copying data around.
398 int StEStructHijing::countGoodTracks() {
399  mnumTracks = 0;
400  int numParticles = mHijing->GetNParticles();
401  for (int i=0;i<numParticles;i++) {
402  if (isTrackGood(i)) {
403  mnumTracks++;
404  }
405  }
406  return mnumTracks;
407 }
408 void StEStructHijing::setHijingReader(THijing* hijing){
409  if (mHijing) delete mHijing;
410  mHijing = hijing;
411 };
412 
413 /**********************************************************************
414  *
415  *********************************************************************/
Generated by   doxygen 1.8.5