StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPeCEvent.cxx
1 //
3 // $Id: StPeCEvent.cxx,v 1.26 2016/02/24 17:14:56 ramdebbe Exp $
4 // $Log: StPeCEvent.cxx,v $
5 // Revision 1.26 2016/02/24 17:14:56 ramdebbe
6 // cleaned up old commented line
7 //
8 // Revision 1.25 2015/07/22 20:11:08 ramdebbe
9 // removed RP container copy to output, needs to come back
10 //
11 // Revision 1.24 2015/02/25 01:34:22 ramdebbe
12 // added copy of the Roman Pot StMuRpsCollection to output tree
13 //
14 // Revision 1.23 2014/12/23 20:44:57 ramdebbe
15 // Modified intantiation of StPeCPair when using MuDst input to get TOF extrapolation in the pPairs branch
16 //
17 // Revision 1.22 2014/06/18 16:30:55 ramdebbe
18 // added more variables to event summary
19 //
20 // Revision 1.21 2013/12/27 16:49:50 ramdebbe
21 // added a set method setTOFgeometry to pass pointer to StBTofGeometry
22 //
23 // Revision 1.20 2013/01/24 15:41:16 ramdebbe
24 // added more flags to choose input or output tracks tof etc.
25 //
26 // Revision 1.19 2012/07/08 00:40:20 ramdebbe
27 // initialized variables per Gene VB advice
28 //
29 // Revision 1.18 2012/06/13 15:44:46 ramdebbe
30 // Added flags to include TOF and Vertex branches in tree
31 //
32 // Revision 1.17 2005/08/24 20:58:00 jeromel
33 // TClone to TObj fix
34 //
35 // Revision 1.16 2003/11/25 01:54:26 meissner
36 // correct several bugs: eta cut for tracks, charge sorting, add counting of FTPC and TPC primary tracks, Add bbc information
37 //
38 // Revision 1.15 2003/09/02 17:58:46 perev
39 // gcc 3.2 updates + WarnOff
40 //
41 // Revision 1.14 2003/03/20 20:10:58 yepes
42 // double counting of tracks corrected
43 //
44 // Revision 1.13 2003/03/18 21:20:28 yepes
45 // correcting problem with bField
46 //
47 // Revision 1.12 2003/02/11 20:45:39 yepes
48 // Events without a vertex in MuDst should vertex parameters to -9999
49 //
50 // Revision 1.11 2003/02/05 17:14:05 yepes
51 // Adding bField and pPairs.psi to tree
52 //
53 // Revision 1.10 2002/12/19 18:09:53 yepes
54 // MuDST input added
55 //
56 // Revision 1.9 2002/03/19 22:23:26 meissner
57 // New variables: zdc unatt., Trigger word, MC tree if Geant Branch, DCA for primary pairs, all tracks for secondary pairs (Test)
58 //
59 // Revision 1.8 2001/04/23 21:44:30 meissner
60 // add dEdx z variable to tree, setFormat(1) for tree, use private BetheBloch (temp solution)
61 //
62 // Revision 1.7 2001/02/21 20:41:58 yepes
63 // Add ctb signals to tree
64 //
65 // Revision 1.6 2001/02/13 17:54:40 yepes
66 // still problems on differnt platforms
67 //
68 // Revision 1.5 2001/02/12 22:33:51 yepes
69 // avoid printing
70 //
71 // Revision 1.4 2001/02/12 21:15:42 yepes
72 // New version of StPeCMaker, lots of changes
73 //
74 // Revision 1.3 2000/04/24 19:15:27 nystrand
75 // Fix of a possible memory leak
76 //
77 // Revision 1.2 2000/04/21 19:10:30 nystrand
78 // Include StPeCPair class
79 //
80 // Revision 1.1 2000/03/24 22:37:06 nystrand
81 // First version of StPeCEvent
82 //
83 // Revision 1.0 2000/03/20 23:28:50 nystrand
84 //
86 // TODO:
87 // Understand what the redifinition of the TClone Array does
88 // -> dereferencing the pointer to get just the array
89 // Place the cut on track->flag() in a unique place (Internal track array?)
90 // Understand the dublication of StPeCTrack (does a tree work with pointers)
91 // ///////////////////////////////////////////////////////////////////
92 #include <Stiostream.h>
93 #include "StPeCEvent.h"
94 #include "StPeCMaker.h"
95 #include "StEventTypes.h"
96 #include "StPeCEnumerations.h"
97 #include "StMuDSTMaker/COMMON/StMuEvent.h"
98 #include "StMuDSTMaker/COMMON/StMuTrack.h"
99 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
100 #include "StMuDSTMaker/COMMON/StMuDst.h"
101 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
102 #include "StMuDSTMaker/COMMON/StMuDebug.h"
103 #include "StPreEclMaker/StPreEclMaker.h"
104 #include "StEmcUtil/geometry/StEmcGeom.h"
105 #include "StEmcUtil/others/emcDetectorName.h"
106 #include "StMuDSTMaker/COMMON/StMuDst2StEventMaker.h"
107 #include "St_db_Maker/St_db_Maker.h"
108 #include "StEpcMaker/StEpcMaker.h"
109 #include "StEmcUtil/projection/StEmcPosition.h"
110 #include "StTofCollection.h"
111 #include "StBTofCollection.h"
112 #include "StMuDSTMaker/COMMON/StMuBTofHitCollection.h"
113 #include "TLine.h"
114 #include "TCanvas.h"
115 #include "THelix.h"
116 #include "TView.h"
117 #include "TPolyMarker3D.h"
118 #include "StMuDSTMaker/COMMON/StMuRpsCollection.h"
119 
120 ClassImp(StPeCEvent)
121 
122  StPeCEvent::StPeCEvent(bool useBemc, bool useTOF, bool useVertex, bool useTracks, bool useRP, bool readStMuDst, bool readStEvent, bool readBothInputs) :
123  muDst(0), eventP(0), treecalo(0), tofHits(0), tofTracks(0), vertices(0) {
124  // StPeCMaxTracks ..... looks more like Fortran......
125  pPairs = new TClonesArray ("StPeCPair", StPeCnMaxTracks);
126  sPairs = new TClonesArray ("StPeCPair", StPeCnMaxTracks);
127  tofHits = new TClonesArray ("StMuBTofHit",5*StPeCnMaxTracks);
128  if((useTOFlocal = useTOF)){
129 
130  tofTracks = new TClonesArray ("StMuTrack",5*StPeCnMaxTracks);
131 
132  LOG_INFO << "StPeCEvent constructor: useTOF ---------- " <<useTOF << endm;
133  }
134  if((useBemcLocal = useBemc)) {
135  treecalo = new TClonesArray ("StEmcCluster",5*StPeCnMaxTracks);
136  LOG_INFO << "StPeCEvent constructor: useBemc ---------- " <<useBemc << endm;
137 
138  }
139  if((useVertexLocal = useVertex)){
140  vertices = new TClonesArray ("StMuPrimaryVertex",5*StPeCnMaxTracks);
141  LOG_INFO << "StPeCEvent constructor: useVertex ---------- " <<useVertex << endm;
142 
143  }
144  if((useTracksLocal = useTracks)){
145 
146  LOG_INFO << "StPeCEvent constructor: useTracks ---------- " <<useTracks << endm;
147 // tracks = new TClonesArray ("StPeCTrack",StPeCnMaxTracks);
148  tracks = new TClonesArray ("StMuTrack", StPeCnMaxTracks);
149  }
150  if((useRPLocal = useRP)){
151 
152 // romanPots = new TClonesArray ("StMuRpsCollection", 3);
153  LOG_INFO << "StPeCEvent constructor: useRP ---------- " <<useRP << endm;
154 
155  }
156  if((readStMuDstLocal = readStMuDst)){
157 
158 // LOG_INFO << "StPeCEvent constructor: readStMuDst ---------- " <<readStMuDst << endm;
159 
160  } if((readStEventLocal = readStEvent)){
161 
162 // LOG_INFO << "StPeCEvent constructor: readStEvent ---------- " <<readStEvent << endm;
163 
164  } if((readBothInputsLocal = readBothInputs)){
165 
166  LOG_INFO << "StPeCEvent constructor: both inputs ---------- " <<readBothInputs << endm;
167 
168  }
169 
170 
171 
172  nPPairs = 0 ;
173  nSPairs = 0 ;
174  nTracks = 0 ;
175  eventP = 0 ;
176  shotCount = 0;
177  eventN = 0;
178  runN = 0;
179  nTot = 0;
180  nPrim = 0;
181  qTot = 0;
182  nGlobalTracks = 0;
183  nPrimaryTracks = 0;
184  nPrimaryTPC = 0;
185  nPrimaryFTPC = 0;
186  bField = 0;
187  pt = 0;
188  xVertex = 0;
189  yVertex = 0;
190  zVertex = 0;
191  rVertex = 0;
192  infoLevel = 0;
193  nTOFhitsSum = 0;
194  nBtofTriggerHitsSum = 0;
195  nTOFtracksSum = 0;
196  zdcEastUASum = 0;
197  zdcWestUASum = 0;
198  zdcCoincidenceRateSum = 0;
199  lastDSM0Sum = 0;
200  lastDSM1Sum = 0;
201 
202 
203  LOG_INFO << "StPeCEvent constructor: leaving constructor ---------- " << endm;
204 
205 }
206 
207 StPeCEvent::~StPeCEvent() {
208  clear() ;
209  delete pPairs ;
210  delete sPairs ;
211  delete tracks ;
212  delete tofHits ;
213  delete tofTracks;
214  delete treecalo;
215  delete vertices;
216 // delete romanPots;
217 }
218 
219 
220 void StPeCEvent::clear ( ) {
221  eventN = 0 ;
222  runN = 0 ;
223  nTot = 0 ;
224  nPrim = 0 ;
225  qTot = 0 ;
226  pt = 0 ;
227  xVertex = 0 ;
228  yVertex = 0 ;
229  zVertex = 0 ;
230  rVertex = 0 ;
231  nPPairs = 0 ;
232  nSPairs = 0 ;
233  nTracks = 0 ;
234 
235  nTOFhitsSum = 0;
236  nBtofTriggerHitsSum = 0;
237  nTOFtracksSum = 0;
238  zdcEastUASum = 0;
239  zdcWestUASum = 0;
240  zdcCoincidenceRateSum = 0;
241  lastDSM0Sum = 0;
242  lastDSM1Sum = 0;
243 
244  pPairs->Clear();
245  sPairs->Clear();
246 
247  if(useTracksLocal)tracks->Clear();
248  tofHits->Clear();
249  if(useTOFlocal) {
250 
251  tofTracks->Clear();
252  }
253  if(useBemcLocal)treecalo->Clear();
254  if(useVertexLocal)vertices->Clear();
255 }
256 
257 // The counting of FTPC tracks is not yet debugged
258 Int_t StPeCEvent::fill ( StEvent *event ) {
259 
260  eventP = event ;
261 
262 
263  // I try to get StEvent
264 
265  // check if there is a collection
266  StEmcCollection *emcStEvent = eventP->emcCollection();
267 
268 
269  if(emcStEvent){
270  StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
271  StEmcDetector* detector = emcStEvent->detector(id);
272  if(detector) {
273  StEmcClusterCollection* coll = detector->cluster();
274  if(coll)
275  {
276 
277  StSPtrVecEmcCluster& clusters = coll->clusters();
278  Int_t n = clusters.size();
279  cout<<" number of clusters in BEMC "<<n<<endl;
280  for(Int_t j = 0;j<n;j++)
281  {
282  StEmcCluster *c =clusters[j];
283  if(c)
284  {
285 
286  new((*treecalo)[j]) StEmcCluster((const StEmcCluster &) *c);
287 
288  cout<<" Number of hits in cluster "<<c->nHits()<<endl;
289  cout<<" Energy of cluster "<<c->energy()<<endl;
290  cout<<" eta and phi of cluster "<<c->eta()<<" "<<c->phi()<<endl;
291  }
292  }
293  }
294 
295  }
296 
297  }
298 
299  // Set Run and Event Number
300  eventN = event->id() ;
301  runN = event->runId();
302  cout << "StEvent Run ID: " << runN << endl;
303  cout << "StEvent ID: " << eventN << endl;
304 
305  bField = event->summary()->magneticField();
306 
307  //get trigger info
308 
309  Int_t NGlobal=0;
310  Int_t NPrimaries=0;
311  Int_t SumQ=0;
312  Float_t SumPx=0.0;
313  Float_t SumPy=0.0;
314  // separate TPC and FTPC
315  nPrimaryTPC=0;
316  nPrimaryFTPC=0;
317 
318  // Get the track nodes
319  StSPtrVecTrackNode& exnode = event->trackNodes();
320  Int_t nnode=exnode.size();
321 
322 
323  // local pointer array for track needed for cuts
324  for( Int_t in=0; in<nnode; in++ ) {
325  UInt_t nprim = exnode[in]->entries(primary);
326  UInt_t nglob = exnode[in]->entries(global);
327 
328  if( nprim>1 || nglob>1 ){
329  cout<<"There could be a problem! nprim= "<<nprim<<" nglob= "<<nglob<<endl;
330 
331  }
332  if( nprim==1 ){
333  StTrack *tp = exnode[in]->track(primary);
334  // DANGER THAT SHOULD BE DONE IN A PROPPER WAY, CHECK 4 PLACES !
335  if (! (tp->flag()>0)) continue;
336  // -----------------------------------------------------------
337  NPrimaries++;
338  float px = tp->geometry()->momentum().x();
339  float py = tp->geometry()->momentum().y();
340  SumPx = SumPx + px; SumPy = SumPy + py;
341  SumQ = SumQ + tp->geometry()->charge();
342  //
343  // need to dereference the Array pointer first
344  //new((*tracks)[nTracks++]) StPeCTrack(1,tp) ;
345 
346  if(fabs(tp->geometry()->momentum().pseudoRapidity())<2.0) {
347  nPrimaryTPC++;
348  } else {
349  nPrimaryFTPC++;
350  }
351  }
352  if( nglob==1 ){
353  StTrack *tnp = exnode[in]->track(global);
354  // DANGER THAT SHOULD BE DONE IN A PROPPER WAY, CHECK 4 PLACES !
355 
356  if (! (tnp->flag()>0)) continue;
357 // cout<<" in StPeCEvent after tnp->flag "<< tnp->detectorInfo()->numberOfPoints() <<endl;
358 // if ( tnp->detectorInfo()->numberOfPoints() < 11 ) continue ;
359 // NGlobal++;
360 // cout<<" in StPeCEvent Nglobal "<<NGlobal<<endl;
361 // Do not store tracks in tree any longer PPY 11/5/02
362 // new((*tracks)[nTracks++]) StPeCTrack(0,tnp) ;
363  }
364  }
365  // nGlobals not a good UPC criteria any more
366 
367  // if ( NGlobal > StPeCnMaxTracks ) return 1 ;
368 
369 
370  if (( nPrimaryTPC > 0 || nPrimaryFTPC>0 ) && // at least one track in either FTPC or TPC
371  ( nPrimaryTPC < StPeCnMaxTracks && nPrimaryFTPC< StPeCnMaxTracks ) &&( nPrimaryFTPC+nPrimaryTPC>=2 ) ) {
372  cout << " analyze event !" << endl;
373  } else {
374  cout << " reject event !" << endl;
375  return 1;
376  }
377 
378  nPrim = NPrimaries ;
379  nTot = NGlobal ;
380  qTot = SumQ ;
381  pt = ::sqrt( SumPx*SumPx + SumPy*SumPy );
382 
383  nGlobalTracks= NPrimaries;
384  nPrimaryTracks= NGlobal;
385 
386 
387 
388  StPrimaryVertex* vtx = event->primaryVertex();
389  if(vtx) {
390  cout << "Vertex flag " << vtx->flag() << endl;
391  xVertex = vtx->position().x();
392  yVertex = vtx->position().y();
393  zVertex = vtx->position().z();
394  rVertex = ::sqrt(xVertex*xVertex + yVertex*yVertex);
395 
396  if ( infoLevel > 1 ) {
397  cout << "StPeCEvent : primary vertex x:" << xVertex << " y: " << yVertex << " z: " << zVertex << " r: " << rVertex <<endl;
398  }
399  } else {
400  xVertex = -9999. ;
401  yVertex = -9999. ;
402  zVertex = -9999. ;
403  rVertex = -9999. ;
404  cout<<"StPeCEvent: There was no primary vertex!"<<endl;
405  }
406 
407  // HERE flag must be tested.
408  StPeCPair* lPair ;
409  nPPairs = 0 ;
410  StTrack *trk1, *trk2 ;
411  for( Int_t i1=0; i1<nnode-1; i1++ ) {
412  if( exnode[i1]->entries(primary) !=1 ) continue ;
413  trk1 = exnode[i1]->track(primary);
414 
415  if ( trk1->detectorInfo()->numberOfPoints() < 11 ) continue ;
416  for( Int_t i2=i1+1; i2<nnode; i2++ ) {
417  if( exnode[i2]->entries(primary) !=1 ) continue ;
418  trk2 = exnode[i2]->track(primary);
419 
420  if ( trk2->detectorInfo()->numberOfPoints() < 11 ) continue ;
421  // DANGER
422  if (! (trk1->flag()>0)) continue;
423  if (! (trk2->flag()>0)) continue;
424  // --------
425  // get pointer to memebr ?
426  // TClonesArray &ppairs = *pPairs;
427  // lPair = new(ppairs[nPPairs++]) StPeCPair(trk1,trk2,1,event) ;
428  lPair = new((*pPairs)[nPPairs++]) StPeCPair(trk1,trk2,1,event) ;
429  cout<<"in StPeCEvent lPair written "<<endl;
430  //#ifdef PECPRINT
431  cout << "StPeCEvent : Primary Pair : "
432  << " sumQ = " << lPair->getSumCharge()
433  << " sumPt = " << lPair->getSumPt()
434  << " mInv = " << lPair->getMInv(pion)
435  << " opening angle = " << lPair->getOpeningAngle()
436  << " cos(theta*) = " << lPair->getCosThetaStar(pion) << endl;
437  //#endif
438  }
439  }
440  //removed secondary pairs:
441 
442  return 0 ;
443 }
444 
445 //===========================================================================================
446 Int_t StPeCEvent::fill(StMuDst *mudst) {
447 
448  TObjArray* muTracks = 0;
449  StMuEvent* event = 0;
450  StMuTrack *tp = 0;
451  Bool_t acceptEvent;
452 
453  //Save the event reference
454  muDst = mudst;
455  event = muDst->event();
456 
457  //Set Run and Event Number
458  eventN = event->eventInfo().id();
459 
460  runN = event->eventInfo().runId();
461  LOG_INFO << "StPeCEvent fill(muDst ): useBemc ---------- " <<useBemcLocal << endm;
462  LOG_INFO << "StPeCEvent fill(muDst ): useTOF ---------- " <<useTOFlocal << endm;
463  LOG_INFO << "StPeCEvent fill(muDst ): useVertex ---------- " <<useVertexLocal << endm;
464  LOG_INFO << "StPeCEvent fill(muDst ): useTracks ---------- " <<useTracksLocal << endm;
465  LOG_INFO << "StPeCEvent fill(muDst ): useRP ---------- " <<useRPLocal << endm;
466  LOG_INFO << "StPeCEvent fill(muDst) : TOF geometry pointer ---------- " << mTOFgeoEv<<endm;
467  LOG_INFO << "StMuEvent Run ID: " << runN << endm;
468  LOG_INFO << "StMuEvent ID: " << eventN << endm;
469 
470  bField = event->eventSummary().magneticField();
471  if (fabs(bField)<0.01) LOG_INFO << "StPeCEvent fill(muDst ): BField off " << endm;;
472  acceptEvent = kFALSE;
473  // RD 11-July 2013 to get BEMC data need to find it out of StMuDst (TO DO)
474  // Get EMC calorimeter clusters from StEvent RD
475 
476  // check if there is a collection
477  StMuEmcCollection *emcStEvent = mudst->muEmcCollection();
478 
479  if(emcStEvent){
480 
481  Int_t n = emcStEvent->getNClusters(1);//enum {bemc=1, bprs=2, bsmde=3, bsmdp=4, eemc=5, eprs=6, esmdu=7, esmdv=8};
482 
483  LOG_INFO<<"StPeCEvent::fill(MuDst) number of clusters in BEMC "<<n<<endm;
484  for(Int_t j = 0;j<n;j++)
485  {
486  StMuEmcCluster *c = emcStEvent->getCluster(j, 1);
487  if(c)
488  {
489  if(useBemcLocal){
490  new((*treecalo)[j]) StMuEmcCluster((const StMuEmcCluster &) *c);
491  }
492 
493  }
494  }
495 
496  } //if emc
497  //
498  //here we transfer TOF information in StMuEvent to UPC ntuple RD
499  //
500 
501  LOG_INFO <<"StPeCEvent::fill number of btof hits "<<mudst->numberOfBTofHit()<<endm;
502  int nMax = mudst->numberOfBTofHit();
503  int globalTrackCounter = 0;
504  for(int i=0;i<nMax;i++) {
505  StMuBTofHit *aHit = (StMuBTofHit *)mudst->btofHit(i);
506  if(aHit)
507  {
508  new((*tofHits)[i]) StMuBTofHit((const StMuBTofHit &) *aHit);
509  if(useTOFlocal) {
510 
511  // global track that matched the hit
512  int trayId = aHit->tray();
513  if(trayId<=120&&trayId>=0) {//TOF
514  StMuTrack *TofGlobalTrack = aHit->globalTrack();
515  if(!TofGlobalTrack) continue;
516  new((*tofTracks)[globalTrackCounter]) StMuTrack((const StMuTrack &) *TofGlobalTrack);
517  globalTrackCounter++;
518  }
519  }
520 
521  }
522  }
523 
524  //
525  // write RP collection to output tree
526  StMuRpsCollection *rp = (StMuRpsCollection*)mudst->RpsCollection();
527  if(!rp) LOG_INFO << "StPeCEvent fill(mudst): No Roman Pot collections " << endm;
528 // if(rp){
529 // new((*romanPots)) StMuRpsCollection((const StMuRpsCollection &) *rp);
530 // romanPots = rp;
531 // }
532  // number of vertices not a good number anymore ! FLK 07/03
533  // if ( event->eventSummary().numberOfVertices() ) {
534 
535  StThreeVectorF vtx = event->primaryVertexPosition();
536  // vertex is set to 0,0,0 of no prim vertex found FLK 07/03
537  if(vtx.x() !=0 && vtx.y()!=0 && vtx.z()!=0) {
538  xVertex = vtx.x();
539  yVertex = vtx.y();
540  zVertex = vtx.z();
541  rVertex = ::sqrt(xVertex*xVertex + yVertex*yVertex);
542  }
543  else {
544  xVertex = -9999.;
545  yVertex = -9999.;
546  zVertex = -9999.;
547  rVertex = -9999.;
548  }
549 
550  // directly from StEvent Summary
551  nGlobalTracks= event->eventSummary().numberOfGoodTracks();
552  nPrimaryTracks= event->eventSummary().numberOfGoodPrimaryTracks();
553  nPrimaryTPC=0;
554  nPrimaryFTPC=0;
555 
556  //Retrieve the primary tracks
557  muTracks = muDst->primaryTracks();
558  // nPrim = muTracks->GetLast() + 1; // this ignores quality check FLK
559  // backwards compatibility, FLK
560  nPrim = nPrimaryTracks;
561  nTot = nGlobalTracks;
562  size_t Nvert = muDst->numberOfPrimaryVertices();
563  LOG_INFO << "StPeCEvent::fill(event mudst) #vertices: " <<Nvert<< endm;
564  zdcCoincidenceRateSum = mudst->event()->runInfo().zdcCoincidenceRate();
565  const StTriggerData * trigData;
566  trigData = const_cast< StTriggerData *> (mudst->event()->triggerData());
567 
568  if(!trigData) {
569  LOG_ERROR << "In StPeCEvent summary: StTriggerData not available in StMuDst "<< endm;
570  return 1; //reject event
571  }
572  lastDSM0Sum = trigData->lastDSM(0);
573  lastDSM1Sum = trigData->lastDSM(1);
574  // attenuated signals
575  zdcWestUASum = trigData->zdcAttenuated(west);
576  zdcEastUASum = trigData->zdcAttenuated(east);
577  //TOF hits as seen by trigger (OR of 8 cells)
578  nBtofTriggerHitsSum = trigData->tofMultiplicity();
579  nTOFtracksSum = globalTrackCounter;
580  //
581  // select tracks that match TOF hits and reconstruct vertex RD
582  //
583  // if(Nvert>=1&&nPrimaryTracks<20)matchTOFhitsToTracks(mudst);
584 
585  //
586  //RD as I found the code, it only reads the tracks related to the last vertex
587  //I do not know yet if the best vertex is placed at the end
588  //I will now read all vertices and try all tracks
589  //26-MAR-2010
590 
591 
592  nPPairs = 0 ;
593  for (size_t verti = 0;verti<Nvert;++verti){
594  //LOG_INFO << "StPeCEvent:: vertex Index: "<<verti<<endm;
595  StMuPrimaryVertex* V= muDst->primaryVertex(verti);
596  //
597  if(useVertexLocal){
598  new((*vertices)[verti]) StMuPrimaryVertex((const StMuPrimaryVertex &) *V);
599  }
600  // assert(V);
601  // Float_t rank = V->ranking();
602  LOG_INFO << "StPeCEvent:: vertex index: "<<verti<<endm;
603  StMuDst::setVertexIndex(verti); //assert vertex; this selects tracks connected to this vertex
604  size_t Ntracks = muDst->primaryTracks()->GetEntries();
605 // LOG_INFO << "StPeCEvent::fill(ev mudst) #track in vertex : " <<Ntracks<< endm;
606  nPrimaryTPC = 0;
607  nPrimaryFTPC = 0;
608  for (size_t trackiter = 0;trackiter<Ntracks;trackiter++){
609  tp = ( StMuTrack*)muDst->primaryTracks(trackiter);
610 
611  if (! (tp->flag()>0)) continue; // Quality check on the track
612 
613  if(fabs(tp->eta())<2.0) {
614  nPrimaryTPC++;
615  } else {
616  nPrimaryFTPC++;
617  }
618 
619  // do not fill track list for now to save space
620  // nTracks should be the same as nPrim afterward
621  if(useTracksLocal){
622 // new((*tracks)[nTracks++])StPeCTrack(0, tp); //old scheme
623 // new((*tracks)[nTracks++]) StPeCTrack((const StPeCTrack &) *tp);
624  new((*tracks)[nTracks++]) StMuTrack((const StMuTrack &) *tp);
625 
626  }
627  } // loop over tracks in vertex
628 
629 
630  if (( nPrimaryTPC > 0 || nPrimaryFTPC>0 ) && // at least one track in either FTPC or TPC
631  ( nPrimaryTPC < StPeCnMaxTracks && nPrimaryFTPC< StPeCnMaxTracks )&& (nPrimaryFTPC+nPrimaryTPC>=2)) {
632  LOG_INFO << " analyze vertex !" << endm;
633  acceptEvent = kTRUE;
634  // nPPairs = 0 ;
635  StPeCPair* lPair ;
636  StMuTrack *trk1, *trk2 ;
637  muTracks = muDst->primaryTracks();
638  for(int i1 = 0; i1 <= muTracks->GetLast(); i1++) {
639  trk1 = (StMuTrack*)muTracks->UncheckedAt(i1);
640  for(int i2 = i1+1; i2 <= muTracks->GetLast(); i2++) {
641  trk2 = (StMuTrack*)muTracks->UncheckedAt(i2);
642 
643 
644  // DANGER
645  if (! (trk1->flag()>0)) continue;
646  if (! (trk2->flag()>0)) continue;
647  //LOG_INFO << " wrote the pair !" << endm;
648 
649  lPair = new((*pPairs)[nPPairs++]) StPeCPair(trk1,trk2,1,event, mTOFgeoEv) ;
650 
651 
652  } //lopp mutracks 2
653  } //loop mutracks 1
654  LOG_INFO << " number of pairs " << nPPairs<<" in vertex index "<<verti<<endm;
655  } else { //accept ev
656  LOG_INFO << " reject vertex !" << endm;
657  // return 1;
658  }
659 
660 
661 #ifdef SPAIRS
662  // HERE flag must be tested.
663  nSPairs = 0 ;
664  StMuTrack *muTrk1, *muTrk2 ;
665  muTracks = muDst->globalTracks();
666 
667  for(int i1 = 0; i1 < nTot; i1++) {
668  muTrk1 = (StMuTrack*)muTracks->UncheckedAt(i1);
669  for(int i2 = i1+1; i2 < nTot; i2++) {
670  muTrk2 = (StMuTrack*)muTracks->UncheckedAt(i2);
671  // DANGER
672  if (! (muTrk1->flag()>0)) continue;
673  if (! (muTrk2->flag()>0)) continue;
674  // --------
675  // get pointer to memebr ?
676  // TClonesArray &ppairs = *pPairs;
677  lPair = new((*sPairs)[nSPairs++]) StPeCPair(muTrk1,muTrk2,0,event) ;
678 
679 #ifdef PECPRINT
680  cout << "StPeCEvent : Primary Pair : "
681  << " sumQ = " << lPair->getSumCharge()
682  << " sumPt = " << lPair->getSumPt()
683  << " mInv = " << lPair->getMInv(pion)
684  << " opening angle = " << lPair->getOpeningAngle()
685  << " cos(theta*) = " << lPair->getCosThetaStar(pion) << endl;
686 #endif
687  }
688  }
689 #endif
690 
691  // return 0; RD 4NOV2013
692  } // loop over vertices
693  if(acceptEvent){
694  return 0;
695  }
696  if(!acceptEvent){
697  return 1;
698  }
699  return 0;
700 }
701 
702  std::vector<int> StPeCEvent::fill(StEvent * eventP, StMuDst *mudst) {
703 
704 
705  TObjArray* muTracks = 0;
706  StThreeVectorD position,momentum;
707  StMuEvent* event = 0;
708  StMuTrack *tp = 0;
709  Bool_t acceptEvent;
710  int countAccVtx = 0;
711  int countVtxWithTOF = 0;
712  int countRejectVtxWithTOF = 0;
713  int saveNumTracks = 0;
714  float saveRejectedVtxZ = 0.;
715  float saveAccepteVtxZ = 0.;
716 
717  std::vector<int> vector(15);
718 
719  LOG_INFO << "StPeCEvent fill(muDst StEvent): useBemc ---------- " <<useBemcLocal << endm;
720  LOG_INFO << "StPeCEvent fill(muDst StEvent): useTOF ---------- " <<useTOFlocal << endm;
721  LOG_INFO << "StPeCEvent fill(muDst StEvent): useVertex ---------- " <<useVertexLocal << endm;
722  LOG_INFO << "StPeCEvent fill(muDst StEvent): useTracks ---------- " <<useTracksLocal << endm;
723  LOG_INFO << "StPeCEvent fill(muDst StEvent): useRP ---------- " <<useRPLocal << endm;
724  LOG_INFO << "StPeCEvent fill(muDst StEvent): TOF geometry pointer ---------- " << mTOFgeoEv<<endm;
725  //
726  // Get Magnetic field from event summary
727  //
728  Double_t bFld;
729  StEventSummary* summary = eventP->summary();
730  if(summary) bFld=summary->magneticField()/10.;
731  if (fabs(bFld)<0.01) LOG_INFO << "StPeCEvent fill(muDst StEvent): BField off " << endm;;
732  acceptEvent = kFALSE;
733 
734  // Get EMC calorimeter clusters from StEvent RD
735 
736  // check if there is a collection
737  StEmcCollection *emcStEvent = eventP->emcCollection();
738 
739  if(emcStEvent){
740  StDetectorId id = static_cast<StDetectorId>(0+kBarrelEmcTowerId);
741  StEmcDetector* detector = emcStEvent->detector(id);
742  if(detector) {
743  StEmcClusterCollection* coll = detector->cluster();
744  if(coll)
745  {
746 
747  StSPtrVecEmcCluster& clusters = coll->clusters();
748  Int_t n = clusters.size();
749 // LOG_INFO<<"StPeCEvent::fill number of clusters in BEMC "<<n<<endm;
750  for(Int_t j = 0;j<n;j++)
751  {
752  StEmcCluster *c =clusters[j];
753  if(c)
754  {
755  if(useBemcLocal){
756  new((*treecalo)[j]) StEmcCluster((const StEmcCluster &) *c);
757  }
758 
759  }
760  }
761  }
762 
763  } // if detector
764  } //if emc
765 
766  //
767  //here we transfer TOF information in StMuEvent to UPC ntuple RD
768  //
769 
770 // LOG_INFO <<"StPeCEvent::fill number of btof hits "<<mudst->numberOfBTofHit()<<endm;
771  nTOFhitsSum = mudst->numberOfBTofHit();
772 
773  int nMax = mudst->numberOfBTofHit();
774  int globalTrackCounter = 0;
775  for(int i=0;i<nMax;i++) {
776  StMuBTofHit *aHit = (StMuBTofHit *)mudst->btofHit(i);
777  if(aHit)
778  {
779  new((*tofHits)[i]) StMuBTofHit((const StMuBTofHit &) *aHit);
780  if(useTOFlocal) {
781 
782  // global track that matched the hit
783  int trayId = aHit->tray();
784  if(trayId<=120&&trayId>=0) {//TOF
785  StMuTrack *TofGlobalTrack = aHit->globalTrack();
786  if(!TofGlobalTrack) continue;
787  new((*tofTracks)[globalTrackCounter]) StMuTrack((const StMuTrack &) *TofGlobalTrack);
788  globalTrackCounter++;
789  }
790  }
791 
792  }
793  }
794  // write RP collection to output tree
795  StMuRpsCollection *rp = (StMuRpsCollection*)mudst->RpsCollection();
796  if(!rp) LOG_INFO << "StPeCEvent fill(mudst StEvent): No Roman Pot collections " << endm;
797 // if(rp){
798 // new((*romanPots)) StMuRpsCollection((const StMuRpsCollection &) *rp);
799 // romanPots = rp;
800 // }
801 
802  //Save the event reference
803  muDst = mudst;
804  event = muDst->event();
805  zdcCoincidenceRateSum = mudst->event()->runInfo().zdcCoincidenceRate();
806  const StTriggerData * trigData;
807  trigData = const_cast< StTriggerData *> (mudst->event()->triggerData());
808 
809  if(!trigData) {
810  LOG_ERROR << "In StPeCEvent summary: StTriggerData not available in StMuDst StEvent"<< endm;
811  return vector;
812  }
813  lastDSM0Sum = trigData->lastDSM(0);
814  lastDSM1Sum = trigData->lastDSM(1);
815  // attenuated signals
816  zdcWestUASum = trigData->zdcAttenuated(west);
817  zdcEastUASum = trigData->zdcAttenuated(east);
818  //TOF hits as seen by trigger (OR of 8 cells)
819  nBtofTriggerHitsSum = trigData->tofMultiplicity();
820  nTOFtracksSum = globalTrackCounter;
821  //Set Run and Event Number
822  eventN = event->eventInfo().id();
823 
824  runN = event->eventInfo().runId();
825  LOG_INFO <<"StEvent StMuEvent Run ID: " << runN << " StMuEvent ID: " << eventN << endm;
826 // LOG_INFO << "StMuEvent ID: " << eventN << endm;
827  LOG_INFO << "Found: " << globalTrackCounter <<" tracks with TOF hit match" << endm;
828 
829  bField = event->eventSummary().magneticField();
830 
831  // number of vertices not a good number anymore ! FLK 07/03
832  // if ( event->eventSummary().numberOfVertices() ) {
833  nVertices = event->eventSummary().numberOfVertices();
834  StThreeVectorF vtx = event->primaryVertexPosition();
835  // vertex is set to 0,0,0 of no prim vertex found FLK 07/03
836  if(vtx.x() !=0 && vtx.y()!=0 && vtx.z()!=0) {
837  xVertex = vtx.x();
838  yVertex = vtx.y();
839  zVertex = vtx.z();
840  rVertex = ::sqrt(xVertex*xVertex + yVertex*yVertex);
841  }
842  else {
843  xVertex = -9999.;
844  yVertex = -9999.;
845  zVertex = -9999.;
846  rVertex = -9999.;
847  }
848 
849  // directly from StEvent Summary
850  nGlobalTracks= event->eventSummary().numberOfGoodTracks();
851  nPrimaryTracks= event->eventSummary().numberOfGoodPrimaryTracks();
852  nPrimaryTPC=0;
853  nPrimaryFTPC=0;
854 
855  //Retrieve the primary tracks
856  muTracks = muDst->primaryTracks();
857  // nPrim = muTracks->GetLast() + 1; // this ignores quality check FLK
858  // backwards compatibility, FLK
859  nPrim = nPrimaryTracks;
860  nTot = nGlobalTracks;
861 
862  size_t Nvert = muDst->numberOfPrimaryVertices();
863  LOG_INFO << "StPeCEvent::fill(event mudst) #vertices: " <<Nvert<< endm;
864  //
865  // select tracks that match TOF hits and reconstruct vertex RD
866  //
867  // if(Nvert>=1&&nPrimaryTracks<20)matchTOFhitsToTracks(mudst);
868 
869  //
870  //RD as I found the code, it only reads the tracks related to the last vertex
871  //I do not know yet if the best vertex is placed at the end
872  //I will now read all vertices and try all tracks
873  //26-MRA-2010
874 // for(int i = 0; i <= muTracks->GetLast(); i++) {
875 // tp = (StMuTrack*)muTracks->UncheckedAt(i);
876  //
877  nPPairs = 0 ;
878  for (size_t verti = 0;verti<Nvert;++verti){
879  //LOG_INFO << "StPeCEvent:: vertex Index: "<<verti<<endm;
880  StMuPrimaryVertex* V= muDst->primaryVertex(verti);
881  //
882  if(useVertexLocal){
883  new((*vertices)[verti]) StMuPrimaryVertex((const StMuPrimaryVertex &) *V);
884  }
885  // assert(V);
886  // Float_t rank = V->ranking();
887 
888  LOG_INFO << "StPeCEvent:: loop vertices #BEMC match: " <<V->nBEMCMatch()<<" num tracks used: "<<V->nTracksUsed()<< endm;
889  StMuDst::setVertexIndex(verti);
890  //assert vertex; this selects tracks connected to this vertex
891  size_t Ntracks = muDst->primaryTracks()->GetEntries();
892  LOG_INFO << "StPeCEvent::fill(ev mudst) #track in vertex : " <<Ntracks<<" vertex index "<<verti<< endm;
893  nPrimaryTPC = 0;
894  nPrimaryFTPC = 0;
895  int countTOF = 0;
896  for (int trackiter = 0;trackiter<Ntracks;trackiter++){
897  tp = ( StMuTrack*)muDst->primaryTracks(trackiter);
898  if(tp->index2BTofHit()>-1){
899  LOG_INFO << "StPeCEvent:: a tracks with TOF match -----index to global " <<tp->index2Global()<<" num BTof hits "<<mudst->numberOfBTofHit()<<endm;
900  countTOF++;
901  }
902  if (! (tp->flag()>0)) continue; // Quality check on the track
903 // LOG_INFO << "StPeCEvent:: a tracks key nPrimary loop "<<tp->id()<<endm;
904  if(fabs(tp->eta())<2.0) {
905  nPrimaryTPC++;
906  } else {
907  nPrimaryFTPC++;
908  }
909 
910  // do not fill track list for now to save space
911  // nTracks should be the same as nPrim afterwards
912  if(useTracksLocal){
913 // new((*tracks)[nTracks++])StPeCTrack(0, tp); //old scheme
914 // new((*tracks)[nTracks++]) StPeCTrack((const StPeCTrack &) *tp);
915  new((*tracks)[nTracks++]) StMuTrack((const StMuTrack &) *tp);
916 
917  }
918  }
919  LOG_INFO << "Number of primary TPC tracks: " << nPrimaryTPC << " FTPC tracks " << nPrimaryFTPC<<endm;
920 // LOG_INFO << "Number of primary tracks event summary: " << nPrimaryTracks << " global tracks " << nGlobalTracks<<endm;
921 
922  //
923  // RD 19 APR 2013 comment this and modify to work with UPCpp to make pairs with tracks that have TOF hit match **START**
924  if (( nPrimaryTPC > 0 || nPrimaryFTPC>0 ) && // at least one track in either FTPC or TPC
925  ( nPrimaryTPC < StPeCnMaxTracks && nPrimaryFTPC< StPeCnMaxTracks )&& (nPrimaryFTPC+nPrimaryTPC>=2)) {
926 
927  acceptEvent = kTRUE;
928  countAccVtx++;
929  if(countTOF>0) {
930  countVtxWithTOF++;
931  saveAccepteVtxZ = V->position().z();
932  }
933  LOG_INFO << " analyze vertex ! " <<countAccVtx<< endm;
934  // nPPairs = 0 ;
935  StPeCPair* lPair ;
936  StMuTrack *trk1, *trk2 ;
937  muTracks = muDst->primaryTracks();
938  for(int i1 = 0; i1 <= muTracks->GetLast(); i1++) {
939  trk1 = (StMuTrack*)muTracks->UncheckedAt(i1);
940 // LOG_INFO << "StPeCEvent:: a tracks key lPair loop "<<trk1->id()<<endm;
941  for(int i2 = i1+1; i2 <= muTracks->GetLast(); i2++) {
942  trk2 = (StMuTrack*)muTracks->UncheckedAt(i2);
943 
944  // DANGER
945  if (! (trk1->flag()>0)) continue;
946  if (! (trk2->flag()>0)) continue;
947 
948 
949  lPair = new((*pPairs)[nPPairs++]) StPeCPair(trk1,trk2,1,event, eventP, mTOFgeoEv) ;
950  LOG_INFO << " wrote the pair ! nPPair: " <<nPPairs<< endm;
951 
952  } //lopp mutracks 2
953  } //loop mutracks 1
954  } else { //accept ev
955  LOG_INFO << " reject vertex !" << endm;
956  if(countTOF>0) {
957  countRejectVtxWithTOF++;
958  saveNumTracks = nPrimaryTPC;
959  saveRejectedVtxZ = V->position().z();
960  }
961  // return 1;
962  } // **END**
963 
964 
965  } //verti loop
966 
967 
968 #ifdef SPAIRS
969  // HERE flag must be tested.
970  nSPairs = 0 ;
971  StMuTrack *muTrk1, *muTrk2 ;
972  muTracks = muDst->globalTracks();
973 
974  for(int i1 = 0; i1 < nTot; i1++) {
975  muTrk1 = (StMuTrack*)muTracks->UncheckedAt(i1);
976  for(int i2 = i1+1; i2 < nTot; i2++) {
977  muTrk2 = (StMuTrack*)muTracks->UncheckedAt(i2);
978  // DANGER
979  if (! (muTrk1->flag()>0)) continue;
980  if (! (muTrk2->flag()>0)) continue;
981  // --------
982  // get pointer to member ?
983  // TClonesArray &ppairs = *pPairs;
984  lPair = new((*sPairs)[nSPairs++]) StPeCPair(muTrk1,muTrk2,0,event) ;
985 
986 #ifdef PECPRINT
987  cout << "StPeCEvent : Primary Pair : "
988  << " sumQ = " << lPair->getSumCharge()
989  << " sumPt = " << lPair->getSumPt()
990  << " mInv = " << lPair->getMInv(pion)
991  << " opening angle = " << lPair->getOpeningAngle()
992  << " cos(theta*) = " << lPair->getCosThetaStar(pion) << endl;
993 #endif
994  }
995  }
996 #endif
997  LOG_INFO << "Exiting StPeCEvent::fill(StEvent *event, StMuDst *mudst) acceptEvent " << countAccVtx<<endm;
998 
999  // transfer counts to returning vector
1000  vector[0] = countAccVtx;
1001  vector[1] = countVtxWithTOF;
1002  vector[2] = countRejectVtxWithTOF;
1003  vector[3] = Nvert;
1004  vector[4] = saveNumTracks;
1005  vector[5] = saveRejectedVtxZ;
1006  vector[6] = saveAccepteVtxZ;
1007 
1008  if(acceptEvent){
1009  matchTOFhitsToTracks(mudst);
1010  return vector;
1011  }
1012  if(!acceptEvent){
1013  // matchTOFhitsToTracks(mudst);
1014  return vector;
1015  }
1016  return vector;
1017 }
1018 
1019 StPeCPair* StPeCEvent::getPriPair ( Int_t i ) {
1020  if ( i < 0 || i >= nPPairs ) return 0 ;
1021  TClonesArray &pairs = *pPairs;
1022  return (StPeCPair *)pairs[i] ;
1023 }
1024 
1025 StPeCPair* StPeCEvent::getSecPair ( Int_t i ) {
1026  if ( i < 0 || i >= nSPairs ) return 0 ;
1027  TClonesArray &pairs = *sPairs;
1028  return (StPeCPair *)pairs[i] ;
1029 }
1030 
1031 void StPeCEvent::reset() {
1032  delete pPairs ;
1033  delete sPairs ;
1034  delete tracks ;
1035  pPairs = 0 ;
1036  sPairs = 0 ;
1037  nPPairs = 0 ;
1038  nSPairs = 0 ;
1039  tracks = 0 ;
1040  nTracks = 0 ;
1041 }
1042 
1043 Long_t StPeCEvent::eventNumber() const{ return eventN; }
1044 Long_t StPeCEvent::runNumber() const{ return runN; }
1045 Int_t StPeCEvent::getNTot() const{ return nTot; }
1046 Int_t StPeCEvent::getNPrim() const{ return nPrim; }
1047 Int_t StPeCEvent::getQTot() const{ return qTot; }
1048 Float_t StPeCEvent::getPt() const{ return pt; }
1049 Float_t StPeCEvent::getXVertex() const{ return xVertex; }
1050 Float_t StPeCEvent::getYVertex() const{ return yVertex; }
1051 Float_t StPeCEvent::getZVertex() const{ return zVertex; }
1052 #ifndef __CINT__
1053 
1054 
1055 StLorentzVectorF StPeCEvent::getEvent4Momentum(StPeCSpecies pid) const{
1056  Float_t mptcle=0.0;
1057  if(pid==pion){
1058  mptcle = pion_plus_mass_c2;
1059  }
1060  if(pid==kaon){
1061  mptcle = 493.677*MeV;
1062  }
1063  if(pid==proton){
1064  mptcle = proton_mass_c2;
1065  }
1066  if(pid==electron){
1067  mptcle = electron_mass_c2;
1068  }
1069  if(pid==muon){
1070  mptcle = 105.6584*MeV;
1071  }
1072  StLorentzVectorF p4event(0.0,0.0,0.0,0.0);
1073 
1074  StSPtrVecTrackNode& exnode = eventP->trackNodes();
1075  if ( !eventP ) {
1076  printf ( "StPeCEvent::getEvent4Momentum eventP null \n" ) ;
1077  return p4event ;
1078  }
1079  Int_t nnode=exnode.size();
1080 
1081  for( Int_t in=0; in<nnode; in++ ) {
1082  if( exnode[in]->entries(global) != 1 ) continue ;
1083  StTrack* trk = exnode[in]->track(primary);
1084  StThreeVectorF p = trk->geometry()->momentum();
1085  Float_t energy = p.massHypothesis(mptcle);
1086  StLorentzVectorF pfour(energy,p);
1087  p4event = p4event + pfour;
1088  }
1089 
1090  return p4event;
1091 }
1092 #endif /*__CINT__*/
1093 
1094 Float_t StPeCEvent::mInv(StPeCSpecies pid) const{
1095 
1096  StLorentzVectorF p4event = getEvent4Momentum(pid);
1097 
1098  return p4event.m();
1099 }
1100 
1101 Float_t StPeCEvent::yRap(StPeCSpecies pid) const{
1102 
1103 
1104  StLorentzVectorF p4event = getEvent4Momentum(pid);
1105 
1106  return p4event.rapidity();
1107 }
1108 
1109 void StPeCEvent::matchTOFhitsToTracks(StMuDst *mudst) {
1110 
1111  TObjArray* priTracks = 0;
1112  TObjArray* gloTracks = 0;
1113 // StMuTrack *trk;
1114  StMuEvent* event = 0;
1115  Int_t snapLimit = 0;
1116  Float_t angleHit;
1117 
1118  if (gROOT->IsBatch()) return;
1119 
1120  event = mudst->event();
1121  //Set Run and Event Number
1122  eventN = event->eventInfo().id();
1123 
1124  runN = event->eventInfo().runId();
1125  LOG_INFO << "-------------- got into matchTOFhitsToTracks number of shots "<<shotCount<< endm;
1126  TDirectory * saveDir = gDirectory;
1127  //
1128  //create a canvas to draw the event and then save to pdf file
1129  //
1130 
1131  TCanvas * cSnap;
1132  TCanvas * cSnap3D;
1133 
1134  int nMax = mudst->numberOfBTofHit();
1135  if(shotCount<snapLimit && snapLimit >0) {
1136  //
1137  //create a canvas to draw the event and then save to pdf file
1138  //
1139 
1140  cSnap = new TCanvas("cSnap", "Snapshot", 600, 600);
1141  cSnap3D = new TCanvas("cSnap3D","Snapshot3D", 600, 600);
1142 
1143  gDirectory->pwd();
1144  gDirectory->cd("test1.tree.root:/snapShots/");
1145  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks after directory: " << endm;
1146  gDirectory->pwd();
1147  // gDirectory->ls();
1148  TH2F * hist = (TH2F*)gDirectory->FindObject(Form("snapShot%d",shotCount));
1149  if(!hist)LOG_INFO << "FindObject failed " << endm;
1150  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks name of hist: " << hist->GetName()<<endm;
1151 
1152  StMuTriggerIdCollection tt=mudst->event()->triggerIdCollection();
1153 
1154  const StTriggerId ttid= tt.nominal();
1155  hist->SetTitle(Form("Run: %d Event number %d Main: %d topo: %d #tof: %d",runN, eventN, ttid.isTrigger(400631), ttid.isTrigger(1), nMax));
1156  cSnap->SetName(Form("snapShot%d",shotCount));
1157  cSnap3D->SetName(Form("snapShot%d_3D",shotCount));
1158  cSnap->Draw();
1159  // cSnap3D->Draw();
1160  cSnap->cd();
1161  hist->Draw();
1162  }
1163 
1164 
1165  priTracks = mudst->primaryTracks();
1166  //LOG_INFO << "StPeCEvent::matchTOFhitsToTracks #primary tracks: " <<priTracks->GetEntries()<< endm;
1167  gloTracks = mudst->globalTracks();
1168 
1169 // int nGlobalTracks = gloTracks->GetEntries();
1170 // int nPrimaryTracks = priTracks->GetEntries();
1171  int globalTrackCounter = 0;
1172  //
1173  //loop on vertex tracks to find those that match TOF hits
1174  //
1175  size_t Nvert = mudst->numberOfPrimaryVertices();
1176 
1177  gFile->cd("test1.tree.root:/snapShots");
1178  //gDirectory->ls();
1179 
1180  TH1F * histV = (TH1F*)gDirectory->FindObject("hNumVtx");
1181 
1182  histV->Fill(Nvert);
1183  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks #vertices: " <<Nvert<< endm;
1184  if(shotCount<snapLimit){
1185  cSnap3D->cd();
1186  Double_t xyz[3]={0.,0.,0.};
1187  Double_t v[3]={0.,0.,0.};
1188  Double_t angularFreq = 0.;
1189  Double_t range[2] = {0., 0.3};
1190 
1191  for (size_t verti = 0;verti<Nvert;++verti){
1192  LOG_INFO << "StPeCEvent:: vertex Index: "<<verti<<endm;
1193  StMuDst::setVertexIndex(verti);
1194  size_t Ntracks = mudst->primaryTracks()->GetEntries();
1195  // LOG_INFO << "StPeCEvent::matchTOFhitsToTracks #track in vertex : " <<Ntracks<< endm;
1196  for (size_t trackiter = 0;trackiter<Ntracks;trackiter++){
1197  StMuTrack* track = mudst->primaryTracks(trackiter);
1198  //LOG_INFO << "Y track end-point : " <<track->lastPoint().y()<< " start point "<< track->firstPoint().y()<<endm;
1199  LOG_INFO << "track start-point x : " <<track->firstPoint().x()<< " start point y "<< track->firstPoint().y()<< " start point z "<< track->firstPoint().z()<<endm;
1200  LOG_INFO << "track charge : " <<track->charge()<< " mag field "<< event->eventSummary().magneticField()<<endm;
1201  LOG_INFO << "track momentum x : " <<track->momentum().x()<< " p y "<< track->momentum().y()<< " p z "<< track->momentum().z()<<endm;
1202  if(shotCount<snapLimit) {
1203  TLine * straightTrack = new TLine(track->firstPoint().z(), track->firstPoint().y(), track->lastPoint().z(), track->lastPoint().y());
1204  straightTrack->SetLineColor(verti+1);
1205 
1206 // THelix * myHelix = new THelix(track->helix().origin().x(),track->helix().origin().y(),
1207 // track->helix().origin().z(),
1208 // track->momentum().x(),track->momentum().y(),
1209 // track->momentum().z(),track->charge()*track->helix().h());
1210 // THelix * myHelix = new THelix(track->firstPoint().x(),track->firstPoint().y(),
1211 // 0.,
1212 // track->momentum().x(),track->momentum().y(),
1213 // track->momentum().z(),track->charge()*event->eventSummary().magneticField());
1214  xyz[0] = track->firstPoint().x();
1215  xyz[1] = track->firstPoint().y();
1216  xyz[2] = track->firstPoint().z();
1217  v[0] = track->momentum().x();
1218  v[1] = track->momentum().y();
1219  v[2] = track->momentum().z();
1220  angularFreq = track->charge()*event->eventSummary().magneticField();
1221  THelix * myHelix = new THelix(&xyz[0],&v[0],angularFreq, &range[0],kHelixZ);
1222  gPad->SetFillColor(37);
1223  // myHelix->SetRange(0.,0.1, kHelixZ);
1224  myHelix->Draw();
1225  }
1226 
1227 
1228 
1229  } // loop over primary tracks
1230 
1231  } // loop over vertices
1232  Float_t radTOF = 0.06;
1233  //test polymarker
1234  TPolyMarker3D *pm = new TPolyMarker3D(360, 20);
1235  pm->SetMarkerSize(0.001);
1236  for(int deg=0;deg<360;deg++){
1237  pm->SetPoint(deg, radTOF*TMath::Cos(deg),radTOF*TMath::Sin(deg),0.);
1238  }
1239  pm->Draw();
1240 
1241  TPolyMarker3D *pmTOF = new TPolyMarker3D(60, 20);
1242  pmTOF->SetMarkerSize(0.004);
1243  pmTOF->SetMarkerColor(kBlue);
1244 
1245  for(int i=0;i<nMax;i++) {
1246  StMuBTofHit *aHit = (StMuBTofHit *)mudst->btofHit(i);
1247  if(aHit)
1248  {
1249  // global track that matched the hit
1250  int trayId = aHit->tray();
1251  if(trayId<=120&&trayId>=0) {//TOF
1252  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks tray id: "<<trayId<<" index i "<<i<<endm;
1253  trayId -= 1;
1254  angleHit = (72-trayId*6)*(3.14/180.);
1255  if(trayId>59) {
1256  trayId -= 61;
1257  angleHit = (108+trayId*6)*(3.14/180.);
1258  }
1259  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks angle: "<<angleHit<<endm;
1260  pmTOF->SetPoint(i, radTOF*TMath::Cos(angleHit),radTOF*TMath::Sin(angleHit),0.);
1261  }
1262  }
1263 
1264  }
1265  pmTOF->Draw();
1266 
1267  if(shotCount<snapLimit){
1268  gPad->GetViewer3D("ogl");
1269  cSnap->cd();
1270  }
1271  //
1272  //repeat loops to display yz view
1273  //
1274 
1275  for (size_t verti = 0;verti<Nvert;++verti){
1276 
1277  StMuDst::setVertexIndex(verti);
1278  size_t Ntracks = mudst->primaryTracks()->GetEntries();
1279  // LOG_INFO << "StPeCEvent::matchTOFhitsToTracks #track in vertex : " <<Ntracks<< endm;
1280  for (size_t trackiter = 0;trackiter<Ntracks;trackiter++){
1281  StMuTrack* track = mudst->primaryTracks(trackiter);
1282 
1283  if(shotCount<snapLimit) {
1284  TLine * straightTrack = new TLine(track->firstPoint().z(), track->firstPoint().y(), track->lastPoint().z(), track->lastPoint().y());
1285  straightTrack->SetLineColor(verti+1);
1286  straightTrack->Draw("same");
1287  }
1288 
1289  //for each track loop over tof hits
1290  for(int i=0;i<nMax;i++) {
1291  StMuBTofHit *aHit = (StMuBTofHit *)mudst->btofHit(i);
1292  if(aHit)
1293  {
1294  int trayId = aHit->tray();
1295  if(trayId<=120&&trayId>=0) {//TOF
1296  StMuTrack *TofGlobalTrack = aHit->globalTrack();
1297  if(!TofGlobalTrack) continue;
1298 
1299  //we loop over global tracks and exit loop as match is made
1300  //
1301  int tofHitTrkId = TofGlobalTrack->id();
1302  if(track->id()==tofHitTrkId) {
1303  LOG_INFO << "StPeCEvent::match GOT MATCH: "<<track->id()<<" vertex index: "<<track->vertexIndex()<<endm;
1304  if(shotCount<snapLimit) {
1305  TLine * straightTof = new TLine(track->firstPoint().z(), track->firstPoint().y(), track->lastPoint().z(), track->lastPoint().y());
1306  straightTof->SetLineWidth(3.3);
1307  straightTof->Draw("same");
1308  }
1309  break;
1310  }
1311  // }
1312  globalTrackCounter++;
1313  }
1314 
1315  }
1316  }
1317  } //
1318  } // loop over vertices yz
1319  } //if shotCount...
1320 
1321 
1322  //
1323  //display global tracks
1324  size_t Gtracks = mudst->globalTracks()->GetEntries();
1325  LOG_INFO << "StPeCEvent::matchTOFhitsToTracks #Global tracks : " <<Gtracks<< endm;
1326  for (size_t Gtrackiter = 0;Gtrackiter<Gtracks;Gtrackiter++){
1327  StMuTrack* track = mudst->globalTracks(Gtrackiter);
1328 
1329  if(shotCount<snapLimit) {
1330  TLine * straightTrack = new TLine(track->firstPoint().z(), track->firstPoint().y(), track->lastPoint().z(), track->lastPoint().y());
1331  straightTrack->SetLineColor(2);
1332  // straightTrack->Draw("same");
1333  }
1334  }
1335 
1336  gDirectory = saveDir;
1337  if(shotCount<snapLimit) {
1338  cSnap->SaveAs(Form("snapShot%d.pdf",shotCount));
1339  cSnap->Close();
1340  cSnap3D->Close();
1341  }
1342  shotCount++;
1343 }
1344 
StMuDst::btofHit
static StMuBTofHit * btofHit(int i)
returns pointer to the i-th muBTofHit
Definition: StMuDst.h:419
StMuDst::globalTracks
static TObjArray * globalTracks()
returns pointer to the global tracks list
Definition: StMuDst.h:303
StMuTrack::vertexIndex
Int_t vertexIndex() const
Returns index of associated primary vertex.
Definition: StMuTrack.cxx:600
StMuTrack::index2Global
Int_t index2Global() const
Returns index of associated global track. If not in order can be set with StMuDst::fixTrackIndeces() ...
Definition: StMuTrack.h:231
StMuTrack::id
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
StMuTrack::flag
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StMuTrack::eta
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
StMuTrack::firstPoint
const StThreeVectorF & firstPoint() const
Returns positions of first measured point.
Definition: StMuTrack.h:261
StMuDst::primaryTracks
static TObjArray * primaryTracks()
returns pointer to a list of tracks belonging to the selected primary vertex
Definition: StMuDst.h:301
StMuDst::muEmcCollection
static StMuEmcCollection * muEmcCollection()
returns pointer to current StMuEmcCollection
Definition: StMuDst.h:389
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMuTrack::momentum
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:260
StMuTrack::lastPoint
const StThreeVectorF & lastPoint() const
Returns positions of last measured point.
Definition: StMuTrack.h:262
StPeCEvent
#include &quot;StEventTypes.h&quot;
Definition: StPeCEvent.h:83
StMuDst::setVertexIndex
static void setVertexIndex(Int_t vtx_id)
Set the index number of the current primary vertex (used by both primaryTracks() functions and for St...
Definition: StMuDst.cxx:273
StMuTriggerIdCollection
Collection of trigger ids as stored in MuDst.
Definition: StMuTriggerIdCollection.h:31
Generated by   doxygen 1.8.5