StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructEStructReader.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructEStructReader.cxx,v 1.1 2013/05/04 23:45:43 prindle Exp $
4  *
5  * Author: Duncan Prindle
6  *
7  **********************************************************************
8  *
9  * Description: event reader class for reading EStruct event format files
10  * Essentially duplicated from StEstructMuDstReader.
11  *
12  ***********************************************************************/
13 #include "StEStructEStructReader.h"
14 
15 #include "StEStructEventCuts.h"
16 #include "StEStructTrackCuts.h"
17 #include "StEStructPool/EventMaker/StEStructEvent.h"
18 #include "StEStructPool/EventMaker/StEStructTrack.h"
19 #include "StEStructPool/EventMaker/StEStructCentrality.h"
20 
21 #include "StDetectorId.h"
22 #include "TFile.h"
23 #include "TList.h"
24 
25 
26 StEStructEStructReader::StEStructEStructReader() {
27  mhasdEdxCuts = 0;
28  mhasToFCuts = 0;
29  mhasPrimaryCuts = 0;
30  mhasVertexRadiusCuts = 0;
31  mFile = 0;
32  mkeyList = 0;
33  miEvent = 0;
34  mnEvents = 0;
35  mECuts = 0;
36  mTCuts = 0;
37  mInChain = false;
38  mAmDone = false;
39  mUseGlobalTracks=false;
40  mPileup = new Pileup();
41 };
42 StEStructEStructReader::StEStructEStructReader(char* estructFormatFile,
43  StEStructEventCuts* ecuts,
44  StEStructTrackCuts* tcuts) {
45  mhasdEdxCuts = 0;
46  mhasToFCuts = 0;
47  mhasPrimaryCuts = 0;
48  mhasVertexRadiusCuts = 0;
49  mFile = new TFile(estructFormatFile);
50  miEvent = 0;
51  mkeyList = mFile->GetListOfKeys();
52  mnEvents = mkeyList->GetEntries();
53  setEventCuts(ecuts);
54  setTrackCuts(tcuts);
55  mInChain = false;
56  mAmDone = false;
57  mUseGlobalTracks=false;
58  mPileup = new Pileup();
59 };
60 
61 //-------------------------------------------------------------------------
62 StEStructEStructReader::~StEStructEStructReader() {
63  delete mPileup;
64  delete mFile;
65 };
66 
67 void StEStructEStructReader::setEventCuts(StEStructEventCuts* ecuts) {
68  mECuts=ecuts;
69  if (!mECuts->doFillHists()) {
70  return;
71  }
72  if ( mECuts->hasToFFractionCut() ) {
73  mhasToFCuts = 1;
74  ToFBefore = new TH2F("dEdxToFNoCut2D","dEdxToFNoCut2D",50,1,1000,50,1,750);
75  ToFAfter = new TH2F("dEdxToFCut2D","dEdxToFCut2D",50,1,1000,50,1,750);
76  mECuts->addCutHists(ToFBefore, ToFAfter, "ToFPlots2D");
77  }
78  if ( mECuts->hasPrimaryFractionCut() ) {
79  mhasPrimaryCuts = 1;
80  PrimaryBefore = new TH2F("PrimVGlobNoCut2D","PrimVGlobNoCut2D",120,1,1200,200,1,2000);
81  PrimaryAfter = new TH2F("PrimVGlobCut2D","PrimVGlobCut2D",120,1,1200,200,1,2000);
82  mECuts->addCutHists(PrimaryBefore, PrimaryAfter, "PrimaryPlots2D");
83  }
84  if ( mECuts->hasPrimaryVertexRadiusCut() ) {
85  mhasVertexRadiusCuts = 1;
86  VRadiusBefore = new TH2F("VertexRadiusNoCut2D","VertexRadiusNoCut2D",50,-10,10,50,-10,10);
87  VRadiusAfter = new TH2F("VertexRadiusCut2D","VertexRadiusCut2D",50,-10,10,50,-10,10);
88  mECuts->addCutHists(VRadiusBefore, VRadiusAfter, "VertexRadiusPlots2D");
89  }
90 };
91 void StEStructEStructReader::setTrackCuts(StEStructTrackCuts* tcuts) {
92  mTCuts=tcuts;
93  if (!mTCuts->doFillHists()) {
94  return;
95  }
96  if ( mTCuts->hasElectronCut() ) {
97  mhasdEdxCuts = 1;
98  dEdxBefore = new TH2F("dEdxNoCut2D","dEdxNoCut2D",150,-1.5,1.5,150,0,0.000015);
99  dEdxAfter = new TH2F("dEdxCut2D","dEdxCut2D",150,-1.5,1.5,150,0,0.000015);
100  mTCuts->addCutHists(dEdxBefore, dEdxAfter, "dEdxPlots2D");
101  }
102  dEdxBetaBefore = new TH3F("dEdxBetaNoCut3D","dEdxBetaNoCut3D",100,-2.0,2.0,100,0.0000015,0.000006,100,0.55,1.05);
103  dEdxBetaAfter = new TH3F("dEdxBetaCut3D","dEdxBetaCut3D",100,-2.0,2.0,100,0.0000015,0.000006,100,0.55,1.05);
104  mTCuts->addCutHists(dEdxBetaBefore, dEdxBetaAfter, "dEdxBetaPlots3D");
105 };
106 
107 //-------------------------------------------------------------------------
108 StEStructEvent* StEStructEStructReader::next() {
109 
110  StEStructEvent* retVal=NULL;
111  if(!mkeyList) return retVal;
112 
113  mEventIn = (StEStructEvent *) mFile->Get(mkeyList->At(miEvent)->GetName());
114  miEvent++;
115  mEventIn->FillChargeCollections();
116 
117  return fillEvent();
118 }
119 StEStructEvent* StEStructEStructReader::fillEvent(){
120 
121  StEStructEvent* retVal = new StEStructEvent();
122  mNumGoodTracks = 0;
123 
124  float x;
125  float y;
126  float z;
127  bool useEvent = true;
128 
129  // Note: Recommended primary vertex cuts are made in StEStructEventCuts.h
130  x = mEventIn->Vx();
131  y = mEventIn->Vy();
132  z = mEventIn->Vz();
133 
134  // In StEStructMuDstReader we invoke mECuts->goodTrigger(muDst).
135  // Comments suggest this is a topological cut. Skip it for now.
136  // >>>>>>Check: If it is really required we need to figure out how to invoke it. I think is has been superseded.
137  int good = 0;
138  if ((fabs(x) < 1e-5) && (fabs(y) < 1e-5) && (fabs(z) < 1e-5)) {
139  useEvent = false;
140  } else if ( !mECuts->goodPrimaryVertexZ(z) ) {
141  useEvent = false;
142  } else {
143  good++;
144  }
145  if (!mECuts->goodPrimaryVertexRadius(sqrt(x*x+y*y))) {
146  useEvent = false;
147  }
148 
149  // In StEStructMuDstReader we cut on distance between VPD Z vertex and
150  // vertex determined by tracking. Didn't store VPD information in EStruct format.
151  // >>>>>>Check: If we really need to tighten this cut we need to store VPD in EStruct.
152 
153  int ndEdx = 0;
154  int nToF = 0;
155  int nTracks = countGoodTracks(&ndEdx, &nToF);
156  if (!mECuts->goodToFFraction(ndEdx, nToF)) {
157  useEvent = false;
158  }
159 
160  // >>>>>>Check: Are you going to redefine primary and global tracks?
161  // >>>>>>Check: If so do that before this cut (probably in countGoodTracks) and implement cut here.
162  //int nGlobal = muDst->globalTracks()->GetEntries();
163  //int nPrimary = muDst->primaryTracks()->GetEntries();
164  //if (!mECuts->goodPrimaryFraction(nPrimary, nGlobal)) {
165  // useEvent = false;
166  //}
167 
168  retVal->SetCentrality((double)nTracks);
169  if(!mECuts->goodCentrality(retVal->Centrality())) useEvent=false;
170 
171  retVal->SetRefMult( mEventIn->RefMult() );
172  retVal->SetctbMult( mEventIn->ctbMult() );
173  retVal->SetNumPrim( mEventIn->NumPrim() );
174 
175  if(useEvent){
176 
177  retVal->SetEventID(mEventIn->EventID());
178  retVal->SetRunID(mEventIn->RunID());
179  retVal->SetVertex(x,y,z);
180  retVal->SetZDCe((float)mEventIn->ZDCe());
181  retVal->SetZDCw((float)mEventIn->ZDCw());
182  // I think this zdc coincidence rate is one number for the event
183  // Don't know how to get rate at time of event.
184  // Actually, in old runs this was true. Now it should be updated every
185  // 10 or 15 seconds throughout the run.
186  retVal->SetZDCCoincidence((float)mEventIn->ZDCCoincidence());
187  retVal->SetBField((float)mEventIn->BField());
188 
189  if (!fillTracks(retVal)) useEvent=false;
190 
191  }
192 
193  mECuts->fillHistogram(mECuts->primaryVertexZName(),z,useEvent);
194  // >>>>>>Check: Don't think we have VPD information in StEStructEvent.
195  //if (btofHeader) {
196  // mECuts->fillHistogram(mECuts->primaryVPDVertexName(),z-vpdZ,useEvent);
197  //}
198  mECuts->fillHistogram(mECuts->centralityName(),(float)nTracks,useEvent);
199  mECuts->fillHistogram(mECuts->goodtoffractionName(),(float)ndEdx,(float)nToF,useEvent);
200 
201  mECuts->fillHistogram("ToFPlots",(double)ndEdx,(double)nToF,useEvent);
202  //mECuts->fillHistogram("PrimaryPlots",(double)nPrimary,(double)nGlobal,useEvent);
203  mECuts->fillHistogram("VertexRadiusPlots",x,y,useEvent);
204 
205  // >>>>>>Check: My pileup finder relies on global tracks. If you really want to tighten up pileup cuts
206  // >>>>>>Check: here we should probably store some information from the pileup code in the StEStructEvent.
207  // Two pileup cuts using my vertex finding.
208  // First we look for vertices which have tracks from only one side of TPC.
209  // Then we look for matching pileup vertices where one of them reconstructs near the primary vertex.
210  //if (useEvent && mECuts->hasZVertMatchCut()) {
211  // double zDist = 0; // Distance from z to nearest pileup half vertex.
212  // int nPile = mPileup->find(muDst);
213  // if ((nPile > 0) && (6 == mPileup->flag(0))) {
214  // double zd1 = z - mPileup->z(0,2);
215  // double zd2 = z - mPileup->z(0,3);
216  // zDist = (fabs(zd1) < fabs(zd2)) ? zd1 : zd2;
217  // if (!mECuts->goodZVertMatch(zDist)) {
218  // useEvent = false;
219  // }
220  // mECuts->fillHistogram(mECuts->zVertMatchName(),zDist,useEvent);
221  // }
222  //}
223 
224  //if (useEvent && mECuts->hasZVertSepCut()) {
225  // double zDist; // Distance from z to nearest pileup vertex.
226  // int mPile; // Multiplicity of nearest pileup vertex.
227  // int nPile; // Total number of pileup vertices found (multiple of 2).
228  // nPile = mPileup->nearest(muDst,z,&zDist,&mPile);
229  // if (nPile > 0 && !mECuts->goodZVertSep(zDist)) {
230  // useEvent = false;
231  // }
232  // mECuts->fillHistogram(mECuts->zVertSepName(),zDist,useEvent);
233  //}
234 
235  if (!useEvent) {
236  delete retVal;
237  retVal=NULL;
238  }
239  if (retVal) retVal->FillChargeCollections(); //creates track list by charge
240 
241  return retVal;
242 }
243 
244 
245 //-------------------------------------------------------------
246 bool StEStructEStructReader::fillTracks(StEStructEvent* estructEvent) {
247 
248  //
249  // create a single EStruct track, check cuts,
250  // fill variables, add to StEStructEvent tracks which
251  // does a copy.
252  //
253 
254  int numTracks = mEventIn->Ntrack();
255  if (0==numTracks) {
256  return false;
257  }
258 
259  for (int ich=0;ich<2;ich++) {
260  StEStructTrackCollection *tc;
261  if (ich==0) {
262  tc=mEventIn->TrackCollectionP();
263  } else {
264  tc=mEventIn->TrackCollectionM();
265  }
266  for (StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
267  StEStructTrack* track = (*iter);
268 
269  bool useTrack=true;
270  float ptot = track->Charge() * track->Ptot();
271  float beta = track->beta();
272 
273  useTrack = isTrackGoodToUse( track ); //this includes track kine cuts
274  mTCuts->fillHistograms(useTrack);
275 
276  mTCuts->fillHistogram("dEdxPlots",ptot,track->Dedx(),useTrack);
277  if (beta != 0) {
278  mTCuts->fillHistogram("dEdxBetaPlots",ptot,track->Dedx(),beta,useTrack);
279  }
280 
281  if (useTrack) {
282  estructEvent->AddTrack(track);
283  }
284  }
285  }
286  return true;
287 }
288 //-------------------------------------------------------------
289 // This method checks all track cuts.
290 // No histogramming or copying data around.
291 bool StEStructEStructReader::isTrackGood(StEStructTrack* track) {
292 
293  bool useTrack=true;
294 
295  // When using global tracks use the outerHelix to calculate eta, phi and dca.
296  // >>>>>> Not sure it makes sense to allow global track cuts for EStruct format events.
297  // >>>>>> Need to look into this if that is required.
298  Float_t mDCA;
299  mEta = track->Eta();
300  mPhi = track->Phi();
301  mDCA = track->DcaGlobal();
302 
303  // Do eta cut first so my ThisCut can use the eta value.
304  useTrack = (mTCuts->goodEta(mEta) && useTrack);
305  useTrack = (mTCuts->goodFlag(track->Flag()) && useTrack);
306  useTrack = (mTCuts->goodCharge(track->Charge()) && useTrack);
307  useTrack = (mTCuts->goodNFitPoints(track->NFitPoints()) && useTrack);
308  useTrack = (mTCuts->goodNFitNMax((float)((float)track->NFitPoints()/(float)track->NMaxPoints())) && useTrack);
309  useTrack = (mTCuts->goodGlobalDCA(mDCA) && useTrack);
310  useTrack = (mTCuts->goodChi2(track->Chi2()) && useTrack);
311  useTrack = (mTCuts->goodPhi(mPhi) && useTrack);
312  if(track->Pt() < 0.15) useTrack = false; // basic pt cut, ranges checked in isTrackGoodToUse
313 
314  //--> But add a quick electron removal... for selected p ranges
315  // Note I only want to do this if I am defining an electron dEdx cut.
316  if (mhasdEdxCuts && mTCuts->goodElectron( track->PIDe_dEdx() ) && useTrack ) {
317  float p = track->Ptot();
318  if( (p>0.2 && p<0.45) || (p>0.7 && p<0.8) ) {
319  useTrack=false;
320  }
321  }
322  //--> end of electron pid
323 
324  return useTrack;
325 };
326 
327 //----includes call to isTrackGood ... and adds some more ....
328 bool StEStructEStructReader::isTrackGoodToUse(StEStructTrack* track){
329 
330  bool useTrack=isTrackGood(track);
331  if(useTrack){
332 
333  float pt = track->Pt();
334  useTrack = (mTCuts->goodPt(pt) && useTrack);
335  float yt = track->Yt();
336  useTrack = (mTCuts->goodYt(yt) && useTrack);
337  float xt = track->Xt();
338  useTrack = (mTCuts->goodXt(xt) && useTrack);
339 
340  }
341  return useTrack;
342 };
343 
344 //-------------------------------------------------------------------------
345 //-------------------------------------------------------------
346 // This method counts all good track.
347 // No histogramming or copying data around.
348 int StEStructEStructReader::countGoodTracks(int *ndEdx, int *nToF) {
349  mNumGoodTracks = 0;
350  int ndedx = 0;
351  int ntof = 0;
352 
353  int numTracks = mEventIn->Ntrack();
354  // >>>>>> If we need to redefine global and primary tracks we need to do more work.
355  if (0==numTracks) {
356  return 0;
357  }
358  // Want to compare pp analysis with |\eta| < 0.5 and |\eta| < 1.0
359  // countGoodTracks is used for centrality, so we change \eta cuts here
360  // (then change back at and)
361  // >>>>>Expect this code to normally be commented out.<<<<<
362 // double etaMin = mTCuts->meta[0];
363 // double etaMax = mTCuts->meta[1];
364 // mTCuts->meta[0] = -1;
365 // mTCuts->meta[1] = +1;
366 
367  for (int ich=0;ich<2;ich++) {
368  StEStructTrackCollection *tc;
369  if (ich==0) {
370  tc=mEventIn->TrackCollectionP();
371  } else {
372  tc=mEventIn->TrackCollectionM();
373  }
374  for (StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
375  StEStructTrack* track = (*iter);
376  if (isTrackGood(track)) {
377  mNumGoodTracks++;
378  if ((fabs(track->PIDpi_dEdx()) < 2) ||
379  (fabs(track->PIDk_dEdx()) < 2) ||
380  (fabs(track->PIDp_dEdx()) < 2)) {
381  ndedx++;
382  }
383  if ((fabs(track->PIDpi_ToF()) < 2) ||
384  (fabs(track->PIDpi_ToF()) < 2) ||
385  (fabs(track->PIDpi_ToF()) < 2)) {
386  ntof++;
387  }
388  }
389  }
390  }
391  *ndEdx = ndedx;
392  *nToF = ntof;
393 // mTCuts->meta[0] = etaMin;
394 // mTCuts->meta[1] = etaMax;
395  return mNumGoodTracks;
396 }
397 
398 //-------------------------------------------------------------------------
399 void StEStructEStructReader::fillEStructTrack(StEStructTrack* eTrack,StEStructTrack* mTrack){
400 
401  eTrack->SetPx(mTrack->Px());
402  eTrack->SetPy(mTrack->Py());
403  eTrack->SetPz(mTrack->Pz());
404 
405  eTrack->SetBx(mTrack->Bx());
406  eTrack->SetBy(mTrack->By());
407  eTrack->SetBz(mTrack->Bz());
408 
409  if (0 == mTrack->BxGlobal() && 0 == mTrack->ByGlobal() && 0 == mTrack->BzGlobal()) {
410  printf("Found track with dca = (0,0,0)\n");
411  }
412  eTrack->SetBxGlobal(mTrack->BxGlobal());
413  eTrack->SetByGlobal(mTrack->ByGlobal());
414  eTrack->SetBzGlobal(mTrack->BzGlobal());
415 
416  eTrack->SetEta(mEta);
417  eTrack->SetPhi(mPhi);
418  eTrack->SetDedx(mTrack->Dedx());
419  eTrack->SetChi2(mTrack->Chi2());
420  eTrack->SetHelix(mTrack->Helix());
421 
422  //
423  // -> note in my analysis I chose nSigma instead of prob.
424  //
425  eTrack->SetPIDe_dEdx(mTrack->PIDe_dEdx());
426  eTrack->SetPIDpi_dEdx(mTrack->PIDpi_dEdx());
427  eTrack->SetPIDp_dEdx(mTrack->PIDp_dEdx());
428  eTrack->SetPIDk_dEdx(mTrack->PIDk_dEdx());
429  eTrack->SetPIDd_dEdx(10000.);
430 
431 
432  eTrack->SetPIDe_ToF(mTrack->PIDe_ToF());
433  eTrack->SetPIDpi_ToF(mTrack->PIDpi_ToF());
434  eTrack->SetPIDp_ToF(mTrack->PIDp_ToF());
435  eTrack->SetPIDk_ToF(mTrack->PIDk_ToF());
436  eTrack->SetPIDd_ToF(999);
437  eTrack->SetBeta(mTrack->beta());
438 
439  eTrack->SetNFitPoints(mTrack->NFitPoints());
440  eTrack->SetNFoundPoints(mTrack->NFoundPoints());
441  eTrack->SetNMaxPoints(mTrack->NMaxPoints());
442 
443  eTrack->SetTopologyMapTPCNHits(mTrack->TopologyMapTPCNHits());
444  eTrack->SetTopologyMapData(0,mTrack->TopologyMapData(0));
445  eTrack->SetTopologyMapData(1,mTrack->TopologyMapData(1));
446 
447  eTrack->SetDetectorID(1); //TPC
448  eTrack->SetFlag(mTrack->Flag());
449  eTrack->SetCharge(mTrack->Charge());
450 };
451 
452 
453 /***********************************************************************
454  *
455  * $Log: StEStructEStructReader.cxx,v $
456  * Revision 1.1 2013/05/04 23:45:43 prindle
457  * Code to read EStruct format files. Modified from StEStructMuDstReader.
458  *
459  *
460  * Revision 1.0 2013/05/02 18:20:32 prindle
461  * Copy of StEStructMuDstReader specialized for reading EStruct format files..
462  *
463  *********************************************************************/
Pileup
Definition: Pileup.h:79
Generated by   doxygen 1.8.5