diff --git a/adaptiveNB.session.py b/adaptiveNB.session.py new file mode 100644 index 0000000..d513b76 --- /dev/null +++ b/adaptiveNB.session.py @@ -0,0 +1,405 @@ +# coding: utf-8 + + +from Pyspectr.pydamm import * +import pandas as pd +from sklearn import naive_bayes as nbca + + +e=Experiment('SignalNoise.his') +e1=Experiment('Unclassified.his') +labels = ['avgBLpre','sigBLpre','avgBLpost','sigBLpost','traceMaxTime','traceMin', +'traceMinTime','factorPSDx10','traceRise','traceFall','traceReturnTime','signalNoise'] +UDF = pd.read_pickle('UnclassifiedDataAll.pkl') +classesu = UDF["signalNoise"] +for i in range(0,65536): + if abs( (UDF.iloc[i])['avgBLpre']-194.3 ) < 194.3/2: + (UDF.iloc[i])['signalNoise']=1 + +classesu = UDF["signalNoise"] +gnb = nbca.GaussianNB() + +gnb.fit( +UDF[[labels[0]]].values, +classesu +) + +gnb.theta_ +gnb.sigma_ + +for lab in labels[1:-1]: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for i in range(0,65536): + if (UDF.iloc[i])['signalNoise']==1 and abs( (UDF.iloc[i])['traceRise']-1.1e06 ) < 90*3: + (UDF.iloc[i])['signalNoise']=2 + else if (UDF.iloc[i])['signalNoise']==0 and abs( (UDF.iloc[i])['traceRise']-361 ) < 12*3: + (UDF.iloc[i])['signalNoise']=3 +for i in range(0,65536): + if (UDF.iloc[i])['signalNoise']==1 and abs( (UDF.iloc[i])['traceRise']-1.1e06 ) < 90*3: + (UDF.iloc[i])['signalNoise']=2 + elif (UDF.iloc[i])['signalNoise']==0 and abs( (UDF.iloc[i])['traceRise']-361 ) < 12*3: + (UDF.iloc[i])['signalNoise']=3 + +used_lab = [] +labels[8] +used_lab.append(labels[1]) +used_lab.append(labels[2]) +used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) + used_lab.append(labels[7]) +used_lab.append(labels[1]) +used_lab.append(labels[2]) +used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +used_lab +for lab in labels[1:-1]: + classesu = UDF["signalNoise"] + gnb.fit( UDF[used_lab], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + + +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for i in range(0,65536): + if (UDF.iloc[i])['signalNoise']==1 and abs( (UDF.iloc[i])['traceRise']-1.1e06 ) < 90: + (UDF.iloc[i])['signalNoise']=2 + elif (UDF.iloc[i])['signalNoise']==0 and abs( (UDF.iloc[i])['traceRise']-361 ) < 12: + (UDF.iloc[i])['signalNoise']=3 + +any(UDF['signalNoise']==3) +any(UDF['signalNoise']==2) +for i in range(0,65536): + if (UDF.iloc[i])['signalNoise']==1 and abs( (UDF.iloc[i])['traceRise']-1.1e06 ) < 90: + (UDF.iloc[i])['signalNoise']=2 + elif (UDF.iloc[i])['signalNoise']==0 and abs( (UDF.iloc[i])['traceRise']-361 ) < 12*5: + (UDF.iloc[i])['signalNoise']=3 + +any(UDF['signalNoise']==3) +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +any(UDF['signalNoise']==2) +any(UDF['signalNoise']==1) +any(UDF['signalNoise']==0) +for i in range(0,65536): + if abs( (UDF.iloc[i])['sigBLpost']-97 ) < 14.5*3: + (UDF.iloc[i])['signalNoise']=4 + elif abs( (UDF.iloc[i])['traceRise']-361 ) < 12*5: + (UDF.iloc[i])['signalNoise']=3 + elif abs( (UDF.iloc[i])['sigBLpost']-1.2e5 ) < 170*3: + (UDF.iloc[i])['signalNoise']=5 + + +used_lab = [] +#used_lab.append(labels[1]) +used_lab.append(labels[2]) +used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +#used_lab.append(labels[1]) +used_lab.append(labels[2]) +#used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +used_lab=[] +#used_lab.append(labels[1]) +used_lab.append(labels[2]) +#used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +used_lab=[] +used_lab.append(labels[1]) +used_lab.append(labels[2]) +#used_lab.append(labels[3]) +used_lab.append(labels[4]) +used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +used_lab=[] +#used_lab.append(labels[1]) +used_lab.append(labels[2]) +#used_lab.append(labels[3]) +used_lab.append(labels[4]) +#used_lab.append(labels[5]) +used_lab.append(labels[6]) +used_lab.append(labels[7]) +used_lab.append(labels[9]) +used_lab.append(labels[10]) +for i in range(0,65536): + if abs( (UDF.iloc[i])['sigBLpre']-27 ) < 79*3: + (UDF.iloc[i])['signalNoise']=6 + elif abs( (UDF.iloc[i])['traceMin']-194 ) < 194/2: + (UDF.iloc[i])['signalNoise']=7 + elif abs( (UDF.iloc[i])['traceMin']-199 ) < 199/2: + (UDF.iloc[i])['signalNoise']=8 + + + +for lab in used_lab: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for lab in labels: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for lab in labels[:-1] + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +for lab in labels[:-1]: + classesu = UDF["signalNoise"] + gnb.fit( UDF[[lab]], classesu) + print(gnb.theta_,gnb.sigma_,lab) + +gnb.fit( +UDF[labels[:-1]].values, +classesu +) +gnb.sigma_ +gnb.theta_ +get_ipython().run_line_magic('save', '"adaptiveNB.session" "0-50"') +for i in range(0,200): + setRegion=[0,0,0,0,0] + if (UDF.iloc[i])['signalNoise']==1 +for i in range(0,200): + setRegion=[0,0,0,0,0] + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +setRegion=[0,0,0,0,0] +plt.cla() +e.gy(911,(14,14)) +for i in range(0,200): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +for i in range(200,1000): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +for i in range(1000,10000): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +for i in range(0,3000): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +setRegion=[0,0,0,0,0] +for i in range(0,3000): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +plt.cla() +e.d(-5,-4,-3,-2,-1) +setRegion=[0,0,0,0,0] +for i in range(2100,6000): + if (UDF.iloc[i])['signalNoise']==1 and setRegion[0]==0: + e1.gy(900,(i,i),clear=False) + print(1,UDF.iloc[i]) + setRegion[0]+=1 + if (UDF.iloc[i])['signalNoise']==2 and setRegion[1]==0: + e1.gy(900,(i,i),clear=False) + print(2,UDF.iloc[i]) + setRegion[1]+=1 + if (UDF.iloc[i])['signalNoise']==4 and setRegion[2]==0: + e1.gy(900,(i,i),clear=False) + print(4,UDF.iloc[i]) + setRegion[2]+=1 + if (UDF.iloc[i])['signalNoise']==6 and setRegion[3]==0: + e1.gy(900,(i,i),clear=False) + print(6,UDF.iloc[i]) + setRegion[3]+=1 + if (UDF.iloc[i])['signalNoise']==7 and setRegion[4]==0: + e1.gy(900,(i,i),clear=False) + print(7,UDF.iloc[i]) + setRegion[4]+=1 + + + + +plt.subplots(2,5) +plt.subplots(5,2) +plt.subplot(5,2,0) +plt.subplot(5,2,1) +e.d(-10,clear=False) +e.d(-10) +plt.subplot(5,2,1) +e.d(-10,clear=False) +gnb.class_count_ +gnb.class_prior__ +gnb.class_prior_ diff --git a/pixie_ldf_c_slim b/pixie_ldf_c_slim index ad3f04e..3c609db 100755 Binary files a/pixie_ldf_c_slim and b/pixie_ldf_c_slim differ diff --git a/source/#map.txt# b/source/#map.txt# deleted file mode 100644 index dc87e0c..0000000 --- a/source/#map.txt# +++ /dev/null @@ -1,50 +0,0 @@ -%mod ch damm type subtype loc trace -0 0 616 mtas C1F 1 0 -0 1 617 mtas C1B 2 0 -0 2 618 mtas C2F 3 0 -0 3 619 mtas C2B 4 0 -0 4 620 mtas C3F 5 0 -0 5 621 mtas C3B 6 0 -0 6 622 mtas C4F 7 0 -0 7 623 mtas C4B 8 0 -0 8 624 mtas C5F 9 0 -0 9 625 mtas C5B 10 0 -0 10 626 mtas C6F 11 0 -0 11 627 mtas C6B 12 0 -0 12 628 mtas I1F 13 0 -0 13 629 mtas I1B 14 0 -0 14 630 mtas I2F 15 0 -0 15 631 mtas I2B 16 0 -1 0 632 mtas I3F 17 0 -1 1 633 mtas I3B 18 0 -1 2 634 mtas I4F 19 0 -1 3 635 mtas I4B 20 0 -1 4 636 mtas I5F 21 0 -1 5 637 mtas I5B 22 0 -1 6 638 mtas I6F 23 0 -1 7 639 mtas I6B 24 0 -1 8 640 mtas M1F 25 0 -1 9 641 mtas M1B 26 0 -1 10 642 mtas M2F 27 0 -1 11 643 mtas M2B 28 0 -1 12 644 mtas M3F 29 0 -1 13 645 mtas M3B 30 0 -1 14 646 mtas M4F 31 0 -1 15 647 mtas M4B 32 0 -2 0 648 mtas M5F 33 0 -2 1 649 mtas M5B 34 0 -2 2 650 mtas M6F 35 0 -2 3 651 mtas M6B 36 0 -2 4 651 mtas O1F 37 0 -2 5 653 mtas O1B 38 0 -2 6 654 mtas O2F 39 0 -2 7 655 mtas O2B 40 0 -2 8 656 mtas O3F 41 0 -2 9 657 mtas O3B 42 0 -2 10 658 mtas O4F 43 0 -2 11 659 mtas O4B 44 0 -2 12 660 mtas O5F 45 0 -2 13 661 mtas O5B 46 0 -2 14 662 mtas O6F 47 0 -2 15 663 mtas O6B 48 1 -3 15 664 ge single 0 0 diff --git a/source/Makefile b/source/Makefile index ed3936c..113744f 100644 --- a/source/Makefile +++ b/source/Makefile @@ -103,7 +103,6 @@ MCPPROCESSORO = McpProcessor.$(ObjSuf) MTCPROCESSORO = MtcProcessor.$(ObjSuf) DSSDPROCESSORO = DssdProcessor.$(ObjSuf) SSDPROCESSORO = SsdProcessor.$(ObjSuf) -MTASMUONPROCESSORO = MtasMuonProcessor.$(ObjSuf) MTASPROCESSORO = MtasProcessor.$(ObjSuf) TRACESUBO = TraceAnalyzer.$(ObjSuf) DETECTORDRIVERO = DetectorDriver.$(ObjSuf) @@ -112,7 +111,6 @@ RAWEVENTO = RawEvent.$(ObjSuf) ROOTPROCESSORO = RootProcessor.$(ObjSuf) RANDOMPOOLO = RandomPool.$(ObjSuf) STATSDATAO = StatsData.$(ObjSuf) -WAVEFORMPROCESSORO = WaveformProcessor.$(ObjSuf) PULSERPROCESSORO = PulserProcessor.$(ObjSuf) #THERMPROCESSORO = ThermometerProcessor.$(ObjSuf) @@ -137,8 +135,7 @@ OBJS = $(READBUFFDATAO) $(SET2CCO) $(DSSDSUBO) $(DETECTORDRIVERO) \ $(HISTOGRAMMERO) $(EVENTPROCESSORO) $(SCINTPROCESSORO) \ $(GEPROCESSORO) $(SPLINEFITPROCESSORO) $(SPLINEPROCESSORO) \ $(DSSDPROCESSORO) $(SSDPROCESSORO) $(RAWEVENTO) $(RANDOMPOOLO) \ - $(MTASPROCESSORO) $(MTASMUONPROCESSORO) $(STATSDATAO) \ - $(WAVEFORMPROCESSORO) $(PULSERPROCESSORO) + $(MTASPROCESSORO) $(STATSDATAO) $(PULSERPROCESSORO) #$(THERMPROCESSORO) #$(VANDLEPROCESSORO) $(PULSERPROCESSORO) diff --git a/source/PixieStd.o b/source/PixieStd.o index df988b5..552072c 100644 Binary files a/source/PixieStd.o and b/source/PixieStd.o differ diff --git a/source/ReadBuffData.RevD.o b/source/ReadBuffData.RevD.o index 2dd6a48..d94f4dd 100644 Binary files a/source/ReadBuffData.RevD.o and b/source/ReadBuffData.RevD.o differ diff --git a/source/WaveformProcessor.o b/source/WaveformProcessor.o deleted file mode 100644 index f7623ab..0000000 Binary files a/source/WaveformProcessor.o and /dev/null differ diff --git a/source/include/MtasMuonProcessor.h b/source/include/MtasMuonProcessor.h deleted file mode 100644 index 032acd5..0000000 --- a/source/include/MtasMuonProcessor.h +++ /dev/null @@ -1,78 +0,0 @@ -/** \file MtasMuonProcessor.h - * - * Header file for analysis of muons in Mtas - */ - -#ifndef __MTAS_MUON_PROCESSOR_H_ -#define __MTAS_MUON_PROCESSOR_H_ - -#include "EventProcessor.h" -#include - -class DetectorSummary; -class RawEvent; -class ChanEvent; - - -//extern std::vector emptyTrace; - -class MtasMuonProcessor : public EventProcessor -{ - private: - DetectorSummary *mtasSummary; - static double measureOnTime; - double firstTime; - public: - MtasMuonProcessor(); - virtual void DeclarePlots(void) const; - virtual bool Process(RawEvent &event); - public: - struct MtasData{ - MtasData() { - st = " "; - E = 0.0; - calE = 0.0; - loc = -1; - sat = 0; - pile = 0; - // tr = emptyTrace; - } - - MtasData(std::string subtype, double energy,double Calenergy, double time, int location, - int saturated, int pileup) { //, std::vector trace) { - st = subtype; - E = energy; - calE = Calenergy; - t = time; - loc = location; - sat = saturated; - pile = pileup; - //tr = trace; - } - - std::string st; - double E; - double calE; - double t; - int loc; - int sat; - int pile; - //std::vector tr; - }; - std::vector > FBEventsMatch_; -/*struct MtasData - { - MtasData(ChanEvent *chan); - - std::string detSubtype; - double energy; - double calEnergy; - double time; - double location; - int saturation; - int pileup; - std::vector trace; - };*/ -}; - -#endif // __MTAS_MUON_PROCESSOR_H_ diff --git a/source/include/WaveformProcessor.h b/source/include/WaveformProcessor.h deleted file mode 100644 index 37a6391..0000000 --- a/source/include/WaveformProcessor.h +++ /dev/null @@ -1,27 +0,0 @@ -/** \file SptFitProcessor.h - * - * Class for handling position-sensitive mcp signals - */ - -#ifndef __WAVEFORMPROCESSOR_H_ -#define __WAVEFORMPROCESSOR_H_ - -#include "EventProcessor.h" - -class WaveformProcessor : public EventProcessor -{ - public: - WaveformProcessor(); // no virtual c'tors - virtual void DeclarePlots(void) const; - virtual bool Process(RawEvent &event); - - struct FitData { - size_t n; - double * y; - double * sigma; - }; - - private: - int counter, TrcCtr, counter_1; -}; -#endif // __WAVEFORMPROCESSOR_H_ diff --git a/source/include/damm_plotids.h b/source/include/damm_plotids.h index 5339d29..6d981fb 100644 --- a/source/include/damm_plotids.h +++ b/source/include/damm_plotids.h @@ -121,10 +121,6 @@ namespace dammIds { const int DD_TDIFFVSCORCORTOF = 2470; const int DD_TQDCAVEVSCORTOF = 2510; }//vandle namespace - namespace waveformprocessor{ - const int D_DISCRIM = 5100; - const int DD_NGVSE = 5101; - } // waveformprocessor namespace namespace pulserprocessor{ //in PulserProcessor.cpp const int D_TIMEDIFF = 5000; const int D_PROBLEMSTUFF = 5001; @@ -176,8 +172,8 @@ namespace dammIds { const int MTAS_EVO_NOLOGIC = 3600; const int MTAS_LIGHTPULSER_EVO = 3700; const int MTAS_REFERENCE_EVO = 3800; - //For MtasMuonProcessor 1D (4000-4299) - const int MTAS_TOTAL = 4000; + //For MtasMuonProcessor 1D (4000-4299) //keeping incase we want to rename things. + /*const int MTAS_TOTAL = 4000; const int MTAS_CENTRAL = 4001; const int MTAS_INNER = 4002; const int MTAS_MIDDLE = 4003; @@ -232,7 +228,7 @@ namespace dammIds { const int MTAS_PREV_TMULT_V_E = 4311; const int MTAS_CURR_IMULT_V_E = 4312; const int MTAS_CURR_MOMULT_V_E = 4313; - const int MTAS_CURR_IDE_V_E = 4314; + const int MTAS_CURR_IDE_V_E = 4314;*/ } // mtas namespace diff --git a/source/src/DeclareHistogram.cpp b/source/src/DeclareHistogram.cpp index 835d948..56428b2 100644 --- a/source/src/DeclareHistogram.cpp +++ b/source/src/DeclareHistogram.cpp @@ -75,7 +75,6 @@ extern "C" void drrsub_(unsigned int& iexist) // unfortunately this function is called before we have read in the map, // and know the EXACT number of channels we'll need const int numberChannels = 96; - const int stripsDSSD = 40; extern DetectorDriver driver; @@ -91,11 +90,6 @@ extern "C" void drrsub_(unsigned int& iexist) DeclareHistogram1D(offsets::D_CAL_ENERGY+i,SE,"Channel Cal. Energy"); } - // not needed for MTAS, Jan '12, DTM - //for (int i=0; i < stripsDSSD; i++) { - //DeclareHistogram1D(offsets::D_FRONT_STRIP + i, SE, "dssd front"); - //DeclareHistogram1D(offsets::D_BACK_STRIP + i, SE, "dssd back"); - //} DeclareHistogram1D(D_HIT_SPECTRUM, S7, "channel hit spectrum"); DeclareHistogram1D(D_SUBEVENT_GAP, SE, "time btwn chan-in event,10ns bin"); diff --git a/source/src/DetectorDriver.cpp b/source/src/DetectorDriver.cpp index 8055215..d5e6513 100644 --- a/source/src/DetectorDriver.cpp +++ b/source/src/DetectorDriver.cpp @@ -17,7 +17,7 @@ SVP - Oct. '10 Added the VandleProcessor for use with VANDLE. Added the PulserProcessor for use with Pulsers. - Added the WaveformProcessor to determine ps time resolutions. + */ /*! @@ -47,12 +47,10 @@ #include "DssdProcessor.h" #include "SsdProcessor.h" #include "MtasProcessor.h" -#include "MtasMuonProcessor.h" #include "GeProcessor.h" #include "McpProcessor.h" #include "MtcProcessor.h" #include "ScintProcessor.h" -#include "WaveformProcessor.h" //#include "VandleProcessor.h" //#include "PulserProcessor.h" @@ -78,7 +76,6 @@ extern RandomPool randoms; */ DetectorDriver::DetectorDriver() { -// vecProcess.push_back(new WaveformProcessor()); // vecProcess.push_back(new ScintProcessor()); // vecProcess.push_back(new GeProcessor()); // vecProcess.push_back(new McpProcessor()); diff --git a/source/src/DetectorDriver_working.cpp b/source/src/DetectorDriver_working.cpp deleted file mode 100644 index d44e2e0..0000000 --- a/source/src/DetectorDriver_working.cpp +++ /dev/null @@ -1,562 +0,0 @@ -/* the detector driver class. - - The main analysis program. A complete event is create in PixieStd - passed into this class. See manual for further details. - - SNL - 7-2-07 - SNL - 7-12-07 - Add root analysis. If the ROOT program has been - detected on the computer system the and the - makefile has the useroot flag declared, ROOT - analysis will be included. - DTM - Oct. '09 - Significant structural/cosmetic changes. Event processing is - now primarily handled by individual event processors which - handle their own DetectorDrivers - - SVP - Oct. '10 - Added the VandleProcessor for use with VANDLE. - Added the PulserProcessor for use with Pulsers. - Added the WaveformProcessor to determine ps time resolutions. -*/ - -/*! - \file DetectorDriver.cpp - - \brief event processing - - In this file are the details for experimental processing of a raw event - created by ScanList() in PixieStd.cpp. Event processing includes things - which do not change from experiment to experiment (such as energy - calibration and raw parameter plotting) and things that do (differences - between MTC and RMS experiment, for example). -*/ - -#include -#include -#include -#include -#include - -#include "DetectorDriver.h" -#include "RandomPool.h" -#include "RawEvent.h" - -#include "damm_plotids.h" - -#include "DssdProcessor.h" -#include "SsdProcessor.h" -#include "MtasProcessor.h" -#include "GeProcessor.h" -#include "McpProcessor.h" -#include "MtcProcessor.h" -#include "ScintProcessor.h" -#include "WaveformProcessor.h" -//#include "VandleProcessor.h" -//#include "PulserProcessor.h" - -#ifdef useroot -#include "RootProcessor.h" -#endif - -using namespace std; - -/* rawevent declared in PixieStd.cpp - * - * driver relies on the respective DetectorSummary addresses remaining constant - */ -extern RawEvent rawev; - -// pool of random numbers declared in RandomPool.cpp -extern RandomPool randoms; - -/*! - detector driver constructor - - Creates instances of all event processors -*/ -DetectorDriver::DetectorDriver() -{ -// vecProcess.push_back(new WaveformProcessor()); -// vecProcess.push_back(new ScintProcessor()); -// vecProcess.push_back(new GeProcessor()); -// vecProcess.push_back(new McpProcessor()); -// vecProcess.push_back(new SsdProcessor()); -// static int i; -// if(i++%100==1) - vecProcess.push_back(new MtasProcessor()); - -// vecProcess.push_back(new DssdProcessor()); -// vecProcess.push_back(new MtcProcessor()); -// vecProcess.push_back(new PulserProcessor()); -// vecProcess.push_back(new VandleProcessor()); - - -#ifdef useroot - // and finally the root processor - vecProcess.push_back(new RootProcessor("tree.root", "tree")); -#endif -} - -/*! - detector driver deconstructor - - frees memory for all event processors - */ -DetectorDriver::~DetectorDriver() -{ - for (vector::iterator it = vecProcess.begin(); - it != vecProcess.end(); it++) { - delete *it; - } - - vecProcess.clear(); -} - -/*! - Retrieves a vector containing all detector types for which an analysis - routine has been defined making it possible to declare this detector type - in the map.txt file. The currently known detector types are in detectorStrings -*/ -const set& DetectorDriver::GetKnownDetectors() -{ - const unsigned int detTypes = 17; - const string detectorStrings[detTypes] = { - "dssd_front", "dssd_back", "idssd_front", "position", "timeclass", - "ge", "si", "scint", "mcp","mtc", "mtas","ssd","generic", "vandle", "pulser", "sili", "logi"}; - - // only call this once - if (!knownDetectors.empty()) - return knownDetectors; - - // this is a list of the detectors that are known to this program. - cout << "constructing the list of known detectors " << endl; - - //? get these from event processors - for (unsigned int i=0; i < detTypes; i++) - knownDetectors.insert(detectorStrings[i]); - - return knownDetectors; -} - -/*! - Called from PixieStd.cpp during initialization. - The calibration file cal.txt is read using the function ReadCal() and - checked to make sure that all channels have a calibration. -*/ - -int DetectorDriver::Init(void) -{ - // initialize the trace analysis routine - traceSub.Init(); - - // initialize processors in the event processing vector - for (vector::iterator it = vecProcess.begin(); - it != vecProcess.end(); it++) { - (*it)->Init(*this); - } - - /* - Read in the calibration parameters from the file cal.txt - */ - //cout << "read in the calibration parameters" << endl; - ReadCal(); - - return 0; -} - -/*! - \brief controls event processing - - The ProcessEvent() function is called from ScanList() in PixieStd.cpp - after an event has been constructed. This function is passed the mode - the analysis is currently in (the options are either "scan" or - "standaloneroot"). The function checks the thresholds for the individual - channels in the event and calibrates their energies. - The raw and calibrated energies are plotted if the appropriate DAMM spectra - have been created. Then experiment specific processing is performed. - Currently, both RMS and MTC processing is available. After all processing - has occured, appropriate plotting routines are called. -*/ -int DetectorDriver::ProcessEvent(const string &mode){ - /* - Begin the event processing looping over all the channels - that fired in this particular event. - */ - plot(dammIds::misc::D_NUMBER_OF_EVENTS, GENERIC_CHANNEL); - - const vector &eventList = rawev.GetEventList(); - for(size_t i=0; i < eventList.size(); i++) { - ChanEvent *chan = eventList[i]; - - PlotRaw(chan); - ThreshAndCal(chan); // check threshold and calibrate - PlotCal(chan); - } //end chan by chan event processing - - // have each processor in the event processing vector handle the event - for (vector::iterator iProc = vecProcess.begin(); - iProc != vecProcess.end(); iProc++) { - if ( (*iProc)->HasEvent() ) { - (*iProc)->Process(rawev); - } - } - - return 0; -} - -const set& DetectorDriver::GetUsedDetectors(void) const -{ - return rawev.GetUsedDetectors(); -} - -// declare plots for all the event processors -void DetectorDriver::DeclarePlots(void) const -{ - for (vector::const_iterator it = vecProcess.begin(); - it != vecProcess.end(); it++) { - (*it)->DeclarePlots(); - } - traceSub.DeclarePlots(); -} - -// sanity check for all our expectations -bool DetectorDriver::SanityCheck(void) const -{ - return true; -} - -/*! - \brief check threshold and calibrate each channel. - - Check the thresholds and calibrate the energy for each channel using the - calibrations contained in the calibration vector filled during ReadCal() -*/ - -int DetectorDriver::ThreshAndCal(ChanEvent *chan) -{ - // retrieve information about the channel - Identifier chanId = chan->GetChanID(); - int id = chan->GetID(); - string type = chanId.GetType(); - string subtype = chanId.GetSubtype(); - - double energy; - - if (type == "ignore" || type == "") { - return 0; - } - /* - If the channel has a trace get it, analyze it and set the energy. - */ - if ( !chan->GetTraceRef().empty() ) { - plot(dammIds::misc::D_HAS_TRACE,id); - vector values; - - traceSub.Analyze(chan->GetTraceRef(), type, subtype); - energy = traceSub.GetEnergy(); - chan->SetEnergy(energy); - } else { - // otherwise, use the Pixie on-board calculated energy - // add a random number to convert an integer value to a - // uniformly distributed floating point - - energy = chan->GetEnergy() + randoms.Get(); - energy /= ChanEvent::pixieEnergyContraction; - } - /* - Set the calibrated energy for this channel - */ - chan->SetCalEnergy( cal[id].Calibrate(energy) ); - - /* - update the detector summary - */ - rawev.GetSummary(type)->AddEvent(chan); - - return 1; -} - -/*! - Plot the raw energies of each channel into the damm spectrum number assigned - to it in the map file with an offset as defined in damm_plotids.h -*/ -int DetectorDriver::PlotRaw(const ChanEvent *chan) const -{ - int id = chan->GetID(); - float energy = chan->GetEnergy(); - - - plot(dammIds::misc::offsets::D_RAW_ENERGY + id, energy); - - return 0; -} - -/*! - Plot the calibrated energies of each channel into the damm spectrum number - assigned to it in the map file with an offset as defined in damm_plotids.h -*/ -int DetectorDriver::PlotCal(const ChanEvent *chan) const -{ - int id = chan->GetID(); - // int dammid = chan->GetChanID().GetDammID(); - float calEnergy = chan->GetCalEnergy(); - - plot(dammIds::misc::offsets::D_CAL_ENERGY + id, calEnergy); - // plot(dammid, calEnergy); - - return 0; -} - -vector DetectorDriver::GetProcessors(const string& type) const -{ - vector retVec; - - for (vector::const_iterator it = vecProcess.begin(); - it != vecProcess.end(); it++) { - if ( (*it)->GetTypes().count(type) > 0 ) - retVec.push_back(*it); - } - - return retVec; -} - -/*! - Read in the calibration for each channel according to the data in cal.txt -*/ -void DetectorDriver::ReadCal() -{ - /* - The file cal.txt contains the calibration for each channel. The first - five variables describe the detector's physical location (strip number, - detector number, ...), the detector type, the detector subtype, the number - of calibrations, and their polynomial order. Using this information, the - rest of a channel's calibration is read in as -- lower threshold for the - current calibration, followed by the polynomial constants in increasing - polynomial order. The lower thresholds and polynomial values are read in - for each distinct calibration specified by the number of calibrations. - */ - - // lookup table for information from map.txt (from PixieStd.cpp) - extern vector modChan; - Identifier lookupID; - - /* - Values used to read in the thresholds and polynomials from cal.txt - The numbers can not be read directly into the vectors - */ - float thresh; - float val; - - /* - The channels module number, channel number, detector location - the number of calibrations, polynomial order, if the detector - should be ignored, the detector type and subtype. - */ - unsigned int detLocation; - string detType, detSubtype; - - const string calFilename("cal.txt"); - - ifstream calFile(calFilename.c_str()); - - // make sure there is a generic calibration for each channel in the map - cal.resize(modChan.size()); - - if (!calFile) { - cout << "Can not open file " << calFilename << endl; - } else { - cout << "reading in calibrations from " << calFilename << endl; - while (calFile) { - /* - While the end of the calibration file has not been reached, - increment the number of lines read and if the first input on a - line is a number, read in the first five parameters for a channel - */ - if ( isdigit(calFile.peek()) ) { - calFile >> detLocation >> detType >> detSubtype; - lookupID.SetLocation(detLocation); - lookupID.SetType(detType); - lookupID.SetSubtype(detSubtype); - // find the identifier in the map - vector::iterator mapIt = - find(modChan.begin(), modChan.end(), lookupID); - if (mapIt == modChan.end()) { - cout << "Can not match detector type " << detType - << " and subtype " << detSubtype - << " with location " << detLocation - << " to a channel in the map." << endl; - exit(EXIT_FAILURE); - } - size_t id = distance(modChan.begin(), mapIt); - Calibration &detCal = cal.at(id); - //? make this a member function of Calibration - detCal.id = id; - calFile >> detCal.numCal >> detCal.polyOrder; - detCal.thresh.clear(); - detCal.val.clear(); - detCal.detLocation = detLocation; - detCal.detType = detType; - detCal.detSubtype = detSubtype; - /* - For the number of calibrations read in the - thresholds and the polynomial values - */ - for (unsigned int i = 0; i < detCal.numCal; i++) { - calFile >> thresh; - detCal.thresh.push_back(thresh); - - for(unsigned int j = 0; j < detCal.polyOrder+1; j++){ - /* - For the calibration order, read in the polynomial - constants in ascending order - */ - calFile >> val; - detCal.val.push_back(val); - } // finish looping on polynomial order - } // finish looping on number of calibrations - - /* - Add the value of MAX_PAR from the param.h file - as the upper limit of all calibrations - */ - detCal.thresh.push_back(MAX_PAR); - - /* - Check the detector type that was read in from cal.txt against - the list of detectors that can be used in the analysis (from - known detectors). If a detector is not found, exit the program - and print a warning. - */ - if (knownDetectors.find(detType) == knownDetectors.end()) { - // This is redundant while this is explicitly matched to the - // map which has identical conditions - cout << endl; - cout << "The detector called '" << detType <<"'"<< endl; - cout << "read in from the file " << calFilename << endl; - cout << "is unknown to this program!. This is a" << endl; - cout << "fatal error. Program execution halted." << endl; - cout << "Please check the " << calFilename - << " file for errors" << endl; - cout << "The currently known detectors include: " << endl; - copy(knownDetectors.begin(), knownDetectors.end(), - ostream_iterator(cout, " ")); - exit(EXIT_FAILURE); - } - } else { - // this is a comment, skip line - calFile.ignore(1000,'\n'); - } - } // end while (!calFile) loop - end reading cal.txt file - } - calFile.close(); - - // check to make sure every channel has a calibration, no default - // calibration is allowed - vector::const_iterator mapIt = modChan.begin(); - vector::const_iterator calIt = cal.begin(); - for (;mapIt != modChan.end(); mapIt++, calIt++) { - string type = mapIt->GetType(); - if (type != "ignore" && type != "" && calIt->detType!= type) { - cout << "Uncalibrated detector found for type " << type - << " at location " << mapIt->GetLocation() - << ". No default calibration is given, please correct." - << endl; - exit(EXIT_FAILURE); - } - } - /* - Print the calibration values that have been read in - */ - //cout << "calibration parameters are: " << cal.size() << endl; - - cout << setw(4) << "mod" - << setw(4) << "ch" - << setw(4) << "loc" - << setw(10) << "type" - << setw(8) << "subtype" - << setw(5) << "cals" - << setw(6) << "order" - << setw(31) << "cal values:thresh, low - high " << endl; - - //? calibration print command? - for(size_t a = 0; a < cal.size(); a++){ - cout << setw(4) << int(a/16) - << setw(4) << (a % 16) - << setw(4) << cal[a].detLocation - << setw(10) << cal[a].detType - << setw(8) << cal[a].detSubtype - << setw(5) << cal[a].numCal - << setw(6) << cal[a].polyOrder; - for(unsigned int b = 0; b < cal[a].numCal; b++){ - cout << setw(6) << cal[a].thresh[b]; - for(unsigned int c = 0; c < cal[a].polyOrder+1; c++){ - cout << setw(7) << setprecision(5) - << cal[a].val[b*(cal[a].polyOrder+1)+c]; - } - cout << setw(6) << cal[a].thresh[b+1]; - } - - cout << endl; - } -} - -/*! - Construct calibration parameters using Zero() method -*/ -Calibration::Calibration() : - id(-1), detType(""), detSubtype(""), detLocation(-1), - numCal(1), polyOrder(1) -{ - thresh.push_back(0); - thresh.push_back(MAX_PAR); - // simple linear calibration - val.push_back(0); // constant coeff - val.push_back(1); // coeff linear in raw energy -} - -double Calibration::Calibrate(double raw) -{ - /* - Make sure we don't have any calibration values below the lowest - calibration theshold or any calibrated energies above the - maximum threshold value set in cal.txt - */ - if(raw < thresh[0]) { - return 0; - } - - if(raw >= thresh[numCal]) { - return thresh[numCal] - 1; - } - - double calVal = 0; - /* - Begin threshold check and calibration, first - loop over the number of calibrations - */ - for(unsigned int a = 0; a < numCal; a++) { - //check to see if energy falls in this calibration range - if (raw >= thresh[a] && raw < thresh[a+1]) { - //loop over the polynomial order - for(unsigned int b = 0; b < polyOrder+1; b++) { - calVal += pow(raw,b) * val[a*(polyOrder+1) + b]; - } - break; - } - } - - return calVal; -} - -/*! - This function is called from the scan program - when scan is either killed or ended. If - ROOT has been enabled, close the ROOT files. - If ROOT is not enabled do nothing. -*/ -extern "C" void detectorend_() -{ - //cout << "ending, no rootfile " << endl; -} - diff --git a/source/src/MtasMuonProcessor.cpp b/source/src/MtasMuonProcessor.cpp deleted file mode 100644 index 6d2a524..0000000 --- a/source/src/MtasMuonProcessor.cpp +++ /dev/null @@ -1,361 +0,0 @@ -/*! \file MtasMuonProcessor.cpp - * - * The MTAS processor handles detectors of type mtas */ - -#include "damm_plotids.h" -#include "param.h" -#include "MtasMuonProcessor.h" -#include "DetectorDriver.h" -#include "RawEvent.h" -#include -#include -#include -#include -#include -#include - - -using std::cout; -using std::endl; -using std::vector; -using std::string; -using std::pair; -MtasMuonProcessor::MtasMuonProcessor() : - EventProcessor(), mtasSummary(NULL) -{ - firstTime = -1.0; - name = "mtas"; - associatedTypes.insert("mtas"); -} - -void MtasMuonProcessor::DeclarePlots(void) const -{ - using namespace dammIds::mtas;//MTAS_MUON - - const int EnergyBins = SE; - //QR:SE=16384,SB=2048,S8=256,S6=64,S3=8 - - //TAS spectras DeclareHistogram1D(, , " "); - DeclareHistogram1D(MTAS_TOTAL,EnergyBins, "Mtas Total "); - DeclareHistogram1D(MTAS_CENTRAL,EnergyBins, "Mtas Central "); - DeclareHistogram1D(MTAS_INNER,EnergyBins, "Mtas Inner "); - DeclareHistogram1D(MTAS_MIDDLE,EnergyBins, "Mtas Middle "); - DeclareHistogram1D(MTAS_OUTER,EnergyBins, "Mtas Outer "); - DeclareHistogram1D(MTAS_TOTAL_CR,EnergyBins, "Mtas Total crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_CENTRAL_CR,EnergyBins, "Mtas Central crunch E/10 (~100 kev/ch) "); - DeclareHistogram1D(MTAS_INNER_CR,EnergyBins, "Mtas Inner crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_MIDDLE_CR,EnergyBins, "Mtas Middle crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_OUTER_CR,EnergyBins, "Mtas Outer crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_TOTAL_CORR,EnergyBins, "Corr. Mtas Total "); - DeclareHistogram1D(MTAS_CENTRAL_CORR,EnergyBins, "Corr. Mtas Central "); - DeclareHistogram1D(MTAS_INNER_CORR,EnergyBins, "Corr. Mtas Inner "); - DeclareHistogram1D(MTAS_MIDDLE_CORR,EnergyBins, "Corr. Mtas Middle "); - DeclareHistogram1D(MTAS_OUTER_CORR,EnergyBins, "Corr. Mtas Outer "); - DeclareHistogram1D(MTAS_TOTAL_CR_CORR,EnergyBins, "Corr. Mtas Total crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_CENTRAL_CR_CORR,EnergyBins, "Corr. Mtas Central crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_INNER_CR_CORR, EnergyBins, "Corr. Mtas Inner crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_MIDDLE_CR_CORR,EnergyBins, "Corr. Mtas Middle crunch E/10 (~100 kev/ch)"); - DeclareHistogram1D(MTAS_OUTER_CR_CORR,EnergyBins, "Corr. Mtas Outer crunch E/10 (~100 kev/ch)"); - - DeclareHistogram1D(MTAS_TOTAL_AFTER,EnergyBins, "Mtas Total After HE "); - DeclareHistogram1D(MTAS_CENTRAL_AFTER,EnergyBins, "Mtas Central After HE"); - DeclareHistogram1D(MTAS_INNER_AFTER,EnergyBins, "Mtas Inner After HE"); - DeclareHistogram1D(MTAS_MIDDLE_AFTER,EnergyBins, "Mtas Middle After HE"); - DeclareHistogram1D(MTAS_OUTER_AFTER,EnergyBins, "Mtas Outer After HE"); - DeclareHistogram1D(MTAS_TOTAL_CR_AFTER,EnergyBins, "Mtas Total After HE (~100 kev/ch) "); - DeclareHistogram1D(MTAS_CENTRAL_CR_AFTER,EnergyBins, "Mtas Central After HE (~100 kev/ch)"); - DeclareHistogram1D(MTAS_INNER_CR_AFTER, EnergyBins, "Mtas Inner After HE (~100 kev/ch)"); - DeclareHistogram1D(MTAS_MIDDLE_CR_AFTER,EnergyBins, "Mtas Middle After HE (~100 kev/ch)"); - DeclareHistogram1D(MTAS_OUTER_CR_AFTER,EnergyBins, "Mtas Outer After HE (~100 kev/ch)"); - - DeclareHistogram1D(MTAS_TOTAL_NEAT,EnergyBins, "Mtas Total Interesting "); - DeclareHistogram1D(MTAS_CENTRAL_NEAT,EnergyBins, "Mtas Central Interesting"); - DeclareHistogram1D(MTAS_INNER_NEAT,EnergyBins, "Mtas Inner Interesting"); - DeclareHistogram1D(MTAS_MIDDLE_NEAT,EnergyBins, "Mtas Middle Interesting"); - DeclareHistogram1D(MTAS_OUTER_NEAT,EnergyBins, "Mtas Outer Interesting"); - DeclareHistogram1D(MTAS_TOTAL_CR_NEAT,EnergyBins, "Mtas Total Interesting (~100 kev/ch) "); - DeclareHistogram1D(MTAS_CENTRAL_CR_NEAT,EnergyBins, "Mtas Central Interesting (~100 kev/ch)"); - DeclareHistogram1D(MTAS_INNER_CR_NEAT, EnergyBins, "Mtas Inner Interesting (~100 kev/ch)"); - DeclareHistogram1D(MTAS_MIDDLE_CR_NEAT,EnergyBins, "Mtas Middle Interesting (~100 kev/ch)"); - DeclareHistogram1D(MTAS_OUTER_CR_NEAT,EnergyBins, "Mtas Outer Interesting (~100 kev/ch)"); - - DeclareHistogram1D(MTAS_MULT,S6,"MTAS Multiplicity"); - - //2D spectras DeclareHistogram2D(, , , " "); - DeclareHistogram2D(MTAS_E_AFTER_V_T,SE,SB, "MTAS Energy v time after HE (10KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_E_AFTER_V_T_CR,SE,SB, "MTAS Energy v time after HE (100KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_MCUT_E_AFTER_V_T,SE,SB, "MTAS Mult.>2 Energy v time after HE (10KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_MCUT_E_AFTER_V_T_CR,SE,SB, "MTAS Mult.>2 Energy v time after HE (100KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_E_NEAT_V_T,SE,SB, "MTAS Energy v time Interesting (10KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_E_NEAT_V_T_CR,SE,SB, "MTAS Energy v time Interesting (100KeV/ch vs ms) "); - DeclareHistogram2D(MTAS_LOC_V_E_AFTER, EnergyBins,S6, " MTAS Energy per Module after HE"); - - DeclareHistogram2D(MTAS_TMULT_V_E,EnergyBins , S5, " MTAS Total Mult vs. Total Energy"); - DeclareHistogram2D(MTAS_PREV_TMULT_V_E,EnergyBins , S5, " MTAS Total Prev. Mult vs. Total Prev. Energy/10"); - DeclareHistogram2D(MTAS_CURR_IMULT_V_E,EnergyBins ,S3 , " MTAS Inner Mult vs. Total Energy"); - DeclareHistogram2D(MTAS_CURR_MOMULT_V_E,EnergyBins ,S4 , " MTAS MO Mult vs. Total Energy "); - DeclareHistogram2D(MTAS_CURR_IDE_V_E,EnergyBins ,SD , " MTAS dE(LR) when 'Muon' IMult=1"); - -} - -bool MtasMuonProcessor::Process(RawEvent &event) -{ - using namespace dammIds::mtas; - if (!EventProcessor::Process(event)) - return false; - - // first time through, grab the according detector summaries - if (mtasSummary == NULL) - mtasSummary = event.GetSummary("mtas"); - vector mtasList= mtasSummary->GetList(); - vector mtasFList; - vector mtasBList; - - - int mtasMult=mtasSummary->GetMult(); - plot(MTAS_MULT, mtasMult); - //double maxLocation =0; - int numCpmts = 0; - double lasttime =-1.; - static double previousHighEtime =-1.0; - static double currentHighEtime =-1.0; - static double previousHighE =-1.0; - static double currentHighE =-1.0; - vector totalMtasEnergy (5,-1); - vector corrMtasEnergy (5,-1); - vector corrMtasChanEnergy (49,-1); - static int prevTMult=0, currTMult=0; - int IsInteresting = 0; - static int prevIMult=0, currIMult=0; - static int prevMOMult=0,currMOMult=0; - double currIdE=0; - FBEventsMatch_.clear(); - //Preprocess events, make detector pairs and sums of rings. - if (currentHighE > 1000. ) - { - previousHighEtime=currentHighEtime; - previousHighE=currentHighE; - prevMOMult=currMOMult; - prevIMult=currIMult; - prevTMult=currTMult; - if ( currIMult>2 && currMOMult==1) //currMOMult<3 && currMOMult>0 - { - IsInteresting=1; - } - } else - { - previousHighEtime=0.; - previousHighE=0.; - prevMOMult=-1; - prevIMult=-1; - prevTMult=-1; - } - - currIMult=0; - currMOMult=0; - currTMult=0; - for(vector::const_iterator it = mtasList.begin(); it != mtasList.end(); it++) - { - string subtype = (*it)->GetChanID().GetSubtype(); - double energy = ((*it)-> GetCalEnergy()); - int loc= (*it)->GetChanID().GetLocation(); - int hasSat=(*it)->IsSaturated(); - const vector tr=(*it)->GetTraceRef(); - /* if (tr.size()!=0) { - cout << tr.size() <IsPileup(); - double time= (*it)->GetTime()* pixie::clockInSeconds; - if (time > lasttime ) - { - lasttime=time; - } - /*static int pilecnt=0; - if (pileup>0) { - pilecnt+=1; - if (pilecnt%100==0||energy>10) { - cout << pilecnt << " " << energy << " dt:" << time - previousHighEtime << endl; - } - }*/ - if (hasSat>0 || energy > 31500) //pileup >0 - { - continue; - } - MtasData ev(subtype, energy, (*it)-> GetCalEnergy(), time, - loc, hasSat, pileup ); - if(subtype[0] =='C') - numCpmts++; - if (energy== 0) - continue; - if (subtype[2]=='F') - mtasFList.push_back(ev); - if (subtype[2]=='B') - mtasBList.push_back(ev); - - if (subtype[0]=='C'){ - totalMtasEnergy.at(1) += energy/12.0; - totalMtasEnergy.at(0) += energy/12.0; - } - if (subtype[0]=='I') { - totalMtasEnergy.at(2) += energy/2.0; - totalMtasEnergy.at(0) += energy/2.0; - } - if (subtype[0]=='M') { - totalMtasEnergy.at(3) += energy/2.0; - totalMtasEnergy.at(0) += energy/2.0; - } - if (subtype[0]=='O') { - totalMtasEnergy.at(4) += energy/2.0; - totalMtasEnergy.at(0) += energy/2.0; - } - } - - plot(MTAS_TOTAL,totalMtasEnergy.at(0)); - plot(MTAS_CENTRAL,totalMtasEnergy.at(1)); - plot(MTAS_INNER,totalMtasEnergy.at(2)); - plot(MTAS_MIDDLE,totalMtasEnergy.at(3)); - plot(MTAS_OUTER,totalMtasEnergy.at(4)); - plot(MTAS_TOTAL_CR,totalMtasEnergy.at(0)/10); - plot(MTAS_CENTRAL_CR,totalMtasEnergy.at(1)/10); - plot(MTAS_INNER_CR,totalMtasEnergy.at(2)/10); - plot(MTAS_MIDDLE_CR,totalMtasEnergy.at(3)/10); - plot(MTAS_OUTER_CR,totalMtasEnergy.at(4)/10); - - //PreProcess: Make a map of pairs of pmts. - for(vector::const_iterator fit = mtasFList.begin(); fit != mtasFList.end(); fit++) - { - string subtypeF = (*fit).st; - string modF (subtypeF, 0 ,2 ); - - for(vector::const_iterator bit = mtasBList.begin(); bit != mtasBList.end(); bit++) - { - string subtypeB = (*bit).st; - string modB (subtypeB, 0 ,2 ); - if (modF==modB) { - corrMtasChanEnergy.at((*fit).loc)=(*fit).calE; - corrMtasChanEnergy.at((*bit).loc)=(*bit).calE; - if (subtypeB[0]=='C'){ - corrMtasEnergy.at(1) += ((*fit).calE + (*bit).calE)/numCpmts; - corrMtasEnergy.at(0) += ((*fit).calE + (*bit).calE)/numCpmts; - } - if (subtypeB[0]=='I') { - corrMtasEnergy.at(2) += ((*fit).calE + (*bit).calE)/2.0; - corrMtasEnergy.at(0) += ((*fit).calE + (*bit).calE)/2.0; - currTMult+=1; - currIMult+=1; - currIdE=((*fit).calE-(*bit).calE); - } - if (subtypeB[0]=='M') { - corrMtasEnergy.at(3) += ((*fit).calE + (*bit).calE)/2.0; - corrMtasEnergy.at(0) += ((*fit).calE + (*bit).calE)/2.0; - currMOMult+=1; - currTMult+=1; - } - if (subtypeB[0]=='O') { - corrMtasEnergy.at(4) += ((*fit).calE + (*bit).calE)/2.0; - corrMtasEnergy.at(0) += ((*fit).calE + (*bit).calE)/2.0; - currMOMult+=1; - currTMult+=1; - } - pair match((*fit),(*bit)); - FBEventsMatch_.push_back(match); - - } - - } - } - if (corrMtasEnergy.at(1)>0.) - { - currTMult+=1; - } - currentHighEtime=lasttime; - currentHighE=corrMtasEnergy.at(0); - - plot(MTAS_TOTAL_CORR,corrMtasEnergy.at(0)); - plot(MTAS_CENTRAL_CORR,corrMtasEnergy.at(1)); - plot(MTAS_INNER_CORR,corrMtasEnergy.at(2)); - plot(MTAS_MIDDLE_CORR,corrMtasEnergy.at(3)); - plot(MTAS_OUTER_CORR,corrMtasEnergy.at(4)); - plot(MTAS_TOTAL_CR_CORR,corrMtasEnergy.at(0)/10); - plot(MTAS_CENTRAL_CR_CORR,corrMtasEnergy.at(1)/10); - plot(MTAS_INNER_CR_CORR,corrMtasEnergy.at(2)/10); - plot(MTAS_MIDDLE_CR_CORR,corrMtasEnergy.at(3)/10); - plot(MTAS_OUTER_CR_CORR,corrMtasEnergy.at(4)/10); - - plot(MTAS_PREV_TMULT_V_E,previousHighE/10,prevTMult); - //Process Loop - /*for(vector< pair >:: iterator mit = FBEventsMatch_.begin(); mit != FBEventsMatch_.end(); mit++) - { - double energy=0.; - char ring = (*mit).first.st[0]; - - if (ring=='C') { - energy=((*mit).first.E+(*mit).second.E)/numCpmts; - } else { - energy=((*mit).first.E+(*mit).second.E)/2.; - } - - - - - }*/ - double dt = currentHighEtime-previousHighEtime; - if (dt>0 && dt < (20.48)*1e-6 ) - { - plot(MTAS_TOTAL_AFTER,corrMtasEnergy.at(0)); - plot(MTAS_CENTRAL_AFTER,corrMtasEnergy.at(1)); - plot(MTAS_INNER_AFTER,corrMtasEnergy.at(2)); - plot(MTAS_MIDDLE_AFTER,corrMtasEnergy.at(3)); - plot(MTAS_OUTER_AFTER,corrMtasEnergy.at(4)); - plot(MTAS_TOTAL_CR_AFTER,corrMtasEnergy.at(0)/10.); - plot(MTAS_CENTRAL_CR_AFTER,corrMtasEnergy.at(1)/10.); - plot(MTAS_INNER_CR_AFTER,corrMtasEnergy.at(2)/10.); - plot(MTAS_MIDDLE_CR_AFTER,corrMtasEnergy.at(3)/10.); - plot(MTAS_OUTER_CR_AFTER,corrMtasEnergy.at(4)/10.); - //2d plotting - plot(MTAS_E_AFTER_V_T, corrMtasEnergy.at(0), dt*1e8); - plot(MTAS_E_AFTER_V_T_CR, corrMtasEnergy.at(0)/10, dt*1e8); - if (currTMult>2) - { - plot(MTAS_MCUT_E_AFTER_V_T, corrMtasEnergy.at(0), dt*1e8); - plot(MTAS_MCUT_E_AFTER_V_T_CR, corrMtasEnergy.at(0)/10, dt*1e8); - } - /*for ( int id=1; id!=49; id++) //Calibrations - { - plot(MTAS_LOC_V_E_AFTER, corrMtasChanEnergy.at(id),id); - }*/ - if (IsInteresting>0 && currTMult>2) - { - plot(MTAS_TOTAL_NEAT,corrMtasEnergy.at(0)); - plot(MTAS_CENTRAL_NEAT,corrMtasEnergy.at(1)); - plot(MTAS_INNER_NEAT,corrMtasEnergy.at(2)); - plot(MTAS_MIDDLE_NEAT,corrMtasEnergy.at(3)); - plot(MTAS_OUTER_NEAT,corrMtasEnergy.at(4)); - - //2d plotting - plot(MTAS_E_NEAT_V_T, corrMtasEnergy.at(0), dt*1e8); - plot(MTAS_E_NEAT_V_T_CR, corrMtasEnergy.at(0)/10, dt*1e8); - - } - if (IsInteresting>0) - { - plot(MTAS_TOTAL_CR_NEAT,corrMtasEnergy.at(0)/10.); - plot(MTAS_CENTRAL_CR_NEAT,corrMtasEnergy.at(1)/10.); - plot(MTAS_INNER_CR_NEAT,corrMtasEnergy.at(2)/10.); - plot(MTAS_MIDDLE_CR_NEAT,corrMtasEnergy.at(3)/10.); - plot(MTAS_OUTER_CR_NEAT,corrMtasEnergy.at(4)/10.); - plot(MTAS_TMULT_V_E,corrMtasEnergy.at(0),currTMult); - plot(MTAS_CURR_IMULT_V_E,corrMtasEnergy.at(0),currIMult); - plot(MTAS_CURR_MOMULT_V_E,corrMtasEnergy.at(0),currMOMult); - - if (currIMult == 1) { - plot(MTAS_CURR_IDE_V_E,corrMtasEnergy.at(0),currIdE+4000); - } - } - } - - EndProcess(); // update the processing time - - return true; -} diff --git a/source/src/MtasMuonProcessor.v1.cpp b/source/src/MtasMuonProcessor.v1.cpp deleted file mode 100644 index b6bd3d8..0000000 --- a/source/src/MtasMuonProcessor.v1.cpp +++ /dev/null @@ -1,71 +0,0 @@ -/*! \file MtasMuonProcessor.cpp - * - * The MTAS processor handles detectors of type mtas */ - -#include "damm_plotids.h" -#include "param.h" -#include "MtasMuonProcessor.h" -#include "DetectorDriver.h" -#include "RawEvent.h" -#include -#include -#include -#include -#include - - -using std::cout; -using std::endl; -using std::vector; -using std::string; - -MtasMuonProcessor::MtasMuonProcessor() : - EventProcessor(), mtasSummary(NULL) -{ - firstTime = -1.0; - name = "mtas"; - associatedTypes.insert("mtas"); -} - -void MtasMuonProcessor::DeclarePlots(void) const -{ - using namespace dammIds::mtas;//MTAS_MUON - - const int EnergyBins = SE; - //QR:SE=16384,SB=2048,S8=256,S6=64,S3=8 - - //TAS spectras DeclareHistogram1D(MTAS_MUON, , " "); - DeclareHistogram1D(MTAS_MULT,S6,"MTAS Multiplicity"); - - //2D spectras DeclareHistogram2D(MTAS_MUON, , , " "); -} - -bool MtasMuonProcessor::Process(RawEvent &event) -{ - using namespace dammIds::mtas; - if (!EventProcessor::Process(event)) - return false; - - // first time through, grab the according detector summaries - if (mtasSummary == NULL) - mtasSummary = event.GetSummary("mtas"); - vector mtasList= mtasSummary->GetList(); - int mtasMult=mtasSummary->GetMult(); - - for(vector::const_iterator it = mtasList.begin(); it != mtasList.end(); it++) - { - double energy = ((*it)-> GetEnergy()); - int location = (*it)->GetChanID().GetLocation(); - string subtype = (*it)->GetChanID().GetSubtype(); - int hasSat=0; - - } - - - - EndProcess(); // update the processing time - - return true; -} - - diff --git a/source/src/PixieStd.cpp b/source/src/PixieStd.cpp index f47baf4..d027fb6 100644 --- a/source/src/PixieStd.cpp +++ b/source/src/PixieStd.cpp @@ -206,10 +206,10 @@ extern "C" void hissub_(unsigned short *sbuf[],unsigned short *nhw) spillValidCount++; bufInSpill = 0; dataWords = 0; lastBuf = -1; } else if (bufNum == 0) { -// cout << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE" -// << " INCOMPLETE BUFFER " << spillInvalidCount++ -// << "\n " << spillValidCount << " valid spills so far." -// << " Starting fresh spill." << endl; +// cout << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE" //NTB this is annoying. + cout << " INCOMPLETE BUFFER " << spillInvalidCount++ + << "\n " << spillValidCount << " valid spills so far." + << " Starting fresh spill." << endl; // throw away previous collected data and start fresh bufInSpill = 0; dataWords = 0; lastBuf = -1; } @@ -256,14 +256,14 @@ extern "C" void hissub_(unsigned short *sbuf[],unsigned short *nhw) /* make sure we retrieved all the chunks of the spill */ if (bufInSpill != totBuf) { -// cout << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE INCOMPLETE BUFFER " -// << spillInvalidCount++ -// << "\n I/B [ " << bufInSpill << " of " << totBuf << " : pos " << totWords -// << " " << spillValidCount << " total spills" -// << "\n| " << hex << buf[0] << " " << buf[1] << " " -// << buf[2] << " " << buf[3] -// << "\n| " << dec << buf[totWords] << " " << buf[totWords+1] << " " -// << buf[totWords+2] << " " << buf[totWords+3] << endl; +// cout << "EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE INCOMPLETE BUFFER " //consider verbose toggle + cout << "Incomplete Buffer " << spillInvalidCount++ + << "\n I/B [ " << bufInSpill << " of " << totBuf << " : pos " << totWords + << " " << spillValidCount << " total spills" + << "\n| " << hex << buf[0] << " " << buf[1] << " " + << buf[2] << " " << buf[3] + << "\n| " << dec << buf[totWords] << " " << buf[totWords+1] << " " + << buf[totWords+2] << " " << buf[totWords+3] << endl; } else { spillValidCount++; MakeModuleData(totData, dataWords); diff --git a/source/src/PulserProcessor.cpp b/source/src/PulserProcessor.cpp index 35c9715..8182f3d 100644 --- a/source/src/PulserProcessor.cpp +++ b/source/src/PulserProcessor.cpp @@ -90,7 +90,7 @@ void PulserProcessor::AnalyzeData(void) map::iterator itStart = pulserMap.find("start"); map::iterator itStop = pulserMap.find("stop"); - int maxPosStart = (*itStart).second.maxPos; + //int maxPosStart = (*itStart).second.maxPos; if((*itStart).second.maxPos ==41) //if(((*itStart).second.trace.at(maxPosStart+1) > 1466) && ((*itStart).second.trace.at(maxPosStart+1) < 1506)) diff --git a/source/src/ReadBuffData.RevD.cpp b/source/src/ReadBuffData.RevD.cpp index c44da42..d75b784 100644 --- a/source/src/ReadBuffData.RevD.cpp +++ b/source/src/ReadBuffData.RevD.cpp @@ -108,8 +108,6 @@ int ReadBuffData(word_t *buf, unsigned long *bufLen, //cout << "buf[0] , mod" << buf[0] << " , " << modNum << endl; - if (modNum != 5) - { if( *bufLen > 0 ) { // check if the buffer has data if (*bufLen == 2) { // this is an empty channel @@ -154,7 +152,7 @@ Remodified (NTB) to use interpret saturated and pileup bits { outfile.open("headerandsums.FA3.txt", std::ofstream::out | std::ofstream::app); outfile << " MOD:CHAN " << modNum << ":" << chanNum; - for(int j=0;j> 16)/32768; //(NTB) //cout << cfdTime << " =? " << ((buf[2] & 0x7FFF0000) >> 31) << endl; //(NTB) //adjusted by xia for 14bit allowance. so our 12bit is written too high. (NTB) - word_t cfdForcedTrig = (buf[2] & 0x80000000) >> 31; //new NTB + //word_t cfdForcedTrig = (buf[2] & 0x80000000) >> 31; //new NTB unused in out implementation word_t energy = buf[3] & 0x00007FFF;//changed NTB word_t traceFlag = (buf[3] & 0x00008000) >> 15; //new @@ -264,7 +262,7 @@ Remodified (NTB) to use interpret saturated and pileup bits eventList.push_back(currentEvt); numEvents++; - } while ( buf < bufStart + *bufLen ); + } while ( buf < bufStart + *bufLen ); //?? NTB } else {// if buffer has data cout << "ERROR BufNData " << *bufLen << endl; cout << "ERROR IN ReadBuffData" << endl; @@ -273,145 +271,4 @@ Remodified (NTB) to use interpret saturated and pileup bits } return numEvents; - } else { - if( *bufLen > 0 ) { // check if the buffer has data - if (*bufLen == 2) { - // this is an empty channel - return 0; - } - do { - ChanEvent *currentEvt = new ChanEvent; - //cout << "detsub" << currentEvt->GetChanID().GetSubtype() << endl; - // decoding event data... see pixie16app.c - // buf points to the start of channel data - word_t chanNum = (buf[0] & 0x0000000F); - word_t slotNum = (buf[0] & 0x000000F0) >> 4; - word_t crateNum = (buf[0] & 0x00000F00) >> 8; - word_t headerLength = (buf[0] & 0x0001F000) >> 12; - - word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; // Event Length now in [30:17] - word_t finishCode = (buf[0] & 0x80000000) >> 31; - - //currentEvt->virtualChannel = ((buf[0] & 0x20000000) != 0); - //currentEvt->saturatedBit = ((buf[0] & 0x40000000) != 0); - currentEvt->pileupBit = (finishCode != 0); //!!!(NTB) should be modified to match Manual ie Pileup to FinishCode throughout but this is a major job. TBD. - -/* // MODIFIED to ignore saturated bit (DTM) - // word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; - word_t eventLength = (buf[0] & 0x3FFE0000) >> 17; - word_t finishCode = (buf[0] & 0x80000000) >> 31; -Remodified (NTB) to use interpret saturated and pileup bits -*/ - // Rev. D header lengths not clearly defined in pixie16app_defs - //! magic numbers here for now - // make some sanity checks - if (headerLength == stats.headerLength) { - // this is a manual statistics block inserted by the poll program - stats.DoStatisticsBlock(&buf[1], modNum); - buf += eventLength; - numEvents = readbuff::STATS; - continue; - } - if (headerLength != 4 && headerLength != 8 && - headerLength != 12 && headerLength != 16) { // (NTB) not sure if this (8,12,16) is appropriate for this version of the pixie cards. - cout << " Unexpected header length: " << headerLength << endl; - cout << " Buffer " << modNum << " of length " << *bufLen << endl; - cout << " CHAN:SLOT:CRATE " - << chanNum << ":" << slotNum << ":" << crateNum << endl; - // advance to next event and continue - // buf += EventLength; - // continue; - - // skip the rest of this buffer (nope, NTB) - return numEvents; - //(NTB) return readbuff::ERROR; - } - - word_t lowTime = buf[1];//ok - - word_t highTime = buf[2] & 0x0000FFFF;//ok - word_t cfdTime = ((buf[2] & 0x7FFF0000) >> 31); //(NTB) - //cout << cfdTime << " =? " << ((buf[2] & 0x7FFF0000) >> 31) << endl; //(NTB) - //adjusted by xia for 14bit allowance. Changed from previous. (NTB) - word_t cfdForcedTrig = (buf[2] & 0x80000000) >> 31; //new NTB - - word_t energy = buf[3] & 0x0000FFFF;//changed NTB - word_t traceLength = (buf[3] & 0x7FFF0000) >> 16; //ok - word_t flag_alt = (buf[3] & 0x40000000) >> 30; - word_t traceFlag = (buf[3] & 0x80000000) >> 31; //new NTB - /*if ( chanNum ==15 && modNum==4) { - //if (traceFlag != 0 && chanNum ==15 && modNum==4) { - cout << " Energy 15bit = " << energy << " 16bit = " << energy_alt << " ." << endl; - cout << " traceFlag = " << traceFlag << "." << endl; - }*/ - // one last sanity check - // - if (traceFlag) { - buf += eventLength; - continue; - } - if ( traceLength / 2 + headerLength != eventLength ) { - cout << " Bad event length (" << eventLength - << ") does not correspond with length of header (" << headerLength - << ") and length of trace (" << traceLength << ")" << endl; - cout << "slot:chan:mod " << slotNum <<":"<< chanNum <<":"<< modNum << " energy " << energy << " Sat. " << traceFlag << " tLen " << traceLength << " hLen " << headerLength << " eLen " << eventLength - << " Pileup " << finishCode << " extra " << flag_alt << endl; - buf += eventLength; - continue; - } - - // handle multiple crates - modNum += 100 * crateNum; - //TBD pass new vars to currentEvt and process as necessary - currentEvt->chanNum = chanNum; - currentEvt->modNum = modNum; - currentEvt->energy = energy; - currentEvt->trigTime = lowTime; - currentEvt->cfdTime = cfdTime; - currentEvt->eventTimeHi = highTime; - currentEvt->eventTimeLo = lowTime; - currentEvt->time = highTime * HIGH_MULT + lowTime;//should also have + cfdTime?? - //Do I need to check CFDForced = 1 -> cfdTime == 0?? - //For checking syncronization - /*if (chanNum == 15) { - cout << "sync check. " << endl; - cout << "chanNum " << chanNum << " modNum " << modNum << " energy " << energy << " time " << highTime * HIGH_MULT + lowTime << endl; - }*/ - /*if (traceLength !=0) - { - cout << "slot:chan:mod " << slotNum <<":"<< chanNum <<":"<< modNum << " hLen " << headerLength << " eLen " << eventLength << " Pileup " << finishCode << " w3 " << buf[3] << endl; -// cout << "hw " << headerWords[0] << " , " << headerWords[1] << " , " << headerWords[2] << " , " << headerWords[3] << " , " << endl; - - }*/ - buf += headerLength; - /* Check if trace data follows the channel header */ - if ( traceLength > 0 ) { - // sbuf points to the beginning of trace data - halfword_t *sbuf = (halfword_t *)buf; - currentEvt->trace.reserve(traceLength); //NTB - - /*if(currentEvt->saturatedBit) - currentEvt->trace.SetTraceInfo("saturation", 1); //NTB*/ - - // Read the trace data (2-bytes per sample, i.e. 2 samples per word) - for(unsigned int k = 0; k < traceLength; k ++) { - currentEvt->trace.push_back(sbuf[k]); - } - buf += traceLength / 2; - } - eventList.push_back(currentEvt); - - numEvents++; - } while ( buf < bufStart + *bufLen ); - } else {// if buffer has data - cout << "ERROR BufNData " << *bufLen << endl; - cout << "ERROR IN ReadBuffData" << endl; - cout << "LIST UNKNOWN" << endl; - return readbuff::ERROR; - } - - return numEvents; - - } -} - +} //end ReadBuffData diff --git a/source/src/ReadBuffData.RevD2018.cpp b/source/src/ReadBuffData.RevD2018.cpp new file mode 100644 index 0000000..c44da42 --- /dev/null +++ b/source/src/ReadBuffData.RevD2018.cpp @@ -0,0 +1,417 @@ +/*! + \file ReadBuffData.RevD.fip32495.cpp + Read Buffer Data program for the following configuration + from pixie16_revd_general_release_08292015 + Jun 15 2015 fippixie16_revdgeneral_r32495.bin + Aug 19 2015 syspixie16_revdgeneral_r33157.bin + and Pixie16RevF_General_072016 pulled my JMA from Igor distr. + Jul 2017 fippixie16_revfgeneral_14b100m_r29406.bin + Jul 2017 syspixie16_revfgeneral_adc100mhz_r33338.bin + Id: + utilities.c 29501 + pixie16app.c 32533 + communication.c 27938 + 2cm24c64.c 15626 + pixie16sys.c 27118 + tools.c 28166 + + \brief retrieve data from raw buffer array ibuf +*/ + +/*---------------------------------------------------------------------- + * Copyright (c) 2005, XIA LLC + * All rights reserved. + * + * Redistribution and use in source and binary forms, + * with or without modification, are permitted provided + * that the following conditions are met: + * + * * Redistributions of source code must retain the above + * copyright notice, this list of conditions and the + * following disclaimer. + * * Redistributions in binary form must reproduce the + * above copyright notice, this list of conditions and the + * following disclaimer in the documentation and/or other + * materials provided with the distribution. + * * Neither the name of XIA LLC nor the names of its + * contributors may be used to endorse or promote + * products derived from this software without + * specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND + * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, + * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR + * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF + * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + *----------------------------------------------------------------------*/ +#include +#include +#include + +#include + +// data related to pixie packet structure +#include "pixie16app_defs.h" + +// our event structure +#include "param.h" +#include "RawEvent.h" + +using pixie::word_t; +using pixie::halfword_t; +using std::cout; +using std::endl; + +extern StatsData stats; + +// define tst bit function from pixie16 files +unsigned long TstBit(unsigned short bit, unsigned long value) +{ + return ((value & (unsigned long)(pow(2.0, (double)bit))) >> bit); +} + +/*! + \brief extract channel information from raw data + + ReadBuffData extracts channel information from the raw data arrays + and places it into a structure called evt. A pointer to each + of the evt objects is placed in the eventlist vector for later time + sorting. +*/ + +vector headerWords(4); +int ReadBuffData(word_t *buf, unsigned long *bufLen, + vector &eventList) +{ + // multiplier for high bits of 48-bit time + static const double HIGH_MULT = pow(2., 32.); + + word_t modNum; + + unsigned long numEvents = 0; + word_t *bufStart = buf; + + /* Determine the number of words in the buffer */ + *bufLen = *buf++; + + /* Read the module number */ + modNum = *buf++; + + + //cout << "buf[0] , mod" << buf[0] << " , " << modNum << endl; + + if (modNum != 5) + { + if( *bufLen > 0 ) { // check if the buffer has data + if (*bufLen == 2) { + // this is an empty channel + return 0; + } + do { + ChanEvent *currentEvt = new ChanEvent; + //cout << "detsub" << currentEvt->GetChanID().GetSubtype() << endl; + // decoding event data... see pixie16app.c + // buf points to the start of channel data + word_t chanNum = (buf[0] & 0x0000000F); + word_t slotNum = (buf[0] & 0x000000F0) >> 4; + word_t crateNum = (buf[0] & 0x00000F00) >> 8; + word_t headerLength = (buf[0] & 0x0001F000) >> 12; + + word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; // Event Length now in [30:17] + word_t finishCode = (buf[0] & 0x80000000) >> 31; + + //currentEvt->virtualChannel = ((buf[0] & 0x20000000) != 0); + //currentEvt->saturatedBit = ((buf[0] & 0x40000000) != 0); + currentEvt->pileupBit = (finishCode != 0); //!!!(NTB) should be modified to match Manual ie Pileup to FinishCode throughout but this is a major job. TBD. + +/* // MODIFIED to ignore saturated bit (DTM) + // word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; + word_t eventLength = (buf[0] & 0x3FFE0000) >> 17; + word_t finishCode = (buf[0] & 0x80000000) >> 31; +Remodified (NTB) to use interpret saturated and pileup bits +*/ + // Rev. D header lengths not clearly defined in pixie16app_defs + //! magic numbers here for now + // make some sanity checks +/* if (modNum == 4) + { + headerWords[0] = buf[0]; + headerWords[1] = buf[1]; + headerWords[2] = buf[2]; + headerWords[3] = buf[3]; + }*/ + + std::ofstream outfile; + if (headerLength >=8) + { + outfile.open("headerandsums.FA3.txt", std::ofstream::out | std::ofstream::app); + outfile << " MOD:CHAN " << modNum << ":" << chanNum; + for(int j=0;j> 16)/32768; //(NTB) + //cout << cfdTime << " =? " << ((buf[2] & 0x7FFF0000) >> 31) << endl; //(NTB) + //adjusted by xia for 14bit allowance. so our 12bit is written too high. (NTB) + word_t cfdForcedTrig = (buf[2] & 0x80000000) >> 31; //new NTB + + word_t energy = buf[3] & 0x00007FFF;//changed NTB + word_t traceFlag = (buf[3] & 0x00008000) >> 15; //new + word_t traceLength = (buf[3] & 0xFFFF0000) >> 16; //ok + //word_t energy_alt = buf[3] & 0x0000FFFF; // for comparison to software documentation + /*if ( chanNum ==15 && modNum==4) { + //if (traceFlag != 0 && chanNum ==15 && modNum==4) { + cout << " Energy 15bit = " << energy << " 16bit = " << energy_alt << " ." << endl; + cout << " traceFlag = " << traceFlag << "." << endl; + }*/ + // one last sanity check + if ( traceLength / 2 + headerLength != eventLength ) { + cout << " Bad event length (" << eventLength + << ") does not correspond with length of header (" << headerLength + << ") and length of trace (" << traceLength << ")" << endl; + buf += eventLength; + continue; + } + + // handle multiple crates + modNum += 100 * crateNum; + //TBD pass new vars to currentEvt and process as necessary + currentEvt->chanNum = chanNum; + currentEvt->modNum = modNum; + currentEvt->energy = energy; + currentEvt->trigTime = lowTime; + currentEvt->cfdTime = cfdTime; + currentEvt->eventTimeHi = highTime; + currentEvt->eventTimeLo = lowTime; + //currentEvt->time = highTime * HIGH_MULT + lowTime;//should also have + cfdTime?? + + if (modNum == 3 )//&& (chanNum == 13 || chanNum == 14 ) ) + { + currentEvt->time = highTime * HIGH_MULT + lowTime - 220;//Make sure CFD time || phase is added to time in 'Detector'Processor! + } else { + currentEvt->time = highTime * HIGH_MULT + lowTime;//Make sure CFD time || phase is added to time in 'Detector'Processor! + } + + currentEvt->saturatedBit = (traceFlag != 0); + //Do I need to check CFDForced = 1 -> cfdTime == 0?? + //For checking syncronization + /*if (chanNum == 15) { + cout << "sync check. " << endl; + cout << "chanNum " << chanNum << " modNum " << modNum << " energy " << energy << " time " << highTime * HIGH_MULT + lowTime << endl; + }*/ + /*if (traceLength !=0) + { + cout << "slot:chan:mod " << slotNum <<":"<< chanNum <<":"<< modNum << " energy " << energy << " Sat. " << traceFlag << " tLen " << traceLength << " hLen " << headerLength << " eLen " << eventLength + << " Pileup " << finishCode << endl; + cout << "t " << buf[1] << " 3: " << buf[3] << endl; + }*/ + // one last sanity check + + buf += headerLength; + /* Check if trace data follows the channel header */ + if ( traceLength > 0 ) { + // sbuf points to the beginning of trace data + halfword_t *sbuf = (halfword_t *)buf; + currentEvt->trace.reserve(traceLength); //NTB + + /*if(currentEvt->saturatedBit) + currentEvt->trace.SetTraceInfo("saturation", 1); //NTB*/ + + // Read the trace data (2-bytes per sample, i.e. 2 samples per word) + for(unsigned int k = 0; k < traceLength; k ++) { + currentEvt->trace.push_back(sbuf[k]); + } + buf += traceLength / 2; + } + eventList.push_back(currentEvt); + + numEvents++; + } while ( buf < bufStart + *bufLen ); + } else {// if buffer has data + cout << "ERROR BufNData " << *bufLen << endl; + cout << "ERROR IN ReadBuffData" << endl; + cout << "LIST UNKNOWN" << endl; + return readbuff::ERROR; + } + + return numEvents; + } else { + if( *bufLen > 0 ) { // check if the buffer has data + if (*bufLen == 2) { + // this is an empty channel + return 0; + } + do { + ChanEvent *currentEvt = new ChanEvent; + //cout << "detsub" << currentEvt->GetChanID().GetSubtype() << endl; + // decoding event data... see pixie16app.c + // buf points to the start of channel data + word_t chanNum = (buf[0] & 0x0000000F); + word_t slotNum = (buf[0] & 0x000000F0) >> 4; + word_t crateNum = (buf[0] & 0x00000F00) >> 8; + word_t headerLength = (buf[0] & 0x0001F000) >> 12; + + word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; // Event Length now in [30:17] + word_t finishCode = (buf[0] & 0x80000000) >> 31; + + //currentEvt->virtualChannel = ((buf[0] & 0x20000000) != 0); + //currentEvt->saturatedBit = ((buf[0] & 0x40000000) != 0); + currentEvt->pileupBit = (finishCode != 0); //!!!(NTB) should be modified to match Manual ie Pileup to FinishCode throughout but this is a major job. TBD. + +/* // MODIFIED to ignore saturated bit (DTM) + // word_t eventLength = (buf[0] & 0x7FFE0000) >> 17; + word_t eventLength = (buf[0] & 0x3FFE0000) >> 17; + word_t finishCode = (buf[0] & 0x80000000) >> 31; +Remodified (NTB) to use interpret saturated and pileup bits +*/ + // Rev. D header lengths not clearly defined in pixie16app_defs + //! magic numbers here for now + // make some sanity checks + if (headerLength == stats.headerLength) { + // this is a manual statistics block inserted by the poll program + stats.DoStatisticsBlock(&buf[1], modNum); + buf += eventLength; + numEvents = readbuff::STATS; + continue; + } + if (headerLength != 4 && headerLength != 8 && + headerLength != 12 && headerLength != 16) { // (NTB) not sure if this (8,12,16) is appropriate for this version of the pixie cards. + cout << " Unexpected header length: " << headerLength << endl; + cout << " Buffer " << modNum << " of length " << *bufLen << endl; + cout << " CHAN:SLOT:CRATE " + << chanNum << ":" << slotNum << ":" << crateNum << endl; + // advance to next event and continue + // buf += EventLength; + // continue; + + // skip the rest of this buffer (nope, NTB) + return numEvents; + //(NTB) return readbuff::ERROR; + } + + word_t lowTime = buf[1];//ok + + word_t highTime = buf[2] & 0x0000FFFF;//ok + word_t cfdTime = ((buf[2] & 0x7FFF0000) >> 31); //(NTB) + //cout << cfdTime << " =? " << ((buf[2] & 0x7FFF0000) >> 31) << endl; //(NTB) + //adjusted by xia for 14bit allowance. Changed from previous. (NTB) + word_t cfdForcedTrig = (buf[2] & 0x80000000) >> 31; //new NTB + + word_t energy = buf[3] & 0x0000FFFF;//changed NTB + word_t traceLength = (buf[3] & 0x7FFF0000) >> 16; //ok + word_t flag_alt = (buf[3] & 0x40000000) >> 30; + word_t traceFlag = (buf[3] & 0x80000000) >> 31; //new NTB + /*if ( chanNum ==15 && modNum==4) { + //if (traceFlag != 0 && chanNum ==15 && modNum==4) { + cout << " Energy 15bit = " << energy << " 16bit = " << energy_alt << " ." << endl; + cout << " traceFlag = " << traceFlag << "." << endl; + }*/ + // one last sanity check + // + if (traceFlag) { + buf += eventLength; + continue; + } + if ( traceLength / 2 + headerLength != eventLength ) { + cout << " Bad event length (" << eventLength + << ") does not correspond with length of header (" << headerLength + << ") and length of trace (" << traceLength << ")" << endl; + cout << "slot:chan:mod " << slotNum <<":"<< chanNum <<":"<< modNum << " energy " << energy << " Sat. " << traceFlag << " tLen " << traceLength << " hLen " << headerLength << " eLen " << eventLength + << " Pileup " << finishCode << " extra " << flag_alt << endl; + buf += eventLength; + continue; + } + + // handle multiple crates + modNum += 100 * crateNum; + //TBD pass new vars to currentEvt and process as necessary + currentEvt->chanNum = chanNum; + currentEvt->modNum = modNum; + currentEvt->energy = energy; + currentEvt->trigTime = lowTime; + currentEvt->cfdTime = cfdTime; + currentEvt->eventTimeHi = highTime; + currentEvt->eventTimeLo = lowTime; + currentEvt->time = highTime * HIGH_MULT + lowTime;//should also have + cfdTime?? + //Do I need to check CFDForced = 1 -> cfdTime == 0?? + //For checking syncronization + /*if (chanNum == 15) { + cout << "sync check. " << endl; + cout << "chanNum " << chanNum << " modNum " << modNum << " energy " << energy << " time " << highTime * HIGH_MULT + lowTime << endl; + }*/ + /*if (traceLength !=0) + { + cout << "slot:chan:mod " << slotNum <<":"<< chanNum <<":"<< modNum << " hLen " << headerLength << " eLen " << eventLength << " Pileup " << finishCode << " w3 " << buf[3] << endl; +// cout << "hw " << headerWords[0] << " , " << headerWords[1] << " , " << headerWords[2] << " , " << headerWords[3] << " , " << endl; + + }*/ + buf += headerLength; + /* Check if trace data follows the channel header */ + if ( traceLength > 0 ) { + // sbuf points to the beginning of trace data + halfword_t *sbuf = (halfword_t *)buf; + currentEvt->trace.reserve(traceLength); //NTB + + /*if(currentEvt->saturatedBit) + currentEvt->trace.SetTraceInfo("saturation", 1); //NTB*/ + + // Read the trace data (2-bytes per sample, i.e. 2 samples per word) + for(unsigned int k = 0; k < traceLength; k ++) { + currentEvt->trace.push_back(sbuf[k]); + } + buf += traceLength / 2; + } + eventList.push_back(currentEvt); + + numEvents++; + } while ( buf < bufStart + *bufLen ); + } else {// if buffer has data + cout << "ERROR BufNData " << *bufLen << endl; + cout << "ERROR IN ReadBuffData" << endl; + cout << "LIST UNKNOWN" << endl; + return readbuff::ERROR; + } + + return numEvents; + + } +} + diff --git a/source/src/SsdProcessor.cpp b/source/src/SsdProcessor.cpp index b1b3cdf..76eb41d 100644 --- a/source/src/SsdProcessor.cpp +++ b/source/src/SsdProcessor.cpp @@ -28,7 +28,7 @@ void SsdProcessor::DeclarePlots(void) const const int EnergyBins = SE; const int positionBins = S4; - const int timeBins = S8; + //const int timeBins = S8; NTB Not used DeclareHistogram2D(SSD1_POSITION_ENERGY, EnergyBins, positionBins, "SSD1 Strip vs E - RF"); @@ -51,7 +51,7 @@ bool SsdProcessor::Process(RawEvent &event) // static Correlator &corr = event.GetCorrelator(); int ssdPos = UINT_MAX; - double ssdEnergy, ssdTime = 0.; + double ssdEnergy;//, ssdTime = 0.; NTB Not used string ssdSubtype=""; bool hasSsd = ssdSummary && (ssdSummary->GetMult() > 0); @@ -68,7 +68,7 @@ bool SsdProcessor::Process(RawEvent &event) ssdPos = ch->GetChanID().GetLocation(); // cout <<" TST " <GetCalEnergy(); - ssdTime = ch->GetTime(); + //ssdTime = ch->GetTime(); NTB Not used } else ssdEnergy = 0.; // plot stuff diff --git a/source/src/TraceAnalyzer.cpp b/source/src/TraceAnalyzer.cpp index c0ba7b8..8b1ae63 100644 --- a/source/src/TraceAnalyzer.cpp +++ b/source/src/TraceAnalyzer.cpp @@ -647,8 +647,8 @@ void TraceAnalyzer::TracePlot(const vector &trace) stdDevBaselinePost, double(traceMaxTime), double(traceMin), double(traceMinTime), avgBaseline, stdDevBaseline, factorPSD, (traceRise+1)*10, (traceFall+1)*10, (traceRise/traceFall +1)*10}; - double meanSignal=0., sigmaSignal=0., ampSignal =0.; - double meanNoise=0., sigmaNoise=0., ampNoise =0.; + double meanSignal=0., sigmaSignal=0.; //, ampSignal =0.; + double meanNoise=0., sigmaNoise=0.; //, ampNoise =0.; double probSignal=.02, probClassifierSignal=1., probNoise=0.98, probClassifierNoise =1.; for (int tcIndex =0; tcIndex < 13; tcIndex++) { @@ -674,7 +674,7 @@ void TraceAnalyzer::TracePlot(const vector &trace) meanSignal /= 1000.; meanSignal -= tcValue; sigmaSignal = signalVector[tcIndex][pIndex*3+1]/1000.; - ampSignal = signalVector[tcIndex][pIndex*3+2]/1000.; + // ampSignal = signalVector[tcIndex][pIndex*3+2]/1000.; double ratioSignal = meanSignal/sigmaSignal; probClassifierSignal += exp(-(ratioSignal*ratioSignal)/2.) / @@ -691,7 +691,7 @@ void TraceAnalyzer::TracePlot(const vector &trace) meanNoise /= 1000.; meanNoise -= tcValue; sigmaNoise = noiseVector[tcIndex][pIndex*3+1]/1000.; - ampNoise = noiseVector[tcIndex][pIndex*3+2]/1000.; + //ampNoise = noiseVector[tcIndex][pIndex*3+2]/1000.; double ratioNoise = meanNoise/sigmaNoise; probClassifierNoise += exp(-(ratioNoise*ratioNoise)/2.) / diff --git a/source/src/WaveformProcessor.cpp b/source/src/WaveformProcessor.cpp deleted file mode 100644 index f569c7e..0000000 --- a/source/src/WaveformProcessor.cpp +++ /dev/null @@ -1,370 +0,0 @@ -/************************************ -This code will obtain the phase of a trace -using either a Chi^2 fitting routine in GSL, -or a Single Point Analysis. - -The data is passed into the raw event so that other -subroutines are able to access it. - - S.V.P. 16 July 2009 - -*************************************/ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -#include "damm_plotids.h" -#include "DetectorDriver.h" -#include "RawEvent.h" -#include "WaveformProcessor.h" -#include "StatsAccumulator.h" - -#ifdef pulsefit -#include -#include -#include -#include -#include -#include -#endif - -//A FEW MAGIC NUMBERS -#define WID 3.9626 //define gaussian width for fit 3.9626 -#define DKAY 3.2334 //define decay constant for fit 3.2334 -#define WAVEFORMLOW 2 //The starting position of the Waveform (referenced from max) -#define WAVEFORMHIGH 12 //The point when the Waveform returns to baseline (referenced from max) - -extern "C" void count1cc_(const int &, const int &, const int &); -extern "C" void set2cc_(const int &, const int &, const int &, const int &); - -void ngdiscrim(const vector &trace, const double &traceQDC, const double &ave_baseline, const int &maxX); - -#ifdef pulsefit -double fitroutine(const vector &trace, const int &maxX, const double &ave_baseline, const double &s_dev_baseline); -int my_f(const gsl_vector *x, void *FitData, gsl_vector *f); -int expb_df(const gsl_vector *x, void *FitData, gsl_matrix *J); -int expb_fdf(const gsl_vector *x, void *FitData, gsl_vector *f, gsl_matrix *J); -#else -double spt_analysis(const vector &trace, const int &maxX, const double &ave_baseline, const double &trcQDC, const int &counter); -#endif - -using namespace dammIds::waveformprocessor; -using namespace std; - -WaveformProcessor::WaveformProcessor(): EventProcessor() -{ - name = "FittingRoutine"; - associatedTypes.insert("scint"); - associatedTypes.insert("vandle"); - associatedTypes.insert("pulser"); - - counter = 0; -} - -void WaveformProcessor::DeclarePlots(void) const -{ -// DeclareHistogram1D(D_DISCRIM, SD, "N-Gamma Discrimination"); -// DeclareHistogram2D(DD_NGVSE, SE, SE,"N-G Discrim vs Energy"); -} - -//process trace -bool WaveformProcessor::Process(RawEvent &event) -{ - if (!EventProcessor::Process(event)) - return(false); - - static const DetectorSummary* vandleEvents = event.GetSummary("vandle"); - static const DetectorSummary* scintEvents = event.GetSummary("scint"); - static const DetectorSummary* pulserEvents = event.GetSummary("pulser"); - - if(vandleEvents && scintEvents) - if(vandleEvents->GetList().empty() || scintEvents->GetList().empty()) - return(false); - - vector allEvents; - if(vandleEvents) - allEvents.insert(allEvents.end(), vandleEvents->GetList().begin(), vandleEvents->GetList().end()); - if(scintEvents) - allEvents.insert(allEvents.end(), scintEvents->GetList().begin(), scintEvents->GetList().end()); - if(pulserEvents) - allEvents.insert(allEvents.end(), pulserEvents->GetList().begin(), pulserEvents->GetList().end()); - - for(vector::const_iterator it = allEvents.begin(); it != allEvents.end(); it++) - { - ChanEvent *chan = *it; - const unsigned int location = chan->GetChanID().GetLocation(); - const string subType = chan->GetChanID().GetSubtype(); - const vector &trace = chan->GetTraceRef(); - - //initalize the variables to be passed to RawEvent - chan->SetTrcQDC(-9999); - chan->SetPhase(-9999); - chan->SetMaxValue(-9999); - chan->SetStdDevBaseline(-9999); - chan->SetAveBaseline(-9999); - chan->SetMaxPos(-9999); - - int satTraceCount = count(trace.begin(),trace.end(),4095); - - if(trace.empty()|| (satTraceCount > 0)) //SKIP IF NO TRACE OR SATURATION - continue; - -/**** LOCATE THE WAVEFORM PEAK ****/ - vector::const_iterator itTrace = max_element(trace.begin()+WAVEFORMLOW, trace.end()-WAVEFORMHIGH); - double max_value = *itTrace; - int max_x = int(itTrace-trace.begin()); - -/**** CALCULATE THE AVERAGE AND S_DEV OF THE BASELINE (REPLACE WITH STATS ACCUMULATOR) ****/ - double sumOfBaseline = 0, baselineOffset = 0; - int numBinsBaseline = 15; - - for (unsigned int v=0; v < numBinsBaseline; v++) - { - unsigned int baselineValue = trace.at(int(v + baselineOffset)); - sumOfBaseline += baselineValue; - } - - double aveBaseline = sumOfBaseline / numBinsBaseline; - double stdDevSum = 0; - - for (unsigned int z=0; z < numBinsBaseline; z++) - { - unsigned int baselineValue = trace.at(z); - stdDevSum += (baselineValue - aveBaseline) * (baselineValue - aveBaseline); - } - double stdDevBaseline = sqrt((1/double(numBinsBaseline))*stdDevSum); - -/**** GET THE QDC VALUE OF THE WAVEFORM ****/ - double traceQDC=0; - for(unsigned int j = (max_x - WAVEFORMLOW); (j < (max_x + WAVEFORMHIGH)) && (j < trace.size()); j++) - traceQDC += (trace.at(j)-aveBaseline); - -/**** N-GAMMA DISCRIMINATION ****/ - if(subType == "liquid") - ngdiscrim(trace, traceQDC, aveBaseline, max_x); - - if((stdDevBaseline <= 3) && ((max_value-aveBaseline) >= 5)) //reject noise condition - { - chan->SetStdDevBaseline(stdDevBaseline); - chan->SetAveBaseline(aveBaseline); -#ifdef pulsefit - chan->SetPhase(fitroutine(trace, max_x, aveBaseline, stdDevBaseline)); -#else - chan->SetPhase(spt_analysis(trace, max_x, aveBaseline, traceQDC, counter)); -#endif - chan->SetTrcQDC(traceQDC); - chan->SetMaxValue(max_value-aveBaseline); - chan->SetMaxPos(max_x); - } - } //END LOOP OVER WHOLE EVENT - - EndProcess(); - return(true); -} - -void ngdiscrim(const vector &trace, const double &traceQDC, const double &aveBaseline, const int &maxX) -{ - double discrim = 0, discrim_norm = 0; - - for(size_t j = maxX+5; (j < maxX+WAVEFORMHIGH) && (j < trace.size()); j++) - discrim += (trace.at(j)-aveBaseline); - - discrim_norm = (discrim/traceQDC)*10000+150; - - plot(D_DISCRIM,int(discrim_norm),1); - plot(DD_NGVSE, int(discrim), int(traceQDC)); -} - -#ifdef pulsefit -double fitroutine(const vector &trace, const int &maxX, const double &aveBaseline, const double &stdDevBaseline) -{ - vector trace_array_fit; - - for(int ll = (maxX-WAVEFORMLOW); ll <= (maxX+WAVEFORMHIGH); ll++) - trace_array_fit.push_back(trace.at(ll) - aveBaseline); - - const int size_fit = trace_array_fit.size(); - - const gsl_multifit_fdfsolver_type *T; - gsl_multifit_fdfsolver *s; - int status; - unsigned int iter = 0; - const size_t numParams = 2; - - gsl_matrix *covar = gsl_matrix_alloc (numParams, numParams); - double y[size_fit], sigma[size_fit]; - struct WaveformProcessor::FitData d = {size_fit, y, sigma}; - gsl_multifit_function_fdf f; - double x_init[numParams] = {2, trace.at(maxX)-aveBaseline}; - gsl_vector_view x = gsl_vector_view_array (x_init, numParams); - - f.f = &my_f; - f.df = &expb_df; - f.fdf = &expb_fdf; - f.n = size_fit; - f.p = numParams; - f.params = &d; - - for(int b = 0; b < size_fit; b++) - { - y[b] = trace_array_fit.at(b); - sigma[b] = stdDevBaseline; - } - - T = gsl_multifit_fdfsolver_lmsder; - s = gsl_multifit_fdfsolver_alloc (T, size_fit, numParams); - gsl_multifit_fdfsolver_set (s, &f, &x.vector); - - do - { - iter++; - - status = gsl_multifit_fdfsolver_iterate(s); - - if (status) - break; - - status = gsl_multifit_test_delta (s->dx, s->x, - 0.001, 0.001); - } - while ((status == GSL_CONTINUE) && (iter < 1000000)); - - gsl_multifit_covar (s->J, 0.0, covar); - - double fit_results[2]; - for(int i=0; i < 2; i++) - fit_results[i] = gsl_vector_get(s->x,i); - - gsl_multifit_fdfsolver_free (s); - gsl_matrix_free (covar); - - return (fit_results[0]+maxX); -} - -int my_f (const gsl_vector * x, void *FitData, gsl_vector * f) -{ - size_t n = ((struct WaveformProcessor::FitData *)FitData)->n; - double *y = ((struct WaveformProcessor::FitData *)FitData)->y; - double *sigma = ((struct WaveformProcessor::FitData *)FitData)->sigma; - - double alpha = gsl_vector_get (x,0); - double beta = gsl_vector_get (x, 1); - - for(size_t i = 0; i < n; i++) - { - double t = i; - double Yi = 0; - - if(t < alpha) - Yi = 0; - else - Yi = beta*(1-exp(-((t-alpha)*(t-alpha))/WID))*exp(-(t-alpha)/DKAY); - - gsl_vector_set (f, i, (Yi - y[i])/sigma[i]); - } - return GSL_SUCCESS; -} - -int expb_df (const gsl_vector * x, void *FitData, gsl_matrix * J) -{ - size_t n = ((struct WaveformProcessor::FitData *)FitData)->n; - double *sigma = ((struct WaveformProcessor::FitData *) FitData)->sigma; - - double alpha = gsl_vector_get (x, 0); - double beta = gsl_vector_get (x, 1); - - double dalf, dbet; - - for (size_t i = 0; i < n; i++) - { - /* Jacobian matrix J(i,j) = dfi / dxj, */ - /* where fi = (Yi - yi)/sigma[i], */ - /* Yi = A * exp(-beta * i) + b */ - /* and the xj are the parameters (A,beta,b) */ - double t = i; - double s = sigma[i]; - - if(t < alpha) - { - dbet = 0; - dalf = 0; - } - else - { - dbet = (1-exp(-(t-alpha)*(t-alpha)/WID))*exp(-(t-alpha)/DKAY); - dalf = beta*exp(-(t-alpha)/DKAY)*((1-exp(-(t-alpha)*(t-alpha)/WID))/DKAY - (2*(t-alpha)/WID)*exp(-(t-alpha)*(t-alpha)/WID)); - } - - gsl_matrix_set (J,i,0, dalf/s); - gsl_matrix_set (J,i,1, dbet/s); - } - return GSL_SUCCESS; -} - -int expb_fdf (const gsl_vector * x, void *FitData, - gsl_vector * f, gsl_matrix * J) -{ - my_f (x, FitData, f); - expb_df (x, FitData, J); - - return GSL_SUCCESS; -} -#else -double spt_analysis(const vector &trace, const int &maxX, const double &ave_baseline, const double &trcQDC, const int &counter) -{ - //Normalize the trace - vector normtrc(trace.size()); - for(unsigned int j = 0; j < trace.size(); j++) - normtrc.at(j) = ((trace.at(j)-ave_baseline)/trcQDC)*6226.55; -// normtrc.at(j) = ((trace.at(j)-ave_baseline)); - - double bassum=0; - for (unsigned int v=0; v < 15; v++) - { - unsigned int baselineValue = normtrc.at(v); - bassum += baselineValue; - } - double ave_bass = bassum / 15; - - double normE=0; - for(unsigned int j = (maxX - WAVEFORMLOW); j < (maxX + WAVEFORMHIGH); j++) - normE += (normtrc.at(j)); - - double delta1=fabs(trace.at(maxX)-trace.at(maxX-1)); - double inv_fcn = 0; - - //pulser settings -// const double sigma1 = 0.102; //solid parameters -// const double amp1 = 1904.293314; - - const double sigma1 = 0.166597; //test parameters - const double amp1 = 1816.27; - - //short bar settings -// const double sigma1 = 0.039185; -// const double amp1 = 1997.157; - - //dev settings - //const double sigma1 = 0.102; //0.097481 0.096199429 - //const double amp1 = 0.3192159*normE-1.361152; //1986.2 - //const double amp1 = 0.3192159*trcQDC-1.361152; //1986.2 - - if(delta1 > ((trace.at(maxX)-ave_baseline)*0.4)) - inv_fcn = (pow(-log((normtrc.at(maxX-1))/(amp1))*pow(sigma1,3),0.25)/sigma1)+(maxX-1); - else if ((delta1 < (trace.at(maxX)-ave_baseline)*0.4) && (normtrc.at(maxX-2) > 0)) - inv_fcn = (pow(-log((normtrc.at(maxX-2))/(amp1))*pow(sigma1,3),0.25)/sigma1)+(maxX-2); - else - inv_fcn = (pow(-log((normtrc.at(maxX-1))/(amp1))*pow(sigma1,3),0.25)/sigma1)+(maxX-1); - - return inv_fcn; -} -#endif //#ifdef pulsefit