00001 #include "StEEmcMixMaker.h"
00002 #include "StEEmcPointMaker.h"
00003 #include "StEEmcA2EMaker.h"
00004
00005 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00006 #include "StMuDSTMaker/COMMON/StMuDst.h"
00007 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00008 #include "StEvent/StTriggerIdCollection.h"
00009 #include "StEvent/StTriggerId.h"
00010
00011 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00012
00013 #include "TRandom.h"
00014 #include "TH1F.h"
00015 #include "TH2F.h"
00016
00017 #define SIMPLE_MC
00018
00019 ClassImp(StEEmcMixMaker);
00020
00021
00022 StEEmcMixMaker::StEEmcMixMaker(const Char_t *name,Int_t s):StMaker(name)
00023 {
00024 mRandom=new TRandom();
00025 mPoolSize=s;
00026
00027 mEEmcTow = new EEmcGeomSimple();
00028
00032 minET(1.0);
00033 maxZ(0.5);
00034 minEpoint(0.8);
00035 range(0.11,0.17);
00036 mTrigMode=0;
00037
00039 StEEmcPointVec_t points;
00040 for ( Int_t i=0; i<40; i++ ) mPool.push_back(points);
00041
00042 mFixedVertex=TVector3(-999.,-999.,-999.);
00043 mSigmaVertex=-999.;
00044
00045 }
00046
00047
00048 Int_t StEEmcMixMaker::Init()
00049 {
00050 mEEpoints=(StEEmcPointMaker *)GetMaker(mPointMakerName);
00051 assert(mEEpoints);
00052 mEEanalysis=(StEEmcA2EMaker *)GetMaker(mAnalysisName);
00053 assert(mEEanalysis);
00054 mMuDstMaker=(StMuDstMaker *)GetMaker(mMuDstMakerName);
00055 assert(mMuDstMaker);
00056
00058 book();
00059
00060 return StMaker::Init();
00061 }
00062
00063
00064 Int_t StEEmcMixMaker::Make()
00065 {
00066
00068 if ( !accept( mMuDstMaker -> muDst() -> event() ) )
00069 return kStOK;
00070
00071 mH1[1]->Fill("accepted",1.0);
00072
00074 mPoints=mEEpoints->points();
00075
00077 if ( mPoints.size() <= 1 ) return kStOK;
00078 mH1[1]->Fill( "2+ points", 1.0);
00079
00081 mixReal();
00082
00084 mixBackground();
00085
00087 fill();
00088
00091 fillPool();
00092
00093 return kStOK;
00094 }
00095
00096
00097 void StEEmcMixMaker::Clear( Option_t *opts )
00098 {
00099 mCandidates.clear();
00100 mBackground.clear();
00101 mPoints.clear();
00102
00103 return;
00104 }
00105
00106
00107 void StEEmcMixMaker::mixReal()
00108 {
00109
00110 if ( !mPoints.size() ) return;
00111
00112 Float_t emax=0.;
00113 Int_t imax=0;
00114 Int_t count=0;
00115
00116
00117 StMuEvent *event = mMuDstMaker -> muDst() -> event();
00118 if ( !event ) return ;
00119 StThreeVectorF v=event->primaryVertexPosition();
00120 TVector3 vertex = TVector3(v.x(),v.y(),v.z());
00121
00122
00124 if ( mFixedVertex.Z()<-500. ) {
00125 if ( v.z()==0. && v.x()==0. && v.y()==0. ) return;
00126 }
00127 else {
00128 vertex=mFixedVertex;
00129 if ( mSigmaVertex > 0. ) vertex[2]=( mFixedVertex.Z() + mRandom->Gaus(0.,mSigmaVertex));
00130 }
00131
00132
00134 for ( UInt_t ipoint=0; ipoint<mPoints.size()-1; ipoint++ )
00135 for ( UInt_t jpoint=ipoint+1; jpoint<mPoints.size(); jpoint++ )
00136 {
00137
00138 StEEmcPoint point1=mPoints[ipoint];
00139 StEEmcPoint point2=mPoints[jpoint];
00140
00142 Bool_t go1 = false;
00143 Bool_t go2 = false;
00144 for ( UInt_t isec=0;isec<mSectorList.size();isec++ )
00145 {
00146 go1 |= point1.sector() == mSectorList[isec];
00147 go2 |= point2.sector() == mSectorList[isec];
00148 }
00149 if ( !(go1&&go2) ) continue;
00150
00151
00153 if ( point1.sector() != point2.sector() ) continue;
00154
00155 count++;
00156 mCandidates.push_back ( StEEmcPair( point1, point2, vertex, vertex ) );
00157 if ( mCandidates.back().energy() > emax ) {
00158 emax = mCandidates.back().energy();
00159 imax=count;
00160 }
00161
00162 }
00163
00164
00165 }
00166
00167 void StEEmcMixMaker::mixBackground()
00168 {
00169
00170 if ( !mPoints.size() ) return ;
00171
00173 StEEmcTower high_tower = mEEanalysis->hightower();
00174
00175 Int_t bin=(Int_t)(high_tower.adc()/100);
00176 if ( bin<0 ) bin=0;
00177 if ( bin>40 ) bin=40;
00178
00181 StEEmcPointVec_t &points = mPool[bin];
00182
00183 StMuEvent *event = mMuDstMaker -> muDst() -> event();
00184 if ( !event ) return;
00185 StThreeVectorF v=event->primaryVertexPosition();
00186 TVector3 vertex = TVector3(v.x(),v.y(),v.z());
00187
00188 if ( mFixedVertex.Z()<-500. ) {
00189 if ( v.z()==0. && v.x()==0. && v.y()==0. ) return;
00190 }
00191 else {
00192 vertex=mFixedVertex;
00193 if ( mSigmaVertex > 0. ) vertex[2]=( mFixedVertex.Z() + mRandom->Gaus(0.,mSigmaVertex));
00194 }
00195
00196
00198 for ( UInt_t i=0; i<mPoints.size(); i++ ) {
00199
00200 StEEmcPoint point1=mPoints[i];
00201
00203 for ( UInt_t j=0; j<points.size(); j++ ) {
00204
00205 StEEmcPoint point2=points[j];
00206
00208 Bool_t go1 = false;
00209 Bool_t go2 = false;
00210 for ( UInt_t isec=0;isec<mSectorList.size();isec++ )
00211 {
00212 go1 |= point1.sector() == mSectorList[isec];
00213 go2 |= point2.sector() == mSectorList[isec];
00214 }
00215 if ( !(go1&&go2) ) continue;
00216
00218 if ( point1.sector() != point2.sector() ) continue;
00219
00221 mBackground.push_back( StEEmcPair( point1, point2, vertex, vertex ) );
00222
00223
00224
00225 }
00226
00227 }
00228
00229 }
00230
00231
00232
00233
00234 Bool_t StEEmcMixMaker::accept( StMuEvent *event )
00235 {
00236 if ( !event ) return false;
00237 StMuTriggerIdCollection tic = event -> triggerIdCollection();
00238 StTriggerId l1trig = tic.l1();
00239
00241 if ( mTriggerList.size() <= 0 ) {
00242 mH1[0]->Fill("no selection",1.0);
00243 return true;
00244 }
00245
00246 Int_t go=0;
00247 std::vector<Int_t>::iterator iter=mTriggerList.begin();
00248 while ( iter != mTriggerList.end() ) {
00249 go = l1trig.isTrigger( (*iter) );
00250 if ( go ) {
00251 go = (*iter);
00252 break;
00253 }
00254 iter++;
00255 }
00256 TString name=go;
00257 mH1[0]->Fill(name,1.0);
00258 return (go!=0);
00259 }
00260
00261
00262
00263 void StEEmcMixMaker::book()
00264 {
00265
00266
00267 mH1.push_back(new TH1F("triggers","Number of triggers fired",1,0.,1.));
00268 mH1[0]->SetBit(TH1::kCanRebin);
00269 mH1.push_back(new TH1F("status","Events processed up to...",1,0.,1.));
00270 mH1[1]->SetBit(TH1::kCanRebin);
00271
00272
00273 mH2.push_back(new TH2F("uvha","<u> vs <v> for higher-energy gamma",288,0.,288.,288,0.,288.));
00274 mH2.push_back(new TH2F("uvla","<u> vs <v> for lower-energy gamma",288,0.,288.,288,0.,288.));
00275 mH2.push_back(new TH2F("uvhc","<u> vs <v> for higher-energy gamma, mass cut",288,0.,288.,288,0.,288.));
00276 mH2.push_back(new TH2F("uvlc","<u> vs <v> for lower-energy gamma, mass cut",288,0.,288.,288,0.,288.));
00277
00278
00279 mH1real.push_back(new TH1F("massR","Invariant mass of photon pairs",360,0.,3.6) );
00280 mH1real.push_back(new TH1F("energyR","Energy of photon pairs",200,0.,40.));
00281 mH1real.push_back(new TH1F("zggR","Energy sharing of photon pairs",50,0.,1.));
00282 mH1real.push_back(new TH1F("phiggR","Opening angle of photon pairs",100,0.,0.1));
00283 mH1real.push_back(new TH1F("ptR","p_{T} of photon pairs",100,0.,10.));
00284 mH1real.push_back(new TH1F("zvertexR","Z_{vertex} [cm]",100,-100.,100.));
00285
00286 mH1mix.push_back(new TH1F("massM","Invariant mass of photon pairs",360,0.,3.6) );
00287 mH1mix.push_back(new TH1F("energyM","Energy of photon pairs",200,0.,40.));
00288 mH1mix.push_back(new TH1F("zggM","Energy sharing of photon pairs",50,0.,1.));
00289 mH1mix.push_back(new TH1F("phiggM","Opening angle of photon pairs",100,0.,0.1));
00290 mH1mix.push_back(new TH1F("ptM","p_{T} of photon pairs",100,0.,10.));
00291 mH1mix.push_back(new TH1F("zvertexM","Z_{vertex} [cm]",100,-100.,100.));
00292
00293 }
00294
00295
00296 void StEEmcMixMaker::fill()
00297 {
00298
00299 StEEmcPairVec_t::iterator ipair=mCandidates.begin();
00300 while ( ipair != mCandidates.end() ) {
00301 fill( mH1real, (*ipair ) );
00302 fillQA( mH2, (*ipair ) );
00303 ipair++;
00304 }
00305
00306 ipair=mBackground.begin();
00307 while ( ipair != mBackground.end() ) {
00308 fill( mH1mix, (*ipair) );
00309 ipair++;
00310 }
00311
00312 }
00313
00314 void StEEmcMixMaker::fillQA( std::vector<TH2F *> &h, StEEmcPair pair )
00315 {
00316
00317 StEEmcPoint p1=pair.point(0);
00318 StEEmcPoint p2=pair.point(1);
00319 StEEmcSmdCluster uh,ul,vh,vl;
00320 if ( p1.energy() > p2.energy() ) {
00321 uh=p1.cluster(0); vh=p1.cluster(1);
00322 ul=p2.cluster(0); vl=p2.cluster(1);
00323 }
00324 else {
00325 ul=p1.cluster(0); vl=p1.cluster(1);
00326 uh=p2.cluster(0); vh=p2.cluster(1);
00327 }
00328 mH2[0]->Fill( uh.mean(), vh.mean() );
00329 mH2[1]->Fill( ul.mean(), vl.mean() );
00330 if ( pair.mass() > mMinMass && pair.mass() < mMaxMass ) {
00331 mH2[2]->Fill( uh.mean(), vh.mean() );
00332 mH2[3]->Fill( ul.mean(), vl.mean() );
00333 }
00334
00335 }
00336
00337 void StEEmcMixMaker::fill( std::vector<TH1F*> &h, StEEmcPair pair )
00338 {
00339
00340 h[0]->Fill( pair.mass() );
00341 if ( pair.mass() > mMinMass && pair.mass() < mMaxMass ) {
00342 h[1]->Fill( pair.energy() );
00343 h[2]->Fill( pair.zgg() );
00344 h[3]->Fill( pair.phigg() );
00345 h[4]->Fill( pair.vertex().Z() );
00346 }
00347
00348 }
00349
00350
00351
00352 void StEEmcMixMaker::fillPool()
00353 {
00354
00355 if ( !mPoints.size() ) return;
00356
00358 StEEmcTower high_tower = mEEanalysis->hightower();
00359
00360 Int_t bin=(Int_t)(high_tower.adc()/100);
00361 if ( bin<0 ) bin=0;
00362 if ( bin>40 ) bin=40;
00363
00366 StEEmcPointVec_t &points = mPool[bin];
00367
00369 std::reverse(points.begin(),points.end());
00370
00372 for ( UInt_t i=0; i<mPoints.size(); i++ ) {
00373
00374 TVector3 d=mPoints[i].position();
00375
00377 Int_t sec,sub,eta;
00378 if ( !mEEmcTow->getTower(d,sec,sub,eta) ) continue;
00379
00381 if ( mTrigMode == 1 ) {
00382 if ( sec==high_tower.sector() &&
00383 sub==high_tower.subsector() &&
00384 eta==high_tower.etabin() ) continue;
00385 }
00386
00388 points.push_back( mPoints[i] );
00389
00390 }
00391
00393 std::reverse(points.begin(),points.end());
00394
00396 if ( points.size() > (UInt_t)mPoolSize ) points.resize(mPoolSize);
00397
00398 }
00399
00400
