Skip to content

Instantly share code, notes, and snippets.

@gadamc

gadamc/test.C

Created August 22, 2011 15:05
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save gadamc/1162588 to your computer and use it in GitHub Desktop.
Save gadamc/1162588 to your computer and use it in GitHub Desktop.
Download ZIP
example of using the FIR filter and filling in a KAmpEvent file.
Raw
test.C
//start root.
//compile this file
// root [0] .L test.C+g
//run the program
// root [1] test();
//Of course, this won't work because you have to point KDataReader and KDataWriter
//to the appropriate files.
#ifndef __CINT__
#include "KDataReader.h"
#include "KDataWriter.h"
#include "KAmpEvent.h"
#include "KRawEvent.h"
#include "KBaselineRemoval.h"
#include "KPatternRemoval.h"
#include "KFIRFilter.h"
#include "KRawMuonVetoSysRecord.h"
#include "KRawBolometerRecord.h"
#include "KAmpBolometerRecord.h"
#include "KRawSambaRecord.h"
#include "KRawBoloPulseRecord.h"
#include "KAmpBoloPulseRecord.h"
#include "KPulseAnalysisRecord.h"
#endif
void test(){
KDataReader f("/Users/adam/analysis/edelweiss/data/kdata/run15/raw/le31b007_007.root");
KDataWriter ff("/Users/adam/analysis/edelweiss/data/kdata/run15/raw/le31b007_007_amp.root","KAmpEvent");
KAmpEvent *ee = (KAmpEvent *)ff.GetEvent();
KRawEvent *e = (KRawEvent *)f.GetEvent();
KBaselineRemoval bas;
KPatternRemoval pat;
KFIRFilter fir;
double coef500[] = {-7.96080145e-04, -1.67277533e-03, -3.96947296e-03, -6.81156559e-03, -9.11271994e-03, 9.89212136e-01, -9.11271994e-03, -6.81156559e-03, -3.96947296e-03, -1.67277533e-03, -7.96080145e-04};
double coef2500[] = {-0.00358845, -0.00782326, -0.0190963 , -0.03342973, -0.04525845, 0.94661957, -0.04525845, -0.03342973, -0.0190963 , -0.00782326, -0.00358845};
double* pulseIon1 = new double[8196];
double* pulseIon2 = new double[8196];
double* pulseHeat1 = new double[512];
double* pulseHeat2 = new double[512];
double *p1,*p2;
unsigned int psize;
for(int i = 0; i < f.GetEntries(); i++){
f.GetEntry(i);
ee->Clear("C");
if(i % 100 == 0) cout << "entry " << i << endl;
KRawMuonVetoSysRecord *muonRaw = (KRawMuonVetoSysRecord *)e->GetMuonVetoSystemRecord();
KRawMuonVetoSysRecord *muonAmp = (KRawMuonVetoSysRecord *)ee->GetMuonVetoSystemRecord();
*muonAmp = *muonRaw;
for(int j = 0; j < e->GetNumBolos(); j++){
if(i % 100 == 0) cout << " bolo " << j << endl;
KRawBolometerRecord *boloRaw = (KRawBolometerRecord *)e->GetBolo(j);
KAmpBolometerRecord *boloAmp = ee->AddBolo();
boloAmp->SetMass(boloRaw->GetMass());
boloAmp->SetDetectorName(boloRaw->GetDetectorName());
KRawSambaRecord *samRaw = (KRawSambaRecord *)boloRaw->GetSambaRecord();
KRawSambaRecord *samAmp = ee->AddSamba();
*samAmp = *samRaw;
boloAmp->SetSambaRecord(samAmp);
for(int k = 0; k < boloRaw->GetNumPulseRecords(); k++){
if(i % 100 == 0) cout << " pulse " << k;
KRawBoloPulseRecord *pRaw = (KRawBoloPulseRecord *)boloRaw->GetPulseRecord(k);
if(pRaw->GetPulseLength() == 0)
continue; //apparently, sometimes there are events where there is no pulse data!
KAmpBoloPulseRecord *pAmp = ee->AddBoloPulse();
boloAmp->AddPulseRecord(pAmp);
pAmp->SetBolometerRecord(boloAmp);
pAmp->SetChannelName(pRaw->GetChannelName());
pAmp->SetPulseTimeWidth(pRaw->GetPulseTimeWidth());
pAmp->SetPretriggerSize(pRaw->GetPretriggerSize());
pAmp->SetFilterSize(pRaw->GetFilterSize());
pAmp->SetHeatPulseStampWidth(pRaw->GetHeatPulseStampWidth());
pAmp->SetCryoPosition(pRaw->GetCryoPosition());
pAmp->SetPolarFet(pRaw->GetPolarFet());
pAmp->SetCorrPied(pRaw->GetCorrPied());
pAmp->SetCompModul(pRaw->GetCompModul());
pAmp->SetCorrTrngl(pRaw->GetCorrTrngl());
pAmp->SetAmplModul(pRaw->GetAmplModul());
pAmp->SetIsHeatPulse(pRaw->GetIsHeatPulse());
KPulseAnalysisRecord *sambarec = ee->AddPulseAnalysisRecord();
sambarec->SetAmp(pRaw->GetAmplitude());
sambarec->SetCalculationType("samba");
sambarec->SetIsBaseline(false);
sambarec->SetBaselineAmplitudeWidth(pRaw->GetAmplitudeBaselineNoise());
sambarec->SetUnit(0);
sambarec->SetBolometerRecord(boloAmp);
sambarec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(sambarec);
sambarec = ee->AddPulseAnalysisRecord();
sambarec->SetAmp(pRaw->GetAmplitudeBaseline());
sambarec->SetCalculationType("samba");
sambarec->SetIsBaseline(true);
sambarec->SetBaselineAmplitudeWidth(pRaw->GetAmplitudeBaselineNoise());
sambarec->SetUnit(0);
sambarec->SetBolometerRecord(boloAmp);
sambarec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(sambarec);
//process the pulse
if(pRaw->GetIsHeatPulse()){
p1= pulseHeat1; p2= pulseHeat2;
psize = 512;
}
else{
p1= pulseIon1; p2= pulseIon2;
psize = 8196;
}
bas.SetInputPulseSize(psize);
bas.SetOutputPulseSize(psize);
pat.SetInputPulseSize(psize);
pat.SetOutputPulseSize(psize);
fir.SetInputPulseSize(psize);
fir.SetOutputPulseSize(psize);
pRaw->GetTrace(p1);
bas.SetInputPulse(p1);
bas.SetOutputPulse(p2);
if(!bas.RunProcess()){
cout << "error removing basline" << endl;
continue;
}
if(!pRaw->GetIsHeatPulse()){
pat.SetPatternLength(pRaw->GetHeatPulseStampWidth());
pat.SetInputPulse(bas.GetOutputPulse());
pat.SetOutputPulse(bas.GetInputPulse());
if(!pat.RunProcess()){
cout << "error removing pattern" << endl;
continue;
}
fir.SetInputPulse(pat.GetOutputPulse());
fir.SetOutputPulse(pat.GetInputPulse());
}
else {
fir.SetInputPulse(bas.GetOutputPulse());
fir.SetOutputPulse(bas.GetInputPulse());
}
//first filter
fir.SetCoefficients(coef500, (unsigned int)11);
if(!fir.RunProcess()) {
continue;
}
int start, end;
if(pRaw->GetIsHeatPulse()){
start = 245;
end = 300;
}
else{
start = 6500;
end = 7000;
}
double maxValue = fir.GetOutputPulse()[start];
double peak = start;
for(int b = start+1; b < end; b++){
if(*(fir.GetOutputPulse()+b) < maxValue){
maxValue = *(fir.GetOutputPulse()+b);
peak = b;
}
}
KPulseAnalysisRecord *rec = ee->AddPulseAnalysisRecord();
rec->SetAmp(maxValue);
rec->SetCalculationType("hp11tapfir500hz");
rec->SetIsBaseline(false);
rec->SetPeakPosition(peak);
rec->SetUnit(0);
rec->SetBolometerRecord(boloAmp);
rec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(rec);
//get the amplitude of the baseline
peak = fir.GetCoefficientsSize()+10;
maxValue = *(fir.GetOutputPulse()+(int)peak);
rec = ee->AddPulseAnalysisRecord();
rec->SetAmp(maxValue);
rec->SetCalculationType("hp11tapfir500hz");
rec->SetIsBaseline(true);
rec->SetPeakPosition(peak);
rec->SetUnit(0);
rec->SetBolometerRecord(boloAmp);
rec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(rec);
//second filter
if(!pRaw->GetIsHeatPulse()) {
fir.SetInputPulse(pat.GetOutputPulse());
fir.SetOutputPulse(pat.GetInputPulse());
}
else {
fir.SetInputPulse(bas.GetOutputPulse());
fir.SetOutputPulse(bas.GetInputPulse());
}
fir.SetCoefficients(coef2500, 11);
if(!fir.RunProcess()) {
cout << "fir 2. error runing fir filter2" << endl;
continue;
}
if(pRaw->GetIsHeatPulse()){
start = 245;
end = 300;
}
else{
start = 6500;
end = 7000;
}
maxValue = fir.GetOutputPulse()[start];
peak = start;
for(int b = start+1; b < end; b++){
if(*(fir.GetOutputPulse()+b) < maxValue){
maxValue = *(fir.GetOutputPulse()+b);
peak = b;
}
}
rec = ee->AddPulseAnalysisRecord();
rec->SetAmp(maxValue);
rec->SetCalculationType("hp11tapfir2500hz");
rec->SetIsBaseline(false);
rec->SetPeakPosition(peak);
rec->SetUnit(0);
rec->SetBolometerRecord(boloAmp);
rec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(rec);
//get the amplitude of the baseline
peak = fir.GetCoefficientsSize()+10;
maxValue = *(fir.GetOutputPulse()+(int)peak);
rec = ee->AddPulseAnalysisRecord();
rec->SetAmp(maxValue);
rec->SetCalculationType("hp11tapfir500hz");
rec->SetIsBaseline(true);
rec->SetPeakPosition(peak);
rec->SetUnit(0);
rec->SetBolometerRecord(boloAmp);
rec->SetBoloPulseRecord(pAmp);
pAmp->AddPulseAnalysisRecord(rec);
}//for pulse loop
if(i % 100 == 0) cout << endl;
}//for bolo loop
ff.Fill();
}//for entry loop
ff.Write();
ff.Close();
f.Close();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment