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

In This Package:

TimeRecoTool.cc
Go to the documentation of this file.
00001 #include "TimeRecoTool.h"
00002 
00003 #include "Event/CalibReadout.h"
00004 #include "Event/CalibReadoutPmtCrate.h"
00005 #include "Event/CalibReadoutPmtChannel.h"
00006 #include "Event/RecTrigger.h"
00007 #include "Event/UserDataHeader.h"
00008 
00009 #include "Conventions/Reconstruction.h"
00010 #include "Conventions/Detectors.h"
00011 
00012 #include "DataSvc/ICableSvc.h"
00013 #include "DataSvc/ICalibDataSvc.h"
00014 #include "DybKernel/IPmtCalibSvc.h"
00015 #include "DetHelpers/IPmtGeomInfoSvc.h"
00016 #include "CLHEP/Vector/LorentzVector.h"
00017 
00018 #include "DybKernel/IChannelQualitySvc.h"
00019 
00020 #include <iostream>
00021 #include <fstream>
00022 #include <math.h>
00023 #include "TROOT.h"
00024 #include "TStyle.h"
00025 #include "TFile.h"
00026 #include "TChain.h"
00027 #include "TTree.h"
00028 #include "TF1.h"
00029 #include "TH2.h"
00030 #include "TMath.h"
00031 #include "TCanvas.h"
00032 #include "TMinuit.h"
00033 
00034 #define NOCQ_ERRORS 10
00035 
00036 //Number of steps for Newton-Raphson Iteration
00037 const Int_t maxstep_t = 20;
00038 
00039 //Number of attempts to broaden PMT selection window for reach convergence
00040 const Int_t maxselect_t = 5;
00041 Int_t nPMT_t[maxselect_t];
00042 
00043 
00044 // Physical boundary for the reconstructed position
00045 Double_t z_phys_min_t = -2400;
00046 Double_t z_phys_max_t = 2400;
00047 Double_t x_phys_max_t = 2300;
00048 Double_t x_phys_min_t = -2300;
00049 Double_t y_phys_max_t = 2300;
00050 Double_t y_phys_min_t = -2300;
00051 Double_t r_phys_max_t = 2300;
00052 
00053 Double_t total_pe_t;
00054 Double_t early_deltime_t;
00055 Double_t median_deltime_t;
00056 Double_t deltime_t[8][24];
00057 
00058 Double_t pmt_charge_t[8][24];
00059 Double_t pmt_time_t[8][24];
00060 Double_t pmt_x_t[8][24];
00061 Double_t pmt_y_t[8][24];
00062 Double_t pmt_z_t[8][24];
00063 Double_t pmt_error_t[8][24];
00064 Int_t pmt_flag_t[8][24];
00065 
00066 Double_t initial_pos_t[4];
00067 Double_t initial_energy_t;
00068 
00069 Int_t totalsteps;
00070 Double_t final_pos_t[4];
00071 
00072 Double_t chi2_sum_t(Double_t x, Double_t y, Double_t z, Double_t t, Int_t iSelect);
00073 Double_t chi2_time_t(Double_t x, Double_t y, Double_t z, Double_t t, Int_t iRing, Int_t iColumn);
00074 void IterativeFit(Int_t n, Int_t maxstep, Double_t x, Double_t y, Double_t z, Int_t iSelect);
00075 void SimpleStart();
00076 Double_t SimpleR(Int_t ring, Int_t column);
00077 void CalculateDeltime(Double_t x, Double_t y, Double_t z);
00078 void SelectPMT();
00079 
00080 Int_t entry;
00081 
00082 using namespace DayaBay;
00083 
00084 TimeRecoTool::TimeRecoTool(const std::string& type,
00085                            const std::string& name, 
00086                            const IInterface* parent)
00087   : GaudiTool(type,name,parent)
00088   , m_cableSvc_t(0)
00089   , m_calibDataSvc_t(0)
00090   , m_pmtGeomSvc_t(0)
00091   , m_chanqualSvc_t(0)
00092 {
00093 
00094 
00095     declareInterface< IReconTool >(this) ;   
00096     
00097     declareProperty("CableSvcName",m_cableSvcName_t="CableSvc",
00098                     "Name of service to map between detector, hardware, and electronic IDs");
00099     declareProperty("CalibDataSvcName", m_calibDataSvcName_t="DybPmtCalibSvc",
00100                     "Name of calibration data service");
00101     declareProperty("PmtGeomSvcName", m_pmtGeomSvcName_t="PmtGeomInfoSvc",
00102                     "Name of PMT geometry data service");
00103     declareProperty("DebugOutput", m_debugOutput_t=false,
00104                     "Output debug informaiton");
00105     declareProperty("ChannelQualitySvcName", m_ChannelQualitySvcName_t="DybChannelQualitySvc",
00106                     "Name of channel quality service");
00107     declareProperty("UseChannelQualitySvc", m_useChannelQualitySvc_t=true,
00108                     "Use ChannelQualityService to check if the PMTs are alive. Turn it off if it's too slow");
00109 
00110 }
00111 
00112 TimeRecoTool::~TimeRecoTool(){}
00113 
00114 StatusCode TimeRecoTool::initialize()
00115 {
00116   
00117   entry = 0;
00118 
00119   // Get Cable Service
00120   m_cableSvc_t = svc<ICableSvc>(m_cableSvcName_t,true);
00121   m_calibDataSvc_t = svc<IPmtCalibSvc>(m_calibDataSvcName_t,true);
00122   m_pmtGeomSvc_t = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName_t,true);
00123 
00124   // channel quality service
00125   m_chanqualSvc_t = svc<IChannelQualitySvc>(m_ChannelQualitySvcName_t,true);
00126   nocqCount=0;
00127 
00128   Debug_t = m_debugOutput_t;
00129   
00130   if (Debug_t) {
00131     
00132     //Debugging tools
00133     TDirectory * dir = gDirectory;
00134     savefile = new TFile("debug.root","RECREATE");
00135 
00136     //fittime1 = new TF1("fittime1","([0]+(1.465/300)*sqrt(pow(([1]*cos([2])+2250*cos(((x-.5)/12)*3.141592)),2)+pow(([1]*sin([2])+2250*sin(((x-.5)/12)*3.141592)),2)+pow(([3]-500.0*(int((x-1)/24)-3.5)),2)))",1,192); 
00137     //fittime0 = new TF1("fittime0","([0]+(1.465/300)*sqrt(pow(([1]*cos([2])+2250*cos(((x-.5)/12)*3.141592)),2)+pow(([1]*sin([2])+2250*sin(((x-.5)/12)*3.141592)),2)+pow(([3]-500.0*(int((x-1)/24)-3.5)),2)))",1,192); 
00138 
00139     fittime1 = new TF1("fittime1","([0]+(1.465/300)*sqrt(pow(([1]+2250*cos(((x-.5)/12)*3.141592)),2)+pow(([2]+2250*sin(((x-.5)/12)*3.141592)),2)+pow(([3]-500.0*(int((x-1)/24)-3.5)),2)))",1,192); 
00140     fittime0 = new TF1("fittime0","([0]+(1.465/300)*sqrt(pow(([1]+2250*cos(((x-.5)/12)*3.141592)),2)+pow(([2]+2250*sin(((x-.5)/12)*3.141592)),2)+pow(([3]-500.0*(int((x-1)/24)-3.5)),2)))",1,192); 
00141     reftime0 = new TF1("reftime0","[0]",1,192); 
00142     reftime1 = new TF1("reftime1","[0]",1,192); 
00143 
00144     hchantime1 = new TH1F("hchantime1","earlytime vs channum",192,0.5,192.5);
00145     hchantime0 = new TH1F("hchantime0","earlytime vs channum",192,0.5,192.5);
00146     hchancharge1 = new TH1F("hchancharge1","charge vs channum",192,0,192);
00147     hchancharge0 = new TH1F("hchancharge0","charge vs channum",192,0,192);
00148     
00149     hdifftime0 = new TH1F("hdifftime0","hdifftime0",300,-50,250);
00150     hchantime0->SetLineColor(3);
00151     hchantime0->SetMarkerStyle(4);
00152 
00153     hchancharge0->SetLineColor(3);
00154 
00155     ccharge = new TCanvas();
00156     ccharge->SetBorderMode(0);
00157     ccharge->SetFrameBorderMode(0);
00158     ccharge->SetFillColor(0);
00159 
00160     ctime = new TCanvas();
00161     ctime->SetBorderMode(0);
00162     ctime->SetFrameBorderMode(0);
00163     ctime->SetFillColor(0);
00164     
00165     hsteps = new TH1F("hsteps","hsteps",maxstep_t+1,0,maxstep_t+1);
00166     hstepsVdDis = new TH2F("hstepsVdDis","hstepsVdDis",maxstep_t+1,0,maxstep_t+1,100,0,1000);
00167     hstepsVchi2 = new TH2F("hstepsVchi2","hstepsVchi2",maxstep_t+1,0,maxstep_t+1,50,0,1);
00168 
00169     dir->cd();
00170     
00171   }
00172  
00173   return StatusCode::SUCCESS;
00174 }
00175 
00176 StatusCode TimeRecoTool::finalize()
00177 {
00178   //Status report for null Channel Quality Service
00179   if (nocqCount > 0)
00180     warning() << "******* Got null Channel Quality objects for "<< nocqCount << " times  *********" << endreq;
00181  
00182  if (Debug_t){
00183      
00184    TDirectory * dir = gDirectory;
00185    savefile->cd();
00186    hsteps->Write();
00187    hstepsVdDis->Write();
00188    hstepsVchi2->Write();
00189    dir->cd();
00190 
00191     savefile->Close();
00192   }
00193 
00194   return StatusCode::SUCCESS;
00195 }
00196 
00197 StatusCode TimeRecoTool::reconstruct(const CalibReadout& readout,
00198                                      RecTrigger& recTrigger)
00199 {
00200   
00201   entry ++;
00202   
00203   if( !readout.detector().isAD() ){
00204     debug() << "Not an AD readout; ignoring detector " 
00205             << readout.detector().detName() << endreq;
00206     recTrigger.setPositionStatus( ReconStatus::kNotProcessed );
00207     return StatusCode::SUCCESS;
00208   }
00209 
00210   const CalibReadoutPmtCrate* pmtReadout 
00211     = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
00212   if(!pmtReadout){
00213     error() << "Incorrect type of readout crate for detector "
00214             << readout.detector().detName() << endreq;
00215     recTrigger.setPositionStatus( ReconStatus::kBadReadout );
00216     return StatusCode::FAILURE;
00217   }
00218 
00219 
00220   // Context for this data
00221   int task = 0;
00222   ServiceMode svcMode(readout.header()->context(), task);
00223 
00224   Int_t nActivePMTs_t = 0;
00225 
00226   if (m_useChannelQualitySvc_t){
00227     const IChannelQuality * cq = m_chanqualSvc_t->channelQuality(svcMode);
00228     if (!cq){
00229       nocqCount ++;
00230       if (nocqCount <= NOCQ_ERRORS){
00231         warning() << "\tGot null channel quality object for " << readout.header()->context().AsString()
00232                   << ". Assuming all channels are good." << endreq;
00233       }
00234       if (nocqCount == NOCQ_ERRORS){
00235         warning() << "\tFurther warning for channel quality service will be suppressed" << endreq;
00236       }
00237 
00238       nActivePMTs_t = 192;
00239       for (Int_t iring = 0; iring < 8; iring++){
00240         for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00241           pmt_flag_t[iring][icolumn] = 1;
00242         }
00243       }
00244 
00245     }else{
00246       
00247       // loop over channels and check if PMTs are active
00248       IChannelQuality::ChannelSet_t chans = cq->channels();
00249       for (int i = 0; i < chans.size(); i++){
00250         DayaBay::FeeChannelId chid = chans[i];
00251         AdPmtSensor pmtId = m_cableSvc_t->adPmtSensor(chid, svcMode);
00252         if (pmtId.is8inch()){
00253           if (cq->good(chid)){
00254             pmt_flag_t[pmtId.ring()-1][pmtId.column()-1] = 1;
00255             nActivePMTs_t ++;
00256           }else{
00257             pmt_flag_t[pmtId.ring()-1][pmtId.column()-1] = -1;
00258             //    std::cout << "\tGod bad one: ring " << pmtId.ring() << " column " << pmtId.column() << std::endl;
00259           
00260           }
00261         }
00262       }
00263     }
00264     //    std::cout << "Number of active PMTs = " << nActivePMTs << std::endl;
00265   }else{ // all PMTs are assumed to be active
00266     nActivePMTs_t = 192;
00267     for (Int_t iring = 0; iring < 8; iring++){
00268       for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00269         pmt_flag_t[iring][icolumn] = 1;
00270       }
00271     }
00272   }
00273  
00274   for (Int_t iring = 0; iring < 8; iring++){
00275     for (Int_t icolumn = 0; icolumn < 24; icolumn++){ 
00276       
00277       pmt_charge_t[iring][icolumn] = 0;
00278       pmt_time_t[iring][icolumn] = 0;
00279       pmt_error_t[iring][icolumn] = 0; 
00280       deltime_t[iring][icolumn] = 0;
00281     
00282     }
00283   }
00284   
00285   //Loop over channels 
00286   //******************
00287   Int_t nHits_t = 0;
00288   total_pe_t = 0;
00289   CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
00290     chanEnd = pmtReadout->channelReadout().end();
00291   for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
00292       chanIter++){
00293     const CalibReadoutPmtChannel& channel = *chanIter;
00294     AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
00295     if( !pmtId.is8inch() ){ continue; }  // Ignore non 8-inch pmt channels
00296     const PmtCalibData* pmtCalib = m_calibDataSvc_t->pmtCalibData(pmtId, svcMode);
00297     if( !pmtCalib ){ 
00298       error() << "No calibration data for pmt ID: " << pmtId << endreq;
00299       recTrigger.setPositionStatus( ReconStatus::kBadReadout );
00300       return StatusCode::FAILURE;
00301     }
00302     const CLHEP::Hep3Vector& pmtPos 
00303       = m_pmtGeomSvc_t->get(pmtId.fullPackedData())->localPosition();
00304     double peakAdc=channel.earliestCharge(-1650,-1250);
00305     double firstTdc = channel.earliestTime();
00306     if(pmt_flag_t[pmtId.ring()-1][pmtId.column()-1] == 1){
00307 
00308       pmt_charge_t[pmtId.ring()-1][pmtId.column()-1] = peakAdc;
00309       pmt_time_t[pmtId.ring()-1][pmtId.column()-1] = firstTdc;
00310 
00311       pmt_x_t[pmtId.ring()-1][pmtId.column()-1] = pmtPos[0];
00312       pmt_y_t[pmtId.ring()-1][pmtId.column()-1] = pmtPos[1];
00313       pmt_z_t[pmtId.ring()-1][pmtId.column()-1] = pmtPos[2];
00314 
00315       if (peakAdc>0){
00316         nHits_t++;
00317         total_pe_t += peakAdc;
00318       }
00319  
00320     }
00321   }
00322   
00323   if(nHits_t<1){
00324     CLHEP::HepLorentzVector position(1e10, 1e10, 1e10, 1e10); // set dummy number in case someone try to use it
00325     recTrigger.setPosition( position );
00326     recTrigger.setPositionQuality( 1e10 );
00327     recTrigger.setPositionStatus( ReconStatus::kNoHits ); //No Hits
00328     return StatusCode::SUCCESS;
00329   }
00330 
00331  // correction of the total pe based on the number of active PMTs
00332   total_pe_t *= 192./ (Double_t)nActivePMTs_t;
00333 
00334   //Initialize final position
00335   final_pos_t[0] = 1e10;
00336   final_pos_t[1] = 1e10;
00337   final_pos_t[2] = 1e10;
00338   final_pos_t[3] = 1e10;
00339 
00340   //Load Prefit Vertex (AdSimple), otherwise use a simple Charge guess
00341   //****************************************************************** 
00342   if(recTrigger.positionStatus()==ReconStatus::kGood){
00343     const CLHEP::HepLorentzVector oldpos = recTrigger.position();
00344     initial_pos_t[0] = oldpos.x();
00345     initial_pos_t[1] = oldpos.y();
00346     initial_pos_t[2] = oldpos.z();
00347     initial_pos_t[3] = oldpos.t();
00348 
00349     //Check if initial r is too close to boundary
00350     Double_t r = sqrt(pow(initial_pos_t[0],2) + pow(initial_pos_t[1],2));
00351     if (r > 1800){
00352       Double_t phi = atan2(initial_pos_t[1],initial_pos_t[0]);
00353       initial_pos_t[0] = 1800*cos(phi);
00354       initial_pos_t[1] = 1800*sin(phi);
00355     }
00356     
00357     //Check if initial z is too close to boundary
00358     if (initial_pos_t[2] > 1700) initial_pos_t[2] =  1700;
00359     else if (initial_pos_t[2] < -1700) initial_pos_t[2] = -1700;
00360     
00361     initial_energy_t = recTrigger.energy();
00362     if (initial_energy_t<0) initial_energy_t = total_pe_t / 163;
00363 
00364   }
00365   else{
00366     SimpleStart();
00367     initial_energy_t = total_pe_t / 163;
00368   }
00369 
00370   Double_t minPE;
00371   if (initial_energy_t > 2.5) minPE = 0.5;
00372   else minPE = initial_energy_t/5;
00373   
00374   Int_t nGoodHits_t = 0;
00375   for (Int_t iring = 0; iring < 8; iring++){
00376     for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00377       if (pmt_flag_t[iring][icolumn] == 1){
00378 
00379         //Eliminate low charge and nohit PMT
00380         if (pmt_charge_t[iring][icolumn] < minPE) 
00381           pmt_flag_t[iring][icolumn] = -1;
00382         
00383         //Set PMT timing error for good PMT     
00384         else{
00385           nGoodHits_t++;
00386           if (pmt_charge_t[iring][icolumn]<41.08)
00387             pmt_error_t[iring][icolumn] =  (0.43 + 31.5 * pow(pmt_charge_t[iring][icolumn]*18.5,-0.67)) / sqrt(5);
00388           else
00389             pmt_error_t[iring][icolumn] =  0.80 / sqrt(5);  
00390         }
00391       
00392       }
00393     }
00394   }
00395  
00396   if(nGoodHits_t<3){
00397     recTrigger.setPositionStatus( ReconStatus::kNotConverged ); //Needs at least 3 good PMT
00398     return StatusCode::SUCCESS;
00399   }  
00400   
00401   //Calculate deltime and select PMTs
00402   //*********************************
00403   CalculateDeltime(initial_pos_t[0],initial_pos_t[1],initial_pos_t[2]);
00404   SelectPMT();
00405   
00406   Int_t iSelect = 0;
00407   Bool_t SuccessFlag = 0;
00408   Int_t nFitPMT_t = 0;
00409 
00410   while (SuccessFlag == 0 && iSelect<maxselect_t){
00411     
00412     if(nPMT_t[iSelect] > 0){
00413       nFitPMT_t += nPMT_t[iSelect];
00414       
00415       if(nFitPMT_t > 2) //Do fit only if there are at least 3 PMTs
00416         IterativeFit(nFitPMT_t,maxstep_t,initial_pos_t[0],initial_pos_t[1],initial_pos_t[2], iSelect);
00417     }
00418 
00419     if (sqrt(pow(final_pos_t[0],2)+pow(final_pos_t[1],2)) < 2300 && final_pos_t[2] < 2200 && final_pos_t[2] > -2200) SuccessFlag = 1; //Good convergence
00420     else iSelect++;
00421     
00422   }
00423   
00424   if(SuccessFlag == 0){ //Check if convergence occurred outsite physical boundary
00425     if (sqrt(pow(final_pos_t[0],2)+pow(final_pos_t[1],2)) < 10*r_phys_max_t && final_pos_t[2] < 10*z_phys_max_t && final_pos_t[2] > 10*z_phys_min_t)
00426       SuccessFlag = 1; //Convergenced
00427   }
00428 
00429   if(SuccessFlag == 0){  //Did not converge
00430     recTrigger.setPositionStatus( ReconStatus::kNotConverged );
00431     return StatusCode::SUCCESS;
00432   }
00433  
00434   CLHEP::HepLorentzVector position(final_pos_t[0], final_pos_t[1], final_pos_t[2], final_pos_t[3]);
00435   recTrigger.setPosition( position );
00436   recTrigger.setPositionQuality( chi2_sum_t(final_pos_t[0],final_pos_t[1],final_pos_t[2],final_pos_t[3],iSelect) );
00437   recTrigger.setPositionStatus( ReconStatus::kGood );
00438   
00439   //Debugging
00440   //*********
00441 
00442   if (Debug_t){
00443     if ((readout.triggerType() & DayaBay::Trigger::kCalib)==268435472){
00444       hsteps->Fill(totalsteps);
00445       hstepsVdDis->Fill(totalsteps,sqrt((final_pos_t[0]-initial_pos_t[0])*(final_pos_t[0]-initial_pos_t[0])+(final_pos_t[1]-initial_pos_t[1])*(final_pos_t[1]-initial_pos_t[1])+(final_pos_t[2]-initial_pos_t[2])*(final_pos_t[2]-initial_pos_t[2])));
00446       hstepsVchi2->Fill(totalsteps,recTrigger.positionQuality());
00447     }
00448 
00449     if ((readout.triggerType() & DayaBay::Trigger::kCalib)==268435472){ // && recTrigger.positionQuality()>0.3){ // && recTrigger.positionStatus() > 2){
00450       hchantime1->Reset();
00451       hchantime0->Reset();
00452       hchancharge1->Reset();
00453       hchancharge0->Reset();
00454       
00455       hdifftime0->Reset();
00456       
00457       for (Int_t iring = 0; iring < 8; iring++){
00458         for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00459           Int_t ichan = (icolumn+1)+ 24*iring; 
00460           
00461           if(pmt_flag_t[iring][icolumn] == 1){
00462             hchantime1->SetBinContent(ichan,pmt_time_t[iring][icolumn] - median_deltime_t);  
00463             hchantime1->SetBinError(ichan,pmt_error_t[iring][icolumn]);
00464             hchancharge1->SetBinContent(ichan,pmt_charge_t[iring][icolumn]);  
00465           }
00466           else if (pmt_flag_t[iring][icolumn] == 0){
00467             hchantime0->SetBinContent(ichan,pmt_time_t[iring][icolumn] - median_deltime_t);  
00468             hchantime0->SetBinError(ichan,pmt_error_t[iring][icolumn]);
00469             hchancharge0->SetBinContent(ichan,pmt_charge_t[iring][icolumn]); 
00470           }
00471           
00472           if(pmt_charge_t[iring][icolumn]>0.2) hdifftime0->Fill(deltime_t[iring][icolumn] - median_deltime_t);
00473         }
00474       }
00475       
00476       fittime1->FixParameter(0,final_pos_t[3] - median_deltime_t); 
00477       fittime1->FixParameter(1,final_pos_t[0]);
00478       fittime1->FixParameter(2,final_pos_t[1]);  
00479       fittime1->FixParameter(3,final_pos_t[2]);
00480       
00481       fittime1->SetLineWidth(1);
00482       fittime1->SetLineColor(2);
00483       hchantime1->Fit(fittime1,"Q");
00484       
00485       fittime0->FixParameter(0,final_pos_t[3] - median_deltime_t ); 
00486       fittime0->FixParameter(1,initial_pos_t[0]);
00487       fittime0->FixParameter(2,initial_pos_t[1]);
00488       fittime0->FixParameter(3,initial_pos_t[2]);
00489       
00490       fittime0->SetLineWidth(1);
00491       fittime0->SetLineColor(4);
00492       fittime0->SetLineStyle(9);
00493       hchantime1->Fit(fittime0,"+Q");
00494       
00495       reftime0->FixParameter(0,early_deltime_t - median_deltime_t);
00496       reftime0->SetLineWidth(1);
00497       reftime0->SetLineColor(2);
00498       reftime0->SetLineStyle(2);
00499       hchantime1->Fit(reftime0,"+Q");
00500       
00501       hchantime1->SetTitle(DayaBay::Trigger::AsString(readout.triggerType()));
00502       
00503       Char_t hname[1024];
00504       sprintf(hname, "%s_%i_hchantime1", DayaBay::Trigger::AsString(readout.triggerType()), entry);
00505       hchantime1->SetName(hname);
00506       ctime->SetName(hname);
00507       
00508       Char_t hname2[1024];
00509       sprintf(hname2, "%s_%i_hchancharge1", DayaBay::Trigger::AsString(readout.triggerType()), entry);
00510       hchancharge1->SetName(hname2);
00511       ccharge->SetName(hname2);
00512       
00513       Char_t hname3[1024];
00514       sprintf(hname3, "%s_%i_hdifftime0", DayaBay::Trigger::AsString(readout.triggerType()), entry);
00515       hdifftime0->SetName(hname3);
00516       
00517       Char_t ytitle[1024];
00518       sprintf(ytitle, "Initial t=%.2f,x=%.2f,y=%.2f,z=%.2f ", 0.0, initial_pos_t[0] ,initial_pos_t[1] ,  initial_pos_t[2] );
00519       hchantime1->GetYaxis()->SetTitle(ytitle);
00520       hchancharge1->GetYaxis()->SetTitle(ytitle);
00521       
00522       Char_t xtitle[1024];
00523       sprintf(xtitle, "(t=%.2f,x=%.2f,y=%.2f,z=%.2f)", final_pos_t[3] - median_deltime_t, final_pos_t[0], final_pos_t[1], final_pos_t[2]);
00524       hchantime1->GetXaxis()->SetTitle(xtitle);
00525       hchancharge1->GetXaxis()->SetTitle(xtitle);
00526       
00527       Char_t early_xtitle[1024];
00528       sprintf(early_xtitle, "(early_deltime_t=%.2f)", early_deltime_t - median_deltime_t);
00529       hdifftime0->GetXaxis()->SetTitle(early_xtitle);
00530       
00531       TDirectory * dir = gDirectory;
00532       
00533       savefile->cd();
00534       
00535       ctime->cd();
00536       hchantime1->Draw();
00537       hchantime0->Draw("same");
00538       ctime->Write();
00539       
00540       ccharge->cd();
00541       hchancharge1->Draw();
00542       hchancharge0->Draw("same");
00543       //htemplatecharge1->Draw("same");
00544       //htemplatecharge0->Draw("same");
00545       ccharge->Write();
00546       
00547       hdifftime0->Write();
00548       
00549       dir->cd();
00550     }
00551   }
00552   
00553   return StatusCode::SUCCESS;
00554 }
00555 
00556 void CalculateDeltime(Double_t x, Double_t y, Double_t z){
00557 
00558   //Calculate deltime = PMT Observed time - Predicted Expected time and store into a vector for chronological sort   
00559   std::vector<double> sorted_deltime;
00560   for (Int_t iring = 0; iring < 8; iring++){
00561     for (Int_t icolumn = 0; icolumn < 24; icolumn++){ 
00562       if (pmt_flag_t[iring][icolumn] != -1){
00563 
00564         deltime_t[iring][icolumn] = pmt_time_t[iring][icolumn] - (1.465/300)*sqrt(pow((x-pmt_x_t[iring][icolumn]),2)+pow((y-pmt_y_t[iring][icolumn]),2)+pow((z-pmt_z_t[iring][icolumn]),2));
00565         sorted_deltime.push_back(deltime_t[iring][icolumn]);
00566       }
00567     }
00568   }
00569 
00570   std::sort(sorted_deltime.begin(),sorted_deltime.end());
00571   Double_t median_t = sorted_deltime[int(sorted_deltime.size()/2)];  
00572   Double_t early_t = 1e10;
00573   
00574   //Find early_deltime by finding a 3ns window containing (nPMT/8) number of PMT
00575   for (Int_t i = 0; i<int(sorted_deltime.size()/2); i++){
00576     if (sorted_deltime[i+int(sorted_deltime.size()/8)] - sorted_deltime[i] < 3.0){
00577       early_t =  sorted_deltime[i];
00578       break;
00579     }
00580   }  
00581  
00582   //If first attempt to find early_deltime_t failed ...
00583   if (early_t > 9e9){
00584     //Find early_deltime by finding a 4ns window containing 5 PMT
00585     for (Int_t i = 0; i<int(sorted_deltime.size()/2); i++){
00586       if (sorted_deltime[i+4] - sorted_deltime[i] < 4.0){
00587         early_t =  sorted_deltime[i];
00588         break;
00589       }
00590     }
00591   }
00592 
00593   //If attempts yield an absurd early_t
00594   if (median_t - early_t > 10 || median_t - early_t < 0){
00595     early_t = median_t - 10.0;
00596   }
00597  
00598   early_deltime_t = early_t;
00599   median_deltime_t = median_t;
00600   
00601 }
00602  
00603  
00604 //Selecting PMTs with cuts
00605 void SelectPMT(){
00606   
00607   for (Int_t i=0; i<maxselect_t; i++){
00608     nPMT_t[i] = 0;
00609   } 
00610 
00611   Double_t min_window_t = median_deltime_t - early_deltime_t;
00612   Double_t max_window_t = 2.25 - 1.45e-3*sqrt(initial_pos_t[0]*initial_pos_t[0] + initial_pos_t[1]*initial_pos_t[1])*sqrt(initial_energy_t/2.2);
00613   
00614   //Ensure that at least half the PMTs are utilized in fit
00615   if (max_window_t <0.0) max_window_t = 0.0;
00616   
00617   for (Int_t iring = 0; iring < 8; iring++){
00618     for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00619       
00620       if (pmt_flag_t[iring][icolumn] != -1){
00621          
00622         pmt_flag_t[iring][icolumn] = -1;
00623         
00624         //Too Early
00625         if (deltime_t[iring][icolumn] < median_deltime_t - min_window_t) continue;
00626         
00627         for (Int_t i=0; i<maxselect_t; i++){
00628           if (deltime_t[iring][icolumn] < median_deltime_t + max_window_t + i*0.5){
00629             pmt_flag_t[iring][icolumn] = i;
00630             nPMT_t[i]++;
00631             break;
00632           }
00633         }
00634                 
00635       } //pmt flag
00636       
00637     } //icolumn
00638   } //iring
00639   
00640 }
00641  
00642  Double_t chi2_time_t(Double_t x, Double_t y, Double_t z, Double_t t, Int_t iRing, Int_t iColumn){
00643   Double_t t_exp0 = t +(1.465/300)*sqrt(pow((x-pmt_x_t[iRing][iColumn]),2)+pow((y-pmt_y_t[iRing][iColumn]),2)+pow((z-pmt_z_t[iRing][iColumn]),2));
00644 
00645   Double_t t_obs = pmt_time_t[iRing][iColumn];
00646   Double_t chi2 = pow((t_exp0 - t_obs)/pmt_error_t[iRing][iColumn],2);
00647   //std::cout << "iring: " << iRing << " iColumn: " << iColumn << " Theory: " << t_exp0 << " Observed: " << t_obs << " Error: " << pmt_error_t[iRing][iColumn] << " Charge: " << pmt_charge_t[iRing][iColumn] << " chi2: " << chi2 << std::endl;        
00648   return chi2;
00649 }
00650 
00651 
00652 Double_t chi2_sum_t(Double_t x, Double_t y, Double_t z, Double_t t, Int_t iSelect){
00653   Double_t chi2_tmp = 0;
00654   Int_t nHits = 0;
00655 
00656   for (Int_t iring = 0; iring < 8; iring++){
00657     for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00658       
00659       
00660       //This is the time portion
00661       if (pmt_flag_t[iring][icolumn] != -1 && pmt_flag_t[iring][icolumn] <= iSelect){
00662         chi2_tmp += chi2_time_t(x, y, z, t, iring, icolumn);
00663         nHits ++;
00664       }
00665         
00666     }
00667   }
00668 
00669   //std::cout << "Sum Chi2: " <<  chi2_tmp << " nHits: " << nHits << std::endl;
00670   if (nHits == 0) chi2_tmp = 1e10;
00671   else chi2_tmp = chi2_tmp / nHits;
00672 
00673   return chi2_tmp;
00674 }
00675 void SimpleStart(){
00676   //If AdSimple did not provide a good vertex, have to produce an initial guess
00677   Double_t Qrow[8] = {0}; //Total charge in a row  
00678   Double_t Qcol[24] = {0}; //Total charge in a column
00679   Double_t Qcol3[24] = {0};//Total charge in 3 adjacent columns
00680   
00681   for (Int_t iring = 0; iring < 8; iring++){
00682     for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00683       Qrow[iring] += pmt_charge_t[iring][icolumn];
00684       Qcol[icolumn] += pmt_charge_t[iring][icolumn];
00685 
00686       Qcol3[icolumn] += pmt_charge_t[iring][icolumn];
00687       Qcol3[(icolumn+1)%24] += pmt_charge_t[iring][icolumn];    
00688       Qcol3[(icolumn+23)%24] += pmt_charge_t[iring][icolumn];   
00689     }
00690   }
00691   
00692   Int_t maxring = 0, maxcolumn = 0;
00693   Double_t maxQrow = 0.0, maxQcol3 = 0.0;
00694   
00695   for (Int_t iring = 0; iring < 8; iring++){
00696     if (Qrow[iring] > maxQrow){
00697       maxQrow = Qrow[iring];
00698       maxring = iring;
00699     }
00700   }
00701   
00702   for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00703     if (Qcol3[icolumn] > maxQcol3){
00704       maxQcol3 = Qcol3[icolumn];
00705       maxcolumn = icolumn;
00706     }
00707   }
00708 
00709   Double_t phiMinus1 = atan2(pmt_y_t[maxring][(maxcolumn+23)%24],pmt_x_t[maxring][(maxcolumn+23)%24]); 
00710   Double_t phi0 = atan2(pmt_y_t[maxring][maxcolumn],pmt_x_t[maxring][maxcolumn]);
00711   Double_t phiPlus1 = atan2(pmt_y_t[maxring][(maxcolumn+1)%24],pmt_x_t[maxring][(maxcolumn+1)%24]);
00712 
00713   Double_t phi = ( Qcol[(maxcolumn+23)%24]*phiMinus1 + Qcol[maxcolumn]*phi0 + Qcol[(maxcolumn+1)%24]*phiPlus1 );
00714   phi = phi / ( Qcol[(maxcolumn+23)%24] + Qcol[maxcolumn] + Qcol[(maxcolumn+1)%24]);
00715 
00716   Double_t r = SimpleR(maxring,maxcolumn);
00717 
00718   
00719   initial_pos_t[0] = r*cos(phi);
00720   initial_pos_t[1] = r*sin(phi);
00721   initial_pos_t[2] = pmt_z_t[maxring][maxcolumn];
00722   initial_pos_t[3] = 1e10;
00723 }
00724 
00725 Double_t SimpleR(Int_t ring, Int_t column){
00726     
00727   Double_t r = 1100.0;
00728   
00729   std::vector<double> cluster_t1;
00730   Double_t median_t1, early_t1;
00731   //Search for usable PMT corresponding to highest charge row & column
00732   for (Int_t i = 0; i < 3; i++){
00733     Int_t di = int(pow(-1,i) * (i+1)/2);
00734     if (pmt_flag_t[ring][(column+24+di)%24] == 1){
00735       cluster_t1.push_back(pmt_time_t[ring][(column+24+di)%24]);
00736     }
00737   }
00738   
00739   if (cluster_t1.size()!=0){ 
00740     std::sort(cluster_t1.begin(),cluster_t1.end());
00741     median_t1 = cluster_t1[int(cluster_t1.size()/2)];
00742     
00743     //Find early_t1 such that it is consistent with the median_t1 to within 2ns
00744     early_t1 = median_t1;
00745     for (Int_t i=0;i<cluster_t1.size();i++){
00746       if (median_t1 - cluster_t1[i] < 2){
00747         early_t1 = cluster_t1[i];
00748         break;
00749       }
00750     }
00751     
00752     //Scan 3x3 cluster of PMT on the opposite side
00753     std::vector<double> cluster_t2;
00754     Double_t median_t2, early_t2;
00755     
00756     for (Int_t dring = -1; dring <= 1; dring++){
00757       if (ring+dring < 0 || ring+dring > 7) continue;
00758       for (Int_t dcolumn = -1; dcolumn <= 1; dcolumn++){
00759         
00760         if (pmt_flag_t[ring+dring][(column+12+dcolumn)%24] == 1){
00761           if (pmt_charge_t[ring+dring][(column+12+dcolumn)%24]>0.5 && fabs(pmt_time_t[ring+dring][(column+12+dcolumn)%24]-early_t1) < 25){
00762             cluster_t2.push_back(pmt_time_t[ring+dring][(column+12+dcolumn)%24]);
00763           }
00764         }
00765       }
00766     }
00767     
00768     
00769     if (cluster_t2.size()!=0){ 
00770       
00771       //Determine the median time of ths cluster
00772       std::sort(cluster_t2.begin(),cluster_t2.end());
00773       if (cluster_t2.size()%2 == 0){ //even
00774         median_t2 = (cluster_t2[int(cluster_t2.size()/2)-1] + cluster_t2[int(cluster_t2.size()/2)]) / 2;
00775       }
00776       else { //odd
00777         median_t2 = cluster_t2[int(cluster_t2.size()/2)];
00778       }
00779       
00780       //Find early_t2 such that it is consistent with the median_t2
00781       early_t2 = median_t2;
00782       for (Int_t i=0;i<cluster_t2.size();i++){
00783         if (median_t2 - cluster_t2[i] < 3){
00784           early_t2 = cluster_t2[i];
00785           break;
00786         }
00787       }
00788       
00789       r = 0.5*(early_t2 - early_t1) * (300/1.465);
00790     }
00791   }
00792 
00793   return r;
00794 }
00795 
00796  void IterativeFit(Int_t n, Int_t maxstep, Double_t x, Double_t y, Double_t z, Int_t iSelect){
00797   
00798   Double_t x0 = x;
00799   Double_t y0 = y;
00800   Double_t z0 = z;
00801   Double_t t0 = 0.0;
00802   Double_t x1 = 1e10;
00803   Double_t y1 = 1e10;
00804   Double_t z1 = 1e10;
00805   Double_t t1 = 1e10;
00806 
00807   Double_t dis[n];
00808   
00809   Double_t tempFunc;
00810   Double_t Func_t, Func_x, Func_y, Func_z;
00811   Double_t tempB, tempC;
00812   Double_t A, Bx, By, Bz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz;
00813   Double_t I00, I01, I02, I03, I11, I22, I33, I12, I13, I23;
00814   Double_t DetH;
00815 
00816   //Initializing
00817   Double_t pmtinfo[n][5]; //{pmt x, pmt y, pmt z, pmt time, 1/pmt error^2}
00818   Int_t i = 0;
00819   
00820   Double_t T1=0;
00821   Double_t T2=0;
00822   Double_t T3=0;
00823   for (Int_t iring = 0; iring < 8; iring++){
00824     for (Int_t icolumn = 0; icolumn < 24; icolumn++){
00825       
00826       if (pmt_flag_t[iring][icolumn] != -1 && pmt_flag_t[iring][icolumn] <= iSelect){
00827         pmtinfo[i][0] = pmt_x_t[iring][icolumn];
00828         pmtinfo[i][1] = pmt_y_t[iring][icolumn];
00829         pmtinfo[i][2] = pmt_z_t[iring][icolumn];
00830         pmtinfo[i][3] = pmt_time_t[iring][icolumn] - median_deltime_t;
00831         pmtinfo[i][4] = 1/(pmt_error_t[iring][icolumn]*pmt_error_t[iring][icolumn]);
00832         
00833         dis[i] = sqrt((x0-pmtinfo[i][0])*(x0-pmtinfo[i][0])+(y0-pmtinfo[i][1])*(y0-pmtinfo[i][1])+(z0-pmtinfo[i][2])*(z0-pmtinfo[i][2]));
00834 
00835         T1 += pmtinfo[i][4];
00836         T2 += pmtinfo[i][3]*pmtinfo[i][4];
00837         T3 += dis[i]*pmtinfo[i][4];
00838         
00839         i++;
00840       }
00841 
00842       if(i==n) break; //Get out of loop when done filling info
00843     }
00844     if(i==n) break; //Get out of loop when done filling info
00845   }
00846 
00847   Int_t step = 0;
00848   //If T1 is 0, meaningless to carry on
00849   if (T1==0){  
00850     t1 = 1e10;
00851     x1 = 1e10;
00852     y1 = 1e10;
00853     z1 = 1e10;
00854   }
00855   else{
00856    
00857     t0 = (T2 - (1.465/300)*T3) / T1;
00858     
00859     A=2*T1; 
00860     
00861     while (step < maxstep){
00862       
00863       Func_t=0;
00864       Func_x=0;
00865       Func_y=0;
00866       Func_z=0;
00867       
00868       Bx=0;
00869       By=0;
00870       Bz=0;
00871       Cxx=0;
00872       Cyy=0;
00873       Czz=0;
00874       Cxy=0;
00875       Cxz=0;
00876       Cyz=0;
00877       
00878       //Since this was already done for step 0
00879       if (step>0){
00880         for (Int_t m = 0; m<n; m++){
00881           dis[m] = sqrt((x0-pmtinfo[m][0])*(x0-pmtinfo[m][0])+(y0-pmtinfo[m][1])*(y0-pmtinfo[m][1])+(z0-pmtinfo[m][2])*(z0-pmtinfo[m][2]));
00882         }
00883       }
00884       
00885       for (Int_t m = 0; m<n; m++){  
00886         
00887         Func_t += 2 * (t0 + (1.465/300)*dis[m] - pmtinfo[m][3]) * pmtinfo[m][4];
00888         
00889         tempFunc = 2 * (1.465/300) * (t0 + (1.465/300)*dis[m] - pmtinfo[m][3])/dis[m] * pmtinfo[m][4];
00890         Func_x += tempFunc * (x0 - pmtinfo[m][0]);
00891         Func_y += tempFunc * (y0 - pmtinfo[m][1]);
00892         Func_z += tempFunc * (z0 - pmtinfo[m][2]);
00893         
00894         tempB = 2 * (1.465/300) * 1/dis[m] * pmtinfo[m][4];
00895         Bx += tempB * (x0 - pmtinfo[m][0]);
00896         By += tempB * (y0 - pmtinfo[m][1]);
00897         Bz += tempB * (z0 - pmtinfo[m][2]);
00898         
00899         tempC = -2 * (1.465/300) * (t0 - pmtinfo[m][3])/(dis[m]*dis[m]*dis[m]) * pmtinfo[m][4];
00900         Cxx += tempFunc + (x0 - pmtinfo[m][0])*(x0 - pmtinfo[m][0])*tempC;
00901         Cyy += tempFunc + (y0 - pmtinfo[m][1])*(y0 - pmtinfo[m][1])*tempC;
00902         Czz += tempFunc + (z0 - pmtinfo[m][2])*(z0 - pmtinfo[m][2])*tempC;
00903         
00904         Cxy += (x0 - pmtinfo[m][0])*(y0 - pmtinfo[m][1])*tempC;
00905         Cxz += (x0 - pmtinfo[m][0])*(z0 - pmtinfo[m][2])*tempC;
00906         Cyz += (y0 - pmtinfo[m][1])*(z0 - pmtinfo[m][2])*tempC;
00907         
00908       }
00909       
00910       I00 = Cxx*Cyy*Czz - Cxx*Cyz*Cyz - Cyy*Cxz*Cxz - Czz*Cxy*Cxy + 2*Cxy*Cyz*Cxz;
00911       I01 = Bx*(Cyz*Cyz - Cyy*Czz) + By*(Cxy*Czz - Cxz*Cyz) + Bz*(Cyy*Cxz - Cxy*Cyz);
00912       I02 = By*(Cxz*Cxz - Cxx*Czz) + Bz*(Cyz*Cxx - Cxy*Cxz) + Bx*(Czz*Cxy - Cyz*Cxz);
00913       I03 = Bz*(Cxy*Cxy - Cxx*Cyy) + Bx*(Cxz*Cyy - Cyz*Cxy) + By*(Cxx*Cyz - Cxz*Cxy);
00914       
00915       I11 = A*(Cyy*Czz - Cyz*Cyz) - By*By*Czz - Bz*Bz*Cyy + 2*By*Bz*Cyz;
00916       I22 = A*(Czz*Cxx - Cxz*Cxz) - Bz*Bz*Cxx - Bx*Bx*Czz + 2*Bz*Bx*Cxz;
00917       I33 = A*(Cxx*Cyy - Cxy*Cxy) - Bx*Bx*Cyy - By*By*Cxx + 2*Bx*By*Cxy;
00918       
00919       I13 = A*(Cxy*Cyz - Cyy*Cxz) + By*By*Cxz - Bx*By*Cyz + Bx*Bz*Cyy - By*Bz*Cxy;
00920       I12 = A*(Cyz*Cxz - Czz*Cxy) + Bz*Bz*Cxy - By*Bz*Cxz + By*Bx*Czz - Bz*Bx*Cyz;
00921       I23 = A*(Cxz*Cxy - Cxx*Cyz) + Bx*Bx*Cyz - Bz*Bx*Cxy + Bz*By*Cxx - Bx*By*Cxz;
00922       
00923       DetH = A*I00 + Bx*I01 + By*I02 + Bz*I03;
00924       
00925       t1 = t0 - (I00*Func_t + I01*Func_x + I02*Func_y + I03*Func_z)/DetH;
00926       x1 = x0 - (I01*Func_t + I11*Func_x + I12*Func_y + I13*Func_z)/DetH;
00927       y1 = y0 - (I02*Func_t + I12*Func_x + I22*Func_y + I23*Func_z)/DetH;
00928       z1 = z0 - (I03*Func_t + I13*Func_x + I23*Func_y + I33*Func_z)/DetH;
00929       
00930       //Convergence reached!
00931       if (fabs(x1-x0)<0.1 && fabs(y1-y0)<0.1 && fabs(z1-z0)<0.1 && fabs(t1-t0)<0.01) break;
00932       
00933       //Setting conditions for next iteration
00934       t0 = t1;
00935       x0 = x1;
00936       y0 = y1;
00937       z0 = z1;
00938       
00939       step ++;
00940     }
00941     
00942     
00943     //If for loop did not break, assume it did not converge
00944     if (step == maxstep){
00945       t1 = 1e10;
00946       x1 = 1e10;
00947       y1 = 1e10;
00948       z1 = 1e10;
00949     }
00950   }
00951   
00952   final_pos_t[0] = x1;
00953   final_pos_t[1] = y1;
00954   final_pos_t[2] = z1;
00955   final_pos_t[3] = t1 + median_deltime_t;
00956    
00957   totalsteps = step;
00958 
00959 }//End of IterativeFit
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

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