StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEEmcIUPi0Analysis.cxx
1 
14 // For real data analyses, you need to switch off #define EFFICIENCY.
15 //#define EFFICIENCY
16 #include "StEEmcIUPi0Analysis.h"
17 #include "StEEmcIUMixMaker.h"
18 #include "StEEmcIUPointMaker.h"
19 #include "TH1F.h"
20 #include "TH2F.h"
21 #include "TH2D.h"
22 #include "TTree.h"
23 #include "TFile.h"
24 #include <stdio.h>
25 #include <iostream>
26 #include "SpinCutsIU.h"
27 #include "StEEmcPool/StEEmcPointMaker/EEmcSectorFit.h"
28 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
29 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
30 #include "StEEmcPool/StRFEmcTrigger/StRFEmcTrigMaker.h"
31 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
32 #include "StMcOutputMaker.h"
33 ClassImp(StEEmcIUPi0Analysis);
34 
35 // ----------------------------------------------------------------------------
36 StEEmcIUPi0Analysis::StEEmcIUPi0Analysis(const Char_t *name):StMaker(name)
37 {
38  mCuts=new SpinCutsIU();
39  mTrigSim = 0;
40  mTrigSimThreshold = 0;
41  mEEgeom=new EEmcGeomSimple();
42  mSpinSort = true;
43 
44 }
45 
46 // ----------------------------------------------------------------------------
47 Int_t StEEmcIUPi0Analysis::Init()
48 {
49 
50  // book histograms for pi0 analysis
51  mHistograms[ kAny ] = new SpinIUHistos("Any", "Pi0 spectra, unsorted");
52  // book histograms for pi0 analysis
53  mBackgrounds[ kAny ] = new SpinIUHistos("AnyB", "combinatoric spectra, unsorted");
54  hPi0Mass=new TH1F("hMass","invariant mass in windows [0.08,0.18]Gev and DeltaR<0.04",120,0.,1.2);
55  hPTvsPiNum = new TH1F("ReconPt", "Recon PT",50,0.,25. );
56  hMcEta = new TH1F("McEta", "Mc Eta",50,0.,2.5 );
57  hEEmcEta = new TH1F("EEmcEta", "EEmc Eta",50,0.,2.5 );
58  hMcPhi = new TH1F("McPhi", "Mc Phi",360,-180.,180. );
59  hMcPt = new TH1F("McPt", "Mc Pt",50,0.,25. );
60  hMcEnergy=new TH1F("McEnergy","Mc energy for recon pi0",80,0.,80.);
61  hREtaPhi=new TH2F("REEmcEtaPhi", "Reconstructed pi0 track; Recon phi / deg;ReconEta",360,-180.,180,20,1.,2.);
62  hReconEta = new TH1F("ReconEta", "Reconstructed Eta",50,0.,2.5 );
63  hReconPhi = new TH1F("ReconPhi", "Recon Phi",360,-180.,180. );
64  hRecoEnergy = new TH1F("ReconEnergy","Reconstructed Pi0 Energy",60,0.,60.);
65  hResoEta = new TH1F("ResoEta", "Resolution of Eta: ReconEta - EEmcEta",10,-1.,1. );
66  hResoPhi = new TH1F("ResoPhi", "Resolution of Phi: ReconPhi - McPhi",20,-20.,20. );
67  hResoPt = new TH1F("ResoPt", "Resolution of Pt",160,-4.,4. );
68  hResoEnergy=new TH1F("ResoEnergy", "Resolution of Energy",200,-10.,10.);
69  hResoPtvsG = new TH2F("ResoPtvsG", "Resolution of Pt; mcPt-ReconPt;mcPt",80,-4.,4.,250,0.,25.);
70  hResoEnergyvsG = new TH2F("ResoEnergyvsG", "Resolution of Energy; mcEnergy-ReconEnergy;mcEnergy",200,-10.,10.,80,0.,80.);
71  hRealMcPD= new TH1F("RealMcdistance","Distance between real and MC pi0",100,0.,0.5);
72  hRealMcPR= new TH1F("RRealMcdistance","Distance between real and MC pi0 in the EEmc detector",100,0.,0.5);
73  mHistograms[1]=new SpinIUHistos("PP","Pi0 spectra, PP, spin4=dec5");
74  mHistograms[2]=new SpinIUHistos("PN","Pi0 spectra, PN, spin4=dec6");
75  mHistograms[3]=new SpinIUHistos("NP","Pi0 spectra, NP, spin4=dec9");
76  mHistograms[4]=new SpinIUHistos("NN","Pi0 spectra, NN, spin4=dec10");
77  EventCounter=new TH1F("EventCounter","Event counter",10,0.,10.);
78  hEventCounter=new TH1F("hEventCounter","Event counts",10,0.,10.);
79  hEventCounter -> GetXaxis() -> SetBinLabel(1,"raw");
80  hEventCounter -> GetXaxis() -> SetBinLabel(2,"minb");
81  hEventCounter -> GetXaxis() -> SetBinLabel(3,"hw trg");
82  hEventCounter -> GetXaxis() -> SetBinLabel(5,"sw trg");
83  hEventCounter -> GetXaxis() -> SetBinLabel(5,"sw trg[3,6]");
84  hEventCounter -> GetXaxis() -> SetBinLabel(6,"vertex");
85  hPairCounter = new TH1F("hPairCounter","Pair counts",10,0.,10.);
86  hPi0Counter = new TH1F("hPi0Counter","Pi0 counts",10,0.,10.);
87  hFillPatternI = new TH1F("hFillPatternI","Intended fill pattern:bunch X-ing @ STAR IP:n runs",120,0.,120.);
88  hSpin4 = new TH1F("hSpin4","Spin 4:spin4 state",11,0.,11.);
89  hSpin4mb = new TH1F("hSpin4mb","Spin 4:spin4 state",11,0.,11.);
90  hBx7 = new TH1F("hBx7","7-bit bunch crossing:bx7",120,0.,120.);
91  hBx48 = new TH1F("hBx48","48-bit bunch crossing:bx48",120,0.,120.);
92  hBx7diffBx48 = new TH2F("hBx7diffBx48","bx1=(bx7-off7)%120, bx2=(bx48-off48)%120;bx7;bx1-bx2",120,0.,120.,21,-10.5,10.5);
93  hBxStar = new TH1F("hBxStar","Beam x-ing at STAR IP;star bunch x-ing",120,0.,120.);
94  hBxStarPi0 = new TH1F("hBxStarPi0","Beam x-ing at STAR IP with pi0 detected;star bunch x-ing",120,0.,120.);
95  hMassBx=new TH2F("hMassBx","Beam x-ing vs Mass;M [GeV];STAR bunch xing",120,0.,1.2,120,0.,120.);
96  hZvertBx=new TH2F("hZvertBx","Beam x-ing vs Z vertex [0.1,0.18]",150,-150.,150.,120,0.,120.);
97  hZggBx=new TH2F("hZggBx","Beam x-ing vs Zgg",100,0.,1.,120,0.,120.);
98  hEtaBx=new TH2F("hEtaBx","Beam x-ing vs eta",120,0.8,2.2,120,0.,120.);
99  DEtaMass=new TH2F("DEtaMass","Delta Eta of the two points vs Mass of the reconstructed pi0",60,0.,0.6,20,0.,0.1);
100  DPhiMass=new TH2F("DPhiMass","Delta Phi of the two points vs Mass of the reconstructed pi0",60,0.,0.6,20,0.,0.1);
101  dUvsdV=new TH2D("dUvsdV","Seed distance between smd clusters; dU; dV",10,0,10,10,0,10);
102  dUvsdVGood=new TH2D("dUvsdVGood","Good Seed distance between smd clusters; dU; dV",10,0,10,10,0,10);
103  dUvsRatio=new TH2F("dUvsRatio","Ucluster ratio along with seed distance;dU;Ratio",20,0.,20.,10,0.,1.);
104  dVvsRatio=new TH2F("dVvsRatio","Vcluster ratio along with seed distance;dV;Ratio",20,0.,20.,10,0.,1.);
105  dUvsRatioGood=new TH2F("dUvsRatioGood","Good Ucluster ratio along with seed distance;dU;Ratio",20,0.,20.,10,0.,1.);
106  dVvsRatioGood=new TH2F("dVvsRatioGood","Good Vcluster ratio along with seed distance;dV;Ratio",20,0.,20.,10,0.,1.);
107  dUVzeroE=new TH1F("dUVzeroE","Pi0 energy with dU or dV=0",60,0.,60);
108  dUVzeroEta=new TH1F("dUVzeroEta","Pi0 Eta with dU or dV=0",20,1.,2.);
109  GoodPiGeo=new TH2F("GoodPiGeo","Pi0 distribution in EEMC with mass range [0.08,0.18]Gev;x;y",120,-240.,240.,120,-240.,240.);
110  GoodPiPt=new TH1F("GoodPiPt", "Recon PT",50,0.,25.);
111  mRealEvent=new StEEmcIUMixEvent();
112  mMixedEvent=new StEEmcIUMixEvent();
113  mRealTree=new TTree("mRealTree","Real Events");
114  mRealTree->Branch("MixEvent","StEEmcIUMixEvent",&mRealEvent);
115  mMixedTree=new TTree("mMixedTree","Mixed Events");
116  mMixedTree->Branch("MixEvent","StEEmcIUMixEvent",&mMixedEvent);
117 
118  AddObj( mRealTree, ".hist" );
119  AddObj( mMixedTree, ".hist" );
120 
121  return StMaker::Init();
122 }
123 // ----------------------------------------------------------------------------
124 Int_t StEEmcIUPi0Analysis::InitRun(Int_t run)
125 {
126 
127  mSpinDb->print(0); // 0=short, 1=huge
128  Info("InitRun","run number = %d",run);
129  const Int_t *spin8bits=mSpinDb->getSpin8bits();
130  for(int bx=0;bx<120;bx++){
131  Bool_t isFilled=(spin8bits[bx] & 0x11)==0x11;
132  if(isFilled) hFillPatternI->Fill(bx);
133  }
134 
135  mRunNumber = run;
136 
137  return StMaker::InitRun(run);
138 
139 }
140 // ----------------------------------------------------------------------------
141 Int_t StEEmcIUPi0Analysis::Make()
142 {
143 
145  EventCounter->Fill(0.,1.);
146  Int_t spin4 = -1;
147  Int_t bxStar = -1;
148  if ( mSpinDb -> isValid() && mSpinDb -> isPolDirLong() )
149  spin4=getSpinState( mMuDst->muDst(), bxStar );
150 
151 
152  hBxStar -> Fill( bxStar );
153 
154  mRealEvent -> setEvent( mMuDst->muDst()->event() );
155  mRealEvent -> setSpin4( spin4 );
156  mMixedEvent -> setEvent( mMuDst->muDst()->event() );
157  mMixedEvent -> setSpin4( spin4 );
158 
159  for ( Int_t ii=0;ii<720;ii++ ) {
160  StEEmcTower tower=mEEanalysis->tower(ii);
161  if ( tower.fail() ) continue;
162  mRealEvent->mADC[ii]=(60./4096.)*tower.adc();
163  mRealEvent->mStat[ii]=tower.stat();
164  }
165 
166 
168 
169  if ( !accept( mMuDst->muDst() ) ) goto END_OF_MAKE;
170  EventCounter->Fill(2.,1.);
171  hPairCounter -> Fill( kEvent );
172 
173 
175  static Int_t mymap[]={0,0,0,0,0,1,2,0,0,3,4};
176 
178  mRealEvent -> mTotalEnergyT = mEEanalysis -> energy(0);
179  mRealEvent -> mTotalEnergyP = mEEanalysis -> energy(1);
180  mRealEvent -> mTotalEnergyQ = mEEanalysis -> energy(2);
181  mRealEvent -> mTotalEnergyR = mEEanalysis -> energy(3);
182  mRealEvent -> mTotalEnergyU = mEEanalysis -> energy(4);
183  mRealEvent -> mTotalEnergyV = mEEanalysis -> energy(5);
184 
186  numcutCan=0;
187  numoacceptedpair=0;
188 
189 
190 
191  for ( Int_t i=0;i<mEEmixer->numberOfCandidates();i++ )
192  {
193 
194  StEEmcIUPair pair = mEEmixer->candidate(i);
195  hPairCounter -> Fill( kPair );
196 
197  if ( !accept( pair ) ) continue;
198 
199  mRealEvent -> addPair( pair );
200  numoacceptedpair++;
201 
202  //std::cout << "spin4=" << spin4 << " ";
203  //pair.print();
204 
205  mHistograms[ kAny ] -> Fill( pair );
206  //assert(mySputMk);
207  StEEmcIUPoint point1=pair.point(0);
208  StEEmcIUPoint point2=pair.point(1);
209  //int sector=point1.sector();
210  TVector3 point1pos=point1.position();
211  TVector3 point2pos=point2.position();
212  float deta=fabs(point1pos.Eta()-point2pos.Eta());
213  float dphi=fabs(point1pos.Phi()-point2pos.Phi());
214  DEtaMass->Fill(pair.mass(),deta);
215  DPhiMass->Fill(pair.mass(),dphi);
216  //calculate detector Eta
217  float Rxy=TMath::Sqrt(pair.vertex().x()*pair.vertex().x()+pair.vertex().y()*pair.vertex().y());
218  float hHeight=pair.pt()*(270.0-pair.vertex().Z())/pair.pz()+Rxy;
219  float etatheta=TMath::ATan(hHeight/270.0);
220  //printf("accept pz=%f\n",pair.pz());
221  float mideta=tan(etatheta/2.0);
222  float eemceta=-log(mideta);
223  hEEmcEta->Fill(eemceta);
224 
225  float EsmallU=(point1.cluster(0).energy()<point2.cluster(0).energy())?point1.cluster(0).energy():point2.cluster(0).energy();
226  float ElargeU=(point1.cluster(0).energy()>point2.cluster(0).energy())?point1.cluster(0).energy():point2.cluster(0).energy();
227  float EsmallV=(point1.cluster(1).energy()<point2.cluster(1).energy())?point1.cluster(1).energy():point2.cluster(1).energy();
228  float ElargeV=(point1.cluster(1).energy()>point2.cluster(1).energy())?point1.cluster(1).energy():point2.cluster(1).energy();
229  float dU=abs(point1.cluster(0).strip(0).index()-point2.cluster(0).strip(0).index());
230  float dV=abs(point1.cluster(1).strip(0).index()-point2.cluster(1).strip(0).index());
231  if(pair.mass()<0.06) {
232  dUvsdV->Fill(dU,dV);
233 
234  if(dU==0||dV==0)
235  {
236  dUVzeroE->Fill(pair.energy());
237  dUVzeroEta->Fill(pair.momentum().Eta());
238  }
239  if(dU==0)
240  {
241  dVvsRatio->Fill(dV,EsmallV/ElargeV);
242  }
243  if(dV==0)
244  {
245  dUvsRatio->Fill(dU,EsmallU/ElargeU);
246  }
247  }//end of mass cut 0.06
248  if(pair.mass()>=0.08 && pair.mass()<=0.18){
249  Float_t e1 = pair.point(0).energy();
250  Float_t e2 = pair.point(1).energy();
251 
252  TVector3 pp = (e1*point1pos + e2*point2pos) * ( 1/(e1+e2) );
253  GoodPiGeo->Fill(pp.X(),pp.Y());
254  GoodPiPt->Fill(pair.pt());
255  dUvsdVGood->Fill(dU,dV);
256 
257  if(dU==0)
258  {
259  dVvsRatioGood->Fill(dV,EsmallV/ElargeV);
260  }
261 
262  if(dV==0)
263  {
264  dUvsRatioGood->Fill(dU,EsmallU/ElargeU);
265  }
266  }//end of good mass cut
267 #ifdef EFFICIENCY
268 
269  //start to calculate pi0 reconstruction efficiency
270  StMcOutputMaker *mkMc=(StMcOutputMaker *)GetMaker("mcRead");
271  assert(mkMc);
272 
273  float rPt=0;
274  float rEta=0;
275  float rPhi=0;
276  float rE=0;
277 
278  if(pair.mass()>=0.08 && pair.mass()<=0.18)
279  {
280 
281  numcutCan++;
282  rPt=pair.pt();
283  rEta=pair.momentum().Eta();
284  rPhi=pair.momentum().Phi();
285  rE=pair.energy();
286  int size=0;
287  size=mkMc->gTr.size();
288  float minD=10000000;
289  float minR=10000000;
290  //printf("ana size=%d\n",size);
291  for ( Int_t j=0;j<size;j++ )
292  {
293  StMcTrack *gt=mkMc->gTr[j];
294  float eemceta=mkMc->geemcEta[j];
295  StLorentzVectorF p4=gt->fourMomentum();
296  float etaMc=p4.pseudoRapidity();
297  float phiMc=p4.phi();
298  //printf("rPhi=%f phiMc=%f rEta=%f etaMc=%f\n",rPhi,phiMc,rEta,etaMc);
299 
300  float deltaphi=fabs(rPhi-phiMc);
301  if(deltaphi>=3.14159265)
302  {deltaphi=6.283-deltaphi;}
303  float deltaeta=rEta-etaMc;
304  float Rdeltaeta=rEta-eemceta;
305  float deltad=TMath::Sqrt(deltaphi*deltaphi+deltaeta*deltaeta);
306  float deltaR=TMath::Sqrt(deltaphi*deltaphi+Rdeltaeta*Rdeltaeta);
307  if(deltad<=minD)
308  {
309  //printf("dphi=%f deta=%f\n",deltaphi,deltaeta);
310  minD=deltad;
311  }
312  if(deltaR<=minR)
313  {
314  minR=deltaR;
315  }
316  }
317  hRealMcPD->Fill(minD);
318  hRealMcPR->Fill(minR);
319  for ( Int_t j=0;j<size;j++ )
320  {
321  StMcTrack *gt=mkMc->gTr[j];
322  float eemceta=mkMc->geemcEta[j];
323  //float eeeta=mkMc->geeEta[j];
324  StLorentzVectorF p4=gt->fourMomentum();
325  float etaMc=p4.pseudoRapidity();
326  //printf("etaMc=%f eemceta=%f\n",etaMc,eemceta);
327  float phiMc=p4.phi();
328  //printf("rPhi=%f phiMc=%f rEta=%f etaMc=%f\n",rPhi,phiMc,rEta,etaMc);
329  float ptMc=p4.perp();
330  float energyMc=p4.e();
331  //printf("energyMc=%f\n",energyMc);
332  float deltaphi=fabs(rPhi-phiMc);
333  if(deltaphi>=3.14159265)
334  {deltaphi=6.283-deltaphi;}
335  float deltaeta=rEta-etaMc;
336  float Rdeltaeta=rEta-eemceta;
337  float deltad=TMath::Sqrt(deltaphi*deltaphi+deltaeta*deltaeta);
338  float deltaR=TMath::Sqrt(deltaphi*deltaphi+Rdeltaeta*Rdeltaeta);
339 
340 
341 
342  //deltad for pythia deltaR for simple MC
343  if(deltaR<=0.04)
344  {
345  hPi0Mass->Fill(pair.mass());
346  hPTvsPiNum->Fill(rPt);
347  hREtaPhi->Fill(rPhi/3.1416*180,rEta);
348  hReconEta->Fill(rEta);
349  hReconPhi->Fill(rPhi/3.1416*180);
350  hRecoEnergy->Fill(rE);
351 
352  hMcEta->Fill(etaMc);
353  hEEmcEta->Fill(eemceta);
354  hMcPhi->Fill(phiMc/3.1416*180);
355  hMcPt->Fill(ptMc);
356  hResoPt->Fill(rPt-ptMc);
357  hResoPtvsG->Fill(rPt-ptMc,ptMc);
358  hResoEta->Fill(rEta-etaMc);
359  hResoPhi->Fill((rPhi-phiMc)/3.1416*180);
360  hMcEnergy->Fill(energyMc);
361  hResoEnergyvsG->Fill(rE-energyMc,energyMc);
362  hResoEnergy->Fill(rE-energyMc);
363  //printf("rE=%f mcE=%f\n",rE,energyMc);
364  //hZgg->Fill(pair.zgg());
365 
366 
367  //printf("zvertex=%f\n",pair.vertex().Z());
368  //McEvsMass->Fill(energyMc,pair.mass());
369  }
370 
371  }
372 
373  }//end of efficiency calculation
374 
375 #endif
376 
378  if ( !mSpinSort ) {
379  std::cout << "Problem detected, spin sorting disabled" << std::endl;
380  continue;
381  }
382 
383  spin4=getSpinState( mMuDst->muDst(),bxStar );
384  if ( spin4>=0 )
385  if ( mymap[spin4] ) {
386  mHistograms[ mymap[spin4] ] -> Fill( pair );
387  }
388 
389  hMassBx->Fill(pair.mass(),bxStar);
390  if ( pair.mass()>0.1 && pair.mass() < 0.18 ) {
391  hBxStarPi0->Fill(bxStar);
392  hZvertBx->Fill(pair.vertex().Z(),bxStar);
393  hZggBx->Fill(pair.zgg(),bxStar);
394  hEtaBx->Fill(pair.momentum().Eta(),bxStar);
395  }
396 
397 
398  }
399  hPi0Counter->Fill(numoacceptedpair);
400  //printf("numcutCan=%d\n",numcutCan);
401  //}//end of mix candidate 6
403  for ( Int_t i=0;i<mEEmixer->numberOfMixedCandidates();i++ )
404  {
405 
406  StEEmcIUPair pair = mEEmixer->mixedCandidate(i);
407  if ( !accept( pair, false ) ) continue;
408  mBackgrounds[ kAny ] -> Fill( pair );
409  mMixedEvent -> addPair( pair );
410 
411  }
412 
413  END_OF_MAKE:
414  mRealTree->Fill();
415  mMixedTree->Fill();
416  return kStOK;
417 }
418 
419 // ----------------------------------------------------------------------------
420 void StEEmcIUPi0Analysis::mixer(const Char_t *name)
421 {
422  mEEmixer=(StEEmcIUMixMaker *)GetMaker(name);
423  assert(mEEmixer);// please specify a valid mix maker
424 }
425 
426 // ----------------------------------------------------------------------------
427 void StEEmcIUPi0Analysis::points(const Char_t *name)
428 {
429  mEEpoints=(StEEmcIUPointMaker *)GetMaker(name);
430  assert(mEEpoints);
431 }
432 
433 // ----------------------------------------------------------------------------
434 void StEEmcIUPi0Analysis::mudst(const Char_t *name)
435 {
436  mMuDst=(StMuDstMaker*)GetMaker(name);
437  assert(mMuDst);
438 }
439 // ----------------------------------------------------------------------------
440 void StEEmcIUPi0Analysis::analysis(const Char_t *name)
441 {
442  mEEanalysis=(StEEmcA2EMaker*)GetMaker(name);
443  assert(mEEanalysis);
444 }
445 // ----------------------------------------------------------------------------
446 void StEEmcIUPi0Analysis::spin(const Char_t *name)
447 {
448  mSpinDb=(StSpinDbMaker*)GetMaker(name);
450 }
451 // ----------------------------------------------------------------------------
452 Bool_t StEEmcIUPi0Analysis::twoBodyCut( StEEmcIUPair &pair )
453 {
454 
455  static const Int_t maxPerCluster = 4;
456  static const Float_t minTowerFrac = 0.;
457 
459 
460  Bool_t towers[720]; for (Int_t i=0;i<720;i++ ) towers[i]=false;
461  StEEmcTower t1 = pair.point(0).tower(0);
462  StEEmcTower t2 = pair.point(1).tower(0);
463 
464  Float_t Etowers=0.;
465  Etowers += t1.energy();
466 
467 
468  towers[ t1.index() ] = true;
469  towers[ t2.index() ] = true;
470 
471  for ( Int_t i=0;i<t1.numberOfNeighbors();i++ ) {
472  StEEmcTower mytow=t1.neighbor(i);
473  towers[ mytow.index() ] = true;
474  Etowers += mytow.energy();
475  }
476  for ( Int_t i=0;i<t2.numberOfNeighbors();i++ ) {
477  StEEmcTower mytow=t2.neighbor(i);
478  if ( !towers[mytow.index()] ) Etowers += mytow.energy();
479  towers[ mytow.index() ] = true;
480  }
481 
482 
483 
485 
488 
489  Int_t count = 0;
490  Float_t pe1=pair.point(0).energy();
491  Float_t pe2=pair.point(1).energy();
492  Float_t min2=(pe1<=pe2)?pe1:pe2;
493  Bool_t Most2E=true;
494 
495  for ( Int_t i=0;i<mEEpoints->numberOfPoints();i++ )
496  {
497  StEEmcIUPoint p=mEEpoints->point(i);
498  Float_t ppe=p.energy();
499  StEEmcTower t=p.tower(0);
500  if ( towers[ t.index() ] ){
501  count++;
502  if(t.index()!=t1.index() && t.index()!=t2.index() && ppe>min2){
503  Most2E=false;
504  }
505  }
506  }
507  //off most energetic two points
508  Most2E=true;
509 
510  Float_t Epoints = pair.energy();
511 
512  return ( count <= maxPerCluster && Most2E && Epoints > minTowerFrac * Etowers );
513 
514 }
515 
516 // ----------------------------------------------------------------------------
517 // mudst based cuts
518 Bool_t StEEmcIUPi0Analysis::accept( StMuDst *mudst )
519 {
520 
521  // verify that the trigger is in the trigger list
522  StMuEvent *event=mudst->event();
523  assert(event);
524 
525  Bool_t good_trigger = false;
526 
527  // -------------------------------------------------------
528  //
529  // Trigger from real data
530  //
531 
532  StMuTriggerIdCollection tic = event -> triggerIdCollection();
533  StTriggerId l1trig = tic.l1();
534  // if no triggers are specified, always return true
535  if ( mTriggerList.size() <=0 ) {
536  good_trigger = true;
537  }
538 
539  Int_t go=0;
540  std::vector<Int_t>::iterator iter=mTriggerList.begin();
541  while ( iter != mTriggerList.end() ) {
542  go = l1trig.isTrigger( (*iter) );
543  good_trigger |= go;
544  iter++;
545  }
546 
547 
548  // -------------------------------------------------------
549  //
550  // Trigger from sim data
551  //
552  if ( mTrigSim )
553  {
554  good_trigger = mTrigSim -> getEEMCtrigHT( mTrigSimThreshold );
555  good_trigger &= mTrigSim -> getBBCtrig();
556  }
557 
558 
559  if ( l1trig.isTrigger( mMinBias ) )
560  hEventCounter -> Fill( kMinBias );
561 
562  // must have a valid trigger to proceed
563  if ( good_trigger )
564  {
565 
566  hEventCounter -> Fill( kTrig );
567  }
568  else
569  {
570 
571  return false;
572  }
573 
574  // loop over all hit towers and verify that one tower
575  // exceeds the software threshold in mCuts
576  StEEmcTower ht;
577  Bool_t sw36=0;
578 
579  for ( Int_t i=0;i<mEEanalysis->numberOfHitTowers(0);i++ )
580  {
581  StEEmcTower t=mEEanalysis->hittower(i,0);
582  if ( !t.fail() && t.et() > ht.et() ) {
583  ht=t;
584  if ( ht.et()>=mCuts->tower_et_cut &&
585  ht.sector()>=2&&ht.sector()<=5 ) sw36=true;
586  }
587  }
588  //printf("ht.et=%f\n",ht.et());
589  if ( ht.et() < mCuts->tower_et_cut ) {
590  good_trigger = false;
591  //printf("false\n");
592  return false;
593  }
594  else
595  {
596  hEventCounter -> Fill( kSoftTrig );
597  }
598  if ( sw36 ) hEventCounter -> Fill( kSoftTrig36 );
599 
600 
601  StThreeVectorF vert=event->primaryVertexPosition();
602  if ( vert.z() < mCuts->z_vertex_min) {
603  // std::cout << "No vertex" << std::endl;
604  return false;
605  }
606  if ( vert.z() > mCuts->z_vertex_max ) {
607  // std::cout << "No vertex" << std::endl;
608  return false;
609  }
610 
611  hEventCounter -> Fill( kVertex );
612 
613 
614  return true;
615 
616 }
617 
618 Bool_t StEEmcIUPi0Analysis::accept( StEEmcIUPair pair, Bool_t fill )
619 {
620 
622  if ( !twoBodyCut(pair) ) return false;
623  if (fill) hPairCounter -> Fill( kTwoBody );
625  StEEmcIUPoint point1=pair.point(0);
626  StEEmcIUPoint point2=pair.point(1);
627  StEEmcTower tower1=pair.point(0).tower(0);
628  StEEmcTower tower2=pair.point(1).tower(0);
629  TVector3 pointpos1=point1.position();
630  TVector3 pointpos2=point2.position();
631  Float_t peta1=pointpos1.Eta();
632  Float_t peta2=pointpos2.Eta();
633  Float_t pphi1=pointpos1.Phi();
634  Float_t pphi2=pointpos2.Phi();
635  Int_t mod_phi1=int(pphi1*180./3.14159265+180.)%6;
636  Int_t mod_phi2=int(pphi2*180./3.14159265+180.)%6;
637  pointpos1=pointpos1.Unit();
638  pointpos2=pointpos2.Unit();
639 
640  Float_t pet1 = (point1.energy()*pointpos1).Perp();
641  Float_t pet2 = (point2.energy()*pointpos2).Perp();
642 
643  //require both points to be above 1.5 GeV because of fluctuation
644  if(pet1<=1.5 || pet2<=1.5) return false;
645 
646  if ( tower1.adc() < mCuts->adc_cut && tower2.adc() < mCuts->adc_cut ) return false;
647 
648  if ( fill) hPairCounter -> Fill( kAdcTow );
649 
650  /*
651  if ( tower1.et() < mCuts->tower_et_cut &&
652  tower2.et() < mCuts->tower_et_cut ) return false;
653  */
654  TVector3 d1 = mEEgeom -> getTowerCenter( (UInt_t)tower1.sector(), (UInt_t)tower1.subsector(), (UInt_t)tower1.etabin() );
655  TVector3 d2 = mEEgeom -> getTowerCenter( (UInt_t)tower2.sector(), (UInt_t)tower2.subsector(), (UInt_t)tower2.etabin() );
656 
657  TVector3 dd1(d1.x(),d1.y(),d1.z()-pair.vertex().z());
658  TVector3 dd2(d2.x(),d2.y(),d2.z()-pair.vertex().z());
659 
660 
661  d1=d1.Unit();
662  d2=d2.Unit();
663  dd1=dd1.Unit();
664  dd2=dd2.Unit();
665  //Float_t et1 = (tower1.energy()*d1).Perp();
666  //Float_t et2 = (tower2.energy()*d2).Perp();
667  Float_t et11 = (tower1.energy()*dd1).Perp();
668  Float_t et22 = (tower2.energy()*dd2).Perp();
669 
670  //place setting bivariate Gaussian tower et cut according to point positions
671  Float_t ueta1,ueta2,tower_et_cut1=3,tower_et_cut2=3;
672  if(peta1>=1.901 && peta1<=2) ueta1=1.932;
673  if(peta1>=1.807 && peta1<1.901) ueta1=1.845;
674  if(peta1>=1.718 && peta1<1.807) ueta1=1.761;
675  if(peta1>=1.633 && peta1<1.718) ueta1=1.664;
676  if(peta1>=1.552 && peta1<1.633) ueta1=1.585;
677  if(peta1>=1.476 && peta1<1.552) ueta1=1.507;
678  if(peta1>=1.403 && peta1<1.476) ueta1=1.4374;
679  if(peta1>=1.334 && peta1<1.403) ueta1=1.36;
680  if(peta1>=1.268 && peta1<1.334) ueta1=1.294;
681  if(peta1>=1.205 && peta1<1.268) ueta1=1.2315;
682  if(peta1>=1.146 && peta1<1.205) ueta1=1.174;
683  if(peta1>=1.089 && peta1<1.146) ueta1=1.1265;
684 
685  if(peta2>=1.901 && peta2<=2) ueta2=1.932;
686  if(peta2>=1.807 && peta2<1.901) ueta2=1.845;
687  if(peta2>=1.718 && peta2<1.807) ueta2=1.761;
688  if(peta2>=1.633 && peta2<1.718) ueta2=1.664;
689  if(peta2>=1.552 && peta2<1.633) ueta2=1.585;
690  if(peta2>=1.476 && peta2<1.552) ueta2=1.507;
691  if(peta2>=1.403 && peta2<1.476) ueta2=1.4374;
692  if(peta2>=1.334 && peta2<1.403) ueta2=1.36;
693  if(peta2>=1.268 && peta2<1.334) ueta2=1.294;
694  if(peta2>=1.205 && peta2<1.268) ueta2=1.2315;
695  if(peta2>=1.146 && peta2<1.205) ueta2=1.174;
696  if(peta2>=1.089 && peta2<1.146) ueta2=1.1265;
697  float toweretcut=mCuts->tower_et_cut;
698  tower_et_cut1=toweretcut*(exp(-0.5*((peta1-ueta1)*(peta1-ueta1)/pow(0.035,2)+(mod_phi1-0.0)*(mod_phi1-0.0)/pow(2.3,2))));
699  tower_et_cut2=toweretcut*(exp(-0.5*((peta2-ueta2)*(peta2-ueta2)/pow(0.035,2)+(mod_phi2-0.0)*(mod_phi2-0.0)/pow(2.3,2))));
700 
701  if ( et11 < tower_et_cut1 && et22 < tower_et_cut2 ) return false;
702 
703 
704 
705  if ( fill ) hPairCounter -> Fill( kEtTow );
706 
708 
709  float Rxy=TMath::Sqrt(pair.vertex().x()*pair.vertex().x()+pair.vertex().y()*pair.vertex().y());
710  float hHeight=pair.pt()*(270.0-pair.vertex().Z())/pair.pz()+Rxy;
711  float etatheta=TMath::ATan(hHeight/270.0);
712  //printf("accept pz=%f\n",pair.pz());
713  float mideta=tan(etatheta/2.0);
714  float eta=-log(mideta);
715  if ( eta < mCuts->eta_min || eta > mCuts->eta_max ) return false;
716  if ( fill ) hPairCounter->Fill( kKinematics );
717 
718 
719 
720 
721 
722 
723  return true;
724 
725 }
726 
727 // ----------------------------------------------------------------------------
728 Int_t StEEmcIUPi0Analysis::getSpinState( StMuDst *mudst, Int_t &bxStar )
729 {
730 
731  StMuEvent *event = mudst -> event();
732  StL0Trigger *trig =&(event->l0Trigger());
733 
734  StMuTriggerIdCollection tic = event -> triggerIdCollection();
735  StTriggerId l1trig = tic.l1();
736 
737 
738  Int_t bx48 = trig->bunchCrossingId();
739  Int_t bx7 = trig->bunchCrossingId7bit( mRunNumber );
740 
741  //bxStar = mSpinDb->BXstarUsingBX48(bx48);
742  bxStar = mSpinDb->BXyellowUsingBX48(bx48);
743 
744  mRealEvent->bx48 = bx48;
745  mRealEvent->bx7 = bx7;
746  mRealEvent->bxStar = bxStar;
747 
748  mMixedEvent->bx48=bx48;
749  mMixedEvent->bx7 =bx7;
750  mMixedEvent->bxStar=bxStar;
751 
754  if ( bx7 == 0 || bx7 == 119 ) return -1; // kick 0 and 119 out of mix
755  if ( bx7 == 14 || bx7 == 54 ) return -1;
757 
758  hBx7->Fill( bx7 );
759  hBx48->Fill( bx48 );
760  hBx7diffBx48->Fill( bx7, mSpinDb->offsetBX48minusBX7(bx48,bx7) );
761 
762  //$$$::cout << "bx7=" << bx7 << " bx48=" << bx48 << std::endl;
763 
764  if ( mSpinDb -> isMaskedUsingBX48(bx48) ) return -1; // return an error flag
765  if ( mSpinDb->offsetBX48minusBX7(bx48,bx7)!=0 ) std::cout << "BUNCH CROSSINGS INCONSISTENT" << std::endl;
766 
767  // assert(mSpinDb->offsetBX48minusBX7(bx48,bx7)==0); // scaler boards were not in sync
768  // mSpinSort = (mSpinDb->offsetBX48minusBX7(bx48,bx7)==0);
769  // disable spin sorting and clear histograms
770  if ( mSpinDb->offsetBX48minusBX7(bx48,bx7)!=0 )
771  {
772  mSpinSort = false;
773  for ( Int_t ii=1;ii<=4;ii++ ) mHistograms[ii]->Clear();
774  }
775 
776 
777  Int_t spin4 = mSpinDb->spin4usingBX48(bx48);
778  hSpin4->Fill(spin4);
779 
780  if ( l1trig.isTrigger(96011) ) hSpin4mb -> Fill(spin4);
781 
782 
783  return spin4;
784 
785 }
786 // ----------------------------------------------------------------------------
787 void StEEmcIUPi0Analysis::Clear(Option_t *opts)
788 {
789  mRealEvent -> Clear();
790  mMixedEvent -> Clear();
791 }
StEEmcIUPi0Analysis::spin
void spin(const Char_t *name)
specifies the name of the spin db maker
Definition: StEEmcIUPi0Analysis.cxx:446
StEEmcIUPi0Analysis::mEEpoints
StEEmcIUPointMaker * mEEpoints
pointer to the point maker
Definition: StEEmcIUPi0Analysis.h:88
StEEmcIUPi0Analysis::points
void points(const Char_t *name)
specifies the name of the point maker
Definition: StEEmcIUPi0Analysis.cxx:427
StEEmcIUMixMaker::candidate
StEEmcIUPair candidate(Int_t c)
Return a specified candidate pair.
Definition: StEEmcIUMixMaker.h:45
StEEmcTower::numberOfNeighbors
Int_t numberOfNeighbors() const
get the number of neighboring towers
Definition: StEEmcTower.h:54
StEEmcA2EMaker
EEmc ADC –&gt; energy maker.
Definition: StEEmcA2EMaker.h:19
StEEmcElement::stat
void stat(unsigned s)
Set a status bit for this element.
Definition: StEEmcElement.h:23
StEEmcIUPi0Analysis::mEEanalysis
StEEmcA2EMaker * mEEanalysis
pointer to analysis maker
Definition: StEEmcIUPi0Analysis.h:92
StEEmcIUMixMaker
A class for mixing pi0 candidates.
Definition: StEEmcIUMixMaker.h:22
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StEEmcIUMixMaker::numberOfMixedCandidates
Int_t numberOfMixedCandidates()
returns the number of mixed-background candidates
Definition: StEEmcIUMixMaker.h:206
StEEmcIUPair::pt
Float_t pt()
Returns pt of pair.
Definition: StEEmcIUPair.h:78
StEEmcA2EMaker::hittower
StEEmcTower & hittower(Int_t hit, Int_t layer)
Definition: StEEmcA2EMaker.h:63
StEEmcIUPi0Analysis::mSpinDb
StSpinDbMaker * mSpinDb
pointer to the spin database
Definition: StEEmcIUPi0Analysis.h:94
StEEmcIUPair::vertex
TVector3 vertex()
Returns vertex of pair.
Definition: StEEmcIUPair.h:81
StEEmcTower::neighbor
void neighbor(StEEmcTower *n)
add a tower to list of neighbors
Definition: StEEmcTower.h:52
StEEmcIUPi0Analysis::mHistograms
SpinIUHistos * mHistograms[5]
Spin-sorted pi0 histograms.
Definition: StEEmcIUPi0Analysis.h:102
StEEmcIUPoint::energy
void energy(Float_t e, Int_t layer=0)
Set the energy of this point.
Definition: StEEmcIUPoint.h:25
StEEmcIUPi0Analysis::analysis
void analysis(const Char_t *name)
specifies the name of the analysis maker
Definition: StEEmcIUPi0Analysis.cxx:440
StEEmcIUPi0Analysis::accept
Bool_t accept(StMuDst *mu)
method to cut events
Definition: StEEmcIUPi0Analysis.cxx:518
StEEmcIUPi0Analysis
A maker for creating pi0 histograms \Weihong He The StEEmcIUPi0Analysis takes as input the list of pi...
Definition: StEEmcIUPi0Analysis.h:35
StEEmcIUPoint::position
void position(TVector3 p)
Set the position of this point at the SMD plane.
Definition: StEEmcIUPoint.h:23
StEEmcIUPi0Analysis::mudst
void mudst(const Char_t *name)
specifies the name of the mudst maker
Definition: StEEmcIUPi0Analysis.cxx:434
StSpinDbMaker::spin4usingBX48
int spin4usingBX48(int bx48)
8bit spin information
Definition: StSpinDbMaker.cxx:271
StEEmcIUPointMaker
Class for building points from smd clusters.
Definition: StEEmcIUPointMaker.h:21
StEEmcIUPi0Analysis::hPairCounter
TH1F * hPairCounter
histogram to keep track of where candidates get cut
Definition: StEEmcIUPi0Analysis.h:178
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StEEmcIUPi0Analysis::Make
Int_t Make()
processes a single event
Definition: StEEmcIUPi0Analysis.cxx:141
StEEmcIUPointMaker::numberOfPoints
Int_t numberOfPoints()
Return number of points.
Definition: StEEmcIUPointMaker.h:183
StEEmcIUPi0Analysis::InitRun
Int_t InitRun(Int_t)
init run
Definition: StEEmcIUPi0Analysis.cxx:124
StEEmcIUMixMaker::numberOfCandidates
Int_t numberOfCandidates()
returns the number of candidates
Definition: StEEmcIUMixMaker.h:205
StEEmcIUPi0Analysis::Clear
void Clear(Option_t *opts="")
clears the maker
Definition: StEEmcIUPi0Analysis.cxx:787
StEEmcA2EMaker::numberOfHitTowers
Int_t numberOfHitTowers(Int_t layer) const
Definition: StEEmcA2EMaker.h:58
StEEmcTower::etabin
Int_t etabin() const
Returns the etabin of this tower, pre- or postshower element.
Definition: StEEmcTower.h:45
StEEmcIUPi0Analysis::getSpinState
Int_t getSpinState(StMuDst *mu, Int_t &bxs)
method to retrieve 4bit spin state
Definition: StEEmcIUPi0Analysis.cxx:728
StEEmcElement::fail
void fail(unsigned f)
Set a fail bit for this element.
Definition: StEEmcElement.h:25
StEEmcIUPi0Analysis::hEventCounter
TH1F * hEventCounter
histogram to keep track of where events get cut
Definition: StEEmcIUPi0Analysis.h:175
StEEmcIUPi0Analysis::mTrigSim
StRFEmcTrigMaker * mTrigSim
Trigger simulation for MC.
Definition: StEEmcIUPi0Analysis.h:191
StEEmcIUPi0Analysis::mixer
void mixer(const Char_t *name)
Definition: StEEmcIUPi0Analysis.cxx:420
StEEmcTower::subsector
Int_t subsector() const
Returns subsector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:43
StEEmcIUPi0Analysis::mMuDst
StMuDstMaker * mMuDst
pointer to mudst
Definition: StEEmcIUPi0Analysis.h:90
StEEmcIUPoint::cluster
void cluster(StEEmcIUSmdCluster c, Int_t plane)
Add an smd cluster to this point.
Definition: StEEmcIUPoint.h:31
StEEmcTower::index
void index(Int_t i)
Definition: StEEmcTower.cxx:76
StEEmcTower
Base class for representing tower, preshower and postshower elements.
Definition: StEEmcTower.h:11
StEEmcIUPair::mass
Float_t mass()
Returns invariant mass of pair.
Definition: StEEmcIUPair.h:74
StSpinDbMaker::print
void print(int level=0)
vs. STAR==yellow bXing
Definition: StSpinDbMaker.cxx:464
StEEmcTower::sector
Int_t sector() const
Returns sector of this tower, pre- or postshower element.
Definition: StEEmcTower.h:41
StEEmcIUPair::energy
Float_t energy()
Returns energy of pair.
Definition: StEEmcIUPair.h:75
StSpinDbMaker::BXyellowUsingBX48
int BXyellowUsingBX48(int bx48)
4bit spin information
Definition: StSpinDbMaker.cxx:341
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
StEEmcIUPointMaker::point
StEEmcIUPoint point(Int_t ipoint)
Return a specified point.
Definition: StEEmcIUPointMaker.h:184
EEmcGeomSimple
EEMC simple geometry.
Definition: EEmcGeomSimple.h:22
kStOK
Definition: Stypes.h:40
StEEmcIUPoint
Base class for representing EEMC points.
Definition: StEEmcIUPoint.h:12
StEEmcIUPi0Analysis::StEEmcIUPi0Analysis
StEEmcIUPi0Analysis(const Char_t *name)
constructor
Definition: StEEmcIUPi0Analysis.cxx:36
StEEmcIUPi0Analysis::mEEmixer
StEEmcIUMixMaker * mEEmixer
Pointer to the pi0 mixer.
Definition: StEEmcIUPi0Analysis.h:86
StEEmcTower::et
void et(Float_t e)
Definition: StEEmcTower.h:34
StEEmcElement::adc
void adc(Float_t a)
Set the pedestal-subtracted ADC for this element.
Definition: StEEmcElement.h:19
StEEmcIUPair::point
StEEmcIUPoint point(Int_t index)
Definition: StEEmcIUPair.h:30
StEEmcIUPoint::tower
void tower(StEEmcTower t, Float_t w=1.)
Add a tower with specified weight to the point.
Definition: StEEmcIUPoint.h:29
StEEmcIUPi0Analysis::mEEgeom
EEmcGeomSimple * mEEgeom
EEMC tower geometry.
Definition: StEEmcIUPi0Analysis.h:199
StEEmcA2EMaker::tower
StEEmcTower & tower(Int_t index, Int_t layer=0)
Definition: StEEmcA2EMaker.h:69
StSpinDbMaker::getSpin8bits
const int * getSpin8bits()
experts only
Definition: StSpinDbMaker.cxx:424
SpinIUHistos
Spin sorted pi0 histograms.
Definition: SpinIUHistos.h:22
StEEmcIUPi0Analysis::twoBodyCut
Bool_t twoBodyCut(StEEmcIUPair &p)
Definition: StEEmcIUPi0Analysis.cxx:452
StEEmcIUMixMaker::mixedCandidate
StEEmcIUPair mixedCandidate(Int_t m)
Returns the specified mixed candidate pair.
Definition: StEEmcIUMixMaker.h:47
StEEmcIUPi0Analysis::Init
Int_t Init()
initializes the maker
Definition: StEEmcIUPi0Analysis.cxx:47
StEEmcElement::energy
void energy(Float_t e)
Set the energy (adc-ped+0.5)/gain for this element.
Definition: StEEmcElement.h:21
StEEmcIUPair::momentum
TVector3 momentum()
Returns momentum of pair.
Definition: StEEmcIUPair.h:80
StMuTriggerIdCollection
Collection of trigger ids as stored in MuDst.
Definition: StMuTriggerIdCollection.h:31
StEEmcIUPair::zgg
Float_t zgg()
Returns energy-sharing of pair.
Definition: StEEmcIUPair.h:76
StEEmcIUPi0Analysis::mTriggerList
std::vector< Int_t > mTriggerList
Definition: StEEmcIUPi0Analysis.h:96
Generated by   doxygen 1.8.5