SD0Tag Class Reference

#include <SD0Tag.h>

Inheritance diagram for SD0Tag:

Public Member Functions
	SD0Tag (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()

Detailed Description

Definition at line 13 of file SD0Tag.h.

Constructor & Destructor Documentation

SD0Tag()

SD0Tag::SD0Tag	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 63 of file SD0Tag.cxx.

                                                              :
  Algorithm(name, pSvcLocator) {
    //Declare the properties
    declareProperty("PID_FLAG",            PID_flag = 1);
    declareProperty("MC_sample",           m_MC_sample = 1);
 
    declareProperty("m_Seperate_Charge",   Seperate_Charge = 1);
    declareProperty("m_Charge_default",    Charge_default = 1);
             
 
  }

Member Function Documentation

execute()

StatusCode SD0Tag::execute

(

)

Definition at line 127 of file SD0Tag.cxx.

                           {
 
  MsgStream log(msgSvc(), name());
 
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  int runNo=eventHeader->runNumber();
  int event=eventHeader->eventNumber();
 
  SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
  //cout<<"//////////////////////////////////////"<<endl;
  nD0++;
 
  //
  // The default energy is 3.773 GeV for psi(3770) data. 
  // Alayizer can add calibrated energy here.
  double Ebeam = 3.773/2.0;
 
  if(runNo >=11414 && runNo <= 23500) {
    for(int i = 0; i < 3467; i++) {
      if(runNo == p_run[i]) {Ebeam = p_Ecm[i]/2.0;}
    }
  }
 
  if(runNo < 0)
  {
    int irun = abs(runNo);
//    Ecmset*  ECMSET = Ecmset::instance();
//    ECMSET->EcmSet(irun);
//    Ebeam =  ECMSET->ECM()/2.0;
  }
 
  if(nD0%1000==0) cout<<"SD0TagAlg-13-11-15pp = "<<nD0<<"  "<<runNo<<"   "<<event<<"   "<<Ebeam*2<<endl;
 
  Hep3Vector xorigin(0,0,0);
  IVertexDbSvc*  vtxsvc;
  Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
  if(vtxsvc->isVertexValid()){
    double* dbv = vtxsvc->PrimaryVertex();
    double*  vv = vtxsvc->SigmaPrimaryVertex();
    xorigin.setX(dbv[0]);
    xorigin.setY(dbv[1]);
    xorigin.setZ(dbv[2]);
  }
 
  /////////////////////////////////////
  Vint iCharge_good; iCharge_good.clear();
  for(int i = 0; i <  evtRecEvent->totalCharged(); i++){
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
    if(!(*itTrk)->isMdcTrackValid()) continue;
    RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
    //%%%%%%%%%%%%%%%%%%%%%%%%
    HepVector a = mdcTrk->helix();
    HepSymMatrix Ea = mdcTrk->err();
    HepPoint3D point0(0.,0.,0.);   // the initial point for MDC recosntruction
    HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
    VFHelix helixip(point0,a,Ea);
    helixip.pivot(IP);
    HepVector vecipa = helixip.a();
    double  Rvxy0=fabs(vecipa[0]);  //the nearest distance to IP in xy plane
    double  Rvz0=vecipa[3];         //the nearest distance to IP in z direction
    double  costheta = cos(mdcTrk->theta());
    //%%%%%%%%%%%%%%%%%%%%%%%%
    if(fabs(Rvxy0) >= 1.0) continue;
    if(fabs(Rvz0)  >= 10.0) continue;
    if(fabs(costheta) >= 0.930) continue;
    iCharge_good.push_back(i);
  }
 
  ///////////////////////////////////////////////////////
  Vint iGam; iGam.clear();
  for(int i =  evtRecEvent->totalCharged(); i<  evtRecEvent->totalTracks(); i++) {
    EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
    if(!(*itTrk)->isEmcShowerValid()) continue;
    RecEmcShower *emcTrk = (*itTrk)->emcShower();
    //
    //   We here remove the hot channels of EMC 
    //
    if(abs(runNo)>=11615&&abs(runNo)<=11617&&emcTrk->cellId()==805438481)  continue;
    if(abs(runNo)>=11615&&abs(runNo)<=11655&&emcTrk->cellId()==805376008)  continue;
    if(abs(runNo)>=13629&&abs(runNo)<=13669&&emcTrk->cellId()==805374754)  continue;
    if(abs(runNo)>=21190&&abs(runNo)<=21231&&emcTrk->cellId()==805375262)  continue;
 
    int emctdc = emcTrk->time();
    if(emctdc>14||emctdc<0) continue;
 
    double eraw = emcTrk->energy();
    double theta = emcTrk->theta();
    int Module = 0;
    if(fabs(cos(theta)) < 0.80 && eraw > 0.025){   Module = 1;   }
    if(fabs(cos(theta)) > 0.86 && fabs(cos(theta)) < 0.92 && eraw > 0.05) {   Module = 2;      }
    if(Module == 0)  continue; 
    iGam.push_back((*itTrk)->trackId());
  }  
  /////////////////////////////////////////////////////////////////
  DTagTool trk;
  bool cosmic_ok = trk.cosmicandleptonVeto();
  /////////////////////////////////////////////////////////////////
 
  int   ntk;
  double tagmass,ebeam,tagmode,ksmass,dlte;
 
  Vint tagtrk;   tagtrk.clear();
  Vint tagGam;   tagGam.clear();
 
  HepLorentzVector  tagp;
 
  double mass_bc_temp, mass_kf_temp, kf_chi2_temp, mks_temp, mpi0_temp, ptag_temp; 
  int Charge_candidate_D = 0;
  double  EGam_max_0 = 0;
 
  for(int i1 = 0; i1 < 2; i1++) {
    if(Seperate_Charge == 2 ) {  Charge_candidate_D = Charge_default; i1 = 2;}
    if(Seperate_Charge == 1 ) {  Charge_candidate_D = pow(-1,i1); }
    if(Seperate_Charge == 0 ) {  Charge_candidate_D = 0; i1 = 2; }
 
    for(int i = 0; i < 15; i++) {
      EGam_max_0 = 0;   
      int  mdslct=pow(2.,i);
      Sing sing;
      sing.Mdset(event,evtRecTrkCol,iCharge_good,iGam,mdslct,Ebeam, PID_flag,Charge_candidate_D);
      bool oktag=sing.Getoktg();
      if(oktag) {
        //
        //  Here analysizer should pick up variables from anti-D0 tags
        //  to do physics analysis
        //
        tagtrk=sing.Gettagtrk1();
        tagmode = sing.Gettagmd();
        //==========================================================================
        if((abs(tagmode) == 11 || abs(tagmode) == 15 || abs(tagmode) == 19)) {
          //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          for(int i =  evtRecEvent->totalCharged(); i <  evtRecEvent->totalTracks(); i++) {
            EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
            int itrk = (*itTrk)->trackId();
            if(!(*itTrk)->isEmcShowerValid()) continue;
            RecEmcShower *emcTrk = (*itTrk)->emcShower();
 
            int emctdc = emcTrk->time();
            if(emctdc>14||emctdc<0) continue;
 
            Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z());
            double dang = 200.;
            for(int j = 0; j <  tagtrk.size(); j++) {
              EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + tagtrk[j];
              if(!(*jtTrk)->isExtTrackValid()) continue;
              RecExtTrack *extTrk = (*jtTrk)->extTrack();
              if(extTrk->emcVolumeNumber() == -1) continue;
              Hep3Vector extpos = extTrk->emcPosition();
              double angd = extpos.angle(emcpos);
              if(angd < dang) dang = angd;
            } 
            dang = dang * 180 / (CLHEP::pi);
            if(fabs(dang)<20.) continue;
 
            double eraw = emcTrk->energy();
            if(eraw > EGam_max_0) EGam_max_0 = eraw;
          }          
        }     
        //==========================================================================
        m_cosmic_ok  = cosmic_ok;
        m_EGam_max_0 = EGam_max_0;
        m_nGood = iCharge_good.size();
        m_nGam = iGam.size(); 
        m_runNo = runNo;
        m_event = event;
 
        m_tagmode = sing.Gettagmd();
        m_mass_bc = sing.Getmass_bc();
        m_delE_tag = sing.GetdelE_tag();
 
        m_tuple1->write();
 
        //      Here one can do double tag analysis based on the variables
        //      from single anti-D0 tag analysis.
        //      To do so, analyzer should make his/her codes
        //
      }
 
    }
  }  // The end of loopping over mdslct
 
  return StatusCode::SUCCESS;
 
}

