StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPeCPair.cxx
1 //
3 // $Id: StPeCPair.cxx,v 1.19 2015/07/22 19:55:20 ramdebbe Exp $
4 // $Log: StPeCPair.cxx,v $
5 // Revision 1.19 2015/07/22 19:55:20 ramdebbe
6 // stopped doing extrapolation to TOF, now copy information from StMuBTofPidTraits, needs to do the same for fill(mudst, stevent)
7 //
8 // Revision 1.18 2014/12/23 20:43:10 ramdebbe
9 // copied the fill functionality of method with both inputs to the one with MuDst input. This gives TOF extrapolation in the pPairs branch
10 //
11 // Revision 1.17 2013/12/27 16:52:20 ramdebbe
12 // added a input argument StBTofGeometry to fill method (StMuDst + StEvent) and x y z coordinates of intercept to TOF cylinder
13 //
14 // Revision 1.16 2012/06/13 15:10:22 ramdebbe
15 // added tof information to both tracks
16 //
17 // Revision 1.15 2004/01/26 23:01:03 perev
18 // WarnOff
19 //
20 // Revision 1.14 2003/11/25 01:54:33 meissner
21 // correct several bugs: eta cut for tracks, charge sorting, add counting of FTPC and TPC primary tracks, Add bbc information
22 //
23 // Revision 1.13 2003/09/02 17:58:46 perev
24 // gcc 3.2 updates + WarnOff
25 //
26 // Revision 1.12 2003/03/19 15:35:47 yepes
27 // *** empty log message ***
28 //
29 // Revision 1.11 2003/03/18 21:20:41 yepes
30 // correcting problem with bField
31 //
32 // Revision 1.10 2003/02/05 17:14:05 yepes
33 // Adding bField and pPairs.psi to tree
34 //
35 // Revision 1.9 2002/12/16 23:05:38 yepes
36 // *** empty log message ***
37 //
38 // Revision 1.8 2002/12/16 23:04:02 yepes
39 // Field comes in KGauss and should be passed to routines in Teslas
40 // problem pointed out by Vladimir
41 //
42 //
43 // Revision 1.7 2002/03/19 22:23:49 meissner
44 // New variables: zdc unatt., Trigger word, MC tree if Geant Branch, DCA for primary pairs, all tracks for secondary pairs (Test)
45 //
46 // Revision 1.6 2001/02/21 20:42:12 yepes
47 // Add ctb signals to tree
48 //
49 // Revision 1.5 2001/02/13 17:54:43 yepes
50 // still problems on differnt platforms
51 //
52 // Revision 1.4 2001/02/12 21:15:59 yepes
53 // New version of StPeCMaker, lots of changes
54 //
55 // Revision 1.1 2000/04/21 19:12:30 nystrand
56 // First Version
57 //
58 // Revision 1.1 2000/03/24 22:37:06 nystrand
59 // First version of StPeCPair
60 //
61 // Revision 1.0 2000/03/20 23:28:50 nystrand
62 //
64 #include <Stiostream.h>
65 #include "StPeCPair.h"
66 #include "StEventTypes.h"
67 #include "StEmcUtil/geometry/StEmcGeom.h"
68 #include "StEmcUtil/projection/StEmcPosition.h"
69 #include "StEmcUtil/filters/StEmcFilter.h"
70 #include "StEmcUtil/geometry/StEmcGeom.h"
71 #include "StEmcUtil/others/emcDetectorName.h"
72 #include "St_geant_Maker/St_geant_Maker.h"
73 
74 
75 
76 ClassImp(StPeCPair)
77 
78 StPeCPair::StPeCPair() {
79 }
80 
81 StPeCPair::~StPeCPair() {
82 }
83 
84 void StPeCPair::Clear(const char *) {
85  pCharge=0;
86  pPt=0.;
87  pPz =0.;
88  pPsi =0.;
89  pAngle =0.;
90  pXyAngle =0.;
91  pPtArm =0.;
92  pAlpha =0.;
93  pPartDca =0.;
94  pV0Dca =0.;
95  rV0 =0.;
96  phiV0 =0.;
97  zV0 =0.;
98 
99  // Not implemented
100  // tr1.Clear();
101  // tr2.Clear();
102  // pionH.Clear();
103  // kaonH.Clear();
104  // protonH.Clear();
105  // electronH.Clear();
106  // muonH.Clear();
107 
108  track1=NULL;
109  track2=NULL;
110  muTrack1=NULL;
111  muTrack2=NULL;
112 }
113 
114 #ifndef __CINT__
115 StPeCPair::StPeCPair ( StTrack* trk1, StTrack* trk2,
116  Bool_t primaryFlag, StEvent* event ) {
117  this->Clear();
118 
119  if(trk1->geometry()->charge() < 0 && trk2->geometry()->charge()>0 ) { // swap to 1+ 2- if different charges
120  track1 = trk2;
121  track2 = trk1;
122  } else { // keep +-, ++, --
123  track1 = trk1;
124  track2 = trk2;
125  }
126  muTrack1=0; // clean, might crash otherwise
127  muTrack2=0;
128  fill ( primaryFlag, event ) ;
129 }
130 
131 StPeCPair::StPeCPair ( StMuTrack* trk1, StMuTrack* trk2,
132  Bool_t primaryFlag, StMuEvent* event, StBTofGeometry * pairTOFgeo ) {
133  this->Clear();
134  //
135  // Get Magnetic field from event summary
136  //
137  Double_t bFld;
138 
139  pairTOFgeoLocal = pairTOFgeo;
140 
141  bFld=event->eventSummary().magneticField()/10.;
142 
143 
144 
145  if (trk1->charge() <0 && trk2->charge() >0 ) { // swap to 1+ 2- if different charges
146  muTrack1 = trk2;
147  muTrack2 = trk1;
148  } else { // keep for +- and ++, --
149  muTrack1 = trk1;
150  muTrack2 = trk2;
151  }
152  track1=0; // clean; might crash otherwise
153  track2=0;
154  fill ( primaryFlag, event ) ;
155 }
156 
157 StPeCPair::StPeCPair ( StMuTrack* trk1, StMuTrack* trk2,
158  Bool_t primaryFlag, StMuEvent* event, StEvent* eventP, StBTofGeometry * pairTOFgeo ) {
159  this->Clear();
160  //
161  // Get Magnetic field from event summary
162  //
163 
164 
165  bFld=event->eventSummary().magneticField()/10.;
166  pairTOFgeoLocal = pairTOFgeo;
167  cout<<" StPeCPair StEvent StMuEvent ************ Mag field from summary: "<<bFld<<" TOF geometry pointer: "<<pairTOFgeo<<endl;
168 // // Get EMC calorimeter clusters from StEvent
169 
170 // // check if there is a collection
171 // StEmcCollection *emcStEvent = eventP->emcCollection();
172 
173 // //
174 // //EMC detector information
175 // //
176 // StEmcDetector* EMCdetector = emcStEvent->detector(kBarrelSmdEtaStripId);
177 // if (!EMCdetector)
178 // {cout<<"There is no kBarrelSmdEtaStripId Detector ---------"<<endl; return;}
179 
180 
181  if (trk1->charge() <0 && trk2->charge() >0 ) { // swap to 1+ 2- if different charges
182  muTrack1 = trk2;
183  muTrack2 = trk1;
184  } else { // keep for +- and ++, --
185  muTrack1 = trk1;
186  muTrack2 = trk2;
187  }
188  track1=0; // clean; might crash otherwise
189  track2=0;
190  fill ( primaryFlag, event, eventP ) ;
191 }
192 
193 
194 
195 void StPeCPair::setTrack1(StTrack* trk) {
196  track1 = trk;
197 }
198 void StPeCPair::setTrack2(StTrack* trk) {
199  track2 = trk;
200 }
201 
202 void StPeCPair::setTrack1(StMuTrack* trk) {
203  muTrack1 = trk;
204 }
205 void StPeCPair::setTrack2(StMuTrack* trk) {
206  muTrack2 = trk;
207 }
208 
209 StTrack* StPeCPair::getTrack1() { return track1; }
210 StTrack* StPeCPair::getTrack2() { return track2; }
211 
212 StMuTrack* StPeCPair::getMuTrack1() { return muTrack1; }
213 StMuTrack* StPeCPair::getMuTrack2() { return muTrack2; }
214 
215 StLorentzVectorF StPeCPair::getPair4Momentum(StPeCSpecies pid) const{
216  StLorentzVectorF p4pair(0.0,0.0,0.0,0.0);
217  if ( pid == 0 ) return pionH.Mom4 ;
218  else if ( pid == 1 ) return kaonH.Mom4 ;
219  else if ( pid == 2 ) return protonH.Mom4 ;
220  else if ( pid == 3 ) return electronH.Mom4 ;
221  else if ( pid == 4 ) return muonH.Mom4 ;
222  else {
223  printf ( "StPeCPair::getPair4Momentum: wrong pid %d \n", pid ) ;
224  return p4pair ;
225  }
226 }
227 #endif /*__CINT__*/
228 
229 
230 
231 Int_t StPeCPair::fill ( Bool_t primaryFlag, StEventSummary* summary,
232  StThreeVectorF& p1, StPhysicalHelixD& h1, short charge1,
233  StThreeVectorF& p2, StPhysicalHelixD& h2, short charge2,
234  StThreeVectorF& primaryVertexPosition ) {
235 //
236 // Check whether tracks are primary or secondary
237 // if they are secondary find point of closest approach
238 // and work with momentum at that point
239 //
240  pPartDca = 0. ;
241  pV0Dca = 0. ;
242  rV0 = 0. ;
243  phiV0 = 0. ;
244  zV0 = 0. ;
245  // Want the DCA Vertex info also for the Primary Pair !!
246  // if ( !primaryFlag ) {
247  pairD dcaLengths ;
248  dcaLengths = h1.pathLengths(h2);
249 
250  Float_t bField ;
251  if ( summary != 0 ) bField = summary->magneticField();
252  else bField = 2.5 ;
253 
254  // The momentum we do not need for the primary pair ....
255  if ( !primaryFlag ) {
256  p1 = h1.momentumAt(dcaLengths.first, tesla*bField*0.1 ) ;
257  p2 = h2.momentumAt(dcaLengths.second, tesla*bField*0.1 ) ;
258  }
259 
260  StThreeVectorD x1 = h1.at(dcaLengths.first);
261  StThreeVectorD x2 = h2.at(dcaLengths.second);
262  StThreeVectorD x = (x1-x2) ;
263  pPartDca = x.mag();
264  //
265  // Construct a helix with very large momentums
266  // to get intersection of V0 with vertex
267  //
268  StThreeVectorD xMean = (x1+x2)/2. ;
269  rV0 = xMean.perp();
270  phiV0 = xMean.phi();
271  zV0 = xMean.z();
272 
273  StThreeVectorD pSum = p1+p2 ;
274  StPhysicalHelixD v0Helix ( pSum, xMean, 0.1*bField/tesla, 100000./GeV ) ;
275  pV0Dca = v0Helix.distance ( primaryVertexPosition ) ;
276  //
277  StThreeVectorF p = p1 + p2 ;
278  pPt = p.perp() ;
279  pPz = p.z() ;
280  pPsi = p.phi();
281  // Opening angle of the pair in the CM (lab)
282  // cos(theta) = (p1"dot"p2)/(abs(p1)*abs(p2))
283  Float_t ScalarProduct = p1*p2;
284  Float_t Denominator = p1.mag()*p2.mag();
285  if(Denominator) {
286  pAngle = acos(ScalarProduct/Denominator);
287  }else{
288  pAngle = -999;
289  }
290  if (p1.perp() * p2.perp()) {
291  pXyAngle = acos((p1.x()*p2.x()+p1.y()*p2.y())/p1.perp()/p2.perp());
292  } else {
293  pXyAngle = -999;
294  }
295 //
296 // Calculate Armenteros variables
297 //
298  Float_t p1AlongPtot = p*p1/p.mag() ;
299  Float_t p2AlongPtot = p*p2/p.mag() ;
300 
301  Float_t pt1Ptot = ::sqrt(p1.mag()*p1.mag()-p1AlongPtot*p1AlongPtot);
302 
303  pPtArm = pt1Ptot ;
304  pAlpha = (p1AlongPtot-p2AlongPtot)/(p1AlongPtot+p2AlongPtot);
305 
306 
307  Float_t mptcle=0.0;
308 //
309 // Loop over different species
310 //
311  Float_t mInv, cosThetaStar ;
312  StLorentzVectorF FourMomentum ;
313  StPeCSpec* species ;
314  for ( int i = 0 ; i < nSpecies ; i++ )
315  {
316  if ( i == pion ) {
317  mptcle = pion_plus_mass_c2;
318  species = &pionH ;
319  }
320  else if ( i == kaon ) {
321  mptcle = 493.677*MeV;
322  species = &kaonH ;
323  }
324  else if ( i == proton ) {
325  mptcle = proton_mass_c2;
326  species = &protonH ;
327  }
328  else if ( i == electron ) {
329  mptcle = electron_mass_c2;
330  species = &electronH ;
331  }
332  else if ( i == muon ) {
333  mptcle = 105.6584*MeV;
334  species = &muonH ;
335  }
336  else {
337  printf ( "StPecPair:calculatePair4Momentum; wrong pid %d \n", i ) ;
338  continue ;
339  }
340 
341  StLorentzVectorF p4pair(0.0,0.0,0.0,0.0);
342  Float_t e1 = p1.massHypothesis(mptcle);
343  Float_t e2 = p2.massHypothesis(mptcle);
344  StLorentzVectorF pf1(e1,p1);
345  StLorentzVectorF pf2(e2,p2);
346  p4pair = pf1 + pf2;
347  FourMomentum = p4pair ;
348  mInv = p4pair.m() ;
349 
350  // ThetaStar is the angle between of one of the daughter tracks
351  // and the Z-axis in the Helicity frame. The Helicity frame is
352  // that rest frame of the parent particle in which the direction of
353  // the scattered nucleus is along the negative Z-axis. See K.Schilling,
354  // P.Seyboth, G. Wolf Nucl. Phys. B 15(1970)397-412.
355  // Since the outgoing nuclei are not tagged, this direction cannot
356  // be determined exactly. Because of the low momentum transfers it is,
357  // however, a reasonable approximation to assume that the nuclei move
358  // parallel to the Z-axis of the lab frame.
359 
360 //
361 // Get the sign right for the boost lab --> parent rest frame
362 // Default is in the other direction
363  StThreeVectorF sp = -1.0*p4pair.vect();
364  p4pair.setVect(sp);
365  pf1 = pf1.boost(p4pair);
366  pf2 = pf2.boost(p4pair);
367  Float_t d1th = pf1.cosTheta();
368  Float_t d2th = pf2.cosTheta();
369  // Define cosThetaStar in the interval 0<cosThetaStar<1
370  if( d1th > 0 ) cosThetaStar = d1th;
371  else cosThetaStar = d2th;
372 
373  species->pid = i ;
374  species->mInv = mInv ;
375  species->yRap = FourMomentum.rapidity() ;
376  species->Mom4 = FourMomentum ;
377  species->cosThetaStar = cosThetaStar ;
378  }
379 
380 //
381 // fill our local Track class; not save if track pointers not properly reset !!!
382 // set does not belong here ! FLK works with pointers which are not arguments of this routine !!
383 // #ifndef __CINT__
384 // Int_t prim = 1 ;
385 // if (track1 && track2)
386 // {
387 // tr1.set(prim,track1);
388 // tr2.set(prim,track2);
389 // }
390 // else if (muTrack1 && muTrack2)
391 // {
392 // tr1.set(prim,muTrack1);
393 // tr2.set(prim,muTrack2);
394 // }
395 // #endif
396 
397 
398  return 0 ;
399 }
400 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
401 //
402 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
403 
404 Int_t StPeCPair::fill ( Bool_t primaryFlag, StMuEvent* event ) {
405 
406  pCharge = muTrack1->charge()+muTrack2->charge();
407 
408  StThreeVectorF p1 ;
409  StThreeVectorF p2 ;
410  StPhysicalHelixD h1, hOuter1 ;
411  StPhysicalHelixD h2, hOuter2 ;
412 
413  short charge1, charge2 ;
414 // Int_t mod,eta,sub;
415  StThreeVectorD position,momentum;
416 //
417 
418  tr1_bemcModule = -999;
419  tr1_bemcEtabin = -999;
420  tr1_bemcEtaValue = -999;
421  tr1_bemcPhiValue = -999;
422  tr1_bemcSub = -999;
423  tr1_bemcEnergy = -999;
424  tr2_bemcModule = -999;
425  tr2_bemcEtabin = -999;
426  tr2_bemcEtaValue = -999;
427  tr2_bemcPhiValue = -999;
428  tr2_bemcSub = -999;
429  tr2_bemcEnergy = -999;
430 
431  tr1_bsmdeModule = -999;
432  tr1_bsmdeEtabin = -999;
433  tr1_bsmdeEtaValue = -999;
434  tr1_bsmdePhiValue = -999;
435  tr1_bsmdeSub = -999;
436  tr1_bsmdeEnergy = -999;
437  tr2_bsmdeModule = -999;
438  tr2_bsmdeEtabin = -999;
439  tr2_bsmdeEtaValue = -999;
440  tr2_bsmdePhiValue = -999;
441  tr2_bsmdeSub = -999;
442  tr2_bsmdeEnergy = -999;
443 
444  tr1_timeOfFlight = -999;
445  tr1_pathLength = -999;
446  tr1_Beta = -999;
447  tr2_timeOfFlight = -999;
448  tr2_pathLength = -999;
449  tr2_Beta = -999;
450 
451  tr1_extrapolatedTOF_mX = -999;
452  tr1_extrapolatedTOF_mY = -999;
453  tr1_extrapolatedTOF_mZ = -999;
454 
455  tr2_extrapolatedTOF_mX = -999;
456  tr2_extrapolatedTOF_mY = -999;
457  tr2_extrapolatedTOF_mZ = -999;
458 
459 
460  p1 = muTrack1->momentum();
461  p2 = muTrack2->momentum();
462  charge1 = muTrack1->charge();
463  charge2 = muTrack2->charge();
464  h1 = muTrack1->helix() ;
465  h2 = muTrack2->helix() ;
466  hOuter1 = muTrack1->outerHelix() ;
467  hOuter2 = muTrack2->outerHelix() ;
468 
469  StThreeVectorF vtx = event->primaryVertexPosition() ;
470 
471 
472  StMuBTofPidTraits mBTofPidTraits_1 = muTrack1->btofPidTraits();
473  StMuBTofPidTraits mBTofPidTraits_2 = muTrack2->btofPidTraits();
474 
475  //skip calorimeter cluster information, if necessary, it should be filled with the method using StEvent and MuDst
476 
477 
478 
479  fill ( primaryFlag, &(event->eventSummary()),
480  p1, h1, charge1, p2, h2, charge2, vtx ) ;
481 
482 #ifndef __CINT__
483  tr1.set(1,muTrack1, event); // 1=primary
484  tr2.set(1,muTrack2, event);
485 #endif
486  tr1_timeOfFlight = mBTofPidTraits_1.timeOfFlight();
487  tr1_pathLength = mBTofPidTraits_1.pathLength();
488  tr1_Beta = mBTofPidTraits_1.beta();
489  const float a = mBTofPidTraits_1.position().x();
490  const float b = mBTofPidTraits_1.position().y();
491  const float c = mBTofPidTraits_1.position().z();
492  tr1_extrapolatedTOF_mX = a;
493  tr1_extrapolatedTOF_mY = b;
494  tr1_extrapolatedTOF_mZ = c;
495 
496  tr2_timeOfFlight = mBTofPidTraits_2.timeOfFlight();
497  tr2_pathLength = mBTofPidTraits_2.pathLength();
498  tr2_Beta = mBTofPidTraits_2.beta();
499  const float a2 = mBTofPidTraits_2.position().x();
500  const float b2 = mBTofPidTraits_2.position().y();
501  const float c2 = mBTofPidTraits_2.position().z();
502  tr2_extrapolatedTOF_mX = a2;
503  tr2_extrapolatedTOF_mY = b2;
504  tr2_extrapolatedTOF_mZ = c2;
505  //
506  //extrapolate tracks to TOF
507  //
508  vector<Int_t> idVec;
509  vector<Double_t> pathVec;
510  PointVec crossVec;
511 
512 
513  //no more private extrapolation; commented out 18-June-2015 RD
514 
515 
516 // if(pairTOFgeoLocal->HelixCrossCellIds(hOuter1,idVec,pathVec,crossVec))
517 // {
518 
519 // Int_t cellId = -999;
520 // Int_t moduleId = -999;
521 // Int_t trayId = -999;
522 // pairTOFgeoLocal->DecodeCellId(idVec[0],cellId,moduleId,trayId);
523 // tr1_extrapolatedTOF_mX = crossVec[0].x();
524 // tr1_extrapolatedTOF_mY = crossVec[0].y();
525 // tr1_extrapolatedTOF_mZ = crossVec[0].z();
526 // }
527 // if(pairTOFgeoLocal->HelixCrossCellIds(hOuter2,idVec,pathVec,crossVec))
528 // {
529 // Int_t cellId = -999;
530 // Int_t moduleId = -999;
531 // Int_t trayId = -999;
532 // pairTOFgeoLocal->DecodeCellId(idVec[0],cellId,moduleId,trayId);
533 // tr2_extrapolatedTOF_mX = crossVec[0].x();
534 // tr2_extrapolatedTOF_mY = crossVec[0].y();
535 // tr2_extrapolatedTOF_mZ = crossVec[0].z();
536 // }
537 
538  return 0 ;
539 
540 }
541 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
542 //
543 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
544 
545 Int_t StPeCPair::fill ( Bool_t primaryFlag, StMuEvent* event, StEvent* eventP ) {
546 
547  pCharge = muTrack1->charge()+muTrack2->charge();
548 
549  StThreeVectorF p1 ;
550  StThreeVectorF p2 ;
551  StPhysicalHelixD h1, hOuter1 ;
552  StPhysicalHelixD h2, hOuter2 ;
553 
554  short charge1, charge2 ;
555  Int_t mod,eta,sub;
556  StThreeVectorD position,momentum;
557 
558 
559 
560  tr1_bemcModule = -999;
561  tr1_bemcEtabin = -999;
562  tr1_bemcEtaValue = -999;
563  tr1_bemcPhiValue = -999;
564  tr1_bemcSub = -999;
565  tr1_bemcEnergy = -999;
566  tr2_bemcModule = -999;
567  tr2_bemcEtabin = -999;
568  tr2_bemcEtaValue = -999;
569  tr2_bemcPhiValue = -999;
570  tr2_bemcSub = -999;
571  tr2_bemcEnergy = -999;
572 
573  tr1_bsmdeModule = -999;
574  tr1_bsmdeEtabin = -999;
575  tr1_bsmdeEtaValue = -999;
576  tr1_bsmdePhiValue = -999;
577  tr1_bsmdeSub = -999;
578  tr1_bsmdeEnergy = -999;
579  tr2_bsmdeModule = -999;
580  tr2_bsmdeEtabin = -999;
581  tr2_bsmdeEtaValue = -999;
582  tr2_bsmdePhiValue = -999;
583  tr2_bsmdeSub = -999;
584  tr2_bsmdeEnergy = -999;
585 
586  tr1_timeOfFlight = -999;
587  tr1_pathLength = -999;
588  tr1_Beta = -999;
589  tr2_timeOfFlight = -999;
590  tr2_pathLength = -999;
591  tr2_Beta = -999;
592 
593  tr1_extrapolatedTOF_mX = -999;
594  tr1_extrapolatedTOF_mY = -999;
595  tr1_extrapolatedTOF_mZ = -999;
596 
597  tr2_extrapolatedTOF_mX = -999;
598  tr2_extrapolatedTOF_mY = -999;
599  tr2_extrapolatedTOF_mZ = -999;
600 
601 
602  p1 = muTrack1->momentum();
603  p2 = muTrack2->momentum();
604  charge1 = muTrack1->charge();
605  charge2 = muTrack2->charge();
606  h1 = muTrack1->helix() ;
607  h2 = muTrack2->helix() ;
608  hOuter1 = muTrack1->outerHelix() ;
609  hOuter2 = muTrack2->outerHelix() ;
610 
611 
612 
613  StMuBTofPidTraits mBTofPidTraits_1 = muTrack1->btofPidTraits();
614  StMuBTofPidTraits mBTofPidTraits_2 = muTrack2->btofPidTraits();
615 
616 
617  StThreeVectorF vtx = event->primaryVertexPosition() ;
618  // Get EMC calorimeter clusters from StEvent
619 
620  // check if there is a collection
621  StEmcCollection *emcStEvent = eventP->emcCollection();
622  //
623  //EMC detector information
624  //
625  cout<<" in StPeCPair fill "<<" TOF geo pointer inside fill: "<<pairTOFgeoLocal<<endl;
626  StEmcDetector* EMCdetector = emcStEvent->detector(kBarrelEmcTowerId); //kBarrelSmdEtaStripId);
627  if (!EMCdetector)
628  {cout<<"There is no kBarrelSmdEtaStripId Detector ---------"<<endl;}
629  //
630  // instantiate object to project track to EMC barrel
631  StEmcPosition * project = new StEmcPosition();
632  StEmcGeom * mGeom2=StEmcGeom::instance("bemc");
633  Bool_t ok=project->trackOnEmc(&position,&momentum,muTrack1,bFld);
634 
635  //if(! ok) cout<<" track projection failed ************************************"<<endl;
636  if(ok) {
637  mGeom2->getBin(position.phi(),position.pseudoRapidity(),mod,eta,sub);
638  //
639  //see if EMC has energy
640  //
641  for(unsigned int i=1;i<=120;i++)
642  {
643  if(fabs(mod)!=i) continue;
644  StEmcModule* module=EMCdetector->module(i);
645  StSPtrVecEmcRawHit& hits=module->hits();
646  for(unsigned int k=0;k<hits.size();k++) if(hits[k]){
647  unsigned int module=hits[k]->module();
648  unsigned int Eta=hits[k]->eta();
649 // float energyT=hits[k]->energy();
650  int s=fabs(hits[k]->sub());
651  int did(0);
652  if (module==fabs(mod) && Eta == fabs(eta)){
653  float energyT1=hits[k]->energy();
654  mGeom2->getId(module,Eta,s,did);
655 
656  tr1_bemcModule = module;
657  tr1_bemcEtabin = Eta;
658  tr1_bemcEtaValue = position.pseudoRapidity();
659  tr1_bemcPhiValue = position.phi();
660  tr1_bemcSub = s;
661  tr1_bemcEnergy = energyT1;
662 
663  }
664  }
665  }
666  }
667  //second track
668  Bool_t ok2=project->trackOnEmc(&position,&momentum,muTrack2,bFld);
669  //if(! ok) cout<<" track projection failed ************************************"<<endl;
670  if(ok2) {
671  mGeom2->getBin(position.phi(),position.pseudoRapidity(),mod,eta,sub);
672  //
673  //see if EMC has energy
674  //
675  for(unsigned int i=1;i<=120;i++)
676  {
677  if(fabs(mod)!=i) continue;
678  StEmcModule* module=EMCdetector->module(i);
679  StSPtrVecEmcRawHit& hits=module->hits();
680  for(unsigned int k=0;k<hits.size();k++) if(hits[k]){
681  unsigned int module=hits[k]->module();
682  unsigned int Eta=hits[k]->eta();
683 // float energyT=hits[k]->energy();
684  int s=fabs(hits[k]->sub());
685  int did(0);
686  if (module==fabs(mod) && Eta == fabs(eta)){
687  float energyT1=hits[k]->energy();
688  mGeom2->getId(module,Eta,s,did);
689  tr2_bemcModule = module;
690  tr2_bemcEtabin = Eta;
691  tr2_bemcEtaValue = position.pseudoRapidity();
692  tr2_bemcPhiValue = position.phi();
693  tr2_bemcSub = s;
694  tr2_bemcEnergy = energyT1;
695 
696  }
697  }
698  }
699  }
700  //
701  //repeat extrapolation, this time to shower max detectors
702  //
703  StEmcDetector* SMDdetector = emcStEvent->detector(kBarrelSmdEtaStripId);
704  if (!SMDdetector)
705  {cout<<"There is no kBarrelSmdEtaStripId Detector ---------"<<endl;}
706  StEmcPosition * projectSmde = new StEmcPosition();
707  StEmcGeom * mGeomSmde=StEmcGeom::instance("bsmde");
708  Bool_t okSmd=projectSmde->trackOnEmc(&position,&momentum,muTrack1,bFld);
709  if(okSmd) {
710  mGeomSmde->getBin(position.phi(),position.pseudoRapidity(),mod,eta,sub);
711  //
712  //see if SMDe has energy
713  //
714  for(unsigned int i=1;i<=120;i++)
715  {
716  if(fabs(mod)!=i) continue;
717  StEmcModule* module=SMDdetector->module(i);
718  StSPtrVecEmcRawHit& hits=module->hits();
719  for(unsigned int k=0;k<hits.size();k++) if(hits[k]){
720  unsigned int module=hits[k]->module();
721  unsigned int Eta=hits[k]->eta();
722 // float energyT=hits[k]->energy();
723  int s=fabs(hits[k]->sub());
724  int did(0);
725  if (module==fabs(mod) && Eta == fabs(eta)){
726  float energyT1=hits[k]->energy();
727  mGeomSmde->getId(module,Eta,s,did);
728 
729  tr1_bsmdeModule = module;
730  tr1_bsmdeEtabin = Eta;
731  tr1_bsmdeEtaValue = position.pseudoRapidity();
732  tr1_bsmdePhiValue = position.phi();
733  tr1_bsmdeSub = s;
734  tr1_bsmdeEnergy = energyT1;
735 
736  }
737  }
738  }
739  }
740  fill ( primaryFlag, &(event->eventSummary()),
741  p1, h1, charge1, p2, h2, charge2, vtx ) ;
742 
743 #ifndef __CINT__
744  tr1.set(1,muTrack1, event); // 1=primary
745  tr2.set(1,muTrack2, event);
746 #endif
747  tr1_timeOfFlight = mBTofPidTraits_1.timeOfFlight();
748  tr1_pathLength = mBTofPidTraits_1.pathLength();
749  tr1_Beta = mBTofPidTraits_1.beta();
750  tr2_timeOfFlight = mBTofPidTraits_2.timeOfFlight();
751  tr2_pathLength = mBTofPidTraits_2.pathLength();
752  tr2_Beta = mBTofPidTraits_2.beta();
753  //
754  //extrapolate tracks to TOF
755  //
756 
757  vector<Int_t> idVec;
758  vector<Double_t> pathVec;
759  PointVec crossVec;
760 
761 
762 
763 
764 
765  if(pairTOFgeoLocal->HelixCrossCellIds(hOuter1,idVec,pathVec,crossVec))
766  {
767 
768  Int_t cellId = -999;
769  Int_t moduleId = -999;
770  Int_t trayId = -999;
771  pairTOFgeoLocal->DecodeCellId(idVec[0],cellId,moduleId,trayId);
772  tr1_extrapolatedTOF_mX = crossVec[0].x();
773  tr1_extrapolatedTOF_mY = crossVec[0].y();
774  tr1_extrapolatedTOF_mZ = crossVec[0].z();
775  cout<<" extrapolated first track"<<endl;
776  }
777  if(pairTOFgeoLocal->HelixCrossCellIds(hOuter2,idVec,pathVec,crossVec))
778  {
779  Int_t cellId = -999;
780  Int_t moduleId = -999;
781  Int_t trayId = -999;
782  pairTOFgeoLocal->DecodeCellId(idVec[0],cellId,moduleId,trayId);
783  tr2_extrapolatedTOF_mX = crossVec[0].x();
784  tr2_extrapolatedTOF_mY = crossVec[0].y();
785  tr2_extrapolatedTOF_mZ = crossVec[0].z();
786  }
787 
788  return 0 ;
789 
790 }
791 Int_t StPeCPair::fill ( Bool_t primaryFlag, StEvent* event ) {
792  //
793  pCharge = track1->geometry()->charge()+track2->geometry()->charge();
794 
795  StThreeVectorF p1 ;
796  StThreeVectorF p2 ;
797  StPhysicalHelixD h1 ;
798  StPhysicalHelixD h2 ;
799  short charge1, charge2 ;
800  //
801  p1 = track1->geometry()->momentum();
802  p2 = track2->geometry()->momentum();
803  charge1 = track1->geometry()->charge();
804  charge2 = track2->geometry()->charge();
805  // if ( !primaryFlag ) {
806 
807  h1 = track1->geometry()->helix() ;
808  h2 = track2->geometry()->helix() ;
809 
810  // }
811 
812  StEventSummary* summary = 0 ;
813  StPrimaryVertex* vtx = 0;
814  vtx = event->primaryVertex();
815  summary = event->summary();
816 
817  StThreeVectorF vtxP ;
818  //
819  // If there is no primary vertex assume (0,0,0)
820  //
821  if ( vtx ) {
822  vtxP = vtx->position() ;
823  } else {
824  vtxP.setX(0.);
825  vtxP.setY(0.);
826  vtxP.setZ(0.);
827  }
828 
829  fill ( primaryFlag, summary, p1, h1, charge1, p2, h2, charge2, vtxP ) ;
830 
831  // fill local track class
832 #ifndef __CINT__
833 
834  tr1.set(1,track1); // 1=primary
835  tr2.set(1,track2);
836 
837 #endif
838  return 0 ;
839 }
840 
841 
842 Int_t StPeCPair::getSumCharge() const{
843  return pCharge ;
844 }
845 
846 Float_t StPeCPair::getSumPt() const{
847  return pPt ;
848 }
849 
850 Float_t StPeCPair::getSumPz() const{
851  return pPz ;
852 }
853 
854 Float_t StPeCPair::getMInv(StPeCSpecies pid) const{
855  if ( pid == 0 ) return pionH.mInv ;
856  else if ( pid == 1 ) return kaonH.mInv ;
857  else if ( pid == 2 ) return protonH.mInv ;
858  else if ( pid == 3 ) return electronH.mInv ;
859  else if ( pid == 4 ) return muonH.mInv ;
860  else {
861  printf ( "StPeCPair::getMInv: wrong pid %d \n", pid ) ;
862  return 0 ;
863  }
864 }
865 
866 Float_t StPeCPair::getOpeningAngle() const{
867  return pAngle ;
868 }
869 Float_t StPeCPair::getCosThetaStar(StPeCSpecies pid) const{
870  if ( pid == 0 ) return pionH.cosThetaStar ;
871  else if ( pid == 1 ) return kaonH.cosThetaStar ;
872  else if ( pid == 2 ) return protonH.cosThetaStar ;
873  else if ( pid == 3 ) return electronH.cosThetaStar ;
874  else if ( pid == 4 ) return muonH.cosThetaStar ;
875  else {
876  printf ( "StPeCPair::getCosThetaStar: wrong pid %d \n", pid ) ;
877  return 0 ;
878  }
879 }
880 
881 
StEmcPosition::trackOnEmc
Bool_t trackOnEmc(StThreeVectorD *position, StThreeVectorD *momentum, const StTrack *const track, Double_t magField, Double_t emcRadius=225.405) const
Track projection utility.
Definition: StEmcPosition.cxx:270
StBTofGeometry::HelixCrossCellIds
Bool_t HelixCrossCellIds(const StHelixD &helix, IntVec &idVec, DoubleVec &pathVec, PointVec &crossVec) const
Definition: StBTofGeometry.cxx:1531
StHelix::pathLengths
pair< double, double > pathLengths(const StHelix &, double minStepSize=10 *micrometer, double minRange=10 *centimeter) const
path lengths at dca between two helices
Definition: StHelix.cc:511
StMuTrack::helix
StPhysicalHelixD helix() const
Returns inner helix (first measured point)
Definition: StMuTrack.cxx:407
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
StMuBTofPidTraits::timeOfFlight
float timeOfFlight() const
timing for PID
Definition: StMuBTofPidTraits.h:119
StEmcGeom::getBin
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Definition: StEmcGeom.h:321
StMuTrack::outerHelix
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
Generated by   doxygen 1.8.5