StRoot: StSpinPool/StWalgoB2009/St2009ZMaker.cxx Source File

00001 // $Id: St2009ZMaker.cxx,v 1.11 2011/09/14 14:23:21 stevens4 Exp $
00002 //
00003 //*-- Author : Ross Corliss, MIT
00004 //  changes Jan Balewski, MIT
00005 // 
00006 
00007 #include "St2009WMaker.h"
00008 #include "WeventDisplay.h"
00009 #include "St2009ZMaker.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 ClassImp(St2009ZMaker)
00017 
00018 //_____________________________________________________________________________
00019 //
00020 St2009ZMaker::St2009ZMaker(const char *name):StMaker(name){
00021   wMK=0;muMK=0;HList=0;
00022 
00023 }
00024 
00025 
00026 //_____________________________________________________________________________
00027 //
00028 Int_t St2009ZMaker::Init(){
00029   assert(wMK);
00030   assert(HList);
00031   initHistos();
00032   return StMaker::Init();
00033 }
00034 
00035 
00036 //_____________________________________________________________________________
00037 //
00038 Int_t St2009ZMaker::InitRun  (int runumber){
00039   LOG_INFO<<Form("::InitRun(%d) done, Z-algo params: nearTotEtFrac=%.2f,  clusterEt=%.1f GeV, delPhi12>%.2f rad, Zmass in[%.1f,%.1f]\n",
00040                  runumber, par_nearTotEtFracZ,par_clusterEtZ,par_delPhi12,par_minMassZ,par_maxMassZ)<<endm;
00041   return 0;
00042 }
00043 
00044 //_____________________________________________________________________________
00045 //
00046 Int_t St2009ZMaker::FinishRun  (int runnumber){
00047 return 0;
00048 }
00049 
00050 //_____________________________________________________________________________
00051 //
00052 Int_t 
00053 St2009ZMaker::Make(){
00054 
00055   //fill various histograms and arrays with event data
00056   find_Z_boson();
00057 
00058   return kStOK;
00059 }
00060 
00061 //============================
00062 void
00063 St2009ZMaker::printJan(WeveEleTrack *T) {
00064   int ibp=kBTow;
00065   WevePointTower poiTw=T->pointTower;
00066   WeveCluster cl=T->cluster;
00067   int id= poiTw.id;
00068   float adc= wMK->wEve.bemc.adcTile[ibp][id-1];
00069   float frac= adc/4096*60 /cl.ET;
00070   printf("Ztower Q=%d pointTw: id=%d ADC=%.0f  2x2ET=%.1f frac=%.2f\n",T->prMuTrack->charge(),id,adc,cl.ET,frac);
00071 }
00072 
00073 
00074 
00075 //_____________________________________________________________________________
00076 //
00077 void
00078 St2009ZMaker::find_Z_boson(){
00079   const float PI=TMath::Pi();
00080   Wevent2009 &wEve=wMK->wEve;
00081   float zdcRate=wMK->mMuDstMaker->muDst()->event()->runInfo().zdcCoincidenceRate();
00082   // printf("========= find_Z_boson() \n");
00083   
00084   hA[31]->Fill(wEve.vertex.size());
00085   hA[0]->Fill("inp",1.);
00086 
00087   // search for  Zs ............
00088   for(uint iv=0;iv<wEve.vertex.size();iv++) {
00089     hA[0]->Fill("vert",1.);
00090     WeveVertex &V=wEve.vertex[iv];
00091     hA[32]->Fill(V.eleTrack.size());
00092     if(V.eleTrack.size()<2) continue;
00093     hA[0]->Fill("TT",1.); // at least 2 tracks exist
00094 
00095     //fill stack plots before doing full Z algo below
00096     for(uint itStack=0;itStack<V.eleTrack.size()-1;itStack++) {
00097       WeveEleTrack &stackT1=V.eleTrack[itStack]; //get first track
00098       
00099       //make plot of mached clusters
00100       if(stackT1.pointTower.id==0) continue;
00101       for (uint itStack2=itStack+1;itStack2<V.eleTrack.size();itStack2++) {
00102         WeveEleTrack &stackT2=V.eleTrack[itStack2]; //get second track
00103         if(stackT2.pointTower.id==0) continue;
00104         // invariant mass of "matched clusters"
00105         if(stackT1.prMuTrack->charge()*stackT2.prMuTrack->charge()<0) hA[40]->Fill(calcMass(stackT1,stackT2)); //opposite-sign 
00106         else hA[50]->Fill(calcMass(stackT1,stackT2));
00107       }
00108     
00109       //make plot of clusters passing deltaR cut
00110       if(stackT1.cluster.ET<par_clusterEtZ) continue;
00111       TVector3 stackD1=stackT1.pointTower.R-stackT1.cluster.position;
00112       if(stackD1.Mag()>par_delR3DZ) continue;
00113       for (uint itStack2=itStack+1;itStack2<V.eleTrack.size();itStack2++) {
00114         WeveEleTrack &stackT2=V.eleTrack[itStack2]; //get second track
00115         if(stackT2.cluster.ET<par_clusterEtZ) continue;
00116         TVector3 stackD2=stackT2.pointTower.R-stackT2.cluster.position;
00117         if(stackD2.Mag()>par_delR3DZ) continue;
00118         // invariant mass of clusters passing deltaR cut
00119         if(stackT1.prMuTrack->charge()*stackT2.prMuTrack->charge()<0) hA[41]->Fill(calcMass(stackT1,stackT2)); //opposite-sign 
00120         else hA[51]->Fill(calcMass(stackT1,stackT2));
00121       }
00122 
00123       //make plot of clusters passing 2x2/4x4 cut
00124       if(stackT1.cluster.ET/stackT1.cl4x4.ET < par_clustFrac24Z) continue;
00125       for (uint itStack2=itStack+1;itStack2<V.eleTrack.size();itStack2++) {
00126         WeveEleTrack &stackT2=V.eleTrack[itStack2]; //get second track
00127         if(stackT2.cluster.ET<par_clusterEtZ) continue;
00128         TVector3 stackD2=stackT2.pointTower.R-stackT2.cluster.position;
00129         if(stackD2.Mag()>par_delR3DZ) continue;
00130         if(stackT2.cluster.ET/stackT2.cl4x4.ET < par_clustFrac24Z) continue;
00131         // invariant mass of clusters passing 2x2/4x4 cut
00132         if(stackT1.prMuTrack->charge()*stackT2.prMuTrack->charge()<0) hA[42]->Fill(calcMass(stackT1,stackT2)); //opposite-sign 
00133         else hA[52]->Fill(calcMass(stackT1,stackT2));
00134       }
00135 
00136       //make plot of clusters passing near cone cut
00137       if(stackT1.cluster.ET/stackT1.nearTotET < par_nearTotEtFracZ) continue;
00138       for (uint itStack2=itStack+1;itStack2<V.eleTrack.size();itStack2++) {
00139         WeveEleTrack &stackT2=V.eleTrack[itStack2]; //get second track
00140         if(stackT2.cluster.ET<par_clusterEtZ) continue;
00141         if(stackT2.cluster.ET/stackT2.cl4x4.ET < par_clustFrac24Z) continue;
00142         TVector3 stackD2=stackT2.pointTower.R-stackT2.cluster.position;
00143         if(stackD2.Mag()>par_delR3DZ) continue;
00144         if(stackT2.cluster.ET/stackT2.nearTotET < par_nearTotEtFracZ) continue;
00145         // invariant mass of clusters passing near cone cut
00146         if(stackT1.prMuTrack->charge()*stackT2.prMuTrack->charge()<0) hA[43]->Fill(calcMass(stackT1,stackT2)); //opposite-sign 
00147         else hA[53]->Fill(calcMass(stackT1,stackT2));
00148       }
00149 
00150     }
00151 
00152     //only one Z can come from a vertex, and it should be the highest-energy object
00153     //hence, the two highest-et clusters should correspond to the z.  Pick those 
00154     //eventually, but for now, just try all of them.
00155     for(uint it=0;it<V.eleTrack.size()-1;it++) { //.....select first track:
00156       WeveEleTrack &T1=V.eleTrack[it];
00157       if(T1.isMatch2Cl==false) continue;
00158       assert(T1.cluster.nTower>0); // internal logical error
00159       assert(T1.nearTotET>0); // internal logical error
00160 
00161       float isoET1=T1.cluster.ET /T1.cl4x4.ET;
00162       hA[29]->Fill(isoET1);
00163 
00164       hA[23]->Fill(T1.cluster.ET);
00165       hA[0]->Fill("tr1",1.);
00166       if(T1.cluster.ET<par_clusterEtZ) continue;
00167       hA[0]->Fill("et1",1.);
00168 
00169       float fracET1=T1.cluster.ET /T1.nearTotET;
00170       hA[24]->Fill(fracET1);
00171       if(fracET1< par_nearTotEtFracZ) continue; 
00172       hA[0]->Fill("con1",1.);
00173 
00174       for (uint it2=it+1;it2<V.eleTrack.size();it2++) { //.....select second track:
00175         WeveEleTrack &T2=V.eleTrack[it2];
00176         if(T2.isMatch2Cl==false) continue;
00177         assert(T2.cluster.nTower>0); // internal logical error
00178         assert(T2.nearTotET>0); // internal logical error
00179 
00180         float isoET2=T2.cluster.ET /T2.cl4x4.ET;
00181         hA[30]->Fill(isoET2);
00182 
00183         hA[25]->Fill(T2.cluster.ET);
00184         hA[0]->Fill("tr2",1.);
00185         if(T2.cluster.ET<par_clusterEtZ) continue;
00186         hA[0]->Fill("et2",1.);
00187 
00188         float fracET2=T2.cluster.ET /T2.nearTotET;
00189         hA[26]->Fill(fracET2);
00190         if(fracET2< par_nearTotEtFracZ) continue; 
00191         hA[0]->Fill("con2",1.);
00192 
00193         float e1=T1.cluster.energy;
00194         float e2=T2.cluster.energy;
00195         TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
00196         TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
00197 
00198         float del_phi=p1.DeltaPhi(p2);
00199         //printf("del Phi=%f\n",del_phi);
00200         float xx=del_phi;
00201         if(xx<-PI+1) xx+=2*PI;
00202         hA[27]->Fill(xx);
00203         if(fabs(del_phi)<par_delPhi12) continue;
00204         hA[0]->Fill("phi12",1.);
00205 
00206         TVector3 psum=p1+p2;
00207         float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
00208         if(mass2<1.) continue; // 9GeV^2) should be param, I'm tired today
00209         hA[0]->Fill("m2",1.);
00210 
00211         float mass=sqrt(mass2);
00212         int Q1Q2=T1.prMuTrack->charge()*T2.prMuTrack->charge();
00213         if (Q1Q2==1) { //..  same sign , can't be Z-> e+ e-
00214           hA[14]->Fill(mass);
00215           hA[54]->Fill(calcMass(T1,T2));
00216           hA[61]->Fill(zdcRate,calcMass(T1,T2));
00217           continue;
00218         }
00219 
00220         //..... now only opposite sign
00221         hA[0]->Fill("QQ",1.);
00222         hA[15]->Fill(mass);
00223         hA[33]->Fill(T1.cluster.ET,T1.prMuTrack->charge()/T1.prMuTrack->pt()); 
00224         hA[33]->Fill(T2.cluster.ET,T2.prMuTrack->charge()/T2.prMuTrack->pt()); 
00225         hA[34]->Fill(T1.pointTower.iEta ,T1.cluster.energy);
00226         hA[34]->Fill(T2.pointTower.iEta ,T2.cluster.energy);
00227         hA[44]->Fill(calcMass(T1,T2));
00228         hA[60]->Fill(zdcRate,calcMass(T1,T2));
00229 
00230 #if 0
00231         printf("RCC:  Found Z w/ invmass=%f\n",mass);
00232         printJan(&T1);
00233         printJan(&T2);
00234 
00235 
00236         if (!wMK->isMC || (wMK->isMC&& wEve.id<500) )
00237           { printf("\n ZZZZZZZZZZZZZZZZZZZ\n");
00238             if(mass<par_minMassZ) 
00239             wMK->wDisaply->exportEvent("Zlow",V,T1);
00240             else
00241               wMK->wDisaply->exportEvent("Zgood",V,T1);
00242             printf("RCC:  Found Z w/ invmass=%f\n",mass);
00243             wEve.print();
00244           }
00245         
00246 #endif
00247 
00248         if (mass<par_minMassZ) continue; //enforce a lower bound
00249         hA[0]->Fill("Zlow",1.);
00250 
00251         if (mass>par_maxMassZ) continue; //enforce an upper bound
00252         hA[0]->Fill("Zhigh",1.);
00253 
00254         // **** I stoped changes here, Jan 
00255 
00256         float fmax1=T1.cluster.ET/T1.cl4x4.ET;
00257         float fmax2=T2.cluster.ET/T2.cl4x4.ET;
00258 
00259         hA[21]->Fill(fmax1,fmax2);
00260         hA[22]->Fill(T1.cluster.ET,T2.cluster.ET);
00261 
00262         hA[1]->Fill(mass);
00263         hA[2]->Fill(T1.prMuTrack->charge(),T2.prMuTrack->charge());
00264         hA[3]->Fill(T1.prMuTrack->charge()*T2.prMuTrack->charge());
00265         hA[4]->Fill(p1.Phi(),p2.Phi());
00266         hA[5]->Fill(del_phi);
00267         hA[6]->Fill(mass,T1.prMuTrack->charge()/T1.primP.Perp()*T2.prMuTrack->charge()/T1.primP.Perp());
00268         hA[7]->Fill(mass,T1.prMuTrack->charge()*T2.prMuTrack->charge());
00269         hA[8]->Fill(T1.cluster.ET);
00270         if (T1.prMuTrack->charge()>0)
00271           {
00272             hA[9]->Fill(p1.Eta(),p1.Phi());
00273             hA[10]->Fill(p2.Eta(),p2.Phi());
00274           }
00275         else
00276           {
00277             hA[10]->Fill(p1.Eta(),p1.Phi());
00278             hA[9]->Fill(p2.Eta(),p2.Phi());
00279           }
00280 
00281         hA[11]->Fill(fmax1,fmax2);
00282         hA[12]->Fill(T1.cluster.ET,T2.cluster.ET);
00283         hA[13]->Fill(mass,del_phi);
00284 
00285       }
00286 
00287 
00288     
00289     }// loop over first track
00290   }// loop over vertices
00291  
00292 }
00293 
00294 
00295 //_____________________________________________________________________________
00296 //
00297 float
00298 St2009ZMaker::calcMass(WeveEleTrack T1, WeveEleTrack T2){
00299   float e1=T1.cluster.energy;
00300   float e2=T2.cluster.energy;
00301   TVector3 p1=T1.primP; p1.SetMag(e1);//cluster.position;
00302   TVector3 p2=T2.primP; p2.SetMag(e2);//cluster.position;
00303   TVector3 psum=p1+p2;
00304   float mass2=(e1+e2)*(e1+e2)-(psum.Dot(psum));
00305   if(mass2<1.) return 0; 
00306   float mass=sqrt(mass2);
00307   
00308   return mass;
00309 }
00310 
00311 // $Log: St2009ZMaker.cxx,v $
00312 // Revision 1.11  2011/09/14 14:23:21  stevens4
00313 // update used for cross section PRD paper
00314 //
00315 // Revision 1.10  2010/05/01 01:31:44  balewski
00316 // added W->JJ code & JES calibration
00317 //
00318 // Revision 1.9  2010/01/10 03:01:37  balewski
00319 // cleanup & nicer histos
00320 //
00321 // Revision 1.8  2010/01/06 14:11:13  balewski
00322 // one Z-plot added
00323 //
00324 // Revision 1.7  2010/01/06 05:21:57  balewski
00325 // cleanup
00326 //
00327 // Revision 1.6  2010/01/06 04:22:15  balewski
00328 // added Q/PT plot for Zs, more cleanup
00329 //
00330 // Revision 1.5  2010/01/05 03:22:55  balewski
00331 // change logic for filling btow status tables, added printout to Z-code
00332 //
00333 // Revision 1.4  2010/01/04 05:12:00  balewski
00334 // added 4x4 cut to Z-algo, cleanup
00335 //
00336 // Revision 1.3  2010/01/03 04:38:24  balewski
00337 // reorganized Z-algo
00338 //
00339 // Revision 1.2  2009/12/11 20:37:09  rcorliss
00340 // Updated Z Maker to use track direction rather than tower direction when computing invariant mass.  Added new histograms.
00341 //
00342 // Revision 1.1  2009/12/07 20:37:56  rcorliss
00343 // Start
00344 //