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

In This Package:

Spallation.cc
Go to the documentation of this file.
00001 #include "Spallation.h"
00002 #include "Event/DaqCrate.h"
00003 #include "Event/DaqPmtCrate.h"
00004 #include "Event/DaqPmtChannel.h"
00005 #include "Event/ReadoutHeader.h"
00006 #include "Event/ReadoutPmtCrate.h"
00007 #include "Event/CalibReadoutHeader.h"
00008 #include "Event/CalibReadoutPmtCrate.h"
00009 #include "StatisticsSvc/IStatisticsSvc.h"
00010 #include "DataSvc/ICableSvc.h"
00011 #include "DataSvc/ICalibDataSvc.h"
00012 #include "DetHelpers/IPmtGeomInfoSvc.h"
00013 #include "DetHelpers/ICoordSysSvc.h"
00014 #include "Context/ServiceMode.h"
00015 #include "TH1F.h"
00016 #include "TF1.h"
00017 #include "TH2.h"
00018 #include "TParameter.h"
00019 #include <math.h>
00020 #include <numeric>
00021 #include "Event/RawEventHeader.h"
00022 #include "Event/RawRom.h"
00023 #include "Event/RawRomLtb.h"
00024 #include "Event/RawPmtChannel.h"
00025 #include "CLHEP/Vector/ThreeVector.h"
00026 #include "CLHEP/Units/SystemOfUnits.h"
00027 #include <fstream>
00028 #include "TMath.h"
00029 
00030 using namespace std;
00031 using namespace DayaBay;
00032 
00033 
00034 Spallation::Spallation(const std::string& type,
00035                                          const std::string& name, 
00036                                          const IInterface* parent)
00037   : GaudiTool(type,name,parent)
00038   , m_cableSvc(0)
00039   , m_statsSvc(0)
00040 {
00041   
00042   declareInterface< IEnergyCalibTool >(this) ;   
00043   declareProperty("MuonCut",m_MuonCut=10000,"PE threshold of a muon");
00044   declareProperty("TimeLow",m_TimeLow=2000,"small time for spallation neutron selection");
00045   declareProperty("TimeHigh",m_TimeHigh=200000,"large time for spallation neutron selection");
00046   declareProperty("CableSvcName",m_cableSvcName="CableSvc",
00047                   "Name of service to map between detector, hardware, and electronic IDs");
00048   declareProperty("OutFilePath",m_outfilepath="/file1/Spallation",
00049                   "File path of registered histograms.");
00050   declareProperty("InFilePath",m_infilepath="/file0/Background",
00051                   "File path of background histograms.");
00052   declareProperty("TDCCutLow",m_tdclow=900,"low tdc cut for good hits");
00053   declareProperty("TDCCutHigh",m_tdchigh=1050,"high tdc cut for good hits");
00054   declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
00055       "Name of Pmt Geometry Information Service");
00056 
00057 }
00058 
00059 Spallation::~Spallation(){}
00060 
00061 StatusCode Spallation::initialize(){
00062          info() << "initialize()" << endreq;
00063   StatusCode sc = this->GaudiTool::initialize();
00064   if( sc != StatusCode::SUCCESS ) return sc;
00065   
00066   //Get PMT geometry service
00067   m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);
00068    if(!m_pmtGeomSvc ){
00069     error() << "Failed to access PMT geometry svc: " << m_pmtGeomSvcName << endreq;
00070     return StatusCode::FAILURE;
00071   }
00072   // Get Cable Service
00073   m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
00074   if(!m_cableSvc){
00075     error() << "Failed to access cable svc: " << m_cableSvcName << endreq;
00076     return StatusCode::FAILURE;
00077   }
00078   // Get the histogram service
00079   m_statsSvc = svc<IStatisticsSvc>("StatisticsSvc",true);
00080   if(!m_statsSvc) {
00081     error() << " No StatisticsSvc available." << endreq;
00082     return StatusCode::FAILURE;
00083   }
00084   first = true;
00085     return StatusCode::SUCCESS;
00086 }
00087 
00088 StatusCode Spallation::finalize()
00089 {
00090   StatusCode sc = this->GaudiTool::finalize();
00091   if( sc != StatusCode::SUCCESS ) return sc;
00092 
00093   return StatusCode::SUCCESS;
00094 }
00095 
00096 StatusCode Spallation::process(const DayaBay::CalibReadoutHeader& calibReadoutHeader,int st, int acu, int zpo){
00097         
00098         CalibReadout* calibReadout = const_cast<CalibReadout*>(calibReadoutHeader.calibReadout());
00099         const DayaBay::Detector& det = calibReadout->detector();
00100         const TimeStamp TriTime = calibReadout->triggerTime();
00101         Context context = calibReadoutHeader.context();
00102         if( !hasStats(det) ){
00103     // New detector; initialize all the detector statistics
00104                 StatusCode sc = prepareStats(context);
00105         if( sc != StatusCode::SUCCESS ) return sc;
00106     }
00107    TH1F* Spectrum = m_statsSvc->getTH1F( this->getOutPath(det) + "/SpallationRawSpectrum");
00108    TH2F* TimeToLastMuon = m_statsSvc->getTH2F( this->getOutPath(det) + "/TimeToLastMuon");
00109    TH2F* NeutronMultiplicity = m_statsSvc->getTH2F( this->getOutPath(det) + "/NeutronMultiplicity");
00110    TH1F* Interval = m_statsSvc->getTH1F( this->getOutPath(det) + "/TimeInterval");
00111         if(first){
00112                 m_vetoMuonTime = TriTime;
00113                 m_AD1MuonTime = TriTime;
00114                 m_AD2MuonTime = TriTime;
00115                 m_AD3MuonTime = TriTime;
00116                 m_AD4MuonTime = TriTime;
00117                 first = false;
00118                 return StatusCode::SUCCESS;
00119         }
00121         if(det.isWaterShield() || det.detectorId() == DetectorId::kRPC){
00122                 m_vetoMuonTime = TriTime;
00123                 return StatusCode::SUCCESS;
00124         }
00126         CalibReadoutPmtCrate* calibCrate = dynamic_cast<CalibReadoutPmtCrate*>(calibReadout);
00127         CalibReadoutPmtCrate::PmtChannelReadouts calibChannels = calibCrate->channelReadout();
00128         double qsum = 0;
00129         CalibReadoutPmtCrate::PmtChannelReadouts::iterator it;
00130         for (it=calibChannels.begin(); it != calibChannels.end(); it++) {
00131                 std::vector<double> charge = it->charge();
00132                 std::vector<double> time = it->time();
00133                 double largecharge=0;
00134         for(int i = 0; i<time.size();i++){
00135             //cout<<time[i]<<endl;      //////more precise cuts research////
00136             if(time[i]<-1*m_tdclow*1.5625&&time[i]>-1*m_tdchigh*1.5625&&charge[i]>largecharge){
00137                 largecharge = charge[i];
00138             }
00139         }
00140                 //double largecharge = it->maxCharge();
00141         qsum+=largecharge;
00142         }
00143         if(qsum > m_MuonCut) {
00144                 if(det.detectorId() == 1) { m_AD1MuonTime = TriTime; m_AD1MuonCharge = qsum;}
00145                 else if(det.detectorId() == 2)  { m_AD2MuonTime = TriTime;m_AD2MuonCharge = qsum;}
00146                 else if(det.detectorId() == 3)  { m_AD3MuonTime = TriTime;m_AD3MuonCharge = qsum;}
00147                 else if(det.detectorId() == 4)  { m_AD4MuonTime = TriTime;m_AD4MuonCharge = qsum;}
00148                 else ;
00149                 return StatusCode::SUCCESS;
00150         }
00151         if(det.detectorId() == 1&&( TriTime.GetNanoSec()-m_AD1MuonTime.GetNanoSec() )< m_TimeHigh && ( TriTime.GetNanoSec()-m_AD1MuonTime.GetNanoSec() )> m_TimeLow){
00152                 Spectrum->Fill(qsum);
00153                 TimeToLastMuon->Fill(TriTime.GetNanoSec()-m_AD1MuonTime.GetNanoSec(), qsum);
00154                 if(qsum>1500&&qsum<1800) Interval->Fill(TriTime.GetNanoSec()-m_AD1MuonTime.GetNanoSec());
00155         }
00156         if(det.detectorId() == 2&&( TriTime.GetNanoSec()-m_AD2MuonTime.GetNanoSec() )< m_TimeHigh && ( TriTime.GetNanoSec()-m_vetoMuonTime.GetNanoSec() )> m_TimeLow){
00157                 Spectrum->Fill(qsum);
00158                 TimeToLastMuon->Fill(TriTime.GetNanoSec()-m_AD2MuonTime.GetNanoSec(), qsum);
00159                  if(qsum>1500&&qsum<1800) Interval->Fill(TriTime.GetNanoSec()-m_AD2MuonTime.GetNanoSec());      
00160         }
00161         if(det.detectorId() == 3&&( TriTime.GetNanoSec()-m_AD3MuonTime.GetNanoSec() )< m_TimeHigh && ( TriTime.GetNanoSec()-m_vetoMuonTime.GetNanoSec() )> m_TimeLow){
00162                 Spectrum->Fill(qsum);
00163                 TimeToLastMuon->Fill(TriTime.GetNanoSec()-m_AD3MuonTime.GetNanoSec(), qsum);
00164                 if(qsum>1500&&qsum<1800) Interval->Fill(TriTime.GetNanoSec()-m_AD3MuonTime.GetNanoSec());
00165         }
00166         if(det.detectorId() == 4&&( TriTime.GetNanoSec()-m_AD4MuonTime.GetNanoSec() )< m_TimeHigh && ( TriTime.GetNanoSec()-m_vetoMuonTime.GetNanoSec() )> m_TimeLow){
00167                 Spectrum->Fill(qsum);
00168                 TimeToLastMuon->Fill(TriTime.GetNanoSec()-m_AD4MuonTime.GetNanoSec(), qsum);
00169                 if(qsum>1500&&qsum<1800)        Interval->Fill(TriTime.GetNanoSec()-m_AD4MuonTime.GetNanoSec());
00170         }
00171         return StatusCode::SUCCESS; 
00172 }
00173 
00174 StatusCode Spallation::calibrate(float TimeInterval,int st, int acu, int zpo,int runN){
00175         std::vector<DayaBay::Detector>::iterator detIter, 
00176    detEnd = m_processedDetectors.end();
00177         for(detIter = m_processedDetectors.begin(); detIter != detEnd; detIter++){
00178                 DayaBay::Detector detector = *detIter;
00179                 if(detector.isAD()){
00180                         TObject* nTime = m_statsSvc->get(this->getOutPath(detector) + "/Time");
00181                         TParameter<float>* nTimePar= dynamic_cast<TParameter<float>*>(nTime);
00182                         nTimePar->SetVal(TimeInterval);
00183                 }
00184         }
00185         return StatusCode::SUCCESS;
00186 }
00187 
00188 bool Spallation::hasStats(const DayaBay::Detector& detector){
00189   // Check if statistics have been initialized for this detector
00190         return (std::find(m_processedDetectors.begin(), m_processedDetectors.end(), detector) != m_processedDetectors.end());
00191 }
00192 
00193 std::string Spallation::getOutPath( const DayaBay::Detector& detector ){
00194         return m_outfilepath + "/" + detector.detName();
00195 }
00196 std::string Spallation::getInPath( const DayaBay::Detector& detector ){
00197         return m_infilepath + "/" + detector.detName();
00198 }
00199 StatusCode Spallation::prepareStats(const Context& context){
00200         DayaBay::Detector detector(context.GetSite(), context.GetDetId());
00202   {
00203     std::ostringstream title, path;
00204     std::string name = "SpallationRawSpectrum";
00205     title << "Spallation neutrons Raw Spectrum in " << detector.detName();
00206     path << this->getOutPath(detector) << "/" << name;
00207     TH1F* RawSpectrum = new TH1F(name.c_str(),title.str().c_str(),1000,0,10000); // to hold water PMT fee board 
00208     RawSpectrum->GetXaxis()->SetTitle("P.E number after calibration");
00209     RawSpectrum->GetYaxis()->SetTitle("Entries");
00210     if( m_statsSvc->put(path.str(),RawSpectrum).isFailure() ) {
00211       error() << "prepareStats(): Could not register " << path << endreq;
00212        delete RawSpectrum;
00213        return StatusCode::FAILURE;
00214     }
00215     }
00216 
00217         // Time latence in the run
00218   
00219     std::string name = "Time";
00220     std::ostringstream path;
00221     path << this->getOutPath(detector) << "/" << name;
00222     TParameter<float>* par = new TParameter<float>(name.c_str(), 
00223                                                0);
00224     if( m_statsSvc->put(path.str(),par).isFailure() ) {
00225       error() << "prepareStats(): Could not register " << name 
00226               << " at " << path << endreq;
00227       delete par;
00228       par = 0;
00229       return StatusCode::FAILURE;
00230     }
00232         {
00233                 std::ostringstream title, path;
00234                 std::string name = "TimeToLastMuon";
00235                 title << "Time to last muon vs energy in " << detector.detName();
00236                 path << this->getOutPath(detector) << "/" << name;
00237                 TH2F* timetolast = new TH2F(name.c_str(),title.str().c_str(),200,0,200000,1000,0.5,20000); 
00238                 timetolast->GetXaxis()->SetTitle("Time to last muon in ns");
00239                 timetolast->GetYaxis()->SetTitle("Total charge of the event");
00240                 if( m_statsSvc->put(path.str(),timetolast).isFailure() ) {
00241                         error() << "prepareStats(): Could not register " << path << endreq;
00242                         delete timetolast;
00243                         return StatusCode::FAILURE;
00244                         }
00245         }
00247         {
00248                 std::ostringstream title, path;
00249                 std::string name = "NeutronMultiplicity";
00250                 title << "Neutron Multiplicity in " << detector.detName();
00251                 path << this->getOutPath(detector) << "/" << name;
00252                 TH2F* multi = new TH2F(name.c_str(),title.str().c_str(),100,0,100,10000,0,6000000); 
00253                 multi->GetXaxis()->SetTitle("neutron multiplicity");
00254                 multi->GetYaxis()->SetTitle("muon total charge");
00255                 if( m_statsSvc->put(path.str(),multi).isFailure() ) {
00256                         error() << "prepareStats(): Could not register " << path << endreq;
00257                         delete multi;
00258                         return StatusCode::FAILURE;
00259                         }
00260         }
00262         {
00263                 std::ostringstream title, path;
00264                 std::string name = "TimeInterval";
00265                 title << "Neutron time interval to last muon in " << detector.detName();
00266                 path << this->getOutPath(detector) << "/" << name;
00267                 TH1F* interval = new TH1F(name.c_str(),title.str().c_str(),200,0,200000); 
00268                 interval->GetXaxis()->SetTitle("Time to last muon in ns");
00269                 interval->GetYaxis()->SetTitle("entries");
00270                 if( m_statsSvc->put(path.str(),interval).isFailure() ) {
00271                         error() << "prepareStats(): Could not register " << path << endreq;
00272                         delete interval;
00273                         return StatusCode::FAILURE;
00274                         }
00275         }
00276   m_processedDetectors.push_back(detector);
00277 
00278   return StatusCode::SUCCESS;
00279 }
00280 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:03:20 for CalibEnergy by doxygen 1.7.4