StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
St2009WMaker.cxx
1 
2 // $Id: St2009WMaker.cxx,v 1.16 2011/09/14 14:23:20 stevens4 Exp $
3 //
4 //*-- Author : Jan Balewski, MIT
5 //*-- Author for Endcap: Justin Stevens, IUCF
6 //*-- Author for JetFinder/JetReader interface: Ilya Selyuzhenkov, IUCF
7 // test 1, after DNP2009 tag, Jan
8 
9 #include <TH1.h>
10 #include <TH2.h>
11 #include <StMessMgr.h>
12 #include <StThreeVectorF.hh>
13 
14 //MuDst
15 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
16 #include <StMuDSTMaker/COMMON/StMuDst.h>
17 #include <StMuDSTMaker/COMMON/StMuEvent.h>
18 
19 #include "StEmcUtil/database/StBemcTables.h"
20 #include "StEmcUtil/geometry/StEmcGeom.h"
21 
22 #include "StEEmcUtil/database/StEEmcDb.h"
23 #include "StEEmcUtil/database/EEmcDbItem.h"
24 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
25 
26 //jets
27 #include "StSpinPool/StJets/StJet.h"
28 #include "StSpinPool/StJets/StJets.h"
29 #include "StJetMaker/StJetMaker.h"
30 #include "StSpinPool/StSpinDbMaker/StSpinDbMaker.h"
31 #include "StJetMaker/StJetReader.h"
32 #include "StJetMaker/StJetSkimEventMaker.h"
33 
34 
35 #include "WeventDisplay.h"
36 #include "St2009WMaker.h"
37 
38 ClassImp(St2009WMaker)
39 
40 //_____________________________________________________________________________
41 //
42 St2009WMaker::St2009WMaker(const char *name):StMaker(name){
43  char muDstMakerName[]="MuDst";
44  mMuDstMaker= (StMuDstMaker*)GetMaker(muDstMakerName); assert(mMuDstMaker);
45  mJetReaderMaker = (StJetReader*)GetMaker("JetReader");
46  if(mJetReaderMaker ==0) {
47  LOG_WARN<<GetName()<<Form("::constructor() NO JETS , W-algo is not working properly, continue")<<endm;
48  }
49 
50  // preset or clear some params
51  par_bht3TrgID= par_l2wTrgID=0;
52  setHList(0);
53  setMC(0);
54  nInpEve= nTrigEve= nAccEve=0;
55 
56  //... MC trigger simulator
57  par_l0emulAdcThresh=31;
58  par_l2emulSeedThresh=5.0;
59  par_l2emulClusterThresh=13.0;
60 
61  //.. vertex
62  par_minPileupVert=3; // to reject events w/o TPC, lower it for MC
63  par_vertexZ=100; // (cm)
64 
65  //... track
66  par_nFitPts=15; // hits on the track
67  par_nHitFrac=0.51;
68  par_trackRin=90; par_trackRout=160; // cm
69  par_trackPt=10.;//GeV
70 
71  //... ETOW
72  par_kSigPed=3; // rawADC-ped cut off
73  par_AdcThres=8; // ADC threshold to avoid correlated noise
74  par_maxADC=200.; // (adc chan) on the highest tower in events
75  par_clustET=15.; // (GeV/c) 2x2 cluster ET
76  par_clustFrac24=0.95; // ET ratio 2x2/4x4 cluster
77  par_nearDeltaR=0.7; //(~rad) near-cone size
78  par_nearTotEtFrac=0.88; // ratio 2x2/near Tot ET
79  par_smallNearDeltaR=0.1; //(~rad) small near-cone size for tpc tracks
80 
81  //... track-cluster
82  par_delR3D=7.; // cm, dist between projected track and center of cluster
83 
84  //... near-,away-side track counter thresholds
85  par_countTrPt=1.; // (GeV), track pT thres
86  par_countTowEt=1.; // (GeV), tower ET thres
87 
88  //... search for W's
89  par_awayDeltaPhi=0.7; // (rad) away-'cone' size
90  par_highET=25.; // (GeV), cut-off for final W-cluster ET
91  par_ptBalance=15.; // (GeV), ele cluster vector + jet sum vector
92  par_leptonEta=1.0; // bracket acceptance
93 
94  par_useEtow=1; // flag for how to use ETOW in algo
95  // flag == 0 -> don't use ETOW
96  // flag >= 1 -> use in event display plots
97  // flag >= 2 -> use in away sum
98  // flag == 3 -> use in near sum
99  setEtowScale(1.0); // for old the Endcap geometr you need ~1.3
100  setBtowScale(1.0);
101 
102  setJetNeutScaleMC(1.0);//vary neutral ET scale for systematic
103  setJetChrgScaleMC(1.0);//vary charged ET scale for systematic
104 
105  mRunNo=0;
106 
107  // irrelevant for W analysis
108  par_DsmThres=28; // only for monitoring
109  par_maxDisplEve=1; // # of displayed selected events
110 
111  use_gains_file = 0;
112 
113 }
114 
115 
116 //_____________________________________________________________________________
117 //
118 Int_t
119 St2009WMaker::Init(){
120  assert(HList);
121  initHistos();
122 
123  mBarrelTables = new StBemcTables();
124  mBtowGeom = StEmcGeom::instance("bemc");
125  mBSmdGeom[kBSE] = StEmcGeom::instance("bsmde");
126  mBSmdGeom[kBSP] = StEmcGeom::instance("bsmdp");
127  wDisaply= new WeventDisplay(this,par_maxDisplEve);
128 
129  mDbE = (StEEmcDb*)GetDataSet("StEEmcDb");
130  assert(mDbE);
131  geomE= new EEmcGeomSimple();
132  mRand = new TRandom3(0);
133 
134  initGeom();
135 
136  if (use_gains_file == 1) {
137  fstream f1; f1.open(gains_file,ios::in);
138  cout << "Opening gains file " << gains_file << endl;
139  char str[200];
140  while (f1 >> str) {
141  int softID = atoi(str);
142  f1 >> str; gains_BTOW[softID] = atof(str);
143  f1 >> str; f1 >> str;
144  }
145  }
146 
147  if(isMC >= 350){ // load vertex reweighting histo
148  TString filename="/star/u/stevens4/wAnalysis/efficXsec/zVertReweight.root";
149  TFile* reweightFile = new TFile(filename);
150  cout<<"Re-weighting vertex z distribution with '"<<nameReweight<<"' histo from file "<<endl<<filename<<endl;
151  hReweight = (TH1F*) reweightFile->Get(nameReweight);
152  }
153 
154  if(isMC) par_minPileupVert=1;
155 
156  // ..... initialization of TPC cuts is run dependent, call it 'hack of the day', should be moved to InitRun() and handle multipl runs per job, after APS, JB
157  for(int isec=0;isec<mxTpcSec;isec++) {
158  int sec=isec+1;
159  float Rin=par_trackRin,Rout=par_trackRout;
160  //.... Rin ..... changes
161  if(sec==4 && par_inpRunNo>=10090089) Rin=125.;
162  if(sec==11 && par_inpRunNo>=10083013) Rin=125.;
163  if(sec==15 && par_inpRunNo>=10088096 && par_inpRunNo<=10090112 ) Rin=125.;
164  //.... Rout ..... changes
165  if(sec==5 && par_inpRunNo>=10098029) Rout=140.;
166  if(sec==6 ) Rout=140.;
167  if(sec==20 && par_inpRunNo>=10095120 && par_inpRunNo<=10099078 ) Rout=140.;
168 
169  mTpcFilter[isec].setCuts(par_nFitPts,par_nHitFrac,Rin,Rout);
170  mTpcFilter[isec].init("sec",sec,HList);
171  }
172 
173  return StMaker::Init();
174 }
175 
176 //________________________________________________
177 //________________________________________________
178 Int_t
179 St2009WMaker::InitRun(int runNo){
180  LOG_INFO<<Form("::InitRun(%d) start",runNo)<<endm;
181  mBarrelTables->loadTables(this );
182 
183  //get run number from MuDst name for embedding (each file goes through InitRun since each has different geneated runNum < 5000)
184  if(isMC>=350){
185  const char *file = mMuDstMaker->GetFile();
186  const char *p1=strstr(file,"adc_");
187  int run=atoi(p1+4);
188  par_inpRunNo=run;
189 
190  LOG_INFO<<Form("::InitRun(%d) reset TPC cuts",par_inpRunNo)<<endm;
191  for(int isec=0;isec<mxTpcSec;isec++) {
192  int sec=isec+1;
193  float Rin=par_trackRin,Rout=par_trackRout;
194  //.... Rin ..... changes
195  if(sec==4 && par_inpRunNo>=10090089) Rin=125.;
196  if(sec==11 && par_inpRunNo>=10083013) Rin=125.;
197  if(sec==15 && par_inpRunNo>=10088096 && par_inpRunNo<=10090112 ) Rin=125.;
198  //.... Rout ..... changes
199  if(sec==5 && par_inpRunNo>=10098029) Rout=140.;
200  if(sec==6 ) Rout=140.;
201  if(sec==20 && par_inpRunNo>=10095120 && par_inpRunNo<=10099078 ) Rout=140.;
202 
203  //reset cuts for TPC filter based on runNum for embedding
204  if(Rin != par_trackRin || Rout!=par_trackRout)
205  LOG_INFO<<"Sector "<<sec<<" cuts aren't default: Rin="<<Rin<<" and Rout="<<Rout<<endl;
206  mTpcFilter[isec].setCuts(par_nFitPts,par_nHitFrac,Rin,Rout);
207  }
208  }
209 
210  mRunNo=runNo;
211  if(runNo>1000000) assert(mRunNo==par_inpRunNo); // what a hack, assurs TPC cuts change w/ time for data, JB. It will crash if multiple runs are analyzed by the same job.
212 
213  LOG_INFO<<Form("::InitRun(%d) done, W-algo params: trigID: bht3=%d L2W=%d isMC=%d\n TPC: nPileupVert>%d, vertex |Z|<%.1fcm, primEleTrack: nFit>%d, hitFrac>%.2f Rin<%.1fcm, Rout>%.1fcm, PT>%.1fGeV/c\n BTOW ADC: kSigPed=%d AdcThr>%d maxAdc>%.0f clustET>%.1f GeV ET2x2/ET4x4>%0.2f ET2x2/nearTotET>%0.2f\n dist(track-clust)<%.1fcm, nearDelR<%.1f\n Counters Thresholds: track>%.1f GeV, tower>%.1f GeV Use ETOW: flag=%d ScaleFact=%.2f\nBtowScaleFacor=%.2f\n W selection highET>%.1f awayDelPhi<%.1frad ptBalance>%.1fGeV |leptonEta|<%.1f",
214  par_inpRunNo,par_bht3TrgID, par_l2wTrgID,isMC,
215  par_minPileupVert,par_vertexZ,
216  par_nFitPts,par_nHitFrac, par_trackRin, par_trackRout, par_trackPt,
217  par_kSigPed, par_AdcThres,par_maxADC,par_clustET,par_clustFrac24,par_nearTotEtFrac,
218  par_delR3D,par_nearDeltaR,
219  par_countTrPt,par_countTowEt,par_useEtow,par_etowScale,
220  par_btowScale,
221  par_highET,par_awayDeltaPhi,par_ptBalance,par_leptonEta
222  )<<endm;
223 
224  return kStOK;
225 }
226 
227 //________________________________________________
228 //________________________________________________
229 Int_t
230 St2009WMaker::FinishRun(int runNo){
231  LOG_INFO<<Form("::FinishRun(%d)",runNo)<<endm;
232 return kStOK;
233 }
234 
235 //________________________________________________
236 //________________________________________________
237 void
238 St2009WMaker::Clear(const Option_t*){
239  wEve.clear();
240 }
241 
242 
243 //--------------------------------------
244 //--------------------------------------
245 Int_t
246 St2009WMaker::Make(){
247  nInpEve++;
248  wEve.id=mMuDstMaker->muDst()->event()->eventId();
249  const char *afile = mMuDstMaker->GetFile();
250  //printf("inpEve=%d eveID=%d daqFile=%s\n",nInpEve, wEve.id,afile);
251  if(nInpEve%2000==1) printf("\n-----in---- %s, nEve: inp=%d trig=%d accpt=%d daqFile=%s\n", GetName(),nInpEve,nTrigEve, nAccEve,afile);
252  hA[0]->Fill("inp",1.);
253 
254  int btowStat=accessBTOW(); //need btow info for MC trg simu
255 
256  if( accessTrig()) return kStOK; //skip event w/o valid trig ID
257  nTrigEve++;
258 
259  accessBSMD();
260  if( accessVertex()) return kStOK; //skip event w/o ~any reasonable vertex
261 
262  if(mJetReaderMaker) {// just QA plots for jets
263  mJets = getJets(mJetTreeBranch); //get input jet info
264  for (int i_jet=0; i_jet< nJets; ++i_jet){
265  StJet* jet = getJet(i_jet);
266  float jet_pt = jet->Pt();
267  float jet_eta = jet->Eta();
268  float jet_phi = jet->Phi();
269  hA[117]->Fill(jet_eta,jet_phi);
270  hA[118]->Fill(jet_pt);
271  }
272  }
273 
274  if( accessTracks()) return kStOK; //skip event w/o ~any highPt track
275 
276  if(wEve.bemc.tileIn[0]==1)
277  hA[0]->Fill("B-in",1.0);
278  if( btowStat ) return kStOK; //skip event w/o energy in BTOW
279  hA[0]->Fill("B200",1.0);
280 
281  accessETOW();// get energy in ETOW
282 
283  if( extendTrack2Barrel()) return kStOK; //skip if not extends to barrel
284  if( matchTrack2Cluster()) return kStOK; //skip if cluster too soft
285 
286  nAccEve++;
287 
288  /* now it starts to get interesting, process every track
289  on the list till the end. */
290  findNearJet();
291  findAwayJet();
292 
293  if(mJetReaderMaker) findPtBalance();
294 
295  hadronicRecoil();
296 
297  if(mJetReaderMaker) tag_Z_boson();
298 
299  find_W_boson();
300  if(nAccEve<2 ||nAccEve%1000==1 ) wEve.print(0x0,isMC);
301 
302  return kStOK;
303 }
304 
305 
306 //--------------------------------------
307 //--------------------------------------
308 void
309 St2009WMaker::initGeom(){
310 
311  //...... BTOW ...........
312  memset(mapBtowIJ2ID,0,sizeof(mapBtowIJ2ID));
313  for(int softID=1; softID<=mxBtow; softID++) {
314  //........... querry BTOW geom
315  int m,e,s;
316  mBtowGeom->getBin(softID,m,e,s);
317  float etaF,phiF;
318  mBtowGeom->getEta(m,e,etaF);
319  mBtowGeom->getPhi(m,s,phiF); // -pi <= phi < pi
320 
321  int iEta, iPhi;
322  assert(L2algoEtaPhi2IJ(etaF,phiF,iEta, iPhi)==0); // tower must be localized at the known position
323  int IJ=iEta+ iPhi*mxBTetaBin;
324  assert(mapBtowIJ2ID[IJ ]==0); // avoid overlaping mapping
325  mapBtowIJ2ID[IJ ]=softID;
326 
327  Float_t x,y,z;
328  assert( mBtowGeom->getXYZ(softID,x,y,z)==0);
329  positionBtow[softID-1]=TVector3(x,y,z);
330  }// end of loop over towers
331 
332 
333  //...... BSMD-E, -P ...........
334  for(int iep=0;iep<mxBSmd;iep++) {
335  for(int softID=1; softID<=mxBStrips; softID++) {
336  Float_t x,y,z;
337  assert( mBSmdGeom[iep]->getXYZ(softID,x,y,z)==0);
338  positionBsmd[iep][softID-1]=TVector3(x,y,z);
339  }// end of loop over towers
340  }
341 
342 #if 0
343  const Float_t* A=mSmdEGeom->EtaB();
344  for(int i=0;i<mxBetaStrMod+1;i++) {
345  float etaF,phiF;
346  mSmdEGeom->getEta(i+1,etaF);
347  printf("i=%d A=%f str=%f del=%f\n",i,A[i],etaF,A[i]-etaF);
348  }
349 #endif
350 
351  //...... ETOW .............
352  for(int isec=0;isec<mxEtowSec;isec++){
353  for(int isub=0;isub<mxEtowSub;isub++){
354  for(int ieta=0;ieta<mxEtowEta;ieta++){
355  positionEtow[isec*mxEtowSub+isub][ieta]=geomE->getTowerCenter(isec,isub,ieta);
356 #if 0 // trash
357  if(isec==0 && isub==0) {//find radius for each eta ring
358  float x=positionEtow[isec*mxEtowSub+isub][ieta].x();
359  float y=positionEtow[isec*mxEtowSub+isub][ieta].y();
360  etowR[ieta] = x*x+y*y;
361  }
362 #endif
363  }
364  }
365  }
366 
367 }
368 
369 
370 //--------------------------------------
371 //--------------------------------------
372 int // returns error code
373 St2009WMaker::L2algoEtaPhi2IJ(float etaF,float phiF,int &iEta, int &iPhi) {
374  if( phiF<0) phiF+=2*C_PI; // I want phi in [0,2Pi]
375  if(fabs(etaF)>=0.99) return -1;
376  int kEta=1+(int)((etaF+1.)/0.05);
377  iPhi=24-(int)( phiF/C_PI*60.);
378  if(iPhi<0) iPhi+=120;
379  // convention: iPhi=[0,119], kEta=[1,40]
380  iEta=kEta-1;
381  //printf("IJ=%d %d\n",iEta,iPhi);
382  if(iEta<0 || iEta>=mxBTetaBin) return -2;
383  if(iPhi<0 || iPhi>=mxBTphiBin) return -3;
384  return 0;
385 }
386 
387 
388 //________________________________________________
389 //________________________________________________
390 TClonesArray*
391 St2009WMaker::getJets(TString branchName){
392  if(mJetReaderMaker ==0) {
393  nJets=-1; return 0;
394  }
395  assert(mJetReaderMaker->getStJets(branchName)->eventId()==wEve.id);
396  assert(mJetReaderMaker->getStJets(branchName)->runId()==mRunNo);
397  nJets = mJetReaderMaker->getStJets(branchName)->nJets();
398  return mJetReaderMaker->getStJets(branchName)->jets();
399 
400 }
401 
402 
403 // $Log: St2009WMaker.cxx,v $
404 // Revision 1.16 2011/09/14 14:23:20 stevens4
405 // update used for cross section PRD paper
406 //
407 // Revision 1.15 2010/05/04 12:14:35 balewski
408 // runs now w/o jet tree
409 //
410 // Revision 1.14 2010/04/27 16:53:44 stevens4
411 // add code to remove events tagged as Zs from W candidates
412 //
413 // Revision 1.13 2010/04/14 22:23:30 balewski
414 // *** empty log message ***
415 //
416 // Revision 1.12 2010/03/23 15:33:55 seelej
417 // Edit to files to allow the use of a text file for the gains instead of using the DB.
418 //
419 // Revision 1.11 2010/03/14 22:50:31 balewski
420 // *** empty log message ***
421 //
422 // Revision 1.10 2010/02/18 22:34:50 stevens4
423 // add tpc effic study and allow energy scaling for data and MC
424 //
425 // Revision 1.9 2010/01/27 22:12:24 balewski
426 // spin code matched to x-section code
427 //
428 // Revision 1.8 2010/01/23 02:35:38 stevens4
429 // add ability to scale jet et and use real btow peds for rcf mc
430 //
431 // Revision 1.7 2010/01/21 00:15:25 balewski
432 // added sector & run dependent TPC cuts on Rin, Rout
433 //
434 // Revision 1.6 2010/01/18 03:26:15 balewski
435 // expanded TPC track filtering, not finished
436 //
437 // Revision 1.5 2010/01/13 03:34:20 stevens4
438 // give trig emulator access to barrel hits
439 //
440 // Revision 1.4 2010/01/09 00:07:16 stevens4
441 // add jet finder
442 //
443 // Revision 1.3 2010/01/06 19:16:47 stevens4
444 // track cuts now on primary component, cleanup
445 //
446 // Revision 1.2 2009/12/30 19:49:58 balewski
447 // test
448 //
449 // Revision 1.1 2009/11/23 23:00:18 balewski
450 // code moved spin-pool
451 //
452 
453 
454 // $Log: St2009WMaker.cxx,v $
455 // Revision 1.16 2011/09/14 14:23:20 stevens4
456 // update used for cross section PRD paper
457 //
458 // Revision 1.15 2010/05/04 12:14:35 balewski
459 // runs now w/o jet tree
460 //
461 // Revision 1.14 2010/04/27 16:53:44 stevens4
462 // add code to remove events tagged as Zs from W candidates
463 //
464 // Revision 1.13 2010/04/14 22:23:30 balewski
465 // *** empty log message ***
466 //
467 // Revision 1.12 2010/03/23 15:33:55 seelej
468 // Edit to files to allow the use of a text file for the gains instead of using the DB.
469 //
470 // Revision 1.11 2010/03/14 22:50:31 balewski
471 // *** empty log message ***
472 //
473 // Revision 1.10 2010/02/18 22:34:50 stevens4
474 // add tpc effic study and allow energy scaling for data and MC
475 //
476 // Revision 1.9 2010/01/27 22:12:24 balewski
477 // spin code matched to x-section code
478 //
479 // Revision 1.8 2010/01/23 02:35:38 stevens4
480 // add ability to scale jet et and use real btow peds for rcf mc
481 //
482 // Revision 1.7 2010/01/21 00:15:25 balewski
483 // added sector & run dependent TPC cuts on Rin, Rout
484 //
485 // Revision 1.6 2010/01/18 03:26:15 balewski
486 // expanded TPC track filtering, not finished
487 //
488 // Revision 1.5 2010/01/13 03:34:20 stevens4
489 // give trig emulator access to barrel hits
490 //
491 // Revision 1.4 2010/01/09 00:07:16 stevens4
492 // add jet finder
493 //
494 // Revision 1.3 2010/01/06 19:16:47 stevens4
495 // track cuts now on primary component, cleanup
496 //
497 // Revision 1.2 2009/12/30 19:49:58 balewski
498 // test
499 //
500 // Revision 1.1 2009/11/23 23:00:18 balewski
501 // code moved spin-pool
502 //
503 // Revision 1.1 2009/11/23 21:11:18 balewski
504 // start
505 //
EEmcGeomSimple::getTowerCenter
TVector3 getTowerCenter(const UInt_t sec, const UInt_t sub, const UInt_t etabin) const
Definition: EEmcGeomSimple.cxx:104
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
St2009WMaker
muDst based extraction of W-signal from pp500 data from 2009
Definition: St2009WMaker.h:50
StMuDstMaker::GetFile
virtual const char * GetFile() const
Returns name of current input or output file, depending on mode (GetFileName does the same...
Definition: StMuDstMaker.cxx:762
StJets::jets
TClonesArray * jets()
Access to the jets in this event.
Definition: StJets.h:42
St2009WMaker::Make
virtual Int_t Make()
Definition: St2009WMaker.cxx:246
StJets::nJets
int nJets() const
The number of jets found in this event.
Definition: StJets.h:39
StJets::eventId
int eventId() const
access to event numbers, used to synchronize with StMuDstMaker for simultaneous reading ...
Definition: StJets.h:59
StBemcTables::loadTables
void loadTables(StMaker *anyMaker)
load tables.
Definition: StBemcTables.cxx:73
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
EEmcGeomSimple
EEMC simple geometry.
Definition: EEmcGeomSimple.h:22
kStOK
Definition: Stypes.h:40
StEmcGeom::getBin
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Definition: StEmcGeom.h:321
StJet
Definition: StJet.h:91
Generated by   doxygen 1.8.5