00001 #include "GaudiKernel/MsgStream.h"
00002 #include "GaudiKernel/AlgFactory.h"
00003 #include "GaudiKernel/ISvcLocator.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/IDataProviderSvc.h"
00006 #include "GaudiKernel/PropertyMgr.h"
00007 #include "EventModel/EventModel.h"
00008 #include "EventModel/Event.h"
00009 #include "EvtRecEvent/EvtRecEvent.h"
00010 #include "EvtRecEvent/EvtRecTrack.h"
00011
00012 #include "DstEvent/TofHitStatus.h"
00013 #include "EventModel/EventHeader.h"
00014
00015
00016 #include "TMath.h"
00017 #include "GaudiKernel/INTupleSvc.h"
00018 #include "GaudiKernel/NTuple.h"
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/IHistogramSvc.h"
00021 #include "CLHEP/Vector/ThreeVector.h"
00022 using CLHEP::Hep3Vector;
00023 #include "CLHEP/Vector/LorentzVector.h"
00024 using CLHEP::HepLorentzVector;
00025 #include "CLHEP/Vector/TwoVector.h"
00026 using CLHEP::Hep2Vector;
00027 #include "CLHEP/Geometry/Point3D.h"
00028 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00029 typedef HepGeom::Point3D<double> HepPoint3D;
00030 #endif
00031 #include "VertexFit/KinematicFit.h"
00032 #include "VertexFit/VertexFit.h"
00033 #include "VertexFit/SecondVertexFit.h"
00034 #include "ParticleID/ParticleID.h"
00035 #include "TwoGammaAlg/TwoGamma.h"
00036 #include <vector>
00037 #include <fstream>
00038 typedef std::vector<int> Vint;
00039 typedef std::vector<HepLorentzVector> Vp4;
00040 const double ecms = 3.686;
00041 const double velc = 299.792458;
00042 const double xmass[5] = {0.000511, 0.105658, 0.139570, 0.493677, 0.938272};
00043 const double pai=3.1415926;
00044
00046 TwoGamma::TwoGamma(const std::string& name, ISvcLocator* pSvcLocator) :
00047 Algorithm(name, pSvcLocator)
00048 {
00049
00050 for(int i = 0; i < 10; i++)
00051 {
00052 m_pass[i] = 0;
00053 }
00054 m_event = 0;
00055
00056 Ndata1 = 0;
00057 Ndata2 = 0;
00058 m_runNo = 0;
00059
00060
00061 declareProperty("CmsEnergy", m_ecms = 3.686);
00062
00063 declareProperty("max1", m_max1 = 1.2);
00064 declareProperty("max2", m_max2 = 0.8);
00065 declareProperty("costheta", m_costheta = 0.8);
00066 declareProperty("dphi1", m_dphi1 = 2.5);
00067 declareProperty("dphi2", m_dphi2 = 5);
00068 declareProperty("eff", m_eff = 0.07);
00069 declareProperty("sec", m_sec = 184.1);
00070 }
00071
00072
00073 StatusCode TwoGamma::initialize()
00074 {
00075 MsgStream log(msgSvc(), name());
00076
00077 log << MSG::INFO << "in initialize()" << endmsg;
00078
00079 StatusCode status;
00080
00081 NTuplePtr nt1(ntupleSvc(), "FILE1/DiGam");
00082 if ( nt1 ) m_tuple1 = nt1;
00083 else
00084 {
00085 m_tuple1 = ntupleSvc()->book ("FILE1/DiGam", CLID_ColumnWiseTuple, "DiGam N-Tuple signal");
00086 if ( m_tuple1 )
00087 {
00088 status = m_tuple1->addItem ("run", m_run );
00089 status = m_tuple1->addItem ("rec", m_rec );
00090 status = m_tuple1->addItem ("time", m_time );
00091
00092 status = m_tuple1->addItem ("ngood", m_ngood);
00093 status = m_tuple1->addItem ("nchrg", m_nchrg);
00094
00095 status = m_tuple1->addItem ("Cms_e1", m_e1);
00096 status = m_tuple1->addItem ("Cms_e2", m_e2);
00097 status = m_tuple1->addItem ("Cms_e", m_e);
00098 status = m_tuple1->addItem ("Cms_costheta1", m_costheta1);
00099 status = m_tuple1->addItem ("Cms_costheta2", m_costheta2);
00100 status = m_tuple1->addItem ("Cms_dltphi", m_dltphi);
00101 status = m_tuple1->addItem ("Cms_dltphi_1", m_dltphi_1);
00102 status = m_tuple1->addItem ("Cms_dlttheta", m_dlttheta);
00103 status = m_tuple1->addItem ("Cms_phi1", m_phi1);
00104 status = m_tuple1->addItem ("Cms_phi2", m_phi2);
00105
00106 status = m_tuple1->addItem ("Lab_e1", m_e1_lab);
00107 status = m_tuple1->addItem ("Lab_e2", m_e2_lab);
00108 status = m_tuple1->addItem ("Lab_e", m_e_lab);
00109 status = m_tuple1->addItem ("Lab_costheta1", m_costheta1_lab);
00110 status = m_tuple1->addItem ("Lab_costheta2", m_costheta2_lab);
00111 status = m_tuple1->addItem ("Lab_dltphi", m_dltphi_lab);
00112 status = m_tuple1->addItem ("Lab_dlttheta", m_dlttheta_lab);
00113 status = m_tuple1->addItem ("Lab_phi1", m_phi1_lab);
00114 status = m_tuple1->addItem ("Lab_phi2", m_phi2_lab);
00115
00116 status = m_tuple1->addItem ("xBoost", m_xBoost);
00117 status = m_tuple1->addItem ("yBoost", m_yBoost);
00118 status = m_tuple1->addItem ("zBoost", m_zBoost);
00119 }
00120 else
00121 {
00122 log << MSG::ERROR << "Cannot book N-tuple1:"<<long(m_tuple1)<<endmsg;
00123 return StatusCode::FAILURE;
00124 }
00125 }
00126 Ndata1 = 0;
00127 Ndata2 = 0;
00128 m_runNo = 0;
00129
00130 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00131 return StatusCode::SUCCESS;
00132 }
00133
00134
00135 StatusCode TwoGamma::execute()
00136 {
00137 StatusCode sc=StatusCode::SUCCESS;
00138 m_event++;
00139
00140 MsgStream log(msgSvc(),name());
00141 log<<MSG::INFO<<"in execute()"<<endreq;
00142 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00143 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(),EventModel::EvtRec::EvtRecEvent);
00144 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(),EventModel::EvtRec::EvtRecTrackCol);
00145
00146 m_run = eventHeader->runNumber();
00147 m_rec = eventHeader->eventNumber();
00148 int runNo = m_run;
00149 int event = m_rec;
00150 int time = eventHeader->time();
00151 if(m_event > 1 && runNo != m_runNo) return sc;
00152 m_runNo = runNo;
00153 m_time = time;
00154
00155 if(m_rec%1000==0)cout<<"Run "<<m_run<<" Event "<<m_rec<<endl;
00156
00157 HepLorentzVector p4psip(0.011*m_ecms, 0.0, 0.005, m_ecms);
00158 double psipBeta = (p4psip.vect()).mag()/(p4psip.e());
00159
00160 m_pass[0]+=1;
00161
00162 int NCharge = evtRecEvent->totalCharged();
00163 if(NCharge != 0)return sc;
00164
00165 m_pass[1]+=1;
00166
00167 double Emax1=0.0;
00168 double Emax2=0.0;
00169 int Imax[2];
00170
00171 if(((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )<2) || ((evtRecEvent->totalTracks() - evtRecEvent->totalCharged() )>15)) return sc;
00172 m_pass[2]+=1;
00173
00174 HepLorentzVector p4Gam_1;
00175 HepLorentzVector p4Gam_2;
00176
00177 for(int i=evtRecEvent->totalCharged(); i<evtRecEvent->totalTracks(); i++)
00178 {
00179 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00180 if(!(*itTrk)->isEmcShowerValid()) continue;
00181 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00182
00183 HepLorentzVector p4Gam;
00184 double Phi = emcTrk->phi();
00185 double Theta = emcTrk->theta();
00186 double Ener=emcTrk->energy();
00187
00188 p4Gam.setPx(Ener * sin(Theta) * cos(Phi));
00189 p4Gam.setPy(Ener * sin(Theta) * sin(Phi));
00190 p4Gam.setPz(Ener * cos(Theta));
00191 p4Gam.setE(Ener);
00192 p4Gam.boost(-0.011, 0 ,-0.005 * 1.0 / 3.686);
00193
00194 if(Ener>Emax2)
00195 {
00196 Emax2=Ener;
00197 Imax[1]=i;
00198 p4Gam_2 = p4Gam;
00199 }
00200 if(Ener>Emax1)
00201 {
00202 Emax2=Emax1;
00203 p4Gam_2 = p4Gam_1;
00204 Imax[1]=Imax[0];
00205 Emax1=Ener;
00206 p4Gam_1 = p4Gam;
00207 Imax[0]=i;
00208 }
00209 }
00210 m_e1_lab = Emax1;
00211 m_e2_lab = Emax2;
00212 m_e_lab = Emax1 + Emax2;
00213 log << MSG::INFO << "Emax1 = " << Emax1 <<"Emax2= "<<Emax2<< endreq;
00214 if(Emax1 < m_max1 || Emax2 < m_max2) return sc;
00215 m_pass[3]+=1;
00216
00217
00218 double emcphi[2],emctht[2];
00219 for(int i = 0;i < 2; i++)
00220 {
00221 if(i>=evtRecTrkCol->size()) break;
00222 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + Imax[i];
00223 if(!(*itTrk)->isEmcShowerValid()) continue;
00224 RecEmcShower *emcTrk = (*itTrk)->emcShower();
00225 emcphi[i]=emcTrk->phi();
00226 emctht[i]=emcTrk->theta();
00227 }
00228 double dltphi_lab = (fabs(emcphi[0] - emcphi[1]) - pai) * 180.0 / pai;
00229 double dlttheta_lab = (fabs(emctht[0] + emctht[1]) - pai) * 180.0 / pai;
00230 m_costheta1_lab = cos(emctht[0]);
00231 m_costheta2_lab = cos(emctht[1]);
00232 m_phi1_lab = emcphi[0] * 180.0 / pai;
00233 m_phi2_lab = emcphi[1] * 180.0 / pai;
00234 m_dlttheta_lab = dlttheta_lab;
00235 m_dltphi_lab = dltphi_lab;
00236
00237 if(fabs(m_costheta1_lab) > m_costheta || fabs(m_costheta2_lab) > m_costheta) return sc;
00238 m_pass[4]+=1;
00239
00240
00241
00242 double px1 = p4Gam_1.px();
00243 double py1 = p4Gam_1.py();
00244 double pz1 = p4Gam_1.pz();
00245 double pxy1 = sqrt(px1 * px1 + py1 * py1);
00246 double e1 = p4Gam_1.e();
00247
00248 double px2 = p4Gam_2.px();
00249 double py2 = p4Gam_2.py();
00250 double pz2 = p4Gam_2.pz();
00251 double pxy2 = sqrt(px2 * px2 + py2 * py2);
00252 double e2 = p4Gam_2.e();
00253
00254
00255 double phi1 = 0;
00256 double phi2 = 0;
00257 if(atan(py1 * 1.0 / px1) > 0){
00258 if(px1 > 0)phi1 = atan(py1 * 1.0 / px1);
00259 else phi1 = -(pai - atan(py1 * 1.0 / px1));
00260 }
00261 else{
00262 if(px1 > 0)phi1 = atan(py1 * 1.0 / px1);
00263 else phi1 = pai + atan(py1 * 1.0 / px1);
00264 }
00265
00266 if(atan(py2 * 1.0 / px2) > 0){
00267 if(px2 > 0)phi2 = atan(py2 * 1.0 / px2);
00268 else phi2 = -(pai - atan(py2 * 1.0 / px2));
00269 }
00270 else{
00271 if(px2 > 0)phi2 = atan(py2 * 1.0 / px2);
00272 else phi2 = pai + atan(py2 * 1.0 / px2);
00273 }
00274
00275 double dltphi = (fabs(phi1 - phi2) - pai) * 180.0 / pai;
00276
00277 double theta1 = 0;
00278 double theta2 = 0;
00279
00280 if(pz1 > 0) theta1 = atan(pxy1 * 1.0 / pz1);
00281 else theta1 = (pai + atan(pxy1 * 1.0 / pz1));
00282
00283 if(pz2 > 0) theta2 = atan(pxy2 * 1.0 / pz2);
00284 else theta2 = (pai + atan(pxy2 * 1.0 / pz2));
00285
00286 double dlttheta = (theta1 + theta2 - pai) * 180.0 / pai;
00287
00288 double xBoost = p4Gam_1.px() + p4Gam_2.px();
00289 double yBoost = p4Gam_1.py() + p4Gam_2.py();
00290 double zBoost = p4Gam_1.pz() + p4Gam_2.pz();
00291 m_xBoost = xBoost;
00292 m_yBoost = yBoost;
00293 m_zBoost = zBoost;
00294
00295 m_costheta1 = cos(theta1);
00296 m_costheta2 = cos(theta2);
00297 m_phi1 = phi1 * 180.0 / pai;
00298 m_phi2 = phi2 * 180.0 / pai;
00299 m_dlttheta = dlttheta;
00300 m_dltphi = dltphi;
00301 m_e1 = e1;
00302 m_e2 = e2;
00303 m_e = e1 + e2;
00304
00305
00306 if(fabs(m_dltphi) < m_dphi1) Ndata1++;
00307 if(fabs(m_dltphi) > m_dphi1 && fabs(m_dltphi) < m_dphi2) Ndata2++;
00308
00309
00310 m_tuple1->write();
00311
00312
00313 return StatusCode::SUCCESS;
00314 }
00315
00316
00317 StatusCode TwoGamma::finalize()
00318 {
00319 char head[100];
00320 char foot[100] = ".txt";
00321 sprintf(head,"OffLineLum_%04d",m_runNo);
00322 strcat(head,foot);
00323 ofstream fout(head,ios::out);
00324 fout<<"==============================================="<<endl;
00325 fout<<"DESIGNED For OffLine Luminosity BY 2Gam Events"<<endl<<endl;
00326 fout<<"MANY THANKS to Prof. C.Z Yuan & Z.Y Wang"<<endl<<endl;
00327 fout<<" 2009/05"<<endl;
00328
00329 fout<<"==============================================="<<endl<<endl<<endl;
00330
00331 double lum = (Ndata1 - Ndata2) * 1.0 / (m_eff * m_sec);
00332 fout<<" Table List "<<endl<<endl;
00333 fout<<"-----------------------------------------------"<<endl;
00334 fout<<"Firstly Energy Gamma "<<m_max1<<" Gev"<<endl;
00335 fout<<"Secondly Energy Gamma "<<m_max2<<" Gev"<<endl;
00336 fout<<"Solid Angle Acceptance "<<m_costheta<<endl;
00337 fout<<"QED XSection "<<m_sec<<" NB"<<endl;
00338 fout<<"Monte Carto Efficiency "<<m_eff * 100<<"%"<<endl<<endl;
00339 fout<<"runNo Luminosity "<<endl<<endl;
00340 fout<<m_runNo<<" "<<lum<<" nb^(-1)"<<endl;
00341 fout<<"-----------------------------------------------"<<endl;
00342 fout.close();
00343
00344 MsgStream log(msgSvc(), name());
00345 cout<< "in finalize()" <<endl;
00346 cout<< "all events " <<m_pass[0]<<endl;
00347 cout<< "NCharged==0 " <<m_pass[1]<<endl;
00348 cout<< "Number of Gam " <<m_pass[2]<<endl;
00349 cout<< "N events " <<m_pass[3]<<endl;
00350 cout<< "N events 0.8 " <<m_pass[4]<<endl;
00351 return StatusCode::SUCCESS;
00352 }
00353
