/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 | Private Attributes
ChargeTemplatePosTool Class Reference

#include <ChargeTemplatePosTool.h>

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

List of all members.

Public Member Functions

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

Static Public Member Functions

static const InterfaceID & interfaceID ()

Private Attributes

std::string m_templateFileName
std::string m_cableSvcName
ICableSvcm_cableSvc
std::string m_pmtGeomSvcName
IPmtGeomInfoSvcm_pmtGeomSvc
std::string m_ChannelQualitySvcName
IChannelQualitySvcm_chanqualSvc
bool m_useChannelQualitySvc
int nocqCount
bool Rsq
bool Debug
bool Debug2
bool Interpol
bool FineFit
bool FullScan
bool m_calculateChi2Only
bool m_interpolateChi2
bool m_minuitMinimization
bool m_debugOutput

Detailed Description

Definition at line 27 of file ChargeTemplatePosTool.h.


Constructor & Destructor Documentation

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

Definition at line 113 of file ChargeTemplatePosTool.cc.

  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_pmtGeomSvc(0)
  , m_chanqualSvc(0)
{
    declareInterface< IReconTool >(this) ;   
    declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
   
    declareProperty("PmtGeomSvcName", m_pmtGeomSvcName="PmtGeomInfoSvc",
                    "Name of PMT geometry data service");

    std::string Path = getenv ("DATASVCROOT");
    //    declareProperty("TemplateFileName", m_templateFileName = "/u/ynakajim/mywork/vertex_recon_template/template_outputs/templates/ChargeTemplate_T1_0001-0298_1-9.root",
    declareProperty("TemplateFileName", m_templateFileName = Path+"/share/ChargeTemplate.root",
                    "File name of the charge distribution template");

    declareProperty("CalculateChi2Only", m_calculateChi2Only=false,
           "Calculate chi2 vaule only for the vertex position given by other tool");

    declareProperty("InterpolateChi2", m_interpolateChi2=true,
           "Make simple Chi2 interpolation to find Minimum point");

    declareProperty("MinuitMinimization", m_minuitMinimization=false,
           "use interpolated templates and Minuit for minimization");

    declareProperty("DebugOutput", m_debugOutput=false,
                    "Output debug informaiton");
    
    declareProperty("ChannelQualitySvcName", m_ChannelQualitySvcName="DybChannelQualitySvc",
                    "Name of channel quality service");

    declareProperty("UseChannelQualitySvc", m_useChannelQualitySvc=true,
                    "Use ChannelQualityService to check if the PMTs are alive. Turn it off if it's too slow");
    

    
}
ChargeTemplatePosTool::~ChargeTemplatePosTool ( ) [virtual]

Definition at line 155 of file ChargeTemplatePosTool.cc.

{}

Member Function Documentation

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

Implements IReconTool.

Definition at line 274 of file ChargeTemplatePosTool.cc.

