6 #include "TMatrixTSym.h"
7 #include "TMatrixDSymEigen.h"
28 vector<vector<string> >
ParseTSV(
const char* s);
29 TMatrixD
GetNewGains(TMatrixD oldmatrix,
string lastfile);
32 int main(
int argc,
char** argv ){
35 if(
string(argv[1]) ==
"help" ||
string(argv[1]) ==
"-help" ||
string(argv[1]) ==
"--help" ||
string(argv[1]) ==
"-h") {
36 cout <<
"This program expects the following argument(s) in order: " << endl;
37 cout <<
"Argument 1: path to root file that contains C,D,L matrices as histos" << endl;
38 cout <<
"Argument 2: path to root file containing eigenvalues and eigenvectors" << endl;
39 cout <<
"Argument 3: lower limit cut on eigenvalues of C matrix" << endl;
40 cout <<
"Argument 4: old gain factor text file (if not first iteration)" << endl;
44 else if(argc != 4 && argc != 5) {
45 cout <<
"Error: unexpected number of arguments" << endl;
53 vector<vector<string> > oldgains;
56 cout <<
"Current iteration is: 1" << endl;
57 nextroot =
"GainFactorIter1.root";
58 nextfilename =
"GainFactorIter1.txt";
62 char* textfilename = argv[4];
63 const char* it = &textfilename[14];
64 int iteration = atoi(it)+1;
68 cout <<
"Error: Something's wrong with GainFactorIter%i.txt file!!" << endl;
75 cout <<
"Current iteration is " << iteration << endl;
76 char iter[50];
char nexty[50];
77 sprintf(iter,
"GainFactorIter%i.txt",iteration);
78 sprintf(nexty,
"GainFactorIter%i.root",iteration);
79 nextfilename = iter; nextroot = nexty;
80 cout <<
"Next file name: " << nextfilename << endl;
81 lastfilename = argv[4];
83 cout <<
"if it didnt work by now it never will :(" << endl;
88 TFile* eigfile =
new TFile(argv[2]);
89 TFile* matrixfile =
new TFile(argv[1]);
91 TH1F* h1D_mD = (TH1F*)matrixfile->Get(
"h1D_mD");
92 TH1F* h1D_nhits = (TH1F*)matrixfile->Get(
"h1D_nhits");
94 TH2F* h2_evecs = (TH2F*)eigfile->Get(
"h2_evecs DO NOT PLOT");
95 TH1F* h1_evals = (TH1F*)eigfile->Get(
"h1_evals");
97 int n_channels = h1D_mD->GetNbinsX();
99 int m_nElements = n_channels;
101 cout <<
"Found historgram h1D_mD, dimension=" << n_channels << endl;
106 TH1F* h_GainFactors =
new TH1F(
"h_GainFactors",
"Gain Factor as Function of Channel Number" ,n_channels, 0., n_channels );
107 TH1F* h_GainFactors_old =
new TH1F(
"h_GainFactors_old",
"Old Gain Factor as Function of Channel Number" ,n_channels, 0., n_channels );
108 TH1F* h_unordered_new =
new TH1F(
"h_unordered_new",
"Gain Factors Unordered" ,100, 0., 2. );
109 TH1F* h_unordered =
new TH1F(
"h_unordered",
"Gain Factors Unordered" ,100, 0., 2. );
110 TH1F* h_unordered_diff =
new TH1F(
"h_unordered_diff",
"Old gain factor - New Gain Factor" ,100, -1., 1. );
111 TH1F* h_old_unordered =
new TH1F(
"h_old_unordered",
"Old Gain Factors Unordered" ,100, 0., 2. );
112 TH1F* h_GainFactorsNew =
new TH1F(
"h_GainFactorsNew",
"Actual gain factors, ordered",n_channels,0,n_channels);
113 TH2D* h_yvsx =
new TH2D(
"h_yvsx",
"Gain factors y vs x",60,0,60,60,0,60);
120 m_mD.ResizeTo(m_nElements,1);
121 m_mL.ResizeTo(m_nElements,1);
122 m_mLt.ResizeTo(1,m_nElements);
127 evec.ResizeTo(n_channels,n_channels);
131 eval.ResizeTo(n_channels);
135 eiginv.ResizeTo(n_channels,n_channels);
139 evecinv.ResizeTo(n_channels,n_channels);
142 double cutoff = atof(argv[3]);
144 cout <<
"Generating Matrices..." << endl;
145 for(
int i = 0; i<n_channels; ++i) {
146 m_mD[i][0]=h1D_mD->GetBinContent(i+1);
147 if(h1_evals->GetBinContent(i+1) > cutoff) eiginv[i][i] = (1./h1_evals->GetBinContent(i+1));
148 for(
int j = 0 ; j<n_channels; ++j) {
149 evec[i][j] = h2_evecs->GetBinContent(i+1,j+1);
153 evecinv.Transpose(evec);
155 cout <<
"Doing some matrix multiplication now (this could take a while...)" << endl;
156 TMatrixD C_inv = evec * eiginv * evecinv;
160 TMatrixD GainFactors = C_inv*m_mD;
166 oldgains =
ParseTSV(lastfilename.c_str());
169 for(
int i = 0; i<n_channels; ++i) {
170 h_GainFactorsNew->SetBinContent(i+1,GainFactors[i][0]);
171 h_unordered_new->Fill(GainFactors[i][0]);
176 GainFactors =
GetNewGains(GainFactors, lastfilename);
181 next_file.open(nextfilename.c_str());
185 for(
int i=0; i<n_channels; i++){
186 if(h1D_nhits->GetBinContent(i+1) < hits_cut) {
188 GainFactors[i][0] = 1;
193 for(
int i=0; i<n_channels; i++){
195 double my_gain = GainFactors[i][0];
196 ostringstream sstream; sstream << my_gain;
197 string smy_gain = sstream.str();
198 ostringstream channel; channel<<i+1;
204 int row = fcalgeom.
row(i);
207 double r =
sqrt(pos.X()*pos.X() + pos.Y()*pos.Y());
208 next_file << smy_gain <<
" channel=" << i <<
" row=" << row <<
" column=" << column <<
" x=" << pos.X() <<
" y=" << pos.Y() <<
" r=" << r << endl;
209 if (my_gain > 0) h_yvsx->Fill(column,row,my_gain);
219 for(
int i=0; i<n_channels; i++){
223 double oldgain = atof(oldgains[i][0].c_str());
224 h_GainFactors_old->SetBinContent(i+1,oldgain);
225 h_old_unordered->Fill(oldgain);
226 h_unordered_diff->Fill(oldgain-GainFactors[i][0]);
229 h_GainFactors->SetBinContent(i+1,GainFactors[i][0]);
230 h_unordered->Fill(GainFactors[i][0]);
234 TFile* m_rootFile =
new TFile( nextroot.c_str(),
"RECREATE" );
236 h_GainFactors->Write();
237 h_GainFactorsNew->Write();
238 h_unordered_new->Write();
239 h_GainFactors_old->Write();
240 h_unordered->Write();
241 h_old_unordered->Write();
242 h_unordered_diff->Write();
248 cout <<
"Total Channels skipped: " << nskipped << endl;
250 cout <<
"done" << endl;
259 typedef vector<vector<string> > Rows;
262 char const row_delim =
'\n';
263 char const field_delim =
'\t';
264 for (
string row; getline(input, row, row_delim); ) {
265 parsed.push_back(Rows::value_type());
266 istringstream ss(row);
267 for (
string field; getline(ss, field, field_delim); ) {
268 parsed.back().push_back(field);
277 int n_channels = 2800;
280 newmatrix.ResizeTo(n_channels,1);
282 vector<vector<string> > oldgains =
ParseTSV(lastfile.c_str());
284 for(
int i = 0; i<n_channels; ++i) {
285 double lastgain = atof(oldgains[i][0].c_str());
286 newmatrix[i][0] = oldmatrix[i][0]*lastgain;
pair< int, int > AbsNumtoXY(int channel)
sprintf(text,"Post KinFit Cut")
DVector2 positionOnFace(int row, int column) const
vector< vector< string > > ParseTSV(const char *s)
TMatrixD GetNewGains(TMatrixD oldmatrix, string lastfile)
int column(int channel) const
int row(int channel) const
int main(int argc, char *argv[])