StRoot: StSpinPool/StWalgoB2009/St2009W

00001 // $Id: St2009W_algo.cxx,v 1.25 2011/09/14 14:23:20 stevens4 Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 //*-- Author for Endcap: Justin Stevens, IUCF
00005 
00006 #include "StEmcUtil/geometry/StEmcGeom.h"
00007 #include "WeventDisplay.h"
00008 #include "StSpinPool/StJets/StJet.h"
00009 #include "StSpinPool/StJets/TowerToJetIndex.h"
00010 
00011 //muDst needed for zdcRate
00012 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00013 #include <StMuDSTMaker/COMMON/StMuDst.h>
00014 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00015 
00016 #include "St2009WMaker.h"
00017 
00018 //________________________________________________
00019 //________________________________________________
00020 void
00021 St2009WMaker::find_W_boson(){
00022   
00023   // printf("========= find_W_boson() \n");
00024   int nGoldW=0;
00025   //remove events tagged as Zs
00026   if(wEve.zTag) return;
00027 
00028   // search for  Ws ............
00029   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00030     WeveVertex &V=wEve.vertex[iv];
00031     for(uint it=0;it<V.eleTrack.size();it++) {
00032       WeveEleTrack &T=V.eleTrack[it];
00033       if(T.isMatch2Cl==false) continue;
00034       assert(T.cluster.nTower>0); // internal logical error
00035       assert(T.nearTotET>0); // internal logical error
00036       
00037       //make cut on lepton |eta| for cross section 
00038       if(fabs(T.primP.Eta()) > par_leptonEta) continue;      
00039       hA[0]->Fill("eta1",1.);
00040 
00041       //signal plots w/o EEMC in veto
00042       if(T.cluster.ET/T.nearTotET_noEEMC>par_nearTotEtFrac){
00043         if(T.awayTotET_noEEMC < 8)//old awayside pt cut
00044           hA[141]->Fill(T.cluster.ET);
00045         if(T.sPtBalance_noEEMC>par_ptBalance ) {//only signed ptBalance cut 
00046           hA[140]->Fill(T.cluster.ET);
00047           if (T.prMuTrack->charge() < 0) {
00048             hA[184+3]->Fill(T.cluster.ET);
00049             hA[200+3]->Fill(T.primP.Eta(),T.cluster.ET);
00050             hA[280+3]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00051           } else if (T.prMuTrack->charge() > 0) {
00052             hA[184+4]->Fill(T.cluster.ET);
00053             hA[200+4]->Fill(T.primP.Eta(),T.cluster.ET);
00054             hA[280+4]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00055           }
00056         }
00057       }
00058 
00059       if(T.cluster.ET /T.nearTotET< par_nearTotEtFrac) continue; // too large nearET
00060 
00061       hA[20]->Fill("noNear",1.);
00062       hA[113]->Fill( T.cluster.ET); // for Joe
00063       hA[50]->Fill(T.awayTpcPT);
00064       hA[51]->Fill(T.awayBtowET);
00065       hA[54]->Fill(T.awayTotET); 
00066       hA[52]->Fill(T.cluster.ET,T.awayTotET);
00067       hA[53]->Fill(T.cluster.ET,T.awayEmcET);
00068       hA[55]->Fill(T.awayEtowET);
00069       hA[60]->Fill(T.cluster.ET,T.awayTpcPT);
00070 
00071       hA[132]->Fill(T.cluster.ET,T.ptBalance.Perp());
00072       hA[133]->Fill(T.awayTotET,T.ptBalance.Perp());
00073       hA[134]->Fill(T.cluster.ET,T.sPtBalance);
00074       hA[135]->Fill(T.awayTotET,T.sPtBalance);
00075 
00076       int i=0;
00077       //background variations for systematic
00078       for (int j=0; j<=20; j++) { //loop over possible sPtBal cuts
00079         float pTBal_cut = 5.+((float) j);
00080         if (T.sPtBalance<pTBal_cut) { //background
00081           if (T.prMuTrack->charge() < 0) {
00082             hA[142+i]->Fill(T.cluster.ET,j);
00083             hA[210+j]->Fill(T.primP.Eta(),T.cluster.ET);
00084             hA[290+j]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00085           } else if (T.prMuTrack->charge() > 0) {
00086             hA[163+i]->Fill(T.cluster.ET,j);
00087             hA[231+j]->Fill(T.primP.Eta(),T.cluster.ET);
00088             hA[311+j]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00089           }
00090         }
00091       }
00092       
00093       //plots for backg sub yield
00094       if(T.sPtBalance>par_ptBalance ) {
00095         hA[136]->Fill(T.cluster.ET);//signal
00096         hA[62]->Fill(T.pointTower.iEta ,T.cluster.energy);
00097         if (T.prMuTrack->charge() < 0) {
00098           hA[184+1]->Fill(T.cluster.ET);
00099           hA[200+1]->Fill(T.primP.Eta(),T.cluster.ET);
00100           hA[280+1]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00101         } else if (T.prMuTrack->charge() > 0) {
00102           hA[184+2]->Fill(T.cluster.ET);
00103           hA[200+2]->Fill(T.primP.Eta(),T.cluster.ET);
00104           hA[280+2]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00105         }
00106       } else {
00107         hA[137]->Fill(T.cluster.ET);//background
00108         if (T.prMuTrack->charge() < 0) {
00109           hA[184+5]->Fill(T.cluster.ET);
00110           hA[200+5]->Fill(T.primP.Eta(),T.cluster.ET);
00111           hA[280+5]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00112         } else if (T.prMuTrack->charge() > 0) {
00113           hA[184+6]->Fill(T.cluster.ET);
00114           hA[200+6]->Fill(T.primP.Eta(),T.cluster.ET);
00115           hA[280+6]->Fill(abs(T.primP.Eta()),T.cluster.ET);
00116         }
00117       }
00118 
00119       //plots for backg sub yield (old awayTot cut DNP)
00120       if(T.awayTotET < 8)
00121         hA[138]->Fill(T.cluster.ET);//old signal
00122       else 
00123         hA[139]->Fill(T.cluster.ET);//old background
00124       
00125       if(0){/***************************/
00126           printf("\n WWWWWWWWWWWWWWWWWWWWW\n");
00127           wDisaply->exportEvent( "W", V, T);
00128           wEve.print();
00129         }/***************************/
00130  
00131         
00132       //put final W cut here
00133       if(T.sPtBalance<par_ptBalance)  continue;
00134       //::::::::::::::::::::::::::::::::::::::::::::::::
00135       //:::::accepted W events for x-section :::::::::::
00136       //::::::::::::::::::::::::::::::::::::::::::::::::
00137       wEve.wTag=true;
00138 
00139       hA[20]->Fill("noAway",1.0);  
00140       hA[114]->Fill( T.cluster.ET);//for Joe
00141 
00142       hA[90]->Fill( T.cluster.ET); 
00143       hA[92]->Fill( T.cluster.ET,T.glMuTrack->dEdx()*1e6); 
00144       hA[93]->Fill( T.cluster.ET,T.glMuTrack->dca().mag());
00145       int k=0; if(T.prMuTrack->charge()<0) k=1;
00146       hA[94+k]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
00147       // h95 used above
00148 
00149  
00150       if(T.cluster.ET<par_highET) continue;  // very likely Ws
00151       hA[91]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.position.Phi());
00152       hA[96]->Fill(V.id);
00153       hA[97]->Fill(V.funnyRank);
00154       hA[98]->Fill(V.z);
00155       hA[99]->Fill( T.prMuTrack->eta());
00156       hA[20]->Fill("goldW",1.);
00157       nGoldW++;
00158 
00159       //some final plots for comparing to embedding
00160       float zdcRate=mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
00161       hA[260]->Fill(zdcRate,wEve.bx7);
00162       hA[261]->Fill(zdcRate,wEve.vertex.size());
00163       hA[262]->Fill(zdcRate,T.prMuTrack->nHitsFit());
00164       float hitFrac=1.*T.prMuTrack->nHitsFit()/T.prMuTrack->nHitsPoss();
00165       hA[263]->Fill(zdcRate,hitFrac);
00166       hA[264]->Fill(zdcRate,T.prMuTrack->globalTrack()->lastPoint().perp());
00167       hA[265]->Fill(zdcRate,T.prMuTrack->globalTrack()->firstPoint().perp());
00168       hA[266]->Fill(zdcRate,T.prMuTrack->pt());
00169       hA[267]->Fill(zdcRate,T.prMuTrack->globalTrack()->chi2());
00170       hA[268]->Fill(zdcRate,1./T.prMuTrack->pt());
00171       hA[269]->Fill(zdcRate,T.glMuTrack->dca().mag());
00172       hA[270]->Fill(zdcRate,T.glMuTrack->dcaD());
00173       hA[271]->Fill(zdcRate,T.glMuTrack->dcaZ());
00174 
00175       if(V.funnyRank<10) {
00176         hA[272]->Fill(zdcRate,T.glMuTrack->pt());
00177         hA[273]->Fill(zdcRate,1./T.glMuTrack->pt());
00178       }
00179 
00180     }// loop over tracks
00181   }// loop over vertices
00182   if(nGoldW>0)
00183     hA[0]->Fill("goldW",1.);
00184 
00185 }
00186 
00187 
00188 //________________________________________________
00189 //________________________________________________
00190 void
00191 St2009WMaker::tag_Z_boson(){
00192 
00193   float par_jetPt=10.;
00194   float lowMass=70.; float highMass=140.;
00195   mJets = getJets("ConeJets12_100"); //select specific jet-type
00196 
00197   //form invariant mass from lepton candidate and jet
00198   for(uint iv=0;iv<wEve.vertex.size();iv++) {//vertex loop
00199     WeveVertex &V=wEve.vertex[iv];
00200     for(uint it=0;it<V.eleTrack.size();it++) {// select track
00201       WeveEleTrack &T1=V.eleTrack[it];
00202       if(T1.isMatch2Cl==false) continue;
00203       assert(T1.cluster.nTower>0); // internal logical error
00204       assert(T1.nearTotET>0); // internal logical error
00205       if(T1.cluster.ET/T1.nearTotET< par_nearTotEtFrac) continue; // too large nearET
00206 
00207       //match lepton candidate with jet
00208       TLorentzVector jetVec;
00209       for (int i_jet=0; i_jet< nJets; i_jet++){//jet loop
00210         jetVec = *((StJet*)mJets->At(i_jet));
00211         if(jetVec.Pt()<par_jetPt) continue;//remove low pt jets
00212         
00213         //electron like cut on jets
00214         StJet* jet = getJet(i_jet);  float maxCluster=0.; 
00215         int totTowers=jet->nBtowers+jet->nEtowers;
00216         for(int itow=0;itow<totTowers;itow++){//loop over towers
00217           if(jet->tower(itow)->detectorId()==13)//drop endcap towers
00218             continue;
00219           
00220           int softId=jet->tower(itow)->towerId();
00221           //find highest 2x2 BTOW cluster in jet
00222           TVector3 pos=positionBtow[softId-1]; int iEta,iPhi;
00223           if( L2algoEtaPhi2IJ(pos.Eta(),pos.Phi(),iEta,iPhi)) 
00224             continue;
00225           float cluster=maxBtow2x2(iEta,iPhi,jet->zVertex).ET;
00226           if(cluster>maxCluster) maxCluster=cluster;
00227         }
00228 
00229         TVector3 jetVec3(jetVec.X(),jetVec.Y(),jetVec.Z());
00230         if(jetVec3.DeltaR(T1.primP)<par_nearDeltaR)
00231           continue;//skip jets in candidate phi isolation'cone'
00232         
00233         //form invM
00234         float e1=T1.cluster.energy;
00235         TVector3 p1=T1.primP; p1.SetMag(e1);
00236         TLorentzVector ele1(p1,e1); //lepton candidate 4- mom
00237         TLorentzVector sum=ele1+jetVec;
00238         float invM=sqrt(sum*sum);
00239         if(maxCluster/jet->jetPt < 0.5) continue;
00240         if(invM > lowMass && invM < highMass)
00241           wEve.zTag=true;
00242       }
00243     }
00244   }
00245   
00246 }
00247 
00248 //________________________________________________
00249 //________________________________________________
00250 void
00251 St2009WMaker::findPtBalance(){
00252 
00253   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00254     WeveVertex &V=wEve.vertex[iv];
00255     for(uint it=0;it<V.eleTrack.size();it++) {
00256       WeveEleTrack &T=V.eleTrack[it];
00257       if(T.isMatch2Cl==false) continue;
00258 
00259       //****loop over branch with EEMC****
00260       mJets = getJets(mJetTreeBranch);
00261       int nJetsWE=nJets;
00262       for (int i_jet=0; i_jet< nJetsWE; i_jet++){//loop over jets
00263         StJet* jet = getJet(i_jet);
00264         TVector3 jetVec; //vector for jet momentum
00265         //vary neutral and charged et in jets for systematic
00266         float neutral=jet->neutralFraction()*jet->Pt();
00267         float charged=jet->chargedFraction()*jet->Pt();
00268         neutral=neutral*par_mcJetNeutScale;
00269         charged=charged*par_mcJetChrgScale;
00270         float sum=neutral+charged;
00271         jetVec.SetPtEtaPhi(sum,jet->Eta(),jet->Phi());
00272         if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
00273               T.ptBalance+=jetVec;
00274       }
00275       TVector3 clustPt(T.primP.X(),T.primP.Y(),0);
00276       clustPt.SetMag(T.cluster.ET);
00277       T.ptBalance+=clustPt;
00278       T.sPtBalance = T.ptBalance.Perp();
00279       if(T.ptBalance.Dot(clustPt)<0) T.sPtBalance *=-1.;
00280       
00281       //****loop over branch without EEMC****
00282       mJets = getJets(mJetTreeBranch_noEEMC);
00283       int nJetsNE=nJets;
00284       TVector3 highJet; //zero out highJet 
00285       for (int i_jet=0; i_jet< nJetsNE; i_jet++){//loop over jets
00286         StJet* jet = getJet(i_jet);
00287         TVector3 jetVec; //vector for jet momentum
00288         //vary neutral and charged et in jets for systematic
00289         float neutral=jet->neutralFraction()*jet->Pt();
00290         float charged=jet->chargedFraction()*jet->Pt();
00291         neutral=neutral*par_mcJetNeutScale;
00292         charged=charged*par_mcJetChrgScale;
00293         float sum=neutral+charged;
00294         jetVec.SetPtEtaPhi(sum,jet->Eta(),jet->Phi());
00295         if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
00296           T.ptBalance_noEEMC+=jetVec;
00297       }
00298       T.ptBalance_noEEMC+=clustPt;
00299       T.sPtBalance_noEEMC = T.ptBalance_noEEMC.Perp();
00300       if(T.ptBalance_noEEMC.Dot(clustPt)<0) T.sPtBalance_noEEMC *=-1.;
00301 
00302     }// end of loop over tracks
00303   }// end of loop over vertices
00304   
00305 }
00306 
00307 
00308 //________________________________________________
00309 //________________________________________________
00310 void
00311 St2009WMaker::findAwayJet(){
00312   // printf("\n******* find AwayJet() nVert=%d\n",wEve.vertex.size());
00313   //wEve.print();
00314   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00315     WeveVertex &V=wEve.vertex[iv];
00316     for(uint it=0;it<V.eleTrack.size();it++) {
00317       WeveEleTrack &T=V.eleTrack[it];
00318       if(T.isMatch2Cl==false) continue;
00319       
00320       // .... sum opposite in phi EMC components
00321       T.awayBtowET=sumBtowCone(V.z,-T.primP,1,T.awayNTow); // '1'= only cut on delta phi
00322       T.awayEmcET=T.awayBtowET;
00323       T.awayEtowET=sumEtowCone(V.z,-T.primP,1,T.awayNTow);
00324       if(par_useEtow >= 2) T.awayEmcET+=T.awayEtowET;
00325 
00326       //..... add TPC ET
00327       T.awayTpcPT=sumTpcCone(V.id,-T.primP,1,T.awayNTr);
00328       T.awayTotET=T.awayEmcET+T.awayTpcPT;
00329       T.awayTotET_noEEMC=T.awayBtowET+T.awayTpcPT;
00330       //printf("\n*** in   awayTpc=%.1f awayEmc=%.1f\n  ",T.awayTpcPT,T.awayEmcET); T.print(); 
00331     }// end of loop over tracks
00332   }// end of loop over vertices
00333 }
00334 
00335 //________________________________________________
00336 //________________________________________________
00337 void
00338 St2009WMaker::findNearJet(){
00339   //printf("\n******* findNearJet() nVert=%d\n",wEve.vertex.size());
00340 
00341   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00342     WeveVertex &V=wEve.vertex[iv];
00343     for(uint it=0;it<V.eleTrack.size();it++) {
00344       WeveEleTrack &T=V.eleTrack[it];
00345       if(T.isMatch2Cl==false) continue;
00346 
00347        // .... sum EMC-jet component
00348       T.nearBtowET=sumBtowCone(V.z,T.primP,2,T.nearNTow); // '2'=2D cone
00349       T.nearEmcET+=T.nearBtowET;
00350       T.nearEtowET=sumEtowCone(V.z,T.primP,2,T.nearNTow);
00351       if(par_useEtow >= 3) T.nearEmcET+=T.nearEtowET;
00352       // .... sum TPC-near component
00353       T.nearTpcPT=sumTpcCone(V.id,T.primP,2,T.nearNTr); // '2'=2D cone
00354       hA[47]->Fill(T.nearTpcPT);
00355       hA[48]->Fill(T.nearEmcET,T.nearTpcPT);
00356       float nearSum=T.nearEmcET+T.nearTpcPT;
00357       
00358       /* correct for double counting of electron track in near cone
00359          rarely primTrPT<10 GeV & globPT>10 - handle this here */
00360       if(T.primP.Pt()>par_trackPt) nearSum-=par_trackPt; 
00361       else  nearSum-=T.primP.Pt();
00362       T.nearTotET=nearSum;
00363       T.nearTotET_noEEMC=nearSum-T.nearEtowET;
00364 
00365       hA[49]->Fill(nearSum);
00366 
00367       // printf("\n*** in nearCone T.print\n  "); T.print();
00368       hA[40]->Fill( T.nearEmcET);
00369       hA[41]->Fill( T.cluster.ET,T.nearEmcET-T.cluster.ET);
00370       float nearTotETfrac=T.cluster.ET/ T.nearTotET;
00371       hA[42]->Fill(nearTotETfrac);
00372       hA[192]->Fill(T.cluster.ET,nearTotETfrac);
00373   
00374     } // end of loop over tracks
00375   }// end of loop over vertices
00376  
00377 
00378 }
00379 
00380 
00381 //________________________________________________
00382 //________________________________________________
00383 int
00384 St2009WMaker::matchTrack2Cluster(){
00385   // printf("******* matchCluster() nVert=%d\n",wEve.vertex.size());
00386   //wEve.print();
00387   int nTr=0;
00388   float Rcylinder= mBtowGeom->Radius();
00389   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00390     WeveVertex &V=wEve.vertex[iv];
00391     float zVert=V.z;
00392     for(uint it=0;it<V.eleTrack.size();it++) {
00393       WeveEleTrack &T=V.eleTrack[it];
00394       T.isMatch2Cl=false; // just in case (was already set in constructor)
00395       if(T.pointTower.id==0) continue;
00396       
00397       float trackPT=T.prMuTrack->momentum().perp();
00398       T.cluster=maxBtow2x2( T.pointTower.iEta, T.pointTower.iPhi,zVert);      
00399       hA[33]->Fill( T.cluster.ET);
00400       hA[34]->Fill(T.cluster.adcSum,trackPT);
00401       hA[110]->Fill( T.cluster.ET);
00402       if (T.cluster.ET <par_clustET) continue; // too low energy
00403       hA[20]->Fill("CL",1.);
00404 
00405       //.. spacial separation (track - cluster)
00406       TVector3 D=T.pointTower.R-T.cluster.position;
00407       hA[43]->Fill( T.cluster.energy,D.Mag());
00408       hA[44]->Fill( T.cluster.position.z(),D.z());
00409       float delPhi=T.pointTower.R.DeltaPhi(T.cluster.position);
00410       //   printf("aaa %f %f %f   phi=%f\n",D.x(),D.y(),D.z(),delPhi);
00411       hA[45]->Fill( T.cluster.energy,Rcylinder*delPhi);// wrong?
00412       hA[46]->Fill( D.Mag());
00413 
00414       if(D.Mag()>par_delR3D) continue; 
00415       hA[20]->Fill("#Delta R",1.);
00416       hA[111]->Fill( T.cluster.ET);
00417 
00418       // ........compute surroinding cluster energy
00419       int iEta=T.cluster.iEta;
00420       int iPhi=T.cluster.iPhi;
00421       T.cl4x4=sumBtowPatch(iEta-1,iPhi-1,4,4,zVert); // needed for lumi monitor
00422       hA[37]->Fill( T.cl4x4.ET);
00423       hA[38]->Fill(T.cluster.energy, T.cl4x4.energy-T.cluster.energy);
00424       
00425       //sum tpc PT  near lepton candidate track
00426       T.smallNearTpcPT=sumTpcCone(V.id,T.primP,3,T.smallNearNTr);
00427       hA[56]->Fill(T.smallNearTpcPT);
00428       T.smallNearTpcPT-=par_trackPt;
00429 
00430       float frac24=T.cluster.ET/(T.cl4x4.ET);
00431       hA[39]->Fill(frac24);
00432       hA[191]->Fill(T.cluster.ET,frac24);
00433       if(frac24<par_clustFrac24) continue;
00434       T.isMatch2Cl=true; // cluster is matched to TPC track
00435 
00436       hA[20]->Fill("fr24",1.);
00437       hA[112]->Fill( T.cluster.ET);
00438       
00439       nTr++;
00440     }// end of one vertex
00441   }// end of vertex loop
00442 
00443   if(nTr<=0) return -1; // cancel event w/o good tracks
00444   hA[0]->Fill("Tr2Cl",1.0);
00445   return 0;
00446 
00447 }
00448 
00449 //________________________________________________
00450 //________________________________________________ 
00451 WeveCluster 
00452 St2009WMaker::maxBtow2x2(int iEta, int iPhi, float zVert){
00453   //printf("   maxBtow2x2  seed iEta=%d iPhi=%d \n",iEta, iPhi);
00454   const int L=2; // size of the summed square 
00455 
00456   WeveCluster maxCL;
00457   // just 4 cases of 2x2 clusters
00458   float maxET=0;
00459   int I0=iEta-1;
00460   int J0=iPhi-1;
00461   for(int I=I0;I<=I0+1;I++){
00462     for(int J=J0;J<=J0+1;J++) {
00463       WeveCluster CL=sumBtowPatch(I,J,L,L,zVert);
00464       if(maxET>CL.ET) continue;
00465       maxET=CL.ET;
00466       maxCL=CL;
00467       // printf("   newMaxETSum=%.1f iEta=%d iPhi=%d \n",maxET, I,J);
00468     }
00469   }// 4 combinations done
00470   //printf(" final inpEve=%d SumET2x2=%.1f \n",nInpEve,maxET);
00471   return maxCL;  
00472 }
00473 
00474 
00475 //________________________________________________
00476 //________________________________________________
00477 WeveCluster
00478 St2009WMaker::sumBtowPatch(int iEta, int iPhi, int Leta,int  Lphi, float zVert){
00479   //printf("  eveID=%d btowSquare seed iEta=%d[+%d] iPhi=%d[+%d] zVert=%.0f \n",wEve.id,iEta,Leta, iPhi,Lphi,zVert);
00480   WeveCluster CL; // object is small, not to much overhead in creating it
00481   CL.iEta=iEta;
00482   CL.iPhi=iPhi;
00483   TVector3 R; 
00484   double sumW=0;
00485   float Rcylinder= mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
00486   for(int i=iEta; i<iEta+Leta;i++){// trim in eta-direction
00487     if(i<0) continue;
00488     if(i>=mxBTetaBin) continue;
00489     for(int j=iPhi;j<iPhi+Lphi;j++) {// wrap up in the phi-direction
00490       int jj=(j+mxBTphiBin)%mxBTphiBin;// keep it always positive
00491       //if(L<5) printf("n=%2d  i=%d jj=%d\n",CL.nTower,i,jj);
00492       int softID= mapBtowIJ2ID[ i+ jj*mxBTetaBin];
00493       float ene= wEve.bemc.eneTile[kBTow][softID-1];
00494       if(ene<=0) continue; // skip towers w/o energy
00495       float adc= wEve.bemc.adcTile[kBTow][softID-1];
00496       float delZ=positionBtow[softID-1].z()-zVert;
00497       float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
00498       float ET=ene*e2et;
00499       float logET=log10(ET+0.5);
00500       CL.nTower++;
00501       CL.energy+=ene;
00502       CL.ET+=ET;
00503       CL.adcSum+=adc;
00504       if(logET>0) {
00505         R+=logET*positionBtow[softID-1];
00506         sumW+=logET;
00507       }
00508       //if(Leta==2) printf("      iEta=%d iPhi=%d  ET=%.1f  ene=%.1f   sum=%.1f logET=%f sumW=%f\n",i,j,ET,ene,CL.energy,logET,sumW); 
00509     }
00510     // printf(" end btowSquare: iEta=%d  nTw=%d, ET=%.1f adc=%.1f\n",i,CL.nTower,CL.ET,CL.adcSum);
00511     if(sumW>0) {
00512       CL.position=1./sumW*R; // weighted cluster position
00513     } else {
00514       CL.position=TVector3(0,0,999);
00515     } 
00516   }
00517   return CL;
00518 }
00519 
00520 //________________________________________________
00521 //________________________________________________
00522 int
00523 St2009WMaker::extendTrack2Barrel(){// return # of extended tracks
00524   //wEve.print();
00525   //printf("******* extendTracks() nVert=%d\n",wEve.vertex.size());
00526   int nTrB=0; 
00527   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00528     WeveVertex &V=wEve.vertex[iv];
00529     for(uint it=0;it<V.eleTrack.size();it++) {
00530       WeveEleTrack &T=V.eleTrack[it];
00531       //.... extrapolate track to the barrel @ R=entrance....
00532       const StPhysicalHelixD TrkHlx=T.prMuTrack->outerHelix(); 
00533       float Rcylinder= mBtowGeom->Radius();
00534       pairD  d2; 
00535       d2 = TrkHlx.pathLength(Rcylinder);
00536       //printf(" R=%.1f path 1=%f, 2=%f, period=%f, R=%f\n",Rctb,d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
00537 
00538       // assert(d2.first<0); // propagate backwards
00539       // assert(d2.second>0); // propagate forwards
00540       if(d2.first>=0 || d2.second<=0) {
00541         LOG_WARN<< Form("MatchTrk , unexpected solution for track crossing CTB\n  d2.firts=%f, second=%f, swap them",  d2.first, d2.second)<<endm;
00542         float xx=d2.first;
00543         d2.first=d2.second;
00544         d2.second=xx;
00545       }
00546 
00547       // extrapolate track to cylinder
00548       StThreeVectorD posR = TrkHlx.at(d2.second);
00549       //printf(" punch2 x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f\n",posCTB.x(),posCTB.y(),posCTB.z(),xmagn);
00550 
00551       float etaF=posR.pseudoRapidity();
00552       float phiF=posR.phi();
00553       int iEta,iPhi;
00554       if( L2algoEtaPhi2IJ(etaF, phiF,iEta,iPhi)) continue;
00555       nTrB++;
00556        hA[20]->Fill("@B",1.);
00557       //printf(" phi=%.0f deg,  eta=%.2f, iEta=%d, iPhi=%d\n",posCTB.phi()/3.1416*180.,posCTB. pseudoRapidity(),iEta, iPhi);
00558       int twID= mapBtowIJ2ID[ iEta+ iPhi*mxBTetaBin];
00559       // printf("hit Tower ID=%d\n",twID);
00560 
00561       T.pointTower.id=twID;
00562       T.pointTower.R=TVector3(posR.x(),posR.y(),posR.z());
00563       T.pointTower.iEta=iEta;
00564       T.pointTower.iPhi=iPhi;
00565     }
00566   }// end of loop over vertices
00567 
00568   if(nTrB<=0) return -1; // cancel event w/o good tracks
00569   hA[0]->Fill("TrB",1.0);
00570   return 0;
00571 
00572 }
00573 
00574 
00575 //________________________________________________
00576 //________________________________________________
00577 //________________________________________________
00578 float
00579 St2009WMaker::sumBtowCone( float zVert,  TVector3 refAxis, int flag, int &nTow){ 
00580   /* flag=1 : only delta phi cut;  flag=2 use 2D cut */
00581   assert(flag==1 || flag==2);
00582   double ptSum=0;
00583   
00584   //.... process BTOW hits
00585   for(int i=0;i< mxBtow;i++) {
00586     float ene=wEve.bemc.eneTile[kBTow][i];
00587     if(ene<=0) continue;
00588     TVector3 primP=positionBtow[i]-TVector3(0,0,zVert);
00589     primP.SetMag(ene); // it is 3D momentum in the event ref frame
00590     if(flag==1) {
00591       float deltaPhi=refAxis.DeltaPhi(primP);
00592       if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
00593     }
00594     if(flag==2) {
00595       float deltaR=refAxis.DeltaR(primP);
00596       if(deltaR> par_nearDeltaR) continue;
00597     }
00598     if(primP.Perp()>par_countTowEt) nTow++;
00599     ptSum+=primP.Perp();
00600   }
00601   
00602   return ptSum;  
00603 }
00604 
00605 
00606 //________________________________________________
00607 //________________________________________________
00608 float
00609 St2009WMaker::sumEtowCone(float zVert, TVector3 refAxis, int flag,int &nTow){
00610   /* flag=1 : only delta phi cut;  flag=2 use 2D cut */
00611   assert(flag==1 || flag==2);
00612   
00613   float ptsum=0;
00614   
00615   //....loop over all phi bins
00616   for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
00617     for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
00618       float ene=wEve.etow.ene[iphi][ieta];
00619       if(ene<=0) continue; //skip towers with no energy
00620       TVector3 primP=positionEtow[iphi][ieta]-TVector3(0,0,zVert);
00621       primP.SetMag(ene); // it is 3D momentum in the event ref frame
00622       if(flag==1) {
00623         float deltaPhi=refAxis.DeltaPhi(primP);
00624         if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
00625       }
00626       if(flag==2) {
00627         float deltaR=refAxis.DeltaR(primP);
00628         if(deltaR> par_nearDeltaR) continue;
00629       }
00630       if(primP.Perp()>par_countTowEt) nTow++;
00631       ptsum+=primP.Perp();
00632     }
00633   }
00634   
00635   return ptsum;
00636 }
00637 
00638 // $Log: St2009W_algo.cxx,v $
00639 // Revision 1.25  2011/09/14 14:23:20  stevens4
00640 // update used for cross section PRD paper
00641 //
00642 // Revision 1.24  2010/04/27 16:53:45  stevens4
00643 // add code to remove events tagged as Zs from W candidates
00644 //
00645 // Revision 1.23  2010/04/16 14:35:32  balewski
00646 // fix borken header
00647 //
00648 // Revision 1.22  2010/03/23 01:31:40  seelej
00649 // Fix to the filling of the histograms for the background systematic.
00650 //
00651 // Revision 1.21  2010/03/22 01:45:58  seelej
00652 // *** empty log message ***
00653 //
00654 // Revision 1.20  2010/03/22 01:33:13  seelej
00655 // Additional change.
00656 //
00657 // Revision 1.19  2010/03/22 00:55:45  seelej
00658 // Change to use the signed pT balance in background subtraction instead of the unsigned pT balance.
00659 //
00660 // Revision 1.18  2010/03/18 18:46:40  balewski
00661 // simplified sPtBalance calculation
00662 //
00663 // Revision 1.17  2010/03/18 16:52:17  balewski
00664 // corrected sPtBalance for no-endcap
00665 //
00666 // Revision 1.16  2010/03/18 15:34:43  seelej
00667 // Changed the definition of sPtBalance - Joe
00668 //
00669 // Revision 1.15  2010/03/14 22:50:31  balewski
00670 // *** empty log message ***
00671 //
00672 // Revision 1.14  2010/03/12 21:08:11  balewski
00673 // simplify logic for filling histos for noE background
00674 //
00675 // Revision 1.13  2010/02/26 21:40:00  seelej
00676 // Joe : Fix to code. Forgot to remove an older piece of code when doing a previous update.
00677 //
00678 // Revision 1.12  2010/02/22 15:49:34  seelej
00679 // Joe : Changes to code for inclusion of background subtraction and systematic studies
00680 //
00681 // Revision 1.11  2010/01/28 03:42:55  balewski
00682 // cleanup
00683 //
00684 // Revision 1.10  2010/01/27 22:12:24  balewski
00685 // spin code matched to x-section code
00686 //
00687 // Revision 1.9  2010/01/26 18:48:11  stevens4
00688 // include |eta|<1 cut for all W yields
00689 //
00690 // Revision 1.8  2010/01/23 02:35:38  stevens4
00691 // add ability to scale jet et and use real btow peds for rcf mc
00692 //
00693 // Revision 1.7  2010/01/10 03:01:37  balewski
00694 // cleanup & nicer histos
00695 //
00696 // Revision 1.6  2010/01/10 01:45:10  stevens4
00697 // fix plots w/o EEMC in veto
00698 //
00699 // Revision 1.5  2010/01/09 00:07:16  stevens4
00700 // add jet finder
00701 //
00702 // Revision 1.4  2010/01/06 19:16:47  stevens4
00703 // track cuts now on primary component, cleanup
00704 //
00705 // Revision 1.3  2010/01/06 14:11:13  balewski
00706 // one Z-plot added
00707 //
00708 // Revision 1.2  2010/01/06 04:22:15  balewski
00709 // added Q/PT plot for Zs, more cleanup
00710 //
00711 // Revision 1.1  2009/11/23 23:00:18  balewski
00712 // code moved spin-pool
00713 //
StSpinPool/StWalgoB2009/St2009W_algo.cxx
