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

In This Package:

Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Static Private Member Functions
MuonRecTool Class Reference

#include <MuonRecTool.h>

Inheritance diagram for MuonRecTool:
Inheritance graph
[legend]
Collaboration diagram for MuonRecTool:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 MuonRecTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~MuonRecTool ()
virtual StatusCode reconstruct (const DayaBay::CalibReadout &readout, DayaBay::RecTrigger &recTrigger)
virtual StatusCode initialize ()
virtual StatusCode finalize ()
int InitialValue ()

Static Public Member Functions

static const InterfaceID & interfaceID ()

Public Attributes

int m_PmtNumber
map< unsigned int, int > m_PmtId
std::string m_cableSvcName
ICableSvcm_cableSvc
std::string m_calibDataSvcName
ICalibDataSvcm_calibDataSvc
ICoordSysSvcmCss
string m_site
IDetectorElement * m_de
DayaBay::CalibReadoutPmtCratem_calibCrate
double xini
double yini
double zini
double par [5]

Static Public Attributes

static map< unsigned int, PmtPropm_PmtPropMap
static const int maxPMTnum = 384
static const double PI = 3.1415926

Static Private Member Functions

static void fcnmu (int &npar, double *gin, double &f, double *par, int iflag)

Detailed Description

Definition at line 47 of file MuonRecTool.h.


Constructor & Destructor Documentation

MuonRecTool::MuonRecTool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 37 of file MuonRecTool.cc.

  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_calibDataSvc(0)
{
  declareInterface< IReconTool >(this) ;   
  declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                  "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("CalibDataSvcName", m_calibDataSvcName="CalibDataSvc",
                  "Name of calibration data service");

  // fixme: this should not be a property.  Instead the site should be
  // pulled out of the CalibReaout's CalibReadoutHeader.
  declareProperty("Site",m_site="DayaBay",
                  "Site name (DayaBay/Far)");
}
MuonRecTool::~MuonRecTool ( ) [virtual]

Definition at line 56 of file MuonRecTool.cc.

{}

Member Function Documentation

StatusCode MuonRecTool::reconstruct ( const DayaBay::CalibReadout readout,
DayaBay::RecTrigger recTrigger 
) [virtual]

check AD or water pool

Loop over channel data

Implements IReconTool.

Definition at line 173 of file MuonRecTool.cc.

