00001
00014
00015
00016 #include "StEEmcIUPi0Analysis.h"
00017 #include "StEEmcIUMixMaker.h"
00018 #include "StEEmcIUPointMaker.h"
00019 #include "TH1F.h"
00020 #include "TH2F.h"
00021 #include "TH2D.h"
00022 #include "TTree.h"
00023 #include "TFile.h"
00024 #include <stdio.h>
00025 #include <iostream>
00026 #include "SpinCutsIU.h"
00027 #include "StEEmcPool/StEEmcPointMaker/EEmcSectorFit.h"
00028 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00029 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00030 #include "StEEmcPool/StRFEmcTrigger/StRFEmcTrigMaker.h"
00031 #include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
00032 #include "StMcOutputMaker.h"
00033 ClassImp(StEEmcIUPi0Analysis);
00034
00035
00036 StEEmcIUPi0Analysis::StEEmcIUPi0Analysis(const Char_t *name):StMaker(name)
00037 {
00038 mCuts=new SpinCutsIU();
00039 mTrigSim = 0;
00040 mTrigSimThreshold = 0;
00041 mEEgeom=new EEmcGeomSimple();
00042 mSpinSort = true;
00043
00044 }
00045
00046
00047 Int_t StEEmcIUPi0Analysis::Init()
00048 {
00049
00050
00051 mHistograms[ kAny ] = new SpinIUHistos("Any", "Pi0 spectra, unsorted");
00052
00053 mBackgrounds[ kAny ] = new SpinIUHistos("AnyB", "combinatoric spectra, unsorted");
00054 hPi0Mass=new TH1F("hMass","invariant mass in windows [0.08,0.18]Gev and DeltaR<0.04",120,0.,1.2);
00055 hPTvsPiNum = new TH1F("ReconPt", "Recon PT",50,0.,25. );
00056 hMcEta = new TH1F("McEta", "Mc Eta",50,0.,2.5 );
00057 hEEmcEta = new TH1F("EEmcEta", "EEmc Eta",50,0.,2.5 );
00058 hMcPhi = new TH1F("McPhi", "Mc Phi",360,-180.,180. );
00059 hMcPt = new TH1F("McPt", "Mc Pt",50,0.,25. );
00060 hMcEnergy=new TH1F("McEnergy","Mc energy for recon pi0",80,0.,80.);
00061 hREtaPhi=new TH2F("REEmcEtaPhi", "Reconstructed pi0 track; Recon phi / deg;ReconEta",360,-180.,180,20,1.,2.);
00062 hReconEta = new TH1F("ReconEta", "Reconstructed Eta",50,0.,2.5 );
00063 hReconPhi = new TH1F("ReconPhi", "Recon Phi",360,-180.,180. );
00064 hRecoEnergy = new TH1F("ReconEnergy","Reconstructed Pi0 Energy",60,0.,60.);
00065 hResoEta = new TH1F("ResoEta", "Resolution of Eta: ReconEta - EEmcEta",10,-1.,1. );
00066 hResoPhi = new TH1F("ResoPhi", "Resolution of Phi: ReconPhi - McPhi",20,-20.,20. );
00067 hResoPt = new TH1F("ResoPt", "Resolution of Pt",160,-4.,4. );
00068 hResoEnergy=new TH1F("ResoEnergy", "Resolution of Energy",200,-10.,10.);
00069 hResoPtvsG = new TH2F("ResoPtvsG", "Resolution of Pt; mcPt-ReconPt;mcPt",80,-4.,4.,250,0.,25.);
00070 hResoEnergyvsG = new TH2F("ResoEnergyvsG", "Resolution of Energy; mcEnergy-ReconEnergy;mcEnergy",200,-10.,10.,80,0.,80.);
00071 hRealMcPD= new TH1F("RealMcdistance","Distance between real and MC pi0",100,0.,0.5);
00072 hRealMcPR= new TH1F("RRealMcdistance","Distance between real and MC pi0 in the EEmc detector",100,0.,0.5);
00073 mHistograms[1]=new SpinIUHistos("PP","Pi0 spectra, PP, spin4=dec5");
00074 mHistograms[2]=new SpinIUHistos("PN","Pi0 spectra, PN, spin4=dec6");
00075 mHistograms[3]=new SpinIUHistos("NP","Pi0 spectra, NP, spin4=dec9");
00076 mHistograms[4]=new SpinIUHistos("NN","Pi0 spectra, NN, spin4=dec10");
00077 EventCounter=new TH1F("EventCounter","Event counter",10,0.,10.);
00078 hEventCounter=new TH1F("hEventCounter","Event counts",10,0.,10.);
00079 hEventCounter -> GetXaxis() -> SetBinLabel(1,"raw");
00080 hEventCounter -> GetXaxis() -> SetBinLabel(2,"minb");
00081 hEventCounter -> GetXaxis() -> SetBinLabel(3,"hw trg");
00082 hEventCounter -> GetXaxis() -> SetBinLabel(5,"sw trg");
00083 hEventCounter -> GetXaxis() -> SetBinLabel(5,"sw trg[3,6]");
00084 hEventCounter -> GetXaxis() -> SetBinLabel(6,"vertex");
00085 hPairCounter = new TH1F("hPairCounter","Pair counts",10,0.,10.);
00086 hPi0Counter = new TH1F("hPi0Counter","Pi0 counts",10,0.,10.);
00087 hFillPatternI = new TH1F("hFillPatternI","Intended fill pattern:bunch X-ing @ STAR IP:n runs",120,0.,120.);
00088 hSpin4 = new TH1F("hSpin4","Spin 4:spin4 state",11,0.,11.);
00089 hSpin4mb = new TH1F("hSpin4mb","Spin 4:spin4 state",11,0.,11.);
00090 hBx7 = new TH1F("hBx7","7-bit bunch crossing:bx7",120,0.,120.);
00091 hBx48 = new TH1F("hBx48","48-bit bunch crossing:bx48",120,0.,120.);
00092 hBx7diffBx48 = new TH2F("hBx7diffBx48","bx1=(bx7-off7)%120, bx2=(bx48-off48)%120;bx7;bx1-bx2",120,0.,120.,21,-10.5,10.5);
00093 hBxStar = new TH1F("hBxStar","Beam x-ing at STAR IP;star bunch x-ing",120,0.,120.);
00094 hBxStarPi0 = new TH1F("hBxStarPi0","Beam x-ing at STAR IP with pi0 detected;star bunch x-ing",120,0.,120.);
00095 hMassBx=new TH2F("hMassBx","Beam x-ing vs Mass;M [GeV];STAR bunch xing",120,0.,1.2,120,0.,120.);
00096 hZvertBx=new TH2F("hZvertBx","Beam x-ing vs Z vertex [0.1,0.18]",150,-150.,150.,120,0.,120.);
00097 hZggBx=new TH2F("hZggBx","Beam x-ing vs Zgg",100,0.,1.,120,0.,120.);
00098 hEtaBx=new TH2F("hEtaBx","Beam x-ing vs eta",120,0.8,2.2,120,0.,120.);
00099 DEtaMass=new TH2F("DEtaMass","Delta Eta of the two points vs Mass of the reconstructed pi0",60,0.,0.6,20,0.,0.1);
00100 DPhiMass=new TH2F("DPhiMass","Delta Phi of the two points vs Mass of the reconstructed pi0",60,0.,0.6,20,0.,0.1);
00101 dUvsdV=new TH2D("dUvsdV","Seed distance between smd clusters; dU; dV",10,0,10,10,0,10);
00102 dUvsdVGood=new TH2D("dUvsdVGood","Good Seed distance between smd clusters; dU; dV",10,0,10,10,0,10);
00103 dUvsRatio=new TH2F("dUvsRatio","Ucluster ratio along with seed distance;dU;Ratio",20,0.,20.,10,0.,1.);
00104 dVvsRatio=new TH2F("dVvsRatio","Vcluster ratio along with seed distance;dV;Ratio",20,0.,20.,10,0.,1.);
00105 dUvsRatioGood=new TH2F("dUvsRatioGood","Good Ucluster ratio along with seed distance;dU;Ratio",20,0.,20.,10,0.,1.);
00106 dVvsRatioGood=new TH2F("dVvsRatioGood","Good Vcluster ratio along with seed distance;dV;Ratio",20,0.,20.,10,0.,1.);
00107 dUVzeroE=new TH1F("dUVzeroE","Pi0 energy with dU or dV=0",60,0.,60);
00108 dUVzeroEta=new TH1F("dUVzeroEta","Pi0 Eta with dU or dV=0",20,1.,2.);
00109 GoodPiGeo=new TH2F("GoodPiGeo","Pi0 distribution in EEMC with mass range [0.08,0.18]Gev;x;y",120,-240.,240.,120,-240.,240.);
00110 GoodPiPt=new TH1F("GoodPiPt", "Recon PT",50,0.,25.);
00111 mRealEvent=new StEEmcIUMixEvent();
00112 mMixedEvent=new StEEmcIUMixEvent();
00113 mRealTree=new TTree("mRealTree","Real Events");
00114 mRealTree->Branch("MixEvent","StEEmcIUMixEvent",&mRealEvent);
00115 mMixedTree=new TTree("mMixedTree","Mixed Events");
00116 mMixedTree->Branch("MixEvent","StEEmcIUMixEvent",&mMixedEvent);
00117
00118 AddObj( mRealTree, ".hist" );
00119 AddObj( mMixedTree, ".hist" );
00120
00121 return StMaker::Init();
00122 }
00123
00124 Int_t StEEmcIUPi0Analysis::InitRun(Int_t run)
00125 {
00126
00127 mSpinDb->print(0);
00128 Info("InitRun","run number = %d",run);
00129 const Int_t *spin8bits=mSpinDb->getSpin8bits();
00130 for(int bx=0;bx<120;bx++){
00131 Bool_t isFilled=(spin8bits[bx] & 0x11)==0x11;
00132 if(isFilled) hFillPatternI->Fill(bx);
00133 }
00134
00135 mRunNumber = run;
00136
00137 return StMaker::InitRun(run);
00138
00139 }
00140
00141 Int_t StEEmcIUPi0Analysis::Make()
00142 {
00143
00145 EventCounter->Fill(0.,1.);
00146 Int_t spin4 = -1;
00147 Int_t bxStar = -1;
00148 if ( mSpinDb -> isValid() && mSpinDb -> isPolDirLong() )
00149 spin4=getSpinState( mMuDst->muDst(), bxStar );
00150
00151
00152 hBxStar -> Fill( bxStar );
00153
00154 mRealEvent -> setEvent( mMuDst->muDst()->event() );
00155 mRealEvent -> setSpin4( spin4 );
00156 mMixedEvent -> setEvent( mMuDst->muDst()->event() );
00157 mMixedEvent -> setSpin4( spin4 );
00158
00159 for ( Int_t ii=0;ii<720;ii++ ) {
00160 StEEmcTower tower=mEEanalysis->tower(ii);
00161 if ( tower.fail() ) continue;
00162 mRealEvent->mADC[ii]=(60./4096.)*tower.adc();
00163 mRealEvent->mStat[ii]=tower.stat();
00164 }
00165
00166
00168
00169 if ( !accept( mMuDst->muDst() ) ) goto END_OF_MAKE;
00170 EventCounter->Fill(2.,1.);
00171 hPairCounter -> Fill( kEvent );
00172
00173
00175 static Int_t mymap[]={0,0,0,0,0,1,2,0,0,3,4};
00176
00178 mRealEvent -> mTotalEnergyT = mEEanalysis -> energy(0);
00179 mRealEvent -> mTotalEnergyP = mEEanalysis -> energy(1);
00180 mRealEvent -> mTotalEnergyQ = mEEanalysis -> energy(2);
00181 mRealEvent -> mTotalEnergyR = mEEanalysis -> energy(3);
00182 mRealEvent -> mTotalEnergyU = mEEanalysis -> energy(4);
00183 mRealEvent -> mTotalEnergyV = mEEanalysis -> energy(5);
00184
00186 numcutCan=0;
00187 numoacceptedpair=0;
00188
00189
00190
00191 for ( Int_t i=0;i<mEEmixer->numberOfCandidates();i++ )
00192 {
00193
00194 StEEmcIUPair pair = mEEmixer->candidate(i);
00195 hPairCounter -> Fill( kPair );
00196
00197 if ( !accept( pair ) ) continue;
00198
00199 mRealEvent -> addPair( pair );
00200 numoacceptedpair++;
00201
00202
00203
00204
00205 mHistograms[ kAny ] -> Fill( pair );
00206
00207 StEEmcIUPoint point1=pair.point(0);
00208 StEEmcIUPoint point2=pair.point(1);
00209
00210 TVector3 point1pos=point1.position();
00211 TVector3 point2pos=point2.position();
00212 float deta=fabs(point1pos.Eta()-point2pos.Eta());
00213 float dphi=fabs(point1pos.Phi()-point2pos.Phi());
00214 DEtaMass->Fill(pair.mass(),deta);
00215 DPhiMass->Fill(pair.mass(),dphi);
00216
00217 float Rxy=TMath::Sqrt(pair.vertex().x()*pair.vertex().x()+pair.vertex().y()*pair.vertex().y());
00218 float hHeight=pair.pt()*(270.0-pair.vertex().Z())/pair.pz()+Rxy;
00219 float etatheta=TMath::ATan(hHeight/270.0);
00220
00221 float mideta=tan(etatheta/2.0);
00222 float eemceta=-log(mideta);
00223 hEEmcEta->Fill(eemceta);
00224
00225 float EsmallU=(point1.cluster(0).energy()<point2.cluster(0).energy())?point1.cluster(0).energy():point2.cluster(0).energy();
00226 float ElargeU=(point1.cluster(0).energy()>point2.cluster(0).energy())?point1.cluster(0).energy():point2.cluster(0).energy();
00227 float EsmallV=(point1.cluster(1).energy()<point2.cluster(1).energy())?point1.cluster(1).energy():point2.cluster(1).energy();
00228 float ElargeV=(point1.cluster(1).energy()>point2.cluster(1).energy())?point1.cluster(1).energy():point2.cluster(1).energy();
00229 float dU=abs(point1.cluster(0).strip(0).index()-point2.cluster(0).strip(0).index());
00230 float dV=abs(point1.cluster(1).strip(0).index()-point2.cluster(1).strip(0).index());
00231 if(pair.mass()<0.06) {
00232 dUvsdV->Fill(dU,dV);
00233
00234 if(dU==0||dV==0)
00235 {
00236 dUVzeroE->Fill(pair.energy());
00237 dUVzeroEta->Fill(pair.momentum().Eta());
00238 }
00239 if(dU==0)
00240 {
00241 dVvsRatio->Fill(dV,EsmallV/ElargeV);
00242 }
00243 if(dV==0)
00244 {
00245 dUvsRatio->Fill(dU,EsmallU/ElargeU);
00246 }
00247 }
00248 if(pair.mass()>=0.08 && pair.mass()<=0.18){
00249 Float_t e1 = pair.point(0).energy();
00250 Float_t e2 = pair.point(1).energy();
00251
00252 TVector3 pp = (e1*point1pos + e2*point2pos) * ( 1/(e1+e2) );
00253 GoodPiGeo->Fill(pp.X(),pp.Y());
00254 GoodPiPt->Fill(pair.pt());
00255 dUvsdVGood->Fill(dU,dV);
00256
00257 if(dU==0)
00258 {
00259 dVvsRatioGood->Fill(dV,EsmallV/ElargeV);
00260 }
00261
00262 if(dV==0)
00263 {
00264 dUvsRatioGood->Fill(dU,EsmallU/ElargeU);
00265 }
00266 }
00267 #ifdef EFFICIENCY
00268
00269
00270 StMcOutputMaker *mkMc=(StMcOutputMaker *)GetMaker("mcRead");
00271 assert(mkMc);
00272
00273 float rPt=0;
00274 float rEta=0;
00275 float rPhi=0;
00276 float rE=0;
00277
00278 if(pair.mass()>=0.08 && pair.mass()<=0.18)
00279 {
00280
00281 numcutCan++;
00282 rPt=pair.pt();
00283 rEta=pair.momentum().Eta();
00284 rPhi=pair.momentum().Phi();
00285 rE=pair.energy();
00286 int size=0;
00287 size=mkMc->gTr.size();
00288 float minD=10000000;
00289 float minR=10000000;
00290
00291 for ( Int_t j=0;j<size;j++ )
00292 {
00293 StMcTrack *gt=mkMc->gTr[j];
00294 float eemceta=mkMc->geemcEta[j];
00295 StLorentzVectorF p4=gt->fourMomentum();
00296 float etaMc=p4.pseudoRapidity();
00297 float phiMc=p4.phi();
00298
00299
00300 float deltaphi=fabs(rPhi-phiMc);
00301 if(deltaphi>=3.14159265)
00302 {deltaphi=6.283-deltaphi;}
00303 float deltaeta=rEta-etaMc;
00304 float Rdeltaeta=rEta-eemceta;
00305 float deltad=TMath::Sqrt(deltaphi*deltaphi+deltaeta*deltaeta);
00306 float deltaR=TMath::Sqrt(deltaphi*deltaphi+Rdeltaeta*Rdeltaeta);
00307 if(deltad<=minD)
00308 {
00309
00310 minD=deltad;
00311 }
00312 if(deltaR<=minR)
00313 {
00314 minR=deltaR;
00315 }
00316 }
00317 hRealMcPD->Fill(minD);
00318 hRealMcPR->Fill(minR);
00319 for ( Int_t j=0;j<size;j++ )
00320 {
00321 StMcTrack *gt=mkMc->gTr[j];
00322 float eemceta=mkMc->geemcEta[j];
00323
00324 StLorentzVectorF p4=gt->fourMomentum();
00325 float etaMc=p4.pseudoRapidity();
00326
00327 float phiMc=p4.phi();
00328
00329 float ptMc=p4.perp();
00330 float energyMc=p4.e();
00331
00332 float deltaphi=fabs(rPhi-phiMc);
00333 if(deltaphi>=3.14159265)
00334 {deltaphi=6.283-deltaphi;}
00335 float deltaeta=rEta-etaMc;
00336 float Rdeltaeta=rEta-eemceta;
00337 float deltad=TMath::Sqrt(deltaphi*deltaphi+deltaeta*deltaeta);
00338 float deltaR=TMath::Sqrt(deltaphi*deltaphi+Rdeltaeta*Rdeltaeta);
00339
00340
00341
00342
00343 if(deltaR<=0.04)
00344 {
00345 hPi0Mass->Fill(pair.mass());
00346 hPTvsPiNum->Fill(rPt);
00347 hREtaPhi->Fill(rPhi/3.1416*180,rEta);
00348 hReconEta->Fill(rEta);
00349 hReconPhi->Fill(rPhi/3.1416*180);
00350 hRecoEnergy->Fill(rE);
00351
00352 hMcEta->Fill(etaMc);
00353 hEEmcEta->Fill(eemceta);
00354 hMcPhi->Fill(phiMc/3.1416*180);
00355 hMcPt->Fill(ptMc);
00356 hResoPt->Fill(rPt-ptMc);
00357 hResoPtvsG->Fill(rPt-ptMc,ptMc);
00358 hResoEta->Fill(rEta-etaMc);
00359 hResoPhi->Fill((rPhi-phiMc)/3.1416*180);
00360 hMcEnergy->Fill(energyMc);
00361 hResoEnergyvsG->Fill(rE-energyMc,energyMc);
00362 hResoEnergy->Fill(rE-energyMc);
00363
00364
00365
00366
00367
00368
00369 }
00370
00371 }
00372
00373 }
00374
00375 #endif
00376
00378 if ( !mSpinSort ) {
00379 std::cout << "Problem detected, spin sorting disabled" << std::endl;
00380 continue;
00381 }
00382
00383 spin4=getSpinState( mMuDst->muDst(),bxStar );
00384 if ( spin4>=0 )
00385 if ( mymap[spin4] ) {
00386 mHistograms[ mymap[spin4] ] -> Fill( pair );
00387 }
00388
00389 hMassBx->Fill(pair.mass(),bxStar);
00390 if ( pair.mass()>0.1 && pair.mass() < 0.18 ) {
00391 hBxStarPi0->Fill(bxStar);
00392 hZvertBx->Fill(pair.vertex().Z(),bxStar);
00393 hZggBx->Fill(pair.zgg(),bxStar);
00394 hEtaBx->Fill(pair.momentum().Eta(),bxStar);
00395 }
00396
00397
00398 }
00399 hPi0Counter->Fill(numoacceptedpair);
00400
00401
00403 for ( Int_t i=0;i<mEEmixer->numberOfMixedCandidates();i++ )
00404 {
00405
00406 StEEmcIUPair pair = mEEmixer->mixedCandidate(i);
00407 if ( !accept( pair, false ) ) continue;
00408 mBackgrounds[ kAny ] -> Fill( pair );
00409 mMixedEvent -> addPair( pair );
00410
00411 }
00412
00413 END_OF_MAKE:
00414 mRealTree->Fill();
00415 mMixedTree->Fill();
00416 return kStOK;
00417 }
00418
00419
00420 void StEEmcIUPi0Analysis::mixer(const Char_t *name)
00421 {
00422 mEEmixer=(StEEmcIUMixMaker *)GetMaker(name);
00423 assert(mEEmixer);
00424 }
00425
00426
00427 void StEEmcIUPi0Analysis::points(const Char_t *name)
00428 {
00429 mEEpoints=(StEEmcIUPointMaker *)GetMaker(name);
00430 assert(mEEpoints);
00431 }
00432
00433
00434 void StEEmcIUPi0Analysis::mudst(const Char_t *name)
00435 {
00436 mMuDst=(StMuDstMaker*)GetMaker(name);
00437 assert(mMuDst);
00438 }
00439
00440 void StEEmcIUPi0Analysis::analysis(const Char_t *name)
00441 {
00442 mEEanalysis=(StEEmcA2EMaker*)GetMaker(name);
00443 assert(mEEanalysis);
00444 }
00445
00446 void StEEmcIUPi0Analysis::spin(const Char_t *name)
00447 {
00448 mSpinDb=(StSpinDbMaker*)GetMaker(name);
00450 }
00451
00452 Bool_t StEEmcIUPi0Analysis::twoBodyCut( StEEmcIUPair &pair )
00453 {
00454
00455 static const Int_t maxPerCluster = 4;
00456 static const Float_t minTowerFrac = 0.;
00457
00459
00460 Bool_t towers[720]; for (Int_t i=0;i<720;i++ ) towers[i]=false;
00461 StEEmcTower t1 = pair.point(0).tower(0);
00462 StEEmcTower t2 = pair.point(1).tower(0);
00463
00464 Float_t Etowers=0.;
00465 Etowers += t1.energy();
00466
00467
00468 towers[ t1.index() ] = true;
00469 towers[ t2.index() ] = true;
00470
00471 for ( Int_t i=0;i<t1.numberOfNeighbors();i++ ) {
00472 StEEmcTower mytow=t1.neighbor(i);
00473 towers[ mytow.index() ] = true;
00474 Etowers += mytow.energy();
00475 }
00476 for ( Int_t i=0;i<t2.numberOfNeighbors();i++ ) {
00477 StEEmcTower mytow=t2.neighbor(i);
00478 if ( !towers[mytow.index()] ) Etowers += mytow.energy();
00479 towers[ mytow.index() ] = true;
00480 }
00481
00482
00483
00485
00488
00489 Int_t count = 0;
00490 Float_t pe1=pair.point(0).energy();
00491 Float_t pe2=pair.point(1).energy();
00492 Float_t min2=(pe1<=pe2)?pe1:pe2;
00493 Bool_t Most2E=true;
00494
00495 for ( Int_t i=0;i<mEEpoints->numberOfPoints();i++ )
00496 {
00497 StEEmcIUPoint p=mEEpoints->point(i);
00498 Float_t ppe=p.energy();
00499 StEEmcTower t=p.tower(0);
00500 if ( towers[ t.index() ] ){
00501 count++;
00502 if(t.index()!=t1.index() && t.index()!=t2.index() && ppe>min2){
00503 Most2E=false;
00504 }
00505 }
00506 }
00507
00508 Most2E=true;
00509
00510 Float_t Epoints = pair.energy();
00511
00512 return ( count <= maxPerCluster && Most2E && Epoints > minTowerFrac * Etowers );
00513
00514 }
00515
00516
00517
00518 Bool_t StEEmcIUPi0Analysis::accept( StMuDst *mudst )
00519 {
00520
00521
00522 StMuEvent *event=mudst->event();
00523 assert(event);
00524
00525 Bool_t good_trigger = false;
00526
00527
00528
00529
00530
00531
00532 StMuTriggerIdCollection tic = event -> triggerIdCollection();
00533 StTriggerId l1trig = tic.l1();
00534
00535 if ( mTriggerList.size() <=0 ) {
00536 good_trigger = true;
00537 }
00538
00539 Int_t go=0;
00540 std::vector<Int_t>::iterator iter=mTriggerList.begin();
00541 while ( iter != mTriggerList.end() ) {
00542 go = l1trig.isTrigger( (*iter) );
00543 good_trigger |= go;
00544 iter++;
00545 }
00546
00547
00548
00549
00550
00551
00552 if ( mTrigSim )
00553 {
00554 good_trigger = mTrigSim -> getEEMCtrigHT( mTrigSimThreshold );
00555 good_trigger &= mTrigSim -> getBBCtrig();
00556 }
00557
00558
00559 if ( l1trig.isTrigger( mMinBias ) )
00560 hEventCounter -> Fill( kMinBias );
00561
00562
00563 if ( good_trigger )
00564 {
00565
00566 hEventCounter -> Fill( kTrig );
00567 }
00568 else
00569 {
00570
00571 return false;
00572 }
00573
00574
00575
00576 StEEmcTower ht;
00577 Bool_t sw36=0;
00578
00579 for ( Int_t i=0;i<mEEanalysis->numberOfHitTowers(0);i++ )
00580 {
00581 StEEmcTower t=mEEanalysis->hittower(i,0);
00582 if ( !t.fail() && t.et() > ht.et() ) {
00583 ht=t;
00584 if ( ht.et()>=mCuts->tower_et_cut &&
00585 ht.sector()>=2&&ht.sector()<=5 ) sw36=true;
00586 }
00587 }
00588
00589 if ( ht.et() < mCuts->tower_et_cut ) {
00590 good_trigger = false;
00591
00592 return false;
00593 }
00594 else
00595 {
00596 hEventCounter -> Fill( kSoftTrig );
00597 }
00598 if ( sw36 ) hEventCounter -> Fill( kSoftTrig36 );
00599
00600
00601 StThreeVectorF vert=event->primaryVertexPosition();
00602 if ( vert.z() < mCuts->z_vertex_min) {
00603
00604 return false;
00605 }
00606 if ( vert.z() > mCuts->z_vertex_max ) {
00607
00608 return false;
00609 }
00610
00611 hEventCounter -> Fill( kVertex );
00612
00613
00614 return true;
00615
00616 }
00617
00618 Bool_t StEEmcIUPi0Analysis::accept( StEEmcIUPair pair, Bool_t fill )
00619 {
00620
00622 if ( !twoBodyCut(pair) ) return false;
00623 if (fill) hPairCounter -> Fill( kTwoBody );
00625 StEEmcIUPoint point1=pair.point(0);
00626 StEEmcIUPoint point2=pair.point(1);
00627 StEEmcTower tower1=pair.point(0).tower(0);
00628 StEEmcTower tower2=pair.point(1).tower(0);
00629 TVector3 pointpos1=point1.position();
00630 TVector3 pointpos2=point2.position();
00631 Float_t peta1=pointpos1.Eta();
00632 Float_t peta2=pointpos2.Eta();
00633 Float_t pphi1=pointpos1.Phi();
00634 Float_t pphi2=pointpos2.Phi();
00635 Int_t mod_phi1=int(pphi1*180./3.14159265+180.)%6;
00636 Int_t mod_phi2=int(pphi2*180./3.14159265+180.)%6;
00637 pointpos1=pointpos1.Unit();
00638 pointpos2=pointpos2.Unit();
00639
00640 Float_t pet1 = (point1.energy()*pointpos1).Perp();
00641 Float_t pet2 = (point2.energy()*pointpos2).Perp();
00642
00643
00644 if(pet1<=1.5 || pet2<=1.5) return false;
00645
00646 if ( tower1.adc() < mCuts->adc_cut && tower2.adc() < mCuts->adc_cut ) return false;
00647
00648 if ( fill) hPairCounter -> Fill( kAdcTow );
00649
00650
00651
00652
00653
00654 TVector3 d1 = mEEgeom -> getTowerCenter( (UInt_t)tower1.sector(), (UInt_t)tower1.subsector(), (UInt_t)tower1.etabin() );
00655 TVector3 d2 = mEEgeom -> getTowerCenter( (UInt_t)tower2.sector(), (UInt_t)tower2.subsector(), (UInt_t)tower2.etabin() );
00656
00657 TVector3 dd1(d1.x(),d1.y(),d1.z()-pair.vertex().z());
00658 TVector3 dd2(d2.x(),d2.y(),d2.z()-pair.vertex().z());
00659
00660
00661 d1=d1.Unit();
00662 d2=d2.Unit();
00663 dd1=dd1.Unit();
00664 dd2=dd2.Unit();
00665
00666
00667 Float_t et11 = (tower1.energy()*dd1).Perp();
00668 Float_t et22 = (tower2.energy()*dd2).Perp();
00669
00670
00671 Float_t ueta1,ueta2,tower_et_cut1=3,tower_et_cut2=3;
00672 if(peta1>=1.901 && peta1<=2) ueta1=1.932;
00673 if(peta1>=1.807 && peta1<1.901) ueta1=1.845;
00674 if(peta1>=1.718 && peta1<1.807) ueta1=1.761;
00675 if(peta1>=1.633 && peta1<1.718) ueta1=1.664;
00676 if(peta1>=1.552 && peta1<1.633) ueta1=1.585;
00677 if(peta1>=1.476 && peta1<1.552) ueta1=1.507;
00678 if(peta1>=1.403 && peta1<1.476) ueta1=1.4374;
00679 if(peta1>=1.334 && peta1<1.403) ueta1=1.36;
00680 if(peta1>=1.268 && peta1<1.334) ueta1=1.294;
00681 if(peta1>=1.205 && peta1<1.268) ueta1=1.2315;
00682 if(peta1>=1.146 && peta1<1.205) ueta1=1.174;
00683 if(peta1>=1.089 && peta1<1.146) ueta1=1.1265;
00684
00685 if(peta2>=1.901 && peta2<=2) ueta2=1.932;
00686 if(peta2>=1.807 && peta2<1.901) ueta2=1.845;
00687 if(peta2>=1.718 && peta2<1.807) ueta2=1.761;
00688 if(peta2>=1.633 && peta2<1.718) ueta2=1.664;
00689 if(peta2>=1.552 && peta2<1.633) ueta2=1.585;
00690 if(peta2>=1.476 && peta2<1.552) ueta2=1.507;
00691 if(peta2>=1.403 && peta2<1.476) ueta2=1.4374;
00692 if(peta2>=1.334 && peta2<1.403) ueta2=1.36;
00693 if(peta2>=1.268 && peta2<1.334) ueta2=1.294;
00694 if(peta2>=1.205 && peta2<1.268) ueta2=1.2315;
00695 if(peta2>=1.146 && peta2<1.205) ueta2=1.174;
00696 if(peta2>=1.089 && peta2<1.146) ueta2=1.1265;
00697 float toweretcut=mCuts->tower_et_cut;
00698 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))));
00699 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))));
00700
00701 if ( et11 < tower_et_cut1 && et22 < tower_et_cut2 ) return false;
00702
00703
00704
00705 if ( fill ) hPairCounter -> Fill( kEtTow );
00706
00708
00709 float Rxy=TMath::Sqrt(pair.vertex().x()*pair.vertex().x()+pair.vertex().y()*pair.vertex().y());
00710 float hHeight=pair.pt()*(270.0-pair.vertex().Z())/pair.pz()+Rxy;
00711 float etatheta=TMath::ATan(hHeight/270.0);
00712
00713 float mideta=tan(etatheta/2.0);
00714 float eta=-log(mideta);
00715 if ( eta < mCuts->eta_min || eta > mCuts->eta_max ) return false;
00716 if ( fill ) hPairCounter->Fill( kKinematics );
00717
00718
00719
00720
00721
00722
00723 return true;
00724
00725 }
00726
00727
00728 Int_t StEEmcIUPi0Analysis::getSpinState( StMuDst *mudst, Int_t &bxStar )
00729 {
00730
00731 StMuEvent *event = mudst -> event();
00732 StL0Trigger *trig =&(event->l0Trigger());
00733
00734 StMuTriggerIdCollection tic = event -> triggerIdCollection();
00735 StTriggerId l1trig = tic.l1();
00736
00737
00738 Int_t bx48 = trig->bunchCrossingId();
00739 Int_t bx7 = trig->bunchCrossingId7bit( mRunNumber );
00740
00741
00742 bxStar = mSpinDb->BXyellowUsingBX48(bx48);
00743
00744 mRealEvent->bx48 = bx48;
00745 mRealEvent->bx7 = bx7;
00746 mRealEvent->bxStar = bxStar;
00747
00748 mMixedEvent->bx48=bx48;
00749 mMixedEvent->bx7 =bx7;
00750 mMixedEvent->bxStar=bxStar;
00751
00754 if ( bx7 == 0 || bx7 == 119 ) return -1;
00755 if ( bx7 == 14 || bx7 == 54 ) return -1;
00757
00758 hBx7->Fill( bx7 );
00759 hBx48->Fill( bx48 );
00760 hBx7diffBx48->Fill( bx7, mSpinDb->offsetBX48minusBX7(bx48,bx7) );
00761
00762
00763
00764 if ( mSpinDb -> isMaskedUsingBX48(bx48) ) return -1;
00765 if ( mSpinDb->offsetBX48minusBX7(bx48,bx7)!=0 ) std::cout << "BUNCH CROSSINGS INCONSISTENT" << std::endl << std::flush;
00766
00767
00768
00769
00770 if ( mSpinDb->offsetBX48minusBX7(bx48,bx7)!=0 )
00771 {
00772 mSpinSort = false;
00773 for ( Int_t ii=1;ii<=4;ii++ ) mHistograms[ii]->Clear();
00774 }
00775
00776
00777 Int_t spin4 = mSpinDb->spin4usingBX48(bx48);
00778 hSpin4->Fill(spin4);
00779
00780 if ( l1trig.isTrigger(96011) ) hSpin4mb -> Fill(spin4);
00781
00782
00783 return spin4;
00784
00785 }
00786
00787 void StEEmcIUPi0Analysis::Clear(Option_t *opts)
00788 {
00789 mRealEvent -> Clear();
00790 mMixedEvent -> Clear();
00791 }
