00001
#include "StEEmcPi0Maker.h"
00002
#include "StMessMgr.h"
00003
00004
#include "StEEmcPair.h"
00005
00006
#include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00007
#include "StMuDSTMaker/COMMON/StMuDst.h"
00008
#include "StMuDSTMaker/COMMON/StMuEvent.h"
00009
#include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00010
#include "StMuDSTMaker/COMMON/StMuTriggerIdCollection.h"
00011
00012
#include "StEvent/StEvent.h"
00013
#include "StEvent/StTriggerIdCollection.h"
00014
#include "StEvent/StTriggerId.h"
00015
00016
#include "StarClassLibrary/StThreeVectorF.hh"
00017
00018
#include "TH1F.h"
00019
#include "TH2F.h"
00020
00021
#include "TFile.h"
00022
#include "TTree.h"
00023
#include "StEEmcMixEvent.h"
00024
00025 ClassImp(StEEmcPi0Maker);
00026
00027
00028
00029
00030
00031 StEEmcPi0Maker::StEEmcPi0Maker(
const Char_t *name,
00032
StEEmcA2EMaker *a2e,
00033
StEEmcGenericClusterMaker *cl,
00034
StEEmcGenericPointMaker *pt ) :
StMaker(name)
00035 {
00036
mEEanalysis=a2e; assert(a2e);
00037 mEEclusters=cl; assert(cl);
00038
mEEpoints=pt; assert(pt);
00039 mCheckTrigger=
false;
00040 setFile(0);
00041 mTree=0;
00042 }
00043
00044
void StEEmcPi0Maker::setFile( TFile *f ){ mFile=f; }
00045 TTree *StEEmcPi0Maker::tree(){
return mTree; }
00046
00047
00048 Int_t StEEmcPi0Maker::Init()
00049 {
00050 hMass =
new TH2F(
"hMass",
"Diphoton invariant mass; M [GeV]",120,0.,1.2,20,0.,20.);
00051 hMass_cut =
new TH2F(
"hMass_cut",
"Diphoton invariant mass; M [GeV]",120,0.,1.2,20,0.,20.);
00052 hMass_split =
new TH2F(
"hMass_split",
"Mass for events w/ a split smd cluster",120,0.,1.2,20,0.,20.);
00053
00054 hPT =
new TH1F(
"hPT",
"Diphoton p_{T}; p_{T} [GeV]",20,0.,20.);
00055 hPT_cut =
new TH1F(
"hPT_cut",
"Diphoton p_{T}; p_{T} [GeV]",20,0.,20.);
00056 hXF =
new TH1F(
"hXF",
"x_{F}=#frac{p_{L}}{#sqrt{s}/2}",50,0.,1.);
00057 hEnergy =
new TH1F(
"hEnergy",
"Diphoton energy; E [GeV]",50,0.,50.);
00058 hEta =
new TH1F(
"hEta",
"Diphoton #eta",24,0.8,2.2);
00059 hPhi =
new TH1F(
"hPhi",
"Diphoton #phi",30,-3.142,3.142);
00060 hZvertex=
new TH1F(
"hZvertex",
"z vertex",150,-150.,150.);
00061 hZgg =
new TH1F(
"hZgg",
"energy sharing",50,0.,1.);
00062 hZgg_cut =
new TH1F(
"hZgg_cut",
"energy sharing",50,0.,1.);
00063
00064 hEChi2 =
new TH1F(
"hEChi2",
"#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00065 hE1Chi2 =
new TH1F(
"hE1Chi2",
"point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00066 hE2Chi2 =
new TH1F(
"hE2Chi2",
"point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00067
00068 hEChi2_low =
new TH1F(
"hEChi2_low",
"#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00069 hE1Chi2_low =
new TH1F(
"hE1Chi2_low",
"point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00070 hE2Chi2_low =
new TH1F(
"hE2Chi2_low",
"point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00071
00072 hEChi2_hi =
new TH1F(
"hEChi2_hi",
"#sum_{points} #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00073 hE1Chi2_hi =
new TH1F(
"hE1Chi2_hi",
"point1 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00074 hE2Chi2_hi =
new TH1F(
"hE2Chi2_hi",
"point2 #frac{(Eu-Ev)^{2}}{Nmip}",50,0.,10.);
00075
00076 hRatio =
new TH1F(
"hRatio",
"E_{smd} / E_{towers}",20,0.,0.05);
00077 hRatio_low =
new TH1F(
"hRatio_low",
"E_{smd} / E_{towers}",20,0.,0.05);
00078 hRatio_hi =
new TH1F(
"hRatio_hi",
"E_{smd} / E_{towers}",20,0.,0.05);
00079
00080 hEvents =
new TH1F(
"hEvents",
"event counter",1,0.,1.); hEvents->SetBit(TH1::kCanRebin);
00081
00082
00083
00084
if ( mFile )
00085 {
00086 mTree=
new TTree(
"mPi0Tree",
"EEmc pi0 events");
00087 mPi0Event=
new StEEmcMixEvent();
00088 mTree->Branch(
"pi0event",&mPi0Event,32000,99);
00089 mTree->SetDirectory( mFile );
00090 }
00091
00092
return StMaker::Init();
00093 }
00094
00095
00096 Int_t StEEmcPi0Maker::Make()
00097 {
00098
00099
StMuDstMaker *
mudst=(
StMuDstMaker*)GetMaker(
"MuDst");
00100 assert(mudst);
00101
00102 hEvents->Fill(
"event",1.0);
00103
00104
if ( !checkTrigger() )
return kStOK;
00105 hEvents->Fill(
"trigger",1.0);
00106
00107 StThreeVectorF pv=mudst->
muDst()->
event()->
primaryVertexPosition();
00108 TVector3 vertex( pv.x(), pv.y(), pv.z() );
00109
00110
if ( pv.z() < -150. || pv.z() > 150.0 || pv.z()==0.0 )
return kStOK;
00111 hEvents->Fill(
"vertex",1.0);
00112
00113
00114 StEEmcPointVec_t
points =
mEEpoints->
points();
00115
00116
if ( points.size() < 2 )
return kStOK;
00117 hEvents->Fill(
"two+ points",1.0);
00118
00119 mPi0Event->setEvent( mudst->
muDst()->
event() );
00120
00121
00122
for ( UInt_t ipoint=0;ipoint<points.size();ipoint++ )
00123 {
00124
StEEmcPoint p1 = points[ipoint];
00125
00126
00127
for ( UInt_t jpoint=0;jpoint<points.size();jpoint++ )
00128 {
00129
if (ipoint>=jpoint)
continue;
00130
StEEmcPoint p2=points[jpoint];
00131
00132
00133 Int_t phibin1 = p1.
tower(0).phibin();
00134 Int_t phibin2 = p2.tower(0).phibin();
00135
00136
00137
00138
00139
00140 Int_t dphi = phibin1-phibin2;
00141
if ( dphi < 0 ) dphi+=60;
00142
if ( dphi > 30 ) dphi=60-dphi;
00143
00144
if ( dphi > 10 )
continue;
00147
StEEmcPair pair(p1,p2,vertex,vertex);
00148
00149
00150
00151
00152
00153
00154
00155
00156 mPairs.push_back( pair );
00157 mPi0Event->addPair( pair );
00158
00159
00160
00161
00162 StEEmcClusterVec_t tower_clusters;
00163
StEEmcCluster c1 = p1.
clusters(0)[0];
00164
StEEmcCluster c2 = p2.clusters(0)[0];
00165 tower_clusters.push_back( c1 );
00166
if ( c1.key() != c2.
key() ) {
00167 tower_clusters.push_back(c2);
00168 }
00169
00170 mPi0Event->mNumberOfTowerClusters.push_back( (Int_t)tower_clusters.size() );
00171
00172 Int_t ntracks = 0;
00173 Int_t npoints = 0;
00174
for ( UInt_t i=0;i<tower_clusters.size();i++ )
00175 {
00176 ntracks += mEEclusters->numberOfTracks( tower_clusters[i] );
00177 npoints +=
mEEpoints->
numberOfPoints( tower_clusters[i] );
00178 }
00179 mPi0Event->mNumberOfTracks.push_back( ntracks );
00180 mPi0Event->mNumberOfPoints.push_back( npoints );
00181
00182
00183
00184
00185
00186
00187
00188
00190
00191
00192 hEvents->Fill(
"pair",1.0);
00193
00194
if ( pair.
momentum().Perp() < 3.0 )
continue;
00195
00196 hEvents->Fill(
"pair pt>3",1.0);
00197
00198
00199
00200
StEEmcCluster cluster1=p1.
clusters(0).at(0);
00201
StEEmcCluster cluster2=p2.clusters(0).at(0);
00202
00203
00204
00205
00206 hEvents->Fill(
"pair in cl",1.0);
00207
00208
00209
00210
if ( cluster1.
numberOfMatchingClusters(4) > 2 ||
00211 cluster1.
numberOfMatchingClusters(5) > 2 )
continue;
00212
00213
StEEmcSmdCluster u1=p1.
cluster(0);
00214
StEEmcSmdCluster u2=p2.cluster(0);
00215
StEEmcSmdCluster v1=p1.
cluster(1);
00216
StEEmcSmdCluster v2=p2.cluster(1);
00217
00218 Bool_t didSplit =
false;
00219
if ( u1.
key() == u2.
key() || v1.
key() == v2.
key() ) didSplit=
true;
00220
00221
00222
00223
00224
00225 Float_t esum_smd = u1.
energy()+u2.
energy()+v1.
energy()+v2.
energy();
00226 Float_t esum_tow = cluster1.
energy();
00227
if ( cluster1.
key() != cluster2.
key() ) esum_tow += cluster2.
energy();
00228
00229
00230
00231 hEvents->Fill(
"pair ncl",1.0);
00232
00233 hMass->Fill( pair.
mass(), pair.
pt() );
00234
if( didSplit ) hMass_split->Fill( pair.
mass(), pair.
pt() );
00235
00236 Float_t ediff1 = 1000.0*(u1.
energy()-v1.
energy());
00237 Float_t ediff2 = 1000.0*(u2.
energy()-v2.
energy());
00238 Float_t nmip1 = 1000.0*(u1.
energy()+v1.
energy())/1.3;
00239 Float_t nmip2 = 1000.0*(u2.
energy()+v2.
energy())/1.3;
00240
if ( nmip1>0. && nmip2>0. )
00241 {
00242 Float_t e1chi2 = ediff1*ediff1/nmip1;
00243 Float_t e2chi2 = ediff2*ediff2/nmip2;
00244 Float_t echi2 = e1chi2+e2chi2;
00245
00246
if ( e1chi2 < 4.0 && e2chi2 < 4.0 ) {
00247 hMass_cut -> Fill( pair.
mass(), pair.
pt() );
00248
if ( pair.
mass()>0.1 && pair.
mass()<0.18 ) hPT_cut->Fill( pair.
pt() );
00249
if ( pair.
mass()>0.1 && pair.
mass()<0.18 ) hZgg_cut->Fill( pair.
zgg() );
00250 }
00251
00252
if ( pair.
mass() < 0.1 ) {
00253 hEChi2_low->Fill(echi2);
00254 hE1Chi2_low->Fill(e1chi2);
00255 hE2Chi2_low->Fill(e2chi2);
00256 hRatio_low -> Fill( esum_smd/esum_tow );
00257 }
00258
if ( pair.
mass() > 0.2 ) {
00259 hEChi2_hi->Fill(echi2);
00260 hE1Chi2_hi->Fill(e1chi2);
00261 hE2Chi2_hi->Fill(e2chi2);
00262 hRatio_hi -> Fill( esum_smd/esum_tow );
00263 }
00264 }
00265
00266
00267
if ( pair.
mass() > 0.1 && pair.
mass() < 0.18 ) {
00268
00269 hPT->Fill( pair.
momentum().Perp() );
00270 hXF->Fill( pair.
momentum().Z()/100.0 );
00271 hEnergy->Fill( pair.
energy() );
00272 hEta->Fill( pair.
momentum().Eta() );
00273 hPhi->Fill( pair.
momentum().Phi() );
00274 hZvertex->Fill( pair.
vertex().Z() );
00275 hZgg->Fill( pair.
zgg() );
00276
00277
if ( nmip1==0. || nmip2==0. )
continue;
00278 Float_t e1chi2 = ediff1*ediff1/nmip1;
00279 Float_t e2chi2 = ediff2*ediff2/nmip2;
00280 Float_t echi2 = e1chi2+e2chi2;
00281
00282 hEChi2->Fill(echi2);
00283 hE1Chi2->Fill(e1chi2);
00284 hE2Chi2->Fill(e2chi2);
00285 hRatio -> Fill( esum_smd/esum_tow );
00286 printf(
"esum_smd/esum_tow = %5.3f\n", esum_smd/esum_tow );
00287
00288 }
00289
00290 }
00291 }
00292
00293
if ( mTree )
00294 mTree->Fill();
00295
00296
return kStOK;
00297 }
00298
00299
00300
00301
00302
void StEEmcPi0Maker::Clear(Option_t *opts)
00303 {
00304 mPairs.clear();
00305 mPi0Event->Clear();
00306 }
00307
00308
00309
00310
00311
void StEEmcPi0Maker::setCheckTrigger(Bool_t t){ mCheckTrigger=t; }
00312
void StEEmcPi0Maker::addTrigger(Int_t t){
mTriggerList.push_back(t); }
00313
00314 Bool_t StEEmcPi0Maker::checkTrigger()
00315 {
00316
00318
if (
mTriggerList.size() == 0 )
return 1;
00319
00320 StTriggerId nominal;
00321
00323
StMuDstMaker *mumk = (
StMuDstMaker*)GetMaker(
"MuDst");
00324
StEvent *event = (
StEvent*)GetInputDS(
"StEvent");
00325
00326
if ( mumk )
00327 {
00328 nominal = mumk->
muDst()->
event()->
triggerIdCollection().
nominal();
00329
goto CHECK_TRIGGER;
00330 }
00331
00333
00334
if ( event )
00335 {
00336 nominal=*event->
triggerIdCollection()->
nominal();
00337
goto CHECK_TRIGGER;
00338 }
00339
00341
goto NO_DATA;
00342
00343
00344 CHECK_TRIGGER:
00345
00346
for ( UInt_t ii=0;ii<
mTriggerList.size();ii++ )
00347 {
00348
if ( nominal.isTrigger( mTriggerList[ii] ) )
return mTriggerList[ii];
00349 }
00350
return 0;
00351
00352 NO_DATA:
00353 assert(2+2==5);
00354
return 0;
00355
00356 }