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