StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2011W_algo.cxx
1 // $Id: St2011W_algo.cxx,v 1.22 2016/01/08 02:08:49 jlzhang Exp $
2 //
3 //*-- Author : Jan Balewski, MIT
4 //*-- Author for Endcap: Justin Stevens, IUCF
5 
6 #include "TF1.h"
7 
8 #include "StEmcUtil/geometry/StEmcGeom.h"
9 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
10 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
11 #include "WeventDisplay.h"
12 
13 //new jet tree format
14 #include "StSpinPool/StJetEvent/StJetEvent.h"
15 #include "StSpinPool/StJetEvent/StJetVertex.h"
16 #include "StSpinPool/StJetEvent/StJetCandidate.h"
17 #include "StSpinPool/StJetEvent/StJetTower.h"
18 
19 #include "St2011WMaker.h"
20 
21 //________________________________________________
22 //________________________________________________
23 void
24 St2011WMaker::find_W_boson(){
25 
26  if(!wEve->l2bitET) return;
27 
28  //printf("========= find_W_boson() \n");
29  int nNoNear=0,nNoAway=0,nEta1=0,nGoldW=0,nGoldWp=0,nGoldWn=0;
30  //remove events tagged as Zs
31  if(wEve->zTag) return;
32 
33  // search for Ws ............
34  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
35  WeveVertex &V=wEve->vertex[iv];
36  for(unsigned int it=0;it<V.eleTrack.size();it++) {
37  WeveEleTrack &T=V.eleTrack[it];
38 
39  getJetEvent();
40  hA[114]->Fill(T.cluster.position.PseudoRapidity(),T.sPtBalance2);
41  hA[115]->Fill(T.cluster.position.Phi(),T.sPtBalance2);
42 
43  if(T.pointTower.id<=0) continue; //skip endcap towers
44  if(T.isMatch2Cl==false) continue;
45  assert(T.cluster.nTower>0); // internal logical error
46  assert(T.nearTotET>0); // internal logical error
47 
48  //make cut on lepton eta
49  if(T.primP.Eta() < par_leptonEtaLow || T.primP.Eta() > par_leptonEtaHigh) continue;
50  hA[20]->Fill("eta1",1.);
51  nEta1++;
52 
53  //define charge separation
54  float q2pt_g = T.glMuTrack->charge()/T.glMuTrack->pt();
55  float q2pt_p = T.prMuTrack->charge()/T.prMuTrack->pt();
56  float hypCorr_g = q2pt_g*(T.cluster.ET);
57  float hypCorr_p = q2pt_p*(T.cluster.ET);
58 
59  //remove ambiguous charges from BG treatment histos
60  if( fabs(hypCorr_p) > par_QET2PTlow && fabs(hypCorr_p) < par_QET2PThigh) {
61  //signal plots w/o EEMC in awayside veto
62  if(T.cluster.ET/T.nearTotET_noEEMC>par_nearTotEtFrac){
63  if(T.sPtBalance_noEEMC2>par_ptBalance ) {//only signed ptBalance cut
64  hA[140]->Fill(T.cluster.ET);
65  hA[240]->Fill(T.prMuTrack->eta(),T.cluster.ET);
66  if (T.prMuTrack->charge() < 0) {
67  hA[184+3]->Fill(T.cluster.ET);
68  } else if (T.prMuTrack->charge() > 0) {
69  hA[184+4]->Fill(T.cluster.ET);
70  }
71  }
72  }
73  }
74 
75  //fill plot for background
76  if(T.cluster.ET > par_highET) {
77  if(T.prMuTrack->charge()>0) hA[251]->Fill(T.cluster.ET/T.nearTotET,T.sPtBalance2);
78  else if(T.prMuTrack->charge()<0) hA[252]->Fill(T.cluster.ET/T.nearTotET,T.sPtBalance2);
79  hA[135]->Fill(T.awayTotET,T.sPtBalance2);
80  }
81 
82  // track matched to cluster plots
83  StThreeVectorF ri=T.glMuTrack->firstPoint();
84  StThreeVectorF ro=T.glMuTrack->lastPoint();
85  int sec = WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity());
86  if((sec < 5 || sec > 7) && sec!=21) { //skip sectors with dead padrows for this
87  hA[63]->Fill(T.prMuTrack->nHitsFit());
88  hA[64]->Fill(1.*T.prMuTrack->nHitsFit()/T.prMuTrack->nHitsPoss());
89  hA[65]->Fill(ri.perp());
90  }
91 
92  if(T.cluster.ET /T.nearTotET< par_nearTotEtFrac) continue; // too large nearET
93 
94  hA[20]->Fill("noNear",1.);
95  nNoNear++;
96  hA[112]->Fill( T.cluster.ET); // for Joe
97  hA[50]->Fill(T.awayTpcPT);
98  hA[51]->Fill(T.awayBtowET);
99  hA[54]->Fill(T.awayTotET);
100  hA[52]->Fill(T.cluster.ET,T.awayTotET);
101  hA[53]->Fill(T.cluster.ET,T.awayEmcET);
102  hA[55]->Fill(T.awayEtowET);
103  hA[60]->Fill(T.cluster.ET,T.awayTpcPT);
104 
105  hA[132]->Fill(T.cluster.ET,T.ptBalance.Perp());
106  hA[133]->Fill(T.awayTotET,T.ptBalance.Perp());
107  hA[134]->Fill(T.cluster.ET,T.sPtBalance);
108  hA[135]->Fill(T.awayTotET,T.sPtBalance);
109  hA[209]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
110  if(T.cluster.ET > par_highET) hA[253]->Fill(T.awayTotET,T.sPtBalance);
111 
112 
113 
114  //if(T.cluster.ET> par_highET && T.awayTotET>par_awayET) continue;
115 
116  // add 2 histograms hA[138],hA[139], the charged separation sPtBalance vs Et 2D plots for background estimation;
117  // Jinlong 12/09/2014
118  if(T.prMuTrack->charge()>0) {
119  hA[130]->Fill(T.cluster.ET,T.sPtBalance);
120  hA[138]->Fill(T.cluster.ET,T.sPtBalance2);
121  } else if(T.prMuTrack->charge()<0) {
122  hA[131]->Fill(T.cluster.ET,T.sPtBalance);
123  hA[139]->Fill(T.cluster.ET,T.sPtBalance2);
124  }
125 
126  //remove ambiguous charges from BG treatment histos
127  if( fabs(hypCorr_p) > par_QET2PTlow && fabs(hypCorr_p) < par_QET2PThigh) {
128  for (int i=0; i<=20; i++) { // i index not currently used
129  for (int j=0; j<=80; j++) {
130  float pTBal_cut = 5.+0.25*((float) j);
131  if (T.sPtBalance2<pTBal_cut) {
132  if (T.prMuTrack->charge() < 0) {
133  hA[142+i]->Fill(T.cluster.ET,j);
134  } else if (T.prMuTrack->charge() > 0) {
135  hA[163+i]->Fill(T.cluster.ET,j);
136  }
137  }
138  }
139  }
140 
141  //plots for backg sub yield
142  if(T.sPtBalance2>par_ptBalance ) {
143  hA[136]->Fill(T.cluster.ET);//signal
144  hA[241]->Fill(T.prMuTrack->eta(),T.cluster.ET);
145  hA[62]->Fill(T.pointTower.iEta ,T.cluster.energy);
146  if (T.prMuTrack->charge() < 0) {
147  hA[184+1]->Fill(T.cluster.ET);
148  } else if (T.prMuTrack->charge() > 0) {
149  hA[184+2]->Fill(T.cluster.ET);
150  }
151  } else {
152  hA[137]->Fill(T.cluster.ET);//background
153  if (T.prMuTrack->charge() < 0) {
154  hA[184+5]->Fill(T.cluster.ET);
155  } else if (T.prMuTrack->charge() > 0) {
156  hA[184+6]->Fill(T.cluster.ET);
157  }
158  hA[202]->Fill(T.cluster.position.PseudoRapidity(),T.prMuTrack->pt());
159  hA[204]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.energy/T.prMuTrack->p().mag());
160  }
161  }
162 
163 
164  if(T.sPtBalance2>par_ptBalance && T.cluster.ET>par_highET){/***************************/
165  printf("\n WWWWWWWWWWWWWWWWWWWWW Barrel \n");
166  wDisaply->exportEvent( "WB", V, T, iv);
167  wEve->print();
168  }/***************************/
169 
170 
171  //put final W cut here
172  if(T.sPtBalance2<par_ptBalance) continue;
173  hA[20]->Fill("noAway",1.0);
174  nNoAway++;
175 
176  //****loop over branch with EEMC**** Feb 23:27:35, Jinlong Zhang
177  StJetVertex* jetVertex = mJetEvent->vertex(V.id);
178  assert(jetVertex->position().z() == V.z); //check that vert-z position match
179  int nJetsWE = jetVertex->numberOfJets();
180  for (int i_jet=0; i_jet< nJetsWE; i_jet++){//loop over jets
181  StJetCandidate* jet = jetVertex->jet(i_jet);
182  TVector3 jetVec; //vector for jet momentum
183  jetVec.SetPtEtaPhi(jet->pt(),jet->eta(),jet->phi());
184  if(jetVec.DeltaR(T.primP) > par_nearDeltaR) {
185  hA[126]->Fill(T.primP.Eta(), jet->eta());
186  hA[127]->Fill(T.primP.Phi(), jet->phi());
187  hA[128]->Fill(jet->pt());
188  hA[129]->Fill(jet->pt(),jetVec.DeltaPhi(T.primP));
189  hA[130]->Fill(jet->pt(),jet->eta());
190  }
191  }
192  hA[125]->Fill(T.jetCount);
193 
194  //::::::::::::::::::::::::::::::::::::::::::::::::
195  //:::::accepted W events for x-section :::::::::::
196  //::::::::::::::::::::::::::::::::::::::::::::::::
197 
198  hA[113]->Fill( T.cluster.ET);//for Joe
199 
200  hA[90]->Fill( T.cluster.ET);
201  hA[92]->Fill( T.cluster.ET,T.glMuTrack->dEdx()*1e6);
202  //hA[93]->Fill( T.cluster.ET,T.glMuTrack->dca(V.id).mag());
203  int k=0; if(T.prMuTrack->charge()<0) k=1;
204  hA[94+k]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
205  // h95 used above
206 
207  //plots to investigate east/west yield diff
208  hA[200]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
209  hA[201]->Fill(T.cluster.position.PseudoRapidity(),T.prMuTrack->pt());
210  hA[203]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.energy/T.prMuTrack->p().mag());
211  hA[205]->Fill(T.prMuTrack->lastPoint().pseudoRapidity(),T.prMuTrack->lastPoint().phi());
212 
213  //Q/pT plot
214  hA[100]->Fill(T.cluster.ET,T.glMuTrack->charge()/T.glMuTrack->pt());
215  hA[101]->Fill(T.cluster.ET,T.prMuTrack->charge()/T.prMuTrack->pt());
216  hA[102]->Fill(T.cluster.ET,hypCorr_g);
217  hA[103]->Fill(T.cluster.ET,hypCorr_p);
218 
219  //for each sector
220  int isec = WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity())-1;
221  hA[260+isec]->Fill(T.cluster.ET,T.glMuTrack->charge()/T.glMuTrack->pt());
222  hA[284+isec]->Fill(T.cluster.ET,T.prMuTrack->charge()/T.prMuTrack->pt());
223  hA[356+isec]->Fill(T.cluster.ET,hypCorr_p);
224  if(k==0) hA[308+isec]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
225  else hA[332+isec]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
226 
227  if(T.cluster.ET<par_highET) continue; // very likely Ws
228  hA[91]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.position.Phi());
229  hA[96]->Fill(V.id);
230  hA[97]->Fill(V.funnyRank);
231  hA[98]->Fill(V.z);
232  hA[99]->Fill(T.prMuTrack->eta());
233  hA[191+k]->Fill(T.prMuTrack->eta(),T.cluster.ET);
234  hA[106]->Fill(wEve->zdcRate);
235 
236  hA[20]->Fill("goldW",1.);
237  nGoldW++;
238  if(T.prMuTrack->charge()>0) nGoldWp++;
239  else if(T.prMuTrack->charge()<0) nGoldWn++;
240  hA[104]->Fill(wEve->time);
241 
242  // free quark search
243  hA[105]->Fill(hypCorr_p,T.glMuTrack->dEdx()*1e6);
244 
245  }// loop over tracks
246  }// loop over vertices
247  if(nNoNear>0) hA[0]->Fill("noNear",1.);
248  if(nNoAway>0) hA[0]->Fill("noAway",1.);
249  if(nEta1>0) hA[0]->Fill("eta1",1.);
250  if(nGoldW>0) hA[0]->Fill("goldW",1.);
251  if(nGoldWp>0) hA[0]->Fill("goldW+",1.);
252  if(nGoldWn>0) hA[0]->Fill("goldW-",1.);
253 
254 }
255 
256 
257 //________________________________________________
258 //________________________________________________
259 void
260 St2011WMaker::tag_Z_boson(){
261 
262  float par_jetPt=10.;
263  float lowMass=70.; float highMass=140.;
264 
265  //form invariant mass from lepton candidate and jet
266  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {//vertex loop
267  WeveVertex &V=wEve->vertex[iv];
268  for(unsigned int it=0;it<V.eleTrack.size();it++) {// select track
269  WeveEleTrack &T1=V.eleTrack[it];
270  if(T1.isMatch2Cl==false) continue;
271  assert(T1.cluster.nTower>0); // internal logical error
272  assert(T1.nearTotET>0); // internal logical error
273 
274  //match lepton candidate with jet
275  TLorentzVector jetVec;
276 
277  getJetEvent(); //check that jet and W event match
278 
279  StJetVertex* jetVertex = mJetEvent->vertex(V.id);
280  int nJets = jetVertex->numberOfJets();
281  for (int i_jet=0; i_jet< nJets; i_jet++){//loop over jets
282  StJetCandidate* jet = jetVertex->jet(i_jet);
283 
284  jetVec = jet->fourMomentum();
285  if(jetVec.Pt()<par_jetPt) continue;//remove low pt jets
286 
287  //electron like cut on jets
288  float maxCluster=0.;
289  int totTowers = jet->numberOfTowers();
290  for(int itow=0;itow<totTowers;itow++){//loop over towers
291  StJetTower *tower = jet->tower(itow);
292  if(tower->detectorId()==13)//drop endcap towers
293  continue;
294 
295  int softId=tower->id();
296  //find highest 2x2 tower cluster in jet
297  TVector3 pos=positionBtow[softId-1]; int iEta,iPhi;
298  if( L2algoEtaPhi2IJ(pos.Eta(),pos.Phi(),iEta,iPhi))
299  continue;
300  float cluster=maxBtow2x2(iEta,iPhi,V.z).ET;
301  if(cluster>maxCluster) maxCluster=cluster;
302  }
303 
304  TVector3 jetVec3(jetVec.X(),jetVec.Y(),jetVec.Z());
305  if(jetVec3.DeltaR(T1.primP)<par_nearDeltaR)
306  continue;//skip jets in candidate phi isolation'cone'
307 
308  //form invM
309  float e1=T1.cluster.energy;
310  TVector3 p1=T1.primP; p1.SetMag(e1);
311  TLorentzVector ele1(p1,e1); //lepton candidate 4- mom
312  TLorentzVector sum=ele1+jetVec;
313  float invM=sqrt(sum*sum);
314  if(maxCluster/jetVec3.Pt() < 0.5) continue;
315  if(invM > lowMass && invM < highMass){
316  wEve->zTag=true;
317  //cout<<"tagged as Z mass = "<<invM<<endl;
318  }
319  }
320  }
321  }
322 
323 }
324 
325 //________________________________________________
326 //________________________________________________
327 void
328 St2011WMaker::findPtBalance(){
329 
330  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
331  WeveVertex &V=wEve->vertex[iv];
332  for(unsigned int it=0;it<V.eleTrack.size();it++) {
333  WeveEleTrack &T=V.eleTrack[it];
334 
335  getJetEvent(); //check that jet and W event match
336 
337  //****loop over branch with EEMC****
338  StJetVertex* jetVertex = mJetEvent->vertex(V.id);
339  assert(jetVertex);
340  assert(jetVertex->position().z() == V.z); //check that vert-z position match
341  int nJetsWE = jetVertex->numberOfJets();
342  int nJetsOutNearCone=0; // add to count jet in signed pT; jinlong 2014/12/19
343  for (int i_jet=0; i_jet< nJetsWE; i_jet++){//loop over jets
344  StJetCandidate* jet = jetVertex->jet(i_jet);
345  TVector3 jetVec; //vector for jet momentum
346  jetVec.SetPtEtaPhi(jet->pt(),jet->eta(),jet->phi());
347  if(jetVec.DeltaR(T.primP) > par_nearDeltaR) {
348  T.ptBalance+=jetVec;
349  nJetsOutNearCone++; // add to count jet in signed pT; jinlong 2014/12/19
350  hA[120]->Fill(T.primP.Eta(), jet->eta());
351  hA[121]->Fill(T.primP.Phi(), jet->phi());
352  hA[122]->Fill(jet->pt());
353  hA[123]->Fill(jet->pt(),jetVec.DeltaPhi(T.primP));
354  hA[124]->Fill(jet->pt(),jet->eta());
355  }
356  }
357 
358  hA[119]->Fill(nJetsOutNearCone);
359  T.jetCount = nJetsOutNearCone; // add to count jet in signed pT; jinlong 2014/12/19
360  TVector3 clustPt(T.primP.X(),T.primP.Y(),0);
361  clustPt.SetMag(T.cluster.ET);
362  T.ptBalance+=clustPt;
363  T.sPtBalance = T.ptBalance.Perp();
364  if(T.ptBalance.Dot(clustPt)<0) T.sPtBalance *=-1.;
365  T.sPtBalance2 = T.ptBalance.Dot(clustPt)/T.cluster.ET; //invariant
366 
367  //****loop over branch without EEMC****
368  StJetVertex* jetVertex_noEEMC = mJetEvent_noEEMC->vertex(V.id);
369  assert(jetVertex_noEEMC->position().z() == V.z); //check that vert-z position match
370  int nJetsNE = jetVertex_noEEMC->numberOfJets();
371  for (int i_jet=0; i_jet< nJetsNE; i_jet++){//loop over jets
372  StJetCandidate* jet = jetVertex_noEEMC->jet(i_jet);
373  TVector3 jetVec; //vector for jet momentum
374  jetVec.SetPtEtaPhi(jet->pt(),jet->eta(),jet->phi());
375  if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
376  T.ptBalance_noEEMC+=jetVec;
377  }
378  T.ptBalance_noEEMC+=clustPt;
379  T.sPtBalance_noEEMC = T.ptBalance_noEEMC.Perp();
380  if(T.ptBalance_noEEMC.Dot(clustPt)<0) T.sPtBalance_noEEMC *=-1.;
381  T.sPtBalance_noEEMC2 = T.ptBalance_noEEMC.Dot(clustPt)/T.cluster.ET; //invariant
382 
383  }// end of loop over tracks
384  }// end of loop over vertices
385 
386 }
387 
388 
389 //________________________________________________
390 //________________________________________________
391 void
392 St2011WMaker::findAwayJet(){
393  // printf("\n******* find AwayJet() nVert=%d\n",wEve->vertex.size());
394  //wEve->print();
395  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
396  WeveVertex &V=wEve->vertex[iv];
397  for(unsigned int it=0;it<V.eleTrack.size();it++) {
398  WeveEleTrack &T=V.eleTrack[it];
399 
400  // .... sum opposite in phi EMC components
401  T.awayBtowET=sumBtowCone(V.z,-T.primP,1); // '1'= only cut on delta phi
402  T.awayEmcET=T.awayBtowET;
403  T.awayEtowET=sumEtowCone(V.z,-T.primP,1);
404  T.awayEmcET+=T.awayEtowET;
405 
406  //..... add TPC ET
407  if(mMuDstMaker) T.awayTpcPT=sumTpcCone(V.id,-T.primP,1,T.pointTower.id);
408  else T.awayTpcPT=sumTpcConeFromTree(iv,-T.primP,1,T.pointTower.id);
409  T.awayTotET=T.awayEmcET+T.awayTpcPT;
410  T.awayTotET_noEEMC=T.awayBtowET+T.awayTpcPT;
411  //printf("\n*** in awayTpc=%.1f awayEmc=%.1f\n ",T.awayTpcPT,T.awayEmcET); T.print();
412 
413  }// end of loop over tracks
414  }// end of loop over vertices
415 }
416 
417 //________________________________________________
418 //________________________________________________
419 void
420 St2011WMaker::findNearJet(){
421  //printf("\n******* findNearJet() nVert=%d\n",wEve->vertex.size());
422 
423  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
424  WeveVertex &V=wEve->vertex[iv];
425  for(unsigned int it=0;it<V.eleTrack.size();it++) {
426  WeveEleTrack &T=V.eleTrack[it];
427 
428  // .... sum EMC-jet component
429  T.nearBtowET=sumBtowCone(V.z,T.primP,2); // '2'=2D cone
430  T.nearEmcET+=T.nearBtowET;
431  T.nearEtowET=sumEtowCone(V.z,T.primP,2);
432  T.nearEmcET+=T.nearEtowET;
433  // .... sum TPC-near component
434  if(mMuDstMaker) T.nearTpcPT=sumTpcCone(V.id,T.primP,2,T.pointTower.id); // '2'=2D cone
435  else T.nearTpcPT=sumTpcConeFromTree(iv,T.primP,2,T.pointTower.id);
436  float nearSum=T.nearEmcET+T.nearTpcPT;
437 
438  //fill histos separately for 2 types of events
439  if(T.pointTower.id>0) { //only barrel towers
440  /* correct for double counting of electron track in near cone rarely primTrPT<10 GeV & globPT>10 - handle this here */
441  if(T.primP.Pt()>par_trackPt) nearSum-=par_trackPt;
442  else nearSum-=T.primP.Pt();
443  T.nearTotET=nearSum;
444  T.nearTotET_noEEMC=nearSum-T.nearEtowET;
445  float nearTotETfrac=T.cluster.ET/ T.nearTotET;
446  float nearTotETfrac_noEEMC=T.cluster.ET/ T.nearTotET_noEEMC;
447 
448  //move requirement here for consistency, but now calc nearCone for all candidates
449  if(T.isMatch2Cl==false) continue;
450 
451  hA[40]->Fill(T.nearEmcET);
452  hA[41]->Fill(T.cluster.ET,T.nearEmcET-T.cluster.ET);
453  hA[42]->Fill(nearTotETfrac);
454  hA[131]->Fill(nearTotETfrac_noEEMC);
455  hA[47]->Fill(T.nearTpcPT);
456  hA[48]->Fill(T.nearEmcET,T.nearTpcPT);
457  hA[49]->Fill(nearSum);
458  hA[250]->Fill(T.cluster.ET,nearTotETfrac);
459 
460  // check east/west yield diff
461  hA[210]->Fill(T.cluster.position.PseudoRapidity(),T.nearEtowET);
462  if(T.cluster.position.PseudoRapidity()>0) hA[211]->Fill(T.cluster.position.Phi(),T.nearEtowET);
463  else hA[212]->Fill(T.cluster.position.Phi(),T.nearEtowET);
464 
465  }
466  else if(T.pointTower.id<0) { //only endcap towers
467  /* correct for double counting of electron track in near cone rarely primTrPT<10 GeV & globPT>10 - handle this here */
468  if(T.prMuTrack->flag()==301){ //short tracks aren't added to nearCone
469  if(T.primP.Pt()>parE_trackPt) nearSum-=parE_trackPt;
470  else nearSum-=T.primP.Pt();
471  }
472  T.nearTotET=nearSum;
473  T.nearTotET_noEEMC=nearSum-T.nearEtowET;
474  float nearTotETfrac=T.cluster.ET/ T.nearTotET;
475 
476  //move requirement here for consistency, but now calc nearCone for all candidates
477  if(T.isMatch2Cl==false) continue;
478 
479  hE[40]->Fill(T.nearEmcET);
480  hE[41]->Fill(T.cluster.ET,T.nearEmcET-T.cluster.ET);
481  hE[42]->Fill(nearTotETfrac);
482  hE[70]->Fill(T.cluster.ET/T.nearEmcET);
483  hE[71]->Fill(T.cluster.ET/T.nearEtowET);
484  hE[47]->Fill(T.nearTpcPT);
485  hE[48]->Fill(T.nearEmcET,T.nearTpcPT);
486  hE[49]->Fill(nearSum);
487  }
488 
489  } // end of loop over tracks
490  }// end of loop over vertices
491 
492 
493 }
494 
495 //________________________________________________
496 //________________________________________________
497 float
498 St2011WMaker::sumBtowCone( float zVert, TVector3 refAxis, int flag){
499  /* flag=1 : only delta phi cut; flag=2 use 2D cut */
500  assert(flag==1 || flag==2);
501  double ptSum=0;
502 
503  //.... process BTOW hits
504  for(int i=0;i< mxBtow;i++) {
505  float ene=wEve->bemc.eneTile[kBTow][i];
506  if(ene<=0) continue;
507  TVector3 primP=positionBtow[i]-TVector3(0,0,zVert);
508  primP.SetMag(ene); // it is 3D momentum in the event ref frame
509  if(flag==1) {
510  float deltaPhi=refAxis.DeltaPhi(primP);
511  if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
512  }
513  if(flag==2) {
514  float deltaR=refAxis.DeltaR(primP);
515  if(deltaR> par_nearDeltaR) continue;
516  }
517  ptSum+=primP.Perp();
518  }
519 
520  return ptSum;
521 }
522 
523 //________________________________________________
524 //________________________________________________
525 float
526 St2011WMaker::sumTpcConeFromTree(int vertID, TVector3 refAxis, int flag,int pointTowId){
527 
528  // flag=2 use 2D cut, 1= only delta phi
529 
530  assert(vertID>=0);
531  assert(vertID<(int)wEve->vertex.size());
532 
533  double ptSum=0;
534  WeveVertex &V=wEve->vertex[vertID];
535  for(unsigned int it=0;it<V.prTrList.size();it++){
536  StMuTrack *prTr=V.prTrList[it];
537  if(prTr->flag()<=0) continue;
538  if(prTr->flag()!=301 && pointTowId>0) continue;// TPC-only regular tracks for barrel candidate
539  if(prTr->flag()!=301 && pointTowId<0) continue;// TPC-only regular tracks for endcap candidate
540  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
541  if(hitFrac<par_nHitFrac) continue;
542  StThreeVectorF prPvect=prTr->p();
543  TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
544  // printf(" prTrID=%4d prTrEta=%.3f prTrPhi/deg=%.1f prPT=%.1f nFitPts=%d\n", prTr->id(),prTr->eta(),prTr->phi()/3.1416*180.,prTr->pt(),prTr->nHitsFit());
545  if(flag==1) {
546  float deltaPhi=refAxis.DeltaPhi(primP);
547  if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
548  }
549  if(flag==2) {
550  float deltaR=refAxis.DeltaR(primP);
551  if(deltaR>par_nearDeltaR) continue;
552 
553  }
554  float pT=prTr->pt();
555  // printf(" passed pt=%.1f\n",pT);
556 
557  //separate quench for barrel and endcap candidates
558  if(pT>par_trackPt && pointTowId>0) ptSum+=par_trackPt;
559  else if(pT>parE_trackPt && pointTowId<0) ptSum+=parE_trackPt;
560  else ptSum+=pT;
561  }
562  return ptSum;
563 }
564 
565 
566 // ************* Barrel Code ************ //
567 // ************************************** //
568 
569 //________________________________________________
570 //________________________________________________
571 int
572 St2011WMaker::extendTrack2Barrel(){// return # of extended tracks
573  //printf("******* extendTracks() nVert=%d\n",wEve->vertex.size());
574  if(!wEve->l2bitET) return 0; //fire barrel trigger
575 
576  int nTrB=0;
577  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
578  WeveVertex &V=wEve->vertex[iv];
579  if(V.rank<0) continue; //remove vertex for endcap algo only
580  for(unsigned int it=0;it<V.eleTrack.size();it++) {
581  WeveEleTrack &T=V.eleTrack[it];
582  if(T.prMuTrack->flag()!=301) continue; //remove track for endcap algo only
583 
584  //do eta sorting at track level (tree analysis)
585  if(T.primP.Eta() < par_leptonEtaLow || T.primP.Eta() > par_leptonEtaHigh) continue;
586 
587  //.... extrapolate track to the barrel @ R=entrance....
588  const StPhysicalHelixD TrkHlx=T.prMuTrack->outerHelix();
589  float Rcylinder= mBtowGeom->Radius();
590  pairD d2;
591  d2 = TrkHlx.pathLength(Rcylinder);
592  //printf(" R=%.1f path 1=%f, 2=%f, period=%f, R=%f\n",Rctb,d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
593 
594  // assert(d2.first<0); // propagate backwards
595  // assert(d2.second>0); // propagate forwards
596  if(d2.first>=0 || d2.second<=0) {
597  LOG_WARN<< Form("MatchTrk , unexpected solution for track crossing CTB\n d2.firts=%f, second=%f, swap them", d2.first, d2.second)<<endm;
598  float xx=d2.first;
599  d2.first=d2.second;
600  d2.second=xx;
601  }
602 
603  // extrapolate track to cylinder
604  StThreeVectorD posR = TrkHlx.at(d2.second);
605  //printf(" punch2 x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f\n",posCTB.x(),posCTB.y(),posCTB.z(),xmagn);
606  float etaF=posR.pseudoRapidity();
607  float phiF=posR.phi();
608  int iEta,iPhi;
609  if( L2algoEtaPhi2IJ(etaF, phiF,iEta,iPhi)) continue;
610  nTrB++;
611  hA[20]->Fill("@B",1.);
612  //printf(" phi=%.0f deg, eta=%.2f, iEta=%d, iPhi=%d\n",posCTB.phi()/3.1416*180.,posCTB. pseudoRapidity(),iEta, iPhi);
613  int twID= mapBtowIJ2ID[ iEta+ iPhi*mxBTetaBin];
614  // printf("hit Tower ID=%d\n",twID);
615 
616  T.pointTower.id=twID;
617  T.pointTower.R=TVector3(posR.x(),posR.y(),posR.z());
618  T.pointTower.iEta=iEta;
619  T.pointTower.iPhi=iPhi;
620  //T.print();
621  }
622  }// end of loop over vertices
623 
624  if(nTrB<=0) return -1;
625  hA[0]->Fill("TrB",1.0);
626  return 0;
627 
628 }
629 
630 //________________________________________________
631 //________________________________________________
632 int
633 St2011WMaker::matchTrack2BtowCluster(){
634  // printf("******* matchCluster() nVert=%d\n",wEve->vertex.size());
635  int nTr=0;
636  float Rcylinder= mBtowGeom->Radius();
637  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
638  WeveVertex &V=wEve->vertex[iv];
639  float zVert=V.z;
640  for(unsigned int it=0;it<V.eleTrack.size();it++) {
641  WeveEleTrack &T=V.eleTrack[it];
642  if(T.pointTower.id<=0) continue; //skip endcap towers
643 
644  float trackPT=T.prMuTrack->momentum().perp();
645  T.cluster=maxBtow2x2( T.pointTower.iEta, T.pointTower.iPhi,zVert);
646  hA[33]->Fill( T.cluster.ET);
647  hA[34]->Fill(T.cluster.adcSum,trackPT);
648  hA[110]->Fill( T.cluster.ET);
649 
650  // ........compute surroinding cluster energy
651  int iEta=T.cluster.iEta;
652  int iPhi=T.cluster.iPhi;
653  T.cl4x4=sumBtowPatch(iEta-1,iPhi-1,4,4,zVert); // needed for lumi monitor
654 
655  if (T.cluster.ET <par_clustET) continue; // too low energy
656  hA[20]->Fill("CL",1.);
657 
658  hA[206]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
659 
660  hA[37]->Fill( T.cl4x4.ET);
661  hA[38]->Fill(T.cluster.energy, T.cl4x4.energy-T.cluster.energy);
662 
663  float frac24=T.cluster.ET/(T.cl4x4.ET);
664  hA[39]->Fill(frac24);
665  if(frac24<par_clustFrac24) continue;
666 
667  hA[20]->Fill("fr24",1.);
668 
669  //.. spacial separation (track - cluster)
670  TVector3 D=T.pointTower.R-T.cluster.position;
671 
672  hA[43]->Fill( T.cluster.energy,D.Mag());
673  hA[44]->Fill( T.cluster.position.z(),D.z());
674  float delPhi=T.pointTower.R.DeltaPhi(T.cluster.position);
675  // printf("aaa %f %f %f phi=%f\n",D.x(),D.y(),D.z(),delPhi);
676  hA[45]->Fill( T.cluster.energy,Rcylinder*delPhi);// wrong?
677  hA[46]->Fill( D.Mag());
678  hA[199]->Fill(T.cluster.position.PseudoRapidity(),D.Mag());
679  hA[207]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
680 
681  if(D.Mag()>par_delR3D) continue;
682  T.isMatch2Cl=true; // cluster is matched to TPC track
683  hA[20]->Fill("#Delta R",1.);
684  hA[111]->Fill( T.cluster.ET);
685 
686  hA[208]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
687 
688  nTr++;
689  }// end of one vertex
690  }// end of vertex loop
691 
692  if(nTr<=0) return -1;
693  hA[0]->Fill("Tr2Cl",1.0);
694  return 0;
695 
696 }
697 
698 //________________________________________________
699 //________________________________________________
700 WeveCluster
701 St2011WMaker::maxBtow2x2(int iEta, int iPhi, float zVert){
702  //printf(" maxBtow2x2 seed iEta=%d iPhi=%d \n",iEta, iPhi);
703  const int L=2; // size of the summed square
704 
705  WeveCluster maxCL;
706  // just 4 cases of 2x2 clusters
707  float maxET=0;
708  int I0=iEta-1;
709  int J0=iPhi-1;
710  for(int I=I0;I<=I0+1;I++){
711  for(int J=J0;J<=J0+1;J++) {
712  WeveCluster CL=sumBtowPatch(I,J,L,L,zVert);
713  if(maxET>CL.ET) continue;
714  maxET=CL.ET;
715  maxCL=CL;
716  // printf(" newMaxETSum=%.1f iEta=%d iPhi=%d \n",maxET, I,J);
717  }
718  }// 4 combinations done
719  //printf(" final inpEve=%d SumET2x2=%.1f \n",nInpEve,maxET);
720  return maxCL;
721 }
722 
723 
724 //________________________________________________
725 //________________________________________________
726 WeveCluster
727 St2011WMaker::sumBtowPatch(int iEta, int iPhi, int Leta,int Lphi, float zVert){
728  //printf(" eveID=%d btowSquare seed iEta=%d[+%d] iPhi=%d[+%d] zVert=%.0f \n",wEve->id,iEta,Leta, iPhi,Lphi,zVert);
729  WeveCluster CL; // object is small, not to much overhead in creating it
730  CL.iEta=iEta;
731  CL.iPhi=iPhi;
732  TVector3 R;
733  double sumW=0;
734  float Rcylinder= mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
735  for(int i=iEta; i<iEta+Leta;i++){// trim in eta-direction
736  if(i<0) continue;
737  if(i>=mxBTetaBin) continue;
738  for(int j=iPhi;j<iPhi+Lphi;j++) {// wrap up in the phi-direction
739  int jj=(j+mxBTphiBin)%mxBTphiBin;// keep it always positive
740  //if(L<5) printf("n=%2d i=%d jj=%d\n",CL.nTower,i,jj);
741  int softID= mapBtowIJ2ID[ i+ jj*mxBTetaBin];
742  float ene= wEve->bemc.eneTile[kBTow][softID-1];
743  if(ene<=0) continue; // skip towers w/o energy
744  float adc= wEve->bemc.adcTile[kBTow][softID-1];
745  float delZ=positionBtow[softID-1].z()-zVert;
746  float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
747  float ET=ene*e2et;
748  float logET=log10(ET+0.5);
749  CL.nTower++;
750  CL.energy+=ene;
751  CL.ET+=ET;
752  CL.adcSum+=adc;
753  if(logET>0) {
754  R+=logET*positionBtow[softID-1];
755  sumW+=logET;
756  }
757  //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);
758  }
759  // printf(" end btowSquare: iEta=%d nTw=%d, ET=%.1f adc=%.1f\n",i,CL.nTower,CL.ET,CL.adcSum);
760  if(sumW>0) {
761  CL.position=1./sumW*R; // weighted cluster position
762  } else {
763  CL.position=TVector3(0,0,999);
764  }
765  }
766  return CL;
767 }
768 
769 
770 // $Log: St2011W_algo.cxx,v $
771 // Revision 1.22 2016/01/08 02:08:49 jlzhang
772 // added couples histograms and fixed a small bug
773 //
774 // Revision 1.21 2013/09/13 19:33:13 stevens4
775 // Updates to code for combined 2011+2012 result presented to spin PWG 9.12.13
776 //
777 // Revision 1.20 2012/10/05 17:53:53 balewski
778 // added correlation plots for reco Q in Z, W algos
779 //
780 // Revision 1.19 2012/09/28 16:00:42 stevens4
781 // add Q*ET/PT requirement to WB histos used for background estimation to be consistent with spin sorting
782 //
783 // Revision 1.18 2012/09/26 14:20:59 stevens4
784 // use PtBal cos(phi) for WB and WE algos and use Q*ET/PT for barrel charge sign
785 //
786 // Revision 1.17 2012/09/18 22:30:18 stevens4
787 // change to new jet tree format with access to all rank>0 vertices
788 //
789 // Revision 1.16 2012/09/18 21:10:07 stevens4
790 // Include all rank>0 vertex again (new jet format coming next), and remove rank<0 endcap vertices.
791 //
792 // Revision 1.15 2012/09/17 03:29:30 stevens4
793 // Updates to Endcap algo and Q*ET/PT charge separation
794 //
795 // Revision 1.14 2012/08/31 20:10:52 stevens4
796 // switch to second EEMC background using both isolation and sPt-Bal (for mirror symmetry (also adjust eta binning)
797 //
798 // Revision 1.13 2012/08/28 14:28:27 stevens4
799 // add histos for barrel and endcap algos
800 //
801 // Revision 1.12 2012/08/21 21:28:22 stevens4
802 // Add spin sorting for endcap Ws
803 //
804 // Revision 1.11 2012/08/21 18:29:16 stevens4
805 // Updates to endcap W selection using ESMD strip ratio
806 //
807 // Revision 1.10 2012/08/07 21:06:38 stevens4
808 // update to tree analysis to produce independent histos in a TDirectory for each eta-bin
809 //
810 // Revision 1.9 2012/07/13 20:53:16 stevens4
811 // Add filling of empty events in W tree
812 // Minor modifications to histograms
813 //
814 // Revision 1.8 2012/07/06 17:47:02 stevens4
815 // Updates for tree reader
816 //
817 // Revision 1.7 2012/07/05 19:03:54 stevens4
818 // preset bin labels for TPC histos to preven merge warning
819 //
820 // Revision 1.6 2012/06/29 21:20:08 stevens4
821 // *** empty log message ***
822 //
823 // Revision 1.5 2012/06/25 20:53:24 stevens4
824 // algo and histo cleanup
825 //
826 // Revision 1.4 2012/06/18 18:28:01 stevens4
827 // Updates for Run 9+11+12 AL analysis
828 //
829 // Revision 1.3 2011/02/25 06:03:47 stevens4
830 // addes some histos and enabled running on MC
831 //
832 // Revision 1.2 2011/02/17 04:16:22 stevens4
833 // move sector dependent track QA cuts before track pt>10 cut and lower par_clustET and par_ptBalance thresholds to 14 GeV
834 //
835 // Revision 1.1 2011/02/10 20:33:23 balewski
836 // start
837 //
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
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
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::p
const StThreeVectorF & p() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:259
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
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
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
StMuTrack::outerHelix
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Definition: StMuTrack.cxx:412
Generated by   doxygen 1.8.5