//************************************************
// Author: Federica Lionetto
// Created on: 11/12/2014
//************************************************
/*
CMNSubStraightLine is one of the CMN subtraction algorithms.
Pedestal subtracted ADC counts ---> Calculate mean and rms ---> Exclude signal (requiring |ADC-mean|<factor*rms) ---> Calculate mean and rms ---> Exclude signal (requiring |ADC-mean|<factor*rms) ---> Calculate the parameters of the straight line that approximates the data ---> Evaluate the expected ADC counts on each Beetle channel and subtract it to the measured ones ---> Pedestal and CMN subtracted ADC counts.
It loops over events, creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex_1evt), calculates the mean and the rms, ignores strips that may potentially contain a signal and creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_1evt), calculates the mean and the rms, ignores strips that may potentially contain a signal and creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig2_1evt), calculates the parameters of the straight line that approximates the data, evaluates the expected ADC counts on each Beetle channel, subtracts it to every strip, and creates the event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex_1evt).
If there is no common mode noise, the pedestal and the pedestal and CMN versions will look the same.
It creates the distribution of the two parameters of the straight line y = a0+a1*x (ha0Beetlex and ha1Beetlex).
It creates the distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex) and the distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex), comparing the two.
It creates the distribution of the RMS of the ADC count distribution for the following cases:
- pedestal subtracted ADC counts;
- pedestal subtracted ADC counts, after the first outlier removal;
- pedestal subtracted ADC counts, after the second outlier removal.
In beam off data the ADC count distribution after CMN subtraction is not symmetric.
This might be explained by muons coming from the beam line and acting as "signal".
To check this hypothesis, it is interesting to look at the distribution of the Beetle channels contributing to the tail on the right (TH1F *hchTailBeetlex) and the distribution of the event numbers where this happens (TH1F *hevTailBeetlex).
The configuration of the CMN algorithm is defined by factor.
The other arguments are:
- the number of one event to display, to check that the algorithm is working properly;
- the string corresponding to the input ROOT file, where to get the pedestal subtracted ADC counts from;
- the pointer to the output ROOT file, where to save the output of this algorithm;
- the string corresponding to the path where to save the pictures.
CMNSubStraightLine will create a directory in the output ROOT file called "CMNSubStraightLine", where the tree (with CMN1, CMN2, and ADCPedCMNSub) and the histograms, graphs, and canvases will be stored.
*/
// Header guard.
#ifndef __CMNSUBSTRAIGHTLINE_C_INCLUDED__
#define __CMNSUBSTRAIGHTLINE_C_INCLUDED__
void CMNSubStraightLine(float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures);
float yStraightLine(float a0, float a1, float x);
void CMNSubStraightLine(float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures)
{
cout << "**************************************************" << endl;
cout << "Computing CMN according to the straight line algorithm..." << endl;
cout << "**************************************************" << endl;
cout << "**************************************************" << endl;
cout << "Configuration" << endl;
cout << "factor: " << factor << endl;
cout << "sample: " << sample << endl;
cout << "**************************************************" << endl;
// Open input ROOT file.
cout << "Open input ROOT file: " << input_ROOT << endl;
TFile *input = TFile::Open(TString(input_ROOT));
input->cd();
TTree *EventInfo = (TTree *)input->Get("EventInfo");
int NEvents = EventInfo->GetEntries();
TString dirCMNSubStraightLine = "CMNSubStraightLine";
output->cd();
output->mkdir(dirCMNSubStraightLine);
output->cd(dirCMNSubStraightLine);
TString alg_path_to_figures = path_to_figures+dirCMNSubStraightLine;
TString alg_path_to_make_figures = "mkdir -p "+alg_path_to_figures;
system(alg_path_to_make_figures.Data());
// Copy the EventInfo tree from the input to the output file and add some more information.
float a0Beetle1 = 0.;
float a1Beetle1 = 0.;
float a0Beetle2 = 0.;
float a1Beetle2 = 0.;
std::vector<float> ADCPedCMNSub; // Raw ADC - pedestal - CMN.
TTree *EventInfoProcessed = EventInfo->CloneTree(0);
std::vector<float> *ADCProcessed = 0;
TBranch *b_ADCProcessed = 0;
EventInfo->SetBranchAddress("ADCProcessed",&ADCProcessed,&b_ADCProcessed);
EventInfoProcessed->Branch("a0Beetle1",&a0Beetle1);
EventInfoProcessed->Branch("a1Beetle1",&a1Beetle1);
EventInfoProcessed->Branch("a0Beetle2",&a0Beetle2);
EventInfoProcessed->Branch("a1Beetle2",&a1Beetle2);
EventInfoProcessed->Branch("ADCPedCMNSub",&ADCPedCMNSub);
// Event-by-event histograms.
// Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex_1evt).
TH1F *hADCPedSubBeetle1_1evt = new TH1F("hADCPedSubBeetle1_1evt","",400,-200,200); // Modify it.
TCanvas *cADCPedSubBeetle1_1evt = new TCanvas("cADCPedSubBeetle1_1evt","",400,300);
InitHist(hADCPedSubBeetle1_1evt,"Pedestal subtracted ADC counts on Beetle 1","ADC counts","");
TH1F *hADCPedSubBeetle2_1evt = new TH1F("hADCPedSubBeetle2_1evt","",400,-200,200); //Modify it.
TCanvas *cADCPedSubBeetle2_1evt = new TCanvas("cADCPedSubBeetle2_1evt","",400,300);
InitHist(hADCPedSubBeetle2_1evt,"Pedestal subtracted ADC counts on Beetle 2","ADC counts","");
// Ignore strips that may potentially contain a signal.
// Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_1evt).
// First iteration.
TH1F *hADCPedSubBeetle1NoSig_1evt = new TH1F("hADCPedSubBeetle1NoSig_1evt","",400,-200,200); // Modify it.
TCanvas *cADCPedSubBeetle1NoSig_1evt = new TCanvas("cADCPedSubBeetle1NoSig_1evt","",400,300);
InitHist(hADCPedSubBeetle1NoSig_1evt,"Pedestal subtracted ADC counts on Beetle 1 - excluding signal","ADC counts","");
TH1F *hADCPedSubBeetle2NoSig_1evt = new TH1F("hADCPedSubBeetle2NoSig_1evt","",400,-200,200); //Modify it.
TCanvas *cADCPedSubBeetle2NoSig_1evt = new TCanvas("cADCPedSubBeetle2NoSig_1evt","",400,300);
InitHist(hADCPedSubBeetle2NoSig_1evt,"Pedestal subtracted ADC counts on Beetle 2 - excluding signal","ADC counts","");
// Second iteration.
TH1F *hADCPedSubBeetle1NoSig2_1evt = new TH1F("hADCPedSubBeetle1NoSig2_1evt","",400,-200,200); // Modify it.
TCanvas *cADCPedSubBeetle1NoSig2_1evt = new TCanvas("cADCPedSubBeetle1NoSig2_1evt","",400,300);
InitHist(hADCPedSubBeetle1NoSig2_1evt,"Pedestal subtracted ADC counts on Beetle 1 - excluding signal","ADC counts","");
TH1F *hADCPedSubBeetle2NoSig2_1evt = new TH1F("hADCPedSubBeetle2NoSig2_1evt","",400,-200,200); //Modify it.
TCanvas *cADCPedSubBeetle2NoSig2_1evt = new TCanvas("cADCPedSubBeetle2NoSig2_1evt","",400,300);
InitHist(hADCPedSubBeetle2NoSig2_1evt,"Pedestal subtracted ADC counts on Beetle 2 - excluding signal","ADC counts","");
// Event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex_1evt).
TH1F *hADCPedCMNSubBeetle1_1evt = new TH1F("hADCPedCMNSubBeetle1_1evt","",400,-200,200); // Modify it.
TCanvas *cADCPedCMNSubBeetle1_1evt = new TCanvas("cADCPedCMNSubBeetle1_1evt","",400,300);
InitHist(hADCPedCMNSubBeetle1_1evt,"Pedestal and CMN subtracted ADC counts on Beetle 1","ADC counts","");
TH1F *hADCPedCMNSubBeetle2_1evt = new TH1F("hADCPedCMNSubBeetle2_1evt","",400,-200,200); //Modify it.
TCanvas *cADCPedCMNSubBeetle2_1evt = new TCanvas("cADCPedCMNSubBeetle2_1evt","",400,300);
InitHist(hADCPedCMNSubBeetle2_1evt,"Pedestal and CMN subtracted ADC counts on Beetle 2","ADC counts","");
// Event-by-event distribution of how the event looks like, to check that the straight line is determined properly.
TH1F *hADCStraightLineBeetle1_1evt = new TH1F("hADCStraightLineBeetle1_1evt","",NBeetle,0,NBeetle);
TCanvas *cADCStraightLineBeetle1_1evt = new TCanvas("cADCStraightLineBeetle1_1evt","",400,300);
InitHist(hADCStraightLineBeetle1_1evt,"Event display of Beetle 1","Beetle channel","ADC counts");
TH1F *hADCStraightLineBeetle2_1evt = new TH1F("hADCStraightLineBeetle2_1evt","",NBeetle,NBeetle,N);
TCanvas *cADCStraightLineBeetle2_1evt = new TCanvas("cADCStraightLineBeetle2_1evt","",400,300);
InitHist(hADCStraightLineBeetle2_1evt,"Event display of Beetle 2","Beetle channel","ADC counts");
// Histograms containing the information on all events.
// Beeetle 1.
TH1F *ha0Beetle1 = new TH1F("ha0Beetle1","",100,-50,50); // Modify it.
TCanvas *ca0Beetle1 = new TCanvas("ca0Beetle1","",400,300);
InitHist(ha0Beetle1,"a0 parameter for Beetle 1","a_{0} (ADC counts)","");
TH1F *ha1Beetle1 = new TH1F("ha1Beetle1","",100,-0.5,0.5); // Modify it.
TCanvas *ca1Beetle1 = new TCanvas("ca1Beetle1","",400,300);
InitHist(ha1Beetle1,"a1 parameter for Beetle 1","a_{1} (ADC counts)","");
// Beetle 2.
TH1F *ha0Beetle2 = new TH1F("ha0Beetle2","",100,-50,50); // Modify it.
TCanvas *ca0Beetle2 = new TCanvas("ca0Beetle2","",400,300);
InitHist(ha0Beetle2,"a0 parameter for Beetle 2","a_{0} (ADC counts)","");
TH1F *ha1Beetle2 = new TH1F("ha1Beetle2","",100,-0.5,0.5); // Modify it.
TCanvas *ca1Beetle2 = new TCanvas("ca1Beetle2","",400,300);
InitHist(ha1Beetle2,"a1 parameter for Beetle 2","a_{1} (ADC counts)","");
// Distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex).
TH1F *hADCPedSubBeetle1 = new TH1F("hADCPedSubBeetle1","",400,-200,200);
InitHist(hADCPedSubBeetle1,"Pedestal subtracted ADC counts on Beetle 1","ADC counts","");
TH1F *hADCPedSubBeetle2 = new TH1F("hADCPedSubBeetle2","",400,-200,200);
InitHist(hADCPedSubBeetle2,"Pedestal subtracted ADC counts on Beetle 2","ADC counts","");
// Distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex).
TH1F *hADCPedCMNSubBeetle1 = new TH1F("hADCPedCMNSubBeetle1","",400,-200,200);
InitHist(hADCPedCMNSubBeetle1,"Pedestal and CMN subtracted ADC counts on Beetle 1","ADC counts","");
TH1F *hADCPedCMNSubBeetle2 = new TH1F("hADCPedCMNSubBeetle2","",400,-200,200);
InitHist(hADCPedCMNSubBeetle2,"Pedestal and CMN subtracted ADC counts on Beetle 2","ADC counts","");
TCanvas *cADCCompBeetle1 = new TCanvas("cADCCompBeetle1","",400,300);
TCanvas *cADCCompBeetle2 = new TCanvas("cADCCompBeetle2","",400,300);
// Distribution of the RMS of the ADC count distribution for the following cases:
// - pedestal subtracted ADC counts;
// - pedestal subtracted ADC counts, after the first outlier removal;
// - pedestal subtracted ADC counts, after the second outlier removal.
TH1F *hRMSPedSubBeetle1 = new TH1F("hRMSPedSubBeetle1","",90,0,30);
TH1F *hRMSPedSubBeetle1NoSig = new TH1F("hRMSPedSubBeetle1NoSig","",90,0,30);
TH1F *hRMSPedSubBeetle1NoSig2 = new TH1F("hRMSPedSubBeetle1NoSig2","",90,0,30);
TH1F *hRMSPedSubBeetle2 = new TH1F("hRMSPedSubBeetle2","",90,0,30);
TH1F *hRMSPedSubBeetle2NoSig = new TH1F("hRMSPedSubBeetle2NoSig","",90,0,30);
TH1F *hRMSPedSubBeetle2NoSig2 = new TH1F("hRMSPedSubBeetle2NoSig2","",90,0,30);
InitHist(hRMSPedSubBeetle1,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
InitHist(hRMSPedSubBeetle1NoSig,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
InitHist(hRMSPedSubBeetle1NoSig2,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
InitHist(hRMSPedSubBeetle2,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");
InitHist(hRMSPedSubBeetle2NoSig,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");
InitHist(hRMSPedSubBeetle2NoSig2,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");
TCanvas *cRMSBeetle1 = new TCanvas("cRMSBeetle1","",400,300);
TCanvas *cRMSBeetle2 = new TCanvas("cRMSBeetle2","",400,300);
// In beam off data the ADC count distribution after CMN subtraction is not symmetric.
// This might be explained by muons coming from the beam line and acting as "signal".
// To check this hypothesis, it is interesting to look at the distribution of the Beetle channels contributing to the tail on the right (TH1F *hchTailBeetlex) and the distribution of the event numbers where this happens (TH1F *hevTailBeetlex).
TH1F *hchTailBeetle1 = new TH1F("hcdTailBeetle1","",NBeetle,0,NBeetle);
TH1F *hchTailBeetle2 = new TH1F("hcdTailBeetle2","",NBeetle,NBeetle,N);
InitHist(hchTailBeetle1,"Asymmetric tail on Beetle 1","Beetle channel","");
InitHist(hchTailBeetle2,"Asymmetric tail on Beetle 2","Beetle channel","");
TCanvas *cchTailBeetle1 = new TCanvas("cchTailBeetle1","",400,300);
TCanvas *cchTailBeetle2 = new TCanvas("cchTailBeetle2","",400,300);
TH1F *hevTailBeetle1 = new TH1F("hevTailBeetle1","",50,0,NEvents+1);
TH1F *hevTailBeetle2 = new TH1F("hevTailBeetle2","",50,0,NEvents+1);
InitHist(hevTailBeetle1,"Asymmetric tail on Beetle 1","Event number","");
InitHist(hevTailBeetle2,"Asymmetric tail on Beetle 2","Event number","");
TCanvas *cevTailBeetle1 = new TCanvas("cevTailBeetle1","",400,300);
TCanvas *cevTailBeetle2 = new TCanvas("cevTailBeetle2","",400,300);
// Add the histograms with the pedestal and CMN subtracted distribution of the ADC counts of each Beetle channel (TH1F *hADCPedCMNSub[N]) to the output file.
TH1F *hADCPedCMNSub[N];
for (int iChannel=0;iChannel<N;++iChannel) {
hADCPedCMNSub[iChannel] = new TH1F(Form("hADCPedCMNSub%d",iChannel),"",1024,0,1024);
InitHist(hADCPedCMNSub[iChannel],Form("Pedestal and CMN subtracted ADC counts on Beetle channel %d",iChannel),"ADC counts","");
}
// Beetle 1.
float meanBeetle1_temp;
float rmsBeetle1_temp;
float xmeanBeetle1 = 0.;
float ymeanBeetle1 = 0.;
float x2meanBeetle1 = 0.;
float xymeanBeetle1 = 0.;
float nAvBeetle1 = 0.;
TF1 *fBeetle1 = new TF1("fBeetle1","pol1",0.,(float)NBeetle);
fBeetle1->SetLineColor(kAzure);
// Beetle2.
float meanBeetle2_temp;
float rmsBeetle2_temp;
float xmeanBeetle2 = 0.;
float ymeanBeetle2 = 0.;
float x2meanBeetle2 = 0.;
float xymeanBeetle2 = 0.;
float nAvBeetle2 = 0.;
TF1 *fBeetle2 = new TF1("fBeetle2","pol1",(float)NBeetle,(float)N);
fBeetle2->SetLineColor(kAzure);
// Fill histograms.
for (int iEvent=0;iEvent<NEvents;++iEvent) {
if (iEvent%1000==0)
cout << "Processing event " << iEvent << endl;
// Reset.
if (iEvent>0)
{
ADCPedCMNSub.clear();
hADCPedSubBeetle1_1evt->Reset();
hADCPedSubBeetle2_1evt->Reset();
hADCPedSubBeetle1NoSig_1evt->Reset();
hADCPedSubBeetle2NoSig_1evt->Reset();
hADCPedSubBeetle1NoSig2_1evt->Reset();
hADCPedSubBeetle2NoSig2_1evt->Reset();
hADCPedCMNSubBeetle1_1evt->Reset();
hADCPedCMNSubBeetle2_1evt->Reset();
hADCStraightLineBeetle1_1evt->Reset();
hADCStraightLineBeetle2_1evt->Reset();
xmeanBeetle1 = 0.;
ymeanBeetle1 = 0.;
x2meanBeetle1 = 0.;
xymeanBeetle1 = 0.;
nAvBeetle1 = 0.;
xmeanBeetle2 = 0.;
ymeanBeetle2 = 0.;
x2meanBeetle2 = 0.;
xymeanBeetle2 = 0.;
nAvBeetle2 = 0.;
}
EventInfo->GetEntry(iEvent);
// Fill histograms for Beetle 1.
for (int iChannel=0;iChannel<NBeetle;++iChannel) {
if (maskBeetle1[iChannel] == 0)
{
hADCPedSubBeetle1->Fill((float)ADCProcessed->at(iChannel));
hADCPedSubBeetle1_1evt->Fill((float)ADCProcessed->at(iChannel));
}
}
// Exclude signal.
meanBeetle1_temp = hADCPedSubBeetle1_1evt->GetMean();
rmsBeetle1_temp = hADCPedSubBeetle1_1evt->GetRMS();
hRMSPedSubBeetle1->Fill(rmsBeetle1_temp);
// Fill histograms again for Beetle 1.
for (int iChannel=0;iChannel<NBeetle;++iChannel) {
if (maskBeetle1[iChannel] == 0)
{
if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1_temp)<factor*rmsBeetle1_temp)
hADCPedSubBeetle1NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
}
}
// Exclude signal.
meanBeetle1_temp = hADCPedSubBeetle1NoSig_1evt->GetMean();
rmsBeetle1_temp = hADCPedSubBeetle1NoSig_1evt->GetRMS();
hRMSPedSubBeetle1NoSig->Fill(rmsBeetle1_temp);
// Fill histograms again for Beetle 1.
for (int iChannel=0;iChannel<NBeetle;++iChannel) {
if (maskBeetle1[iChannel] == 0)
{
if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1_temp)<factor*rmsBeetle1_temp)
{
hADCPedSubBeetle1NoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
xmeanBeetle1 = xmeanBeetle1+(float)iChannel;
ymeanBeetle1 = ymeanBeetle1+(float)ADCProcessed->at(iChannel);
x2meanBeetle1 = x2meanBeetle1+((float)iChannel)*((float)iChannel);
xymeanBeetle1 = xymeanBeetle1+(float)iChannel*(float)ADCProcessed->at(iChannel);
nAvBeetle1++;
hADCStraightLineBeetle1_1evt->Fill((float)iChannel,(float)ADCProcessed->at(iChannel));
}
}
}
// Divide by nAvBeetle1 to get the averaged quantities.
xmeanBeetle1 = xmeanBeetle1/nAvBeetle1;
ymeanBeetle1 = ymeanBeetle1/nAvBeetle1;
x2meanBeetle1 = x2meanBeetle1/nAvBeetle1;
xymeanBeetle1 = xymeanBeetle1/nAvBeetle1;
// Calculate CMN1.
a0Beetle1 = (x2meanBeetle1*ymeanBeetle1-xmeanBeetle1*xymeanBeetle1)/(x2meanBeetle1-xmeanBeetle1*xmeanBeetle1);
a1Beetle1 = (xymeanBeetle1-xmeanBeetle1*ymeanBeetle1)/(x2meanBeetle1-xmeanBeetle1*xmeanBeetle1);
ha0Beetle1->Fill(a0Beetle1);
ha1Beetle1->Fill(a1Beetle1);
// Draw CMN1.
fBeetle1->SetParameter(0,a0Beetle1);
fBeetle1->SetParameter(1,a1Beetle1);
if (iEvent == sample)
{
cout << "Beetle 1" << endl;
cout << "Value of a0: " << a0Beetle1 << endl;
cout << "Value of a1: " << a1Beetle1 << endl;
cout << "Value of a0: " << fBeetle1->GetParameter(0) << endl;
cout << "Value of a1: " << fBeetle1->GetParameter(1) << endl;
}
rmsBeetle1_temp = hADCPedSubBeetle1NoSig2_1evt->GetRMS();
hRMSPedSubBeetle1NoSig2->Fill(rmsBeetle1_temp);
// Fill histograms for Beetle 2.
for (int iChannel=NBeetle;iChannel<N;++iChannel) {
if (maskBeetle2[iChannel-NBeetle] == 0)
{
hADCPedSubBeetle2->Fill((float)ADCProcessed->at(iChannel));
hADCPedSubBeetle2_1evt->Fill((float)ADCProcessed->at(iChannel));
}
}
// Exclude signal.
meanBeetle2_temp = hADCPedSubBeetle2_1evt->GetMean();
rmsBeetle2_temp = hADCPedSubBeetle2_1evt->GetRMS();
hRMSPedSubBeetle2->Fill(rmsBeetle2_temp);
// Fill histograms again for Beetle 2.
for (int iChannel=NBeetle;iChannel<N;++iChannel) {
if (maskBeetle2[iChannel-NBeetle] == 0)
{
if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp)<factor*rmsBeetle2_temp)
hADCPedSubBeetle2NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
}
}
// Exclude signal.
meanBeetle2_temp = hADCPedSubBeetle2NoSig_1evt->GetMean();
rmsBeetle2_temp = hADCPedSubBeetle2NoSig_1evt->GetRMS();
hRMSPedSubBeetle2NoSig->Fill(rmsBeetle2_temp);
// Fill histograms again for Beetle 2.
for (int iChannel=NBeetle;iChannel<N;++iChannel) {
if (maskBeetle2[iChannel-NBeetle] == 0)
{
if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp)<factor*rmsBeetle2_temp)
{
hADCPedSubBeetle2NoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
xmeanBeetle2 = xmeanBeetle2+(float)iChannel;
ymeanBeetle2 = ymeanBeetle2+(float)ADCProcessed->at(iChannel);
x2meanBeetle2 = x2meanBeetle2+((float)iChannel)*((float)iChannel);
xymeanBeetle2 = xymeanBeetle2+(float)iChannel*(float)ADCProcessed->at(iChannel);
nAvBeetle2++;
hADCStraightLineBeetle2_1evt->Fill((float)iChannel,(float)ADCProcessed->at(iChannel));
}
}
}
// Divide by nAvBeetle2 to get the averaged quantities.
xmeanBeetle2 = xmeanBeetle2/nAvBeetle2;
ymeanBeetle2 = ymeanBeetle2/nAvBeetle2;
x2meanBeetle2 = x2meanBeetle2/nAvBeetle2;
xymeanBeetle2 = xymeanBeetle2/nAvBeetle2;
// Calculate CMN2.
a0Beetle2 = (x2meanBeetle2*ymeanBeetle2-xmeanBeetle2*xymeanBeetle2)/(x2meanBeetle2-xmeanBeetle2*xmeanBeetle2);
a1Beetle2 = (xymeanBeetle2-xmeanBeetle2*ymeanBeetle2)/(x2meanBeetle2-xmeanBeetle2*xmeanBeetle2);
ha0Beetle2->Fill(a0Beetle2);
ha1Beetle2->Fill(a1Beetle2);
// Draw CMN2.
fBeetle2->SetParameter(0,a0Beetle2);
fBeetle2->SetParameter(1,a1Beetle2);
if (iEvent == sample)
{
cout << "Beetle 2" << endl;
cout << "Value of a0: " << a0Beetle2 << endl;
cout << "Value of a1: " << a1Beetle2 << endl;
cout << "Value of a0: " << fBeetle2->GetParameter(0) << endl;
cout << "Value of a1: " << fBeetle2->GetParameter(1) << endl;
}
rmsBeetle2_temp = hADCPedSubBeetle2NoSig2_1evt->GetRMS();
hRMSPedSubBeetle2NoSig2->Fill(rmsBeetle2_temp);
// Subtract CMN1.
for (int iChannel=0;iChannel<NBeetle;++iChannel) {
ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle1,a1Beetle1,(float)iChannel));
hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle1,a1Beetle1,(float)iChannel));
if (maskBeetle1[iChannel] == 0)
{
hADCPedCMNSubBeetle1->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle1,a1Beetle1,(float)iChannel));
hADCPedCMNSubBeetle1_1evt->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle1,a1Beetle1,(float)iChannel));
if (((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle1,a1Beetle1,(float)iChannel)) > 30.)
{
cout << "The event #" << iEvent << " has at least one Beetle channels (" << iChannel << ") contributing to the tail" << endl;
hchTailBeetle1->Fill(iChannel);
hevTailBeetle1->Fill(iEvent);
}
}
}
// Subtract CMN2.
for (int iChannel=NBeetle;iChannel<N;++iChannel) {
ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle2,a1Beetle2,(float)iChannel));
hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle2,a1Beetle2,(float)iChannel));
if (maskBeetle2[iChannel-NBeetle] == 0)
{
hADCPedCMNSubBeetle2->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle2,a1Beetle2,(float)iChannel));
hADCPedCMNSubBeetle2_1evt->Fill((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle2,a1Beetle2,(float)iChannel));
if (((float)ADCProcessed->at(iChannel)-yStraightLine(a0Beetle2,a1Beetle2,(float)iChannel)) > 30.)
{
cout << "The event #" << iEvent << " has at least one Beetle channels (" << iChannel << ") contributing to the tail" << endl;
hchTailBeetle2->Fill(iChannel);
hevTailBeetle2->Fill(iEvent);
}
}
}
EventInfoProcessed->Fill();
if (iEvent==sample) {
DrawHist(cADCPedSubBeetle1_1evt,hADCPedSubBeetle1_1evt,"",alg_path_to_figures);
DrawHist(cADCPedSubBeetle2_1evt,hADCPedSubBeetle2_1evt,"",alg_path_to_figures);
DrawHist(cADCPedSubBeetle1NoSig_1evt,hADCPedSubBeetle1NoSig_1evt,"",alg_path_to_figures);
DrawHist(cADCPedSubBeetle2NoSig_1evt,hADCPedSubBeetle2NoSig_1evt,"",alg_path_to_figures);
DrawHist(cADCPedSubBeetle1NoSig2_1evt,hADCPedSubBeetle1NoSig2_1evt,"",alg_path_to_figures);
DrawHist(cADCPedSubBeetle2NoSig2_1evt,hADCPedSubBeetle2NoSig2_1evt,"",alg_path_to_figures);
DrawHist(cADCPedCMNSubBeetle1_1evt,hADCPedCMNSubBeetle1_1evt,"",alg_path_to_figures);
DrawHist(cADCPedCMNSubBeetle2_1evt,hADCPedCMNSubBeetle2_1evt,"",alg_path_to_figures);
hADCStraightLineBeetle1_1evt->SetStats(0);
hADCStraightLineBeetle2_1evt->SetStats(0);
DrawHistFunc(cADCStraightLineBeetle1_1evt,hADCStraightLineBeetle1_1evt,fBeetle1,alg_path_to_figures);
DrawHistFunc(cADCStraightLineBeetle2_1evt,hADCStraightLineBeetle2_1evt,fBeetle2,alg_path_to_figures);
}
}
DrawHist(ca0Beetle1,ha0Beetle1,"",alg_path_to_figures);
DrawHist(ca1Beetle1,ha1Beetle1,"",alg_path_to_figures);
DrawHist(ca0Beetle2,ha0Beetle2,"",alg_path_to_figures);
DrawHist(ca1Beetle2,ha1Beetle2,"",alg_path_to_figures);
TLegend *legADCCompBeetle1 = CreateLegend2(hADCPedSubBeetle1,"ped sub",hADCPedCMNSubBeetle1,"ped + CMN sub");
TLegend *legADCCompBeetle2 = CreateLegend2(hADCPedSubBeetle2,"ped sub",hADCPedCMNSubBeetle2,"ped + CMN sub");
cADCCompBeetle1->SetLogy();
cADCCompBeetle2->SetLogy();
DrawHistCompare(cADCCompBeetle1,hADCPedSubBeetle1,hADCPedCMNSubBeetle1,legADCCompBeetle1,alg_path_to_figures);
DrawHistCompare(cADCCompBeetle2,hADCPedSubBeetle2,hADCPedCMNSubBeetle2,legADCCompBeetle2,alg_path_to_figures);
TLegend *legRMSBeetle1 = CreateLegend3(hRMSPedSubBeetle1,"ped sub",hRMSPedSubBeetle1NoSig,"after 1st outlier removal",hRMSPedSubBeetle1NoSig2,"after 2nd outlier removal","lpfw",0.42,0.59,0.94,0.89);
TLegend *legRMSBeetle2 = CreateLegend3(hRMSPedSubBeetle2,"ped sub",hRMSPedSubBeetle2NoSig,"after 1st outlier removal",hRMSPedSubBeetle2NoSig2,"after 2nd outlier removal","lpfw",0.42,0.59,0.94,0.89);
DrawHistCompare3(cRMSBeetle1,hRMSPedSubBeetle1,hRMSPedSubBeetle1NoSig,hRMSPedSubBeetle1NoSig2,legRMSBeetle1,alg_path_to_figures);
DrawHistCompare3(cRMSBeetle2,hRMSPedSubBeetle2,hRMSPedSubBeetle2NoSig,hRMSPedSubBeetle2NoSig2,legRMSBeetle2,alg_path_to_figures);
DrawHist(cchTailBeetle1,hchTailBeetle1,"E",alg_path_to_figures);
DrawHist(cchTailBeetle2,hchTailBeetle2,"E",alg_path_to_figures);
DrawHist(cevTailBeetle1,hevTailBeetle1,"E",alg_path_to_figures);
DrawHist(cevTailBeetle2,hevTailBeetle2,"E",alg_path_to_figures);
output->cd();
// Close input ROOT file.
input->Close();
return;
}
float yStraightLine(float a0, float a1, float x)
{
return a0+a1*x;
}
#endif
Federica Lionetto