/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

MuonSimpleAlg.cc
Go to the documentation of this file.
00001 #include "MuonSimpleAlg.h"
00002 
00003 #include "Event/Readout.h"
00004 #include "Event/ReadoutPmtCrate.h"
00005 #include "Event/ReadoutPmtChannel.h"
00006 #include "Event/CalibReadout.h"
00007 #include "Event/CalibReadoutPmtCrate.h"
00008 #include "Event/CalibReadoutPmtChannel.h"
00009 #include "DataSvc/ICableSvc.h"
00010 #include "Event/RawEventHeader.h"
00011 #include "DetDesc/IGeometryInfo.h"
00012 #include "DetDesc/DetectorElement.h"
00013 #include "DetHelpers/IPmtGeomInfoSvc.h"
00014 #include "CLHEP/Vector/LorentzVector.h"
00015 #include <math.h>
00016 #include <sstream>
00017 
00018 #include "Event/RecHeader.h"
00019 #include "Event/RecTrigger.h"
00020 #include "Event/CalibReadoutHeader.h"
00021 #include "Conventions/Reconstruction.h"
00022 #include "Conventions/Detectors.h"
00023 #include <vector>
00024 
00025 using namespace std;
00026 using namespace DayaBay;
00027 
00028 MuonSimpleAlg::MuonSimpleAlg(const std::string& name, ISvcLocator* pSvcLocator)
00029   : DybAlgorithm<DayaBay::UserDataHeader>(name,pSvcLocator)
00030 {
00031   declareProperty("Site",m_site="DayaBay","Site name (DayaBay/Far)");
00032   declareProperty("m_timeCutLow", m_timeCutLow=-1700., "Lower limit of the TDC cut to get rid of noise");
00033   declareProperty("m_timeCutHigh", m_timeCutHigh=-1400., "Higher limit of the TDC cut to get rid of noise");
00034 }
00035 
00036 MuonSimpleAlg::~MuonSimpleAlg(){}
00037 
00038 StatusCode MuonSimpleAlg::initialize()
00039 {
00040   StatusCode sc = this->GaudiAlgorithm::initialize();
00041   if(sc.isFailure()) return sc;
00042 
00044   sc = service("StaticCableSvc",m_cableSvc);
00045   if( sc.isFailure() ) {
00046     error() << "Can't get StaticCableSvc" << endreq;
00047     return sc;
00048   }
00049 
00050   sc = service("StaticCalibDataSvc",m_calibDataSvc);
00051   if( sc.isFailure() ) {
00052     error() << "Can't get StaticCalibDataSvc" << endreq;
00053     return sc;
00054   }
00055 
00056   sc = service("CoordSysSvc",mCss);
00057   if(sc.isFailure()) {
00058     error()<<"Can't get CoordSysSvc"<<endreq;
00059     return sc;
00060   }
00061 
00062   // Make use of svc<>
00063   IMessageSvc* msg = svc<IMessageSvc>("MessageSvc");
00064   MsgStream log(msg, name());
00065   log << MSG::DEBUG << "intializing ......" << endreq;
00066   
00067   // Get PmtGeoSvc
00068   IPmtGeomInfoSvc* pmtGeoSvc;
00069   sc = service("PmtGeomInfoSvc",pmtGeoSvc);
00070   if(sc.isFailure()) {
00071     error()<<"Can't get PmtGeomInfoSvc"<<endreq;
00072     return sc;
00073   }
00074   
00075   //choose Site
00076   //Site :: Site_t sitestr = Site :: FromString("DayaBay");
00077   Site :: Site_t sitestr = Site :: FromString(m_site.c_str());
00078   //choose DetectorId
00079   DetectorId::DetectorId_t idetid = DetectorId :: FromString("IWS");
00080   DetectorId::DetectorId_t odetid = DetectorId :: FromString("OWS");
00081   //for the inward_facing or outward_facing
00082   bool inward_facing = true;
00083   bool outward_facing = false;
00084   CLHEP::Hep3Vector global;
00085   Gaudi::XYZPoint pos;
00086 
00087   string searchDet="/dd/Structure/Pool/";
00088   if(m_site=="Far")     searchDet.append("far-ows");
00089   if(m_site=="DayaBay") searchDet.append("db-ows");
00090   m_de = getDet<IDetectorElement>(searchDet);
00091 
00092   int PMTcount = 0;
00093   //IWS
00094   for(int iwall_number=1;iwall_number<10;++iwall_number){
00095     for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
00096       PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,idetid);
00097       unsigned int pmtid = poolPmtId.fullPackedData();
00098       bool p = pmtGeoSvc->get(pmtid); if(!p){break;}
00099 
00100       global = pmtGeoSvc->get(pmtid)->globalPosition();
00101       Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
00102       m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
00103 
00104       m_PmtId[pmtid] = PMTcount;
00105       PMTcount++;
00106     }
00107   }
00108   //Outer Pool OWS facing outward
00109   for(int iwall_number=1;iwall_number<9;++iwall_number){
00110     for(int iwall_spot=1;iwall_spot<17;++iwall_spot){
00111       PoolPmtSensor poolPmtId(iwall_number,iwall_spot,outward_facing,sitestr,odetid);
00112       unsigned int pmtid = poolPmtId.fullPackedData();
00113       bool p = pmtGeoSvc->get(pmtid); if(!p){break;}
00114 
00115       global = pmtGeoSvc->get(pmtid)->globalPosition();
00116       Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
00117       m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
00118 
00119       m_PmtId[pmtid] = PMTcount;
00120       PMTcount++;
00121     }      
00122   }
00123   //Outer Pool OWS facing inward
00124   for(int iwall_number=1;iwall_number<10;++iwall_number){
00125     for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
00126       PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,odetid);
00127       unsigned int pmtid = poolPmtId.fullPackedData();
00128       bool p = pmtGeoSvc->get(pmtid); if(!p){break;}
00129 
00130       global = pmtGeoSvc->get(pmtid)->globalPosition();
00131       Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
00132       m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
00133 
00134       m_PmtId[pmtid] = PMTcount;
00135       PMTcount++;
00136     }
00137   }
00138 
00139   return StatusCode::SUCCESS;
00140 }
00141 
00142 StatusCode MuonSimpleAlg::finalize()
00143 {
00144   return this->GaudiAlgorithm::finalize();
00145 }
00146 
00147 StatusCode MuonSimpleAlg::execute()
00148 {
00149   DayaBay::HeaderObject *muonTag = 0;
00150 
00151   if (exist<DayaBay::HeaderObject>(evtSvc(), "/Event/Tag/Muon/MuonAll")) {
00152     muonTag = get<DayaBay::HeaderObject>("/Event/Tag/Muon/MuonAll");
00153     //    calibHdrs = muonTag->findHeaders("/Event/CalibReadout/CalibReadoutHeader");
00154     std::vector<const DayaBay::IHeader*> calibHdrs = muonTag->findHeaders(CalibReadoutHeader::classID());
00155 
00156     //info() << "Calib header size:" << calibHdrs.size() << endreq;
00157 
00158     Gaudi::XYZPoint p1(0,0,0), p2(0,0,0);
00159     Gaudi::XYZVector vect;
00160     int ntrigger = 0;
00161     vector<const DayaBay::IHeader*>::iterator it;
00162     for(it=calibHdrs.begin(); it!=calibHdrs.end(); it++){
00163       const CalibReadoutHeader* calibHeader = dynamic_cast<const CalibReadoutHeader*>(*it);
00164       
00165       const CalibReadout* calibReadout = calibHeader->calibReadout();
00166       if(!calibReadout) {
00167         info() << "No CalibReadout retrieved this cycle." << endreq;
00168         continue;
00169       } 
00170       if(!calibReadout->detector().isWaterShield()) continue;
00171 
00172       ntrigger++;
00173       //reset values
00174       map<unsigned int/*PmtId*/,PmtProp>::iterator pmtit,idend=m_PmtPropMap.end();
00175       for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
00176         PmtProp& PoolPmt = pmtit->second;
00177         PoolPmt.Charge = 0;
00178         PoolPmt.FirstHit = 999999;
00179       }
00180 
00181       const CalibReadoutPmtCrate* pCalibPmtCrate = 0;    
00182       pCalibPmtCrate = dynamic_cast< const DayaBay::CalibReadoutPmtCrate* >(calibReadout);
00183       if(!pCalibPmtCrate) {
00184         error()<<"Failed to get CalibReadoutPmtCrate"<<endreq;
00185         return StatusCode::FAILURE;
00186       }
00187       
00188       Detector det = pCalibPmtCrate->detector();
00189       //ServiceMode svcMode( pCalibReadoutHdr->context(),0 );
00190     
00192       CalibReadoutPmtCrate::PmtChannelReadouts channels = pCalibPmtCrate->channelReadout();
00193       CalibReadoutPmtCrate::PmtChannelReadouts::iterator ci, ci_end = channels.end();
00194       for(ci = channels.begin(); ci!=ci_end; ci++) {
00195         unsigned int multi = ci->size();
00196         
00197         double charge = 0;
00198         for(unsigned int i=0; i<multi; i++){
00199           if(ci->time(i)>m_timeCutLow && ci->time(i)<m_timeCutHigh){
00200             if(charge<ci->charge(i)) charge = ci->charge(i);
00201             //if(hittime>calibChannel->time(i)) hittime = calibChannel->time(i);
00202           }
00203         } 
00204 
00205         //int sensor = ci->pmtSensorId().sensorId();
00206         //if(sensor==0)  continue;//skip fake PMT hits
00207         //PoolPmtSensor PmtId = PoolPmtSensor(sensor);
00208         //unsigned int pmtid = PmtId.fullPackedData();
00209         unsigned int pmtid = ci->pmtSensorId().fullPackedData();
00210         
00211         m_PmtPropMap[pmtid].Charge = charge;
00212         m_PmtPropMap[pmtid].FirstHit = ci->earliestTime();
00213       }
00214       
00215       //find maximum PMT
00216       double max = 0;
00217       double maxid = 0;
00218       for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
00219         PmtProp& PoolPmt = pmtit->second;
00220         if(max<PoolPmt.Charge){
00221           max = PoolPmt.Charge;
00222           maxid = pmtit->first;
00223         }
00224       }
00225       
00226       //if(ntrigger==1) p1 = m_PmtPropMap[maxid].Position;
00227       //if(ntrigger==2) p2 = m_PmtPropMap[maxid].Position;
00228       if(calibReadout->detector().detectorId()==5) p1 = m_PmtPropMap[maxid].Position;
00229       if(calibReadout->detector().detectorId()==6) p2 = m_PmtPropMap[maxid].Position;
00230     }
00231     
00232     int PositionStatus = 0;
00233     float PositionX = 0;
00234     float PositionY = 0;
00235     float PositionZ = 0;
00236     int DirectionStatus = 0;
00237     float DirectionX = 0;
00238     float DirectionY = 0;
00239     float DirectionZ = 0;
00240     
00241     //when both IWS and OWS triggers exist
00242     if(p1.Mag2()>0 && p2.Mag2()>0){
00243       //info() << "simple muon reconstruction using IWS and OWS" << endreq;
00244 
00245       vect = p1-p2;
00246       PositionStatus = 1;
00247       PositionX = p1.X();
00248       PositionY = p1.Y();
00249       PositionZ = p1.Z();
00250       DirectionStatus = 1;    
00251       DirectionX = vect.X();
00252       DirectionY = vect.Y();
00253       DirectionZ = vect.Z();
00254       
00255     }
00256     
00257     // Create output user data header
00258     UserDataHeader* muonSimple = MakeHeaderObject();
00259     muonSimple->set("PositionStatus",PositionStatus);
00260     muonSimple->set("PositionX",PositionX);
00261     muonSimple->set("PositionY",PositionY);
00262     muonSimple->set("PositionZ",PositionZ);
00263     muonSimple->set("DirectionStatus",DirectionStatus);
00264     muonSimple->set("DirectionX",DirectionX);
00265     muonSimple->set("DirectionY",DirectionY);
00266     muonSimple->set("DirectionZ",DirectionZ);
00267     
00268   }
00269   return StatusCode::SUCCESS;
00270 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:11:18 for MuonSimple by doxygen 1.7.4