TFile *f1=new TFile("fall.root");
TCanvas *c1=new TCanvas;
gPad->SetLogz();
tree->Draw("cathod:Iteration$>>hc(4096,0,4096,2500,4000,14000)","","colz");
hc->SetStats(0);
c1->Draw();
tree->Draw("anode:Iteration$>>ha(4096,0,4096,2500,3000,11000)","","colz");
ha->SetStats(0);
c1->Draw();
float baseline(const float *v, float *s, int npts)
{
float base=0;
float mloss=0;
for(int j=0;j<1900;j++) {
base+=v[j]/1900.;
}
for(int j=0;j<4096;j++) {
s[j]=v[j]-base;
if(j<1900)
mloss+=s[j]*s[j]/1900;
}
return mloss;
}
void medianfilter(const float *v, float *s,int npts, int nM)
{
for(int j=0;j<npts;j++) {
vector<float> vs;
for(int k=-nM/2;k<=nM/2;k++) {
int ch=0;
if(j+k>=0 && j+k<npts) ch=j+k;
if(j+k>=npts) ch=npts-1;
vs.push_back(v[ch]);
}
sort(vs.begin(),vs.end());
s[j]=vs[nM/2];
}
}
void averagefilter(const float *v, float *s,int npts, int nA)
{
for(int j=0;j<npts;j++) {
s[j]=0;
for(int k=-nA/2;k<=nA/2;k++) {
int ch=0;
if(j+k>=0 && j+k<npts) ch=j+k;
if(j+k>=npts) ch=npts-1;
s[j]+=v[ch]/(nA+1);
//cout<<s[j]<<" "<<v[ch]<<endl;
}
}
}
void ana(int run,int nM,int nA)
{
float cathod[4096],anode[4096];
float cwave[4096],awave[4096];
float mcwave[4096],mawave[4096];
float acwave[4096],aawave[4096];
float ea,ec;//energy
float closs,aloss;
//TFile *fin=new TFile(Form("f%02d.root",run));
TFile *fin=new TFile("fall.root");
TTree *tin= (TTree*)fin->Get("tree");
tin->SetBranchAddress("cathod",&cathod);
tin->SetBranchAddress("anode",&anode);
//TFile *fout=new TFile(Form("ana%02d.root",run),"recreate");
TFile *fout=new TFile("fana.root","recreate");
fout->cd();
TTree *tout=new TTree("tree","tree");
tout->Branch("cwave",&cwave,"cwave[4096]/F");
tout->Branch("awave",&awave,"awave[4096]/F");
tout->Branch("closs",&closs,"closs/F");
tout->Branch("aloss",&aloss,"aloss/F");
tout->Branch("mcwave",&mcwave,"mcwave[4096]/F");
tout->Branch("mawave",&mawave,"mawave[4096]/F");
tout->Branch("acwave",&acwave,"acwave[4096]/F");
tout->Branch("aawave",&aawave,"aawave[4096]/F");
tout->Branch("ea",&ea,"ea/F");
tout->Branch("ec",&ec,"ec/F");
Long64_t nentries=tin->GetEntriesFast();
for(Long64_t jentry=0;jentry<nentries;jentry++) {
tin->GetEntry(jentry);
closs=baseline(cathod,cwave,4096);//baseline correction,closs-cathod noise fluctuation
aloss=baseline(anode,awave,4096);//baseline correction, aloss-anode noise fluctuation
medianfilter(cwave,mcwave,4096,nM);//median filter
averagefilter(mcwave,acwave,4096,nA);//moving average
medianfilter(awave,mawave,4096,nM);
averagefilter(mawave,aawave,4096,nA);
ea=TMath::MinElement(4096,aawave);
ec=TMath::MaxElement(4096,acwave);
tout->Fill();
if(jentry%50==0) cout<<jentry<< " ";
}
cout<<endl;
tout->Write();
fout->Close();
}
//ana(1,100,20); create new root file for analysis
TFile *f=new TFile("fana.root");
tree->Draw("acwave:Iteration$>>(2500,1500,4000,2500,-1000,8000)","","colz");//cathod after median->moving average filter
c1->Draw();
tree->Draw("aawave:Iteration$>>(2500,1500,4000,2500,-7000,1000)","","colz");
c1->Draw();
tree->Draw("-ea:ec>>(500,1000,6000,500,0,8000)","","colz");
c1->Draw();