StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtGstarTxtReader.cxx
1 /***************************************************************************
2  *
3  *
4  *
5  * Author: Mike Lisa, Ohio State, lisa@mps.ohio-state.edu
6  ***************************************************************************
7  *
8  * Description: part of STAR HBT Framework: StHbtMaker package
9  * This is the HbtEventReader class to be used when reading
10  * event-generator (e.g. mevsim) files in GSTAR text format
11  *
12  ***************************************************************************
13  *
14  *
15  **************************************************************************/
16 #define HBT_BFIELD 0.5*tesla
17 #include <Stiostream.h>
18 #include "Stiostream.h"
19 #include "StHbtMaker/Reader/StHbtGstarTxtReader.h"
20 #include "StHbtMaker/Base/StHbtEventCut.h"
21 #include "StHbtMaker/Base/StHbtTrackCut.h"
22 #include "StHbtMaker/Base/StHbtV0Cut.h"
23 #include "StHbtMaker/Base/StHbtKinkCut.h"
24 // hbt stuff
25 #include "SystemOfUnits.h" // has "tesla" in it
26 #include "StHbtMaker/Infrastructure/StHbtTrackCollection.hh"
27 #include "StHbtMaker/Infrastructure/StHbtV0Collection.hh"
28 #include "phys_constants.h"
29 #include "StPhysicalHelix.hh"
30 
31 #ifdef __ROOT__
32 ClassImp(StHbtGstarTxtReader)
33 #endif
34 
35 #if !(ST_NO_NAMESPACES)
36  using namespace units;
37 #endif
38 
39 /* =========== some useful functions in parsing the strings ============== */
40 int get_next_int(string strline,int start_at,int& wordends){
41  int wordstarts = strline.find_first_not_of(" ",start_at);
42  wordends = strline.find_first_of(" ",wordstarts);
43  int ians = atoi((strline.substr(wordstarts,wordends-1)).c_str());
44  return ians;
45 }
46 float get_next_float(string strline,int start_at, int& wordends){
47  int wordstarts = strline.find_first_not_of(" ",start_at);
48  wordends = strline.find_first_of(" ",wordstarts);
49  float ans = atof((strline.substr(wordstarts,wordends-1)).c_str());
50  return ans;
51 }
52 double dedxMean_geantTxt(double mass, double momentum){
53  double dedxMean;
54  double tpcDedxGain = 0.174325e-06;
55  double tpcDedxOffset = -2.71889;
56  double tpcDedxRise = 776.626;
57 
58  double gamma = ::sqrt(::pow(momentum/mass,2)+1.);
59  double beta = ::sqrt(1. - 1./::pow(gamma,2));
60  double rise = tpcDedxRise*::pow(beta*gamma,2);
61  if ( beta > 0)
62  dedxMean = tpcDedxGain/::pow(beta,2) * (0.5*::log(rise)-::pow(beta,2)- tpcDedxOffset);
63  else
64  dedxMean = 1000.;
65  return dedxMean;
66 }
67 /* ======================================================================== */
68 
69 
70 //_______________________________
71 StHbtGstarTxtReader::StHbtGstarTxtReader() : mInputStream(0){
72  mFileName = "GstarTextFile"; // default name
73  mReaderStatus = 0; // means "good"
74 }
75 
76 //_______________________________
77 StHbtGstarTxtReader::StHbtGstarTxtReader(char* file) : mInputStream(0), mFileName(file)
78 {
79  mReaderStatus = 0; // means "good"
80 }
81 
82 //_______________________________
83 StHbtGstarTxtReader::~StHbtGstarTxtReader(){
84  if (!mInputStream){
85  delete mInputStream;
86  mInputStream = 0;
87  }
88 }
89 
90 //_______________________________
91 StHbtEvent* StHbtGstarTxtReader::ReturnHbtEvent(){
92  if (!mInputStream){
93  cout << "StHbtGstarTxtReader::ReturnHbtEvent() - there is no input stream!";
94  mReaderStatus = 1; // 0 means "good"
95  return (0);
96  }
97  if (!(*mInputStream)){
98  cout << "StHbtGstarTxtReader::ReturnHbtEvent() - input stream in bad state!" << endl;
99  cout << "State is " << mInputStream->rdstate() << endl;
100  mReaderStatus = 1; // 0 means "good"
101  return (0);
102  }
103  // ok find the next EVENT keyword...
104  string strline;
105  char line[100];
106  string keyword;
107  int stringposition;
108 
109  //vvvvvvvvvvvvvvvvvvvv find and decode EVENT line vvvvvvvvvvvvvvvvvvv
110  keyword = "EVENT:";
111 
112  cout << "StHbtGstarTxtReader::ReturnHbtEvent() -- find and decode EVENT line..." << endl;
113  while (strline.substr(0,keyword.size()) != keyword){
114  (*mInputStream).getline(line,100);
115  strline = line;
116  if (!(mInputStream->good())){
117  cout << "StHbtGstarTxtReader::ReturnHbtEvent() - input stream in bad state!" <<endl;
118  cout << "State is " << mInputStream->rdstate() << endl;
119  mReaderStatus = 1; // 0 means "good"
120  return (0);
121  //mInputStream->clear();
122  }
123  }
124 
125 
126  // event number...
127  int ievent = get_next_int(strline,keyword.size()+1,stringposition);
128  int ntracks = get_next_int(strline,stringposition,stringposition);
129  int nvertices = get_next_int(strline,stringposition,stringposition);
130  cout << "Event number is " << ievent << endl;
131  cout << "Number of tracks is " << ntracks << endl;
132  cout << "Number of vertices is " << nvertices << endl;
133 
134  if (ntracks <= 0){ //
135  cout << "StHbtGstarTxtReader::ReturnHbtEvent() - hit end-of-file " << endl;
136  cout << "State is " << mInputStream->rdstate() << endl;
137  mReaderStatus = 1; // 0 means "good"
138  return (0);
139  }
140  //^^^^^^^^^^^^^^^^^^^^^^^^ done finding and decoding EVENT line ^^^^^^^^^^^^^^^^
141 
142  StHbtEvent* event = new StHbtEvent;
143  const StHbtThreeVector vertexPos(0.0,0.0,0.0);
144 
145  event->SetEventNumber(ievent);
146  event->SetCtbMult(0);
147  event->SetZdcAdcEast(0);
148  event->SetZdcAdcWest(0);
149  event->SetNumberOfTpcHits(0);
150  event->SetNumberOfTracks(ntracks);
151  event->SetNumberOfGoodTracks(ntracks); // same for now
152  event->SetReactionPlane(0.);
153  event->SetReactionPlaneSubEventDifference(0.);
154  event->SetPrimVertPos(vertexPos);
155 
156 
157  // now decode TRACK lines...
158  int pid;
159  float px,py,pz;
160  StHbtThreeVector p;
161  keyword = "TRACK:";
162  for (int itrack=0; itrack<ntracks; itrack++){
163  strline = " ";
164  while (strline.substr(0,keyword.size()) != keyword){
165  (*mInputStream).getline(line,100);
166  strline = line;
167  if (!(mInputStream->good())){
168  cout << "StHbtGstarTxtReader::ReturnHbtEvent() - input stream in bad state!" <<endl;
169  cout << "State is " << mInputStream->rdstate() << endl;
170  mReaderStatus = 1; // 0 means "good"
171  return (0);
172  //mInputStream->clear();
173  }
174  }
175 
176 
177  pid = get_next_int(strline,keyword.size()+1,stringposition);
178  px = get_next_float(strline,stringposition,stringposition);
179  py = get_next_float(strline,stringposition,stringposition);
180  pz = get_next_float(strline,stringposition,stringposition);
181 
182  int charge = -1;
183  float mass=0.0;
184 
185  StHbtTrack* hbtTrack = new StHbtTrack;
186 
187  // cout << "Pid is " << pid;
188  hbtTrack->SetTrackId(itrack);
189 
190  switch (pid){
191  case 2: // intentional fall-through
192  charge = 1;
193  case 3: // pion
194  hbtTrack->SetNSigmaElectron(0.);
195  hbtTrack->SetNSigmaPion(999.);
196  hbtTrack->SetNSigmaKaon(-999.);
197  hbtTrack->SetNSigmaProton(-999.);
198  mass = 5.1099905E-04;
199  break;
200  case 8: // intentional fall-through
201  charge = 1;
202  case 9: // pion
203  hbtTrack->SetNSigmaElectron(999.);
204  hbtTrack->SetNSigmaPion(0.);
205  hbtTrack->SetNSigmaKaon(-999.);
206  hbtTrack->SetNSigmaProton(-999.);
207  mass = 0.139567;
208  break;
209  case 11: // intensional fall-thru
210  charge = 1;
211  case 12:
212  hbtTrack->SetNSigmaElectron(999.);
213  hbtTrack->SetNSigmaPion(999.0);
214  hbtTrack->SetNSigmaKaon(0.);
215  hbtTrack->SetNSigmaProton(-999.);
216  mass = 0.493667;
217  break;
218  case 14: // proton
219  charge = 1;
220  hbtTrack->SetNSigmaElectron(999.);
221  hbtTrack->SetNSigmaPion(999.);
222  hbtTrack->SetNSigmaKaon(999.);
223  hbtTrack->SetNSigmaProton(0.);
224  mass = 0.93828;
225  break;
226  default:
227  // cout << "Non-recognized track -- pid = " << pid << endl;
228  charge = -99;
229  mass = 0.0;
230  break;
231  }
232 
233  if (charge == -99){
234  delete hbtTrack;
235  // cout << "Removing StHbtTrack" << endl;
236  continue; // do NOT make a StHbtTrack-- continue with loop...
237  }
238 
239 
240  StHbtThreeVector tmpP;
241  tmpP.setX(px);
242  tmpP.setY(py);
243  tmpP.setZ(pz);
244  hbtTrack->SetP( tmpP );
245 
246  // cout << "Pid is " << pid << " and momentum is " << tmpP << endl;
247 
248  // place-holder
249  hbtTrack->SetNHits(45);
250  hbtTrack->SetNHitsPossible(45);
251  hbtTrack->SetdEdx( dedxMean_geantTxt( mass, tmpP.mag() ) );
252  hbtTrack->SetPt( tmpP.perp());
253  hbtTrack->SetCharge(charge);
254 
255  StPhysicalHelixD helix = StPhysicalHelixD( hbtTrack->P(), vertexPos, HBT_BFIELD, hbtTrack->Charge() );
256  hbtTrack->SetHelix(helix);
257 
258  hbtTrack->SetDCAxy(0.001);
259  hbtTrack->SetDCAz(0.001);
260  hbtTrack->SetChiSquaredXY( 0.);
261  hbtTrack->SetChiSquaredZ( 0.);
262 
263  event->TrackCollection()->push_back(hbtTrack);
264  }
265  cout << event->TrackCollection()->size() << " tracks pushed to collection" << endl;
266 
267  return event;
268 }
269 
270 //_______________________________
271 StHbtString StHbtGstarTxtReader::Report(){
272  StHbtString temp = "\n This is the StHbtGstarTxtReader - no Early Cuts applied\n";
273  temp += " *** NOTE I am kinda stupid-- I do NOT handle vertices, and I ONLY\n";
274  temp += " *** know about pions protons and kaons\n";
275  return temp;
276 }
277 
278 
279 //_______________________________
280 int StHbtGstarTxtReader::Init(const char* ReadWrite, StHbtString& Message){
281  cout << " *\n *\n *\n StHbtGstarTxtReader::Init() being called*\n *\n";
282  mReaderStatus = 0; // means "good"
283  // if ((ReadWrite=="r")|| (ReadWrite=="R")){ // this object will be a reader
284  if (((*ReadWrite)=='r')|| ((*ReadWrite)=='R')){ // this object will be a reader
285  mInputStream = new ifstream;
286  mInputStream->open(mFileName);
287  if (!(*mInputStream)){
288  cout << "StHbtGstarTxtReader::Init() - Cannot open input file! " << endl;
289  return (1);
290  }
291  cout << "StHbtGstarTxtReader::Init() - being configured as a Reader - " << ReadWrite << endl;
292  }
293  else{ // this object will be a writer
294  cout << " CANNOT BE CONFIGURED AS A WRITER" << endl;
295  return (1);
296  }
297  return (0);
298 }
299 
300 //_______________________________
301 void StHbtGstarTxtReader::Finish(){
302  if (mInputStream) mInputStream->close();
303 }
304 
Generated by   doxygen 1.8.5