StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2011pubMcMaker.cxx
1 // $Id: St2011pubMcMaker.cxx,v 1.3 2013/09/13 19:33:13 stevens4 Exp $
2 //
3 //*-- Author : Justin Stevens, IUCF
4 //
5 
6 #include "St2011WMaker.h"
7 #include "St2011pubMcMaker.h"
8 #include "StEmcUtil/geometry/StEmcGeom.h"
9 
10 //need these to get MC record
11 #include "tables/St_g2t_tpc_hit_Table.h"
12 #include "StMcEventMaker/StMcEventMaker.h"
13 #include "StMcEvent/StMcEvent.hh"
14 #include "StMcEvent/StMcVertex.hh"
15 #include "StMcEvent/StMcTrack.hh"
16 
17 ClassImp(St2011pubMcMaker)
18 
19 //_______________________________________________________
20 //_______________________________________________________
21 St2011pubMcMaker::St2011pubMcMaker(const char *name):StMaker(name){
22  wMK=0;HList=0;
23 
24 }
25 
26 
27 //_______________________________________________________
28 //_______________________________________________________
29 St2011pubMcMaker::~St2011pubMcMaker(){
30  //
31 }
32 
33 
34 //_______________________________________________________
35 //_______________________________________________________
36 Int_t St2011pubMcMaker::Init(){
37  assert(wMK);
38  assert(HList);
39  initHistos();
40  return StMaker::Init();
41 }
42 
43 
44 //_______________________________________________________
45 //_______________________________________________________
46 Int_t
47 St2011pubMcMaker::Make(){
48  //printf("-----------in %s\n", GetName());
49  //only get geant particle info for W MC
50  if(wMK->isMC==350){
51  if(doMCanalysis()){
52  doWanalysis();
53  doWefficiency();
54  }
55  }
56  return kStOK;
57 }
58 
59 //_______________________________________________________
60 //_______________________________________________________
61 void
62 St2011pubMcMaker::doWanalysis(){
63  //has access to whole W-algo-maker data via pointer 'wMK'
64 
65  // run through W cuts to fill other histos............
66  for(unsigned int iv=0;iv<wMK->wEve->vertex.size();iv++) {
67  WeveVertex &V=wMK->wEve->vertex[iv];
68  for(unsigned int it=0;it<V.eleTrack.size();it++) {
69  WeveEleTrack &T=V.eleTrack[it];
70  if(T.isMatch2Cl==false) continue;
71  assert(T.cluster.nTower>0); // internal logical error
72  assert(T.nearTotET>0); // internal logical error
73 
74  if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
75  if(T.awayTotET> 30.) continue; // too large awayET , Jan
76  //Full W cuts applied at this point
77 
78  //W info from pythia record
79  hA[1]->Fill(mWP.Perp());
80  hA[2]->Fill(mWP.z());
81 
82  hA[24]->Fill(mWP.Eta());
83  hA[25]->Fill(T.hadronicRecoil.Eta());
84  hA[26]->Fill(mWP.Perp(),T.hadronicRecoil.Eta());
85  if(mWP.Eta()>0) {
86  hA[27]->Fill(T.hadronicRecoil.Eta());
87  }
88  else {
89  hA[28]->Fill(T.hadronicRecoil.Eta());
90  }
91 
92  //hadronic recoil and correlations with W from pythia
93  TVector3 hadronicPt(T.hadronicRecoil.X(),T.hadronicRecoil.Y(),0); //transverse momentum vector
94  hA[3]->Fill(mWP.Perp()-hadronicPt.Perp());
95  hA[4]->Fill(mWP.Perp(),hadronicPt.Perp());
96  hA[5]->Fill(hadronicPt.Perp());
97 
98  float delPhi=mWP.DeltaPhi(-hadronicPt);
99  hA[6]->Fill(mWP.Perp(),delPhi);
100  hA[7]->Fill(mWP.Perp(),mWP.Perp()-hadronicPt.Perp());
101  hA[8]->Fill(T.hadronicRecoil.Perp(),delPhi);
102 
103  //electron reco and from geant record
104  hA[9]->Fill(mElectronP.Perp());
105  hA[10]->Fill(T.cluster.ET);
106  hA[11]->Fill(mElectronP.Perp(),T.cluster.ET);
107  hA[12]->Fill(mElectronP.Perp(),mElectronP.Perp()-T.cluster.ET);
108 
109  TVector3 electronPt(T.cluster.position.X(),T.cluster.position.Y(),0); //transvers momentum vector
110  electronPt.SetMag(T.cluster.ET);
111 
112  //neutrino reco and from geant record
113  TVector3 neutrinoPt=-1*(hadronicPt+electronPt);
114  hA[13]->Fill(neutrinoPt.Perp());
115  hA[14]->Fill(mNeutrinoP.Perp());
116  hA[15]->Fill(mNeutrinoP.Perp(),neutrinoPt.Perp());
117  hA[16]->Fill(mNeutrinoP.Perp()-neutrinoPt.Perp());
118 
119  hA[17]->Fill(mNeutrinoP.Perp(),mElectronP.Perp());
120 
121  float recoDeltaPhi=neutrinoPt.DeltaPhi(electronPt);
122  float geantDeltaPhi=mNeutrinoP.DeltaPhi(mElectronP);
123 
124  hA[18]->Fill(geantDeltaPhi,recoDeltaPhi);
125  hA[19]->Fill(cos(geantDeltaPhi)-cos(recoDeltaPhi));
126 
127  float Mt=sqrt(2*T.cluster.ET*neutrinoPt.Perp()*(1-cos(recoDeltaPhi))); //real data
128  float gMt=sqrt(2*mElectronP.Perp()*mNeutrinoP.Perp()*(1-cos(geantDeltaPhi)));
129 
130  hA[20]->Fill(Mt);
131  hA[21]->Fill(gMt);
132 
133  hA[22]->Fill(Mt,T.cluster.ET);
134  hA[23]->Fill(gMt-Mt);
135 
136  //Kinematics
137  //reconstruct W pL from reconstructed quantities
138  float trueWpL=mWP.z();
139  float eleTheta=T.primP.Theta();
140  float ratioE=T.cluster.energy/40.0;
141  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
142  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
143  hA[29]->Fill(pLRecoPlus);
144  hA[30]->Fill(trueWpL-pLRecoPlus);
145  hA[31]->Fill(trueWpL,pLRecoPlus);
146 
147  hA[32]->Fill(pLRecoMinus);
148  hA[33]->Fill(trueWpL-pLRecoMinus);
149  hA[34]->Fill(trueWpL,pLRecoMinus);
150 
151  const StMuTrack *prTr=T.prMuTrack; assert(prTr);
152  float p_chrg=prTr->charge();
153  if(p_chrg > 0) continue;
154 
155  float eleEta=T.primP.Eta();
156  //sort 2 solutions by electron eta
157  if(eleEta<-0.8) {
158  hA[35]->Fill(trueWpL,pLRecoMinus);
159  hA[37]->Fill(trueWpL,pLRecoPlus);
160  }
161  else if(eleEta>0.8) {
162  hA[36]->Fill(trueWpL,pLRecoMinus);
163  hA[38]->Fill(trueWpL,pLRecoPlus);
164  }
165 
166  if(T.cluster.ET < 30) continue; //only W's we find in data
167  //Correlate W pL with electron E in 3 electron eta ranges
168  if(eleEta < -0.8) {
169  hA[39]->Fill(mWP.z(),T.cluster.energy);
170  hA[42]->Fill(T.cluster.energy);
171  }
172  if(eleEta > 0.8){
173  hA[40]->Fill(mWP.z(),T.cluster.energy);
174  hA[43]->Fill(T.cluster.energy);
175  }
176  if(eleEta > -0.1 && eleEta < 0.1) {
177  hA[41]->Fill(mWP.z(),T.cluster.energy);
178  hA[44]->Fill(T.cluster.energy);
179  }
180 
181  }
182  }
183 }
184 
185 
186 //_______________________________________________________
187 //_______________________________________________________
188 void
189 St2011pubMcMaker::doWefficiency(){
190 
191  //only count leptons in our eta range
192  if(fabs(mElectronP.Eta()) > 1.) return;
193  //ele has |eta| < 1
194  hA[50]->Fill(mElectronP.Perp());
195  hA[54]->Fill(mElectronP.Eta());
196  hA[58]->Fill(mVertex.Z());
197  hA[62]->Fill(mElectronP.Phi());
198 
199  hA[68]->Fill(mElectronP.Perp()); //forJoe
200 
201  //plot for Scott
202  TVector3 detEle; //where lepton would hit BEMC
203  float Rcylinder= wMK->mBtowGeom->Radius();
204  detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
205  detEle.SetZ(detEle.Z()+mVertex.Z());
206  hA[66]->Fill(detEle.Eta(),mElectronP.Perp());
207 
208  //trigger efficiency
209  if(!wMK->wEve->l2bitET) return;
210  //good trig
211  hA[51]->Fill(mElectronP.Perp());
212  hA[55]->Fill(mElectronP.Eta());
213  hA[59]->Fill(mVertex.Z());
214  hA[63]->Fill(mElectronP.Phi());
215  hA[69]->Fill(mElectronP.Perp());//forJoe
216 
217  hA[67]->Fill(detEle.Eta(),mElectronP.Perp());
218 
219  //vertex efficiency
220  if(wMK->wEve->vertex.size()<=0) return;
221  //vertex rank>0 and |z|<100
222  hA[52]->Fill(mElectronP.Perp());
223  hA[56]->Fill(mElectronP.Eta());
224  hA[60]->Fill(mVertex.Z());
225  hA[64]->Fill(mElectronP.Phi());
226 
227  hA[70]->Fill(mElectronP.Perp());//forJoe
228 
229  //float diff=wMK->wEve->vertex[0].z - mVertex.Z();
230  //cout<<"diff="<<diff<<endl;
231 
232  //reco efficiency
233  for(unsigned int iv=0;iv<wMK->wEve->vertex.size();iv++) {
234  WeveVertex &V=wMK->wEve->vertex[iv];
235  for(unsigned int it=0;it<V.eleTrack.size();it++) {
236  WeveEleTrack &T=V.eleTrack[it];
237  if(T.isMatch2Cl==false) continue;
238  assert(T.cluster.nTower>0); // internal logical error
239  assert(T.nearTotET>0); // internal logical error
240 
241  if(T.cluster.ET/T.nearTotET< wMK->par_nearTotEtFrac)
242  continue; // too large nearET
243  if(T.ptBalance.Perp() < wMK->par_ptBalance || T.awayTotET > 30.) //Jan
244  continue;
245 
246  //pass all W cuts
247  hA[53]->Fill(mElectronP.Perp());
248  hA[57]->Fill(mElectronP.Eta());
249  hA[61]->Fill(mVertex.Z());
250  hA[65]->Fill(mElectronP.Phi());
251 
252  hA[71]->Fill(mElectronP.Perp());//forJoe
253  }
254  }
255 
256 }
257 
258 
259 //________________________________________________
260 //________________________________________________
261 bool
262 St2011pubMcMaker::doMCanalysis(){
263  StMcEvent* mMcEvent = 0;
264  mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
265  assert(mMcEvent);
266 
267  //initialize momentum vectors
268  StThreeVectorF pW; float eW;
269  StThreeVectorF pNeutrino; //float eNeutrino;
270  StThreeVectorF pElectron; //float eElectron;
271 
272  StMcVertex *V=mMcEvent->primaryVertex();
273  mVertex=TVector3(V->position().x(),V->position().y(),V->position().z());
274 
275  unsigned int i=1;
276  int found=0;
277  while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
278  StMcTrack* mcTrack = mMcEvent->tracks()[i];
279  int pdgId=mcTrack->pdgId();
280  //float pt=mcTrack->pt();
281  //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
282  if(pdgId==11 || pdgId==-11){ //select e+ and e-
283  if(abs(mcTrack->parent()->pdgId()) == 24 ){
284  pElectron=mcTrack->momentum();
285  //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
286  pW=mcTrack->parent()->momentum();
287  eW=mcTrack->parent()->energy();
288  //LOG_INFO<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
289  found++;
290  }
291  }
292  if(pdgId==12 || pdgId==-12){ //select neutrino
293  if(abs(mcTrack->parent()->pdgId()) == 24 ){
294  pNeutrino=mcTrack->momentum();
295  //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
296  pW=mcTrack->parent()->momentum();
297  eW=mcTrack->parent()->energy();
298  //LOG_INFO<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
299  found++;
300  }
301  }
302  i++;
303  }
304 
305  if(found!=2) return false;
306 
307  mWP=TVector3(pW.x(),pW.y(),pW.z());
308  mNeutrinoP=TVector3(pNeutrino.x(),pNeutrino.y(),pNeutrino.z());
309  mElectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
310  TVector3 diff=mWP-mNeutrinoP-mElectronP;
311  if(diff.Mag() > 0.0001) //should get exactly right
312  LOG_INFO<<"\n \n W+e+nu vector sum ="<<diff.Mag()<<endm;
313  if(mElectronP.Mag() < 0.0001)
314  LOG_INFO<<"\n \n no lepton track ="<<endm;
315 
316  //calculate x1 and x2 from W rapidity
317  float rapW = 0.5*log((eW+mWP.Z())/(eW-mWP.Z()));
318  float mw2sqs=80.4/500.;
319  float x1=mw2sqs*exp(rapW);
320  float x2=mw2sqs*exp(-rapW);
321 
322  hA[72]->Fill(rapW);
323  if(fabs(mElectronP.Eta())<1){ //require midrapidity e
324  hA[73]->Fill(x1);
325  hA[74]->Fill(x2);
326  hA[75]->Fill(x1,x2);
327  hA[76]->Fill(x1-x2);
328  }
329 
330  return true;
331 
332 }
333 
334 // $Log: St2011pubMcMaker.cxx,v $
335 // Revision 1.3 2013/09/13 19:33:13 stevens4
336 // Updates to code for combined 2011+2012 result presented to spin PWG 9.12.13
337 //
338 // Revision 1.2 2013/07/03 16:53:07 stevens4
339 // Update for efficiency studies with embedding
340 //
341 // Revision 1.1 2011/02/10 20:33:25 balewski
342 // start
343 //
St2011pubMcMaker::Make
virtual Int_t Make()
Definition: St2011pubMcMaker.cxx:47
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
St2011pubMcMaker
maker to retrieve info from geant.root files for comparison with reco quantities from MC ...
Definition: St2011pubMcMaker.h:30
kStOK
Definition: Stypes.h:40
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