finalize()

StatusCode SD0Tag::finalize

(

)

Definition at line 318 of file SD0Tag.cxx.

                            {
  MsgStream log(msgSvc(), name());
  cout<<"nD0         =    "<<nD0<<endl;
  cout<<"n1         =    "<<n1<<endl;
  cout<<"n2         =    "<<n2<<endl;
  cout<<"ND0     ====   "<<ND0<<endl;
  cout<<"NDp     ====   "<<NDp<<endl;
  log << MSG::INFO << "in finalize()" << endmsg;
  return StatusCode::SUCCESS;
}

initialize()

StatusCode SD0Tag::initialize

(

)

Definition at line 77 of file SD0Tag.cxx.

                             {
 
 
  MsgStream log(msgSvc(), name());
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  StatusCode status;
 
  ifstream readrunEcm;
//
// Read in the energy calibration constant for run-by-run data
  readrunEcm.open("../share/psipp_ecms_calrunbyrun_runno_3648runs.dat");
  cout<<"readrunEcm = "<<readrunEcm.is_open()<<endl;
  for(int i = 0; i < 3467; i++) {
    readrunEcm>>p_run[i]; readrunEcm>>p_Ecm[i];
  }
 
 
  NTuplePtr nt1(ntupleSvc(), "FILE1/tag");
  if ( nt1 ) m_tuple1 = nt1;
  else {
    m_tuple1 = ntupleSvc()->book ("FILE1/tag", CLID_ColumnWiseTuple, "tag N-Tuple example");
    if ( m_tuple1 )  {
      status = m_tuple1->addItem ("tagmode",          m_tagmode);
      status = m_tuple1->addItem ("mass_bc",          m_mass_bc);
      status = m_tuple1->addItem ("delE_tag",         m_delE_tag);
 
      status = m_tuple1->addItem ("cosmic_ok",        m_cosmic_ok);    
      status = m_tuple1->addItem ("EGam_max_0",       m_EGam_max_0);
      status = m_tuple1->addItem ("nGood",            m_nGood);      
      status = m_tuple1->addItem ("nGam",             m_nGam);      
      status = m_tuple1->addItem ("runNo",            m_runNo);      
      status = m_tuple1->addItem ("event",            m_event);      
    }
    else  {
      log << MSG::ERROR << "    Cannot book N-tuple:" << long(m_tuple1) <<endmsg;
      return StatusCode::FAILURE;
    }
  }  
 
  log << MSG::INFO << "successfully return from initialize()" <<endmsg;
  return StatusCode::SUCCESS;
}

---

The documentation for this class was generated from the following files:

• SD0Tag.h
• SD0Tag.cxx