AlcapDAQ: rootana/src/TAP_generators/TemplateFitAPGenerator.cpp Source File

00001 #include "TAPGeneratorFactory.h"
00002 #include "TemplateFitAPGenerator.h"
00003 #include "TPulseIsland.h"
00004 #include "TAnalysedPulse.h"
00005 #include "TTemplateFitAnalysedPulse.h"
00006 #include "InterpolatePulse.h"
00007 #include "EventNavigator.h"
00008 #include "debug_tools.h"
00009 #include <iostream>
00010 using std::cout;
00011 using std::endl;
00012 
00013 TemplateArchive* TemplateFitAPGenerator::fTemplateArchive=NULL;
00014 TemplateArchive* TemplateFitAPGenerator::fTemplateArchive2=NULL;
00015 
00016 TemplateFitAPGenerator::TemplateFitAPGenerator(TAPGeneratorOptions* opts):
00017    TVAnalysedPulseGenerator("TemplateFitAPGenerator",opts),
00018     fInitPedestal(EventNavigator::Instance().GetSetupRecord().GetPedestal(GetChannel())),
00019     fTemplate2(NULL),
00020     fIntegralRatio(NULL),
00021     fMaxBin(fInitPedestal, EventNavigator::Instance().GetSetupRecord().GetPolarity(GetChannel()))
00022 {
00023    // get the templates from the archive
00024    if(!fTemplateArchive){
00025      fTemplateArchive=new TemplateArchive(opts->GetString("template_archive").c_str(),"READ");
00026    }
00027    if(opts->HasOption("template_archive_2")){
00028      std::string template_2_src=opts->GetString("template_archive_2");
00029      if (template_2_src!=fTemplateArchive->GetName() && !fTemplateArchive2)
00030         fTemplateArchive2=new TemplateArchive(opts->GetString("template_archive_2").c_str(),"READ");
00031      else if(!fTemplateArchive2){
00032         fTemplateArchive2=fTemplateArchive;
00033      }
00034    }
00035    fTemplate=fTemplateArchive->GetTemplate(GetChannel());
00036    if(fTemplateArchive2) fTemplate2=fTemplateArchive2->GetTemplate(GetChannel());
00037 
00038    // make the fitters
00039    fFitter = new TemplateMultiFitter(GetChannel().str());
00040    fFitter->AddTemplate(fTemplate);
00041    fFitter->Init();
00042    if(fTemplate2){
00043      fDoubleFitter = new TemplateMultiFitter(GetChannel().str());
00044      fDoubleFitter->AddTemplate(fTemplate);
00045      fDoubleFitter->AddTemplate(fTemplate2);
00046      fDoubleFitter->Init();
00047    }
00048 
00049    // prepare the integral ratio cuts
00050    if(opts->GetFlag("use_IR_cut")){
00051       fIntegralRatio=new Algorithm::IntegralRatio(
00052                       opts->GetInt("start_integral",10),
00053                       opts->GetInt("start_tail",60),
00054                       opts->GetInt("stop_integral",0),
00055                       EventNavigator::Instance().GetSetupRecord().GetPolarity(GetChannel()));
00056       fIntegralMax=opts->GetDouble("max_integral");
00057       fIntegralMin=opts->GetDouble("min_integral");
00058       fIntegralRatioMax=opts->GetDouble("max_ratio");
00059       fIntegralRatioMin=opts->GetDouble("min_ratio");
00060   }
00061 
00062   // cache the template parameters so we dont need to look them up
00063   fTemplateAmp=fTemplate->GetAmplitude();
00064   fTemplateTime=fTemplate->GetTime();
00065   fTemplatePed=fTemplate->GetPedestal();
00066 
00067   if(fTemplate2){
00068      fTemplate2Amp=fTemplate2->GetAmplitude();
00069      fTemplate2Time=fTemplate2->GetTime();
00070      fTemplate2Ped=fTemplate2->GetPedestal();
00071   }
00072 
00073   // chi2 to determine when we refit with 2 pulses
00074   fChi2MinToRefit=opts->GetDouble("min_chi2_to_refit",2e5);
00075 }
00076 
00077 TemplateFitAPGenerator::~TemplateFitAPGenerator(){
00078   delete fFitter;
00079     
00080   delete fDoubleFitter;
00081   if(fIntegralRatio) delete fIntegralRatio;
00082   if(fTemplateArchive) {
00083     delete fTemplateArchive;
00084     fTemplateArchive=NULL;
00085   }
00086   if(fTemplateArchive2) {
00087     delete fTemplateArchive2;
00088     fTemplateArchive2=NULL;
00089   }
00090 }
00091 
00092 int TemplateFitAPGenerator::ProcessPulses( 
00093       const PulseIslandList& pulseList,
00094       AnalysedPulseList& analysedList){
00095 
00096   // Loop over all the TPIs given to us
00097   TTemplateFitAnalysedPulse* tap;
00098   double integral, ratio;
00099   for (PulseIslandList::const_iterator tpi=pulseList.begin();
00100        tpi!=pulseList.end(); tpi++){
00101 
00102     if(fIntegralRatio && !PassesIntegralRatio(*tpi, integral,ratio)){
00103        continue;
00104     }
00105     double init_amp=fMaxBin(*tpi)/fTemplateAmp;
00106     double init_time=  fMaxBin.GetTime() * fTemplate->GetRefineFactor()- fTemplateTime ;
00107     fFitter->SetPedestal(fInitPedestal);
00108     fFitter->SetPulseEstimates(0, init_amp, init_time);
00109 
00110     // Make histo for TPI
00111     TH1D* hPulseToFit=InterpolatePulse(*tpi,"hPulseToFit","hPulseToFit",false,fTemplate->GetRefineFactor());
00112     
00113     int fit_status = fFitter->FitWithOneTimeFree(0, hPulseToFit);
00114 
00115     // Now that we've found the information we were looking for make a TAP to
00116     // hold it.  This method makes a TAP and sets the parent TPI info.  It needs
00117     // the index of the parent TPI in the container as an argument
00118     tap = MakeNewTAP<TTemplateFitAnalysedPulse>(tpi-pulseList.begin());
00119     tap->SetTemplate(fTemplate);
00120     tap->SetPedestal(fFitter->GetPedestal());
00121     tap->SetAmplitude(fFitter->GetAmplitude(0));
00122     tap->SetTime(-fFitter->GetTime(0));
00123     tap->SetChi2(fFitter->GetChi2());
00124     tap->SetNDoF(fFitter->GetNDoF());
00125     tap->SetFitStatus(fit_status);
00126     tap->SetIntegral(integral);
00127     tap->SetIntegralRatio(ratio);
00128 
00129     if(fTemplate2 && fFitter->GetChi2() > fChi2MinToRefit){
00130         TTemplateFitAnalysedPulse* second_tap=NULL;
00131         bool better_with_2 = RefitWithTwo(hPulseToFit, tap, second_tap);
00132         if(better_with_2){
00133            analysedList.push_back(tap);
00134            analysedList.push_back(second_tap);
00135         } else{
00136            analysedList.push_back(tap);
00137            if(second_tap) delete second_tap;
00138         }
00139     } else {
00140        analysedList.push_back(tap);
00141     }
00142 
00143    }
00144 
00145    // returning 0 tells the caller that we were successful, return non-zero if not
00146    return 0;
00147 }
00148 
00149 bool TemplateFitAPGenerator::PassesIntegralRatio(const TPulseIsland* pulse, double& integral, double& ratio)const{
00150    if(!fIntegralRatio) return true;
00151    try{
00152      fIntegralRatio->SetPedestalToMinimum(pulse);
00153      (*fIntegralRatio)(pulse);
00154    }catch(std::out_of_range& e){
00155      return false;
00156    }
00157    integral=fIntegralRatio->GetTotal();
00158    ratio=fIntegralRatio->GetRatio();
00159    if( fIntegralMax < integral || fIntegralMin > integral
00160        || fIntegralRatioMax < ratio || fIntegralRatioMin > ratio) {
00161        return false;
00162    }
00163    return true;
00164 }
00165 
00166 bool TemplateFitAPGenerator::RefitWithTwo(TH1D* tpi, TTemplateFitAnalysedPulse*& tap_one, TTemplateFitAnalysedPulse*& tap_two)const{
00167 
00168   // find the max bin after the current fitted tap
00169   int second_time=-9999999;
00170   double second_scale=-9999999;
00171   InitializeSecondPulse(tpi,tap_one,second_time,second_scale);
00172 
00173   // fit the second template around the second peak
00174   fDoubleFitter->SetPedestal(tap_one->GetPedestal());
00175   fDoubleFitter->SetPulseEstimates( 0, tap_one->GetAmplitude(), tap_one->GetTime());
00176   fDoubleFitter->SetPulseEstimates( 1, second_scale, second_time);
00177   int fit_status = fDoubleFitter->FitWithOneTimeFree(1, tpi);
00178 
00179   // fit the first pulse again
00180   fit_status = fDoubleFitter->FitWithOneTimeFree(0, tpi);
00181 
00182   // fit the second pulse again
00183   fit_status = fDoubleFitter->FitWithAllTimesFixed(tpi);
00184 
00185   if(fit_status==0){
00186     // Has the double fit improved the chi-2 per NDoF
00187     double old_chi2=tap_one->GetChi2()/tap_one->GetNDoF();
00188     double new_chi2=fDoubleFitter->GetChi2()/fDoubleFitter->GetNDoF();
00189     if( (old_chi2 - new_chi2) < 0){
00190       // refill tap_one's values
00191       tap_one->SetPedestal(fDoubleFitter->GetPedestal());
00192       tap_one->SetAmplitude(fDoubleFitter->GetAmplitude(0));
00193       tap_one->SetTime(fDoubleFitter->GetTime(0));
00194       tap_one->SetChi2(fDoubleFitter->GetChi2());
00195       tap_one->SetNDoF(fDoubleFitter->GetNDoF());
00196       tap_one->SetFitStatus(fit_status);
00197 
00198       // make the second tap
00199       tap_two = MakeNewTAP<TTemplateFitAnalysedPulse>(tap_one->GetParentID());
00200       tap_two->SetTemplate(fTemplate2);
00201       tap_two->SetPedestal(fDoubleFitter->GetPedestal());
00202       tap_two->SetAmplitude(fDoubleFitter->GetAmplitude(1));
00203       tap_two->SetTime(fDoubleFitter->GetTime(1));
00204       tap_two->SetChi2(fDoubleFitter->GetChi2());
00205       tap_two->SetNDoF(fDoubleFitter->GetNDoF());
00206       tap_two->SetFitStatus(fit_status);
00207 
00208       // link the two pulses together
00209       tap_one->SetOtherPulse(tap_two);
00210       tap_two->SetOtherPulse(tap_one);
00211       tap_two->IsPileUpPulse();
00212 
00213       return true;
00214     }
00215   } 
00216 
00217   return false;
00218 }
00219 
00220 void TemplateFitAPGenerator::InitializeSecondPulse(
00221           TH1D* tpi, const TTemplateFitAnalysedPulse* tap_one,
00222           int& second_time, double& second_scale)const{
00223     // scan from the time of the first tap and find the bin with the largest residual (within some window)
00224     int limit = tpi->GetNbinsX();
00225     int template_limit=tap_one->GetTime() - fTemplate->GetTime() + fTemplate->GetHisto()->GetNbinsX()-10;
00226     if(template_limit <  limit ) limit= template_limit;
00227     double residual, max_residual=-1e6;
00228     int max_bin=0;
00229     for( int bin=tap_one->GetTime(); bin < limit; ++bin){
00230       residual=tpi->GetBinContent(bin) - tap_one->GetBinContent(bin);
00231       if(residual>max_residual){
00232         max_residual=residual;
00233         max_bin=bin;
00234       }
00235     }
00236     second_time=max_bin - fTemplate2->GetTime();
00237     second_scale=max_residual / fTemplate2->GetAmplitude();
00238 }
00239 
00240 ALCAP_TAP_GENERATOR(TemplateFit,template_archive,use_IR_cut, min_integral, max_integral, min_ratio, max_ratio, template_archive_2,min_chi2_to_refit);