#include <DiGam.h>
Public Member Functions | |
| DiGam (const std::string &name, ISvcLocator *pSvcLocator) | |
| DiGam (const std::string &name, ISvcLocator *pSvcLocator) | |
| StatusCode | execute () |
| StatusCode | execute () |
| StatusCode | finalize () |
| StatusCode | finalize () |
| StatusCode | initialize () |
| StatusCode | initialize () |
Private Attributes | |
| int | boost_signal |
| int | boost_tot |
| double | boostMinEmax |
| double | boostMinEmin |
| double | boostPhimax |
| double | boostPhimin |
| double | CrossSection |
| TH1F * | h_lum |
| TH1F * | h_lum |
| double | jpsiCrossSection |
| double | jpsiMCEff |
| double | jpsiMCEffBoost |
| NTuple::Item< double > | m_angle |
| NTuple::Item< double > | m_angle |
| NTuple::Item< double > | m_boostAngle |
| NTuple::Item< double > | m_boostAngle |
| NTuple::Item< double > | m_boostDelPhi |
| NTuple::Item< double > | m_boostDelPhi |
| NTuple::Item< double > | m_boostDelTheta |
| NTuple::Item< double > | m_boostDelTheta |
| NTuple::Item< double > | m_boostIM |
| NTuple::Item< double > | m_boostIM |
| NTuple::Item< double > | m_boostM |
| NTuple::Item< double > | m_boostM |
| NTuple::Item< double > | m_boostMaxE |
| NTuple::Item< double > | m_boostMaxE |
| NTuple::Item< double > | m_boostMinE |
| NTuple::Item< double > | m_boostMinE |
| NTuple::Item< double > | m_delPhi |
| NTuple::Item< double > | m_delPhi |
| NTuple::Item< double > | m_delTheta |
| NTuple::Item< double > | m_delTheta |
| NTuple::Item< long > | m_idxmc |
| NTuple::Item< long > | m_idxmc |
| NTuple::Item< double > | m_IM |
| NTuple::Item< double > | m_IM |
| NTuple::Item< double > | m_m |
| NTuple::Item< double > | m_m |
| NTuple::Item< double > | m_maxE |
| NTuple::Item< double > | m_maxE |
| NTuple::Item< double > | m_maxPhi |
| NTuple::Item< double > | m_maxPhi |
| NTuple::Item< double > | m_maxTheta |
| NTuple::Item< double > | m_maxTheta |
| NTuple::Item< double > | m_minE |
| NTuple::Item< double > | m_minE |
| NTuple::Item< double > | m_minPhi |
| NTuple::Item< double > | m_minPhi |
| NTuple::Item< double > | m_minTheta |
| NTuple::Item< double > | m_minTheta |
| NTuple::Array< long > | m_motheridx |
| NTuple::Array< long > | m_motheridx |
| NTuple::Array< long > | m_pdgid |
| NTuple::Array< long > | m_pdgid |
| NTuple::Item< long > | m_rec |
| NTuple::Item< long > | m_rec |
| NTuple::Item< long > | m_run |
| NTuple::Item< long > | m_run |
| ITHistSvc * | m_thsvc |
| ITHistSvc * | m_thsvc |
| NTuple::Item< double > | m_tot |
| NTuple::Item< double > | m_tot |
| NTuple::Tuple * | m_tuple2 |
| NTuple::Tuple * | m_tuple2 |
| double | MCEff |
| double | MCEffBoost |
| double | psi2sCrossSection |
| double | psi2sMCEff |
| double | psi2sMCEffBoost |
| double | psi3770CrossSection |
| double | psi3770MCEff |
| double | psi3770MCEffBoost |
| int | RunModel |
| int | runNo |
| int | signal |
| int | tot |
|
||||||||||||
|
00038 :
00039 Algorithm(name, pSvcLocator) {
00040 declareProperty("jpsiCrossSection",jpsiCrossSection);
00041 declareProperty("jpsiMCEff",jpsiMCEff);
00042 declareProperty("jpsiMCEffBoost",jpsiMCEffBoost);
00043
00044 declareProperty("psi2sCrossSection",psi2sCrossSection);
00045 declareProperty("psi2sMCEff",psi2sMCEff);
00046 declareProperty("psi2sMCEffBoost",psi2sMCEffBoost);
00047
00048 declareProperty("psi3770CrossSection",psi3770CrossSection);
00049 declareProperty("psi3770MCEff",psi3770MCEff);
00050 declareProperty("psi3770MCEffBoost",psi3770MCEffBoost);
00051
00052 declareProperty("RunModel",RunModel=1);
00053 }
//***************************************************************************************
|
|
||||||||||||
|
|
|
|
|
|
|
00139 {
00140 MsgStream log(msgSvc(), name());
00141 log << MSG::INFO << "in execute()" << endreq;
00142 N0++;
00143 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00144 runNo=eventHeader->runNumber();
00145 int event=eventHeader->eventNumber();
00146 log<<MSG::INFO<<"run,evtnum="<<runNo<<","<<event<<endreq;
00147 m_run = eventHeader->runNumber();
00148 m_rec = eventHeader->eventNumber();
00149
00150 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00151 log<<MSG::INFO<<"ncharg,nneu,tottks="
00152 <<evtRecEvent->totalCharged()<<","
00153 <<evtRecEvent->totalNeutral()<<","
00154 <<evtRecEvent->totalTracks()<<endreq;
00155 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00156 std::vector<int> iGam;
00157 iGam.clear();
00158 N1++;
00159 for(int i=evtRecEvent->totalCharged();i< evtRecEvent->totalTracks();i++){
00160 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00161 if(!(*itTrk)->isEmcShowerValid()) continue;
00162 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00163 if(emcTrk->energy()<0.6)continue;
00164 if(fabs(cos(emcTrk->theta()))>0.83)continue;
00165 iGam.push_back((*itTrk)->trackId());
00166 }
00167 if(iGam.size()<2)return StatusCode::SUCCESS;
00168 N2++;
00169 double energy1=0.5;
00170 double energy2=0.5;
00171 int gam1=-9999;
00172 int gam2=-9999;
00173 for(int i=0;i<iGam.size();i++){
00174 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + iGam[i];
00175 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00176 //std::cout<<emcTrk->energy()<<std::endl;
00177 if(emcTrk->energy()>energy1){ //max emc
00178 energy2=energy1;
00179 gam2=gam1;
00180 energy1=emcTrk->energy();
00181 gam1=iGam[i];
00182 }
00183 else if(emcTrk->energy()>energy2){ //second max emc
00184 energy2=emcTrk->energy();
00185 gam2=iGam[i];
00186 }
00187 }
00188 if(gam1==-9999 || gam2==-9999)return StatusCode::SUCCESS;
00189 N3++;
00190 m_tot=energy1+energy2;
00191 RecEmcShower* maxEmc=(*(evtRecTrkCol->begin()+gam1))->emcShower();
00192 RecEmcShower* minEmc=(*(evtRecTrkCol->begin()+gam2))->emcShower();
00193 m_maxE=maxEmc->energy();
00194 m_minE=minEmc->energy();
00195 m_maxTheta=maxEmc->theta();
00196 m_minTheta=minEmc->theta();
00197 m_maxPhi=maxEmc->phi();
00198 m_minPhi=minEmc->phi();
00199 m_delPhi=(fabs(m_maxPhi-m_minPhi)-pai)*180./pai;
00200 m_delTheta=(fabs(m_maxTheta+m_minTheta)-pai)*180./pai;
00201
00202 HepLorentzVector max,min;
00203 max.setPx(m_maxE*sin(m_maxTheta)*cos(m_maxPhi));
00204 max.setPy(m_maxE*sin(m_maxTheta)*sin(m_maxPhi));
00205 max.setPz(m_maxE*cos(m_maxTheta));
00206 max.setE(m_maxE);
00207 min.setPx(m_minE*sin(m_minTheta)*cos(m_minPhi));
00208 min.setPy(m_minE*sin(m_minTheta)*sin(m_minPhi));
00209 min.setPz(m_minE*cos(m_minTheta));
00210 min.setE(m_minE);
00211 m_angle=max.angle(min.vect())*180./pai;
00212 m_m=(max+min).m();
00213 m_IM=max.invariantMass(min);
00214 //select signal and sideband to get lum;
00215 if(m_minE>=boostMinEmin && m_delPhi>-7 && m_delPhi<5 && m_minE<=boostMinEmax)tot++;
00216 if(m_minE>=boostMinEmin && m_delPhi>-4 && m_delPhi<2 && m_minE<=boostMinEmax)signal++;
00217
00218 //boost
00219 //HepLorentzVector p_cms(0.034067,0.0,0.0,3.097);
00220 HepLorentzVector boost1=max.boost(-0.011,0,0);
00221 HepLorentzVector boost2=min.boost(-0.011,0,0);
00222 m_boostAngle=boost1.angle(boost2.vect())*180./pai;
00223 m_boostMaxE=boost1.e();
00224 m_boostMinE=boost2.e();
00225 m_boostDelPhi=(fabs(boost1.phi()-boost2.phi())-pai)*180./pai;
00226 m_boostDelTheta=(fabs(boost1.theta()+boost2.theta())-pai)*180./pai;
00227 m_boostM=(boost1+boost2).m();
00228 m_boostIM=boost1.invariantMass(boost2);
00229 log << MSG::INFO << "num Good Photon " << iGam.size()<<endreq;
00230 if(m_boostMinE>=boostMinEmin && m_boostMinE<=boostMinEmax && fabs(m_boostDelPhi)<=boostPhimin)boost_signal++;
00231 if(m_boostMinE>=boostMinEmin && m_boostMinE<=boostMinEmax && fabs(m_boostDelPhi)<=boostPhimax)boost_tot++;
00232 //mcTruth
00233 if (eventHeader->runNumber()<0)
00234 {
00235 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(), "/Event/MC/McParticleCol");
00236 int m_numParticle = 0;
00237 if (!mcParticleCol)
00238 {
00239 //std::cout << "Could not retrieve McParticelCol" << std::endl;
00240 return StatusCode::FAILURE;
00241 }
00242 else
00243 {
00244 bool jpsipDecay = false;
00245 int rootIndex = -1;
00246 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
00247 for (; iter_mc != mcParticleCol->end(); iter_mc++)
00248 {
00249 if ((*iter_mc)->primaryParticle()) continue;
00250 if (!(*iter_mc)->decayFromGenerator()) continue;
00251
00252 if ((*iter_mc)->particleProperty()==443)
00253 {
00254 jpsipDecay = true;
00255 rootIndex = (*iter_mc)->trackIndex();
00256 }
00257 if (!jpsipDecay) continue;
00258 int mcidx = ((*iter_mc)->mother()).trackIndex() - rootIndex;
00259 int pdgid = (*iter_mc)->particleProperty();
00260 m_pdgid[m_numParticle] = pdgid;
00261 m_motheridx[m_numParticle] = mcidx;
00262 m_numParticle += 1;
00263 }
00264 }
00265 m_idxmc = m_numParticle;
00266 }
00267 N4++;
00268 m_tuple2->write();
00269 return StatusCode::SUCCESS;
00270 }
|
|
|
|
|
|
00272 {
00273 MsgStream log(msgSvc(), name());
00274 log << MSG::INFO << "in finalize()" << endmsg;
00275 //get intLumEndcapEE
00276 double ssg=0.;
00277 double lum=0.;
00278 double boost_ssg=0.;
00279 double boost_lum=0.;
00280 if(CrossSection*MCEff==0||CrossSection*MCEffBoost==0){
00281 log << MSG::FATAL<<"DiGam Error: cross_section or MC_eff is not correct!"<<endreq;
00282 exit(1);
00283 }
00284 ssg=(signal-(tot-signal))*1.0;
00285 //boss 650 mc:7.135%
00286 lum=(ssg)/(CrossSection*MCEff);
00287 //boost LUM
00288 boost_ssg=(boost_signal-(boost_tot-boost_signal))*1.0;
00289 //boost 650 mc:7.097%
00290 boost_lum=(boost_ssg)/(CrossSection*MCEffBoost);
00291 std::cout<<"minE,phi"<<boostMinEmin<<","<<boostMinEmax<<","<<boostPhimin<<","<<boostPhimax<<std::endl;
00292 std::cout<<"zhaohs is:"<<CrossSection<<","<<MCEffBoost<<std::endl;
00293 std::cout<<"zhaohs generall is:"<<CrossSection<<","<<MCEff<<std::endl;
00294 std::cout<<"ssg lum boost boost_lum is:"<<ssg<<","<<lum<<","<<boost_ssg<<","<<boost_lum<<std::endl;
00295 h_lum->SetBinContent(1,lum);
00296 h_lum->SetBinContent(2,ssg);
00297 h_lum->SetBinContent(3,boost_lum);
00298 h_lum->SetBinContent(4,boost_ssg);
00299 if(m_thsvc->regHist("/DQAHist/zhsLUM/lum", h_lum).isFailure()){
00300 log << MSG::FATAL<<"DiGam Error:can't write data to LUM!"<<endreq;
00301 exit(1);
00302 }
00303 return StatusCode::SUCCESS;
00304 }
|
|
|
|
|
|
00055 {
00056 MsgStream log(msgSvc(), name());
00057 log << MSG::INFO << "in initialize()" << endmsg;
00058
00059 if(RunModel==1){//jpsi
00060 CrossSection=jpsiCrossSection;
00061 MCEff=jpsiMCEff;
00062 MCEffBoost=jpsiMCEffBoost;
00063 boostPhimin=2.5;
00064 boostPhimax=5;
00065 boostMinEmin=1.2;
00066 boostMinEmax=1.6;
00067 }
00068 else if(RunModel==2){//psi2s
00069 CrossSection=psi2sCrossSection;
00070 MCEff=psi2sMCEff;
00071 MCEffBoost=psi2sMCEffBoost;
00072 boostPhimin=2.5;
00073 boostPhimax=5;
00074 boostMinEmin=1.4;
00075 boostMinEmax=1.9;
00076 }
00077 else if(RunModel==3){//psi3770
00078 CrossSection=psi3770CrossSection;
00079 MCEff=psi3770MCEff;
00080 MCEffBoost=psi3770MCEffBoost;
00081 boostPhimin=2.5;
00082 boostPhimax=5;
00083 boostMinEmin=1.4;
00084 boostMinEmax=2;
00085 }
00086
00087 tot=0;
00088 signal=0;
00089 boost_signal=0;
00090 boost_tot=0;
00091 StatusCode ssc=service("THistSvc", m_thsvc);
00092 if(ssc.isFailure()) {
00093 log << MSG::FATAL << "DiGamAlg:Couldn't get THistSvc" << endreq;
00094 exit(1);
00095 }
00096 h_lum=new TH1F("lum","lum",4,0.5,4.5);
00097 StatusCode status;
00098 NTuplePtr nt2(ntupleSvc(),"DQAFILE/zhsDiGam");
00099 if(nt2)m_tuple2=nt2;
00100 else{
00101 m_tuple2=ntupleSvc()->book("DQAFILE/zhsDiGam",CLID_ColumnWiseTuple,"ks N-Tuple example");
00102 if(m_tuple2){
00103 status=m_tuple2->addItem("ETol",m_tot);
00104 status=m_tuple2->addItem("maxE",m_maxE);
00105 status=m_tuple2->addItem("minE",m_minE);
00106 status=m_tuple2->addItem("maxTheta",m_maxTheta);
00107 status=m_tuple2->addItem("minTheta",m_minTheta);
00108 status=m_tuple2->addItem("maxPhi",m_maxPhi);
00109 status=m_tuple2->addItem("minPhi",m_minPhi);
00110 status=m_tuple2->addItem("delTheta",m_delTheta);
00111 status=m_tuple2->addItem("delPhi",m_delPhi);
00112 status=m_tuple2->addItem("angle",m_angle);
00113 status=m_tuple2->addItem("boostAngle",m_boostAngle);
00114 status=m_tuple2->addItem("boostMaxE",m_boostMaxE);
00115 status=m_tuple2->addItem("boostMinE",m_boostMinE);
00116 status=m_tuple2->addItem("boostDelPhi",m_boostDelPhi);
00117 status=m_tuple2->addItem("boostDelTheta",m_boostDelTheta);
00118 status=m_tuple2->addItem("boostM",m_boostM);
00119 status=m_tuple2->addItem("boostIM",m_boostIM);
00120 status=m_tuple2->addItem("m",m_m);
00121 status=m_tuple2->addItem("IM",m_IM);
00122
00123 status=m_tuple2->addItem ("run",m_run );
00124 status=m_tuple2->addItem ("rec",m_rec );
00125 status=m_tuple2->addItem("indexmc", m_idxmc, 0, 100);
00126 status=m_tuple2->addIndexedItem("pdgid", m_idxmc, m_pdgid);
00127 status=m_tuple2->addIndexedItem("motheridx", m_idxmc, m_motheridx);
00128
00129 }
00130 else {
00131 log<<MSG::FATAL<<"Cannot book N-tuple:"<<long(m_tuple2)<<endmsg;
00132 return StatusCode::FAILURE;
00133 }
00134 }
00135
00136 return StatusCode::SUCCESS;
00137 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1