{
  if( !(readout.detector().isWaterShield()) ) {
    return StatusCode::SUCCESS;
   }

  //reset values
  m_PmtNumber = 0;
  map<unsigned int/*PmtId*/,PmtProp>::iterator pmtit,idend=m_PmtPropMap.end();
  for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
    PmtProp& PoolPmt = pmtit->second;
    PoolPmt.Charge = 0;
    PoolPmt.FirstHit = 999999;
    PoolPmt.Used = 0;
    
    m_PmtNumber++;
  }
  
  const CalibReadoutPmtCrate* pCalibPmtCrate = 0;

  pCalibPmtCrate = dynamic_cast< const DayaBay::CalibReadoutPmtCrate* > (&readout);
  if(!pCalibPmtCrate) {
    error()<<"Failed to get CalibReadoutPmtCrate"<<endreq;
    return StatusCode::FAILURE;
  }
  
  Detector det = pCalibPmtCrate->detector();
  //ServiceMode svcMode( pCalibReadoutHdr->context(),0 );
  double   tempCharge=0, tempfirstHit;
  CalibReadoutPmtCrate::PmtChannelReadouts channels = pCalibPmtCrate->channelReadout();
  CalibReadoutPmtCrate::PmtChannelReadouts::iterator ci, ci_end = channels.end();
  for(ci = channels.begin(); ci!=ci_end; ci++) {
    unsigned int multi = ci->size();

    double charge = 0;
    for(unsigned int nr = 0; nr<multi; ++nr) {
      charge += ci->charge(nr);
    }

    unsigned int pmtid = ci->pmtSensorId().fullPackedData();
    
    m_PmtPropMap[pmtid].Charge = charge;
    double firsthit = ci->earliestTime();
    m_PmtPropMap[pmtid].FirstHit = firsthit;
    if(charge>5) m_PmtPropMap[pmtid].Used = 1;
  
   //find the max. charge PMT 
    if(charge>tempCharge)  {
       tempCharge = charge;
       xini = m_PmtPropMap[pmtid].Position.x() ;
       yini = m_PmtPropMap[pmtid].Position.y();
       zini = m_PmtPropMap[pmtid].Position.z();
       tempfirstHit = m_PmtPropMap[pmtid].FirstHit = firsthit; 
     } 
  //  cout<<"PMT id"<<pmtid<<" ; charge ="<<charge<<" ;  firsthit ="<<firsthit<<endl;
  }
  // Using the minuit for  reconstruction--------------------------------------------------
    // Initialize TMinuit via generic fitter interface with a maximum of 5 params
   int  Npar =5; 
   TVirtualFitter *minuit = TVirtualFitter::Fitter(0, Npar);
   
   minuit->SetFCN(fcnmu);
   Double_t arglist[100];
   Int_t ierflg = 0;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
   minuit->SetParameter(0, "x0",  xini,   10, xini-500, xini+500);
   minuit->SetParameter(1, "y0b", yini,   10, yini-500, yini+500);
   minuit->SetParameter(2, "z0",  zini,   10, zini-500, zini+500 );
   minuit->SetParameter(3, "theta0",  10,         0.1  , 0,          90);
   minuit->SetParameter(4, "phi0",    0  ,         0.1  ,-180,       180);
   
   arglist[0] = 1000;
   arglist[1] = 0.01;
   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg); //MINI, SIMPLEX
  // gMinuit->mnexcm("MINI", arglist ,2,ierflg); //MINI, SIMPLEX
  // gMinuit->mnexcm("SIMPLEX", arglist ,2,ierflg); //MINI, SIMPLEX

 // Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(-1,amin);
   // error matrix; mnemat()
   double emat[Npar][Npar];
   CLHEP::HepMatrix matrix(Npar, Npar);
   gMinuit->mnemat(&emat[0][0], Npar);
  for(int row=0; row<Npar; row++) {
    for(int col=0; col<Npar; col++) {
      matrix[row][col] = emat[row][col];
    }
  }

 
   double re_x0, re_y0, re_z0, re_theta0, re_phi0, re_px, re_py,re_pz;
   // get the restrcution parameters
   re_x0 = minuit->GetParameter(0);
   re_y0 = minuit->GetParameter(1);
   re_z0 = minuit->GetParameter(2);
   re_theta0 = minuit->GetParameter(3);
   re_phi0 = minuit->GetParameter(4);
   re_px = sin(re_theta0*PI/180)*cos(re_phi0*PI/180);
   re_py = sin(re_theta0*PI/180)*sin(re_phi0*PI/180);
   re_pz = (-1)*cos(re_theta0*PI/180);

  // cout<<"icstat ="<<icstat<<endl;
  // cout<<"px1="<<re_px <<",py1="<<re_py<<",pz1="<<re_pz<<endl;
   CLHEP::HepLorentzVector pos(re_x0,re_y0,re_z0,0);
   CLHEP::HepLorentzVector dir(re_px,re_py,re_pz,0);
   int fitstat=0;
   if(icstat<3)  
      fitstat =0;
   else if(icstat==3)
      fitstat =1;
   // set the position and direction to the recTrigger 
   recTrigger.setPosition(pos);
   recTrigger.setPositionStatus(fitstat);
   recTrigger.setDirection(dir);  
   recTrigger.setDirectionStatus(fitstat);  
   recTrigger.setErrorMatrix(matrix);
  
   return StatusCode::SUCCESS;
}
StatusCode MuonRecTool::initialize ( ) [virtual]

Definition at line 58 of file MuonRecTool.cc.