{
  if( !readout.detector().isAD() ){
    debug() << "Not an AD readout; ignoring detector " 
            << readout.detector().detName() << endreq;
    recTrigger.setPositionStatus( ReconStatus::kNotProcessed );
    return StatusCode::SUCCESS;
  }

  const CalibReadoutPmtCrate* pmtReadout 
    = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
  if(!pmtReadout){
    error() << "Incorrect type of readout crate for detector "
            << readout.detector().detName() << endreq;
    recTrigger.setPositionStatus( ReconStatus::kBadReadout );
    return StatusCode::FAILURE;
  }

  // Context for this data
  int task = 0;
  ServiceMode svcMode(readout.header()->context(), task);


  Int_t nActivePMTs = 0;

  if (m_useChannelQualitySvc){
    const IChannelQuality * cq = m_chanqualSvc->channelQuality(svcMode);
    if (!cq){
      nocqCount ++;
      if (nocqCount <= NOCQ_ERRORS){
        warning() << "\tGot null channel quality object for " << readout.header()->context().AsString()
                  << ". Assuming all channels are good." << endreq;
      }
      if (nocqCount == NOCQ_ERRORS){
        warning() << "\tFurther warning for channel quality service will be suppressed" << endreq;
      }
      
      nActivePMTs = 192;
      for (Int_t iring = 0; iring < 8; iring++){
        for (Int_t icolumn = 0; icolumn < 24; icolumn++){
          pmt_active_flag[iring][icolumn] = 1;
        }
      }
      
    }else{
      
      // loop over channels and check if PMTs are active
      IChannelQuality::ChannelSet_t chans = cq->channels();
      for (int i = 0; i < chans.size(); i++){
        DayaBay::FeeChannelId chid = chans[i];
        AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chid, svcMode);
        if (pmtId.is8inch()){
          if (cq->good(chid)){
            pmt_active_flag[pmtId.ring()-1][pmtId.column()-1] = 1;
            nActivePMTs ++;
          }else{
            pmt_active_flag[pmtId.ring()-1][pmtId.column()-1] = 0;
            //    std::cout << "\tGod bad one: ring " << pmtId.ring() << " column " << pmtId.column() << std::endl;
          
          }
        }
      }
    }
    //    std::cout << "Number of active PMTs = " << nActivePMTs << std::endl;
  }else{ // all PMTs are assumed to be active
    nActivePMTs = 192;
    for (Int_t iring = 0; iring < 8; iring++){
      for (Int_t icolumn = 0; icolumn < 24; icolumn++){
        pmt_active_flag[iring][icolumn] = 1;
      }
    }
  }

  
  
  // Get vertex position reconstructed by AdSimple (or whatevetr other methods)

  double old_pos[3] = {0,0,0};
  double old_time = 0;
  if(recTrigger.positionStatus()==ReconStatus::kGood){
    const CLHEP::HepLorentzVector oldpos = recTrigger.position();
    old_pos[0] = oldpos.x();
    old_pos[1] = oldpos.y();
    old_pos[2] = oldpos.z();
    old_time = oldpos.t();

    // In rare occasion, the old reconstructed vertex contains NaN
    // This crashes reconstruction.
    if (old_pos[0] != old_pos[0] || old_pos[1] != old_pos[1] || old_pos[2] != old_pos[2]){
      warning() << "\tInitial vertex location conains NaN. Resetting to the origin...." << endreq;
      std::cout << "\t"<< old_pos[0] << " " << old_pos[1] << " " << old_pos[2] << std::endl;
      old_pos[0] = 0;
      old_pos[1] = 0;
      old_pos[2] = 0;
    }
  }
  


  Double_t rsq_start = old_pos[0]*old_pos[0] + old_pos[1]*old_pos[1];
  Double_t z_start = old_pos[2];
  Double_t phi_start = 180.0 * TMath::ATan2(old_pos[1],old_pos[0])  / TMath::Pi();
  //        Double_t rsq_start = 0;
  //        Double_t z_start = 0;
  //        Double_t phi_start = 0;
  if (phi_start < 0) phi_start += 360.;

  
  Double_t chi2_ini = chi2_interpol(rsq_start,z_start,phi_start);
  
  if (m_calculateChi2Only){ // just update the position quality based on the old vertex position, and finish
    recTrigger.setPositionQuality( chi2_ini);
    return StatusCode::SUCCESS;
  }

  total_pe = 0;
  
  for (Int_t iring = 0; iring < 8; iring++){
    for (Int_t icolumn = 0; icolumn < 24; icolumn++){ 
      pe_obs[iring][icolumn] = 0;
    }
  }

  for (Int_t iphi = 0; iphi < phibins; iphi++){
    for (Int_t ir = 0; ir < rbins+1; ir++){
      for (Int_t iz = 0; iz < zbins; iz++){
        chi2_flag[iphi][ir][iz] = 0;
      }
    }
  }
       
  int nHits = 0;
  CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
    chanEnd = pmtReadout->channelReadout().end();
  for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
      chanIter++){
    const CalibReadoutPmtChannel& channel = *chanIter;
    AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
    if( !pmtId.is8inch() ){ continue; }  // Ignore non 8-inch pmt channels
   
    double peakAdc=channel.earliestCharge(-1650,-1250);
    double firstTdc = channel.earliestTime();
    pe_obs[pmtId.ring()-1][pmtId.column()-1] = peakAdc;
    total_pe += peakAdc;
    if (peakAdc > 0){
      nHits++;
    }
  }


  
  if(nHits<1){
    CLHEP::HepLorentzVector position(1e10, 1e10, 1e10, 1e10); // set dummy number in case someone try to use it
    recTrigger.setPosition( position );
    recTrigger.setPositionQuality( 1e10 );
    recTrigger.setPositionStatus( ReconStatus::kNoHits );
    return StatusCode::SUCCESS;
  }


  
  // new reconstrcution

  // correction of the total pe based on the number of active PMTs
  total_pe *= 192./ (Double_t)nActivePMTs;

  
  Double_t chi2_min =1e10;
  
  Int_t ir_min = -1;
  Int_t iz_min = -1;
  Int_t iphi_min = -1;
  
  const Int_t max_step = 100;
  Double_t chi2_min_step[max_step];
  Int_t ir_min_step[max_step];
  Int_t iz_min_step[max_step];
  Int_t iphi_min_step[max_step];
  Int_t nstep = 0;

  if (FullScan){
    for (Int_t ir = 0; ir < rbins+1; ir++){
      for (Int_t iz = 0; iz < zbins; iz++){
        for (Int_t iphi = 0; iphi < phibins; iphi++){
          Double_t chi2_tmp = chi2_grid(iphi, ir, iz);
          if (chi2_tmp > 0 && chi2_min > chi2_tmp) {
            chi2_min = chi2_tmp;
            ir_min = ir;
            iz_min = iz;
            iphi_min = iphi;
          }
        }
        //         std::cout << "---------" << std::endl;
        //         std::cout << ir << "  " << iz << "  " << chi2[ir][iz] << std::endl;
      }
    }
  }else{
       
    Int_t irsq_grid[3];
    Int_t iz_grid[3];
    Int_t iphi_grid[3];
       
    while (nstep < max_step){
      if (nstep == 0){ // first step
        irsq_grid[0] = (Int_t)(rsq_start/rsqbinw+0.5);
        if (irsq_grid[0] > rsqbins -1)
          irsq_grid[0] = rsqbins - 1;
        irsq_grid[1] = irsq_grid[0] + 1;

        iz_grid[0] = (Int_t)((z_start-zmin)/zbinw-0.5);
        if (iz_grid[0] < 0)
          iz_grid[0] = 0;
        if (iz_grid[0] > zbins-2)
          iz_grid[0] = zbins-2;
        iz_grid[1] = iz_grid[0] + 1;
           
        iphi_grid[0] = (Int_t)(phi_start/phibinw - 0.5);
        iphi_grid[1] = iphi_grid[0] + 1;
        iphi_grid[0] = (24 + iphi_grid[0])%24;
        iphi_grid[1] = (24 + iphi_grid[1])%24;

        for (Int_t irsq = 0; irsq < 2; irsq ++){
          for (Int_t iz = 0; iz < 2; iz ++){
            for (Int_t iphi = 0; iphi < 2; iphi ++){
              if (iphi_grid[iphi] < 0 || iphi_grid[iphi] >=  phibins || irsq_grid[irsq] < 0 || irsq_grid[irsq] > rbins || iz_grid[iz] < 0 ||
                  iz_grid[iz] >= zbins){
                std::cout << "something wrong 2! " << iphi_grid[iphi] << " " << irsq_grid[irsq]<< " " << iz_grid[iz] << std::endl;
                std::cout << rsq_start << " " << phi_start << " " << z_start << std::endl;
              }

              Double_t chi2_tmp = chi2_grid(iphi_grid[iphi], irsq_grid[irsq], iz_grid[iz]);
              if (chi2_tmp > 0 && chi2_min > chi2_tmp) {
                chi2_min = chi2_tmp;
                chi2_min_step[nstep] = chi2_tmp;
                ir_min_step[nstep] = irsq_grid[irsq];
                iz_min_step[nstep] = iz_grid[iz];
                iphi_min_step[nstep] = iphi_grid[iphi];
              }
            }
          }
        }

        if (Debug){
          std::cout << " minimization step " << nstep+1
               << " " << chi2_min_step[nstep]
               << " " << ir_min_step[nstep]
               << " " << iz_min_step[nstep]
               << " " << iphi_min_step[nstep]
               << std::endl;
        }
        nstep ++;


      }else{
        for (Int_t irsq = 0; irsq < 3; irsq ++){
          irsq_grid[irsq] = ir_min_step[nstep-1] - 1 + irsq;
          //      if (irsq_grid[irsq]  >= 0 && irsq_grid[irsq] < rsqbins+1){
          if (irsq_grid[irsq] < rsqbins+1){
            for (Int_t iz = 0; iz < 3; iz ++){
              iz_grid[iz] = iz_min_step[nstep-1] - 1 + iz;
              if (iz_grid[iz]  >= 0 && iz_grid[iz] < zbins){
                for (Int_t iphi = 0; iphi < 3; iphi ++){
                  iphi_grid[iphi] = iphi_min_step[nstep-1] - 1 + iphi;
                  if (irsq_grid[irsq] >= 0){
                    iphi_grid[iphi] = (24 + iphi_grid[iphi])%24;
                  }else{
                    iphi_grid[iphi] = (12 + iphi_grid[iphi])%24; // rotate 180 degrees
                    irsq_grid[irsq] = -1 * irsq_grid[irsq];
                  }
                  
                  //                      if (iphi_grid[iphi] < 0 || iphi_grid[iphi] >=  phibins || irsq_grid[irsq] < 0 || irsq_grid[irsq] > rbins || iz_grid[iz] < 0 ||
                  //                          iz_grid[iz] >= zbins){
                  //                        std::cout << "something wrong 3! " << iphi_grid[iphi] << " " << irsq_grid[irsq]<< " " << iz_grid[iz] << std::endl;
                  //                        std::cout << rsq_start << " " << phi_start << " " << z_start << std::endl;
                  //                      }
                  Double_t chi2_tmp = chi2_grid(iphi_grid[iphi], irsq_grid[irsq], iz_grid[iz]);
                  if (chi2_tmp > 0 && chi2_min > chi2_tmp) {
                    chi2_min = chi2_tmp;
                    chi2_min_step[nstep] = chi2_tmp;
                    ir_min_step[nstep] = irsq_grid[irsq];
                    iz_min_step[nstep] = iz_grid[iz];
                    iphi_min_step[nstep] = iphi_grid[iphi];
                  }
                }
              }
            }
          }
        }
           
        if (Debug){
          std::cout << " minimization step " << nstep+1
               << " " << chi2_min_step[nstep]
               << " " << ir_min_step[nstep]
               << " " << iz_min_step[nstep]
               << " " << iphi_min_step[nstep]
               << std::endl;
        }

        if (chi2_min_step[nstep-1] - chi2_min < 1e-10){
          ir_min = ir_min_step[nstep-1];
          iz_min = iz_min_step[nstep-1];
          iphi_min = iphi_min_step[nstep-1];
          if (iz_min < 0 ||iz_min >= zbins){
            std::cout << "something wrong 4! "<< nstep << std::endl;
          }

          break;
        }
        nstep ++;

      }
    }
  }


  Double_t r0;
  if (ir_min == 0){
    r0 = 0;
    // special case for ir_min == 0 to resolve phi_min
    Double_t chi2_second_min = 1e10;
    for (Int_t iPhi = 0; iPhi < phibins; iPhi++){
      Double_t chi2_tmp = chi2_grid(iPhi, ir_min+1, iz_min);
      if (chi2_second_min > chi2_tmp){
        chi2_second_min = chi2_tmp;
        iphi_min = iPhi;
      }
    }
  }else{
    if (Rsq)
      r0 = (ir_min-0.5)*rsqbinw;
    else
      r0 = (ir_min-0.5)*rbinw;
  }
     
  Double_t z0 = (iz_min+0.5)*zbinw +zmin;
  Double_t phi0 = (iphi_min+0.5)*phibinw;

  Double_t RecR = TMath::Sqrt(r0);
  Double_t RecRSQ = r0;
  Double_t RecZ = z0;
  Double_t RecPhi = phi0;

  Double_t chi2_0 = chi2_grid(iphi_min,ir_min,iz_min);
  Double_t chi2_z1;
  Double_t chi2_z2;
  Double_t z1;
  Double_t z2;
       
  Double_t chi2_r1;
  Double_t chi2_r2;
  Double_t r1;
  Double_t r2;
       
  Double_t chi2_phi1;
  Double_t chi2_phi2;
  Double_t phi1;
  Double_t phi2;


  Double_t chi2_recz;
  Double_t chi2_recr;
  Double_t chi2_recphi;
  
     
  if (Interpol){
    
    Double_t A;
    Double_t B;
    
    if (iz_min == 0){
      chi2_z1 = chi2_interpol(r0, z_far_min, phi0);// used extrapolated templated for z_far_min
      chi2_z2 = chi2_grid(iphi_min,ir_min,iz_min+1);
      z1 = z_far_min;
      z2 = (iz_min+1+0.5)*zbinw +zmin;
    }else if (iz_min == zbins-1){
      chi2_z1 = chi2_grid(iphi_min,ir_min,iz_min-1);
      chi2_z2 = chi2_interpol(r0, z_far_max, phi0); // used extrapolated templated for z_far_max
      z1 = (iz_min-1+0.5)*zbinw +zmin;
      z2 = z_far_max;
    }else{
      chi2_z1 = chi2_grid(iphi_min,ir_min,iz_min-1);
      chi2_z2 = chi2_grid(iphi_min,ir_min,iz_min+1);
      z1 = (iz_min-1+0.5)*zbinw +zmin; 
      z2 = (iz_min+1+0.5)*zbinw +zmin; 
    } 
       
    Double_t chi2_z10 = chi2_z1 - chi2_0;
    Double_t chi2_z20 = chi2_z2 - chi2_0;
       

    A = 2 * ((z2-z0)/chi2_z20 - (z1-z0)/chi2_z10);
    B =  (z2*z2-z0*z0)/chi2_z20 - (z1*z1-z0*z0)/chi2_z10;
    RecZ = B/A;

    //     if (RecZ > z_phys_max){
    //       RecZ = z_phys_max;
    //     }
    
    //     if (RecZ < z_phys_min){
    //       RecZ = z_phys_min;
    //     }
    
       
    if (Rsq){
      if (ir_min == 0){
        chi2_r1 = chi2_grid((iphi_min + phibins/2)%phibins,ir_min+1,iz_min); // take negative r
        chi2_r2 = chi2_grid(iphi_min,ir_min+1,iz_min);
        r1 = -1 * (ir_min+1-0.5)*rsqbinw;
        r2 = (ir_min+1-0.5)*rsqbinw;
      }else if (ir_min == rbins){
        chi2_r1 = chi2_grid(iphi_min,ir_min-1,iz_min); 
        chi2_r2 = chi2_interpol(rsq_far_max, z0, phi0); // used extrapolated templated for r_far_max
        r1 = (ir_min-1-0.5)*rsqbinw;
        r2 = rsq_far_max;
      }else{
        chi2_r1 = chi2_grid(iphi_min,ir_min-1,iz_min);
        chi2_r2 = chi2_grid(iphi_min,ir_min+1,iz_min);
        if (ir_min == 1)
          r1 = 0;
        else
          r1 = (ir_min-1-0.5)*rsqbinw;
        r2 = (ir_min+1-0.5)*rsqbinw;
      }
      // Dec. 22, 2011
      // Use liner r to interpolate!
      // r0 is used for z interpolation. So make sure to run this part after z interpolation.
      r0 = TMath::Sqrt(r0);
      if (r1 < 0)
        r1 = - TMath::Sqrt(-r1);
      else
        r1 = TMath::Sqrt(r1);

      r2 = TMath::Sqrt(r2);
      
      
    }else{
      if (ir_min == 0){
        chi2_r1 = chi2_grid((iphi_min + phibins/2)%phibins,ir_min+1,iz_min); // take negative r
        chi2_r2 = chi2_grid(iphi_min,ir_min+1,iz_min);
        r1 = -1 * (ir_min+1-0.5)*rbinw;
        r2 = (ir_min+1-0.5)*rbinw;
      }else if (ir_min == rbins){
        chi2_r1 = chi2_grid(iphi_min,ir_min-1,iz_min);
        chi2_r2 = chi2_interpol(r_far_max, z0, phi0); // used extrapolated templated for r_far_max
        r1 = (ir_min-1-0.5)*rbinw;
        r2 = r_far_max;
      }else{
        chi2_r1 = chi2_grid(iphi_min,ir_min-1,iz_min);
        chi2_r2 = chi2_grid(iphi_min,ir_min+1,iz_min);
        r1 = (ir_min-1-0.5)*rbinw;
        r2 = (ir_min+1-0.5)*rbinw;
      }
    }
       
    Double_t chi2_r10 = chi2_r1 - chi2_0;
    Double_t chi2_r20 = chi2_r2 - chi2_0;
       
    A = 2 * ((r2-r0)/chi2_r20 - (r1-r0)/chi2_r10);
    B =  (r2*r2-r0*r0)/chi2_r20 - (r1*r1-r0*r0)/chi2_r10;
       
//     if (Rsq){
//       RecRSQ = B/A;
//       if (RecRSQ < 0){
//      RecRSQ = - RecRSQ;
//      RecR = -TMath::Sqrt(RecRSQ);
//       }else{
//      RecR = TMath::Sqrt(RecRSQ);
//       }
//     }else{
//       RecR =  B/A;
//       RecRSQ =RecR*RecR;
//     }

    RecR =  B/A;
    RecRSQ =RecR*RecR;

//     if (RecRSQ > rsq_phys_max){
//       RecRSQ = rsq_phys_max;
//       if (RecR > 0)
//      RecR = r_phys_max;
//       else
//      RecR = - r_phys_max;
//     }
    

    Int_t ir_min_for_interpol = ir_min;
    Double_t  chi2_0_for_interpol = chi2_0;
    if (ir_min == 0){
      // because there is no phi dependence for ir_min = 0, and thus interpolation dose not work.
      ir_min_for_interpol = 1; 
      chi2_0_for_interpol = chi2_grid(iphi_min,ir_min_for_interpol,iz_min);
    }
       
    if (iphi_min == 0){
      chi2_phi1 = chi2_grid(iphi_min+23,ir_min_for_interpol,iz_min);
      chi2_phi2 = chi2_grid(iphi_min+1,ir_min_for_interpol,iz_min);
      phi1 = (iphi_min-1+0.5)*phibinw;
      phi2 = (iphi_min+1+0.5)*phibinw;
    }else if (iphi_min == phibins-1){
      chi2_phi1 = chi2_grid(iphi_min-1,ir_min_for_interpol,iz_min);
      chi2_phi2 = chi2_grid(iphi_min-23,ir_min_for_interpol,iz_min);
      phi1 = (iphi_min-1+0.5)*phibinw;
      phi2 = (iphi_min+1+0.5)*phibinw;
    }else{
      chi2_phi1 = chi2_grid(iphi_min-1,ir_min_for_interpol,iz_min);
      chi2_phi2 = chi2_grid(iphi_min+1,ir_min_for_interpol,iz_min);
      phi1 = (iphi_min-1+0.5)*phibinw;
      phi2 = (iphi_min+1+0.5)*phibinw;
    } 
       
    Double_t chi2_phi10 = chi2_phi1 - chi2_0_for_interpol;
    Double_t chi2_phi20 = chi2_phi2 - chi2_0_for_interpol;
         

    A = 2 * ((phi2-phi0)/chi2_phi20 - (phi1-phi0)/chi2_phi10);
    B =  (phi2*phi2-phi0*phi0)/chi2_phi20 - (phi1*phi1-phi0*phi0)/chi2_phi10;
    RecPhi = B/A;
    if (RecPhi < 0) RecPhi += 360.;
    //       if (RecPhi > 0) RecPhi -= 360.;


    // now estimate Minimum Chi2;

    chi2_recz = (pow((z2-RecZ),2)*chi2_z1 - pow((z1-RecZ),2)*chi2_z2)
      /(pow((z2-RecZ),2) - pow((z1-RecZ),2));
//     if(Rsq){
//       chi2_recr = (pow((r2-RecRSQ),2)*chi2_r1 - pow((r1-RecRSQ),2)*chi2_r2)
//      /(pow((r2-RecRSQ),2) - pow((r1-RecRSQ),2));
//     }else{
//       chi2_recr = (pow((r2-RecR),2)*chi2_r1 - pow((r1-RecR),2)*chi2_r2)
//      /(pow((r2-RecR),2) - pow((r1-RecR),2));
//     }
    chi2_recr = (pow((r2-RecR),2)*chi2_r1 - pow((r1-RecR),2)*chi2_r2)
      /(pow((r2-RecR),2) - pow((r1-RecR),2));
    chi2_recphi = (pow((phi2-RecPhi),2)*chi2_phi1 - pow((phi1-RecPhi),2)*chi2_phi2)
      /(pow((phi2-RecPhi),2) - pow((phi1-RecPhi),2));

    if (chi2_min > chi2_recz) chi2_min = chi2_recz;
    if (chi2_min > chi2_recr) chi2_min = chi2_recr;
    if (chi2_min > chi2_recphi) chi2_min = chi2_recphi;
    
    
    
  }
     
     
  Double_t chi2_fine_min =1e10;
  Int_t ir_fine_min = -1;
  Int_t iz_fine_min = -1;
  Int_t iphi_fine_min = -1;
  Double_t rsq_min = 0;
  Double_t r_min = 0;
  Double_t z_min = 0;
  Double_t phi_min = 0;
  if (FineFit){
    if (FullScan){
      for (Int_t ir = 0; ir < rbins_fine; ir++){
        Double_t rsq = (ir+0.5)*rsqbinw_fine;

        for (Int_t iz = 0; iz < zbins_fine; iz++){
          Double_t z = (iz+0.5)*zbinw_fine +zmin;
             
          for (Int_t iphi = 0; iphi < phibins_fine; iphi++){
            Double_t phi = (iphi+0.5)*phibinw_fine;
            Double_t chi2_tmp =  chi2_interpol(rsq, z, phi);

            if (chi2_fine_min > chi2_tmp) {
              chi2_fine_min = chi2_tmp;
              ir_fine_min = ir;
              iz_fine_min = iz;
              iphi_fine_min = iphi;
              rsq_min = rsq;
              z_min = z;
              phi_min = phi;
            }
          }
        }
      }
    }else{

      // use TMinuit minimization here!
      Int_t ierflag;
      Double_t arglist[NPars];
      arglist[0] = 1.0;
         
      minu->mnexcm("SET ERR",arglist,1,ierflag);
         
      //   minu->mnparm(iPar,
      //                parname[iPar].c_str(),
      //                ini_value[iPar],
      //                ini_step[iPar],
      //                0,100,ierflag);
      minu->mnparm(0,
                   "Rsq",
                   rsq_start,
                   10*10,
                   0,3000*3000,ierflag);
      minu->mnparm(1,
                   "Z",
                   z_start,
                   10,
                   -30000,3000,ierflag);
      minu->mnparm(2,
                   "Phi",
                   phi_start,
                   3,
                   -30,390,ierflag);
         
      arglist[0] = 10000;
      arglist[1] = 1.0;
      //   minu->mnexcm("SIMPLEX", arglist ,2,ierflag);
      //   std::cout << "================ SIMPLEX finished with error flag = " << ierflag  << " ============================" << std::endl;
      //   if (ierflag != 0) return;
      arglist[1] = 1.0;
      minu->mnexcm("MIGRAD", arglist ,2,ierflag);
      //  std::cout << "================ MIGRAD finished with error flag = " << ierflag  << " ============================" << std::endl;

      Double_t fpar,ferr;
      minu->GetParameter(0,fpar,ferr);
      rsq_min = fpar;
      minu->GetParameter(1,fpar,ferr);
      z_min = fpar;
      minu->GetParameter(2,fpar,ferr);
      phi_min = fpar;

      chi2_min = chi2_interpol(rsq_min, z_min, phi_min);
         
    }
    if (rsq_min < 0){
      rsq_min = - rsq_min;
      r_min = -TMath::Sqrt(rsq_min);
    }else{
      r_min = TMath::Sqrt(rsq_min);
    }
  }
     
  if (Debug){
    std::cout << r0 << " " << r1 << " " << r2 << " " << chi2_0 << " " << chi2_r1 << " " << chi2_r2 << std::endl;
    std::cout << "test "  << RecRSQ << std::endl;
    std::cout << "test2 " << RecZ<< std::endl;
    std::cout << "test_interpol "<< rsq_min << " " << z_min << " " << phi_min << std::endl;
    for (Int_t ir = 0; ir < rbins+1; ir++){
      for (Int_t iz = 0; iz < zbins; iz++){
        chi2_map->SetBinContent(ir+1,iz+1,chi2_grid(iphi_min,ir,iz));
      }
    }
    chi2_map->Draw("zcol");
    c1->Update();
    c1->Print("chi2_map.ps","Portrait");
    TString dummy;
    std::cin >> dummy;
    if (dummy.Contains("q")) return StatusCode::FAILURE;
    if (FineFit){
      for (Int_t ir = 0; ir < rbins_fine; ir++){
        Double_t rsq = (ir+0.5)*rsqbinw_fine;
        for (Int_t iz = 0; iz < zbins_fine; iz++){
          Double_t z = (iz+0.5)*zbinw_fine +zmin;
          Double_t chi2_tmp = chi2_interpol(rsq, z, phi_min);
          chi2_map_fine->SetBinContent(ir+1,iz+1,chi2_tmp);
        }
      }
      chi2_map_fine->Draw("zcol");
      c1->Update();
      c1->Print("chi2_map.ps","Portrait");
      std::cin >> dummy;
      if (dummy.Contains("q")) return StatusCode::FAILURE;
    }
  }

  
  double pos[3];
  if (FineFit){
    phi_min = TMath::Pi() * phi_min  / 180.;
    pos[0] = r_min * cos(phi_min);
    pos[1] = r_min * sin(phi_min);
    pos[2] = z_min;
  }else{
    RecPhi = TMath::Pi() * RecPhi  / 180.;
    pos[0] = RecR * cos(RecPhi);
    pos[1] = RecR * sin(RecPhi);
    pos[2] = RecZ;
  }


  // Set quality to charge-weighted standard deviation of hit positions
  recTrigger.setPositionQuality( chi2_min);
  // adding protection for having NaN in reconstructed vertex...
  if (pos[0] != pos[0] || pos[1] != pos[1] || pos[2] != pos[2]){
    warning() << "\tReconstructed vertex location conains NaN. Resetting to the origin, and setting flag ReconStatus::kNotConverged." << endreq;
    std::cout << "\t"<< pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
    pos[0] = 0;
    pos[1] = 0;
    pos[2] = 0;
    recTrigger.setPositionStatus( ReconStatus::kNotConverged );
  }else{
    recTrigger.setPositionStatus( ReconStatus::kGood );
  }

  CLHEP::HepLorentzVector position(pos[0], pos[1], pos[2], old_time);
  recTrigger.setPosition( position );
  
  return StatusCode::SUCCESS;
}
StatusCode ChargeTemplatePosTool::initialize ( ) [virtual]

