StRoot: StSpinPool/StWalgo2011/St2011W

00001 // $Id: St2011W_algo.cxx,v 1.3 2011/02/25 06:03:47 stevens4 Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 //*-- Author for Endcap: Justin Stevens, IUCF
00005 
00006 #include "TF1.h"
00007 
00008 #include "StEmcUtil/geometry/StEmcGeom.h"
00009 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00010 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00011 #include "WeventDisplay.h"
00012 #include "StSpinPool/StJets/StJet.h"
00013 #include "StSpinPool/StJets/TowerToJetIndex.h"
00014 
00015 #include "St2011WMaker.h"
00016 
00017 //________________________________________________
00018 //________________________________________________
00019 void
00020 St2011WMaker::find_W_boson(){
00021   
00022   if(!wEve->l2bitET) return;
00023 
00024   //printf("========= find_W_boson() \n");
00025   int nNoNear=0,nNoAway=0,nEta1=0,nGoldW=0;
00026   //remove events tagged as Zs
00027   if(wEve->zTag) return;
00028 
00029   // search for  Ws ............
00030   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00031     WeveVertex &V=wEve->vertex[iv];
00032     for(uint it=0;it<V.eleTrack.size();it++) {
00033       WeveEleTrack &T=V.eleTrack[it];
00034       if(T.pointTower.id<=0) continue; //skip endcap towers
00035       if(T.isMatch2Cl==false) continue;
00036       assert(T.cluster.nTower>0); // internal logical error
00037       assert(T.nearTotET>0); // internal logical error
00038       
00039       //make cut on lepton |eta|  
00040       if(fabs(T.primP.Eta()) > par_leptonEta) continue;      
00041       hA[20]->Fill("eta1",1.);
00042       nEta1++;
00043 
00044       //signal plots w/o EEMC in veto
00045       if(T.cluster.ET/T.nearTotET_noEEMC>par_nearTotEtFrac){
00046         if(T.sPtBalance_noEEMC>par_ptBalance ) {//only signed ptBalance cut 
00047           hA[140]->Fill(T.cluster.ET);
00048           hA[240]->Fill(T.prMuTrack->eta(),T.cluster.ET);
00049           if (T.prMuTrack->charge() < 0) {
00050             hA[184+3]->Fill(T.cluster.ET);
00051           } else if (T.prMuTrack->charge() > 0) {
00052             hA[184+4]->Fill(T.cluster.ET);
00053           }
00054         }
00055       }
00056 
00057       if(T.cluster.ET /T.nearTotET< par_nearTotEtFrac) continue; // too large nearET
00058 
00059       hA[20]->Fill("noNear",1.);
00060       nNoNear++;
00061       hA[112]->Fill( T.cluster.ET); // for Joe
00062       hA[50]->Fill(T.awayTpcPT);
00063       hA[51]->Fill(T.awayBtowET);
00064       hA[54]->Fill(T.awayTotET); 
00065       hA[52]->Fill(T.cluster.ET,T.awayTotET);
00066       hA[53]->Fill(T.cluster.ET,T.awayEmcET);
00067       hA[55]->Fill(T.awayEtowET);
00068       hA[60]->Fill(T.cluster.ET,T.awayTpcPT);
00069 
00070       hA[132]->Fill(T.cluster.ET,T.ptBalance.Perp());
00071       hA[133]->Fill(T.awayTotET,T.ptBalance.Perp());
00072       hA[134]->Fill(T.cluster.ET,T.sPtBalance);
00073       hA[135]->Fill(T.awayTotET,T.sPtBalance);
00074       hA[209]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
00075 
00076       for (int i=0; i<=20; i++) {
00077         //  float awayTot_cut = 10.+2.*((float) i);
00078         for (int j=0; j<=20; j++) {
00079           float pTBal_cut = 5.+((float) j);
00080           if (T.sPtBalance<pTBal_cut) { 
00081             if (T.prMuTrack->charge() < 0) {
00082               hA[142+i]->Fill(T.cluster.ET,j);
00083             } else if (T.prMuTrack->charge() > 0) {
00084               hA[163+i]->Fill(T.cluster.ET,j);
00085             }
00086           }
00087         }
00088       }
00089 
00090       //plots for backg sub yield
00091       if(T.sPtBalance>par_ptBalance ) {
00092         hA[136]->Fill(T.cluster.ET);//signal
00093         hA[241]->Fill(T.prMuTrack->eta(),T.cluster.ET);
00094         hA[62]->Fill(T.pointTower.iEta ,T.cluster.energy);
00095         if (T.prMuTrack->charge() < 0) {
00096           hA[184+1]->Fill(T.cluster.ET);
00097         } else if (T.prMuTrack->charge() > 0) {
00098           hA[184+2]->Fill(T.cluster.ET);
00099         }
00100       } else {
00101         hA[137]->Fill(T.cluster.ET);//background
00102         if (T.prMuTrack->charge() < 0) {
00103           hA[184+5]->Fill(T.cluster.ET);
00104         } else if (T.prMuTrack->charge() > 0) {
00105           hA[184+6]->Fill(T.cluster.ET);
00106         }
00107         hA[202]->Fill(T.cluster.position.PseudoRapidity(),T.prMuTrack->pt());
00108         hA[204]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.energy/T.prMuTrack->p().mag());
00109       }
00110 
00111       if(0){/***************************/
00112         printf("\n WWWWWWWWWWWWWWWWWWWWW  Barrel \n");
00113         wDisaply->exportEvent( "WB", V, T);
00114         wEve->print();
00115       }/***************************/
00116  
00117       
00118       //put final W cut here
00119       if(T.sPtBalance<par_ptBalance)  continue;
00120       hA[20]->Fill("noAway",1.0);
00121       nNoAway++;
00122 
00123       //::::::::::::::::::::::::::::::::::::::::::::::::
00124       //:::::accepted W events for x-section :::::::::::
00125       //::::::::::::::::::::::::::::::::::::::::::::::::
00126 
00127       hA[113]->Fill( T.cluster.ET);//for Joe
00128 
00129       hA[90]->Fill( T.cluster.ET); 
00130       hA[92]->Fill( T.cluster.ET,T.glMuTrack->dEdx()*1e6); 
00131       hA[93]->Fill( T.cluster.ET,T.glMuTrack->dca(V.id).mag());
00132       int k=0; if(T.prMuTrack->charge()<0) k=1;
00133       hA[94+k]->Fill( T.cluster.ET,T.glMuTrack->dcaD());
00134       // h95 used above
00135 
00136       //plots to investigate east/west yield diff
00137       hA[200]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
00138       hA[201]->Fill(T.cluster.position.PseudoRapidity(),T.prMuTrack->pt());
00139       hA[203]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.energy/T.prMuTrack->p().mag());
00140       hA[205]->Fill(T.prMuTrack->lastPoint().pseudoRapidity(),T.prMuTrack->lastPoint().phi());
00141       
00142 
00143       if(T.cluster.ET<par_highET) continue;  // very likely Ws
00144       hA[91]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.position.Phi());
00145       hA[96]->Fill(V.id);
00146       hA[97]->Fill(V.funnyRank);
00147       hA[98]->Fill(V.z);
00148       hA[99]->Fill(T.prMuTrack->eta());
00149       hA[190+k]->Fill(T.prMuTrack->eta(),T.cluster.ET);
00150       
00151       hA[20]->Fill("goldW",1.);
00152       nGoldW++;
00153     
00154     }// loop over tracks
00155   }// loop over vertices
00156   if(nNoNear>0) hA[0]->Fill("noNear",1.);
00157   if(nNoAway>0) hA[0]->Fill("noAway",1.);
00158   if(nEta1>0) hA[0]->Fill("eta1",1.);
00159   if(nGoldW>0) hA[0]->Fill("goldW",1.);
00160 
00161 }
00162 
00163 
00164 //________________________________________________
00165 //________________________________________________
00166 void
00167 St2011WMaker::tag_Z_boson(){
00168 
00169   float par_jetPt=10.;
00170   float lowMass=70.; float highMass=140.;
00171   //select specific jet-type
00172   mJets = getJets(mJetTreeBranch); 
00173   if(mJetTreeChain) mJets = getJetsTreeAnalysis(mJetTreeBranch);
00174 
00175   //form invariant mass from lepton candidate and jet
00176   for(uint iv=0;iv<wEve->vertex.size();iv++) {//vertex loop
00177     WeveVertex &V=wEve->vertex[iv];
00178     for(uint it=0;it<V.eleTrack.size();it++) {// select track
00179       WeveEleTrack &T1=V.eleTrack[it];
00180       if(T1.isMatch2Cl==false) continue;
00181       assert(T1.cluster.nTower>0); // internal logical error
00182       assert(T1.nearTotET>0); // internal logical error
00183       if(T1.cluster.ET/T1.nearTotET< par_nearTotEtFrac) continue; // too large nearET
00184 
00185       //match lepton candidate with jet
00186       TLorentzVector jetVec;
00187       for (int i_jet=0; i_jet< nJets; i_jet++){//jet loop
00188         jetVec = *((StJet*)mJets->At(i_jet));
00189         if(jetVec.Pt()<par_jetPt) continue;//remove low pt jets
00190         
00191         //electron like cut on jets
00192         StJet* jet = getJet(i_jet);  float maxCluster=0.; 
00193         int totTowers=jet->nBtowers+jet->nEtowers;
00194         for(int itow=0;itow<totTowers;itow++){//loop over towers
00195           if(jet->tower(itow)->detectorId()==13)//drop endcap towers
00196             continue;
00197           
00198           int softId=jet->tower(itow)->towerId();
00199           //find highest 2x2 BTOW cluster in jet
00200           TVector3 pos=positionBtow[softId-1]; int iEta,iPhi;
00201           if( L2algoEtaPhi2IJ(pos.Eta(),pos.Phi(),iEta,iPhi)) 
00202             continue;
00203           float cluster=maxBtow2x2(iEta,iPhi,jet->zVertex).ET;
00204           if(cluster>maxCluster) maxCluster=cluster;
00205         }
00206 
00207         TVector3 jetVec3(jetVec.X(),jetVec.Y(),jetVec.Z());
00208         if(jetVec3.DeltaR(T1.primP)<par_nearDeltaR)
00209           continue;//skip jets in candidate phi isolation'cone'
00210         
00211         //form invM
00212         float e1=T1.cluster.energy;
00213         TVector3 p1=T1.primP; p1.SetMag(e1);
00214         TLorentzVector ele1(p1,e1); //lepton candidate 4- mom
00215         TLorentzVector sum=ele1+jetVec;
00216         float invM=sqrt(sum*sum);
00217         if(maxCluster/jet->jetPt < 0.5) continue;
00218         if(invM > lowMass && invM < highMass)
00219           wEve->zTag=true;
00220       }
00221     }
00222   }
00223   
00224 }
00225 
00226 //________________________________________________
00227 //________________________________________________
00228 void
00229 St2011WMaker::findPtBalance(){
00230 
00231   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00232     WeveVertex &V=wEve->vertex[iv];
00233     for(uint it=0;it<V.eleTrack.size();it++) {
00234       WeveEleTrack &T=V.eleTrack[it];
00235       if(T.isMatch2Cl==false) continue;
00236 
00237       //****loop over branch with EEMC****
00238       mJets = getJets(mJetTreeBranch); 
00239       if(mJetTreeChain) mJets = getJetsTreeAnalysis(mJetTreeBranch); 
00240       int nJetsWE=nJets;
00241       for (int i_jet=0; i_jet< nJetsWE; i_jet++){//loop over jets
00242         StJet* jet = getJet(i_jet);
00243         TVector3 jetVec; //vector for jet momentum
00244         jetVec.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
00245         if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
00246               T.ptBalance+=jetVec;
00247       }
00248       TVector3 clustPt(T.primP.X(),T.primP.Y(),0);
00249       clustPt.SetMag(T.cluster.ET);
00250       T.ptBalance+=clustPt;
00251       T.sPtBalance = T.ptBalance.Perp();
00252       if(T.ptBalance.Dot(clustPt)<0) T.sPtBalance *=-1.;
00253       
00254       //****loop over branch without EEMC****
00255       mJets = getJets(mJetTreeBranch_noEEMC); 
00256       if(mJetTreeChain) mJets = getJetsTreeAnalysis(mJetTreeBranch_noEEMC);
00257       int nJetsNE=nJets;
00258       for (int i_jet=0; i_jet< nJetsNE; i_jet++){//loop over jets
00259         StJet* jet = getJet(i_jet);
00260         TVector3 jetVec; //vector for jet momentum
00261         jetVec.SetPtEtaPhi(jet->Pt(),jet->Eta(),jet->Phi());
00262         if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
00263           T.ptBalance_noEEMC+=jetVec;
00264       }
00265       T.ptBalance_noEEMC+=clustPt;
00266       T.sPtBalance_noEEMC = T.ptBalance_noEEMC.Perp();
00267       if(T.ptBalance_noEEMC.Dot(clustPt)<0) T.sPtBalance_noEEMC *=-1.;
00268 
00269     }// end of loop over tracks
00270   }// end of loop over vertices
00271   
00272 }
00273 
00274 
00275 //________________________________________________
00276 //________________________________________________
00277 void
00278 St2011WMaker::findAwayJet(){
00279   // printf("\n******* find AwayJet() nVert=%d\n",wEve->vertex.size());
00280   //wEve->print();
00281   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00282     WeveVertex &V=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       
00287       // .... sum opposite in phi EMC components
00288       T.awayBtowET=sumBtowCone(V.z,-T.primP,1); // '1'= only cut on delta phi
00289       T.awayEmcET=T.awayBtowET;
00290       T.awayEtowET=sumEtowCone(V.z,-T.primP,1);
00291       T.awayEmcET+=T.awayEtowET; 
00292 
00293       //..... add TPC ET
00294       if(mMuDstMaker) T.awayTpcPT=sumTpcCone(V.id,-T.primP,1,T.pointTower.id);
00295       else T.awayTpcPT=sumTpcConeFromTree(V.id,-T.primP,1,T.pointTower.id);
00296       T.awayTotET=T.awayEmcET+T.awayTpcPT;
00297       T.awayTotET_noEEMC=T.awayBtowET+T.awayTpcPT;
00298       //printf("\n*** in   awayTpc=%.1f awayEmc=%.1f\n  ",T.awayTpcPT,T.awayEmcET); T.print(); 
00299     
00300     }// end of loop over tracks
00301   }// end of loop over vertices
00302 }
00303 
00304 //________________________________________________
00305 //________________________________________________
00306 void
00307 St2011WMaker::findNearJet(){
00308   //printf("\n******* findNearJet() nVert=%d\n",wEve->vertex.size());
00309 
00310   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00311     WeveVertex &V=wEve->vertex[iv];
00312     for(uint it=0;it<V.eleTrack.size();it++) {
00313       WeveEleTrack &T=V.eleTrack[it];
00314       if(T.isMatch2Cl==false) continue;
00315 
00316        // .... sum EMC-jet component
00317       T.nearBtowET=sumBtowCone(V.z,T.primP,2); // '2'=2D cone
00318       T.nearEmcET+=T.nearBtowET;
00319       T.nearEtowET=sumEtowCone(V.z,T.primP,2);
00320       T.nearEmcET+=T.nearEtowET; 
00321       // .... sum TPC-near component
00322       if(mMuDstMaker) T.nearTpcPT=sumTpcCone(V.id,T.primP,2,T.pointTower.id); // '2'=2D cone
00323       else T.nearTpcPT=sumTpcConeFromTree(V.id,T.primP,2,T.pointTower.id);
00324       float nearSum=T.nearEmcET+T.nearTpcPT;
00325       
00326       //fill histos separately for 2 types of events
00327       if(T.pointTower.id>0) { //only barrel towers
00328         /* correct for double counting of electron track in near cone rarely primTrPT<10 GeV & globPT>10 - handle this here */
00329         if(T.primP.Pt()>par_trackPt) nearSum-=par_trackPt; 
00330         else  nearSum-=T.primP.Pt();
00331         T.nearTotET=nearSum;
00332         T.nearTotET_noEEMC=nearSum-T.nearEtowET; 
00333         float nearTotETfrac=T.cluster.ET/ T.nearTotET;
00334 
00335         hA[40]->Fill(T.nearEmcET);
00336         hA[41]->Fill(T.cluster.ET,T.nearEmcET-T.cluster.ET);
00337         hA[42]->Fill(nearTotETfrac);
00338         hA[47]->Fill(T.nearTpcPT);
00339         hA[48]->Fill(T.nearEmcET,T.nearTpcPT);
00340         hA[49]->Fill(nearSum);
00341 
00342         // check east/west yield diff
00343         hA[210]->Fill(T.cluster.position.PseudoRapidity(),T.nearEtowET);
00344         if(T.cluster.position.PseudoRapidity()>0) hA[211]->Fill(T.cluster.position.Phi(),T.nearEtowET);
00345         else hA[212]->Fill(T.cluster.position.Phi(),T.nearEtowET);
00346 
00347       }
00348       else if(T.pointTower.id<0) { //only endcap towers
00349         /* correct for double counting of electron track in near cone rarely primTrPT<10 GeV & globPT>10 - handle this here */
00350         if(T.primP.Pt()>parE_trackPt) nearSum-=parE_trackPt; 
00351         else  nearSum-=T.primP.Pt();
00352         T.nearTotET=nearSum;
00353         T.nearTotET_noEEMC=nearSum-T.nearEtowET;
00354         float nearTotETfrac=T.cluster.ET/ T.nearTotET;
00355 
00356         hE[40]->Fill(T.nearEmcET);
00357         hE[41]->Fill(T.cluster.ET,T.nearEmcET-T.cluster.ET);
00358         hE[42]->Fill(nearTotETfrac);
00359         hE[47]->Fill(T.nearTpcPT);
00360         hE[48]->Fill(T.nearEmcET,T.nearTpcPT);
00361         hE[49]->Fill(nearSum);
00362       }  
00363 
00364     } // end of loop over tracks
00365   }// end of loop over vertices
00366  
00367 
00368 }
00369 
00370 //________________________________________________
00371 //________________________________________________
00372 float
00373 St2011WMaker::sumBtowCone( float zVert,  TVector3 refAxis, int flag){
00374   /* flag=1 : only delta phi cut;  flag=2 use 2D cut */
00375   assert(flag==1 || flag==2);
00376   double ptSum=0;
00377 
00378   //.... process BTOW hits
00379   for(int i=0;i< mxBtow;i++) {
00380     float ene=wEve->bemc.eneTile[kBTow][i];
00381     if(ene<=0) continue;
00382     TVector3 primP=positionBtow[i]-TVector3(0,0,zVert);
00383     primP.SetMag(ene); // it is 3D momentum in the event ref frame
00384     if(flag==1) {
00385       float deltaPhi=refAxis.DeltaPhi(primP);
00386       if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
00387     }
00388     if(flag==2) {
00389       float deltaR=refAxis.DeltaR(primP);
00390       if(deltaR> par_nearDeltaR) continue;
00391     }
00392     ptSum+=primP.Perp();
00393   }
00394 
00395   return ptSum;
00396 }
00397 
00398 //________________________________________________
00399 //________________________________________________
00400 float
00401 St2011WMaker::sumTpcConeFromTree(int vertID, TVector3 refAxis, int flag,int pointTowId){
00402 
00403   // flag=2 use 2D cut, 1= only delta phi
00404 
00405   assert(vertID>=0);
00406   assert(vertID<(int)wEve->vertex.size());
00407   
00408   double ptSum=0;
00409   WeveVertex &V=wEve->vertex[vertID];
00410   for(uint it=0;it<V.prTrList.size();it++){
00411     StMuTrack *prTr=V.prTrList[it];
00412     if(prTr->flag()<=0) continue;
00413     if(prTr->flag()!=301 && pointTowId>0) continue;// TPC-only regular tracks for barrel candidate
00414     if(prTr->flag()!=301 && prTr->flag()!=311 && pointTowId<0) continue;// TPC regular and short EEMC tracks for endcap candidate
00415     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00416     if(hitFrac<par_nHitFrac) continue;
00417     StThreeVectorF prPvect=prTr->p();
00418     TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
00419     // printf(" prTrID=%4d  prTrEta=%.3f prTrPhi/deg=%.1f prPT=%.1f  nFitPts=%d\n", prTr->id(),prTr->eta(),prTr->phi()/3.1416*180.,prTr->pt(),prTr->nHitsFit());
00420     if(flag==1) {
00421       float deltaPhi=refAxis.DeltaPhi(primP);
00422       if(fabs(deltaPhi)> par_awayDeltaPhi) continue;
00423     }
00424     if(flag==2) {
00425       float deltaR=refAxis.DeltaR(primP);
00426       if(deltaR>par_nearDeltaR) continue;
00427 
00428     }
00429     float pT=prTr->pt();
00430     //    printf(" passed pt=%.1f\n",pT);
00431     
00432     //separate quench for barrel and endcap candidates
00433     if(pT>par_trackPt && pointTowId>0) ptSum+=par_trackPt;
00434     else if(pT>parE_trackPt && pointTowId<0) ptSum+=parE_trackPt;
00435     else  ptSum+=pT;
00436   }
00437   return ptSum;
00438 }
00439 
00440 
00441 // ************* Barrel Code ************ //
00442 // ************************************** //
00443 
00444 //________________________________________________
00445 //________________________________________________
00446 int
00447 St2011WMaker::extendTrack2Barrel(){// return # of extended tracks
00448   //printf("******* extendTracks() nVert=%d\n",wEve->vertex.size());
00449   if(!wEve->l2bitET) return 0; //fire barrel trigger
00450 
00451   int nTrB=0;
00452   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00453     WeveVertex &V=wEve->vertex[iv];
00454     if(V.rank<0) continue; //remove vertex for endcap algo only
00455     for(uint it=0;it<V.eleTrack.size();it++) {
00456       WeveEleTrack &T=V.eleTrack[it];
00457       if(T.prMuTrack->flag()!=301) continue; //remove track for endcap algo only
00458 
00459       //.... extrapolate track to the barrel @ R=entrance....
00460       const StPhysicalHelixD TrkHlx=T.prMuTrack->outerHelix();
00461       float Rcylinder= mBtowGeom->Radius();
00462       pairD  d2;
00463       d2 = TrkHlx.pathLength(Rcylinder);
00464       //printf(" R=%.1f path 1=%f, 2=%f, period=%f, R=%f\n",Rctb,d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
00465 
00466       // assert(d2.first<0); // propagate backwards
00467       // assert(d2.second>0); // propagate forwards
00468       if(d2.first>=0 || d2.second<=0) {
00469         LOG_WARN<< Form("MatchTrk , unexpected solution for track crossing CTB\n  d2.firts=%f, second=%f, swap them",  d2.first, d2.second)<<endm;
00470         float xx=d2.first;
00471         d2.first=d2.second;
00472         d2.second=xx;
00473       }
00474 
00475       // extrapolate track to cylinder
00476       StThreeVectorD posR = TrkHlx.at(d2.second);
00477       //printf(" punch2 x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f\n",posCTB.x(),posCTB.y(),posCTB.z(),xmagn);
00478       float etaF=posR.pseudoRapidity();
00479       float phiF=posR.phi();
00480       int iEta,iPhi;
00481       if( L2algoEtaPhi2IJ(etaF, phiF,iEta,iPhi)) continue;
00482       nTrB++;
00483        hA[20]->Fill("@B",1.);
00484       //printf(" phi=%.0f deg,  eta=%.2f, iEta=%d, iPhi=%d\n",posCTB.phi()/3.1416*180.,posCTB. pseudoRapidity(),iEta, iPhi);
00485       int twID= mapBtowIJ2ID[ iEta+ iPhi*mxBTetaBin];
00486       // printf("hit Tower ID=%d\n",twID);
00487 
00488       T.pointTower.id=twID;
00489       T.pointTower.R=TVector3(posR.x(),posR.y(),posR.z());
00490       T.pointTower.iEta=iEta;
00491       T.pointTower.iPhi=iPhi;
00492       //T.print();
00493     }
00494   }// end of loop over vertices
00495 
00496   if(nTrB<=0) return -1;
00497   hA[0]->Fill("TrB",1.0);
00498   return 0;
00499  
00500 }
00501 
00502 //________________________________________________
00503 //________________________________________________
00504 int
00505 St2011WMaker::matchTrack2BtowCluster(){
00506   // printf("******* matchCluster() nVert=%d\n",wEve->vertex.size());
00507   int nTr=0;
00508   float Rcylinder= mBtowGeom->Radius();
00509   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00510     WeveVertex &V=wEve->vertex[iv];
00511     float zVert=V.z;
00512     for(uint it=0;it<V.eleTrack.size();it++) {
00513       WeveEleTrack &T=V.eleTrack[it];
00514       if(T.pointTower.id<=0) continue; //skip endcap towers
00515       
00516       float trackPT=T.prMuTrack->momentum().perp();
00517       T.cluster=maxBtow2x2( T.pointTower.iEta, T.pointTower.iPhi,zVert);      
00518       hA[33]->Fill( T.cluster.ET);
00519       hA[34]->Fill(T.cluster.adcSum,trackPT);
00520       hA[110]->Fill( T.cluster.ET);
00521 
00522       // ........compute surroinding cluster energy
00523       int iEta=T.cluster.iEta;
00524       int iPhi=T.cluster.iPhi;
00525       T.cl4x4=sumBtowPatch(iEta-1,iPhi-1,4,4,zVert); // needed for lumi monitor
00526     
00527       if (T.cluster.ET <par_clustET) continue; // too low energy
00528       hA[20]->Fill("CL",1.);
00529 
00530       hA[206]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
00531 
00532       hA[37]->Fill( T.cl4x4.ET);
00533       hA[38]->Fill(T.cluster.energy, T.cl4x4.energy-T.cluster.energy);
00534       
00535       float frac24=T.cluster.ET/(T.cl4x4.ET);
00536       hA[39]->Fill(frac24);
00537       if(frac24<par_clustFrac24) continue;
00538   
00539       hA[20]->Fill("fr24",1.);
00540  
00541       //.. spacial separation (track - cluster)
00542       TVector3 D=T.pointTower.R-T.cluster.position;
00543 
00544       hA[43]->Fill( T.cluster.energy,D.Mag());
00545       hA[44]->Fill( T.cluster.position.z(),D.z());
00546       float delPhi=T.pointTower.R.DeltaPhi(T.cluster.position);
00547       //   printf("aaa %f %f %f   phi=%f\n",D.x(),D.y(),D.z(),delPhi);
00548       hA[45]->Fill( T.cluster.energy,Rcylinder*delPhi);// wrong?
00549       hA[46]->Fill( D.Mag());
00550       hA[199]->Fill(T.cluster.position.PseudoRapidity(),D.Mag());
00551       hA[207]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
00552 
00553       if(D.Mag()>par_delR3D) continue; 
00554       T.isMatch2Cl=true; // cluster is matched to TPC track
00555       hA[20]->Fill("#Delta R",1.);
00556       hA[111]->Fill( T.cluster.ET);
00557       
00558       hA[208]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.ET);
00559 
00560       nTr++;
00561     }// end of one vertex
00562   }// end of vertex loop
00563 
00564   if(nTr<=0) return -1; 
00565   hA[0]->Fill("Tr2Cl",1.0);
00566   return 0;
00567 
00568 }
00569 
00570 //________________________________________________
00571 //________________________________________________ 
00572 WeveCluster 
00573 St2011WMaker::maxBtow2x2(int iEta, int iPhi, float zVert){
00574   //printf("   maxBtow2x2  seed iEta=%d iPhi=%d \n",iEta, iPhi);
00575   const int L=2; // size of the summed square 
00576 
00577   WeveCluster maxCL;
00578   // just 4 cases of 2x2 clusters
00579   float maxET=0;
00580   int I0=iEta-1;
00581   int J0=iPhi-1;
00582   for(int I=I0;I<=I0+1;I++){
00583     for(int J=J0;J<=J0+1;J++) {
00584       WeveCluster CL=sumBtowPatch(I,J,L,L,zVert);
00585       if(maxET>CL.ET) continue;
00586       maxET=CL.ET;
00587       maxCL=CL;
00588       // printf("   newMaxETSum=%.1f iEta=%d iPhi=%d \n",maxET, I,J);
00589     }
00590   }// 4 combinations done
00591   //printf(" final inpEve=%d SumET2x2=%.1f \n",nInpEve,maxET);
00592   return maxCL;  
00593 }
00594 
00595 
00596 //________________________________________________
00597 //________________________________________________
00598 WeveCluster
00599 St2011WMaker::sumBtowPatch(int iEta, int iPhi, int Leta,int  Lphi, float zVert){
00600   //printf("  eveID=%d btowSquare seed iEta=%d[+%d] iPhi=%d[+%d] zVert=%.0f \n",wEve->id,iEta,Leta, iPhi,Lphi,zVert);
00601   WeveCluster CL; // object is small, not to much overhead in creating it
00602   CL.iEta=iEta;
00603   CL.iPhi=iPhi;
00604   TVector3 R; 
00605   double sumW=0;
00606   float Rcylinder= mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
00607   for(int i=iEta; i<iEta+Leta;i++){// trim in eta-direction
00608     if(i<0) continue;
00609     if(i>=mxBTetaBin) continue;
00610     for(int j=iPhi;j<iPhi+Lphi;j++) {// wrap up in the phi-direction
00611       int jj=(j+mxBTphiBin)%mxBTphiBin;// keep it always positive
00612       //if(L<5) printf("n=%2d  i=%d jj=%d\n",CL.nTower,i,jj);
00613       int softID= mapBtowIJ2ID[ i+ jj*mxBTetaBin];
00614       float ene= wEve->bemc.eneTile[kBTow][softID-1];
00615       if(ene<=0) continue; // skip towers w/o energy
00616       float adc= wEve->bemc.adcTile[kBTow][softID-1];
00617       float delZ=positionBtow[softID-1].z()-zVert;
00618       float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
00619       float ET=ene*e2et;
00620       float logET=log10(ET+0.5);
00621       CL.nTower++;
00622       CL.energy+=ene;
00623       CL.ET+=ET;
00624       CL.adcSum+=adc;
00625       if(logET>0) {
00626         R+=logET*positionBtow[softID-1];
00627         sumW+=logET;
00628       }
00629       //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); 
00630     }
00631     // printf(" end btowSquare: iEta=%d  nTw=%d, ET=%.1f adc=%.1f\n",i,CL.nTower,CL.ET,CL.adcSum);
00632     if(sumW>0) {
00633       CL.position=1./sumW*R; // weighted cluster position
00634     } else {
00635       CL.position=TVector3(0,0,999);
00636     } 
00637   }
00638   return CL;
00639 }
00640 
00641 
00642 // $Log: St2011W_algo.cxx,v $
00643 // Revision 1.3  2011/02/25 06:03:47  stevens4
00644 // addes some histos and enabled running on MC
00645 //
00646 // Revision 1.2  2011/02/17 04:16:22  stevens4
00647 // move sector dependent track QA cuts before track pt>10 cut and lower par_clustET and par_ptBalance thresholds to 14 GeV
00648 //
00649 // Revision 1.1  2011/02/10 20:33:23  balewski
00650 // start
00651 //
StSpinPool/StWalgo2011/St2011W_algo.cxx
