StRoot: StSpinPool/StWalgo2011/St2011pubMcMaker.cxx Source File

00001 // $Id: St2011pubMcMaker.cxx,v 1.1 2011/02/10 20:33:25 balewski Exp $
00002 //
00003 //*-- Author : Justin Stevens, IUCF
00004 // 
00005 
00006 #include "St2011WMaker.h"
00007 #include "St2011pubMcMaker.h"
00008 #include "StEmcUtil/geometry/StEmcGeom.h"
00009 
00010 //need these to get MC record
00011 #include "tables/St_g2t_tpc_hit_Table.h"
00012 #include "StMcEventMaker/StMcEventMaker.h"
00013 #include "StMcEvent/StMcEvent.hh"
00014 #include "StMcEvent/StMcVertex.hh"
00015 #include "StMcEvent/StMcTrack.hh"
00016 
00017 ClassImp(St2011pubMcMaker)
00018 
00019 //_______________________________________________________
00020 //_______________________________________________________
00021 St2011pubMcMaker::St2011pubMcMaker(const char *name):StMaker(name){
00022   wMK=0;HList=0;
00023     
00024 }
00025 
00026 
00027 //_______________________________________________________
00028 //_______________________________________________________
00029 St2011pubMcMaker::~St2011pubMcMaker(){
00030   //
00031 }
00032 
00033 
00034 //_______________________________________________________
00035 //_______________________________________________________
00036 Int_t St2011pubMcMaker::Init(){
00037   assert(wMK);
00038   assert(HList);
00039   initHistos();
00040   return StMaker::Init();
00041 }
00042 
00043 
00044 //_______________________________________________________
00045 //_______________________________________________________
00046 Int_t 
00047 St2011pubMcMaker::Make(){
00048   //printf("-----------in %s\n", GetName());
00049   //only get geant particle info for W MC
00050   if(wMK->isMC==8 || wMK->isMC==30 || wMK->isMC==20){
00051     if(doMCanalysis()){
00052       doWanalysis();
00053       doWefficiency();
00054     }
00055   }
00056   return kStOK;
00057 }
00058 
00059 //_______________________________________________________
00060 //_______________________________________________________
00061 void 
00062 St2011pubMcMaker::doWanalysis(){
00063   //has access to whole W-algo-maker data via pointer 'wMK'
00064  
00065   // run through W cuts to fill other histos............
00066   for(uint iv=0;iv<wMK->wEve->vertex.size();iv++) {
00067     WeveVertex &V=wMK->wEve->vertex[iv];
00068     for(uint it=0;it<V.eleTrack.size();it++) {
00069       WeveEleTrack &T=V.eleTrack[it];
00070       if(T.isMatch2Cl==false) continue;
00071       assert(T.cluster.nTower>0); // internal logical error
00072       assert(T.nearTotET>0); // internal logical error
00073       
00074       if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
00075       if(T.awayTotET> 30.) continue; // too large awayET , Jan
00076       //Full W cuts applied at this point
00077 
00078       //W info from pythia record
00079       hA[1]->Fill(mWP.Perp());
00080       hA[2]->Fill(mWP.z());
00081       
00082       hA[24]->Fill(mWP.Eta());
00083       hA[25]->Fill(T.hadronicRecoil.Eta());
00084       hA[26]->Fill(mWP.Perp(),T.hadronicRecoil.Eta());
00085       if(mWP.Eta()>0) {
00086         hA[27]->Fill(T.hadronicRecoil.Eta());
00087       }
00088       else {
00089         hA[28]->Fill(T.hadronicRecoil.Eta());
00090       }
00091       
00092       //hadronic recoil and correlations with W from pythia
00093       TVector3 hadronicPt(T.hadronicRecoil.X(),T.hadronicRecoil.Y(),0); //transverse momentum vector
00094       hA[3]->Fill(mWP.Perp()-hadronicPt.Perp());
00095       hA[4]->Fill(mWP.Perp(),hadronicPt.Perp());
00096       hA[5]->Fill(hadronicPt.Perp()); 
00097 
00098       float delPhi=mWP.DeltaPhi(-hadronicPt);
00099       hA[6]->Fill(mWP.Perp(),delPhi);
00100       hA[7]->Fill(mWP.Perp(),mWP.Perp()-hadronicPt.Perp());
00101       hA[8]->Fill(T.hadronicRecoil.Perp(),delPhi);
00102       
00103       //electron reco and from geant record
00104       hA[9]->Fill(mElectronP.Perp());
00105       hA[10]->Fill(T.cluster.ET);      
00106       hA[11]->Fill(mElectronP.Perp(),T.cluster.ET);
00107       hA[12]->Fill(mElectronP.Perp(),mElectronP.Perp()-T.cluster.ET);
00108 
00109       TVector3 electronPt(T.cluster.position.X(),T.cluster.position.Y(),0); //transvers momentum vector
00110       electronPt.SetMag(T.cluster.ET);
00111       
00112       //neutrino reco and from geant record
00113       TVector3 neutrinoPt=-1*(hadronicPt+electronPt);
00114       hA[13]->Fill(neutrinoPt.Perp());  
00115       hA[14]->Fill(mNeutrinoP.Perp());
00116       hA[15]->Fill(mNeutrinoP.Perp(),neutrinoPt.Perp());
00117       hA[16]->Fill(mNeutrinoP.Perp()-neutrinoPt.Perp());
00118 
00119       hA[17]->Fill(mNeutrinoP.Perp(),mElectronP.Perp());
00120 
00121       float recoDeltaPhi=neutrinoPt.DeltaPhi(electronPt);
00122       float geantDeltaPhi=mNeutrinoP.DeltaPhi(mElectronP);
00123 
00124       hA[18]->Fill(geantDeltaPhi,recoDeltaPhi);
00125       hA[19]->Fill(cos(geantDeltaPhi)-cos(recoDeltaPhi));
00126 
00127       float Mt=sqrt(2*T.cluster.ET*neutrinoPt.Perp()*(1-cos(recoDeltaPhi)));  //real data
00128       float gMt=sqrt(2*mElectronP.Perp()*mNeutrinoP.Perp()*(1-cos(geantDeltaPhi)));
00129       
00130       hA[20]->Fill(Mt);  
00131       hA[21]->Fill(gMt);
00132 
00133       hA[22]->Fill(Mt,T.cluster.ET);
00134       hA[23]->Fill(gMt-Mt);
00135 
00136       //Kinematics 
00137       //reconstruct W pL from reconstructed quantities
00138       float trueWpL=mWP.z();
00139       float eleTheta=T.primP.Theta();
00140       float ratioE=T.cluster.energy/40.0;  
00141       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
00142       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
00143       hA[29]->Fill(pLRecoPlus);
00144       hA[30]->Fill(trueWpL-pLRecoPlus);
00145       hA[31]->Fill(trueWpL,pLRecoPlus);
00146       
00147       hA[32]->Fill(pLRecoMinus);
00148       hA[33]->Fill(trueWpL-pLRecoMinus);
00149       hA[34]->Fill(trueWpL,pLRecoMinus);
00150 
00151       const StMuTrack *prTr=T.prMuTrack; assert(prTr);
00152       float p_chrg=prTr->charge();
00153       if(p_chrg > 0) continue;
00154 
00155       float eleEta=T.primP.Eta();
00156       //sort 2 solutions by electron eta
00157       if(eleEta<-0.8) {
00158         hA[35]->Fill(trueWpL,pLRecoMinus);
00159         hA[37]->Fill(trueWpL,pLRecoPlus);
00160       }
00161       else if(eleEta>0.8) {
00162         hA[36]->Fill(trueWpL,pLRecoMinus);
00163         hA[38]->Fill(trueWpL,pLRecoPlus);
00164       }
00165       
00166       if(T.cluster.ET < 30) continue; //only W's we find in data
00167       //Correlate W pL with electron E in 3 electron eta ranges
00168       if(eleEta < -0.8) {
00169         hA[39]->Fill(mWP.z(),T.cluster.energy);
00170         hA[42]->Fill(T.cluster.energy);
00171       }
00172       if(eleEta > 0.8){
00173         hA[40]->Fill(mWP.z(),T.cluster.energy);
00174         hA[43]->Fill(T.cluster.energy);
00175       }      
00176       if(eleEta > -0.1 && eleEta < 0.1) {
00177         hA[41]->Fill(mWP.z(),T.cluster.energy);
00178         hA[44]->Fill(T.cluster.energy); 
00179       }
00180 
00181     }
00182   }
00183 }
00184 
00185 
00186 //_______________________________________________________
00187 //_______________________________________________________
00188 void 
00189 St2011pubMcMaker::doWefficiency(){
00190 
00191   //only count leptons in our eta range
00192   if(fabs(mElectronP.Eta()) > 1.) return;
00193   //ele has |eta| < 1
00194   hA[50]->Fill(mElectronP.Perp());
00195   hA[54]->Fill(mElectronP.Eta());
00196   hA[58]->Fill(mVertex.Z());
00197   hA[62]->Fill(mElectronP.Phi());
00198 
00199   hA[68]->Fill(mElectronP.Perp()); //forJoe
00200 
00201   //plot for Scott
00202   TVector3 detEle; //where lepton would hit BEMC
00203   float Rcylinder= wMK->mBtowGeom->Radius();
00204   detEle.SetPtEtaPhi(Rcylinder,mElectronP.Eta(),mElectronP.Phi());
00205   detEle.SetZ(detEle.Z()+mVertex.Z());
00206   hA[66]->Fill(detEle.Eta(),mElectronP.Perp());
00207 
00208   //trigger efficiency
00209   if(!wMK->wEve->l2bitET) return;
00210   //good trig
00211   hA[51]->Fill(mElectronP.Perp());
00212   hA[55]->Fill(mElectronP.Eta());
00213   hA[59]->Fill(mVertex.Z());
00214   hA[63]->Fill(mElectronP.Phi());
00215   hA[69]->Fill(mElectronP.Perp());//forJoe
00216 
00217   hA[67]->Fill(detEle.Eta(),mElectronP.Perp());
00218 
00219   //vertex efficiency
00220   if(wMK->wEve->vertex.size()<=0) return;
00221   //vertex rank>0 and |z|<100
00222   hA[52]->Fill(mElectronP.Perp()); 
00223   hA[56]->Fill(mElectronP.Eta());
00224   hA[60]->Fill(mVertex.Z());
00225   hA[64]->Fill(mElectronP.Phi());
00226 
00227   hA[70]->Fill(mElectronP.Perp());//forJoe
00228 
00229   //float diff=wMK->wEve->vertex[0].z - mVertex.Z();
00230   //cout<<"diff="<<diff<<endl;
00231 
00232   //reco efficiency
00233   for(uint iv=0;iv<wMK->wEve->vertex.size();iv++) {
00234     WeveVertex &V=wMK->wEve->vertex[iv];
00235     for(uint it=0;it<V.eleTrack.size();it++) {
00236       WeveEleTrack &T=V.eleTrack[it];
00237       if(T.isMatch2Cl==false) continue;
00238       assert(T.cluster.nTower>0); // internal logical error
00239       assert(T.nearTotET>0); // internal logical error
00240       
00241       if(T.cluster.ET/T.nearTotET< wMK->par_nearTotEtFrac) 
00242         continue; // too large nearET
00243       if(T.ptBalance.Perp() < wMK->par_ptBalance || T.awayTotET > 30.) //Jan
00244         continue;
00245       
00246       //pass all W cuts 
00247       hA[53]->Fill(mElectronP.Perp());
00248       hA[57]->Fill(mElectronP.Eta());
00249       hA[61]->Fill(mVertex.Z());
00250       hA[65]->Fill(mElectronP.Phi());
00251 
00252       hA[71]->Fill(mElectronP.Perp());//forJoe
00253     }
00254   }
00255 
00256 }
00257 
00258 
00259 //________________________________________________
00260 //________________________________________________
00261 bool
00262 St2011pubMcMaker::doMCanalysis(){
00263   StMcEvent* mMcEvent = 0;  
00264   mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
00265   assert(mMcEvent);
00266 
00267   //initialize momentum vectors
00268   StThreeVectorF pW;        float eW;
00269   StThreeVectorF pNeutrino; //float eNeutrino;
00270   StThreeVectorF pElectron; //float eElectron;
00271 
00272   StMcVertex *V=mMcEvent->primaryVertex(); 
00273   mVertex=TVector3(V->position().x(),V->position().y(),V->position().z());
00274   
00275   uint i=1;
00276   int found=0;
00277   while(found<2 && i<mMcEvent->tracks().size()){//loop tracks
00278     StMcTrack* mcTrack = mMcEvent->tracks()[i];
00279     int pdgId=mcTrack->pdgId();
00280     //float pt=mcTrack->pt();
00281     //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00282     if(pdgId==11 || pdgId==-11){ //select e+ and e-
00283       if(abs(mcTrack->parent()->pdgId()) == 24 ){ 
00284         pElectron=mcTrack->momentum();
00285         //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00286         pW=mcTrack->parent()->momentum();
00287         eW=mcTrack->parent()->energy();
00288         //LOG_INFO<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00289         found++;
00290       }
00291     }
00292     if(pdgId==12 || pdgId==-12){ //select neutrino
00293       if(abs(mcTrack->parent()->pdgId()) == 24 ){ 
00294         pNeutrino=mcTrack->momentum();
00295         //LOG_INFO<<"pdgId "<<pdgId<<" pt "<<pt<<" pz "<<mcTrack->momentum().z()<<endm;
00296         pW=mcTrack->parent()->momentum();
00297         eW=mcTrack->parent()->energy();
00298         //LOG_INFO<<"pdgId "<<mcTrack->parent()->pdgId()<<" pt "<<mcTrack->parent()->pt()<<" pz "<<mcTrack->parent()->momentum().z()<<endm;
00299         found++;
00300       }
00301     }
00302     i++;
00303   }
00304  
00305   if(found!=2) return false;
00306  
00307   mWP=TVector3(pW.x(),pW.y(),pW.z());
00308   mNeutrinoP=TVector3(pNeutrino.x(),pNeutrino.y(),pNeutrino.z());
00309   mElectronP=TVector3(pElectron.x(),pElectron.y(),pElectron.z());
00310   TVector3 diff=mWP-mNeutrinoP-mElectronP;
00311   if(diff.Mag() > 0.0001) //should get exactly right
00312     LOG_INFO<<"\n \n W+e+nu vector sum ="<<diff.Mag()<<endm;
00313   if(mElectronP.Mag() < 0.0001)
00314     LOG_INFO<<"\n \n no lepton track ="<<endm;
00315 
00316   //calculate x1 and x2 from W rapidity
00317   float rapW = 0.5*log((eW+mWP.Z())/(eW-mWP.Z()));
00318   float mw2sqs=80.4/500.;
00319   float x1=mw2sqs*exp(rapW);
00320   float x2=mw2sqs*exp(-rapW);
00321 
00322   hA[72]->Fill(rapW);
00323   if(fabs(mElectronP.Eta())<1){ //require midrapidity e
00324     hA[73]->Fill(x1);
00325     hA[74]->Fill(x2);
00326     hA[75]->Fill(x1,x2);
00327     hA[76]->Fill(x1-x2);
00328   }
00329 
00330   return true;
00331 
00332 }
00333 
00334 // $Log: St2011pubMcMaker.cxx,v $
00335 // Revision 1.1  2011/02/10 20:33:25  balewski
00336 // start
00337 //