Definition at line 157 of file ChargeTemplatePosTool.cc.

{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName,true);

  // channel quality service
  m_chanqualSvc = svc<IChannelQualitySvc>(m_ChannelQualitySvcName,true);
  nocqCount=0;
  

  Rsq = true;
  Debug = m_debugOutput;
  Debug2 = false;
  Interpol = m_interpolateChi2;
  FineFit = m_minuitMinimization;
  FullScan = false;
  
  
  zmin = -2000;
  zmax = 2000;
  xmin = -2000;
  xmax = 2000;
  
  xbinw = (xmax - xmin)/(Double_t)xbins;
  ybinw = (xmax - xmin)/(Double_t)ybins;
  zbinw = (zmax - zmin)/(Double_t)zbins;
  rbinw = (xmax - 0)/(Double_t)rbins;
  rsqbinw = (xmax*xmax - 0)/(Double_t)rsqbins;
  phibinw = 360./(Double_t)phibins;
  
  
  xbinw_fine = (xmax - xmin)/(Double_t)xbins_fine;
  ybinw_fine = (xmax - xmin)/(Double_t)ybins_fine;
  zbinw_fine = (zmax - zmin)/(Double_t)zbins_fine;
  rbinw_fine = (xmax - 0)/(Double_t)rbins_fine;
  rsqbinw_fine = (xmax*xmax - 0)/(Double_t)rsqbins_fine;
  phibinw_fine = 360./(Double_t)phibins_fine;
  
  
  
  
  info() << "Reading vertex template teble from " << m_templateFileName.c_str() << endreq;
  TFile * f_template = new TFile(m_templateFileName.c_str());
  
  for (Int_t ir = 0; ir < rbins+1; ir++){
    for (Int_t iz = 0; iz < zbins; iz++){
      rz_temp_sum[ir][iz] = 0;
      TH1F * h;
      if (Rsq)
        h = (TH1F*)f_template->Get(Form("/stats/templates/adc_template_rsq%d_z%d",ir,iz));
      else
        h = (TH1F*)f_template->Get(Form("/stats/templates/adc_template_r%d_z%d",ir,iz));
      for (Int_t iring = 0; iring < 8; iring++){
        for (Int_t icolumn = 0; icolumn < 24; icolumn++){
          //          Int_t ich = iring*24 + icolumn;
          Int_t ich = iring*24 + (icolumn + 12)%24;  //temporally fix for 180-degree rotation of the corrdinates
          rz_temp[ir][iz][iring][icolumn] = h->GetBinContent(ich+1);
          rz_temp_sum[ir][iz] += h->GetBinContent(ich+1);
        }
      }
    }
  }

  // normalize
  for (Int_t ir = 0; ir < rbins+1; ir++){
    for (Int_t iz = 0; iz < zbins; iz++){
      for (Int_t iring = 0; iring < 8; iring++){
        for (Int_t icolumn = 0; icolumn < 24; icolumn++){
          rz_temp[ir][iz][iring][icolumn] /= rz_temp_sum[ir][iz];
        }
      }
    }
  
  }

  if (Debug)  c1 = new TCanvas("c1","c1",600,600);

  if(Rsq){
    chi2_map = new TH2F("chi2_map","chi2_map",rsqbins,0,xmax*xmax,zbins,zmin,zmax);
  }else{
    chi2_map = new TH2F("chi2_map","chi2_map",rbins,0,xmax,zbins,zmin,zmax);
  }
  
  if(Rsq){
    chi2_map_fine = new TH2F("chi2_map_fine","chi2_map_fine",rsqbins_fine,0,xmax*xmax,zbins_fine,zmin,zmax);
  }else{
    chi2_map_fine = new TH2F("chi2_map_fine","chi2_map_fine",rbins_fine,0,xmax,zbins_fine,zmin,zmax);
  }
  
  if (Debug){
    c1->Print("chi2_map.ps[","Portrait");
  }


  

  minu = new TMinuit(NPars);
  if (!Debug)
    minu->SetPrintLevel(-1);
  
  minu->SetFCN(minuit_fcn);

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

