StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StTofpMcAnalysisMaker.cxx
1 /*************************************************
2  *
3  * $Id: StTofpMcAnalysisMaker.cxx,v 1.6 2018/02/26 23:26:51 smirnovd Exp $
4  *
5  *************************************************/
6 
7 #include <Stiostream.h>
8 #include "TH1.h"
9 #include "TH2.h"
10 #include "StThreeVectorF.hh"
11 #include "StThreeVectorD.hh"
12 #include "StEventTypes.h"
13 #include "StMcEventTypes.hh"
14 #include "StAssociationMaker/StAssociationMaker.h"
15 #include "StAssociationMaker/StTrackPairInfo.hh"
16 #include "StEventMaker/StEventMaker.h"
17 #include "StTofUtil/StTofGeometry.h"
18 #include "StTofUtil/tofPathLength.hh"
19 #include "StMessMgr.h"
20 #include "StTofpMcAnalysisMaker.h"
21 
22 
23 
24 const Int_t StTofpMcAnalysisMaker::mPtBin = 50;
25 const Int_t StTofpMcAnalysisMaker::mYBin = 14;
26 const Float_t StTofpMcAnalysisMaker::mPtMin = 0.;
27 const Float_t StTofpMcAnalysisMaker::mPtMax = 5.;
28 const Float_t StTofpMcAnalysisMaker::mYMin = -1.2;
29 const Float_t StTofpMcAnalysisMaker::mYMax = 0.2;
30 
31 
32 
33 //----------------------------------------------------------------
34 const StGlobalTrack* partnerTrack(mcTrackMapType* map, StMcTrack* mT) {
35  mcTrackMapIter i = map->find(mT);
36  const StGlobalTrack* rT = 0;
37  if (i != map->end()) rT = (*i).second->partnerTrack();
38  return rT;
39 }
40 
41 
42 
43 //_________________________________________________
44 StTofpMcAnalysisMaker::StTofpMcAnalysisMaker(const char *name, const char *title):StMaker(name,title)
45 {
46  // StTofpMcAnalysisMaker Constructor
47 
48  mOuterTrackGeometry = true;
49  mMinHitsPerTrack = 0;
50 
51 
52  // - zero all pointers defined in the header file
53  hMcPionPlus = 0;
54  hMcPionMin = 0;
55  hMcKaonPlus = 0;
56  hMcKaonMin = 0;
57  hMcProton = 0;
58  hMcAntiProton = 0;
59  hMcElectron = 0;
60  hRcPionPlus = 0;
61  hRcPionMin = 0;
62  hRcKaonPlus = 0;
63  hRcKaonMin = 0;
64  hRcProton = 0;
65  hRcAntiProton = 0;
66  hRcElectron = 0;
67  hMatchPionPlus = 0;
68  hMatchPionMin = 0;
69  hMatchKaonPlus = 0;
70  hMatchKaonMin = 0;
71  hMatchProton = 0;
72  hMatchAntiProton = 0;
73  hMatchElectron = 0;
74  //-
75  hMomResPtPion = 0;
76  hMultRef = 0;
77  //-
78  hTofpHitMap1 = 0;
79  hTofpHitMap2 = 0;
80  hTofpHitMap3 = 0;
81  hTofpHitMap4 = 0;
82  hTofpSlatIdA0 = 0;
83  hTofpSlatIdA1 = 0;
84  hTofpSlatIdB1 = 0;
85  hTofpSlatIdD1 = 0;
86  hTofpSlatIdD2 = 0;
87  hTofpSlatIdE1 = 0;
88  hTofpSlatIdE2 = 0;
89  hTofpSlatIdE3 = 0;
90  hTofpSlatIdE4 = 0;
91  hTofpSlatIdE5 = 0;
92  hTofpSlatIdF1 = 0;
93  hTofpSlatHitVecSize = 0;
94 }
95 
96 //_________________________________________________
97 StTofpMcAnalysisMaker::~StTofpMcAnalysisMaker()
98 {
99  // StTofpMcAnalysisMaker Destructor
100  if (Debug()) LOG_INFO << "Inside StTofpMcAnalysisMaker Destructor" << endm;
101 }
102 
103 //_____________________________________________________________________________
104 
105 void StTofpMcAnalysisMaker::Clear(const char*)
106 {
107  StMaker::Clear();
108 }
109 
110 //_________________________________________________
112 {
113  return StMaker::Finish();
114 }
115 
116 
117 //_________________________________________________
118 Int_t StTofpMcAnalysisMaker::Init()
119 {
120  // StTofpMcAnalysisMaker - Init
121 
122  hMomResPtPion = new TH2F("momResMomTPion","(|pmc| - |p|)/|pmc| as a Function of Pt Pion",mPtBin,mPtMin,mPtMax,200,-1.,1.);
123  hMomResPtPion->SetXTitle("Pt (GeV/c)");
124  hMomResPtPion->SetYTitle("(|pmc| - |p|)/|pmc| (GeV/c)");
125 
126  hMultRef = new TH1F("multRef","Mult Ref",100,0.,1000.);
127  hMultRef->SetXTitle("uncorrected negative primary tracks");
128 
129  //-
130  hMcElectron = new TH2F("mcElectron","Mc Electron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
131  hMcElectron->SetXTitle("p_{T} (GeV/c)");
132  hMcElectron->SetYTitle("pseudorapidity");
133  hMcElectron->Sumw2();
134  hRcElectron = new TH2F("rcElectron","RcElectron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
135  hRcElectron->SetXTitle("p_{T} (GeV/c)");
136  hRcElectron->SetYTitle("pseudorapidity");
137  hRcElectron->Sumw2();
138  hMatchElectron = new TH2F("mElectron","TOF matched electron",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
139  hMatchElectron->Sumw2();
140  hMatchElectron->SetXTitle("p_{T} (GeV/c)");
141  hMatchElectron->SetYTitle("pseudorapidity");
142 
143  //-
144  hMcPionPlus = new TH2F("mcPionPlus","Mc Pion Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
145  hMcPionPlus->SetXTitle("p_{T} (GeV/c)");
146  hMcPionPlus->SetYTitle("pseudorapidity");
147  hMcPionPlus->Sumw2();
148  hRcPionPlus = new TH2F("rcPionPlus","Rc Pion Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
149  hRcPionPlus->SetXTitle("p_{T} (GeV/c)");
150  hRcPionPlus->SetYTitle("pseudorapidity");
151  hRcPionPlus->Sumw2();
152  hMatchPionPlus = new TH2F("mPionPlus","TOF matched pi+",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
153  hMatchPionPlus->Sumw2();
154  hMatchPionPlus->SetXTitle("p_{T} (GeV/c)");
155  hMatchPionPlus->SetYTitle("pseudorapidity");
156 
157  //-
158  hMcPionMin = new TH2F("mcPionMin","Mc Pion Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
159  hMcPionMin->SetXTitle("p_{T} (GeV/c)");
160  hMcPionMin->SetYTitle("pseudorapidity");
161  hMcPionMin->Sumw2();
162  hRcPionMin = new TH2F("rcPionMin","Rc Pion Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
163  hRcPionMin->SetXTitle("p_{T} (GeV/c)");
164  hRcPionMin->SetYTitle("pseudorapidity");
165  hRcPionMin->Sumw2();
166  hMatchPionMin = new TH2F("mPionMin","TOF matched pi-",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
167  hMatchPionMin->Sumw2();
168  hMatchPionMin->SetXTitle("p_{T} (GeV/c)");
169  hMatchPionMin->SetYTitle("pseudorapidity");
170 
171  //-
172  hMcKaonPlus = new TH2F("mcKaonPlus","Mc Kaon Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
173  hMcKaonPlus->SetXTitle("p_{T} (GeV/c)");
174  hMcKaonPlus->SetYTitle("pseudorapidity");
175  hMcKaonPlus->Sumw2();
176  hRcKaonPlus = new TH2F("rcKaonPlus","Rc Kaon Plus",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
177  hRcKaonPlus->SetXTitle("p_{T} (GeV/c)");
178  hRcKaonPlus->SetYTitle("pseudorapidity");
179  hRcKaonPlus->Sumw2();
180  hMatchKaonPlus = new TH2F("mKaonPlus","TOF matched K+",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
181  hMatchKaonPlus->Sumw2();
182  hMatchKaonPlus->SetXTitle("p_{T} (GeV/c)");
183  hMatchKaonPlus->SetYTitle("pseudorapidity");
184 
185  //-
186  hMcKaonMin = new TH2F("mcKaonMin","Mc Kaon Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
187  hMcKaonMin->SetXTitle("p_{T} (GeV/c)");
188  hMcKaonMin->SetYTitle("pseudorapidity");
189  hMcKaonMin->Sumw2();
190  hRcKaonMin = new TH2F("rcKaonMin","Rc Kaon Min",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
191  hRcKaonMin->SetXTitle("p_{T} (GeV/c)");
192  hRcKaonMin->SetYTitle("pseudorapidity");
193  hRcKaonMin->Sumw2();
194  hMatchKaonMin = new TH2F("mKaonMin","TOF matched K-",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
195  hMatchKaonMin->Sumw2();
196  hMatchKaonMin->SetXTitle("p_{T} (GeV/c)");
197  hMatchKaonMin->SetYTitle("pseudorapidity");
198 
199  //-
200  hMcProton = new TH2F("mcProton","Mc Proton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
201  hMcProton->SetXTitle("p_{T} (GeV/c)");
202  hMcProton->SetYTitle("pseudorapidity");
203  hMcProton->Sumw2();
204  hRcProton = new TH2F("rcProton","RcProton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
205  hRcProton->SetXTitle("p_{T} (GeV/c)");
206  hRcProton->SetYTitle("pseudorapidity");
207  hRcProton->Sumw2();
208  hMatchProton = new TH2F("mProton","TOF matched proton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
209  hMatchProton->Sumw2();
210  hMatchProton->SetXTitle("p_{T} (GeV/c)");
211  hMatchProton->SetYTitle("pseudorapidity");
212 
213  //-
214  hMcAntiProton = new TH2F("mcAntiProton","Mc Antiproton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
215  hMcAntiProton->SetXTitle("p_{T} (GeV/c)");
216  hMcAntiProton->SetYTitle("pseudorapidity");
217  hMcAntiProton->Sumw2();
218  hRcAntiProton = new TH2F("rcAntiProton","Rc Antiproton",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
219  hRcAntiProton->SetXTitle("p_{T} (GeV/c)");
220  hRcAntiProton->SetYTitle("pseudorapidity");
221  hRcAntiProton->Sumw2();
222  hMatchAntiProton = new TH2F("mAntiProton","TOF matched pbar",mPtBin,mPtMin,mPtMax,mYBin,mYMin,mYMax);
223  hMatchAntiProton->Sumw2();
224  hMatchAntiProton->SetXTitle("p_{T} (GeV/c)");
225  hMatchAntiProton->SetYTitle("pseudorapidity");
226 
227  // Match QA Histograms
228  hTofpSlatIdA0 = new TH1D("SlatIdA0","events per slat",41,0.5,41.5);
229  hTofpSlatIdA1 = new TH1D("SlatIdA1","valid slat",41,0.5,41.5);
230  hTofpSlatIdB1 = new TH1D("SlatIdB1","#tracks match valid slat",41,0.5,41.5);
231  hTofpSlatIdD1 = new TH1D("SlatIdD1","track match per valid slat",41,0.5,41.5);
232  hTofpSlatIdD2 = new TH1D("SlatIdD2","single track match per slat",41,0.5,41.5);
233  hTofpSlatIdE1 = new TH1D("SlatIdE1","one slat for one track match",41,0.5,41.5);
234  hTofpSlatIdE2 = new TH1D("SlatIdE2","recovered from hitprof-weight",41,0.5,41.5);
235  hTofpSlatIdE3 = new TH1D("SlatIdE3","recovered from ss",41,0.5,41.5);
236  hTofpSlatIdE4 = new TH1D("SlatIdE4","recovered from closest hitplane",41,0.5,41.5);
237  hTofpSlatIdE5 = new TH1D("SlatIdE5","total recovered slat per track match",41,0.5,41.5);
238  hTofpSlatIdF1 = new TH1D("SlatIdF1","primary track match per slat",41,0.5,41.5);
239  hTofpSlatHitVecSize = new TH1D("SlatMult","Slat Mult per Track",10,0,10);
240  hTofpHitMap1 = new TH2D("tofpHitMap1","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
241  hTofpHitMap2 = new TH2D("tofpHitMap2","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
242  hTofpHitMap3 = new TH2D("tofpHitMap3","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
243  hTofpHitMap4 = new TH2D("tofpHitMap4","valid hit positions", 500,-250,0, 120,-1.23, -1.08);
244 
245  return StMaker::Init();
246 }
247 
248 
249 Int_t StTofpMcAnalysisMaker::InitRun(int runnumber){
250  gMessMgr->Info("StTofpMcAnalysisMaker -- reinitializing TofGeometry","OS");
251  mTofGeom = new StTofGeometry();
252  mTofGeom->init(this);
253  return kStOk;
254 }
255 
256 Int_t StTofpMcAnalysisMaker::FinishRun(int runnumber){
257  gMessMgr->Info("StTofpMcAnalysisMaker -- cleaning up geometry","OS" );
258  if (mTofGeom) delete mTofGeom;
259  mTofGeom=0;
260  return kStOk;
261 }
262 
263 //_________________________________________________
265 
266  // Get the pointers we need, we have to use the titles we gave them in the
267  // macro. I just used the defaults.
268 
269  // StEvent
270  StEvent* rEvent = 0;
271  rEvent = (StEvent*) GetInputDS("StEvent");
272 
273  // StMcEvent
274  StMcEvent* mEvent = (StMcEvent*) GetDataSet("StMcEvent");
275 
276  // StAssociationMaker
277  StAssociationMaker* assoc = 0;
278  assoc = (StAssociationMaker*) GetMaker("StAssociationMaker");
279 
280  // the Multimaps...
281  rcTpcHitMapType* theHitMap = 0;
282  theHitMap = assoc->rcTpcHitMap();
283  rcTrackMapType* theTrackMap = 0;
284  theTrackMap = assoc->rcTrackMap();
285  mcV0MapType* theMcV0Map = 0;
286  theMcV0Map = assoc->mcV0Map();
287 
288  if (!theHitMap) {
289  gMessMgr->Warning() << "----------WARNING----------\n"
290  << "No Hit Map found for this event!" << endm;
291  return kStWarn;
292  }
293 
294  mcTrackMapType* mcTrackMap = 0;
295  mcTrackMap = assoc->mcTrackMap();
296  if (!mcTrackMap){
297  LOG_INFO << "No mcTrackMap!!! " << endm;
298  return 0;
299  }
300 
301 
302  // Example: look at the position of the primary vertex
303  // Map is not needed for this, but it's a good check,
304  // tracking will not be good if primary vertex was not well placed.
305 
306 
307  // First check whether the Primary Vertex is there at all.
308  StThreeVectorD VertexPos(0,0,0);
309  if (rEvent->primaryVertex()) {
310  VertexPos = rEvent->primaryVertex()->position();
311  LOG_INFO << "Position of Primary Vertex from StEvent:" << endm;
312  LOG_INFO << VertexPos << endm;
313  }
314  else {
315  LOG_INFO << "----------- WARNING ------------" << endm;
316  LOG_INFO << "No primary vertex from StEvent!!" << endm;
317  LOG_INFO << "Assume it is at (0,0,0) " << endm;
318  }
319  LOG_INFO << "Position of Primary Vertex from StMcEvent:" << endm;
320  LOG_INFO << mEvent->primaryVertex()->position() << endm;
321 
322  if (!theTrackMap) {
323  gMessMgr->Warning() << "----------WARNING----------\n"
324  << "No Track Map found for this event!" << endm;
325  return kStWarn;
326  }
327 
328  // -----------------------------------------
329  // --- Stage I: MATCH RC-TRACKS to SLATS ---
330  // -----------------------------------------
331 
332  // Note: why not replaced Stage I w/ the actual TOFpMatchMaker and
333  // use the StTofSlat in StEvent as the reference for Stage IIb?
334 
335 
336  StEvent *event = rEvent;
337  bool mHisto(true);
338  int NTOFP(41);
339  //.........................................................................
340  // A. build vector of candidate slats with valid ADC signals
341  idVector validSlatIdVec;
342  for (int i=0;i<NTOFP;i++){
343  unsigned short slatId = mTofGeom->daqToSlatId(i);
344  hTofpSlatIdA0->Fill(i+1);
345 
346  validSlatIdVec.push_back(slatId);
347  hTofpSlatIdA1->Fill(i+1);
348  }
349  gMessMgr->Info("","OST") << "A: #valid slats: " << validSlatIdVec.size() << endm;
350  // end of Sect.A
351 
352 
353  //.........................................................................
354  // B. loop over global tracks and determine all slat-track matches
355  //
356  tofSlatHitVector allSlatsHitVec;
357  allSlatsHitVec.clear();
358  StSPtrVecTrackNode& nodes = event->trackNodes();
359  int nAllTracks=0;
360  for (unsigned int iNode=0; iNode<nodes.size(); iNode++){
361  tofSlatHitVector slatHitVec;
362  slatHitVec.clear();
363  StTrack *theTrack = nodes[iNode]->track(global);
364 
365  // make sure we have a track, a miniDST might have removed it...
366  if (validTrack(theTrack)){
367  nAllTracks++;
368  StPhysicalHelixD theHelix = trackGeometry(theTrack)->helix();
369  slatHitVec = mTofGeom->tofHelixToArray(theHelix, validSlatIdVec);
370  if (slatHitVec.size()>0 && mHisto) hTofpSlatHitVecSize->Fill(slatHitVec.size());
371 
372  for (size_t ii = 0; ii < slatHitVec.size(); ii++) {
373  slatHitVec[ii].trackIdVec.push_back(iNode);
374  allSlatsHitVec.push_back(slatHitVec[ii]);
375  if (mHisto){
376  int id = mTofGeom->slatIdToDaq(slatHitVec[ii].slatIndex);
377  float xhit = slatHitVec[ii].hitPosition.x();
378  float yhit = slatHitVec[ii].hitPosition.y();
379  float zhit = slatHitVec[ii].hitPosition.z();
380  float phihit = atan2(yhit,xhit);
381  hTofpSlatIdB1->Fill(id);
382  hTofpHitMap1->Fill(zhit,phihit);
383  }
384 
385  if (Debug()){
386  LOG_INFO << "B: trackid=";
387  idVectorIter ij = slatHitVec[ii].trackIdVec.begin();
388  while (ij != slatHitVec[ii].trackIdVec.end()) {LOG_INFO << " " << *ij; ij++;}
389  LOG_INFO << "\tind=" << mTofGeom->slatIdToDaq(slatHitVec[ii].slatIndex)
390  << "\thitprof="<< slatHitVec[ii].hitProfile
391  << "\ts="<<slatHitVec[ii].s << "\tthxy="<<slatHitVec[ii].theta_xy
392  << "\tthzr="<<slatHitVec[ii].theta_zr;
393  if (slatHitVec.size()>1) LOG_INFO << " M" << endm;
394  else LOG_INFO << endm;
395  }
396  }
397  } // existing global track
398  } // loop over nodes
399 
400  gMessMgr->Info("","OST") << "B: #matched/#avail/#total tracknodes: "
401  <<allSlatsHitVec.size() << "/" <<nAllTracks
402  << "/" << nodes.size() << endm;
403  // end of Sect.B
404 
405 
406  //.........................................................................
407  // C Neighbours -- identify crosstalk, geometry shifts (only fill histograms)
408  //
409  // [ ... ]
410  // end of Sect.C
411 
412 
413  //.........................................................................
414  // D. sort hit vectors and deal with (discard) slats matched by multiple tracks
415  //
416  int nSingleHitSlats(0);
417  tofSlatHitVector singleHitSlatsVec;
418  StructSlatHit slatHit;
419 
420  tofSlatHitVector tempVec = allSlatsHitVec;
421  tofSlatHitVector erasedVec = tempVec;
422  while (tempVec.size() != 0) {
423  int nTracks = 0;
424  vector<StThreeVectorD> vPosition;
425  vector<vector<StThreeVectorD> > vLayerHitPositions;
426  vector<Int_t> vHitProfile;
427  vector<Float_t> vS, vTheta_xy, vTheta_zr;
428  idVector trackIdVec;
429 
430  tofSlatHitVectorIter tempIter=tempVec.begin();
431  tofSlatHitVectorIter erasedIter=erasedVec.begin();
432  while(erasedIter!= erasedVec.end()) {
433  if(tempIter->slatIndex == erasedIter->slatIndex) {
434  nTracks++;
435  // save all hit data in temporary vectors
436  trackIdVec.push_back(erasedIter->trackIdVec.back());
437  vPosition.push_back(erasedIter->hitPosition);
438  vLayerHitPositions.push_back(erasedIter->layerHitPositions);
439  vHitProfile.push_back(erasedIter->hitProfile);
440  vS.push_back(erasedIter->s);
441  vTheta_xy.push_back(erasedIter->theta_xy);
442  vTheta_zr.push_back(erasedIter->theta_zr);
443 
444  if (mHisto){
445  float xhit = erasedIter->hitPosition.x();
446  float yhit = erasedIter->hitPosition.y();
447  float zhit = erasedIter->hitPosition.z();
448  float phihit = atan2(yhit,xhit);
449  hTofpHitMap2->Fill(zhit,phihit);
450  }
451 
452  erasedVec.erase(erasedIter);
453  erasedIter--;
454  }
455  erasedIter++;
456  }
457 
458  if (mHisto){
459  int daqId = mTofGeom->slatIdToDaq(tempIter->slatIndex);
460  hTofpSlatIdD1->Fill(daqId);
461  }
462 
463  if (nTracks==1){
464  nSingleHitSlats++;
465  // for singly hit slat, copy data in singleHitSlatsVec
466  slatHit.slatIndex = tempIter->slatIndex;
467  slatHit.hitPosition = vPosition[0];
468  slatHit.layerHitPositions = vLayerHitPositions[0];
469  slatHit.trackIdVec = trackIdVec;
470  slatHit.hitProfile = vHitProfile[0];
471  slatHit.s = vS[0];
472  slatHit.theta_xy = vTheta_xy[0];
473  slatHit.theta_zr = vTheta_zr[0];
474 
475  singleHitSlatsVec.push_back(slatHit);
476 
477  if (mHisto){
478  int daqId = mTofGeom->slatIdToDaq(tempIter->slatIndex);
479  float xhit = slatHit.hitPosition.x();
480  float yhit = slatHit.hitPosition.y();
481  float zhit = slatHit.hitPosition.z();
482  float phihit = atan2(yhit,xhit);
483  hTofpSlatIdD2->Fill(daqId);
484  hTofpHitMap3->Fill(zhit,phihit);
485  }
486 
487  // debugging output
488  if (Debug()) {
489  LOG_INFO << "D: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
490  << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
491  << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:";
492  idVectorIter ij=trackIdVec.begin();
493  while (ij != trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
494  LOG_INFO <<endm;
495  }
496  }
497  else if (nTracks>1){
498  // for multiple hit slats either discard (yes) or
499  // find the most likely candidate.
500  } else
501  gMessMgr->Warning("","OST") << "D: no tracks extrapolate to matched slat ... should not happen!" << endm;
502 
503  tempVec = erasedVec;
504  }
505  gMessMgr->Info("","OST") << "D: #before/#after: " << allSlatsHitVec.size()
506  << "/" << singleHitSlatsVec.size() << endm;
507  //end of Sect.D
508 
509 
510 
511  //.........................................................................
512  // E. sort and deal singleHitSlatsVector for multiple slats associated to single tracks
513  //
514  tofSlatHitVector allMatchedSlatsVec;
515  tempVec = singleHitSlatsVec;
516  erasedVec = tempVec;
517  while (tempVec.size() != 0) {
518  StructSlatHit slatHit;
519  int nSlats = 0;
520  vector<StThreeVectorD> vPosition;
521  vector< vector<StThreeVectorD> > vLayerHitPositions;
522  vector<Int_t> vHitProfile;
523  vector<Float_t> vS, vTheta_xy, vTheta_zr;
524  idVector vTrackId;
525  vector<Int_t> slatIndex;
526 
527  tofSlatHitVectorIter tempIter=tempVec.begin();
528  tofSlatHitVectorIter erasedIter=erasedVec.begin();
529  while(erasedIter!= erasedVec.end()) {
530  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back()) {
531  nSlats++;
532  // save all hit data in temporary vectors
533  slatIndex.push_back(erasedIter->slatIndex);
534  vTrackId.push_back(erasedIter->trackIdVec.back());
535  vPosition.push_back(erasedIter->hitPosition);
536  vLayerHitPositions.push_back(erasedIter->layerHitPositions);
537  vHitProfile.push_back(erasedIter->hitProfile);
538  vS.push_back(erasedIter->s);
539  vTheta_xy.push_back(erasedIter->theta_xy);
540  vTheta_zr.push_back(erasedIter->theta_zr);
541 
542  erasedVec.erase(erasedIter);
543  erasedIter--;
544  }
545  erasedIter++;
546  }
547 
548 
549  if (nSlats==1){
550  // for singly hit slat, copy data in singleHitSlatsVec
551  slatHit.slatIndex = slatIndex[0];
552  slatHit.hitPosition = vPosition[0];
553  slatHit.layerHitPositions = vLayerHitPositions[0];
554  slatHit.trackIdVec.push_back(vTrackId[0]);
555  slatHit.hitProfile = vHitProfile[0];
556  slatHit.s = vS[0];
557  slatHit.theta_xy = vTheta_xy[0];
558  slatHit.theta_zr = vTheta_zr[0];
559  slatHit.matchFlag = 0;
560 
561  allMatchedSlatsVec.push_back(slatHit);
562 
563  if (mHisto){
564  int daqId = mTofGeom->slatIdToDaq(slatIndex[0]);
565  hTofpSlatIdE1->Fill(daqId);
566  }
567 
568  // debugging output
569  if (Debug()) {
570  LOG_INFO << "E: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
571  << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
572  << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:";
573  idVectorIter ij=vTrackId.begin();
574  while (ij != vTrackId.end()) { LOG_INFO << " " << *ij; ij++; }
575  LOG_INFO <<endm;
576  }
577  }
578  else if (nSlats>1){ // for multiple hit slats find the most likely candidate.
579  int thiscandidate(-99);
580  int thisMatchFlag(0);
581 
582  // 1. sort on hitprofile weight
583  int weight(0);
584  vector<int> weightCandidates;
585  thisMatchFlag = 1;
586  if (Debug()) LOG_INFO << "E: find ... weight ";
587  for (int i=0;i<nSlats;i++){
588  int hitWeight = vLayerHitPositions[i].size();
589  if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[i]) << "("<<hitWeight<<")"<<" ";
590  if (hitWeight>weight) {
591  weight=hitWeight;
592  weightCandidates.clear();
593  weightCandidates.push_back(i);
594  } else if (hitWeight == weight)
595  weightCandidates.push_back(i);
596  }
597  if (weightCandidates.size()==1){
598  thiscandidate = weightCandidates[0];
599  int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
600  if (mHisto) hTofpSlatIdE2->Fill(daqId);
601  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
602  }
603 
604  // 2. if still undecided check on ss
605  if (weightCandidates.size()>1){
606  Float_t ss(0);
607  vector<int> ssCandidates;
608  thisMatchFlag = 2;
609  if (Debug()) LOG_INFO << " ss ";
610  for (unsigned int i=0;i<weightCandidates.size();i++){
611  int ii=weightCandidates[i];
612  if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
613  if (vS[ii]>ss){
614  ss = vS[ii];
615  ssCandidates.clear();
616  ssCandidates.push_back(ii);
617  }else if (vS[ii]==ss)
618  ssCandidates.push_back(ii);
619  }
620  if (ssCandidates.size()==1){
621  thiscandidate = ssCandidates[0];
622  int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
623  if (mHisto) hTofpSlatIdE3->Fill(daqId);
624  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
625  }
626 
627  // 3. if still undecided go for closest/first hit
628  if (ssCandidates.size()>1){
629  Int_t hitprof(0);
630  vector<int> profileCandidates;
631  thisMatchFlag = 3;
632  if (Debug()) LOG_INFO << " hprof ";
633  for (unsigned int i=0;i<ssCandidates.size();i++){
634  int ii=ssCandidates[i];
635  if (Debug()) LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
636  if (vHitProfile[ii]>hitprof){
637  hitprof = vHitProfile[ii];
638  profileCandidates.clear();
639  profileCandidates.push_back(ii);
640  }else if (vHitProfile[ii]==hitprof)
641  profileCandidates.push_back(ii);
642  }
643  if (profileCandidates.size()==1){
644  thiscandidate = profileCandidates[0];
645  int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
646  if (mHisto) hTofpSlatIdE4->Fill(daqId);
647  if (Debug()) LOG_INFO << "candidate =" << daqId << endm;
648  }
649  else
650  if (Debug()) LOG_INFO << "none" << endm;
651  }
652 
653 
654  // forget it, and let user know of the non-decision
655  if (thiscandidate == -99 && Debug()){
656  LOG_INFO << "E: ind=";
657  for (unsigned int ii=0;ii<slatIndex.size();ii++)
658  LOG_INFO << mTofGeom->slatIdToDaq(slatIndex[ii]) << " ";
659  LOG_INFO << "\ttrkid:" << vTrackId[0] << " Unable to decide. ";
660  LOG_INFO << "(hitprofs:";
661  for (unsigned int ii=0;ii<slatIndex.size();ii++)
662  LOG_INFO << vHitProfile[ii] << " ";
663  LOG_INFO << " ss:";
664  for (unsigned int ii=0;ii<slatIndex.size();ii++)
665  LOG_INFO << vS[ii] << " ";
666  LOG_INFO << ")" << endm;
667  }
668 
669  }
670 
671 
672  if (thiscandidate>=0){
673  slatHit.slatIndex = slatIndex[thiscandidate];
674  slatHit.hitPosition = vPosition[thiscandidate];
675  slatHit.layerHitPositions = vLayerHitPositions[thiscandidate];
676  //slatHit.trackIdVec.clear();
677  slatHit.trackIdVec.push_back(vTrackId[thiscandidate]);
678  slatHit.hitProfile = vHitProfile[thiscandidate];
679  slatHit.s = vS[thiscandidate];
680  slatHit.theta_xy = vTheta_xy[thiscandidate];
681  slatHit.theta_zr = vTheta_zr[thiscandidate];
682  slatHit.matchFlag = thisMatchFlag;
683 
684  allMatchedSlatsVec.push_back(slatHit);
685 
686  if (mHisto){
687  int daqId = mTofGeom->slatIdToDaq(slatIndex[thiscandidate]);
688  hTofpSlatIdE5->Fill(daqId);
689  }
690 
691  // debugging output
692  if (Debug()) {
693  LOG_INFO << "E: ind=" << mTofGeom->slatIdToDaq(slatHit.slatIndex)
694  << "\thitprof="<< slatHit.hitProfile << "\ts="<<slatHit.s
695  << "\tthxy="<<slatHit.theta_xy << "\tthzr="<<slatHit.theta_zr << "\ttrackid:"
696  << vTrackId[thiscandidate] << endm;
697  }
698  }
699  } else
700  gMessMgr->Warning("","OS") << "E: no slats belong to this track ... should not happen!" << endm;
701 
702  tempVec = erasedVec;
703  }
704  gMessMgr->Info("","OST") << "E: #before/#after: " << singleHitSlatsVec.size()
705  << "/" << allMatchedSlatsVec.size() << endm;
706  // end of Sect.E
707 
708 
709  // --------------------------------------------------------------------//
710  // Note: section F is not required: the primary-track requirement will //
711  // be enforced lateron in the MC-RC track association. //
712  // --------------------------------------------------------------------//
713 
714 
715  //.........................................................................
716  // F. perform further selection and
717  // fill valid track histograms and SlatCollection
718  //
719  //StTofSlatCollection *mSlatCollection = new StTofSlatCollection;
720  int nValidSinglePrimHitSlats(0); //nValidSingleHitSlats(0)
721 
722  for (size_t ii=0; ii < allMatchedSlatsVec.size(); ii++){
723  int daqId = mTofGeom->slatIdToDaq(allMatchedSlatsVec[ii].slatIndex);
724  //int jj = daqId-1;
725 
726  if (allMatchedSlatsVec[ii].trackIdVec.size()!=1)
727  gMessMgr->Warning("","OST") << "F: WHAT!?! mult.matched slat in single slat list " << daqId << endm;
728 
729  // 1. fill valid single track AND valid tdc histograms
730  //if (validTdc(mTofpTdc[jj])) nValidSingleHitSlats++;
731 
732  // get track-id from slat hit vector
733  unsigned int trackNode = allMatchedSlatsVec[ii].trackIdVec[0];
734  StTrack *theTrack = nodes[trackNode]->track(primary);
735 
736  // 2. continue only if the (primary) track exists
737  if (validTofTrack(theTrack)){
738  nValidSinglePrimHitSlats++;
739  if (mHisto){
740  float xhit = allMatchedSlatsVec[ii].hitPosition.x();
741  float yhit = allMatchedSlatsVec[ii].hitPosition.y();
742  float zhit = allMatchedSlatsVec[ii].hitPosition.z();
743  float phihit = atan2(yhit,xhit);
744  hTofpHitMap4->Fill(zhit,phihit);
745  hTofpSlatIdF1->Fill(daqId);
746  }
747 
748  //--- store number of hits per track
749  int nHitsPerTrack = theTrack->topologyMap().numberOfHits(kTpcId);
750  //if (mHisto) hTofpNumberOfTrackHits->Fill(nHitsPerTrack);
751 
752  // select the apropriate track geometry
753  const StTrackGeometry *theTrackGeometry = trackGeometry(theTrack);
754 
755  //--- get momentum from track
756  const StThreeVectorF momentum = theTrackGeometry->momentum();
757  //if (mHisto) hTofpPtTrack->Fill(momentum.perp());
758 
759  //--- calculate flight path
760  double pathLength = tofPathLength(&event->primaryVertex()->position(),
761  &allMatchedSlatsVec[ii].hitPosition,
762  theTrackGeometry->helix().curvature());
763 
764 
765  //--- calculate local hit position on slat based first, last and middle plane
766  // (middle plane is the average of first and last plane, which is mathematically
767  // the same as SlatHitVec.hitPosition ... )
768  StThreeVectorD *pInnerLayer, *pOuterLayer;
769  pInnerLayer = &(*(allMatchedSlatsVec[ii].layerHitPositions.begin()));
770  pOuterLayer = &(*(allMatchedSlatsVec[ii].layerHitPositions.end() - 1));
771 
772  //--- dig out from the dedx and rich pid traits (only for Debug mode)
773  float dedx(0.), cherang(0);
774  int dedx_np(0), cherang_nph(0);
775  if (Debug()){
776  StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
777  for (unsigned int it=0;it<traits.size();it++){
778  if (traits[it]->detector() == kTpcId){
779  StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
780  if (pid && pid->method() ==kTruncatedMeanId){
781  dedx = pid->mean();
782  dedx_np = pid->numberOfPoints();
783  }
784  } else if (traits[it]->detector() == kRichId){
785  StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
786  if (pid){
787  StRichSpectra* pidinfo = pid->getRichSpectra();
788  if (pidinfo && pidinfo->getCherenkovPhotons()>2){
789  cherang = pidinfo->getCherenkovAngle();
790  cherang_nph = pidinfo->getCherenkovPhotons();
791  }
792  }
793  }
794  }
795  }
796 
797  //--- calculate local hit position on slat based on average hitposition
798  //float localHitPos = mTofGeom->slatHitPosition(&allMatchedSlatsVec[ii].hitPosition);
799 
800  // Fill TOF Slat Collection
801  //StTofSlat *tofSlat = new StTofSlat(jj+1,(int)mTofpAdc[jj],(int)mTofpTdc[jj],theTrack);
802  //StTofSlat *tofSlat = new StTofSlat(jj+1,(int)mTofpAdc[jj],(int)mTofpTdc[jj],theTrack,
803  // localHitPos, allMatchedSlatsVec[ii].hitProfile,
804  // allMatchedSlatsVec[ii].matchFlag);
805  //tofSlat->setPosition(allMatchedSlatsVec[ii].hitPosition);
806  //mSlatCollection->push_back(tofSlat);
807 
808  // dump debug data
809  if (Debug()){
810  LOG_INFO << "F: ind=" << mTofGeom->slatIdToDaq(allMatchedSlatsVec[ii].slatIndex)
811  << "\ttrackid:";
812  idVectorIter ij=allMatchedSlatsVec[ii].trackIdVec.begin();
813  while (ij != allMatchedSlatsVec[ii].trackIdVec.end()) { LOG_INFO << " " << *ij; ij++; }
814  LOG_INFO << "\tR=" << 1/(theTrackGeometry->helix().curvature())
815  << "\tpT=" << momentum.perp() << "\tp=" << momentum.mag()
816  << "\thits="<< nHitsPerTrack << "\ts="<< pathLength
817  << "\t#fitp=" <<theTrack->fitTraits().numberOfFitPoints(kTpcId)
818  << "\t#trkp=" <<theTrack->detectorInfo()->numberOfPoints(kTpcId)
819  << " \tdedx=" << dedx
820  << " \tdca="<< theTrack->geometry()->helix().distance(event->primaryVertex()->position());
821  if (cherang!=0) LOG_INFO << " \trich="<< cherang << " (" << cherang_nph << ")";
822  LOG_INFO << endm;
823  }
824 
825  } // track exists
826  }
827 
828 
829 
830  // ---------------------------------------------
831  // --- Stage II: FIND RC-MC pairs (TPC eff.) ---
832  // ---------------------------------------------
833 
834  // ----
835  // IIa. loop over mcTracks and check for associated rcTrack (TPC tracking eff.)
836  // IIb. check matching results for rcTrack (TOFp matching eff.)
837  // ----
838 
839 
840  const StSPtrVecMcTrack& mcTracks = mEvent->primaryVertex()->daughters();
841  if (Debug()) LOG_INFO << "mcTrack Loop: "<< mcTracks.size() << " entries" << endm;
842  StMcTrack* mcTrack;
843  const StTrack* rcTrack = 0;
844 
845  int nRecon(0), nMatch(0);
846 
847  //--- Loop over MC tracks
848  for (StMcTrackConstIterator mcItr = mcTracks.begin(); mcItr != mcTracks.end(); mcItr++) {
849  mcTrack = *mcItr;
850  if (!mcTrack) {
851  gMessMgr->Warning() << "No McTrack?!? Should not happen!"<<endm;
852  continue;
853  }
854 
855 
856  //-- Fill histograms with McTrack data
857  StThreeVectorF mcP = mcTrack->momentum();
858  if (mcTrack->particleDefinition()->pdgEncoding()== 11) hMcElectron->Fill(mcP.perp(),mcP.pseudoRapidity());
859  else if (mcTrack->particleDefinition()->pdgEncoding()== 211) hMcPionPlus->Fill(mcP.perp(),mcP.pseudoRapidity());
860  else if (mcTrack->particleDefinition()->pdgEncoding()==-211) hMcPionMin->Fill(mcP.perp(),mcP.pseudoRapidity());
861  else if (mcTrack->particleDefinition()->pdgEncoding()== 321) hMcKaonPlus->Fill(mcP.perp(),mcP.pseudoRapidity());
862  else if (mcTrack->particleDefinition()->pdgEncoding()==-321) hMcKaonMin->Fill(mcP.perp(),mcP.pseudoRapidity());
863  else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hMcProton->Fill(mcP.perp(),mcP.pseudoRapidity());
864  else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hMcAntiProton->Fill(mcP.perp(),mcP.pseudoRapidity());
865  else LOG_INFO << "mcTracks has wrong pdgEncoding: "<< mcTrack->particleDefinition()->pdgEncoding() << endm;
866 
867  //-- get associated reconstructed rctrack
868  rcTrack = partnerTrack(mcTrackMap,mcTrack);
869  if (!rcTrack) continue;
870  if (rcTrack->flag()<0) continue;
871 
872  //-- require primary track
873  rcTrack = rcTrack->node()->track(primary);
874  if (!rcTrack) continue;
875 
876  nRecon++;
877 
878  //-- Fill histograms for the TPC-reconstructed tracks
879  // (select the approriate trackgometry)
880 
881  //StThreeVectorF rcP = rcTrack->geometry()->momentum();
882  StThreeVectorF rcP = trackGeometry(rcTrack)->momentum();
883 
884  if (mcTrack->particleDefinition()->pdgEncoding()== 11) hRcElectron->Fill(rcP.perp(),rcP.pseudoRapidity());
885  else if (mcTrack->particleDefinition()->pdgEncoding()== 211) hRcPionPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
886  else if (mcTrack->particleDefinition()->pdgEncoding()== -211) hRcPionMin->Fill(rcP.perp(),rcP.pseudoRapidity());
887  else if (mcTrack->particleDefinition()->pdgEncoding()== 321) hRcKaonPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
888  else if (mcTrack->particleDefinition()->pdgEncoding()== -321) hRcKaonMin->Fill(rcP.perp(),rcP.pseudoRapidity());
889  else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hRcProton->Fill(rcP.perp(),rcP.pseudoRapidity());
890  else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hRcAntiProton->Fill(rcP.perp(),rcP.pseudoRapidity());
891  else LOG_INFO << "cannot fill rc histo" << endm;
892  float diff = (mcP.mag()-rcP.mag()) / mcP.mag();
893  if (mcTrack->particleDefinition()->pdgEncoding()==-211) hMomResPtPion->Fill(mcP.perp(),diff);
894 
895 
896 
897 
898  //-- Loop over matched slats from Stage 1 (allMatchedSlatsVec) and look for similar trackKey
899  //
900 
901  if (Debug()) LOG_INFO << "RC/MC track Id="<< rcTrack->key() << "/" << mcTrack->key() << endm;
902 
903  // Note on trackId vs. StTrack->key():
904  // these are *not* the same. trackId (based on the StEvent->trackNodes
905  // array index) is used internally to quickly locate tracks in the trackNodes
906  // array. StTrack->key(), however, should *always* be used to identify tracks
907  // outside the scope of this local code.
908  // To get the proper Key() from the trackId use this:
909  // trackKey = StEvent->trackNodes[trackId]->track()->key();
910 
911  for (tofSlatHitVectorIter matchedSlat=allMatchedSlatsVec.begin();
912  matchedSlat != allMatchedSlatsVec.end(); matchedSlat++){
913 
914  //- All matches from Stage I should have one and only one associated track
915  if (matchedSlat->trackIdVec.size()!=1){
916  gMessMgr->Warning() << "Slat with " << matchedSlat->trackIdVec.size()
917  << "track matches? Should not happen!" << endm;
918  continue;
919  }
920 
921  //- Map trackId to trackKey
922  unsigned int thisTrackId = matchedSlat->trackIdVec[0];
923  unsigned int thisKey = nodes[thisTrackId]->track(global)->key();
924 
925  if (Debug()) LOG_INFO << " match trackKey (trackId): " << thisKey
926  << " (" << thisTrackId << ")";
927 
928  if (thisKey == rcTrack->key()){
929  if (Debug()) LOG_INFO << "RC-MC-TOFp match";
930 
931  nMatch++;
932 
933  //- Fill histograms for the reconstructed and matched tracks
934  if (mcTrack->particleDefinition()->pdgEncoding()== 11) hMatchElectron->Fill(rcP.perp(),rcP.pseudoRapidity());
935  else if (mcTrack->particleDefinition()->pdgEncoding()== 211) hMatchPionPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
936  else if (mcTrack->particleDefinition()->pdgEncoding()== -211) hMatchPionMin->Fill(rcP.perp(),rcP.pseudoRapidity());
937  else if (mcTrack->particleDefinition()->pdgEncoding()== 321) hMatchKaonPlus->Fill(rcP.perp(),rcP.pseudoRapidity());
938  else if (mcTrack->particleDefinition()->pdgEncoding()== -321) hMatchKaonMin->Fill(rcP.perp(),rcP.pseudoRapidity());
939  else if (mcTrack->particleDefinition()->pdgEncoding()== 2212) hMatchProton->Fill(rcP.perp(),rcP.pseudoRapidity());
940  else if (mcTrack->particleDefinition()->pdgEncoding()==-2212) hMatchAntiProton->Fill(rcP.perp(),rcP.pseudoRapidity());
941  else LOG_INFO << "cannot fill match histo" << endm;
942  }
943  if (Debug()) LOG_INFO << endm;
944 
945  }
946 
947 
948 
949  // Alternative check ... reextrapolate the track and compare slatIds
950 
951 // // A. make sure the track makes it to the ctb
952 // StPhysicalHelixD thisHelix = rcTrack->geometry()->helix();
953 // //pairD pairL = thisHelix.pathLength(mTofGeom->tofParam().r + mTofGeom->tofParam().counter_thickness);
954 // //if(fabs(pairL.first)>=500.0 && fabs(pairL.second)>=500.0) continue;
955 //
956 // // B. only consider slat extrapolations (no Nslat or Ntrack/slat cuts)
957 // tofSlatHitVector hitVec = mTofGeom->tofHelixToArray(thisHelix, validSlatIdVec);
958 // if (hitVec.size()==0) continue;
959 //
960 // {
961 // LOG_INFO << "TOF HitVec ("<< hitVec.size() << ") slatIn reads:";
962 // for (tofSlatHitVectorIter ii=hitVec.begin();ii!=hitVec.end();ii++)
963 // LOG_INFO << " " << ii->slatIndex;
964 // LOG_INFO << endm;
965 // }
966 //
967 //
968 //
969 // // C i. loop over AllHits and select tracks w/ one or more single hit slats
970 // // ii. select on hit profile
971 // bool singleHit(false),profiledSingleHit(false);
972 // tofSlatHitVectorIter thisHit=hitVec.begin();
973 // while (thisHit!=hitVec.end()){
974 // tofSlatHitVectorIter allHits = allMatchedSlatsVec.begin();
975 // while (allHits!=allMatchedSlatsVec.end()){
976 // if (thisHit->slatIndex == allHits->slatIndex){
977 // LOG_INFO << "TOF found it in allHits "
978 // << allHits->trackIdVec.size() << endm;
979 // if (allHits->trackIdVec.size()==1){
980 // singleHit=true;
981 // if (thisHit->hitProfile == 31) profiledSingleHit=true;
982 // }
983 // }
984 // allHits++;
985 // }
986 // thisHit++;
987 // }
988 // if (!singleHit) {LOG_INFO << "TOF multiple hit"<<endm;continue;}
989 // LOG_INFO << "TOF single hit"<< endm;
990 // //if (!profiledSingleHit) continue;
991 // //LOG_INFO << "TOF profiled single hit"<< endm;
992 
993 
994  }
995  gMessMgr->Info("","OST") << "Embedding: #mcTrack/#recon/#match: " << mcTracks.size()
996  <<"/"<< nRecon << "/" << nMatch << endm;
997 
998  return kStOK;
999 }
1000 
1001 
1002 //---------------------------------------------------------------------------
1003 bool StTofpMcAnalysisMaker::validTrack(StTrack *track){
1004  // 1. no track, no go.
1005  if (!track) return false;
1006 
1007  // 2. track quality flag, should be >0
1008  if (track->flag()<=0) return false;
1009 
1010  // 3. minimum #hits per track
1011  if (track->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return false;
1012 
1013  // 4. rlast (not implemented, yet)
1014 
1015  return true;
1016 }
1017 
1018 
1019 bool StTofpMcAnalysisMaker::validTofTrack(StTrack *track){
1020  // select valid tracks for time-of-flight calculations
1021 
1022  // 1. track must exist
1023  if (!track) return false;
1024 
1025  // 2. track must be a primary track
1026  if (!dynamic_cast<StPrimaryTrack*>(track)) return false;
1027 
1028  // 3. DCA cut (obsolete?)
1029  double DCA= track->impactParameter();
1030  //int charge = track->geometry()->charge();
1031  //if (mHisto) hTofpDCATrackprimVertex->Fill(DCA*charge);
1032  double mMaxDCA = 999.;
1033  if (DCA > mMaxDCA) {
1034  gMessMgr->Info("","OST") << "dca>max:" << DCA<< endm;
1035  return false;
1036  }
1037 
1038  return true;
1039 }
1040 
1041 
1042 const StTrackGeometry* StTofpMcAnalysisMaker::trackGeometry(const StTrack* track){
1043  // returns apropriate StTrackGeometry (standard or outerGeometry)
1044  if (!track) return 0;
1045  const StTrackGeometry *thisTrackGeometry;
1046  if (mOuterTrackGeometry)
1047  thisTrackGeometry = track->outerGeometry();
1048  else
1049  thisTrackGeometry = track->geometry();
1050  return thisTrackGeometry;
1051 }
1052 
1053 //--------------------------------------------------------------------------
1054 /* $Log: StTofpMcAnalysisMaker.cxx,v $
1055  * Revision 1.6 2018/02/26 23:26:51 smirnovd
1056  * StTof: Remove outdated ClassImp macro
1057  *
1058 /* Revision 1.5 2018/02/26 23:13:20 smirnovd
1059 /* Move embedded CVS log messages to the end of file
1060 /*
1061  * Revision 1.4 2011/04/03 15:52:57 fisyak
1062  * Fix effect of constness in StAssociationMaker
1063  *
1064  * Revision 1.3 2007/04/17 23:11:04 dongx
1065  * replaced with standard STAR Loggers
1066  *
1067  * Revision 1.2 2005/10/06 19:58:15 fisyak
1068  * Adjust for persistent StMcEvent
1069  *
1070  * Revision 1.1 2004/03/16 04:58:55 geurts
1071  * *** empty log message ***
1072  */
TH1D * hTofpSlatIdE4
recovered from ss
TH1D * hTofpSlatIdA0
2D tray hit positions (F)
TH1D * hTofpSlatIdB1
valid slat
TH2F * hRcElectron
TPC reconstructed pbar.
TH1F * hMultRef
MC/RC momentum resolution.
TH2F * hRcKaonMin
TPC reconstructed K+.
TH2F * hRcAntiProton
TPC reconstructed p+.
TH2F * hMatchKaonMin
TOF matched K+.
TH1D * hTofpSlatIdE2
one slat for one track match
TH1D * hTofpSlatIdE5
recovered from closest hitplane
StTofGeometry * mTofGeom
slat mult per StTrack
TH2D * hTofpHitMap4
2D tray hit positions (D)
TH1D * hTofpSlatIdA1
events per slat
TH2F * hMcKaonMin
montecarlo K+
TH2D * hTofpHitMap3
2D tray hit positions (pre-D)
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
TH2F * hMatchPionMin
TOF matched pi+.
TH1D * hTofpSlatHitVecSize
primary track match per slat
TH2F * hMcKaonPlus
montecarlo pi-
TH1D * hTofpSlatIdE1
single track match per slat
TH2F * hMcPionMin
montecarlo pi+
TH1D * hTofpSlatIdD1
#tracks match valid slat
TH2F * hMatchProton
TOF matched K-.
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
TH2F * hMatchPionPlus
TPC reconstructed e-.
tofSlatHitVector tofHelixToArray(const StPhysicalHelixD &helix, idVector slatIdVec)
finds slats in an array of trays which are crossed by a track-helix.
TH2F * hMatchKaonPlus
TOF matched pi-.
TH2F * hRcProton
TPC reconstructed K-.
TH2F * hMatchElectron
TOF matched pbar.
Time-of-Flight Geometry Utilities.
TH2F * hMcAntiProton
montecarlo p+
TH2D * hTofpHitMap1
multiplicity reference
TH1D * hTofpSlatIdE3
recovered from hitprof-weight
TH2F * hMomResPtPion
TOF matched e-.
TH2F * hRcPionPlus
montecarlo e-
TH2F * hRcPionMin
TPC reconstructed pi+.
TH2F * hMcElectron
montecarlo pbar
TH1D * hTofpSlatIdF1
total recovered slat per track match
Definition: Stypes.h:42
Definition: Stypes.h:40
TH1D * hTofpSlatIdD2
track match per valid slat
TH2F * hRcKaonPlus
TPC reconstructed pi-.
TH2F * hMatchAntiProton
TOF matched p+.
virtual Int_t Finish()
Definition: StMaker.cxx:776
rcTpcHitMapType * rcTpcHitMap()
Diff btw global r and phi coords of FTPC hits.
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Definition: Stypes.h:41
TH2D * hTofpHitMap2
2D tray hit positions (B)
TH2F * hMcProton
montecarlo K-
void init()
initialize geometry class from XDF file and set-up DAQ/Slat mappings