StRoot: StSpinPool/StWalgoB2009/St2009WjjMaker.cxx Source File

00001 // $Id: St2009WjjMaker.cxx,v 1.5 2010/06/25 15:42:09 balewski Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 // 
00005 #include <math.h>
00006 
00007 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00008 #include <StMuDSTMaker/COMMON/StMuDst.h>
00009 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00010 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
00011 
00012 #include <StSpinPool/StSpinDbMaker/StSpinDbMaker.h>
00013 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
00014 
00015 #include "St2009WMaker.h"
00016 
00017 #include "St2009WjjMaker.h"
00018 #include <StSpinPool/StJets/StJet.h>
00019 
00020 ClassImp(St2009WjjMaker)
00021 
00022 //_____________________________________________________________________________
00023 //
00024 St2009WjjMaker::St2009WjjMaker(const char *name):StMaker(name){
00025   wMK=0;HList=0;
00026   core=name;
00027  
00028   par_jetPtLow=5.; par_jetPtHigh=60.;// GeV/c
00029   par_jetEtaLow=1.0;  par_jetEtaHigh=1.8;
00030 
00031   par_djPtLow=1.;  par_djPtHigh=40;// GeV/c
00032   par_djPzLow=3; par_djPzHigh=90; // GeV/c
00033   par_djEtaMin=0.1;
00034   par_etaSumLow=0.01; par_etaSumHigh=2.5;
00035 
00036   par_spinSort=false;
00037   par_vertexZ=100;// cm
00038   isMC=0;
00039   spinDb=0;
00040   par_corLevel=0; // default no corrections
00041   mJEScorrFile="fixMe"; memset(mJEScorrH,0,sizeof(mJEScorrH));
00042  }
00043 
00044 
00045 //_____________________________________________________________________________
00046 //
00047 Int_t 
00048 St2009WjjMaker::Init(){
00049   assert(wMK);
00050   assert(HList);
00051   initHistos();
00052   nRun=0;
00053   
00054   LOG_INFO<<GetName()<<Form("::Init JES corr=%d, input=%s",par_corLevel,mJEScorrFile.Data())<<endl;;
00055   if(par_corLevel) {
00056     TFile *fd=new TFile(mJEScorrFile); assert(fd);
00057     for(int i=0;i<mxJESeta;i++) {
00058       TString tit= "jesCorr_iEta"; tit+=i;
00059       mJEScorrH[i]=(TH1F*)fd->Get(tit); assert(mJEScorrH[i]);
00060       mJEScorrH[i]->Print();
00061 
00062       // add derivate for interpolation    
00063       TH1F *h=mJEScorrH[i];      
00064       TAxis *ax=h->GetXaxis();      int nb=ax->GetNbins(); 
00065       for(int k=1;k<=nb-1;k++) {
00066         float x1=ax->GetBinCenter(k);   
00067         float x2=ax->GetBinCenter(k+1); 
00068         float y1=h->GetBinContent(k);   
00069         float y2=h->GetBinContent(k+1); 
00070         float tg=(y2-y1)/(x2-x1);       
00071         h->SetBinError(k,tg);     
00072       }    
00073     } // end of loop over eta bins
00074   }// end of JES initialialization
00075 
00076   return StMaker::Init();
00077 }
00078   
00079 //_____________________________________________________________________________
00080 //
00081 TLorentzVector  
00082 St2009WjjMaker::trueJet( TLorentzVector rJ) {
00083   if (par_corLevel==0) return rJ;
00084   int iEta=(rJ.Eta()+1.2)/0.4;
00085   if(iEta<0) iEta=0;
00086   if(iEta>=mxJESeta ) iEta=mxJESeta-1;
00087   TH1F* h=mJEScorrH[iEta];
00088   float rPt=rJ.Pt();
00089   int bin=h->FindBin(rPt);
00090   float x1=h->GetBinCenter(bin);
00091 
00092   if(x1>rPt && bin>1) { // use lower pT bin for interpolation
00093     bin--;
00094     x1=h->GetBinCenter(bin);
00095   }
00096   float y1=h->GetBinContent(bin);
00097   float tg=h->GetBinError(bin);
00098   
00099   float truePt=y1+ tg*(rPt-x1);
00100   
00101   float fac=truePt/rJ.Pt();
00102   //printf(" rPt=%.2f  tPt=%.2f  ratio=%.3f\n",rJ.Pt(),truePt ,fac);
00103 
00104   
00105   return fac*rJ;
00106 }
00107 
00108 
00109 //_____________________________________________________________________________
00110 //
00111 Int_t 
00112 St2009WjjMaker::FinishRun  (int runNo){
00113   char txt[1000];
00114 
00115   if(par_spinSort) {
00116     sprintf(txt,"events T= %d %d",Tfirst,Tlast);
00117     printf("Finish run=%d , events time range %s\n",runNo,txt);
00118     hbxIdeal->GetYaxis()->SetTitle(txt);
00119   }
00120   return kStOK;
00121 }
00122 
00123 //_____________________________________________________________________________
00124 //
00125 Int_t 
00126 St2009WjjMaker::InitRun  (int runNo){
00127 
00128   if(par_spinSort) {
00129     assert(spinDb);
00130     assert(runNo>= 10081007); // F10407, first pp500 long fill with defined pol pattern
00131     assert(runNo<=10103046); // F10536, last pp500 fill in run9
00132     
00133     char txt[1000],txt0[100];
00134     sprintf(txt0,"bxIdeal%d",nRun);
00135     sprintf(txt,"intended fill pattern  R%d-%d vs. bXing; %s", runNo,nRun,spinDb->getV124comment());
00136     nRun++;
00137     Tfirst=int(2e9); Tlast=-Tfirst;
00138     hbxIdeal=new TH1F(core+txt0,txt,128,-0.5,127.5);
00139     hbxIdeal->SetFillColor(kYellow);
00140     HList->Add(hbxIdeal);
00141   
00142     spinDb->print(0); // 0=short, 1=huge
00143     for(int bx=0;bx<120;bx++){
00144       if(spinDb->isBXfilledUsingInternalBX(bx))  hbxIdeal->Fill(bx);   
00145     }
00146     
00147     sprintf(txt,"bXing= bx7+off=%d",spinDb->BX7offset());
00148     hA[4]->GetXaxis()->SetTitle(txt);
00149   } // end of spin sort
00150     
00151 
00152   LOG_INFO<<GetName()<<Form("::InitRun(%d) done, W->jet+jet sorting  params: doSpinSort=%d |vertZ|<%.0f cm,\n   jetPt=[%.1f,%.1f] GeV/c, jetEta=[%.1f,%.1f]\n   DJ: pT=[%.1f,%.1f] GeV/c, |Pz|=[%.1f,%.1f] GeV/c, eta1+2=[%.1f,%.1f]",
00153                  runNo,par_spinSort,par_vertexZ ,par_jetPtLow,par_jetPtHigh,par_jetEtaLow,par_jetEtaHigh,
00154                  par_djPtLow,par_djPtHigh,par_djPzLow,par_djPzHigh,par_etaSumLow, par_etaSumHigh
00155                  )<<endm;        
00156   return kStOK;
00157 }
00158 
00159 //_____________________________________________________________________________
00160 //
00161 Int_t 
00162 St2009WjjMaker::Make(){
00163   int T=wMK->mMuDstMaker->muDst()->event()->eventInfo().time();
00164   if(Tlast<T) Tlast=T;
00165   if(Tfirst>T) Tfirst=T;
00166 
00167   bXingSort();
00168   return kStOK;
00169 }
00170 
00171 //_____________________________________________________________________________
00172 //
00173 void 
00174 St2009WjjMaker::bXingSort(){
00175   //has access to whole W-algo-maker data via pointer 'wMK'
00176   
00177   hA[0]->Fill("inp",1.);
00178 
00179   //...... trigger .........
00180 
00181   if(!isMC ){ // fixed n the middle of processing
00182     StMuEvent* muEve = wMK->mMuDstMaker->muDst()->event();
00183     StMuTriggerIdCollection *tic=&(muEve->triggerIdCollection());
00184     assert(tic);
00185     const StTriggerId &l1=tic->l1();
00186     vector<unsigned int> idL=l1.triggerIds();
00187     //  printf("nTrig=%d, trigID: ",idL.size());
00188     int trgOK=0;
00189     for(unsigned int i=0;i<idL.size(); i++){
00190       if(idL[i]==230420) trgOK+=1; // AJP
00191       if(idL[i]==230411) trgOK+=2; // JP2
00192     }
00193     if(trgOK)    hA[0]->Fill("trig",1.);
00194     // tmp, do not abort events here, now chain is trigger filtering
00195   }
00196   
00197   //....... find vertex .....
00198   int nInpPrimV=wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices();
00199   int nVer=0;
00200   for(int iv=0;iv<nInpPrimV;iv++) {
00201     if(nVer) break; // tmp - consider only 1st vertex, since jets are made only for the 1st vertex
00202     StMuPrimaryVertex* V= wMK->mMuDstMaker->muDst()->primaryVertex(iv);
00203     assert(V);
00204     wMK->mMuDstMaker->muDst()->setVertexIndex(iv);
00205     float rank=V->ranking();
00206     if (rank<=0) continue;
00207     const StThreeVectorF &r=V->position();
00208     if(fabs(r.z()) > par_vertexZ) continue;
00209     nVer++; // count valid vertices
00210   }
00211 
00212   if(nVer<=0) return; 
00213   hA[0]->Fill("vert",1.);
00214 
00215   //......... require:  vertex is reasonable .......
00216   
00217 
00218   //****loop over jet branch with B+EEMC****
00219   if(GetMaker("JetReader")==0) return; // most likely jet-maker was not in the chain
00220   TClonesArray* jets = wMK->getJets("ConeJets12_100"); // select specific jet-tree type 
00221   if(jets==0) return; // most likely jet-maker was not in the chain
00222   int nJets= wMK->nJets;
00223   //printf("nJets=%d  n2J=%.0f\n",nJets, hA[0]->GetBinContent(1)); 
00224   if(nJets<1) return;
00225   hA[0]->Fill("anyJ",1.);
00226   if(nJets<2) return;
00227   hA[0]->Fill("mulJ",1.);
00228 
00229   const int mxJ=2;
00230   TLorentzVector jet[mxJ]; 
00231   int nJ=0;
00232   for (int i_jet=0; i_jet< nJets; i_jet++){// try to find 2 jets passing cuts
00233     TLorentzVector Jreco = *((StJet*)jets->At(i_jet));
00234 
00235     //...... compute correction .......
00236     TLorentzVector J=Jreco; //1:1
00237     J=trueJet(J);
00238 
00239     if(J.Pt()<par_jetPtLow) continue;
00240 
00241     if(nJ>=mxJ) { // 3rd jet, kill for now 
00242       hA[0]->Fill("J3",1.);
00243       hA[6]->Fill(fabs(jet[0].DeltaPhi(J)),fabs(jet[1].DeltaPhi(J)));
00244       return;
00245     }
00246     jet[nJ++]=J;
00247     if(nJ==1)hA[0]->Fill("J1",1.);
00248   }
00249   
00250   if(nJ<mxJ) return;
00251   hA[0]->Fill("J2",1.);
00252 
00253   // 2 jet-event, process it
00254   assert(nJ<=mxJ);
00255 
00256   // decide which is jet #1
00257   if( fabs(jet[0].Eta()) > fabs(jet[1].Eta()) ) {// swap
00258     TLorentzVector J=jet[0];
00259     jet[0]=jet[1]; jet[1]=J;
00260   }
00261 
00262   //  now jet1 has smaller |eta|
00263   hA[10]->Fill(jet[0].Et(),jet[0].Eta());
00264   hA[11]->Fill(jet[1].Et(),jet[1].Eta());
00265   hA[17]->Fill(jet[0].E(),jet[1].E());
00266   hA[18]->Fill(jet[0].Pt(),jet[1].Pt());
00267   hA[19]->Fill(jet[0].Eta(),jet[1].Eta());
00268   float phiCDdeg=jet[0].DeltaPhi(jet[1])/3.1416*180.;
00269   if(phiCDdeg<-90) phiCDdeg+=360; // choose different phi range
00270   hA[13]->Fill(phiCDdeg); 
00271 
00272   
00273   TLorentzVector diJet=jet[0]+jet[1];
00274   float invM=sqrt(diJet*diJet);
00275  
00276   hA[14]->Fill(invM);
00277   hA[16]->Fill(invM, diJet.Pt());
00278   hA[15]->Fill( diJet.Z(),diJet.Pt());
00279  
00280 
00281   // for Pavel N.
00282   if(invM<60) {
00283       hA[21]->Fill(diJet.Pt());
00284       hA[22]->Fill(fabs(jet[0].DeltaPhi(jet[1])));
00285       hA[23]->Fill(jet[0].Eta()-jet[1].Eta());
00286       hA[24]->Fill(jet[0].Eta(),jet[1].Eta());
00287   }
00288 
00289   if(par_spinSort) {//........ do spin sorting
00290     
00291     assert(spinDb->isValid());  // all 3 DB records exist 
00292     assert(spinDb->isPolDirLong());  // you do not want mix Long & Trans by accident
00293     int bx48=wMK->wEve.bx48;
00294     int bx7=wMK->wEve.bx7;
00295     if(spinDb->offsetBX48minusBX7(bx48,bx7)) {
00296       printf("BAD bx7=%d bx48=%d del=%d\n",bx7,bx48,spinDb->offsetBX48minusBX7(bx48,bx7));
00297       hA[0]->Fill("badBx48",1.);
00298       return; // both counters must be in sync
00299     }
00300     hA[2]->Fill(bx7);
00301     int bxStar7=spinDb->BXstarUsingBX7(bx7);
00302     hA[4]->Fill(bxStar7);
00303     
00304     int spin4=spinDb->spin4usingBX48(bx48); 
00305     hA[5]->Fill(bxStar7,spin4);
00306     
00307     hA[20]->Fill(invM,spin4);
00308     
00309   } // end of spin sorting
00310 
00311   /*   TLorentzVector (px,py,pz,E). 
00312        TLorentzVector v4(TVector3(1., 2., 3.),4.);
00313        v.SetVect(TVector3(1,2,3)); 
00314        v.SetXYZT(x,y,z,t); 
00315        v.SetPxPyPzE(px,py,pz,e); 
00316        v.SetXYZM(x,y,z,m);   //   ->  v=(x,y,z,e=Sqrt(x*x+y*y+z*z+m*m))
00317   */
00318 
00319 
00320 #if 0
00321   for (int i_jet=1; i_jet< nJetsWE; i_jet++){//loop over jets
00322     StJet* jet = wMK->getJet(i_jet);
00323     TVector3 jetVec; //vector for jet momentum
00324         //vary neutral and charged et in jets for systematic
00325         float neutral=jet->neutralFraction()*jet->Pt();
00326         float charged=jet->chargedFraction()*jet->Pt();
00327         neutral=neutral*par_mcJetNeutScale;
00328         charged=charged*par_mcJetChrgScale;
00329         float sum=neutral+charged;
00330         jetVec.SetPtEtaPhi(sum,jet->Eta(),jet->Phi());
00331         if(jetVec.DeltaR(T.primP) > par_nearDeltaR)
00332               T.ptBalance+=jetVec;
00333       }
00334       TVector3 clustPt(T.primP.X(),T.primP.Y(),0);
00335       clustPt.SetMag(T.cluster.ET);
00336       T.ptBalance+=clustPt;
00337       T.sPtBalance = T.ptBalance.Perp();
00338       if(T.ptBalance.Dot(clustPt)<0) T.sPtBalance *=-1.;
00339       
00340 
00341 #endif
00342 
00343   if(par_spinSort) {
00344     assert(spinDb->isValid());  // all 3 DB records exist 
00345     assert(spinDb->isPolDirLong());  // you do not want mix Long & Trans by accident  
00346     int bx48=wMK->wEve.bx48;
00347     int bx7=wMK->wEve.bx7;
00348     if(spinDb->offsetBX48minusBX7(bx48,bx7)) {
00349       printf("BAD bx7=%d bx48=%d del=%d\n",bx7,bx48,spinDb->offsetBX48minusBX7(bx48,bx7));
00350       hA[0]->Fill("badBx48",1.);  
00351       return;   
00352     }
00353         
00354     hA[2]->Fill(bx7);
00355     
00356     int bxStar7=spinDb->BXstarUsingBX7(bx7);
00357     
00358     hA[4]->Fill(bxStar7);
00359     
00360     int spin4=spinDb->spin4usingBX48(bx48); 
00361     hA[5]->Fill(bxStar7,spin4);
00362   }// end of spin processing  
00363     
00364 #if 0
00365   StMuEvent* muEve = wMK->mMuDstMaker->muDst()->event();
00366   int max=0;
00367     for (int m=0;m<300;m++)  {    // access L0-HT data
00368       int val=muEve->emcTriggerDetector().highTower(m);
00369       if(max<val) max=val;
00370     }
00371 #endif
00372 
00373 
00374       //::::::::::::::::::::::::::::::::::::::::::::::::
00375       //:::::accepted W events for x-section :::::::::::
00376       //::::::::::::::::::::::::::::::::::::::::::::::::
00377 
00378 }
00379 
00380 
00381 // $Log: St2009WjjMaker.cxx,v $
00382 // Revision 1.5  2010/06/25 15:42:09  balewski
00383 // more plots for Pavel
00384 //
00385 // Revision 1.4  2010/05/04 12:14:35  balewski
00386 // runs now w/o jet tree
00387 //
00388 // Revision 1.3  2010/05/03 17:24:37  balewski
00389 // added spin sorting of di-jets
00390 //
00391 // Revision 1.2  2010/05/01 01:31:44  balewski
00392 // added W->JJ code & JES calibration
00393 //
00394 // Revision 1.1  2010/04/16 01:04:43  balewski
00395 // start
00396 //
00397