Commit 00a99df5 authored by dauncey's avatar dauncey
Browse files

Add Hack.hh for Mac ROOT problem and Gaussian PU mean

parent 89394d5a
No preview for this file type
......@@ -196,6 +196,7 @@ i_x \\ i_y
\end{pmatrix}
\end{equation}
\newpage
\section{Optimisation of layer weights}
Each event $e$ gives some values $d_{e,i}$ of the deposited energy in
layer $i$; these
......
......@@ -18,7 +18,7 @@ public:
const std::string &sRoot="",
const std::string &sOut="") // ctor
: fTwoPi(2.0*acos(-1.0)), fName(sName),
fTFileHandler((sOut==""?"":sOut+"/")+fName+(sRoot==""?sRoot:"_"+sRoot)) {
fTFileHandler((sOut==""?"":sOut+"/")+fName+(sRoot==""?"/Base":"/"+sRoot)) {
fEventSelect=0;
fPrintLevel=0;
......
This diff is collapsed.
......@@ -46,7 +46,10 @@ public:
hPuGenEvent=new TH1F((fName+"_PuGenEvent").c_str(),
";Number of PU events",
400,0.0,400.0);
100,0.0,500.0);
hParticlesVsPuGenEvent=new TH2F((fName+"_ParticlesVsPuGenEvent").c_str(),
";Number of PU events;Number of particles",
100,0.0,500.0,100,0.0,500.0);
hSimHitEventNumber[0]=new TH1F((fName+"_SimHitEventNumberSig").c_str(),
";Signal SimHit event number",
......@@ -57,28 +60,28 @@ public:
hMipsVsLayer[0]=new TH2F((fName+"_"+"MipsVsLayer0").c_str(),
";Layer;MIPs",
52,0.0,52.0,100,0.0,10.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,100,0.0,10.0);
hMipsVsLayerP[0]=new TProfile((fName+"_"+"MipsVsLayer0P").c_str(),
";Layer;MIPs",
52,0.0,52.0,0.0,10.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,0.0,10.0);
hMipsVsLayer[1]=new TH2F((fName+"_"+"MipsVsLayer1").c_str(),
";Layer;MIPs",
52,0.0,52.0,100,0.0,100.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,100,0.0,100.0);
hMipsVsLayerP[1]=new TProfile((fName+"_"+"MipsVsLayer1P").c_str(),
";Layer;MIPs",
52,0.0,52.0,0.0,100.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,0.0,100.0);
hMipsVsLayer[2]=new TH2F((fName+"_"+"MipsVsLayer2").c_str(),
";Layer;MIPs",
52,0.0,52.0,100,0.0,1000.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,100,0.0,1000.0);
hMipsVsLayerP[2]=new TProfile((fName+"_"+"MipsVsLayer2P").c_str(),
";Layer;MIPs",
52,0.0,52.0,0.0,1000.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,0.0,1000.0);
hMipsVsLayer[3]=new TH2F((fName+"_"+"MipsVsLayer3").c_str(),
";Layer;MIPs",
52,0.0,52.0,100,0.0,10000.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,100,0.0,10000.0);
hMipsVsLayerP[3]=new TProfile((fName+"_"+"MipsVsLayer3P").c_str(),
";Layer;MIPs",
52,0.0,52.0,0.0,10000.0);
Geometry::kNumberOfLayers+2,0.0,Geometry::kNumberOfLayers+2,0.0,10000.0);
hPhotonEta=new TH1R((fName+"_"+"PhotonEta").c_str(),
";Truth photon |#eta|",
......@@ -381,11 +384,14 @@ public:
unsigned simHitNum[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
double simHitSum[Geometry::kNumberOfEndcaps][Geometry::kNumberOfLayers];
unsigned nParticles(0);
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
std::vector<Particle> &vL1Track(event.l1Tracks(e));
std::cout << fName << " Event " << fEventNumber
<< " numbers of particles in endcap " << e << " = "
<< vL1Track.size() << std::endl;
nParticles+=vL1Track.size();
std::vector<Particle> &vPhoton(event.photons(e));
std::cout << fName << " Event " << fEventNumber
......@@ -603,14 +609,14 @@ public:
hHitNumber[l]->Fill(nSimHitLayer);
hHitNumber[Geometry::kNumberOfLayers]->Fill(nSimHit);
hMipsVsLayer[0]->Fill(l,simHitMips[l]);
hMipsVsLayerP[0]->Fill(l,simHitMips[l]);
hMipsVsLayer[1]->Fill(l,simHitMips[l]);
hMipsVsLayerP[1]->Fill(l,simHitMips[l]);
hMipsVsLayer[2]->Fill(l,simHitMips[l]);
hMipsVsLayerP[2]->Fill(l,simHitMips[l]);
hMipsVsLayer[3]->Fill(l,simHitMips[l]);
hMipsVsLayerP[3]->Fill(l,simHitMips[l]);
hMipsVsLayer[0]->Fill(l+1,simHitMips[l]);
hMipsVsLayerP[0]->Fill(l+1,simHitMips[l]);
hMipsVsLayer[1]->Fill(l+1,simHitMips[l]);
hMipsVsLayerP[1]->Fill(l+1,simHitMips[l]);
hMipsVsLayer[2]->Fill(l+1,simHitMips[l]);
hMipsVsLayerP[2]->Fill(l+1,simHitMips[l]);
hMipsVsLayer[3]->Fill(l+1,simHitMips[l]);
hMipsVsLayerP[3]->Fill(l+1,simHitMips[l]);
}
if(fabs(simHitMips[0]-event.simMips(e,0))<1.0e-5*event.simMips(e,0)) {
......@@ -700,6 +706,9 @@ public:
*/
}
hParticlesVsPuGenEvent->Fill(event.puGenEvent().size(),nParticles);
bool noHits(true);
for(unsigned l(0);l<Geometry::kNumberOfLayers;l++) {
if(simHitNum[0][l]>0 || simHitNum[1][l]>0) noHits=false;
......@@ -745,6 +754,7 @@ protected:
TH1F *hGenTrackEta;
TH1F *hPuGenEvent;
TH2F *hParticlesVsPuGenEvent;
TH1F *hSimHitEventNumber[2];
......
......@@ -79,6 +79,15 @@ public:
else fDetId&=0xfffffdff;
}
bool pbxSelected() const {
return (fDetId&0x00001000)!=0;
}
void pbxSelected(bool f) {
if(f) fDetId|=0x00001000;
else fDetId&=0xffffefff;
}
bool tpgSelected() const {
return (fDetId&0x00000400)!=0;
}
......
......@@ -5,25 +5,26 @@
#include <cstdint>
#include <vector>
#include "AdcReading.hh"
#include "AdcReadingMultiple.hh"
#include "DigHitDetId.hh"
#include "Random.hh"
#include "SimHit.hh"
#include "DigNoise.hh"
class DigHitMultiple : public DigHitDetId, public AdcReading {
class DigHitMultiple : public DigHitDetId, public AdcReadingMultiple {
public:
/* Now in AdcReadingMultiple
enum {
kNumberOfBx=6
};
*/
DigHitMultiple()
//: RecDetId(0) {
: DigHitDetId(0) {
fSimHitDepositSignal=0.0;
for(unsigned bx(0);bx<kNumberOfBx;bx++) {
for(unsigned bx(0);bx<AdcReadingMultiple::kNumberOfBx;bx++) {
fSimHitDepositPileup[bx]=0.0;
}
}
......@@ -92,14 +93,14 @@ public:
void zeroSimHitDeposits() {
fSimHitDepositSignal=0.0;
fSimHitDepositPileup[fBxNumber%kNumberOfBx]=0.0;
fSimHitDepositPileup[AdcReadingMultiple::currentBx()]=0.0;
}
void addSimHitDeposit(float d, bool signal=false) {
if(invalid()) return;
if(signal) fSimHitDepositSignal+=d;
else fSimHitDepositPileup[fBxNumber%kNumberOfBx]+=d;
else fSimHitDepositPileup[AdcReadingMultiple::currentBx()]+=d;
}
bool addSimHitVector(std::vector<SimHit> &s) {
......@@ -154,32 +155,32 @@ public:
// Result
double charge() const {
return chargeResult();
double charge(unsigned dBx=0) const {
return chargeResult(dBx);
}
double sigmas() const {
return charge()/sigma();
double sigmas(unsigned dBx=0) const {
return charge(dBx)/sigma();
}
double deposit() const {
return charge()/Constants::fFcPerGev;
double deposit(unsigned dBx=0) const {
return charge(dBx)/Constants::fFcPerGev;
}
double mips() const {
return charge()*Constants::fMipPerFc[Geometry::numberOfSublayers(layer(),wafer())-1];
double mips(unsigned dBx=0) const {
return charge(dBx)*Constants::fMipPerFc[Geometry::numberOfSublayers(layer(),wafer())-1];
}
double transverseMips() const {
return mips()*Geometry::point(*this).sinTheta();
double transverseMips(unsigned dBx=0) const {
return mips(dBx)*Geometry::point(*this).sinTheta();
}
double energy() const {
return mips()*Calibration::layerWeight(layer());
double energy(unsigned dBx=0) const {
return mips(dBx)*Calibration::layerWeight(layer());
}
double transverseEnergy() const {
return energy()*Geometry::point(*this).sinTheta();
double transverseEnergy(unsigned dBx=0) const {
return energy(dBx)*Geometry::point(*this).sinTheta();
}
// Access to input charge (pre-digitisation) information
......@@ -204,18 +205,22 @@ public:
// Access to SimHit information
unsigned numberOfSimHits() const {
return (fSimHitDepositSignal+fSimHitDepositPileup[fBxNumber%kNumberOfBx])>0.0?1:0;
return (fSimHitDepositSignal+fSimHitDepositPileup[/*AdcReadingMultiple::*/currentBx()])>0.0?1:0;
}
float simDeposit(bool signal=false) const {
//double w[kNumberOfBx]={0.199,0.406,0.736,1.000};
//double w[kNumberOfBx]={0.0,0.0,0.0,1.000};
//double w[AdcReadingMultiple::kNumberOfBx]={0.199,0.406,0.736,1.000};
//double w[AdcReadingMultiple::kNumberOfBx]={0.0,0.0,0.0,1.000};
double d(fSimHitDepositSignal);
if(!signal) {
for(unsigned bx(0);bx<kNumberOfBx;bx++) {
//unsigned bx(kNumberOfBx) {
d+=fWeightBx[bx]*fSimHitDepositPileup[(fBxNumber+bx+1)%kNumberOfBx];
/*
for(unsigned bx(0);bx<AdcReadingMultiple::kNumberOfBx;bx++) {
d+=fWeightBx[bx]*fSimHitDepositPileup[(AdcReadingMultiple::fBxNumber+bx+1)%AdcReadingMultiple::kNumberOfBx];
}
*/
for(unsigned dBx(0);dBx<AdcReadingMultiple::kNumberOfBx;dBx++) {
d+=fWeightBx[AdcReadingMultiple::kNumberOfBx-dBx-1]*fSimHitDepositPileup[previousBx(dBx)];
}
}
return d;
......@@ -333,7 +338,7 @@ public:
//if(!RecDetId::read(fin)) return false;
if(!DigHitDetId::read(fin)) return false;
if(!AdcReading::read(fin)) return false;
if(!AdcReadingMultiple::read(fin)) return false;
fin >> std::dec;
fin >> fGaussian;
......@@ -348,7 +353,7 @@ public:
if(!DigHitDetId::write(fout)) return false;
fout << " ";
if(!AdcReading::write(fout)) return false;
if(!AdcReadingMultiple::write(fout)) return false;
fout << std::dec;
fout << " " << fGaussian;
......@@ -361,7 +366,7 @@ public:
std::cout << "DigHitMultiple::";
//RecDetId::print();
DigHitDetId::print();
std::cout << " ";AdcReading::print();std::cout << std::endl;
std::cout << " ";AdcReadingMultiple::print();std::cout << std::endl;
std::cout << " Deposit " << deposit() << std::endl;
std::cout << " Charge " << charge() << std::endl;
std::cout << " MIPs " << mips() << std::endl;
......@@ -390,22 +395,29 @@ public:
// CMSSW values
if(fPeakingTime<0.0 || fRcPower==0) {
/*
fWeightBx[0]=0.0;
fWeightBx[1]=0.0;
fWeightBx[2]=0.0;
fWeightBx[3]=0.004;
fWeightBx[4]=0.200;
fWeightBx[5]=1.000;
*/
fWeightBx[0]=0.0;
fWeightBx[1]=0.004;
fWeightBx[2]=0.200;
fWeightBx[3]=1.000;
return;
}
// CR-RC shaper
for(unsigned bx(0);bx<kNumberOfBx;bx++) {
for(unsigned bx(0);bx<AdcReadingMultiple::kNumberOfBx;bx++) {
if(fPeakingTime>0.0) {
unsigned dbx(kNumberOfBx-bx-1);
unsigned dbx(AdcReadingMultiple::kNumberOfBx-bx-1);
fWeightBx[bx]=(1.0+dbx/fPeakingTime)*exp(-(dbx/fPeakingTime));
} else {
fWeightBx[bx]=(bx==(kNumberOfBx-1)?1.0:0.0);
fWeightBx[bx]=(bx==(AdcReadingMultiple::kNumberOfBx-1)?1.0:0.0);
}
fWeightBx[bx]=pow(fWeightBx[bx],fRcPower);
......@@ -424,36 +436,36 @@ public:
static void printStatics() {
std::cout << "DigHitMultiple::printStatics()" << std::endl
<< " Number of BX saved = " << kNumberOfBx << std::endl
<< " Number of BX saved = " << AdcReadingMultiple::kNumberOfBx << std::endl
<< " Peaking time in units of BX = " << fPeakingTime
<< ", in units of ns = " << 25.0*fPeakingTime << std::endl
<< " CR-RC^N with N = " << fRcPower << std::endl
<< " Weights per BX =";
for(unsigned bx(0);bx<kNumberOfBx;bx++) {
for(unsigned bx(0);bx<AdcReadingMultiple::kNumberOfBx;bx++) {
std::cout << " " << fWeightBx[bx];
}
std::cout << std::endl << std::endl;
}
static unsigned fBxNumber;
//static unsigned fBxNumber;
private:
static double fWeightBx[kNumberOfBx];
static double fWeightBx[AdcReadingMultiple::kNumberOfBx];
static double fPeakingTime;
static unsigned fRcPower;
float fGaussian;
float fSimHitDepositSignal;
float fSimHitDepositPileup[kNumberOfBx];
float fSimHitDepositPileup[AdcReadingMultiple::kNumberOfBx];
};
unsigned DigHitMultiple::fBxNumber=0;
//unsigned DigHitMultiple::fBxNumber=0;
double DigHitMultiple::fPeakingTime=1.0; // 25 ns in units of BX
unsigned DigHitMultiple::fRcPower=1;
double DigHitMultiple::fWeightBx[kNumberOfBx]={
(1+5.0/fPeakingTime)*exp(-5.0/fPeakingTime),
(1+4.0/fPeakingTime)*exp(-4.0/fPeakingTime),
double DigHitMultiple::fWeightBx[AdcReadingMultiple::kNumberOfBx]={
//(1+5.0/fPeakingTime)*exp(-5.0/fPeakingTime),
//(1+4.0/fPeakingTime)*exp(-4.0/fPeakingTime),
(1+3.0/fPeakingTime)*exp(-3.0/fPeakingTime),
(1+2.0/fPeakingTime)*exp(-2.0/fPeakingTime),
(1+1.0/fPeakingTime)*exp(-1.0/fPeakingTime),
......
......@@ -367,7 +367,7 @@ public:
private:
//double fGaussian;
float fGaussian;
std::vector<SimHit> *vSimHit;
//std::vector<SimHit> *vSimHit;
float fSimHitDeposit[2];
};
......
......@@ -91,6 +91,11 @@ public:
// Access SimHits
unsigned minBiasEvents() const {
//return fMinBiasEvents;
return vPuGenEvent.size();
}
std::vector<SimHit>& simHits(unsigned e, unsigned l, unsigned w, unsigned c) {
BACKTRACE_AND_ASSERT(e<Geometry::kNumberOfEndcaps);
BACKTRACE_AND_ASSERT(l<Geometry::kNumberOfLayers);
......@@ -778,6 +783,7 @@ public:
for(unsigned c(0);c<Geometry::numberOfSensorCellsPerWafer(l,w);c++) {
if(!vDigHit[e][l][w][c].invalid()) {
vDigHit[e][l][w][c].daqSelected(vDigHit[e][l][w][c].mips()>=fDigHitDaqSelectorMipsCut[sl]);
vDigHit[e][l][w][c].pbxSelected(vDigHit[e][l][w][c].mips(1)>=fDigHitPbxSelectorMipsCut[sl]);
vDigHit[e][l][w][c].tpgSelected(vDigHit[e][l][w][c].sigmas()>=fDigHitTpgSelectorSigmasCut);
} else {
vDigHit[e][l][w][c].daqSelected(false);
......@@ -1371,6 +1377,9 @@ static void printStatics() {
<< " fDigHitDaqSelectorMipsCut 100mu = " << fDigHitDaqSelectorMipsCut[0] << std::endl
<< " fDigHitDaqSelectorMipsCut 200mu = " << fDigHitDaqSelectorMipsCut[1] << std::endl
<< " fDigHitDaqSelectorMipsCut 300mu = " << fDigHitDaqSelectorMipsCut[2] << std::endl
<< " fDigHitPbxSelectorMipsCut 100mu = " << fDigHitPbxSelectorMipsCut[0] << std::endl
<< " fDigHitPbxSelectorMipsCut 200mu = " << fDigHitPbxSelectorMipsCut[1] << std::endl
<< " fDigHitPbxSelectorMipsCut 300mu = " << fDigHitPbxSelectorMipsCut[2] << std::endl
<< " fDigHitTpgSelectorSigmasCut = " << fDigHitTpgSelectorSigmasCut << std::endl
<< " fTrgHitSelectorMethod = " << fTrgHitSelectorMethod << std::endl
<< " fTrgHitSelectorMipsCut = " << fTrgHitSelectorMipsCut << std::endl
......@@ -1381,6 +1390,7 @@ static void printStatics() {
static double fDigHitDaqSelectorMipsCut[3];
static double fDigHitPbxSelectorMipsCut[3];
static double fDigHitTpgSelectorSigmasCut;
static unsigned fTrgHitSelectorMethod;
......@@ -1455,6 +1465,7 @@ TH2F *hXY3d;
//double Event::fDigHitDaqSelectorMipsCut[3]={0.5376,0.5230,0.5196}; // JBS: from CMSSW
double Event::fDigHitDaqSelectorMipsCut[3]={0.5,0.5,0.5};
double Event::fDigHitPbxSelectorMipsCut[3]={2.5,2.5,2.5};
double Event::fDigHitTpgSelectorSigmasCut=3.0;
unsigned Event::fTrgHitSelectorMethod=2;
......
......@@ -11,8 +11,7 @@
class GeometryRawV1toAnaV10 {
public:
GeometryRawV1toAnaV10() {
fGcm.read("GeometryCellMapRawV1toAnaV10.txt");
//fGcm.read("GeometryCellMapRawV1toAnaV10_FH.txt");
assert(fGcm.read("GeometryCellMapRawV1toAnaV10.txt"));
}
~GeometryRawV1toAnaV10() {
......
......@@ -66,7 +66,7 @@ public:
f >> d;
if(d!=fFileEvent+1) {
std::cout << "!!! d = " << d << ", fFileEvent = " << fFileEvent << std::endl;
//std::cout << "!!! d = " << d << ", fFileEvent = " << fFileEvent << std::endl;
//return false;
}
......
......@@ -21,6 +21,12 @@ public:
virtual void normalise() = 0;
static void incrementNormalisationAll(double n=1.0) {
for(unsigned i(0);i<vTHN.size();i++) {
vTHN[i]->incrementNormalisation(n);
}
}
static void normaliseAll() {
for(unsigned i(0);i<vTHN.size();i++) {
vTHN[i]->normalise();
......
......@@ -2,6 +2,10 @@
g++ -Wall -std=c++0x -I. code/AdcReadingTest.cc -o AdcReadingTest.exe `root-config --libs --cflags`
*/
//#include "RConfigure.h"
//#define R__HAS_STD_STRING_VIEW
#include "Hack.hh"
//#include <unistd.h>
#include <signal.h>
......@@ -15,7 +19,7 @@
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
//#include "TMath.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TFileHandler.hh"
......
//#include <unistd.h>
#include "Hack.hh"
#include <signal.h>
#include <cassert>
......
......@@ -2,7 +2,8 @@
g++ -Wall -std=c++0x -I. code/CalibrationTest.cc -o CalibrationTest.exe `root-config --libs --cflags`
*/
//#include <unistd.h>
#include "Hack.hh"
#include <signal.h>
#include <cassert>
......
#include "Hack.hh"
#include <iostream>
#include <sstream>
#include <vector>
......
#include "Hack.hh"
#include <sstream>
#include <fstream>
#include <cassert>
......@@ -51,6 +53,7 @@ int main(int argc, const char **argv) {
TFileHandler tfh("DataCompression");
TH1D *hRaw=new TH1D("Raw",";ADC value;Number",2048,0.0,2048.0);
TH1D *hRawInt=new TH1D("RawInt",";ADC value;Integral",2048,0.0,2048.0);
TH1D *hRawData=new TH1D("RawData",";ADC value;Bit-weighted number",2048,0.0,2048.0);
TH1D *hOne=new TH1D("One",";Division;Average bits",10,0.0,10.0);
......@@ -75,7 +78,7 @@ int main(int argc, const char **argv) {
assert(fin);
unsigned n,nBins(0);
double d[2048],dTotal(0.0);;
double d[2048],dTotal(0.0);
fin >> nBins;
assert(nBins==2048);
......@@ -93,6 +96,12 @@ int main(int argc, const char **argv) {
std::cout << "Total of data = " << dTotal << std::endl;
for(unsigned i(0);i<2048;i++) {
for(unsigned j(i);j<2048;j++) {
hRawInt->Fill(j+0.5,d[i]/dTotal);
}
}
double aMin(100.0);
for(unsigned j(0);j<11;j++) {
unsigned b(1<<j);
......@@ -390,5 +399,9 @@ int main(int argc, const char **argv) {
hHuffData->SetBinContent(i+1,bAdd*d[i]);
}
std::cout << "Huffman average = " << bSum/dTotal << std::endl;
std::cout << "RawInt[33] = " << hRawInt->GetBinContent(34) << std::endl;
return 0;
}
......@@ -167,6 +167,12 @@ int main(int argc, const char **argv) {
ec.daqSelected(false);
assert(!ec.daqSelected());
assert(!ec.pbxSelected());
ec.pbxSelected(true);
assert(ec.pbxSelected());
ec.pbxSelected(false);
assert(!ec.pbxSelected());
assert(!ec.tpgSelected());
ec.tpgSelected(true);
assert(ec.tpgSelected());
......@@ -178,6 +184,20 @@ int main(int argc, const char **argv) {
assert(ec.invalid());
ec.invalid(false);
assert(!ec.invalid());
ec.setSensorCell(c);
assert(ec.endcap()==e);
assert(ec.layer()==l);
assert(ec.wafer()==w);
assert(ec.sensorCell()==c);
assert(ec.reserved()==0);
ec.setSensorCell(511-c);
assert(ec.endcap()==e);
assert(ec.layer()==l);
assert(ec.wafer()==w);
assert(ec.sensorCell()==511-c);
assert(ec.reserved()==0);
}
}
}
......
#include "Hack.hh"
#include <iostream>
#include <sstream>
......
#include "Hack.hh"
#include <iostream>
#include <sstream>
......