StRoot: StSpinPool/StWalgo2011/St2011ZMaker.cxx Source File

00001 // $Id: St2011ZMaker.cxx,v 1.1 2011/02/10 20:33:24 balewski Exp $
00002 //
00003 //*-- Author : Ross Corliss, MIT
00004 //  changes Jan Balewski, MIT
00005 // 
00006 
00007 #include "St2011WMaker.h"
00008 #include "WeventDisplay.h"
00009 #include "St2011ZMaker.h"
00010 
00011 ClassImp(St2011ZMaker)
00012 const float PI=TMath::Pi();
00013 
00014 //_____________________________________________________________________________
00015 //
00016 St2011ZMaker::St2011ZMaker(const char *name):StMaker(name){
00017   wMK=0;muMK=0;HList=0;
00018 
00019 }
00020 
00021 
00022 //_____________________________________________________________________________
00023 //
00024 Int_t St2011ZMaker::Init(){
00025   assert(wMK);
00026   assert(HList);
00027   initHistos();
00028   return StMaker::Init();
00029 }
00030 
00031 
00032 //_____________________________________________________________________________
00033 //
00034 Int_t St2011ZMaker::InitRun  (int runumber){
00035   LOG_INFO<<Form("::InitRun(%d) done, Z-algo params: nearTotEtFrac=%.2f,  clusterEt=%.1f GeV, delPhi12>%.2f rad, Zmass in[%.1f,%.1f]\n",
00036                  runumber, par_nearTotEtFracZ,par_clusterEtZ,par_delPhi12,par_minMassZ,par_maxMassZ)<<endm;
00037   return 0;
00038 }
00039 
00040 //_____________________________________________________________________________
00041 //
00042 Int_t St2011ZMaker::FinishRun  (int runnumber){
00043 return 0;
00044 }
00045 
00046 //_____________________________________________________________________________
00047 //
00048 Int_t 
00049 St2011ZMaker::Make(){
00050 
00051   //fill various histograms and arrays with event data
00052   find_Z_boson();
00053   findEndcap_Z_boson();
00054 
00055   return kStOK;
00056 }
00057 
00058 //============================
00059 void
00060 St2011ZMaker::printJan(WeveEleTrack *T) {
00061   int ibp=kBTow;
00062   WevePointTower poiTw=T->pointTower;
00063   WeveCluster cl=T->cluster;
00064   int id= poiTw.id;
00065   float adc= wMK->wEve->bemc.adcTile[ibp][id-1];
00066   float frac= adc/4096*60 /cl.ET;
00067   printf("Ztower Q=%d pointTw: id=%d ADC=%.0f  2x2ET=%.1f frac=%.2f\n",T->prMuTrack->charge(),id,adc,cl.ET,frac);
00068 }
00069 
00070 
00071 //_____________________________________________________________________________
00072 //
00073 void
00074 St2011ZMaker::findEndcap_Z_boson(){
00075   Wevent2011 *wEve=wMK->wEve;
00076   // printf("========= findEndcap_Z_boson() \n");
00077   
00078   hA[50]->Fill("inp",1.);
00079 
00080   // search for  Zs ............
00081   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00082     hA[50]->Fill("vert",1.);
00083     WeveVertex &V=wEve->vertex[iv];
00084     for(uint it=0;it<V.eleTrack.size();it++) {
00085       WeveEleTrack &T=V.eleTrack[it];
00086       if(T.isMatch2Cl==false) continue;
00087       assert(T.cluster.nTower>0); // internal logical error
00088       assert(T.nearTotET>0); // internal logical error
00089    
00090       //place cuts on reco track first (both barrel and endcap tracks allowed)
00091       float isoET1=T.cluster.ET /T.cl4x4.ET;
00092       hA[51]->Fill(isoET1);
00093       hA[52]->Fill(T.cluster.ET);
00094       hA[50]->Fill("tr1",1.);
00095       if(T.cluster.ET<par_clusterEtZ) continue;
00096       hA[50]->Fill("et1",1.);
00097 
00098       float fracET1=T.cluster.ET /T.nearTotET;
00099       hA[53]->Fill(fracET1);
00100       if(fracET1<par_nearTotEtFracZ) continue; 
00101       hA[50]->Fill("con1",1.);
00102 
00103       //find highest ET endcap cluster
00104       float maxET=0; 
00105       WeveCluster maxCluster;
00106       for(int iEta=0; iEta<12; iEta++){ //loop over eta bins
00107         for(int iPhi=0; iPhi<60; iPhi++){ //loop over phi bins
00108           
00109           WeveCluster eclust=wMK->maxEtow2x1(iEta,iPhi,V.z);
00110           if(eclust.ET < par_clusterEtZ) continue;
00111           if(maxET > eclust.ET) continue;
00112           else {
00113             maxET=eclust.ET;
00114             maxCluster=eclust;
00115           }
00116         
00117         }
00118       }
00119       if(maxCluster.ET <= 1.0) continue;//remove low E clusters
00120 
00121       //apply cuts to ETOW cluster
00122       WeveCluster cl4x4=wMK->sumEtowPatch(maxCluster.iEta-1,maxCluster.iPhi-1,4,4,V.z);
00123       hA[54]->Fill(maxCluster.ET/cl4x4.ET);
00124       if(maxCluster.ET/cl4x4.ET<wMK->parE_clustFrac24) continue;
00125       hA[55]->Fill(maxCluster.ET);
00126       hA[50]->Fill("tr2",1.);
00127       if(maxCluster.ET < par_clusterEtZ) continue;
00128       hA[50]->Fill("et2",1.);
00129       
00130       //assume poor tracking effic so only towers in nearCone
00131       float nearBtow=wMK->sumBtowCone(V.z,maxCluster.position,2);
00132       float nearEtow=wMK->sumEtowCone(V.z,maxCluster.position,2);
00133       float nearSum=nearBtow; nearSum+=nearEtow;
00134       hA[56]->Fill(maxCluster.ET/nearSum);
00135       if(maxCluster.ET/nearSum<wMK->parE_nearTotEtFrac) continue;
00136       hA[50]->Fill("con2",1.);
00137       
00138       //add plots of good candidates
00139       float e1=T.cluster.energy;
00140       float e2=maxCluster.energy;
00141       TVector3 p1=T.primP; p1.SetMag(e1);
00142       TVector3 p2=maxCluster.position; p2.SetMag(e2);
00143       
00144       float del_phi=p1.DeltaPhi(p2);
00145       float xx=del_phi;
00146       if(xx<-PI+1) xx+=2*PI;
00147       hA[57]->Fill(xx);
00148       if(fabs(del_phi)<par_delPhi12) continue;
00149       hA[50]->Fill("phi12",1.);
00150       
00151       TVector3 psum=p1+p2;
00152       float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
00153       if(mass2<1.) continue; 
00154       hA[50]->Fill("m2",1.);
00155       
00156       float mass=sqrt(mass2);
00157       hA[58]->Fill(mass);
00158 
00159     } //track loop
00160   } //vertex loop
00161 }
00162 
00163 
00164 //_____________________________________________________________________________
00165 //
00166 void
00167 St2011ZMaker::find_Z_boson(){
00168   Wevent2011 *wEve=wMK->wEve;
00169   // printf("========= find_Z_boson() \n");
00170   
00171   hA[31]->Fill(wEve->vertex.size());
00172   hA[0]->Fill("inp",1.);
00173 
00174   // search for  Zs ............
00175   for(uint iv=0;iv<wEve->vertex.size();iv++) {
00176     hA[0]->Fill("vert",1.);
00177     WeveVertex &V=wEve->vertex[iv];
00178     hA[32]->Fill(V.eleTrack.size());
00179     if(V.eleTrack.size()<2) continue;
00180     hA[0]->Fill("TT",1.); // at least 2 isolated tracks exist
00181 
00182     //only one Z can come from a vertex, and it should be the highest-energy object
00183     //hence, the two highest-et clusters should correspond to the z.  Pick those 
00184     //eventually, but for now, just try all of them.
00185     for(uint it=0;it<V.eleTrack.size()-1;it++) { //.....select first track:
00186       WeveEleTrack &T1=V.eleTrack[it];
00187       if(T1.isMatch2Cl==false) continue;
00188       assert(T1.cluster.nTower>0); // internal logical error
00189       assert(T1.nearTotET>0); // internal logical error
00190 
00191       float isoET1=T1.cluster.ET /T1.cl4x4.ET;
00192       hA[29]->Fill(isoET1);
00193 
00194       hA[23]->Fill(T1.cluster.ET);
00195       hA[0]->Fill("tr1",1.);
00196       if(T1.cluster.ET<par_clusterEtZ) continue;
00197       hA[0]->Fill("et1",1.);
00198 
00199       float fracET1=T1.cluster.ET /T1.nearTotET;
00200       hA[24]->Fill(fracET1);
00201       if(fracET1< par_nearTotEtFracZ) continue; 
00202       hA[0]->Fill("con1",1.);
00203 
00204       for (uint it2=it+1;it2<V.eleTrack.size();it2++) { //.....select second track:
00205         WeveEleTrack &T2=V.eleTrack[it2];
00206         if(T2.isMatch2Cl==false) continue;
00207         assert(T2.cluster.nTower>0); // internal logical error
00208         assert(T2.nearTotET>0); // internal logical error
00209 
00210         float isoET2=T2.cluster.ET /T2.cl4x4.ET;
00211         hA[30]->Fill(isoET2);
00212 
00213         hA[25]->Fill(T2.cluster.ET);
00214         hA[0]->Fill("tr2",1.);
00215         if(T2.cluster.ET<par_clusterEtZ) continue;
00216         hA[0]->Fill("et2",1.);
00217 
00218         float fracET2=T2.cluster.ET /T2.nearTotET;
00219         hA[26]->Fill(fracET2);
00220         if(fracET2< par_nearTotEtFracZ) continue; 
00221         hA[0]->Fill("con2",1.);
00222 
00223         float e1=T1.cluster.energy;
00224         float e2=T2.cluster.energy;
00225         TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
00226         TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
00227 
00228         float del_phi=p1.DeltaPhi(p2);
00229         //printf("del Phi=%f\n",del_phi);
00230         float xx=del_phi;
00231         if(xx<-PI+1) xx+=2*PI;
00232         hA[27]->Fill(xx);
00233         if(fabs(del_phi)<par_delPhi12) continue;
00234         hA[0]->Fill("phi12",1.);
00235 
00236         TVector3 psum=p1+p2;
00237         float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
00238         if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
00239         hA[0]->Fill("m2",1.);
00240 
00241         float mass=sqrt(mass2);
00242         int Q1Q2=T1.prMuTrack->charge()*T2.prMuTrack->charge();
00243         if (Q1Q2==1) { //..  same sign , can't be Z-> e+ e-
00244           hA[14]->Fill(mass);
00245           continue;
00246         }
00247 
00248         //..... now only opposite sign
00249         hA[0]->Fill("QQ",1.);
00250         hA[15]->Fill(mass);
00251         hA[33]->Fill(T1.cluster.ET,T1.prMuTrack->charge()/T1.prMuTrack->pt()); 
00252         hA[33]->Fill(T2.cluster.ET,T2.prMuTrack->charge()/T2.prMuTrack->pt()); 
00253         hA[34]->Fill(T1.pointTower.iEta ,T1.cluster.energy);
00254         hA[34]->Fill(T2.pointTower.iEta ,T2.cluster.energy);
00255 #if 0
00256         printf("RCC:  Found Z w/ invmass=%f\n",mass);
00257         printJan(&T1);
00258         printJan(&T2);
00259 
00260 
00261         if (!wMK->isMC || (wMK->isMC&& wEve->id<500) )
00262           { printf("\n ZZZZZZZZZZZZZZZZZZZ\n");
00263             if(mass<par_minMassZ) 
00264             wMK->wDisaply->exportEvent("Zlow",V,T1);
00265             else
00266               wMK->wDisaply->exportEvent("Zgood",V,T1);
00267             printf("RCC:  Found Z w/ invmass=%f\n",mass);
00268             wEve->print();
00269           }
00270         
00271 #endif
00272 
00273         if (mass<par_minMassZ) continue; //enforce a lower bound
00274         hA[0]->Fill("Zlow",1.);
00275 
00276         if (mass>par_maxMassZ) continue; //enforce an upper bound
00277         hA[0]->Fill("Zhigh",1.);
00278 
00279         // **** I stoped changes here, Jan 
00280 
00281         float fmax1=T1.cluster.ET/T1.cl4x4.ET;
00282         float fmax2=T2.cluster.ET/T2.cl4x4.ET;
00283 
00284         hA[21]->Fill(fmax1,fmax2);
00285         hA[22]->Fill(T1.cluster.ET,T2.cluster.ET);
00286 
00287         hA[1]->Fill(mass);
00288         hA[2]->Fill(T1.prMuTrack->charge(),T2.prMuTrack->charge());
00289         hA[3]->Fill(T1.prMuTrack->charge()*T2.prMuTrack->charge());
00290         hA[4]->Fill(p1.Phi(),p2.Phi());
00291         hA[5]->Fill(del_phi);
00292         hA[6]->Fill(mass,T1.prMuTrack->charge()/T1.primP.Perp()*T2.prMuTrack->charge()/T1.primP.Perp());
00293         hA[7]->Fill(mass,T1.prMuTrack->charge()*T2.prMuTrack->charge());
00294         hA[8]->Fill(T1.cluster.ET);
00295         if (T1.prMuTrack->charge()>0)
00296           {
00297             hA[9]->Fill(p1.Eta(),p1.Phi());
00298             hA[10]->Fill(p2.Eta(),p2.Phi());
00299           }
00300         else
00301           {
00302             hA[10]->Fill(p1.Eta(),p1.Phi());
00303             hA[9]->Fill(p2.Eta(),p2.Phi());
00304           }
00305 
00306         hA[11]->Fill(fmax1,fmax2);
00307         hA[12]->Fill(T1.cluster.ET,T2.cluster.ET);
00308         hA[13]->Fill(mass,del_phi);
00309 
00310       }
00311 
00312 
00313     
00314     }// loop over first track
00315   }// loop over vertices
00316  
00317 }
00318 
00319 
00320 // $Log: St2011ZMaker.cxx,v $
00321 // Revision 1.1  2011/02/10 20:33:24  balewski
00322 // start
00323 //