Commit 94cfffcb authored by dauncey's avatar dauncey

Rewriting AnalysisPhoton

parent 1ea3f838
No preview for this file type
......@@ -75,8 +75,9 @@ The resulting RMS using the best fit values is given by
\end{eqnarray*}
But since the solution is defined by $Ma=v$, then $a^T M a = a^T v$. Hence
\begin{equation}
\mathrm{RMS}^2_\mathrm{min} = \frac{1}{N} \sum_e T_e^2 - a^T M a
= \frac{1}{N} \sum_e T_e^2 - v^T M^{-1} v
\mathrm{RMS}^2_\mathrm{min} = \frac{1}{N} \left(\sum_e T_e^2\right) - a^T M a
= \frac{1}{N} \left(\sum_e T_e^2\right) - v^T M^{-1} v
= \frac{1}{N} \left(\sum_e T_e^2\right) - a^T v
\end{equation}
The above can be extended slightly, which may improve the energy response
linearity as well as the RMS. The energy estimation for the event (i.e. the
......
......@@ -10,12 +10,15 @@
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TVectorD.h"
#include "TMatrixDSym.h"
#include "Random.hh"
#include "TE1F.hh"
#include "TE2F.hh"
#include "AnalysisBase.hh"
#include "SubAnalysisPhotonSimHit.hh"
#include "SubAnalysisPhotonTrgHit.hh"
#include "Event.hh"
#include "ErrorMatrix.hh"
......@@ -26,7 +29,113 @@ public:
enum {
kNumberOfTests=AdcReading::fMaximumNumberOfMethods+26
};
/*
class Coefficients {
public:
enum {
kNumberOfCoefficients=Geometry::kNumberOfLayers
};
Coefficients() : fSum(kNumberOfCoefficients), fSSq(kNumberOfCoefficients) {
fNum=0.0;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fSum(i)=0.0;
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fSSq(i,j)=0.0;
}
}
}
void fill(double e, const double *p) {
fNum++;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fSum(i)+=e*p[i];
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fSSq(i,j)+=p[i]*p[j];
}
}
}
void normalise() {
assert(fNum>0.5);
if(fNum==1.0) return;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fSum(i)/=fNum;
}
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fSSq(i,j)/=fNum;
}
}
fNum=1.0;
}
const TMatrixDSym& coefficientMatrix() const {
return fSSq;
}
TVectorD coefficients() {
normalise();
TMatrixDSym im(fSSq);
im.Invert();
return im*fSum;
}
void read(const std::string &fName) {
std::ifstream fin(fName.c_str());
if(!fin) return;
unsigned n;
fin >> n;
assert(n==kNumberOfCoefficients);
fin >> fNum;
double d;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fin >> d;
fSum(i)=d;
}
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fin >> d;
fSSq(i,j)=d;
}
}
}
void write(const std::string &fName) {
std::ofstream fout(fName.c_str());
if(!fout) return;
normalise();
fout << kNumberOfCoefficients;
fout << fNum << std::endl;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fout << " " << std::setprecision(12) << fSum(i);
}
fout << std::endl;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fout << " " << std::setprecision(12) << fSSq(i,j);
}
fout << std::endl;
}
}
double fNum;
TVectorD fSum;
TMatrixDSym fSSq;
};
*/
/*
class ErrorMatrix {
public:
......@@ -134,7 +243,18 @@ public:
*/
AnalysisPhoton(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("AnalysisPhoton",sRoot,sOut) {
: AnalysisBase("AnalysisPhoton",sRoot,sOut),
fSAPSimHit(sRoot,sOut), fSAPTrgHit(sRoot,sOut) {
// What to do?
fSimEvent=true;
fDigEvent=false;
fTrgEvent=true;
fC2DEvent=false;
fC3DEvent=false;
hSelectSinThetaAll=new TH1F((fName+"_SelectSinThetaAll").c_str(),
";sin(#theta)",100,0.0,1.0);
......@@ -2387,8 +2507,8 @@ public:
}
virtual ~AnalysisPhoton() {
fErrorMatrix.write("AnalysisPhoton_ErrorMatrix.out");
fErrorMatrix.normalise();
//fErrorMatrix.write("AnalysisPhoton_ErrorMatrix.out");
//fErrorMatrix.normalise();
//TMatrixDSym fm;//(fErrorMatrix.fitMatrix());
/*
......@@ -2405,6 +2525,8 @@ public:
<< ";" << std::endl;
}
*/
/*
for(unsigned i(0);i<40;i++) {
hErrMean->Fill(i,fErrorMatrix.fSum[i]);
std::cout << "fErrMean[" << i << "]=" << fErrorMatrix.fSum[i]
......@@ -2425,7 +2547,8 @@ public:
hErrMatrix->Fill(i,j,em[i][j]);//fErrorMatrix.fSSq[i][j]);
hCovMatrix->Fill(i,j,em[i][j]/sqrt(em[i][i]*em[j][j]));//fErrorMatrix.fSSq[i][j]/sqrt(fErrorMatrix.fSSq[i][i]*fErrorMatrix.fSSq[j][j]));
}
}
}
*/
}
bool fitC3dShape(const double *a, double &et) {
......@@ -2621,10 +2744,15 @@ public:
//double stInv(1.0/proj.position().sinTheta());
/*
if(proj.position().rho()> 30.0 &&
proj.position().rho()<145.0
//&& stInv>=2.5 && stInv<3.0
) good=true;
*/
if(fabs(proj.position().eta())>=1.7 &&
fabs(proj.position().eta())< 2.7) good=true;
}
}
......@@ -2709,8 +2837,14 @@ public:
double tanMax[2];
double phiMax[2];
//simEvent(event);
//return true;
if(fSimEvent) fSAPSimHit.event(event);
//event->digHitNoiseCreator();
//event->digHitProcessor();
//event->trgHitProcessor();
if(fTrgEvent) fSAPTrgHit.event(event);
return true;
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
......@@ -3232,7 +3366,17 @@ public:
}
protected:
bool fSimEvent;
bool fDigEvent;
bool fTrgEvent;
bool fC2DEvent;
bool fC3DEvent;
Calibration cal;
SubAnalysisPhotonSimHit fSAPSimHit;
SubAnalysisPhotonTrgHit fSAPTrgHit;
TH1F *hSelectSinThetaAll,*hSelectSinThInvAll;
TH1F *hSelectSinThetaSel,*hSelectSinThInvSel;
......
#ifndef Coefficients_HH
#define Coefficients_HH
#include "TVectorD.h"
#include "TMatrixDSym.h"
class Coefficients {
public:
enum {
kNumberOfCoefficients=Geometry::kNumberOfHGCLayers+Geometry::kNumberOfBHLayers
};
Coefficients(bool b=false) : fPlainRms(b),
fSum(kNumberOfCoefficients),
fCoe(kNumberOfCoefficients),
fSSq(kNumberOfCoefficients) {
fNum=0.0;
fErg=0.0;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fSum(i)=0.0;
fCoe(i)=0.0;
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fSSq(i,j)=0.0;
}
fSSq(i,i)=1.0e-12;
}
}
void fill(double e, const double *p) {
TVectorD v(kNumberOfCoefficients);
for(unsigned i(0);i<kNumberOfCoefficients;i++) v(i)=p[i];
fill(e,v);
}
void fill(double e, const TVectorD &p) {
assert(p.GetNrows()==kNumberOfCoefficients);
fNum++;
if(fPlainRms) fErg+=e*e;
else fErg+=e;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
if(fPlainRms) fSum(i)+=e*p(i);
else fSum(i)+=p(i);
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
if(fPlainRms) fSSq(i,j)+=p(i)*p(j);
else fSSq(i,j)+=p(i)*p(j)/e;
}
}
}
double solve(uint64_t lm=0xffffffffffffffff) {
assert(fNum>0.5);
assert(lm!=0);
uint64_t one(1);
TVectorD sum(fSum);
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
bool ignore((lm&(one<<i))==0);
if(ignore) sum(i)=0.0;
}
TMatrixDSym ssq(fSSq);
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
bool ignore((lm&(one<<i))==0);
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
bool jgnore((lm&(one<<j))==0);
if(ignore || jgnore) ssq(i,j)=0.0;
}
if(ignore) ssq(i,i)=1.0; // Allow inverse
}
ssq.Invert();
fCoe=ssq*sum;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
bool ignore((lm&(one<<i))==0);
if(ignore) assert(fCoe(i)==0.0);
}
return (fErg-sum*fCoe)/fNum;
}
const TMatrixDSym& coefficientMatrix() const {
return fSSq;
}
const TVectorD& coefficients() {
return fCoe;
}
bool read(const std::string &fName) {
std::ifstream fin(fName.c_str());
if(!fin) return false;
unsigned n;
fin >> n;
if(n!=kNumberOfCoefficients) return false;
fin >> fNum >> fErg;
double d;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fin >> d;
fSum(i)=d;
}
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fin >> d;
fSSq(i,j)=d;
}
}
if(!fin) return false;
return true;
}
bool write(const std::string &fName) {
std::ofstream fout(fName.c_str());
if(!fout) return false;
fout << kNumberOfCoefficients << " "
<< std::setprecision(12) << fNum << " "
<< std::setprecision(12) << fErg
<< std::endl;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
fout << " " << std::setprecision(12) << fSum(i);
}
fout << std::endl;
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
fout << " " << std::setprecision(12) << fSSq(i,j);
}
fout << std::endl;
}
if(!fout) return false;
fout.close();
return true;
}
void print() {
std::cout << "Coefficients::print() Number of entries = "
<< fNum << ", energy^2 " << fErg << std::endl;
std::cout << "Sum vector =";
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
std::cout << " " << fSum(i);
}
std::cout << std::endl;
std::cout << "Sum matrix =";
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
for(unsigned j(0);j<kNumberOfCoefficients;j++) {
std::cout << " " << fSSq(i,j);
}
std::cout << std::endl;
}
std::cout << " Coefficients =";
for(unsigned i(0);i<kNumberOfCoefficients;i++) {
std::cout << " " << fCoe(i);
}
std::cout << std::endl;
}
const bool fPlainRms;
double fNum,fErg;
TVectorD fSum,fCoe;
TMatrixDSym fSSq;
};
#endif
......@@ -602,6 +602,8 @@ public:
}
void trgHitProcessor() {
std::cout << "Event::trgHitProcessor called" << std::endl;
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
trgHitProcessor(e);
}
......
......@@ -49,11 +49,13 @@ class ReadSimHitTTree {
fTTree->SetBranchAddress( "SimHit_E", &_SimHit_E );
fTTree->SetBranchAddress( "SimHit_Id_BH", &_SimHit_Id_BH );
fTTree->SetBranchAddress( "SimHit_E_BH", &_SimHit_E_BH );
fTTree->SetBranchAddress( "StringTestVector", &_SimHit_str );
_SimHit_Id=0;
_SimHit_E=0;
_SimHit_Id_BH=0;
_SimHit_E_BH=0;
_SimHit_str=0;
if(PrintLevel::print(2)) std::cout << "ReadSimHitTTree opening file "
<< fFileName <<" successful" << std::endl;
......@@ -115,7 +117,15 @@ class ReadSimHitTTree {
vHbh.push_back(s);
}
}
if(_SimHit_str!=0) {
std::cout << "_SimHit_str size = " << _SimHit_str->size() << std::endl;
std::vector<std::string> &vStr(*_SimHit_str);
for(unsigned j(0);j<_SimHit_str->size() && j<100;j++){
if(vStr[j]!="abc0") std::cout << "Line " << j << " = " << vStr[j] << std::endl;
}
}
return true;
}
......@@ -163,6 +173,8 @@ class ReadSimHitTTree {
std::vector<UInt_t> *_SimHit_Id_BH;
std::vector<Double_t> *_SimHit_E_BH;
std::vector<std::string> *_SimHit_str;
};
#endif
......@@ -47,6 +47,12 @@ public:
return fEvent;
}
// Position
Point point() const {
return Geometry::point(*this);
}
// Deposited energy in GeV
double deposit() const {
return fScale*fDeposit;
......@@ -64,6 +70,11 @@ public:
return Constants::fMipPerGev[Geometry::numberOfSublayers(layer(),wafer())-1]*deposit();
}
// Equivalent transverse number of MIPs
double transverseMips() const {
return mips()*point().sinTheta();
}
// Equivalent incoming energy in GeV
double energy() const {
return mips()*Calibration::layerWeight(layer());
......@@ -71,7 +82,7 @@ public:
// Equivalent incoming transverse energy in GeV
double transverseEnergy() const {
return energy()*Geometry::point(*this).sinTheta();
return energy()*point().sinTheta();
}
bool read(std::istream &fin) {
......
#ifndef SubAnalysisPhotonSimHit_HH
#define SubAnalysisPhotonSimHit_HH
#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TH1A.hh"
#include "Average.hh"
#include "Coefficients.hh"
#include "utilities.cc"
class SubAnalysisPhotonSimHit {
public:
SubAnalysisPhotonSimHit(const std::string &sRoot="",
const std::string &sOut="")
: fName("SubAnalysisPhotonSimHit"),
fFile((sOut==""?"":sOut+"/")+fName+(sRoot==""?sRoot:"_"+sRoot)) {
nEvt=0;
hEvent=new TH1F((fName+"_Event").c_str(),
";Layer;Number of trgHits",100,0.0,1000.0);
for(unsigned l(0);l<Geometry::kNumberOfLayers;l++) {
hSimHitNumber[l]=new TH1F((fName+"_Layer"+printNum(l,2)+"SimHitNumber").c_str(),
";Number",
100,0.0,1000.0);
}
hCoefficients[0]=new TH1F((fName+"_Coefficients0").c_str(),
";Layer;Coefficient (GeV/MIP)",
100,0.0,100.0);
hCoefficients[1]=new TH1F((fName+"_Coefficients1").c_str(),
";Layer;Coefficient (GeV/MIP)",
100,0.0,100.0);
hCoeffResult[0]=new TH2F((fName+"_CoeffResult0").c_str(),
";Truth E_{T} (GeV);Reco E_{T} (GeV)",
50,0.0,100.0,50,0.0,100.0);
hCoeffResultP[0]=new TProfile((fName+"_CoeffResultP0").c_str(),
";Truth E_{T} (GeV);Reco E_{T} (GeV)",
50,0.0,100.0);
hCoeffResultA[0]=new TH1A((fName+"_CoeffResultA0").c_str(),
";Truth E_{T} (GeV);Reco-Truth E_{T} (GeV)",
50,0.0,100.0);
hCoeffResult[1]=new TH2F((fName+"_CoeffResult1").c_str(),
";Truth E_{T} (GeV);Reco E_{T} (GeV)",
50,0.0,100.0,50,0.0,100.0);
hCoeffResultP[1]=new TProfile((fName+"_CoeffResultP1").c_str(),
";Truth E_{T} (GeV);Reco E_{T} (GeV)",
50,0.0,100.0);
hCoeffResultA[1]=new TH1A((fName+"_CoeffResultA1").c_str(),
";Truth E_{T} (GeV);Reco-Truth E_{T} (GeV)",
50,0.0,100.0);
}
virtual ~SubAnalysisPhotonSimHit() {
Average avg[2][Coefficients::kNumberOfCoefficients];
for(unsigned k(0);k<2;k++) {
vCoefficients[k][4].print();
std::ostringstream sout;
sout << fFile << "_" << k << ".coe";
assert(vCoefficients[k][4].write(sout.str()));
/*
double rmsSqMin(1.0e100);
for(uint64_t lm(0x000ffffffaaa0000);lm<0x0010000000000000;lm++) {
double rmsSq(vCoefficients[k][4].solve());
if(rmsSqMin>rmsSq) {
rmsSqMin=rmsSq;
std::cout << "New minimum " << rmsSqMin << " found for lm = "
<< printBin(lm) << std::endl;
}
}
*/
uint64_t lm(0x000ffffffaaaaaaa);
std::cout << "k = " << k << " RMS^2 for trigger coeffs = "
<< vCoefficients[k][4].solve(lm) << std::endl;
vCoefficients[k][4].print();
std::cout << "k = " << k << " RMS^2 for all coeffs = "
<< vCoefficients[k][4].solve() << std::endl;
vCoefficients[k][4].print();
TVectorD c4(vCoefficients[k][4].coefficients());
assert(vTVectorD[k].size()==vTruthE.size());
double rmsSq(0.0),rmsSqOverE(0.0);
for(unsigned i(0);i<vTVectorD[k].size();i++) {
double eReco(c4*vTVectorD[k][i]);
//std::cout << "k = " << k << ", event " << i
// << ", truth E = " << vTruthE[k][i]
// << ", reco E = " << eReco << std::endl;
rmsSq+=(eReco-vTruthE[i])*(eReco-vTruthE[i]);
rmsSqOverE+=(eReco-vTruthE[i])*(eReco-vTruthE[i])/vTruthE[i];
hCoeffResult[k]->Fill(vTruthE[i],eReco);
hCoeffResultP[k]->Fill(vTruthE[i],eReco);
hCoeffResultA[k]->Fill(vTruthE[i],eReco-vTruthE[i]);
}