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