StRoot: StSpinPool/StWalgoB2009/StMcJetCalibMaker.cxx Source File

00001 // $Id: StMcJetCalibMaker.cxx,v 1.1 2010/05/01 01:31:45 balewski Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 // 
00005 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00006 #include <StMuDSTMaker/COMMON/StMuDst.h>
00007 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00008 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.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 
00018 #include "St2009WMaker.h"
00019 #include "St2009WjjMaker.h" // for jetEcorr
00020 #include "StMcJetCalibMaker.h"
00021 #include <StSpinPool/StJets/StJet.h>
00022 
00023 ClassImp(StMcJetCalibMaker)
00024  
00025 //________________________________________________________
00026 //
00027 StMcJetCalibMaker::StMcJetCalibMaker(const char *name):StMaker(name){
00028   wMK=0; wjjMK=0;HList=0;
00029   core=name;
00030   par_vertexZ=70;// cm
00031   par_verZerr=1;//cm
00032   par_delRmax=0.1;// match gen & reco jets
00033   par_corLevel=0; // level of jet energy correction
00034   par_jetEtaLow=-1.5; par_jetEtaHigh=2; 
00035   isMC=1; 
00036 }
00037 
00038 
00039 //_____________________________________________________________________________
00040 //
00041 Int_t 
00042 StMcJetCalibMaker::Init(){
00043   assert(wMK);
00044   if( par_corLevel)  assert(wjjMK);
00045   assert(HList);
00046   initHistos();
00047 
00048   return StMaker::Init();
00049 }
00050 
00051 
00052 //_____________________________________________________________________________
00053 //
00054 Int_t 
00055 StMcJetCalibMaker::FinishRun  (int runNo){
00056   return kStOK;
00057 }
00058 
00059 //_____________________________________________________________________________
00060 //
00061 Int_t 
00062 StMcJetCalibMaker::InitRun  (int runNo){
00063 
00064 
00065   LOG_INFO<<GetName()<<Form("::InitRun(%d),calib params: |geant vertZ|<%.0f cm  |zG-zR|<%.1f cm, jet: eta=[%.1f,%1f],  delR(reco-gen)<%.2f \n  JES corr level=%d",
00066                             runNo,par_vertexZ ,par_verZerr,par_jetEtaLow,par_jetEtaHigh, 
00067                            par_delRmax,
00068                             par_corLevel
00069                  )<<endm;        
00070   assert(isMC);
00071   return kStOK;
00072 }
00073 
00074 //_____________________________________________________________________________
00075 //
00076 Int_t 
00077 StMcJetCalibMaker::Make(){
00078   calibrate();
00079 
00080   return kStOK;
00081 }
00082 
00083 //_____________________________________________________________________________
00084 //
00085 void 
00086 StMcJetCalibMaker::calibrate(){
00087   StMcEvent* mMcEvent = 0;
00088   mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent");
00089   if(mMcEvent==0) return; // STARNGE - fix it
00090   assert(mMcEvent);
00091   
00092 
00093   //........ restrict M-C vertex .........
00094 
00095   StMcVertex *V=mMcEvent->primaryVertex();
00096   TVector3 mcV=TVector3(V->position().x(),V->position().y(),V->position().z());
00097   //printf("MC ver=%f %f %f\n", mcV.x(), mcV.y(), mcV.z());
00098   
00099   hA[0]->Fill("inp",1.);
00100   hA[1]->Fill(mcV.z());
00101   if(fabs(mcV.z()) > par_vertexZ) return;
00102   hA[0]->Fill("verG",1.);
00103 
00104 #if 0 
00105   for(uint i=0;i<mMcEvent->tracks().size();i++){//loop tracks
00106     StMcTrack* mcTrack = mMcEvent->tracks()[i];
00107     int pdgId=mcTrack->pdgId();
00108     const StLorentzVectorF &p4=mcTrack-> fourMomentum();
00109     printf("i=%d pdgId=%d  p=%f %f %f  m=%f\n",i,pdgId,p4.x(),p4.y(),p4.z(),p4.m());
00110     if(i>15) break;
00111   }
00112 #endif
00113  
00114    //........ select Pythia particles .........
00115   int type=0; // QCD process
00116   int pdgId=mMcEvent->tracks()[6]->pdgId(); // boson?
00117   if(pdgId==23) type=1; //Z
00118   if( pdgId==-24 || pdgId==24 ) type=2;// W+, W-
00119 
00120   // A+B --> C+D  or  A+B-->W --> C+D
00121  
00122   if(mMcEvent->tracks().size()<10) return;
00123   // access 2:2 kinematics for W, Z production or QCD scattering
00124   //  StMcTrack* mctA= mMcEvent->tracks()[4]; // q or qbar 1
00125   // StMcTrack* mctB= mMcEvent->tracks()[5]; //  q or qbar 2
00126 
00127   StMcTrack* mctW=0, *mctC=0, *mctD=0;
00128   if(type) {// Z,W+, or W-
00129     mctW=mMcEvent->tracks()[6]; // boson
00130     mctC= mMcEvent->tracks()[7];// decay1 (lepton or q or qbar)
00131     mctD= mMcEvent->tracks()[8];// decay 2 product
00132     hA[0]->Fill("W,Z",1.);
00133   } else { // QCD process
00134     mctC= mMcEvent->tracks()[6];// paron1 
00135     mctD= mMcEvent->tracks()[7];// parton 2
00136   }
00137 
00138   //... define matching array for gen vs. thrown jets
00139   const int mxGJ=2, mxRJ=4;
00140   int nGJ=0,nRJ=0 , mrJI[mxGJ];;
00141   TLorentzVector gJA[mxGJ], rJA[mxRJ];
00142 
00143   float delRA[mxGJ][mxRJ],delRAc[mxGJ][mxRJ];
00144   memset(mrJI,-1,sizeof(mrJI)); // clear all gen-jets as not-matched
00145 
00146   for(int k=0;k<2;k++) { // find parton in the barrel
00147     const StLorentzVectorF &ttC=mctC-> fourMomentum();
00148     const StLorentzVectorF &ttD=mctD-> fourMomentum();
00149     TLorentzVector gJ;
00150     if(k==0) gJ=TLorentzVector(ttC.px(),ttC.py(),ttC.pz(),ttC.e());
00151     if(k==1) gJ=TLorentzVector(ttD.px(),ttD.py(),ttD.pz(),ttD.e());
00152     hA[11]->Fill(gJ.Eta(),gJ.Phi());
00153     if(gJ.Eta()<par_jetEtaLow )  continue;
00154     if(gJ.Eta()>par_jetEtaHigh ) continue;
00155     assert(nGJ<mxGJ);
00156     gJA[nGJ++]=gJ;
00157     hA[4]->Fill(gJ.Pt(),gJ.Eta());
00158   }
00159   
00160   if(!nGJ) return;
00161   hA[0]->Fill("etaG",1.);
00162 
00163   if(type) {
00164     const StLorentzVectorF &pW=mctW-> fourMomentum();
00165     hA[10]->Fill(pW.m());
00166   }
00167 
00168 
00169   //....... find  reco 1st vertex .....
00170 
00171   int nInpPrimV=wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices();
00172   int nVer=0;
00173   for(int iv=0;iv<nInpPrimV;iv++) {
00174     if(nVer) break; // consider only 1st valid vertex
00175     StMuPrimaryVertex* V= wMK->mMuDstMaker->muDst()->primaryVertex(iv);
00176     assert(V);
00177     wMK->mMuDstMaker->muDst()->setVertexIndex(iv);
00178     float rank=V->ranking();
00179     if (rank<=0) continue;
00180     const StThreeVectorF &r=V->position();
00181     hA[2]->Fill(r.z()-mcV.z());
00182     if(fabs(r.z()-mcV.z()) > par_verZerr) continue;
00183     nVer++; // count valid vertices
00184     //printf("  reco ver=%f %f %f\n", r.x(), r.y(), r.z());
00185   }
00186 
00187   if(nVer<=0) return; 
00188   hA[0]->Fill("verR",1.);
00189 
00190 
00191   //....... select reco jets  ........
00192 
00193   TClonesArray* jets = wMK->getJets("ConeJets12_100"); // select jet-tree type 
00194   int nJets= wMK->nJets;
00195   //printf("nJets=%d  n2J=%.0f\n",nJets, hA[0]->GetBinContent(1)); 
00196 
00197   if(nJets<1) return;
00198   hA[0]->Fill("anyJ",1.);
00199   if(nJets>=2)  hA[0]->Fill("mulJ",1.);
00200   
00201   // ..... compute R-distance for all possible pairs (gen_i,reco_j) ...
00202   float par_recJetPtMin=5.; 
00203 
00204   for (int ir=0; ir< nJets; ir++){
00205      StJet* jet = (StJet*)jets->At(ir);
00206      TLorentzVector rJ = *jet;
00207      if(par_corLevel) rJ=wjjMK->trueJet(rJ);// apply JES correction
00208 
00209      if(rJ.Pt()<par_recJetPtMin) continue;
00210      if(nRJ>=mxRJ) continue; // skip 
00211      for(int ig=0;ig<nGJ;ig++)  // compute R
00212        delRA[ig][nRJ]=gJA[ig].DrEtaPhi(rJ);
00213      rJA[nRJ++]=rJ; // store jet
00214   }
00215   //  printf("aa nRJ=%d\n",nRJ);
00216   
00217 
00218   memcpy(delRAc,delRA,sizeof(delRA));// just for QA & printouts
00219 
00220   //...... select best matching starting from lowest minimum in R
00221   int nMJ=0;
00222   while(nMJ<nRJ) {
00223     float minR=1e37;
00224     int irm=-1,igm=-1;
00225     for (int ir=0; ir< nRJ; ir++)      
00226       for(int ig=0;ig<nGJ;ig++) {
00227         if(delRA[ig][ir]>minR) continue;
00228         minR= delRA[ig][ir];
00229         igm=ig; irm=ir;
00230       }
00231     mrJI[igm]=irm;
00232     nMJ++;
00233     // disable this  this pair
00234     for (int ir=0; ir< nRJ; ir++)  delRA[igm][ir]=99999; 
00235     for(int ig=0;ig<nGJ;ig++) delRA[ig][irm]=99999;
00236 
00237   }// end of while
00238 
00239   //.... work only with matched jets
00240   for(int ig=0;ig<nGJ;ig++) {
00241     if(mrJI[ig]<0)  continue;
00242     int ir=mrJI[ig];
00243     float delR=delRAc[ig][ir];
00244     TLorentzVector gJ=gJA[ig], rJ=rJA[ ir];
00245     
00246     hA[3]->Fill(delR);
00247     if(delR>par_delRmax) continue;
00248     hA[0]->Fill("delR",1.);
00249 
00250  
00251 
00252     hA[15]->Fill(gJ.Pt());
00253     hA[12]->Fill(gJ.Pt(),gJ.Eta());
00254     hA[21]->Fill(gJ.E(),rJ.E());
00255     hA[22]->Fill(gJ.Pt(),rJ.Pt());
00256 
00257 
00258     int iEta=(rJ.Eta()+1.2)/0.4;
00259     if(iEta<0) iEta=0;
00260     if(iEta>=mxEta) iEta=mxEta-1;
00261     hA[40+iEta]->Fill(gJ.Pt(),rJ.Pt());
00262 
00263     float ratioPT=rJ.Pt()/gJ.Pt();
00264     //float neutFrac=->neutralFraction();
00265     //float chargFrac=jet->chargedFraction();
00266     hA[16]->Fill(ratioPT);
00267 
00268     //... neutral....
00269     //hA[26]->Fill(ratio*neutFrac);
00270  
00271     //.... charged....
00272     //hA[36]->Fill(ratio*chargFrac);
00273 
00274   }
00275 
00276 #if 0  // pppppppppp test printing
00277   TLorentzVector pG;
00278     printf("gen type=%d  pdgId=%d  eveID=%d\n",type,pdgId,wMK->mMuDstMaker->muDst()->event()->eventId());
00279     if(type && par_corLevel==0) {
00280       const StLorentzVectorF &ttW=mctW-> fourMomentum();
00281       pG=TLorentzVector(ttW.px(),ttW.py(),ttW.pz(),ttW.e());
00282       printf("gen W,  P=%.1f %.1f %.1f   E=%.1f  M=%.1f  PT=%.1f eta=%.2f  phi/deg=%.0f\n",pG.X(),pG.Y(),pG.Z(),pG.E(), pG.M(),pG.Pt(),pG.Eta(), pG.Phi()/3.14*180);
00283     }
00284     for(int k=0;k<2;k++) { // find parton in the barrel
00285       const StLorentzVectorF &ttC=mctC-> fourMomentum();
00286       const StLorentzVectorF &ttD=mctD-> fourMomentum();
00287       if(k==0) pG=TLorentzVector(ttC.px(),ttC.py(),ttC.pz(),ttC.e());
00288       if(k==1) pG=TLorentzVector(ttD.px(),ttD.py(),ttD.pz(),ttD.e());
00289       printf("gen q=%c,  P=%.1f %.1f %.1f   E=%.1f  M=%.1f  PT=%.1f eta=%.2f  phi/deg=%.0f\n",'C'+k,pG.X(),pG.Y(),pG.Z(),pG.E(), pG.M(),pG.Pt(),pG.Eta(), pG.Phi()/3.14*180);
00290     }
00291  
00292 
00293     printf("reco inp nJ=%d eveID=%d: selected nGJ=%d nRJ=%d\n",nJets,wMK->mMuDstMaker->muDst()->event()->eventId(),nGJ,nRJ);
00294     for (int i_jet=0; i_jet< nJets; i_jet++){
00295       StJet* jet = (StJet*)jets->At(i_jet);
00296       TLorentzVector J = *jet;
00297       float neutPt=jet->neutralFraction()*J.Pt();
00298       if(J.Pt()<5) continue;
00299       printf("recJ=%d,  P=%.1f %.1f %.1f   E=%.1f  M=%.1f  PT=%.1f neutPt=%.1f  eta=%.2f  phi/deg=%.0f\n",i_jet,J.X(),J.Y(),J.Z(),J.E(), J.M(),J.Pt(),neutPt,J.Eta(), J.Phi()/3.14*180);
00300     }
00301 
00302     printf("matching gen-->reco :\n");
00303     for(int ig=0;ig<nGJ;ig++) {
00304       printf("gen J PT=%.1f : delR(ir)= ",gJA[ig].Pt());
00305       for (int ir=0; ir< nRJ; ir++){
00306         if( mrJI[ig]==ir) printf("*");
00307         else printf(" ");
00308         printf("%.2f, ", delRAc[ig][ir]);
00309       }
00310       if(mrJI[ig]<0)    printf("no match");
00311       else printf(" rec PT=%.1f", rJA[mrJI[ig]].Pt());
00312       printf("\n");
00313     }
00314     printf("\n---------\n");
00315 
00316 #endif
00317 
00318 
00319 
00320 #if 0
00321 
00322 
00323 
00324 
00325     //...... compute correction .......
00326     TLorentzVector pJ=pJreco; //1:1
00327     float cor1=1./St2009WjjMaker::jetEcorr(pJreco.E());
00328 
00329     // apply all corrections
00330     if(par_corLevel) pJ=cor1*pJ;
00331 
00332     float neutFrac=jet->neutralFraction();
00333     float chargFrac=jet->chargedFraction();
00334     
00335     float frac=pJ.E()/pG.E();
00336     if(fabs(frac-1.)> par_ptMargin) continue;
00337     
00338     if(delR>par_delRmargin) continue;
00339     hA[3]->Fill(delR);
00340     if(delR>par_delRmax) continue;
00341     hA[0]->Fill("delR",1.);
00342     hA[12]->Fill(pG.Pt(),pG.Eta());
00343     
00344     //..... study energy resolution.....
00345     //printf("%d %f %f %f \n",par_corLevel,frac,pJ.E(),pG.E());
00346 
00347     //...... full jet energy .....
00348     hA[16]->Fill(frac);
00349     hA[17]->Fill(frac,pG.Eta());
00350     hA[18]->Fill(pG.Eta()-pJ.Eta(),pG.Eta());
00351     hA[19]->Fill(pG.DeltaPhi(pJ)/3.1416*180 ,pG.Eta());
00352     hA[21]->Fill(pG.E(),pJ.E());
00353     hA[22]->Fill(pG.Pt(),pJ.Pt());
00354 
00355     int iEta=(pJ.Eta()+1.)/0.5;
00356     if(iEta<0) iEta=0;
00357     if(iEta>=mxEta) iEta=mxEta-1;
00358     hA[40+iEta]->Fill(pG.Pt(),pJ.Pt());
00359 
00360     //... neutral....
00361     hA[26]->Fill(frac*neutFrac);
00362     hA[27]->Fill(pJ.Eta(),pJ.E()*neutFrac);
00363  
00364     //.... charged....
00365     hA[36]->Fill(frac*chargFrac);
00366     hA[37]->Fill(pJ.Eta(),pJ.E()*chargFrac);
00367     break;
00368   }
00369 
00370 #endif
00371 
00372   /*   TLorentzVector (px,py,pz,E). 
00373        TLorentzVector v4(TVector3(1., 2., 3.),4.);
00374        v.SetVect(TVector3(1,2,3)); 
00375        v.SetXYZT(x,y,z,t); 
00376        v.SetPxPyPzE(px,py,pz,e); 
00377        v.SetXYZM(x,y,z,m);   //   ->  v=(x,y,z,e=Sqrt(x*x+y*y+z*z+m*m))
00378   */
00379   
00380   
00381   
00382   //::::::::::::::::::::::::::::::::::::::::::::::::
00383   //:::::accepted W events for x-section :::::::::::
00384   //::::::::::::::::::::::::::::::::::::::::::::::::
00385   
00386 }
00387 
00388 
00389 // $Log: StMcJetCalibMaker.cxx,v $
00390 // Revision 1.1  2010/05/01 01:31:45  balewski
00391 // added W->JJ code & JES calibration
00392 //
00393 // Revision 1.1  2010/04/16 01:04:43  balewski
00394 // start
00395 //
00396 
00397