{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  m_calibDataSvc = svc<ICalibDataSvc>(m_calibDataSvcName,true);

  //------------------
  StatusCode sc;
  sc = service("CoordSysSvc",mCss);
  if(sc.isFailure()) {
    error()<<"Can't get CoordSysSvc"<<endreq;
    return sc;
  }

  // Make use of svc<>
  IMessageSvc* msg = svc<IMessageSvc>("MessageSvc");
  MsgStream log(msg, name());
  log << MSG::DEBUG << "intializing ......" << endreq;
  
  // Get PmtGeoSvc
  IPmtGeomInfoSvc* pmtGeoSvc;
  sc = service("PmtGeomInfoSvc",pmtGeoSvc);
  if(sc.isFailure()) {
    error()<<"Can't get PmtGeomInfoSvc"<<endreq;
    return sc;
  }
 
  
  //choose Site
  //Site :: Site_t sitestr = Site :: FromString("DayaBay");
  Site :: Site_t sitestr = Site :: FromString(m_site.c_str());
  //choose DetectorId
  DetectorId::DetectorId_t idetid = DetectorId :: FromString("IWS");
  DetectorId::DetectorId_t odetid = DetectorId :: FromString("OWS");
  //for the inward_facing or outward_facing
  bool inward_facing = true;
  bool outward_facing = false;
  CLHEP::Hep3Vector global;
  Gaudi::XYZPoint pos;

  string searchDet="/dd/Structure/Pool/";
  if(m_site=="Far")     searchDet.append("far-ows");
  if(m_site=="DayaBay") searchDet.append("db-ows");
  m_de = getDet<IDetectorElement>(searchDet);

  int PMTcount = 0;
  //IWS
  for(int iwall_number=1;iwall_number<10;++iwall_number){
    for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,idetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}

      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
      m_PmtId[pmtid] = PMTcount;
      PMTcount++;
      //cout<<"IWS-------wall ="<<iwall_number<<"; spot= "<<iwall_spot<<"---------------------------------"<<endl;
     // cout<<"PMT x ="<<m_PmtPropMap[pmtid].Position.x()<<";  y ="<<m_PmtPropMap[pmtid].Position.y()<<"; z ="<<m_PmtPropMap[pmtid].Position.z()<<endl;  
     }
  }
  //Outer Pool OWS facing outward
  for(int iwall_number=1;iwall_number<9;++iwall_number){
    for(int iwall_spot=1;iwall_spot<17;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,outward_facing,sitestr,odetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}
      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
      PMTcount++;
      //   cout<<"OWS faceout-------wall ="<<iwall_number<<"; spot= "<<iwall_spot<<"---------------------------------"<<endl;
     // cout<<"PMT x ="<<m_PmtPropMap[pmtid].Position.x()<<";  y ="<<m_PmtPropMap[pmtid].Position.y()<<"; z ="<<m_PmtPropMap[pmtid].Position.z()<<endl;  
    }      
  }
  //Outer Pool OWS facing inward
  for(int iwall_number=1;iwall_number<10;++iwall_number){
    for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,odetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}
//      pos = pmtGeoSvc->get(pmtid)->localPosition();
//      cout<< "PMT[" << count << "]=new TVector3(" << pos.x() <<"," 
//        <<pos.y()<<","<<pos.z()<<");"<<endl;

      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);
      m_PmtId[pmtid] = PMTcount;
      PMTcount++;
      //cout<<"OWS facein-------wall ="<<iwall_number<<"; spot= "<<iwall_spot<<"---------------------------------"<<endl;
     //  cout<<"PMT x ="<<m_PmtPropMap[pmtid].Position.x()<<";  y ="<<m_PmtPropMap[pmtid].Position.y()<<"; z ="<<m_PmtPropMap[pmtid].Position.z()<<endl;  
    }
  }

  //-----print out pmt positions 
//  map<unsigned int/*PmtId*/,PmtProp>::iterator pmtit,idend=m_PmtPropMap.end();
//  for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
//    PmtProp& PoolPmt = pmtit->second;
//    cout << pmtit->first << " ("
//       << PoolPmt.Position.X() << "," 
//       << PoolPmt.Position.Y() << "," 
//       << PoolPmt.Position.Z() << ")" 
//       << endl;
//  }

  return StatusCode::SUCCESS;
}
StatusCode MuonRecTool::finalize ( ) [virtual]

Definition at line 168 of file MuonRecTool.cc.

{
  return StatusCode::SUCCESS;
}
int MuonRecTool::InitialValue ( )
void MuonRecTool::fcnmu ( int &  npar,
double *  gin,
double &  f,
double *  par,
int  iflag 
) [static, private]

