StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2011ZMaker.cxx
1 // $Id: St2011ZMaker.cxx,v 1.13 2016/01/08 02:08:49 jlzhang Exp $
2 //
3 //*-- Author : Ross Corliss, MIT
4 // changes Jan Balewski, MIT
5 // changes Justin Stevens, MIT
6 
7 #include <TMath.h>
8 
9 #include "St2011WMaker.h"
10 #include "WeventDisplay.h"
11 #include "St2011ZMaker.h"
12 
13 ClassImp(St2011ZMaker)
14 const float PI=TMath::Pi();
15 
16 //_____________________________________________________________________________
17 //
18 St2011ZMaker::St2011ZMaker(const char *name):StMaker(name){
19  wMK=0;muMK=0;HList=0;
20 
21 }
22 
23 
24 //_____________________________________________________________________________
25 //
26 Int_t St2011ZMaker::Init(){
27  assert(wMK);
28  assert(HList);
29  initHistos();
30  return StMaker::Init();
31 }
32 
33 
34 //_____________________________________________________________________________
35 //
36 Int_t St2011ZMaker::InitRun (int runumber){
37  LOG_INFO<<Form("::InitRun(%d) done, Z-algo params: nearTotEtFrac=%.2f, clusterEt=%.1f GeV, delPhi12>%.2f rad, Zmass in[%.1f,%.1f]\n",
38  runumber, par_nearTotEtFracZ,par_clusterEtZ,par_delPhi12,par_minMassZ,par_maxMassZ)<<endm;
39  return 0;
40 }
41 
42 //_____________________________________________________________________________
43 //
44 Int_t St2011ZMaker::FinishRun (int runnumber){
45 return 0;
46 }
47 
48 //_____________________________________________________________________________
49 //
50 Int_t
51 St2011ZMaker::Make(){
52 
53  //fill various histograms and arrays with event data
54  find_Z_boson();
55  findEndcap_Z_boson();
56 
57  return kStOK;
58 }
59 
60 //============================
61 void
62 St2011ZMaker::printJan(WeveEleTrack *T) {
63  int ibp=kBTow;
64  WevePointTower poiTw=T->pointTower;
65  WeveCluster cl=T->cluster;
66  int id= poiTw.id;
67  float adc= wMK->wEve->bemc.adcTile[ibp][id-1];
68  float frac= adc/4096*60 /cl.ET;
69  printf("Ztower Q=%d pointTw: id=%d ADC=%.0f 2x2ET=%.1f frac=%.2f\n",T->prMuTrack->charge(),id,adc,cl.ET,frac);
70 }
71 
72 
73 //_____________________________________________________________________________
74 //
75 void
76 St2011ZMaker::findEndcap_Z_boson(){
77  Wevent2011 *wEve=wMK->wEve;
78  // printf("========= findEndcap_Z_boson() \n");
79 
80  hA[50]->Fill("inp",1.);
81 
82  // search for Zs ............
83  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
84  hA[50]->Fill("vert",1.);
85  WeveVertex &V=wEve->vertex[iv];
86 
87  //first loop over good barrel tracks
88  for(unsigned int it=0;it<V.eleTrack.size();it++) {
89  WeveEleTrack &TB=V.eleTrack[it];
90  if(TB.pointTower.id<=0) continue; //skip endcap towers
91  if(TB.isMatch2Cl==false) continue;
92  assert(TB.cluster.nTower>0); // internal logical error
93  assert(TB.nearTotET>0); // internal logical error
94 
95  //place cuts on reco track first (both barrel and endcap tracks allowed)
96  float isoET1=TB.cluster.ET /TB.cl4x4.ET;
97  hA[51]->Fill(isoET1);
98  hA[52]->Fill(TB.cluster.ET);
99  hA[50]->Fill("trB",1.);
100  if(TB.cluster.ET<par_clusterEtZ) continue;
101  hA[50]->Fill("etB",1.);
102 
103  float fracET1=TB.cluster.ET /TB.nearTotET;
104  hA[53]->Fill(fracET1);
105  if(fracET1<par_nearTotEtFracZ) continue;
106  hA[50]->Fill("conB",1.);
107 
108  // 1) try to find candidate track in the endcap
109  for(unsigned int it=0;it<V.eleTrack.size();it++) {
110  WeveEleTrack &TE=V.eleTrack[it];
111  if(TE.pointTower.id>=0) continue; //skip barrel towers
112  if(TE.isMatch2Cl==false) continue;
113  assert(TE.cluster.nTower>0); // internal logical error
114  assert(TE.nearTotET>0); // internal logical error
115 
116  float isoET2=TE.cluster.ET/TE.cl4x4.ET;
117  hA[71]->Fill(isoET2);
118  hA[72]->Fill(TE.cluster.ET);
119  hA[70]->Fill("trE",1.);
120  if(TE.cluster.ET<par_clusterEtZ) continue;
121  hA[70]->Fill("etE",1.);
122 
123  float fracET2=TE.cluster.ET/TE.nearTotET;
124  hA[73]->Fill(fracET2);
125  if(fracET2<par_nearTotEtFracZ) continue;
126  hA[70]->Fill("conE",1.);
127 
128  float e1=TB.cluster.energy;
129  float e2=TE.cluster.energy;
130  TVector3 p1=TB.primP; p1.SetMag(e1);//cluster.position;
131  TVector3 p2=TE.primP; p2.SetMag(e2);//cluster.position;
132 
133  float del_phi=p1.DeltaPhi(p2);
134  //printf("del Phi=%f\n",del_phi);
135  float xx=del_phi;
136  if(xx<-PI+1) xx+=2*PI;
137  hA[74]->Fill(xx);
138  if(fabs(del_phi)<par_delPhi12) continue;
139  hA[70]->Fill("phi12",1.);
140 
141  TVector3 psum=p1+p2;
142  float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
143  if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
144  hA[70]->Fill("m2",1.);
145  hA[77]->Fill(p1.Eta(),p2.Eta());
146  hA[78]->Fill(psum.Eta());
147  hA[79]->Fill(psum.Pt());
148 
149  float mass=sqrt(mass2);
150  int Q1Q2=TB.prMuTrack->charge()*TE.prMuTrack->charge();
151  if (Q1Q2==1) { //.. same sign , can't be Z-> e+ e-
152  hA[76]->Fill(mass);
153  hA[80]->Fill(TE.cluster.ET,TE.prMuTrack->charge()/TE.prMuTrack->pt()); continue;
154  }
155 
156  //..... now only opposite sign
157  hA[70]->Fill("QQ",1.);
158  hA[75]->Fill(mass);
159  hA[81]->Fill(TE.cluster.ET,TE.prMuTrack->charge()/TE.prMuTrack->pt());
160  }
161 
162  // 2) use highest ET endcap cluster with no track requirement
163  float maxET=0;
164  WeveCluster maxCluster;
165  for(int iEta=0; iEta<12; iEta++){ //loop over eta bins
166  for(int iPhi=0; iPhi<60; iPhi++){ //loop over phi bins
167 
168  WeveCluster eclust=wMK->maxEtow2x2(iEta,iPhi,V.z);
169  if(eclust.ET < par_clusterEtZ) continue;
170  if(maxET > eclust.ET) continue;
171  else {
172  maxET=eclust.ET;
173  maxCluster=eclust;
174  }
175  }
176  }
177  if(maxCluster.ET <= 1.0) continue;//remove low E clusters
178 
179  //apply cuts to max ETOW cluster and isolation sums
180  WeveCluster cl4x4=wMK->sumEtowPatch(maxCluster.iEta-1,maxCluster.iPhi-1,4,4,V.z);
181  hA[54]->Fill(maxCluster.ET/cl4x4.ET);
182  if(maxCluster.ET/cl4x4.ET<wMK->parE_clustFrac24) continue;
183  hA[55]->Fill(maxCluster.ET);
184  hA[50]->Fill("trE",1.);
185  if(maxCluster.ET < par_clusterEtZ) continue;
186  hA[50]->Fill("etE",1.);
187 
188  //assume poor tracking effic so only towers in nearCone
189  float nearBtow=wMK->sumBtowCone(V.z,maxCluster.position,2);
190  float nearEtow=wMK->sumEtowCone(V.z,maxCluster.position,2);
191  float nearSum=nearBtow; nearSum+=nearEtow;
192  hA[56]->Fill(maxCluster.ET/nearSum);
193  if(maxCluster.ET/nearSum<wMK->parE_nearTotEtFrac) continue;
194  hA[50]->Fill("conE",1.);
195 
196  //add plots of good candidates
197  float e1=TB.cluster.energy;
198  float e2=maxCluster.energy;
199  TVector3 p1=TB.primP; p1.SetMag(e1);
200  TVector3 p2=maxCluster.position; p2.SetMag(e2);
201  float del_phi=p1.DeltaPhi(p2);
202  float xx=del_phi;
203  if(xx<-PI+1) xx+=2*PI;
204  hA[57]->Fill(xx);
205  if(fabs(del_phi)<par_delPhi12) continue;
206  hA[50]->Fill("phi12",1.);
207  TVector3 psum=p1+p2;
208  float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
209  if(mass2<1.) continue;
210  hA[50]->Fill("m2",1.);
211  float mass=sqrt(mass2);
212  hA[58]->Fill(mass);
213  hA[59]->Fill(p1.Eta(),p2.Eta());
214  hA[60]->Fill(psum.Eta());
215  hA[61]->Fill(psum.Pt());
216 
217  } //track loop
218  } //vertex loop
219 }
220 
221 
222 //_____________________________________________________________________________
223 //
224 void
225 St2011ZMaker::find_Z_boson(){
226  Wevent2011 *wEve=wMK->wEve;
227  // printf("========= find_Z_boson() \n");
228 
229  hA[31]->Fill(wEve->vertex.size());
230  hA[0]->Fill("inp",1.);
231 
232  // search for Zs ............
233  for(unsigned int iv=0;iv<wEve->vertex.size();iv++) {
234  hA[0]->Fill("vert",1.);
235  WeveVertex &V=wEve->vertex[iv];
236  hA[32]->Fill(V.eleTrack.size());
237  if(V.eleTrack.size()<2) continue;
238  hA[0]->Fill("TT",1.); // at least 2 isolated tracks exist
239 
240  //only one Z can come from a vertex, and it should be the highest-energy object
241  //hence, the two highest-et clusters should correspond to the z. Pick those
242  //eventually, but for now, just try all of them.
243  for(unsigned int it=0;it<V.eleTrack.size()-1;it++) { //.....select first track:
244  WeveEleTrack &T1=V.eleTrack[it];
245  if(T1.pointTower.id<=0) continue; //skip endcap towers
246  if(T1.isMatch2Cl==false) continue;
247  assert(T1.cluster.nTower>0); // internal logical error
248  assert(T1.nearTotET>0); // internal logical error
249 
250  float isoET1=T1.cluster.ET /T1.cl4x4.ET;
251  hA[29]->Fill(isoET1);
252 
253  hA[23]->Fill(T1.cluster.ET);
254  hA[0]->Fill("tr1",1.);
255  if(T1.cluster.ET<par_clusterEtZ) continue;
256  hA[0]->Fill("et1",1.);
257 
258  float fracET1=T1.cluster.ET /T1.nearTotET;
259  hA[24]->Fill(fracET1);
260  if(fracET1< par_nearTotEtFracZ) continue;
261  hA[0]->Fill("con1",1.);
262 
263  for (unsigned int it2=it+1;it2<V.eleTrack.size();it2++) { //.....select second track:
264  WeveEleTrack &T2=V.eleTrack[it2];
265  if(T2.pointTower.id<=0) continue; //skip endcap towers
266  if(T2.isMatch2Cl==false) continue;
267  assert(T2.cluster.nTower>0); // internal logical error
268  assert(T2.nearTotET>0); // internal logical error
269 
270  float isoET2=T2.cluster.ET /T2.cl4x4.ET;
271  hA[30]->Fill(isoET2);
272 
273  hA[25]->Fill(T2.cluster.ET);
274  hA[0]->Fill("tr2",1.);
275  if(T2.cluster.ET<par_clusterEtZ) continue;
276  hA[0]->Fill("et2",1.);
277 
278  float fracET2=T2.cluster.ET /T2.nearTotET;
279  hA[26]->Fill(fracET2);
280  if(fracET2< par_nearTotEtFracZ) continue;
281  hA[0]->Fill("con2",1.);
282 
283  float e1=T1.cluster.energy;
284  float e2=T2.cluster.energy;
285  TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
286  TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
287 
288  float del_phi=p1.DeltaPhi(p2);
289  //printf("del Phi=%f\n",del_phi);
290  float xx=del_phi;
291  if(xx<-PI+1) xx+=2*PI;
292  hA[27]->Fill(xx);
293  if(fabs(del_phi)<par_delPhi12) continue;
294  hA[0]->Fill("phi12",1.);
295 
296  TVector3 psum=p1+p2;
297  float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
298  float yZ=0.5*log((e1+e2+psum.Z())/(e1+e2-psum.Z()));
299  if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
300  hA[0]->Fill("m2",1.);
301 
302  float mass=sqrt(mass2);
303  int Q1Q2=T1.prMuTrack->charge()*T2.prMuTrack->charge();
304  if (Q1Q2==1) { //.. same sign , can't be Z-> e+ e-
305  hA[14]->Fill(mass);
306  continue;
307  }
308 
309  //..... now only opposite sign
310  hA[0]->Fill("QQ",1.);
311  hA[15]->Fill(mass);
312  hA[33]->Fill(T1.cluster.ET,T1.prMuTrack->charge()/T1.prMuTrack->pt());
313  hA[33]->Fill(T2.cluster.ET,T2.prMuTrack->charge()/T2.prMuTrack->pt());
314  hA[34]->Fill(T1.pointTower.iEta ,T1.cluster.energy);
315  hA[34]->Fill(T2.pointTower.iEta ,T2.cluster.energy);
316  hA[35]->Fill(p1.Eta(),p2.Eta());
317 
318  if(T1.prMuTrack->charge()>0) hA[42]->Fill(p1.Phi(),p2.Phi());
319  else hA[42]->Fill(p2.Phi(),p1.Phi());
320  hA[43]->Fill(T1.cluster.ET,T1.prMuTrack->charge()*T1.cluster.ET/T1.prMuTrack->pt());
321  hA[43]->Fill(T2.cluster.ET,T2.prMuTrack->charge()*T2.cluster.ET/T2.prMuTrack->pt());
322 
323 #if 0
324  printf("RCC: Found Z w/ invmass=%f\n",mass);
325  printJan(&T1);
326  printJan(&T2);
327 
328 
329  if (!wMK->isMC || (wMK->isMC&& wEve->id<500) )
330  { printf("\n ZZZZZZZZZZZZZZZZZZZ\n");
331  if(mass<par_minMassZ)
332  wMK->wDisaply->exportEvent("Zlow",V,T1);
333  else
334  wMK->wDisaply->exportEvent("Zgood",V,T1);
335  printf("RCC: Found Z w/ invmass=%f\n",mass);
336  wEve->print();
337  }
338 
339 #endif
340 
341  if (mass<par_minMassZ) continue; //enforce a lower bound
342  hA[0]->Fill("Zlow",1.);
343 
344  if (mass>par_maxMassZ) continue; //enforce an upper bound
345  hA[0]->Fill("Zhigh",1.);
346 
347  hA[36]->Fill(yZ);//psum.Eta());
348  hA[37]->Fill(psum.Pt());
349 
350  int bxStar7=wMK->wEve->bxStar7;
351  int bxStar48=wMK->wEve->bxStar48;
352  if(bxStar48!=bxStar7) {
353  hA[0]->Fill("badBx48",1.);
354  return; // both counters must be in sync
355  }
356 
357  int spin4=wMK->wEve->spin4;
358  hA[38]->Fill(spin4);
359  if(yZ<0) hA[39]->Fill(spin4);
360  else if(yZ>0) hA[40]->Fill(spin4);
361 
362  // L0 x1 and x2 computation
363  float mZ = 91.188;
364  float roots = 510.;
365  float x1 = mZ/roots * TMath::Exp(yZ);
366  float x2 = mZ/roots * TMath::Exp(-1.*yZ);
367  hA[44]->Fill(x1,x2);
368  hA[45]->Fill(x1*mass/mZ,x2*mass/mZ);
369 
370  // free quark search
371  if(T1.prMuTrack->charge()>0)
372  hA[46]->Fill( T1.prMuTrack->charge()*T1.cluster.ET/T1.prMuTrack->pt(),
373  T2.prMuTrack->charge()*T2.cluster.ET/T2.prMuTrack->pt());
374  else
375  hA[46]->Fill( T2.prMuTrack->charge()*T2.cluster.ET/T2.prMuTrack->pt(),
376  T1.prMuTrack->charge()*T1.cluster.ET/T1.prMuTrack->pt());
377 
378  // **** I stoped changes here, Jan
379 
380  float fmax1=T1.cluster.ET/T1.cl4x4.ET;
381  float fmax2=T2.cluster.ET/T2.cl4x4.ET;
382 
383  hA[21]->Fill(fmax1,fmax2);
384  hA[22]->Fill(T1.cluster.ET,T2.cluster.ET);
385 
386  hA[1]->Fill(mass);
387  hA[2]->Fill(T1.prMuTrack->charge(),T2.prMuTrack->charge());
388  hA[3]->Fill(T1.prMuTrack->charge()*T2.prMuTrack->charge());
389  hA[4]->Fill(p1.Phi(),p2.Phi());
390  hA[5]->Fill(del_phi);
391  hA[6]->Fill(mass,T1.prMuTrack->charge()/T1.primP.Perp()*T2.prMuTrack->charge()/T1.primP.Perp());
392  hA[7]->Fill(mass,T1.prMuTrack->charge()*T2.prMuTrack->charge());
393  hA[8]->Fill(T1.cluster.ET);
394  if (T1.prMuTrack->charge()>0)
395  {
396  hA[9]->Fill(p1.Eta(),p1.Phi());
397  hA[10]->Fill(p2.Eta(),p2.Phi());
398  }
399  else
400  {
401  hA[10]->Fill(p1.Eta(),p1.Phi());
402  hA[9]->Fill(p2.Eta(),p2.Phi());
403  }
404 
405  hA[11]->Fill(fmax1,fmax2);
406  hA[12]->Fill(T1.cluster.ET,T2.cluster.ET);
407  hA[13]->Fill(mass,del_phi);
408 
409  }
410 
411 
412 
413  }// loop over first track
414  }// loop over vertices
415 
416 }
417 
418 
419 // $Log: St2011ZMaker.cxx,v $
420 // Revision 1.13 2016/01/08 02:08:49 jlzhang
421 // added couples histograms and fixed a small bug
422 //
423 // Revision 1.12 2012/10/05 17:53:53 balewski
424 // added correlation plots for reco Q in Z, W algos
425 //
426 // Revision 1.11 2012/10/05 16:44:31 stevens4
427 // final z plots: update x1-x2 correlation in z maker
428 //
429 // Revision 1.10 2012/10/01 19:48:20 stevens4
430 // add plots for Z result and move esmd cross point calculation outside plane loop
431 //
432 // Revision 1.9 2012/09/26 14:20:59 stevens4
433 // use PtBal cos(phi) for WB and WE algos and use Q*ET/PT for barrel charge sign
434 //
435 // Revision 1.8 2012/09/24 19:28:15 balewski
436 // moved eta & pT hitso to be filled after invM cut, now they are 'golden Z'
437 //
438 // Revision 1.7 2012/09/18 19:34:22 balewski
439 // fill eta-dependent spin sort
440 //
441 // Revision 1.6 2012/09/14 21:02:29 balewski
442 // *lumi-maker re-written to accumulate alternative rel lumi monitors,
443 // * added spin sorting to Zs
444 //
445 // Revision 1.5 2012/08/21 17:40:09 stevens4
446 // Revert to previous version
447 //
448 // Revision 1.3 2012/06/29 21:20:08 stevens4
449 // *** empty log message ***
450 //
451 // Revision 1.2 2012/06/26 20:30:23 stevens4
452 // Updates ZMaker for mixing barrel and endcap arms
453 //
454 // Revision 1.1 2011/02/10 20:33:24 balewski
455 // start
456 //
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
St2011ZMaker::Make
virtual Int_t Make()
Definition: St2011ZMaker.cxx:51
kStOK
Definition: Stypes.h:40
St2011ZMaker
uses tree from W-algo to find Zs
Definition: St2011ZMaker.h:26
Generated by   doxygen 1.8.5