StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2009W_algo.cxx
1 // $Id: St2009W_algo.cxx,v 1.25 2011/09/14 14:23:20 stevens4 Exp $
2 //
3 //*-- Author : Jan Balewski, MIT
4 //*-- Author for Endcap: Justin Stevens, IUCF
5 
6 #include "StEmcUtil/geometry/StEmcGeom.h"
7 #include "WeventDisplay.h"
8 #include "StSpinPool/StJets/StJet.h"
9 #include "StSpinPool/StJets/TowerToJetIndex.h"
10 
11 //muDst needed for zdcRate
12 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
13 #include <StMuDSTMaker/COMMON/StMuDst.h>
14 #include <StMuDSTMaker/COMMON/StMuEvent.h>
15 
16 #include "St2009WMaker.h"
17 
18 //________________________________________________
19 //________________________________________________
20 void
21 St2009WMaker::find_W_boson(){
22 
23  // printf("========= find_W_boson() \n");
24  int nGoldW=0;
25  //remove events tagged as Zs
26  if(wEve.zTag) return;
27 
28  // search for Ws ............
29  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
30  WeveVertex &V=wEve.vertex[iv];
31  for(unsigned int it=0;it<V.eleTrack.size();it++) {
32  WeveEleTrack &T=V.eleTrack[it];
33  if(T.isMatch2Cl==false) continue;
34  assert(T.cluster.nTower>0); // internal logical error
35  assert(T.nearTotET>0); // internal logical error
36 
37  //make cut on lepton |eta| for cross section
38  if(fabs(T.primP.Eta()) > par_leptonEta) continue;
39  hA[0]->Fill("eta1",1.);
40 
41  //signal plots w/o EEMC in veto
42  if(T.cluster.ET/T.nearTotET_noEEMC>par_nearTotEtFrac){
43  if(T.awayTotET_noEEMC < 8)//old awayside pt cut
44  hA[141]->Fill(T.cluster.ET);
45  if(T.sPtBalance_noEEMC>par_ptBalance ) {//only signed ptBalance cut
46  hA[140]->Fill(T.cluster.ET);
47  if (T.prMuTrack->charge() < 0) {
48  hA[184+3]->Fill(T.cluster.ET);
49  hA[200+3]->Fill(T.primP.Eta(),T.cluster.ET);
50  hA[280+3]->Fill(abs(T.primP.Eta()),T.cluster.ET);
51  } else if (T.prMuTrack->charge() > 0) {
52  hA[184+4]->Fill(T.cluster.ET);
53  hA[200+4]->Fill(T.primP.Eta(),T.cluster.ET);
54  hA[280+4]->Fill(abs(T.primP.Eta()),T.cluster.ET);
55  }
56  }
57  }
58 
59  if(T.cluster.ET /T.nearTotET< par_nearTotEtFrac) continue; // too large nearET
60 
61  hA[20]->Fill("noNear",1.);
62  hA[113]->Fill( T.cluster.ET); // for Joe
63  hA[50]->Fill(T.awayTpcPT);
64  hA[51]->Fill(T.awayBtowET);
65  hA[54]->Fill(T.awayTotET);
66  hA[52]->Fill(T.cluster.ET,T.awayTotET);
67  hA[53]->Fill(T.cluster.ET,T.awayEmcET);
68  hA[55]->Fill(T.awayEtowET);
69  hA[60]->Fill(T.cluster.ET,T.awayTpcPT);
70 
71  hA[132]->Fill(T.cluster.ET,T.ptBalance.Perp());
72  hA[133]->Fill(T.awayTotET,T.ptBalance.Perp());
73  hA[134]->Fill(T.cluster.ET,T.sPtBalance);
74  hA[135]->Fill(T.awayTotET,T.sPtBalance);
75 
76  int i=0;
77  //background variations for systematic
78  for (int j=0; j<=20; j++) { //loop over possible sPtBal cuts
79  float pTBal_cut = 5.+((float) j);
80  if (T.sPtBalance<pTBal_cut) { //background
81  if (T.prMuTrack->charge() < 0) {
82  hA[142+i]->Fill(T.cluster.ET,j);
83  hA[210+j]->Fill(T.primP.Eta(),T.cluster.ET);
84  hA[290+j]->Fill(abs(T.primP.Eta()),T.cluster.ET);
85  } else if (T.prMuTrack->charge() > 0) {
86  hA[163+i]->Fill(T.cluster.ET,j);
87  hA[231+j]->Fill(T.primP.Eta(),T.cluster.ET);
88  hA[311+j]->Fill(abs(T.primP.Eta()),T.cluster.ET);
89  }
90  }
91  }
92 
93  //plots for backg sub yield
94  if(T.sPtBalance>par_ptBalance ) {
95  hA[136]->Fill(T.cluster.ET);//signal
96  hA[62]->Fill(T.pointTower.iEta ,T.cluster.energy);
97  if (T.prMuTrack->charge() < 0) {
98  hA[184+1]->Fill(T.cluster.ET);
99  hA[200+1]->Fill(T.primP.Eta(),T.cluster.ET);
100  hA[280+1]->Fill(abs(T.primP.Eta()),T.cluster.ET);
101  } else if (T.prMuTrack->charge() > 0) {
102  hA[184+2]->Fill(T.cluster.ET);
103  hA[200+2]->Fill(T.primP.Eta(),T.cluster.ET);
104  hA[280+2]->Fill(abs(T.primP.Eta()),T.cluster.ET);
105  }
106  } else {
107  hA[137]->Fill(T.cluster.ET);//background
108  if (T.prMuTrack->charge() < 0) {
109  hA[184+5]->Fill(T.cluster.ET);
110  hA[200+5]->Fill(T.primP.Eta(),T.cluster.ET);
111  hA[280+5]->Fill(abs(T.primP.Eta()),T.cluster.ET);
112  } else if (T.prMuTrack->charge() > 0) {
113  hA[184+6]->Fill(T.cluster.ET);
114  hA[200+6]->Fill(T.primP.Eta(),T.cluster.ET);
115  hA[280+6]->Fill(abs(T.primP.Eta()),T.cluster.ET);
116  }
117  }
118 
119  //plots for backg sub yield (old awayTot cut DNP)
120  if(T.awayTotET < 8)
121  hA[138]->Fill(T.cluster.ET);//old signal
122  else
123  hA[139]->Fill(T.cluster.ET);//old background
124 
125  if(0){/***************************/
126  printf("\n WWWWWWWWWWWWWWWWWWWWW\n");
127  wDisaply->exportEvent( "W", V, T);
128  wEve.print();
129  }/***************************/
130 
131 
132  //put final W cut here
133  if(T.sPtBalance<par_ptBalance) continue;
134  //::::::::::::::::::::::::::::::::::::::::::::::::
135  //:::::accepted W events for x-section :::::::::::
136  //::::::::::::::::::::::::::::::::::::::::::::::::
137  wEve.wTag=true;
138 
139  hA[20]->Fill("noAway",1.0);
140  hA[114]->Fill( T.cluster.ET);//for Joe
141 
142  hA[90]->Fill( T.cluster.ET);
143  hA[92]->Fill( T.cluster.ET,T.glMuTrack->dEdx()*1e6);
144  hA[93]->Fill( T.cluster.ET,T.glMuTrack->dca().mag());
145  int k=0; if(T.prMuTrack->charge()<0) k=1;
146  hA[94+k]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
147  // h95 used above
148 
149 
150  if(T.cluster.ET<par_highET) continue; // very likely Ws
151  hA[91]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.position.Phi());
152  hA[96]->Fill(V.id);
153  hA[97]->Fill(V.funnyRank);
154  hA[98]->Fill(V.z);
155  hA[99]->Fill( T.prMuTrack->eta());
156  hA[20]->Fill("goldW",1.);
157  nGoldW++;
158 
159  //some final plots for comparing to embedding
160  float zdcRate=mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
161  hA[260]->Fill(zdcRate,wEve.bx7);
162  hA[261]->Fill(zdcRate,wEve.vertex.size());
163  hA[262]->Fill(zdcRate,T.prMuTrack->nHitsFit());
164  float hitFrac=1.*T.prMuTrack->nHitsFit()/T.prMuTrack->nHitsPoss();
165  hA[263]->Fill(zdcRate,hitFrac);
166  hA[264]->Fill(zdcRate,T.prMuTrack->globalTrack()->lastPoint().perp());
167  hA[265]->Fill(zdcRate,T.prMuTrack->globalTrack()->firstPoint().perp());
168  hA[266]->Fill(zdcRate,T.prMuTrack->pt());
169  hA[267]->Fill(zdcRate,T.prMuTrack->globalTrack()->chi2());
170  hA[268]->Fill(zdcRate,1./T.prMuTrack->pt());
171  hA[269]->Fill(zdcRate,T.glMuTrack->dca().mag());
172  hA[270]->Fill(zdcRate,T.glMuTrack->dcaD());
173  hA[271]->Fill(zdcRate,T.glMuTrack->dcaZ());
174 
175  if(V.funnyRank<10) {
176  hA[272]->Fill(zdcRate,T.glMuTrack->pt());
177  hA[273]->Fill(zdcRate,1./T.glMuTrack->pt());
178  }
179 
180  }// loop over tracks
181  }// loop over vertices
182  if(nGoldW>0)
183  hA[0]->Fill("goldW",1.);
184 
185 }
186 
187 
188 //________________________________________________
189 //________________________________________________
190 void
191 St2009WMaker::tag_Z_boson(){
192 
193  float par_jetPt=10.;
194  float lowMass=70.; float highMass=140.;
195  mJets = getJets("ConeJets12_100"); //select specific jet-type
196 
197  //form invariant mass from lepton candidate and jet
198  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {//vertex loop
199  WeveVertex &V=wEve.vertex[iv];
200  for(unsigned int it=0;it<V.eleTrack.size();it++) {// select track
201  WeveEleTrack &T1=V.eleTrack[it];
202  if(T1.isMatch2Cl==false) continue;
203  assert(T1.cluster.nTower>0); // internal logical error
204  assert(T1.nearTotET>0); // internal logical error
205  if(T1.cluster.ET/T1.nearTotET< par_nearTotEtFrac) continue; // too large nearET
206 
207  //match lepton candidate with jet
208  TLorentzVector jetVec;
209  for (int i_jet=0; i_jet< nJets; i_jet++){//jet loop
210  jetVec = *((StJet*)mJets->At(i_jet));
211  if(jetVec.Pt()<par_jetPt) continue;//remove low pt jets
212 
213  //electron like cut on jets
214  StJet* jet = getJet(i_jet); float maxCluster=0.;
215  int totTowers=jet->nBtowers+jet->nEtowers;
216  for(int itow=0;itow<totTowers;itow++){//loop over towers
217  if(jet->tower(itow)->detectorId()==13)//drop endcap towers
218  continue;
219 
220  int softId=jet->tower(itow)->towerId();
221  //find highest 2x2 BTOW cluster in jet
222  TVector3 pos=positionBtow[softId-1]; int iEta,iPhi;
223  if( L2algoEtaPhi2IJ(pos.Eta(),pos.Phi(),iEta,iPhi))
224  continue;
225  float cluster=maxBtow2x2(iEta,iPhi,jet->zVertex).ET;
226  if(cluster>maxCluster) maxCluster=cluster;
227  }
228 
229  TVector3 jetVec3(jetVec.X(),jetVec.Y(),jetVec.Z());
230  if(jetVec3.DeltaR(T1.primP)<par_nearDeltaR)
231  continue;//skip jets in candidate phi isolation'cone'
232 
233  //form invM
234  float e1=T1.cluster.energy;
235  TVector3 p1=T1.primP; p1.SetMag(e1);
236  TLorentzVector ele1(p1,e1); //lepton candidate 4- mom
237  TLorentzVector sum=ele1+jetVec;
238  float invM=sqrt(sum*sum);
239  if(maxCluster/jet->jetPt < 0.5) continue;
240  if(invM > lowMass && invM < highMass)
241  wEve.zTag=true;
242  }
243  }
244  }
245 
246 }
247 
248 //________________________________________________
249 //________________________________________________
250 void
251 St2009WMaker::findPtBalance(){
252 
253  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
254  WeveVertex &V=wEve.vertex[iv];
255  for(unsigned int it=0;it<V.eleTrack.size();it++) {
256  WeveEleTrack &T=V.eleTrack[it];
257  if(T.isMatch2Cl==false) continue;
258 
259  //****loop over branch with EEMC****
260  mJets = getJets(mJetTreeBranch);
261  int nJetsWE=nJets;
262  for (int i_jet=0; i_jet< nJetsWE; i_jet++){//loop over jets
263  StJet* jet = getJet(i_jet);
264  TVector3 jetVec; //vector for jet momentum
265  //vary neutral and charged et in jets for systematic
266  float neutral=jet->neutralFraction()*jet->Pt();
267  float charged=jet->chargedFraction()*jet->Pt();
268  neutral=neutral*par_mcJetNeutScale;
269  charged=charged*par_mcJetChrgScale;
270  float sum=neutral+charged;
271  jetVec.SetPtEtaPhi(sum,jet->Eta(),jet->Phi());
272  if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
273  T.ptBalance+=jetVec;
274  }
275  TVector3 clustPt(T.primP.X(),T.primP.Y(),0);
276  clustPt.SetMag(T.cluster.ET);
277  T.ptBalance+=clustPt;
278  T.sPtBalance = T.ptBalance.Perp();
279  if(T.ptBalance.Dot(clustPt)<0) T.sPtBalance *=-1.;
280 
281  //****loop over branch without EEMC****
282  mJets = getJets(mJetTreeBranch_noEEMC);
283  int nJetsNE=nJets;
284  TVector3 highJet; //zero out highJet
285  for (int i_jet=0; i_jet< nJetsNE; i_jet++){//loop over jets
286  StJet* jet = getJet(i_jet);
287  TVector3 jetVec; //vector for jet momentum
288  //vary neutral and charged et in jets for systematic
289  float neutral=jet->neutralFraction()*jet->Pt();
290  float charged=jet->chargedFraction()*jet->Pt();
291  neutral=neutral*par_mcJetNeutScale;
292  charged=charged*par_mcJetChrgScale;
293  float sum=neutral+charged;
294  jetVec.SetPtEtaPhi(sum,jet->Eta(),jet->Phi());
295  if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
296  T.ptBalance_noEEMC+=jetVec;
297  }
298  T.ptBalance_noEEMC+=clustPt;
299  T.sPtBalance_noEEMC = T.ptBalance_noEEMC.Perp();
300  if(T.ptBalance_noEEMC.Dot(clustPt)<0) T.sPtBalance_noEEMC *=-1.;
301 
302  }// end of loop over tracks
303  }// end of loop over vertices
304 
305 }
306 
307 
308 //________________________________________________
309 //________________________________________________
310 void
311 St2009WMaker::findAwayJet(){
312  // printf("\n******* find AwayJet() nVert=%d\n",wEve.vertex.size());
313  //wEve.print();
314  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
315  WeveVertex &V=wEve.vertex[iv];
316  for(unsigned int it=0;it<V.eleTrack.size();it++) {
317  WeveEleTrack &T=V.eleTrack[it];
318  if(T.isMatch2Cl==false) continue;
319 
320  // .... sum opposite in phi EMC components
321  T.awayBtowET=sumBtowCone(V.z,-T.primP,1,T.awayNTow); // '1'= only cut on delta phi
322  T.awayEmcET=T.awayBtowET;
323  T.awayEtowET=sumEtowCone(V.z,-T.primP,1,T.awayNTow);
324  if(par_useEtow >= 2) T.awayEmcET+=T.awayEtowET;
325 
326  //..... add TPC ET
327  T.awayTpcPT=sumTpcCone(V.id,-T.primP,1,T.awayNTr);
328  T.awayTotET=T.awayEmcET+T.awayTpcPT;
329  T.awayTotET_noEEMC=T.awayBtowET+T.awayTpcPT;
330  //printf("\n*** in awayTpc=%.1f awayEmc=%.1f\n ",T.awayTpcPT,T.awayEmcET); T.print();
331  }// end of loop over tracks
332  }// end of loop over vertices
333 }
334 
335 //________________________________________________
336 //________________________________________________
337 void
338 St2009WMaker::findNearJet(){
339  //printf("\n******* findNearJet() nVert=%d\n",wEve.vertex.size());
340 
341  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
342  WeveVertex &V=wEve.vertex[iv];
343  for(unsigned int it=0;it<V.eleTrack.size();it++) {
344  WeveEleTrack &T=V.eleTrack[it];
345  if(T.isMatch2Cl==false) continue;
346 
347  // .... sum EMC-jet component
348  T.nearBtowET=sumBtowCone(V.z,T.primP,2,T.nearNTow); // '2'=2D cone
349  T.nearEmcET+=T.nearBtowET;
350  T.nearEtowET=sumEtowCone(V.z,T.primP,2,T.nearNTow);
351  if(par_useEtow >= 3) T.nearEmcET+=T.nearEtowET;
352  // .... sum TPC-near component
353  T.nearTpcPT=sumTpcCone(V.id,T.primP,2,T.nearNTr); // '2'=2D cone
354  hA[47]->Fill(T.nearTpcPT);
355  hA[48]->Fill(T.nearEmcET,T.nearTpcPT);
356  float nearSum=T.nearEmcET+T.nearTpcPT;
357 
358  /* correct for double counting of electron track in near cone
359  rarely primTrPT<10 GeV & globPT>10 - handle this here */
360  if(T.primP.Pt()>par_trackPt) nearSum-=par_trackPt;
361  else nearSum-=T.primP.Pt();
362  T.nearTotET=nearSum;
363  T.nearTotET_noEEMC=nearSum-T.nearEtowET;
364 
365  hA[49]->Fill(nearSum);
366 
367  // printf("\n*** in nearCone T.print\n "); T.print();
368  hA[40]->Fill( T.nearEmcET);
369  hA[41]->Fill( T.cluster.ET,T.nearEmcET-T.cluster.ET);
370  float nearTotETfrac=T.cluster.ET/ T.nearTotET;
371  hA[42]->Fill(nearTotETfrac);
372  hA[192]->Fill(T.cluster.ET,nearTotETfrac);
373 
374  } // end of loop over tracks
375  }// end of loop over vertices
376 
377 
378 }
379 
380 
381 //________________________________________________
382 //________________________________________________
383 int
384 St2009WMaker::matchTrack2Cluster(){
385  // printf("******* matchCluster() nVert=%d\n",wEve.vertex.size());
386  //wEve.print();
387  int nTr=0;
388  float Rcylinder= mBtowGeom->Radius();
389  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
390  WeveVertex &V=wEve.vertex[iv];
391  float zVert=V.z;
392  for(unsigned int it=0;it<V.eleTrack.size();it++) {
393  WeveEleTrack &T=V.eleTrack[it];
394  T.isMatch2Cl=false; // just in case (was already set in constructor)
395  if(T.pointTower.id==0) continue;
396 
397  float trackPT=T.prMuTrack->momentum().perp();
398  T.cluster=maxBtow2x2( T.pointTower.iEta, T.pointTower.iPhi,zVert);
399  hA[33]->Fill( T.cluster.ET);
400  hA[34]->Fill(T.cluster.adcSum,trackPT);
401  hA[110]->Fill( T.cluster.ET);
402  if (T.cluster.ET <par_clustET) continue; // too low energy
403  hA[20]->Fill("CL",1.);
404 
405  //.. spacial separation (track - cluster)
406  TVector3 D=T.pointTower.R-T.cluster.position;
407  hA[43]->Fill( T.cluster.energy,D.Mag());
408  hA[44]->Fill( T.cluster.position.z(),D.z());
409  float delPhi=T.pointTower.R.DeltaPhi(T.cluster.position);
410  // printf("aaa %f %f %f phi=%f\n",D.x(),D.y(),D.z(),delPhi);
411  hA[45]->Fill( T.cluster.energy,Rcylinder*delPhi);// wrong?
412  hA[46]->Fill( D.Mag());
413 
414  if(D.Mag()>par_delR3D) continue;
415  hA[20]->Fill("#Delta R",1.);
416  hA[111]->Fill( T.cluster.ET);
417 
418  // ........compute surroinding cluster energy
419  int iEta=T.cluster.iEta;
420  int iPhi=T.cluster.iPhi;
421  T.cl4x4=sumBtowPatch(iEta-1,iPhi-1,4,4,zVert); // needed for lumi monitor
422  hA[37]->Fill( T.cl4x4.ET);
423  hA[38]->Fill(T.cluster.energy, T.cl4x4.energy-T.cluster.energy);
424 
425  //sum tpc PT near lepton candidate track
426  T.smallNearTpcPT=sumTpcCone(V.id,T.primP,3,T.smallNearNTr);
427  hA[56]->Fill(T.smallNearTpcPT);
428  T.smallNearTpcPT-=par_trackPt;
429 
430  float frac24=T.cluster.ET/(T.cl4x4.ET);
431  hA[39]->Fill(frac24);
432  hA[191]->Fill(T.cluster.ET,frac24);
433  if(frac24<par_clustFrac24) continue;
434  T.isMatch2Cl=true; // cluster is matched to TPC track
435 
436  hA[20]->Fill("fr24",1.);
437  hA[112]->Fill( T.cluster.ET);
438 
439  nTr++;
440  }// end of one vertex
441  }// end of vertex loop
442 
443  if(nTr<=0) return -1; // cancel event w/o good tracks
444  hA[0]->Fill("Tr2Cl",1.0);
445  return 0;
446 
447 }
448 
449 //________________________________________________
450 //________________________________________________
451 WeveCluster
452 St2009WMaker::maxBtow2x2(int iEta, int iPhi, float zVert){
453  //printf(" maxBtow2x2 seed iEta=%d iPhi=%d \n",iEta, iPhi);
454  const int L=2; // size of the summed square
455 
456  WeveCluster maxCL;
457  // just 4 cases of 2x2 clusters
458  float maxET=0;
459  int I0=iEta-1;
460  int J0=iPhi-1;
461  for(int I=I0;I<=I0+1;I++){
462  for(int J=J0;J<=J0+1;J++) {
463  WeveCluster CL=sumBtowPatch(I,J,L,L,zVert);
464  if(maxET>CL.ET) continue;
465  maxET=CL.ET;
466  maxCL=CL;
467  // printf(" newMaxETSum=%.1f iEta=%d iPhi=%d \n",maxET, I,J);
468  }
469  }// 4 combinations done
470  //printf(" final inpEve=%d SumET2x2=%.1f \n",nInpEve,maxET);
471  return maxCL;
472 }
473 
474 
475 //________________________________________________
476 //________________________________________________
477 WeveCluster
478 St2009WMaker::sumBtowPatch(int iEta, int iPhi, int Leta,int Lphi, float zVert){
479  //printf(" eveID=%d btowSquare seed iEta=%d[+%d] iPhi=%d[+%d] zVert=%.0f \n",wEve.id,iEta,Leta, iPhi,Lphi,zVert);
480  WeveCluster CL; // object is small, not to much overhead in creating it
481  CL.iEta=iEta;
482  CL.iPhi=iPhi;
483  TVector3 R;
484  double sumW=0;
485  float Rcylinder= mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
486  for(int i=iEta; i<iEta+Leta;i++){// trim in eta-direction
487  if(i<0) continue;
488  if(i>=mxBTetaBin) continue;
489  for(int j=iPhi;j<iPhi+Lphi;j++) {// wrap up in the phi-direction
490  int jj=(j+mxBTphiBin)%mxBTphiBin;// keep it always positive
491  //if(L<5) printf("n=%2d i=%d jj=%d\n",CL.nTower,i,jj);
492  int softID= mapBtowIJ2ID[ i+ jj*mxBTetaBin];
493  float ene= wEve.bemc.eneTile[kBTow][softID-1];
494  if(ene<=0) continue; // skip towers w/o energy
495  float adc= wEve.bemc.adcTile[kBTow][softID-1];
496  float delZ=positionBtow[softID-1].z()-zVert;
497  float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
498  float ET=ene*e2et;
499  float logET=log10(ET+0.5);
500  CL.nTower++;
501  CL.energy+=ene;
502  CL.ET+=ET;
503  CL.adcSum+=adc;
504  if(logET>0) {
505  R+=logET*positionBtow[softID-1];
506  sumW+=logET;
507  }
508  //if(Leta==2) printf(" iEta=%d iPhi=%d ET=%.1f ene=%.1f sum=%.1f logET=%f sumW=%f\n",i,j,ET,ene,CL.energy,logET,sumW);
509  }
510  // printf(" end btowSquare: iEta=%d nTw=%d, ET=%.1f adc=%.1f\n",i,CL.nTower,CL.ET,CL.adcSum);
511  if(sumW>0) {
512  CL.position=1./sumW*R; // weighted cluster position
513  } else {
514  CL.position=TVector3(0,0,999);
515  }
516  }
517  return CL;
518 }
519 
520 //________________________________________________
521 //________________________________________________
522 int
523 St2009WMaker::extendTrack2Barrel(){// return # of extended tracks
524  //wEve.print();
525  //printf("******* extendTracks() nVert=%d\n",wEve.vertex.size());
526  int nTrB=0;
527  for(unsigned int iv=0;iv<wEve.vertex.size();iv++) {
528  WeveVertex &V=wEve.vertex[iv];
529  for(unsigned int it=0;it<V.eleTrack.size();it++) {
530  WeveEleTrack &T=V.eleTrack[it];
531  //.... extrapolate track to the barrel @ R=entrance....
532  const StPhysicalHelixD TrkHlx=T.prMuTrack->outerHelix();
533  float Rcylinder= mBtowGeom->Radius();
534  pairD d2;
535  d2 = TrkHlx.pathLength(Rcylinder);
536  //printf(" R=%.1f path 1=%f, 2=%f, period=%f, R=%f\n",Rctb,d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
537 
538  // assert(d2.first<0); // propagate backwards
539  // assert(d2.second>0); // propagate forwards
540  if(d2.first>=0 || d2.second<=0) {
541  LOG_WARN<< Form("MatchTrk , unexpected solution for track crossing CTB\n d2.firts=%f, second=%f, swap them", d2.first, d2.second)<<endm;
542  float xx=d2.first;
543  d2.first=d2.second;
544  d2.second=xx;
545  }
546 
547  // extrapolate track to cylinder
548  StThreeVectorD posR = TrkHlx.at(d2.second);
549  //printf(" punch2 x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f\n",posCTB.x(),posCTB.y(),posCTB.z(),xmagn);
550 
551  float etaF=posR.pseudoRapidity();
552  float phiF=posR.phi();
553  int iEta,iPhi;
554  if( L2algoEtaPhi2IJ(etaF, phiF,iEta,iPhi)) continue;
555  nTrB++;
556  hA[20]->Fill("@B",1.);
557  //printf(" phi=%.0f deg, eta=%.2f, iEta=%d, iPhi=%d\n",posCTB.phi()/3.1416*180.,posCTB. pseudoRapidity(),iEta, iPhi);
558  int twID= mapBtowIJ2ID[ iEta+ iPhi*mxBTetaBin];
559  // printf("hit Tower ID=%d\n",twID);
560 
561  T.pointTower.id=twID;
562  T.pointTower.R=TVector3(posR.x(),posR.y(),posR.z());
563  T.pointTower.iEta=iEta;
564  T.pointTower.iPhi=iPhi;
565  }
566  }// end of loop over vertices
567 
568  if(nTrB<=0) return -1; // cancel event w/o good tracks
569  hA[0]->Fill("TrB",1.0);
570  return 0;
571 
572 }
573 
574 
575 //________________________________________________
576 //________________________________________________
577 //________________________________________________
578 float
579 St2009WMaker::sumBtowCone( float zVert, TVector3 refAxis, int flag, int &nTow){
580  /* flag=1 : only delta phi cut; flag=2 use 2D cut */
581  assert(flag==1 || flag==2);
582  double ptSum=0;
583 
584  //.... process BTOW hits
585  for(int i=0;i< mxBtow;i++) {
586  float ene=wEve.bemc.eneTile[kBTow][i];
587  if(ene<=0) continue;
588  TVector3 primP=positionBtow[i]-TVector3(0,0,zVert);
589  primP.SetMag(ene); // it is 3D momentum in the event ref frame
590  if(flag==1) {
591  float deltaPhi=refAxis.DeltaPhi(primP);
592  if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
593  }
594  if(flag==2) {
595  float deltaR=refAxis.DeltaR(primP);
596  if(deltaR> par_nearDeltaR) continue;
597  }
598  if(primP.Perp()>par_countTowEt) nTow++;
599  ptSum+=primP.Perp();
600  }
601 
602  return ptSum;
603 }
604 
605 
606 //________________________________________________
607 //________________________________________________
608 float
609 St2009WMaker::sumEtowCone(float zVert, TVector3 refAxis, int flag,int &nTow){
610  /* flag=1 : only delta phi cut; flag=2 use 2D cut */
611  assert(flag==1 || flag==2);
612 
613  float ptsum=0;
614 
615  //....loop over all phi bins
616  for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
617  for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
618  float ene=wEve.etow.ene[iphi][ieta];
619  if(ene<=0) continue; //skip towers with no energy
620  TVector3 primP=positionEtow[iphi][ieta]-TVector3(0,0,zVert);
621  primP.SetMag(ene); // it is 3D momentum in the event ref frame
622  if(flag==1) {
623  float deltaPhi=refAxis.DeltaPhi(primP);
624  if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
625  }
626  if(flag==2) {
627  float deltaR=refAxis.DeltaR(primP);
628  if(deltaR> par_nearDeltaR) continue;
629  }
630  if(primP.Perp()>par_countTowEt) nTow++;
631  ptsum+=primP.Perp();
632  }
633  }
634 
635  return ptsum;
636 }
637 
638 // $Log: St2009W_algo.cxx,v $
639 // Revision 1.25 2011/09/14 14:23:20 stevens4
640 // update used for cross section PRD paper
641 //
642 // Revision 1.24 2010/04/27 16:53:45 stevens4
643 // add code to remove events tagged as Zs from W candidates
644 //
645 // Revision 1.23 2010/04/16 14:35:32 balewski
646 // fix borken header
647 //
648 // Revision 1.22 2010/03/23 01:31:40 seelej
649 // Fix to the filling of the histograms for the background systematic.
650 //
651 // Revision 1.21 2010/03/22 01:45:58 seelej
652 // *** empty log message ***
653 //
654 // Revision 1.20 2010/03/22 01:33:13 seelej
655 // Additional change.
656 //
657 // Revision 1.19 2010/03/22 00:55:45 seelej
658 // Change to use the signed pT balance in background subtraction instead of the unsigned pT balance.
659 //
660 // Revision 1.18 2010/03/18 18:46:40 balewski
661 // simplified sPtBalance calculation
662 //
663 // Revision 1.17 2010/03/18 16:52:17 balewski
664 // corrected sPtBalance for no-endcap
665 //
666 // Revision 1.16 2010/03/18 15:34:43 seelej
667 // Changed the definition of sPtBalance - Joe
668 //
669 // Revision 1.15 2010/03/14 22:50:31 balewski
670 // *** empty log message ***
671 //
672 // Revision 1.14 2010/03/12 21:08:11 balewski
673 // simplify logic for filling histos for noE background
674 //
675 // Revision 1.13 2010/02/26 21:40:00 seelej
676 // Joe : Fix to code. Forgot to remove an older piece of code when doing a previous update.
677 //
678 // Revision 1.12 2010/02/22 15:49:34 seelej
679 // Joe : Changes to code for inclusion of background subtraction and systematic studies
680 //
681 // Revision 1.11 2010/01/28 03:42:55 balewski
682 // cleanup
683 //
684 // Revision 1.10 2010/01/27 22:12:24 balewski
685 // spin code matched to x-section code
686 //
687 // Revision 1.9 2010/01/26 18:48:11 stevens4
688 // include |eta|<1 cut for all W yields
689 //
690 // Revision 1.8 2010/01/23 02:35:38 stevens4
691 // add ability to scale jet et and use real btow peds for rcf mc
692 //
693 // Revision 1.7 2010/01/10 03:01:37 balewski
694 // cleanup & nicer histos
695 //
696 // Revision 1.6 2010/01/10 01:45:10 stevens4
697 // fix plots w/o EEMC in veto
698 //
699 // Revision 1.5 2010/01/09 00:07:16 stevens4
700 // add jet finder
701 //
702 // Revision 1.4 2010/01/06 19:16:47 stevens4
703 // track cuts now on primary component, cleanup
704 //
705 // Revision 1.3 2010/01/06 14:11:13 balewski
706 // one Z-plot added
707 //
708 // Revision 1.2 2010/01/06 04:22:15 balewski
709 // added Q/PT plot for Zs, more cleanup
710 //
711 // Revision 1.1 2009/11/23 23:00:18 balewski
712 // code moved spin-pool
713 //
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StMuTrack::dcaZ
Float_t dcaZ(Int_t vtx_id=-1) const
Z component of global DCA.
Definition: StMuTrack.cxx:351
StMuTrack::chi2
Double_t chi2() const
Returns chi2 of fit.
Definition: StMuTrack.h:251
StJet::jetPt
float jetPt
Pt (stored for convenience when drawing TTree)
Definition: StJet.h:125
StJet::zVertex
float zVertex
position of vertex used to reconstruct jet
Definition: StJet.h:134
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StMuTrack::nHitsFit
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
StJet::nBtowers
int nBtowers
The number of Barrel towers in this jet.
Definition: StJet.h:107
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
StMuTrack::nHitsPoss
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
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::dca
StThreeVectorF dca(Int_t vtx_id=-1) const
Returns 3D distance of closest approach to primary vertex.
Definition: StMuTrack.cxx:359
StMuTrack::lastPoint
const StThreeVectorF & lastPoint() const
Returns positions of last measured point.
Definition: StMuTrack.h:262
StMuTrack::dEdx
Double_t dEdx() const
Returns measured dE/dx value.
Definition: StMuTrack.h:248
StMuTrack::dcaD
Float_t dcaD(Int_t vtx_id=-1) const
Signed radial component of global DCA (projected)
Definition: StMuTrack.cxx:332
StJet::nEtowers
int nEtowers
The number of Endcap towers in this jet.
Definition: StJet.h:110
StMuTrack::globalTrack
const StMuTrack * globalTrack() const
Returns pointer to associated global track. Null pointer if no global track available.
Definition: StMuTrack.h:272
StMuTrack::outerHelix
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
StJet
Definition: StJet.h:91
Generated by   doxygen 1.8.5