StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2009pubMcMaker.cxx
1 // $Id: St2009pubMcMaker.cxx,v 1.7 2011/09/14 14:23:21 stevens4 Exp $
2 //
3 //*-- Author : Justin Stevens, IUCF
4 //
5 
6 #include "St2009WMaker.h"
7 #include "St2009ZMaker.h"
8 #include "St2009pubMcMaker.h"
9 #include "StEmcUtil/geometry/StEmcGeom.h"
10 #include <TMath.h>
11 
12 //need these to get MC record
13 #include "tables/St_g2t_tpc_hit_Table.h"
14 #include "StMcEventMaker/StMcEventMaker.h"
15 #include "StMcEvent/StMcEvent.hh"
16 #include "StMcEvent/StMcVertex.hh"
17 #include "StMcEvent/StMcTrack.hh"
18 
19 //muDst needed for zdcRate
20 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
21 #include <StMuDSTMaker/COMMON/StMuDst.h>
22 #include <StMuDSTMaker/COMMON/StMuEvent.h>
23 #include <StMuDSTMaker/COMMON/StMuMcTrack.h>
24 
25 ClassImp(St2009pubMcMaker)
26 
27 //_______________________________________________________
28 //_______________________________________________________
29 St2009pubMcMaker::St2009pubMcMaker(const char *name):StMaker(name){
30  wMK=0;HList=0;
31 
32 }
33 
34 
35 //_______________________________________________________
36 //_______________________________________________________
37 St2009pubMcMaker::~St2009pubMcMaker(){
38  //
39 }
40 
41 
42 //_______________________________________________________
43 //_______________________________________________________
44 Int_t St2009pubMcMaker::Init(){
45  assert(wMK);
46  assert(HList);
47  initHistos();
48  return StMaker::Init();
49 }
50 
51 
52 //_______________________________________________________
53 //_______________________________________________________
54 Int_t
55 St2009pubMcMaker::Make(){
56  //printf("-----------in %s\n", GetName());
57  //only get geant particle info for W MC
58  if(wMK->isMC==350){
59  if(doWMCanalysis()){
60  doWefficiency();
61  doWanalysis();
62  }
63  }
64  if(wMK->isMC==352){
65  if(doZMCanalysis()){
66  doZefficiency();
67  }
68  }
69  return kStOK;
70 }
71 
72 
73 //_______________________________________________________
74 //_______________________________________________________
75 void
76 St2009pubMcMaker::doWefficiency(){
77 
78  //initialize some kinematic variables for later use
79 
80  //wrap barrel onto 2 modules to see boundaries
81  const float PI=TMath::Pi();
82  float elePhi=mElectronP.Phi();
83  if(elePhi<0) elePhi+=2*PI;
84  elePhi=elePhi*180/PI;
85  elePhi+=3.0;
86  float modulePair=12.0;
87  int resid = (int) (elePhi/modulePair);
88  float phiMod= elePhi - resid*modulePair;
89  // 'detector' eta
90  TVector3 detEle; //where lepton would hit BEMC
91  float Rcylinder= wMK->mBtowGeom->Radius();
92  detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
93  detEle.SetZ(detEle.Z()+mVertex.Z());
94 
95  float w=1; //determine zVertex weighting for event
96  if(wMK->isMC ==350) w=wMK->hReweight->GetBinContent(wMK->hReweight->FindBin(mVertex.Z()));
97  float zdcRate=wMK->mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
98 
99  //only count leptons in our eta range
100  if(fabs(mElectronP.Eta()) > 1.0) return;
101  //only count leptons in our ET range
102  if(fabs(mElectronP.Perp()) < 25.) return;
103 
104  int bx7=wMK->wEve.bx7;
105  if( (bx7>30 && bx7<40) || (bx7>110) ) return; //skip abort gaps
106 
107  //ele has |eta| < 1 and ET > 25
108  hA[50]->Fill(mElectronP.Perp(),w);
109  hA[54]->Fill(mElectronP.Eta(),w);
110  hA[58]->Fill(mVertex.Z(),w);
111  hA[62]->Fill(mElectronP.Phi(),w);
112  hB[83]->Fill(mElectronP.Eta(),mElectronP.Perp(),w);
113  hA[89]->Fill(zdcRate,w);
114 
115  hA[68]->Fill(mElectronP.Perp(),w); //forJoe
116  hA[77]->Fill(phiMod,w); //wrap barrel onto 2 modules
117  hB[66]->Fill(detEle.Eta(),mElectronP.Perp(),w);
118  hA[79]->Fill(detEle.Eta(),w);
119  hA[122]->Fill(zdcRate,wMK->wEve.bx7);
120 
121  //check reconstructed vs thrown ET
122  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
123  WeveVertex &V=wMK->wEve.vertex[iv];
124  for(unsigned int it=0;it<V.eleTrack.size();it++) {
125  WeveEleTrack &T=V.eleTrack[it];
126  if(T.isMatch2Cl==false) continue;
127  assert(T.cluster.nTower>0); // internal logical error
128  assert(T.nearTotET>0); // internal logical error
129 
130  hB[68]->Fill(mElectronP.Perp(),T.cluster.ET,w);
131  }
132  }
133 
134  //trigger efficiency
135  if(!wMK->wEve.l2bitET) return;
136  //good trig
137  hA[51]->Fill(mElectronP.Perp(),w);
138  hA[55]->Fill(mElectronP.Eta(),w);
139  hA[59]->Fill(mVertex.Z(),w);
140  hA[63]->Fill(mElectronP.Phi(),w);
141  hB[84]->Fill(mElectronP.Eta(),mElectronP.Perp(),w);
142  hA[90]->Fill(zdcRate,w);
143 
144  hA[69]->Fill(mElectronP.Perp(),w);
145  hA[78]->Fill(phiMod,w);//wrap barrel onto 2 modules
146  hB[67]->Fill(detEle.Eta(),mElectronP.Perp(),w);
147  hA[80]->Fill(detEle.Eta(),w);
148 
149  //vertex efficiency
150  //require one vertex w/ rank>0 and |z|<100
151  if(wMK->wEve.vertex.size()<=0) return;
152  //require geant vert match reco vert (w/in 5 cm)
153  bool matchVert=false;
154  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
155  WeveVertex &V=wMK->wEve.vertex[iv];
156  //cout<<endl<<"vertex diff "<<fabs(V.z - mVertex.Z())<<endl;
157  if(fabs(V.z - mVertex.Z()) < 5.0) matchVert=true;
158  }
159  if(!matchVert) return;
160  //cout<<"after matched vert"<<endl;
161  hA[52]->Fill(mElectronP.Perp(),w);
162  hA[56]->Fill(mElectronP.Eta(),w);
163  hA[60]->Fill(mVertex.Z(),w);
164  hA[64]->Fill(mElectronP.Phi(),w);
165  hA[91]->Fill(zdcRate,w);
166 
167  hA[70]->Fill(mElectronP.Perp(),w);//forJoe
168 
169  //get MC track info
170  int eleTrId=999; TVector3 eleTrP;
171  int size=wMK->mMuDstMaker->muDst()->mcArray(1)->GetEntries();
172  cout<<"size="<<size<<endl;
173  for(int i=0; i<size; i++){
174  StMuMcTrack* muMcTrack=(StMuMcTrack*)wMK->mMuDstMaker->muDst()->mcArray(1)->UncheckedAt(i);
175  int geantId=muMcTrack->GePid();
176  int id=muMcTrack->Id();
177  if(geantId==2 || geantId==3){
178  cout<<i<<" id="<<id<<" geantId="<<geantId<<endl;
179  eleTrId=id;
180  eleTrP=TVector3(muMcTrack->Pxyz().x(),muMcTrack->Pxyz().y(),muMcTrack->Pxyz().z());
181  cout<<"phiDiff="<<eleTrP.Phi()-mElectronP.Phi()<<" etaDiff="<<eleTrP.Eta()-mElectronP.Eta()<<endl;
182  break;
183  }
184  }
185 
186  int nTrackAssoc=0; float recoPtAssoc=0;
187  int flagAssoc=0;
188  int nHitAssoc=0; float hitFracAssoc=0.;
189  float rInAssoc=0.; float rOutAssoc=0.;
190  //check tracking efficiency before any track cuts
191  for(unsigned int iv=0;iv<wMK->wEve.vertex.size(); iv++) {
192  unsigned int vertID=wMK->wEve.vertex[iv].id;
193  wMK->mMuDstMaker->muDst()->setVertexIndex(vertID);
194  Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
195  for(int itr=0;itr<nPrimTrAll;itr++) {
196  StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
197  if(prTr->flag()<=0) continue;
198  const StMuTrack *glTr=prTr->globalTrack();
199  if(glTr==0) continue; // see the reason in _accessMuDst
200  if(prTr->flag()!=301) continue;// TPC+prim vertex tracks
201  if(prTr->pt()<10.) continue;
202 
203  StThreeVectorF ro=glTr->lastPoint();
204  int secID=WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity());
205  if(secID==20) continue;
206 
207  //only use with private copy of StMuTrack which makes these functions public (currently in dev they're private)
208  //int idTruth=prTr->idTruth();
209  //int qaTruth=prTr->qaTruth();
210 
211  int idTruth=0;
212  int qaTruth=0;
213 
214  float delPhi=prTr->phi()-mElectronP.Phi();
215  float delEta=prTr->eta()-mElectronP.Eta();
216  if( eleTrId==idTruth ) { //reco track matches to electron
217 
218  flagAssoc=prTr->flag();
219  nHitAssoc=prTr->nHitsFit();
220  hitFracAssoc=1.*prTr->nHitsFit()/prTr->nHitsPoss();
221  rInAssoc=glTr->firstPoint().perp();
222  rOutAssoc=glTr->lastPoint().perp();
223 
224  cout<<"found electron, delR="<<sqrt(delPhi*delPhi+delEta*delEta)<<" nHits="<<prTr->nHitsFit()<<" fitFrac="<<hitFracAssoc<<endl;
225  cout<<"idTruth="<<idTruth<<" qaTruth="<<qaTruth<<endl;
226  recoPtAssoc=prTr->pt();
227  nTrackAssoc++;
228  }
229  }
230  }
231  if(nTrackAssoc==1) {
232  hA[120]->Fill(recoPtAssoc,recoPtAssoc-mElectronP.Perp());
233  hA[124]->Fill(recoPtAssoc,w);
234  hA[126]->Fill(1./recoPtAssoc,w);
235  hA[128]->Fill(recoPtAssoc,nHitAssoc);
236  hA[129]->Fill(recoPtAssoc,hitFracAssoc);
237  hA[130]->Fill(recoPtAssoc,rInAssoc);
238  hA[131]->Fill(recoPtAssoc,rOutAssoc);
239  hA[132]->Fill(1./recoPtAssoc,nHitAssoc);
240  hA[133]->Fill(1./recoPtAssoc,hitFracAssoc);
241  hA[134]->Fill(1./recoPtAssoc,rInAssoc);
242  hA[135]->Fill(1./recoPtAssoc,rOutAssoc);
243  hA[136]->Fill(flagAssoc);
244  }
245  cout<<"number of associated tracks = "<<nTrackAssoc<<endl;
246 
247  // track efficiency
248  bool foundTrack=false; int nTrk=0; float recoPt=0;
249  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
250  WeveVertex &V=wMK->wEve.vertex[iv];
251  if(V.eleTrack.size()>0){ //track with pt > 10
252  hA[85]->Fill(mElectronP.Perp(),w);
253  hA[86]->Fill(mElectronP.Eta(),w);
254  hA[87]->Fill(mVertex.Z(),w);
255  hA[88]->Fill(mElectronP.Phi(),w);
256  hA[92]->Fill(zdcRate,w);
257  foundTrack=true;
258 
259  //find track
260  for(unsigned int it=0;it<V.eleTrack.size();it++) {
261  WeveEleTrack &T=V.eleTrack[it];
262  int idTruth=0; //T.prMuTrack->idTruth(); //again need private StMuTrack
263  if( eleTrId==idTruth ) { //reco track matches to electron
264  recoPt=T.primP.Perp();
265  nTrk++;
266  }
267  }
268  break; //only count event once
269  }
270  }
271  if(!foundTrack) return;
272 
273  if(nTrackAssoc==1) { //accepted tracks
274  hA[125]->Fill(recoPtAssoc,w);
275  hA[127]->Fill(1./recoPtAssoc,w);
276  }
277 
278  if(wMK->wEve.wTag) hA[140]->Fill(mElectronP.Perp(),w);
279 
280  //reco efficiency
281  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
282  WeveVertex &V=wMK->wEve.vertex[iv];
283  for(unsigned int it=0;it<V.eleTrack.size();it++) {
284  WeveEleTrack &T=V.eleTrack[it];
285  if(T.isMatch2Cl==false) continue;
286  assert(T.cluster.nTower>0); // internal logical error
287  assert(T.nearTotET>0); // internal logical error
288 
289  if(wMK->wEve.zTag)
290  continue; //event tagged as Z
291  if(fabs(T.primP.Eta()) > wMK->par_leptonEta)
292  continue; // reconstructed eta too large
293 
294  if(T.cluster.ET/T.nearTotET < wMK->par_nearTotEtFrac)
295  continue; // too large nearET
296  if(T.sPtBalance < wMK->par_ptBalance)
297  continue; // too large signed pt balance
298 
299  if(T.cluster.ET < 25.)
300  continue; // too small reco ET (use same cut as data)
301 
302  //pass all W cuts
303  hA[53]->Fill(mElectronP.Perp(),w);
304  hA[57]->Fill(mElectronP.Eta(),w);
305  hA[61]->Fill(mVertex.Z(),w);
306  hA[65]->Fill(mElectronP.Phi(),w);
307  hA[93]->Fill(zdcRate,w);
308 
309  hA[71]->Fill(mElectronP.Perp(),w);//forJoe
310  hA[121]->Fill(T.primP.Phi()-mElectronP.Phi(),T.primP.Eta()-mElectronP.Eta());
311  hA[123]->Fill(zdcRate,wMK->wEve.bx7);
312 
313  //fill thrown ele ET after all W cuts
314  hA[111]->Fill(detEle.Eta(),mElectronP.Perp());
315  hA[113]->Fill(detEle.Eta(),mElectronSmearP.Perp());
316  hA[115]->Fill(T.cluster.ET/mElectronP.Perp());
317  hA[116]->Fill(T.cluster.energy/mElectronP.Mag());
318 
319  //template function stuff
320  for(int ires=0; ires<10; ires++){
321  for(int iscale=0; iscale<101; iscale++){
322  float scaleFact = iscale*0.004 + 0.8;
323  hB[120+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
324  hB[150+ires]->Fill(scaleFact,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
325  if(fabs(detEle.Eta()) < 0.5) hB[130+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
326  else if(fabs(detEle.Eta()) < 1.0) hB[140+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
327  }
328  }
329 
330  break; //only count event once
331  }
332  }
333 
334 }
335 
336 
337 //_______________________________________________________
338 //_______________________________________________________
339 void
340 St2009pubMcMaker::doZefficiency(){
341 
342  //Acceptance cuts applied
343  //only count leptons in our eta range
344  if(fabs(mZelectronP.Eta()) > 1.) return;
345  if(fabs(mZpositronP.Eta()) > 1.) return;
346  //only count leptons in our ET range
347  if(mZelectronP.Perp() < 15.) return;
348  if(mZpositronP.Perp() < 15.) return;
349 
350  TVector3 sumP=mZelectronP+mZpositronP;
351  float sumE=mZelectronP.Mag()+mZpositronP.Mag();
352  float geantZmass=sqrt(sumE*sumE-sumP.Dot(sumP));
353  cout<<"geant mass"<<geantZmass<<endl;
354  float zdcRate=wMK->mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
355 
356  float w=1; //determine zVertex weighting for event
357  if(wMK->isMC==352) w=wMK->hReweight->GetBinContent(wMK->hReweight->FindBin(mVertex.Z()));
358 
359  //fill all events in acceptance
360  hA[100]->Fill(geantZmass,w);
361  if(geantZmass>70. && geantZmass<110.) hA[105]->Fill(zdcRate,w);
362 
363  //trigger efficiency
364  if(!wMK->wEve.l2bitET) return;
365  //fill good trig events
366  hA[101]->Fill(geantZmass,w);
367  if(geantZmass>70. && geantZmass<110.) hA[106]->Fill(zdcRate,w);
368 
369  //vertex efficiency
370  if(wMK->wEve.vertex.size()<=0) return;
371  //fill vertex rank>0 and |z|<100 events
372  bool matchVert=false;
373  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
374  WeveVertex &V=wMK->wEve.vertex[iv];
375  //cout<<endl<<"vertex diff "<<fabs(V.z - mVertex.Z())<<endl;
376  if(fabs(V.z - mZvertex.Z()) < 5.0) matchVert=true;
377  }
378  if(!matchVert) return;
379  hA[102]->Fill(geantZmass,w);
380  if(geantZmass>70. && geantZmass<110.) hA[107]->Fill(zdcRate,w);
381 
382  //track efficiency
383  bool gt1track=false;
384  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
385  WeveVertex &V=wMK->wEve.vertex[iv];
386  if(V.eleTrack.size()>1){ //fill 2 track with pt > 10 events
387  hA[103]->Fill(geantZmass,w);
388  if(geantZmass>70. && geantZmass<110.) hA[108]->Fill(zdcRate,w);
389 
390  gt1track=true;
391  break; //only count event once
392  }
393  }
394  if(!gt1track) return;
395 
396  //reco efficiency
397  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
398  WeveVertex &V=wMK->wEve.vertex[iv];
399  if(V.eleTrack.size()<2) continue;
400  for(unsigned int it=0;it<V.eleTrack.size()-1;it++) {//.....select first track:
401 
402  WeveEleTrack &T1=V.eleTrack[it];
403  if(T1.isMatch2Cl==false) continue;
404  assert(T1.cluster.nTower>0); // internal logical error
405  assert(T1.nearTotET>0); // internal logical error
406 
407  if(T1.cluster.ET<zMK->par_clusterEtZ) continue;
408  float fracET1=T1.cluster.ET /T1.nearTotET;
409  if(fracET1<zMK->par_nearTotEtFracZ) continue;
410 
411  for (unsigned int it2=it+1;it2<V.eleTrack.size();it2++) { //.....select second track:
412 
413  WeveEleTrack &T2=V.eleTrack[it2];
414  if(T2.isMatch2Cl==false) continue;
415  assert(T2.cluster.nTower>0); // internal logical error
416  assert(T2.nearTotET>0); // internal logical error
417 
418  if(T2.cluster.ET<zMK->par_clusterEtZ) continue;
419  float fracET2=T2.cluster.ET /T2.nearTotET;
420  if(fracET2<zMK->par_nearTotEtFracZ) continue;
421 
422  float e1=T1.cluster.energy;
423  float e2=T2.cluster.energy;
424  TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
425  TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
426 
427  float del_phi=p1.DeltaPhi(p2);
428  if(fabs(del_phi)<zMK->par_delPhi12) continue;
429  TVector3 psum=p1+p2;
430  float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
431  if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
432  //float mass=sqrt(mass2);
433 
434  int Q1Q2=T1.prMuTrack->charge()*T2.prMuTrack->charge();
435  if (Q1Q2==1) continue;
436 
437  //fill reco Z histos
438  hA[104]->Fill(geantZmass,w);
439  if(geantZmass>70. && geantZmass<110.) hA[109]->Fill(zdcRate,w);
440 
441  } //second track
442  } //first track
443  } //vertex
444 
445  return;
446 }
447 
448 
449 //________________________________________________
450 //________________________________________________
451 bool
452 St2009pubMcMaker::doWMCanalysis(){
453  StMcEvent* mMcEvent = 0;
454  mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
455  //assert(mMcEvent);
456 
457  if(!mMcEvent) return false;
458 
459  //initialize momentum vectors
460  StThreeVectorF pW; float eW;
461  StThreeVectorF pNeutrino; //float eNeutrino;
462  StThreeVectorF pElectron; //float eElectron;
463 
464  StMcVertex *V=mMcEvent->primaryVertex();
465  mVertex=TVector3(V->position().x(),V->position().y(),V->position().z());
466 
467  //find geant tracks for W and decay daughters
468  unsigned int i=1;
469  int found=0;
470  while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
471  StMcTrack* mcTrack = mMcEvent->tracks()[i];
472  int pdgId=mcTrack->pdgId();
473  float pt=mcTrack->pt();
474  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
475  if(pdgId==11 || pdgId==-11){ //select e+ and e-
476  if(abs(mcTrack->parent()->pdgId()) == 24 ){
477  pElectron=mcTrack->momentum();
478  //mKeyElectron=mcTrack->key();
479  //mEveGenElectron=mcTrack->eventGenLabel();
480  //mIdElectron=i;
481  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
482  pW=mcTrack->parent()->momentum();
483  eW=mcTrack->parent()->energy();
484  LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
485  found++;
486  }
487  }
488  if(pdgId==12 || pdgId==-12){ //select neutrino
489  if(abs(mcTrack->parent()->pdgId()) == 24 ){
490  pNeutrino=mcTrack->momentum();
491  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
492  pW=mcTrack->parent()->momentum();
493  eW=mcTrack->parent()->energy();
494  LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
495  found++;
496  }
497  }
498  i++;
499  }
500 
501  i=1;
502  while(i<mMcEvent->tracks().size()){
503  StMcTrack* mcTrack = mMcEvent->tracks()[i];
504  int pdgId=mcTrack->pdgId();
505  float pt=mcTrack->pt();
506  i++;
507  if(abs(pdgId)==11 || abs(pdgId)==12 || abs(pdgId)==24 || abs(pdgId)==21 || abs(pdgId) < 10 || abs(pdgId)==92 || pt<10.0) continue;
508  cout<<"high pt MC track: pdgId="<<pdgId<<" pt="<<pt<<endl;
509 
510  }
511 
512  if(found!=2) return false;
513 
514  mWP=TVector3(pW.x(),pW.y(),pW.z());
515  mNeutrinoP=TVector3(pNeutrino.x(),pNeutrino.y(),pNeutrino.z());
516  mElectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
517  TVector3 diff=mWP-mNeutrinoP-mElectronP;
518  if(diff.Mag() > 0.0001) //should get exactly right
519  LOG_INFO<<"\n \n W+e+nu vector sum ="<<diff.Mag()<<endm;
520  if(mElectronP.Mag() < 0.0001)
521  LOG_INFO<<"\n \n no lepton track ="<<endm;
522 
523  //calculate x1 and x2 from W rapidity
524  float rapW = 0.5*log((eW+mWP.Z())/(eW-mWP.Z()));
525  float mw2sqs=80.4/500.;
526  float x1=mw2sqs*exp(rapW);
527  float x2=mw2sqs*exp(-rapW);
528 
529  hA[72]->Fill(rapW);
530  if(fabs(mElectronP.Eta())<1){ //require midrapidity e
531  hA[73]->Fill(x1);
532  hA[74]->Fill(x2);
533  hA[75]->Fill(x1,x2);
534  hA[76]->Fill(x1-x2);
535  }
536 
537  //original electron Pt spectrum
538  if(fabs(mElectronP.Eta()) < 1.)
539  hA[117]->Fill(mElectronP.Pt());
540 
541  //define weighting for rhicbos
542  wRB = 0.0; //init weight at 0
543  float weight[120] = {0,0,0,0,0,0,0,0,0,0,0,0.0179866,0,0,0,1.02642,0.597259,3.00388,0.88189,1.68209,0.877086,0.787812,0.920501,0.737976,0.735082,0.977448,1.21798,0.745058,0.693323,1.02895,1.69843,1.17129,0.945612,0.917748,0.977399,0.87226,1.0391,0.86036,0.769629,1.03845,0.923467,1.06595,0.843156,1.05421,0.859137,0.896914,0.941134,0.826568,0.857237,0.915694,0.892815,0.901463,1.07729,1.00146,0.935025,1.01036,1.03116,0.873262,1.04329,0.914255,0.976182,0.88009,1.02751,1.12671,1.07694,1.10049,1.04973,1.11188,1.10395,0.914941,1.08072,1.11682,1.0614,1.11588,0.960717,0.992878,1.05434,1.0802,0.967019,1.02817,1.04221,1.01508,1.0033,1.2565,0.987695,1.17665,1.05089,1.00209,0.97153,1.07144,0.824123,1.62538,0.883172,0.860085,0.860167,1.26625,0.58652,0.925282,0.998732,1.06802,0.620375,0.686918,0.845012,0.863964,0.966573,0.52009,0.974503,0.603462,0.305834,0.220393,0.347469,0.849899,0.674756,0.630642,0.613394,0.173667,0.152909,0.325029,0.200083,0.316456};
544 
545  float et = mElectronP.Pt();
546  for(int j=0; j<120; j++){
547  if(et > j*0.5 && et < (j+1)*0.5 && fabs(mElectronP.Eta()) < 1.){
548  hA[118]->Fill(et,weight[j]);
549  wRB = weight[j];
550  break;
551  }
552  }
553 
554  wRB=1.0; // drop weighting for now
555 
556  //smear electron energy by BTOW resolution
557  float initE=mElectronP.Mag();
558  float width=sqrt(0.14*0.14/initE + 0.015*0.015); //resolution factor from NIM
559  float smear = gRandom->Gaus(0, width);
560  float finalE=mElectronP.Mag()+smear;
561  cout<<"initE="<<initE<<" width="<<width<<" smear="<<smear<<" finalE="<<finalE<<endl;
562  mElectronSmearP=mElectronP; mElectronSmearP.SetMag(finalE);
563 
564  //where lepton would hit BEMC
565  TVector3 detEle;
566  float Rcylinder= wMK->mBtowGeom->Radius();
567  detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
568  detEle.SetZ(detEle.Z()+mVertex.Z());
569 
570  //thrown ET for jacob peak comp
571  hA[110]->Fill(detEle.Eta(),mElectronP.Perp());
572  //smeared ET
573  hA[112]->Fill(detEle.Eta(),mElectronSmearP.Perp());
574 
575  //template function stuff
576  for(int ires=0; ires<10; ires++){
577  float widthTemp = ires/100.;
578  float smearTemp = gRandom->Gaus(0, widthTemp);
579  float finalETemp = mElectronP.Mag()+smearTemp*mElectronP.Mag();
580  //cout<<"initE="<<initE<<" width="<<width<<" smear="<<smear<<" finalE="<<finalE<<endl;
581  mElectronSmearTempP[ires]=mElectronP; mElectronSmearTempP[ires].SetMag(finalETemp);
582  for(int iscale=0; iscale<101; iscale++){
583  float scaleFact = iscale*0.004 +0.9;
584  hB[110+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact);
585  }
586  }
587 
588  return true;
589 
590 }
591 
592 
593 
594 //________________________________________________
595 //________________________________________________
596 bool
597 St2009pubMcMaker::doZMCanalysis(){
598  StMcEvent* mMcEvent = 0;
599  mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
600  //assert(mMcEvent);
601 
602  if(!mMcEvent) return false;
603 
604  //initialize momentum vectors
605  StThreeVectorF pZ; float eZ;
606  StThreeVectorF pPositron; //float eNeutrino;
607  StThreeVectorF pElectron; //float eElectron;
608 
609  StMcVertex *V=mMcEvent->primaryVertex();
610  mZvertex=TVector3(V->position().x(),V->position().y(),V->position().z());
611 
612  //find geant tracks for Z and e+,e-
613  unsigned int i=1;
614  int found=0;
615  while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
616  StMcTrack* mcTrack = mMcEvent->tracks()[i];
617  int pdgId=mcTrack->pdgId();
618  float pt=mcTrack->pt();
619  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
620  if(pdgId==11){ //select e-
621  if(abs(mcTrack->parent()->pdgId()) == 23 ){
622  pElectron=mcTrack->momentum();
623  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
624  pZ=mcTrack->parent()->momentum();
625  eZ=mcTrack->parent()->energy();
626  LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
627  found++;
628  }
629  }
630  if(pdgId==-11){ //select e+
631  if(abs(mcTrack->parent()->pdgId()) == 23 ){
632  pPositron=mcTrack->momentum();
633  LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
634  pZ=mcTrack->parent()->momentum();
635  eZ=mcTrack->parent()->energy();
636  LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
637  found++;
638  }
639  }
640  i++;
641  }
642 
643  if(found!=2) return false;
644 
645  mZP=TVector3(pZ.x(),pZ.y(),pZ.z());
646  mZpositronP=TVector3(pPositron.x(),pPositron.y(),pPositron.z());
647  mZelectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
648  TVector3 diff=mZP-mZpositronP-mZelectronP;
649  if(diff.Mag() > 0.0001) //should get exactly right
650  LOG_INFO<<"\n \n Z+e+e vector sum ="<<diff.Mag()<<endm;
651  if(mZelectronP.Mag() < 0.0001)
652  LOG_INFO<<"\n \n no lepton track ="<<endm;
653 
654  return true;
655 }
656 
657 
658 //_______________________________________________________
659 //_______________________________________________________
660 void
661 St2009pubMcMaker::doWanalysis(){
662  //has access to whole W-algo-maker data via pointer 'wMK'
663 
664  // run through W cuts to fill other histos............
665  for(unsigned int iv=0;iv<wMK->wEve.vertex.size();iv++) {
666  WeveVertex &V=wMK->wEve.vertex[iv];
667  for(unsigned int it=0;it<V.eleTrack.size();it++) {
668  WeveEleTrack &T=V.eleTrack[it];
669  if(T.isMatch2Cl==false) continue;
670  assert(T.cluster.nTower>0); // internal logical error
671  assert(T.nearTotET>0); // internal logical error
672 
673  if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
674  if(T.awayTotET> 30.) continue; // too large awayET , Jan
675  //Full W cuts applied at this point
676 
677  //W info from pythia record
678  hA[1]->Fill(mWP.Perp());
679  hA[2]->Fill(mWP.z());
680 
681  hA[24]->Fill(mWP.Eta());
682  hA[25]->Fill(T.hadronicRecoil.Eta());
683  hA[26]->Fill(mWP.Perp(),T.hadronicRecoil.Eta());
684  if(mWP.Eta()>0) {
685  hA[27]->Fill(T.hadronicRecoil.Eta());
686  }
687  else {
688  hA[28]->Fill(T.hadronicRecoil.Eta());
689  }
690 
691  //hadronic recoil and correlations with W from pythia
692  TVector3 hadronicPt(T.hadronicRecoil.X(),T.hadronicRecoil.Y(),0); //transverse momentum vector
693  hA[3]->Fill(mWP.Perp()-hadronicPt.Perp());
694  hA[4]->Fill(mWP.Perp(),hadronicPt.Perp());
695  hA[5]->Fill(hadronicPt.Perp());
696 
697  float delPhi=mWP.DeltaPhi(-hadronicPt);
698  hA[6]->Fill(mWP.Perp(),delPhi);
699  hA[7]->Fill(mWP.Perp(),mWP.Perp()-hadronicPt.Perp());
700  hA[8]->Fill(T.hadronicRecoil.Perp(),delPhi);
701 
702  //electron reco and from geant record
703  hA[9]->Fill(mElectronP.Perp());
704  hA[10]->Fill(T.cluster.ET);
705  hA[11]->Fill(mElectronP.Perp(),T.cluster.ET);
706  hA[12]->Fill(mElectronP.Perp(),mElectronP.Perp()-T.cluster.ET);
707 
708  TVector3 electronPt(T.cluster.position.X(),T.cluster.position.Y(),0); //transvers momentum vector
709  electronPt.SetMag(T.cluster.ET);
710 
711  //neutrino reco and from geant record
712  TVector3 neutrinoPt=-1*(hadronicPt+electronPt);
713  hA[13]->Fill(neutrinoPt.Perp());
714  hA[14]->Fill(mNeutrinoP.Perp());
715  hA[15]->Fill(mNeutrinoP.Perp(),neutrinoPt.Perp());
716  hA[16]->Fill(mNeutrinoP.Perp()-neutrinoPt.Perp());
717 
718  hA[17]->Fill(mNeutrinoP.Perp(),mElectronP.Perp());
719 
720  float recoDeltaPhi=neutrinoPt.DeltaPhi(electronPt);
721  float geantDeltaPhi=mNeutrinoP.DeltaPhi(mElectronP);
722  hA[18]->Fill(geantDeltaPhi,recoDeltaPhi);
723  hA[19]->Fill(cos(geantDeltaPhi)-cos(recoDeltaPhi));
724 
725  float Mt=sqrt(2*T.cluster.ET*neutrinoPt.Perp()*(1-cos(recoDeltaPhi))); //real data
726  float gMt=sqrt(2*mElectronP.Perp()*mNeutrinoP.Perp()*(1-cos(geantDeltaPhi)));
727 
728  hA[20]->Fill(Mt);
729  hA[21]->Fill(gMt);
730 
731  hA[22]->Fill(Mt,T.cluster.ET);
732  hA[23]->Fill(gMt-Mt);
733 
734  //Kinematics
735  //reconstruct W pL from reconstructed quantities
736  float trueWpL=mWP.z();
737  float eleTheta=T.primP.Theta();
738  float ratioE=T.cluster.energy/40.0;
739  float pLRecoPlus=80.0*(ratioE)*((cos(eleTheta))+sqrt(cos(eleTheta)*cos(eleTheta)+sin(eleTheta)*sin(eleTheta)*(1-ratioE*ratioE)))/(ratioE*ratioE*sin(eleTheta)*sin(eleTheta));//+ sqrt solution
740  float pLRecoMinus=80.0*(ratioE)*((cos(eleTheta))-sqrt(cos(eleTheta)*cos(eleTheta)+sin(eleTheta)*sin(eleTheta)*(1-ratioE*ratioE)))/(ratioE*ratioE*sin(eleTheta)*sin(eleTheta));//- sqrt solution
741  hA[29]->Fill(pLRecoPlus);
742  hA[30]->Fill(trueWpL-pLRecoPlus);
743  hA[31]->Fill(trueWpL,pLRecoPlus);
744 
745  hA[32]->Fill(pLRecoMinus);
746  hA[33]->Fill(trueWpL-pLRecoMinus);
747  hA[34]->Fill(trueWpL,pLRecoMinus);
748 
749  const StMuTrack *prTr=T.prMuTrack; assert(prTr);
750  float p_chrg=prTr->charge();
751  if(p_chrg > 0) continue;
752 
753  float eleEta=T.primP.Eta();
754  //sort 2 solutions by electron eta
755  if(eleEta<-0.8) {
756  hA[35]->Fill(trueWpL,pLRecoMinus);
757  hA[37]->Fill(trueWpL,pLRecoPlus);
758  }
759  else if(eleEta>0.8) {
760  hA[36]->Fill(trueWpL,pLRecoMinus);
761  hA[38]->Fill(trueWpL,pLRecoPlus);
762  }
763 
764  if(T.cluster.ET < 30) continue; //only W's we find in data
765  //Correlate W pL with electron E in 3 electron eta ranges
766  if(eleEta < -0.8) {
767  hA[39]->Fill(mWP.z(),T.cluster.energy);
768  hA[42]->Fill(T.cluster.energy);
769  }
770  if(eleEta > 0.8){
771  hA[40]->Fill(mWP.z(),T.cluster.energy);
772  hA[43]->Fill(T.cluster.energy);
773  }
774  if(eleEta > -0.1 && eleEta < 0.1) {
775  hA[41]->Fill(mWP.z(),T.cluster.energy);
776  hA[44]->Fill(T.cluster.energy);
777  }
778 
779  }
780  }
781 }
782 
783 
784 // $Log: St2009pubMcMaker.cxx,v $
785 // Revision 1.7 2011/09/14 14:23:21 stevens4
786 // update used for cross section PRD paper
787 //
788 // Revision 1.6 2010/05/03 19:54:35 stevens4
789 // only try to calc effic if W->e+nu is found in McEvent
790 //
791 // Revision 1.5 2010/03/14 22:50:31 balewski
792 // *** empty log message ***
793 //
794 // Revision 1.4 2010/02/19 21:04:10 stevens4
795 // cleanup unused variables
796 //
797 // Revision 1.3 2010/02/18 19:48:10 stevens4
798 // add more effic histograms and cleanup
799 //
800 // Revision 1.2 2010/01/21 17:54:31 stevens4
801 // add effic histos and charge seperated background plots
802 //
803 // Revision 1.1 2009/11/23 23:00:18 balewski
804 // code moved spin-pool
805 //
806 // Revision 1.1 2009/11/23 21:11:18 balewski
807 // start
808 //
St2009pubMcMaker
maker to retrieve info from geant.root files for comparison with reco quantities from MC ...
Definition: St2009pubMcMaker.h:31
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
St2009pubMcMaker::Make
virtual Int_t Make()
Definition: St2009pubMcMaker.cxx:55
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
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::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
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
kStOK
Definition: Stypes.h:40
StMuTrack::lastPoint
const StThreeVectorF & lastPoint() const
Returns positions of last measured point.
Definition: StMuTrack.h:262
StMuTrack::phi
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
StMuTrack::globalTrack
const StMuTrack * globalTrack() const
Returns pointer to associated global track. Null pointer if no global track available.
Definition: StMuTrack.h:272
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
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Generated by   doxygen 1.8.5