Definition at line 299 of file MuonRecTool.cc.

                                                                                {
    
     double  nwater =1.35; // using the water index n =1.35
     double  t_sigma = 1.2; //assumed the jetter is 1.2 ns
     double  clight  =299.792458; //unit of  mm/ns, the light speed
     double  t0 = 1561; 
     double  Vcmuon= clight; // unit mm/ns, muon speed in water
     double  Vcwater = clight/nwater;
     double  thetac = TMath::ACos(1/nwater); // the Cerenkov degree, in unit of rad
   
     double  dist[maxPMTnum],trackvalue[maxPMTnum],expect_time[maxPMTnum];
     double  pmtUsed_posx[maxPMTnum],pmtUsed_posy[maxPMTnum],pmtUsed_posz[maxPMTnum],pmtUsed_time[maxPMTnum];
     int useNum =0,ht =0;
     double  t_Trans_mean =0, t_Trans=0,chisq0 = 0, delta0;
     double px,py,pz;
     map<unsigned int,PmtProp>::iterator inum, hend=m_PmtPropMap.end();
    
     for(inum=m_PmtPropMap.begin();inum!=hend;inum++) {
        PmtProp& PoolPmt = inum->second;
       if(PoolPmt.Used==1){
          useNum++;
          pmtUsed_posx[ht] = PoolPmt.Position.x();
          pmtUsed_posy[ht] = PoolPmt.Position.y();
          pmtUsed_posz[ht] = PoolPmt.Position.z();
          pmtUsed_time[ht]  = PoolPmt.FirstHit+ t0;
          px= sin(par[3]*PI/180)*cos(par[4]*PI/180);
          py= sin(par[3]*PI/180)*sin(par[4]*PI/180);
          pz= (-1)*cos(par[3]*PI/180);
          dist[ht]=TMath::Sqrt((py*(par[2]-pmtUsed_posz[ht])-pz*(par[1]-pmtUsed_posy[ht]))*(py*(par[2]-pmtUsed_posz[ht])-pz*(par[1]-pmtUsed_posy[ht]))+(pz*(par[0]-pmtUsed_posx[ht])-px*(par[2]-pmtUsed_posz[ht]))*(pz*(par[0]-pmtUsed_posx[ht])-px*(par[2]-pmtUsed_posz[ht]))+(px*(par[1]-pmtUsed_posy[ht])-py*(par[0]-pmtUsed_posx[ht]))*(px*(par[1]-pmtUsed_posy[ht])-py*(par[0]-pmtUsed_posx[ht])))/TMath::Sqrt(px*px+py*py+pz*pz);
          trackvalue[ht]=px*(pmtUsed_posx[ht]-par[0])+py*(pmtUsed_posy[ht]-par[1])+pz*(pmtUsed_posz[ht]-par[2]);
        
          //get the PMT expect time ............................................
          expect_time[ht]=(trackvalue[ht]-dist[ht]/tan(thetac))/Vcmuon+dist[ht]/(sin(thetac))/Vcwater;
          t_Trans =  expect_time[ht]-pmtUsed_time[ht];
          t_Trans_mean =t_Trans_mean + t_Trans;
          ht++;
       }//if PMT used
     }// for ht 
    t_Trans_mean =  t_Trans_mean/useNum ;
   // cout<<" t_Trans_mean="<< t_Trans_mean<<endl;
  // cout<<"useNum ="<<useNum<<endl;
    useNum =0;
    // get the chi^2...
   for(inum=m_PmtPropMap.begin();inum!=hend;inum++){
        PmtProp& PoolPmt = inum->second;
        if(PoolPmt.Used==1) {
            useNum++;
            delta0 =(pmtUsed_time[ht]+t_Trans_mean-expect_time[ht])/t_sigma;
            chisq0 += delta0*delta0;
          }
       }
         f = chisq0;
  }

Member Data Documentation

map< unsigned int, PmtProp > MuonRecTool::m_PmtPropMap [static]

Definition at line 63 of file MuonRecTool.h.

Definition at line 65 of file MuonRecTool.h.

map<unsigned int,int> MuonRecTool::m_PmtId

Definition at line 66 of file MuonRecTool.h.

Definition at line 70 of file MuonRecTool.h.

Definition at line 73 of file MuonRecTool.h.

Definition at line 76 of file MuonRecTool.h.

Definition at line 79 of file MuonRecTool.h.

Definition at line 83 of file MuonRecTool.h.

Definition at line 85 of file MuonRecTool.h.

IDetectorElement* MuonRecTool::m_de

Definition at line 86 of file MuonRecTool.h.

Definition at line 87 of file MuonRecTool.h.

Definition at line 88 of file MuonRecTool.h.

Definition at line 88 of file MuonRecTool.h.

Definition at line 88 of file MuonRecTool.h.

double MuonRecTool::par[5]

Definition at line 89 of file MuonRecTool.h.

const int MuonRecTool::maxPMTnum = 384 [static]

Definition at line 90 of file MuonRecTool.h.

const double MuonRecTool::PI = 3.1415926 [static]

Definition at line 91 of file MuonRecTool.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:12:29 for MuonRec by doxygen 1.7.4