StRoot: StSpinPool/StWalgoB2009/St2009pubMcMaker.cxx Source File

00001 // $Id: St2009pubMcMaker.cxx,v 1.7 2011/09/14 14:23:21 stevens4 Exp $
00002 //
00003 //*-- Author : Justin Stevens, IUCF
00004 // 
00005 
00006 #include "St2009WMaker.h"
00007 #include "St2009ZMaker.h"
00008 #include "St2009pubMcMaker.h"
00009 #include "StEmcUtil/geometry/StEmcGeom.h"
00010 #include <TMath.h>
00011 
00012 //need these to get MC record
00013 #include "tables/St_g2t_tpc_hit_Table.h"
00014 #include "StMcEventMaker/StMcEventMaker.h"
00015 #include "StMcEvent/StMcEvent.hh"
00016 #include "StMcEvent/StMcVertex.hh"
00017 #include "StMcEvent/StMcTrack.hh"
00018 
00019 //muDst needed for zdcRate
00020 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00021 #include <StMuDSTMaker/COMMON/StMuDst.h>
00022 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00023 #include <StMuDSTMaker/COMMON/StMuMcTrack.h>
00024 
00025 ClassImp(St2009pubMcMaker)
00026 
00027 //_______________________________________________________
00028 //_______________________________________________________
00029 St2009pubMcMaker::St2009pubMcMaker(const char *name):StMaker(name){
00030   wMK=0;HList=0;
00031     
00032 }
00033 
00034 
00035 //_______________________________________________________
00036 //_______________________________________________________
00037 St2009pubMcMaker::~St2009pubMcMaker(){
00038   //
00039 }
00040 
00041 
00042 //_______________________________________________________
00043 //_______________________________________________________
00044 Int_t St2009pubMcMaker::Init(){
00045   assert(wMK);
00046   assert(HList);
00047   initHistos();
00048   return StMaker::Init();
00049 }
00050 
00051 
00052 //_______________________________________________________
00053 //_______________________________________________________
00054 Int_t 
00055 St2009pubMcMaker::Make(){
00056   //printf("-----------in %s\n", GetName());
00057   //only get geant particle info for W MC
00058   if(wMK->isMC==350){
00059     if(doWMCanalysis()){
00060       doWefficiency();
00061       doWanalysis();
00062     }
00063   }
00064   if(wMK->isMC==352){
00065     if(doZMCanalysis()){
00066       doZefficiency();
00067     }
00068   }
00069   return kStOK;
00070 }
00071 
00072      
00073 //_______________________________________________________
00074 //_______________________________________________________
00075 void 
00076 St2009pubMcMaker::doWefficiency(){
00077 
00078   //initialize some kinematic variables for later use 
00079   
00080   //wrap barrel onto 2 modules to see boundaries
00081   const float PI=TMath::Pi();
00082   float elePhi=mElectronP.Phi();
00083   if(elePhi<0) elePhi+=2*PI;
00084   elePhi=elePhi*180/PI;
00085   elePhi+=3.0;
00086   float modulePair=12.0;
00087   int resid = (int) (elePhi/modulePair);
00088   float phiMod= elePhi - resid*modulePair;
00089   // 'detector' eta
00090   TVector3 detEle; //where lepton would hit BEMC
00091   float Rcylinder= wMK->mBtowGeom->Radius();
00092   detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
00093   detEle.SetZ(detEle.Z()+mVertex.Z());
00094 
00095   float w=1; //determine zVertex weighting for event
00096   if(wMK->isMC ==350) w=wMK->hReweight->GetBinContent(wMK->hReweight->FindBin(mVertex.Z()));
00097   float zdcRate=wMK->mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
00098     
00099   //only count leptons in our eta range
00100   if(fabs(mElectronP.Eta()) > 1.0) return;
00101   //only count leptons in our ET range
00102   if(fabs(mElectronP.Perp()) < 25.) return;
00103 
00104   int bx7=wMK->wEve.bx7;
00105   if( (bx7>30 && bx7<40) || (bx7>110) ) return; //skip abort gaps
00106 
00107   //ele has |eta| < 1 and ET > 25
00108   hA[50]->Fill(mElectronP.Perp(),w);
00109   hA[54]->Fill(mElectronP.Eta(),w);
00110   hA[58]->Fill(mVertex.Z(),w);
00111   hA[62]->Fill(mElectronP.Phi(),w);
00112   hB[83]->Fill(mElectronP.Eta(),mElectronP.Perp(),w);
00113   hA[89]->Fill(zdcRate,w);
00114 
00115   hA[68]->Fill(mElectronP.Perp(),w); //forJoe
00116   hA[77]->Fill(phiMod,w); //wrap barrel onto 2 modules
00117   hB[66]->Fill(detEle.Eta(),mElectronP.Perp(),w);
00118   hA[79]->Fill(detEle.Eta(),w);
00119   hA[122]->Fill(zdcRate,wMK->wEve.bx7);
00120 
00121   //check reconstructed vs thrown ET
00122   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00123     WeveVertex &V=wMK->wEve.vertex[iv];
00124     for(uint it=0;it<V.eleTrack.size();it++) {
00125       WeveEleTrack &T=V.eleTrack[it];
00126       if(T.isMatch2Cl==false) continue;
00127       assert(T.cluster.nTower>0); // internal logical error
00128       assert(T.nearTotET>0); // internal logical error
00129       
00130       hB[68]->Fill(mElectronP.Perp(),T.cluster.ET,w);
00131     }
00132   }
00133 
00134   //trigger efficiency
00135   if(!wMK->wEve.l2bitET) return;
00136   //good trig
00137   hA[51]->Fill(mElectronP.Perp(),w);
00138   hA[55]->Fill(mElectronP.Eta(),w);
00139   hA[59]->Fill(mVertex.Z(),w);
00140   hA[63]->Fill(mElectronP.Phi(),w);
00141   hB[84]->Fill(mElectronP.Eta(),mElectronP.Perp(),w);
00142   hA[90]->Fill(zdcRate,w);
00143 
00144   hA[69]->Fill(mElectronP.Perp(),w);
00145   hA[78]->Fill(phiMod,w);//wrap barrel onto 2 modules
00146   hB[67]->Fill(detEle.Eta(),mElectronP.Perp(),w);
00147   hA[80]->Fill(detEle.Eta(),w);
00148 
00149   //vertex efficiency
00150   //require one vertex w/ rank>0 and |z|<100
00151   if(wMK->wEve.vertex.size()<=0) return; 
00152   //require geant vert match reco vert (w/in 5 cm)
00153   bool matchVert=false;
00154   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00155     WeveVertex &V=wMK->wEve.vertex[iv];
00156     //cout<<endl<<"vertex diff "<<fabs(V.z - mVertex.Z())<<endl;
00157     if(fabs(V.z - mVertex.Z()) < 5.0) matchVert=true;
00158   }
00159   if(!matchVert) return;
00160   //cout<<"after matched vert"<<endl;
00161   hA[52]->Fill(mElectronP.Perp(),w); 
00162   hA[56]->Fill(mElectronP.Eta(),w);
00163   hA[60]->Fill(mVertex.Z(),w);
00164   hA[64]->Fill(mElectronP.Phi(),w);
00165   hA[91]->Fill(zdcRate,w);
00166 
00167   hA[70]->Fill(mElectronP.Perp(),w);//forJoe
00168 
00169   //get MC track info
00170   int eleTrId=999; TVector3 eleTrP;
00171   int size=wMK->mMuDstMaker->muDst()->mcArray(1)->GetEntries();
00172   cout<<"size="<<size<<endl;
00173   for(int i=0; i<size; i++){
00174     StMuMcTrack* muMcTrack=(StMuMcTrack*)wMK->mMuDstMaker->muDst()->mcArray(1)->UncheckedAt(i);
00175     int geantId=muMcTrack->GePid();
00176     int id=muMcTrack->Id();
00177     if(geantId==2 || geantId==3){
00178       cout<<i<<" id="<<id<<" geantId="<<geantId<<endl;
00179       eleTrId=id;
00180       eleTrP=TVector3(muMcTrack->Pxyz().x(),muMcTrack->Pxyz().y(),muMcTrack->Pxyz().z());
00181       cout<<"phiDiff="<<eleTrP.Phi()-mElectronP.Phi()<<" etaDiff="<<eleTrP.Eta()-mElectronP.Eta()<<endl;
00182       break;
00183     }
00184   }
00185 
00186   int nTrackAssoc=0; float recoPtAssoc=0;
00187   int flagAssoc=0;
00188   int nHitAssoc=0; float hitFracAssoc=0.;
00189   float rInAssoc=0.; float rOutAssoc=0.;
00190   //check tracking efficiency before any track cuts
00191   for(uint iv=0;iv<wMK->wEve.vertex.size(); iv++) {
00192     uint vertID=wMK->wEve.vertex[iv].id;
00193     wMK->mMuDstMaker->muDst()->setVertexIndex(vertID);
00194     Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
00195     for(int itr=0;itr<nPrimTrAll;itr++) {
00196       StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
00197       if(prTr->flag()<=0) continue;
00198       const StMuTrack *glTr=prTr->globalTrack();
00199       if(glTr==0) continue; // see the reason in _accessMuDst
00200       if(prTr->flag()!=301) continue;// TPC+prim vertex tracks
00201       if(prTr->pt()<10.) continue;
00202       
00203       StThreeVectorF ro=glTr->lastPoint();
00204       int secID=WtpcFilter::getTpcSec(ro.phi(),ro.pseudoRapidity());
00205       if(secID==20) continue;
00206       
00207       //only use with private copy of StMuTrack which makes these functions public (currently in dev they're private)
00208       //int idTruth=prTr->idTruth();
00209       //int qaTruth=prTr->qaTruth();
00210       
00211       int idTruth=0;
00212       int qaTruth=0;
00213 
00214       float delPhi=prTr->phi()-mElectronP.Phi();
00215       float delEta=prTr->eta()-mElectronP.Eta();
00216       if( eleTrId==idTruth ) { //reco track matches to electron
00217         
00218         flagAssoc=prTr->flag();
00219         nHitAssoc=prTr->nHitsFit();
00220         hitFracAssoc=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00221         rInAssoc=glTr->firstPoint().perp();
00222         rOutAssoc=glTr->lastPoint().perp();
00223         
00224         cout<<"found electron, delR="<<sqrt(delPhi*delPhi+delEta*delEta)<<" nHits="<<prTr->nHitsFit()<<" fitFrac="<<hitFracAssoc<<endl;
00225         cout<<"idTruth="<<idTruth<<" qaTruth="<<qaTruth<<endl;          
00226         recoPtAssoc=prTr->pt();
00227         nTrackAssoc++;
00228       }
00229     }
00230   }
00231   if(nTrackAssoc==1) {
00232     hA[120]->Fill(recoPtAssoc,recoPtAssoc-mElectronP.Perp());
00233     hA[124]->Fill(recoPtAssoc,w);
00234     hA[126]->Fill(1./recoPtAssoc,w);
00235     hA[128]->Fill(recoPtAssoc,nHitAssoc);
00236     hA[129]->Fill(recoPtAssoc,hitFracAssoc);
00237     hA[130]->Fill(recoPtAssoc,rInAssoc);
00238     hA[131]->Fill(recoPtAssoc,rOutAssoc);
00239     hA[132]->Fill(1./recoPtAssoc,nHitAssoc);
00240     hA[133]->Fill(1./recoPtAssoc,hitFracAssoc);
00241     hA[134]->Fill(1./recoPtAssoc,rInAssoc);
00242     hA[135]->Fill(1./recoPtAssoc,rOutAssoc);
00243     hA[136]->Fill(flagAssoc);
00244   }
00245   cout<<"number of associated tracks = "<<nTrackAssoc<<endl;
00246 
00247   // track efficiency
00248   bool foundTrack=false; int nTrk=0; float recoPt=0; 
00249   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00250     WeveVertex &V=wMK->wEve.vertex[iv];
00251     if(V.eleTrack.size()>0){ //track with pt > 10
00252       hA[85]->Fill(mElectronP.Perp(),w);
00253       hA[86]->Fill(mElectronP.Eta(),w);
00254       hA[87]->Fill(mVertex.Z(),w);
00255       hA[88]->Fill(mElectronP.Phi(),w);
00256       hA[92]->Fill(zdcRate,w);
00257       foundTrack=true;
00258       
00259       //find track 
00260       for(uint it=0;it<V.eleTrack.size();it++) {
00261         WeveEleTrack &T=V.eleTrack[it];
00262         int idTruth=0; //T.prMuTrack->idTruth(); //again need private StMuTrack
00263         if( eleTrId==idTruth ) { //reco track matches to electron
00264           recoPt=T.primP.Perp();
00265           nTrk++;
00266         }
00267       }
00268       break; //only count event once
00269     }
00270   }
00271   if(!foundTrack) return;
00272 
00273   if(nTrackAssoc==1) { //accepted tracks
00274     hA[125]->Fill(recoPtAssoc,w);
00275     hA[127]->Fill(1./recoPtAssoc,w);
00276   }
00277 
00278   if(wMK->wEve.wTag) hA[140]->Fill(mElectronP.Perp(),w);
00279 
00280   //reco efficiency
00281   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00282     WeveVertex &V=wMK->wEve.vertex[iv];
00283     for(uint it=0;it<V.eleTrack.size();it++) {
00284       WeveEleTrack &T=V.eleTrack[it];
00285       if(T.isMatch2Cl==false) continue;
00286       assert(T.cluster.nTower>0); // internal logical error
00287       assert(T.nearTotET>0); // internal logical error
00288       
00289       if(wMK->wEve.zTag)
00290         continue; //event tagged as Z
00291       if(fabs(T.primP.Eta()) > wMK->par_leptonEta) 
00292         continue; // reconstructed eta too large
00293 
00294       if(T.cluster.ET/T.nearTotET < wMK->par_nearTotEtFrac) 
00295         continue; // too large nearET
00296       if(T.sPtBalance < wMK->par_ptBalance) 
00297         continue; // too large signed pt balance
00298       
00299       if(T.cluster.ET < 25.) 
00300         continue; // too small reco ET (use same cut as data)
00301 
00302       //pass all W cuts 
00303       hA[53]->Fill(mElectronP.Perp(),w);
00304       hA[57]->Fill(mElectronP.Eta(),w);
00305       hA[61]->Fill(mVertex.Z(),w);
00306       hA[65]->Fill(mElectronP.Phi(),w);
00307       hA[93]->Fill(zdcRate,w);
00308 
00309       hA[71]->Fill(mElectronP.Perp(),w);//forJoe
00310       hA[121]->Fill(T.primP.Phi()-mElectronP.Phi(),T.primP.Eta()-mElectronP.Eta());
00311       hA[123]->Fill(zdcRate,wMK->wEve.bx7);
00312 
00313       //fill thrown ele ET after all W cuts
00314       hA[111]->Fill(detEle.Eta(),mElectronP.Perp());
00315       hA[113]->Fill(detEle.Eta(),mElectronSmearP.Perp());
00316       hA[115]->Fill(T.cluster.ET/mElectronP.Perp());
00317       hA[116]->Fill(T.cluster.energy/mElectronP.Mag());
00318 
00319       //template function stuff
00320       for(int ires=0; ires<10; ires++){
00321         for(int iscale=0; iscale<101; iscale++){
00322           float scaleFact = iscale*0.004 + 0.8;
00323           hB[120+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
00324           hB[150+ires]->Fill(scaleFact,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
00325           if(fabs(detEle.Eta()) < 0.5) hB[130+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
00326           else if(fabs(detEle.Eta()) < 1.0) hB[140+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact,wRB);
00327         }
00328       }
00329       
00330       break; //only count event once
00331     }
00332   }
00333 
00334 }
00335 
00336 
00337 //_______________________________________________________
00338 //_______________________________________________________
00339 void 
00340 St2009pubMcMaker::doZefficiency(){
00341 
00342   //Acceptance cuts applied
00343   //only count leptons in our eta range
00344   if(fabs(mZelectronP.Eta()) > 1.) return;
00345   if(fabs(mZpositronP.Eta()) > 1.) return;
00346   //only count leptons in our ET range
00347   if(mZelectronP.Perp() < 15.) return;
00348   if(mZpositronP.Perp() < 15.) return;
00349 
00350   TVector3 sumP=mZelectronP+mZpositronP;
00351   float sumE=mZelectronP.Mag()+mZpositronP.Mag();
00352   float geantZmass=sqrt(sumE*sumE-sumP.Dot(sumP));
00353   cout<<"geant mass"<<geantZmass<<endl;
00354   float zdcRate=wMK->mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
00355 
00356   float w=1; //determine zVertex weighting for event
00357   if(wMK->isMC==352) w=wMK->hReweight->GetBinContent(wMK->hReweight->FindBin(mVertex.Z()));
00358 
00359   //fill all events in acceptance
00360   hA[100]->Fill(geantZmass,w);
00361   if(geantZmass>70. && geantZmass<110.) hA[105]->Fill(zdcRate,w);
00362   
00363   //trigger efficiency
00364   if(!wMK->wEve.l2bitET) return;
00365   //fill good trig events
00366   hA[101]->Fill(geantZmass,w);
00367   if(geantZmass>70. && geantZmass<110.) hA[106]->Fill(zdcRate,w);
00368   
00369   //vertex efficiency
00370   if(wMK->wEve.vertex.size()<=0) return;
00371   //fill vertex rank>0 and |z|<100 events
00372   bool matchVert=false;
00373   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00374     WeveVertex &V=wMK->wEve.vertex[iv];
00375     //cout<<endl<<"vertex diff "<<fabs(V.z - mVertex.Z())<<endl;
00376     if(fabs(V.z - mZvertex.Z()) < 5.0) matchVert=true;
00377   }
00378   if(!matchVert) return;
00379   hA[102]->Fill(geantZmass,w);
00380   if(geantZmass>70. && geantZmass<110.) hA[107]->Fill(zdcRate,w);
00381   
00382   //track efficiency
00383   bool gt1track=false;
00384   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00385     WeveVertex &V=wMK->wEve.vertex[iv];
00386     if(V.eleTrack.size()>1){ //fill 2 track with pt > 10 events
00387       hA[103]->Fill(geantZmass,w);
00388       if(geantZmass>70. && geantZmass<110.) hA[108]->Fill(zdcRate,w);
00389       
00390       gt1track=true;
00391       break; //only count event once
00392     }
00393   }
00394   if(!gt1track) return;
00395 
00396   //reco efficiency
00397   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00398     WeveVertex &V=wMK->wEve.vertex[iv];
00399     if(V.eleTrack.size()<2) continue;
00400     for(uint it=0;it<V.eleTrack.size()-1;it++) {//.....select first track:
00401 
00402       WeveEleTrack &T1=V.eleTrack[it];
00403       if(T1.isMatch2Cl==false) continue;
00404       assert(T1.cluster.nTower>0); // internal logical error
00405       assert(T1.nearTotET>0); // internal logical error
00406 
00407       if(T1.cluster.ET<zMK->par_clusterEtZ) continue;
00408       float fracET1=T1.cluster.ET /T1.nearTotET;
00409       if(fracET1<zMK->par_nearTotEtFracZ) continue;
00410       
00411       for (uint it2=it+1;it2<V.eleTrack.size();it2++) { //.....select second track:
00412 
00413         WeveEleTrack &T2=V.eleTrack[it2];
00414         if(T2.isMatch2Cl==false) continue;
00415         assert(T2.cluster.nTower>0); // internal logical error
00416         assert(T2.nearTotET>0); // internal logical error
00417 
00418         if(T2.cluster.ET<zMK->par_clusterEtZ) continue;
00419         float fracET2=T2.cluster.ET /T2.nearTotET;
00420         if(fracET2<zMK->par_nearTotEtFracZ) continue;
00421         
00422         float e1=T1.cluster.energy;
00423         float e2=T2.cluster.energy;
00424         TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
00425         TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
00426 
00427         float del_phi=p1.DeltaPhi(p2);
00428         if(fabs(del_phi)<zMK->par_delPhi12) continue;
00429         TVector3 psum=p1+p2;
00430         float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
00431         if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
00432         //float mass=sqrt(mass2); 
00433         
00434         int Q1Q2=T1.prMuTrack->charge()*T2.prMuTrack->charge();
00435         if (Q1Q2==1) continue;
00436         
00437         //fill reco Z histos
00438         hA[104]->Fill(geantZmass,w); 
00439         if(geantZmass>70. && geantZmass<110.) hA[109]->Fill(zdcRate,w);
00440         
00441       } //second track
00442     } //first track
00443   } //vertex
00444 
00445   return;
00446 }
00447 
00448 
00449 //________________________________________________
00450 //________________________________________________
00451 bool
00452 St2009pubMcMaker::doWMCanalysis(){
00453   StMcEvent* mMcEvent = 0;  
00454   mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
00455   //assert(mMcEvent);
00456 
00457   if(!mMcEvent) return false;
00458 
00459   //initialize momentum vectors
00460   StThreeVectorF pW;        float eW;
00461   StThreeVectorF pNeutrino; //float eNeutrino;
00462   StThreeVectorF pElectron; //float eElectron;
00463 
00464   StMcVertex *V=mMcEvent->primaryVertex(); 
00465   mVertex=TVector3(V->position().x(),V->position().y(),V->position().z());
00466   
00467   //find geant tracks for W and decay daughters
00468   uint i=1;
00469   int found=0;
00470   while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
00471     StMcTrack* mcTrack = mMcEvent->tracks()[i];
00472     int pdgId=mcTrack->pdgId();
00473     float pt=mcTrack->pt();
00474     LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00475     if(pdgId==11 || pdgId==-11){ //select e+ and e-
00476       if(abs(mcTrack->parent()->pdgId()) == 24 ){ 
00477         pElectron=mcTrack->momentum();
00478         //mKeyElectron=mcTrack->key();
00479         //mEveGenElectron=mcTrack->eventGenLabel();
00480         //mIdElectron=i;
00481         LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00482         pW=mcTrack->parent()->momentum();
00483         eW=mcTrack->parent()->energy();
00484         LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00485         found++;
00486       }
00487     }
00488     if(pdgId==12 || pdgId==-12){ //select neutrino
00489       if(abs(mcTrack->parent()->pdgId()) == 24 ){ 
00490         pNeutrino=mcTrack->momentum();
00491         LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00492         pW=mcTrack->parent()->momentum();
00493         eW=mcTrack->parent()->energy();
00494         LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00495         found++;
00496       }
00497     }
00498     i++;
00499   }
00500 
00501   i=1;
00502   while(i<mMcEvent->tracks().size()){
00503     StMcTrack* mcTrack = mMcEvent->tracks()[i];
00504     int pdgId=mcTrack->pdgId();
00505     float pt=mcTrack->pt();
00506     i++;
00507     if(abs(pdgId)==11 || abs(pdgId)==12 || abs(pdgId)==24 || abs(pdgId)==21 || abs(pdgId) < 10 || abs(pdgId)==92 || pt<10.0) continue;
00508     cout<<"high pt MC track: pdgId="<<pdgId<<" pt="<<pt<<endl;
00509     
00510   }
00511  
00512   if(found!=2) return false;
00513   
00514   mWP=TVector3(pW.x(),pW.y(),pW.z());
00515   mNeutrinoP=TVector3(pNeutrino.x(),pNeutrino.y(),pNeutrino.z());
00516   mElectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
00517   TVector3 diff=mWP-mNeutrinoP-mElectronP;
00518   if(diff.Mag() > 0.0001) //should get exactly right
00519     LOG_INFO<<"\n \n W+e+nu vector sum ="<<diff.Mag()<<endm;
00520   if(mElectronP.Mag() < 0.0001)
00521     LOG_INFO<<"\n \n no lepton track ="<<endm;
00522 
00523   //calculate x1 and x2 from W rapidity
00524   float rapW = 0.5*log((eW+mWP.Z())/(eW-mWP.Z()));
00525   float mw2sqs=80.4/500.;
00526   float x1=mw2sqs*exp(rapW);
00527   float x2=mw2sqs*exp(-rapW);
00528 
00529   hA[72]->Fill(rapW);
00530   if(fabs(mElectronP.Eta())<1){ //require midrapidity e
00531     hA[73]->Fill(x1);
00532     hA[74]->Fill(x2);
00533     hA[75]->Fill(x1,x2);
00534     hA[76]->Fill(x1-x2);
00535   }
00536 
00537   //original electron Pt spectrum
00538   if(fabs(mElectronP.Eta()) < 1.)
00539     hA[117]->Fill(mElectronP.Pt());
00540 
00541   //define weighting for rhicbos
00542   wRB = 0.0; //init weight at 0
00543   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};
00544 
00545   float et = mElectronP.Pt();
00546   for(int j=0; j<120; j++){
00547     if(et > j*0.5 && et < (j+1)*0.5 && fabs(mElectronP.Eta()) < 1.){
00548       hA[118]->Fill(et,weight[j]);
00549       wRB = weight[j];
00550       break;
00551     }
00552   }
00553 
00554   wRB=1.0; // drop weighting for now
00555 
00556   //smear electron energy by BTOW resolution
00557   float initE=mElectronP.Mag();
00558   float width=sqrt(0.14*0.14/initE + 0.015*0.015); //resolution factor from NIM
00559   float smear = gRandom->Gaus(0, width);
00560   float finalE=mElectronP.Mag()+smear;
00561   cout<<"initE="<<initE<<" width="<<width<<" smear="<<smear<<" finalE="<<finalE<<endl;
00562   mElectronSmearP=mElectronP; mElectronSmearP.SetMag(finalE);
00563 
00564   //where lepton would hit BEMC
00565   TVector3 detEle;
00566   float Rcylinder= wMK->mBtowGeom->Radius();
00567   detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
00568   detEle.SetZ(detEle.Z()+mVertex.Z());
00569 
00570   //thrown ET for jacob peak comp
00571   hA[110]->Fill(detEle.Eta(),mElectronP.Perp());
00572   //smeared ET
00573   hA[112]->Fill(detEle.Eta(),mElectronSmearP.Perp());
00574 
00575   //template function stuff
00576   for(int ires=0; ires<10; ires++){
00577     float widthTemp = ires/100.;
00578     float smearTemp = gRandom->Gaus(0, widthTemp);
00579     float finalETemp = mElectronP.Mag()+smearTemp*mElectronP.Mag();
00580     //cout<<"initE="<<initE<<" width="<<width<<" smear="<<smear<<" finalE="<<finalE<<endl;
00581     mElectronSmearTempP[ires]=mElectronP; mElectronSmearTempP[ires].SetMag(finalETemp);
00582     for(int iscale=0; iscale<101; iscale++){
00583       float scaleFact = iscale*0.004 +0.9;
00584       hB[110+ires]->Fill(iscale,mElectronSmearTempP[ires].Perp()*scaleFact);
00585     }
00586   }
00587 
00588   return true;
00589 
00590 }
00591 
00592 
00593 
00594 //________________________________________________
00595 //________________________________________________
00596 bool
00597 St2009pubMcMaker::doZMCanalysis(){
00598   StMcEvent* mMcEvent = 0;  
00599   mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
00600   //assert(mMcEvent);
00601 
00602   if(!mMcEvent) return false;
00603 
00604   //initialize momentum vectors
00605   StThreeVectorF pZ;        float eZ;
00606   StThreeVectorF pPositron; //float eNeutrino;
00607   StThreeVectorF pElectron; //float eElectron;
00608 
00609   StMcVertex *V=mMcEvent->primaryVertex(); 
00610   mZvertex=TVector3(V->position().x(),V->position().y(),V->position().z());
00611   
00612   //find geant tracks for Z and e+,e-
00613   uint i=1;
00614   int found=0;
00615   while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
00616     StMcTrack* mcTrack = mMcEvent->tracks()[i];
00617     int pdgId=mcTrack->pdgId();
00618     float pt=mcTrack->pt();
00619     LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00620     if(pdgId==11){ //select e-
00621       if(abs(mcTrack->parent()->pdgId()) == 23 ){ 
00622         pElectron=mcTrack->momentum();
00623         LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00624         pZ=mcTrack->parent()->momentum();
00625         eZ=mcTrack->parent()->energy();
00626         LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00627         found++;
00628       }
00629     }
00630     if(pdgId==-11){ //select e+
00631       if(abs(mcTrack->parent()->pdgId()) == 23 ){ 
00632         pPositron=mcTrack->momentum();
00633         LOG_DEBUG<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00634         pZ=mcTrack->parent()->momentum();
00635         eZ=mcTrack->parent()->energy();
00636         LOG_DEBUG<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00637         found++;
00638       }
00639     }
00640     i++;
00641   }
00642  
00643   if(found!=2) return false; 
00644   
00645   mZP=TVector3(pZ.x(),pZ.y(),pZ.z());
00646   mZpositronP=TVector3(pPositron.x(),pPositron.y(),pPositron.z());
00647   mZelectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
00648   TVector3 diff=mZP-mZpositronP-mZelectronP;
00649   if(diff.Mag() > 0.0001) //should get exactly right
00650     LOG_INFO<<"\n \n Z+e+e vector sum ="<<diff.Mag()<<endm;
00651   if(mZelectronP.Mag() < 0.0001)
00652     LOG_INFO<<"\n \n no lepton track ="<<endm;
00653 
00654   return true;
00655 }
00656 
00657 
00658 //_______________________________________________________
00659 //_______________________________________________________
00660 void 
00661 St2009pubMcMaker::doWanalysis(){
00662   //has access to whole W-algo-maker data via pointer 'wMK'
00663  
00664   // run through W cuts to fill other histos............
00665   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00666     WeveVertex &V=wMK->wEve.vertex[iv];
00667     for(uint it=0;it<V.eleTrack.size();it++) {
00668       WeveEleTrack &T=V.eleTrack[it];
00669       if(T.isMatch2Cl==false) continue;
00670       assert(T.cluster.nTower>0); // internal logical error
00671       assert(T.nearTotET>0); // internal logical error
00672       
00673       if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
00674       if(T.awayTotET> 30.) continue; // too large awayET , Jan
00675       //Full W cuts applied at this point
00676 
00677       //W info from pythia record
00678       hA[1]->Fill(mWP.Perp());
00679       hA[2]->Fill(mWP.z());
00680       
00681       hA[24]->Fill(mWP.Eta());
00682       hA[25]->Fill(T.hadronicRecoil.Eta());
00683       hA[26]->Fill(mWP.Perp(),T.hadronicRecoil.Eta());
00684       if(mWP.Eta()>0) {
00685         hA[27]->Fill(T.hadronicRecoil.Eta());
00686       }
00687       else {
00688         hA[28]->Fill(T.hadronicRecoil.Eta());
00689       }
00690       
00691       //hadronic recoil and correlations with W from pythia
00692       TVector3 hadronicPt(T.hadronicRecoil.X(),T.hadronicRecoil.Y(),0); //transverse momentum vector
00693       hA[3]->Fill(mWP.Perp()-hadronicPt.Perp());
00694       hA[4]->Fill(mWP.Perp(),hadronicPt.Perp());
00695       hA[5]->Fill(hadronicPt.Perp()); 
00696 
00697       float delPhi=mWP.DeltaPhi(-hadronicPt);
00698       hA[6]->Fill(mWP.Perp(),delPhi);
00699       hA[7]->Fill(mWP.Perp(),mWP.Perp()-hadronicPt.Perp());
00700       hA[8]->Fill(T.hadronicRecoil.Perp(),delPhi);
00701       
00702       //electron reco and from geant record
00703       hA[9]->Fill(mElectronP.Perp());
00704       hA[10]->Fill(T.cluster.ET);      
00705       hA[11]->Fill(mElectronP.Perp(),T.cluster.ET);
00706       hA[12]->Fill(mElectronP.Perp(),mElectronP.Perp()-T.cluster.ET);
00707 
00708       TVector3 electronPt(T.cluster.position.X(),T.cluster.position.Y(),0); //transvers momentum vector
00709       electronPt.SetMag(T.cluster.ET);
00710       
00711       //neutrino reco and from geant record
00712       TVector3 neutrinoPt=-1*(hadronicPt+electronPt);
00713       hA[13]->Fill(neutrinoPt.Perp());  
00714       hA[14]->Fill(mNeutrinoP.Perp());
00715       hA[15]->Fill(mNeutrinoP.Perp(),neutrinoPt.Perp());
00716       hA[16]->Fill(mNeutrinoP.Perp()-neutrinoPt.Perp());
00717 
00718       hA[17]->Fill(mNeutrinoP.Perp(),mElectronP.Perp());
00719 
00720       float recoDeltaPhi=neutrinoPt.DeltaPhi(electronPt);
00721       float geantDeltaPhi=mNeutrinoP.DeltaPhi(mElectronP);
00722       hA[18]->Fill(geantDeltaPhi,recoDeltaPhi);
00723       hA[19]->Fill(cos(geantDeltaPhi)-cos(recoDeltaPhi));
00724 
00725       float Mt=sqrt(2*T.cluster.ET*neutrinoPt.Perp()*(1-cos(recoDeltaPhi)));  //real data
00726       float gMt=sqrt(2*mElectronP.Perp()*mNeutrinoP.Perp()*(1-cos(geantDeltaPhi)));
00727       
00728       hA[20]->Fill(Mt);  
00729       hA[21]->Fill(gMt);
00730 
00731       hA[22]->Fill(Mt,T.cluster.ET);
00732       hA[23]->Fill(gMt-Mt);
00733 
00734       //Kinematics 
00735       //reconstruct W pL from reconstructed quantities
00736       float trueWpL=mWP.z();
00737       float eleTheta=T.primP.Theta();
00738       float ratioE=T.cluster.energy/40.0;  
00739       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
00740       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
00741       hA[29]->Fill(pLRecoPlus);
00742       hA[30]->Fill(trueWpL-pLRecoPlus);
00743       hA[31]->Fill(trueWpL,pLRecoPlus);
00744       
00745       hA[32]->Fill(pLRecoMinus);
00746       hA[33]->Fill(trueWpL-pLRecoMinus);
00747       hA[34]->Fill(trueWpL,pLRecoMinus);
00748 
00749       const StMuTrack *prTr=T.prMuTrack; assert(prTr);
00750       float p_chrg=prTr->charge();
00751       if(p_chrg > 0) continue;
00752 
00753       float eleEta=T.primP.Eta();
00754       //sort 2 solutions by electron eta
00755       if(eleEta<-0.8) {
00756         hA[35]->Fill(trueWpL,pLRecoMinus);
00757         hA[37]->Fill(trueWpL,pLRecoPlus);
00758       }
00759       else if(eleEta>0.8) {
00760         hA[36]->Fill(trueWpL,pLRecoMinus);
00761         hA[38]->Fill(trueWpL,pLRecoPlus);
00762       }
00763       
00764       if(T.cluster.ET < 30) continue; //only W's we find in data
00765       //Correlate W pL with electron E in 3 electron eta ranges
00766       if(eleEta < -0.8) {
00767         hA[39]->Fill(mWP.z(),T.cluster.energy);
00768         hA[42]->Fill(T.cluster.energy);
00769       }
00770       if(eleEta > 0.8){
00771         hA[40]->Fill(mWP.z(),T.cluster.energy);
00772         hA[43]->Fill(T.cluster.energy);
00773       }      
00774       if(eleEta > -0.1 && eleEta < 0.1) {
00775         hA[41]->Fill(mWP.z(),T.cluster.energy);
00776         hA[44]->Fill(T.cluster.energy); 
00777       }
00778 
00779     }
00780   }
00781 }
00782 
00783 
00784 // $Log: St2009pubMcMaker.cxx,v $
00785 // Revision 1.7  2011/09/14 14:23:21  stevens4
00786 // update used for cross section PRD paper
00787 //
00788 // Revision 1.6  2010/05/03 19:54:35  stevens4
00789 // only try to calc effic if W->e+nu is found in McEvent
00790 //
00791 // Revision 1.5  2010/03/14 22:50:31  balewski
00792 // *** empty log message ***
00793 //
00794 // Revision 1.4  2010/02/19 21:04:10  stevens4
00795 // cleanup unused variables
00796 //
00797 // Revision 1.3  2010/02/18 19:48:10  stevens4
00798 // add more effic histograms and cleanup
00799 //
00800 // Revision 1.2  2010/01/21 17:54:31  stevens4
00801 // add effic histos and charge seperated background plots
00802 //
00803 // Revision 1.1  2009/11/23 23:00:18  balewski
00804 // code moved spin-pool
00805 //
00806 // Revision 1.1  2009/11/23 21:11:18  balewski
00807 // start
00808 //