pedestalsCorrection.cxx 2.19 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
#include <iostream>

using namespace std;

void pedestalsCorrection()
{
	TString outPosFile = "parameters/correctedPedsPositions.par";

	TString histName = "fhAdcAllWoBaselineNonLinear";

	TString location = "../data/WPT_tests";

	TString fIn;
	fIn.Form("%s/680pF_7okt/raw_data/TestTriggerMode_680pF_10minutes_7okt_0000.root", location.Data());

	TFile *fr = new TFile(fIn.Data());
	if (fr->IsZombie()) {
		Error("pedestalsCorrection.cxx", "Files %s was not open.", fIn.Data());
		return;
	}

	//output files
	ofstream calPosFile;
	calPosFile.open(outPosFile.Data());
	if (calPosFile.is_open()) {
		Info("calTransition47.cxx", "Positions of 0.5V peak will be saved in %s.", outPosFile.Data());
	} else {
		Error("calTransition47.cxx", "File %s was not open.", outPosFile.Data());
		return;
	}

	TH2D *hA = (TH2D*)fr->Get(histName.Data());

	TCanvas *c1 = new TCanvas();
	c1->Divide(2,2);

	c1->cd(1);
	hA->GetXaxis()->SetRangeUser(0, 64);
	hA->Draw();

	c1->cd(2);

	TString cName;
	TString cTitle;

	TCanvas *cWork[11];
	for (Int_t i = 0; i < 11; i++) {
		cWork[i] = new TCanvas();
		cName.Form("can%d", i);
		cTitle.Form("canvas %d", i);
		cWork[i]->SetName(cName.Data());
		cWork[i]->SetTitle(cTitle.Data());
		cWork[i]->Divide(2, 3);
		cWork[i]->cd(1);
	}

	TH1D *hch[64];

	Double_t x, y;
	Double_t relHight = 0.5;

	for (Int_t i = 0; i < 64; i++) {
		hch[i] = extractHist(hA, i);

		hch[i]->GetXaxis()->SetRangeUser(-100, 500);
		Int_t cNumber = i/6;
		cWork[cNumber]->cd(i-6*cNumber+1);
		hch[i]->Draw();

		hch[i]->ShowPeaks(2., "", relHight);

		TList *functions = hch[i]->GetListOfFunctions();
		TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");

		const Int_t noPeaks = pm->GetN();
		if (noPeaks != 1) {
			Error("calTransition47.cxx", "%d peaks found in hist %s", noPeaks, hch[i]->GetName());
			continue;
		}

		Double_t *energy = pm->GetX();
		calPosFile << i << "\t" << energy[0] << endl;

	}

	calPosFile.close();

	return;

}

TH1D* extractHist(TH2D *hInput, Int_t channel) {

	TH1D *hTemp;
	TString hName;
	TString hTitle;

	hName.Form("ch%d", channel);
	hTitle.Form("channel %d", channel);
	hTemp = hInput->ProjectionY(hName.Data(), channel+1, channel+1);
	hTemp->SetTitle(hTitle.Data());

	return hTemp;

}