Experimental Data Analysis for GIC

Rawdata-waveform of cathod/anode

In [1]:
TFile *f1=new TFile("fall.root");
TCanvas *c1=new TCanvas;
gPad->SetLogz();

cathod waveform

In [2]:
tree->Draw("cathod:Iteration$>>hc(4096,0,4096,2500,4000,14000)","","colz");
hc->SetStats(0);
c1->Draw();

anode waveform

In [3]:
tree->Draw("anode:Iteration$>>ha(4096,0,4096,2500,3000,11000)","","colz");
ha->SetStats(0);
c1->Draw();

Energy:

  • median filtering (100points as a unit) is used firstly to eliminate sharp noise and then moving average filtering (20 points as a unit) is adopted.
  • J.Liu etc., Nuclear Inst. and Methods in Physics Research, A 1004 (2021) 165363

Baseline subtraction

In [4]:
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;
}

median filter

In [5]:
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];
    }
  
}

moving average filter for noise reduction

In [6]:
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;
        }
    }
}

Analysis Routine

In [7]:
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();
}
In [8]:
//ana(1,100,20); create new root file for analysis

Pulses after filters

In [9]:
TFile *f=new TFile("fana.root");
In [10]:
tree->Draw("acwave:Iteration$>>(2500,1500,4000,2500,-1000,8000)","","colz");//cathod after median->moving average filter
c1->Draw();
In [11]:
tree->Draw("aawave:Iteration$>>(2500,1500,4000,2500,-7000,1000)","","colz");
c1->Draw();

anode energy vs. cathod energy

In [12]:
tree->Draw("-ea:ec>>(500,1000,6000,500,0,8000)","","colz");
c1->Draw();