StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMtdMatchMaker.cxx
1 /*******************************************************************
2  * $Id: StMtdMatchMaker.cxx,v 1.42 2018/12/07 13:16:48 marr Exp $
3  * Author: Bingchu Huang
4  *****************************************************************
5  *
6  * Description: MTD Match Maker to do the matching between the
7  * fired cells and TPC tracks
8  *
9  *****************************************************************
10  *
11  * $Log: StMtdMatchMaker.cxx,v $
12  * Revision 1.42 2018/12/07 13:16:48 marr
13  * Bug fix: check the existence of vertex when analyzing StEvent
14  *
15  * Revision 1.41 2018/12/06 18:11:13 marr
16  * Improvement: extrapolate tracks to the proper primary vertex when available. This eliminates large negative dTof values
17  *
18  * Revision 1.40 2018/09/04 19:29:14 marr
19  * Use the pairD definition in StHelixD.hh
20  *
21  * Revision 1.39 2017/07/31 14:19:14 marr
22  * Add protection for BL9, installed in 2017 for test, since the geometry file
23  * does not include this backleg.
24  *
25  * Revision 1.38 2017/03/08 20:48:54 marr
26  * 1) Add a new data member mYear to indicate run year
27  * 2) Invoke appropriate functions in StMtdGeometry class to calculate local y
28  * to make the class backward compatible
29  *
30  * Revision 1.37 2017/02/13 02:57:10 marr
31  * From 2017, do not move BL 8&24 along y direction by hand since this is already
32  * done in the geometry file. Calibration, production and analysis should use
33  * the new version consistently.
34  *
35  * Revision 1.36 2016/08/05 16:12:24 marr
36  * Add MTD hit IdTruth to avoid applying dy shift for BL 8 and 24 for MC hits
37  *
38  * Revision 1.35 2016/07/28 14:31:23 marr
39  * Fix coverity check: initialization of data member
40  *
41  * Revision 1.34 2016/07/27 15:46:34 marr
42  * Fix coverity check: initialization of data members
43  *
44  * Revision 1.33 2015/10/16 19:04:55 marr
45  * Remove filling trees
46  *
47  * Revision 1.32 2015/07/24 16:21:24 marr
48  * Create geometry in InitRun(), which is needed to use StMtdGeometry class.
49  *
50  * Revision 1.31 2015/07/24 16:03:25 marr
51  * *** empty log message ***
52  *
53  * Revision 1.30 2015/07/10 16:07:40 marr
54  * Add the distance along radius to the calculation of the distance between
55  * a MTD hit and a projected track position
56  *
57  * Revision 1.29 2015/05/01 01:55:02 marr
58  * Fix the geometry of shifted backleg 8 and 24
59  *
60  * Revision 1.28 2015/04/24 19:55:16 marr
61  * Add a member function cleanUpMtdPidTraits() to clean up the MTD pidTraits for
62  * all global and primary tracks before the matching process. This is needed when
63  * running MuDst in afterburner mode.
64  *
65  * Revision 1.27 2015/04/10 18:30:42 marr
66  * Remove lines that are commented out
67  *
68  * Revision 1.26 2015/04/10 18:25:59 marr
69  * 1. Use global coordinates, instead of local ones, to calculate the distances
70  * between projected track positions and MTD hits. Using local coordinates causes
71  * a unphysical shoulder at large DeltaZ.
72  * 2. Fix the following bug: if a track can be projected to two adjacent modules,
73  * the position in the module with a smaller module number gets dropped before
74  * matching. This caused an asymmetry in the local z distribution of matched tracks.
75  *
76  * Revision 1.25 2015/04/07 16:25:16 marr
77  * 1. Make use the constants defined in StMtdConstants
78  * 2. Disable the print-out when running in MuDst mode
79  * 3. Cleaning up
80  *
81  * Revision 1.24 2014/11/12 22:27:23 marr
82  * Clean up the matching information when running on StEvent in afterburner mode
83  *
84  * Revision 1.23 2014/09/20 04:25:20 marr
85  * Assign matching information to all the primary tracks
86  *
87  * Revision 1.22 2014/09/18 22:03:01 marr
88  * Do not set default geometry tag
89  *
90  * Revision 1.21 2014/09/09 14:00:39 marr
91  * Fill the expected time-of-flight calculated via track extrapolation
92  *
93  * Revision 1.20 2014/08/18 17:44:45 marr
94  * Add assert() statement to abort if the loaded geometry is wrong
95  *
96  * Revision 1.19 2014/07/29 19:23:47 marr
97  * Remove the dependency on "StarGenerator/StarLight/starlightconstants.h" as it is not needed anymore.
98  *
99  * Revision 1.18 2014/07/24 02:53:04 marr
100  * 1) Add log info of the matched track-hit pair
101  * 2) Set DeltaY and DeltaZ in PidTraits
102  *
103  * Revision 1.17 2014/07/18 15:52:00 marr
104  * Initialize trgTime[2]
105  *
106  * Revision 1.16 2014/07/16 15:43:53 huangbc
107  * use mMtdGeom->SetLockBField();
108  *
109  * Revision 1.15 2014/07/15 02:01:04 huangbc
110  * Implement multi-tracks to 1 hit matching algo. Set neighbor module matching and 3 extra cells as default.
111  *
112  * Revision 1.14 2014/07/10 20:50:35 huangbc
113  * Use new MTD geometry class. Load geometry volume from geant.
114  * Choose closest one for multi-tracks which associated with same hit.
115  * Remove dca cut
116  *
117  * Revision 1.13 2014/07/08 03:09:07 huangbc
118  * Change dca<10 to dca2Beam_R<10.
119  *
120  * Revision 1.12 2014/06/19 19:16:27 huangbc
121  * Fixed an issue of reading SL12d production data. Add expTof for MTD pidtraits.
122  *
123  * Revision 1.11 2014/04/16 02:23:39 huangbc
124  * 1. fixed a bug of un-initialized variable nDedxPts in MtdTrack construction function.
125  * 2. reoriganized project2Mtd function. Made it more readable.
126  * 3. save pathlengths of extrapolated tracks.
127  * 4. tot selection < 40 ns. drop off maximum tot selection.
128  * 5. add hits <-> track index association in mMuDstIn=true mode.
129  *
130  * Revision 1.10 2014/03/11 22:18:11 geurts
131  * corrected pvtx retrieval in StEvent environment
132  *
133  * Revision 1.9 2014/03/11 02:15:53 geurts
134  * Protect against potentially non-existing primary vertex in StEvent [Bingchu]
135  *
136  * Revision 1.8 2013/12/09 22:53:25 geurts
137  * update: enable filling of MTD Pid traits and include a few more protections against zero-pointers [Bingchu]
138  *
139  * Revision 1.7 2013/11/25 16:10:43 geurts
140  * Remove AddHist for uninitialized histogram [Jason Webb]
141  *
142  * Revision 1.6 2013/11/19 22:30:30 jeromel
143  * Added name
144  *
145  * Revision 1.5 2013/11/19 00:17:16 geurts
146  * include protection against zero triggerIdCollection() pointers. This is relevant when using simulated data as input.
147 *
148 * Revision 1.4 2013/04/25 14:52:13 geurts
149 * Minor adjustments that address SL44 compiler warnings
150 *
151 * Revision 1.3 2013/04/18 21:01:10 geurts
152 * Bugfix (RT#2575): protection against events that have tracks, but no vertex.
153 * - Warning messages
154 * - Consistently reset variables that depend on a vertex to -9999.
155 *
156 *
157 *
158 *******************************************************************/
159 #include <iostream>
160 #include <fstream>
161 #include <stdlib.h>
162 #include <string>
163 #include <vector>
164 #include <math.h>
165 #include <cmath>
166 
167 #include "TH1.h"
168 #include "TH2.h"
169 #include "TProfile.h"
170 #include "TChain.h"
171 #include "TSystem.h"
172 #include "TTree.h"
173 #include "TBranch.h"
174 #include "TGeoManager.h"
175 
176 #include "StEvent.h"
177 #include "StTrackNode.h"
178 #include "StContainers.h"
179 #include "StPrimaryVertex.h"
180 #include "StVertex.h"
181 #include "StMeasuredPoint.h"
182 #include "StDedxPidTraits.h"
183 #include "StTrackPidTraits.h"
184 #include "StTrackGeometry.h"
185 #include "StTrackDetectorInfo.h"
186 #include "StGlobalTrack.h"
187 #include "StParticleTypes.hh"
188 #include "StThreeVector.hh"
189 #include "StThreeVectorF.hh"
190 #include "StThreeVectorD.hh"
191 #include "StLorentzVectorD.hh"
192 #include "StGlobals.hh" // for PR()
193 #include "StGetConfigValue.hh"
194 #include "StTimer.hh"
195 #include "StMemoryInfo.hh"
196 #include "StMessMgr.h"
197 #include "StChain.h"
198 #include "St_DataSet.h"
199 #include "St_DataSetIter.h"
200 #include "StEventTypes.h"
201 #include "StTpcDedxPidAlgorithm.h"
202 #include "StParticleDefinition.hh"
203 #include "StPhysicalHelix.hh"
204 #include "StHelixD.hh"
205 
206 #include "PhysicalConstants.h"
207 #include "SystemOfUnits.h" // has "tesla" in it
208 
209 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
210 #include "StMuDSTMaker/COMMON/StMuEvent.h"
211 #include "StMuDSTMaker/COMMON/StMuDst.h"
212 #include "StMuDSTMaker/COMMON/StMuTrack.h"
213 #include "StMuDSTMaker/COMMON/StMuMtdCollection.h"
214 #include "StMuDSTMaker/COMMON/StMuMtdRawHit.h"
215 #include "StMuDSTMaker/COMMON/StMuMtdHeader.h"
216 #include "StMuDSTMaker/COMMON/StMuMtdHit.h"
217 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
218 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
219 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
220 #include "StMuDSTMaker/COMMON/StMuMtdPidTraits.h"
221 
222 #include "StEventMaker/StEventMaker.h"
223 #include "StAssociationMaker/StAssociationMaker.h"
224 #include "StMcEventMaker/StMcEventMaker.h"
225 #include "StDetectorDbMaker/St_MagFactorC.h"
226 #include "tables/St_vertexSeed_Table.h"
227 
228 #include "StBTofPidTraits.h"
229 #include "StBTofCollection.h"
230 #include "StBTofHit.h"
231 #include "StBTofRawHit.h"
232 #include "StBTofHeader.h"
233 #include "StBTofUtil/tofPathLength.hh"
234 #include "StBTofUtil/StBTofGeometry.h"
235 #include "StBTofUtil/StBTofDaqMap.h"
236 #include "StBTofUtil/StBTofHitCollection.h"
237 #include "StMtdPidTraits.h"
238 
239 #include "StMtdUtil/StMtdGeometry.h"
240 #include "StMtdMatchMaker.h"
241 
242 using namespace std;
243 
244 
246 StMtdMatchMaker::StMtdMatchMaker(const Char_t *name): StMaker(name)
247 {
248  mBeamHelix = NULL;
249  doPrintMemoryInfo = kFALSE;
250  doPrintCpuInfo = kFALSE;
251  mCosmicFlag=kFALSE;
252 
253  mMinFitPointsPerTrack=15;
256  mMinEta=-0.8;
257  mMaxEta=0.8;
258  mMinPt = 1.0;
259 
260  mMuDstIn = kFALSE;
261  mHisto = kFALSE;
262  memset(mVDrift,0,sizeof(mVDrift));
263  mnNeighbors = kTRUE;
264  mNExtraCells = 3;
265  index2Primary.clear();
266 
267  mELossFlag=kTRUE;
268  mLockBField = kFALSE;
269  mMtdGeom = 0;
270 
271  ngTracks = 0;
272  mEvent = NULL;
273  mMuDst = NULL;
274  mYear = -1;
275  mGeomTag = "";
276 
277  fZReso = new TF1("fZReso","sqrt([0]/x/x+[1])",0,100);
278  fZReso->SetParameters(148.7,1.654); //cm
279  fPhiReso = new TF1("fPhiReso","sqrt([0]/x/x+[1])",0,100);
280  fPhiReso->SetParameters(9.514e-4,7.458e-6); //rad
281 
282  mEventCounterHisto = NULL;
283  mCellsMultInEvent = NULL;
284  mHitsMultInEvent = NULL;
285  mHitsPrimaryInEvent = NULL;
286  mHitsGlobalInEvent = NULL;
287  mHitsMultPerTrack = NULL;
288  mHitsPosition = NULL;
289  for(int i=0; i<mNBacklegs; i++)
290  {
291  mDaqOccupancy[i] = NULL;
292  mDaqOccupancyProj[i] = NULL;
293  mHitCorr[i] = NULL;
294  mHitCorrModule[i] = NULL;
295  mDeltaHitFinal[i] = NULL;
296  }
297  mTrackPtEta = NULL;
298  mTrackPtPhi = NULL;
299  mTrackNFitPts = NULL;
300  mTrackdEdxvsp = NULL;
301  mNSigmaPivsPt = NULL;
302  mCellsPerEventMatch1 = NULL;
303  mHitsPerEventMatch1 = NULL;
304  mCellsPerTrackMatch1 = NULL;
305  mTracksPerCellMatch1 = NULL;
306  mDaqOccupancyMatch1 = NULL;
307  mDeltaHitMatch1 = NULL;
308  mCellsPerEventMatch2 = NULL;
309  mHitsPerEventMatch2 = NULL;
310  mCellsPerTrackMatch2 = NULL;
311  mTracksPerCellMatch2 = NULL;
312  mDaqOccupancyMatch2 = NULL;
313  mDeltaHitMatch2 = NULL;
314  mCellsPerEventMatch3 = NULL;
315  mHitsPerEventMatch3 = NULL;
316  mCellsPerTrackMatch3 = NULL;
317  mTracksPerCellMatch3 = NULL;
318  mDaqOccupancyMatch3 = NULL;
319  mDeltaHitMatch3 = NULL;
320  mCellsPrimaryPerEventMatch3 = NULL;
321  hphivsz = NULL;
322  hTofPhivsProj = NULL;
323  hTofZvsProj = NULL;
324  hMtdZvsProj = NULL;
325  hMtdPhivsProj = NULL;
326  hMtddPhivsBackleg = NULL;
327  hMtddZvsBackleg = NULL;
328 }
329 
330 StMtdMatchMaker::~StMtdMatchMaker(){
331  return;
332 }
333 
334 
335 void StMtdMatchMaker::Clear(const char*){
336  StMaker::Clear();
337 }
338 
341 
342  for(int i=0;i<mNAllTrays;i++){
343  for(int j=0;j<mNStrips;j++){
344  mVDrift[i][j] = 60.;
345  }
346  }
347  if(mHisto) bookHistograms();
348 
349  return StMaker::Init();
350 }
351 
353 void StMtdMatchMaker::bookHistograms(){
354 
355  mEventCounterHisto = new TH1D("eventCounter","eventCounter",20,0,20);
356  AddHist(mEventCounterHisto);
357 
358  mCellsMultInEvent = new TH1D("cellsPerEvent","cellsPerEvent",1000,0,1000);
359  AddHist(mCellsMultInEvent);
360 
361  mHitsMultInEvent = new TH1D("hitsPerEvent","hitsPerEvent",1000,0,1000);
362  AddHist(mHitsMultInEvent);
363 
364  mHitsPrimaryInEvent = new TH1D("hitsPrimaryPerEvent","hitsPrimaryPerEvent",1000,0,1000);
365  AddHist(mHitsPrimaryInEvent);
366 
367  mHitsMultPerTrack = new TH1D("hitsPerTrack","hitsPerTrack",10,0,10);
368  AddHist(mHitsMultPerTrack);
369 
370  mHitsPosition = new TH2D("hitsPosition","hitsPositions",1000,-500.,500.,1000,-500.,500);
371  AddHist(mHitsPosition);
372 
373  // histogram not created, memory location invalid
374  // AddHist(mHitsGlobalInEvent);
375 
376  // occupancy
377  for(int i=0;i<mNBacklegs;i++) {
378  char hisname[100];
379  sprintf(hisname,"Occupancy_Backleg_%d",i+1);
380  mDaqOccupancy[i] = new TH1D(hisname,";fired cell(module*12+stripId)",60,0,60);
381  sprintf(hisname,"OccupancyProj_Backleg_%d",i+1);
382  mDaqOccupancyProj[i] = new TH1D(hisname,";projected cell",60,0,60);
383 
384  AddHist(mDaqOccupancy[i]);
385  AddHist(mDaqOccupancyProj[i]);
386  }
387 
388  // correlation
389  for(int i=0;i<mNBacklegs;i++) {
390  char hisname[100];
391  sprintf(hisname,"Corr_Backleg_%d",i+1);
392  mHitCorr[i] = new TH2D(hisname,";project stripId;fired stripId",60,0,60,60,0,60);
393  sprintf(hisname,"Corr_Backleg_%d_module",i+1);
394  mHitCorrModule[i] = new TH2D(hisname,";project moduleId;fired moduleId",6,0,6,6,0,6);
395  AddHist(mHitCorr[i]);
396  AddHist(mHitCorrModule[i]);
397  }
398 
399  // project hit position
400  for(int i=0;i<mNBacklegs;i++) {
401  char hisname[100];
402  sprintf(hisname,"LocalYZ_Backleg_%d",i+1);
403  mDeltaHitFinal[i] = new TH2D(hisname,";localY;localZ",300,-15.,15.,320,-80.,80);
404  AddHist(mDeltaHitFinal[i]);
405  }
406 
407  mTrackPtEta = new TH2D("trackPtEta",";p_{T} (GeV/c);#eta",500,0.,10.,60,-1.5,1.5);
408  mTrackPtPhi = new TH2D("trackPtPhi",";p_{T} (GeV/c);#phi",500,0.,10.,120,0.,2.*M_PI);
409  mTrackNFitPts = new TH1D("trackNFitPts",";nHitsFit",50,0.,50.);
410  mTrackdEdxvsp = new TH2D("trackdEdxvsp",";p_{T} (GeV/c);dE/dx",500,0.,10.,1000,0.,10.);
411  mNSigmaPivsPt = new TH2D("nSigmaPivsPt",";p_{T} (GeV/c);n#sigma_{#pi}",500,0.,10.,1000,-10.,10.);
412 
414  mCellsPerEventMatch1 = new TH1D("cellsPerEventMatch1","cellPerEventMatch1;# hits matched with track(s) per event",10,0,10);
415  mHitsPerEventMatch1 = new TH1D("hitsPerEventMatch1","hitsPerEventMatch1;# tracks matched with hits per event",10,0,10);
416  mCellsPerTrackMatch1 = new TH1D("cellsPerTrackMatch1","cellsPerTrackMatch1;# tracks matched with same hit",10,0,10);
417  mTracksPerCellMatch1 = new TH1D("tracksPerCellMatch1","tracksPerCellMatch1;# hits matched with same track",10,0,10);
418  mDaqOccupancyMatch1 = new TH1D("daqPerCellMatch1","daqPerCellMatch1;stripId",60,0,60);
419  mDeltaHitMatch1 = new TH2D("deltaHitMatch1","deltaHitMatch1;localY;localZ",300,-15,15,800,-80.,80);
420 
422  mCellsPerEventMatch2 = new TH1D("cellsPerEventMatch2","cellPerEventMatch2;# hits matched with track(s) per event",10,0,10);
423  mHitsPerEventMatch2 = new TH1D("hitsPerEventMatch2","hitsPerEventMatch2;# tracks matched with hits per event",10,0,10);
424  mCellsPerTrackMatch2 = new TH1D("cellsPerTrackMatch2","cellsPerTrackMatch2;# tracks matched with same hit",10,0,10);
425  mTracksPerCellMatch2 = new TH1D("tracksPerCellMatch2","tracksPerCellMatch2;# hits matched with same track",10,0,10);
426  mDaqOccupancyMatch2 = new TH1D("daqPerCellMatch2","daqPerCellMatch2;stripId",60,0,60);
427  mDeltaHitMatch2 = new TH2D("deltaHitMatch2","deltaHitMatch2;localY;localZ",300,-15,15,800,-80.,80);
428 
430  mCellsPerEventMatch3 = new TH1D("cellsPerEventMatch3","cellPerEventMatch3;# hits matched with track(s) per event",10,0,10);
431  mHitsPerEventMatch3 = new TH1D("hitsPerEventMatch3","hitsPerEventMatch3;# tracks matched with hits per event",10,0,10);
432  mCellsPerTrackMatch3 = new TH1D("cellsPerTrackMatch3","cellsPerTrackMatch3;# tracks matched with same hit",10,0,10);
433  mTracksPerCellMatch3 = new TH1D("tracksPerCellMatch3","tracksPerCellMatch3;# hits matched with same track",10,0,10);
434  mDaqOccupancyMatch3 = new TH1D("daqPerCellMatch3","daqPerCellMatch3;stripId",60,0,60);
435  mDeltaHitMatch3 = new TH2D("deltaHitMatch3","deltaHitMatch3;localY;localZ",300,-15,15,800,-80.,80);
436 
437  mCellsPrimaryPerEventMatch3 = new TH1D("cellsPrimaryPerEventMatch3","cellsPrimaryPerEventMatch3",10,0,10);
438 
439  hphivsz =new TH2F("hphivsz","hphivsz",500,0,TMath::Pi()*2,500,-500,500);
440  hTofPhivsProj=new TH2F("hTofPhivsProj","hTofPhivsProj",100,0,TMath::Pi()*2,100,0,TMath::Pi()*2);
441  hTofZvsProj=new TH2F("hTofZvsProj","hTofzvsProj",600,-300,300,600,-300,300);
442  hMtdPhivsProj=new TH2F("hMtdPhivsProj","hMtdPhivsProj;projected #phi; fired #phi",360,-TMath::Pi(),TMath::Pi(),180,0,2.*TMath::Pi());
443  hMtddPhivsBackleg=new TH2F("hMtddPhivsBackleg","hMtddPhivsBackleg;backleg; d#phi",30,0,30,1000,-2.*TMath::Pi(),2.*TMath::Pi());
444  hMtddZvsBackleg =new TH2F("hMtddZvsBackleg","hMtddZvsBackleg;backleg; dz",30,0,30,400,-200,200);
445  hMtdZvsProj=new TH2F("hMtdZvsProj","hMtdzvsProj;projected z(cm); fired z(cm)",600,-300,300,600,-300,300);
446 
447  AddHist(mTrackPtEta);
448  AddHist(mTrackPtPhi);
449  AddHist(mTrackNFitPts);
450  AddHist(mTrackdEdxvsp);
451  AddHist(mNSigmaPivsPt);
452 
453  AddHist(mCellsPerEventMatch1);
454  AddHist(mHitsPerEventMatch1);
455  AddHist(mCellsPerTrackMatch1);
456  AddHist(mTracksPerCellMatch1);
457  AddHist(mDaqOccupancyMatch1);
458  AddHist(mDeltaHitMatch1);
459 
460  AddHist(mCellsPerEventMatch2);
461  AddHist(mHitsPerEventMatch2);
462  AddHist(mCellsPerTrackMatch2);
463  AddHist(mTracksPerCellMatch2);
464  AddHist(mDaqOccupancyMatch2);
465  AddHist(mDeltaHitMatch2);
466 
467  AddHist(mCellsPerEventMatch3);
468  AddHist(mHitsPerEventMatch3);
469  AddHist(mCellsPerTrackMatch3);
470  AddHist(mTracksPerCellMatch3);
471  AddHist(mDaqOccupancyMatch3);
472  AddHist(mDeltaHitMatch3);
473 
474  AddHist(mCellsPrimaryPerEventMatch3);
475 
476  AddHist(hphivsz);
477  AddHist(hTofPhivsProj);
478  AddHist(hTofZvsProj);
479  AddHist(hMtdZvsProj);
480  AddHist(hMtdPhivsProj);
481  AddHist(hMtddPhivsBackleg);
482  AddHist(hMtddZvsBackleg);
483 
484  return;
485 }
486 
488 Int_t StMtdMatchMaker::InitRun(int runnumber) {
489 
490  // Get run year
491  mYear= (Int_t)(runnumber/1000000) + 1999;
492  LOG_INFO << "Run year = " << mYear << endm;
493 
494  //=======================================================//
495  // MTD Geometry initialization
496  //=======================================================//
497  LOG_INFO << "Initializing MTD Geometry:" << endm;
498  if(gGeoManager)
499  {
500  LOG_INFO << "TGeoManager (" << gGeoManager->GetName() << "," << gGeoManager->GetTitle() << ") exists" << endm;
501  }
502  else
503  {
504  GetDataBase("VmcGeometry");
505  if(!gGeoManager)
506  {
507  int year = GetDateTime().GetYear();
508  if(mGeomTag.Length()==0)
509  {
510  mGeomTag = Form("y%da",year);
511  if(year<2012) mGeomTag = "y2014a";
512  LOG_INFO << "Load default geometry " << mGeomTag.Data() << " for year " << year << endm;
513  }
514  else
515  {
516  LOG_INFO << "Load input geometry " << mGeomTag.Data() << " for year " << year << endm;
517  }
518 
519  TString ts = Form("$STAR/StarVMC/Geometry/macros/loadStarGeometry.C(\"%s\",1)",mGeomTag.Data());
520  int ierr=0;
521  gROOT->Macro(ts.Data(),&ierr);
522  assert(!ierr);
523  }
524  }
525  assert(gGeoManager);
526 
527  if(!gMtdGeometry)
528  {
529  mMtdGeom = new StMtdGeometry("mtdGeometry","mtdGeometry in StMtdMatchMaker",gGeoManager);
530  mMtdGeom->SetCosmicFlag(mCosmicFlag);
531  mMtdGeom->SetELossFlag(mELossFlag);
532  mMtdGeom->SetNExtraCells(mNExtraCells);
533  if(mLockBField) mMtdGeom->SetLockBField(mLockBField);
534  mMtdGeom->Init(this);
535  LOG_INFO<<" Created a new mtdGeometry ..."<<endm;
536  }
537  assert(gMtdGeometry);
538  mMtdGeom = gMtdGeometry;
539 
540  //========== Set Beam Line =====================
541  double x0 = 0.;
542  double y0 = 0.;
543  double dxdz = 0.;
544  double dydz = 0.;
545 
546  //Get Current Beam Line Constraint from database
547  TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic");
548 
549  if (dbDataSet) {
550  vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
551 
552  x0 = vSeed->x0;
553  y0 = vSeed->y0;
554  dxdz = vSeed->dxdz;
555  dydz = vSeed->dydz;
556  }
557  else {
558  LOG_INFO << "StMtdMatchMaker -- No Database for beamline" << endm;
559  }
560 
561  LOG_INFO << "BeamLine Constraint (StMtdMatchMaker): " << endm;
562  LOG_INFO << "x(z) = " << x0 << " + " << dxdz << " * z" << endm;
563  LOG_INFO << "y(z) = " << y0 << " + " << dydz << " * z" << endm;
564 
565  StThreeVectorD origin(x0,y0,0.0);
566  double pt = 88889999;
567  double nxy=::sqrt(dxdz*dxdz + dydz*dydz);
568  if(nxy<1.e-5){ // beam line _MUST_ be tilted
569  LOG_WARN << "StMtdMatchMaker:: Beam line must be tilted!" << endm;
570  nxy=dxdz=1.e-5;
571  }
572  double p0=pt/nxy;
573  double px = p0*dxdz;
574  double py = p0*dydz;
575  double pz = p0; // approximation: nx,ny<<0
576  StThreeVectorD MomFstPt(px*GeV, py*GeV, pz*GeV);
577  mBeamHelix = new StPhysicalHelixD(MomFstPt,origin,0.5*tesla,1.);
578  return kStOK;
579 }
580 
581 
583 Int_t StMtdMatchMaker::FinishRun(int runnumber)
584 {
585 
586  delete mBeamHelix;
587 
588  return kStOK;
589 }
590 
591 
594 
595  return kStOK;
596 }
597 
598 
601 
602  mEvent=(StEvent *) GetInputDS("StEvent");
603  if(mEvent){
604  mMuDstIn = kFALSE;
605  }else{
606  mMuDst = (StMuDst *)GetInputDS("MuDst");
607  mMuDstIn = kTRUE;
608  }
609 
610  if(mHisto) mEventCounterHisto->Fill(0);
611  StTimer timer;
612  if(doPrintCpuInfo) timer.start();
613  if(doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
614 
615  // clean up mtdPidTraits in MuDst
616  if(mMuDstIn) cleanUpMtdPidTraits();
617 
618 
619  // read data from StMtdHit
621  //
622  mtdCellHitVector daqCellsHitVec;
623  idVector validModuleVec;
624  if(!readMtdHits(daqCellsHitVec,validModuleVec)) return kStOK;
625  //end of Sect.A
626  if(Debug()){
627  LOG_INFO<<" Sect.A: =============================="<<endm;
628  LOG_INFO <<" total # of cells =" << daqCellsHitVec.size() << endm;
629  for(size_t iv=0;iv<validModuleVec.size();iv++){
630  LOG_INFO << " module# "<< validModuleVec[iv] <<" Valid! "<<endm;
631  }
632  LOG_INFO<<" daqCellsHitVec:"<<endm;
633  mtdCellHitVectorIter daqIter = daqCellsHitVec.begin();
634  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
635  int daqChn = (daqIter->module-1)*12 + (daqIter->cell);
636  int daqBackleg = daqIter->backleg;
637  LOG_INFO<<" daq backleg:"<<daqBackleg<<" chn:"<<daqChn<<endm;
638  }
639  }
640 
641  if(mHisto) {
642  mCellsMultInEvent->Fill(daqCellsHitVec.size());
643  if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
644  }
645  if(doPrintCpuInfo) {
646  timer.stop();
647  LOG_INFO <<"CPU time after Step A - loading hits: "
648  << timer.elapsedTime() <<" sec"<<endm;
649  timer.start();
650  }
651  //...................................................................................
652  //
653  mtdCellHitVector allCellsHitVec;
654  Int_t nPrimaryHits(0);
656  project2Mtd(daqCellsHitVec,allCellsHitVec,nPrimaryHits);
657 
658  if(mHisto) {
659  mHitsMultInEvent->Fill(allCellsHitVec.size());
660  mHitsPrimaryInEvent->Fill(nPrimaryHits);
661  if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
662  }
663 
664  if(Debug()){
665  LOG_INFO<<" Sect.B: =============================="<<endm;
666  LOG_INFO<<" allCellsHitVec:"<<endm;
667  mtdCellHitVectorIter proIter = allCellsHitVec.begin();
668  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
669  int proChn = (proIter->module-1)*12 + (proIter->cell);
670  int proBackleg = proIter->backleg;
671  LOG_INFO<<" proj backleg:"<<proBackleg<<" chn:"<<proChn<<endm;
672  }
673  }
674  // end of Sect.B
675  if(doPrintCpuInfo) {
676  timer.stop();
677  LOG_INFO << "CPU time after Step B - projection:"
678  << timer.elapsedTime() << "sec" <<endm;
679  timer.start();
680  }
681  //....................................................................................
682  //
684  mtdCellHitVector matchHitCellsVec;
685  matchMtdHits(daqCellsHitVec,allCellsHitVec,matchHitCellsVec);
686 
687  if(Debug()){
688  LOG_INFO<<" Sect.C: =============================="<<endm;
689  LOG_INFO<<" matchCellsHitVec:"<<endm;
690  mtdCellHitVectorIter matIter = matchHitCellsVec.begin();
691  for(unsigned int imat=0;imat<matchHitCellsVec.size();imat++, matIter++) {
692  int matChn = (matIter->module-1)*12 + (matIter->cell);
693  int matBackleg = matIter->backleg;
694  LOG_INFO<<" match backleg:"<<matBackleg<<" chn:"<<matChn<<endm;
695  }
696  }
697 
698  //end of Sect.C
699  if(Debug()){ LOG_INFO <<"C:before/after:"<< allCellsHitVec.size()<< "/"<<matchHitCellsVec.size()<<endm;}
700  if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
701  if(doPrintCpuInfo){
702  timer.stop();
703  LOG_INFO << "CPU time after Step C - matching:"
704  <<timer.elapsedTime()<< "sec" <<endm;
705  timer.start();
706  }
707  //....................................................................................
709  //
710  mtdCellHitVector singleHitCellsVec;
711  mtdCellHitVector multiHitsCellsVec;
712  sortSingleAndMultiHits(matchHitCellsVec,singleHitCellsVec,multiHitsCellsVec);
713 
714  if(Debug()){
715  LOG_INFO<<" Sect.D: =============================="<<endm;
716  LOG_INFO<<" singleCellsHitVec:"<<endm;
717  mtdCellHitVectorIter singleIter = singleHitCellsVec.begin();
718  for(unsigned int isingle=0;isingle<singleHitCellsVec.size();isingle++, singleIter++) {
719  int singleChn = (singleIter->module-1)*12 + (singleIter->cell);
720  int singleBackleg = singleIter->backleg;
721  LOG_INFO<<" single backleg:"<<singleBackleg<<" chn:"<<singleChn<<endm;
722  LOG_INFO<<" LeTimeWest:"<<singleIter->leadingEdgeTime.first<<" LeTimeEast:"<<singleIter->leadingEdgeTime.second<<endm;
723  }
724  LOG_INFO<<" multiCellsHitVec:"<<endm;
725  mtdCellHitVectorIter multiIter = multiHitsCellsVec.begin();
726  for(unsigned int imulti=0;imulti<multiHitsCellsVec.size();imulti++, multiIter++) {
727  int multiChn = (multiIter->module-1)*12 + (multiIter->cell);
728  int multiBackleg = multiIter->backleg;
729  LOG_INFO<<" multi backleg:"<<multiBackleg<<" chn:"<<multiChn<<endm;
730  }
731  }
732  //end of Sect.D
733  if(doPrintCpuInfo){
734  timer.stop();
735  LOG_INFO <<" CPU time after Step D - sorting:"
736  <<timer.elapsedTime()<< " sec"<<endm;
737  timer.start();
738  }
739  //......................................................................................
741  //
742  mtdCellHitVector finalMatchedCellsVec;
743  finalMatchedMtdHits(singleHitCellsVec,finalMatchedCellsVec);
744 
745  if(Debug()){
746  LOG_INFO<<" Sect.E: =============================="<<endm;
747  LOG_INFO<<" finalMatchedCellsHitVec:"<<endm;
748  mtdCellHitVectorIter finalIter = finalMatchedCellsVec.begin();
749  for(unsigned int ifinal=0;ifinal<finalMatchedCellsVec.size();ifinal++, finalIter++) {
750  int finalChn = (finalIter->module-1)*12 + (finalIter->cell);
751  int finalBackleg = finalIter->backleg;
752  LOG_INFO<<" final backleg:"<<finalBackleg<<" chn:"<<finalChn<<endm;
753  }
754  }
755  //end of Sect.E
756  if(Debug()){ LOG_INFO <<"E: before/after:"<< singleHitCellsVec.size() <<"/" <<finalMatchedCellsVec.size()<<endm;}
757  if(doPrintCpuInfo){
758  LOG_INFO <<"CPU time after Step E - final matched:"
759  <<timer.elapsedTime() <<" sec"<<endm;
760  timer.start();
761  }
762  //.......................................................................
764  //
765  if(mHisto) {
766  if(finalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
767  mCellsPerEventMatch3->Fill(finalMatchedCellsVec.size());
768  }
769  Int_t nValidSingleHitCells(0),nValidSinglePrimHitCells(0);
770  fillPidTraits(finalMatchedCellsVec,nValidSingleHitCells,nValidSinglePrimHitCells);
771  //end of Sect.F
772  if(Debug()){
773  LOG_INFO<<" Sect.F: =============================="<<endm;
774  LOG_INFO<<" nValidSingleHitCells:"<<nValidSingleHitCells<<endm;
775 
776  LOG_INFO <<"#(daq):" <<daqCellsHitVec.size()
777  << "#(proj):"<<allCellsHitVec.size()
778  << "#(prim proj):" <<nPrimaryHits
779  << "#(single valid):"<<nValidSingleHitCells <<endm;
780  LOG_INFO <<"F: before/after"<< finalMatchedCellsVec.size()<< "/"<<nValidSingleHitCells <<endm;
781  }
782  if(doPrintCpuInfo){
783  timer.stop();
784  LOG_INFO << "CPU time after Step F - final :"
785  << timer.elapsedTime()<< "sec"<<endm;
786  timer.start();
787  }
789  if(mMuDstIn){
790  StMuMtdCollection *theMtd = mMuDst->MtdCollection();
791  if(theMtd){
792  if(theMtd->hitsPresent()){
793  LOG_DEBUG <<" MtdCollection: hit container present. "<<endm;
794  if(Debug()){
795  LOG_INFO << "# of hits in this event:" << theMtd->hitsPresent() <<endm;
796  for(Int_t i=0;i<theMtd->hitsPresent();i++){
797  StMuMtdHit* aHit = theMtd->MtdHit(i) ;
798  if(!aHit) continue;
799  LOG_INFO <<"backleg:"<<aHit->backleg()
800  <<" module:"<<aHit->backleg()
801  <<" cell:"<<aHit->cell()
802  <<" tof:"<<aHit->tof()<<endm;
803  }
804  }
805  }
806  }
807  }else{
808  StMtdCollection *theMtd = mEvent->mtdCollection();
809  if(theMtd){
810  if(theMtd->hitsPresent()){
811  LOG_DEBUG <<" MtdCollection: hit container present. "<<endm;
812  if(Debug()){
813  StSPtrVecMtdHit& tmpCellMtdVec = theMtd->mtdHits();
814  LOG_INFO << "# of hits in this event:" << tmpCellMtdVec.size()<<endm;
815  for(size_t i=0;i<tmpCellMtdVec.size();i++){
816  StMtdHit* p=tmpCellMtdVec[i];
817  LOG_INFO <<(*p)<<endm;
818  }
819  }
820  }
821  }
822  }
823  //-- end check
824  if(doPrintMemoryInfo){
825  StMemoryInfo::instance()->snapshot();
826  StMemoryInfo::instance()->print();
827  }
828  if(doPrintCpuInfo){
829  timer.stop();
830  LOG_INFO<< "CPU time for StMtdMatchMaker::Make():"
831  <<timer.elapsedTime()<< "sec"<<endm;
832  }
833  if(Debug()){
834  LOG_INFO<<"StMtdMatchMaker -- bye--bye"<<endm;
835  }
836 
837  return kStOK;
838 }
839 
840 //---------------------------------------------------------------------------
842 void StMtdMatchMaker::cleanUpMtdPidTraits()
843 {
844  index2Primary.clear();
845  for(Int_t ii=0;ii<(Int_t)mMuDst->array(muPrimary)->GetEntries();ii++)
846  {
847  StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(ii);
848  if(!pTrack) continue;
849  Int_t index2Global = pTrack->index2Global();
850  if(index2Global<0) continue;
851  index2Primary[index2Global] = ii;
852  }
853 
854  UInt_t Nnodes = mMuDst->numberOfGlobalTracks();
855  for(UInt_t iNode=0;iNode<Nnodes;iNode++)
856  {
857  StMuTrack *theTrack = mMuDst->globalTracks(iNode);
858  if(!theTrack) continue;
859  if(theTrack->index2MtdHit()<0) continue;
860 
861  //clean up any association done before
862  StMuMtdPidTraits pidMtd;
863  theTrack->setMtdPidTraits(pidMtd);
864  theTrack->setIndex2MtdHit(-999);
865 
866  Int_t pIndex = -999;
867  map<Int_t, Int_t>::iterator it = index2Primary.find(iNode);
868  if(it!=index2Primary.end()){
869  pIndex = it->second;
870  }
871  if(pIndex>=0){
872  StMuTrack *thePrimaryTrack= (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(pIndex);
873  if(thePrimaryTrack){
874  thePrimaryTrack->setMtdPidTraits(pidMtd);
875  thePrimaryTrack->setIndex2MtdHit(-999);
876  }
877  }
878  }
879 }
880 
881 //---------------------------------------------------------------------------
883 Bool_t StMtdMatchMaker::readMtdHits(mtdCellHitVector& daqCellsHitVec,idVector& validModuleVec){
884 
885  // StMuMtdCollection can't save Hits yet. Pending MuDst Mode.
886  if(mMuDstIn){
887  if(!mMuDst){
888  LOG_INFO << "No Mudst ... bye-bye" <<endm;
889  return kFALSE;
890  }
892  int nMtdHits = mMuDst->numberOfMTDHit();
893  if(nMtdHits<=0) return kFALSE;
894  for(Int_t i=0;i<nMtdHits;i++){
895  StMuMtdHit* aHit = mMuDst->mtdHit(i) ;
896  if(!aHit) continue;
897  if(aHit->backleg()<=0||aHit->backleg()>mNBacklegs) continue; // barrel BackLeg hits
898 
899  //clean up any association done before
900  aHit->setIndex2Primary(-1);
901  aHit->setIndex2Global(-1);
902  aHit->setAssociatedTrackKey(-1);
903  int backlegId = aHit->backleg();
904  int moduleId = aHit->module();
905  int cellId = aHit->cell();
906  if(Debug()) {LOG_INFO <<"A: fired hit in " << "backleg=" << backlegId <<" module="<<moduleId<<" cell="<<cellId<<endm;}
907  StructCellHit aDaqCellHit;
908  aDaqCellHit.backleg = backlegId;
909  aDaqCellHit.module= moduleId;
910  aDaqCellHit.cell=cellId;
911  aDaqCellHit.tot=aHit->tot();
912  aDaqCellHit.leadingEdgeTime=aHit->leadingEdgeTime();
913  aDaqCellHit.index2MtdHit=i;
914  aDaqCellHit.idTruth = aHit->idTruth();
915  daqCellsHitVec.push_back(aDaqCellHit);
916 
917  //additional valid number configuration
918  int id=backlegId*100+moduleId;
919  if(find(validModuleVec.begin(),validModuleVec.end(),id) == validModuleVec.end())
920  validModuleVec.push_back(id);
921  int hisIndex = backlegId - 1;
922  if(mHisto) {
923  mDaqOccupancy[hisIndex]->Fill((moduleId-1)*12+cellId);
924  }
925  }
926  }else{
927  if(!mEvent||!(mEvent->mtdCollection())||!(mEvent->mtdCollection()->hitsPresent())){
928  if(!mEvent){LOG_INFO << "no StEvent" <<endm; return kFALSE;}
929  else if(!(mEvent->mtdCollection())) {
930  LOG_INFO << "no MTD Collection" <<endm;
931  return kFALSE;
932  }
933  else if(!(mEvent->mtdCollection()->hitsPresent())){LOG_INFO << "no MTD hits present" <<endm; return kFALSE;}
934  return kFALSE;
935  }
936 
937 
939  //.........................................................................
941  StMtdCollection *theMtd = mEvent->mtdCollection();
942 
943  //.........................................................................
945  StSPtrVecMtdHit& mtdHits= theMtd->mtdHits();
946 
947  for(size_t i=0;i<mtdHits.size();i++){
948  StMtdHit* aHit = mtdHits[i];
949  if(!aHit) continue;
950  if(aHit->backleg()<=0||aHit->backleg()>mNBacklegs) continue; // barrel BackLeg hits
951 
952  //clean up any association done before
953  aHit->setAssociatedTrack(0);
954 
955  int backlegId = aHit->backleg();
956  int moduleId = aHit->module();
957  int cellId = aHit->cell();
958 
959  if(Debug()) {LOG_INFO <<"A: fired hit in " << "backleg=" << backlegId <<" module="<<moduleId<<" cell="<<cellId<<endm;}
960  LOG_DEBUG <<"A: fired hit in " << "backleg=" << backlegId <<" module="<<moduleId<<" cell="<<cellId<<" leadingWest="<<aHit->leadingEdgeTime().first<<" leadingEast="<<aHit->leadingEdgeTime().second<<endm;
961  StructCellHit aDaqCellHit;
962  aDaqCellHit.backleg = backlegId;
963  aDaqCellHit.module= moduleId;
964  aDaqCellHit.cell=cellId;
965  aDaqCellHit.tot=aHit->tot();
966  aDaqCellHit.leadingEdgeTime=aHit->leadingEdgeTime();
967  aDaqCellHit.index2MtdHit=i;
968  aDaqCellHit.idTruth = aHit->idTruth();
969  daqCellsHitVec.push_back(aDaqCellHit);
970 
971  //additional valid number configuration
972  int id=backlegId*100+moduleId;
973  if(find(validModuleVec.begin(),validModuleVec.end(),id) == validModuleVec.end())
974  validModuleVec.push_back(id);
975 
976  int hisIndex = backlegId - 1;
977  if(mHisto) {
978  mDaqOccupancy[hisIndex]->Fill((moduleId-1)*12+cellId);
979  }
980  }
981  }
982  return kTRUE;
983 }
984 
985 
987 void StMtdMatchMaker::project2Mtd(mtdCellHitVector daqCellsHitVec,mtdCellHitVector& allCellsHitVec,Int_t& nPrimaryHits){
988  int nAllTracks=0;
989  ngTracks = 0;
990 
991  StThreeVectorD pVtx(0,0,0);
992  UInt_t Nnodes = 0;
993  if(mMuDstIn){
994  Nnodes = mMuDst->numberOfGlobalTracks();
995  for(UInt_t iNode=0;iNode<Nnodes;iNode++){
996  StMuTrack *theTrack=mMuDst->globalTracks(iNode);
997  if(!theTrack) continue;
998  if(!validTrack(theTrack)) continue;
999 
1000  bool isPrimary=kFALSE;
1001  Int_t pIndex = -999;
1002  map<Int_t, Int_t>::iterator it = index2Primary.find(iNode);
1003  if(it!=index2Primary.end())
1004  pIndex = it->second;
1005  if(pIndex>=0) isPrimary=kTRUE;
1006 
1007  pVtx = mMuDst->event()->primaryVertexPosition();
1008  if(isPrimary)
1009  {
1010  int vtxIndex = ((StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(pIndex))->vertexIndex();
1011  if(vtxIndex>-1)
1012  {
1013  StMuPrimaryVertex *vertex = mMuDst->primaryVertex(vtxIndex);
1014  if(vertex) pVtx = vertex->position();
1015  }
1016  }
1017  if(matchTrack2Mtd(daqCellsHitVec,theTrack->outerHelix(),theTrack->charge(),allCellsHitVec,iNode,pVtx)){
1018  nAllTracks++;
1019  if(isPrimary) nPrimaryHits++;
1020  }
1021  ngTracks++;
1022  }
1023  }else{
1024  Nnodes = mEvent->trackNodes().size();
1025  for(UInt_t iNode=0;iNode<Nnodes;iNode++){
1026  StSPtrVecTrackNode& nodes=mEvent->trackNodes();
1027  StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
1028  if(!theTrack) continue;
1029 
1030  //clean up any association done before
1031  StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
1032  for (StSPtrVecTrackPidTraitsIterator it = traits.begin(); it != traits.end(); it++)
1033  {
1034  if( (*it)->detector() == kMtdId )
1035  {
1036  traits.erase(it);
1037  break;
1038  }
1039  }
1040 
1041  if (mEvent->primaryVertex()) pVtx = mEvent->primaryVertex()->position();
1042  else pVtx.set(0,0,0);
1043  bool isPrimary =kFALSE;
1044  StPrimaryTrack *pTrack =dynamic_cast<StPrimaryTrack*>(theTrack->node()->track(primary));
1045  if(pTrack)
1046  {
1047  isPrimary = kTRUE;
1048  if(pTrack->vertex()) pVtx = pTrack->vertex()->position();
1049 
1050  //clean up any association done before
1051  StSPtrVecTrackPidTraits& ptraits = pTrack->pidTraits();
1052  for (StSPtrVecTrackPidTraitsIterator it = ptraits.begin(); it != ptraits.end(); it++)
1053  {
1054  if( (*it)->detector() == kMtdId )
1055  {
1056  ptraits.erase(it);
1057  break;
1058  }
1059  }
1060  }
1061  if(!validTrack(theTrack)) continue;
1062  if(matchTrack2Mtd(daqCellsHitVec,theTrack->outerGeometry()->helix(),theTrack->geometry()->charge(),allCellsHitVec,iNode,pVtx)){
1063  nAllTracks++;
1064  if(isPrimary) nPrimaryHits++;
1065  }
1066  ngTracks++;
1067  } // end for
1068  } // end if (mMuDstIn)
1069 
1070 }
1071 
1072 
1074 bool StMtdMatchMaker::matchTrack2Mtd(mtdCellHitVector daqCellsHitVec,const StPhysicalHelixD &helix, Int_t gq, mtdCellHitVector& allCellsHitVec, unsigned int iNode, StThreeVectorD pVtx){
1075  float mField = 0;
1076  if(mMuDstIn) mField = mMuDst->event()->runInfo().magneticField();
1077  else mField = mEvent->runInfo()->magneticField();
1078 
1079  StThreeVector<double> dcaPos = helix.at(helix.pathLength(pVtx));
1080  StThreeVector<double> dca = dcaPos - pVtx;
1081  //LOG_INFO<<"******************* B Field from data = "<<mField/10.<<" charge = "<<gq<<" **********************"<<endm;
1082 
1083  float mFieldFromGeom = mMtdGeom->GetFieldZ(100,100,0.);
1084  if(fabs(mFieldFromGeom-mField/10.)>0.2){
1085  LOG_WARN<<"Wrong magnetc field mField = "<<mField/10.<<" mFieldFromGeom = "<<mFieldFromGeom<<" check the magF input!!!"<<endm;
1086  }
1087  IntVec idVec;
1088  DoubleVec pathVec;
1089  DoubleVec tofVec;
1090  PointVec crossVec;
1091  if(mCosmicFlag){
1092  mMtdGeom->HelixCrossCellIds(helix,idVec,pathVec,crossVec,tofVec);
1093  }else{
1094  mMtdGeom->HelixCrossCellIds(helix,pVtx,idVec,pathVec,crossVec,tofVec);
1095  }
1096 
1097  if((idVec.size()!=pathVec.size()) || idVec.size()!=crossVec.size() || idVec.size()!=tofVec.size()){
1098  if(Debug()){
1099  LOG_INFO<<"Inconsistent size idVec = "<<idVec.size()<<" pathVec = "<<pathVec.size()<<" crossVec = "<<crossVec.size()<<" tofVec = "<<tofVec.size()<<endm;
1100  }
1101  return kFALSE;
1102  }
1103 
1104  int nCells = 0;
1105  if(Debug()){
1106  LOG_INFO<<"dca:"<<dca.mag()<<endm;
1107  LOG_INFO<<"idVec:"<<idVec.size()<<endm;
1108  LOG_INFO<<"iNode:"<<iNode<<endm;
1109  if(idVec.size()>1) LOG_INFO<<"MARK: one track cross two modules!"<<endm;
1110  }
1111  for (UInt_t i = 0; i < idVec.size(); i++) {
1112  Int_t iBL = -9999;
1113  Int_t iMod = -9999;
1114  Int_t iCell = -9999;
1115  mMtdGeom->DecodeCellId(idVec[i],iBL,iMod,iCell);
1116  StMtdGeoModule *mMtdGeoModule = mMtdGeom->GetGeoModule(iBL,iMod);
1117  Double_t local[3]={0,0,0};
1118  Double_t global[3]={crossVec[i].x(),crossVec[i].y(),crossVec[i].z()};
1119  if(mMtdGeoModule) mMtdGeoModule->MasterToLocal(global,local);
1120  StThreeVectorD ol(local[0],local[1],local[2]);
1121 
1122  StructCellHit cellHit;
1123  cellHit.backleg = iBL;
1124  cellHit.module = iMod;
1125  cellHit.cell = iCell;
1126  cellHit.trackIdVec.push_back((Int_t)iNode);
1127  cellHit.hitPosition = crossVec[i];
1128  cellHit.zhit = ol.z();
1129  cellHit.yhit = ol.y();
1130  cellHit.theta = ol.theta();
1131  cellHit.pathLength = pathVec[i];
1132  cellHit.expTof2MTD = tofVec[i];
1133  allCellsHitVec.push_back(cellHit);
1134  nCells++;
1135 
1136  if(Debug()){
1137  if(idVec.size()>1) LOG_INFO<<"iBL:iMod:iCell="<<iBL<<" "<<iMod<<" "<<iCell<<endm;
1138  }
1139  }
1140  return kTRUE;
1141 }
1142 
1144 void StMtdMatchMaker::matchMtdHits(mtdCellHitVector& daqCellsHitVec,mtdCellHitVector& allCellsHitVec,mtdCellHitVector& matchHitCellsVec){
1145  StructCellHit cellHit;
1146  mtdCellHitVectorIter daqIter = daqCellsHitVec.begin();
1147 
1148  for(unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
1149  mtdCellHitVectorIter proIter = allCellsHitVec.begin();
1150  StMtdGeoModule *geoModule = mMtdGeom->GetGeoModule(daqIter->backleg,daqIter->module);
1151  if(daqIter->backleg==9)
1152  {
1153  if(!geoModule)
1154  {
1155  LOG_WARN << "Geometry module not available for BL = " << daqIter->backleg << " Mod = " << daqIter->module << endm;
1156  continue;
1157  }
1158  }
1159  else
1160  assert(geoModule);
1161  for(unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
1162 
1163  int daqIndex = (daqIter->module-1)*12 + (daqIter->cell);
1164  int proIndex = (proIter->module-1)*12 + (proIter->cell);
1165  int hisIndex = daqIter->backleg - 1;
1166  if (mHisto) {
1167 
1168  double stripPhiCen = 0.;
1169  int channel = daqIter->cell;
1170  stripPhiCen = geoModule->GetCellPhiCenter(channel);
1171 
1172  double mLeTimeWest = daqIter->leadingEdgeTime.first;
1173  double mLeTimeEast = daqIter->leadingEdgeTime.second;
1174  StThreeVectorD modCen = geoModule->GetNodePoint();
1175  double stripZCen = modCen.z() - (mLeTimeWest-mLeTimeEast)/2./gMtdCellDriftV*1000.;
1176 
1177  double daqphi = stripPhiCen;
1178  double daqz = stripZCen;
1179 
1180  hMtdZvsProj->Fill(proIter->hitPosition.z(),daqz);
1181  hMtdPhivsProj->Fill(proIter->hitPosition.phi(),daqphi);
1182  hMtddPhivsBackleg->Fill(hisIndex,proIter->hitPosition.phi()-daqphi);
1183  hMtddZvsBackleg->Fill(hisIndex,proIter->hitPosition.z()-daqz);
1184 
1185  }
1186 
1187  if(daqIter->backleg==proIter->backleg) {
1188  if (mHisto) {
1189  if(hisIndex>=0&&hisIndex<30) {
1190  mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
1191  mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
1192  } else {
1193  LOG_WARN << " weird tray # " << daqIter->module<< endm;
1194  }
1195  }
1196  }
1197 
1198  /*
1199  int iNode = proIter->trackIdVec[0];
1200  if(mMuDstIn){
1201  StMuTrack *theTrack=mMuDst->globalTracks(iNode);
1202  }else{
1203  StSPtrVecTrackNode& nodes=mEvent->trackNodes();
1204  StGlobalTrack *theTrack = dynamic_cast<StGlobalTrack*>(nodes[iNode]->track(global));
1205  }
1206  */
1207  StThreeVectorD modCen = geoModule->GetNodePoint();
1208 
1209  Int_t ibackleg = daqIter->backleg;
1210  Int_t imodule = daqIter->module;
1211  Int_t icell = daqIter->cell;
1212 
1213  double zdaq = (daqIter->leadingEdgeTime.second-daqIter->leadingEdgeTime.first)/2./mVDrift[(ibackleg-1)*gMtdNModules+imodule-1][icell]*1e3;
1214  bool isMatch = false;
1215 
1216 
1217  if(daqIter->backleg==proIter->backleg){
1218  if(daqIter->module==proIter->module) isMatch = true;
1219  }
1220 
1221  if(mnNeighbors){
1222  if(daqIter->backleg==proIter->backleg && daqIter->module!=proIter->module){
1223  if(zdaq<0){
1224  if((daqIter->module-1>0)&&daqIter->module-1==proIter->module) isMatch = true;
1225  }else if(zdaq==0){
1226  if(abs(daqIter->module-proIter->module)<=1) isMatch = true;
1227  }else{
1228  if((daqIter->module+1<6)&&daqIter->module+1==proIter->module) isMatch = true;
1229  }
1230  }
1231  }
1232 
1233 
1234  if(isMatch){
1235  cellHit.backleg = daqIter->backleg;
1236  cellHit.module = daqIter->module;
1237  cellHit.cell = daqIter->cell;
1238  cellHit.hitPosition = proIter->hitPosition;
1239  cellHit.trackIdVec = proIter->trackIdVec;
1240  cellHit.zhit = proIter->zhit;
1241  cellHit.yhit = proIter->yhit;
1242  cellHit.tot = daqIter->tot;
1243  cellHit.leadingEdgeTime = daqIter->leadingEdgeTime;
1244  cellHit.index2MtdHit = daqIter->index2MtdHit;
1245  cellHit.theta = proIter->theta;
1246  cellHit.pathLength = proIter->pathLength;
1247  cellHit.expTof2MTD = proIter->expTof2MTD;
1248  cellHit.idTruth = daqIter->idTruth;
1249  matchHitCellsVec.push_back(cellHit);
1250  }
1251  }
1252  } //end {sec. C}'
1253 }
1254 
1255 
1256 
1258 void StMtdMatchMaker::sortSingleAndMultiHits(mtdCellHitVector& matchHitCellsVec,mtdCellHitVector& singleHitCellsVec,mtdCellHitVector& multiHitsCellsVec)
1259 {
1260  StructCellHit cellHit;
1261  StructCellHit cellHit_candidate;
1262  mtdCellHitVector tempVec = matchHitCellsVec;
1263  mtdCellHitVector erasedVec = tempVec;
1264  mtdCellHitVector multiHitsCellsVec_temp;
1265  while (tempVec.size() != 0)
1266  {
1267  Int_t nTracks = 0;
1268  idVector trackIdVec;
1269 
1270  mtdCellHitVectorIter tempIter=tempVec.begin();
1271  mtdCellHitVectorIter erasedIter=erasedVec.begin();
1272  while(erasedIter!= erasedVec.end())
1273  {
1274  if(tempIter->backleg == erasedIter->backleg &&
1275  tempIter->module == erasedIter->module &&
1276  tempIter->cell == erasedIter->cell)
1277  {
1278  //
1279  if(nTracks>0)
1280  {
1281  if(nTracks==1)
1282  {
1283  cellHit.cell = tempIter->cell;
1284  cellHit.module = tempIter->module;
1285  cellHit.backleg = tempIter->backleg;
1286  cellHit.hitPosition = tempIter->hitPosition;
1287  cellHit.trackIdVec.push_back(tempIter->trackIdVec.back());
1288  cellHit.zhit = tempIter->zhit;
1289  cellHit.matchFlag = -999;
1290  cellHit.yhit = tempIter->yhit;
1291  cellHit.tot = tempIter->tot;
1292  cellHit.leadingEdgeTime = tempIter->leadingEdgeTime;
1293  cellHit.index2MtdHit = tempIter->index2MtdHit;
1294  cellHit.theta = tempIter->theta;
1295  cellHit.pathLength = tempIter->pathLength;
1296  cellHit.expTof2MTD = tempIter->expTof2MTD;
1297  cellHit.idTruth = tempIter->idTruth;
1298  if(Debug())
1299  {
1300  LOG_INFO<<"track Info: "<<cellHit.zhit<<" "<<cellHit.yhit<<endm;
1301  LOG_INFO<<"hit Info: "<<cellHit.backleg<<" "<<cellHit.module<<" "<<cellHit.cell<<endm;
1302  }
1303  multiHitsCellsVec_temp.push_back(cellHit);
1304  }
1305 
1306  cellHit_candidate.cell = erasedIter->cell;
1307  cellHit_candidate.module = erasedIter->module;
1308  cellHit_candidate.backleg = erasedIter->backleg;
1309  cellHit_candidate.hitPosition = erasedIter->hitPosition;
1310  cellHit_candidate.trackIdVec.push_back(erasedIter->trackIdVec.back());
1311  cellHit_candidate.zhit = erasedIter->zhit;
1312  cellHit_candidate.matchFlag = -999;
1313  cellHit_candidate.yhit = erasedIter->yhit;
1314  cellHit_candidate.tot = erasedIter->tot;
1315  cellHit_candidate.leadingEdgeTime = erasedIter->leadingEdgeTime;
1316  cellHit_candidate.index2MtdHit = erasedIter->index2MtdHit;
1317  cellHit_candidate.theta = erasedIter->theta;
1318  cellHit_candidate.pathLength = erasedIter->pathLength;
1319  cellHit_candidate.expTof2MTD = erasedIter->expTof2MTD;
1320  cellHit_candidate.idTruth = erasedIter->idTruth;
1321  multiHitsCellsVec_temp.push_back(cellHit_candidate);
1322  if(Debug())
1323  {
1324  LOG_INFO<<"track Info: "<<cellHit_candidate.zhit<<" "<<cellHit_candidate.yhit<<endm;
1325  LOG_INFO<<"hit Info: "<<cellHit_candidate.backleg<<" "<<cellHit_candidate.module<<" "<<cellHit_candidate.cell<<endm;
1326  }
1327  }
1328 
1329  nTracks++;
1330  trackIdVec.push_back(erasedIter->trackIdVec.back()); // merge
1331  erasedVec.erase(erasedIter);
1332  --erasedIter;
1333  }
1334  ++erasedIter;
1335  }
1336 
1337  cellHit.cell = tempIter->cell;
1338  cellHit.module = tempIter->module;
1339  cellHit.backleg = tempIter->backleg;
1340  cellHit.hitPosition = tempIter->hitPosition;
1341  cellHit.trackIdVec = trackIdVec;
1342  cellHit.zhit = tempIter->zhit;
1343  cellHit.matchFlag = -999;
1344  cellHit.yhit = tempIter->yhit;
1345  cellHit.tot = tempIter->tot;
1346  cellHit.leadingEdgeTime = tempIter->leadingEdgeTime;
1347  cellHit.index2MtdHit = tempIter->index2MtdHit;
1348  cellHit.theta = tempIter->theta;
1349  cellHit.pathLength = tempIter->pathLength;
1350  cellHit.expTof2MTD = tempIter->expTof2MTD;
1351  cellHit.idTruth = tempIter->idTruth;
1352 
1353  if (nTracks==1)
1354  {
1355  // hit is matched one track
1356  cellHit.matchFlag = 1;
1357  singleHitCellsVec.push_back(cellHit);
1358  }
1359  else if (nTracks>1)
1360  {
1361  // hit is matched multiple tracks
1362  bool isMatchToSingleTrack = kFALSE;
1363  idVector tmpIdVec = trackIdVec;
1364  sort(tmpIdVec.begin(),tmpIdVec.end());
1365  tmpIdVec.erase(unique(tmpIdVec.begin(),tmpIdVec.end()),tmpIdVec.end());
1366  if(tmpIdVec.size()==1)
1367  {
1368  isMatchToSingleTrack = kTRUE;
1369  if(Debug())
1370  {
1371  LOG_INFO<<"A hit is matched to one track that is projected to two modules!"<<endm;
1372  }
1373  }
1374  for(int iTrack = 0; iTrack < nTracks; iTrack++)
1375  {
1376  if(isMatchToSingleTrack) multiHitsCellsVec_temp[iTrack].matchFlag = 1;
1377  multiHitsCellsVec.push_back(multiHitsCellsVec_temp[iTrack]);
1378  }
1379  }
1380  else
1381  {
1382  LOG_WARN << "D: no tracks extrapolate to matched cell ... should not happen!" << endm;
1383  }
1384 
1385  if(Debug())
1386  {
1387  LOG_DEBUG << "D: backleg=" << cellHit.backleg << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1388  idVectorIter ij=trackIdVec.begin();
1389  while (ij != trackIdVec.end()) { LOG_DEBUG<< " " << *ij<<endm; ij++; }
1390  }
1391  tempVec = erasedVec;
1392  multiHitsCellsVec_temp.clear();
1393  }
1394 
1395  if(Debug())
1396  {
1397  LOG_INFO<<"multiHitsCellsVec size="<<multiHitsCellsVec.size()<<endm;
1398  }
1399 
1402  mtdCellHitVector tempVec_2Trck = multiHitsCellsVec;
1403  mtdCellHitVector erasedVec_2Trck = tempVec_2Trck;
1404  while(tempVec_2Trck.size()!=0)
1405  {
1406  StructCellHit Cellhit;
1407  Int_t ntracks = 0;
1408  idVector vTrackId;
1409  vector<StThreeVectorD> vPosition;
1410  vector<Int_t> vchannel, vbackleg, vmodule, vcell;
1411  vector<Int_t> vflag;
1412  vector<Float_t> vzhit, vyhit;
1413  vector<pairD> vtot;
1414  vector<pairD> vtdc;
1415  vector<Double_t> vtheta;
1416  vector<Double_t> vpathLength;
1417  vector<Double_t> vexpTof2MTD;
1418  vector<Int_t> vindex2MtdHit;
1419  vector<Int_t> vidTruth;
1420 
1421  mtdCellHitVectorIter temp_Iter=tempVec_2Trck.begin();
1422  mtdCellHitVectorIter erased_Iter=erasedVec_2Trck.begin();
1423 
1424  while(erased_Iter!= erasedVec_2Trck.end())
1425  {
1426  if(temp_Iter->backleg == erased_Iter->backleg &&
1427  temp_Iter->module == erased_Iter->module &&
1428  temp_Iter->cell == erased_Iter->cell)
1429  {
1430  ntracks++;
1431  vbackleg.push_back(erased_Iter->backleg);
1432  vmodule.push_back(erased_Iter->module);
1433  vcell.push_back(erased_Iter->cell);
1434  vflag.push_back(erased_Iter->matchFlag);
1435  vPosition.push_back(erased_Iter->hitPosition);
1436  vTrackId.push_back(erased_Iter->trackIdVec.back());
1437  vzhit.push_back(erased_Iter->zhit);
1438  vyhit.push_back(erased_Iter->yhit);
1439  vtot.push_back(erased_Iter->tot);
1440  vtdc.push_back(erased_Iter->leadingEdgeTime);
1441  vindex2MtdHit.push_back(erased_Iter->index2MtdHit);
1442  vtheta.push_back(erased_Iter->theta);
1443  vpathLength.push_back(erased_Iter->pathLength);
1444  vexpTof2MTD.push_back(erased_Iter->expTof2MTD);
1445  vidTruth.push_back(erased_Iter->idTruth);
1446  if(Debug())
1447  {
1448  LOG_INFO<<"ntracks ="<<ntracks<<endm;
1449  LOG_INFO<<"tack INFO: "<<vzhit[ntracks-1]<<" "<<vyhit[ntracks-1]<<endm;
1450  LOG_INFO<<"hits INFO: "<<vbackleg[ntracks-1]<<" "<<vmodule[ntracks-1]<<" "<<vcell[ntracks-1]<<endm;
1451  }
1452  erasedVec_2Trck.erase(erased_Iter);
1453  erased_Iter--;
1454  }
1455  erased_Iter++;
1456  }// end of the multitrack finding loop
1457 
1458  if (ntracks<2)
1459  {
1460  LOG_INFO<<"What? one hit just match with one track in multimatch!!!!!!! "<<endm;
1461  }
1462 
1463  if(ntracks>1)
1464  {
1465  // assign the closest track as the match
1466 
1467  Int_t thiscandidate = -99;
1468  Int_t thisMatchFlag = 0;
1469  Float_t maxR = 9999.;
1470 
1471  for(int j=0; j<ntracks; j++)
1472  {
1473  Int_t ibackleg = vbackleg[j];
1474  Int_t imodule = vmodule[j];
1475  Int_t icell = vcell[j];
1476  Int_t iidTruth = vidTruth[j];
1477 
1478  Float_t trkLocalY = vyhit[j];
1479  Float_t hitLocalY = -999.;
1480  if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell,ibackleg,iidTruth);
1481  else hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell);
1482  Float_t dy = fabs(trkLocalY-hitLocalY);
1483 
1484  Float_t trkGlobalZ = vPosition[j].z();
1485  Float_t hitGlobalZ = getMtdHitGlobalZ(vtdc[j].first, vtdc[j].second, imodule);
1486  Float_t dz = fabs(trkGlobalZ-hitGlobalZ);
1487 
1488  Int_t projMod = getProjModule(vzhit[j],vPosition[j].z());
1489  Float_t dr = abs(imodule-projMod)*gMtdRadiusDiff;
1490 
1491  Float_t rr = 9999.;
1492  if(mCosmicFlag) rr = dy;
1493  else rr = sqrt(dy*dy+dz*dz+dr*dr);
1494  if(Debug())
1495  {
1496  LOG_INFO<<"hit information::"<<" "<<ibackleg<<" "<<imodule<<" "<<icell<<endm;
1497  LOG_INFO<<"track information::"<<" "<<vTrackId[j] << " " << trkGlobalZ<<" "<<trkLocalY<<endm;
1498  LOG_INFO<<"distance::"<<" "<<dz<<" "<<dy<<" "<<rr<<endm;
1499  }
1500 
1501  if(rr<maxR)
1502  {
1503  maxR = rr;
1504  thiscandidate = j;
1505  thisMatchFlag = (vflag[j]==1? 1 : 7);
1506  }
1507  }
1508 
1509  if (thiscandidate>=0)
1510  {
1511  Cellhit.backleg = vbackleg[thiscandidate];
1512  Cellhit.module = vmodule[thiscandidate];
1513  Cellhit.cell = vcell[thiscandidate];
1514  Cellhit.trackIdVec.push_back(vTrackId[thiscandidate]);
1515  Cellhit.hitPosition = vPosition[thiscandidate];
1516  Cellhit.matchFlag = thisMatchFlag;
1517  Cellhit.zhit = vzhit[thiscandidate];
1518  Cellhit.yhit = vyhit[thiscandidate];
1519  Cellhit.tot = vtot[thiscandidate];
1520  Cellhit.leadingEdgeTime = vtdc[thiscandidate];
1521  Cellhit.index2MtdHit = vindex2MtdHit[thiscandidate];
1522  Cellhit.theta = vtheta[thiscandidate];
1523  Cellhit.pathLength = vpathLength[thiscandidate];
1524  Cellhit.expTof2MTD = vexpTof2MTD[thiscandidate];
1525  Cellhit.idTruth = vidTruth[thiscandidate];
1526  singleHitCellsVec.push_back(Cellhit);
1527  if(Debug())
1528  {
1529  LOG_INFO<<"flag = "<<thisMatchFlag<<endm;
1530  LOG_INFO<<"Pick track " << vTrackId[thiscandidate] << " with local z" << vzhit[thiscandidate] << endm;
1531  }
1532  }
1533  }
1534  tempVec_2Trck = erasedVec_2Trck;
1535  }// end of the multihitscellsvec loop
1536 }
1537 
1538 
1540 void StMtdMatchMaker::finalMatchedMtdHits(mtdCellHitVector& singleHitCellsVec,mtdCellHitVector& finalMatchedCellsVec){
1541  mtdCellHitVector tempVec = singleHitCellsVec;
1542  mtdCellHitVector erasedVec = tempVec;
1543 
1544  while (tempVec.size() != 0)
1545  {
1546  StructCellHit cellHit;
1547  Int_t nCells = 0;
1548  idVector vTrackId;
1549  vector<StThreeVectorD> vPosition;
1550  vector<Int_t> vchannel, vbackleg, vmodule, vcell;
1551  vector<Float_t> vzhit, vyhit;
1552  vector<pairD> vtot;
1553  vector<pairD> vtdc;
1554  vector<Double_t> vtheta;
1555  vector<Double_t> vpathLength;
1556  vector<Double_t> vexpTof2MTD;
1557  vector<Int_t> vindex2MtdHit;
1558  vector<Int_t> vflag;
1559  vector<Int_t> vidTruth;
1560  mtdCellHitVectorIter tempIter=tempVec.begin();
1561  mtdCellHitVectorIter erasedIter=erasedVec.begin();
1562  while(erasedIter!= erasedVec.end())
1563  {
1564  if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back())
1565  {
1566  nCells++;
1567  vbackleg.push_back(erasedIter->backleg);
1568  vmodule.push_back(erasedIter->module);
1569  vcell.push_back(erasedIter->cell);
1570  vPosition.push_back(erasedIter->hitPosition);
1571  vTrackId.push_back(erasedIter->trackIdVec.back());
1572  vzhit.push_back(erasedIter->zhit);
1573  vyhit.push_back(erasedIter->yhit);
1574  vtot.push_back(erasedIter->tot);
1575  vtdc.push_back(erasedIter->leadingEdgeTime);
1576  vindex2MtdHit.push_back(erasedIter->index2MtdHit);
1577  vtheta.push_back(erasedIter->theta);
1578  vflag.push_back(erasedIter->matchFlag);
1579  vpathLength.push_back(erasedIter->pathLength);
1580  vexpTof2MTD.push_back(erasedIter->expTof2MTD);
1581  vidTruth.push_back(erasedIter->idTruth);
1582  if(Debug())
1583  {
1584  LOG_INFO<<"flag 1 ::"<<" "<<nCells<<" "<<erasedIter->matchFlag<<endm;
1585  }
1586 
1587  erasedVec.erase(erasedIter);
1588  --erasedIter;
1589  }
1590  ++erasedIter;
1591  }
1592 
1593  if (nCells==1)
1594  {
1595  // for singly hit cell, copy data in singleHitCellsVec
1596  cellHit.backleg = vbackleg[0];
1597  cellHit.module = vmodule[0];
1598  cellHit.cell = vcell[0];
1599  cellHit.trackIdVec.push_back(vTrackId[0]);
1600  cellHit.hitPosition = vPosition[0];
1601  cellHit.matchFlag = vflag[0];
1602  cellHit.zhit = vzhit[0];
1603  cellHit.yhit = vyhit[0];
1604  cellHit.tot = vtot[0];
1605  cellHit.leadingEdgeTime = vtdc[0];
1606  cellHit.index2MtdHit = vindex2MtdHit[0];
1607  cellHit.theta = vtheta[0];
1608  cellHit.pathLength = vpathLength[0];
1609  cellHit.expTof2MTD = vexpTof2MTD[0];
1610  cellHit.idTruth = vidTruth[0];
1611  finalMatchedCellsVec.push_back(cellHit);
1612 
1613  if(Debug())
1614  {
1615  LOG_INFO<<"flag 2::"<<vflag[0]<<endm;
1616  }
1617 
1618  // debugging output
1619  if(Debug())
1620  {
1621  LOG_DEBUG << "E: ibackleg=" << cellHit.backleg << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:";
1622  idVectorIter ij=vTrackId.begin();
1623  while (ij != vTrackId.end()) { LOG_DEBUG << " " << *ij; ij++; }
1624  LOG_DEBUG << endm;
1625  }
1626  }
1627  else if (nCells>1)
1628  {
1629  // for multiple hit cells find the most likely candidate.
1630  Int_t thiscandidate = -99;
1631  Int_t thisMatchFlag = 0;
1632 
1633  // sort on tot
1634  vector<Int_t> ttCandidates;
1635  for (Int_t i=0;i<nCells;i++)
1636  {
1637  pair<Double_t,Double_t> tt = vtot[i];
1638  if(tt.first<40.) ttCandidates.push_back(i);
1639  }
1640  if (ttCandidates.size()==1)
1641  {
1642  if(Debug())
1643  {
1644  LOG_INFO<<"flag 3 ::"<<vflag[ttCandidates[0]]<<endm;
1645  }
1646  thiscandidate = ttCandidates[0];
1647  thisMatchFlag = vflag[ttCandidates[0]] + 1;
1648  }
1649  else if (ttCandidates.size()>1)
1650  { // sort on hitposition
1651  Float_t ss(9999.);
1652  vector<Int_t> ssCandidates;
1653  for(size_t j=0;j<ttCandidates.size();j++)
1654  {
1655  Int_t ibackleg = vbackleg[ttCandidates[j]];
1656  Int_t imodule = vmodule[ttCandidates[j]];
1657  Int_t icell = vcell[ttCandidates[j]];
1658  Int_t iidTruth = vidTruth[ttCandidates[j]];
1659 
1660  Float_t trkLocalY = vyhit[ttCandidates[j]];
1661  Float_t hitLocalY = -999.;
1662  if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell,ibackleg,iidTruth);
1663  else hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell);
1664  Float_t dy = fabs(trkLocalY-hitLocalY);
1665 
1666  Float_t trkGlobalZ = vPosition[ttCandidates[j]].z();
1667  Float_t hitGlobalZ = getMtdHitGlobalZ(vtdc[ttCandidates[j]].first, vtdc[ttCandidates[j]].second, imodule);
1668  Float_t dz = fabs(trkGlobalZ-hitGlobalZ);
1669 
1670  Int_t projMod = getProjModule(vzhit[ttCandidates[j]],vPosition[ttCandidates[j]].z());
1671  Float_t dr = abs(imodule-projMod)*gMtdRadiusDiff;
1672 
1673  Float_t rr = 9999.;
1674  if(mCosmicFlag) rr = dy;
1675  else rr = sqrt(dy*dy+dz*dz+dr*dr);
1676  if(rr<ss)
1677  {
1678  ss = rr;
1679  ssCandidates.clear();
1680  ssCandidates.push_back(ttCandidates[j]);
1681  }
1682  else if(rr==ss)
1683  ssCandidates.push_back(ttCandidates[j]);
1684  }
1685  if (ssCandidates.size()>=1)
1686  {
1687  if(Debug())
1688  {
1689  LOG_INFO<<"flag 4::"<<ssCandidates.size()<<" "<<vflag[ttCandidates[0]]<<endm;
1690  }
1691  thiscandidate = ssCandidates[0];
1692  thisMatchFlag = vflag[ssCandidates[0]]+2;
1693  }
1694  }
1695 
1696  if (thiscandidate>=0)
1697  {
1698  cellHit.backleg = vbackleg[thiscandidate];
1699  cellHit.module = vmodule[thiscandidate];
1700  cellHit.cell = vcell[thiscandidate];
1701  cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
1702  cellHit.hitPosition = vPosition[thiscandidate];
1703  cellHit.matchFlag = thisMatchFlag;
1704  cellHit.zhit = vzhit[thiscandidate];
1705  cellHit.yhit = vyhit[thiscandidate];
1706  cellHit.tot = vtot[thiscandidate];
1707  cellHit.leadingEdgeTime = vtdc[thiscandidate];
1708  cellHit.index2MtdHit = vindex2MtdHit[thiscandidate];
1709  cellHit.theta = vtheta[thiscandidate];
1710  cellHit.pathLength = vpathLength[thiscandidate];
1711  cellHit.expTof2MTD = vexpTof2MTD[thiscandidate];
1712  cellHit.idTruth = vidTruth[thiscandidate];
1713 
1714  finalMatchedCellsVec.push_back(cellHit);
1715 
1716  // debugging output
1717  if(Debug()) { LOG_DEBUG << "E: ibackleg=" << cellHit.backleg << " imodule=" << cellHit.module << " icell=" << cellHit.cell << "\ttrackid:" << vTrackId[thiscandidate] << endm; }
1718  }
1719  }
1720  else
1721  {
1722  LOG_WARN << "E: no cells belong to this track ... should not happen!" << endm;
1723  }
1724 
1725  tempVec = erasedVec;
1726  } // end of Sect.E
1727 }
1728 
1729 
1731 void StMtdMatchMaker::fillPidTraits(mtdCellHitVector& finalMatchedCellsVec,Int_t& nValidSingleHitCells,Int_t& nValidSinglePrimHitCells){
1732 
1733  for (size_t ii=0; ii < finalMatchedCellsVec.size(); ii++){
1734  Int_t backleg = finalMatchedCellsVec[ii].backleg;
1735  Int_t module = finalMatchedCellsVec[ii].module;
1736  Int_t cell = finalMatchedCellsVec[ii].cell;
1737 
1738  if (finalMatchedCellsVec[ii].trackIdVec.size()!=1)
1739  LOG_WARN << "F: WHAT!?! mult.matched cell in single cell list " << backleg << " " << module << " " << cell << endm;
1740 
1741  Float_t trkLocalY = finalMatchedCellsVec[ii].yhit;
1742  Float_t trkLocalZ = finalMatchedCellsVec[ii].zhit;
1743  Float_t trkGlobalZ = finalMatchedCellsVec[ii].hitPosition.z();
1744 
1745  Int_t hitIdTruth = finalMatchedCellsVec[ii].idTruth;
1746  Float_t hitLocalY = -999.;
1747  if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(backleg,module)->GetCellLocalYCenter(cell,backleg,hitIdTruth);
1748  else hitLocalY = mMtdGeom->GetGeoModule(backleg,module)->GetCellLocalYCenter(cell);
1749  Float_t LeTimeWest = finalMatchedCellsVec[ii].leadingEdgeTime.first;
1750  Float_t LeTimeEast = finalMatchedCellsVec[ii].leadingEdgeTime.second;
1751  Float_t hitGlobalZ = getMtdHitGlobalZ(LeTimeWest, LeTimeEast, module);
1752 
1753  Float_t dy = trkLocalY - hitLocalY;
1754  Float_t dz = trkGlobalZ - hitGlobalZ;
1755 
1756 
1757  if(mHisto) {
1758  mTracksPerCellMatch3->Fill(finalMatchedCellsVec[ii].trackIdVec.size());
1759  mDeltaHitMatch3->Fill(dy, dz);
1760  int hisIndex = backleg-1;
1761  mDeltaHitFinal[hisIndex]->Fill(dy,dz);
1762  }
1763 
1764  // get track-id from cell hit vector
1765  Int_t trackNode = finalMatchedCellsVec[ii].trackIdVec[0];
1766  if(mMuDstIn){
1767 
1768  LOG_DEBUG<<"In StMuDst mode: mtd hit matched with track successfully : track nodeId:"<<finalMatchedCellsVec[ii].trackIdVec[0]<<" mtd hitId:"<<finalMatchedCellsVec[ii].index2MtdHit<<endm;
1769 
1770  StMuTrack *gTrack = mMuDst->globalTracks(trackNode);
1771  if(!gTrack) {
1772  LOG_WARN << "Wrong global track!" << endm;
1773  continue;
1774  }
1775  StMuMtdHit *mtdHit = mMuDst->mtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1776  if(mtdHit->backleg()!=backleg || mtdHit->module()!=module || mtdHit->cell()!=cell) {
1777  LOG_WARN << "Wrong hit in the MtdHitCollection!" << endm;
1778  continue;
1779  }
1780  nValidSingleHitCells++;
1781 
1782  mtdHit->setAssociatedTrackKey(gTrack->id());
1783  mtdHit->setIndex2Global(trackNode);
1784  gTrack->setIndex2MtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1785 
1786  StMuMtdPidTraits pidMtd = gTrack->mtdPidTraits();
1787  pidMtd.setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1788  pidMtd.setYLocal(trkLocalY);
1789  pidMtd.setZLocal(trkLocalZ);
1790  pidMtd.setDeltaY(dy);
1791  pidMtd.setDeltaZ(dz);
1792  pidMtd.setThetaLocal(finalMatchedCellsVec[ii].theta);
1793  pidMtd.setPosition(finalMatchedCellsVec[ii].hitPosition);
1794  pidMtd.setPathLength(finalMatchedCellsVec[ii].pathLength);
1795  pidMtd.setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1796  gTrack->setMtdPidTraits(pidMtd);
1797 
1798  Int_t pIndex = -999;
1799  map<Int_t, Int_t>::iterator it = index2Primary.find(trackNode);
1800  if(it!=index2Primary.end()){
1801  pIndex = it->second;
1802  }
1803  StMuTrack *pTrack = (StMuTrack *)mMuDst->array(muPrimary)->UncheckedAt(pIndex);
1804  if(pTrack && pIndex>-1){
1805  mtdHit->setIndex2Primary(pIndex);
1806  pTrack->setIndex2MtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1807  StMuMtdPidTraits ppidMtd = pTrack->mtdPidTraits();
1808  ppidMtd.setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1809  ppidMtd.setYLocal(trkLocalY);
1810  ppidMtd.setZLocal(trkLocalZ);
1811  ppidMtd.setDeltaY(dy);
1812  ppidMtd.setDeltaZ(dz);
1813  ppidMtd.setThetaLocal(finalMatchedCellsVec[ii].theta);
1814  ppidMtd.setPosition(finalMatchedCellsVec[ii].hitPosition);
1815  ppidMtd.setPathLength(finalMatchedCellsVec[ii].pathLength);
1816  ppidMtd.setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1817  pTrack->setMtdPidTraits(ppidMtd);
1818  }
1819 
1820  }else{
1821 
1822  LOG_INFO<<"In StEvent mode: mtd hit matched with track successfully : track nodeId:"<<finalMatchedCellsVec[ii].trackIdVec[0]<<" mtd hitId:"<<finalMatchedCellsVec[ii].index2MtdHit<<endm;
1823 
1824  // get track-id from cell hit vector
1825  StSPtrVecTrackNode &nodes = mEvent->trackNodes();
1826  StGlobalTrack *globalTrack = dynamic_cast<StGlobalTrack*>(nodes[trackNode]->track(global));
1827  if(!globalTrack) {
1828  LOG_WARN << "Wrong global track!" << endm;
1829  continue;
1830  }
1831  StPrimaryTrack *primaryTrack =dynamic_cast<StPrimaryTrack*>(globalTrack->node()->track(primary));
1832 
1833  // Fill association in MTD Hit Collection
1834  StMtdCollection *theMtd = mEvent->mtdCollection();
1835  StSPtrVecMtdHit& mtdHits = theMtd->mtdHits();
1836  StMtdHit *mtdHit = mtdHits[finalMatchedCellsVec[ii].index2MtdHit];
1837  if(mtdHit->backleg()!=backleg || mtdHit->module()!=module || mtdHit->cell()!=cell) {
1838  LOG_WARN << "Wrong hit in the MtdHitCollection!" << endm;
1839  continue;
1840  }
1841  nValidSingleHitCells++;
1842 
1843  mtdHit->setAssociatedTrack(globalTrack);
1844 
1845  StMtdPidTraits *pidMtd = new StMtdPidTraits();
1846  pidMtd->setMtdHit(mtdHit);
1847  pidMtd->setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1848  pidMtd->setYLocal(trkLocalY);
1849  pidMtd->setZLocal(trkLocalZ);
1850  pidMtd->setDeltaY(dy);
1851  pidMtd->setDeltaZ(dz);
1852  pidMtd->setThetaLocal(finalMatchedCellsVec[ii].theta);
1853  pidMtd->setPosition(finalMatchedCellsVec[ii].hitPosition);
1854  pidMtd->setPathLength(finalMatchedCellsVec[ii].pathLength);
1855  pidMtd->setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1856  globalTrack->addPidTraits(pidMtd);
1857 
1858  if(primaryTrack){
1859  StMtdPidTraits *ppidMtd = new StMtdPidTraits();
1860  ppidMtd->setMtdHit(mtdHit);
1861  ppidMtd->setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1862  ppidMtd->setYLocal(trkLocalY);
1863  ppidMtd->setZLocal(trkLocalZ);
1864  ppidMtd->setDeltaY(dy);
1865  ppidMtd->setDeltaZ(dz);
1866  ppidMtd->setThetaLocal(finalMatchedCellsVec[ii].theta);
1867  ppidMtd->setPosition(finalMatchedCellsVec[ii].hitPosition);
1868  ppidMtd->setPathLength(finalMatchedCellsVec[ii].pathLength);
1869  ppidMtd->setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1870  primaryTrack->addPidTraits(ppidMtd);
1871  }
1872  }
1873  }
1874 }
1875 
1878  return validTrack(MtdTrack(track));
1879 }
1880 
1883  return validTrack(MtdTrack(track));
1884 }
1885 
1889  float gpt = mtt.pt;
1890  float geta = mtt.eta;
1891  int gnFtPts = mtt.nFtPts;
1892  int gnDedxPts = mtt.nDedxPts;
1893  if(gnDedxPts<mMindEdxFitPoints) return kFALSE;
1894  if(geta<mMinEta||geta>mMaxEta) return kFALSE;
1895  if(gpt<mMinPt) return kFALSE;
1897  if (mtt.flag<=0 || mtt.flag>=1000) return kFALSE;
1898 
1900  // if (mtt->topologyMap().numberOfHits(kTpcId) < mMinHitsPerTrack) return kFALSE;
1902  if (gnFtPts < mMinFitPointsPerTrack) return kFALSE;
1904  //fg float ratio = (1.0*mtt->fitTraits().numberOfFitPoints(kTpcId)) / (1.0*mtt->numberOfPossiblePoints(kTpcId));
1905  float ratio = 99.;
1906  if(mtt.nHitsPoss!=0) ratio = gnFtPts / (1.0*mtt.nHitsPoss);
1907  if (ratio < mMinFitPointsOverMax) return kFALSE;
1908 
1909  return kTRUE;
1910 }
1911 
1913 Float_t StMtdMatchMaker::getMtdHitGlobalZ(Float_t leadingWestTime, Float_t leadingEastTime, Int_t module)
1914 {
1915  return (module-3)*gMtdCellLength + (leadingEastTime - leadingWestTime)/2./gMtdCellDriftV*1e3;
1916 }
1917 
1919 Int_t StMtdMatchMaker::getProjModule(Float_t local_z, Float_t global_z)
1920 {
1921  int module = 0;
1922  for(int i=0; i<5; i++)
1923  {
1924  if(abs(global_z-(i-2)*gMtdCellLength-local_z)<5)
1925  {
1926  module = i+1;
1927  break;
1928  }
1929  }
1930  return module;
1931 }
1932 
1933 //___________________________________________________
1934 
1935 MtdTrack::MtdTrack(StTrack *stt){
1936 
1937  pt = -999.; eta = -999.; nFtPts = 0;
1938  nDedxPts = 0; flag = 0; nHitsPoss = 999;
1939  if(stt){
1940  pt = stt->geometry()->momentum().perp();
1941  eta = stt->geometry()->momentum().pseudoRapidity();
1942  nFtPts = stt->fitTraits().numberOfFitPoints(kTpcId);
1943  static StTpcDedxPidAlgorithm PidAlgorithm;
1944  const StParticleDefinition* pd=stt->pidTraits(PidAlgorithm);
1945  if(pd){
1946  if(PidAlgorithm.traits()){
1947  nDedxPts=PidAlgorithm.traits()->numberOfPoints();
1948  }
1949  }
1950  flag = stt->flag();
1951  nHitsPoss = stt->numberOfPossiblePoints(kTpcId);
1952  }
1953 }
1954 //------------------------------------------------
1955 MtdTrack::MtdTrack(StMuTrack *mut){
1956 
1957  pt = -999.; eta = -999.; nFtPts = 0;
1958  nDedxPts = 0; flag = 0; nHitsPoss = 999;
1959  if(mut){
1960  pt = mut->momentum().perp();
1961  eta = mut->momentum().pseudoRapidity();
1962  nFtPts = mut->nHitsFit(kTpcId);
1963  nDedxPts = mut->nHitsDedx();
1964  flag = mut->flag();
1965  nHitsPoss = mut->nHitsPoss(kTpcId);
1966  }
1967 }
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StThreeVectorF primaryVertexPosition(int vtx_id=-1) const
The StMuDst is supposed to be structured in &#39;physical events&#39;. Therefore there is only 1 primary vert...
Definition: StMuEvent.cxx:221
Int_t index2Global() const
Returns index of associated global track. If not in order can be set with StMuDst::fixTrackIndeces() ...
Definition: StMuTrack.h:231
short id() const
Returns the track id(or key), is unique for a track node, i.e. global and primary tracks have the sam...
Definition: StMuTrack.h:228
virtual Int_t Init()
initialize drifting velocity and histograms.
Float_t mMinFitPointsOverMax
minimum dE/dx fit points per track
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
Definition: StMuTrack.h:238
MTD track class.
virtual void Clear(Option_t *option="")
User defined functions.
Definition: StMaker.cxx:634
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
bool validTrack(StTrack *track)
check track quality
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
Float_t mMinEta
minimum ratio
StMtdMatchMaker(const char *name="MtdMatch")
Default constructor.
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
Definition: Stypes.h:40
Int_t getProjModule(Float_t local_z, Float_t global_z)
calculate module of projected position
Float_t mMinPt
maximum pseudorapidity
static TClonesArray * array(int type)
returns pointer to the n-th TClonesArray
Definition: StMuDst.h:262
Float_t mMaxEta
minimum pseudorapidity
virtual Int_t Make()
associate tracks with mtd hits
Int_t InitRun(int runnumber)
InitRun: initialize geometries (retrieve beam line constraint from database)
virtual Int_t Finish()
write QA histograms
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
Int_t mMindEdxFitPoints
minimum fit points per track
Float_t getMtdHitGlobalZ(Float_t leadingWestTime, Float_t leadingEastTime, Int_t module)
calculate global z of the MTD hits based on its timing information
Int_t FinishRun(int runnumber)
FinishRun: clean up BeamHelix (will be reinitialized at the next initRun)