Definition at line 265 of file ChargeTemplatePosTool.cc.

{
  //Status report for null Channel Quality Service
  if (nocqCount > 0)
    warning() << "******* Got null Channel Quality objects for "<< nocqCount << " times  *********" << endreq;

  return StatusCode::SUCCESS;
}

Member Data Documentation

Definition at line 46 of file ChargeTemplatePosTool.h.

Definition at line 49 of file ChargeTemplatePosTool.h.

Definition at line 51 of file ChargeTemplatePosTool.h.

Definition at line 54 of file ChargeTemplatePosTool.h.

Definition at line 57 of file ChargeTemplatePosTool.h.

Definition at line 61 of file ChargeTemplatePosTool.h.

Definition at line 63 of file ChargeTemplatePosTool.h.

Definition at line 65 of file ChargeTemplatePosTool.h.

Definition at line 67 of file ChargeTemplatePosTool.h.

Definition at line 69 of file ChargeTemplatePosTool.h.

Definition at line 70 of file ChargeTemplatePosTool.h.

Definition at line 71 of file ChargeTemplatePosTool.h.

Definition at line 72 of file ChargeTemplatePosTool.h.

Definition at line 73 of file ChargeTemplatePosTool.h.

Definition at line 74 of file ChargeTemplatePosTool.h.

Definition at line 77 of file ChargeTemplatePosTool.h.

Definition at line 78 of file ChargeTemplatePosTool.h.

Definition at line 79 of file ChargeTemplatePosTool.h.

Definition at line 80 of file ChargeTemplatePosTool.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:19 for ChargeTemplatePos by doxygen 1.7.4