StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtMcEventReader.cxx
1 #define HBT_BFIELD 0.25*tesla
2 #define DIFF_CUT_OFF 1.
3 
4 #include "StHbtMaker/Reader/StHbtMcEventReader.h"
5 #include "StChain.h"
6 
7 // StEvent stuff
8 #include "StEventTypes.h"
9 #include "StEventMaker/StEventMaker.h"
10 // StMcEvent stuff
11 #include "StMcEventTypes.hh"
12 #include "StMcEventMaker/StMcEventMaker.h"
13 
14 // c++ stuff
15 //#include <typeinfo>
16 #include <math.h>
17 
18 // hbt stuff
19 #include "SystemOfUnits.h" // has "tesla" in it
20 #include "StHbtMaker/Infrastructure/StHbtTrackCollection.hh"
21 #include "StHbtMaker/Infrastructure/StHbtV0Collection.hh"
22 
23 #include <Stiostream.h>
24 #include <stdlib.h>
25 #include <string>
26 #include <vector>
27 
28 #include "PhysicalConstants.h"
29 #include "SystemOfUnits.h"
30 #include "StThreeVector.hh"
31 #include "phys_constants.h"
32 
33 #include "StChain.h"
34 #include "St_DataSet.h"
35 #include "St_DataSetIter.h"
36 
37 #include "g2t_event.h"
38 #include "g2t_ftp_hit.h"
39 #include "g2t_svt_hit.h"
40 #include "g2t_tpc_hit.h"
41 #include "g2t_track.h"
42 #include "g2t_vertex.h"
43 
44 #include "tables/St_g2t_event_Table.h"
45 #include "tables/St_g2t_ftp_hit_Table.h"
46 #include "tables/St_g2t_svt_hit_Table.h"
47 #include "tables/St_g2t_tpc_hit_Table.h"
48 #include "tables/St_g2t_track_Table.h"
49 #include "tables/St_g2t_vertex_Table.h"
50 
51 #include "StMcEvent.hh"
52 #include "StMcTrack.hh"
53 #include "StMcTpcHit.hh"
54 #include "StMcFtpcHit.hh"
55 #include "StMcSvtHit.hh"
56 #include "StMcVertex.hh"
57 
58 #include "StParticleDefinition.hh"
59 #include "StKaonZeroShort.hh"
60 #include "StLambda.hh"
61 #include "StAntiLambda.hh"
62 #include "StPhysicalHelix.hh"
63 
64 
65 
66 #ifdef __ROOT__
67 ClassImp(StHbtMcEventReader)
68 #endif
69 
70 #ifndef ST_NO_NAMESPACES
71  using namespace units;
72 #endif
73 
74 #define __2POWER16__ 65536
75 
76 struct vertexFlag {
77  StMcVertex* vtx;
78  int primaryFlag; };
79 
80 double dedxMean(double mass, double momentum){
81  double dedxMean;
82  double tpcDedxGain = 0.174325e-06;
83  double tpcDedxOffset = -2.71889;
84  double tpcDedxRise = 776.626;
85 
86  double gamma = ::sqrt(::pow(momentum/mass,2)+1.);
87  double beta = ::sqrt(1. - 1./::pow(gamma,2));
88  double rise = tpcDedxRise*::pow(beta*gamma,2);
89  if ( beta > 0)
90  dedxMean = tpcDedxGain/::pow(beta,2) * (0.5*::log(rise)-::pow(beta,2)- tpcDedxOffset);
91  else
92  dedxMean = 1000.;
93  return dedxMean;
94 }
95 
96 //__________________
97 StHbtMcEventReader::StHbtMcEventReader(){
98  mEventCut=0;
99  mTrackCut=0;
100  mV0Cut=0;
101  mReaderStatus = 0; // "good"
102  mV0=0;
103 
104  mMotherMinvYPt = new StHbt3DHisto("Mother_MinvYPt","Mother_MinvYPt", 100,0.,5., 80,-2.,2., 80,0.,4.);
105  mMotherMinvYMt = new StHbt3DHisto("Mother_MinvYMt","Mother_MinvYMt", 100,0.,5., 80,-2.,2., 80,0.,4.);
106  mMotherMinvEtaPt = new StHbt3DHisto("Mother_MinvEtaPt","Mother_MinvEtaPt",100,0.,5., 80,-2.,2., 80,0.,4.);
107 }
108 //__________________
109 StHbtMcEventReader::~StHbtMcEventReader(){
110  if (mEventCut) delete mEventCut;
111  if (mTrackCut) delete mTrackCut;
112  if (mV0Cut) delete mV0Cut;
113 }
114 //__________________
115 StHbtString StHbtMcEventReader::Report(){
116  StHbtString temp = "\n This is the StHbtMcEventReader\n";
117  temp += "---> EventCuts in Reader: ";
118  if (mEventCut) {
119  temp += mEventCut->Report();
120  }
121  else {
122  temp += "NONE";
123  }
124  temp += "\n---> TrackCuts in Reader: ";
125  if (mTrackCut) {
126  temp += mTrackCut->Report();
127  }
128  else {
129  temp += "NONE";
130  }
131  temp += "\n---> V0Cuts in Reader: ";
132  if (mV0Cut) {
133  temp += mV0Cut->Report();
134  }
135  else {
136  temp += "NONE";
137  }
138  temp += "\n";
139  return temp;
140 }
141 //__________________
142 StHbtEvent* StHbtMcEventReader::ReturnHbtEvent(){
143  cout << " **************************************************************************************" << endl;
144  cout << " StHbtMcEventReader::ReturnHbtEvent() : Seconds elapsed since last call : " << difftime( time(0), timeStamp ) << endl;
145  cout << " **************************************************************************************" << endl;
146  timeStamp = time(0);
147 
148  // **************************************
149  // get pointer to mcEventMaker, Event *
150  // **************************************
151  StMcEvent* mcEvent = 0;
152 
153  // StMcEventMaker* mTempMaker = (StMcEventMaker*) mTheMcEventMaker;
154  // mcEvent = mTempMaker->currentMcEvent();
155  // 28sep2005 - mike lisa replaces the two lines above with the
156  // following one as per http://www.star.bnl.gov/HyperNews-star/protected/get/starsoft/
157  mcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
158 
159  if (!mcEvent){
160  cout << "StHbtMcEventReader - No StMcEvent!!! " << endl;
161  return 0;
162  }
163 
164  // ******************
165  // Event properties
166  // ******************
167  //unsigned long RunNumber = mcEvent->runNumber();
168  //cout << " MC : run: #" << RunNumber;
169  unsigned long EventNumber = mcEvent->eventNumber();
170  //cout << " * event: #" << EventNumber;
171  int Mult = mcEvent->tracks().size();
172  //cout << " * mult = " << Mult;
173  //cout << "StHbtMcEventReader::ReturnHbtEvent - We have " << Mult << " tracks " << endl;
174 
175  StHbtThreeVector VertexPosition = mcEvent->primaryVertex()->position(); // using templates, everything is fine
176  //cout << " * primary Vertex = " << VertexPosition << endl;
177 
178  if (VertexPosition.x() == VertexPosition.y() &&
179  VertexPosition.y() == VertexPosition.z() ){//x==y==z --> mark for bad events from embedding
180  cout << "StHbtMcEventReader - bad vertex !!! " << endl;
181  return 0;
182  }
183 
184 
185  StHbtEvent* hbtEvent = new StHbtEvent;
186  hbtEvent->SetEventNumber(EventNumber);
187  hbtEvent->SetCtbMult(0);
188  hbtEvent->SetZdcAdcEast(0);
189  hbtEvent->SetZdcAdcWest(0);
190  hbtEvent->SetNumberOfTpcHits(0);
191  hbtEvent->SetNumberOfTracks(Mult);
192  hbtEvent->SetNumberOfGoodTracks(Mult); // same for now
193  hbtEvent->SetReactionPlane(0.);
194  hbtEvent->SetReactionPlaneSubEventDifference(0.);
195  hbtEvent->SetPrimVertPos(VertexPosition);
196 
197  // By now, all event-wise information has been extracted and stored in hbtEvent
198  // see if it passes any front-loaded event cut
199 // if (mEventCut){
200 // cout << " performing event cut " << endl;
201 // if (!(mEventCut->Pass(hbtEvent))){ // event failed! - return null pointer (but leave Reader status flag as "good")
202 // delete hbtEvent;
203 // return 0;
204 // }
205 // }
206 
207  StHbtThreeVector p;
208 
209  for (StMcTrackIterator iter=mcEvent->tracks().begin(); iter!=mcEvent->tracks().end(); iter++){
210  StMcTrack* track = *iter;
211 
212  // y-pt for mother particles
213  if (track->particleDefinition()) {
214  if (CheckPdgIdList(mAcceptedMothers,track->particleDefinition()->pdgEncoding())) {
215  mMotherMinvYPt->Fill(track->fourMomentum().m(),track->rapidity(),track->pt());
216  mMotherMinvYMt->Fill(track->fourMomentum().m(),track->rapidity(),track->fourMomentum().mt());
217  mMotherMinvEtaPt->Fill(track->fourMomentum().m(),track->pseudoRapidity(),track->pt());
218  }
219  }
220 
221 
222  // check Pdg Id of the StMcTrack and its mc-mother and mc-daughters
223  int pdgCode = 0;
224  int motherPdgCode = 0;
225  int daughterPdgCode = 0;
226  int motherTrackId = 0;
227  if (CheckPdgIdLists()) {
228  int check=0;
229  if (!track->particleDefinition()) {
230  cout << " track has no particle definiton " << endl;
231  continue;
232  }
233  pdgCode = track->particleDefinition()->pdgEncoding();
234  //this was just stupid... malisa 14nov01 motherPdgCode = track->parent(); // 0 if no start vertex
235  if (track->parent()) {
236  motherPdgCode = track->parent()->pdgId();
237  motherTrackId = track->parent()->key();
238  }
239  if (motherPdgCode==333) cout << " phi 333" << endl;
240 
241  if ( track->stopVertex() == 0 ) {
242  check += CheckPdgIdList(pdgCode,motherPdgCode,0);
243  // cout << " no stop vertex " << endl;
244  }
245  else {
246  for (unsigned int iDaughter=0; iDaughter < track->stopVertex()->daughters().size()-1; iDaughter++) {
247  daughterPdgCode = track->stopVertex()->daughters()[iDaughter]->pdgId();
248  check += CheckPdgIdList(pdgCode,motherPdgCode,daughterPdgCode);
249  // cout << " daughterPdgCode " << daughterPdgCode;
250  }
251  }
252  // cout << " check " << check << endl;
253  if ( !(check) ) {
254  continue; // particle failed, continue with next track
255  }
256  }
257 
258  int nTpcHits = track->tpcHits().size();
259 
260  if (nTpcHits==0) {
261  //cout << "No hits in TPC-- skipping track " << endl;
262  continue;
263  }
264 
265  //cout << "nTpcHits\t" << nTpcHits << endl;
266 
267  //cout << "Getting readty to instantiate new StHbtTrack " << endl;
268 
269  StHbtTrack* hbtTrack = new StHbtTrack;
270  //cout << "StHbtTrack instantiated " << endl;
271 
272 #ifdef TheWorldIsNice
273  hbtTrack->SetP( track->momentum() ); // set momentum
274 #else
275  StHbtThreeVector tmpP;
276  tmpP.setX( track->momentum().x() );
277  tmpP.setY( track->momentum().y() );
278  tmpP.setZ( track->momentum().z() );
279  hbtTrack->SetP( tmpP ); // set momentum
280 #endif
281  //cout << " P " << hbtTrack->P() << endl;
282 
283  hbtTrack->SetTrackId(track->key()+motherTrackId*__2POWER16__);
284 
285  hbtTrack->SetNHits( nTpcHits ); // hits in Tpc
286  hbtTrack->SetNHitsPossible(nTpcHits ); // hits in Tpc
287  //cout << " NHits " << hbtTrack->NHits() << endl;
288 
289  // **************************************************************************************
290  // set NSigma's, derivation in dedx distribustion from being a poin, kaon, proton [sigma]
291  // **************************************************************************************
292  int geantPid = track->particleDefinition()->pdgEncoding();
293  //cout << geantPid << " " << track->particleDefinition()->mass() << " " << track->particleDefinition()->charge() << endl;
294 
295  switch (geantPid) {
296  case 11: // intentional fall-through
297  case -11: // gid=211,-211 is pion
298  hbtTrack->SetNSigmaElectron(0.);
299  hbtTrack->SetNSigmaPion(-999.);
300  hbtTrack->SetNSigmaKaon(-999.);
301  hbtTrack->SetNSigmaProton(-999.);
302  break;
303  case 211: // intentional fall-through
304  case -211: // gid=211,-211 is pion
305  hbtTrack->SetNSigmaElectron(999.);
306  hbtTrack->SetNSigmaPion(0.);
307  hbtTrack->SetNSigmaKaon(-999.);
308  hbtTrack->SetNSigmaProton(-999.);
309  break;
310  case 321: // intentional fall-through
311  case -321: // gid=321,-321 is kaon
312  hbtTrack->SetNSigmaElectron(999.);
313  hbtTrack->SetNSigmaPion(999.0);
314  hbtTrack->SetNSigmaKaon(0.);
315  hbtTrack->SetNSigmaProton(-999.);
316  break;
317  case 2212: // intentional fall-through
318  case -2212: // gid=2212,-2212 is proton
319  hbtTrack->SetNSigmaElectron(999.);
320  hbtTrack->SetNSigmaPion(999.);
321  hbtTrack->SetNSigmaKaon(999.);
322  hbtTrack->SetNSigmaProton(0.);
323  break;
324  default:
325  hbtTrack->SetNSigmaElectron(999.);
326  hbtTrack->SetNSigmaPion(999.);
327  hbtTrack->SetNSigmaKaon(999.);
328  hbtTrack->SetNSigmaProton(999.);
329  break;
330  }
331 
332  //cout << "Nsig electron,pion,kaon,proton : " << hbtTrack->NSigmaElectron() << " ";
333  //cout << hbtTrack->NSigmaPion() << " " hbtTrack->NSigmaKaon() << " " << hbtTrack->NSigmaProton() << " PDG code " << geantPId << endl;
334 
335  hbtTrack->SetdEdx( dedxMean( track->particleDefinition()->mass(), track->momentum().mag() ) ); // not in mike's
336  //cout << " dedx " << hbtTrack->Dedx() << endl;
337 
338  hbtTrack->SetPt( track->momentum().perp() );
339  //cout << " Pt " << hbtTrack->Pt() << endl;
340 
341  hbtTrack->SetCharge( (int)(track->particleDefinition()->charge()) );
342  //cout << " choarge " << hbtTrack->Charge() << endl;
343 
344 #ifdef TheWorldIsNice
345  StPhysicalHelix helix = StPhysicalHelix( hbtTrack->P(), track->startVertex()->position(), HBT_BFIELD, hbtTrack->Charge() );
346 #else
347  StHbtThreeVector tmpSV;
348  tmpSV.setX( track->startVertex()->position().x() );
349  tmpSV.setY( track->startVertex()->position().y() );
350  tmpSV.setZ( track->startVertex()->position().z() );
351  StPhysicalHelixD helix = StPhysicalHelixD( hbtTrack->P(), tmpSV, HBT_BFIELD, hbtTrack->Charge() );
352 #endif
353 
354  hbtTrack->SetHelix(helix);
355 
356 
357  double pathlength = helix.pathLength(VertexPosition);
358  //cout << "pathlength\t" << pathlength << endl;
359  StHbtThreeVector DCAxyz = helix.at(pathlength)-VertexPosition;
360  //cout << "DCA\t\t" << DCAxyz << " " << DCAxyz.perp() << endl;
361 
362  hbtTrack->SetDCAxy( DCAxyz.perp() ); // in xy-plane
363  hbtTrack->SetDCAz( DCAxyz.z() ); // in z direction
364 
365  hbtTrack->SetChiSquaredXY( 0.);
366  hbtTrack->SetChiSquaredZ( 0.);
367 
368 
369  //cout << "pushing..." <<endl;
370  // By now, all track-wise information has been extracted and stored in hbtTrack
371  // see if it passes any front-loaded event cut
372  if (mTrackCut){
373  if (!(mTrackCut->Pass(hbtTrack))){ // track failed - delete it and skip the push_back
374  delete hbtTrack;
375  continue;
376  }
377  }
378  hbtEvent->TrackCollection()->push_back(hbtTrack);
379  }
380  cout << hbtEvent->TrackCollection()->size() << " tracks pushed to collection" << endl;
381 
382  // ******************
383  // fill v0 collection
384  // ******************
385  StKaonZeroShort* k0Short = StKaonZeroShort::instance();
386  StLambda* lambda = StLambda::instance();
387  StAntiLambda* antiLambda = StAntiLambda::instance();
388  for (StMcVertexIterator vIter=mcEvent->vertices().begin(); vIter!=mcEvent->vertices().end(); vIter++){
389  StMcVertex* vertex = *vIter;
390  StMcTrack* parent = (StMcTrack*)vertex->parent(); // get parent
391  if (parent) {
392  if ( parent->particleDefinition() == k0Short ||
393  parent->particleDefinition() == lambda ||
394  parent->particleDefinition() == antiLambda &&
395  vertex->numberOfDaughters() == 2) {
396 #ifdef STHBTDEBUG
397  cout << " v0 Id : " << parent->particleDefinition()->name() << endl;
398 #endif
399  StMcTrack* pos;
400  StMcTrack* neg;
401  if ( (*(vertex->daughters().begin()))->particleDefinition()->charge() > 0 ) {
402  pos = *(vertex->daughters().begin()); // positive daughter
403  neg = *(vertex->daughters().end()-1); // negative daughter
404  }
405  else if ( (*(vertex->daughters().begin()))->particleDefinition()->charge() < 0 ) {
406  neg = *(vertex->daughters().begin()); // negative daughter
407  pos = *(vertex->daughters().end()-1); // positive daughter
408  }
409  else continue;
410 
411  // fill the StHbtV0 structure
412  StHbtV0* hbtv0 = new StHbtV0();
413  hbtv0->SetdecayLengthV0( abs(vertex->position()-VertexPosition) );
414  hbtv0->SetdecayVertexV0( vertex->position() );
415 
416  StPhysicalHelixD posHelix = StPhysicalHelixD( pos->momentum(), vertex->position(),
417  HBT_BFIELD, pos->particleDefinition()->charge() );
418  StPhysicalHelixD negHelix = StPhysicalHelixD( neg->momentum(), vertex->position(),
419  HBT_BFIELD, neg->particleDefinition()->charge() );
420  StPhysicalHelixD v0Helix = StPhysicalHelixD( pos->momentum()+neg->momentum(), vertex->position(),
421  HBT_BFIELD, 0 );
422 
423  double posPathLength = posHelix.pathLength( vertex->position() );
424  double negPathLength = negHelix.pathLength( vertex->position() );
425 
426  hbtv0->SetdcaV0Daughters( abs(posHelix.at(posPathLength)-negHelix.at(negPathLength)) );
427  //cout << " dcaV0Daughters " << hbtv0->dcaV0Daughters() << endl;
428 
429  hbtv0->SetdcaV0ToPrimVertex( v0Helix.distance( VertexPosition ) );
430  //cout << " dcaV0ToPrimVertex " << hbtv0->dcaV0ToPrimVertex() << endl;
431 
432  hbtv0->SetdcaPosToPrimVertex( posHelix.distance( VertexPosition ) ); // VertexPosition = prim vert pos
433  hbtv0->SetdcaNegToPrimVertex( negHelix.distance( VertexPosition ) ); // VertexPosition = prim vert pos
434 
435  hbtv0->SetmomPos( pos->momentum() );
436  hbtv0->SetmomNeg( neg->momentum() );
437 
438  hbtv0->SettpcHitsPos( pos->tpcHits().size() );
439  hbtv0->SettpcHitsNeg( neg->tpcHits().size() );
440 
441  //hbtv0->SetTrackTopologyMapPos(0,pos->topologyMap().data(0));
442  //hbtv0->SetTrackTopologyMapPos(1,pos->topologyMap().data(1));
443  //hbtv0->SetTrackTopologyMapNeg(0,neg->topologyMap().data(0));
444  //hbtv0->SetTrackTopologyMapNeg(1,neg->topologyMap().data(1));
445 
446  hbtv0->SetkeyPos(pos->key());
447  hbtv0->SetkeyNeg(neg->key());
448 
449  hbtv0->UpdateV0();
450 
451  // apply v0 cut
452  if (mV0Cut){
453  if (!(mV0Cut->Pass(hbtv0))){ // track failed - delete it and skip the push_back
454  delete hbtv0;
455  continue;
456  }
457  }
458 
459  hbtEvent->V0Collection()->push_back(hbtv0);
460  }
461  }
462  }
463  cout << hbtEvent->V0Collection()->size() << " v0s pushed to collection" << endl;
464 
465  return hbtEvent;
466 }
467 
468 
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StHelix::distance
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Generated by   doxygen 1.8.5