Commit 42dce392 authored by Luca Mastrolorenzo's avatar Luca Mastrolorenzo
Browse files

addition of the classes allowing the calibration procedure be performed all at...

addition of the classes allowing the calibration procedure be performed all at once + trgCalib.cpp main dedicated only to run these classes + addition of analysis histos in CoefficientDetClu3D.hh
parent 64510062
#ifndef CalibResults_HH
#define CalibResults_HH
#include <iostream>
#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <algorithm>
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "Random.hh"
#include "TE1F.hh"
#include "TE2F.hh"
#include "AnalysisBase.hh"
#include "Event.hh"
#include "TMath.h"
class CalibResults : public AnalysisBase {
public:
enum {
kCutBins=50
};
CalibResults(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("CalibResults",sRoot,sOut) {
//AnalysisBase(sName), fTwoPi(2.0*acos(-1.0)) {
h_events = new TH1F("h_events", "h_events",2,0,2);
h_Matrix_T = new TH2F("h_Matrix_T","d_e_i * d_e_j", 53, 0, 53, 53, 0, 53);
h_N2 = new TH2F("h_N2","N2", 53, 0, 53, 53, 0, 53);
h_Vector_T = new TH1F("h_Vector_T","vector_T", 53, 0, 53);
h_N = new TH1F("h_N","N", 53, 0, 53);
h_C3dPtResponse = new TH1F("h_C3dPtResponse","h_C3dPtResponse",100,0,2);
h_C3dPtResponse2D = new TH2F("h_C3dPtResponse2D","h_C3dPtResponse2D",600,0,300,600,0,300);
double deltaEta(double eta1, double eta2);
double deltaPhi(double phi1, double phi2);
double dist2D(double x1, double y1, double x2, double y2);
for(unsigned e(0); e<Geometry::kNumberOfEndcaps; e++)
for(unsigned l(0); l<52; l++)
for(unsigned w(0); w<700; w++)
MatrixEvent[e][l][w] = -1;
calibFile.open((sOut+"/calibrationCoefficient.txt"));
}
virtual ~CalibResults() {
}
virtual bool event(Event &event) {
double Eta_min = 1.7;
double Eta_max = 2.7;
unsigned minC2dInC3d = 3;
TString ParticleType = "photon";
bool debug_ = false;
if(fEventNumber <= 1){
for(int i=0; i<40; i++){
if(debug_) std::cout<< "layer weight = " << Calibration::layerWeight(i) << std::endl;
}
}
if(true) std::cout << fName << " Event " << fEventNumber << std::endl;
h_events->Fill(1);
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
double t_deposits_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
double t_weightDep_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
for(int i(0); i< Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; i++){
t_deposits_T[i]=0.;
t_weightDep_T[i]=0.;
}
std::vector<Particle> vSig;
//(event.photons(e));
if(ParticleType=="pion"){
std::cout << "Using Pions! \n";
vSig=event.pions(e);
}else if(ParticleType=="photon"){
std::cout << "Using Photons! \n";
vSig=event.photons(e);
}else if(ParticleType=="electron"){
std::cout << "Using Electrons! \n";
vSig=event.electrons(e);
}
//select event with signal particle in the geometrical and kinematical acceptance
int goodLayer = 0;
for(unsigned l(0); l<Geometry::kNumberOfLayers; l++) {
std::vector<Particle> vSigProjTest;
for(unsigned i(0);i<vSig.size();i++){
vSigProjTest.push_back(vSig[i].project(Geometry::layerZ(e,l)));
}
for(unsigned i(0);i<vSigProjTest.size();i++){
if( fabs( vSigProjTest[i].position().eta() ) > Eta_min && fabs( vSigProjTest[i].position().eta() ) < Eta_max ) goodLayer = goodLayer + 1 ;
}
}
//consider only the true particles that are goemetrically allowed to deposits in at least 40 layer (EE+FH?)
if(goodLayer<40) continue;
/** LOOP 3d clusters **/
std::vector<TrgC3d> vTrgC3d( event.trgC3ds(e) );
std::sort( vTrgC3d.begin(), vTrgC3d.end(), eneTranOrderC3d );
if(vTrgC3d.size()<=0) continue;
//select the best 3D-cluster per given endcap and assign the correspondent vector of 2D-cluster
std::vector<TrgC2d*> &vTrgC2d_fromMaxC3d( vTrgC3d[0].trgC2ds() );
//select only those maximum clusters that have more then 3 C2d inside
if(vTrgC2d_fromMaxC3d.size()<=minC2dInC3d) continue;
//loop over the C2d clusters collection of the max-C3d
for( unsigned iC2d = 0; iC2d < vTrgC2d_fromMaxC3d.size(); iC2d++ ) {
if(vTrgC2d_fromMaxC3d[iC2d]->transverseMips()>0 ){
t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = vTrgC2d_fromMaxC3d[iC2d]->transverseMips();
t_weightDep_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = vTrgC2d_fromMaxC3d[iC2d]->transverseMips() * vSig[0].FourP().Et();
if(debug_) std::cout<< "--------->>>>>>>>>> " << t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] << std::endl;
}else{
t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = 0;
t_weightDep_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = 0;
}
//assign 0 to proper layers without energy deposits
for(int ll=0; ll<52; ll++){
if( t_deposits_T[ll]<0.0001 || !(t_deposits_T[ll]==0 || t_deposits_T[ll]!=0)){
t_deposits_T[ll] = 0.;
t_weightDep_T[ll] = 0.;
if( isnan( t_weightDep_T[ll] )) if(debug_) std::cout<< "IS NAN " << std::endl;
}
}
}// end loop over C2d from MaxC3d
t_deposits_T[52] = 1.;
t_weightDep_T[52] = 1. * vSig[0].FourP().Et();
for(int hh=0; hh<52; hh++){
if( isnan( t_deposits_T[hh] )){
t_deposits_T[hh] = 0.;
if(debug_) std::cout<< "IS NAN " << std::endl;
}
if(debug_) std::cout<< "cross check deposit vec = "<< hh << ", " << t_deposits_T[hh] << std::endl;
}
double E_CluCalibMIPt_26trgLayers = 0.;
std::string content;
int i=0;
while(calibFile >> content) {
cout << "READING CALIB FILE : " << content << std::endl;
E_CluCalibMIPt_26trgLayers += t_deposits_T[i] * atof(content.c_str());
i++;
}
std::cout << "CALIBRATED ENERGY = "<< E_CluCalibMIPt_26trgLayers << " GeV" << std::endl;
if(vTrgC3d.size()>0){
//if(debug_) std::cout<< "energy particle mc << " << vSig[0].FourP().Et() << " energy after " << E_CluCalibMIPt_26trgLayers<< std::endl;
h_C3dPtResponse->Fill(E_CluCalibMIPt_26trgLayers/vSig[0].FourP().Et());
h_C3dPtResponse2D->Fill(vSig[0].FourP().Et(), E_CluCalibMIPt_26trgLayers);
}
for(int i=0; i<Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; i++){
h_N->Fill(i,1);
if( t_weightDep_T[i]>0 ){
if(debug_) std::cout<< " =============== fill transverse vector = " << i << " , " << t_weightDep_T[i] << std::endl;
h_Vector_T->Fill(i,t_weightDep_T[i]);
}
for(int j=0; j<Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; j++){
h_N2->Fill(i, j, 1);
if(t_deposits_T[i]>0 && t_deposits_T[j]>0){
h_Matrix_T->Fill( i, j, t_deposits_T[i] * t_deposits_T[j] );
}
}
}
}//end loop ever endcap
return true;
}
private:
// const double fTwoPi;
TH1F *h_events;
TH2F *h_Matrix_T;
TH1F *h_Vector_T;
TH2F *h_N2;
TH1F *h_N;
TH1F *h_C3dPtResponse;
TH2F *h_C3dPtResponse2D;
int MatrixEvent[2][52][700];
double deltaPhi( double phi1, double phi2) {
double dPhi(phi1-phi2);
double pi(acos(-1.0));
if (dPhi<=-pi) dPhi+=2.0*pi;
else if(dPhi> pi) dPhi-=2.0*pi;
return dPhi;
}
double deltaEta(double eta1, double eta2){
double dEta = (eta1-eta2);
return dEta;
}
double deltaR(double eta1, double eta2, double phi1, double phi2) {
double dEta = deltaEta(eta1, eta2);
double dPhi = deltaPhi(phi1, phi2);
return sqrt(dEta*dEta+dPhi*dPhi);
}
double dist2D(double x1, double y1, double x2, double y2){
double dx = ( x1 - x2 );
double dy = ( y1 - y2 );
return sqrt( dx*dx + dy*dy );
}
ifstream calibFile;
};
#endif
This diff is collapsed.
#ifndef FillMatrixForCalib_HH
#define FillMatrixForCalib_HH
#include <iostream>
#include <cassert>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include <algorithm>
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "Random.hh"
#include "TE1F.hh"
#include "TE2F.hh"
#include "AnalysisBase.hh"
#include "Event.hh"
#include "TMath.h"
class FillMatrixForCalib : public AnalysisBase {
public:
enum {
kCutBins=50
};
FillMatrixForCalib(const std::string &sRoot="", const std::string &sOut="")
: AnalysisBase("FillMatrixForCalib",sRoot,sOut) {
//AnalysisBase(sName), fTwoPi(2.0*acos(-1.0)) {
h_events = new TH1F("h_events", "h_events",2,0,2);
h_Matrix_T = new TH2F("h_Matrix_T","d_e_i * d_e_j", 53, 0, 53, 53, 0, 53);
h_N2 = new TH2F("h_N2","N2", 53, 0, 53, 53, 0, 53);
h_Vector_T = new TH1F("h_Vector_T","vector_T", 53, 0, 53);
h_N = new TH1F("h_N","N", 53, 0, 53);
double deltaEta(double eta1, double eta2);
double deltaPhi(double phi1, double phi2);
double dist2D(double x1, double y1, double x2, double y2);
for(unsigned e(0); e<Geometry::kNumberOfEndcaps; e++)
for(unsigned l(0); l<52; l++)
for(unsigned w(0); w<700; w++)
MatrixEvent[e][l][w] = -1;
}
virtual ~FillMatrixForCalib() {
}
virtual bool event(Event &event) {
double Eta_min = 1.7;
double Eta_max = 2.7;
unsigned minC2dInC3d = 3;
TString ParticleType = "photon";
bool debug_ = false;
if(fEventNumber <= 1){
for(int i=0; i<40; i++){
if(debug_) std::cout<< "layer weight = " << Calibration::layerWeight(i) << std::endl;
}
}
if(true) std::cout << fName << " Event " << fEventNumber << std::endl;
h_events->Fill(1);
for(unsigned e(0);e<Geometry::kNumberOfEndcaps;e++) {
double t_deposits_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
double t_weightDep_T[Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1];
for(int i(0); i< Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; i++){
t_deposits_T[i]=0.;
t_weightDep_T[i]=0.;
}
std::vector<Particle> vSig;
//(event.photons(e));
if(ParticleType=="pion"){
std::cout << "Using Pions! \n";
vSig=event.pions(e);
}else if(ParticleType=="photon"){
std::cout << "Using Photons! \n";
vSig=event.photons(e);
}else if(ParticleType=="electron"){
std::cout << "Using Electrons! \n";
vSig=event.electrons(e);
}
//select event with signal particle in the geometrical and kinematical acceptance
int goodLayer = 0;
for(unsigned l(0); l<Geometry::kNumberOfLayers; l++) {
std::vector<Particle> vSigProjTest;
for(unsigned i(0);i<vSig.size();i++){
vSigProjTest.push_back(vSig[i].project(Geometry::layerZ(e,l)));
}
for(unsigned i(0);i<vSigProjTest.size();i++){
if( fabs( vSigProjTest[i].position().eta() ) > Eta_min && fabs( vSigProjTest[i].position().eta() ) < Eta_max ) goodLayer = goodLayer + 1 ;
}
}
//consider only the true particles that are goemetrically allowed to deposits in at least 40 layer (EE+FH?)
if(goodLayer<40) continue;
/** LOOP 3d clusters **/
std::vector<TrgC3d> vTrgC3d( event.trgC3ds(e) );
std::sort( vTrgC3d.begin(), vTrgC3d.end(), eneTranOrderC3d );
if(vTrgC3d.size()<=0) continue;
//select the best 3D-cluster per given endcap and assign the correspondent vector of 2D-cluster
std::vector<TrgC2d*> &vTrgC2d_fromMaxC3d( vTrgC3d[0].trgC2ds() );
//select only those maximum clusters that have more then 3 C2d inside
if(vTrgC2d_fromMaxC3d.size()<=minC2dInC3d) continue;
//loop over the C2d clusters collection of the max-C3d
for( unsigned iC2d = 0; iC2d < vTrgC2d_fromMaxC3d.size(); iC2d++ ) {
if(vTrgC2d_fromMaxC3d[iC2d]->transverseMips()>0 ){
t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = vTrgC2d_fromMaxC3d[iC2d]->transverseMips();
t_weightDep_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = vTrgC2d_fromMaxC3d[iC2d]->transverseMips() * vSig[0].FourP().Et();
if(debug_) std::cout<< "--------->>>>>>>>>> " << t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] << std::endl;
}else{
t_deposits_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = 0;
t_weightDep_T[vTrgC2d_fromMaxC3d[iC2d]->layer()] = 0;
}
//assign 0 to proper layers without energy deposits
for(int ll=0; ll<52; ll++){
if( t_deposits_T[ll]<0.0001 || !(t_deposits_T[ll]==0 || t_deposits_T[ll]!=0)){
t_deposits_T[ll] = 0.;
t_weightDep_T[ll] = 0.;
if( isnan( t_weightDep_T[ll] )) if(debug_) std::cout<< "IS NAN " << std::endl;
}
}
}// end loop over C2d from MaxC3d
t_deposits_T[52] = 1.;
t_weightDep_T[52] = 1. * vSig[0].FourP().Et();
for(int hh=0; hh<52; hh++){
if( isnan( t_deposits_T[hh] )){
t_deposits_T[hh] = 0.;
if(debug_) std::cout<< "IS NAN " << std::endl;
}
if(debug_) std::cout<< "cross check deposit vec = "<< hh << ", " << t_deposits_T[hh] << std::endl;
}
for(int i=0; i<Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; i++){
h_N->Fill(i,1);
if( t_weightDep_T[i]>0 ){
if(debug_) std::cout<< " =============== fill transverse vector = " << i << " , " << t_weightDep_T[i] << std::endl;
h_Vector_T->Fill(i,t_weightDep_T[i]);
}
for(int j=0; j<Geometry::kNumberOfLayers + Geometry::kNumberOfBHLayers + 1; j++){
h_N2->Fill(i, j, 1);
if(t_deposits_T[i]>0 && t_deposits_T[j]>0){
h_Matrix_T->Fill( i, j, t_deposits_T[i] * t_deposits_T[j] );
}
}
}
}//end loop ever endcap
return true;
}
private:
// const double fTwoPi;
TH2F *h_Matrix_T;
TH1F *h_Vector_T;
TH2F *h_N2;
TH1F *h_N;
TH1F *h_events;
int MatrixEvent[2][52][700];
double deltaPhi( double phi1, double phi2) {
double dPhi(phi1-phi2);
double pi(acos(-1.0));
if (dPhi<=-pi) dPhi+=2.0*pi;
else if(dPhi> pi) dPhi-=2.0*pi;
return dPhi;
}
double deltaEta(double eta1, double eta2){
double dEta = (eta1-eta2);
return dEta;
}
double deltaR(double eta1, double eta2, double phi1, double phi2) {
double dEta = deltaEta(eta1, eta2);
double dPhi = deltaPhi(phi1, phi2);
return sqrt(dEta*dEta+dPhi*dPhi);
}
double dist2D(double x1, double y1, double x2, double y2){
double dx = ( x1 - x2 );
double dy = ( y1 - y2 );
return sqrt( dx*dx + dy*dy );
}
};
#endif
#ifndef MatrixInversion_HH
#define MatrixInversion_HH
#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <cstdlib>
#include <TString.h>
#include <algorithm>
#include <map>
#include <sstream>
#include "TMath.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TString.h"
#include "TFile.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TClonesArray.h"
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include "FillMatrixForCalib.hh"
class MatrixInversion{
public:
MatrixInversion(){}
~MatrixInversion(){}
// std::map<int,double> CalibCoefficients( FillMatrixForCalib *analysis ){
std::map<int,double> CalibCoefficients( TString fileName ){
const int Nbin = 26;
TFile*fIn = new TFile(fileName, "READ");
fIn->ls();
//TH2F* h_de_i_de_j = analysis->h_Matrix_T;
//TH2F* h_N2 = analysis->h_N2;
//TH1F* h_de_j_Te = analysis->h_Vector_T;
//TH1F* h_N = analysis->h_N;
TH2F* h_de_i_de_j = (TH2F*)fIn->Get("h_Matrix_T");
TH2F* h_N2 = (TH2F*)fIn->Get("h_N2");
TH1F* h_de_j_Te = (TH1F*)fIn->Get("h_Vector_T");
TH1F* h_N = (TH1F*)fIn->Get("h_N");
std::cout << "Number of Bin "<< Nbin << std::endl;
TH2F* h_M = new TH2F("h_M","h_M",Nbin,0,Nbin,Nbin,0,Nbin);
TH2F* h_Minv = new TH2F("h_Minv","h_Minv",Nbin,0,Nbin,Nbin,0,Nbin);
TH2F* h_I = new TH2F("h_I","h_I",Nbin,0,Nbin,Nbin,0,Nbin);
TH1F* h_a = new TH1F("h_a","h_a",40,0,40);
TH1F* h_v = new TH1F("h_v","h_v",Nbin,0,Nbin);
std::cout << "in calib 1" << std::endl;
TMatrixDSym m(Nbin);
TMatrixDSym m_original(Nbin);
TVectorD V(Nbin);
std::map<int, double> CoefficientsMap;
int k[Nbin] = {1,3,5,7,9,11,13,15,17,19,21,23,25,27,28,29,30,31,32,33,34,35,36,37,38,39};
int r[Nbin] = {1,3,5,7,9,11,13,15,17,19,21,23,25,27,28,29,30,31,32,33,34,35,36,37,38,39};
for(int i=0; i<Nbin; i++){
//std::cout << "trigger layer = " << k[i] << " index: " << i << std::endl;
V(i)= h_de_j_Te->GetBinContent(k[i]+1) / h_N->GetBinContent(k[i]+1) ;
h_v->SetBinContent(i+1, V(i) );
for(int j=0; j<Nbin; j++){
//std::cout << "trigger layer = " << r[j] << " index: " << j << std::endl;
std::cout << "j,r = "<< k[i] << ", " << r[j] << std::endl;
m(i,j) = h_de_i_de_j->GetBinContent(k[i]+1,r[j]+1)/h_N2->GetBinContent(k[i]+1,r[j]+1);
m_original(i,j) = h_de_i_de_j->GetBinContent(k[i]+1,r[j]+1)/h_N2->GetBinContent(k[i]+1,r[j]+1);
h_M->SetBinContent(i+1, j+1, m(i,j));
}
}
TMatrixD M = m.Invert();
TMatrixD I = m_original*M;
for(int i=0; i<Nbin; i++){
std::cout << "V_" << i << " " << V(i) << std::endl;
}
for(int i=0; i<Nbin; i++){
for(int j=0; j<Nbin; j++){
h_Minv->SetBinContent(i+1, j+1, M(i,j));
//if(i==j) std::cout << I(i,j)<< std::endl;
h_I->SetBinContent(i+1, j+1, I(i,j));
}
}
TVectorD a(Nbin);
a = M * V;
std::cout << "from vector \n";
for(int i=0; i<Nbin; i++){
std::cout << i << " " << a(i) << std::endl;
}
std::cout << "Map content" << std::endl;
CoefficientsMap[0] = 0.;
for(int i=1; i<40; i++){
if(i<=27){
if(i%2){
CoefficientsMap[i] = a((int)((i-1)/2));
std::cout << i <<" " << (int)((i-1)/2) << " " << a((int)((i-1)/2))<< " " << std::endl ;
h_a->SetBinContent(i+1, a((int)((i-1)/2)));
}else if(!(i%2)){
CoefficientsMap[i] = 0;
std::cout << i << " 0. " << std::endl;
h_a->SetBinContent(i+1, 0);
}
}else if(i>27){
std::cout << " coefficient in the FH \n";
CoefficientsMap[i] = a(i-14);
std::cout << i << " " << i-14 << " " << a(i-14) << std::endl ;
}
}
return CoefficientsMap;
}
};
#endif
...@@ -172,6 +172,19 @@ public: ...@@ -172,6 +172,19 @@ public:
return e; return e;
} }
double trgC2dMaxVsAllHitRatio() const {
double e(0.0);
double max(-1.0);
double frac(0.0);
for(unsigned i(0);i<vTrgHit.size();i++) {
e+=vTrgHit[i]->transverseEnergy();
if(vTrgHit[i]->transverseEnergy() > max){
max = vTrgHit[i]->transverseEnergy();
}
}
frac = max / e;
return frac;
}
// Access to SimHit information // Access to SimHit information
unsigned numberOfSimHits() const { unsigned numberOfSimHits() const {
......
...@@ -5,7 +5,7 @@ ...@@ -5,7 +5,7 @@
#include <cstdint> #include <cstdint>
#include <cassert> #include <cassert>
#include <vector>